Bayesian Estimation of the Axial Position in Astigmatism-Based Three-Dimensional Particle Tracking

Accurate estimation of the axial position of a molecule using a single lateral image remains a challenge in fluorescent single particle tracking. Here, a principled algorithm for the Bayesian estimation of the axial position of a molecule in three-dimensional astigmatism-based particle tracking is proposed. This technique uses the data from a calibration image set to derive the position without assuming a functional form for the abberated defocusing curve. Using a calibration image set from forty 57 nm beads, the axial position is calculated, and the error associated with position estimation is discussed. This method is compared to previously published algorithms.


Introduction
Light microscopy using fluorescent probes has led to numerous discoveries in biology and has become a standard tool for biological scientists.In the last several decades, techniques that monitor the position of a single fluorescent molecule, herein referred to as single particle tracking (SPT), have found applications from quantifying the motion of single molecular motors to analyzing the spatial distribution of membrane associated complexes in cells; see, for example, [1,2].More recently, SPT has been combined with fluorescent photoswitching to achieve a localization-based imaging technique that allows for the recording of microscopy images at better than diffraction limited resolution [3][4][5].
While SPT in two dimensions is somewhat straightforward, three-dimensional (3D) localization of a fluorescent molecule from a single lateral-section image has remained difficult.Some of these methods use information from two different focal planes [6][7][8].The majority of these algorithms, however, perform axial localization from a single image and use the change in lateral width of the fluorescent image with axial position.These include techniques that estimate position using maximum-likelihood methods [9].However, due to the approximate axial symmetry of the 3D point spread function (PSF), it is often difficult to distinguish particles above the object plane from those below it.As a result, several methods have been proposed to break this symmetry [10][11][12][13].
One class of popular techniques uses a cylindrical optic to separate the axial focus of the two lateral dimensions creating an elliptical PSF.These astigmatism-based techniques are straightforward to implement by adding a single optical element in the imaging path.However, the estimation of the axial position from these techniques is not straightforward.Indeed, a number of different techniques have been proposed for estimation of the axial position.All of these techniques use a user-defined parameterization of the defocusing curve which is then fit to a calibration data set [10,11].These functions, which are at best based on an idealized model of the imaging system, often include ad hoc terms to account for spherical and chromatic aberrations that are not present in the model.While these methods often yield satisfactory localization results, they are prone to systematic errors associated with the choice of model and aberration parameterization and make estimation of the localization error difficult.
An alternative methodology is to use the calibration data as representative of the population of all possible International Journal of Optics molecular images from different heights.This data set can then be thought of as an estimate of the probability distribution that relates parameters of the measured image to the axial position.The approach presented here uses Bayesian estimation to calculate the expectation value of the axial position given this probability distribution and the measured lateral widths of an elliptical image.This procedure generates a likelihood function of the axial position given the data, which allows for a statistical understanding of the uncertainty in the localization.Information theoretical has been previously used to discuss the limits of lateral particle localization [14].

Materials and Methods
Experimental measurements were carried out on a Nikon (Tokyo, Japan) TE2000 inverted microscope.Fluorescence excitation was achieved by focusing a 473 nm diode laser (Dragon Lasers (ChangChun, China)) onto the back aperture of a 100x magnification, 1.49 N.A. objective lens (Nikon (Tokyo, Japan)).Emission light was projected onto an EMCCD camera (Andor Technology, Belfast, Northern Ireland) using a 2.5x projection lens, yielding a pixel size of approximately 80 nm in the recorded image.To create an astigmatism, a 1 m focal length cylindrical lens was placed between the projection lens and the camera, which resulted in a focal separation of ∼100 nm between the x and y axes.The sample was mounted on a two-axis piezo stage while the objective lens was moved via a separate piezo drive (Mad City Labs, Madison, Wis, USA).Data acquisition software to control the stages, camera, and excitation laser was written in LabView (National Instruments, Austin, Tex, USA).All data and image analyses were performed in Igor Pro (Wavemetrics, Lake Oswego, Ore, USA).
For each field of fluorescent molecules to be analyzed, a series of images was recorded at different axial positions.For these experiments, 57 nm Dragon green fluorescent beads (Bangs Labs) were adhered to a coverglass surface in a fluid cell, similar to that described previously [15].Between each image in the series, the objective was moved by 10 nm.This resulted in an axial separation of the images of 7.5 nm due to the focal shift inherent in the imaging system.This value is measured using the optical trap as an interferometer as described in [15,16].All axial positions in this paper have been corrected for the focal shift.A total of 40 individual beads were analyzed.
To calculate the width of a fluorescent spot in the two lateral dimensions, each bead image was fit to a two-dimensional elliptical Gaussian using a Levenberg-Marquardt algorithm to minimize the chi-square: where I 0 is the background offset, A is the peak intensity, (x 0 , y 0 ) is the lateral center of the peak, and (σ x , σ y ) are the lateral widths of the peak.The probability function p(σ x , σ y | z), see Section 3, was created by binning the data at each z-value in 0.04 pixel bins in (σ x , σ y ) space.To ensure a continuous probability function and to compensate for a finite data set, we smoothed p(σ x , σ y | z) for each z with a Gaussian kernel 15 bins in width, yielding a smoothed function with 0.6 pixel resolution.

Bayesian Estimation of the Axial Position
We wish to estimate the most likely z-position of a molecule given a pair of image widths, (σ x , σ y ).This procedure amounts to calculating the posterior probability, p(z | σ x , σ y ), given a calibration image set taken at known zpositions.Bayes' theorem allows us to relate the posterior probability distribution to the calibration data via where p(σ x , σ y | z) is the probability distribution measured from the calibration image set.For this work, we will assume a flat prior probability distribution, p(z).The normalization term, p(σ x , σ y ) = z p(σ x , σ y | z), represents the probability of measuring a pair of widths (σ x , σ y ) summed over all heights.Using a Bayesian approach, we estimate the expectation value of z, herein referred to as z B , using the posterior probability Because this approach is derived from knowledge of the probability distribution, we can estimate the error on z B in several ways.One method is to calculate the standard deviation (SD) via Another method of reporting the error is to use credibility intervals (CIs) which specify a posterior probability interval.For example, an 80% credibility interval for the parameter z of a-b means that the posterior probability that z lies in the interval from a to b is 0.8.

Results
Figure 1(a) shows a typical bead image at 5 different zpositions and demonstrates the astigmatism.As expected, the x-width is minimized at a different z-position than the y-width due to the presence of the cylindrical lens in the imaging path (see Figure 1(b)).Our estimation of the axial position of a particle uses the two lateral widths and attempts to infer the z-position.It is therefore convenient to display the data with the two lateral widths as the independent variables and the z-position as the dependent variable.Figure 1(c) shows the data of 1b plotted in this way.By binning the data in Figure 1(c), we can estimate the probability function p(σ x , σ y | z) as shown in Figure 2(a).In addition, we can calculate the total probability distribution of measuring a pair of widths, p(σ x , σ y ) (see Figure 2 Color scale represents the z-position.Letters denote the location of points in the (σ x , σ y ) plane analyzed and tabulated in Table 1.The solid line is data from the same fit as for (b).function is useful for evaluating the validity of a measured pair of widths in light of the calibration data set and can be used to identify spurious data during an experiment.We chose four points in the σ x , σ y plane to demonstrate this algorithm (labeled A-D in Figure 1(c)).Points A-C generally represent positions below, at, and above the focal plane, respectively.Point D does not lie within the range spanned by the calibration data set and is an example of a spurious data point.Figure 3 shows the posterior probability distributions for each of these points.The expectation value, 80% credibility interval, the standard deviation, and p(σ x , σ y ) of points A-C are shown in Table 1.Credibility intervals are calculated symmetrically about the expectation value with respect to the posterior probability.
While the posterior distributions for points A-C are easily interpreted (see Figure 3(a)), the distribution for point D is difficult to interpret on its own because this pair of lateral widths is not well represented by the calibration data set (see Figure 3(b)).However, because p(σ x , σ y ) for point D is 0.27 × 10 −4 , more than an order of magnitude lower than the values from points A-C, this point can be readily identified as an outlier.

Comparison with Previous Techniques for Axial Estimation
Several techniques have been used previously for axial estimation which rely on a parameterized fit of the microscope defocusing curve [10][11][12].One of the most recent examples was presented by Huang et al. [12] for use in localizationbased microscopy.Those authors use a fit to the change in lateral width as a function of z-position and then calculate, for a given (σ x , σ y ) pair, the closest point within the fit in the (σ x , σ y )-plane.Figure 4: Distance plots using the method of Huang et al. [12] for points A-D using the calibration data displayed in Figure 1(c).
To compare the Bayesian estimation described here to this algorithm, we fit our data with the same function used in [12] (see Figure 1, solid lines) ( For a given (σ x , σ y ), a modified geometric distance in the (σ x , σ y )-plane from this point to the fit value for each z-position is calculated using The z-position that corresponds to the minimum of this distance function is chosen as the true position, referred to herein as z d .The distance functions for points A-D are shown in Figure 4, and the calculated z d for each point is listed in Table 1.This fit-based technique has several drawbacks.First, the defocusing curved used in the fit is chosen by hand, and corrections for aberrations are added ad hoc as high-order polynomial terms.Systematic deviations of the data from this curve will give rise to systematic errors in calculating z d .Additionally, while the value of the minimum of the distance function and the width of this minimum are related to the localization error, they do not give a probabilistic estimate of the error that corresponds to the distribution of widths measured at different z-positions.Finally, the evaluation of spurious data using the distance function, such as Figure 4(b), is difficult.
In order to compare these two methods, we used two different metrics.Figure 5(a) shows a jackknife analysis of the two methods performed on the entire calibration data set.Briefly, for each of the 40 calibration beads, the difference between the estimated and actual z-positions was calculated using the other 39 beads for calibration of each of the two methods.The average and standard deviations of these errors are plotted as a function of z-position.By this measure, the two techniques do equally as well.As expected, the standard deviations calculated from the probability function (∼50 nm, see Table 1) are similar in magnitude to those derived from the jackknife analysis.These errors are distinct from those measured by localizing a single particle many times as they include both measurement uncertainty and the uncertainty caused by the distribution of widths from different particles.
The Bayesian method does slightly worse at the very lowest and highest z-positions where the posterior distribution is not fully sampled.A more complex algorithm that fits the shape of the posterior probability distribution removes this error and recovers a similar residual to the distance method.In general, the calibration data at the edges is under represented for both methods, and therefore it is best to limit experimental samples to a smaller range than the calibrated one.
We additionally examined the difference between z B and z d for a range of points in the (σ x , σ y ) plane (see Figure 5(b)).The distance method deviates from the Bayesian value by up to 20 nm, with regions of systematic difference between the two techniques.

Conclusion
Our Bayesian method of axial position estimation offers two main advantages over previous techniques.First, it allows a probabilistic estimate of the localization error not afforded by previously described techniques.Second, it facilitates the identification of spurious data points that might otherwise be included in data analysis.
While the algorithm presented here describes a Bayesian method for estimating axial particle position, it does not attempt to correct for the effect of optical aberrations on particle localization.We find that using a high-numericalaperture objective, a large systematic error is induced in axial localization above the focal plane [17].Further studies will be necessary to develop a robust method to correct for these errors.
Single-particle tracking has become an important technique in numerous fields of biology.The recent implementation of these techniques for super-resolution 3D imaging has the potential to extend our knowledge of cellular processes on the nanometer level.The axial-localization technique presented here should prove useful as these imaging methods become more quantitative and less descriptive.

Figure 1 :
Figure1(a) shows a typical bead image at 5 different zpositions and demonstrates the astigmatism.As expected, the x-width is minimized at a different z-position than the y-width due to the presence of the cylindrical lens in the imaging path (see Figure1(b)).Our estimation of the axial position of a particle uses the two lateral widths and attempts to infer the z-position.It is therefore convenient to display the data with the two lateral widths as the independent variables and the z-position as the dependent variable.Figure1(c) shows the data of 1b plotted in this way.By binning the data in Figure1(c), we can estimate the probability function p(σ x , σ y | z) as shown in Figure2(a).In addition, we can calculate the total probability distribution of measuring a pair of widths, p(σ x , σ y ) (see Figure2(b)).This

Figure 2 :
Figure 2: Calibration probability distributions.(a) 2D plots of p(σ x , σ y | z) at three different axial positions.(b) Total probability function p(σ x , σ y ) of measuring a pair of widths (σ x , σ y ).The lateral size of each image is 5 pixels.

Figure 3 :
Figure 3: Posterior probability distributions for points A-D from Figure 1(c).

Figure 5 :
Figure 5: Comparison of the Bayesian and distance estimation methods.(a) Residual errors calculated from a jackknife analysis of the calibration data (see text, Bayesian method: red circles, distance method: black squares).Error bars are standard deviations.(b) z B − z d for different points in the (σ x , σ y ) plane.Differences are in nm (color scale).

Table 1 :
[12]lization results for points A-D from Figure1(c).z B : position estimate using the Bayesian algorithm, z d : position estimate using the method of Huang et al.[12], CI: credibility interval, SD: standard deviation.