Localization of a fluorescent source without numerical fitting

We present an algebraic solution to the problem of localizin g a single fluorescent particle with sub-diffraction-limit a ccuracy. The algorithm is derived and its performance studied experimentall y. Isolated 20 nm fluorescent beads were imaged using a wide-field microscop e at two different positions separated by 100 nm and at a range of sign al-to-noise ratios (SNR). The data were analyzed using both the new algor ithm and the standard approach of fitting the data to a Gaussian profile. Re sults indicate that the proposed approach is nearly as accurate as Gaussian fitting across a wide range of SNR while executing over 200 times faster. In ad dition, the new algorithm is able to localize at lower SNR than the fitting method. © 2008 Optical Society of America OCIS codes: (180.2520) Fluorescence microscopy; (300.6280); Spectro scopy, fluorescence and luminescence; (100.2960) Image analysis . References and links 1. M. Lakadamyali, M. J. Rust, H. P. Babcock, and X. Zhuang, “V isualizing infection of individual influenza viruses,” P. Natl. Acad. Sci USA100, 9280–9285 (2003). 2. C. Kural, H. Balci, and P. R. Selvin, “Molecular motors one at a time: FIONA to the rescue,” J. Phys. Condens. Matt. 27, S3979–S3995 (2005). 3. C. K. Payne, “Imaging gene delivery with fluorescence micr os opy,” Nanomedicine2, 847–860 (2007). 4. W. E. Moerner, “New directions in single-molecule imagin g and analysis,” P. Natl. Acad. Sci USA 104, 12,596– 12,602 (2007). 5. C. Joo, H. Balci, Y. Ishitsuka, C. Buranachai, and T. Ha, “A dvances in single-molecule fluorescence methods for molecular biology,” Annu. Rev. Biochem. 77, 51–76 (2008). 6. M. K. Cheezum, W. F. Walker, and W. H. Guilford, “Quantitat ive comparison of algorithms for tracking single fluorescent particles,” Biophys. J. 81, 2378–2388 (2001). 7. A. J. Berglund, M. D. McMahon, J. J. McClelland, and J. A. Li ddle, “Fast, bias-free algorithm for tracking single particles with variable size and shape,” Opt. Express 16, 14,064–14,075 (2008). 8. R. E. Thompson, D. R. Larson, and W. W. Webb, “Precise nanom eter localization analysis for individual fluorescent probes,” Biophys. J. 82, 2775–2783 (2002). 9. R. J. Ober, S. Ram, and E. S. Ward, “Localization accuracy i n single-molecule microscopy,” Biophys. J. 86, 1185–1200 (2007). 10. S. Ram, E. S. Ward, and R. J. Ober, “A stochastic analysis o f performance limits for optical microscopes,” Multdim. Syst. Sign. P. 17, 27–57 (2006). 11. A. J. Berglund and H. Mabuchi, “Tracking-FCS: Fluoresce nce correlation spectroscopy of individual particles,” Opt. Express13, 8069–8082 (2005). 12. S. B. Andersson, “Tracking a single fluorescent molecule with a confocal microscope,” Appl. Phys. B 80, 809– 816 (2005). 13. V. Levi, Q. Ruan, and E. Gratton, “3-D particle tracking i a two-photon microscope: application to the study of molecular dynamics in cells,” Biophys. J. 88, 2919–2928 (2005). 14. S. Bancroft, “An algebraic solution of the GPS pseudoran ge equations,” IEEE T. Aero. Elec. Sys. AES-21, 56–59 (1985). #100791 $15.00 USD Received 29 Aug 2008; revised 9 Oct 2008; accepted 24 Oct 2008; published 29 Oct 2008 (C) 2008 OSA 10 November 2008 / Vol. 16, No. 23 / OPTICS EXPRESS 18714 15. T. Sun and S. B. Andersson, “Precise 3-D localization of fl uorescent probes without numerical fitting,” in Proceedings of the International Conference of IEEE Engineeri ng in Medicine and Biology Society (IEEE, 2007) pp. 4181–4184. 16. S. B. Andersson, “Precise localization of fluorescent pr obes without numerical fitting,” inProceedings of IEEE International Symposium on Biomedical Imaging (ISBI) (IEEE, 2007), pp. 252–255. 17. S. B. Andersson, “Position estimation of fluorescent pro bes in a confocal microscope,” in Proceedings of IEEE Conference on Decision and Control (IEEE, 2007), pp. 4950–4955. 18. R. Juškaitis, “Measuring the real point spread functio of high numerical aperture microscope objective lenses,” in Handbook of Biological Confocal Microscopy , J. B. Pawley, ed. (Springer, 2006), pp. 239-250. 19. A. Ben-Israel and T. N. E. Greville, Generalized Inverses: Theory and Applications , Pure and Applied Mathematics (John Wiley and Sons, 1974). 20. A. Papoulis,Probability, Random Variables, and Stochastic Processes (McGraw-Hill, 1991).


Introduction
The ability to determine the position of an isolated fluorescent particle with sub-diffraction limit accuracy is an invaluable tool in fluorescence microscopy.It has been used in the study of a wide variety of molecular processes including the motion of influenza viruses during infection [1], the motion of molecular motors such as myosin and kinesin [2], and the delivery of genes to the nucleus using synthetic delivery vectors [3].The interested reader is directed to the review articles of Moerner [4] and Joo, et.al. [5] for further references.
Many different techniques for estimating the position of the particle have been proposed (see, e.g., [6,7]).The most commonly used ones typically fit the measured intensity data to a Gaussian profile.The location corresponding to the maximum of the Gaussian is taken to be the position of the source.Theoretical formulas for the accuracy of this approach have been derived [8] and, with sufficiently high signal-to-noise ratio (SNR), accuracy on the order of one nanometer can be achieved [2].Fundamental limits on the accuracy of estimation have been determined and algorithms based on optimal estimation theory proposed that nearly achieve this limit [9,10].
Despite the myriad of successes to date, these methods suffer from drawbacks arising from the numerical nature of the schemes.These include the need for a large set of data and the time it takes to perform a numerical fit of the data.Neither of these is an issue when imaging is performed in wide-field using a CCD camera and when analysis is done offline.In applications where real-time position information is desired, however, it is important to produce estimates as fast as possible.Moreover, when images are acquired using confocal or multi-photon techniques, pixels are obtained sequentially rather than in parallel.The amount of data required to determine accurately the position of a fluorescent particle then has a direct effect on the temporal resolution of any real-time tracking algorithm utilizing position estimation.These issues are particularly relevant given the increasing interest in using confocal and multi-photon setups to track the motion of single particles [11,12,13].
In this paper we describe a novel algebraic algorithm inspired by a solution to the position estimation problem in the Global Positioning System (GPS) known as Bancroft's algorithm [14].The algorithm, termed fluoroBancroft, relies on the fact that the fluorescence intensity depends only on the distance between the location of the measurement (the center of a pixel when a CCD array is used) and the position of the source particle.We focus here on localization in the plane, although the algorithm has been generalized to 3-D [15].Simulation results to date indicate that the method has accuracy similar to Gaussian fitting [16] when large amounts of data are available and performs significantly better than Gaussian fitting when only a few measurements are available [17].In fact as few as three measurements are needed to estimate the position of a source particle in the plane.
The remainder of the paper is organized as follows.In Sec. 2 we describe the fundamental idea of the algorithm in terms of range-based localization and then give the mathematical derivation of the estimation formula in Sec. 3. In Sec. 4 we describe a set of experiments in which 20 nm fluorescent microspheres were imaged at two successive locations separated by 100 nm.A collection of images at different SNR were taken and the data were used to study the performance of the fluoroBancroft algorithm with respect to the Gaussian fitting approach.These results are described in Sec. 5.

FluoroBancroft -the concept
As will be shown in Sec. 3 below, a measurement of the fluorescence intensity can be converted to a measurement of the distance between the fluorescent particle and the position at which the measurement was taken.If a CCD camera is being used, each pixel acts as an individual detector with its position given by the center of the pixel.In confocal and multi-photon applications, the center of the detection volume of the microscope determines the position of the measurement.As illustrated in Fig. 1, a measurement of the range yields a circle of possible locations of the source (or, if estimation is being done in three dimensions, a sphere of possible locations).A second measurement from a different location in space yields a second circle.Under the assumption that the fluorescent particle is fixed, it must be located at the intersection of the two circles.With two measurements there are thus two possible positions for the source.Including a third measurement from a third location produces yet another circle.In the absence of noise, there is a single point of intersection of the three circles and thus a unique solution for the position of the source particle (in three dimensions, a fourth measurement is needed to get a unique solution).Noise in the intensity measurements, however, is propagated into noise in the measurement of the range.In general, then, the circles will not intersect at a common point and an approximate, rather than an exact, solution to the estimation problem must be found.In Sec. 3 below, we develop a linear system in terms of these range estimates (c.f.(10)) and present an analytical solution using the Moore-Penrose generalized inverse.Due to the properties of this inverse, the solution satisfies a least-squares criterion in terms of the error (see (12)).Fig. 1.Range-based estimation in the plane.Each measurement yields an estimate of the range to the source (yellow star) and thus a circle of possible locations for that source.With two measurements from two different locations, the source must lie at the intersection of the two circles, leading to two possible solutions.In the absence of noise a third measurement produces a unique solution given by the single common point of intersection of the circles.

FluoroBancroft -the derivation
We model the intensity at (x, y) of a fluorescent particle located at (x o , y o ) as a Gaussian profile given by where m is an unknown scaling factor determined by the photon emission rate of the fluorescent particle and the integration time of the measurement, σ x and σ y describe the width of the point spread function in the two directions, η shot is a random variable describing the shot noise, and η B is a random variable with mean N B that captures the background intensity noise and the intrinsic noise (dark current) of the detector.The value of N B is assumed known (through measurement).The random variable defined by is assumed to be Poisson with mean value where • represents expectation.The parameters σ x and σ y are determined by setting the e −1 point of the intensity model to be equal to the Rayleigh radius.In the absence of aberrations, the values will be equal and given by where N.A. is the numerical aperture of the objective lens and λ is the wavelength of the emitted light.If σ x is not equal to σ y then we scale the y−coordinate by defining ỹ = (σ y /σ x ).We then have x + η B + η shot (5) where we have defined the distance r between the point (x, y) and the position of the molecule Note that this model is an approximation of the true point spread function of a diffraction limited spot and the noise in the measurement (see, e.g., [18]).Deviations of the measured intensity from this model will introduce error in the estimated position of the source particle.
Taking the expected value of ( 5) and rearranging, we obtain an equation for the range to the source particle: Note that the value of m, related to the true intensity of the fluorophores, is not typically known.
We now assume we have a collection of intensity measurements, denoted as I i , and taken from the (known) positions (x i , y i ).Each measurement yields a range measurement of the form (7). To rearrange this equation, define b = 2σ 2 x ln(m), Then ( 7) can be written as the following scalar equation: Stacking together N such measurements into a vector equation and rearranging yields where are all known.The N × 3 matrix B captures the locations at which each of the measurements was taken.Under the mild assumption that all the points do not lie on a single line, this matrix is full rank.
The linear system defined by ( 10) is overdetermined and thus in general does not have an exact solution.An approximate solution that minimizes the Euclidean norm of the residual error, δ , defined by is given by the Moore-Penrose pseudo-inverse of B [19].To obtain this least-squares solution, pre-multiply ( 10) by B T and rearrange to yield where The unknown position of the fluorescent particle appears on both sides of this equation, explicitly on the left side and implicitly on the right side through Λ.To solve for the position of the fluorescent source, we utilize the following proposition.
Proposition 3.1.Let e be the N × 1 vector of all ones as defined in (11) and let A be any N × m matrix such that the matrix obtained by appending e to A, B △ = (A e), is full rank.Then Proof.Let I denote the identity matrix.We have Thus B † e is equal to the last column of the identity matrix as claimed.
Prop. 3.1 allows us to solve for the position of the fluorescent particle as follows.Define the matrix Q as Pre-multiplying ( 13) by Q and applying the proposition yields the fluoroBancroft solution: Given a set of N measurements, the fluoroBancroft algorithm is applied by first building the vector α and the matrix B defined in (11) , calculating the Moore-Penrose inverse B † , and then performing the matrix multiplication in (16).Although the matrix B T B must be inverted, this matrix is always 3×3, independent of the number of measurements used.

Experimental methods
All experiments were carried out on a Zeiss Axiovert 200 microscope using a 63x, 1.2 N.A. water immersion objective lens (Carl Zeiss, Inc., Thornwood, NY, USA).Actuation was achieved using a three-dimensional nanopositioning stage (Nano-PDQ; Mad City Labs, Inc., Madison, WI, USA) with a positioning accuracy of 0.1 nm.An Xcite EXFO illumination source (Mississauga, Ontario, Canada) was used to generate the excitation light.This light was passed through a dichroic filter set (Chroma Technology, Rockingham, VT, USA) to produce a narrow-band excitation centered at 532 nm.Fluorescence emission collected by the objective was passed through the filter set that blocked any reflected excitation light.The image was captured using a CCD camera (Retiga EXi; QImaging, Inc., Surrey, BC, Canada).Each pixel on the CCD sensor array was 6.45 µm × 6.45 µm; this translates to 102 nm × 102 nm in the image plane with the 63x objective.
20-nm diameter carboxylate-modified microspheres embedded with "Nile Red" fluorophores were purchased (Invitrogen, Inc, Carlsbad, CA).The microspheres were diluted in de-ionized water and dried onto a glass coverslip.The coverslip was then placed on a bead of water on a glass slide and sealed.A region containing four visible and isolated beads was selected.A cropped image of one of the beads taken with an exposure time of 25 ms is shown in Fig. 2. Images were taken with exposure times of from 1 to 51 ms in steps of 2 ms.At each exposure time, forty sets of two images were collected.In between each image the nanostage was displaced 100 nm in the x−direction.The background statistics were determined by selecting a 30 × 30 pixel region that did not contain a microsphere and calculating the mean and variance.The SNR for a given exposure time was calculated from where I 0 denotes the mean maximum signal level above the background, σ 2 BG the variance of the background intensity values, and σ 2 I 0 the variance in the maximum intensities.Means and variances were calculated across the 80 images acquired at a given exposure time.The resulting SNRs ranged from 3.28 to 18.5.
The position of a particle in an image was determined using both the Gaussian fitting and the fluoroBancroft algorithms.This was done by selecting a 16×16 pixel array surrounding each particle and using the 256 data points in the two position estimation algorithms.Once the positions in each image were determined, the displacements between every pair of displaced images at a given exposure time were calculated by simply calculating the length of the displacement vector, leading to 1600 measurements for each microsphere.The values for σ x and σ y in the model ( 1) used in the fluoroBancroft algorithm were calcu- lated from (4) to be 203 nm for the given imaging parameters.
Gaussian fitting was performed using a least-squares minimization (using the built-in Matlab routine lsqnonlin) to fit a model of the form The position of the particle, (x o , y o ) as well as the values for σ x , σ y , and m in the model were varied.The algorithm was initialized using the coordinates of the center of the pixel with the largest intensity as the position of the particle, the maximum intensity value for m, and with σ x and σ y as in (4).All analysis was performed using the Matlab software package (MathWorks, Natick, MA, USA) on an Apple iMac (2.66 GHz Core 2 Duo) desktop computer.

Results and discussion
Typical results of the analysis for the displacement of one microsphere at a fixed SNR of 8.49 are shown in Fig. 3.The mean error in the measured displacement when the fluoroBancroft algorithm was used was 44.5 nm with a standard deviation of 24.5 nm.The performance was nearly as accurate as the Gaussian fitting algorithm which yielded a mean error of 37.0 nm and a standard deviation of 25.5 nm.Despite the similar performance, the fluoroBancroft algorithm executed on average over 200 times faster than the Gaussian fit approach (see also Fig. 7).

Biasedness of the fluoroBancroft estimator
The absolute position of the beads is unknown and thus one cannot determine the bias in the position estimates quantitatively since the position error cannot be calculated.Using the error in the displacements, however, we can investigate the bias qualitatively as follows.If the position estimates were Gaussian distributed with zero mean and variance σ 2 e , then the error in the translation measurement, d, would be Rayleigh distributed with the probability distribution function where U(d) is the unit step function (see, e.g.[20]).The mean and the variance for this distribution are given by Therefore, if the position estimates produced by the fluoroBancroft estimator are unbiased (zero-mean error) then, from (20), the mean and variance of the translation estimates would be related according to In Fig. 5.1 we show the difference between standard deviation of the error in the translation estimates produced by the fluoroBancroft algorithm and the mean of those errors, scaled according to (21).These results indicate that the algorithm exhibits some bias, particularly at low SNR.

Variance as a function of SNR
Although the error in the position estimates cannot be calculated, the accuracy of the algorithms in terms of their standard deviations in the estimates can be compared.In Fig. 5 we show the standard deviations in the position estimates in the x-direction (Fig. 5 (a)) and y-direction (Fig. 5 (b)).The two algorithms are very similar in the x-direction, although at low SNR the Gaussian fit fails to localize the particle.The mean difference in the standard deviation in the x-direction between the two algorithms across all SNR above 6 is only 0.086 nm.In the ydirection, however, the fluoroBancroft algorithm is somewhat worse than Gaussian fit with an average difference of 13.6 nm.In single particle tracking, the displacement measurements between frames of an image sequence reveal information about the motion of the particle.The errors in the position estimation propagate into the displacement measurement and ultimately determine its accuracy.The accuracy of the displacement measurement is shown in Fig. 6 (a) in which we show the standard deviation in the displacement error from the fluoroBancroft algorithm (solid blue line) and from the Gaussian fit algorithm (dotted red line).Note that the Gaussian fit algorithm was unable to consistently localize the particle at SNRs below 6 while the accuracy at these SNRs of the displacement error using the fluoroBancroft algorithm was on the order of 30 nm.The difference between the standard deviation of the Gaussian fit algorithm and the fluoroBancroft algorithm is shown in Fig. 6 (b).At SNRs below 10 the fluoroBancroft algorithm is consistently more accurate than the Gaussian fit with an average difference of 2.15 nm (omitting the results below an SNR of 6).At SNRs above 10 the Gaussian fit is more accurate with an average difference of -4.5 nm.

Computational speed
Exact execution times for the algorithms depend of course on the hardware platform and the software environment.The actual values of the execution times, therefore, have little meaning in and of themselves.The ratio of the two, however, is an indicator of the relative performance of the two algorithms.In Fig. 7  algorithm to the run time of the fluoroBancroft algorithm as a function of SNR.At SNRs above 6, the fluoroBancroft is on average 261 times faster than the numerical approach while at lower SNR the difference is even more dramatic.Note also that this algorithm is algebraic and thus its execution time is a function only of the amount of data (since that determines the size of the matrices in ( 16)).The variations in this ratio, therefore, arise entirely from the Gaussian fit approach.

Conclusions
This paper describes a novel algorithm for determining the position with nanometer-scale accuracy of a sub-diffraction limit fluorescent particle.The algorithm uses a range-based localization scheme to create an algebraic formula for the position of the particle.The performance of the algorithm was demonstrated experimentally and shown to have accuracy similar to the

Fig. 2 .
Fig. 2. Fluorescence image of a 20-nm diameter microsphere embedded with "Nile Red" fluorophores.Axis labels are pixel indices and grayscale represents fluorescence intensity.

Fig. 3 .
Fig.3.Results for images acquired with an SNR of 8.49.The mean of the error is indicated with a black dotted line and the standard deviation with a red dashed line.(a) Error in displacement when the fluoroBancroft algorithm was used for position estimation.The mean was 44.5 nm and the standard deviation was 24.5 nm.(b) Error in displacement when a Gaussian fit was used for position estimation.The mean of 37.0 nm was better than that obtained by fluoroBancroft while the standard deviation of 25.5 nm was slightly worse.

Fig. 4 .
Fig.4.Biasedness of the fluoroBancroft estimator.If the estimator were unbiased, then the mean and variance of the displacement estimates would be related by (21).That the difference, shown here, is not zero indicates the algorithm exhibits some bias in the position estimates.

Fig. 5 .
Fig. 5. Comparison of the standard deviation of the position estimate using fluoroBancroft (blue solid line) and Gaussian fit (red dashed line).(a) Standard deviation of the x−position estimate.(b) Standard deviation of the y−position estimate.In the x-direction, the two estimators are very similar while in the y− direction the Gaussian fit exhibits consistently better performance at SNRs above 6.

Fig. 6 .
Fig. 6.Comparison of the standard deviation of the errors in the calculated displacement.(a) Standard deviation of the error for each estimator.(b) Difference between Gaussian fit and fluoroBancroft results.