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
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 . 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)).
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–14]. Other kinds of noise include dark signal and clock induced charges . 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 . 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 . One simple approach is the calculation of the center of mass of the signal [5, 6]. 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,8], 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 .
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 asEq. (1).
Over the past years, this scheme and derivations have been used for localization in a number of studies, e.g. [9–12]. In  Ober et al. formulated the problem and calculated an absolute limit for accuracy in the presence of additional readout noise. In  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, 12].
The ML approach has recently been applied with a different noise model including excess noise caused by EMCCD amplification . We will refer to the ML location estimator based on this model as excess noise aware ML estimator. The statistics of the amplification and its consequences regarding ML localization have been thoroughly discussed by Chao et al. . We will base our considerations on the simplified model suggested by Basden et al. , which was recommended for high amplification gains in  and used in .
During EMCCD amplification, photo-electrons are shifted through a multiplication register. According to Basden et al.  the number of output electrons Si given a certain number of input electrons ni approximately follows a discrete version of an Erlang distribution:16]. The effects of readout noise become insignificant with high EMCCD gains and we therefore omit them [3, 10]. 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.
When multiple emitters in close proximity are to be localized, the search for the ML estimate (excess noise aware as well as unaware) of their positions is a challenging high dimensional optimization task, because it has to be done jointly for all emitters. While  and  discuss only localization of single emitters, the localization tools from [10–12] achieve parallelization only by prior division of the image into subregions, which are then processed independently. It has not been recognized that the ML localization problem has a structure that allows for the use of a divide and conquer (DC) strategy to achieve higher stability and parallelization on a different level.
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.  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) . We compare the results and required computation time of the NMLA with a well established general purpose optimization method  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.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, 20], but was not utilized for localization. A comparable view has been applied to spectra analysis .
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 . 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: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 . 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:
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  and is a conventional choice used among others in [10–12].
2.2. Localization with the inner loop
We will now describe the inner loop of our algorithm which is capable of performing ML localization on multiple emitters. By utilizing Eq. (3), it is possible to independently acquire ML estimates for the model parameters
- The inner E-step: Regard the current parameters θ⃗t of the mixture distribution as correct and use it to estimate how much of the intensity of each pixel belongs to each light source and how much to the background.
- The inner M-step: Regard the estimates of the E-step as correct and use them to independently determine the parameters and of each light source.
It has been shown that the estimates obtained by the EM scheme improve with each iteration until convergence . 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:
- The inner M-step. The estimates for the a priori probabilities are acquired in the following way: 19]. 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:
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:
Accordingly, the flux of a specific light source can be estimated up to the factor c as . 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 . 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:
We again apply the EM scheme to find the solution. Starting again with an initial estimate by the user or from a previous frame (θ⃗0, E0), the algorithm creates an iteratively improving sequence of tuples that converge to an optimum. As in the optimization of the pure shot noise model, the new estimates are in each iteration calculated by executing two steps:
- The outer E-step: Regard the current parameters (θ⃗t, Et) as correct and use them to estimate each pixel’s expected photon count , given the pixel’s observed signal Si.
- The outer M-step: Regard the results of the E-step as correct and use them to re-estimate the parameters θ⃗t+1 as well as the total image flux Et+1 and the flux for each light source .
- The outer E-step. Based on the current parameters (θ⃗t, Et) we calculate the expected number of photons for each pixel i, given the observed signal Si. It is is defined as: Eq. (15) efficiently as Appendix.
- The outer M-step. Even though the estimates for the model and the mean photon count have to be performed together, in each iteration of the algorithm their new estimates can be found independently of one another. The tasks are essentially equivalent to tasks described in Eq. (5) and Eq. (6). We will first describe the optimization of the emitter parameters θ⃗t+1
Optimizing the estimate of the emitter parameters: Using the expected values calculated in the previous step, the next parameter estimate is defined as:Eq. (5). Indeed the new parameter estimates θ⃗t+1 can be acquired by employing the inner loop described above. The only modification necessary is a substitution of the image intensities Ii by the current expected photon count . This constitutes a nested EM algorithm in which the M-step consists of an EM algorithm itself.
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):
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)) . 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  that the approximation of the CRLB for excess noise, obtained by multiplying the CRLB for a pure shot noise setting by a factor of is not appropriate for low photon counts. Surprisingly, this approximation from  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).
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 . 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  to the log likelihood as described in  (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.
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.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–12], likely because the theoretical understanding of the problem has only recently been achieved [13,14]. 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  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:
First, our algorithm keeps computation time low even as the number of overlapping emitters is increased (Fig. 4(a)). Second, since the algorithm allows independent re-estimation of emitter locations during the inner M-step, it allows for easy coarse level parallelization. This could facilitate GPU implementations. While previous methods [10–12] rely on the division of the image into subregions in a preprocessing step, NMLA becomes faster through parallelization even when all emitters are in close proximity and have to be localized jointly (Fig. 4(a)). Third, our algorithm proved to be significantly less sensitive to random initialization than the general purpose optimizer (Fig. 4(b)).
4.3. Conceptual advances
We achieve the level of robustness and scalability shown in Fig. 4 by reformulating the optimization task in a way that allows for the use of EM. The task of localization in a shot noise setting is very closely related to the problem of Gaussian mixture estimation and EM is well established for this problem and demonstrates solid convergence properties . In addition, it does not require the calculation of derivatives, in contrast to many general purpose optimization methods.
It should also be noted, that our novel formulation of the optimization task defined in Eq. (5) is of lower dimensionality than the previous formulation which can be found in [9–12]. Our approach still yields the same optimum. In all previous formulations, apart from the emitter locations, m + 1 parameters namely the flux of each emitter and the flux of the background had to be determined. Apart from the emitter locations, our formulation requires estimation of only m parameters, namely each emitter’s percentage of light pk. This is due to the fact, that p0 can be calculated directly from the other values pk.
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.
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 , 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 , but it could be included by calculating pSi|ni (Si|ni; g) as convolution of Eq. (2) and a Gaussian distribution . 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 . Dark signal is pixel dependent  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 . According to  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–12] 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.
ML estimation in EMCCD images and its potential in low-light microscopy has accrued considerable theoretical understanding. . 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, 28].
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 .
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 . 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:
Calculating the CRLB
We calculated the CRLB as the inverse of the Fisher information matrix . For the pure shot noise setting the elements of the Fisher information matrix were calculated as described in . For the case with excess noise, the elements of the Fisher information matrix can be computed asEq. (22) can be found in . As done in  we calculate q(Si; Ei) efficiently as
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:Eq. (15): Eq. (15): Eq. (15) with Eq. (28) and Eq. (32) we get: 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
To calculate the standard deviations shown in Fig. 2(a) we generated a series of 250 simulated 15×15 pixel images for each value of flux. For Fig. 2(b) a series of 1000 simulated 15×15 pixel images was created for each value of background flux. In both cases a single emitter was placed in the center of each image. The standard deviation of the Gaussian used as PSF was set to be 1.3 pixels. Photon images were generated by randomly drawing each pixel’s photon count ni from a Poisson distribution parametrized by the pixel’s flux Ei. We calculated the flux as combination of the Gaussian integral over the pixel’s area and the per pixel background flux. The algorithms were initialized in the center of the first frame. The initial value of pk in the first frame was 10% for all emitters. All algorithms subsequently used the results of one frame as initialization for the next one.
After the inner loop was applied to these photon images, the EMCCD amplification was simulated. An EMCCD gain of 300 was used. As in , 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 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.
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 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: , where 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  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  implemented in the Apache Commons Mathematics Library.
For the experiments regarding computation time and robustness we used the center of mass of all 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 , . 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.
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]
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]
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]
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]
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]