## Abstract

One of the ongoing challenges in single particle fluorescence microscopy resides in estimating the axial position of particles with sub-resolution precision. Due to the complexity of the diffraction patterns generated by such particles, the standard fitting methods used to estimate a particle’s lateral position are not applicable. A new approach for axial localization is proposed: it consists of a maximum-likelihood estimator based on a theoretical image formation model that incorporates noise. The fundamental theoretical limits on localization are studied, using Cramér-Rao bounds. These indicate that the proposed approach can be used to localize particles with nanometer-scale precision. Using phantom data generated according to the image formation model, it is then shown that the precision of the proposed estimator reaches the fundamental limits. Moreover, the approach is tested on experimental data, and sub-resolution localization at the 10 *nm* scale is demonstrated.

© 2005 Optical Society of America

## 1. Introduction

Luminescent markers such as fluorescent proteins and quantum dots have become an invaluable tool in biology, where they enable studies of molecular dynamics and interactions in living cells and organisms. Such studies are usually performed with a fluorescence microscope configured to acquire a time-series of two or three-dimensional data, the latter generally in the form of a stack of images taken at different focal distances. These are then processed using particle tracking techniques, which aim to establish the particles’ trajectories from their position in individual acquisitions. Determining these positions with a high level of precision is essential for obtaining biologically significant results; the typical sizes of the commonly employed fluorophores are of the order of 10 *nm*, which is significantly smaller than the optical resolution of the system.

This implies that a fluorescent particle can be assimilated to a point source, and thus, that its image corresponds to a section of the microscope’s three-dimensional point spread function (PSF), degraded by various types of noise. In essence, the localization task then amounts to determining the position of the particle by fitting a model of the PSF to such an image.

In the lateral directions, for particles that are in focus, this is a relatively straightforward task for which several methods have been proposed. Axial localization is more challenging, however, since even when the specimen can be optically sectioned with an arbitrarily fine step [1], localization is made difficult by the poor axial optical resolution, the fact that the PSF of a microscope is non-stationary along the optical axis [2], and the presence of noise [3].

These factors are not only limiting in the case of particle localization, but in any 3D imaging application in microscopy. Consequently, various approaches for improving the resolution of optical microscopes, such as I^{5}M [4], 4-Pi Microscopy [5], and STED [6], have been proposed in recent years, showing that Abbe’s resolution limit can be broken. Alternatively, a system specifically destined for particle tracking localization was introduced by Kao *et al*., who proposed the use of cylindrical optics to encode a particle’s axial position [7], and reached axial resolutions down to 12 *nm*. The downside of these methods is that they require customized hardware, which currently still limits their widespread applicability. In this work, we show that via computational means, particles can be localized with a precision that is clearly beyond the limit traditionally imposed by optical resolution. The method is destined for widefield fluorescence microscopy, which makes it widely applicable.

#### 1.1. Review of computational approaches

The model-based methods proposed for lateral localization typically rely on a simplified diffraction model, or some Gaussian approximation of the PSF (see, e.g. [8, 9]). Notably, Thompson et al. proposed an iterative method based on the minimization of the least-squares difference between an image of a particle and a Gaussian model of the PSF [10], and Ober *et al*. studied the theoretical limits of lateral localization [11]. By computing the Cramér-Rao bound (CRB) for the lateral position, they confirmed that, although the images of single particles are limited by the microscope’s resolution, it is possible to estimate the lateral position with sub-resolution accuracy. In some cases, nanometer-scale localization can be achieved.

To date, only few studies have dealt explicitly with the issue of sub-resolution localization in the axial direction. This can be partly attributed to the scarcity of simple but accurate PSF models for optical microscopes. Several attempts to circumvent the use of a PSF model have been made. Van Oijen *et al*. [12] proposed a method involving a high-resolution z-stack acquisition of the particle (i.e., a series of images taken at different focal distances with regular intervals). It is based on identifying the slice for which the radial size of a Gaussian fit to the diffraction-limited spot is minimal within the z-stack. There are several limitations to this approach, however. It can only work properly if the movement of the particle during the acquisition process is sufficiently slow, and localization is limited by the size and resolution of the z-stack. Additionally, the section of the PSF whose radial size is minimal does not necessarily correspond to the situation where the particle is in the focal plane (we will emphasize this in the following section). Also, for the estimation algorithm to be as fast as possible, it is desirable to localize particles from only one or few acquisitions, without needing to process entire high-resolution z-stacks.

Speidel *et al*. demonstrated the feasibility of sub-resolution localization by experimentally showing that the axial position of a particle can be determined from a single defocused image of the particle [13]. When the particle is sufficiently out of focus, it gives rise to diffraction rings. These authors empirically established a linear relation between the radius of the outermost diffraction ring and the axial position of the particle, which allows them to estimate its position with nanometer precision. This is especially attractive since the estimation becomes possible from a single acquisition. The downside of the approach is again related to the non-stationarity of the PSF, meaning that the linear relationship may vary as a function of the particle’s depth within the specimen. It is also constrained to the localization of particles that are sufficiently out-of-focus such that rings are present in their diffraction patterns.

In principle, it is possible to obtain an analytical solution to axial localization by using a theoretical PSF model to estimate a particle’s position from one or few out-of-focus acquisitions (the diffraction pattern increases in complexity as a function of defocus, thus containing more “information”, but also less signal). In a preliminary report, we have investigated the viability of such an approach by establishing the fundamental theoretical limits with respect to the precision that can be expected in the estimation [14].

#### 1.2. Organization of the paper

In this paper we present a refined particle localization method, built upon a non-stationary theoretical PSF model. We first introduce an image formation model, which also includes the effect of noise. Next we establish the Cramér-Rao bound on axial localization, which gives us the fundamental precision that can be achieved with such an image formation model, independently of the estimator used. We hereby extend the methodology presented by Ober *et al*. [11] to three dimensions. Subsequently, we derive a maximum-likelihood estimator for the axial position, and show that, under ideal circumstances, it reaches the precision predicted by the theoretical bound. In the final part of the paper, we incorporate lateral localization into the maximum-likelihood estimator, and show the validity of our approach by demonstrating the axial localization with sub-resolution precision of fluorescent beads. We also discuss the possibility of optimizing acquisition parameters based on the CRB.

## 2. Materials and methods

#### 2.1. Notations and conventions

To formulate the mathematical expressions throughout this paper, we use an absolute coordinate system placed in the object space of the microscope. We make the hypothesis that a standard microscope setup is used, meaning that the sample consists of a specimen mounted between a microscope slide and a coverslip. We define the origin of our system at the interface between the coverslip and specimen layer (see Fig. 1). The optical axis (z-axis) points from the objective towards the sample, such that distances into the specimen are positive. We denote the position of a particle by (*x*_{p}*,y*_{p}*,z*_{p}
), and a point of observation (corresponding to a point on the focal plane) by (*x,y,z*). When multiple acquisitions at different focal positions are considered, (*x,y,z*_{n}
) corresponds to a point on the nth acquisition. In the first part, where we concentrate on axial localization, we assume that the particle is located in (0,0,*z*_{p}
) for the sole purpose of making expressions as simple as possible. For the sake of consistency, we also express the pixel coordinates of acquisitions in object space (imagine acquisitions being demagnified and projected onto the focal plane). This results in a direct link between the PSF and the image generated on the CCD. Finally, all figures showing xz-sections of PSFs are logarithmically adjusted for intensity in order to emphasize details in the diffraction pattern.

#### 2.2. Simulation parameters and experimental setup

The implementation and simulation of the algorithms were performed using the Matlab programming environment (The Mathworks, Natick, MA). Experimental measurements were carried out on a Zeiss Axioplan 2 microscope system. Both theoretical and experimental results were computed for a 63× magnification, 1.4 NA Zeiss Plan-Apochromat oil-immersion objective. For experimental validation, we prepared samples using fluorescent nanobeads by drying dilutions of TetraSpeck fluorescent microspheres (Molecular Probes, Eugene, OR) onto a slide, and subsequently embedding them under a coverslip using a solid mounting medium of refractive index 1.46. The excitation and emission peaks of these beads are 365 *nm* and 430 *nm*, respectively. In conjunction with this, we used a DAPI beamsplitter corresponding to an excitation wavelength of 365 *nm* and an emission wavelength of 450 *nm*. The physical pixel width of the AxioCam CCD mounted on the microscope is 6.45 *μm*.

In order to verify the estimated position of the particles in our sample during experimental validation, we used a Leica TCS SP2 AOBS confocal microscope, configured to record 50% of the light reflected by the specimen. In this way, the beads along with the coverslip-specimen and specimen-slide interfaces are visible in acquisitions. The microscope has the capability of scanning a single line through the sample, which results in an xz-acquisition. The latter confirmed that all beads were adjacent to the microscope slide. Distance measures on such acquisitions are accurate within 100 *nm*, which is sufficient to indicate the validity of our experimental results.

#### 2.3. Image formation model

We now briefly describe the theoretical PSF model, and put forward an image formation model that incorporates noise.

### 2.3.1. PSF model

The dominant source of aberrations in modern optical microscopes originates from a mismatch between the refractive index of the specimen and those of the immersion and coverslip layers. Objectives are designed for use with a specific immersion medium and coverslip, but cannot compensate for the wide variety of specimen types occurring in practice. In fact, they only produce aberration-free images for sources that are positioned at the coverslip-specimen layer interface. For sources at an arbitrary depth within the specimen, the optical path of light rays differs from the path for a source located at the aforementioned interface. This optical path difference (OPD) then generates spherical aberrations in the images produced by the system. Most importantly, the amount of aberration depends on the distance between the source and the coverslip, implying that the PSF is non-stationary along the optical axis.

In practice, most biological samples feature refractive indices closer to that of water than that of immersion oil (which is required for high NA objectives). Even for objects that are located only a few micrometers below the coverslip, the aberrations induced by the mismatch of refractive indices become non-negligible. Much effort has gone into establishing suitable theoretical models that account for these aberrations. The most accurate ones use vectorial computations based on either the Huygens-Fresnel principle [2] or the Debye approximation [15]. It was recently shown that these two approaches are equivalent in the case of an infinitely high Fresnel number [16], which is a reasonable assumption in biological microscopy. However, when evaluated at the resolutions provided by the sensor grid of currently available CCDs, these models do not yield significant improvements over scalar formulations of the OPD, especially when considering the computational advantages of the latter. This aspect will be further justified in the discussion section of the paper.

Gibson and Lanni [17] proposed a scalar PSF model that is reasonably simple and has the advantage of depending only on the standard parameters of the objective and the optical properties of the specimen, both of which can be determined with adequate accuracy. According to this model (formulated in object space), the response to a point source located in (*x*_{p}*,y*_{p}*,z*_{p}
) is given by

where *W*(*ρ,z*|*z*_{p}
) = *k* · OPD, with *k* being the wave number of the emitted light, NA is the numerical aperture of the objective, *ρ* is the radius of the microscope’s limiting aperture in the microscope’s back focal plane, and *A* is a constant complex amplitude. Due to the hypothesis of spatial invariance in planes orthogonal to the optical axis, the PSF is radially symmetric and can be expressed as a function of the coordinate $r=\sqrt{{\left(x-{x}_{p}\right)}^{2}+{\left(y-{y}_{p}\right)}^{2}}$. We can thus rewrite Eq. (1) as

The detailed expression for the OPD is given in the appendix (Eq. (16)). Note that when imaging a source located at the interface between the coverslip and specimen layers, the PSF corresponds to the standard defocus model [18], where *W*(*ρ,z*|*z*_{p}
) is proportional to -*z*. When, in addition to this, the system is in focus, *W*(*ρ*,0|0) = 0 and *PSF*(*r*,0|0) becomes the familiar Airy function.

### 2.3.2. Noise model

In fluorescence microscopy, noise from a variety of sources contributes to the recorded images of a specimen, depending on the nature of the specimen and the type of image acquisition setup used. The three main sources of noise occurring in CCD devices are photon noise (also called shot noise), dark noise, and read-out noise. For high-performance cameras the latter two can be considered negligible. Photon noise results from statistical variation in the arrival rate of photons incident on the CCD. As a result of the nature of this variation, the recorded signal at a given pixel on the CCD follows a Poisson distribution. Note that the effect of photon noise is particularly important when the energy of the photon-emitting source is low, implying a lower photon flux. We thus define an image formation model where the photon count at a given pixel on the CCD follows a Poisson distribution whose mean is proportional to the intensity predicted by the PSF model. We characterize the ratio between the expected photon count and the predicted intensity by introducing the conversion factor *c*, defined as the amount of photons corresponding to a unitary increase in measured intensity. This factor, along with the constant amplitude |*A*|^{2}, depends on various properties of the experimental setup used, such as the energy of the fluorescent particle, the sensitivity of the CCD sensor (a fixed property of the camera), and the exposure time. Let *q̄* denote the expected number of photons corresponding to the measured intensity due to a point source located at (0,0,*z*_{p}
). Clearly

where the PSF is given by Eq. (1). The probability of observing q photons emitted by a particle located in (0,0,*z*_{p}
) at a point (*x,y*) in the focal plane positioned at *z*_{n}
is then given by

which constitutes the basis for our image acquisition model. Thus, the probability of observing a given spatial distribution *q*(*x,y,z*_{n}
) of photons due to a particle (located in (0,0,*z*_{p}
)) is then given by the joint probability

where *𝓢* is the set of points in object space corresponding to pixels in the CCD array, and *N* corresponds to the number of acquisitions of the particle. In order to simplify the notation, we shall refer to the photon counts *q* and *q̄* without explicitly writing their arguments (*x,y,z*_{n}
|*z*_{p}
).

### 2.3.3. Model-based particle localization

Localization consists of estimating the particle’s position (*x*_{p}*,y*_{p}*,z*_{p}
) from the aforementioned distribution of photons. The estimation is done by fitting a theoretical model to the acquisition(s) of a particle. As opposed to conventional approaches where a generic model such as a Gaussian is used, we perform the localization by fitting our image formation model to the acquisitions.

#### 2.4. Theoretical Bounds

Having formulated the image formation model for a single fluorescent particle, we now proceed with an investigation of the feasibility of axial localization. The aim is to establish the maximal theoretical precision that axial localization can achieve. To determine this maximal precision, we compute the Cramér-Rao bound, which is the theoretical lower bound on the variance of any unbiased estimator. Based on the image formation model, it yields a lower bound on the precision that can be reached in estimating the particle’s axial position *z*_{p}
. Mathematically, the bound states that

where *ẑ*_{p}
is an unbiased estimator of the particle’s position *z*_{p}
. By substituting Eq. (4) into this result and simplifying, we obtain

where the expression for $\frac{\partial \stackrel{\u0304}{q}}{\partial {z}_{p}}$ is given in Eq. (17). The practical relevance of this fundamental result becomes more readily apparent when applying the bound to particular cases, and studying its relationship with the PSF. A simple example is given in Fig. 2, for a source located at the interface between the coverslip and specimen layers. Note the singular behavior of the CRB around the origin, which is related to the depth of field of the microscope. The PSF varies little within the center of that region, and localization becomes less precise. Mathematically speaking, the singularity at the origin is due to the derivative of the PSF, which in this particular case is zero at the in-focus position (since *z*_{p}
= 0).

As indicated by Eq. (7), the shape of the CRB is solely determined by the PSF model, whereas the scale depends on the amount of noise present in the acquisitions. In fact, the amplitude of the CRB is proportional to (*c*|*A*|^{2})^{-1}. As mentioned above, besides exposure time, the energy of the particle is the determining factor for the signal-to-noise ratio. A low-energy particle emits fewer photons, which results in a higher variability in the photon counting process on the CCD. The recorded image will thus be noisier than for a higher energy particle; consequently the CRB will be proportionally higher (meaning that the precision of any estimator will be lower). In our model, the energy of the particle is implicitly related to the amplitude *A*. Another parameter that influences the CRB—but to a lesser extent—is the size of the support within which the particle is observed. At high defocus distances, the support needs to be sufficiently large in order to include the outermost rings of the diffraction pattern.

A more complete illustration of the CRB’s behavior in relation to the PSF and the particle’s depth is given in Fig. 3. For very small changes of *z*_{p}
(not shown, i.e., in our ongoing example, up to ~ 100 *nm*), the PSF can be assumed to be locally stationary. However, as the change of *z*_{p}
increases (up to ~ 1 *μm*), although the shape of the PSF remains essentially the same, a non-negligible axial shift of the PSF with respect to the particle’s position occurs. This phenomenon is accentuated for larger changes of *z*_{p}
, where the “focal shift” increases as a function of the particle’s depth. Incidentally, while the CRB also reflects this shift, it depends much more on the complexity of the diffraction pattern. For sources deeper within the specimen, the diffraction patterns become more complex and the CRB gets lower accordingly. Thus, the bound is much higher for sections of the PSF that resemble a blurred spot, which is not surprising. For a given configuration (i.e., set of acquisition parameters), the value of the CRB changes as a function of the amount of defocus alone. It is minimal only for a specific interval situated approximately between the in-focus position and the positions where the first diffraction rings appear. From this, it is readily apparent that taking out-of-focus acquisitions will lead to a better precision in the estimation.

Having established the fundamental limits on sub-resolution particle localization we now proceed with the development of an estimator whose precision reaches this lower bound.

#### 2.5. A maximum likelihood estimator for axial localization

An optimal maximum-likelihood (ML) estimator for the axial position of a particle is obtained by maximizing the likelihood of our image formation model—in other words of Eq. (5)—with respect to the particle’s position *z*_{p}
:

The ML estimator for the axial position is then obtained by solving for *z*_{p}
in the above expression. Since it is not possible to obtain a closed form solution, we deploy a Newton optimization scheme by linearizing the maximum-likelihood around an estimate of the position. Using the Taylor expansion of the model, we obtain the following first order approximation of Eq. (8):

$$\phantom{\rule{10.2em}{0ex}}+\sum _{n=1}^{N}\sum _{x,y\in \U0001d4e2}\left(\frac{{\partial}^{2}\stackrel{\u0304}{q}}{\partial {z}_{p}^{2}}\left(\frac{q}{\stackrel{\u0304}{q}}-1\right)-{\left(\frac{\partial \stackrel{\u0304}{q}}{\partial {z}_{p}}\right)}^{2}\frac{q}{{\stackrel{\u0304}{q}}^{2}}\right)\left({z}_{p}-{\hat{z}}_{p}\right)\equiv 0,$$

where *ẑ*_{p}
is an initial estimate of the axial position. It is then obvious that the linearization can be performed around the new estimate *z*_{p}
, which implicitly leads to the following iterative expression:

where *m* denotes the *m*th iteration. An adequate initialization for the algorithm is crucial, since the linearization of the likelihood holds only locally. If the initial estimate is too remote from the correct position, convergence of the algorithm is not guaranteed. An efficient way of obtaining an adequate initial estimate ${\widehat{z}}_{p}^{\left(0\right)}$ is to evaluate the normalized cross-correlation between the acquisitions and a number of sections pre-computed from the 3D PSFs corresponding to a range of possible particle positions:

where *μ*_{q}
and *μ*_{q̄}
are the mean values of pixels in the acquisitions and model, respectively. An appropriate stopping criterion for the algorithm can be defined based on the absolute value of the update step. If the latter is smaller than the CRB by an order of magnitude, further refining the estimation is statistically meaningless and the algorithm can thus be stopped.

#### 2.6. Localization in three dimensions

In practice, localizing a particle along the optical axis is not possible without determining its position in the acquisition plane. To this end, an ML estimator for the xy-position can be obtained by making the same developments as for the axial estimator. Since the aim of this work is to demonstrate a new approach for axial localization, we do not state the resulting expressions here. Note, however, that the experimental results presented below were obtained by using an ML estimator for all three dimensions.

## 3. Results

Prior to testing our estimation algorithm on experimental data, we verified its performance in simulation. We generated phantom data by applying Poisson noise to computed sections of the PSF corresponding to a particle at an arbitrary depth within the specimen. The estimation algorithm was then run with these simulated acquisitions, generated for a particle situated at *z*_{p}
= 5 *μm*, using an initial estimate that differed by 0.1 *μm* from the actual value. The process was repeated for various focal distances, using different realizations of the Poisson noise. We then compared the standard deviation of these estimations with the CRB. Fig. 4 shows this result for one particular set of parameters; from other simulations, we have strong evidence that our algorithm achieves the theoretical precision for any configuration.

In their analysis of lateral localization, Ober *et al*. [11] discussed the theoretical limits on estimation precision and used a maximum-likelihood estimator based on their approximative two-dimensional image formation model to show that these limits can be reached. Here, we have presented an analytical expression for a maximum-likelihood estimator based on a complete, three-dimensional formulation of the image formation process, and shown that it reaches the theoretical limits in axial localization. Although they have not been specifically shown here, theoretical bounds on lateral localization can be established for our model as well, and the estimator can be shown to reach these bounds.

#### 3.1. Calibration and experimental setup

Before presenting the results obtained with our localization method, it is necessary to mention how the constant complex amplitude *A* (cf. Eq. (1)) and the photon quantization factor *c* can be estimated in practice. The former can be easily approximated by fitting the model to the data after an initial estimate for the particle’s position has been determined using normalized cross-correlation (which is used precisely because it is independent of *A*). Using a least-squares fit, a sufficiently precise value of *A* is obtained. While an approximation of *c* is not required by the estimation algorithm, it is needed when computing the CRB for experimental data. Assuming that the measures follow a Poisson distribution, the mean is equal to the variance in every point. Since we only have a single realization per point, the mean can be computed using the PSF model (with the estimate of *A*), and the standard deviation approximated with the difference between the model and measure. We obtain an estimate of *c* by computing the sum of pointwise ratios of standard deviation over mean.

Because our method relies on a non-stationary PSF, it requires the knowledge of the focal positions *z*_{n}
in order to estimate *z*_{p}
. In practice, there are two possibilities to obtain these values. The first is to prepare samples in a way such that the focus can be calibrated to the coverslip-specimen interface. This is possible, for example, by drying reference beads onto the coverslip, in order to mark the interface. Since focusing is done using a piezo actuator, the position of the acquisition plane with respect to the interface is then known. The other approach is to include reference beads that are visible over a large axial range in the specimen. By acquiring z-stacks with known relative displacements, and performing a global fit of the PSF model to these stacks, we can determine the position of the acquisition planes together with the locations of the calibration beads. The precision of this calibration increases with the number of views used in the fit (cf. discussion on CRB). In our case, we used the latter approach with 30 views. The initialization step also gave us very precise estimates of the position of the reference beads, which could then be used as gold standards for our experiments.

Experimental data were obtained by acquiring z-stacks with an axial step of 100 *nm* of the fluorescent nano-bead preparations described in the materials and methods section. An xz-section of such a stack is shown in Fig. 5(a). The corresponding section of the theoretical PSF (Fig. 5(b)) shows that the model fits the data well, even for relatively difficult conditions (high NA, depth of the particle). As apparent in the xz-section, but more evidently so in the slices of the z-stack shown in Fig. 6, a non-negligible amount of background noise is present in the acquisitions. This needs to be taken into account for the estimation algorithm and thus requires an extension of our image formation model.

#### 3.2. Extension of the statistical noise model

The mean and variance of the background noise can be estimated to reasonable accuracy from sections of acquisitions that are exempt of fluorescent sources. For the sample discussed at the end of this section, the estimated values for the mean and variance are, respectively, *μ*_{b}
= 514.74 and ${\mathrm{\sigma}}_{b}^{2}$ = 177.82. From these measures it is obvious that the background noise does not follow a Poisson distribution, which suggests that the background in our experiment is due to readout noise (especially since a cooled camera was used). In principle, it is possible to extend our statistical model (Eq. (4)) such as to incorporate background noise, which is typically described as white Gaussian in the literature. To facilitate this extension, we investigate the possibility of approximating background noise with Poisson statistics.

A fundamental property of the Poisson distribution is that it rapidly converges towards a Gaussian with equal mean and variance, given that the latter is large enough, which is usually considered the case when *μ* > 10. Since the variance of the background noise is significantly higher than this value, we make the approximation by splitting the Gaussian distribution into a Poisson distribution and a fixed offset (equal to *μ*_{b}
- ${\mathrm{\sigma}}_{b}^{2}$), which leaves us with the convolution between two Poisson distributions. The convolution of two Poisson distributions yields another Poisson distribution, whose mean is equal to the sum of means from the original distributions. We thus obtain the following extension of our image formation model:

Consequently, the expression for the CRB becomes

with the iterative estimator given by:

To illustrate the appropriateness of this model, we compare in Fig. 6 a few slices of a measured z-stack with the simulated acquisitions obtained using the extended model. When rings are present in the diffraction model, there is an intensity peak at the center of the pattern. If the source is aligned with the xy-grid, this peak is recorded by a single pixel on the CCD. If, however, the source is slightly shifted, the peak’s intensity is distributed across four pixels, like it appears in these examples. Localization in three dimensions was used to determine the particle’s position in the measured z-stack. This estimated position was then used to generate the phantom data.

#### 3.3. Validation with real data

In our acquisitions of the nano-bead sample described in the materials and methods section, several dozen beads were visible. Among these, we chose five beads whose diffraction patterns were well separated, such that estimation errors due to overlap from the patterns of different beads were not possible.

For our setup, the CRB shown in Fig. 3 indicates that localization is much more precise when positioning the focus below the bead, such that diffraction rings appear in the acquisition. The xz-section of the bead confirms this; it is indeed much harder to differentiate two acquisitions that depict a blurred spot of light than two acquisitions that present clearly disparate diffraction patterns. In order to illustrate the performance of our estimator, we thus apply it to acquisitions that feature diffraction rings. Initial values for the particle’s position were obtained using normalized cross-correlation with a series of slices of the PSF model computed with the same axial spacing (100 *nm*) as the experimental acquisitions.

To demonstrate the localization for acquisitions taken at various defocus distances, the estimation was performed using pairs of acquisitions spaced by 200 *nm*, for all such pairs within the acquired z-stacks. Independently of the amount of defocus, the algorithm converges rapidly, requiring 5 iterations on average. Fig. 7 demonstrates the result of the estimation for three beads over a range of 2.5 *μm*. In the best case, localization precision (i.e., standard deviation of the estimation with respect to the reference) of 12.8 *nm* is achieved. The worst result obtained with the selected beads was a precision of 23.8 *nm*. The averages of the estimated positions for the three beads shown in Fig. 7 are 22.046 *μm*, 22.069 *μm*, and 22.085 *μm*, respectively. These values are also in perfect agreement with the reference positions of the beads (22.050 *μm*, 22.073 *μm*, and 22.081 *μm*, respectively), which are obtained using a global fit where all measurements are included. To further confirm our results, we compared our estimates with those obtained using the Leica TCS SP2 AOBS confocal microscope. This acquisition showed that the beads were located approximately between 22.0 *μm* and 22.1 *μm* within the specimen, which is strong evidence for the soundness of our estimations. In Fig. 8 we show the CRB for the shot noise-only image formation model, the CRB for the extended image formation model and the average value of the precision estimation achieved with the beads.

#### 3.4. Optimal acquisition settings

Beyond its theoretical applications, the CRB can also be used to determine optimal acquisition settings that may serve as guidelines to experimentalists. As the evaluation of the CRB showed, it is advantageous to take acquisitions slightly defocused with respect to the particle’s actual position. In practice, however, particles can be situated anywhere within the specimen, and it is therefore not possible to adequately position the focus with respect to an individual particle. Still, the study of particles is usually confined to a predetermined section of the specimen. In such cases, under the hypothesis that the particle’s axial position follows an uniform distribution within the section, optimal focal positions leading to the lowest average CRB can be determined. This optimization is non-trivial, but can be performed by solving the following cost function:

where *a* and *b* are the bounds of the region of interest. In Fig. 9, we show the results of this optimization for a variety of settings. It is immediately clear that the optimal settings are non-trivial. The estimation precision is significantly higher when acquisitions are taken with an optimal focus, especially for particles that are deeper within the specimen. At the same time, these results also show the effect on the CRB of increasing the number of acquisitions. Notice how the CRB decreases as the number of acquisitions is augmented. This is expected, since increasing the amount of “information” on the particle should implicitly lead to a better estimation precision. This property is especially useful in highly noisy acquisition conditions.

## 4. Discussion

By investigating the fundamental theoretical limits of axial localization from defocused acquisitions of sub-resolution fluorescent particles, we have shown that nanometer precision can be achieved. The maximum-likelihood estimator proposed in this work reaches the theoretical limit provided that the image formation model is accurate, which we have experimentally shown to be the case. The use of a non-stationary PSF model makes the localization applicable to any configuration of microscope objectives and specimen preparation; it is especially powerful for localizing particles at any depth within the specimen. Usually, the non-stationary nature of the PSF along the optical axis requires approximative models that suppose stationarity to hold for small layers of the specimen (see, e.g., [19]). Here, we developed an approach based directly on the analytical expression of the PSF, thus guaranteeing convergence within the precision of the theoretical limits.

In our experimental tests we have shown that an axial localization precision below 15 *nm* can be reached. These results confirm the practical applicability of the proposed approach, and demonstrate sub-resolution localization. They also confirm the findings of Speidel *et al*. [13], who were the first to show that nanometer-precision axial localization from defocused acquisitions is possible in widefield fluorescence microscopy. While most localization and tracking approaches that claim such a precision along the optical axis are limited to one or two particles (see, e.g., [20]), the method proposed here can be applied to any number of particles detected within a z-stack. For such a multiple-particle scenario, our model could be extended to account for overlap between the diffraction patterns of particles.

An efficient method to model the combined effect of various sources of noise was introduced, rendering the estimation possible for a wide range of configurations. In particular, incorporating additional sources does not increase the complexity of the model.

An important observation is that the localization algorithm performs significantly better for acquisitions that are taken by placing the focus on the side of the particle where the diffraction pattern is more detailed (in cases where *n*_{s}
< *n*_{i}
, such as in our example, this corresponds to *z*_{n}
> *z*_{p}
). The lesser performance of the estimation on the other side is consistent with the higher value of the CRB (see Fig. 3); we also suspect that it may be partly due to slight discordancies between the PSF model and the experimental observations (see Fig. 5 in the range of -2 to 0 *μm*).

#### 4.1. Influence of the PSF model

We briefly justify our choice of a scalar PSF model for the proposed localization method. Our experiments with the vectorial model proposed by Török *et al*. [15] and Hell *et al*. [2, 16] indicate that the differences with respect to the results obtained using the scalar model are not significant in the context of our work. Studies of the CRB for the vectorial formulation show that in some cases, it is slightly lower than its scalar equivalent (see Fig. 10). However, this is only apparent for strongly out-of-focus acquisitions where the signal intensity is weak and generally undetectable, mainly due to background noise. Also, this effect is most noticeable for less aberrated cases; as one penetrates deeper into the specimen, the CRBs for the two models become virtually equivalent.

Moreover, the scalar model has a clear computational advantage. A vectorial model requires three integrals instead of one for the scalar case, not to mention the fact that the integrands are much more involved. Since the localization algorithm also requires the second derivative of the PSF, the difference in computational cost is considerable.

We note that our methodology is generic enough to accommodate for other theoretical models as further progress is made in this field (see, e.g., [21, 22]).

#### 4.2. Shortcomings and possible extensions of the method

In practice, when a thick section of specimen is considered, a z-stack with sufficient axial resolution (i.e., low spacing between acquisitions) is required to guarantee that all particles present in the specimen are recorded. As a consequence, each particle is visible in multiple slices, which can then be used in the localization. The analytical expression for the CRB can be used to derive the optimal acquisition positions with respect to a particular experiment, in order to maximize the performance of the localization.

A parameter not explicitly taken into account is the temporal resolution of the acquisitions; its determining factor is the movement of the particle during the acquisition of the z-stack. For fast-moving particles, it is still possible to perform the localization, however, by limiting the number of acquisitions. The volume (i.e., depth) of observation is then reduced, and as a consequence localization becomes less precise. Another element that can hinder the efficiency of localization is the diffusion of light occurring within the specimen. Although our approach permits the localization of particles at any depth, it is in this respect limited by a factor that affects any localization method.

The first-order approximation made in the development of the ML estimator holds only locally, meaning that the estimator is very sensitive to the initial estimate. Precision in the latter can be increased, if necessary, by computing the normalized cross-correlation with a finer sampling step (see Eq. (11)). Another possible improvement in this direction might be obtained by using a higher order (e.g. quadratic) approximation of the likelihood function.

These limitations aside, the methodology presented in this work is promising, showing that with a standard widefield fluorescence microscope, particles can be localized with nanometer-scale precision. Our experimental results confirm that the localization precision is comparable to that of specialized hardware such as the setup proposed by Kao *et al*. [7].

## Appendix

In this appendix we recall the essential aspects of the PSF model [17], and provide the full expressions of the derivatives used in the computations of the CRB and the maximum-likelihood estimation.

The optical path difference, noted OPD (= ABCD - PQRS) is approximated as follows:

$$\phantom{\rule{10.2em}{0ex}}+{z}_{p}\left\{\sqrt{{n}_{s}^{2}-{\mathrm{NA}}^{2}{\rho}^{2}}-\frac{{n}_{i}}{{n}_{s}}\sqrt{{n}_{i}^{2}-{\mathrm{NA}}^{2}{\rho}^{2}}\right\}$$

$$\phantom{\rule{10.2em}{0ex}}+{t}_{g}\left\{\sqrt{{n}_{g}^{2}-{\mathrm{NA}}^{2}{\rho}^{2}}-\frac{{n}_{i}}{{n}_{g}}\sqrt{{n}_{i}^{2}-{\mathrm{NA}}^{2}{\rho}^{2}}\right\}$$

$$\phantom{\rule{10.2em}{0ex}}{-t}_{{g}^{*}}\left\{\sqrt{{{n}_{{g}^{*}}}^{2}{-\mathrm{NA}}^{2}{\rho}^{2}}-\frac{{n}_{i}}{{n}_{{g}^{*}}}\sqrt{{n}_{i}^{2}-{\mathrm{NA}}^{2}{\rho}^{2}}\right\}$$

$$\phantom{\rule{10.2em}{0ex}}{-t}_{{i}^{*}}\left\{\sqrt{{{n}_{{i}^{*}}}^{2}-{\mathrm{NA}}^{2}{\rho}^{2}}-\frac{{n}_{i}}{{n}_{{i}^{*}}}\sqrt{{n}_{i}^{2}-{\mathrm{NA}}^{2}{\rho}^{2}}\right\}.$$

As shown in Fig. 11, a ray emitted at point *P* in the design system *PQRS* first passes through a coverslip of thickness ${t}_{{g}^{\mathit{*}}}$
and refractive index *n _{g*}*, then through an immersion medium of thickness ${t}_{{i}^{\mathit{*}}}$
and refractive index ${n}_{{i}^{\mathit{*}}}$
, before entering the front lens of the objective. Under non-design conditions (path

*ABCD*), a ray is emitted at a point

*A*within the specimen. Thus it first traverses a specimen layer of thickness

*t*

_{s}, which in our notation corresponds to

*z*

_{p}and refractive index

*n*

_{s}, before passing through a coverslip of thickness

*t*

_{g}and refractive index

*n*

_{g}. Then it passes the immersion medium of thickness

*t*

_{i}and refractive index

*n*

_{i}before entering the objective. The first and second derivatives of the PSF with respect to the axial position

*z*

_{p}, required in Eqs. (7)–(10), are as follows:

$$\phantom{\rule{6.2em}{0ex}}\xb7{\int}_{0}^{1}\mathrm{cos}\left(W(\rho ,{z}_{n}\mid {z}_{p})\right){J}_{0}\left(\mathit{kr}\mathrm{NA}\rho \right)g\left(\rho \right)\rho \phantom{\rule{.2em}{0ex}}d\rho $$

$${\phantom{\rule{5.2em}{0ex}}-2k\mid A\mid}^{2}{\int}_{0}^{1}\mathrm{cos}\left(W(\rho ,{z}_{n}\mid {z}_{p})\right){J}_{0}\left(\mathit{kr}\mathrm{NA}\rho \right)\rho \phantom{\rule{.2em}{0ex}}d\rho ,$$

$$\phantom{\rule{6.2em}{0ex}}\xb7{\int}_{0}^{1}\mathrm{sin}\left(W(\rho ,{z}_{n}\mid {z}_{p})\right){J}_{0}\left(\mathit{kr}\mathrm{NA}\rho \right)g\left(\rho \right)\rho \phantom{\rule{.2em}{0ex}}d\rho ,$$

$$\phantom{\rule{5.2em}{0ex}}-2{k}^{2}{\mid A\mid}^{2}{\int}_{0}^{1}\mathrm{sin}\left(W(\rho ,{z}_{n}\mid {z}_{p})\right){J}_{0}\left(\mathit{kr}\mathrm{NA}\rho \right)\rho \phantom{\rule{.2em}{0ex}}d\rho $$

$$\phantom{\rule{6.2em}{0ex}}\xb7{\int}_{0}^{1}\mathrm{sin}\left(W(\rho ,{z}_{n}\mid {z}_{p})\right){J}_{0}\left(\mathit{kr}\mathrm{NA}\rho \right)g{\left(\rho \right)}^{2}\rho \phantom{\rule{.2em}{0ex}}\mathrm{d\rho}$$

$$\phantom{\rule{5.2em}{0ex}}+2{k}^{2}{\mid A\mid}^{2}{\left({\int}_{0}^{1}\mathrm{cos}\left(W(\rho ,{z}_{n}\mid {z}_{p})\right){J}_{0}\left(\mathit{kr}\mathrm{NA}\rho \right)g\left(\rho \right)\rho \phantom{\rule{.2em}{0ex}}d\rho \right)}^{2}$$

$$\phantom{\rule{5.2em}{0ex}}-2{k}^{2}{\mid A\mid}^{2}{\int}_{0}^{1}\mathrm{cos}\left(W(\rho ,{z}_{n}\mid {z}_{p})\right){J}_{0}\left(\mathit{kr}\mathrm{NA}\rho \right)\rho \phantom{\rule{.2em}{0ex}}d\rho $$

$$\phantom{\rule{6.2em}{0ex}}\xb7{\int}_{0}^{1}\mathrm{cos}\left(W(\rho ,{z}_{n}\mid {z}_{p})\right){J}_{0}\left(\mathit{kr}\mathrm{NA}\rho \right)g{\left(\rho \right)}^{2}\rho \phantom{\rule{.2em}{0ex}}d\rho ,$$

where

Note that Haeberlé [23] proposed an easy way of evaluating the vectorial model by Török et al. and Hell *et al*., based on the parameters and formulation of the OPD of the Gibson and Lanni model.

## Acknowledgments

We would like to thank the reviewer for providing us with insightful comments that led to the clarification of some essential aspects of the paper. Many thanks also to Nathalie Garin for her help with the preparation of the samples used in the experiments.

## References and links

**
1
. **
In practice, the number of possible optical sections is constrained by the exposure time, the dynamics of the biological process under study, and photobleaching of the fluorescent labels.

**
2
. **
S.
Hell
,
G.
Reiner
,
C.
Cremer
, and
E. H. K.
Stelzer
, “
Aberrations in confocal fluorescence microscopy induced by mismatches in refractive index
,”
J. Microsc.
**
169
**
,
391
–
405
(
1993
). [CrossRef]

**
3
. **
U.
Kubitscheck
, “
Single Protein Molecules Visualized and Tracked in the Interior of Eukaryotic Cells
,”
Single Molecules
**
3
**
,
267
–
274
(
2002
). [CrossRef]

**
4
. **
M. G. L.
Gustafsson
,
D. A.
Agard
, and
J. W.
Sedat
, “
I5M: 3D widefield light microscopy with better than 100nm axial resolution
,”
J. Microsc.
**
195
**
,
10
–
16
(
1999
). [CrossRef] [PubMed]

**
5
. **
S.
Hell
and
E. H. K.
Stelzer
, “
Fundamental improvement of resolution with a 4Pi-confocal fluorescence microscope using two-photon excitation
,”
Opt. Commun.
**
93
**
,
277
–
282
(
1992
). [CrossRef]

**
6
. **
S.
Hell
and
J.
Wichmann
, “
Breaking the diffraction resolution limit by stimulated emission: stimulated-emission-depletion fluorescence microscopy
,”
Opt. Lett.
**
19
**
(11),
780
–
782
(
1994
). [CrossRef] [PubMed]

**
7
. **
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
**
,
1291
–
1300
(
1994
). [CrossRef] [PubMed]

**
8
. **
N.
Bobroff
, “
Position measurement with a resolution and noise-limited instrument
,”
Rev. Sci. Instrum.
**
57
**
,
1152
–
1157
(
1986
). [CrossRef]

**
9
. **
M. K.
Cheezum
,
W. F.
Walker
, and
W. H.
Guilford
, “
Quantitative Comparison of Algorithms for Tracking Single Fluorescent Particles
,”
Biophys J.
**
81
**
,
2378
–
2388
(
2001
). [CrossRef] [PubMed]

**
10
. **
R. E.
Thompson
,
D. R.
Larson
, and
W. W.
Webb
, “
Precise Nanometer Localization Analysis for Individual Fluorescent Probes
,”
Biophys J.
**
82
**
,
2775
–
2783
(
2002
). [CrossRef] [PubMed]

**
11
. **
R. J.
Ober
,
S.
Ram
, and
S.
Ward
, “
Localization Accuracy in Single-Molecule Microscopy
,”
Biophys J.
**
86
**
,
1185
–
1200
(
2004
). [CrossRef] [PubMed]

**
12
. **
A.
Van Oijen
,
J.
Köhler
,
J.
Schmidt
,
M.
Müller
, and
G.
Brakenhoff
, “
3-Dimensional super-resolution by spectrally selective imaging
,”
Chem. Phys. Lett.
**
292
**
,
183
–
187
(
1998
). [CrossRef]

**
13
. **
M.
Speidel
,
A.
Jonas
, and
E.-L.
Florin
, “
Three-dimensional tracking of fluorescent nanoparticles with sub-nanometer precision by use of off-focus imaging
,”
Opt. Lett.
**
28
**
(2),
69
–
71
(
2003
). [CrossRef] [PubMed]

**
14
. **
N.
Subotic
,
D. Van
De Ville
, and
M.
Unser
, “
On the Feasibility of Axial Tracking of a Fluorescent Nano-Particle Using a Defocusing Model
,” in
*
Proceedings of the Second IEEE International Symposium on Biomedical Imaging: From Nano to Macro (ISBI’04)
*
, pp.
1231
–
1234
(
Arlington VA, USA
,
2004
).

**
15
. **
P.
Török
,
R.
Varga
,
Z.
Laczik
, and
R.
Booker
, “
Electromagnetic diffraction of light focused through a planar interface between materials of mismatched refractive indices: an integral representation
,”
J. Opt. Soc. Am. A
**
12
**
(2),
325
–
332
(
1995
). [CrossRef]

**
16
. **
A.
Egner
and
S. W.
Hell
, “
Equivalence of the Huygens-Fresnel and Debye approach for the calculation of high aperture point-spread functions in the presence of refractive index mismatch
,”
J. Microsc.
**
193
**
,
244
–
249
(
1999
). [CrossRef]

**
17
. **
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
**
8
**
(10),
1601
–
1613
(
1991
). [CrossRef]

**
18
. **
M.
Gu
,
*
Advanced Optical Imaging Theory
*
(
Springer
,
2000
).

**
19
. **
C.
Preza
and
J.-A.
Conchello
, “
Depth-variant maximum-likelihood restoration for three-dimensional fluorescence microscopy
,”
J. Opt. Soc. Am. A
**
21
**
(9),
1593
–
1601
(
2004
). [CrossRef]

**
20
. **
V.
Levi
,
Q.
Ruan
, and
E.
Gratton
, “
3-D Particle Tracking in a Two-Photon Microscope: Application to the Study of Molecular Dynamics in Cells
,”
Biophys J.
**
88
**
,
2919
–
2928
(
2005
). [CrossRef] [PubMed]

**
21
. **
P.
Török
and
R.
Varga
, “
Electromagnetic diffraction of light focused through a stratified medium
,”
Applied Optics
**
36
**
(11),
2305
–
2312
(
1997
). [CrossRef] [PubMed]

**
22
. **
A. S.
Marathay
and
J. F.
McCalmont
, “
On the usual approximation used in the Rayleigh-Sommerfeld diffraction theory
,”
J. Opt. Soc. Am. A
**
21
**
(4),
510
–
516
(
2004
). [CrossRef]

**
23
. **
O.
Haeberlé
, “
Focusing of light through a stratified medium: a practical approach for computing microscope point spread functions. Part I: Conventional microscopy
,”
Opt. Commun.
**
216
**
,
55
–
63
(
2002
). [CrossRef]