Determination of localization accuracy based on experimentally acquired image sets : applications to single molecule microscopy

Fluorescence microscopy is a photon-limited imaging modal ity that allows the study of subcellular objects and processes w ith high specificity. The best possible accuracy (standard deviation) wit h hich an object of interest can be localized when imaged using a fluorescence microscope is typically calculated using the Cram ér-Rao lower bound, that is, the inverse of the Fisher information. However, the current approach fo r the calculation of the best possible localization accuracy relies on an anal ytic expression for the image of the object. This can pose practical challeng es since it is often difficult to find appropriate analytical models for the images of general objects. In this study, we instead develop an approach that d irectly uses an experimentally collected image set to calculate the best po ssible localization accuracy for a general subcellular object. In this approach , we fit splines, i.e. smoothly connected piecewise polynomials, to the expe rimentally collected image set to provide a continuous model of the obje ct, which can then be used for the calculation of the best possible localiz ation accuracy. Due to its practical importance, we investigate in detail th e application of the proposed approach in single molecule fluorescence micro scopy. In this case, the object of interest is a point source and, therefore , the acquired image set pertains to an experimental point spread function . © 2015 Optical Society of America OCIS codes: (100.2000) Digital image processing; (180.2520) Fluoresce nce microscopy; (030.5290) Photon statistics; (180.6900) Three-dimension al microscopy. References and links 1. W. E. Moerner and D. P. Fromm, “Methods of single-molecule fluo rescence spectroscopy and microscopy,” Rev. Sci. Instrum.74, 3597–3619 (2003). 2. A. Small and S. Stahlheber, “Fluorophore localization alg orithms for super-resolution microscopy,” Nat. Methods 11(3), 267–279 (2014). 3. S. Ram, P. Prabhat, J. Chao, E. S. Ward, and R. J. Ober, “High a ccuracy 3D quantum dot tracking with multifocal plane microscopy for the study of fast intracellular dynamics in live cells,” Biophys. J.95, 6025–6043 (2008). 4. R. J. Ober, A. Tahmasbi, S. Ram, Z. Lin, and E. S. Ward, “Quanti tative aspects of single molecule microscopy: Information-theoretic analysis of single-molecule data,” I EEE Signal Process. Mag. 32(1), 58–69 (2015). 5. R. J. Ober, S. Ram, and E. S. Ward, “Localization accuracy in single-molecule microscopy,” Biophys. J. 86, 1185–1200 (2004). 6. A. V. Abraham, S. Ram, J. Chao, E. S. Ward, and R. J. Ober, “Quan titative study of single molecule location estimation techniques,” Opt. Express 17, 23352–23373 (2009). 7. A. Tahmasbi, S. Ram, J. Chao, A. V. Abraham, F. W. Tang, E. S. War d, and R. J. Ober, “Designing the focal plane spacing for multifocal plane microscopy,” Opt. Express 22(14), 16706–16721 (2014). 8. S. M. Kay,Fundamentals of Statistical Signal Processing, Volume I: E stimation Theory(Prentice Hall, Upper Saddle River, NJ, USA, 1993). 9. D. L. Snyder and M. I. Miller,Random Point Processes in Time and Space (Springer Verlag, NY, USA, 1991), 2nd ed. 10. S. Ram, E. S. Ward, and R. J. Ober, “A stochastic analysis of performance limits for optical microscopes,” Multidim. Sys. Sig. Proc. 17, 27–57 (2006). 11. C. S. Smith, N. Joseph, B. Rieger, and K. A. Lidke, “Fast, si ngle-molecule localization that achieves theoretically minimum uncertainty,” Nat. Methods 7, 373–375 (2010). 12. M. Badieirostami, M. D. Lew, M. A. Thompson, and W. E. Moerne r, “Three-dimensional localization precision of the double-helix point spread function versus astigmatis m and biplane,” Appl. Phys. Lett. 97, 161103 (2010). 13. Y. Shechtman, S. J. Sahl, A. S. Backer, and W. E. Moerner, “O ptimal point spread function design for 3D imaging,” Phys. Rev. Lett. 113, 133902 (2014). 14. F. Aguet, S. Geissbuhler, I. Marki, T. Lasser, and M. Unse r, “Super-resolution orientation estimation and localization of fluorescent dipoles using 3-D steerable filters,” Opt. Express17, 6829–6848 (2009). 15. X. Michalet and A. J. Berglund, “Optimal diffusion coeffic ient estimation in single-particle tracking,” Phys. Rev. E 85, 061916 (2012). 16. M. Born and E. Wolf,Principles of Optics(Cambridge University, Cambridge, UK, 2002), 7th ed. 17. S. F. Gibson and F. Lanni, “Experimental test of an analyti cal model of aberration in an oil-immersion objective lens used in three-dimensional light microscopy,” J. Opt. Soc . Am. A 9, 154–166 (1992). 18. B. Zhang, J. Zerubia, and J. -C. Olivo-Marin, “Gaussian a pproximations of fluorescence microscope point-spread function models,” Appl. Opt. 46, 1819–1829 (2007). 19. P. Torok and F. J. Kao, Optical Imaging and Microscopy (Springer Verlag, NY, USA, 2003). 20. M. Unser, A. Aldroubi, and M. Eden, “B-spline signal proc essing: part I–theory,” IEEE Trans. Signal Process., 41, 821–832 (1993). 21. M. Unser, “Splines: a perfect fit for signal and image proce ssing,” IEEE Signal Process. Mag. 16, 22–38 (1999). 22. C. De Boor,A Practical Guide to Splines (Springer Verlag, NY, USA, 2001), Rev. ed. 23. I. J. Schoenberg, “Spline functions and the problem of gr aduation,” Proc. Natl. Acad. Sci. USA 52, 947–950 (1964). 24. I. J. Schoenberg, “Cardinal interpolation and spline fu nctions,” J. Approx. Theory2, 167–206 (1969). 25. M. Arigovindan, M. Suhling, P. Hunziker, and M. Unser, “V ariational image reconstruction from arbitrarily spaced samples: a fast multiresolution spline solution,” IEE E Trans. Image Process. 14, 450–460 (2005). 26. X. Lai, Z. Lin, E. S. Ward, and R. J. Ober, “Noise suppressi on of point spread functions and its influence on deconvolution of three-dimensional fluorescence microscopy image sets,” J. Microsc. 217, 93–108 (2005). 27. C. Claxton and R. Staunton, “Measurement of the point-spr ead function of a noisy imaging system,” J. Opt. Soc. Am. A 25, 159–170 (2008). 28. M. Baranski, S. Perrin, N. Passilly, L. Froehly, J. Alber o, S. Bargiel, and C. Gorecki, “A simple method for quality evaluation of micro-optical components based on 3D IP SF measurement,” Opt. Express 22, 13202– 13212 (2014). 29. S. R. P. Pavani and R. Piestun, “Three dimensional trackin g of fluorescent microparticles using a photon-limited double-helix response system,” Opt. Express 16, 22048–22057 (2008). 30. S. Ram, P. Prabhat, E. S. Ward, and R. J. Ober, “Improved sing le particle localization accuracy with dual objective multifocal plane microscopy,” Opt. Express 17, 6881–6898 (2009). 31. R. E. Campbell, O. Tour, A. E. Palmer, P. A. Steinbach, G. S. B aird, D. A. Zacharias, and R. Y. Tsien, “A monomeric red fluorescent protein,” Proc. Natl. Acad. Sci. USA 99, 7877–7882 (2002). 32. J. Chao, E. S. Ward, and R. J. Ober, “Fisher information matr ix fo branching processes with application to electron-multiplying charge-coupled devices,” Multidim. S ys. Sig. Proc. 23, 349–379 (2012). 33. S. Liu, E. Kromann, W. Krueger, J. Bewersdorf, and K. Lidke , “Three dimensional single molecule localization using a phase retrieved pupil function,” Opt. Express 21, 29462–29487 (2013). 34. S. Quirin, S. R. P. Pavani, and R. Piestun, “Optimal 3D sing le-molecule localization for superresolution microscopy with aberrations and engineered point spread func tions,” Proc. Natl. Acad. Sci. USA109(3), 675–679 (2012). 35. S. Proppert, S. Wolter, T. Holm, T. Klein, S. van de Linde, a nd M. Sauer, “Cubic B-spline calibration for 3D super-resolution measurements using astigmatic imaging,” Opt . Ex ress22(9), 10304–10316 (2014). 36. S. Ram, E. S. Ward, and R. J. Ober, “How accurately can a sing le molecule be localized in three dimensions using a fluorescence microscope?” Proc. SPIE 5699, 426–435 (2005).


Introduction
Fluorescence microscopy, a light microscopy technique that enables the detection of specifically-labeled objects, is extensively used to study subcellular structures, proteins and dynamics [1][2][3][4].An important question in fluorescence microscopy concerns the best possible accuracy, in terms of standard deviation, with which an object of interest can be localized.This is particularly important for applications such as localization-based superresolution microscopy in which the spatial resolution is closely related to the localization accuracy [2,4].This problem has been addressed by introducing the practical localization accuracy measure (PLAM) [5][6][7].
The PLAM provides a lower bound on the accuracy with which an unknown parameter, e.g. the location of a subcellular object or a single molecule, can be estimated when imaged using a pixelated detector [5,6].The PLAM is calculated using a well-established statistical tool, the Cramér-Rao lower bound, that is, the inverse of the Fisher information [8,9].The latter represents the amount of information the data provides about an unknown parameter [8].Lower bounds on the accuracy of estimates are useful in assessing the quantitative performance of systems as they can predict how well a system can perform an estimation task without specifying a particular estimator [5,8,10].This feature has turned the PLAM into a reliable measure of accuracy that helps to optimize the design of fluorescence microscopy experiments [4,[10][11][12][13][14][15].
However, the calculation of the PLAM to date relied on an analytical expression for the image of the object, which we refer to as the image function [5,10,16].In practice, this can be problematic owing to the fact that often no accurate analytical image function is available [16][17][18].Even if an appropriate analytical model is available for the image function, the lack of knowledge about the precise values of imaging parameters might also impose difficulties in the calculation of the PLAM, as analytical image functions are typically strongly tied to the parameters of the imaging setup [3,10].For instance, as shown in [19], the experimentally achievable value for the numerical aperture of an objective lens might considerably differ from its nominal value, especially for a high numerical aperture objective lens.
Here, we address the above concerns by developing a new approach that directly makes use of an experimental image set to calculate the PLAM for a general object.A continuously differentiable representation of the experimental image set is necessary to obtain certain derivatives that are required for the calculation of the PLAM.To achieve such a continuously differentiable model we fit splines, i.e. smoothly connected piecewise polynomials, to the experimentally collected image set [20][21][22][23][24]. We use splines since they are well-established in image processing [20,21,25] and have a number of useful properties.Importantly, their derivatives can be obtained analytically [21].Our proposed method only requires the experimental image set and provides the best possible accuracy with which a general subcellular object can be localized.Knowledge of imaging parameters such as the numerical aperture of the objective lens and the refractive index of the immersion oil is not required as these parameters are already encoded in the experimental image set.
Single molecule microscopy is a well-known application of fluorescence microscopy which allows the detection of individual molecules [1,2,5].Due to its practical importance, we study the application of our approach to single molecule microscopy in more detail.In this case, the object of interest is a single molecule which is typically modeled as a point source [2,5] and, as such, the acquired image set pertains to an experimental point spread function (PSF).An experimental PSF is a (3D) PSF obtained by (z-stack) imaging a point source, e.g. a bead [26][27][28][29].We verify our approach using simulations in the presence of extraneous noise sources and give practical examples.We also give non-point source examples.
Although in this study we focus on fluorescence microscopy applications, our proposed approach is generally applicable to imaging problems, and can therefore also be utilized in other applications, e.g.astronomy.

Acquisition of the experimental PSF
A bead sample was prepared using 0.1 µm Tetraspeck™ beads (Life Technologies Corporation, Grand Island, NY) as described in [30].The imaging of the beads was carried out using a Zeiss Axiovert 200 inverted microscope with a Zeiss Plan-apochromat 63x, NA 1.45 oil immersion objective lens (Carl Zeiss, Oberkochen, Germany) at a refractive index of 1.515.The measurement of the 3D PSF was performed by z-stack imaging with a Piezo Flexure Objective Scanner (Physik Instrumente, Karlsruhe, Germany) with a step size of 25 nm.The sample was illuminated by a laser with a wavelength of 488 nm (Toptica Photonics, Munich, Germany) and the emission light (at a wavelength of 650 nm) was passed through a quad-band filter set (Semrock, Inc., Rochester, NY).An electron multiplying charge coupled device (EMCCD) camera iXon DU897-BV (Andor Technologies, South Windsor, CT) with conventional readout and a pixel size of 16 µm × 16 µm was used to acquire the data.

Acquisition of the lysosome images
A human prostate carcinoma epithelial cell line, 22Rv1, was obtained from the ATCC (Manassas, VA) and was maintained in RPMI-1640 (Lonza, Basel, Switzerland) supplemented with 10 % FCS (HyClone Laboratories, Inc., Logan, UT).For imaging studies, the culture medium was replaced with phenol-red free RPMI-1640 (Invitrogen, Carlsbad, CA) supplemented with 10 % FCS (HyClone Laboratories, Inc.).The 22Rv1 cells were transfected with expression plasmids encoding LAMP-1 with C-terminally linked red fluorescent protein (mRFP; [31]) using an Amaxa Nucleofector™ (Lonza) instrument and program X-001.36 hours following transfection, the cells were imaged as live cells using a Zeiss Axio Observer A1 inverted microscope with a Zeiss Plan-apochromat 100x, NA 1.4 oil immersion objective lens (Carl Zeiss) at a refractive index of 1.515.The sample was illuminated by a 543 nm laser (Opto Engine LLC, Midvale, UT).A charge coupled device (CCD) camera Orca ER (Hamamatsu, Bridgewater, NJ) with a pixel size of 6.45 µm × 6.45 µm was used to acquire the data.

Computations and software
All of the computations were carried out in a custom-written software package developed in the MATLAB environment (The MathWorks Inc., Natick, MA).This tool is capable of calculating the Fisher information matrix and the PLAM for both 2D and 3D experimental image sets.

Fisher information matrix and problem formulation
In this section, we briefly explain the theory for determining the best possible localization accuracy in single molecule microscopy.For a 3D localization problem, we denote the location of the object of interest in the object space by the parameter vector θ := (x 0 , y 0 , z 0 ) ∈ Θ, where Θ ⊆ R 3 , is an open parameter space.For a 2D localization problem, the location parameter vector is obviously reduced to θ := (x 0 , y 0 ) ∈ Θ ⊆ R 2 .The best possible accuracy with which the location of the object can be estimated, observing its pixelated image, is given by the practical localization accuracy measure (PLAM) [5][6][7].The PLAM is determined using the Cramér-Rao lower bound (CRLB) [5,10].According to the Cramér-Rao inequality [8,9], the covariance matrix of any unbiased estimator θ of a parameter vector θ ∈ Θ is always greater than or equal to the inverse Fisher information matrix (FIM), i.e. cov( θ) ≥ I −1 (θ).The main diagonal elements of the inverse FIM provide lower bounds on the variance of the estimates of the unknown parameters, whereas we are interested in the estimation accuracy in terms of the standard devi-ation.Hence, the PLAM vector is defined as the element-wise square root of the main diagonal entries of the inverse FIM [6,7].
We next express the FIM for the single molecule microscopy problem.Let {C 1 , . . .,C K pix } be a pixelated detector, where C k ⊆ R 2 denotes the area occupied by the k th pixel and K pix is the total number of pixels.The pixels are assumed to be disjoint.It has been shown that the photon counts detected by the pixels of the detector due to the object of interest are the realizations of independent Poisson random variables with expected values [5,9] where r := (x, y) ∈ R 2 , N is the expected number of photons that impact the infinite detector plane (i.e.R 2 ) due to the object, M is the lateral magnification of the objective lens and q z 0 is the image function [5,10].The image function is a bivariate probability density function (pdf) that describes the image of a stationary object on the detector plane at unit lateral magnification when it is located on the optical axis at position z 0 ∈ R [4,10].The image function for a 2D localization problem is simply given by setting z 0 = 0.In case that the object of interest is a point source, the image function is identical to the PSF of the optical system.For example, considering the standard Born and Wolf 3D PSF model [16], the image function is given by [3] where A is a normalization constant, n oil denotes the refractive index of the immersion medium, e.g.oil, and J 0 is the zeroth order Bessel function of the first kind [16].
It has been shown that in the presence of extraneous noise, the expression of the FIM is given by [5,10] where v θ (k) := µ θ (k) + b k with b k , k = 1, . . ., K pix , denoting the photon count due to the background signal at pixel C k .If the parameter vector is independent of the background level b k , we have The term α(k), k = 1, . . ., K pix , is known as the noise coefficient that depends on the extraneous noise sources and the detector type.In the absence of readout noise, α(k) = 1 for all k = 1, . . ., K pix [5].In the presence of readout noise and when using a CCD camera or a complementary metal oxide semiconductor (CMOS) detector the noise coefficient is given by [10] α(k where η k and σ 2 k denote the mean and the variance of the readout noise at pixel C k , respectively.The expression of the noise coefficient in the presence of readout noise and stochastic signal amplification, i.e. when using an EMCCD camera, can be found in [32]. Supposing that the object of interest is a point source and that an analytical expression is available for the PSF (e.g. the Airy profile [16] or a bivariate Gaussian profile [18] assuming a 2D case, and the Born and Wolf model [16], i.e.Eq. ( 2), assuming a 3D case), Eq. ( 3) can be used to calculate the PLAM for a single molecule microscopy experiment.However, as mentioned earlier, the lack of appropriate analytical models for PSFs and the lack of knowledge about the precise values of imaging parameters often cause major problems in the calculation of the PLAM.Additionally, it is often important to calculate the PLAM for a general experimental object as opposed to a point source.To overcome these problems, we propose an alternative approach by directly making use of an experimental image set for the calculation of the PLAM.
Other approaches are reported in the literature to address the model mismatch issue.For instance, in [33] a phase-retrieved pupil function was used to generate a more accurate model for the PSF of the optical system.This more accurate PSF model was then used for the calculation of the PLAM.In [34], a similar approach was used to model engineered PSFs.Such techniques, however, are limited to point-like objects and depend on a variety of imaging parameters, such as the numerical aperture of the objective lens (see e.g.[33]).

Experimental image sets and experimental PSFs
In this section, we develop notation for an experimental image set which will be useful for our later discussions.A 3D experimental image set is a set of pixelated images of an object acquired at different defocus levels [27,28], which are corrupted by extraneous noise sources, such as background and readout noise, during the measurement process [5].In addition, due to the stochastic nature of light, the acquired images are also inherently stochastic [9,10].Let z p ∈ R, p = 1, . . ., K stk , denote the defocus level in the object space, where K stk is the total number of levels.We define an acquired 3D experimental image set as a realization where N c > 0 is the expected photon count, * denotes the convolution operator, b c k,p ≥ 0 is the background level at pixel C k , k = 1, . . ., K pix , at defocus level z p , p = 1, . . ., K stk , and N (0, σ 2,c k ) denotes a zero-mean Gaussian distribution with variance σ 2,c k associated with the readout noise.If the microscope system is spatially-invariant along the z-axis, we have q z 0 ,z p := q z p −z 0 .The above notation can also be used for a 2D localization problem simply by assuming K stk = 1 and z p = z 0 .In this case, the experimental image set contains only a single image of the object.If the object of interest is a point source, the experimental image set pertains to an experimental PSF which can be collected by imaging a bead sample (see Section 2.1).
After acquiring an experimental image set, e.g. an experimental PSF, the next step is to estimate the image function q z 0 .Once we estimate the image function, we can substitute it into Eq.( 1) to obtain an analytical expression for µ θ (k), k = 1, . . ., K pix , which can then be used to analytically calculate the partial derivatives required in the FIM equation (i.e.Eq. ( 3)).This will be the topic of subsequent sections.

Piecewise polynomial fitting
Splines are piecewise polynomials with pieces that are smoothly connected together [21].They have been extensively used in multidimensional data fitting (e.g.surface fitting) and interpolation problems due to their useful properties [20,22,35].One of the important characteristics of a spline is that it can be represented in the form of a linear combination of basis functions known as B-splines [20].B-splines have a number of important properties, namely affine invariance, local support and positivity [24], which make them of interest for our application.We therefore take advantage of splines to estimate the image function.In particular, we next explain how to fit a volume spline to a 3D experimental image set.Fitting a surface spline to a 2D experimental image set is a special case of the 3D fitting by simply setting K stk = 1 and z p = z 0 .Denote by ∆x > 0, and ∆y > 0, the physical pixel size in the image space in the x and y directions, respectively.Let ∆x 0 := ∆x/M, and ∆y 0 := ∆y/M, be the effective pixel size in the object space in the x and y directions, respectively, where M is the lateral magnification of the microscope optics.Let ∆z 0 > 0 be the step size in the z-direction in the object space.A volume spline of degree d ∈ N 0 with element spacing (∆x 0 , ∆y 0 , ∆z 0 ) in the object space is given by [24] s d a (x, y, z) := where {a m,n,p | m = 1, . . ., K row , n = 1, . . ., K col , p = 1, . . ., K stk } are called the B-spline coefficients, K row and K col denote the number of rows and columns of the image, respectively, such that K row × K col = K pix , K stk denotes the total number of defocus levels, and β d denotes the symmetrical B-spline of degree d given by (see Fig. 1) where Given the noisy measurements h k,p at pixels C k , k = 1, . . ., K pix , and at defocus levels z p , p = 1, . . ., K stk , our problem is to find a volume spline s d a (x, y, z) for (x, y, z) ∈ R 3 , such that where b c k,p denotes the background level at pixel C k and at defocus level z p , and is assumed to be known or can be estimated [3].It is important to note that the above problem is an interpolation problem with the exception that, for each defocus level z p , the data points are calculated by integrating the continuous surface spline over the pixels instead of evaluating it at the centers of the pixels/intervals (see e.g.[20,25]).This integral sampling is the appropriate way of modeling the photon detection process in fluorescence microscopy [5,9,10].Our problem is, in fact, to find a set of B-spline coefficients that minimizes the cost function To introduce a concise matrix notation for the above cost function we define and where K := K pix × K stk is the total number of data points.We also define S ∈ R K×K such that where r = (x, y) ∈ R 2 .Using this matrix notation the cost function is given by ε(a where • denotes the Euclidean (ℓ 2 ) norm.
Minimizing the cost function in Eq. ( 8) leads to an exact spline fit for the experimental image set (i.e.zero error for the cost function).However, in practice and as described in Section 3.2, experimental image sets are inherently stochastic and are typically corrupted by extraneous noise.As such, an exact spline fit does not necessarily provide the best continuous approximation.Hence, we regularize the optimization problem using an additional term that intends to suppress the noise [20,25] where D l is the vector of all possible partial derivatives of order l.For instance, for l = 1 we have We next express the regularization cost function, i.e.Eq. ( 9), in terms of the expansion of the B-spline coefficients.This will help to introduce a matrix notation for the regularization term.From Eq. ( 9), for a ∈ R K , we have Substituting Eq. ( 5) into the above equation and a little manipulation yields (for details see Appendix A) where Using this notation, Eq. ( 10) can be expressed in matrix form as follows (see also [25]) where for conciseness the superscript d and subscript l are dropped.
To estimate the B-spline coefficients in the presence of stochasticity and noise, by making use of the matrix notation introduced above, we solve the following optimization problem which is a regularized least-squares problem [20,21].The first term measures the error between the data and the model in the least squares sense whereas the second term imposes a smoothness constraint on the solution.The regularization (smoothing) factor γ ≥ 0 controls the trade-off between fidelity to the data and the smoothness of the estimate.Using vector differentiation [25], it is easy to verify that the minimizer to Eq. ( 11) is given by the solution of the following equation which can be solved efficiently using Gaussian elimination or singular value decomposition.We note that the solution of the above equation can be derived given a specific choice of the smoothing factor γ, the derivative order l and the B-spline degree d.The smoothing factor can be chosen based on a priori information, e.g. the variance of the measurement noise.By setting γ = 0, the optimization problem in Eq. ( 11) reduces to a standard least squares problem [20].
The typical choice for the derivative order in modern statistics literature is l = 2, although other orders can also be easily used [22].Given the order of derivatives, an appropriate degree for the B-splines can be chosen as d = 2l − 1 [20].The rationale for this choice is Schoenberg's work [23] in which it is demonstrated for a 1D problem that the solution that minimizes the error in Eq. ( 11) is a spline of degree d = 2l − 1 with simple knots at the data points and some natural end conditions.For instance, cubic spline interpolation (i.e. using a spline of degree d = 3) is the appropriate choice when using the derivatives of order l = 2.

Calculation of the Fisher information matrix
Once we estimate the B-spline coefficients â through Eq. ( 12), we can substitute them into Eq.( 5) and find the spline fit to the experimental image set h.This volume fit ŝd a after normalization can be used to obtain an estimate of the image function.For conciseness, define ∑ m,n,p := ∑ We define the normalization factor where r = (x, y) ∈ R 2 and we applied the B-spline property R β d (x)dx = 1 for d ∈ N 0 (see Appendix B for details).The estimated image function is given by where (x, y) ∈ R 2 , and ãz 0 m,n,p := âm,n,p /C(z 0 ), m = 1, . . ., K row , n = 1, . . ., K col , p = 1, . . ., K stk , are termed the normalized B-spline coefficients.
We now have an estimate of the image function that can be used to calculate the PLAM.Substituting the estimated image function into Eq.( 1), for k = 1, . . ., K pix , we have It is important to note that, assuming that the pixel size, the magnification and the location of the object are unchanged, the integral in the above expression is constant and therefore it can be precalculated once and used in different experiments.The next step is to calculate the partial derivatives of µ θ (k) with respect to the unknown parameters.An interesting feature of B-splines is that their first derivatives can be obtained analytically through the following expression [21] Using this identity and taking the partial derivatives of both sides of Eq. ( 14) with respect to x 0 for k = 1, . . ., K pix , we have (for details see Appendix C) where r = (x, y) ∈ R 2 and we assumed ãz 0 m,0,p = ãz 0 m,K col +1,p = 0, m = 1, . . ., K row , p = 1, . . ., K stk .

Limit of the accuracy for estimating other parameters
In the previous section, we were primarily concerned with the calculation of the best possible localization accuracy directly from an experimental image set.However, using the general statistical framework described in Section 3.1 we can also calculate the best possible accuracy with which other parameters can be estimated, such as the photon count and the background level.
To this end, we simply define an extended parameter vector θ := (x 0 , y 0 , z 0 , N, b) ∈ Θ ⊆ R 5 , where b := b k , k = 1, . . ., K pix , denotes the constant background level.Since the photon count is independent of the estimated image function, from Eq. ( 14) it is straightforward to verify that Similarly, for the background level we have By substituting the above equations into Eq.( 3), we can obtain the best possible accuracy for estimating the photon count and the background level.

Verification of the approach in the absence of noise
We have developed an approach for the calculation of the best possible accuracy with which a general object can be localized, i.e. the PLAM, directly from 2D and 3D experimental image sets.We defer to Sections 4.4 and 4.5 examples concerning the PLAM for general experimental objects.Here we primarily focus on point-like objects (e.g. a single molecule) and, as such, the experimental image set pertains to an experimental PSF.We refer to the PLAM deduced from an experimental PSF using the aforementioned approach as the experimental PLAM, whereas the PLAM calculated from an analytical PSF is referred to as the analytical PLAM.We further refer to the limit of the localization accuracy for the x, y and z coordinates of the single molecule as x 0 -PLAM, y 0 -PLAM and z 0 -PLAM, respectively.For a 2D PSF, the z 0 -PLAM is not relevant.
It is important to verify the performance of the developed approach in terms of the deviation of the experimental PLAM from the analytical PLAM for a simulated experimental PSF in idealized imaging conditions.We dedicate this section to this verification.
In particular, we assume idealized imaging conditions where the experimental PSF is devoid of stochasticity, due to the photon statistics, and extraneous noise, such as Poisson-distributed background noise and Gaussian-distributed readout noise.Considering these conditions, we simulate 2D and 3D experimental PSFs using analytical PSFs, such as the Airy PSF model and the Born and Wolf PSF model, respectively [16].We then calculate the experimental PLAM using the proposed approach given the simulated experimental PSFs.We investigate the performance of our approach by comparing the results in two different steps.The first step concerns the comparison of the image profiles of the spline fit to the simulated experimental PSF and its associated analytical model.The second step compares the deduced experimental PLAM with its corresponding analytical PLAM.
We first study a 2D case.Figures 2(a focus point source, i.e. the Airy profile [16], which is used as the 2D analytical PSF model, and the simulated 2D experimental PSF using this model, respectively.A bicubic spline is then fit to the simulated experimental PSF (see Fig. 2(c)).A Comparison of the cross sections of the fit and the analytical PSF suggests that in idealized imaging conditions the spline fit provides an appropriate estimate of the image function (see Fig. 2(d), the (absolute difference) error is consistently less than 4 × 10 −14 photons/pixel).We next compare the deduced experimental PLAM with its corresponding analytical PLAM.In idealized imaging conditions, altering the effective pixel size, i.e. the physical pixel size divided by the lateral magnification, can potentially vary the deduced PLAM since it changes the spatial sampling of the PSF.Therefore, we calculate the experimental and analytical PLAMs as functions of the effective pixel size in Fig. 2(e).For effective pixel sizes ranging from 0.05 µm to 0.2 µm, the experimental x 0 -PLAM is very close to the analytical x 0 -PLAM.For this range the absolute difference error is consistently below 0.6 nm (i.e. the relative error is below 7 %).However, as the effective pixel size increases beyond 0.2 µm, the experimental x 0 -PLAM deviates from the analytical x 0 -PLAM.For example, the relative error for the effective pixel size of 0.25 µm is approximately 35 %.This result is not surprising considering the fact that a large effective pixel size implies a coarser spatial sampling of the analytical PSF which, in turn, leads to a less accurate spline fit.We note that effective pixel sizes beyond 0.2 µm, however, are rare in practice.
For the verification of the 3D case, an experimental PSF is simulated using the Born and Wolf analytical 3D PSF model [16] (see Figs. 3(a) and 3(b) for xy-and xz-projections).A cubic volume spline is then fit to the simulated 3D experimental PSF (see Figs. 3(c) and 3(d) for xy and xz projections).The negligible error between the spline fit and the analytical PSF, evaluated at the pixels and at the z-steps, suggests that in idealized imaging conditions the spline fit provides an accurate estimate of the image function (see Fig. 3(e) for the xz projection of the error, where the absolute error is consistently less than 10 −12 photons/pixel/z-step).We next compare the deduced experimental PLAM with its corresponding analytical PLAM as a function of the z 0 position of the single molecule on the optical axis.This is of significant practical importance as it allows us to verify whether the deduced experimental PLAM remains valid as the single molecule moves along the z-axis.Figure 3(f) shows the deduced experimental x 0 -PLAM, its corresponding analytical x 0 -PLAM and the absolute deviation error over the zrange of [0, 1] µm (since the PLAMs are axially symmetric, the results for the z-range of [-1, 0] µm are omitted).The error is consistently smaller than 0.9 nm over the z-range of [-1, 1] µm and the average percentage error over this range is approximately 2 %.Similar results are observed for the deduced experimental z 0 -PLAM (see Fig. 3(g)), where the error is consistently smaller than 4.3 % over the z-range of [-1, 1] µm, with an average of 3.9 %.We would like to note that these results are for a pixel size of 13 µm × 13 µm and that the error can be further decreased by decreasing the pixel size.

Effects of stochasticity and noise in the experimental PSF on the deduced PLAM
In practice, experimental PSFs are inherently stochastic, due to the Poisson distribution of the collected photons, and are typically corrupted by extraneous noise [5,26,28].Hence, in this section we investigate the effects of stochasticity and noise in the experimental PSF on the deduced experimental PLAM.Figures 4(a), 4(a') and 4(b), 4(b') show an out of focus xyprojection of the Born and Wolf 3D PSF model [16], and the corresponding xy-projection of a simulated 3D experimental PSF using this model, respectively.For the simulation of the experimental PSF, we considered both the stochasticity of the signal and the extraneous noise sources.As mentioned earlier, in the presence of noise, an exact spline fit is not appropriate.Therefore, we fit a volume cubic smoothing spline, with a small smoothing factor γ = 0.01, to the simulated experimental PSF (see Figs. 4(c) and 4(c')).Comparing the error between the simulated experimental PSF and the analytical PSF model (Fig. 4(d)) and between the spline fit and the analytical PSF model (Fig. 4(d')), evaluated at the pixels, suggests that the smoothing spline fitting can suppress the stochasticity and noise (the error range is reduced from [-2, 2] photons/pixel in Fig. 4(d) to [-1, 1] photon/pixel in Fig. 4(d')).This can also be obviously observed in Fig. 4(c'), which shows a similar pattern to Fig. 4(a').
Figure 4(e) shows the deduced experimental x 0 -PLAM, its corresponding analytical x 0 -PLAM and the absolute deviation error over the z-range of [0, 0.8] µm.Although we assumed a significant amount of readout noise (σ c = 10 e − /pixel), the error is relatively small over the zrange of [0, 0.8] µm.Specifically, it is consistently smaller than 1.03 nm over this range (i.e. the average error is 3.2 %).Similar results are observed for the deduced experimental z 0 -PLAM in the presence of noise (the average error over the z-range of [0, 0.8] µm is 11.7 %, see Fig. 4(f)).
However, the error between the experimental and analytical PLAMs is a function of the noise level, the B-spline degree and the smoothing factor.We next investigate this important dependence.For this purpose, we define the root mean square percentage error (RMSPE) as follows where PLAM E and PLAM A denote the experimental and analytical PLAMs, respectively, and z i ∈ [0.2, 0.9] µm, for i = 1, . . ., P. This integral error proves to be useful for studying the dependence of the experimental PLAM on the noise level, since it combines the errors in x 0 -, y 0and z 0 -PLAMs in a single measure.This error is plotted in Figs.5(a) and 5(b) as a function of the standard deviation of the readout noise and the background level, respectively, for different B-spline degrees.Not surprisingly, the error is on average monotonically increasing with the noise level regardless of the degree of the B-spline.At small noise levels, high B-spline degrees yield smaller errors than small B-spline degrees whereas this is reversed at high noise levels.d), the B-spline degree is set to 3. Each data point is the average of the errors for x 0 -PLAM, y 0 -PLAM and z 0 -PLAM calculated over a z-range of [0.2, 0.9] µm.For the calculation of the PLAMs, we assumed a background level of b = 10 photons/pixel and a photon count of N = 500 photons.Other parameters are the same as those used in Fig. 4.
For instance, the error levels for the B-splines of degree d = 1 and degree d = 7 are 11.2 % and 9.3 %, respectively, when the standard deviation of the readout noise is σ c = 0 e − /pixel (see Fig. 5(a)).On the other hand, at σ c = 12 e − /pixel, the aforementioned error levels change to 12 % and 13.82 %, respectively.Importantly, for noise levels that are typically observed in practice, i.e. readout noise with standard deviations 0 to 12 e − /pixel and background levels 0 to 200 photons/pixel, a B-spline of degree 3 appears to provide the smallest amount of error.Specifically, for a cubic B-spline the error remains in the range of 8 % to 11 % for the mentioned noise levels (see Figs. 5(a) and 5(b)).
Figures 5(c) and 5(d) show the error as a function of the standard deviation of the readout noise and the background level, respectively, for different smoothing factors.When the smoothing factor is zero, the error significantly increases with increasing noise level (e.g. it increases from 15 % to 35 % as the background level increases from 0 to 200 photons/pixel).As we increase the smoothing factor (e.g. to 0.001), the numerical values of the error decrease and the error becomes less dependent to the changes in the noise level.A relatively small smoothing factor of 0.01, for example, yields a relatively robust curve to the noise (with this smoothing factor, the error increases only from 8.9 % to 14.3 % as the background level increases from 0 to 200 photons/pixel).On the other hand, a large smoothing factor (e.g.γ = 1) leads to an oversmoothed spline fit.The loss of information due to over-smoothing, in turn, results in a large error (see Figs. 5(c) and 5(d)).Consequently, a relatively small smoothing factor (e.g.0.005 to 0.05) can be an appropriate choice for typical noise levels in practice.

Experimental PSF example
In this section, we provide an example to investigate the performance of the proposed approach in practice.In particular, we collect the 3D experimental PSF of a microscopy setup using the procedure described in Section 2.1.We have deliberately used a setup with an aberrated PSF as it is a good example to illustrate the practical performance of the proposed approach.Figures 6(a) and 6(b) show the yz-and xy-projections of the acquired experimental PSF, respectively.To suppress the stochasticity and noise in the collected experimental PSF, based on the analyses reported in the previous section, we fit a volume smoothing spline of appropriate degree and smoothing factor to the experimental PSF.The yz-and xy-projections of the smoothing spline fit are shown in Figs.6(a') and 6(b'), respectively, where we see a substantial suppression of the extraneous noise.We further calculate the experimental x 0 -PLAM and z 0 -PLAM along the z-axis, which are shown in Figs.6(c) and 6(d), respectively (the experimental y 0 -PLAM is analogous to x 0 -PLAM and is not shown).The experimental x 0 -PLAM has smaller numerical values at or close to the plane of focus and increases as the particle moves away.This is an expected result for typical 3D PSFs (e.g. the Born and Wolf PSF) [36].A subtle point in the behavior of the experimental x 0 -PLAM is that it is not symmetric with respect to the plane of focus.For example, the numerical value of the x 0 -PLAM is 30 nm at z 0 = −0.6 µm, whereas it is approximately 26 nm at z 0 = 0.6 µm.This is not surprising since any mismatch between the refractive indices of the sample and immersion medium contributes to an axially asymmetric PSF [17].
Additionally, the experimental z 0 -PLAM is large near or at the focal plane, e.g. it is 144 nm at the focal plane, and decreases as the point source moves away from the focal plane (see Fig. 6(d)).The large numerical value of the experimental z 0 -PLAM at the focal plane is sometimes referred to as the depth discrimination problem and is expected (see e.g.[3,36]).The axial asymmetry of the z 0 -PLAM can also be explained by the axial asymmetry of the PSF caused by the mismatch between the refractive indices of the sample and immersion oil [17].

Spherical shell example
As mentioned earlier, the proposed approach allows the calculation of the PLAM for general experimental objects.This section provides an example for the calculation of the experimental PLAM for such general objects.In particular, we simulated the 3D image set for a spherical shell by convolving the simulated 3D object with the 3D PSF of the optical system which was assumed to be the Born and Wolf model, and by considering the stochasticity due to Poisson statistics and extraneous noise sources (see Figs. 7(a)-7(c)).We then fit a volume smoothing spline to the experimental image set (see Figs. 7(a')-7(c')).We supposed that at z 0 = 0 the equator of the spherical shell is in focus and calculated the experimental PLAM for this 3D object as it moves along the z-axis.
As shown in see Fig. 7(d), the experimental x 0 -and y 0 -PLAM are small when the object is in focus and increase as the object moves away from the plane of focus.This is an expected result and is analogous to the behavior of the x 0 -and y 0 -PLAM for a point source.Interestingly, the experimental z 0 -PLAM is large when the object is at or near the plane of focus which implies that the z 0 -position of the object cannot be estimated accurately (see Fig. 7(e)).This behavior is analogous to the depth discrimination problem for point sources [3,36].For instance, the experimental z 0 -PLAM is 108 nm when the object is 150 nm away from the plane of focus.This is mainly due to the fact that the image profiles of slices of a spherical shell near its equator appear similar (e.g.Fig. 7(a')).By increasing the z 0 , the experimental z 0 -PLAM reduces significantly.For example, at z 0 = 700 nm, the experimental z 0 -PLAM is 54 nm.Further increasing the z 0 , gradually worsens the experimental z 0 -PLAM.The reason for this behavior is that at very large z 0 -positions, the image profiles of the object become spread out over many pixels and, as such, the photon count per pixel will be negligible compared to the noise level.

Experimental non-point source example
Following the previous section, we now provide an experimental example for the calculation of the PLAM for non-point-like objects.Specifically, we imaged lysosomal compartments in live cells as described in Section 2.2 (see Fig. 8(a)).Figures 8(b) and 8(c) show two individual lysosomal compartments.We fitted surface smoothing splines to theses images (see Figs. 8(b') and 8(c')) and calculated the experimental PLAM.Using the proposed algorithm, the experimental PLAM can be calculated for any arbitrary photon count and noise level.Here we assumed a background level of 20 photons/pixel and a photon count of 5000 photons.The results are shown in Fig. 8(d).For instance, for the lysosome shown in Fig. 8(b) the experimental x 0 -and y 0 -PLAMs are 11.52 nm and 12.86 nm, respectively.

Conclusions
In this paper, we proposed an approach to determining the limit of the accuracy with which a general subcellular object, imaged using a fluorescence microscope, can be localized directly from an experimental image set.This technique, unlike traditional methods, does not rely on an analytical expression for the image of the object and therefore avoids potential model and parameter mismatch issues.We studied in detail a special case where the object of interest is a point source and, as such, the experimental image set reduces to an experimental PSF.We verified our approach using simulations and reported practical and non-point source examples.These developments can help to optimize the design of fluorescence microscopy experiments.

Fig. 1 .
Fig.1.A plot of the symmetrical B-spline of degree d for d = 0, 1, 2, 3, with unit element spacing, which will lead to nearest neighborhood, linear, quadratic and cubic interpolation of the experimental image set, respectively.The vertical dotted lines show the pixel boundaries and the B-splines are located at the center of pixel 3.

Fig. 2 .
Fig. 2. Verification of the approach in idealized imaging conditions for a 2D PSF.(a) The image of a point source simulated using the Airy profile for a 100x, NA 1.4 objective lens with an emission wavelength of 690 nm.(b) An experimental PSF simulated using the model image in (a) in the absence of stochasticity and noise.The pixel size and the detector size are 7 µm × 7 µm and 21 × pixels, respectively.(c) The bicubic spline fit of the simulated experimental PSF in (b).(d)The line profiles that pass through the peaks of the analytical PSF and the spline fit, and the error between the line profiles.(e) The deduced experimental x 0 -PLAM and its corresponding analytical PLAM as a function of the effective pixel size in the object space (the results for y 0 -PLAM are similar due to the radial symmetry of the PSF and are omitted).The absolute difference error is also shown.For the calculation of the PLAM, we assumed a background level of b = 10 photons/pixel and a photon count of N = 500 photons.

Fig. 3 .
Fig.3.Verification of the approach in idealized imaging conditions for a 3D PSF.(a), (b) xy-(at z 0 = 0.8 µm) and xz-projections of an experimental PSF simulated using the Born and Wolf 3D PSF model for a 100x, NA 1.4 objective lens and an emission wavelength of 690 nm in the absence of stochasticity and noise.The refractive index of the immersion oil n oil is set to 1.515.The pixel size, the z-step size and the detector size are 13 µm × 13 µm, 50 nm and 13 × 13 pixels, respectively.(c), (d) xy and xz projections of the cubic volume spline fit of the simulated experimental PSF.(e) The error between the analytical PSF and the spline fit evaluated at the pixels and at the z-steps.(f) The deduced experimental x 0 -PLAM and its corresponding analytical PLAM as a function of the z 0 position of the single molecule in the object space (the results for y 0 -PLAM are similar and are omitted).The error is also shown.(g) The deduced experimental z 0 and its corresponding analytical PLAM as a function of z 0 , and the error between them.For the calculation of the PLAM, we assumed a background level of b = 10 photons/pixel and a photon count of N = 500 photons.

Fig. 4 .
Fig. 4. The performance of the approach in the presence of stochasticity and noise.(a), (a') An xy-projection (at z 0 = 0.9 µm) of the Born and Wolf 3D PSF model for a 100x, NA 1.4 objective lens and an emission wavelength of 690 nm.The refractive index of the immersion oil n oil is set to 1.515.(b), (b') The corresponding xy-projection of the experimental PSF simulated using the model in panels (a) and (a') in the presence of stochasticity and noise, where the standard deviation σ c is 10 e − /pixel and the photon count N c is 10000 photons.The pixel size, the z-step size and the detector size are the same as those in Fig. 3. (c), (c') The xy-projection of the cubic volume spline fit of the simulated experimental PSF, where the smoothing factor γ is set to 0.01.(d), (d') The error between xy-projections of the analytical PSF model and the simulated experimental PSF (i.e.data) and between the analytical PSF model and the spline fit, respectively, evaluated at the pixels.(e), (f) The deduced experimental x 0 -PLAM and z 0 -PLAM and their corresponding analytical PLAMs as functions of the z 0 position of the single molecule in the object space.The errors are also shown.For the calculation of the PLAM, we assumed a background level of b = 10 photons/pixel and a photon count of N = 500 photons.

Fig. 5 .
Fig. 5.The effect of the B-spline degree and the smoothing factor on the error between analytical and experimental PLAMs in the presence of noise.(a), (b) The root mean square percentage error (see Eq. (17) for the definition) between the experimental and analytical PLAMs as functions of the standard deviation σ c of the readout noise and background level b c , respectively, for different B-spline degrees d. (c), (d) The same for different smoothing factors γ.In (a) and (b) the smoothing factor is set to 0.01 and in (c) and (d), the B-spline degree is set to 3. Each data point is the average of the errors for x 0 -PLAM, y 0 -PLAM and z 0 -PLAM calculated over a z-range of [0.2, 0.9] µm.For the calculation of the PLAMs, we assumed a background level of b = 10 photons/pixel and a photon count of N = 500 photons.Other parameters are the same as those used in Fig.4.

Fig. 6 .
Fig. 6.A practical example.(a), (b) The yz-projection and the xy-projection (at z 0 = 2.6 µm) of a deliberately aberrated experimentally collected PSF from a practical microscopy setup, respectively, where the ROI size is 33 × 33 pixels (for information regarding other parameters see Section 2.1).(a'), (b') The corresponding yz-and xy-projections of the cubic volume spline fit with a smoothing factor of γ = 0.01, which is evaluated on a finer grid (color scale bars are in photons).The vertical dashed lines show the location of the plane of focus and the size bars are 1.5 µm (panels (a) and (a') are stretched in the zdirection for better visualization while their scale in the y-direction is the same as panels (b) and (b')).The estimated photon count and background level of the bead sample are approximately N c = 4500 photons and b c = 16 photons/pixel, respectively.(c), (d) The experimental x 0 -PLAM and z 0 -PLAM, respectively, along the z-axis (the reported results are the average of the results for multiple beads).For the calculation of the experimental PLAMs we assumed N = 500 photons and b = 10 photons/pixel.

Fig. 7 .
Fig. 7.The experimental PLAM for a spherical shell.(a) The xy-projection at z 0 = 0 µm, (b) the yz-projection (z 0 ∈ [-0.1, 2.2] µm) and (c) the xy-projection at z 0 = 1.8 µm of the simulated 3D image set of a spherical shell with internal and external radii of 1.5 µm and 1.8 µm, respectively, where the ROI size is 49 × 49 pixels.The image set was obtained by convolving the simulated object with the Born and Wolf 3D PSF model.We considered Poisson statistics, background and readout noise, where the standard deviation σ c is 4 e − /pixel, the background level b c is 10 photons/pixel and the photon count N c is 60000 photons.We assumed a 100x, NA 1.4 objective lens with n oil = 1.515.The pixel size, the zstep size and the emission wavelength are 13 µm × 13 µm, 50 nm and 690 nm, respectively.(a'), (b') and (c') The corresponding xy-, yz-and xy-projections of the cubic volume spline fit with a smoothing factor of γ = 0.015, which is evaluated on a finer grid (color scale bars are in photons).The focal plane is located at 0 µm and the size bars are 1 µm (panels (b) and (b') are stretched in the z-direction for better visualization while their scale in the y-direction is the same as panel (a)).(d), (e) The experimental x 0 -PLAM (y 0 -PLAM) and z 0 -PLAM, respectively, along the z-axis.For the calculation of the experimental PLAMs, we assumed a background level of b = 20 photons/pixel and a photon count of N = 5000 photons.

Fig. 8 .
Fig. 8.The experimental PLAM for lysosomes.(a) The image of a 22Rv1 cell transfected with mRFP-LAMP-1 which was acquired as described in Section 2.2.(b) and (c) The images of two individual lysosomal compartments marked by arrows in panel (a).All of the imaging parameters are reported in Section 2.2.(b') and (c') The corresponding cubic surface spline fit of the lysosomes with a smoothing factor of γ = 0.01, which is evaluated on a finer grid (color scale bars are in photons).The size bar in panel (a) is 3.8 µm and other size bars are 0.645 µm.(d) The experimental x 0 -PLAM and y 0 -PLAM for the lysosomes.For the calculation of the experimental PLAMs, we assumed a background level of b = 20 photons/pixel and a photon count of N = 5000 photons.
the above equation, it follows