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

Authors: Joshua W. Shaevitz
Journal: International Journal of Optics
 (2009)

---

## Abstract

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.


---

## Body

## 1. 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–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–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–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 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].


## 2. 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 thex 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:(1)I(x,y)=I0+Ae−((x−x0)/2σx2)−((y−y0)/2σy2),where I0 is the
background offset, A is the peak
intensity, (x0,y0) is the lateral
center of the peak, and (σx,σy) are the lateral
widths of the peak.The probability functionp(σ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.


## 3. Bayesian Estimation of the Axial Position

We wish to estimate the most likelyz-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 z-positions.
Bayes' theorem allows us to relate the posterior probability distribution to
the calibration data via(2)p(z∣σx,σy)=p(σx,σy∣z)p(z)p(σx,σy),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)=∑zp(σ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 zB, using the posterior probability(3)zB=∑zz×p(z∣σx,σy).Because this approach is derived from knowledge of the
probability distribution, we can estimate the error onzB in several
ways. One method is to calculate the standard
deviation (SD) via(4)SD=∑z(z−zB)2×p(z∣σx,σy).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.


## 4. Results

Figure1(a) shows a typical bead image at 5 different z-positions 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.Variation in the lateral bead profile with axial
position. (a) Five images of the same 57 nm beam taken at different axial
positions demonstrate the astigmatism. (b) Thex- and y-widths (x: red +, y: blue ×) for each
height calculated by fitting the bead image to (1). Solid lines are a fit of
the width data to the defocusing curve described in (5). (c) Two-dimensional
width plot of the same data plotted in (b). 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).
(a)(b)(c)By binning the data in Figure1(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(b)). This 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.Calibration probability distributions. (a) 2D
plots ofp(σ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.
(a)(b)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.Table 1
Localization results for points A–D from Figure1(c). zB: position estimate using the Bayesian algorithm, zd: position estimate using the method of Huang et al.
[12], CI: credibility
interval, SD: standard deviation.
PointzB (nm)zd (nm)80% CI (nm)SD (nm)p(σx,σy)×10−4A441450388–474312.9B243236176– 3436019C10611362–176417.0Posterior probability distributions for points A–D
from Figure1(c).
(a)(b)While the posterior distributions for points A–C are
easily interpreted (see Figure3(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.


## 5. 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–12]. One of the most recent
examples was presented by Huang et al. [12] for use in localization-based 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.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)(5)σ(z)=σ01+(z−z0d)2+A(z−z0d)3+B(z−z0d)4.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(6)D=(σx1/2−σx,fit1/2)2+(σy1/2−σy,fit1/2)2.The z-position that
corresponds to the minimum of this distance function is chosen as the true
position, referred to herein as zd. The distance functions for points A–D are shown in
Figure 4, and the calculated zd for each point
is listed in Table 1.Distance plots using the method of Huang et al.
[12] for points A–D
using the calibration data displayed in Figure 1(c).
(a)(b)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
calculatingzd. 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. Figure5(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.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)zB−zd for different
points in the (σx,σy) plane.
Differences are in nm (color scale).
(a)(b)The Bayesian method does slightly worse at the very
lowest and highestz-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 betweenzB and zd 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.


## 6. 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-numerical-aperture 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.



---

*Source: 896208-2009-03-10.xml*