February 1, 2012 / Vol. 37, No. 3 / OPTICS LETTERS
413
Fast maximum likelihood algorithm for localization of fluorescent molecules Rebecca Starr1, Shane Stahlheber2,3, and Alex Small3,* 1
Kellogg Honors College and Department of Mathematics and Statistics, California State Polytechnic University, 3801 West Temple Avenue, Pomona, California 91768, USA 2
3
Department of Computer Science, California State Polytechnic University, 3801 West Temple Avenue, Pomona, California 91768, USA
Department of Physics and Astronomy, California State Polytechnic University, 3801 West Temple Avenue, Pomona, California 91768, USA *Corresponding author:
[email protected]
Received September 27, 2011; revised November 4, 2011; accepted December 17, 2011; posted December 20, 2011 (Doc. ID 155289); published February 1, 2012 A common task in microscopy is to fit an image of a fluorescent probe to a point spread function (PSF) in order to estimate the position of the probe. The PSF is often approximated as a Gaussian for mathematical simplicity. We show that the separable property of the Gaussian PSF enables a reduction of computational time from OL2 to OL, where L is the width (in pixels) of the image. When tested on realistic simulated data, our algorithm is able to localize the probes with precision close to the Cramér–Rao lower bound. © 2012 Optical Society of America OCIS codes: 180.2520, 100.6640.
Estimating the position of a fluorescent molecule from a microscope image is a key task in single-molecule biological experiments [1] and superresolution microscopy [2]. Although diffraction causes the width of the image to be of the order of the wavelength of light or larger, the coordinates of the image’s center (corresponding to the coordinates of the molecule) can be estimated with high precision, limited only by the number of photons detected [3]. In the ideal case, when pixelation effects are negligible and the only noise is shot noise, if the same molecule is imaged repeatedly (each acquisition lasting the same time and the molecule emitting photons at a steady rate) the maximum precision (i.e., smallest possible standard deviationpof the repeated position estimates) is δx ≥ λ∕2πNA N , where λ is the wavelength of the light emitted by the molecule, NA is the numerical aperture of the imaging system, and N is the number of photons detected from the molecule [3]. The theorem of Cramér and Rao states that the maximum precision is achievable if the position is estimated by fitting the system’s PSF to the image via maximum likelihood estimation (MLE) [3], with molecule coordinates (x0 ,y0 ) as fit parameters. Using MLE requires a model for the PSF (e.g., 2D Gaussian) and the noise in light detection (e.g., shot noise). Via the model we can calculate the likelihood of obtaining the data, and vary the parameters to maximize the likelihood of the data. MLE has been implemented on graphics processing units (GPUs) to estimate molecular positions very quickly [4], and the algorithm that we describe here is based on similar principles and thus easily extendable to GPU implementation. Other popular approaches, not requiring detailed noise models, include least-squares fits [5], a Gaussian mask algorithm based on least-squares [6], fluoroBancroft [7], and iterated centroid methods to mitigate background effects [8]. In each method the computational time is proportional to the number of pixels in the image. The proportionality between time and number of pixels depends on the complexity of the model and convergence criteria, but the number of floating point op0146-9592/12/030413-03$15.00/0
erations in an iteration of MLE or least-squares (or even the single iteration of fluoroBancroft) is proportional to the number of pixels. Although the PSF is (often) closer to an Airy function 2J 1 kr∕kr2 , where k 2πNA∕λ, fitting to a Gaussian is a common practice [4,6,9] that has been shown to give accurate and precise estimates with performance close to the theoretical limit [9]. Here we will perform faster MLE fits by using the separable property of the Gaussian 2 2 2 2 2 2 2 function: e−αk x−x0 y−y0 e−αk x−x0 e−αk y−y0 , where α is a numerical factor, assumed to be 0.28 [5]. Although the separable property of the Gaussian PSF has been used in deconvolution [10] and estimates of the Cramér– Rao lower bound [11], we are unaware of any published studies of the performance of maximum likelihood position estimators using our approach. We can model the signal S on a pixel centered at (x,y) as a Poisson random variable, and fit to it an expression μx; y equal to the sum of a PSF and a background: ZZ μx; y b I 0
2
0
2
0
2
e−αk x −x0 y −y0 dx0 dy0 ;
(1)
pixel area
where b is the average background count per pixel and I 0 is proportional to the integration time and the photon emission rate of the molecule. We assume that the background is a combination of out-of-focus and scattered fluorescence, and is hence a Poisson random variable with mean b. Fits are done on a square region of interest (ROI) of L × L pixels. If the PSF is negligible outside of the ROI, and if the “dead space” between pixels is a negligible fraction of the detector’s surface area, we can approximate the sum of all of the photon counts in the 2 0 2 column centered at x by integrating e−αk y −y0 over (−∞,∞). The integral over the pixel width is a difference of error functions, so the summed signal S s in a column centered at x has mean μs x given by © 2012 Optical Society of America
OPTICS LETTERS / Vol. 37, No. 3 / February 1, 2012
(2)
where a is the pixel width. Equation (2) depends only on the x coordinate of the column; we can thus estimate x0 and y0 by fitting Eq. (2) to L column sums and L row sums rather than L2 pixels. Consequently, we evaluate functions 2L times per iteration rather than L2 times. To estimate x0 via MLE, we must vary a vector of parameters θ x0 ; b; I 0 to maximize the logarithm of the likelihood LS s jθ, i.e., the conditional probability of obtaining the vector Ss of signals from each pixel, given the parameters θ. If we assume Poisson noise in photon detection, the log-likelihood is given by [3,4] log LSs jθ
X
S s x log μs x; θ − μs x; θ − log S s x!;
x
(3) where x is the pixel coordinate and the expected photon count is written μs x; θ to indicate that it depends on the fit parameters as well as the coordinate of the column. We can efficiently find a maximum by varying elements of θ via the Newton–Raphson method [12]. (For details, see the supplemental materials of [4].) Although MLE is not sensitive to small errors in the value of α [9], α may also be varied through a straightforward extension of our approach, if desired, e.g., to obtain depth information through the degree of defocus. To begin the first iteration, we will estimate b by setting the smallest column sum S s;min equal to bL, and estimate x0 with a background-corrected center of mass: xcm
X x
xS s x − bL∕
X
S s x − bL.
(4)
x
To estimate I 0 , we will set the largest column sum S s;max equal to the right side of Eq. (2), using our estimates of bL and x0 and the x coordinate of the appropriate column. The Newton–Raphson method is then used to iterate from this initial point in parameter space to the maximum of L. It is important to note that our initial x0 estimate should typically be very close to the maximum likelihood estimate: In many cases of practical interest, the photon count is high enough that the shot noise is approximately Gaussian, we are approximating the PSF as Gaussian, and we have a reasonable first-pass background correction. For Gaussian noise, a Gaussian PSF, and small pixels, the center of mass estimator is the MLE for the molecular position [3]. Also, when the initial estimate is close to the maximum, Newton–Raphson is known to require only a few iterations to converge [12]. We typically found that going beyond 2 − 3 iterations did not significantly reduce the variance of the position estimates, but six iterations were used for caution. We tested the performance of our algorithm by estimating x0 from simulated images generated with known molecular positions. The performance of MLE in fluorescent probe localization has been extensively studied as a function of photon count and background noise [9]; we will focus on cases where our assumptions (negligible
gaps between pixels, PSF magnitude negligible outside ROI) are of marginal validity. We assume shot noise in light detection and additive background fluorescence, and use an Airy PSF integrated (numerically) over the pixel area to generate photon counts. The wavelength of light is 550 nm, the NA is 1, the center-to-center pixel spacing is d 110 nm, and the number of pixels L along the width of the ROI is varied, as well as the pixel width a ≤ d. We generated the images so that 1000 photons are detected (in the absence of background) when the molecule is located at the center of the ROI in the best case for each figure; shrinking the ROI or displacing the molecule by (x0 ,y0 ) from the center of the ROI causes a smaller number of photons to be detected because some photons now fall outside the ROI. (However, we do not reduce the number of photons or the background if the pixel area is decreased, assuming that integration time increases commensurately, with proportional effects on signal and background.) We compared the backgroundfree case with a background of four photons/pixel (256 or 400 background photons total over the image, i.e., 25% or 40% of the actual signal). For each value of x0 in Fig. 1, we generated 1600 images to calculate the mean and standard deviation of the x0 estimates, thereby determining the bias (if any) and precision of the estimates. All calculations in Fig. 1 were done in Matlab, except that standard deviations in a “best case” were compared with the “practical localization accuracy measure” (PLAM) calculated with the program FandPLimitTool [9]. FandPLimitTool assumed that fitting was done on a 2D image rather than column sums, that there is no gap between pixels, that y0 0 (i.e., the ROI is correctly centered), and that the ROI is 12 × 12 (i.e., 2.4λ across, large enough to capture the entire central maximum of the PSF). In other calculations, we considered worse cases: a smaller ROI (8 × 8 or 10 × 10, corresponding to widths of 1.6λ or 2λ), an off-center ROI (y0 d or 2d), and substantial gaps between pixel edges (either 10% or 20% of the center-to-center spacing d). In both graphs, the PLAM for the ideal case is nearly constant: small displacements do not significantly reduce Localization Precision (nm)
p πI 0 L Erf αkx a∕2 − x0 2 2αk p − Erf αkx − a∕2 − x0 ;
μs x bL
(a)
6.5
Background: 0 photons/pixel
6 5.5 5
PLAM: a=d,12x12 ROI, y =0 0
4.5
a=d, 12x12 ROI, y =0 0
4 3.5
a=0.8d, 8x8 ROI, y =d 0
0
100 x (nm)
200
5.2
(b)
a=0.8d, 8x8 ROI, y =2d 0
a=0.8d, 10x10 ROI, y =d
0
Localization Precision (nm)
414
0
Background: 4 photons/pixel
a=0.8d, 10x10 ROI, y =2d 0
5
a=0.9d, 8x8 ROI, y =d 0
a=0.9d, 8x8 ROI, y0=2d
4.8
a=0.9d, 10x10 ROI, y =d 0
4.6
a=0.9d, 10x10 ROI, y =2d
4.4 4.2
0
100
200
x (nm) 0
Fig. 1. (Color online) Localization precision (standard deviation of position estimate) vsersus fluorophore coordinates (x0 ,y0 ) (relative to ROI center) for ideal and nonideal cases. a pixel edges width, d center-to-center pixel spacing.
February 1, 2012 / Vol. 37, No. 3 / OPTICS LETTERS
the number of photons in the ROI. The precision (standard deviation of estimates) is generally within 16% of the PLAM (comparable to or better than least-squares [9]), even if the molecule is displaced from the center by two pixels (i.e., a poorly drawn ROI). In all cases shown here the bias in the x0 estimate is smaller than the standard deviation. In a few cases the standard deviation is about 2% lower than the PLAM; we attribute this to the fact that when estimating the variance of a distribution, the estimate of the standard deviation will have a fracp tional uncertainty of order 1∕ n 2.5% for a sample of size n 1600. These discrepancies only occur when the sample standard deviation is close to the PLAM. When the sample standard deviation is (on average) substantially above the PLAM, the fluctuations in the plots never go below the PLAM. Furthermore, in the b 4 case the only noticeable changes seem to be fluctuations of order 4%, suggesting that small displacements do not affect localization precision. When we compare the background-free case [Fig. 1(a)] and the b 4 case [Fig. 1(b)], we see that in the absence of background the standard deviation starts at a lower value (as expected) and becomes worse (especially for a large ROI), while in the case with background the standard deviation starts at a higher value but rises more slowly. We attribute this to model mismatch: We are making the common approximation of fitting a Gaussian to data generated with an Airy PSF. This approximation is least accurate in the tails of the PSF, and a very far tail of the PSF will be found on the left side of the image for large x0 and a large ROI. Fitting this poorly matched model to our data will affect the log-likelihood score. This suggests that while our assumption of negligible PSF magnitude outside the ROI is less valid for a small ROI, the issue is mitigated by avoidance of other model mismatch effects. Also, the problem can be corrected in most situations by recentering the ROI if the estimated x0 is more than 1–2 pixels off-center, increasing computational time by at most a factor of 2 while substantially improving precision. Finally, we have created an ImageJ plugin that implements our localization algorithm (available athttp://code .google.com/p/molecule‑localization‑plugin/). In Fig. 2, we show the results of repeated application of this algorithm to fluorescent molecules with different subpixel positions. Red dots indicate actual molecular positions, centers of circles indicate mean position estimates, and radii are equal to 2× the standard deviation of the estimates. Background photon count is 196 per 7 × 7 ROI, and the number of photons emitted by the molecule are 350, 4 × 350 1400, and 4 × 1400 5600. The numbers of circles that fit into each pixel (32 , 72 , and 152 ) show that the standard deviation drops by a factor of ≈2 for each 4-fold increase, consistent with theory [2]. The uniformity of the radii, and the correspondence of
415
Fig. 2. (Color online) Positions of molecules in simulated images (red dots), mean estimated positions (centers of circles), and standard deviations of estimates (radius 2× standard deviation). 104 images were used for each circle.
circle centers with molecular positions, confirm that our algorithm is unbiased and robust against pixelation effects. In conclusion, we have shown that Gaussian fitting via MLE can be sped up from OL2 to OL by taking advantage of the separable property of the Gaussian PSF. The algorithm yields performance comparable to standard MLE, even for realistic cases in which key assumptions of the algorithm (no gaps between pixels, and a PSF completely contained within the ROI) are violated. This work was supported by a grant from the California State University Program for Education and Research in Biotechnology (CSUPERB). The manuscript was written at the Kavli Institute for Theoretical Physics, supported in part by National Science Foundation (NSF) grant PHY05-51164. We thank Raimund Ober and Anish Abraham for help with FandPLimitTool. References 1. W. Moerner, Proc. Natl. Acad. Sci. USA 104, 12596 (2007). 2. X. Zhuang, Nat. Photon. 3, 365 (2009). 3. R. J. Ober, S. Ram, and E. S. Ward, Biophys. J. 86, 1185 (2004). 4. C. S. Smith, N. Joseph, B. Rieger, and K. A. Lidke, Nature Methods 7, 373 (2010). 5. B. Zhang, J. Zerubia, and J. C. Olivo-Marin, Appl. Opt. 46, 1819 (2007). 6. R. E. Thompson, D. R. Larson, and W. W. Webb, Biophys. J. 82, 2775 (2002). 7. S. Andersson, Opt. Express 16, 18714 (2008). 8. A. J. Berglund, M. D. McMahon, J. J. McClelland, and J. A. Liddle, Opt. Express 16, 14064 (2008). 9. A. V. Abraham, S. Ram, J. Chao, E. S. Ward, and R. J. Ober, Opt. Express 17, 23352 (2009). 10. L. Lucy, Astron. J. 104, 1260 (1992). 11. S. Ram, E. Ward, and R. Ober, in Biomedical Engineering: Nano to Macro (IEEE, 2006), pp. 770. 12. W. H. Press, S. A. Teukolsky, W. T. Vetterling, and B. P. Flannery, Numerical Recipes: The Art of Scientific Computing, 3rd ed. (Cambridge University, 2007).
June 1, 2012 / Vol. 37, No. 11 / OPTICS LETTERS
1967
Fast maximum likelihood algorithm for localization of fluorescent molecules: erratum Rebecca Starr,1 Shane Stahlheber,2,3 and Alex Small3,* 1
Kellogg Honors College and Department of Mathematics and Statistics, California State Polytechnic University, 3801 West Temple Avenue, Pomona, California 91768, USA 2
3
Department of Computer Science, California State Polytechnic University, 3801 West Temple Avenue, Pomona, California 91768, USA
Department of Physics and Astronomy, California State Polytechnic University, 3801 West Temple Avenue, Pomona, California 91768, USA *Corresponding author:
[email protected] Received March 26, 2012; posted March 26, 2012 (Doc. ID 165202); published May 29, 2012
The code used to generate figures in our previous Letter [Opt. Lett. 37, 413 (2012)] contained an error. The error has been corrected and revised figures are presented here. The results are unchanged. © 2012 Optical Society of America OCIS codes: 180.2520, 100.6640.
Localization Precision (nm)
In our recently published Letter [1], the code used to perform the maximum likelihood estimation (MLE) contained an error. Our code did not properly vary all the fitting parameters (x0 , I 0 , and b) simultaneously. We have corrected the code to fix this error and have replicated the results of both figures with the corrected code. The results are shown in Fig. 1 and Fig. 2. The corrected code is available online [2]. Although a few cases have slightly higher or lower precision than in the original version, all cases remain clustered in the same range, i.e., within about 16% of the theoretical limit (labeled PLAM, for Practical Localization Accuracy Measure, calculated with code from the Ober Lab [3]).
6
(a) Background: 0 photons/pixel
5.5 5 4.5
PLAM: a=d, 12x12 ROI, y =0 0
a=d,12x12 ROI, y =0
4
0
a=0.8d, 8x8 ROI, y0=d
3.5
0
50
100
150
200
250
a=0.8d, 8x8 ROI, y =2d 0
a=0.8d, 10x10 ROI, y0=d
x (nm) 0
a=0.8d, 10x10 ROI, y =2d
Fig. 2. (Color online) Revised version of Fig. 2 from our original paper.
Although the corrected and original codes give essentially identical results in all cases studied in the original Letter, we further tested the code by reducing the signal to noise ratio, in order to identify cases in which the errors in the code might be problematic and verify that the corrected code is in fact superior. When we considered cases analogous to Fig. 2 but with a substantially lower photon count per molecule (88 photons), we verified that the original code gave poor reconstruction, while the corrected code gave spurious results that remained close to the theoretical limit.
Standard Deviation (nm)
0
5.4
(b) Background: 4 photons/pixel
a=0.9d, 8x8 ROI, y =d 0
5.2
a=0.9d, 8x8 ROI, y =2d 0
a=0.9d, 10x10 ROI, y =d
5
We thank S. Wolter for helping us identify the error in the code.
0
a=0.9d, 10x10 ROI, y =2d 0
4.8
References
4.6 4.4 4.2 0
50
100
150
200
250
x (nm) 0
Fig. 1. (Color online) Revised version of Fig. 1 from our original paper.
0146-9592/12/111967-01$15.00/0
1. R. A. Starr, S. Stahlheber, and A. Small, Opt. Lett. 37, 413 (2012). 2. http://code.google.com/p/molecule‑localization‑plugin/. 3. A. V. Abraham, S. Ram, J. Chao, E. S. Ward, and R. J. Ober, Opt. Express 17, 23352 (2009).
© 2012 Optical Society of America
