## Abstract

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

## 1. 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 continuous-space 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–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.

## 2. 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*) ∈ ℂ^{p}^{×}* ^{m}*,

*p*,

*m*∈ ℕ,

*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.

**Lemma 1**. *Let*
$\tilde{X}\left(k\right)\in \mathbb{R}$, *k* = 1, 2, …, *N*, *be a finite 1D sequence. For n* = 1, 2, …, *N, let*
$X\left(n\right):=\left(IDFT\left(\tilde{X}\right)\right)\left(n\right)=\frac{1}{N}{\displaystyle {\sum}_{k=1}^{N}\tilde{X}\left(k\right){e}^{i2\pi kn/N}}$ *be the inverse discrete Fourier transform (inverse DFT, or IDFT) of*
$\tilde{X}$. *Then, there exists a minimal and asymptotically stable system* (*A*, *B*, *C*), *such that*

*Moreover*,

*where*$\tilde{A}:=A$, $\tilde{B}:=\left(I-{A}^{N}\right)B$, $\tilde{C}=C$.

*If A*= 0,

^{N}*then*$\left(\tilde{A},\tilde{B},\tilde{C}\right)=\left(A,B,C\right)$.

*Proof.* Let
$\tilde{X}\left(k\right)\in \mathbb{R}$, *k* = 1, 2, …, *N*, be a finite 1D sequence. Let

*A*,

*B*,

*C*), such that

According to Eqs. (4) and (5), we then have, for *k* = 1, 2, …, *N*,

For a square matrix *T* ∈ ℂ^{m}^{×}* ^{m}*,

*m*∈ ℕ, where the number 1 is not an eigenvalue of

*T*, we have the identity ${\sum}_{n=0}^{N-1}{T}^{n}={\left(I-T\right)}^{-1}\left(I-{T}^{N}\right)$. 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*

^{−}

^{i}^{2}

^{πk}^{/}

*,*

^{N}*k*= 1, …,

*N*(or equivalently,

*I*−

*Ae*

^{−}

^{i}^{2}

^{πk}^{/}

*,*

^{N}*k*= 1, …,

*N*, is invertible), and ${\sum}_{n=0}^{N-1}{\left(A{e}^{-i2\pi k/N}\right)}^{n}={\left(I-A{e}^{-i2\pi k/N}\right)}^{-1}\left(I-{A}^{N}\right)$. Substituting this expression into Eq. (6), we have, for

*k*= 1, 2, …,

*N*,

*A*= 0, then $\left(\tilde{A},\tilde{B},\tilde{C}\right)=\left(A,B,C\right)$. □

^{N}In the following theorem, we extend the results obtained for 1D sequences to 2D sequences.

**Theorem 1.** *Let*
$\tilde{X}\left({k}_{1},{k}_{2}\right)\in \mathbb{R}$, *k*_{i} = 1, 2, …, *N _{i}*,

*i*= 1, 2,

*be a finite 2D sequence. For n*= 1, 2, …,

_{i}*N*,

_{i}*i*= 1, 2,

*let*

*be the inverse 2D DFT of*$\tilde{X}$.

*Then, there exist minimal and asymptotically stable systems*(

*A*,

_{i}*B*,

_{i}*C*),

_{i}*i*= 1, 2,

*such that*

*where, for i*= 1, 2,

*Moreover*,

*where, for k*= 1, 2, …,

_{j}*N*,

_{j}*j*= 1, 2,

*where*${\tilde{A}}_{j}:={A}_{j}$, ${\tilde{B}}_{j}:=\left(I-{A}_{j}^{{N}_{j}}\right){B}_{j}$, ${\tilde{C}}_{j}:={C}_{j}$.

*For j*= 1, 2,

*if*${A}_{j}^{{N}_{j}}=0$,

*then*$\left({\tilde{A}}_{j},{\tilde{B}}_{j},{\tilde{C}}_{j}\right)=\left({A}_{j},{B}_{j},{C}_{j}\right)$.

*Proof.* Let
$\tilde{X}\left({k}_{1},{k}_{2}\right)\in \mathbb{R}$, *k _{i}* = 1, 2, …,

*N*,

_{i}*i*= 1, 2, be a finite 2D sequence. For

*n*= 1, …,

_{i}*N*,

_{i}*i*= 1, 2, let

*X*to form a matrix

*Q*as

Decompose *Q* via SVD as *Q* = *U*Σ*V*, where for *r* ∈ ℕ,
$U\in {\u2102}^{{N}_{1}\times r}$, Σ ∈ ℂ^{r}^{×}* ^{r}* and
$V\in {\u2102}^{r\times {N}_{2}}$.

For *n _{i}* = 1, 2, …,

*N*,

_{i}*i*= 1, 2, define

*X*

_{1}(

*n*

_{1}) ∈ ℂ

^{1×}

*and*

^{r}*X*

_{2}(

*n*

_{2}) ∈ ℂ

^{r}^{×1}, such that

Then

Moreover, according to Proposition 1, there exist minimal and asymptotically stable systems (*A _{i}*,

*B*,

_{i}*C*),

_{i}*i*= 1, 2, such that, for

*i*= 1, 2,

According to Eqs. (13) and (16),

*k*= 1, 2, …,

_{i}*N*,

_{i}*i*= 1, 2. Then, according to Lemma 1, for

*k*= 1, 2,, …,

_{j}*N*,

_{j}*j*= 1, 2,

*j*= 1, 2, if ${A}_{j}^{{N}_{j}}=0$, then $({\tilde{A}}_{j},{\tilde{B}}_{j},{\tilde{C}}_{j})=({A}_{j},{B}_{j},{C}_{j})$.

## 3. Location estimation

So far, we have shown the existence of minimal and asymptotically stable systems (*A _{i}*,

*B*,

_{i}*C*),

_{i}*i*= 1, 2, that realize a finite 2D sequence $\tilde{X}\left({k}_{1},{k}_{2}\right)\in \mathbb{R}$,

*k*= 1, 2, …,

_{i}*N*,

_{i}*i*= 1, 2, in the frequency domain. Here, we summarize our overall localization approach. Given that $\tilde{X}$ is a single molecule image with multiple peaks of intensity, we determine the locations of the molecules by calculating the pole locations of (

*A*,

_{i}*B*,

_{i}*C*),

_{i}*i*= 1, 2.

In Theorem 1, we have shown that there exist minimal and asymptotically stable systems (*A _{i}*,

*B*,

_{i}*C*),

_{i}*i*= 1, 2, such that

*k*= 1, 2, …,

_{j}*N*,

_{j}*j*= 1, 2,

*j*= 1, 2, if ${A}_{j}^{{N}_{j}}=0$, then $({\tilde{A}}_{j},{\tilde{B}}_{j},{\tilde{C}}_{j})=({A}_{j},{B}_{j},{C}_{j})$. If we diagonalize

*A*,

_{i}*i*= 1, 2, then the diagonal elements of the resulting matrix ${\overline{A}}_{i}$ give the poles of the system. In the following, we use matrix diagonalization in Eq. (20) to express $\tilde{X}$ in terms of the poles of the system.

For *s*_{1}, *s*_{2} ∈ ℕ and
${A}_{i}\in {\u2102}^{{s}_{i}\times {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}\in {\u2102}^{{s}_{i}\times {s}_{i}}$, we have the diagonal matrix ${\overline{A}}_{i}:={T}_{i}{A}_{i}{T}_{i}^{-1}=\text{diag}({a}_{1}^{i},\cdots ,{a}_{{s}_{i}}^{i})$, ${a}_{{t}_{i}}^{i}\in \u2102$, then with ${\overline{B}}_{i}:={T}_{i}{B}_{i}={\left[{b}_{1}^{i},\cdots ,{b}_{{s}_{i}}^{i}\right]}^{T}$, ${\overline{C}}_{i}:={C}_{i}{T}_{i}^{-1}=\left[{c}_{1}^{i},\cdots ,{c}_{{s}_{i}}^{i}\right]$,

*i*= 1, 2, where ${b}_{{t}_{1}}^{1}\in {\u2102}^{1\times r}$, ${b}_{{t}_{2}}^{2}\in \u2102$, ${c}_{{t}_{1}}^{1}\in \u2102$, ${c}_{{t}_{2}}^{2}\in {\u2102}^{r\times 1}$,

*t*= 1, 2, …,

_{i}*s*= 1, 2, for

_{i}i*k*= 1. 2, …,

_{j}*N*,

_{j}*j*= 1, 2 we can write $\tilde{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
$\tilde{X}$ occur at
$({a}_{{t}_{1}}^{1},{a}_{{t}_{2}}^{2})$, *t _{i}* = 1, …,

*s*,

_{i}*i*= 1, 2. (Note that for

*t*= 1, …,

_{i}*s*,

_{i}*i*= 1, 2, if the row ${b}_{{t}_{1}}^{1}$ of ${\overline{B}}_{1}$ and the column ${c}_{{t}_{2}}^{2}$ of ${\overline{C}}_{2}$ are perpendicular, i.e., ${b}_{{t}_{1}}^{1}{c}_{{t}_{2}}^{2}=0$, then the residual of the pole $({a}_{{t}_{1}}^{1},{a}_{{t}_{2}}^{2})$ is equal to zero, and $({a}_{{t}_{1}}^{1},{a}_{{t}_{2}}^{2})$ 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
$\tilde{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
${a}_{{t}_{j}}^{j}=\left|{a}_{{t}_{j}}^{j}\right|{e}^{i{w}_{{t}_{j}}^{j}}$,
$0\le {w}_{{t}_{j}}^{j}\le 2\pi $, *t _{j}* = 1, …,

*s*,

_{j}*j*= 1, 2. 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 $\left\{\left({x}_{{t}_{2}},{y}_{{t}_{1}}\right):{c}_{{t}_{1}}^{1}{b}_{{t}_{1}}^{1}{c}_{{t}_{2}}^{2}{b}_{{t}_{2}}^{2}\ne 0,{t}_{i}=1,\dots ,{s}_{i},i=1,2\right\}$ where

*M*> 0 denotes the lateral magnification of the microscope system.

## 4. Algorithm

We now explain our proposed approach in more detail. In Section 3, we have calculated the poles of a 2D single molecule image
$\tilde{X}$ in terms of the elements of minimal and asymptotically stable systems (*A _{i}*,

*B*,

_{i}*C*),

_{i}*i*= 1, 2. Here, using the balanced state space realization algorithm introduced by Maciejowski [19], we propose a step-by-step algorithm to calculate systems (

*A*,

_{i}*B*,

_{i}*C*),

_{i}*i*= 1, 2, that realize $\tilde{X}$, and to determine the locations of the single molecules using the realization.

**Algorithm.** *Let*
$\tilde{X}\left({k}_{1},{k}_{2}\right)\in \mathbb{R}$, *k _{i}* = 1, 2, …,

*N*,

_{i}*i*= 1, 2,

*represent the acquired image data.*

*Subtract an estimated background level*$\widehat{\beta}$,*e.g., the average of the data points near the boundary of the image data*$\tilde{X}$,*from the image data*$\tilde{X}$,*and define the background-subtracted image*${\tilde{X}}_{bs}$*as**Arrange the entries of X to form a matrix Q as*$$Q:=\left[\begin{array}{cccc}X\left(1,1\right)& X\left(1,2\right)& \cdots & X\left(1,{N}_{2}\right)\\ X\left(2,1\right)& X\left(2,2\right)& \cdots & X\left(2,{N}_{2}\right)\\ \vdots & \vdots & \ddots & \vdots \\ X\left({N}_{1},1\right)& X\left({N}_{1},2\right)& \cdots & X\left({N}_{1},{N}_{2}\right)\end{array}\right].$$*Decompose Q via SVD as Q*=*U*Σ*V. Let the positive integer r*≤*K*,*K*= min(*N*_{1},*N*_{2}),*denote the number of retained singular values (see Section 4.1). Partition*$\mathrm{\Sigma}=diag(\widehat{\mathrm{\Sigma}},\widehat{\widehat{\mathrm{\Sigma}}})$, $\widehat{\mathrm{\Sigma}}\in {\u2102}^{r\times r}$, $U=\left[\widehat{U}\phantom{\rule{1em}{0ex}}\widehat{\widehat{U}}\right]$, $\widehat{U}\in {\u2102}^{{N}_{1}\times r}$,*and*$V={\left[\widehat{V}\phantom{\rule{1em}{0ex}}\widehat{\widehat{V}}\right]}^{T}$, $\widehat{V}\in {\u2102}^{r\times {N}_{2}}$.*For n*= 1, 2, …,_{i}*N*,_{i}*i*= 1, 2,*define*${X}_{1}^{r}\left({n}_{1}\right)\in {\u2102}^{1\times r}$*and*${X}_{2}^{r}\left({n}_{2}\right)\in {\u2102}^{r\times 1}$,*such that*$$\left[\begin{array}{c}{X}_{1}^{r}\left(1\right)\\ {X}_{1}^{r}\left(2\right)\\ \vdots \\ {X}_{1}^{r}\left({N}_{1}\right)\end{array}\right]:=\widehat{U}{\widehat{\mathrm{\Sigma}}}^{1/2},\phantom{\rule{1em}{0ex}}\left[{X}_{2}^{r}\left(1\right)\phantom{\rule{1em}{0ex}}{X}_{2}^{r}\left(2\right)\phantom{\rule{1em}{0ex}}\cdots \phantom{\rule{1em}{0ex}}{X}_{2}^{r}\left({N}_{2}\right)\right]:{\widehat{\mathrm{\Sigma}}}^{1/2}\widehat{V}.$$*Construct the Hankel matrices*${H}_{1}\in {\u2102}^{\left({N}_{1}+1\right)\times \left({N}_{1}+1\right)r}$, ${H}_{2}\in {\u2102}^{\left({N}_{2}+1\right)\times \left({N}_{2}+1\right)}$*as*$${H}_{i}:=\left[\begin{array}{cccccc}{X}_{i}^{r}\left(1\right)& {X}_{i}^{r}\left(2\right)& \cdots & {X}_{i}^{r}\left({N}_{i}-1\right)& {X}_{i}^{r}\left({N}_{i}\right)& 0\\ {X}_{i}^{r}\left(2\right)& {X}_{i}^{r}\left(3\right)& \cdots & {X}_{i}^{r}\left({N}_{i}\right)& 0& 0\\ \vdots & \vdots & \ddots & \vdots & \vdots & \vdots \\ {X}_{i}^{r}\left({N}_{i}\right)& 0& \cdots & 0& 0& 0\\ 0& 0& \cdots & 0& 0& 0\end{array}\right],\phantom{\rule{1em}{0ex}}i=1,2,$$*where*0*denotes a block of zeros of the corresponding size. For i*= 1, 2,*decompose H*=_{i}via SVD as H_{i}*U*Σ_{i}≤_{i}V_{i}. Let the positive integers s_{i}*N*,_{i}*i*= 1, 2,*denote the numbers of retained singular values in the respective SVDs (see Section 4.2). For i*= 1, 2,*partition*${\mathrm{\Sigma}}_{i}=diag({\widehat{\mathrm{\Sigma}}}_{i},{\widehat{\widehat{\mathrm{\Sigma}}}}_{i})$, ${\widehat{\mathrm{\Sigma}}}_{i}\in {\u2102}^{{s}_{i}\times {s}_{i}}$, ${U}_{i}=\left[{\widehat{U}}_{i}\phantom{\rule{1em}{0ex}}{\widehat{\widehat{U}}}_{i}\right]$, ${\widehat{U}}_{1}\in {\u2102}^{\left({N}_{1}+1\right)\times {s}_{1}}$, ${\widehat{U}}_{2}\in {\u2102}^{\left({N}_{2}+1\right)r\times {s}_{2}}$,*and*${V}_{i}={\left[{\widehat{V}}_{i}\phantom{\rule{1em}{0ex}}{\widehat{\widehat{V}}}_{i}\right]}^{T}$, ${\widehat{V}}_{1}\in {\u2102}^{{s}_{1}\times \left({N}_{1}+1\right)r}$, ${\widehat{V}}_{2}\in {\u2102}^{{s}_{2}\times \left({N}_{2}+1\right)}$,*conformally*.*Let*${C}_{1}^{r;{s}_{1}}\in {\u2102}^{1\times {s}_{1}}$*and*${C}_{2}^{r;{s}_{2}}\in {\u2102}^{1\times {s}_{2}}$*be the first row of*${\widehat{U}}_{1}{\widehat{\mathrm{\Sigma}}}_{1}^{1/2}$*and the first r of*${\widehat{U}}_{2}{\widehat{\mathrm{\Sigma}}}_{2}^{1/2}$,*respectively. Also, let*${B}_{1}^{r;{s}_{1}}\in {\u2102}^{{s}_{1}\times r}$*and*${B}_{2}^{r;{s}_{2}}\in {\u2102}^{{s}_{2}\times 1}$*be the first r columns of*${\widehat{\mathrm{\Sigma}}}_{1}^{1/2}{\widehat{V}}_{1}$,*respectively. Assuming*$${\widehat{U}}_{i}=\left[\begin{array}{c}{\overline{U}}_{1}^{i}\\ \vdots \\ {\overline{U}}_{{N}_{i}}^{i}\\ {\overline{U}}_{{N}_{i}+1}^{i}\end{array}\right],\phantom{\rule{1em}{0ex}}i=1,2,$$*where*${\overline{U}}_{{n}_{1}}^{1}\in {\u2102}^{1\times {s}_{1}}$, ${\overline{U}}_{{n}_{2}}^{2}\in {\u2102}^{r\times {s}_{2}}$,*n*= 1, …,_{i}*N*+ 1,_{i}*i*= 1, 2,*define*$${\widehat{U}}_{i}^{\uparrow}:=\left[\begin{array}{c}{\overline{U}}_{2}^{i}\\ \vdots \\ {\overline{U}}_{{N}_{i}+1}^{i}\end{array}\right],\phantom{\rule{1em}{0ex}}{\widehat{U}}_{i}^{\downarrow}:=\left[\begin{array}{c}{\overline{U}}_{1}^{i}\\ \vdots \\ {\overline{U}}_{{N}_{i}}^{i}\end{array}\right],\phantom{\rule{1em}{0ex}}i=1,2.$$Then, let ${A}_{i}^{r;{s}_{i}}={\widehat{\mathrm{\Sigma}}}_{i}^{-1/2}{\widehat{U}}_{i}^{\downarrow *}{\widehat{U}}_{i}^{\uparrow}{\widehat{\mathrm{\Sigma}}}_{i}^{1/2}\in {\u2102}^{{s}_{i}\times {s}_{i}}$, i = 1, 2.*V. Diagonalize*${A}_{j}^{r;{s}_{j}}\in {\u2102}^{{s}_{j}\times {s}_{j}}$,*j*= 1; 2,*i.e., for t*= 1; 2; …;_{j}*s*;_{j}*j*= 1; 2,*and some invertible*${T}_{j}\in {\u2102}^{{s}_{j}\times {s}_{j}}$,*let*${\overline{A}}_{j}^{r;{s}_{j}}:={T}_{j}{A}_{j}^{r;{s}_{j}}{T}_{j}^{-1}=diag\left({a}_{1}^{j},\cdots ,{a}_{{s}_{j}}^{j}\right)$, ${a}_{{t}_{j}}^{j}=\left|{a}_{{t}_{j}}^{j}\right|{e}^{i{w}_{{t}_{j}}^{j}}\in \u2102$, $0\le {w}_{{t}_{j}}^{j}\le 2\pi $,*be a corresponding diagonal matrix for*${A}_{j}^{r;{s}_{j}}$.*Also, let*${\overline{B}}_{j}^{r;{s}_{j}}:={T}_{j}{B}_{j}^{r;{s}_{j}}={\left[{b}_{1}^{j},\cdots ,{b}_{{s}_{j}}^{j}\right]}^{T}$, ${\overline{C}}_{j}^{r;{s}_{j}}:={C}_{j}^{r;{s}_{j}}{T}_{j}^{-1}=\left[{c}_{1}^{j},\cdots ,{c}_{{s}_{j}}^{j}\right]$,*where*${b}_{{t}_{1}}^{1}\in {\u2102}^{1\times r}$, ${b}_{{t}_{2}}^{2}\in \u2102$, ${c}_{{t}_{1}}^{1}\in \u2102$, ${c}_{{t}_{2}}^{2}\in {\u2102}^{r\times 1}$,*t*= 1; 2; …;_{j}*s*;_{j}*j*= 1; 2.*Note that in theory, there is a possibility that*${A}_{1}^{r;{s}_{1}}$*and/or*${A}_{2}^{r;{s}_{2}}$*are not diagonalizable. In practice, however, because the diagonalization is numerically computed*, ${A}_{1}^{r;{s}_{1}}$*and/or*${A}_{2}^{r;{s}_{2}}$*can be expected to be diagonalizable. In the unlikely scenario where*${A}_{1}^{r;{s}_{1}}$*and/or*${A}_{2}^{r;{s}_{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.**For h*= min(*s*_{1};*s*_{2}),*calculate, in the object space, the estimated peak locations*(*x*;_{k}*y*); ${x}_{k}\in \left\{{\widehat{x}}_{1},\dots ,{\widehat{x}}_{{s}_{2}}\right\}$, ${y}_{k}\in \left\{{\widehat{y}}_{1},\dots ,{\widehat{y}}_{{s}_{1}}\right\}$,_{k}*k*= 1, …,*h*, where$${\widehat{x}}_{{t}_{2}}:=\frac{\mathrm{\Delta}x{w}_{{t}_{2}}^{2}{N}_{1}}{2M\pi}+\frac{\mathrm{\Delta}x}{2M},\phantom{\rule{1em}{0ex}}{\widehat{y}}_{{t}_{1}}:=\frac{\mathrm{\Delta}y{w}_{{t}_{1}}^{1}{N}_{2}}{2M\pi}+\frac{\mathrm{\Delta}y}{2M},\phantom{\rule{1em}{0ex}}{t}_{i}=1,2,\dots ;{s}_{i},\phantom{\rule{1em}{0ex}}i=1,2,$$*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.

#### 4.1. 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}:={\displaystyle {\sum}_{i=1}^{r}{\sigma}_{i}^{2}}$ be the energy of the sequence

*σ*,

_{i}*i*= 1, …,

*r*. Estimate the optimal number of retained singular values

*r*in the first SVD as

#### 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

Let
${\sigma}_{1}^{i}\ge \dots \ge {\sigma}_{{N}_{i}}^{i}\ge 0$, *i* = 1, 2, be the singular values in the second and third SVDs, respectively. For *l _{i}* = 1, …,

*N*,

_{i}*i*= 1, 2, let

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
${s}_{i}=1,\dots ,{\widehat{l}}_{i},i=1,2$, let
${\tilde{X}}^{r;{s}_{1},{s}_{2}}\left({k}_{1},{k}_{2}\right)={\displaystyle {\sum}_{l=1}^{{s}_{1}}{\displaystyle {\sum}_{j=1}^{{s}_{2}}\frac{{c}_{l}^{1}{b}_{l}^{1}{c}_{j}^{2}{b}_{j}^{2}}{\left({e}^{i2\pi {k}_{1}/{N}_{1}}-{a}_{l}^{1}\right)\left({e}^{i2\pi {k}_{2}/{N}_{2}}-{a}_{j}^{2}\right)}}}$, *k _{t}* = 1, …,

*N*,

_{t}*t*= 1, 2, 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, ${\tilde{X}}^{r;{s}_{1},{s}_{2}}$ is the reconstructed image of Eq. (22) after discarding the singular values corresponding to noise. Denoting the poles of ${\tilde{X}}^{r;{s}_{1},{s}_{2}}$ by $\left({\overline{a}}_{k}^{1},{\overline{a}}_{k}^{2}\right)$, ${\overline{a}}_{k}^{t}\in \left\{{a}_{1}^{t},\cdots ,{a}_{{s}_{t}}^{t}\right\}$, ${\overline{a}}_{k}^{t}:=\left|{\overline{a}}_{k}^{t}\right|{e}^{i{\overline{w}}_{k}^{t}}$, $0\le {\overline{w}}_{k}^{t}\le 2\pi $,

*k*= 1, …,

*s*

_{1}

*s*

_{2},

*t*= 1, 2, and their corresponding product of coefficients in the numerator by.

*p*∈ ℂ, k = 1, …,

_{k}*s*

_{1}

*s*

_{2}, assume the peak magnitudes to be $\left|{p}_{1}\right|\ge \cdots \ge \left|{p}_{{s}_{1}{s}_{2}}\right|\ge 0$.

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
${\widehat{\theta}}^{h}:=\left({\widehat{\theta}}_{1},\dots ,{\widehat{\theta}}_{h}\right)\in {\mathbb{R}}^{2h}$,
${\widehat{\theta}}_{n}:=\left({\widehat{x}}_{n},{\widehat{y}}_{n}\right)\in {\mathbb{R}}^{2}$, *n* = 1, …, *h*, such that

*h*peaks with the largest magnitudes. In general, we consider all possible

*h*-combinations of the poles of ${\tilde{X}}^{r;{s}_{1},{s}_{2}}$, but in most cases, the single molecules are associated with the peaks with the largest magnitudes. Let $\left\{{z}_{1},\dots ,{z}_{{N}_{pix}}\right\}$ denote our acquired data, where

*N*denotes the number of pixels in the image. Then, the estimated number of single molecules $\widehat{h}$ is given by

_{pix}*k*= 1, …,

*N*, the mean number of photons detected in the

_{pix}*k*

^{th}pixel given by [17]

*N*is the expected number of photons due to the

_{p,n}*n*molecule that impact the detector plane during the image exposure,

^{th}*C*⊂ ℝ

_{k}^{2}denotes the region in the detector plane occupied by the

*k*pixel, and

^{th}*q*is the 2D PSF of the optical system. If the PSF is the Airy profile, then

*q*is given by

*n*denotes the numerical aperture of the objective lens,

_{a}*λ*denotes the emission wavelength of the molecule, and

*J*

_{1}denotes the first order Bessel function of the first kind.

#### 4.3. 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 ℝ* ^{n}*. The maximum likelihood estimate
${\widehat{\theta}}_{mle}$ of

*θ*∈ Θ, for data incorporating Gaussian readout noise, is given by

*N*pixels and

_{pix}*L*is the log-likelihood function given by [17]

In Eq. (39), in the case of one molecule, *µ _{θ}*(

*k*) is the mean photon count in the

*k*pixel due to the molecule and is given by

^{th}*θ*= (

*x*

_{0},

*y*

_{0}) ∈ ℝ

^{2}denotes the location of the molecule in the object space,

*N*is the expected number of photons from the molecule that are detected over the detector plane, and

_{p}*q*is the Airy profile given by Eq. (37). Also,

*β*is the background level in the

_{k}*k*pixel, and

^{th}*η*and

_{k}*σ*denote the mean and standard deviation of the Gaussian readout noise in the

_{k}*k*pixel, respectively.

^{th}## 5. Methods

#### 5.1. 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.

#### 5.2. Imaging experiments

### 5.2.1. 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.

### 5.2.2. 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.

### 5.2.3. 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.

## 6. 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.

#### 6.1. 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.

#### 6.1.1. 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 *TP*, the number of false negatives by *FN*, and the number of false positives by *FP*, we define the precision (*PRE*) and recall (*REC*) measures as [20]

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%).

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

*pixel in both the*

^{th}*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 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 *x*- and *y*-estimates for the correctly detected (i.e., true positive) molecules and the corresponding true *x*-and *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 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 *x*- and *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].

### 6.1.2. 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 distance of

*d*= 100 nm, but differ by the mean photon count per molecule, which ranges from 500 to 4500. More specifically, each data set comprises 200 images, where each image contains molecules with the same given mean photon count. In another case, the mean photon count is the same for all data sets at 2500 photons/molecule, but the data sets differ by

_{min}*d*, which ranges from 100 nm to 500 nm. Here, each data set comprises 200 images, where each image contains molecules separated by the same given minimum distance

_{min}*d*. Figure 6 shows the precision and recall measures for the different data sets (similar to the one-molecule case, the Hungarian algorithm with a search area of radius 100 nm was used to pair the localized molecules with the ground truth molecules). Specifically, the figure shows that the recall for all data sets is more than 95%. The precision is likewise quite good, as even when the mean number of photons is as low as 500 photons/molecule, or the minimum distance between each pair of molecules is as small as 100 nm, a relatively high percentage (around 85%) of the detected molecules are true positives. We also analyzed the detection rate for data sets with more molecules per image (5, 7, and 9 molecules per image), and obtained similar results. To give examples of the reconstructed image calculated from our algorithm, we simulated images of size 60 × 60 pixels containing two or three closely spaced molecules separated by a distance

_{min}*d*of 50 nm, 250 nm, and 300 nm, with a mean photon count of 1500 photons/molecule. We then reconstructed each image by applying our algorithm. As shown in Fig. 7, in all cases we were able to distinguish the closely spaced molecules from each other.

We next analyze the bias and accuracy of the algorithm when applied to data sets consisting of repeat images of multiple molecules. Unlike the one-molecule case, bias is observed here. We first characterize this bias by demonstrating its dependence on the distances between the molecules relative to the image size. In the following, we focus on data sets comprising images of two molecules, though we also analyzed data sets with three and five molecules and obtained similar results. We simulated data sets comprising 15 × 15-pixel, 20 × 20-pixel, and 40 × 40-pixel images in order to evaluate different combinations of the distance *d* between the two molecules and the image size. For a given data set, we simulated 500 images with a mean photon count of 2500 photons/molecule, a pixel size of 6.5 *µ*m × 6.5 *µ*m, and a lateral magnification of 100. For these settings, the area occupied by an *N* × *N*-pixel image in the object space is an *s* × *s* square region, where *s* = 65*N* nm (e.g., for N=20, *s* = 1300 nm). For each data set, the difference between the average of the estimated *x*-locations for the correctly detected molecules and the corresponding true *x*-coordinate is plotted in Fig. 8. For each image size considered, the figure shows that as *d* approaches *s*/2 (e.g., for N=20, *s*/2 = 650 nm), i.e., as the difference between the phases of the poles of the second-order system resulting from the algorithm approaches the maximum of *π* rad on the unit circle, the effect of the poles on each other decreases and the estimated bias for 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.

### 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 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.

#### 6.2. Results for experimental data

Here, we present the results obtained by applying the proposed localization algorithm to experimental single molecule image data acquired as described in Section 5.2. Both the acquired image and the reconstructed image are shown in Fig. 10, demonstrating that we were able to recover the locations of the significant peaks in the original image that are associated with the locations of individual Alexa Fluor 647 dye molecules.

We next applied the algorithm to a relatively small 41 × 41-pixel region of interest (ROI) [Figs. 11(a) and 11(d)] in the acquired image so that a better visual comparison can be made between the reconstructed image obtained with the algorithm and the actual image. In addition to the image reconstructed in terms of the parameters of the multi-order system [Figs. 11(b) and 11(e)], an image is reconstructed using Eq. (36), in which the single molecule locations estimated using our algorithm are used in the computation of the Airy profile *q* [Figs. 11(c) and 11(f)]. For this purpose, we separately estimated the mean photon counts *N _{p,n}* in Eq. (36) and the parameter
$\alpha :=\frac{2\pi {n}_{a}}{\lambda}$ in Eq. (37) using a maximum likelihood estimator. The reconstruction using the Airy profile provides a better visual comparison with the actual image by showing peaks with more comparable intensities.

## Appendix A: Frequency response of a multi-order system as the PSF model

In this section, we present the results of the proposed algorithm applied to images simulated using the frequency response of a multi-order system. In this case, instead of Eq. (36), ${\mu}_{{\widehat{\theta}}^{h}}$ is given as

*N*

_{p}_{;1}for the molecule. For each mean photon count, the data set consists of 1000 repeat images of size 20 × 20 pixels. In Figs. 12(a) and 12(b), an example of an image with a mean photon count of

*N*

_{p}_{;1}= 1000 is shown. To assess the bias of the algorithm, we calculated the differences between the averages of the

*x*- and

*y*-estimates and the corresponding true x- and y-coordinates. Similar to the case of images simulated with the Airy profile [Fig. 3], the evenness of the spread of the estimated bias about 0 nm for both coordinates [Fig. 12(c)] demonstrates that there is no systematic bias associated with our proposed algorithm when there is only one molecule per image.

Also, we calculated the standard deviation of the estimates for nine of the data sets and compared the results with the limit of the localization accuracy calculated using the approach for experimental PSFs presented in [23]. It can be seen in Fig. 13 that when the image of the molecule is simulated as the frequency response of a first-order system, the accuracy of the algorithm comes close to the limit of the localization accuracy.

## Funding

This work was supported in part by the National Institutes of Health (R01 GM085575).

## References and links

**1. **A. Small and S. Stahlheber, “Fluorophore localization algorithms for super-resolution microscopy,” Nat. Methods **11**(3), 267–279 (2014). [CrossRef] [PubMed]

**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). [CrossRef] [PubMed]

**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). [CrossRef] [PubMed]

**4. **T. A. Laurence and B. A. Chromy, “Efficient maximum likelihood estimator fitting of histograms,” Nat. Methods **7**(5), 338–339 (2010). [CrossRef] [PubMed]

**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). [CrossRef] [PubMed]

**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). [CrossRef] [PubMed]

**7. **D. L. Donoho, “Compressed sensing,” IEEE Trans. Inf. Theory **52**(4), 1289–1306 (2006). [CrossRef]

**8. **L. Zhu, W. Zhang, D. Elnatan, and B. Huang, “Faster STORM using compressed sensing,” Nat. Methods **9**(7), 721–723 (2012). [CrossRef] [PubMed]

**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). [CrossRef] [PubMed]

**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). [CrossRef] [PubMed]

**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). [CrossRef] [PubMed]

**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). [CrossRef]

**13. **Y. Hua, “Estimating two-dimensional frequencies by matrix enhancement and matrix pencil,” IEEE Trans. Signal Process. **40**(9), 2267–2280 (1992). [CrossRef]

**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 (IEEE2004), 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). [CrossRef]

**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). [CrossRef]

**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). [CrossRef]

**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). [CrossRef]

**19. **J. M. Maciejowski, “Guaranteed stability with subspace methods,” Systems & Control Letters **26**(2), 153–156 (1995). [CrossRef]

**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). [CrossRef] [PubMed]

**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). [CrossRef]

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

**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). [CrossRef] [PubMed]