State space approach to single molecule localization in fluorescence microscopy

Single molecule super-resolution microscopy enables imaging at sub-diffraction-limit resolution by producing images of subsets of stochastically photoactivated fluorophores over a sequence of frames. In each frame of the sequence, the fluorophores are accurately localized, and the estimated locations are used to construct a high-resolution image of the cellular structures labeled by the fluorophores. Many methods have been developed for localizing fluorophores from the images. The majority of these methods comprise two separate steps: detection and estimation. In the detection step, fluorophores are identified. In the estimation step, the locations of the identified fluorophores are estimated through an iterative approach. Here, we propose a non-iterative state space-based localization method which combines the detection and estimation steps. We demonstrate that the estimated locations obtained from the proposed method can be used as initial conditions in an estimation routine to potentially obtain improved location estimates. The proposed method models the given image as the frequency response of a multi-order system obtained with a balanced state space realization algorithm based on the singular value decomposition of a Hankel matrix. The locations of the poles of the resulting system determine the peak locations in the frequency domain, and the locations of the most significant peaks correspond to the single molecule locations in the original image. The performance of the method is validated using both simulated and experimental data. © 2017 Optical Society of America OCIS codes: (170.2520) Fluorescence microscopy; (100.2960) Image analysis; (110.3010) Image reconstruction techniques; (000.5490) Probability theory, stochastic processes, and statistics; (100.6640 ) Superresolution. References and links 1. A. Small and S. Stahlheber, “Fluorophore localization algorithms for super-resolution microscopy,” Nat. Methods 11(3), 267–279 (2014). 2. 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(5), 1377–1393 (2011). 3. S. J. Holden, S. Uphoff, and A. N. Kapanidis, “DAOSTORM: an algorithm for high-density super-resolution microscopy,” Nat. Methods 8, 279–280 (2011). 4. T. A. Laurence and B. A. Chromy, “Efficient maximum likelihood estimator fitting of histograms,” Nat. Methods 7(5), 338–339 (2010). 5. R. Henriques, M. Lelek, E. F. Fornasiero, F. Valtorta, C. Zimmer, and M. M. Mhlanga, “QuickPALM: 3D real-time photoactivation nanoscopy image processing in ImageJ,” Nat. Methods 7(5), 339–340 (2010). 6. A. J. Berglund, M. D. McMahon, J. J. McClelland, and J. A. Liddle, “Fast, bias-free algorithm for tracking single particles with variable size and shape,” Opt. Express 16(18), 14064–14075 (2008). 7. D. L. Donoho, “Compressed sensing,” IEEE Trans. Inf. Theory 52(4), 1289–1306 (2006). 8. L. Zhu, W. Zhang, D. Elnatan, and B. Huang, “Faster STORM using compressed sensing,” Nat. Methods 9(7), 721–723 (2012). 9. T. Quan, H. Zhu, X. Liu, Y. Liu, J. Ding, S. Zeng, and Z. L. Huang, “High-density localization of active molecules using Structured Sparse Model and Bayesian Information Criterion,” Opt. Express 19(18), 16963–16974 (2011). 10. J. Min, C. Vonesch, H. Kirshner, L. Carlini, N. Olivier, S. Holden, S. Manley, J. C. Ye, and M. Unser, “FALCON: fast and unbiased reconstruction of high-density super-resolution microscopy data,” Sci. Rep. 4, 4577 (2014). Vol. 8, No. 3 | 1 Mar 2017 | BIOMEDICAL OPTICS EXPRESS 1332 #278546 Journal © 2017 https://doi.org/10.1364/BOE.8.001332 Received 13 Oct 2016; revised 14 Jan 2017; accepted 30 Jan 2017; published 6 Feb 2017 11. J. Huang, K. Gumpper, Y. Chi, M. Sun, and J. Ma, “Fast two-dimensional super-resolution image reconstruction algorithm for ultra-high emitter density,” Opt. Lett. 40(13), 2989–2992 (2015). 12. Y. Chi, L. L. Scharf, A. Pezeshki, and A. R. Calderbank, “Sensitivity to basis mismatch in compressed sensing,” IEEE Trans. Signal Process. 59(5), 2182–2195 (2011). 13. Y. Hua, “Estimating two-dimensional frequencies by matrix enhancement and matrix pencil,” IEEE Trans. Signal Process. 40(9), 2267–2280 (1992). 14. R. J. Ober, X. Lai, Z. Lin, and E. S. Ward, “A state space approach to noise reduction of 3D fluorescent microscopy images,” in Proc. IEEE International Conference on Image Processing (IEEE 2004), pp. 1153–1156. 15. X. Lai, E. S. Ward, Z. Lin, and R. J. Ober, “Three-dimensional state space realization algorithm: noise suppression of fluorescence microscopy images and point spread functions,” Proc. SPIE 5701, 53–60 (2005). 16. R. J. Ober, X. Lai, Z. Lin, and E. S. Ward, “State space realization of a three-dimensional image set with application to noise reduction of fluorescent microscopy images of cells,” Multidim. Syst. Sign. P. 16(1), 7–47 (2005). 17. S. Ram, E. S. Ward, and R. J. Ober, “A stochastic analysis of performance limits for optical microscopes,” Multidim. Syst. Sign. P. 17(1), 27–57 (2006). 18. T. McKelvey, H. Akçay, and L. Ljung “Subspace-based multivariable system identification from frequency response data,” IEEE Trans. Automatic Control 41(7), 960–979 (1996). 19. J. M. Maciejowski, “Guaranteed stability with subspace methods,” Systems & Control Letters 26(2), 153–156 (1995). 20. D. Sage, H. Kirshner, T. Pengo, N. Stuurman, J. Min, S. Manley, and M. Unser, “Quantitative evaluation of software packages for single-molecule localization microscopy,” Nat. Methods 12(8), 717–724 (2015). 21. A. V. Abraham, S. Ram, J. Chao, E. S. Ward, and R. J. Ober, “Quantitative study of single molecule location estimation techniques,” Opt. Express 17(26), 23352–23373 (2009). 22. R. J. Ober, S. Ram, and E. S. Ward, “Localization accuracy in single-molecule microscopy,” Biophys. J. 86(2), 1185–1200 (2004). 23. A. Tahmasbi, E. S. Ward, and R. J. Ober, “Determination of localization accuracy based on experimentally acquired image sets: applications to single molecule microscopy,” Opt. Express 23(6), 7630–7652 (2015).


Introduction
Single molecule super-resolution methods have been successful at achieving sub-diffraction-limit resolution based on two key innovations: photoactivable fluorophores and powerful fluorophore localization algorithms [1].In these methods, a fluorescently labeled cellular structure is imaged over a sequence of frames.In each frame, only a small number of stochastically photoactivated fluorophores are detected.In order to construct a high-resolution image of the cellular structure, the locations of individual emitting fluorophores are estimated with sub-pixel precision from each frame and used to re-render the structure.Many fluorophore localization methods are available, and they typically comprise the following separate steps: a detection step that identifies the regions in the image that contain emitting fluorophores, and an estimation step that determines the locations of these fluorophores accurately.In recent years, several methods have been developed that use fitting-based algorithms to solve the estimation problem.The basis of most of these methods is to fit a point spread function (PSF) model to the acquired data and estimate the parameters of the model by minimizing the difference between the data and the model through an iterative approach.For example, in [2], a method was proposed that uses the maximum likelihood estimator to localize multiple emitters simultaneously within a two-dimensional (2D) fitting subregion.Similarly, the DAOSTORM algorithm [3] fits multiple PSFs in a recursive approach by analyzing pixel clusters in the residual image.In this algorithm, the fluorophore locations are determined by minimizing a least-squares criterion.
Fitting-based algorithms are not the only approaches used to solve the estimation problem for multi-emitter images.Many other localization methods have been developed that use non-fitting algorithms for the estimation problem.These methods are preferable when accurate PSF and noise models are not available.As an example, the QuickPALM software uses the simple centroid method [4,5].In this method, the fluorophore location is estimated as the average photon location, or centroid.However, the image background causes a systematic deviation in centroid-based methods.To solve the background bias problem in centroid-based methods, the virtual window center of mass (VWCM) method has been demonstrated to be a good background-corrected centroid estimator [6].Although centroid methods are fast and computationally simple, their accuracy is not comparable to that of good fitting-based methods.Another important class of non-fitting algorithms have been developed based on sparse support recovery methods [7].The compressive-sensing-based method CSSTORM [8], structured sparse model and Bayesian information criterion (SSC-BIC) [9], and fast localization algorithm based on a continuousspace formulation (FALCON) [10] are well-known examples of such algorithms.Among them, CSSTORM has been shown to achieve accurate localization for emitter densities as high as 10 emitters/µm 2 [11].In this method, a large-scale convex optimization problem needs to be solved in an iterative approach [12].Huang et al. [11] have proposed a non-fitting algorithm by transferring the molecule localization problem to the frequency domain.Their proposed algorithm is based on a 2D spectrum-estimation method called matrix enhancement and matrix pencil (MEMP) [13].MEMP can provide a significant speed advantage over CSSTORM while retaining the same level of accuracy (it is 100 times faster than CSSTORM with l 1 -homotopy [11]).Huang et al., however, assume that the PSF can be approximated by a Gaussian function, which can be problematic in practice due to the fact that the Gaussian model is often not an accurate analytical PSF.
Here, we propose a non-fitting state space algorithm for single molecule localization which combines the detection and estimation steps in a non-iterative approach.Generally, a single molecule fluorescence image contains multiple peaks of intensity that correspond to emitting fluorophores.Our solution is to model such an image by the frequency response of a multi-order system, as the locations of the poles of such a system determine the peak locations in the frequency domain.To realize this localization algorithm, we take advantage of the balanced state space realization algorithm used in [14][15][16] for the reduction of noise in fluorescence microscopy images.This realization algorithm is based on the singular value decomposition (SVD) of a Hankel matrix.To associate the peak locations in the image with the poles of the underlying system, we apply this realization algorithm to the inverse Fourier transform of the image rather than to the image itself.In our algorithm, the number of emitting fluorophores, which correspond to the most significant peaks in the image, is ultimately determined using a procedure that utilizes a least-squares criterion.Our algorithm also allows us to derive a theoretical reconstruction of the image.A reconstructed image is an image that looks similar to the original image, but is specified analytically in terms of the state space parameters of the system calculated using our proposed localization algorithm.
Note that the realization algorithm from [16] was developed for the reduction of noise in a three-dimensional (3D) data set comprising a z-stack of microscopy images.Here, we apply a 2D version of that realization algorithm because the super-resolution microscopy data sets for which we develop our localization algorithm consist of 2D images that are analyzed independently of one another.Our localization algorithm, however, can be extended to the 3D localization of fluorescent emitters from a z-stack by simply applying the 3D version of the realization algorithm as presented in [16].
To assess the performance of the proposed localization algorithm, we apply the algorithm to both simulated and experimental data comprising images of closely spaced molecules.In the case of simulated data, we evaluate the detection rate of the algorithm for molecules with different mean photon counts and for different distances between the molecules.We also analyze the bias of the algorithm.The bias is evaluated as the average of the deviations of the estimated molecule locations from the ground truth.In the case where there is only one molecule per image, our results suggest that there is no systematic bias associated with the algorithm.In the case of data sets consisting of repeat images of multiple molecules, however, the results show that bias exists which we have found to be dependent on the distances between the molecules relative to the image size.Also, the accuracy of the algorithm, determined by how far the estimates are spread out from the ground truth, is assessed by looking at the square root of the average of the squared deviations from the ground truth.In the case that we have repeat images of the same molecules, we look at the squared deviations of the estimates from their average, i.e., the accuracy is given by the standard deviation of the estimates.Importantly, for data sets comprising repeat images of one molecule, the standard deviation of the estimates is compared with the limit of the localization accuracy, a theoretical accuracy benchmark given by the square root of the Cramér-Rao lower bound (CRLB) [17].The results show that the accuracy of the algorithm is reasonable, but the difference between the accuracy and the limit of accuracy is nevertheless around twice the limit of accuracy.We show, however, that by using the obtained location estimates as the initial conditions for a maximum likelihood estimator, we can decrease the standard deviation of the estimates and approach the limit of accuracy, as is usually possible with the maximum likelihood estimator for standard single molecule estimation problems.We further demonstrate with experimental data that the algorithm can recover the locations of the significant peaks in the original image that correspond to the locations of individual Alexa Fluor 647 dye molecules.
This paper is organized as follows.In Section 2, we show the existence of minimal and asymptotically stable systems that realize a 2D image in the frequency domain.Section 3 is devoted to summarizing our overall localization approach, developed based on the systems introduced in Section 2. In Section 4, we explain the overall approach proposed in Section 3 in more detail, and develop a state space-based localization algorithm based on the SVD of a Hankel matrix.In Section 5, we specify the parameters used to generate simulated image data, and describe the sample preparation procedure and the microscopy setup used to produce experimental image data.Lastly, the results of applying our proposed algorithm to simulated and experimental single molecule images are reported and discussed in Section 6.

System identification using frequency measurements
In this section, we show the existence of minimal and asymptotically stable systems that realize a finite 2D sequence in the frequency domain.We begin by demonstrating the existence of minimal and asymptotically stable systems for finite one-dimensional (1D) sequences in Lemma 1, using a subspace-based method similar to that described in [18], and then extend the result to two dimensions in Theorem 1.The basis of Lemma 1 is given by Proposition 1, which states that a finite 1D data set can be expressed as the impulse response of a minimal and asymptotically stable system [16].Note that the stability of subspace methods was analyzed previously by Maciejowski in [19], where similar results were reported.Proposition 1.For positive integer N, let X(n) ∈ C p×m , p, m ∈ N, n = 1, 2, ..., N, be a 1D matrix-valued sequence.Then, there exists a minimal and asymptotically stable system (A, B, C), such that Proposition 1 enables us to write the following lemma, which shows the existence of a minimal and asymptotically stable system that realizes a finite 1D sequence in the frequency domain.
be the inverse discrete Fourier transform (inverse DFT, or IDFT) of X.Then, there exists a minimal and asymptotically stable system Moreover, where be the IDFT of X.Then, according to Proposition 1, there exists a minimal and asymptotically stable system (A, B, C), such that According to Eqs. ( 4) and ( 5), we then have, for k = 1, 2, ..., N, For a square matrix T ∈ C m×m , m ∈ N, where the number 1 is not an eigenvalue of T, we have the identity N −1 n=0 T n = (I − T) −1 (I − T N ).Then, since the realization (A, B, C) is asymptotically stable, i.e., |λ(A)| < 1 holds for any eigenvalue λ(A) of A, the number 1 is not an eigenvalue of Ae −i2πk/N , k = 1, ..., N (or equivalently, I − Ae −i2πk/N , k = 1, ..., N, is invertible), and Substituting this expression into Eq.( 6), we have, for k = 1, 2, ..., N, where Ã := A, B := In the following theorem, we extend the results obtained for 1D sequences to 2D sequences.
be the inverse 2D DFT of X.Then, there exist minimal and asymptotically stable systems where, for i = 1, 2, Moreover, where, for k j = 1, 2, ..., N j , j = 1, 2, where Ãj := A j , Bj := (I − A be the inverse 2D DFT of X. Arrange the entries of X to form a matrix Q as Decompose Q via SVD as Q = UΣV, where for r ∈ N, U ∈ C N 1 ×r , Σ ∈ C r×r and V ∈ C r×N 2 . For Then Moreover, according to Proposition 1, there exist minimal and asymptotically stable systems According to Eqs. ( 13) and ( 16), where Xi where Ãj := A j , Bj :

Location estimation
So far, we have shown the existence of minimal and asymptotically stable systems in the frequency domain.Here, we summarize our overall localization approach.Given that X is a single molecule image with multiple peaks of intensity, we determine the locations of the molecules by calculating the pole locations of In Theorem 1, we have shown that there exist minimal and asymptotically stable systems where, for where Ãj := A j , Bj : then the diagonal elements of the resulting matrix Āi give the poles of the system.In the following, we use matrix diagonalization in Eq. ( 20) to express X in terms of the poles of the system.
For s 1 , s 2 ∈ N and A i ∈ C s i ×s i , i = 1, 2, which are diagonalizable, i.e., for t i = 1, 2, ..., s i , i = 1, 2, and some invertible T i ∈ C s i ×s i , we have the diagonal matrix Āi := we can write X in terms of the poles of the system as Equation ( 22) provides an analytical expression for the reconstructed image, in which the poles of X occur at (a is not associated with a peak in the corresponding 2D image.)We next obtain the locations of the molecules in the object space in terms of the phase of the calculated poles.Let the 2D sequence X denote the pixel intensities of our N 1 × N 2 image with pixel width ∆x and pixel height ∆y, obtained by sampling the image at the center of each pixel.Assume Then, by linearly mapping a 2π × 2π square region in the frequency domain to the region with area N 1 × N 2 pixels in the image space (between the center of the first pixel and the center of the last pixel) and converting from image space units to object space units, the set containing the peak locations (i.e., the molecule locations) in the object space is given by (x t 2 , y t 1 ) : , where and M > 0 denotes the lateral magnification of the microscope system.

Algorithm
We now explain our proposed approach in more detail.In Section 3, we have calculated the poles of a 2D single molecule image X in terms of the elements of minimal and asymptotically stable systems Here, using the balanced state space realization algorithm introduced by Maciejowski [19], we propose a step-by-step algorithm to calculate systems , that realize X, and to determine the locations of the single molecules using the realization.
, represent the acquired image data.I. Subtract an estimated background level β, e.g., the average of the data points near the boundary of the image data X, from the image data X, and define the background-subtracted image Xbs as II.Let X be the 2D IDFT of Xbs , i.e., III. Arrange the entries of X to form a matrix Q as , denote the number of retained singular values (see Section 4.1).
IV. Construct the Hankel matrices ) as where 0 denotes a block of zeros of the corresponding size.For i = 1, 2, decompose H i via SVD as )r×s 2 , and where ∈ C s j ×s j , j = 1, 2, i.e., for t j = 1, 2, ..., s j , j = 1, 2, and some invertible Note that in theory, there is a possibility that A r;s 1 1 and/or A r;s 2 2 are not diagonalizable.In practice, however, because the diagonalization is numerically computed, A r;s 1 1 and/or A r;s 2 2 can be expected to be diagonalizable.In the unlikely scenario where A r;s 1 1 and/or A r;s 2 2 are not diagonalizable, very small perturbations of the data can be introduced to alter slightly their eigenvalues and make them diagonalizable.A perturbation can be achieved, for example, by simply adding a very small value to a pixel of the image.
VI.For h = min(s 1 , s 2 ), calculate, in the object space, the estimated peak locations (x k , y k ), x k ∈ x1 , ..., xs 2 , y k ∈ ŷ1 , ..., ŷs 1 , k = 1, ..., h, where where ∆x and ∆y are the width and height of each pixel of the image, respectively, and M > 0 denotes the lateral magnification of the microscope system.
The proposed algorithm crucially depends on SVD.Most of the singular values resulting from an SVD are relatively small and are considered to correspond to noise [16].Here, an important question is how many singular values are associated with noise and should be discarded in each SVD?In the following subsections, we describe the determination of the number of retained singular values in the three SVDs of the algorithm, and importantly, the number of single molecules in the given image.In addition, we give a description of the maximum likelihood estimator with which we will demonstrate the use of the results of the algorithm as the initial conditions for an estimation routine.

Determination of the number of retained singular values in the first SVD
Let σ 1 ≥ ... ≥ σ K ≥ 0, K = min(N 1 , N 2 ), denote the singular values in the first SVD.For r = 1, ..., K, let E r := r i=1 σ 2 i be the energy of the sequence σ i , i = 1, ..., r.Estimate the optimal number of retained singular values r in the first SVD as r = min r=1,...,K r : where E K := K i=1 σ 2 i is the energy of the sequence of all singular values and τ ∈ R denotes a threshold value typically chosen in the range [0.8, 0.9] [11].In Section 6.1.3,we examine the effect of different threshold values on the detection rate of the algorithm.
, be the singular values in the second and third SVDs, respectively.For where 2 are the energies of the sequences σ i 1 , ..., σ i l i and σ i 1 , ..., σ i N i , respectively, and τ i ∈ R denotes a threshold value which is again typically chosen in the range [0.8, 0.9] (see Section 6.1.3).The estimates li , i = 1, 2, thus denote the number of singular values that remain after discarding those that are considered to obviously correspond to noise.
We next try to reduce further the number of singular values to retain using an optimization approach that minimizes the difference between the original image and the reconstructed image obtained by the estimated locations of the peaks of the image.For , be the estimated data calculated via the algorithm by retaining r singular values in the first SVD and s 1 and s 2 singular values in the second and third SVDs, respectively.In other words, Xr;s 1 ,s 2 is the reconstructed image of Eq. ( 22) after discarding the singular values corresponding to noise.Denoting the poles of Xr;s 1 ,s 2 by ( ā1 and their corresponding product of coefficients in the numerator by p k ∈ C, k = 1, ..., s 1 s 2 , assume the peak magnitudes to be Let h = min(s 1 , s 2 ) denote the number of single molecules, assuming that we retain s 1 and s 2 singular values in the second and third SVDs, respectively.In the following, we estimate the optimal number of single molecules.Let θh := ( θ1 , ..., θh are the estimated locations of the h peaks with the largest magnitudes.In general, we consider all possible h-combinations of the poles of Xr;s 1 ,s 2 , but in most cases, the single molecules are associated with the peaks with the largest magnitudes.Let z 1 , ..., z N pi x denote our acquired data, where N pix denotes the number of pixels in the image.Then, the estimated number of single molecules ĥ is given by ĥ = arg min where, in the case that the single molecule image is modeled with a 2D PSF, we have, for k = 1, ..., N pix , the mean number of photons detected in the k th pixel given by [17] µ θh (k where N p,n is the expected number of photons due to the n th molecule that impact the detector plane during the image exposure, C k ⊂ R 2 denotes the region in the detector plane occupied by the k th pixel, and q is the 2D PSF of the optical system.If the PSF is the Airy profile, then q is given by q(x, y) := where n a denotes the numerical aperture of the objective lens, λ denotes the emission wavelength of the molecule, and J 1 denotes the first order Bessel function of the first kind.

Fitting single molecule images using the maximum likelihood estimator
The molecule locations estimated with the proposed algorithm can be used as the initial conditions in any estimation routine.In this paper, we demonstrate this using the maximum likelihood estimation routine [21,22].In the following, we briefly explain the basis of the maximum likelihood estimation.Let Θ denote the parameter space that is an open subset of R n .The maximum likelihood estimate θmle of θ ∈ Θ, for data incorporating Gaussian readout noise, is given by θmle = arg min where z 1 , ..., z N pi x denotes an image with N pix pixels and L is the log-likelihood function given by [17] L(θ|z 1 , ..., In Eq. (39), in the case of one molecule, µ θ (k) is the mean photon count in the k th pixel due to the molecule and is given by where θ = (x 0 , y 0 ) ∈ R 2 denotes the location of the molecule in the object space, N p is the expected number of photons from the molecule that are detected over the detector plane, and q is the Airy profile given by Eq. (37).Also, β k is the background level in the k th pixel, and η k and σ k denote the mean and standard deviation of the Gaussian readout noise in the k th pixel, respectively.

Simulation parameters
To analyze the performance of the proposed algorithm, we simulated different data sets using parameters commonly used in single molecule experiments.Some data sets comprise repeat images of one molecule, and some comprise repeat images of more than one molecule.Also, some data sets are such that each image contains a different set of molecules whose locations are randomly chosen based on uniform distributions that place the molecules within different spatial intervals inside the image.Regardless of the data set, the image of a molecule was generated with the Airy profile of Eq. (37) with a numerical aperture of n a = 1.4 and an emission wavelength of λ = 485 nm.Furthermore, a lateral magnification of M = 100, a detector pixel size of 6.5 µm × 6.5 µm, and a zero-mean Gaussian readout noise with standard deviation σ = 6 e − per pixel, were assumed.Also, we assumed that all simulated images are background-subtracted.

Sample preparation
High-performance Zeiss coverslips (#1.5) were prepared as follows: coverslips were sonicated with the following solutions in succession (each for 20 minutes): 50% HPLC-grade ethanol, 1mM HCl with 50% HPLC-grade ethanol, 1M KOH with 50% HPLC-grade ethanol, and 50% HPLC-grade ethanol.The cleaned coverslips were attached to MatTek dishes.200 µl of Poly-L-lysine (PLL) solution (Sigma-Aldrich) were added to the glass bottom area of the dishes for 10 minutes at room temperature.PPL was removed and 250-pM Alexa Fluor 647 fluorescent dye (Invitrogen) in 200 µl of phosphate-buffered saline (PBS) was applied for 10 minutes at room temperature.The sample was then washed with PBS twice at room temperature followed by the addition of 1 ml of PBS to the sample.

Microscopy setup
Custom laser excitation optics were installed for a Zeiss Axio Observer.A1 microscope.The laser optics were configured with 635-nm and 405-nm diode lasers (OptoEngine) for the excitation and photoactivation, respectively, of Alexa Fluor 647.The excitation light was reflected using a dichroic filter (Di01-R405/488/561/635-25x36; Semrock) and focused on the back focal plane of a 63×, 1.46 NA Zeiss objective lens.The emission light from the Alexa Fluor 647 dye was collected by the objective lens and filtered using a single bandpass filter (FF01-676/29-25; Semrock).The images were recorded using an electron-multiplying charge-coupled device camera (iXon DU897-BV; Andor) in conventional readout mode.The camera pixel size was 16 µm × 16 µm.
All components, including lasers, shutters and cameras, were controlled and synchronized using custom software written in the C programming language.

Super-resolution imaging
We first removed PBS from the single molecule sample prepared in Section 5.2.1 and added the imaging buffer (50-mM beta-mercaptoethylamine (MEA), 0.5-mg/ml glucose oxidase, 40-µg/ml catalase in PBS, pH 7.4, with 10% glucose).The sample was sealed with a coverslip and then positioned on the sample stage of the microscope for 5 to 10 minutes for temperature equilibration and the oxygen scavenging process.Images were subsequently acquired at a rate of 20 frames per second.The sample was illuminated with the 635-nm and 405-nm diode lasers alternately with photoactivation by the 405-nm laser every third frame.The frames with 405-nm laser illumination were excluded from data analysis.

Results and discussion
In this section, we present and discuss the results of the proposed algorithm when applied to both simulated and experimental images of single molecules.

Results for simulated data
Using simulated data sets, we first examine the performance of the algorithm in terms of the detection rate.We then analyze the bias and accuracy of the algorithm.The bias is assessed by the average of the deviations of the estimated molecule locations from the ground truth.The accuracy is assessed by looking at the square root of the average of the squared deviations from the ground truth.For repeat images of the same molecules, however, we look instead at the standard deviation of the estimates.In particular, for data sets containing repeat images of one molecule, we compare the standard deviation of the estimates with the limit of the localization accuracy given by the square root of the CRLB.Besides these analyses, we examine the effect on the detection rate when different threshold values are used in the first, second, and third SVDs of the algorithm.

One molecule
Here, to evaluate the detection rate of the algorithm, we first simulated data sets in which each image contains one molecule, whose location was randomly chosen based on a uniform distribution that places it within the image.For a given data set, the mean photon count is the same for the molecule in every image.Different data sets differ by this mean photon count, which ranges from 500 to 4500.For each mean photon count, we simulated 200 images.To establish statistical measures of the detection rate, we needed to pair the molecules localized by the algorithm with the molecules from the ground truth.For this purpose, we used the Hungarian algorithm with a search area of radius 100 nm [20].Then, we categorized the localized molecules which were successfully paired with ground truth molecules as true positives.The ground truth molecules that were not paired with any localized molecule and the localized molecules which were not paired with any ground truth molecule were categorized as false negatives and false positives, respectively.Denoting the number of true positives by T P, the number of false negatives by F N, and the number of false positives by FP, we define the precision (PRE) and recall (REC) measures as [20] PRE := T P FP + T P , REC := Figure 1 shows the results of the measures of the detection rate for data sets consisting of images containing one molecule each.It can be seen that for all mean photon counts considered, there are no false negatives and the recall is 1.Also, the figure shows that by increasing the mean photon count, the precision increases.However, it is important to note that even when the mean number of photons is as low as 500, a relatively large number of detected molecules are true positives (about 86%).For each mean photon count, 200 images of size 30 × 30 pixels were simulated using the parameters given in Section 5.1.The Hungarian algorithm with a search area of radius 100 nm is used to pair the localized molecules with the ground truth molecules.
We next examine the bias of the algorithm for a data set in which each frame contains one molecule whose location in the image is chosen randomly.For this purpose, we simulated 1000 15 × 15-pixel images, each containing one molecule with a mean photon count of 1500 photons.In each image, the location of the molecule was drawn from a uniform probability distribution that places it between the 2 nd and 14 th pixel in both the x and y dimensions.(We assumed that no molecule was located near the edges of the 15 × 15-pixel image.)As shown in Fig. 2, the deviations of both the x and y location estimates from the ground truth are, overall, centered around 0 nm.Therefore, the results suggest that, in the case where there is only one molecule per Fig. 2. Analysis of the error of location estimates obtained from a data set in which each frame contains one molecule whose location in the image is chosen randomly.Shown in the left and right plots are the differences between the x-estimates and the true x-values, and the differences between the y-estimates and the true y-values, respectively, for the true positives obtained with the algorithm.The data set consists of 1000 15 × 15-pixel images, each of a molecule with a mean photon count of 1500 photons whose location is randomly chosen from a uniform probability distribution that places the molecule between the 2 nd and 14 th pixel in both the x and y dimensions.The images were simulated using the parameters given in Section 5.1.
image, there is no systematic bias associated with the algorithm in this case (the average of x and y deviations are 0.321 nm and 0.335 nm, respectively).Also, the square root of the average of the squares of the x and y deviations are 9.123 nm and 9.467 nm, respectively, which are close to the standard deviations of the estimated locations obtained for a data set consisting of repeat images of one molecule with the same mean photon count of 1500 photons (analysis of data sets with repeat images is presented next).This suggests that the variation of the deviations about the ground truth is reasonable.
To examine further the bias of the algorithm, we simulated data sets containing repeat images of one molecule.The data sets differ by the mean photon count of the molecule, which we assume does not vary from frame to frame in a given data set.This mean photon count ranges from 500 to 4500 for the different data sets.For each data set, we simulated 1000 repeat images.Figure 3 shows, as a function of the mean photon count, the differences between the averages of the xand y-estimates for the correctly detected (i.e., true positive) molecules and the corresponding true xand y-coordinates.Similar to the case of data sets with non-repeat images [Fig.2], the evenness of the spread of the estimated bias about 0 nm for both coordinates suggests that when there is only one molecule per image, there is no systematic bias associated with our proposed algorithm.
We next evaluate the performance of the algorithm in terms of the standard deviation of the estimates for the sets of repeat images.For nine of the data sets from Fig. 3, we calculated the standard deviations of the x-estimates and y-estimates for the correctly detected (i.e., true positive) molecules.The percentage differences between the standard deviations and the CRLB-based limits of the x-localization accuracy and y-localization accuracy [17] are shown in Fig. 4. The percentage difference is the absolute difference between the standard deviation of the estimates and the corresponding limit of accuracy, expressed as a percentage of the limit of accuracy.As shown in Fig. 4, when the mean number of photons increases, the standard deviation of the estimates decreases.Also, as can be seen in the second row of Fig. 4, the differences between the standard deviations of the estimates and their respective limits of the localization accuracy are around twice (i.e., around 200% of) the limits of accuracy.This difference likely arises from the fact that our algorithm approximates an Airy profile with the frequency response of a first-order system, and there is a difference between the shape of the peak of an Airy profile and that of the The percentage difference between the standard deviation of the x-estimates and the limit of the x-localization accuracy, and the percentage difference between the standard deviation of the y-estimates and the limit of the y-localization accuracy.The percentage difference is the absolute difference between the standard deviation of the estimates and the corresponding limit of accuracy, expressed as a percentage of the limit of accuracy.first-order system in the frequency domain.In Appendix A, we applied our algorithm to images simulated using the frequency response of a first-order system rather than an Airy profile, and in that case, the standard deviations of the xand y-estimates came close to their respective limits of accuracy.Also, we used the location estimates obtained with the algorithm from a set of repeat images as initial conditions for the maximum likelihood estimation of the location of the molecule from those same images.This maximum likelihood estimator fits an Airy photon distribution profile to the image data, and the equations that describe how the maximum likelihood estimates are calculated are given in Section 4. 3 [Eqs. (38) and (39)].We calculated the standard deviations of the resulting x-estimates and y-estimates, and the percentage differences between them and the limits of the x-localization accuracy and y-localization accuracy [Fig.5].We only considered those estimates for which the estimated locations were within the image.As can be seen in Fig. 5, the standard deviations are substantially smaller compared to those obtained with the algorithm [Fig.4], and come close to the limits of accuracy, consistent with the results in [21] and [22].

Multiple molecules
So far, we have evaluated the performance of the algorithm in the case where we have only one molecule in any given image.Here, we analyze the results obtained when the algorithm was used to simultaneously localize molecules from images that contain multiple closely spaced molecules.
As before, we first analyze the detection rate of the algorithm.For this purpose, we simulated data sets containing images of either two or three molecules.For each image, the location of each molecule is randomly chosen from a uniform probability distribution that places the molecule inside the image, with the constraint that the distance between each pair of molecules is not less than a minimum distance d min .In one case, all data sets are simulated with the same minimum the location of each molecule approaches 0 nm.The results for d = s/2 are shown with filled symbols in Fig. 8. Also, we calculated the results for the y-estimates and obtained similar results.
Having characterized the nature of the bias, we next analyze, as we did in the case of one molecule, the bias and accuracy of the algorithm as a function of the mean photon count per molecule.For this purpose, we simulated data sets which contain repeat images of two molecules.These data sets again differ by the mean photon count per molecule, which we assume does not vary from frame to frame.This mean photon count ranges from 500 to 4500 for the different data sets.We simulated 500 20 × 20-pixel images per data set.In one case, the distance d between the two molecules is 650 nm (which corresponds to the filled circle in Fig. 8), and in another case, d = 487.5 nm (which corresponds to the first open circle to the left of the filled circle in Fig. 8).Looking at the two distances allows us to verify the effect of different distances between the molecules relative to the image size.To assess the bias of the algorithm, for each molecule in a given data set we calculated the difference between the average of the estimated x-locations and the corresponding true x-coordinate.As can be seen in Fig. 9, when d = 650 nm, the estimated bias is around 0 nm for both molecules.On the other hand, when d = 487.5 nm, the estimated bias levels are around 3.5 nm and -3.5 nm for molecules 1 and 2, respectively.These bias results are consistent with the illustration of bias in Fig. 8.Note that we also analyzed the y-estimates and obtained similar results.For each distance d, we calculated the standard deviations of the estimated x-locations for nine of the data sets.As shown in Fig. 9, for both distances, as the mean number of photons per molecule increases, the standard deviation of the estimates decreases.Also, even when the mean photon count is as low as 500 photons/molecule, the plots show that our algorithm can still localize the molecules with relatively high accuracy (the standard deviations of the x-estimates are around 30 nm for both molecules when d = 487.5 nm, and around 27 nm when d = 650 nm).Similar results were obtained for the y-estimates.Results are shown for three simulated data sets that differ by the mean photon count per molecule per image.Each data set consists of 100 images in which there are two molecules per image.The location of each molecule is randomly chosen from a uniform distribution that places the molecule inside the image, and is subject to the constraint that the distance between the two molecules is not less than 400 nm.show that in the case that we have a relatively large number of photons per molecule, the detection rate of the algorithm is not very sensitive to the threshold values.Provided that a relatively high threshold value (e.g., in the range [0.7, 0.9]) is used to ensure that singular values corresponding to signal are not discarded, it appears that the optimization procedure of Section 4.2 is able to remove the singular values corresponding to noise and yield the correct result.On the other hand, in the case of a low number of photons per molecule, the differences between the singular values that correspond to noise and the singular values associated with signal are often small, and there is no straightforward guideline to choose the threshold values.
In our analysis, we consider three simulated data sets in which there are two molecules per image.The three data sets differ by the mean photon count per molecule per image, which we chose to be 500, 1000, and 2500.As can be seen in Table 1, when the mean photon count is 1000 or 2500 per molecule, the recall and precision remain unchanged for threshold values of 0.7, 0.8, and 0.9.However, in the case of the low mean photon count of 500 per molecule, use of a high threshold value, such as 0.8 or 0.9, for the first SVD, leads to the retention of more noise singular values and results in a nontrivial reduction of the precision.

4. 2 .
Determination of the number of retained singular values in the second and third SVDs, and the number of single molecules in the image

Fig. 1 .
Fig.1.Analysis of the detection rate of the algorithm, applied to data sets in which each image contains one molecule, whose location in the image is chosen randomly according to a uniform probability distribution.For a given data set, the same mean photon count is used to simulate the molecule in each image.Different data sets differ by this mean photon count.For each mean photon count, 200 images of size 30 × 30 pixels were simulated using the parameters given in Section 5.1.The Hungarian algorithm with a search area of radius 100 nm is used to pair the localized molecules with the ground truth molecules.

Fig. 3 .Fig. 4 .
Fig.3.Analysis of the average of location estimates obtained from repeat images of one molecule.Shown in the left and right plots are the difference between the average of the x-estimates and the true x-value, and the difference between the average of the y-estimates and the true y-value, respectively, for data sets that differ by the mean photon count assumed for the molecule per image.For each mean photon count, the data set consists of 1000 images of size 15 × 15 pixels, simulated using the parameters given in Section 5.1.

Fig. 5 .
Fig.5.Analysis of the standard deviation of location estimates produced by the maximum likelihood estimator when the location estimates obtained with the algorithm are used as the initial conditions.(a) The standard deviations of the maximum likelihood xand y-estimates for the same data sets as in Fig.4, which comprise repeat images of one molecule.(b) The percentage difference between the standard deviation of the x-estimates and the limit of the x-localization accuracy, and the percentage difference between the standard deviation of the y-estimates and the limit of the y-localization accuracy.

Fig. 8 .Fig. 9 .
Fig.8.Analysis of the average of the location estimates obtained from sets of repeat images of two molecules.Shown in the left and right plots are the differences between the average of the x-estimates and the true x-value for the first and second molecules, respectively, for data sets comprising 15 × 15-pixel, 20 × 20-pixel, and 40 × 40-pixel images.For each image size, distances d between the two molecules are chosen around half of the side length of the square region occupied by the image in the object space.For a given data set, we simulated 500 images with a mean photon count of 2500 photons/molecule and the parameters given in Section 5.1.The results for d = s/2, where s = 65N nm is the side length of the square region occupied by an N × N-pixel image in the object space, are shown with filled symbols.

Fig. 10 .
Fig. 10.Result of the algorithm applied to an experimental super-resolution image.(a) Image of individual Alexa Fluor 647 molecules acquired using the microscopy setup described in Section 5.2.The pixel size and image size are 16 µm × 16 µm and 192 × 192 pixels, respectively.(b) The magnitude of the reconstructed image obtained with the algorithm.

6. 1 . 3 .
Analysis of the effect of threshold values on the detection rate of the algorithm In Sections 4.1 and 4.2, typical threshold values in the range [0.8, 0.9] are suggested for the first, second, and third SVDs of the algorithm.Here, we carry out a more in-depth analysis on the effect of the threshold values on the detection rate of the algorithm.The results of our analysis

Fig. 11 .
Fig. 11.Results of the algorithm applied to an ROI from an experimental super-resolution image.(a) A 41 × 41-pixel ROI of the super-resolution image shown in Fig. 10.(b) The magnitude of the reconstructed image (algorithm result).(c) The image reconstructed using Eq.(36), in which the single molecule locations estimated using our algorithm are used in the computation of the Airy profile q in Eq. (37), and in which the mean photon counts N p,n and the parameter α := 2πn a λ are separately estimated with a maximum likelihood estimator.(d), (e), and (f) show the mesh plots of the images in (a), (b), and (c), respectively.

Table 1 .
Detection rate of the algorithm as a function of the threshold values used for the retention of singular values in the first, second and third SVDs.