OSA's Digital Library

Virtual Journal for Biomedical Optics

Virtual Journal for Biomedical Optics

| EXPLORING THE INTERFACE OF LIGHT AND BIOMEDICINE

  • Editors: Andrew Dunn and Anthony Durkin
  • Vol. 9, Iss. 3 — Mar. 6, 2014
« Show journal navigation

A divide and conquer strategy for the maximum likelihood localization of low intensity objects

Alexander Krull, André Steinborn, Vaishnavi Ananthanarayanan, Damien Ramunno-Johnson, Uwe Petersohn, and Iva M. Tolić-Nørrelykke  »View Author Affiliations


Optics Express, Vol. 22, Issue 1, pp. 210-228 (2014)
http://dx.doi.org/10.1364/OE.22.000210


View Full Text Article

Acrobat PDF (1524 KB)





Browse Journals / Lookup Meetings

Browse by Journal and Year


   


Lookup Conference Papers

Close Browse Journals / Lookup Meetings

Article Tools

Share
Citations

Abstract

In cell biology and other fields the automatic accurate localization of sub-resolution objects in images is an important tool. The signal is often corrupted by multiple forms of noise, including excess noise resulting from the amplification by an electron multiplying charge-coupled device (EMCCD). Here we present our novel Nested Maximum Likelihood Algorithm (NMLA), which solves the problem of localizing multiple overlapping emitters in a setting affected by excess noise, by repeatedly solving the task of independent localization for single emitters in an excess noise-free system. NMLA dramatically improves scalability and robustness, when compared to a general purpose optimization technique. Our method was successfully applied for in vivo localization of fluorescent proteins.

© 2014 Optical Society of America

1. Introduction

With the advent of high time resolution microscopy and advances in the biological toolkit for microscopy, it has become possible to visualize the movement of cells and of single molecules inside a live cell [1

1. M. Coelho, N. Maghelli, and I. M. Tolić-Nørrelykke, “Single-molecule imaging in vivo: the dancing building blocks of the cell,” Integr. Biol. (Camb) 5, 748–758 (2013). [CrossRef]

]. For extracting information from these time-lapse images, it becomes essential to have the tools to track the motion of biological entities with high accuracy and efficiency. In cell biology, the objects of interest are often below the resolution limit of the microscope, appearing in the form of the system’s point spread function (PSF). Furthermore low-light conditions are employed so as to avoid photo damage and bleaching, making the task even more challenging (Figs. 1(a) and 1(b)).

Fig. 1 The NMLA tracks proteins in vivo (a) Image of a fission yeast cell expressing dynein heavy chain tagged with 3 GFPs (Dhc1-3GFP). The white line marks the cell outline, the green arrowhead marks a dynein molecule. The brightest spot is the spindle pole body. (b) The movement of dyneins is visualized in consecutive maximum intensity projections onto the y-axis (see color bar). (c) Schematic representation of the NMLA. The outer loop, which addresses the excess noise, is depicted in blue while the inner loop, which performs the localization, is depicted in gray. (d) Scheme of cell with the 2D trace of a dynein (green, corresponding to the green arrowhead in a) tracked using the NMLA. (e) Projection onto the y-axis of the traces obtained using the NMLA. The green trace is the path of the molecule marked in a.

Image generation in low-light microscopy is affected by several types of noise, which aggravate the localization task. We will focus on shot noise, readout noise, and excess noise, which are generally considered to be the most relevant with regard to localization problems [9

9. R. J. Ober, S. Ram, and E. S. Ward, “Localization accuracy in single-molecule microscopy,” Biophys. J. 86, 1185–1200 (2004). [CrossRef] [PubMed]

14

14. J. Chao, E. S. Ward, and R. J. Ober, “Fisher information matrix for branching processes with application to electron-multiplying charge-coupled devices,” Multidimens. Syst. Signal Process. 23, 349–379 (2012). [CrossRef] [PubMed]

]. Other kinds of noise include dark signal and clock induced charges [3

3. D. J. Denvir and E. Conroy, “Electron-multiplying CCD: the new ICCD, in Low-Light-Level and Real-Time Imaging Systems, Components, and Applications”, Proc. SPIE 4796, pp. 164–174 (2003). [CrossRef]

]. The latter two will be a topic in the Discussion section.

The first source of noise to be considered lies at the foundation of image generation: The number of photons detected in a pixel during an exposure is nondeterministic. The resulting effect, known as shot noise, is especially severe under low-light conditions. In the commonly used cameras with a charge-coupled device (CCD), the signal is additionally corrupted by additive Gaussian readout noise. It originates from the circuitry used to amplify the analog signal and finally convert it into a digital signal [2

2. K. Irie, A. E. McKinnon, K. Unsworth, and I. M. Woodhead, “A model for measurement of noise in CCD digital-video cameras,” Meas. Sci. Technol. 19, 045207 (2008). [CrossRef]

]. Readout noise poses a severe problem when little light is available. Cameras with electron multiplying CCDs (EMCCDs) reduce the effects of readout noise by amplifying the charge before it is read out. However, the nondeterministic amplification mechanism introduces yet another kind of noise, referred to as multiplicative or excess noise.

A variety of methods have been suggested to address the localization problem in noise affected microscopy images [4

4. B. C. Carter, G. T. Shubeita, and S. P. Gross, “Tracking single particles: a user-friendly quantitative evaluation,” Phys. Biol. 2, 60–72 (2005). [CrossRef] [PubMed]

]. One simple approach is the calculation of the center of mass of the signal [5

5. G. M. Lee, A. Ishihara, and K. A. Jacobson, “Direct observation of brownian motion of lipids in a membrane,” Proc. Natl. Acad. Sci. USA 88, 6274–6278 (1991). [CrossRef] [PubMed]

, 6

6. R. N. Ghosh and W. W. Webb, “Automated detection and tracking of individual and clustered cell surface low density lipoprotein receptor molecules.” Biophys. J. 66, 1301–1318 (1994). [CrossRef] [PubMed]

]. However this approach is not useful for images with multiple emitters and high background. Another widely used approach is based on least squares fitting of an approximated PSF to the image [7

7. C. M. Anderson, G. N. Georgiou, I. E. Morrison, G. V. Stevenson, and R. J. Cherry, “Tracking of cell surface receptors by fluorescence digital imaging microscopy using a charge-coupled device camera. low-density lipoprotein and influenza virus receptor mobility at 4 degrees c,” J. Cell Sci. 101 (Pt 2), 415–425 (1992). [PubMed]

,8

8. G. J. Schutz, H. Schindler, and T. Schmidt, “Single-molecule microscopy on model membranes reveals anomalous diffusion.” Biophys. J. 73, 1073–1080 (1997). [CrossRef] [PubMed]

], but this method is suboptimal. For low-light settings, it has been shown that the theoretically motivated maximum likelihood (ML) approach to localization has superior precision [9

9. R. J. Ober, S. Ram, and E. S. Ward, “Localization accuracy in single-molecule microscopy,” Biophys. J. 86, 1185–1200 (2004). [CrossRef] [PubMed]

].

The ML approach is based on a statistical model of the imaging process with its various sources of noise. ML localization methods address the problem of shot noise by viewing the number of photons in each pixel as random variable. In general, the number of photons ni in each pixel i is assumed to follow a Poisson distribution pni (ni; Ei). The mean of this distribution Ei is the flux at the pixel i, i.e. the expected number of photons detected at the pixel during an exposure. Ei depends on the locations of the objects which are to be localized. We will refer to these objects as emitters and to their manifestations in the image as blobs. By assuming that the number of photons ni at a pixel can be derived from the recorded intensity Ii at the pixel, the probability for the observation of an image is formulated as
i=1Mpni(ni;Ei),
(1)
where M is the number of pixels in the image. Different emitter locations change the expected photon counts Ei and thus result in different probabilities. The ML estimate is defined as the emitter locations which maximize Eq. (1).

Over the past years, this scheme and derivations have been used for localization in a number of studies, e.g. [9

9. R. J. Ober, S. Ram, and E. S. Ward, “Localization accuracy in single-molecule microscopy,” Biophys. J. 86, 1185–1200 (2004). [CrossRef] [PubMed]

12

12. Y. Wang, T. Quan, S. Zeng, and Z.-L. Huang, “PALMER: a method capable of parallel localization of multiple emitters for high-density localization microscopy,” Opt. Express 20, 16039–16049 (2012). [CrossRef] [PubMed]

]. In [9

9. R. J. Ober, S. Ram, and E. S. Ward, “Localization accuracy in single-molecule microscopy,” Biophys. J. 86, 1185–1200 (2004). [CrossRef] [PubMed]

] Ober et al. formulated the problem and calculated an absolute limit for accuracy in the presence of additional readout noise. In [10

10. C. S. Smith, N. Joseph, B. Rieger, and K. A. Lidke, “Fast, single-molecule localization that achieves theoretically minimum uncertainty.” Nat. Methods 7, 373–375 (2010). [CrossRef] [PubMed]

] the ML approach was implemented in a fast graphics processing unit (GPU) based method. Later, similar techniques were used to simultaneously localize multiple emitters [11

11. F. Huang, S. L. Schwartz, J. M. Byars, and K. A. Lidke, “Simultaneous multiple-emitter fitting for single molecule super-resolution imaging,” Biomed. Opt. Express 2, 1377–1393 (2011). [CrossRef] [PubMed]

, 12

12. Y. Wang, T. Quan, S. Zeng, and Z.-L. Huang, “PALMER: a method capable of parallel localization of multiple emitters for high-density localization microscopy,” Opt. Express 20, 16039–16049 (2012). [CrossRef] [PubMed]

].

During EMCCD amplification, photo-electrons are shifted through a multiplication register. According to Basden et al. [15

15. C. H. A. G. Basden, “Photon counting strategies with low light level CCDs,” Monthly Notices of the Royal Astron. Soc. 345, 985–991 (2003). [CrossRef]

] the number of output electrons Si given a certain number of input electrons ni approximately follows a discrete version of an Erlang distribution:
pSi|ni(Si|ni;g)=Sni1exp(Sig)gni(ni1)!,
(2)
where g is the mean amplification gain, which we will refer to as EMCCD gain. In settings with low photon counts, such as in biological applications, a high EMCCD gain is beneficial as it leads to better localization quality [16

16. J. Chao, E. S. Ward, and R. J. Ober, “Localization accuracy in single molecule microscopy using electron-multiplying charge-coupled device cameras, in Three-Dimensional and Multidimensional Microscopy: Image Acquisition and Processing XIX,” Proc. SPIE 8227,p. 82271P (2012). [CrossRef]

]. The effects of readout noise become insignificant with high EMCCD gains and we therefore omit them [3

3. D. J. Denvir and E. Conroy, “Electron-multiplying CCD: the new ICCD, in Low-Light-Level and Real-Time Imaging Systems, Components, and Applications”, Proc. SPIE 4796, pp. 164–174 (2003). [CrossRef]

, 10

10. C. S. Smith, N. Joseph, B. Rieger, and K. A. Lidke, “Fast, single-molecule localization that achieves theoretically minimum uncertainty.” Nat. Methods 7, 373–375 (2010). [CrossRef] [PubMed]

]. Furthermore, it will be assumed that the observed signal in a pixel is identical to the number of output electrons Si. We will use the term signal as opposed to intensity in the context of this model, where there is no deterministic mapping from observable pixel output to the number of photons in a pixel.

Here, we reformulate the optimization task to utilize this property in a novel nested expectation maximization (EM) algorithm that we refer to as Nested Maximum Likelihood Algorithm (NMLA), see Fig. 1(c). The EM meta algorithm was introduced by Dempster et al. [17

17. A. P. Dempster, N. M. Laird, and D. B. Rubin, “Maximum likelihood from incomplete data via the EM algorithm,” J. Royal Stat. Soc., B 39, 1–38 (1977).

] and is an established scheme for ML estimation in a setting with hidden variables. We show in simulation and experiments that our NMLA closely approaches the best theoretically possible precision, the Cramer-Rao lower bound (CRLB) [18

18. C. Rao, Linear Statistical Inference and its Applications (Wiley, 1975).

]. We compare the results and required computation time of the NMLA with a well established general purpose optimization method [19

19. M. J. D. Powell, “An efficient method for finding the minimum of a function of several variables without calculating derivatives,” Comput. J. 7, 155–162 (1964). [CrossRef]

] and find that NMLA yields dramatic improvements in scalability and robustness. In addition, the application of NMLA in a biological setting is demonstrated (Figs. 1(d) and 1(e)).

We will now describe the general outline of our method. Detailed theory and calculations will be discussed in the Method section below. The NMLA consists of an inner and an outer loop (Fig. 1(c)). Starting with the shot- and excess noise corrupted image, the outer loop addresses excess noise by calculating an expected photon count in each pixel. Based on these results, the inner loop estimates the location of the emitters. In a pure shot noise setting, the inner loop on its own is able to perform ML localization. The EM scheme of the inner loop implements a DC strategy by enabling the refinement of parameter estimates to be done independently for each emitter in each iteration. Our nested approach allows us to use this strategy even in a setting affected by excess noise.

We will first describe an alternate view of the pure shot noise model followed by the description of our algorithm’s inner loop and its application in a pure shot noise setting. Finally we will describe how the full NMLA implements an excess noise aware ML estimator.

2. Method

2.1. An alternate view of the shot noise model

Our shot noise model is equivalent to the one formulated in Eq. (1). In order to motivate the inner loop of the NMLA, we suggest an alternative view of the image generation process, focusing on the distribution of photons in the image instead of the number of photons in each pixel. A similar view has been described for image simulation in [9

9. R. J. Ober, S. Ram, and E. S. Ward, “Localization accuracy in single-molecule microscopy,” Biophys. J. 86, 1185–1200 (2004). [CrossRef] [PubMed]

, 20

20. J. Chao, S. Ram, E. S. Ward, and R. J. Ober, “Ultrahigh accuracy imaging modality for super-localization microscopy,” Nat. Methods 10, 335–338 (2013). [CrossRef] [PubMed]

], but was not utilized for localization. A comparable view has been applied to spectra analysis [21

21. A. Steinborn, S. Taut, V. Brendler, G. Geipel, and B. Flach, “TRLFS: analysing spectra with an expectation-maximization (EM) algorithm.” Spectrochim Acta A Mol Biomol Spectrosc. 71, 1425–1432 (2008). [CrossRef] [PubMed]

].

Image generation is seen as a two step process: First, a number N is randomly determined from a Poisson distribution pN(N; E) which is parametrized by the total flux of the image E=i=1MEi. Next, N photons are randomly distributed on the pixel grid of the image I⃗ = (I1,..., IM) according to a probability distribution pi(i; θ⃗). Each photon increases the intensity Ii at the pixel i that it hits by the constant c. A pixel’s intensity is thus Ii = nic. The vector θ⃗ = (θ⃗k, pk, p0|k = 1,...,m) summarizes the parameters of m emitters and additional stray background light. While the vector θ⃗k holds the location of emitter k, the term pk stands for the percentage of light from this emitter. It can also be seen as the a priori probability of a photon to have the origin k. The term p0 denotes the percentage of background light. The numbers pk are under the constraint that their sum is equal to one.

Our perspective of image generation allows us to write the probability for the generation of an image I⃗ as the product:
pI(I;θ,E)=pN(N;E)pI|N(I|N;θ).
(3)
Please refer to the Appendix for a detailed deduction. The term pI⃗|N(I⃗|N; θ⃗) stands for the conditional probability of the image given that it contains a total of N photons. It is a multi-nomial distribution and can be written as pI|N(I|N;θ)=c^i=1Mpi(i;θ)Iic. The constant ĉ is the accordant multinomial coefficient.

Let us now examine the probability distribution pi(i; θ⃗) according to which the photons are distributed on the pixel grid. We consider it to be a mixture of m + 1 components:
pi(i;θ)=k=0mpi|k(i|k;θk)pk.
(4)

The term pi|k(i|k; θ⃗k) denotes the probability for a photon to be recorded at a specific pixel i, given it is of origin k. For photons originating from the background (k = 0) it is considered a uniform distribution. The distribution for photons coming from one of the emitters (k > 0) is derived by integrating the PSF over the pixel’s area and normalizing it with the PSF’s integral over the image area. We use a Gaussian approximation of the PSF for this purpose, which was evaluated to be suitable for most applications in [13

13. K. I. Mortensen, L. S. Churchman, J. A. Spudich, and H. Flyvbjerg, “Optimized localization analysis for single-molecule tracking and super-resolution microscopy,” Nat. Methods 7, 377–381 (2010). [CrossRef] [PubMed]

] and is a conventional choice used among others in [10

10. C. S. Smith, N. Joseph, B. Rieger, and K. A. Lidke, “Fast, single-molecule localization that achieves theoretically minimum uncertainty.” Nat. Methods 7, 373–375 (2010). [CrossRef] [PubMed]

12

12. Y. Wang, T. Quan, S. Zeng, and Z.-L. Huang, “PALMER: a method capable of parallel localization of multiple emitters for high-density localization microscopy,” Opt. Express 20, 16039–16049 (2012). [CrossRef] [PubMed]

].

2.2. Localization with the inner loop

It has been shown that the estimates obtained by the EM scheme improve with each iteration until convergence [22

22. M. I. Schlesinger and V. Hlavac, Ten Lectures on Statistical and Structural Pattern Recognition, (Kluwer Academic Publishers, 2002) vol. 24 of Computational Imaging and Vision. [CrossRef]

]. We will now describe the two steps in detail.

  • The inner E-step. Using Bayes’ theorem, it is possible to calculate for each pixel i, the probability pk|i(k|i; θ⃗t) that a photon measured at this position has a specific origin k:
    pk|i(k|i;θt)=pi|k(i|k;θt)pkk=0mpi|k(i|k;θt)pk
    (7)
    The E-step uses these probabilities to calculate a set of numbers Ii,kt.
    Ii,kt=Iipk|i(k|i;θt)
    (8)
    This can be seen as separating the intensity Ii at each pixel i into fractions Ii,kt, assumed to originate from the different sources k.
  • The inner M-step. The estimates for the a priori probabilities are acquired in the following way:
    pkt+1=i=1MIi,kti=1MIi
    (9)
    The a priori probability of the background can be calculated as p0t+1=1k=1mpkt+1. Based on the results of the E-step the new parameters for each emitter θkt+1 of the new model are independently estimated as:
    θkt+1=argmaxθki=1M(Ii,ktlnpi|k(i|k;θk))
    (10)
    In order to maximize the expression we use the general purpose optimization method known as Powell’s method [19

    19. M. J. D. Powell, “An efficient method for finding the minimum of a function of several variables without calculating derivatives,” Comput. J. 7, 155–162 (1964). [CrossRef]

    ]. Optimization can be performed independently for each emitter.

Estimating the flux. For some applications, it may be desired to estimate not only the above mentioned parameters of each light source but also the amount of light that is emitted by it. The average number of photons Êk recorded per exposure, originating from a source k, will be called the flux of the source. It can be calculated from the total flux of the image E:
E^k=Epk
(11)

While the inner loop described above provides an estimate for the a priori probability pk of each light source, it is straight-forward to calculate an ML estimate for the total image flux E. It can be done independently of the rest of the parameters θ⃗. The result is proportional to the sum of all pixel intensities:
E*=i=1MIic
(12)

Accordingly, the flux of a specific light source can be estimated up to the factor c as E^k*=pkE*. The gain factor c of the system is required in order to compute absolute estimates of the photon counts. If it is unknown, the results obtained still allow for comparison between different light sources.

2.3. Localization with the complete NMLA

Let us now examine the complete model including the probabilistic EMCCD amplification as described in [15

15. C. H. A. G. Basden, “Photon counting strategies with low light level CCDs,” Monthly Notices of the Royal Astron. Soc. 345, 985–991 (2003). [CrossRef]

]. We observe an image of EMCCD signals S⃗ = (S1,..., SM) generated from an unknown photon image n⃗ = (n1,..., nM). During amplification, each pixel’s signal Si is taken from the distribution described in Eq. (2). The new task is now to concurrently find the model parameters θ⃗* and the flux of the image E*, which maximize the probability to observe the image S⃗ of signals. In this model, it is not possible to find the estimates independently:
(θ*,E*)=argmax(θ,E)pS(S;θ,E;g)
(13)
=argmax(θ,E)nnsetpS,n(S,n;θ,E;g),
(14)
where the function pS⃗,n⃗ (S⃗,n⃗; θ⃗, E;g) denotes the joint probability to acquire a photon image n⃗ and observe S⃗. The term n⃗set denotes the set of possible photon images of the same size.

Optimizing the estimate for the total image flux: The maximum likelihood estimation of the total image flux can be identified with the task from Eq. (6). It has an analytical solution analogous to Eq. (12):
E*=i=1Mn^it
(19)
As in the pure shot noise model, the flux for the different light sources E^k* can be calculated accordingly as E^k*=pk*E*.

3. Experiments

3.1. Evaluation with simulated images

To evaluate our method and determine what additional gain is to be expected from the extension of the pure shot noise model, the inner loop as well as the complete NMLA were applied to simulated images with and without excess noise. We compared the standard deviation of the localization results with theoretically predicted bounds for precision (Fig. 2). As expected, the inner loop of our algorithm attains the CRLB in a pure shot noise case for a range of different values of flux (Fig. 2(a)) [10

10. C. S. Smith, N. Joseph, B. Rieger, and K. A. Lidke, “Fast, single-molecule localization that achieves theoretically minimum uncertainty.” Nat. Methods 7, 373–375 (2010). [CrossRef] [PubMed]

]. After deriving the CRLB for the excess noise aware ML estimator, we find that it is attained by the NMLA. We confirm the finding from [20

20. J. Chao, S. Ram, E. S. Ward, and R. J. Ober, “Ultrahigh accuracy imaging modality for super-localization microscopy,” Nat. Methods 10, 335–338 (2013). [CrossRef] [PubMed]

] that the approximation of the CRLB for excess noise, obtained by multiplying the CRLB for a pure shot noise setting by a factor of 2 is not appropriate for low photon counts. Surprisingly, this approximation from [13

13. K. I. Mortensen, L. S. Churchman, J. A. Spudich, and H. Flyvbjerg, “Optimized localization analysis for single-molecule tracking and super-resolution microscopy,” Nat. Methods 7, 377–381 (2010). [CrossRef] [PubMed]

] closely corresponds to the results attained by the ML estimator, which is based on a pure shot noise model and implemented by our inner loop (Figs. 2,3).

Fig. 2 The NMLA and its inner loop closely approach the theoretical bounds in simulation. (a) The standard deviation of the localization results in one dimension as a function of the flux of the blob. The complete NMLA (red squares) and the inner loop (red triangles) were applied to the simulated images. The inner loop was also applied to the same image series without excess noise (blue triangles). The results are compared with theoretical bounds: the Cramer-Rao lower bound (CRLB) for a scenario without excess noise (blue line), the CRLB considering excess noise (red line) and the approximated bound suggested in [13] (black line). The NMLA yields the highest possible precision, given by the CRLB, in settings with excess noise. (b) The standard deviation of the localization results as a function of background flux. The flux of the blob was set to 25 photons. The remaining parameters and the algorithms employed were as in a. Note that the deviation between the theoretical bounds and the results of the estimators (also seen in [23]) is a general feature of the ML localization. In a and b, simulated images (without excess noise outlined in blue, and the same image with excess noise outlined in red) correspond to different values of flux/background flux indicated by the dashed lines under the images. In the images without excess noise the colors indicate the number of photons in each pixel. In the images with excess noise colors indicate the estimated number of photons as described in the Appendix.
Fig. 3 The standard deviation of the localization results as a function of the measured flux of fluorescent bead images. Three 3000-frame-long movies were obtained using Total Internal Reflection Fluorescence (TIRF) microscopy with different laser powers and EMCCD gain of 300. The pixel size was 106 nm (see Appendix for details). The beads were localized with the NMLA (squares) and its inner loop (triangles). Each pair of a triangle and a square at a certain value of flux corresponds to a single bead; displayed are images of three tracked beads connected to their resulting data point. The results are compared with the CRLB considering excess noise. The background flux used in the calculation of the bound was derived from the average intensities in manually selected empty areas of the movies (see Appendix for details). The inset is an enlarged view of the section corresponding to the low flux region (dashed lines) including the theoretical bound from [13].

The same general picture can be observed when different levels of Poisson background noise are applied (Fig. 2(b)). With increasing background, however, a deviation between the results of all three ML estimators and the corresponding theoretical bounds is evident. This apparent general feature of the ML localization has been observed independently in [23

23. S. Stallinga and B. Rieger, “The effect of background on localization uncertainty in single emitter imaging,” Proceedings of 9th IEEE International Symposium on Biomedical Imaging 2012, 988–991(2012).

]. Since heavy background is a common feature in biological applications this effect might be worth further examination.

3.2. Evaluation on images of fluorescent beads

To evaluate our technique experimentally, we imaged fluorescent beads under various light conditions and tracked them using the NMLA and its inner loop (Fig. 3). As in simulation, the results are in agreement with the theoretical bound. The NMLA yielded a lower standard deviation than its inner loop for each of the observed tracks. As in simulation the difference is larger for low photon counts.

3.3. Evaluation of computation time and robustness

We compared the performance of the NMLA and a direct application of Powell’s method [19

19. M. J. D. Powell, “An efficient method for finding the minimum of a function of several variables without calculating derivatives,” Comput. J. 7, 155–162 (1964). [CrossRef]

] to the log likelihood as described in [13

13. K. I. Mortensen, L. S. Churchman, J. A. Spudich, and H. Flyvbjerg, “Optimized localization analysis for single-molecule tracking and super-resolution microscopy,” Nat. Methods 7, 377–381 (2010). [CrossRef] [PubMed]

] (supplementary equation 138). To assess the scalability we applied both optimizers to a grid of overlapping emitters and measured the required computation time for different emitter numbers (Fig. 4(a)). To compare the robustness of the methods we initialized them at varying positions and measured localization success in a single emitter setting (Fig. 4(b)). In both cases our NMLA outperforms Powell’s method.

Fig. 4 Th NMLA outperforms Powell’s method with regard to scalability and robustness. (a) Required computation time as a function of the number of emitters. Different numbers of overlapping emitters were tracked with three different implementations of the excess noise aware ML estimator: The computation time required by Powell’s method (circles) grows significantly faster than the time required by the NMLA (squares). A naively parallelized version of the NMLA (diamonds) yields a further improvement. (b) Number of inliers in localization results as a function of the standard deviation of the random initialization. In each frame the algorithms were initialized at a Gaussian distributed random location around the true position of the emitter. Every location estimate closer than 1.3 pixels to the correct location was considered an inlier. The NMLA achieves 100% inliers for higher standard deviations and does not loose as many as Powell’s method.

Note that the set of parameters Powell’s method is employed on differs from our θ⃗ since it does not include the a priori probabilities pk. Instead the flux of each emitter Êk = pkE, ∀k = 1,...,m as well as the background flux Ê0 = p0E are used. When performing our experiments, we used the implementation of Powell’s method in the Apache Commons Mathematics Library. We set the absolute and relative threshold used as stopping criteria to 0.001, which produced results of similar quality as the NMLA.

4. Discussion

4.1. General findings regarding the excess noise aware ML estimator

The ML estimator based on the pure shot noise model is frequently used on EMCCD images, e.g. in [10

10. C. S. Smith, N. Joseph, B. Rieger, and K. A. Lidke, “Fast, single-molecule localization that achieves theoretically minimum uncertainty.” Nat. Methods 7, 373–375 (2010). [CrossRef] [PubMed]

12

12. Y. Wang, T. Quan, S. Zeng, and Z.-L. Huang, “PALMER: a method capable of parallel localization of multiple emitters for high-density localization microscopy,” Opt. Express 20, 16039–16049 (2012). [CrossRef] [PubMed]

], likely because the theoretical understanding of the problem has only recently been achieved [13

13. K. I. Mortensen, L. S. Churchman, J. A. Spudich, and H. Flyvbjerg, “Optimized localization analysis for single-molecule tracking and super-resolution microscopy,” Nat. Methods 7, 377–381 (2010). [CrossRef] [PubMed]

,14

14. J. Chao, E. S. Ward, and R. J. Ober, “Fisher information matrix for branching processes with application to electron-multiplying charge-coupled devices,” Multidimens. Syst. Signal Process. 23, 349–379 (2012). [CrossRef] [PubMed]

]. We suggest to use the ML estimator based on the pure shot noise model as an approximate solution only if enough light is available. One should also consider the findings from [14

14. J. Chao, E. S. Ward, and R. J. Ober, “Fisher information matrix for branching processes with application to electron-multiplying charge-coupled devices,” Multidimens. Syst. Signal Process. 23, 349–379 (2012). [CrossRef] [PubMed]

] and determine whether EMCCD amplification is the tool of choice in such a setting. On the other hand, for the low photon count setting which is relevant in biology and astronomy, we conclude from theory, simulation, and experiments that significantly better results can be obtained when the statistical model accounts for excess noise in addition to shot noise. The NMLA is a powerful implementation of the excess noise aware ML estimator.

4.2. Practical advantages of NMLA

We see three major advantages of our method with regard to practical application:

4.3. Conceptual advances

By wrapping the inner loop in another EM algorithm we are able to effectively adapt these properties into a setting involving excess noise. To our knowledge, we are the first to explore the excess noise aware ML estimator for multiple emitters and from an algorithmic point of view.

4.4. Perspective

Our approach allows us to reduce the problem of ML localization in an excess noise affected setting to repeated localization in a pure shot noise setting. Therefore, it is possible to apply any new algorithmic finding for the pure shot noise setting in the presence of excess noise. An example could be the algorithm in [25

25. R. Starr, S. Stahlheber, and A. Small, “Fast maximum likelihood algorithm for localization of fluorescent molecules,” Opt. Lett. 37, 413–415 (2012). [CrossRef] [PubMed]

], which radically reduces the complexity for the localization of a single emitter in a shot noise setting. Our framework makes it possible to incorporate such methods in the inner M-step, thereby making them applicable with excess noise.

Even though only two types of noise, shot noise and excess noise were considered, our general scheme is not limited to this noise model. Other components could be included without structural change if they prove to be beneficial. Readout noise has little effect in our setting [3

3. D. J. Denvir and E. Conroy, “Electron-multiplying CCD: the new ICCD, in Low-Light-Level and Real-Time Imaging Systems, Components, and Applications”, Proc. SPIE 4796, pp. 164–174 (2003). [CrossRef]

], but it could be included by calculating pSi|ni (Si|ni; g) as convolution of Eq. (2) and a Gaussian distribution [13

13. K. I. Mortensen, L. S. Churchman, J. A. Spudich, and H. Flyvbjerg, “Optimized localization analysis for single-molecule tracking and super-resolution microscopy,” Nat. Methods 7, 377–381 (2010). [CrossRef] [PubMed]

]. The other parts of our method would remain unaffected. A different source of noise called dark signal is a Poisson distributed charge arising from imperfections in the production process of the chip [2

2. K. Irie, A. E. McKinnon, K. Unsworth, and I. M. Woodhead, “A model for measurement of noise in CCD digital-video cameras,” Meas. Sci. Technol. 19, 045207 (2008). [CrossRef]

]. Dark signal is pixel dependent [2

2. K. Irie, A. E. McKinnon, K. Unsworth, and I. M. Woodhead, “A model for measurement of noise in CCD digital-video cameras,” Meas. Sci. Technol. 19, 045207 (2008). [CrossRef]

] and its intensity pattern on the sensor can be measured. Consequently, one could extend our noise model by adding a measured pixel dependent term to each pixel’s flux Ei. All other calculations could be performed as described above. Finally, our algorithm could include clock induced charges, which are generated by shifting electrons through the CCD [26

26. M. Hirsch, R. J. Wareham, M. L. Martin-Fernandez, M. P. Hobson, and D. J. Rolfe, “A stochastic model for electron multiplication charge-coupled devices – from theory to practice,” PLoS ONE 8, e53671 (2013). [CrossRef]

]. According to [26

26. M. Hirsch, R. J. Wareham, M. L. Martin-Fernandez, M. P. Hobson, and D. J. Rolfe, “A stochastic model for electron multiplication charge-coupled devices – from theory to practice,” PLoS ONE 8, e53671 (2013). [CrossRef]

] these charges should be considered in the same way as dark signal. It would thus be possible to incorporate clock induced charges into our algorithm in the same manner as dark signal.

It should be noted that since we were aiming at a qualitative demonstration, the NMLA currently relies on user initialization and our implementation has not yet been optimized with regard to performance. While an interactive rather than an automatic behavior is desirable for some applications we see no reason why our general algorithmic scheme could not be combined with initialization mechanisms and highly optimized GPU methods as used in [10

10. C. S. Smith, N. Joseph, B. Rieger, and K. A. Lidke, “Fast, single-molecule localization that achieves theoretically minimum uncertainty.” Nat. Methods 7, 373–375 (2010). [CrossRef] [PubMed]

12

12. Y. Wang, T. Quan, S. Zeng, and Z.-L. Huang, “PALMER: a method capable of parallel localization of multiple emitters for high-density localization microscopy,” Opt. Express 20, 16039–16049 (2012). [CrossRef] [PubMed]

] to produce a fully automatic system.

Implementations of the NMLA are invaluable for studying biological phenomena such as membrane diffusion and cellular kinetics of molecules. While conventional methods of studying particle properties inside the cell such as FRAP (Fluorescent Recovery After Photobleaching) or FCS (Fluorescent Correlation Spectroscopy) provide bulk measurements of a population of molecules, single particle tracking enabled by the NMLA implementation offers kinetic data at the single-molecule level. We can thus study the temporal evolution of each particle, as well as detect the existence of different populations of particles. In addition, our NMLA offers the possibility of simultaneous tracking of multiple molecules imaged with low-light in a crowded environment such as inside the cell, while at the same time ensuring a lower computation time but higher precision than other methods.

5. Conclusion

ML estimation in EMCCD images and its potential in low-light microscopy has accrued considerable theoretical understanding. [14

14. J. Chao, E. S. Ward, and R. J. Ober, “Fisher information matrix for branching processes with application to electron-multiplying charge-coupled devices,” Multidimens. Syst. Signal Process. 23, 349–379 (2012). [CrossRef] [PubMed]

]. We have designed and tested a new tool for ML estimation that can help to utilize these findings in cell biology. A variant of this tool has already been successfully employed by us in [27

27. I. Kalinina, A. Nandi, P. Delivani, M. R. Chacón, A. H. Klemm, D. Ramunno-Johnson, A. Krull, B. Lindner, N. Pavin, and I. M. Tolić-Nørrelykke, “Pivoting of microtubules around the spindle pole accelerates kinetochore capture,” Nat. Cell Biol. 15, 82–87 (2013). [CrossRef]

, 28

28. V. Ananthanarayanan, M. Schattat, S. K. Vogel, A. Krull, N. Pavin, and I. M. Tolić-Nørrelykke, “Dynein motion switches from diffusive to directed upon cortical anchoring,” Cell 153, 1526–1536 (2013). [CrossRef] [PubMed]

].

An implementation of the NMLA is available as Open Source software under the GNU General Public License. It can be installed under the name Low Light Tracking Tool (http://fiji.sc/Low_Light_Tracking_Tool) via the image processing software Fiji [29

29. J. Schindelin, I. Arganda-Carreras, E. Frise, V. Kaynig, M. Longair, T. Pietzsch, S. Preibisch, C. Rueden, S. Saalfeld, B. Schmid, J.-Y. Tinevez, D. J. White, V. Hartenstein, K. Eliceiri, P. Tomancak, and A. Cardona, “Fiji: an open-source platform for biological-image analysis,” Nat. Methods 9, 676–682 (2012). [CrossRef] [PubMed]

].

Appendix Imaging using the Total Internal Reflection Fluorescence Microscope

An inverted stand, manual XY stage Olympus IX71 microscope (Olympus, Tokyo, Japan) with custom-built TIRF condenser and manual TIRF angle adjustment was employed for imaging beads (Fig. 3) and dyniens inside the cytoplasm of a fission yeast zygote (Figs. 1(a) and 1(b)). Imaging was performed using an Olympus UApo 150× 1.45 Oil TIRFM inf/0.13–0.21 corr (Olympus, Tokyo, Japan) objective with diode-pumped solid state 491nm laser for GFP excitation and 560 nm laser for mCherry/tdTomato excitation (75mW; Cobolt, Solna, Sweden). Laser intensity was controlled using the acousto-optic tunable filter in the Andor Revolution laser combiner (ALC, Andor Technology plc., Belfast, UK). The wavelength filter used was the BL HC 525/30 (Semrock Inc., Rochester, NY, USA). The microscope was equipped with an Andor iXon EM+ DU-897 BV back illuminated EMCCD (Andor Technology plc., Belfast, UK) with pixel size of EMCCD chip being 16μm and image pixel size being 0.106μm with the 150× objective. The system was controlled using the Andor iQ software version 1.9.1 (Andor Technology plc., Belfast, UK). For imaging dynein in the cytoplasm (Fig. 1), the imaging conditions used were: excitation with 80% power (18mW) of 491nm laser, exposure time of 5–9ms, with 2000 continuous repetitions. The dyneins in the strain used in this paper are labeled with triple GFP [30

30. S. K. Vogel, N. Pavin, N. Maghelli, F. Jülicher, and I. M. Tolić-Nørrelykke, “Self-organization of dynein motors generates meiotic nuclear oscillations,” PLoS Biol. 7, 0918–0928 (2009). [CrossRef]

]. The beads (Estapor 0.103 μm fluorescent beads, Merck KGaA, Darmstadt, Germany) were imaged with 2%, 3% and 5% power (0.45mW, 0.675mW and 1.125mW respectively) of 491nm laser, exposure time of 5ms, with 3000 continuous repetitions. The EMCCD gain used was 300 and the pre-amplification gain was 4.8.

The image probability as product

In Eq. (3) we formulate the probability for the observation of a specific image I⃗ created by a total of N photons as the product of the probability pN(N; E) that an image contains exactly N photons and the probability pI⃗|N(I⃗|N; θ⃗). that the image is observed under the condition, that it contains exactly N photons. We will now explain why this is possible. According to Bayesian theory we can write:
pI(I;θ,E)=N^=0(pN(N^;E)pI|N(I|N;θ)).
(20)
Since it is impossible to generate the image I⃗ from any number of photons other than the exact number of photons N that created it, the probability pI⃗|N(I⃗|N̂; θ⃗) will be 0 in all summands except for the single case = N. We can thus write:
pI(I;θ,E)=pN(N;E)pI|N(I|N;θ).
(21)

Calculating the CRLB

We calculated the CRLB as the inverse of the Fisher information matrix [18

18. C. Rao, Linear Statistical Inference and its Applications (Wiley, 1975).

]. For the pure shot noise setting the elements of the Fisher information matrix were calculated as described in [10

10. C. S. Smith, N. Joseph, B. Rieger, and K. A. Lidke, “Fast, single-molecule localization that achieves theoretically minimum uncertainty.” Nat. Methods 7, 373–375 (2010). [CrossRef] [PubMed]

]. For the case with excess noise, the elements of the Fisher information matrix can be computed as
FIMab=i=1MESi;θ,E[lnq(Si;Ei)alnq(Si;Ei)b]
(22)
=i=1MSiq(Si;Ei)lnq(Si;Ei)alnq(Si;Ei)b,
(23)
where a, bθ⃗ ∪ {E} are the parameters and E′Si;θ⃗, E stands for the expected value with regard to the signal Si of the pixel i and q(Si; Ei) describes the probability distribution of the signal Si, which is parametrized by the pixel’s flux Ei = pi (i; θ⃗) · E. A deduction of Eq. (22) can be found in [14

14. J. Chao, E. S. Ward, and R. J. Ober, “Fisher information matrix for branching processes with application to electron-multiplying charge-coupled devices,” Multidimens. Syst. Signal Process. 23, 349–379 (2012). [CrossRef] [PubMed]

]. As done in [13

13. K. I. Mortensen, L. S. Churchman, J. A. Spudich, and H. Flyvbjerg, “Optimized localization analysis for single-molecule tracking and super-resolution microscopy,” Nat. Methods 7, 377–381 (2010). [CrossRef] [PubMed]

] we calculate q(Si; Ei) efficiently as
q(Si;Ei)=n=0pSi|ni(Si|n;g)pni(n;Ei)
(24)
=exp(Ei)δ(Si)+exp(Si/gEi)g1EiEiSi/gI1(2EiSi/g),
(25)
with I1 being the modified Bessel function of the first kind of order one.

Deduction for the use of the Bessel function

In this section the step from the inconvenient Eq. (15) to the convenient Eq. (17) is derived. Eq. (15) is inconvenient with regard to computation, because of the infinite number of summands. Eq. (17) is convenient because the Bessel functions are very well investigated and fast implementations are available.

Reiterating, Eq. (2) describes the probability that a signal Si is observed in a pixel given a photon count of ni. Since the factorial of a negative number as well as the division by 0 is not defined we assume that this probability pSi|ni (Si|ni = 0; g) = 0 for all Si > 0 and that pSi|ni (Si = 0|ni = 0; g) = 1. We furthermore assume that pSi|ni (Si = 0|ni; g) = 0 for all ni > 0. These are reasonable assumptions when considering the underlying model of electron amplification. At least one electron has to be present to generate an output of Si > 0. Let us furthermore remember the assumption that the number of photons in each pixel is Poisson distributed:
pn(n;Ei)=Einn!exp(Ei).
(26)
We will now examine the denominator in Eq. (15):
n=0pSi|ni(Si|n;g)pn(n;Ei)=exp(Ei)δ(Si)+n=1pSi|ni(Si|n;g)pn(n;Ei)=exp(Ei)δ(Si)+n=1Si(n1)exp(Si/g)gn(n1)!Einn!exp(Ei)=exp(Ei)δ(Si)+exp(Si/gEi)n=1Si(n1)Eingnn!(n1)!
(27)
substituting n = l + 1 we have
exp(Ei)δ(Si)+exp(Si/gEi)l=0SilEiEilggll!(l+1)!=exp(Ei)δ(Si)+exp(Si/gEi)g1Eil=01l!(l+1)!(SiEi/g)l=exp(Ei)δ(Si)+exp(Si/gEi)g1Eil=01l!(l+1)!(SiEi/g)lEiSi/gEiSi/g=exp(Ei)δ(Si)+exp(Si/gEi)g1EiEiSi/gl=01l!(l+1)!(SiEi/g)lEiSi/g=exp(Ei)δ(Si)+exp(Si/gEi)g1EiEiSi/gl=01l!(l+1)!(2EiSi/g2)2l+1I1(2EiSi/g)=exp(Ei)δ(Si)+exp(Si/gEi)g1.EiEiSi/gI1(2EiSi/g)
(28)
where
I1(2EiSi/g)=l=01l!(l+1)!(2EiSi/g2)2l+1
(29)
is the modified Bessel function of the first kind of order one:
I1(z)=l=01l!(l+1)!(z2)2l+1
(30)
Let us now consider the numerator from Eq. (15):
n=1pSi|ni(Si|n;g)pn(n;Ei)n=n=1Si(n1)gn(n1)!exp(Si/g)exp(Ei)Einn!n=exp(Si/gEi)n=1Si(n1)Eingn((n1)!)2
(31)
substituting again n = l + 1 we get
exp(Si/gEi)l=0SilEiEilggl(l!)2=exp(Si/gEi)g1Eil=0(SiEi/g)l(l!)2=exp(αSiEi)g1Eil=01(l!)2(2SiEi/g2)2lI0(2SiEi/g)=exp(Si/gEi)g1EiI0(2SiEi/g)
(32)
where
I0(2SiEi/g)=l=01(l!)2(2SiEi/g2)2l
(33)
is the modified Bessel function of the first kind of order zero:
I0(z)=l=01(l!)2(z2)2l
(34)
Combining Eq. (15) with Eq. (28) and Eq. (32) we get:
n^i=n=0pSi|ni(Si|n;g)pn(n;Ei)nn=0pSi|ni(Si|n;g)pn(n;Ei)=exp(Si/gEi)g1EiI0(2SiEi/g)exp(Ei)δ(Si)+exp(Si/gEi)g1EiEiSi/gI1(2EiSi/g)={0ifSi=0SiEi/gI0(2SiEi/g)I1(2SiEi/g)ifSi>0
(35)
It follows from the considerations above, that Eq. (15) can be calculated by equation Eq. (17) and that it should be interpreted as zero in the special case where Si = 0.

Simulation of images for precision measurement

After the inner loop was applied to these photon images, the EMCCD amplification was simulated. An EMCCD gain of 300 was used. As in [13

13. K. I. Mortensen, L. S. Churchman, J. A. Spudich, and H. Flyvbjerg, “Optimized localization analysis for single-molecule tracking and super-resolution microscopy,” Nat. Methods 7, 377–381 (2010). [CrossRef] [PubMed]

], we drew random numbers from a continuous gamma distribution, approximating for the discrete distribution described by Eq. (2). The results were then rounded to determine the pixel’s final signal Si. We use the Apache Commons Mathematics Library implementation of the Poisson distribution and the Colt Library implementation of the gamma distribution.

Measuring computation time

The simulated images used to produce Fig. 4(a) were generated in the same way as the images used for precision measurement, but using a grid of emitters instead of a single one. The image sequences had a size of 11×11 pixels (1 emitter), 11×16 pixels (2 emitters), 16×16 pixels (4 emitters), 21×21 pixels (9 emitters), 21× 26 pixels (12 emitters), 26×26 pixels (16 emitters) and 26×31 pixels (20 emitters). The first emitter was set at pixel location (5,5) with all other emitters shifted in x and y direction by multiples of 5 pixels. The Gaussians used to model the PSF of the emitters was set to have a standard deviation of 1.3 pixels. Each emitter was set to have a flux of 100 photons. Background flux was set to 1 photon per pixel. A sequence of 10 frames was generated for each number of emitters. The images were post processed to account for EMCCD excess noise as was done for precision measurement.

To perform tracking, the appropriate number of emitters was initialized in the first frame on their true position and the initial pk0 set to 5%. It should be noted that this initialization is problematic when applied on the grid with 20 emitters, because background is implicitly set to 0. We found however, that localization was performed correctly in practice.

The localization results from one frame were used as initialization for the next. Computation time was averaged over all frames. The calculations were performed on a Intel® Core™ i7-3820 Processor. The program was run on a virtual machine with 4 simulated cores (using VMware Player). The guest system was Ubuntu 12.04 hosted by Windows 7.

Measuring robustness

Sequences of 1000 simulated 15×15 pixel images with were generated and post processed as for precision measurement. The flux of the single emitter was set to 100 photons. The background flux was set to 0.25 photons per pixel. For each frame the initial pk0 for the estimators was set to 5%. The initialization of the location was randomly determined in each frame using a normal distribution centered around the true position. Different standard deviations were used for each sequence.

Estimation of photon counts per pixel

The estimated photon counts per pixel indicated by color in Figs. 1 and 3 were computed according to the manufacturer’s information. First, the bias of 100 counts was subtracted from the signal (resulting values below zero were interpreted as zero). Then, the resulting values were divided by the EMCCD gain of 300 and finally multiplied by a sensitivity factor (given in the Andor iXon EM+ DU-897 performance sheet) of 12.17 electrons per analog/digital count.

To calculate the background flux used in the calculation of the theoretical bounds in Fig. 3 an area without beads was manually selected in each movie. The pixel value was then averaged over all frames and pixels in the area. Finally, as for the photon count estimation described above, the bias was subtracted, the value was divided by the EMCCD gain and finally multiplied by the sensitivity factor.

The photon count per pixel for the simulated EMCCD images in Fig. 2 was calculated by dividing the pixel’s signal by the EMCCD gain of 300.

Tracking of fluorescent beads

Isolated beads were selected from each movies first frame. The beads were tracked independently, starting from a manual initialization in the first frame. The initial value of pk in the first frame was 10% for all emitters. The standard deviation of the Gaussian to be used as approximate PSF was experimentally determined and set to be 1.3 pixels for all emitters. The cameras bias of 100 counts in each pixel was subtracted for the complete NMLA as well as for the inner loop alone. In case of the complete NMLA, each pixel’s signal was subsequently divided by the sensitivity factor (given in the Andor iXon EM+ DU-897 performance sheet) of 12.17 electrons per analog/digital count. This factor can be ignored when the inner loop alone is applied.

In order to correct the drift of the beads under the microscope during the 3000 frames long movies, the standard deviations displayed in Fig. 3 were calculated in the following way: (13000l=13000(f(l)xestl)2)12, where xestl is the estimated x-position in frame l and f (l) = a · l + b is a linear least squares fit to all estimated x-positions and their frame number.

Implementation details for NMLA

Each iteration of the NMLA was performed in a window around the last estimated position. The width of the square window was 6σ + 6 pixels, where σ is the standard deviation of the Gaussian approximating the emitter’s PSF. After each iteration, the window was shifted according to the new position estimate. The iteration was continued until the improvement of the log likelihood (as described in [13

13. K. I. Mortensen, L. S. Churchman, J. A. Spudich, and H. Flyvbjerg, “Optimized localization analysis for single-molecule tracking and super-resolution microscopy,” Nat. Methods 7, 377–381 (2010). [CrossRef] [PubMed]

] in supplementary equation 138) between two outer iterations dropped below a threshold of 0.001 or the maximum number of 30 iterations was reached.

The inner loop was employed in a window like the complete NMLA. In order to solve the optimization task in the M-step of the inner loop, we employed Powell’s method [19

19. M. J. D. Powell, “An efficient method for finding the minimum of a function of several variables without calculating derivatives,” Comput. J. 7, 155–162 (1964). [CrossRef]

] implemented in the Apache Commons Mathematics Library.

For the experiments regarding computation time and robustness we used the center of mass of all Ik,it as initialization in every M-step, provided it yielded an improvement regarding the current optimization task. After this initialization, Powell’s method was employed.

After each iteration, we calculated the change of parameters as max(|pktpkt1|, θktθkt1). The iterations in the inner loop were continued until the change of parameters dropped below a threshold of 0.001 or the maximum of 1000 iterations was reached.

Acknowledgments

We thank Britta Schroth-Diez, Daniel White and Nicola Maghelli for the discussion and help regarding EMCCD microscopy; Stephan Saalfeld, Stefan Gumhold, Ivo Sbalzarini and Pavel Tomancak for comments on the manuscript; Thomas Kühn, Tobias Pietzsch and Johannes Schindelin for advice regarding the software development with Imglib and Fiji; and Ivana Šarić for help with the preparation of the figures.

References and links

1.

M. Coelho, N. Maghelli, and I. M. Tolić-Nørrelykke, “Single-molecule imaging in vivo: the dancing building blocks of the cell,” Integr. Biol. (Camb) 5, 748–758 (2013). [CrossRef]

2.

K. Irie, A. E. McKinnon, K. Unsworth, and I. M. Woodhead, “A model for measurement of noise in CCD digital-video cameras,” Meas. Sci. Technol. 19, 045207 (2008). [CrossRef]

3.

D. J. Denvir and E. Conroy, “Electron-multiplying CCD: the new ICCD, in Low-Light-Level and Real-Time Imaging Systems, Components, and Applications”, Proc. SPIE 4796, pp. 164–174 (2003). [CrossRef]

4.

B. C. Carter, G. T. Shubeita, and S. P. Gross, “Tracking single particles: a user-friendly quantitative evaluation,” Phys. Biol. 2, 60–72 (2005). [CrossRef] [PubMed]

5.

G. M. Lee, A. Ishihara, and K. A. Jacobson, “Direct observation of brownian motion of lipids in a membrane,” Proc. Natl. Acad. Sci. USA 88, 6274–6278 (1991). [CrossRef] [PubMed]

6.

R. N. Ghosh and W. W. Webb, “Automated detection and tracking of individual and clustered cell surface low density lipoprotein receptor molecules.” Biophys. J. 66, 1301–1318 (1994). [CrossRef] [PubMed]

7.

C. M. Anderson, G. N. Georgiou, I. E. Morrison, G. V. Stevenson, and R. J. Cherry, “Tracking of cell surface receptors by fluorescence digital imaging microscopy using a charge-coupled device camera. low-density lipoprotein and influenza virus receptor mobility at 4 degrees c,” J. Cell Sci. 101 (Pt 2), 415–425 (1992). [PubMed]

8.

G. J. Schutz, H. Schindler, and T. Schmidt, “Single-molecule microscopy on model membranes reveals anomalous diffusion.” Biophys. J. 73, 1073–1080 (1997). [CrossRef] [PubMed]

9.

R. J. Ober, S. Ram, and E. S. Ward, “Localization accuracy in single-molecule microscopy,” Biophys. J. 86, 1185–1200 (2004). [CrossRef] [PubMed]

10.

C. S. Smith, N. Joseph, B. Rieger, and K. A. Lidke, “Fast, single-molecule localization that achieves theoretically minimum uncertainty.” Nat. Methods 7, 373–375 (2010). [CrossRef] [PubMed]

11.

F. Huang, S. L. Schwartz, J. M. Byars, and K. A. Lidke, “Simultaneous multiple-emitter fitting for single molecule super-resolution imaging,” Biomed. Opt. Express 2, 1377–1393 (2011). [CrossRef] [PubMed]

12.

Y. Wang, T. Quan, S. Zeng, and Z.-L. Huang, “PALMER: a method capable of parallel localization of multiple emitters for high-density localization microscopy,” Opt. Express 20, 16039–16049 (2012). [CrossRef] [PubMed]

13.

K. I. Mortensen, L. S. Churchman, J. A. Spudich, and H. Flyvbjerg, “Optimized localization analysis for single-molecule tracking and super-resolution microscopy,” Nat. Methods 7, 377–381 (2010). [CrossRef] [PubMed]

14.

J. Chao, E. S. Ward, and R. J. Ober, “Fisher information matrix for branching processes with application to electron-multiplying charge-coupled devices,” Multidimens. Syst. Signal Process. 23, 349–379 (2012). [CrossRef] [PubMed]

15.

C. H. A. G. Basden, “Photon counting strategies with low light level CCDs,” Monthly Notices of the Royal Astron. Soc. 345, 985–991 (2003). [CrossRef]

16.

J. Chao, E. S. Ward, and R. J. Ober, “Localization accuracy in single molecule microscopy using electron-multiplying charge-coupled device cameras, in Three-Dimensional and Multidimensional Microscopy: Image Acquisition and Processing XIX,” Proc. SPIE 8227,p. 82271P (2012). [CrossRef]

17.

A. P. Dempster, N. M. Laird, and D. B. Rubin, “Maximum likelihood from incomplete data via the EM algorithm,” J. Royal Stat. Soc., B 39, 1–38 (1977).

18.

C. Rao, Linear Statistical Inference and its Applications (Wiley, 1975).

19.

M. J. D. Powell, “An efficient method for finding the minimum of a function of several variables without calculating derivatives,” Comput. J. 7, 155–162 (1964). [CrossRef]

20.

J. Chao, S. Ram, E. S. Ward, and R. J. Ober, “Ultrahigh accuracy imaging modality for super-localization microscopy,” Nat. Methods 10, 335–338 (2013). [CrossRef] [PubMed]

21.

A. Steinborn, S. Taut, V. Brendler, G. Geipel, and B. Flach, “TRLFS: analysing spectra with an expectation-maximization (EM) algorithm.” Spectrochim Acta A Mol Biomol Spectrosc. 71, 1425–1432 (2008). [CrossRef] [PubMed]

22.

M. I. Schlesinger and V. Hlavac, Ten Lectures on Statistical and Structural Pattern Recognition, (Kluwer Academic Publishers, 2002) vol. 24 of Computational Imaging and Vision. [CrossRef]

23.

S. Stallinga and B. Rieger, “The effect of background on localization uncertainty in single emitter imaging,” Proceedings of 9th IEEE International Symposium on Biomedical Imaging 2012, 988–991(2012).

24.

L. Xu and M. I. Jordan, “On convergence properties of the EM algorithm for gaussian mixtures,” Neural Comput. 8, 129–151 (1995). [CrossRef]

25.

R. Starr, S. Stahlheber, and A. Small, “Fast maximum likelihood algorithm for localization of fluorescent molecules,” Opt. Lett. 37, 413–415 (2012). [CrossRef] [PubMed]

26.

M. Hirsch, R. J. Wareham, M. L. Martin-Fernandez, M. P. Hobson, and D. J. Rolfe, “A stochastic model for electron multiplication charge-coupled devices – from theory to practice,” PLoS ONE 8, e53671 (2013). [CrossRef]

27.

I. Kalinina, A. Nandi, P. Delivani, M. R. Chacón, A. H. Klemm, D. Ramunno-Johnson, A. Krull, B. Lindner, N. Pavin, and I. M. Tolić-Nørrelykke, “Pivoting of microtubules around the spindle pole accelerates kinetochore capture,” Nat. Cell Biol. 15, 82–87 (2013). [CrossRef]

28.

V. Ananthanarayanan, M. Schattat, S. K. Vogel, A. Krull, N. Pavin, and I. M. Tolić-Nørrelykke, “Dynein motion switches from diffusive to directed upon cortical anchoring,” Cell 153, 1526–1536 (2013). [CrossRef] [PubMed]

29.

J. Schindelin, I. Arganda-Carreras, E. Frise, V. Kaynig, M. Longair, T. Pietzsch, S. Preibisch, C. Rueden, S. Saalfeld, B. Schmid, J.-Y. Tinevez, D. J. White, V. Hartenstein, K. Eliceiri, P. Tomancak, and A. Cardona, “Fiji: an open-source platform for biological-image analysis,” Nat. Methods 9, 676–682 (2012). [CrossRef] [PubMed]

30.

S. K. Vogel, N. Pavin, N. Maghelli, F. Jülicher, and I. M. Tolić-Nørrelykke, “Self-organization of dynein motors generates meiotic nuclear oscillations,” PLoS Biol. 7, 0918–0928 (2009). [CrossRef]

OCIS Codes
(040.3780) Detectors : Low light level
(100.2960) Image processing : Image analysis
(110.0180) Imaging systems : Microscopy
(170.2520) Medical optics and biotechnology : Fluorescence microscopy
(150.1135) Machine vision : Algorithms

ToC Category:
Image Processing

History
Original Manuscript: September 9, 2013
Revised Manuscript: November 15, 2013
Manuscript Accepted: November 16, 2013
Published: January 2, 2014

Virtual Issues
Vol. 9, Iss. 3 Virtual Journal for Biomedical Optics

Citation
Alexander Krull, André Steinborn, Vaishnavi Ananthanarayanan, Damien Ramunno-Johnson, Uwe Petersohn, and Iva M. Tolić-Nørrelykke, "A divide and conquer strategy for the maximum likelihood localization of low intensity objects," Opt. Express 22, 210-228 (2014)
http://www.opticsinfobase.org/vjbo/abstract.cfm?URI=oe-22-1-210


Sort:  Author  |  Year  |  Journal  |  Reset  

References

  1. M. Coelho, N. Maghelli, I. M. Tolić-Nørrelykke, “Single-molecule imaging in vivo: the dancing building blocks of the cell,” Integr. Biol. (Camb) 5, 748–758 (2013). [CrossRef]
  2. K. Irie, A. E. McKinnon, K. Unsworth, I. M. Woodhead, “A model for measurement of noise in CCD digital-video cameras,” Meas. Sci. Technol. 19, 045207 (2008). [CrossRef]
  3. D. J. Denvir, E. Conroy, “Electron-multiplying CCD: the new ICCD, in Low-Light-Level and Real-Time Imaging Systems, Components, and Applications”, Proc. SPIE 4796, pp. 164–174 (2003). [CrossRef]
  4. B. C. Carter, G. T. Shubeita, S. P. Gross, “Tracking single particles: a user-friendly quantitative evaluation,” Phys. Biol. 2, 60–72 (2005). [CrossRef] [PubMed]
  5. G. M. Lee, A. Ishihara, K. A. Jacobson, “Direct observation of brownian motion of lipids in a membrane,” Proc. Natl. Acad. Sci. USA 88, 6274–6278 (1991). [CrossRef] [PubMed]
  6. R. N. Ghosh, W. W. Webb, “Automated detection and tracking of individual and clustered cell surface low density lipoprotein receptor molecules.” Biophys. J. 66, 1301–1318 (1994). [CrossRef] [PubMed]
  7. C. M. Anderson, G. N. Georgiou, I. E. Morrison, G. V. Stevenson, R. J. Cherry, “Tracking of cell surface receptors by fluorescence digital imaging microscopy using a charge-coupled device camera. low-density lipoprotein and influenza virus receptor mobility at 4 degrees c,” J. Cell Sci. 101 (Pt 2), 415–425 (1992). [PubMed]
  8. G. J. Schutz, H. Schindler, T. Schmidt, “Single-molecule microscopy on model membranes reveals anomalous diffusion.” Biophys. J. 73, 1073–1080 (1997). [CrossRef] [PubMed]
  9. R. J. Ober, S. Ram, E. S. Ward, “Localization accuracy in single-molecule microscopy,” Biophys. J. 86, 1185–1200 (2004). [CrossRef] [PubMed]
  10. C. S. Smith, N. Joseph, B. Rieger, K. A. Lidke, “Fast, single-molecule localization that achieves theoretically minimum uncertainty.” Nat. Methods 7, 373–375 (2010). [CrossRef] [PubMed]
  11. F. Huang, S. L. Schwartz, J. M. Byars, K. A. Lidke, “Simultaneous multiple-emitter fitting for single molecule super-resolution imaging,” Biomed. Opt. Express 2, 1377–1393 (2011). [CrossRef] [PubMed]
  12. Y. Wang, T. Quan, S. Zeng, Z.-L. Huang, “PALMER: a method capable of parallel localization of multiple emitters for high-density localization microscopy,” Opt. Express 20, 16039–16049 (2012). [CrossRef] [PubMed]
  13. K. I. Mortensen, L. S. Churchman, J. A. Spudich, H. Flyvbjerg, “Optimized localization analysis for single-molecule tracking and super-resolution microscopy,” Nat. Methods 7, 377–381 (2010). [CrossRef] [PubMed]
  14. J. Chao, E. S. Ward, R. J. Ober, “Fisher information matrix for branching processes with application to electron-multiplying charge-coupled devices,” Multidimens. Syst. Signal Process. 23, 349–379 (2012). [CrossRef] [PubMed]
  15. C. H. A. G. Basden, “Photon counting strategies with low light level CCDs,” Monthly Notices of the Royal Astron. Soc. 345, 985–991 (2003). [CrossRef]
  16. J. Chao, E. S. Ward, R. J. Ober, “Localization accuracy in single molecule microscopy using electron-multiplying charge-coupled device cameras, in Three-Dimensional and Multidimensional Microscopy: Image Acquisition and Processing XIX,” Proc. SPIE 8227,p. 82271P (2012). [CrossRef]
  17. A. P. Dempster, N. M. Laird, D. B. Rubin, “Maximum likelihood from incomplete data via the EM algorithm,” J. Royal Stat. Soc., B 39, 1–38 (1977).
  18. C. Rao, Linear Statistical Inference and its Applications (Wiley, 1975).
  19. M. J. D. Powell, “An efficient method for finding the minimum of a function of several variables without calculating derivatives,” Comput. J. 7, 155–162 (1964). [CrossRef]
  20. J. Chao, S. Ram, E. S. Ward, R. J. Ober, “Ultrahigh accuracy imaging modality for super-localization microscopy,” Nat. Methods 10, 335–338 (2013). [CrossRef] [PubMed]
  21. A. Steinborn, S. Taut, V. Brendler, G. Geipel, B. Flach, “TRLFS: analysing spectra with an expectation-maximization (EM) algorithm.” Spectrochim Acta A Mol Biomol Spectrosc. 71, 1425–1432 (2008). [CrossRef] [PubMed]
  22. M. I. Schlesinger, V. Hlavac, Ten Lectures on Statistical and Structural Pattern Recognition, (Kluwer Academic Publishers, 2002) vol. 24 of Computational Imaging and Vision. [CrossRef]
  23. S. Stallinga, B. Rieger, “The effect of background on localization uncertainty in single emitter imaging,” Proceedings of 9th IEEE International Symposium on Biomedical Imaging 2012, 988–991(2012).
  24. L. Xu, M. I. Jordan, “On convergence properties of the EM algorithm for gaussian mixtures,” Neural Comput. 8, 129–151 (1995). [CrossRef]
  25. R. Starr, S. Stahlheber, A. Small, “Fast maximum likelihood algorithm for localization of fluorescent molecules,” Opt. Lett. 37, 413–415 (2012). [CrossRef] [PubMed]
  26. M. Hirsch, R. J. Wareham, M. L. Martin-Fernandez, M. P. Hobson, D. J. Rolfe, “A stochastic model for electron multiplication charge-coupled devices – from theory to practice,” PLoS ONE 8, e53671 (2013). [CrossRef]
  27. I. Kalinina, A. Nandi, P. Delivani, M. R. Chacón, A. H. Klemm, D. Ramunno-Johnson, A. Krull, B. Lindner, N. Pavin, I. M. Tolić-Nørrelykke, “Pivoting of microtubules around the spindle pole accelerates kinetochore capture,” Nat. Cell Biol. 15, 82–87 (2013). [CrossRef]
  28. V. Ananthanarayanan, M. Schattat, S. K. Vogel, A. Krull, N. Pavin, I. M. Tolić-Nørrelykke, “Dynein motion switches from diffusive to directed upon cortical anchoring,” Cell 153, 1526–1536 (2013). [CrossRef] [PubMed]
  29. J. Schindelin, I. Arganda-Carreras, E. Frise, V. Kaynig, M. Longair, T. Pietzsch, S. Preibisch, C. Rueden, S. Saalfeld, B. Schmid, J.-Y. Tinevez, D. J. White, V. Hartenstein, K. Eliceiri, P. Tomancak, A. Cardona, “Fiji: an open-source platform for biological-image analysis,” Nat. Methods 9, 676–682 (2012). [CrossRef] [PubMed]
  30. S. K. Vogel, N. Pavin, N. Maghelli, F. Jülicher, I. M. Tolić-Nørrelykke, “Self-organization of dynein motors generates meiotic nuclear oscillations,” PLoS Biol. 7, 0918–0928 (2009). [CrossRef]

Cited By

Alert me when this paper is cited

OSA is able to provide readers links to articles that cite this paper by participating in CrossRef's Cited-By Linking service. CrossRef includes content from more than 3000 publishers and societies. In addition to listing OSA journal articles that cite this paper, citing articles from other participating publishers will also be listed.

Figures

Fig. 1 Fig. 2 Fig. 3
 
Fig. 4
 

« Previous Article  |  Next Article »

OSA is a member of CrossRef.

CrossCheck Deposited