## Abstract

3D single-molecule localization microscopy relies on fitting the shape of point-spread-functions (PSFs) recorded on a wide-field detector. However, optical aberrations distort those shapes, which compromises the accuracy and precision of single-molecule localization microscopy. Here, we employ a computational phase retrieval based on a vectorial PSF model to quantify the spatial variance of optical aberrations in a two-channel ultrawide-field single-molecule localization microscope. The use of a spatially variant PSF model enables accurate and precise emitter localization in *x*-, *y*- and *z*-directions throughout the entire field of view.

© 2019 Optical Society of America under the terms of the OSA Open Access Publishing Agreement

## 1. Introduction

Super-resolution microscopy represents a family of techniques employed to surpass the diffraction-limited spatial resolution of light microscopy through optical or computational means [1–3]. Single-molecule localization microscopy (SMLM) is one such technique, which has developed into a powerful tool in molecular and cellular biology [4–6]. In camera-based single-molecule localization microscopy, molecular positions are usually estimated by fitting a theoretical model to the experimentally measured point-spread-functions (PSFs) of single fluorescent emitters. Determination of molecular positions with a precision of tens of nanometers enables the sub-cellular localization and co-localization of fluorescently labeled biomolecules, as well as single-molecule tracking [7]. The accuracy and precision of molecular position estimation is of prime importance, because these metrics affect all subsequent steps in data analysis, such as estimation of cluster formation in pointillistic single-molecule localizations [8,9], or molecular displacement estimation in single-molecule tracking [10–13]. The accuracy of the estimated position, i.e. its systematic deviation from the true position, is primarily dependent on the choice of the theoretical PSF model used for fitting the measured fluorescence intensity data [14]. The precision of the estimated localization, i.e. its random statistical variability upon repeated measurements, is primarily determined by the signal-to-background ratio of the fluorescence intensity data [4–6]. Efforts to improve localization accuracy and precision require minimizing systematic and random errors during data acquisition and computational post-processing [15].

A major advancement of single-molecule localization microscopy has been the ability to use PSF engineering to encode the 3D position of an emitter and/or its emission wavelength [16–18] or its molecular orientation [19–21]. To achieve emitter localization below the diffraction limit in both the lateral (*x*,*y*) and axial (*z*) dimensions, PSFs have been engineered to change shape rapidly along the z-dimension, either through stretching or through rotation [22–30]. For example, the double-helix point-spread-function (DHPSF) was engineered to produce two high-intensity lobes that revolve around the emitter’s *x*,*y*-position in a *z*-displacement dependent manner [24,25]. Careful calibration measurements have revealed that the fitted *x*,*y*-position, obtained using analytical PSF models, deviates slightly from the true lateral position of an emitter as a function of *z*-displacement. This ‘wobbling’ behavior of double-helix and astigmatic PSFs can be calibrated and corrected to obtain more accurate estimates of an emitter’s lateral position [31].

The large chip size of modern sCMOS camera has enabled SMLM in ultrawide fields-of-view (diameters of up to 250 µm) [32–34]. Ultrawide-field imaging combined with fast frame rates (up to 100 Hz for full chip readouts), has helped increase the overall data acquision throughput [35]. However, the high-NA objective lenses used for SMLM cannot offer full aberration correction everywhere in the field-of-view. It is well established that optical aberrations become more pronounced as emitters are localized far away from the central optical axis [32,36–39]. Thus, the prevalent assumption of a spatially invariant PSF in the imaging system does not hold and the PSF models used in high-precision, high-accuracy SMLM need to be treated as spatially variant.

Previous work extracted emitter positions from experimental PSFs using simple analytical, interpolated, or scalar/vectorial PSF models. Simple analytical models typically use Gaussian functions to approximate the experimentally measured fluorescence intensity profiles [4–6]. More sophisticated models based on the Gibson-Lanni model [40] or Hermite functions [41] have also been employed. However, simple analytical functions are not flexible enough to accurately model the complex shapes of aberrated PSFs encountered in real microscopes. To overcome these limitations, piecewise polynomial functions have been used to interpolate experimentally acquired PSFs [42–45]. Such interpolated model functions contain larger numbers of adjustable parameters compared to the simple analytical functions, which makes them more flexible to approximate the PSF intensity profiles. Finally, models based on scalar or vectorial diffraction theory explicitly account for the propagation of light through high numerical aperture imaging systems [16,18,46]. Augmented vectorial models also include the glass-water refractive index boundary at the microscope coverslip to enable estimation of PSF shapes for emitters located within a refractive index mismatched medium [47,48]. However, neither interpolated or diffraction-theory-based models have addressed the spatial variance of experimental PSFs in SMLM.

Here, we use a vectorial PSF model and refine it using computational phase retrieval to quantify the spatially variant optical aberrations in two different color channels of a single-molecule localization microscope. We show that the spatially variant vectorial PSF model enables accurate and precise 3D emitter localization throughout an ultrawide field-of-view. We further show that the same emitters can be co-localized in two different color channels with an accuracy less than 25 nm, which matches localization precisions typically achieved in live-cell SMLM. Our work provides a quantitative framework for processing multi-channel 3D SMLM data acquired with arbitrarily engineered PSFs in an ultrawide field of view.

## 2. Theory and methods

#### 2.1. Optical setup

Experiments were performed on a custom-built dual-color inverted fluorescence microscope based on the RM21 platform (Mad City Labs, Inc) [Fig. 1]. Immersion oil was placed between the objective lens (UPLSAPO 60X 1.4 NA, Olympus) and the glass cover slip (#1.5, 22 mm x 22 mm, Schott). The sample was excited with a 514 nm laser (Genesis MX514 MTM, Coherent) and a 561nm laser (Genesis MX561 MTM, Coherent). All excitation lasers were circularly polarized by zero order quarter wave plates (WPQ05M-514 and WPQ05M-561, Thorlabs), and a bandpass filter filtered the spectral profile of the 514 nm laser (ET510/10bp, Chroma).

Fluorescence emission was passed through a shared filter set of LP02-514RU-25 (Semrock), NF03-561E-25 (Semrock), and ET700SP-2P8 (Chroma). A dichroic beam splitter (T560lpxr-uf3, Chroma) was then used to split the emission pathway into ‘red’ and ‘green’ channels, corresponding to the transmitted and reflected light, respectively. The ‘red’ transmitted channel contained an additional 561nm notch filter (ZET561NF, Chroma) to block scattered laser light. Each emission path contained a wavelength-specific dielectric phase mask (Double Helix, LLC) placed in the Fourier plane (FP) of the microscope and mounted on high-precision XYZ translation stages (9062-XYZ, Newport) to generate a double-helix point-spread-function [25,49]. The fluorescence signals in both channels were detected on two separate sCMOS cameras (ORCA-Flash 4.0 V2, Hamamatsu).

#### 2.2. Sample preparation and imaging

For samples consisting of calibration beads, we diluted 100 nm tetraspeck fluorospheres (Invitrogen) diluted in PBS at a ratio of 1:3000 (v/v) and then spin-coated them onto a microscope coverglass. To maintain a refractive index similar to water, we covered the beads with a 1.5% (w/v) low-melting point agarose (Fisher Scientific) pad made with PBS buffer to mimic experimental conditions used for live cell imaging [12,50,51]. To generate a sample of beads immobilized away from the microscope coverslip, we suspended fluorescent beads into molten 1.5% agarose and spotted a 10 µl of mixture onto the coverglass which was previously spin-coated with fluorescent beads.

The phase mask lateral position was optimized prior to each experiment by scanning the nominal focal plane through a single fluorescent bead immobilized on the microscopy coverlip. The resulting PSF sections were analyzed using the Easy-DHPSF software [52]. The phase mask lateral position was adjusted to ensure even lobe intensities and minimal lateral shifts of the fitted positions at different positions of the nominal focal plane (wobble effect).

Calibration images to determine the spatially variant aberrations profile were acquired by scanning the nominal focal plane through fields-of-view containing many immobilized fluorescent beads. Specifically, we stepped the distance *d* in 50 nm increments through a total range of 3 µm and recorded 10 PSF images at each step. All images were acquired with an exposure times of 30 ms and illumination intensities of 3-50 W/cm^{2} on the sample.

#### 2.3. Modelling an isotropic emitter

We employ a vectorial model of light propagation to simulate the PSF of an emitter in a refractive index mismatched medium [47,48]. Briefly, we simulate polarized emission of a dipole emitter and propagate the light through the refractive index boundary. The refracted light waves are then collected by a high numerical aperture (NA) objective lens. This vectorial PSF model provides the amplitudes and phases of the light waves in the back focal plane (BFP) of the objective lens, i.e. the pupil function [Fig. 1]. As detailed further below, the phases of these waves are dependent on the position of the emitter, whereas the amplitudes are not. The light waves are then further propagated through the imaging system (one tube lens and two 4f lenses) using Fourier optics to ultimately yield the final intensity profile in the image plane. To model the PSF of an isotropic emitter, the PSFs of three specific dipole orientations are summed in the image plane [47,53]. We use this model to compute the PSFs of isotropic emitters at arbitrary *x*,*y*- and *z*-positions above the refractive index boundary [Fig. 1, inset]. The mathematical details of the model are described in the following.

The refraction and transmission of the emitted light at the refractive index boundary is described, respectively, by Snell’s law, *n _{1} sinθ_{1} = n_{2} sinθ_{2}*, and the Fresnel coefficients. The Fresnel coefficients for transmission of

*s*- and

*p*-polarized light are given by

*n*=

_{1}*n*≈

_{glass}*n*= 1.52,

_{oil}*n*=

_{2}*n*= 1.33. The variables

_{water}*θ*and

*ϕ*correspond to the polar and azimuthal angles of emitted light rays. The polar angle

*θ*maps to the BFP polar coordinate

*ρ*through

*ρ*= sin

*θ*=

_{1}*n*sin

_{2}*θ*, whereas the azimuthal angle

_{2}/ n_{1}*ϕ*maps one-to-one to the BFP azimuthal coordinate. The numerical aperture of the objective lens defines the range of angle

*θ*that can be collected and thus the maximum radius

_{1}*ρ*in the BFP according toThe emitter’s dipole orientation can be defined in both Cartesian and spherical coordinates according to

_{max}*x-*and

*y-*polarized light, ${E}_{x}^{A}$and ${E}_{y}^{A}$, are calculated using Eqs. (5) and (6), respectively.

The phases of the electric field in the Fourier plane (FP) (the location of the phase mask, see Fig. 1) depend on several factors, namely the *z*-position of the emitter above the refractive index boundary [47]

*d*between the nominal focal plane and the refractive index boundary (modified from [46])

*x*,

*y*-displacement of the emitter from the central optical axis (modified from [46])

*λ*is the emission wavelength and

*M*is the magnification of the imaging system.

The overall phase in the FP is thus given by

and the complex amplitudes in the FP are then given byandTo compute the PSF in the image plane, we propagate the light wave through the rest of the imaging system using Fourier optics

*P(u,v)*specifies the phase manipulation, which, in our case, alters the standard PSF to the DHPSF.

The model described so far produces DHPSFs assuming no optical aberrations. To account for optical aberrations, Eq. (13) is modified to

*j*),

*a*is the spatially variant aberration coefficient, and Z

_{j}(x,y)*is the Zernike polynomial normalized to π.*

_{j}In the vectorial PSF model, we assume that the fluorescence emission from each emitter in our case is isotropic. Isotropic emitters are modelled by summing the PSFs of three anisotropic emission dipoles with orientations angle Φ and Θ of (0, π/2), (π/2, π/2), and (π/2, 0) (See Eq. (4)). Additionally, a Gaussian blur with a radius of 0.58 pixels (62.6 nm) is added to the simulated image to account for the finite spherical volume of fluorescent beads. This blur radius was determined by minimizing the difference between a simulated 100 nm sphere filled with isotropic emitters and a blurred isotropic point emitter in the absence of noise and aberration.

#### 2.4. Data analysis

We use two different methods to fit the DHPSFs in order to extract the *x*,*y*,*z*-positions of individual emitters. The first method is the analytical model built into the Easy-DHPSF fitting routine, namely a model consisting of two laterally displaced Gaussian functions [52]. This routine was modified so that both least squares (LSQ) and/or maximum likelihood estimation (MLE) can be performed [35]. The other method uses the spatially variant vectorial PSF model described above and MLE for PSF fitting. To determine the *x,y,z*-positions of individual emitters accurately using the vectorial PSF model, the distance *d* of current frame and the aberration coefficients *a _{j}(x,y)* (

*j*= 4, … 15) of each emitter must be determined as described in the following.

## 3. Results and discusson

#### 3.1. Quantifying spatially dependent aberrations in large fields-of-view

To determine the spatial dependence of optical aberrations in a large field of view, we immobilized fluorescent beads on glass coverslips and scanned the nominal focal plane of the objective lens over the beads in equally sized (50 nm) steps. This procedure produced a set of 120 averaged images, each acquired at different defocus throughout the microscopes’ field of view. We then fit a 2-dimensional Gaussian function to the rotation angles of varying DHPSFs in each image and defined the apex of this Gaussian function as the position of the optical axis. The bead closest to the optical axis, which is assumed to be the least aberrated, was then used to estimate the nominal focal plane position *d* in each image, as well as the phase mask’s rotation and lateral displacement from the optical axis. MLE was used to minimize the difference between the experimental PSFs and unaberrated vectorial PSF model. To estimate the position of the nominal focal plane relative to the refractive index boundary, the *z*-position relative to the refractive index boundary of the bead was assumed to be 50 nm (the diameter of the beads in our calibration measurements is 100 ± 6 nm, according to the manufacturer). The *x*- and *y*-postion of the abovementioned bead was also used as fitting parameters in this step.

Next, we quantified optical aberrations by estimating the coefficients *a _{j}(x,y)* of the Zernike polynomials with

*j*= 2,…,15 (Noll indices) at all bead locations by minimizing the difference between experimental PSFs and simulated PSFs, again using MLE [Fig. 2(a)], as described previously [18,46]. After initial fitting, we did not consider the Zernike polynomials with Noll indices 2 and 3 further, because they correspond to tip and tilt, which do not change the shape of the PSF. To estimate the Zernike coefficients

*a*with

_{j}(x,y)*j*= 4,…,15 at any location in the field of view, we fit a two-dimensional polynomial function of low order (usually 2-4, not to be confused with the Zernike polynomials) to the experimentally determined coefficients [Fig. 2(b)] [54]. The order of polynomial surfaces was increased, if necessary, to minimize any spatial dependence of residuals [Figs. 2(c) and 2(d)], specifically order 2 for

*j*= 4-6, order 3 for

*j*= 7-10 and order 4 for

*j*= 11-15. The distributions of residuals was analyzed to identify outlier beads with coefficients outside the mean ± 3 standard deviations range. Such beads were iteratively removed from the fit until all remaining beads fell within the above range. Excluded beads were removed from the surface fits of all other Zernike polynomials as well. We determined that outliers were due to fluorescent beads wiggling in the agarose gel, beads immobilized in the agarose instead of on the coverslip, overlapping DHPSFs, or intensity bleeding from neighboring beads. Also, outliers were evenly distributed throughout the field of view, ruling out any dependence on the signal-to-background ratio. The final polynomial surfaces can then each be queried at any location of the field of view to obtain the local aberration coefficients. The combined aberration coefficients provide the overall spatially dependent wavefront aberration phase $W(x,y;\rho ,\phi )$, which allows us to compute spatially variant aberrated PSFs.To verify the accuracy of the polynomial fits, we performed model cross-validation, compared estimates with different signal-to-background ratios, and performed simulations with beads at different

*z*-positions. First, we randomly chose a subset of the acquired fields-of-view and performed the same polynomial surface fitting. The differences between surfaces created from all and subset of fields-of-view were not significant, i.e. they were within one standard deviations of the distribution of residuals [Fig. 3(a)]. Second, we performed the same estimation to the same field-of-view using experimental data with different signal and background levels. Again, we found that the resulting surfaces were the same within errors [Fig. 3(b)]. Third, we considered the effect of bead size heterogeneity on coefficients estimation via simulation. Aberrated images were simulated using the spatially dependent aberration profiles determined above and emitter

*z*-positions of either 20 nm or 80 nm. The aberration coefficients were then estimated by assuming the

*z*-position of 50 nm as described for experimental data. The difference of estimated coefficients relative to the ground truth was within one standard deviation of the coefficient residuals for every considered aberration except for defocus [Fig. 3(c)]. The defocus and primary spherical aberrations are strongly coupled with the distance

*d*and emitter

*z*-position, which explains the larger differences in the aberration coefficients

*a*

_{4}and

*a*

_{11}. However, these simulations overestimated the bead size variation, which is 6 nm in diameter for the unstained microspheres according to manufacturer. When we repeated the same simulations with 6 nm depth variations, all Zernike coefficient differences were withing one standard deviation of the coefficient residuals. We also investigated to what extent the precision of coefficient estimation limits the 3D localization accuracy. When individual coefficients were computationally shifted one standard deviation of the residuals away from their original values, the deviation, defined as the differences between localized positions and the ground truth, was less than 5 nm at

*x*,

*y*-dimensions in all cases. The deviation in the

*z*-dimension was less than 10 nanometers except for defocus (

*j*= 4) or primary spherical aberration (

*j*= 11) [Fig. 3(d)]. Based on these experimental and simulated measurements, we conclude that polynomial surface fitting assuming a bead

*z*-position of 50 nm is robust and yields accurate estimation of the spatially variant aberration coefficient profiles through a large field of view.

Finally, we performed the same estimation and fitting procedure in the second color channel of our microscope. The shapes of the corresponding surfaces between the two color channels are qualitatively similar, differing only in overall magnitude [Fig. 4]. We conclude that aberrations that are described by low order Zernike polynomials originate predominantly from optical components that are shared by both color channels of our microscope.

#### 3.2. Estimating the z-positions of emitters on flat coverslips

To estimate the *x*,*y-* and *z*-positions of all beads in the field of view, the position of the nominal focal plane has to be determined in every frame. We estimated the distance *d* between the refractive index boundary and nominal focal plane using a reference marker with known *z*-position, such as a fluorescent bead of known size immobilized on the glass coverslip (*z* = 50 nm for a 100 nm fluorescent bead). Using the same reference marker in each frame, we then tracked changes in *d* over time due to stage drift [Fig. 5(a)]. Knowledge of the nominal focal plane position allowed us to estimate the *x*,*y-* and *z*-positions of all other beads in the field of view in every frame. We found that the estimated bead depths covered a wider range (stdev = 68 nm, range = [2, 343] nm) than expected, given the <10 nm bead size variation and the statistical localization precisions, defined as the standard deviations in the estimated *x*,*y-* and *z*-coordinates of each bead [blue localizations in top panel of Fig. 5(b) and Fig. 5(c)]. We noticed that several beads exhibited large localization precisions [colored in red in bottom panel of Fig. 5(b) and Fig. 5(c)]. These beads were localized at larger depths and/or showed noticeable translational motion from frame to frame. Not considering these insufficiently immobilized beads, the averaged localization precisions were *σ _{x}* = 2.4 nm,

*σ*= 2.9 nm, and

_{y}*σ*= 3.6 nm.

_{z}The major contributor to the observed bead depth variance is the fact that some beads can get displaced into the agarose during sample preparation. We speculate that a second contributing factor is the thickness variations of the glass coverslip. Microscope coverslips are not polished to optical flatness, which for high quality optics can reach λ/20~25 nm. A third contribution could come from the use of inaccurate aberration coefficients. To rule out this possibility, we displaced the field of view by 10 µm and estimated the depths of the same beads again at their new positions in the field of view [yellow localizations in top panel of Fig. 5(b)]. The average difference between bead depths before and after lateral displacement was 0.6 nm [Fig. 5(b), bottom panel]. We conclude that our method is robust in estimating emitter depths accurately throughout the field of view. However, random thickness variations of the coverslip, which appear to be on the order of tens of nanometers, remain uncorrected.

#### 3.3. Estimating the z-positions of emitters in agarose and correlation between two color channels

Imaging different components of the same structure in two different color channels provides information about their spatial proximity. However, due to different amount of aberrations present in different color channels, the images from two color channels are distorted differentially, limiting co-localization accuracy. Previous work has established the use of 3D transformation functions to correct the distortion and register two 3D super-resolution data sets acquired in different color channels [50]. In order to achieve 3D registration errors of less than 10 nm, transformation functions have to be computed based on large numbers of control points distributed throughout the 3D image volume. Such an approach leads to increased computational cost when imaging control points of similarly high density in an ultrawide field of view.

We hypothesized that the use of spatially variant vectorial PSF models, with aberrations in two separate color channels calibrated, would yield depth accuracies sufficient to justify the use of 2D transformation functions to register 3D localization data. To test this hypothesis, we embedded fluorescent beads in agarose and compared their depths estimated in two color channels of our microscope [top panel of Fig. 6(a), red and green dots, the same bead used as reference marker is colored in black]. The beads were distributed widely across the field of view [Fig. 6(b)], and the depths of the same beads matched up well with a mean difference of 0.8 nm and a standard deviation of 23.3 nm [Fig. 6(c)]. We conclude that systematic *z*-colocalization errors can be eliminated and *z*-colocalization accuracies of about 20 nm can be achieved using the spatially variant vectorial PSF model. We also note that the differences of more than 20 nm still occur for four of the beads considered here [bottom panel of Figs. 6(a) and 6(c)]. This observation may result from the existence of higher-order aberrations or sample-induced aberrations that are not accounted for by our model. To reduce *z*-colocalization errors into sub-10 nanometer regime, a locally calibrated 3D transformation function could be employed, possibly with lower control point densities [50], but we did not explore this option further.

#### 3.4. Comparison to simple analytical models

The above measurements suggest that localization precision of a few nanometers can be readily achieved at high signal-to-background ratios. To assess the improvements in accuracy and precision of position estimates made with the spatially variant vectorial PSF model, we compared its performance to that of the double Gaussian model used in the Easy-DHPSF algorithm [52].

For experimental data of fluorescent beads immobilized on coverslips, we found that the Easy-DHPSF algorithm yielded bead *z*-positions that mapped out a bowl shape [Fig. 7(a)]. Emitters at the periphery were localized up to 800 nm higher than emitters near the center of the field of view because spatially variant aberrations, especially defocus, were not taken into account. We recapitulated this bowl shape using the Easy-DHPSF algorithm to analyze simulated images of bright emitters in the same *z*-plane (depth = 50 nm) imaged through an aberrated imaging system by using the experimentally determined aberration coefficient profiles [Fig. 7(b)]. To simulate these images, isotropic emitters were placed adjacent to the refractive index boundary (*z* = 50 nm) and the nominal focal plane was made coincident with the refractive index boundary (*d* = 0).When using a global look-up table, spatially variant aberrations were not taken into account and the *z*-error, defined as the average z-positions of *N* = 5 individual estimates relative to the ground truth, mapped out a bowl shape similar to the experimental measurement. When we refit the same data using many local look-up tables [36], the *z*-errors were reduced by roughly one order of magnitude, but they still displayed a spatial dependence [Fig. 7(c)]. However, this improvement was dependent on the emitter positions relative to the nominal focal plane. When the nominal focal plane was coincident with the refractive index boundary, most of the estimated z-errors were within the [-60, 80] nm range.

By contrast, when the nominal focal plane was placed 500 nm above the coverslip (*d* = −500 nm), as would be done when localizing emitters inside a 1000 nm sized bacterial cell, the *z*-error range increased to [0, 200] nm [Fig. 7(d)]. Importantly, in both the *d* = 0 nm and *d* = −500 nm cases, the *z*-error remained spatially variant. By comparison, the spatially variant vectorial PSF model estimated the *z*-positions without any spatial dependence (absolute mean z-error = 1.4 nm), because the *z*-errors were determined solely by the signal-to-background-limited localization precision. When the separation of the emitters and the nominal focal plane was random within a range of ± 500 nm, the benefit of using many local lookup table was not as pronounced as for a constant separation distance [Figs. 7(c) and 7(d)]. Such a situation is encountered when single-molecules randomly dispersed in a 1000 nm sized bacterial cells are localized relative to the nominal focal plane stably positioned at *d* = −500 nm [12]. The *z*-errors obtained by double Gaussian fitting using local look-up tables showed spatially dependent *z*-errors in the ranges of [0, 300] nm [Fig. 7(e)]. We found that the Easy-DHPSF algorithm often converged to incorrect estimates, especially when the simulated PSFs were severely aberrated. Incorrect estimates obtain frequently even when the double Gaussian model was heavily constrained and independent of whether global or local look-up tables were used. While grossly incorrect localizations, such as those outside the calibrated DHPSF *z*-range, can be filtered out during data post-processing, such filtering steps have the undesired outcome of reducing data acquisition throughput.

Finally, we evaluated the accuracy and precision of different methods as a function of the signal-to-background ratio. Noised images of a single bead (*z* = 50 nm; *d* = 0) were simulated with different signal-to-noise ratios at an off-center position in the field of view (70 µm lateral shift) and analyzed to extract the emitter *x*,*y-* and *z-*positions. The background was held constant at 10 photons per pixel. At different signal-to-background ratios, the accuracy of the double Gaussian model can be improved by using a local look-up table, but the precison remained unchanged. For all signal-to-background ratios evaluated, the vectorial PSF model yielded better accuracy and precision than the double Guassian model [Fig. 7(f)]. We conclude that the use spatially variant PSF models is the optimal solution for single-molecule imaging in ultrawide fields-of-view, where optical aberration can vary substantially between different locations.

## 4. Conclusions

In this work, we introduce a spatially variant vectorial PSF model to accurately and precisely determine the *x*, *y*- and *z*-positions of individual emitters in an ultrawide field of view. The spatially variant aberrations are calibrated using bead images obtained at different positions in the field of view and at different nominal focal plane positions. This procedure quantifies the optical aberrations of the imaging system as linear combinations of Zernike polynomials. The coefficients of these linear combinations vary throughout the field of view and enable the generation of realistic PSF images for experimental data fitting.

We show that current implementation of the spatially variant vectorial PSF model achieves better precision and accuracy than simple analytical models. Careful measurements in two color channels also showed that *z*-colocalization accuracies of less than 25 nm can be achieved. The spatially variant vectorial PSF model presented here can be applied to model any (engineered)

PSF corresponding to isotropic emitters, such as beads or rotating single-molecule fluorophores, or anisotropic emitters, such as fixed dipoles, in a refractive index mismatched medium.

Limitations of the approach include the inability to calibrate the thickness variation of microscope coverslips *in situ*. Unknown undulations in the refractive index boundary limit the accuracy in Zernike coefficients estimation and thus the absolute *z*-localization accuracy to tens of nanometers, as quantified in Figs. 3(c) and 3(d). Possible improvements of the vectorial PSF model could include performing more extensive calibration measurements, incorporating additional Zernike modes beyond *j* = 15, and modeling refractive index variations in the specimen. Incorporating the finite bandwidth of the fluorescence emission spectrum into the PSF model may also help to further improve the agreement between simulated and experimental PSFs. However, such improvements add to the already substantial computational cost of our model (~2 min/localization when run on a single core of a workstation computer). Data processing time could be shortened by using Graphics Processing Units (GPUs), as recently demonstrated for interpolated PSF models [44]. Alternatively, spatially variant interpolated PSFs could be generated based on the method presented here. Such an approach would result in improved localization accuracy, because locally aberrated PSF models could be used for interpolation by computationally moving the emitter to different *z*-positions instead of experimentally scanning the nominal focal plane through emitters at constant (but unknown) *z*-position.

## Funding

University of Virginia; Hartwell Foundation.

## Acknowledgments

This work was supported by startup funds from the University of Virginia. C. R. was supported by a postdoctoral fellowship from the Hartwell Foundation.

## References

**1. **S. J. Sahl, S. W. Hell, and S. Jakobs, “Fluorescence nanoscopy in cell biology,” Nat. Rev. Mol. Cell Biol. **18**(11), 685–701 (2017). [CrossRef] [PubMed]

**2. **W. E. Moerner, “Microscopy beyond the diffraction limit using actively controlled single molecules,” J. Microsc. **246**(3), 213–220 (2012). [CrossRef] [PubMed]

**3. **B. Huang, M. Bates, and X. Zhuang, “Super-resolution fluorescence microscopy,” Annu. Rev. Biochem. **78**(1), 993–1016 (2009). [CrossRef] [PubMed]

**4. **E. Betzig, G. H. Patterson, R. Sougrat, O. W. Lindwasser, S. Olenych, J. S. Bonifacino, M. W. Davidson, J. Lippincott-Schwartz, and H. F. Hess, “Imaging intracellular fluorescent proteins at nanometer resolution,” Science **313**(5793), 1642–1645 (2006). [CrossRef] [PubMed]

**5. **S. T. Hess, T. P. Girirajan, and M. D. Mason, “Ultra-high resolution imaging by fluorescence photoactivation localization microscopy,” Biophys. J. **91**(11), 4258–4272 (2006). [CrossRef] [PubMed]

**6. **M. J. Rust, M. Bates, and X. Zhuang, “Sub-diffraction-limit imaging by stochastic optical reconstruction microscopy (STORM),” Nat. Methods **3**(10), 793–796 (2006). [CrossRef] [PubMed]

**7. **A. Gahlmann and W. E. Moerner, “Exploring bacterial cell biology with single-molecule tracking and super-resolution imaging,” Nat. Rev. Microbiol. **12**(1), 9–22 (2014). [CrossRef] [PubMed]

**8. **P. Rubin-Delanchy, G. L. Burn, J. Griffié, D. J. Williamson, N. A. Heard, A. P. Cope, and D. M. Owen, “Bayesian cluster identification in single-molecule localization microscopy data,” Nat. Methods **12**(11), 1072–1076 (2015). [CrossRef] [PubMed]

**9. **P. Sengupta, T. Jovanovic-Talisman, D. Skoko, M. Renz, S. L. Veatch, and J. Lippincott-Schwartz, “Probing protein heterogeneity in the plasma membrane using PALM and pair correlation analysis,” Nat. Methods **8**(11), 969–975 (2011). [CrossRef] [PubMed]

**10. **S. Manley, J. M. Gillette, G. H. Patterson, H. Shroff, H. F. Hess, E. Betzig, and J. Lippincott-Schwartz, “High-density mapping of single-molecule trajectories with photoactivated localization microscopy,” Nat. Methods **5**(2), 155–157 (2008). [CrossRef] [PubMed]

**11. **M. Mustafi and J. C. Weisshaar, “Simultaneous Binding of Multiple EF-Tu Copies to Translating Ribosomes in Live Escherichia coli,” MBio **9**(1), e02143 (2018). [CrossRef] [PubMed]

**12. **J. M. Rocha, C. J. Richardson, M. Zhang, C. M. Darch, E. Cai, A. Diepold, and A. Gahlmann, “Single-molecule tracking in live Yersinia enterocolitica reveals distinct cytosolic complexes of injectisome subunits,” Integr. Biol. **10**(9), 502–515 (2018). [CrossRef] [PubMed]

**13. **A. G. Santiago, T. Y. Chen, L. A. Genova, W. Jung, A. M. George Thompson, M. M. McEvoy, and P. Chen, “Adaptor protein mediates dynamic pump assembly for bacterial metal efflux,” Proc. Natl. Acad. Sci. U.S.A. **114**(26), 6694–6699 (2017). [CrossRef] [PubMed]

**14. **J. Engelhardt, J. Keller, P. Hoyer, M. Reuss, T. Staudt, and S. W. Hell, “Molecular orientation affects localization accuracy in superresolution far-field fluorescence microscopy,” Nano Lett. **11**(1), 209–213 (2011). [CrossRef] [PubMed]

**15. **A. Small and S. Stahlheber, “Fluorophore localization algorithms for super-resolution microscopy,” Nat. Methods **11**(3), 267–279 (2014). [CrossRef] [PubMed]

**16. **C. Smith, M. Huisman, M. Siemons, D. Grünwald, and S. Stallinga, “Simultaneous measurement of emission color and 3D position of single molecules,” Opt. Express **24**(5), 4996–5013 (2016). [CrossRef] [PubMed]

**17. **Y. Shechtman, L. E. Weiss, A. S. Backer, M. Y. Lee, and W. E. Moerner, “Multicolour localization microscopy by point-spread-function engineering,” Nat. Photonics **10**(9), 590–594 (2016). [CrossRef] [PubMed]

**18. **M. Siemons, C. N. Hulleman, R. O. Thorsen, C. S. Smith, and S. Stallinga, “High precision wavefront control in point spread function engineering for single emitter localization,” Opt. Express **26**(7), 8397–8416 (2018). [CrossRef] [PubMed]

**19. **M. P. Backlund, M. D. Lew, A. S. Backer, S. J. Sahl, G. Grover, A. Agrawal, R. Piestun, and W. E. Moerner, “Simultaneous, accurate measurement of the 3D position and orientation of single molecules,” Proc. Natl. Acad. Sci. U.S.A. **109**(47), 19087–19092 (2012). [CrossRef] [PubMed]

**20. **A. S. Backer, M. Y. Lee, and W. E. Moerner, “Enhanced DNA imaging using super-resolution microscopy and simultaneous single-molecule orientation measurements,” Optica **3**(6), 659–666 (2016). [CrossRef] [PubMed]

**21. **O. Zhang, J. Lu, T. Ding, and M. D. Lew, “Imaging the three-dimensional orientation and rotational mobility of fluorescent emitters using the Tri-spot point spread function,” Appl. Phys. Lett. **113**(3), 031103 (2018). [CrossRef] [PubMed]

**22. **H. P. Kao and A. S. Verkman, “Tracking of single fluorescent particles in three dimensions: use of cylindrical optics to encode particle position,” Biophys. J. **67**(3), 1291–1300 (1994). [CrossRef] [PubMed]

**23. **B. Huang, W. Wang, M. Bates, and X. Zhuang, “Three-dimensional super-resolution imaging by stochastic optical reconstruction microscopy,” Science **319**(5864), 810–813 (2008). [CrossRef] [PubMed]

**24. **S. R. P. Pavani and R. Piestun, “Three dimensional tracking of fluorescent microparticles using a photon-limited double-helix response system,” Opt. Express **16**(26), 22048–22057 (2008). [CrossRef] [PubMed]

**25. **S. R. Pavani, M. A. Thompson, J. S. Biteen, S. J. Lord, N. Liu, R. J. Twieg, R. Piestun, and W. E. Moerner, “Three-dimensional, single-molecule fluorescence imaging beyond the diffraction limit by using a double-helix point spread function,” Proc. Natl. Acad. Sci. U.S.A. **106**(9), 2995–2999 (2009). [CrossRef] [PubMed]

**26. **M. D. Lew, S. F. Lee, M. Badieirostami, and W. E. Moerner, “Corkscrew point spread function for far-field three-dimensional nanoscale localization of pointlike objects,” Opt. Lett. **36**(2), 202–204 (2011). [CrossRef] [PubMed]

**27. **D. Baddeley, M. B. Cannell, and C. Soeller, “Three-dimensional sub-100 nm super-resolution imaging of biological samples using a phase ramp in the objective pupil,” Nano Res. **4**(6), 589–598 (2011). [CrossRef]

**28. **S. Jia, J. C. Vaughan, and X. Zhuang, “Isotropic 3D Super-resolution Imaging with a Self-bending Point Spread Function,” Nat. Photonics **8**, 302–306 (2014). [CrossRef] [PubMed]

**29. **Y. Shechtman, S. J. Sahl, A. S. Backer, and W. E. Moerner, “Optimal point spread function design for 3D imaging,” Phys. Rev. Lett. **113**(13), 133902 (2014). [CrossRef] [PubMed]

**30. **Y. Shechtman, L. E. Weiss, A. S. Backer, S. J. Sahl, and W. E. Moerner, “Precise Three-Dimensional Scan-Free Multiple-Particle Tracking over Large Axial Ranges with Tetrapod Point Spread Functions,” Nano Lett. **15**(6), 4194–4199 (2015). [CrossRef] [PubMed]

**31. **L. Carlini, S. J. Holden, K. M. Douglass, and S. Manley, “Correction of a Depth-Dependent Lateral Distortion in 3D Super-Resolution Imaging,” PLoS One **10**(11), e0142949 (2015). [CrossRef] [PubMed]

**32. **P. Almada, S. Culley, and R. Henriques, “PALM and STORM: Into large fields and high-throughput microscopy with sCMOS detectors,” Methods **88**, 109–121 (2015). [CrossRef] [PubMed]

**33. **K. M. Douglass, C. Sieben, A. Archetti, A. Lambert, and S. Manley, “Super-resolution imaging of multiple cells by optimised flat-field epi-illumination,” Nat. Photonics **10**(11), 705–708 (2016). [CrossRef] [PubMed]

**34. **Z. Zhao, B. Xin, L. Li, and Z. L. Huang, “High-power homogeneous illumination for super-resolution localization microscopy with large field-of-view,” Opt. Express **25**(12), 13382–13395 (2017). [CrossRef] [PubMed]

**35. **F. Huang, T. M. Hartwich, F. E. Rivera-Molina, Y. Lin, W. C. Duim, J. J. Long, P. D. Uchil, J. R. Myers, M. A. Baird, W. Mothes, M. W. Davidson, D. Toomre, and J. Bewersdorf, “Video-rate nanoscopy using sCMOS camera-specific single-molecule localization algorithms,” Nat. Methods **10**(7), 653–658 (2013). [CrossRef] [PubMed]

**36. **A. von Diezmann, M. Y. Lee, M. D. Lew, and W. E. Moerner, “Correcting field-dependent aberrations with nanoscale accuracy in three-dimensional single-molecule localization microscopy,” Optica **2**(11), 985–993 (2015). [CrossRef] [PubMed]

**37. **A. W. Lohmann and D. P. Paris, “Space-variant image formation,” J. Opt. Soc. Am. **55**(8), 1007–1013 (1965). [CrossRef]

**38. **M. Řeřábek and P. Patá, “The space variant PSF for deconvolution of wide-field astronomical images,” Proc. SPIE **7015**, 70152G (2008). [CrossRef]

**39. **M. Řeřábek, P. Patá, and K. Fliegel, “Enhancement of the accuracy of the astronomical measurements carried on the wide-field astronomical image data,” Proc. SPIE **8135**, 81351M (2011).

**40. **S. F. Gibson and F. Lanni, “Experimental test of an analytical model of aberration in an oil-immersion objective lens used in three-dimensional light microscopy,” J. Opt. Soc. Am. A **9**(1), 154–166 (1992). [CrossRef] [PubMed]

**41. **S. Stallinga and B. Rieger, “Position and orientation estimation of fixed dipole emitters using an effective Hermite point spread function model,” Opt. Express **20**(6), 5896–5921 (2012). [CrossRef] [PubMed]

**42. **A. Tahmasbi, E. S. Ward, and R. J. Ober, “Determination of localization accuracy based on experimentally acquired image sets: applications to single molecule microscopy,” Opt. Express **23**(6), 7630–7652 (2015). [CrossRef] [PubMed]

**43. **H. P. Babcock and X. Zhuang, “Analyzing Single Molecule Localization Microscopy Data Using Cubic Splines,” Sci. Rep. **7**(1), 552 (2017). [CrossRef] [PubMed]

**44. **Y. Li, M. Mund, P. Hoess, J. Deschamps, U. Matti, B. Nijmeijer, V. J. Sabinina, J. Ellenberg, I. Schoen, and J. Ries, “Real-time 3D single-molecule localization using experimental point spread functions,” Nat. Methods **15**(5), 367–369 (2018). [CrossRef] [PubMed]

**45. **H. Kirshner, C. Vonesch, and M. Unser, “Can localization microscopy benefit from approximation theory?” in 2013 IEEE 10th International Symposium on Biomedical Imaging (IEEE, 2013), pp. 588–591. [CrossRef]

**46. **P. N. Petrov, Y. Shechtman, and W. E. Moerner, “Measurement-based estimation of global pupil functions in 3D localization microscopy,” Opt. Express **25**(7), 7945–7959 (2017). [CrossRef] [PubMed]

**47. **A. S. Backer and W. E. Moerner, “Extending single-molecule microscopy using optical Fourier processing,” J. Phys. Chem. B **118**(28), 8313–8329 (2014). [CrossRef] [PubMed]

**48. **D. Axelrod, “Fluorescence excitation and imaging of single molecules near dielectric-coated and bare surfaces: a theoretical study,” J. Microsc. **247**(2), 147–160 (2012). [CrossRef] [PubMed]

**49. **M. A. Thompson, M. D. Lew, M. Badieirostami, and W. E. Moerner, “Localizing and tracking single nanoscale emitters in three dimensions with high spatiotemporal resolution using a double-helix point spread function,” Nano Lett. **10**(1), 211–218 (2010). [CrossRef] [PubMed]

**50. **A. Gahlmann, J. L. Ptacin, G. Grover, S. Quirin, A. R. von Diezmann, M. K. Lee, M. P. Backlund, L. Shapiro, R. Piestun, and W. E. Moerner, “Quantitative multicolor subdiffraction imaging of bacterial protein ultrastructures in three dimensions,” Nano Lett. **13**(3), 987–993 (2013). [CrossRef] [PubMed]

**51. **J. L. Ptacin, A. Gahlmann, G. R. Bowman, A. M. Perez, A. R. von Diezmann, M. R. Eckart, W. E. Moerner, and L. Shapiro, “Bacterial scaffold directs pole-specific centromere segregation,” Proc. Natl. Acad. Sci. U.S.A. **111**(19), E2046–E2055 (2014). [CrossRef] [PubMed]

**52. **M. D. Lew, A. R. S. von Diezmann, and W. E. Moerner, “Easy-DHPSF open-source software for three-dimensional localization of single molecules with precision beyond the optical diffraction limit,” Protoc. Exch. (2013), doi:. [CrossRef]

**53. **B. Richards and E. Wolf, “Electromagnetic diffraction in optical systems, II. Structure of the image field in an aplanatic system,” Proc. R. Soc. Lond. A Math. Phys. Sci. **253**, 358–379 (1959).

**54. **G. Zheng, X. Ou, R. Horstmeyer, and C. Yang, “Characterization of spatially varying aberrations for wide field-of-view microscopy,” Opt. Express **21**(13), 15131–15143 (2013). [CrossRef] [PubMed]