Three-dimensional localization microscopy in live flowing cells

Capturing the dynamics of live cell populations with nanoscale resolution poses a significant challenge, primarily owing to the speed-resolution trade-off of existing microscopy techniques. Flow cytometry would offer sufficient throughput, but lacks subsample detail. Here we show that imaging flow cytometry, in which the point detectors of flow cytometry are replaced with a camera to record 2D images, is compatible with 3D localization microscopy through point-spread-function engineering, which encodes the depth of the emitter into the emission pattern captured by the camera. The extraction of 3D positions from sub-cellular objects of interest is achieved by calibrating the depth-dependent response of the imaging system using fluorescent beads mixed with the sample buffer. This approach enables 4D imaging of up to tens of thousands of objects per minute and can be applied to characterize chromatin dynamics and the uptake and spatial distribution of nanoparticles in live cancer cells. A high-throughput nanoscale 3D microscopy technique, combining imaging flow cytometry with a point-spread function, allows 3D localization measurements in live flowing cells.

F low cytometry enables the rapid characterization of large populations of cells and particles. Increasing reliability and decreasing costs have helped the technology emerge at the forefront of medical diagnostics, drug-efficacy screens, cell sorting, population characterization and more [1][2][3] . While fluorescence intensity, light scatter and absorbance measurements are now standard capabilities of commercially available instruments, subsample characterization is typically not possible when the signal is relayed onto point detectors, thus the interpretability of acquired data is constrained to relative measurements and signature recognition using standards 4 .
Imaging flow cytometry (IFC) replaces point detectors with a camera, enabling widefield microscopy at high-throughput [5][6][7] . Recently, instruments have used multiple spectral channels 8,9 and other microscope modalities; for example, phase imaging 10 and optical sectioning 11,12 . Generally, because the operating principles of IFC and widefield microscopy are the same, tools developed for one can be applied to the other.
One widefield technique used to attain nanoscale information is localization microscopy 13 . The central premise is that spreading the emission light of a point source over several camera pixels enables the emitter's position to be extracted with a precision far below the pixel size. This approach has two main requirements: (1) a low emitter density in each image such that individual emitters can be recognized and (2) knowledge of the optical system's pointspread function (PSF), that is the expected distribution of light on the detector that emanates from a point source.
The low-emitter-density requirement can be satisfied, to some extent, by flow cytometry, where the microfluidics are designed to pass objects through the imaging region one at a time. Knowledge of the PSF is used to extract the underlying parameters from an image, namely an emitter's position. Note that dense yet small volumes of emitters, such as small beads or punctate spots in cells, can be treated collectively as single point sources so long as they occupy a subdiffraction volume.
In two-dimensional (2D) localization microscopy, the in-focus PSF is often approximated as a symmetric, 2D-Gaussian function, which is used to fit the measured image. This model can account for some defocus, which manifests as broadening (or blur) in the image. This approach is poorly suited for three-dimensional (3D) microscopy because there is a redundant widening in the PSF as the object is displaced in-front-of or behind the focal plane 14 . This ambiguity is resolved in iPALM and biplane microscopy by imaging on multiple detectors [15][16][17] . Alternatively, by adding an aberration into the imaging system by inserting an additional optical element, the PSF itself can be modified to remove the ambiguity. This technique is broadly defined as PSF-engineering microscopy and is advantageous in that it can be implemented in a standard, single-camera imaging system 18 .
A simple implementation of PSF engineering is achieved by adding a cylindrical lens into the imaging path, which encodes the depth into the directional stretching of the PSF shape 14,19 . The result can be fit with an asymmetric 2D-Gaussian function, then processed via a shape-to-depth calibration curve, which is normally generated by imaging a stationary object at various defocus levels 19,20 . More complicated PSFs, realized by adding a phase-shaping optical element in the microscope's back focal plane, can also be used to encode 3D spatial information 18 . In a related yet disparate application, the incorporation of diffractive optical elements for enhanced-depthof-field 21 imaging has been applied to IFC, modifying the PSF of the microscope to make it depth independent; 22 this effectively produces a 2D projection of the imaged volume.
Here, we demonstrate high-throughput, nanoscale, multicolour 3D localization by combining IFC with PSF engineering. Furthermore, we show this method can be implemented directly into existing instruments with a minor, removable addition. The technique is composed of four steps: (1) inserting an optical element (a phase mask or cylindrical lens) into the imaging path of an imaging flow cytometer; (2) capturing images of flowing objects, for example, cells, along with a calibration sample that follows a known or measurable depth-distribution function, that is, fluorescent beads; (3) generating a calibration curve by decoding the PSF response for the calibration sample and (4) applying the calibration curve to obtain the 3D positions of the sample images.

Three-dimensional localization microscopy in live flowing cells
We validate the approach in vitro using fluorescently labelled DNA-origami nanorulers, and in vivo by measuring chromosomal compaction states inside live yeast cells. Next, we demonstrate the application to nanomedicine by quantifying the uptake and spatial distributions of liposome nanoparticles in human T lymphocytes. Last, we show how PSF engineering can extend the useful depth range in IFC using the Tetrapod PSF 23 . Notably, the flow-based experiments presented here were based on a commercially available imaging flow cytometer, the Amnis ImageStreamX. The instrument's data output was processed using custom MATLAB scripts, thereby demonstrating the immediate and broad applicability of our approach to existing apparatus.

Device design and in vitro validation
An illustration of the device is shown in Fig. 1a. Loaded samples flow through a laser-illuminated imaging volume. Light from the sample is captured by an objective lens, relayed through a series of lenses, divided into separate spectral channels and directed onto a camera. The camera's readout is synchronized with the flow speed, that is operating in time-delay-integration mode (see Supplementary  Information). This synchronization circumvents motion blur, enabling a longer integration time; that is, higher signal. Inserting a cylindrical lens between two relay lenses (Methods) induces an astigmatism that stretches the PSF: objects closer to the objective relative to the focal plane appear vertically extended, but horizontally narrow and vice versa (Fig. 1b,c). The extent of the astigmatism contains information about the emitter's depth, but is not a direct measurement of the z position, thus a calibration step is needed.
Unlike a standard microscope, whose PSF can be characterized by z scanning over a stationary fiducial 19 , in a flow cell, it is not possible to scan over a point source. Our solution for attaining an accurate experimental-PSF calibration is to exploit the statistical properties of the flowing beads. In brief, since the probability distribution of object depths (z positions) in flow is well defined and measurable (Fig. 1d), a sufficiently large sample size will closely resemble the average distribution. A dataset of calibration images can therefore be recorded and algorithmically sorted by relative depth, for example by the extent of the astigmatism, producing a shape-distribution curve (Fig. 1e). This shape distribution is then parameterized by the depth percentile, where each position along the curve is assigned the statistically most likely z position (Fig. 1f). Increasing the number of calibration measurements decreases the position errors ( Supplementary Fig. 1). New images can be compared to the reference curve to find their depth irrespective of their probability distribution in the flow. In other words, after calibration, any point-like object can be localized in 3D.
While determining the 3D positions in a single colour channel can be used to measure the flow profile and characterize the PSF 24 , of more practical interest in imaging cytometry is measuring the relation between two or more objects, that is colocalization. For emitters separated by distances below the diffraction limit (~300 nm), this can be accomplished by capturing the emission of each object in a different image, for example by spectrally separating the objects. To assess our ability to use the six spectral channels collected by the ImageStreamX and characterize our method's performance, we employed 200 nm Tetraspeck fluorescent beads due to their broad emission spectrum (Methods). We collected a six-channel dataset of 50,000 beads in less than 3 min. Next, we divided the bead data into two independent groups: one for calibration and one for evaluating the performance of our approach; namely, the colocalization distance measured between spectral channels. Images of an example bead are shown in Fig. 2a, where the localized positions are shown in Fig. 2b,c. The interchannel residual distances are plotted as 2D histograms in Fig. 2d, and were used to compute the geometric 3D distances per channel pair. The cross-channel error was defined as the median 3D distance (Fig. 2e), and ranged from 25 to 53 nm for various channel pairs in optimized conditions (that is, cylindrical lens placement, laser power and flow stability, Supplementary Fig. 2). objects is captured by an objective lens, relayed with an extended optical system through a cylindrical lens (red box), split into spectral channels, then imaged on a camera. NA, numerical aperture. b,c, Objects at different depths in the flow cell (b) exhibit depth-dependent image shapes (c). Scale bar, 1 µm. d, The relative fraction of in-focus objects measured at various objective positions. e, The 2D distribution of astigmatic shapes from localized objects, where the density is shown as the number of emitters per 10 nm 2 in terms of the Gaussian shape parameters σ X and σ Y . f, The shape-to-depth calibration curve generated from d,e, where z = 0 µm corresponds to the centre of the object distribution.
To demonstrate 3D distance measurements between pairs of objects, we employed a DNA-origami nanoruler. The two ends were labelled with green (Atto488) and red (Atto647N) fluorophores, with a designed separation of 180 nm (ref. 25 ; Fig. 3a,b). The sample was diluted and mixed with fluorescent beads before imaging. A subset of fluorescent beads was used for calibration and the remainder for comparison to the nanorulers. The calculated interchannel distances in each axis (Fig. 3d) was used to compute the geometric 3D distance between the probes per object; that is, to calculate the object length. We found the results obtained in flow, 164 ± 94 nm (mean ± s.d., see Fig. 3e), were similar to that observed by conventional 2D microscopy where the nanorulers were fixed to a glass coverslip, 151 ± 66 (mean ± s.d., Supplementary Fig. 3).

chromosomal compaction in Saccharomyces cerevisiae
Populations of live cells often exhibit broad heterogeneity with variable responses to stimuli that necessitates high-throughput imaging 26 . For example, chromosomal repair after ultraviolet-induced DNA damage 27 , responses to changes in the chemical environment 28 and gene regulation 26,29 have all been shown to involve dynamic chromosomal reconfiguration in 3D.
To demonstrate the applicability of our approach to such measurements, we imaged fluorescently tagged DNA loci of live yeast cells to investigate a proposed mechanism for gene regulation: chromosomal regions containing suppressed genes are spatially compacted 30 . In brief, two arrays of operator sites were integrated into the chromosome flanking the Gal locus 28 , a region encoding several genes involved in galactose metabolism 31 . Fluorescently labelled operator-binding proteins attach to each operator array forming two fluorescent spots (Fig. 4a), where the interspot distance is affected by the growth conditions (Fig. 4b). To account for the cellular autofluorescence and unbound fluorescent proteins, our localization algorithm was modified to include an additional Gaussian-background term (Fig. 4c,d). Each cell image was then analysed in the two colour channels separately to attain absolute 3D positions (Fig. 4e,g), which were used to calculate the relative displacement in each dimension separately and in 3D (Fig. 4f,h,i). Similar to previously measurements 32,33 , the mean interloci distance  was observed to be dependent on the growth condition (mean ± s.d., 339 ± 174 nm in dextrose and 434 ± 193 nm in galactose). Crucially, while previous 3D datasets were recorded over hours at roughly one image per second 33 , the current method acquires hundreds of images per second, producing a large library of suitable cell images in only a few minutes (N = 8,996 for dextrose and 11,074 for galactose, recorded in <5 and <10 min, respectively). A similar change in the interloci distance was observable for cells that were fixed before imaging (Supplementary Fig. 4 and Methods). The 3D colocalization precision was estimated to be 160 ± 100 nm (mean ± s.d.), found by colocalizing the mCherry fluorophore signal captured in spectral channels 4 and 5.
A fundamental challenge in imaging heterogeneous and temporally dynamic populations is the necessity to gather large statistics over short time windows. This is where our high throughput approach provides a unique capability for 3D imaging. To investigate the dynamics of chromatin compaction, we applied our technique in a continuous measurement of the same system (Fig. 5). Cells were grown to log phase in either raffinose or galactose, corresponding to a closed or open Gal locus, respectively. Concentrated galactose or dextrose was mixed with the prepared cells and immediately loaded into the ImageStreamX for imaging (t = 0). The mean 3D distance is immediately affected by the switching of the medium (Fig. 5a), which corresponds to realtime changes in the distribution of chromosomal conformation states (Fig. 5b).

Nanoparticle penetration into human T lymphocytes
Liposome nanoparticles have received attention as targetable drug carriers, where the chemical composition can be optimized to deliver therapeutic payloads into cells, for example, small-molecule drugs, proteins and nucleic acids [34][35][36] . Of particular interest is how therapeutic nanoparticles are taken up by cells, including lymphocytes and other non-adherent cells. The measurement necessitates 3D imaging because one-dimensional (non-IFC) and 2D measurements only characterize colocalization of particles with cells, but cannot directly quantify the degree of penetration, that is, the intracellular coordinates on the single-cell level. In addition, nonadherent cells pose a challenge for 3D microscopy techniques that require samples to be stationary. Cells can be attached to a surface before imaging 37 , although this preprocessing step can affect viability and morphology, and may preclude further analysis. Flow-based techniques avoid this step, capturing data from suspended live cells as they move through the imaging volume.
To demonstrate how to use 3D IFC to characterize the uptake of liposomes into non-adherent cells, we exposed human T lymphocyte (Jurkat) cells to fluorescently labelled 150 nm liposome nanoparticles (Fig. 6a, see Methods). To localize the positions of the liposomes with respect to the cell, extending our method beyond colocalization of spots, we used a brightfield image to find the x and y centre of the cell and the 3D positions of two or more cell-surfaceadhered fluorescent beads to identify the cell's centre in z and radius. This was possible due to the spherical shape of these non-adherent cells. Cells detected to be aspherical, without two or more surface beads, or lacking liposomes were excluded from analysis. Example images are shown in Fig. 6b. Normalized to the cell size, the radial position of the liposomes, 0.80 ± 0.23 (mean ± s.d.) was significantly smaller than that of the fluorescent polystyrene beads, 0.99 ± 0.12, indicating liposome uptake into the cells (N cells = 323, P < 1 × 10 −6 , using a one-sided, two-sample t-test; Fig. 6c). Of significant interest for drug-delivery applications is the fraction of liposome that enter the nucleus; that is, the centre region of the cell 34   (Hoechst) and imaged in 2D ( Supplementary Fig. 5). The average diameter of the nucleus (7 µm) was 67 ± 9% of the cell (10.5 µm), where ~30% of the liposomes taken up by the cells were localized to the estimated nuclear region (Methods).

Extended imaging depth by PSF engineering
To demonstrate the generality of our approach to other PSFs, we inserted another phase-shaping optical element into the imaging path, replacing the cylindrical lens (Fig. 7a). The Tetrapod phase mask was designed to optimize the Fisher information for 3D localization 38 , and can be used for 3D measurements over extended depths (Fig. 7c).
Extending the depth of field can also have important applications in IFC; for example, where a sample needs to be maintained at a specific concentration, precluding the use of sheath fluid. Up to this point, we have used the extent of the astigmatism to order images by depth; however, this solution is not applicable to all 3D PSFs. A more general solution is generating a calibration curve with no assumed knowledge of the PSF. First, we collected 5,000 bead measurements in 23 separate datasets with narrowdepth distribution (5 µm), where the objective lens was moved forward by 2 µm between each recording. The images in each dataset were aligned by a centre of mass calculation on the basis of their intensities and averaged to produce an image template. Next, a representative image was identified in each dataset by calculating the 2D correlation of each image with the template. These selected images composed an intermediate calibration dataset, as they were already ordered by z. We then increased the flow diameter to 50 µm (Fig. 7b), imaging and ranking 13,562 PSFs by their similarity to the images in the calibration dataset. These ranked images were assigned a probable depth from the known flow dimensions (Fig. 7d) thus creating our shape-to-depth calibration. As a verification procedure, the experimental results were compared to numerical simulations of the optical system (Fig. 7e, see Methods).

conclusions
In this work, we have extended the capabilities of IFC to enable three-dimensional localization. This was made possible by using the high-throughput of IFC to derive a 3D model of the PSF from the underlying flow distribution. Furthermore, we have shown that our approach is fully compatible with existing commercial apparatus and requires only a small hardware addition to the device: the placement of a cylindrical lens or phase mask into the imaging path. Using this approach, we have increased the throughput of 3D-distance measurements of DNA compaction in live yeast cells by orders of magnitude over previously used traditional-microscopy methods 33 , and demonstrated the application of our method to nanomedicine; namely, quantification of nanoparticle uptake in lymphocyte cells.
Generally, the experiments that would benefit most from our method include those that: (1) require the full distribution of a heterogeneous sample (Fig. 4), (2) measure the dynamics of a population's response to stimuli (Fig. 5), (3) are affected by surface attachment, for example non-adherent cells (Fig. 6) or (4) aim to observe and quantify rare events. Specific applications include: evaluating drug delivery and uptake 34 , monitoring metabolicallyengineered pathways, correlative chromosome conformation capture measurements 26 and more.
Like all localization-based methods, a suitably low density of emitters is required so that the PSFs do not significantly overlap in the image. While certain types of objects, such as fluorescent beads, are ideal for this application, we made use of multiple spectral channels to colocalize closely spaced, selectively labelled regions of DNA. Using spectral separation of fluorescence emission, however, is not a fundamental requirement of this method. The recent application of machine learning, namely neural networks, to 3D localization microscopy has been shown by our group 39 and others 40 to be adept at extracting the positions of densely spaced emitters, and could represent one path to enable more complex samples. Favourably, neural nets also have the advantage of rapid-data processing, and have recently been applied to online analysis controlling cell sorting during IFC experiments in a custom instrument 41 .
In addition to new applications directly related to IFC, the approach described here may be used to calibrate any microscope's 3D PSF by using a flow system with well-known depth-distribution properties to produce a shape-to-depth calibration that could be applied to static samples. This approach would solve the longstanding challenge stemming from using calibration curves generated with surface-adhered objects while imaging is subsequently performed in media with a different refractive index 24,42 .
Related to IFC, the work presented here does not exhaust the full potential of combining PSF engineering to flow imaging. Future applications could use custom phase masks designs to optimally match the flow profile (z positions) in a particular experiment 23 , and make the depth determination more robust to variations in flowrates; for example, by encoding the depth in a one-dimensional stretch of the PSF 43,44 orthogonal to the flow, thereby decoupling the depth-estimation error from flow-rate induced errors. PSF engineering could also be used to enhance the autofocusing and coresize characterization by adding a dedicated 3D-imaging channel inside an instrument. Finally, it would be useful to incorporate fluorescence-activated cell sorting on the basis of subcellular 3D colocalizations through online analysis in the IFC operating software, creating a new approach for using selection-marker technologies.
In the context of nanomedicine, our technology enables an unprecedented combination of throughput and precision for the intracellular localization of nanoparticles. A high-throughput system that measures the exact position of therapeutic nanoparticles on cellular uptake would be critical for realizing the promise of directing nanomedicines to specific intracellular targets 36 .

Online content
Any methods, additional references, Nature Research reporting summaries, source data, extended data, supplementary information, acknowledgements, peer review information; details of author  contributions and competing interests; and statements of data and code availability are available at https://doi.org/10.1038/s41565-020-0662-0.

Methods
Sample preparation for flow cytometry. Before loading, samples were diluted to <2 × 10 5 objects per µl and mixed with fluorescent beads. This density balanced a high acquisition rate (≥100 objects per s) with a low probability of imaging overlapping particles. Fields of view containing two or more objects were either cropped or discarded.
To prepare fluorescent beads for surface adhesion to the human T lymphocyte cells (Jurkat, ATCC 77441), Concanavalin A (ConA; Sigma, catalogue no. C2010) was prepared as described previously 46 . Briefly, a stock solution was made in PBS (pH 6.5) containing 0.06 M CaCl 2 and 0.06 M MnCl 2 . Stocks were stored at −80. Before use, samples were thawed on ice and diluted with 200 µl PBS (pH 6.5) and 600 µl H 2 O. Fluorescent, 47 nm polystyrene beads (ThermoFisher, F10720) were added 1:1,000 to the ConA solutions. The mixture was then mixed with cells 1:10 for 10 min before imaging.
The Jurkat cells were cultured as described previously 47 , 37 °C, 5%/ CO 2 . Cells were incubated for 30-60 min with fluorescent liposomes (900 µl of 10 6 cells ml −1 with 100 µl of liposomes at 2.3 × 10 7 nanoparticles per ml). The supernatant, containing unbound liposomes, was removed by 3×, two inutes centrifugation steps at 1,000 relative centrifugal force, resuspending in Dulbecco's phosphate buffer saline. 10 min before imaging, cells were incubated with ConA-treated beads to promote surface attachment. Cells were selected for analysis based on the presence of beads and nanoparticles, as well as a round cell morphology, that is a circular brightfield image, to ensure the shape could be extrapolated from several 3D bead measurements. To assess the effectiveness of the extrapolated 3D surface profile, cells with more than two beads were used to assess the extracted radius. This resulted in an estimated radial error of 660 nm, defined as the standard deviation of (r bead -r fit ). This corresponds to ~13% error of the cell diameter (10.5 ± 1.4 µm).
All buffers and media used for cell growth and imaging was sterilized by filtering before use (0.22 µm).

Inserting the cylindrical lens and phase mask into the imaging flow cytometer.
A cylindrical lens (f = 1 m, Thorlabs, LJ1516RM-A) was mounted in a 2D micrometre translational mount (Thorlabs, LM1XY) and inserted into the imaging path of the ImageStreamX using a detachable magnetic base (Thorlabs, KB2X2). The position of the lens was adjusted until the astigmatic PSFs appeared symmetrical and oriented vertically/horizontally across all channels. The position along the optical axis affects the nature of the astigmatism shape, and was optimized iteratively.
The Tetrapod phase mask was mounted in a custom adapter designed to fit in the filter wheel located at a plane conjugate to the back focal plane of the objective. The mask position was refined in x and y until the PSF was symmetrical.
IFC settings. The default ImageStreamX parameters were used for data collection with the exceptions described next. The instrument was used with the ×60 magnification objective and a core diameter 7 µm (default, see Supplementary Information) for beads, DNA nanorulers and yeast measurements. For imaging T cells, the ×40 magnification objective was used with a core diameter of 10 µm (default). In the Tetrapod experiment, the core diameter was adjusted where noted and speed tracking was set to OFF. Autofocus was set to OFF (default was ON), and the ideal focal position was found after inserting the cylindrical lens. This was done by imaging short acquisitions ~2 min of fluorescent beads and minimizing the distance between the average focal plane (where the shape of the object is symmetric) and the centre of the flow so that the number of objects on either side of focus was approximately even. This value was typically in the range of (−2.5, 2.5) µm and appeared robust between measurements on the same day. Next, after enabling the Advanced Configuration setting, ObjectMaskCoeffParam was set to 0 (default was 0.8), and Keep Clipped Objects was set to ON (default was OFF) to retain the largest field of view possible; note that these two settings were only needed for the Tetrapod PSF, but were used in all measurements described for consistency.
For beads and cell data, the sheath buffer was 1× PBS; however, for the DNA nanoruler experiment, the sheath buffer was exchanged with 1× TAE containing 10 mM MgCl 2 to avoid possible changes in origami structure due to changes in ionic strength.
IFC data preparation for 3D analysis. There were two steps of data preparation before 3D localization: object classification and assignment of offset xy position.
Bead and non-bead objects were classified using the data in the exported feature file generated by the IDEAS companion software supplied with the instrument. The feature Image_MC<channel number> provided sufficient classification. According to the IDEAS manual, this feature is 'the sum of the pixel intensities in the mask, background subtracted' . The described classification scheme worked well due to the differences in fluorescence intensity and ratio between some spectral channels of the fluorescent beads relative to both the cells and the nanorulers. More general classification schemes might require use of the full image data rather than relying on exported feature parameters. In cases of similar calibration and test samples, it may be necessary to measure the PSF calibration separately; where we did observe some subtle differences in the PSF calibration after stopping and restarting the flow, this may be due to hysteresis in the objective position movement that resets on changing samples. We therefore recommend collecting calibration and sample data concurrently.
The absolute xy position of an emitter with respect to the relevant chargecoupled device (CCD) channel is the sum of its position in the image frame and the offset of the image frame with respect to the CCD channel. The offsets in x and y are the exported feature parameters Raw_Centroid_X and Raw_Centroid_Y. Note that retaining the absolute x and y position for each object is necessary to correct for any systematic bias that occurs while fitting the astigmatic PSF when the cylindrical lens is not perfectly aligned.
Software. IFC datasets (.cif files) were generated using INSPIRE software (Amnis). Postexperiment, feature data for all objects was exported to.txt format using IDEAS software (Amnis). The feature data was useful for rapid classification of objects as either calibration beads or cells/DNA nanorulers.
MATLAB (Mathworks, v.2018b) with Bioformats package 48 was used for analysis of all .cif data from cell and nanoruler samples.
The FIJI distribution of ImageJ 49 was used with the Thunderstorm plugin 50 for 2D image analysis of the DNA nanorulers.

Localization of objects.
For data acquired with the cylindrical lens, intensities I(x,y) of fluorescent objects were fit using nonlinear least-squares (MATLAB's lsqnonlin function) to the following 2D-Gaussian model: where B is a constant background for bead and in vitro DNA data, A is the maximum intensity, and (x 0 ,y 0 ) is the xy position of the emitter. The widths, 2σ 1 and 2σ 2 , of the 2D Gaussian and its rotation θ from the CCD axis are related to a, b and c as follows: We defined the axes of the 2D Gaussian as follows: the positive σ 2 axis is rotated 90° clockwise from the positive σ 1 axis. The positive axis of σ 1 may be rotated in the range of (-45°, 45°) with respect to the positive y axis of the CCD. This results in the positive axis of σ 1 aligned closer to the y axis of the CCD and the positive axis of σ 2 aligned closer to the x axis of the CCD.
Emitters in cells were fitted similarly, except the background was also a 2D-Gaussian function, whose shape parameters were restricted to being larger than that of the localized emitters inside the cell: Parameterizing the Gaussian fits and localizing in 3D. To parametrize the PSF shape (σ 1, σ 2 ) as a function of the emitter depth, one or more channels were used for calibration. After fitting with a freely rotating astigmatic Gaussian, the precise orientation angles of the Gaussian were found per channel by taking the mean angle of calibration beads whose shape parameters were between 1 and 3.5 pixels with a finite ratio of 1.1/2 pixels. The average angle for each side of focus was determined separately to account for aberrations in the imaging system that distort the PSF. Images were relocalized with these orientations.
To create a preliminary ordering of the calibration objects images by depth, a 16-parameter function was used to generate a 2D image to match the 2D histogram of shape distributions, Fig. 1e. Once a curve that approximated the range of object shapes was found, the data was projected to the closet point along the curve to find the relative order of the objects. If more than one channel was used for calibration, the average object position was used. In each channel analysed, objects were binned by relative depth to find the average shape parameters for data percentiles. Finally, a spline function was employed to create a continuous shape as a function of percentage along the flow (0-100%), Fig. 1f. To convert the percentile in the flow to absolute z positions, we used the MATLAB function norminv that outputs a position as a function of percentile for a normal distribution of a given mean (in our case, 0) and standard deviation (~1/4 of the flow diameter).
Validating the depth distribution within the stream. The distribution of particle density as a function of position within the core was determined as follows: first, the core diameter was set to either 5 or 7 µm. Next, fluorescent bead data was collected in the absence of an astigmatic lens at different objective positions, over a range of 10 µm. The fraction of objects in focus at each distance from the objective f was defined as: where N total is the total number of objects collected, and N focus is the number of objects for which ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi σ 2 1 þ σ 2 2 p < 330 nm I . The in-focus threshold was chosen to be small enough so that P f pðf Þ < 1 I , indicating that objects were not defined to be in focus for more than one focal position f, and large enough to include a large number of objects and thus ensure a statistically small error for p(f). The resulting p(f) values were fitted to a Gaussian distribution (Fig. 1d), resulting in mean μ fit and standard deviation σ fit . We found σ fit corresponded well to a quarter of the core diameter setting in the INSPIRE software, and thus used the core diameter setting for all measurements (Supplementary Fig. 6).

Colour-channel registration.
To compare localizations between different spectral channels, that is different parts of the detector, it is necessary to perform subpixel, cross-channel registration. Due to field-dependent aberrations that might affect the channels differently, we divided the imaging volume into 8 × 8 × 8 voxels and corrected for the average displacement between two channels. This data was then interpolated in 3D using the MATLAB function interp3, to create a function that corrects each localization by the estimated bias at that point in the field, x 0 y 0 z 0 ¼ f ðx; y; zÞ I ( Supplementary Fig. 7).
Removing outlier localizations and low signal-to-noise ratio objects. For calibration, fluorescent objects were required to be localizable in each channel used to determine the z distribution, and had to have an amplitude at least three times larger than the standard deviation of the background.
For in vitro DNA experiments, objects with a measured 3D distances >0.5 µm, comprising <2.7% of the total analysed objects measured in flow and <0.3% of objects on the coverslip, and were excluded when calculating the mean and standard deviations of the distributions.
For yeast measurements, a minimum amplitude of the Gaussian fit was required to ensure that the spot was localizable (11) digital units and chosen because the distribution broadened at lower thresholds but did not change significantly at higher thresholds.
For lymphocyte measurements, some fluorescence images did not match the 3D calibration curve, which represents the PSF of a single, well-separated object. Objects with localized shape parameters that had a geometric distance to the calibration curve greater than 225 nm, ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi σ 2 1 þ σ 2 2 p < 225 nm I (0.45 pixels), were removed. This could reflect significantly defocused objects, multiple closely spaced objects or periods of flow instability. Additionally, objects that were determined to be significantly outside the cell volume, >125% cell radius, were excluded from further analysis. This was made possible by our 3D analysis.

Intracellular-coordinate determination.
For intracellular localization of nanoparticles, brightfield images (Channel 1) were thresholded and evaluated using the MATLAB function Regionprops to extract the cell's 2D centre position (x and y), circularity, eccentricity and ellipsoidal radii. Objects with an average radius of <4 µm (8 pixels), an eccentricity of >0.6 or circularity of <0.90 were removed from further analysis. These images usually corresponding to debris, multiple cells or non-circular cells. Fluorescence images were then analysed to determine the estimated number of fluorescent surface-bound beads (Channel 2) and liposomes (Channel 5). Cells with two or more beads and one or more liposomes were then localized by fitting cell and image-offset parameters together with N asymmetric-Gaussian fits, where N was the number of objects identified. The Gaussian shape parameters were then compared to a calibration curve to extract the z position as described previously.
The 3D cellular coordinates were determined using the localized positions of the surface beads. Specifically, minimizing the differences of the squared 3D radial displacements using the MATLAB function lsqnonlin to identify the cell's z centre and radius. Note the x and y centre positions were found using the brightfield image. Object radial positions were normalized between the cell's centre position and surface. Note that the surface-bead radii displayed in the histogram were taken from cells with three-or-more localized beads.
Nuclear-volume estimation. Nuclear-stained Jurkat cells were analysed to estimate the nuclear volume. Brightfield images (Channel 6, Supplementary  Fig. 5a) were used to find the 2D projection. The nuclear area was stained by a DNA stain (Channel 1, Supplementary Fig. 5b). The cytoplasm was defined as the cellular area not covered by the nucleus; Supplementary Fig. 5c. The relative nuclear volume was defined as the fraction of nuclear measurements at a given radial bin (Supplementary Fig. 5d). This metric was validated by 2D and 3D simulations of randomly positioned spherical nuclei in cellular spheres. To estimate the fraction of nanoparticles contained within the nucleus, the relative nuclear-volume probability curve was multiplied by the normalized histogram of nanoparticle positions.

Simulation of the Tetrapod PSF.
To simulate the Tetrapod PSF, a numerical model of the imaging system was constructed on the basis of a scalarapproximation model described previously 23 . The simulated emission path consisted of an objective lens (60×, numerical aperture 0.9), a tube lens and relay lenses which form the diffraction-limited image on the detector. The phase mask was inserted in a conjugate back focal plane to the objective. The physical phase mask was simulated as a phase-modulating element combined with an iris described by: where r,ϕ are the polar coordinates in the back focal plane and Φ is the phase mask design 23 (Fig. 7c). Due to the refractive index mismatch between the sample medium (similar to water, n = 1.33) and objective's immersion medium (air, n = 1), a depth change of the emitter manifests as a stronger defocus in phase than an equidistant change in objective position. The distance of the objective relative to the flow cell was therefore determined empirically.
(2) DNA nanorulers were designed and made by GATTAquant GmbH.
(3) KWY4069 yeast cells, provided by K. Weis, were not authenticated and were not tested for mycoplasma. (4) Jurkat cells, provided by D. Yablonski, were authenticated to be precursors of Jurkat Clone E6.1 (Cellosaurus no. CVCL_0367). Authentication was performed using the Promega GenePrint 24 System to determine short tandem repeat profile of 23 loci plus Amelogenin for gender determination (X or XY) and analysed using the 3500xl Genetic Analyzer (Life Technologies) and GeneMapper IDX software. The cells were tested to be mycoplasma free.
Reporting Summary. Further information on research design is available in the Nature Research Reporting Summary linked to this article.

Data availability
The data that supports the plots within this paper and other findings of this study are available from the corresponding authors on reasonable request.

code availability
The analysis scripts for image categorization, calibration, localization and 3D distance measurements were written in MATLAB 2018b (Mathworks) and are available from the corresponding authors on reasonable request.  For a reference copy of the document with all sections, see nature.com/documents/nr-reporting-summary-flat.pdf

Life sciences study design
All studies must disclose on these points even when the disclosure is negative.

Sample size
An estimated sample-size requirement calculation was performed by simulation and is shown in Supplemental Fig. 1. The result dictated that 1000s-10,000s of calibration objects are needed in a measurement to accurately determine the position of a sample object in 3D.
The main objective of our work is to complete as many measurements as possible as quickly as possible. This is limited in several ways: Samples must be dilute enough both not to overlap in images and not to cause aggregation. The acquisition itself is limited by the instrument settings, namely a maximum of 200 microlitres of a sample can be loaded at a time. In certain cases, the ImageStreamX also re-engages the 'autofocus' function, specifically if more sheath fluid needs to be loaded. Since the 3D calibration is focus-position dependent, re-engaging the autofocus represents a termination of experiments. This typically limited experiments to ~10-30 minutes depending on the requirements. Longer experiments were possible, but at the increased risk of equipment malfunction. For the data presented in Figs. 1 & 2, fluorescent beads can be measured with extremely high thoughput and the total number of measured objects was preselected to be 50,000, which is enough to not limit our method's precision.
For the data presented in Fig. 3, the number of samples was determined by the requirement to perform several repeats, and thus the acquisition was limited by the amount of materials available.
For the data presented in Fig. 4, cells were imaged for ~10 minutes or until the instrument disengaged autofocus and the measurement was halted.
For the data presented in Fig. 5, cell treatments were measured periodically to determine the approximate duration of the experiment. Cell densities were determined by the growth rates in the respective medias.
For the data presented in Fig. 6, the measurement was limited by the maximum instrument runtime. For the data presented in Fig. 7, the enlarged flow width, the vastly increased flow rate limited measurements to only several minutes.
Data exclusions Three types of exclusions were enforced (1) Pre-selection of object attributes for acquisition, (2) discarding data due to instrument malfunction, (3) post-acquisition selection. 1) During imaging, exclusions were made for object profiles that were either (a) not sufficiently bright to be localized in all required imaging channels, or (b) appeared too large to be single objects. These criteria were preselected based on visual inspection of imaged samples that were not used for subsequent analysis. These pre-selected exclusions were necessary to contain the filesizes to reasonable values <4GB.
2) Due to the stringent focus requirements of our measurements, any instrument disruption, usually caused by re-engagement of the autofocus or sample clogging, required acquisition to be immediately terminated. In cases where (a) insufficient numbers of calibration beads were recorded, or (b) the instrument disruption was not immediately realized, the entire dataset was discarded.
3) Post acquisition, exclusions were made for images containing more than one object and/or cell, insufficient signal, and poor fits. Details are provided in the manuscript. The decision to exclude closely-spaced or overlapping objects was pre-established. The signal levels that proved insufficient for analysis as well as identification of poor-localizations (dissimilar to calibration objects) were established after recording and analyzing a subset of the data.

Replication
Due to the inherit high throughput of our method, all objects and cell measurements were repeated 100s-1000s of times. Independent recording replicates were also used to verify the robustness of the experimental data in the manuscript. For the data presented in Figs. 1 & 2: All replicates of the Tetraspeck bead analysis succeeded with similar but slightly varying levels of precision per run based on the optimization parameters described in the manuscript. No absolute replication failures were identified. For the data presented in Fig. 3: Replicates successfully produced qualitatively similar results; albeit with lower measurement precision in the direction of flow. Using the correct buffer (see Methods section) proved to be critical for overcoming nanorod aggregation. For the data presented in Fig. 4: Replicates successfully produced qualitatively similar results; With regards to imaging buffer: water and unoptimized growth medias with high fluorescence backgrounds proved to be deterious to this experiment. For the data presented in Fig. 5: All replicates succeeded in producing qualitatively similar results. For the data presented in Fig. 6: All replicates succeeded in producing qualitatively similar results. For the data presented in Fig. 7: Replicates successfully produced qualitatively similar results; importantly, the "track flow-speed" instrument option must be disabled.
Randomization For calibration and comparison to fluorescent beads in each of the measurements, data was assigned to a control or calibration group using a random number generator in Matlab 2018b (Figure 2-6) or using alternating measurements to assign groups ( Figure 1). No selection, and therefore no randomization was performed in the data described by Figure 7.

Blinding
Blinding was not relevant for the data described in Figs. 1, 2, 3, 6 and 7, which represent single-sample types.
For the data described in Figs. 4 & 5, the analysis algorithm was performed automatically without regard to sample type; All group allocations were performed randomly after data was acquired and all data comparisons made were analyzed using identical analysis algorithms, with the exception that cells and beads shown in Fig. 4 and Supplementary Fig. 4, had a different background term to account for cellular autofluorescence as described in the main text; however, this was not the primary comparison under consideration.
Reporting for specific materials, systems and methods We require information from authors about some types of materials, experimental systems and methods used in many studies. Here, indicate whether each material, system or method listed is relevant to your study. If you are not sure if a list item applies to your research, read the appropriate section before selecting a response.

Authentication
The Jurkat cells used were authenticated to be those from which the Jurkat Clone E6.1 was originated (Cellosaurus no. CVCL_0367). Authentication was performed using the Promega GenePrint 24 System in order to determine short tandem repeat (STR) profile of 23 loci plus Amelogenin for gender determination (X or XY) and analyzed using the 3500xl Genetic Analyzer (Life Technologies) and GeneMapper IDX software.

Mycoplasma contamination
Cell lines tested were confirmed not to contain mycoplasma contamination.
Commonly misidentified lines (See ICLAC register) No commonly misidentified cell lines were used.