In recent years three-dimensional (3D) super-resolution fluorescence imaging by single-molecule localization (localization microscopy) has gained considerable interest because of its simple implementation and high optical resolution. Astigmatic and biplane imaging are experimentally simple methods to engineer a 3D-specific point spread function (PSF), but existing evaluation methods have proven problematic in practical application. Here we introduce the use of cubic B-splines to model the relationship of axial position and PSF width in the above mentioned approaches and compare the performance with existing methods. We show that cubic B-splines are the first method that can combine precision, accuracy and simplicity.
© 2014 Optical Society of America
Since the introduction of localization microscopy methods such as photoactivated localization microscopy (PALM) [1, 2], fluorescence photoactivation localization microscopy (FPALM) , stochastic optical reconstruction microscopy (STORM) [4,5] and direct STORM (dSTORM) [6, 7], several techniques have been developed to enable precise 3D localization of single emitters. Most commonly, 3D information is obtained by imaging the lateral plane and inferring the axial coordinate from the lateral width of the point spread function (PSF). Unfortunately, the natural PSF of a wide-field microscope yields little axial information due to two reasons: It is axially symmetric, and the PSF is insensitive to axial changes close to the focus. For both reasons, the PSF has to be modified to contain axial information. With a modified PSF, resolutions of 15–20 nm laterally and 40–80 nm axially are achieved [8,9] (values expressed as full width at half maximum (FWHM)). Using interferometric PALM, axial resolutions of the order of 1–10 nm can be achieved, even though with restrictions to smaller axial ranges [10,11].
Three methods have been used to engineer the PSF to be more responsive to axial variations: Lobe splitting, biplane imaging and astigmatic imaging. For lobe splitting approaches like the double-helix PSF  and phase ramp imaging localization microscopy (PRILM) , a PSF is split into two lobes whose relative position encodes the axial position. This method, however, is limited to very sparse signals because the algorithms used require well separated emitters for accurate 3D localization.
For the biplane approach , the emission beam is split equally and projected onto two cameras or two distinct areas of one camera chip. The focal planes in each path are calibrated to be a few hundred nanometers apart, and recorded simultaneously. In this way, two sections of the PSF with known axial offsets are provided. The second plane both breaks the symmetry of the PSF and ensures that at least one of the planes images an unfocused signal.
The astigmatic approach [8,14] inserts a cylindrical lens into the optical detection path, which leads to slightly different focal planes in the x and y directions resulting in elliptically shaped PSFs whose ellipticity changes with respect to z.
Biplane and astigmatic 3D imaging have the distinct advantage over other PSF engineering approaches, that the Gaussian approximation to the PSF (Eq. (1)) holds for a considerable axial range. Existing 2D evaluation programs and their extensive theory [15–18] can be used. If the Gaussian approximation is dropped, model-free approaches [13, 19] with considerable computational cost and complexity must be used. For single-molecule localization of organic fluorophores, the Gaussian approximation has been used with good results in 3D . Computationally, biplane 3D and astigmatic 3D are very similar, because in biplane 3D the width of the two round spots in both planes takes the role of the two widths of the elliptical spot in astigmatic 3D. The PSF width of astigmatic data tends to be harder to model because of the additional lens. Therefore, we investigate only astigmatic 3D in this article, but expect all conclusions to be applicable also to biplane 3D.
Using the Gaussian approximation, the lateral PSF widths σx and σy encode the axial position z. We know of three algorithms to decode z: Sigma-difference look-up, quadratic approximation and quartic approximation. Henriques et al.  estimated σx and σy from a calibration data set with known axial positions, and created a lookup table for z indexed by σx − σy. The z position of an unknown fluorophore was determined by fitting its PSF with a Gaussian with free widths, and then looking up the observed σx − σy in the table.
Holtzer et al.  used a physically derived model, assuming that the PSF broadening is parabolic. Equation (2) resembles this model, with σ0,(x,y) giving the standard deviation of the PSF in the planes of sharpest x and y, z0,(x,y) giving the focal point for each plane, and a(x,y) being an optics-dependent parameter. The expressions for σx and σy can be inserted into the PSF, and the z position becomes a regular parameter of G and can be fitted like the other parameters.
In practice, the quadratic model does not fit well to typical data and thus Huang et al.  extended the model with terms up to the fourth order, arguing that those are needed to account for the imperfection of optics (Eq. (3)). Δσi,(x,y) represents the effective focus depth for the polynomial term (cf. ). Since the parabolic model is a special case of the quartic model with Δσi,(x,y) = ∞ for all i ≠ 2, we refer to both models as the polynomial model.
From a computational point of view, all three approaches attempt an approximation of the functional relationship between z and σx,y from noisy calibration data. The sigma-difference method only partially achieves that goal, because it models σx − σy instead of the individual σ components and cannot be used to replace the σ terms in the PSF with functions of z. This forces the researcher to fit σx and σy individually and afterwards apply the sigma-difference algorithm. The additional and unphysical degree of freedom in the model has been shown in 2D to cause a loss in precision . The quadratic and quartic methods, on the other hand, use polynomial functions for the approximation, which are bad approximators because all their parameters are global. Small errors in the calibration data perturb the entire approximated function. Such errors are common at the limit of the calibration interval, where the PSF is smeared out over many pixels and measurements of σ are highly uncertain.
Cubic basis splines  (cubic B-splines) are a well-established approximation technique. It avoids these problems and is computationally simpler than polynomial fitting. In this article, we use our rapid, accurate program implementing direct stochastic optical reconstruction microscopy (rapidSTORM)  software to evaluate and compare the fitting performances of the above mentioned techniques and will show that the evaluation with cubic B-splines performs equally well or better concerning localization precision.
2. Methods and instrumentation
2.1. Optical setup
The presented data was acquired on two similar setups which are schematically depicted in Fig. 1. For setup 1, the laser light of a CUBE 640 nm 100 mW (Coherent, USA) enters the IX71 (Olympus, Japan) inverted microscope via the back entrance and is reflected by a quad band dichroic mirror (FF410-504-582-669-Di01, Semrock, USA) to a 60× NA = 1.45 oil-immersion objective (PLAPON 60XOTIRFM, Olympus, Japan). It is thereby focused by a two-lens-system on the back focal plane of the objective in order to obtain wide-field illumination. The fluorescence signal passes the dichroic mirror and is imaged by lenses L1 (f = 80 mm, Qioptic, USA) and L2 (f = 200 mm, Qioptic, USA) on an EMCCD (DU-897 (16 μm p. pixel), Andor, UK) with mounted fluorescence emission filter (ET700/75, Chroma, USA). A cylindrical lens CL with focal length f = 10 m (SCX-25.4-5000.0-C, CVI Melles Griot, USA) is placed in a section of the detection path where fluorescence light is parallel and causes distinct focal planes for x- and y-direction (FPx and FPy). We chose the long focal length to best fit to the other components of the optical system.
Setup 2 is built around a Zeiss Axio Observer.D1 (Carl Zeiss Microscopy, Germany) equipped with a 63× NA = 1.15 water-immersion objective (LD C-Apochromat 63x/1.15 W Corr M27, Carl Zeiss Microscopy, Germany). Since the detection consists of a sCMOS camera (Neo 5.5 (6.5 μm p. pixel), Andor, UK) with smaller pixels than on an EMCCD chip, no further magnification is needed in order to obtain a final pixel-size of about 100 nm/px and the fluorescence is directly imaged by the tube-lens onto the chip. A f = 100 mm cylindrical lens (LJ1567RM-A, Thorlabs, USA), placed ∼25 mm in front of the camera front plate, introduces astigmatism (cf. ). In this position the fluorescence light is convergent and the amount of astigmatic distortion can be tuned by the distance between the lens and the (undistorted) focal plane.
2.2. Theoretical background for cubic B-spline evaluation
Cubic B-splines are computed by dividing the data into segments between knots (black circles in Fig. 2) and recursively constructing a basis function over each segment . Each basis function si covers four intervals. The cubic B-spline S is the weighted sum of all basis functions si with weights ci (Eq. (4)). Rewriting the basis functions in matrix form (Eq. (5)) with its transposed XT and inverse X−1, a least-squares fit (Eq. (6)) extracts the weights from the calibration data set of measured width σi at known z positions zi. For a more detailed description, please refer to the rapidSTORM manual .
2.3. Measurement procedure
In order to record a reference sample we used an objective z piezo scanner (PIfoc, Physik Instrumente, Germany) and 100 nm TetraSpeck fluorescent microspheres (T-7279, life technologies, USA) sparsely dispersed either on a cover slip or in a hydrogel-layer with a refractive index of about 1.34 (Matrigel, BD Biosciences, USA). The setup was chosen to closely match the calibration sample to later experiments, i.e., having the specimen close to the cover slip or in aqueous environment. Measurements using oil immersion far away from the cover slip are prone to aberrations. Thus for deep volume imaging the use of water immersion is advantageous and we consider these two configurations the most reasonable. The specimen is scanned over several micrometers by applying a triangular voltage to the piezo. The fluorescence is imaged by an EMCCD or sCMOS camera.
With the measured data we computed a ground truth z position for each frame by choosing z0 = 0 so that σx(z0) = σy(z0) and computing the true z position for further frames from the piezo parameters. The calibration bead’s image was localized in each frame with the standard rapidSTORM fitting procedure , i.e., fitting a Gaussian PSF model with free PSF width parameters to the data using Levenberg-Marquardt least squares fitting. We used rapidSTORM version 3.3.
The fitting procedure yields z, σx and σy values for each bead image. The points make up two data sets, one with σx as a function of z and one with σy as a function of z. We then used different functional approximations of the data points to gain smooth, continuous functions describing the data points. In the following, we will refer to the set of fitted pairs of z to σx and σy as “sigma plot”, and to its functional approximation as “sigma curve”.
2.3.1. Polynomial method
For the polynomial method, we cut an axial region that could be approximated with a sigma curve with positive curvature, and removed inconsistent points, i.e., outliers that obviously do not resemble the true PSF. Since Eq. (3) can be expanded to a polynomial of fourth order, we used linear least-squares fitting of the squared PSF widths to determine its coefficients. The linear least-squares fit gives five parameters bi for a normalized polynomial around the arbitrarily chosen origin of the z coordinates. We determined the parameter z0 first by determining the global minimum of the fitted polynomial via the roots of the derivative. Afterwards, we subtracted z0 from all z coordinates, fitted a normalized polynomial again, and used simple parameter matching and algebra to compute and σ0 = b0. Since we chose the global minimum for z0, the linear term b1 was always 0.
If the sigma plot and sigma curve were obviously deviating, we varied the choice of the axial region. The fitted parameter values were used as input parameters to rapidSTORM’s polynomial PSF model.
2.3.2. Nelder-Mead optimization of the polynomial parameters
For optimizing the result of the polynomial method with Nelder-Mead optimization, we used the GNU Scientific Library  implementation of the Nelder-Mead simplex method . The objective function was the squared distance between the piezo-determined and the rapidSTORM-fitted z coordinates for all beads in a 1.5 μm wide axial section centered on the equifocal z coordinate. We chose the initial simplex size relative to the obtained parameter values, with a factor between 10% and 50%.
When the objective function failed to improve, or we were unsatisfied with the obtained accuracy, we varied the simplex size factor or recombined the parameters gained from different runs. When this also failed, we used the polynomial method on smaller axial regions and used them as inputs for the Nelder-Mead optimization on larger regions. Finally, we chose the parameter set with the best value for the objective function.
2.3.3. Sigma-difference method
We applied the sigma difference method  by computing σx − σy. We smoothed the points in the sigma plot with a Gaussian low-pass filter of 50 nm width. Then, we cut out the axial section around the equifocal point where the values were monotonous. Afterwards, we looked up the nearest value of σx − σy in the smoothed sigma plot for every point in the unsmoothed sigma plot, with linear interpolation, and used the smoothed point’s z coordinate as the localized z coordinate.
2.3.4. Cubic B-spline method
For the cubic B-spline method, we cut out a suitable axial region as described for the polynomial model. We interpolated the PSF widths with a cubic basis spline with ten knots using the GNU Scientific Library implementation  instead of a polynomial model. We stored the knot coordinates (z and σ) in a data file. When fitting unknown data, we recovered the spline from the data file using a rapidSTORM-internal implementation of splines following McKinley and Levine . Since a cubic B-spline is a smooth analytic function, we substituted σx and σy in the PSF by the spline functions and used the z coordinate as a PSF input parameter.
2.4. Comparison with the Nikon N-STORM
Nikon offers 3D super-resolution imaging with their N-STORM setup using a built in auto-calibration routine (NIS-Elements 4.10 software with N-STORM plugin). The main optical components are an EMCCD (DU-897 (16 μm p. pixel), Andor, UK) and a 100× NA = 1.49 oil-immersion objective (Apo TIRF 100X 1.49 Oil, Nikon, Japan). We were interested in the performance of that black-box routine, which supposedly follows the polynomial method, and wanted to compare it with rapidSTORM. Therefore we took calibration data with the N-STORM, ran the autocalibration and also fed the data to the cubic B-spline-algorithm. In order to rearrange the N-STORM data to a format that rapidSTORM can interpret, we used the custom software “N-Storm to RapidStorm converter v. 0.4”, kindly provided by S. Malkusch. Again we computed the offset between known and fitted z position. We assumed that the N-STORM routine first takes twenty one frames at the initial plane, then jumps to z = −800 nm, and performs a z ramp with 10 nm / frame up to z = +800 nm. Finally, the piezo jumps back to the initial position for another twenty frames.
2.5. Quantifying computational robustness
The computational robustness measures the likelihood of erroneous application of a 3D method. Since the computational robustness is not absolutely quantifiable, we ranked the critical steps, i.e., steps with error sources, of all methods and took the sum of the ranks as a complexity score.
All employed methods start their calibration with a sigma sample, i.e., a calibration data set that gives σx and σy for known z positions. We therefore ignore the common critical steps in the procedure that obtains the sigma sample.
3.1. Axial localization performance
Cubic B-splines accurately model z-σ relationships for both water- and oil-immersed objectives (Fig. 3). Real measurements, especially on oil objectives with TetraSpecks near the phase boundary, commonly show asymmetric and irregular behavior that is hard to model with polynomial models. We found that computing cubic B-splines on different axial ranges resulted in very similar splines, while variations of the z-range could radically impact polynomial parameters.
The axial re-localization performance of the different width determination methods for the Gaussian PSF is shown in Fig. 4. The less bright Sample 2 showed photon counts comparable to single switching events of organic fluorophores. It has an overall lower localization precision due to a poor signal-to-noise ratio. The apparent higher variation in precision and accuracy is due to a lower axial sampling density, i.e., fewer localizations are averaged for each interval.
The polynomial method shows high precision, i.e., a low axial FWHM of the localization distribution, and overall reasonable accuracy, i.e., a close agreement of the true and the average measured z position. However, we found the calibration process for Sample 1 to be highly fragile. Initially, we chose to include an additional 100 nm of data points on the left, steeper slope. The inclusion triggered strong oscillations in the fitted polynomial and prevented accurate localization. While the oscillations are counterintuitive for a functional approximation of a smooth curve, they are a known shortcoming of polynomial approximations .
The fit accuracy was improved slightly by optimizing the polynomial model parameters with the Nelder-Mead algorithm. However, computing time on the scale of hours was consumed due to the repetitive re-fitting of the whole calibration stack with slightly modified parameter values, and several manual interventions were necessary. We suspect that the complex z-σ relationship of Sample 1 can not be accurately represented by a polynomial.
The sigma-difference method reached considerably better accuracy. At the same time, precision suffers noticeably because of the overparameterization of the fit problem with two free parameters (σx and σy) for one physical degree of freedom (z). This overparameterization is known to cause lateral localization imprecision  and we reason that the same effect applies to axial localization.
The cubic B-spline approach combines both the precision of the polynomial model and the accuracy of the sigma-difference lookup table. In our view, the precision stems from the explicit modeling of the z coordinate in the PSF model, like in the polynomial model. The accuracy is due to the model-free description of the z-σ relationship, like in the sigma-difference method. Furthermore, the cubic B-spline routine proves to be easy and robust in use. Therefore, it is ideally suited also for non-experts with relatively little knowledge about 3D fitting.
Figure 5 shows the complexity sources we identified in the 3D inference methods. Since a high complexity precludes the practical application of a method in the lab, complexity is a fundamental benchmark. As discussed before, choosing the fit range for the polynomial method proved surprisingly critical, and too large and too small choices of the axial fit region were within 200 nm. We reason that because all polynomial parameters are influenced by all points in the fit region, too large fit regions put undue weight on the unreliable points far away from the focus. On oil objectives, we saw asymmetric z-σ relationships, and we had to choose asymmetric fit regions. The polynomial method often failed to fit the steeper side, most likely because more points were available on the gentler slope. For too small fit regions, the higher-order terms showed very high uncertainties, and extrapolation behavior was poor.
3.3. Comparison with a commercially available setup
For better classification of our results we performed 3D localization experiments on a Nikon N-STORM system, a commercially available 3D localization microscope. We have no a priori knowledge about the software algorithms used for 3D localization and wanted to assess the performance of this out-of-the-box setup compared to the rapidSTORM cubic B-spline interpolation. We were especially interested in whether our routine proves to be an easy to use and robust method to check the output of Nikon’s highly automated approach.
In order to resemble realistic conditions, we did not try to tweak the calibration data acquisition and only took a not too dense region of TetraSpeck fluorescent microspheres deposited on a cover slip. We did not account actively for drift and did not apply filters for multi-bead conglomerates as for comparability we must assume that those were as well not filtered and accounted for by the N-STORM algorithm.
As depicted in Fig. 6, the z-localizations resulting from a rapidSTORM re-localization job as used throughout this work show a nice plateau around the focal plane. In this axial range of about 1 μm, the cubic B-spline accurately found the correct positions, while NIS-Elements, despite showing a small plateau around z = 0 nm, seems to be failing to find the correct coordinates. The localization precision obtained by rapidSTORM appears also slightly higher than that of the Nikon software.
We reason that indeed our method is able to compete with commercially available algorithms and, more importantly, is capable of helping the user to assess the calibration. Using rapidSTORM’s open-source cubic B-spline calibration, experimentalists can now evaluate the performance of their setup prior to running a real experiment.
We have described and tested a cubic B-spline approach for determining the relationship between the axial emitter position and the PSF width in the Gaussian PSF model. The cubic B-spline approach achieved both high precision and high accuracy by combining the strengths of two approaches: It incorporates the axial emitter position as a true fit parameter, like the highly precise polynomial width model [8, 21], and does not require an a priori model of the widths, like the highly accurate sigma-difference algorithm used in QuickPALM . Calibration for the cubic B-spline algorithm is trivial and well-understood, and a free and open reference implementation is available with the rapidSTORM software on our website at http://super-resolution.de.
We want to thank S. Malkusch for supplying us with his “N-Storm to RapidStorm converter” and P. Zessin for help with the N-STORM setup. This work was supported by the Biophotonics Initiative of the German Ministry of Research and Education (BMBF) (grant 13N11019). This publication was funded by the German Research Foundation (DFG) in the funding programme Open Access Publishing.
References and links
1. 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, 1642–1645 (2006). [CrossRef] [PubMed]
6. M. Heilemann, S. van de Linde, M. Schüttpelz, R. Kasper, B. Seefeldt, A. Mukherjee, P. Tinnefeld, and M. Sauer, “Subdiffraction-resolution fluorescence imaging with conventionalfluorescent probes,” Angew. Chem. Int. Ed. 47, 6172–6176 (2008). [CrossRef]
7. S. van de Linde, A. Löschberger, T. Klein, M. Heidbreder, S. Wolter, M. Heilemann, and M. Sauer, “Direct stochastic optical reconstruction microscopy with standard fluorescent probes,” Nat. Protoc. 6, 991–1009 (2011). [CrossRef] [PubMed]
9. M. F. Juette, T. J. Gould, M. D. Lessard, M. J. Mlodzianoski, B. S. Nagpure, B. T. Bennett, S. T. Hess, and J. Bewersdorf, “Three-dimensional sub-100 nm resolution fluorescence microscopy of thick samples,” Nat. Methods 5, 527–529 (2008). [CrossRef] [PubMed]
10. G. Shtengel, J. A. Galbraith, C. G. Galbraith, J. Lippincott-Schwartz, J. M. Gillette, S. Manley, R. Sougrat, C. M. Waterman, P. Kanchanawong, M. W. Davidson, R. D. Fetter, and H. F. Hess, “Interferometric fluorescent super-resolution microscopy resolves 3D cellular ultrastructure,” Proc. Natl. Acad. Sci. U. S. A. 106, 3125–3130 (2009). [CrossRef] [PubMed]
11. P. Kanchanawong, G. Shtengel, A. M. Pasapera, E. B. Ramko, M. W. Davidson, H. F. Hess, and C. M. Waterman, “Nanoscale architecture of integrin-based cell adhesions,” Nature 468, 580–584 (2010). [CrossRef] [PubMed]
12. S. R. P. 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, 2995–2999 (2009). [CrossRef] [PubMed]
13. D. Baddeley, M. 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, 589–598 (2011). [CrossRef]
17. K. I. Mortensen, L. S. Churchman, J. A. Spudich, and H. Flyvbjerg, “Optimized localization analysis for single-molecule tracking and super-resolution microscopy,” Nat. Methods 7, 377–381 (2010). [CrossRef] [PubMed]
19. A. G. York, A. Ghitani, A. Vaziri, M. W. Davidson, and H. Shroff, “Confined activation and subdiffractive localization enables whole-cell palm with genetically expressed probes,” Nat. Methods 8, 327–333 (2011). [CrossRef] [PubMed]
20. R. Henriques, M. Lelek, E. F. Fornasiero, F. Valtorta, C. Zimmer, and M. M. Mhlanga, “QuickPALM: 3D real-time photoactivation nanoscopy image processing in ImageJ,” Nat. Methods 7, 339–340 (2010). [CrossRef] [PubMed]
21. L. Holtzer, T. Meckel, and T. Schmidt, “Nanometric three-dimensional tracking of individual quantum dots in cells,” Appl. Phys. Lett. 90, 053902 (2007). [CrossRef]
22. S. Wolter, S. Proppert, S. Aufmkolk, A. Lampe, and T. Klein, “rapid STORM manual,” http://www.super-resolution.de/home/rapidstorm/ (2012).
23. S. Wolter, M. Schüttpelz, M. Tscherepanow, S. van de Linde, M. Heilemann, and M. Sauer, “Real-time computation of subdiffraction-resolution fluorescence images,” J. Microsc. 237, 12–22 (2010). [CrossRef] [PubMed]
24. C. de Boor, A Practical Guide to Splines, Applied Mathematical Sciences (Springer, 1978). [CrossRef]
25. S. Wolter, A. Löschberger, T. Holm, S. Aufmkolk, M.-C. Dabauvalle, S. van de Linde, and M. Sauer, “rapid STORM: accurate, fast open-source software for localization microscopy,” Nat. Methods 9, 1040–1041 (2012). [CrossRef] [PubMed]
26. M. Galassi, J. Davies, J. Theiler, B. Gough, G. Jungman, M. Booth, and F. Rossi, Gnu Scientific Library: Reference Manual (Network Theory Ltd., 2003).
27. J. A. Nelder and R. Mead, “A simplex method for function minimization,” Comput. J. 7, 308–313 (1965). [CrossRef]
28. S. McKinley and M. Levine, Cubic Spline Interpolation (College of the Redwoods, 1999).