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² [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₁-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 ℝ$, 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,

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

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 ${\displaystyle {\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^−i2πk/N, k = 1, …, N (or equivalently, I − Ae^−i2πk/N, k = 1, …, N, is invertible), and ${\displaystyle {\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,

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 ℝ$, k_i = 1, 2, …, N_i, i = 1, 2, be a finite 2D sequence. For n_i = 1, 2, …, N_i, i = 1, 2, let

Moreover,

Proof. Let $\tilde{X}\left(k_1,k_2\right)\in ℝ$, k_i = 1, 2, …, N_i, i = 1, 2, be a finite 2D sequence. For n_i = 1, …, N_i, i = 1, 2, let

Decompose Q via SVD as Q = UΣV, where for r ∈ ℕ, $U\in {ℂ}^{N_1\times r}$, Σ ∈ ℂ^r×r and $V\in {ℂ}^{r\times N_2}$.

For n_i = 1, 2, …, N_i, i = 1, 2, define X₁(n₁) ∈ ℂ^1×r and X₂(n₂) ∈ ℂ^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),

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 ℝ$, k_i = 1, 2, …, 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

For s₁, s₂ ∈ ℕ and $A_i\in {ℂ}^{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 {ℂ}^{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 ℂ$, 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 {ℂ}^{1\times r}$, $b_{t_2}^2\in ℂ$, $c_{t_1}^1\in ℂ$, $c_{t_2}^2\in {ℂ}^{r\times 1}$, t_i = 1, 2, …, s_i i = 1, 2, for k_j = 1. 2, …, 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_i = 1, …, 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₁ × N₂ 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₁ × N₂ 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

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 ℝ$, k_i = 1, 2, …, N_i, i = 1, 2, represent the acquired image data.

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 σ₁ ≥ … ≥ σ_K ≥ 0, K = min(N₁, N₂), 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₁ and s₂ 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₁ s₂, t = 1, 2, and their corresponding product of coefficients in the numerator by. p_k ∈ ℂ, k = 1, …, s₁ s₂, 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₁, s₂) denote the number of single molecules, assuming that we retain s₁ and s₂ 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 {ℝ}^{2h}$, ${\widehat{\theta }}_n:=\left({\widehat{x}}_n,{\widehat{y}}_n\right)\in {ℝ}^2$, n = 1, …, h, such that

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 ℝⁿ. The maximum likelihood estimate ${\widehat{\theta }}_{mle}$ of θ ∈ Θ, for data incorporating Gaussian readout noise, is given by