1. Introduction and background

Spatial resolution in optical imaging was considered for more than a century to be limited by the wavelength of light, which is commonly known as Abbe's diffraction limit [1]. Today, and specifically in life sciences, modern microscopy techniques enable the recovery of subwavelength features, by utilizing photo-activated fluorescent molecules and assuming prior knowledge on the recovered signal. More specifically, the techniques known as STED [2], PALM [3], STORM [4] and dSTORM [5] exploit fluorescence to achieve a dramatic spatial resolution improvement, which has enabled imaging of subcellular features and organelles within biological cells with unprecedented resolution. Furthermore, structured illumination microscopy [6, 7] is an additional widely used super-resolution imaging modality due to low photo-toxicity and relatively fast acquisition rates. However, its ability to enhance the spatial resolution over the diffraction limit only by a factor of two a modest increase compared with alternative florescence imaging techniques. Although all fluorescence imaging techniques rely on attaching fluorescent molecules to regions of interest within the specimen, the specific imaging mechanisms differ from one another. PALM and STORM rely on the stochastic activation and deactivation of sequential sets of fluorescent molecules, and a subsequent localization procedure to find the positions of the fluorescing molecules with spatial resolution as good as 20nm [8]. In each cycle, only a small fraction of the molecules is activated, such that with high probability they are well isolated from one another, and therefore their centroids can be found with high precision. The final image is produced from smoothing a grainy image composed of the super-position of all localizations. STED on the other hand, is based on the depletion of spontaneous emissions in a deep subwavelength region and subsequent scanning of the sample, point by point. The STED microscope has an effective point spread function (PSF) smaller than the diffraction limited PSF, so that the spatial resolution of the final image is on the order of 50-60nm [9, 10]. Although these techniques enable imaging of sub-cellular organelles with unprecedented resolution and clarity, they share a common drawback in terms of temporal resolution: the time it takes to acquire all the data needed to recover a full image. STED requires scanning of the specimen, and although it recently achieved video rate speeds with spatial resolution of 60nm [10], its temporal resolution depends on the field of view of the imaged specimen. PALM/STORM techniques require tens of thousands of exposures to ensure that eventually all fluorophores are activated with high probability, leading to a long acquisition cycle, typically on the order of minutes [4].

To circumvent the long acquisition periods of PALM/STORM, an alternative technique - called super-resolution optical fluctuation imaging [11] (SOFI) - has emerged in recent years. SOFI uses higher simultaneously excited fluorophore densities to enable faster capturing of frames, thus achieving increased temporal resolution, or short integration time. Due to the high fluorophore density, the molecules overlap in each frame, making it difficult to apply a super-localization technique. Nonetheless, information is contained in the pixel-wise correlation between consecutive frames. Using a statistical measure known as cumulants [12], this correlative information is exploited to reduce the PSF size in the SOFI image. This resolution improvement scales with the square-root of the order P of the cumulant, effectively making the PSF narrower, such that features spaced more than $16X16$ of the diffraction limit may now be resolved. This reduction in PSF size can be further improved, by using PSF estimation from pixel cross-cumulants, leading to an overall resolution improvement which eventually scales linearly with the order of the cumulant [13]. In practice, however, although high fluorophore density reduces the number of required frames, thus far the spatial resolution offered by SOFI has not reached the super-resolution obtained by other fluorescent imaging techniques. This is because the use of statistical orders higher than two requires long data sequences to retain the same SNR level as in correlation based SOFI, thereby limiting its practical resolution improvement to a factor of two over the diffraction limit. In addition, as the statistical order increases, a phenomenon known as dynamic range expansion becomes more evident, in which fluorescence of strong emitters masks weak emitters [11], leading to contrast issues in high order SOFI.

With the understanding that SOFI fundamentally cannot reach the resolution of PALM/STORM, other methods have recently been proposed, employing multi-emitter fitting algorithms to reduce the number of required acquisition cycles [14–17]. Methods such as FALCON [15] and 3B [16] have shown their ability to achieve sub-diffraction resolution of $60nm$ and $50nm$ on endoplasmic reticulum and podosome structures, while handling data sets of overlapping emitters, thus demonstrating a temporal resolution of 2.5 and 4 seconds, respectively. Alternatively, Deschout et al. integrated both PALM and SOFI to image focal adhesions with 90-100nm resolution [17]. Two more recent high density methods are SRRF [18] and MUSICAL [19], which achieve spatial resolution of $50-100nm$ from movies which are on the order of 100 frames, and in some specific cases, even fewer. While FALCON is a frame-by-frame sparsity based method, MUSICAL relies on correlations estimated from a movie stack by implementing the multiple signal classification algorithm [20] (MUSIC). But, both methods do not exploit the fact that the emitters are uncorrelated, as in SOFI. On the other hand, SRRF estimates for each frame separately in the movie stack the local degree of gradient convergence on a sub-pixel grid. This new image exhibits an improvement in the spatial resolution compared with the original frame. It then uses temporal averaging, maximum projection or high order cumulants estimation with the inexplicit assumption of uncorrelated blinking, as in SOFI, to further increase the spatial resolution, but does not exploit sparsity.

Here, we present SPARCOM: sparsity-based super-resolution correlation microscopy. Our technique, for which we presented preliminary results in [21], displays short integration time with spatial resolution similar to PALM and STORM. It utilizes sparsity in the correlation domain, while assuming that the blinking emitters are uncorrelated over time and space, which can highly increase the number of detected locations of the sources, compared with sparse recovery performed on the signal itself (as in FALCON, for example). Namely, quoting from a recent review on sparse recovery using correlation information [22]: if existing algorithms can recover support of size m, then using such correlation information, the guaranteed size of recoverable support can be increased to O(m²). In a similar vein, exploiting sparsity in the correlation domain facilitates additional advantages such as noise reduction [22, 23]. Indeed, we demonstrate through numerical simulations, experimental data, and using Fourier ring correlation analysis, that SPARCOM achieves spatial resolution comparable to single molecule localization with ~50 frames only, thereby outperforming existing methods for fast acquisition. For example, using Fourier ring correlation analysis, we demonstrate reconstructions with an overall spatial resolution of ~41nm and ~45nm from movies with exposure duration 33 and 200 times shorter than the single molecule localization reconstruction, respectively. In addition to considerably increasing spatial resolution of the recovery, our technique has several other benefits. First, SPARCOM exhibits optical sectioning, rejecting out-of-focus emissions better than other techniques (such as SOFI), by relying on correlations only. Second, SPARCOM avoids the problem of dynamic range expansion since it relies only on second-order statistics, without resorting to high-order statistics for obtaining increased resolution. We demonstrate our technique on both simulated and experimental data. In the next section we provide a detailed mathematical description of SPARCOM, and demonstrate how we exploit sparsity in the correlation domain to achieve super-resolution.

2. Mathematical formulation

2.1 Compressed sensing

We begin with a short introduction to the concepts underlying sparse representations of signals [21, 24] and compressed-sensing [25–27] (CS). The Shannon-Nyquist theorem [20] dictates a minimal sampling rate to accurately describe a band-limited signal. In practice, many signals may be represented compactly in an appropriate basis, implying that fewer samples are needed for their recovery. CS theory states that a finite-dimensional signal which has a sparse representation in some basis, can be recovered stably from a small set of linear measurements chosen in a non-adaptive manner [20]. Many CS applications deal with signal recovery from measurements taken below the Nyquist rate of the signal, generally termed as sub-Nyquist sampling [20]. Such applications are found in various fields such as holography [28, 29], ultrasound [30, 31], communications [32–34], radar [35, 36], magnetic resonance imaging [37, 38] and more. Furthermore, sparsity can be used to recover high frequency content lost in the measurement process, generally termed super-resolution [39]. In optics, this idea was pioneered in the context of sub-wavelength optical imaging [40] and has since been expanded and demonstrated both in simulations and experiments in Ankylography [41], coherent diffraction imaging (CDI) and lensless imaging [42–45], phase retrieval [46–54], point-detector compressive fluorescence microscopy (CFM) [55] and super-resolution in quantum state tomography [56–58]. In this vein, here we rely on sparsity to achieve subwavelength imaging from correlative fluorescence measurements.

In certain applications, the unknown signal is not the quantity of interest, but rather we are interested in the locations of its non-zero values, or the support. Recent works in the context of array processing [22] and cognitive radio [59] have shown that, taking the measurements in the correlation domain (i.e. empirically estimating the auto-correlation matrix of a sequence of consecutive measurements) and recovering the second order statistics of the signal, facilitates support recovery of denser signals than when trying to recover the signal itself, or support recovery with less measurements. We rely on this result to facilitate super-resolution. We do that by using the measurements from a given fluctuations movie in the (sampled) correlation domain, and improve support recovery by estimating only the variance of each emitter, assuming that the emitters are uncorrelated.

2.2 Sparse recovery in the correlation domain

To better understand the advantages of performing sparse recovery in the correlation domain, consider Fig. 1. It is well known [20] that from a measurement of length m (denoted by $y_1$v in panel a of Fig. 1), one can recover vectors $x_1$ that are at most $m/2$ sparse, meaning, signals that have at most $m/2$ nonzero values. This is known as the single measurement vector (SMV) model. By relying on several measurements (e.g. consecutive frames, as denoted by $y_1,\dots ,y_T$ in panel b of Fig. 1), a set of T signals $x_1,\dots ,x_T$ with similar sparsity pattern of 0(m) nonzero entities can be recovered [21,22], using the same measurement matrix, A. In this case, also known as the multiple measurement vector (MMV) model, all recovered vectors have nonzero values at the same locations. Next, assume that in addition the non-zero entries of $x_1,\dots ,x_T$ are uncorrelated random variables, and so their cross-correlation matrix $R_x(0)$ (panel c in Fig. 1), is diagonal. The auto-correlation matrix of the measurements $R_y(0)$ is of size $mXm$ and is estimated from the measurements as in (b), while the auto-correlation matrix of the recovered signal is of size $TXT$ and is diagonal. Due to the uncorrelated entries assumption, $R_x(0)$ has the same number of unknowns as (a) in the form of the variances of $x_1,\dots ,x_T$, but we have now quadrupled the number of measurements in the form of $R_y(0)$, and can recover up to $O(m^2)$ entries from the diagonal of $R_x(0)$ [22]. Since we recover variances, the sparsity pattern of (c) is as that of (a) and (b). This directly translates to the detection of more sources (point emitters) than other algorithms can identify, from the same number of measurements. We illustrate this increase in emitter detection in Fig. 17 in the Appendix, where SPARCOM achieves higher hit-rates than the single-molecule localization procedure. Next, we give a full description of SPACROM and how super-resolution can be achieved by performing sparse recovery in the correlation domain.

 

Fig. 1 The difference between SMV, MMV and correlation based MMV. All models (a)-(c) have the same measurement matrix A. (a) Single measurement vector (SMV) model. For measurement $y_1$of length m, a sparse signal $x_1$with up to $m/2$ non-zero entries can be recovered. (b) Multiple measurement vector (MMV) model. Assume that $x_1,\dots ,x_T$ share a common support, denoted by their colored rectangles. From T measurements, denoted by $y_1,\dots ,y_T$, each of length m, $x_1,\dots ,x_T$ with sparsity level of at most $O(m^2)$ can be recovered. (c) MMV with additional assumption of uncorrelated entries of $x_1,\dots ,x_T$. In this model, support recovery of the emitters up to $O(m^2)$ can be achieved.

Download Full Size | PDF

2.3 The SPARCOM model and method

Before deriving the mathematics of SPARCOM, we first outline its main assumptions. As in SOFI, we assume that emissions by different emitters are uncorrelated over time, which provides further prior information to exploit. This assumption is reasonable since no Förster [60] or Dexter [61] homo-energy transfer is expected for sparsely located emitters (e.g. fluorophore distance > 10nm). In addition (unlike SOFI), we further rely on the fact that the emitters are sparsely distributed. These two assumptions are natural in fluorescence imaging.

Since we expect to recover fine features, smaller than the size of a single camera pixel, we recover our images on a much denser grid than the grid of the captured movie. As we show below, the discretization process does not limit the resolution of the recovered image, and leads to a clear separation of sub-wavelength features. To this end, we first introduce a Cartesian sampling grid with spacing ${\Delta }_L$, which represents the low-resolution sampling grid of the camera in Cartesian coordinates. Thus, each sampled frame $f$ of size $M X M$ pixels from the captured movie can be written as a sum of L active emitters

where $u[m{\Delta }_L,n{\Delta }_L],m,n=[0,\dots ,M-1]$ is the low-resolution sampled PSF, $s_k(t)$ represents the temporal fluctuations of the recorded intensity of emitter k and $[m_k,n_k]$ represent the location of the emitter on the 2D grid. We assume that the PSF of the imaging system, $u(\cdot ),$ is known (In practice, the PSF can be estimated from the data [13] or calculated based on the microscope's acquisition parameters as we show below in an experimental demonstration (Fig. 4).

Next, we introduce a high-resolution grid with spacing ${\Delta }_H$, on which we recover the variance values of the emitters, such that $\left[m_k,n_k\right]=[i_k,j_k]{\Delta }_H$ for some integers $i_k,j_k=0,\dots ,N-1$, where $N$ can be considerably larger than $M$. With this high-resolution grid, rewriting (1) in Cartesian form yields

For simplicity, we choose ${\Delta }_L=P{\Delta }_H$ for $P\ge 1$. Thus, the low-resolution movie is captured over a grid of $MXM$ pixels, while the high-resolution reconstruction is performed over an $NXN$ grid, with $N\gg M$. Our goal is to find (in a super-resolved image of $N^2$ pixels) the locations of the actual L emitters. If the algorithm detects an emitter at position $[i_k,j_k]$, then it places its estimated variance value at that position, while if no emitter is detected, then a value of zero is placed instead.

Rewriting (2) more compactly,

In (3) we use a double sum to describe the possible coordinates of the emitters on the high-resolution grid where only L such emitters are actually present. To achieve a numerically efficient algorithm, we consider (3) in the discrete Fourier domain. Thus, matrix-vector multiplication operations are based on the fast Fourier transform (FFT), which has a computational complexity of r log(r), where r is the length of the vector over which the FFT is applied.

Next, we apply an M -point spatial (two dimensional) discrete Fourier transform (DFT) on (3)

where $W_N=e^{-j\frac{2\pi }{N}}$, which yields

with $U\left[k_m,k_n\right],k_m,k_n=0,\dots ,M-1$, the $MXM$ DFT of the discretized PSF. Denote the column-wise stacking of each frame $F[k_m,k_n,t]$ as a length- $M^2$ vector $y(t)$ and let $s(t)$ be a length- $N^2$ vector stacking of $s_{ij}(t)$. With this notation, (5) becomes

where $A=H(F_M\otimes F_M)$, $\otimes$ is the Kronecker multiplication operator, $F_M$ is a partial $MXN$ DFT matrix which consists of the first M rows of the $NXN$ DFT matrix and $H$ is the $M^{2}X{M}^2$ diagonal matrix

The vector $s(t)$ is an L sparse vector, that is only L elements are non-zero.

Using (6), the auto-correlation matrix of $y(t)$ is given by

for some time-lag $\tau$. Since the emitters are uncorrelated, $R_s(\tau )$ is a diagonal matrix with an L -sparse diagonal $r_s(\tau )$. Taking $\tau =0$ for simplicity, our goal is the recovery of this diagonal, which represents the variance values of each of the emitters on a high-resolution grid. It is important to note that the support of $r_{\text{s}}(\text{τ})$ is equivalent to the support of $s(t)$ which corresponds to the positions of the emitters on the high-resolution grid.

In practice, we estimate $R_y(\tau )$ using the empirical correlation,

with ${(\cdot )}^H$ denoting the complex-conjugate operation and $\overline{y}={\displaystyle \sum }_{\text{t}=1}^{\text{T}}y(\text{t})$.

Since $R_s(\tau )$ is diagonal, it holds that $R_{\text{y}}\left(0\right)={\displaystyle \sum }_{\text{l}=1}^{{\text{N}}^2}{a}_{\text{l}}{a}_{\text{l}}^{\text{H}}{r}_{{\text{s}}_{\text{l}}}(0)$ with $r_{s_l}(0)$ being the $l$ th entry of $r_{\text{s}}(0)$. Denoting $x=r_{\text{s}}(0)$, we suggest to recover $x$ by formulating the following convex optimization problem [62],

with

Here $\lambda \ge 0$ is a regularization parameter, $\left|\right|\cdot |{|}_F$ denotes the Frobenius norm, $a_l$ is the l th column of $A$ and $x_l$ is the l th entry of $x$. The $\left|\right|\cdot |{|}_1$ -regularizer promotes sparsity of $x$ and $f(x)$ enforces consistency of the equations in (7). The vector $x$ can also be reweighted [63] every few iterations (prechosen by the user) by a diagonal weighting matrix $W$ (starting from $W$ being the identity matrix). Each entry in the diagonal of $W$ equals the inverse value of the absolute value of the corresponding entry in $x$ plus a small non-negative value. We solve (8) iteratively using the FISTA approach [64–66], leading to Algorithm 1.

The FISTA algorithm requires knowledge of the gradient of (9). Differentiating (9) with respect to x results in the following gradient

where $M=|A^{\text{H}}A{|}^2$ and $v={[a_1^{\text{H}}{R}_{\text{y}}(0)a_1,\dots ,a_{{\text{N}}^2}^{\text{H}}{R}_{\text{y}}(0)a_{{\text{N}}^2}]}^{\text{T}}$ and we choose the zero time-lag for variance estimation. The application of the absolute value-squared is performed element-wise and ${(\cdot )}^T$ denotes the transpose operation. In practice, (10) is implemented efficiently using FFT and inverse FFT operations only [67]. In Algorithm 1, $L_f$ is the Lipschitz constant of (10), and is equal to the maximum eigenvalue of $M$. To better understand the main processing steps of SPARCOM, a block diagram of our technique is given in Fig. 6 in the Appendix.

Algorithm 1: SPARCOM via FISTA for minimizing (9)

View Table | View all tables in this article

As a final note on the difference between SPARCOM and SOFI, consider the auto-correlation $R_y(\tau )$ of the acquired movie. For a prechosen value of $\tau$, auto-correlation and cross-correlation SOFI estimate only parts of $R_{\text{y}}\left(\text{τ}\right)$, and from it SOFI constructs the final image (e.g. the diagonal of $R_{\text{y}}\left(\text{τ}\right)$, in auto-correlation SOFI). In doing that, SOFI does not estimate the entire auto-correlation matrix of the input movie, and disregards additional information that is present in the off-diagonal elements of $R_{\text{y}}\left(\text{τ}\right)$. In particular, SOFI does not make explicit use of the fact that the emitters are uncorrelated. On the other hand, SPARCOM explicitly uses this structural knowledge of the temporal fluctuations of the emitters, effectively reducing the number of degrees of freedom in the reconstruction problem. This leads to improved reconstruction results in terms of spatial resolution, as we present in the next section, on simulated and experimental results. After examining many frames and cross-sections from different movies, we can reliably say that we expect SPARCOM to perform as good as low density single-molecule localization in terms of spatial resolution, but with much fewer frames and hence a much faster acquisition time. Actually, in a typical PALM/STORM experiment, the user specifically adjusts the acquisition process to achieve well-isolated diffraction-limited spots by controlling the excitation laser intensity, which imposes a time restriction on the maximum frame-rate of the imaging cycle. In contrast, SPARCOM does not impose such a restriction, hence the temporal resolution of SPARCOM is only limited by the frame-rate of the camera. In the next section we provide several examples of our SPARCOM technique, on data from simulations and from actual experiments.

3. Results

3.1 Simulation results

We begin with simulation results given in the panels of Fig. 2. The panels show the simulation and reconstruction of microtubules from a freely available movie [68,69] of 361 high density frames. Figure 2(a) shows the simulated ground truth of the image with subwavelength features, while Fig. 2(b) illustrates the diffraction-limited image, obtained by summing all the 361 frames in the movie. Figures 2(c)-2(d) show SRRF reconstructions from the 361 frames movie and a higher density movie of 60 frames (summing every 6 frames from the original movie), respectively. The emitter density in the 361 frames movie is $~35emitters/\mu m^2$, corresponding to an average of 224-225 active emitters per frame (field of view is $6.4\mu m^2$), while the emitter density of the 60 frames movie is 6 times denser. Figures 2(e)-2(f) depict MUSICAL reconstructions from the 361 and 60 frames movies, respectively, while Figs. 2(g)-2(h) illustrate SPARCOM reconstructions from the same movies, respectively. Red boxes indicate zoomed-in areas; the lower left red box in each panel is zoomed-in, in the adjacent lower right boxes. All MUSICAL and SPARCOM images are smoothed with the same known $64X64$ pixels Gaussian PSF, to achieve a fair comparison, which should not be affected by the smoothing process. Its parameters are listed in the methods section in the Appendix.

 

Fig. 2 Simulation results. Upper row: (a) Ground truth: high resolution image of simulated microtubules. (b) Diffraction-limited image. (c)-(d) SRRF reconstructions from high and very-high density movies of 361 and 60 frames (6 times denser), respectively. Lower row: (e)-(f) MUSICAL reconstructions from high and very-high density movies of 361 and 60 frames, respectively. (g)-(h) SPARCOM reconstructions from high and very-high density movies of 361 and 60 frames, respectively. Comparing the reconstructions for the same number of frames (e.g., the magnified regions in the red boxes), clearly the ability of SPARCOM to separate between closely adjacent subwavelength features, outperforms both SRRF and MUSICAL, even with as few as 60 frames. These separations closely match those in the ground truth.

Download Full Size | PDF

Visually comparing all panels, we observe that all three methods achieve super-resolution, similar to the ground truth in panel a, and recover features which are missing in the diffraction limited image of panel b. Closer examination, especially features which are magnified in the red boxes, reveals that SPARCOM achieves better reconstruction in terms of spatial resolution over SRRF and MUSICAL, both for the movies of 361 and 60 frames. For example, the white arrows in the enlarged red box in panel c demonstrate that SRRF is unable to separate the borders of the microtubules (also indicated by the red arrows in panels c and d). These borders are vague in the MUSICAL reconstructions as well, but are clearly visible in the SPARCOM images.

Although SPARCOM exhibits better separation of adjacent features, some artifacts appear in the reconstruction as discontinuity of the microtubules. Considering panels g (361 frames) and h (60 frames of increased density) of Fig. 2, we can clearly see that as we reduce the number of frames, but simultaneously increase the emitters density, some parts of the filaments seem broken, as indicated for example by the white arrows in Fig. 2(h). This distinction emphasizes a clear trade-off in SPARCOM between the emitter density/number of frames to the creation of these artifacts. Even for the movie of 50 frames with increased density, SPARCOM exhibits better spatial resolution and clear separation of the microtubules, but at a cost of broken filaments. Thus, the user may wish to control the density of emitters, depending on the application at hand. Using a few hundreds of frames will give a clear reconstruction, while reducing the number of frames (and increasing the emitter density to shorten integration time) might result in some discontinuity, but will still manage to separate the fine features of the imaged specimen. Alternatively, it is possible to perform SPARCOM recovery with high density to reveal the fine structures of the images specimen in short integration time, but complement SPARCOM with reconstructions such as SRRF or MUSICAL to fill the areas of discontinuity.

We now turn to a quantitative comparison based on the panels presented in Fig. 2, by considering selected intensity cross-sections taken along the yellow lines in several panels of Fig. 2. The upper (lower) panel of Fig. 3 shows the selected intensity profile marked by the solid (dashed) yellow line in Fig. 2. The lines scans in Fig. 3 where chosen for two reasons. The dashed line scan was chosen at a region were the filaments are closest to one another, to stress that in the case of high density and minimal number of frames, SPARCOM manages to resolve the filaments, and the other methods fail. The second line scan (solid) was chosen over a larger area arbitrarily to show that in general and without any “fixing” of the data – SPARCOM exhibits very good reconstruction.

 

Fig. 3 Intensity profiles comparing the performance of SRRF, MUSICAL and SPARCOM. Normalized cross-sections taken along the solid yellow line (upper) and the dashed yellow line (lower) of Fig. 2, comparing the ground truth (dash-dot green, Fig. 2(a)), diffraction-limited image (dot blue, Fig. 1(b)), SRRF image recovered from a movie of 60 frames (solid yellow, Fig. 2(c)), MUSICAL image recovered from a movie of 60 frames (dash red, Fig. 2(f)), and SPARCOM image recovered from the high density movie of 60 frames (purple circle-head, Fig. 2(h)). These panels substantiate the visual conclusions drawn from Fig. 2: SPARCOM recovers the profiles of sub-diffraction spaced microtubules from a high emitters' density movie clearly and with a good agreement to the ground truth, while SRRF and MUSICAL do not achieve such high resolution between close features.

Download Full Size | PDF

Comparing these cross-sections substantiates the conclusions drawn from the visual comparison: SPARCOM, taken with 60 frames only, displays resolution similar to the ground truth, while improving upon SRRF and MUSICAL. For example, the peak to peak distance of the SPARCOM reconstruction in the lower panel corresponds to $162.7nm$ (ground truth peak to peak distance is $122.4nm$). Observing the upper panel of Fig. 3, one can see a small shift of 23nm in the two inner filaments in the SPARCOM reconstruction, compared with the ground truth (SRRF exhibits a shift of 61nm in the second filament from the left). This shift is marginal and clearly unobservable in the panels of Fig. 2. Although SPARCOM exhibits some shift compared with MUSICAL, it also shows improved separation and resolution. This bias was unobservable for the same comparison with the SPARCOM reconstruction from the movie of 361 frames. Figures 8, 9, 10, and 11 in the Appendix provide additional examples and quantitative comparisons to FALCON and SOFI.

3.2 Experimental results

We next test SPARCOM on actual experimental data, which inevitably contains noise and other non-ideal aspects that always accompany experiments (e.g. accuracy of the PSF knowledge). Ideally, we would like to compare high-density SPARCOM to low-density PALM/STORM reconstruction and show similar spatial resolution, but at a much higher frame-rate. However, such a comparison requires performing two different experiments, often with two different labeling mechanisms and fluorophores on the same specimen. Such a process inevitably results in discrepancies between the spatial distributions of the emitters

within the specimen, which in turn results in recoveries of non-overlapping areas of the same specimen, as was demonstrated in Fig. 1 of the 3B method [16] (a comparison to 3B is given in the Appendix in Fig. 14). We therefore demonstrate the advantages of SPARCOM using an experimental data set of intermediate density, showing regions of low density and regions of high density of fluorescing emitters.

Figure 4(a) depicts the diffraction limited image, created by summing all the frames in the movie, while Fig. 4(b) shows a single diffraction limited image from a 500 frames movie. The density of simultaneously activated emitters is such that some areas contain many overlapping emitters, while in other areas well isolated emitters can be found. Figures 4(c)-4(d) show SRRF reconstructions from the movie of 500 frames and from a movie 10-times denser of 50 frames (by summing every 10 frames from the 500 frames movie), respectively. Figures 4(e)-4(f) illustrate MUSICAL reconstructions for the same movies, respectively, while Figs. 4(g)-4(h) depict SPARCOM reconstructions for the same movies, respectively. The enlarged areas in the red boxes at the lower left corners of each panel correspond to the areas marked by the small red boxes in the upper left corner of each panel. All MUSICAL and SPARCOM images were smoothed with the same $64X64$ pixels Gaussian PSF which was used in the SPARCOM recovery process, for consistency. Its parameters are listed in the methods section in the Appendix.

 

Fig. 4 Experimental results. Panels (a)-(b) illustrate a diffraction-limited image of microtubules and a single noisy frame (from a movie of 500 frames). (c)-(d) SRRF reconstructions from a movie of 500 frames and a 10-times denser movie of 500 frames, respectively. (e)-(f) MUSICAL reconstructions for the same movies. (g)-(h) SPARCOM reconstructions for the same movies. Red boxes in the lower left corners of each panel are enlarged regions of the small red boxes in the upper left corner of each panel, respectively. Judging visually, both SPARCOM and MUSICAL exhibit clear separation of the microtubules within the enlarged boxes, while the SRRF reconstructions do not, for both movies. (i, j) Intensity cross-sections (normalized) taken along the solid (i) and dashed (j) yellow lines. Panel (i) shows that in regions of high emitters' density, SPARCOM clearly achieves better spatial resolution than both SRRF and MUSICAL, depicting a clear bifurcation, which is absent in MUSICAL and SRRF. On the other hand, panel (j) shows that in areas of low emitters' density, the spatial resolution of SPARCOM is similar to SRRF and better than MUSICAL.

Download Full Size | PDF

Visual inspection of Fig. 4 reveals several observations. First, both MUSICAL and SPARCOM achieve clear separation of the microtubules depicted in the enlarged boxes, while SRRF does not, for both movies of increasing emitter density. But, for the case of 500 frames, the separation of SPARCOM (panel g) seems better than that of MUSICAL (panel e). This observation substantiates the ability of SPARCOM to better resolve sub-diffraction features than the other methods considered here. Second, closer inspection reveals that SRRF contains residual noise between the microtubules, which is absent in the MUSICAL and SPARCOM images.

We now quantify the separation ability of SPAROCM compared with SRRF and MUSICAL by examining intensity profiles along the solid (Fig. 4(i)) and dashed lines (Fig. 4(j)) for the recoveries taken from the movie of 50 frames. Panel 4(i) corresponds to a scenario of close by microtubules, while panel 4(j) corresponds to a scenario of a single, isolated microtubule. In the first scenario, SPARCOM (purple circle-head) clearly shows separation of the two microtubules (peak to peak distance of $109nm$). This clear separation is absent in both the MUSICAL (dash red) and SRRF recoveries (solid yellow), though it seems there is an agreement between the SPARCOM and MUSICAL reconstructions. In the case of the isolated microtubule (Fig. 4(j)), there is a clear agreement between SPARCOM and SRRF regarding the width and location of the microtubule, while the MUSICAL recovery seems wider and shifted.

As in Fig. 2, we would like to address the appearance of artifacts in the recovery. A similar trade-off between the emitter density and appearance of artifacts (as broken filaments) is evident in panels g and h. Furthermore, a narrowing of the filaments is evident as the emitter density increases, but this phenomenon is common to all the methods considered here (visually comparing between the pairs of panels c-d, e-f and g-h). Yet, as we demonstrate visually and using the intensity cross-sections, SPARCOM exhibits better separation compared to state-of-the-art methods. Additional inspection shows that a similar trade-off exists also in the SRRF images, where increasing emitter density results in discontinuous filaments. Although MUSICAL recoveries seem smooth, a clear separation between very close filaments was not achieved, as in SPARCOM (Fig. 4(i)). Additional comparisons to FALCON and SOFI are given in Fig. 12 in the Appendix.

3.3 Fourier ring correlation analysis

To substantiate the ability of SPARCOM to achieve a similar spatial resolution as single-molecule localization, but with a much shorter integration time, we carry out a comparative study using Fourier ring correlation analysis (FRC) [70], given in Fig. 5. FRC analysis is considered a reliable method for estimating the spatial resolution of different reconstructions, and has been applied successfully in electron microscopy and super-resolution localization microscopy. The key concept in FRC is to correlate between the Fourier transforms of two reconstructions of the same object over concentric rings of constant spatial frequencies. The overall spatial resolution is determined when the correlation drops below a predefined threshold. First, we perform FRC analysis for the ThunderSTORM, FALCON, MUSICAL, SRRF and SPARCOM recoveries by comparing them with the ground truth object given in Fig. 2(a). In this way, a high-resolution reconstruction will have a high degree of correlation with the ground truth for a wide range of spatial frequencies, until at a certain point it will drop below a predefined threshold which represents its overall resolution. As customary, the resolution criteria was chosen as a fixed 1/7 drop in the value of the normalized correlation (represented by the horizontal black line in the panels of Fig. 5).

 

Fig. 5 Fourier Ring Correlation (FRC) analysis. (a)-(b) FRC analysis from the simulated movies of 361 frames and 60 frames of increased emitter density from Fig. 2, respectively. (c) FRC analysis from the movie of 500 frames of the experimental data set in Fig. 4. Colored curves are the FRC curves and the horizontal black lines represent the fixed 1/7 resolution threshold. Resolution corresponds to the intersection points of the curves with the horizontal threshold lines. In both (a) and (b), the resolution presented by SPARCOM matches that of ThunderSTORM (perfect reconstruction), even when the number of frames decreases to 60, and is higher than the other methods. For example, the spatial resolution of FALCON in both cases is lower by a factor of ~2 compared with SPARCOM. In (c), both SPARCOM and SRRF demonstrate similar resolution, while MUSICAL shows decreased resolution.

Download Full Size | PDF

Panels a and b of Fig. 5 show FRC curves for two cases. Panel a compares the FRC curves for ThunderSTORM (12,000 frames, red), SPARCOM (purple), FALCON (cyan), SRRF (blue) and MUSICAL (orange) for the simulated movie of 361 frames, while panel b shows the same comparison for the simulated movie with increased density of 60 frames. For the quantitative comparison, each curve is a smoothed FRC with a Savitzky-Golay filter [71] with a window of 19 points and a polynomial degree of 1. Inspecting both panels a and b, we see that the resolution of SPARCOM matches that of single molecule ThunderSTORM with 12000 frames, and is higher or matched to other state-of-the-art methods (the resolution gain for SPARCOM and MUSICAL for the 60 frames movie is comparable). Table 1 summarizes the resolution gains for each method in both movies and shows for example that the resolution gain of SPARCOM is ~2 times better than FALCON and ~2.5 better than SRRF.

Table 1. FRC resolution gains for panels a and b of Fig. 5, showing the gain at the intersection points of the FRC lines with the fixed 1/7 threshold, depicted as a black vertical line.

View Table | View all tables in this article

The second FRC analysis was performed over the experimental data set presented in Fig. 4. We compare in panel c between MUSICAL, SRRF and SPARCOM by taking 400 frames from the original movie and dividing them into two movies of 200 frames each. We perform recovery with each method over each of the two movies and calculate the FRC curves. Observing panel c, we can see that SPARCOM and SRRF achieve comparable resolution gain (though SRRF has a slight advantage), and both are better than MUSICAL. Table 2 substantiates this observation (23.1nmdifference between SPARCOM and MUSICAL).

Table 2. FRC resolution gains for panel c of Fig. 5, showing the gain at the intersection points of the FRC lines with the fixed 1/7 threshold, depicted as a black vertical line.

View Table | View all tables in this article

Finally, Fig. 5 demonstrates that the spatial resolution of SPARCOM is preferable over other state-of-the-art methods, while achieving similar values to low density ThunderSTORM, with only a fraction of the frames.

In summary, Figs. 2-5, and additional examples in the Appendix, provide direct numerical and experimental evidence that SPARCOM presents several clear advantages over both single-molecule localization, SOFI, and multi-emitter fitting algorithms (Fig. 13 in the Appendix). SPARCOM displays spatial resolution similar to that of single-molecule localization imaging, while offering up to 200 times reduction in the total acquisition duration. Figures 15, 16, 17, 18, and 19 in the Appendix provide additional quantitative metrics to objectively assess the performance of SPARCOM against the other methods considered here, as a function of increasing molecular density and minimal separation distance. These metrics help substantiate the conclusions drawn here and discussed above.

4. Conclusions and discussion

To conclude, we presented a sparsity based method, called SPARCOM, for super-resolution in fluorescence microscopy. Using SPARCOM, we were able to reconstruct super-resolved structures with high emitter densities and many overlapping emitters in each exposure. SPARCOM enjoys short integration time which in turn leads to high temporal resolution, while achieving spatial resolution comparable to PALM and STORM. We compared the performance of our technique on simulated and experimental data sets against state-of-the-art high density methods, and demonstrated the recovery of subwavelength features that were missing in them. Additionally, we quantified the performance of SPARCOM and single molecule localization in a series of scenarios with increasing emitters’ density. As the density increases, SPARCOM exhibits better support recovery compared with single molecule localization. The ability of SPARCOM to recover the correct support of emitters with high emitters' density leads to a dramatic reduction in the number of required exposures to produce a single super-resolved image. At very high emitter densities, SPARCOM exhibits some reconstruction artifacts, such as occasional filament discontinuity. This may be the result of the empirical auto-correlation estimation from a few number of movie frames. When lower densities are concerned, these artifacts are alleviated, thus a trade-off exists between the overall integration time and the formation of artifacts. Altogether, SPARCOM displays greatly improved temporal resolution while retaining spatial resolution comparable to PALM/STORM.

As a final note, we would like to point out the inherent difference between the sparsity assumptions and the way they are being used in PALM/STORM (or any single-molecule fitting algorithm) and in SPARCOM. Consider a Region Of Interest (ROI) defined by a single diffraction-limited spot, as viewed in the camera. PALM/STORM employs a strong assumption that within this ROI, there is only a single emitter, that is, a signal with sparsity one. The measurement procedure in PALM/STORM is thus specifically adjusted to fulfill this condition, by forcing only a small subset of (well isolated) emitters to fluoresce. Of course, in reality this condition can only be fulfilled partially, and inevitably the measurements contain areas of overlapping emitters as well as isolated ones, due to the stochastic nature of the activation process of the emitters. On the other hand, SPARCOM exploits sparsity in a different way, which allows for the use of higher densities of excited emitters, and therefore leads to improved temporal resolution. First, as in correlation-based SOFI, SPARCOM uses only second-order statistics and relies on the fact that the emission of light from different emitters is uncorrelated. This in itself leads to an overall resolution improvement by a factor of 2 over the diffraction limit. Second, to achieve an even greater resolution improvement, we exploit the sparse nature of the emitters, namely, we assume that every diffraction-limited spot arises from a small number of emitters – not necessarily a single one. Using the sparse distribution of the emitters together with their uncorrelated emissions leads to sparsity in the correlation domain. This, in turn, increases the number of possible detected emitters, even more than by exploiting sparsity without the statistical prior of uncorrelated emissions [8,50]. In the Appendix we provide additional examples which show that exploiting sparsity in the correlation domain leads to much better results than performing sparse recovery on each image and combining all single-frame reconstructions together (Fig. 17).

Appendix

1. SPARCOM block diagram

The following block diagram presents the main building blocks of SPARCOM. Each movie frame is transformed into the frequency domain, followed by the estimation of the auto-correlation matrix of the vectorized frames in the movie for a pre-chosen time-lag (e.g. for the zero time-lag, the covariance matrix of all the pixels in the movie is estimated). This matrix serves as the measurement input to the sparse solver, alongside an estimate of microscope’s PSF. The sparse solver is a specifically tailored FISTA implementation. The output of the algorithm is the variance estimation (if the zero time-lag is chosen) for the emitters on a high-resolution grid, as detailed in Algorithm 1 in the main paper.

 

Fig. 6 Block diagram of the main elements of SPARCOM. Each frame from the diffraction limited movie is converted into the Fourier domain using the fast Fourier transform (FFT2 in Matlab). Afterwards, the auto-correlation matrix (for a pre-chosen time-leg) of the entire Fourier transformed data set is estimated empirically (displayed here in logarithmic scale). Using an estimate of the PSF, a sparse solver (based on a tailored FISTA implementation) recovers the super-resolved image on a denser grid. All SPARCOM reconstructions provided in this work are generated on a high-resolution grid 8 times denser than the low-resolution grid of the input movies.

Download Full Size | PDF

2. Endoplasmic Reticulum example

The following example highlights the ability of SPARCOM to yield reconstructions of sub-cellular organelles, other than microtubules, from fast acquisition movies with nanometric resolution. This example demonstrates SPARCOM imaging of the endoplasmic reticulum (ER) protein (reticulon-4), fused to tdEos in a U2OS cell. The dataset is freely available, including the sample preparation protocol. The imaging wavelength is $561nm$, with NA of 1.3 and each pixel is $100nm$. The input sequence is a 160 frames movie (acquisition time of 2.5 seconds) of $512X512$ pixels each, and the resulting SPARCOM image is of size $512X512$ pixels. Here, we recover the SPARCOM image by assuming sparsity in the wavelet domain, using a Daubechies wavelet decomposition of order 2 with 8 taps. Figure 7 shows the diffraction limited image (a), a single frame (b) and the SPARCOM image (c). Observing panels (a) and (c) reveals that the ER network has many interconnections and an elaborate structure, demonstrating SPARCOM over additional structures than the ones considered in the main paper. White markers in panel c indicate a FWHM of ~130nm.

 

Fig. 7 Experimental reconstruction of ER in a U2OS cell. Panels (a)-(c) illustrate a diffraction-limited image, a single noisy frame and the SPARCOM recovery, under sparsity assumption in the wavelet domain. SPARCOM achieves a clear and smooth depiction of the ER, while rejecting most of the background noise in the raw movie (b).

Download Full Size | PDF

3. Additional simulation results and comparisons

3.1 Comparison with ThunderSTORM, SOFI and FALCON

We begin with additional comparisons to the dataset given in Figs. 2 and 3 in the main manuscript. Our main purpose here is to perform additional comparisons to single molecule localization reconstruction (implemented using ThunderSTORM) with 12000 frames (considered here as the gold standard) and to FALCON. FALCON is a sparsity based super-resolution fluorescence microscopy method which recovers the super-resolved structure frame-by-frame, and does not rely on correlations, as SPARCOM. The simulation setup is the same as described for Figs. 2 and 3. The emitter density in the 12000 frames movie is $~0.94emitters/\mu m^2$, corresponding to an average of 6-7 active emitters per frame (field of view is $6.4\mu m^2$).

Figure 8(a) shows the simulated ground truth, while Fig. 8(b) illustrates the diffraction-limited image, obtained by summing all the 361 frames in the movie. Figure 8(c) shows a smoothed single molecule localization reconstruction from a low-density movie of 12000 frames of the same simulated microtubules and the same number of emitters (the image is constructed using the ThunderSTORM plugin for the Image-J software). Panel c illustrates the best possible recovery of single molecule localization, when temporal constraints are irrelevant. Figure 8(d) illustrates 4th order SOFI recovery (absolute value, zero time-lag) from the 361 frames movie. On the other hand, Figs. 8(e) and 8(f) correspond to FALCON recoveries from the movies of 361 and 60 frames, respectively (FALCON reconstructions are done using a Gsigma value of 1.5, up-sampling ratio of 8, sparsity level of 2 and the built in PSF model of FALCON). Lastly, SPARCOM recoveries for the same movies are given in Figs. 8(g) and 8(h), respectively. All ThunderSTORM, FALCON and SPARCOM images are smoothed with the same known PSF, to achieve a fair comparison, which should not be affected by the smoothing process.

 

Fig. 8 Additional simulation results. Upper row: (a) Ground truth: high resolution image of simulated microtubules. (b) Diffraction-limited image. (c) Smoothed ThunderSTORM recovery from a low-density movie of 12000 frames. (d) 4th SOFI recovery (absolute value, zero time-lag) from a high-density movie of 361 frames. Lower row: (e)-(f) FALCON reconstructions from movies of 361 frames and 60 frames of increased densities, respectively. (g)-(h) SPARCOM reconstructions for the 361 and 60 frames movies, respectively. Comparing the reconstructions for the same number of frames (e.g., the magnified regions in the red boxes), clearly the ability of SPARCOM to separate between closely adjacent subwavelength features is superior to FALCON and SOFI. In this example, SPARCOM has an acquisition rate 33 (361 frames) and 200 (60 frames) times faster than PALM/STORM.

Download Full Size | PDF

Comparing Figs. 8(a) and 8(c), it is evident that the ThunderSTORM (panel c) recovery clearly detects all the microtubules correctly, and represents very accurate recovery, using 12000 frames of low emitters’ density. In contrast, the SOFI reconstruction (panel d) does not present a major resolution improvement when constructed from the corresponding 361 high emitters’ density movie. In other words, SOFI works much faster than single molecule localization but its spatial resolution is insufficient. Comparing these to SPARCOM (panels g-h) we observe that the temporal resolution of SPARCOM is far better than that of the low-density, single molecule localization reconstruction, while achieving similar spatial resolution: in this example, SPARCOM performs well with only 361 and 60 frames, compared with the 12000 frames of ThunderSTORM, needed to produce reliable recovery. This dramatic reduction in the number of required frames from 12000 to 361 and even 60 corresponds to a 33 and 200 times faster acquisition rate of SPARCOM compared with single molecule localization.

Next, comparing FALCON and SPARCOM recoveries from a 361 and 60 frame movies (panels S3e-S3h) leads to several conclusions. When the filaments are well separated from each other, both methods manage to recover them. However, when considering the zoomed-in areas outlined by the red frames, it is evident that for very closely situated sub-diffraction features, SPARCOM performs considerably better than FALCON. That is, the SPARCOM images for both 361 and 60 frames are clearer than the FALCON recoveries and very closely resemble the ground truth image and the ThunderSTORM recovery. This comparison suggests that SPARCOM achieves better resolution than all other contemporary sparsity-based recovery methods.

To obtain a more quantitative comparison between SPARCOM, SOFI and single molecule localization, we compare intensity cross-sections, taken along the yellow lines in several panels of Fig. 8. Figure 9 depicts two selected intensity profiles. The upper (lower) panel shows the selected intensity profile marked by the solid (dashed) yellow lines in Fig. 8, respectively. In both panels the ground truth is given by the dash-dot green line, (Fig. 3(a)), diffraction limited profile by dot blue line (Fig. 3(b)), smoothed ThunderSTORM from a 12000 frames movie by the solid yellow line (Fig. 3(c)), FALCON image recovered from the high density movie of 60 frames by the dash red line (Fig. 3(f)) and SPARCOM image recovered from the high density movie of 60 frames by the purple circle-head line (Fig. 3(h)).

 

Fig. 9 Intensity profiles comparing the performance of ThunderSTORM, FALCON and SPARCOM. Normalized cross-sections taken along the solid yellow line (upper) and the dashed yellow line (lower) of Fig. 8, comparing the ground truth (dash-dot green, Fig. 8(a)), diffraction-limited image (dot blue, Fig. 8(b)), smoothed ThunderSTORM recovered from the low density movie of 12000 frames (solid yellow, Fig. 8(c)), FALCON image recovered from the high density movie of 60 frames (dash red, Fig. 8(f)), and SPARCOM image recovered from the high density movie of 60 frames (purple circle head, Fig. 8(h)). These panels substantiate the visual conclusions drawn from Fig. 8: SPARCOM recovers the profiles of sub-diffraction spaced microtubules from a high emitters' density movie clearly and with a good agreement to the ground truth and ThunderSTORM, while FALCON does not achieve such high resolution.

Download Full Size | PDF

Comparing these cross-sections substantiates the conclusions drawn from the visual comparison: our sparsity-based method, SPARCOM, taken with 60 frames only, displays resolution similar to the ground truth and to the ThunderSTORM reconstruction from 12000 low density frames. Clearly, SPARCOM achieves considerably better resolution than FALCON: the FALCON lines show spurious peaks between the filaments and the width of each true peak appears considerably wider than in the ground truth, and in the SPARCOM reconstruction. To achieve this high spatial resolution, ThunderSTORM requires a clear separation between the far field emissions of the emitters, which implies that each frame in the movie needs to contain a much smaller number of well isolated diffraction limited spots. This always comes at a cost of degraded temporal resolution. In contrast, with only 60 frames and increased emitters' density, SPARCOM reliably recovers the true support of the filaments.

3.2 Additional example

In the following, we present additional simulated results comparing SPARCOM with SOFI and single molecule localization (implemented using the ThunderSTORM ImageJ plug-in), further supporting our conclusions of the main paper.

Figures 10 and 11 provide the simulation and reconstruction of microtubules from a high-density movie containing 1000 noisy frames with an additional out-of-focus layer, corresponding to a distance of $1\mu m$ from the focal plane (represented by the thick lines in Fig. 5(c)). Each movie frame is $64X64$ pixels and the pixel size is ${160}^{2}nm$. The imaging wavelength is $800nm$ and the numerical aperture (NA) is 1.4, to simulate acquisition with a high NA microscope. Since this is a simulated dataset, the PSF is a-priori known and does not need to be estimated. Figure 10(a) shows the simulated ground truth of the image with subwavelength features. The image is of size 512 X 512 pixels. Figure 10(b) presents the diffraction limited image, created by summing all of the frames in the movie and Fig. 10(c) show a single diffraction limited image from the high-density movie. It is evident that - in each frame - most of the emitters are active, and their corresponding PSFs overlap significantly.

 

Fig. 10 Simulation results comparing between SPARCOM, single molecule ThunderSTORM, SOFI and the ground truth in a movie with high emitters' density. Upper row: (a) Ground truth: high resolution image of simulated subwavelength structure. (b) Diffraction limited image created by summing the 1000 frames of high the emitters' density movie. (c) A single diffraction limited frame from the raw data of the high density movie. Clearly most emitters are active simultaneously and their far-fields are overlapping. (d) Smoothed ThunderSTORM reconstruction from a low emitters' density movie with 4000 frames. Comparing visually, the ThunderSTORM image is very similar to the ground truth. Middle row: (e) Smoothed ThunderSTORM, recovered from the movie of 1000 frames of high emitters' density. (f) Correlation SOFI (zero time-lag), recovered from the movie of 1000 frames of high emitters' density. (g) 4th order SOFI (in absolute value, zero time-lag), recovered from the movie of 1000 frames of high emitters' density. (h) SPARCOM, recovered from 200 frames of the high density movie. Lower row: SPARCOM reconstructions from varying number of frames. All reconstructions were performed from the high emitters' density movie. (i) SPARCOM, recovered from 100 frames. (j) SPARCOM, recovered from 500 frames. (k) SPARCOM, recovered from 600 frames. (l) SPARCOM, recovered from 1000 frames. All SPARCOM reconstructions were smoothed with the same kernel as in (d) and (e). Comparing panels (d) with (h)-(l) visually, it can be seen that the spatial resolution of SPARCOM is similar to that of single molecule ThunderSTORM recovered from a low density movie, in which clear separation of emitters in each frame is possible. The SPARCOM resolution is comparable to that of the ground truth image, performed only with 200 frames instead of 4000 frames as in the low density ThunderSTORM recovery, leading to a 20 times faster acquisition rate. Additionally, panel (b) shows the diffraction limited features with out-of-focus features (thicker objects). Correlation SOFI (panel (f)) did not eliminate this interference, while all SPARCOM and ThunderSTORM images successfully reject these features.

Download Full Size | PDF

 

Fig. 11 Comparative simulation intensity profiles taken from Fig. 10. Normalized cross-sections along the solid line (upper) and the dashed line (lower) of Fig. 10, comparing the ground truth (dash blue, Fig. 10(a)), diffraction-limited image (dash-dot green, Fig. 10(b)), smoothed ThunderSTORM reconstructed from the low density movie of 4000 frames (solid dot purple, Fig. 10(d)), smoothed ThunderSTORM reconstructed from the high density movie of 1000 frames (solid dot cyan, Fig. 10(e)), 4th order cumulant SOFI reconstructed from the high density movie of 1000 frames, absolute value, zero time-lag (dot black, Fig. 10(g)), and SPARCOM (solid red, Fig. 10(h)), reconstructed from 200 frames of the high density movie. We omit cross-sections from SPARCOM reconstructions performed with more frames, since they are similar to the SPARCOM cross-section from 200 frames. Considering both panels, it is clear that SPARCOM manages to recover the two adjacent sub-diffraction features with comparable width to the ground truth and the ThunderSTORM reconstruction from the low-density movie, while SOFI and ThunderSTORM recovered from the high density movie do not recover these features correctly.

Download Full Size | PDF

Figure 10(d) corresponds to a smoothed ThunderSTORM reconstruction from a low-density movie of 4000 frames, with the same number of emitters. On the other hand, Fig. 10(e) shows a smoothed ThunderSTORM reconstruction from the high-density movie, while Figs. 10(f) and 10(g) correspond to the cross-correlation and 4th order SOFI for the zero time-lag. ThunderSTORM images were constructed using the ThunderSTORM plugin for the Image-J software and smoothed with the known PSF (the parameters of ThunderSTORM are given in the Methods section in the main paper). Figure 10(h) shows the SPARCOM reconstruction over a 512 X 512 grid when used with only 200 frames and Figs 10(i)-10(l) illustrate SPARCOM reconstructions with 100, 500, 600 and 1000 frames, respectively. The images are smoothed with the same kernel as in Figs. 10(d) and 10(e).

Considering Fig. 10, we arrive at several conclusions. First, Figs. 10(i)-10(l) show that as we increase the number of frames used, the SPARCOM reconstruction improves, since the correlation estimation improves. Second, the high-density ThunderSTORM image (Fig. 10(e)) seems grainy and inconsistent in the sense that it cannot resolve the interior of the sub-diffraction features. Additionally, it is clear that single molecule ThunderSTORM requires many more frames of low density emissions to fully recover the sub-wavelength features (Fig. 10(d)), compared with SPARCOM. In this example, SPARCOM improves the temporal resolution over single molecule localization by requiring 20 times less (high density) frames (Fig. 10(h)).

We also performed an intensity cross-sections comparison to further support our conclusions, as can be seen in Fig. 11. The upper panel and the lower panel correspond to the solid and dashed yellow lines drawn in the panels of Fig. 10, respectively. The dash dot green lines in Fig. 11 correspond to the diffraction limit profiles, the dashed blue lines to the ground truth, solid dot purple to the low density smoothed ThunderSTORM image, solid dot cyan to the high density smoothed ThunderSTORM image, dot black to 4th order SOFI (zero time-lag, absolute value, 1000 frames) and last, solid red to the SPARCOM reconstruction (200 frames). Again, it is clear that SPARCOM correctly detects the true support of the filaments, while high density single molecule localization fails at detecting the sub-wavelength features. The SOFI reconstruction detects the filaments, but with degraded resolution compared to SPARCOM and single molecule localization.

4. Additional experimental results and comparisons

4.1 Comparison with ThunderSTORM, SOFI and FALCON

In this section we provide additional comparisons between single molecule ThunderSTORM, SOFI, FALCON and SPARCOM on the experimental dataset presented in Fig. 4 in the main paper.

Figure 12(a) depicts the diffraction limited image, created by summing all the frames in the movie, while Fig. 12(b) shows a single diffraction limited image from the 500 frames movie. The density of simultaneously activated emitters is such that some areas contain many overlapping emitters, while in other areas well isolated emitters can be found. Figure 12(c) corresponds to 4th order SOFI reconstruction (zero time-lag, absolute value, 500 frames), while Fig. 12(d) illustrates single molecule ThunderSTORM recovery (500 frames). Figures 12(e)-12(f) show FALCON images for the 500 and 50 frames movies, respectively. Finally, Figs. 12(g)-12(h) present SPARCOM recoveries from the 500 and 50 frame movies, respectively. All ThunderSTORM, FALCON and SPARCOM images were smoothed with the same kernel, for consistency.

 

Fig. 12 Experimental results. Panels (a)-(c) illustrate a diffraction-limited microtubule image, a single noisy frame and 4th order SOFI (absolute value, zero time-lag) recovery from a 500 frames movie, respectively. Judging visually, FALCON (500 frames, panel e and 50 frames, panel f) and SPARCOM (500 frames panel g, and 50 frames, panel h) depict a clearer image than ThunderSTORM (500 frames, panel d). Zoomed in red boxes show that FALCON requires hundreds of frames, while SPARCOM requires only tens of frames to resolve the filaments. (i, j) Intensity cross-sections (normalized) taken along the solid (i) and dashed (j) yellow lines. Panel (i) shows that in regions of high emitters' density, SPARCOM achieves better spatial resolution than the other methods, depicting the filaments' bifurcation clearly (FALCON missed this bifurcation). On the other hand, panel (j) shows that in areas of low emitters' density, the spatial resolution of SPARCOM is similar to that of single molecule localization and FALCON.

Download Full Size | PDF

Considering panel c, it is clear that the spatial resolution of SOFI does not match the resolution of the other reconstruction methods (panels e-h). Panel d (ThunderSTORM) demonstrates that in areas of closely spaced microtubules, the reconstruction is unclear, with many false localizations. In these regions, single molecule ThunderSTORM does not produce a clear image (because the density of emitters is too high), while SPARCOM reconstructions seem clearer. The zoomed area marked by the red box illustrates a region of crossing microtubules. Comparing this region in all the recoveries, we can see that single molecule ThunderSTORM does not manage to resolve the filaments, while FALCON requires several hundreds of frame for good detection. On the other hand, SPARCOM resolves the filaments with only several hundreds of frames. This observation substantiates the ability of SPARCOM to better resolve sub-diffraction features than the other methods considered here. The reconstruction parameters of SPARCOM and ThunderSTORM are as described in the Methods section in the main paper, while the FALCON parameters are as described previously.

To further support the conclusions drawn from visual examination of the panels of Fig. 12, we compare intensity profiles (Fig. 12(i), 12(j)) along the solid and dashed lines in panels a-h. Considering Fig. 12(i), in areas of high emitters' density, SPARCOM clearly depicts the contours of two filaments (peak distance of $127nm$), while SOFI and ThunderSTORM show a degraded resolution profile. In this case, FALCON did not manage to separate the two microtubules, and mistakenly detected a single one between the two actual microtubules. On the other hand, Fig. 12(j) shows that when considering low emitters' density (a well isolated filament), the spatial resolution of SPARCOM is similar to that of ThunderSTORM and FALCON, and in both cases, it is better than the spatial resolution of SOFI.

4.2 Comparison with ThunderSTORM, SOFI and FALCON

In this section we provide additional comparisons to the multi-emitter fitting algorithm, present in the ThunderSTORM Image-J plug-in. All ThunderSTORM reconstructions were made using the multi-emitter fitting option (maximum of 10 molecules per fitting region) and using a maximum likelihood estimation (the remaining parameters are as described in the Methods section). All reconstructions were smoothed with the same previously used Gaussian kernel. The first three panels of Fig. 13 correspond to the experimental dataset presented in the main paper (Fig. 4), while the three lower panels correspond to the simulated dataset presented in the main paper (Figs. 2 and 3).

 

Fig. 13 Multi-emitter fitting comparisons on both simulated and experimental data sets. (a), (d) Diffraction limited images of the experimental and simulated movies, respectively. (b), (e) Multi-emitter ThunderSTORM fitting using movies of 500 and 361 frames of experimental and simulated data sets, respectively. (c), (f) Multi-emitter ThunderSTORM fitting using movies of 50 and 60 frames of experimental and simulated data sets, respectively. Red boxes indicated enlarged areas in lower right corner of each panel. As the number of frames is reduced, while preserving the total number of emitters, reconstruction quality degrades.

Download Full Size | PDF

Observing Fig. 13, we deduce the following. As the number of frames decreases while preserving the total number of emitters, the reconstruction quality of the maximum likelihood, multi-emitter ThunderSTORM fit degrades. For example, considering the enlarged areas in the red boxes at the lower right corner of panel f, we observe that the filaments are now broken completely, and no clear depiction of them is evident. Furthermore, even when using lower density (more frames) as is evident from panels b and e, we observe many false positive detections, which result in a noisy reconstruction which does not manage to clearly resolve the filaments (for example, the enlarged box in panel e). In contrast to Multi-emitter ThunderSTORM, these filaments are clearly separated using SPARCOM, as shown in Figs. 2-4 in the main paper.

All of the reconstructions presented in Figs. 7-13 support the same conclusions drawn in the main paper. When temporal resolution is critical, increasing the emitters' density is beneficial and leads to a shorter acquisition time than the time required for a single molecule localization experiment. In the case of high emitter density and for the same measurements, SPARCOM is preferable over state-of-the-art methods and achieves a similar spatial resolution as low density single molecule localization, while requiring much fewer frames, thus improving the spatio-temporal resolution of the super-resolved image.

4.3 Comparison with 3B

We also perform a qualitative comparison between SPARCOM and 3B, which is considered as one of the leading methods for high density super-resolution microscopy. We consider a very limited field of view, due to the highly demanding computational cost of the 3B method (implemented as a stand-alone Image-J plug-in). Figure 14 presents reconstruction results for the simulated and experimental datasets given in Figs. 2 and 4 in the main paper, respectively. Upper row show shows the diffraction limited image (Fig. 14(a)) and 3B and SPARCOM reconstructions from the movie of 361 high density frames (Fig. 14(b) and Fig. 14(c), respectively). We also perform the same comparison for the even denser movie of 60 frames (Fig. 14(d) for 3B and Fig. 14(e) for SPARCOM). Lower row presents respective results for a selected region from the experimental dataset of Fig. 4 of the main paper, where Fig. 14(f) shows the diffraction limited image and Figs. 14(g) and S9h show 3B and SPARCOM reconstructions from a movie of 500 frames, respectively. Lastly, Figs. 14(i) and 14(j) show 3B and SPARCOM reconstructions from a denser movie of 50 frames, respectively.

 

Fig. 14 Qualitative comparison between SPARCOM and 3B for high emitters’ density data sets. (a), (f) are diffraction limited images of a selected area from the simulated and experimental data sets presented in Figs. 2 and 4 in the main paper, respectively. Panels (b) and (c) show 3B and SPARCOM reconstructions using 361 frames, respectively, while panels. (d) and (e) show 3B and SPARCOM reconstructions for 60 frames, respectively (by summing every 6 consecutive frames in the 361 frames movie). Panels (g) and (h) present 3B and SPARCOM reconstructions for 500 frames, while panels (i) and (j) compare 3B and SPARCOM reconstructions of the same object from a movie of 50 frames (by summing every 10 consecutive frames in the 500 frames movie). The 3B reconstructions seem very aggregated and broken, compared with SPARCOM, in all cases.

Download Full Size | PDF

Observing the panels of Fig. 14, we see that the 3B reconstruction is affected dramatically by the reduction in the number of frames (and increment in the emitters’ density). For example, considering the highlighted regions by the white rectangles in Figs. 14(e) and 14(j), SPARCOM reconstructions manage to separate clearly the overlapping microtubules, while the 3B reconstructions of the same areas seem aggregated and inconsistent (panels b and g), especially when using only tens of frames (panels d and i). Furthermore, panel d shows that for 60 frames, 3B does not separate the filaments well, while in panel e, SPARCOM clearly performs better.

We can conclude therefore, that SPARCOM performs better than 3B in terms of super-resolution image reconstruction, with only tens of frames – leading to an improvement in the temporal resolution and total acquisition time.

5. Metrics

5.1 Minimal separation distance

In the following, we present additional metrics by which we objectively assess the performance of SPARCOM, compared with the other methods considered in the paper. Since our main goal is to achieve high spatial resolution with a minimal number of frames, we focus on resolution criteria as a function of molecular density and minimal separation distance.

Figure 15 demonstrates the ability of SPARCOM to separate close by emitters. To show that, two simulated datasets of 50 frames each are generated. The first corresponds to two emitters and the second to four (diffraction limited images are given in panel (a) and (c) respectively). Corresponding SPARCOM reconstructions alongside the emitters' locations are given in panels (b) and (d). The emission wavelength of this simulation is ${\lambda }_{em}=723nm$, NA is 1.4 and each pixel size is $100nm$ on the $64X64$ pixels movie. Clearly, SPARCOM exhibits a 3.44-times resolution increase over the diffraction limit ($258.11nm$), by separating the two emitters situated ${\lambda }_{em}/10$, or $75nm$ apart (panel (b)). When four emitters are simulated as placed on a perfect diamond shape (panel (d)), one can observe that the vertical and horizontal distances are $100nm$, but the closest distance is $70.7nm$.

 

Fig. 15 Simulation of the minimal separation distance for SPARCOM. (a) and (c) show diffraction limited spots created by summing 50 frames of two and four simulated fluctuating emitters, respectively. Purple dots represent true positions of the emitters. (b) and (d) illustrate SPARCOM recoveries of the two and four emitters, respectively, alongside their true locations, marked by the purple dots. Clear separation is achieved even when the emitters are $75nm$ apart, as shown in panel (b).

Download Full Size | PDF

5.2 Hit-rate and false-alarm rate vs. labeling density

To further quantify the ability of SPAROCM in scenarios of increasing labeling density, we perform the following experiment of hit-rate and false-alarm rate vs. labeling density. We also compare SPARCOM against frame-by-frame sparsity-based reconstruction using the FISTA algorithm to show the advantage of further exploiting the uncorrelated blinking of the emitters. Namely, we perform sparse recovery for each diffraction limited frame in a movie and then accumulate all reconstructed images without exploiting statistical information. We perform this test for two datasets of synthetic data. Figures 16(a) and 16(b) show an overlay of diffraction limited spots (grey), SPARCOM (purple) and frame-by-frame sparsity-based reconstruction (green) reconstructions for two different emitters' densities (1000 frames in each movie, $64X64$ pixels and reconstruction is on a $512X512$ pixels grid). Emission wavelength is $800nm$. It is evident that the SPARCOM reconstructions, though having some artifacts, are much clearer than the frame-by-frame formulation and capture far more emitters. These reconstructions support our conclusion that using sparse recovery in the correlations domain, while using the uncorrelated blinking of the emitters is much more beneficial than using sparse recovery alone on each image.

 

Fig. 16 Reconstruction comparison between SPARCOM and a sparsity-based frame-by-frame recovery. Panels (a) and (b) show an overlay of diffraction limited spots (grey), SPARCOM (purple) and frame-by-frame sparsity-based reconstruction (green) reconstructions for two different emitters' densities. Each movie contains 1000 frames. Clearly, SPARCOM exhibits more details and better localization of the emitters than by performing a sparsity-based frame-by-frame recovery. SPARCOM recoveries were performed using $\lambda =0.0001$, while frame-by-frame recoveries were performed using $\lambda =0.1$, such that each method produces the best possible result.

Download Full Size | PDF

To additionally ascertain the ability of SPARCOM to handle highly dense datasets, we study the resolution of SPARCOM as a function of the density of excited emitters, and compare it to single molecule localization (implemented using single molecule ThunderSTORM with parameters as indicated in the Methods section) under the same densities. To this end, we carry out a comparative study, in a series of simulations under varying emitter densities, as shown in Fig. 17. We perform this test to empirically estimate SPARCOM performance in terms of correct support recovery of emitters (that is, correct identification of their locations), by recovering images made up of point emitters in a series of scenarios ranging from very low to high emitter densities. We simulate a series of movies of randomly distributed emitters over a field of view (FOV) of $10.24\mu m^2$, represented by a grid of $512X512$ pixels and imaged with a numerical aperture of 1.4 and emission wavelength of $800nm$. Each movie contains 1000 frames and each frame is of size $64X64$ pixels (8-fold compression). The different movies represent different densities of emitters, and we perform SPARCOM recovery for each of them, under a parametric sweep of different values of the regularization parameter $\lambda$. For comparison, we also consider ThunderSTORM reconstructions for the same values of emitters’ density. Figures 17(a)- 17(c) show three examples of an overlay of diffraction limited spots (grey), SPARCOM (green dots) and ThunderSTORM (purple dots) reconstructions for different emitters' densities. In all cases, no smoothing is performed since the images are composed of random point emitters (i.e. they do not describe a specific feature).

 

Fig. 17 Qualitative comparison between the ability of SPARCOM and ThunderSTORM to correctly identify positions of emitters with increasing emitters' density. (a)-(c) Overlay of diffraction limited spots (grey), SPARCOM (green) and ThunderSTORM (purple) reconstructions for different emitters' densities. White dots indicate that both reconstructions overlap. All images are presented in logarithmic scale, for clarity. Comparisons are made on a linear scale. (d)-(e) Correct support identification (“hit rate”) and erroneous support identification (“false alarm rate”), respectively. (d)-(e) show that when high emitters' density is used, SPARCOM identifies the correct emitters' locations substantially better than single molecule fit, while preserving a low false-alarm rate.

Download Full Size | PDF

Figures 17(d) and 17(e) show a comparative study for the different recovery methods. Figure 17(d) depicts the recovery performance of SPARCOM and ThunderSTORM (single emitter fit) as the emitters' density is increased. The performance is measured by the “hit rate”, that is, out of the total number of emitters present in the data, we measure the percentage of correctly identified emitters in terms of their location (correct identification is considered if a recovered emitter is found near an emitter in the ground truth up to a shift of one pixel in the $512X512$ grid, in any direction). Figure 17(e) illustrates the number of erroneously identified locations both by SPARCOM and ThunderSTORM under increasing emitters' density, referred to as “false alarm rate”. From these panels, we deduce that as the density increases (many emitters fluoresce simultaneously), SPARCOM performs much better than ThunderSTORM. For very low emitter densities, these panels show that the performance of SPARCOM is similar to ThunderSTORM in terms of support identification.

Furthermore, The analysis of the reconstruction performance is fairly robust to the choice of the parameter $\lambda$ ($\lambda =0.0001, 0.0005$ and $\lambda =0.001$, which are different by an order of magnitude). Of course, as we increase the density, at some point all methods fail, but SPARCOM fails at much higher densities than single emitter fit. For example, considering Fig. 17(c), the hit-rate of SPARCOM (blue line) for emitters' densities of $1\mu m^{-2}, 3\mu m^{-2}$, and $7\mu m^{-2}$ are $95.4\text{\%},92.1\text{\%}$ and $85.9\text{\%}$ for SPARCOM and $74.8\text{\%},73.5\text{\%}$ and $34.4\text{\%}$ for single molecule ThunderSTORM, respectively. In addition, the false recovery rate of all methods is roughly similar. We conclude that SPARCOM considerably outperforms single molecule fitting when high emitters' density is used. This implies that a great reduction in the number of required exposures is possible for the recovery of a super-resolved image – leading to a substantial improvement in the temporal resolution of the recovery process – while maintaining the highest resolution of single molecule localization (which is obtained at very low densities of excited emitters).

5.3 Fourier ring correlation resolution enhancement vs. labeling density

Next we compare Fourier Ring Correlation (FRC) curves against varying emitter density and SNR. The first is a comparison between state-of-the-art methods MUSICAL, SRRF and SPARCOM, over the same numerical experiment of increasing labeling densities. Each curve in Fig. 18 represents the resolution enhancement for each movie of increasing labeling density, for each of the methods. That is, each point on each curve represents the FRC enhancement, calculated between the ground truth image and the corresponding reconstruction, for the different movies of increasing labeling densities. Observing Fig. 18, the SPARCOM curve represents the highest resolution (lower is better), especially as the emitters' density increases. This comparison confirms that SPAROM achieves better performance in terms of spatial resolution, especially when increasing the emitters' density, which in turn leads to an improvement in the overall temporal resolution of the reconstruction.

 

Fig. 18 Resolution enhancement versus increasing labeling density. Each curve represents the resolution enhancement achieved by performing FRC analysis over movies of increasing labeling densities for MUSICAL, SRRF and SPARCOM. Each point on each curve is the result of FRC analysis between the ground truth (simulated) image and the reconstruction for each method. This figure clearly shows that the resolution gain of SPARCOM is better than SRRF and MUSICAL, especially as the labeling density increases.

Download Full Size | PDF

5.4 Fourier ring correlation resolution enhancement vs. SNR

The second FRC curve tests the robustness of SPARCOM against noise. We use the simulation of two fluctuating emitters (100 frames in the movie) presented in Fig. 15 with varying SNR values, where the SNR is computed as

$$SN{R}_{dB}=10{\log }_{10}\frac{Powe{r}_{signal}}{Powe{r}_{noise}}.$$

Here, $Powe{r}_{signal}$ is the average variance of each frame in the input movie prior to the addition of noise, and $Powe{r}_{noise}$ is the variance of the additive white Gaussian noise we added to the movie. For each value of SNR, we calculated the FRC resolution enhancement between the reconstruction and the ground truth image of the two close-by emitters (separated by a distance of $75nm$). Figure 19 shows the corresponding FRC vs. SNR curves for SPARCOM, SRRF, MUSICAL and ThunderSTORM (multi-emitter fit, maximum likelihood estimation). Clearly, we can observe that SPARCOM achieves higher resolution for all SNR values (lower is better, corresponding to higher spatial resolution), better than all considered methods. All methods exhibit a resolution increase around $-25dB$.

 

Fig. 19 Resolution enhancement versus SNR. Each curve represents the resolution gain achieved by performing FRC analysis over movies of increasing SNR values (the lower the curve, the higher the spatial resolution). Each point on a curve is the result of FRC analysis between the ground truth (simulated) image and the reconstruction for each method. This figure clearly shows that the resolution gain of SPARCOM is better than the other methods, for all SNR values considered here.

Download Full Size | PDF

6. Methods

6.1 Simulation setup of Figs. 2, 3 and 5(a),5(b)

The results presented in the aforementioned figures are based on a simulation of microtubules from a freely available movie containing 361 noisy frames, each of size $64X64$ pixels and a pixel size of 100²nm. The imaging wavelength is 723nm and the numerical aperture (NA) is 1.4. Additionally, we create a very high-density movie of 60 frames by summing every six frames from the 361 frames movie. Since this is a simulated dataset, the PSF is a-priori known and does not need to be estimated. In all the following examples (and also in the examples provided in the Appendix), the value of the regularization parameter λ in the SPARCOM algorithm is determined empirically. The SPARCOM implementation is written in MATLAB (The MathWorks, Inc.). In this simulation, SPARCOM was performed with $\lambda =0.12,0.13$ for the movies of 60 and 361 frames, respectively, using 400 iterations and on patches of size $16X16$ pixels. Reconstruction was performed on a grid of size $512X512$ pixels. SRRF was performed for both movies with a ring radius of 0.5, radiality magnification of 8 ($512X512$ ) pixels and 6 axes in the ring, using the TRAC2 (correlations) method (using the TRAC4 method resulted in a grainy image, similar to SOFI images due to the dynamic range expansion phenomena). MUSICAL was performed with thresholds of −1 and 0 for the 60 frames and 361 frames movies, respectively and $\alpha =4$ in both cases. Reconstruction was performed on a grid of $512X512$ pixels. Lastly, ThunderSTORM was performed for both movies with a B-spline wavelet filter of order 3 and scale 2, approximate localization with 8 neighbors, subpixel localization with an integrated Gaussian using the weighted least squares method and a fit radius and initial of 3 and 1.6, respectively (additional ThunderSTORM recoveries with maximum likelihood estimation are given in Fig. 13). Image was reconstructed on a grid of $512X512$ pixels using the average shifted histograms rendering method.

6.2 Experimental setup of Figs. 4 and 5(c)

We used a freely available experimental dataset of microtubules over 500 frames, taken with a frame-rate of 25Hz. We also create a highly dense movie of 50 frames, by summing every 10 consecutive frames. Each frame is of size $64X64$ pixels with a pixel size of 100nm. The numerical aperture is 1.3 and the imaging wavelength is 690nm. Based on these parameters, we construct the PSF using the freely available PSF generator. SPARCOM was performed with $\lambda =250$ for both movies, using 2000 and 200 iterations for the 500 frames and 50 frames movies, respectively and on patches of size $8X8$ pixels. Reconstruction was performed on a grid of size $512X512$ pixels. SRRF was performed for both movies with a ring radius of 0.5, radiality magnification of 8 ($512X512$) pixels and 6 axes in the ring, using the TRAC2 (correlations) method. MUSICAL was performed with thresholds of −1 and 0 for the 50 frames and 500 frames movies, respectively and $\alpha =4$ in both cases. Reconstruction was performed on a grid of $512X512$ pixels. Lastly, ThunderSTORM was performed for both movies with a B-spline wavelet filter of order 3 and scale 2, approximate localization with 8 neighbours, subpixel localization with an integrated Gaussian using the weighted least squares method and a fit radius and initial of 3 and 1.6, respectively. Image was reconstructed on a grid of $512X512$ pixels using the average shifted histograms rendering method.

All SPARCOM reconstructions were performed using an 8Gbyte RAM PC, using a single core implementation (Intel I5-4430 CPU, 3 GHz clock). The typical runtime for a SPARCOM reconstruction of a patch of $8X8$ pixels to a grid of $64X64$ (8 times denser) pixels is 4-5 seconds for 2000 Iterations and 1000 frames (typically ~0.3-0.4 milliseconds per a single iteration). For a patch of $32X32$ pixels to a grid of $256X256$ (8 times denser) pixels is 27-28 seconds for 2000 Iterations and 1000 frames (typically ~3-4 milliseconds per a single iteration). Larger image sizes can be analyzed as a single image, or in a patch-based manner (which runs faster). We implemented SPARCOM as a single core algorithm, but a patch-based analysis can be efficiently parallelized.

Funding

European Union's Horizon 2020 research and innovation program (646804-ERC-COG-BNYQ).

Acknowledgements

We thank Prof. Shimon Weiss and Xiyu Yi for fruitful discussions on super-resolution fluorescence microscopy and SOFI.

References and links

1. M. Born and E. Wolf, Principles of Optics: Electromagnetic Theory of Propagation, Interference and Diffraction of Light (Elsevier, 2013).

2. S. W. Hell and J. Wichmann, “Breaking the diffraction resolution limit by stimulated emission: stimulated-emission-depletion fluorescence microscopy,” Opt. Lett. 19(11), 780–782 (1994). [CrossRef]   [PubMed]  

3. E. Betzig, G. H. Patterson, R. Sougrat, O. W. Lindwasser, S. Olenych, J. S. Bonifacino, M. W. Davidson, J. Lippincott-Schwartz, and H. F. Hess, “Imaging intracellular fluorescent proteins at nanometer resolution,” Science 313, 1642–1645 (2006).

4. M. J. Rust, M. Bates, and X. Zhuang, “Sub-diffraction-limit imaging by stochastic optical reconstruction microscopy (STORM),” Nat. Methods 3(10), 793–796 (2006). [CrossRef]   [PubMed]  

5. M. Heilemann, S. van de Linde, M. Schüttpelz, R. Kasper, B. Seefeldt, A. Mukherjee, P. Tinnefeld, and M. Sauer, “Subdiffraction-resolution fluorescence imaging with conventional fluorescent probes,” Angew. Chem. Int. Ed. Engl. 47(33), 6172–6176 (2008). [CrossRef]   [PubMed]  

6. M. G. Gustafsson, “Surpassing the lateral resolution limit by a factor of two using structured illumination microscopy,” J. Microsc. 198(2), 82–87 (2000). [CrossRef]   [PubMed]  

7. P. Kner, B. B. Chhun, E. R. Griffis, L. Winoto, and M. G. L. Gustafsson, “Super-resolution video microscopy of live cells by structured illumination,” Nat. Methods 6(5), 339–342 (2009). [CrossRef]   [PubMed]  

8. L. Schermelleh, R. Heintzmann, and H. Leonhardt, “A guide to super-resolution fluorescence microscopy,” J. Cell Biol. 190(2), 165–175 (2010). [CrossRef]   [PubMed]  

9. B. R. Rankin, G. Moneron, C. A. Wurm, J. C. Nelson, A. Walter, D. Schwarzer, J. Schroeder, D. A. Colón-Ramos, and S. W. Hell, “Nanoscopy in a living multicellular organism expressing GFP,” Biophys. J. 100(12), L63–L65 (2011). [CrossRef]   [PubMed]  

10. V. Westphal, S. O. Rizzoli, M. A. Lauterbach, D. Kamin, R. Jahn, and S. W. Hell, “Video-Rate Far-Field Optical Nanoscopy Dissects Synaptic Vesicle Movement,” Science 320(5873), 246–249 (2008). [CrossRef]   [PubMed]  

11. T. Dertinger, R. Colyer, G. Iyer, S. Weiss, and J. Enderlein, “Fast, background-free, 3D super-resolution optical fluctuation imaging (SOFI),” Proc. Natl. Acad. Sci. U.S.A. 106(52), 22287–22292 (2009). [CrossRef]   [PubMed]  

12. J. M. Mendel, “Tutorial on higher-order statistics (spectra) in signal processing and system theory: theoretical results and some applications,” Proc. IEEE 79(3), 278–305 (1991). [CrossRef]  

13. T. Dertinger, R. Colyer, R. Vogel, J. Enderlein, and S. Weiss, “Achieving increased resolution and more pixels with Superresolution Optical Fluctuation Imaging (SOFI),” Opt. Express 18(18), 18875–18885 (2010). [CrossRef]   [PubMed]  

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

15. 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(1), 4577 (2015). [CrossRef]   [PubMed]  

16. S. Cox, E. Rosten, J. Monypenny, T. Jovanovic-Talisman, D. T. Burnette, J. Lippincott-Schwartz, G. E. Jones, and R. Heintzmann, “Bayesian localization microscopy reveals nanoscale podosome dynamics,” Nat. Methods 9(2), 195–200 (2012). [CrossRef]   [PubMed]  

17. H. Deschout, T. Lukes, A. Sharipov, D. Szlag, L. Feletti, W. Vandenberg, P. Dedecker, J. Hofkens, M. Leutenegger, T. Lasser, and A. Radenovic, “Complementarity of PALM and SOFI for super-resolution live-cell imaging of focal adhesions,” Nat. Commun. 7, 13693 (2016). [CrossRef]   [PubMed]  

18. N. Gustafsson, S. Culley, G. Ashdown, D. M. Owen, P. M. Pereira, and R. Henriques, “Fast live-cell conventional fluorophore nanoscopy with ImageJ through super-resolution radial fluctuations,” Nat. Commun. 7, 12471 (2016). [CrossRef]   [PubMed]  

19. K. Agarwal and R. Macháň, “Multiple signal classification algorithm for super-resolution fluorescence microscopy,” Nat. Commun. 7, 13752 (2016). [CrossRef]   [PubMed]  

20. Y. C. Eldar, Sampling Theory: Beyond Bandlimited Systems (Cambridge University, 2015).

21. O. Solomon, M. Mutzafi, X. Yi, S. Weiss, Y. C. Eldar, and M. Segev, “Sparsity-based super-resolution optical fluctuation imaging,” in CLEO: Applications and Technology, 2016.

22. P. Pal and P. P. Vaidyanathan, “Pushing the Limits of Sparse Support Recovery Using Correlation Information,” IEEE Trans. Signal Process. 63(3), 711–726 (2015). [CrossRef]  

23. S. W. Ellingson, “Sensitivity of antenna arrays for long-wavelength radio astronomy,” IEEE Trans. Antenn. Propag. 59(6), 1855–1863 (2011). [CrossRef]  

24. Y. C. Eldar and G. Kutyniok, Compressed Sensing: Theory and Applications (Cambridge University, 2012).

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

26. E. Candès, J. Romberg, and T. Tao, “Robust uncertainty principles: Exact signal reconstruction from highly incomplete frequency information,” IEEE Trans. Inf. Theory 52(2), 489–509 (2006). [CrossRef]  

27. E. J. Candès, J. K. Romberg, and T. Tao, “Stable signal recovery from incomplete and inaccurate measurements,” Commun. Pure Appl. Math. 59(8), 1207–1223 (2006). [CrossRef]  

28. Y. Rivenson, A. Stern, and J. Rosen, “Compressive multiple view projection incoherent holography,” Opt. Express 19(7), 6109–6118 (2011). [CrossRef]   [PubMed]  

29. Y. Rivenson, A. Stern, and B. Javidi, “Overview of compressive sensing techniques applied in holography [Invited],” Appl. Opt. 52(1), A423–A432 (2013). [CrossRef]   [PubMed]  

30. N. Wagner, Y. C. Eldar, and Z. Friedman, “Compressed beamforming in ultrasound imaging,” IEEE Trans. Signal Process. 60(9), 4643–4657 (2012). [CrossRef]  

31. T. Chernyakova and Y. Eldar, “Fourier-domain beamforming: the path to compressed ultrasound imaging,” IEEE Trans. Ultrason. Ferroelectr. Freq. Control 61(8), 1252–1267 (2014). [CrossRef]   [PubMed]  

32. M. Mishali and Y. C. Eldar, “From theory to practice: Sub-Nyquist sampling of sparse wideband analog signals,” IEEE J. Sel. Top. Signal Process. 4(2), 375–391 (2010). [CrossRef]  

33. M. Mishali and Y. C. Eldar, “Blind multiband signal reconstruction: Compressed sensing for analog signals,” IEEE Trans. Signal Process. 57(3), 993–1009 (2009). [CrossRef]  

34. Z. Tian and G. B. Giannakis, “Compressed sensing for wideband cognitive radios,” ICASSP, IEEE Int. Conf. Acoust. Speech Signal Process. Proc. 4, 1357–1360, 2007. [CrossRef]  

35. R. Baraniuk and P. Steeghs, “Compressive radar imaging,” in IEEE Natl. Radar Conf. Proc.128–133, 2007.

36. O. Bar-Ilan and Y. C. Eldar, “Sub-Nyquist Radar via Doppler Focusing,” IEEE Trans. Signal Process. 62(7), 1796–1811 (2014). [CrossRef]  

37. M. Lustig, D. Donoho, and J. M. Pauly, “Sparse MRI: The application of compressed sensing for rapid MR imaging,” Magn. Reson. Med. 58(6), 1182–1195 (2007). [CrossRef]   [PubMed]  

38. L. Weizman, Y. C. Eldar, and D. Ben Bashat, “Reference-based MRI,” Med. Phys. 43(10), 5357–5369 (2016). [CrossRef]   [PubMed]  

39. V. I. Morgenshtern and E. J. Candès, “Stable Super-Resolution of Positive Sources : the Discrete Setup,” SIAM J. Imaging Sci. Imaging Sci. 9(2), 1–35 (2015).

40. S. Gazit, A. Szameit, Y. C. Eldar, and M. Segev, “Super-resolution and reconstruction of sparse sub-wavelength images,” Opt. Express 17(26), 23920–23946 (2009). [CrossRef]   [PubMed]  

41. M. Mutzafi, Y. Shechtman, Y. C. Eldar, O. Cohen, and M. Segev, “Sparsity-based Ankylography for Recovering 3D molecular structures from single-shot 2D scattered light intensity,” Nat. Commun. 6(1), 7950 (2015). [CrossRef]   [PubMed]  

42. A. Szameit, Y. Shechtman, E. Osherovich, E. Bullkich, P. Sidorenko, H. Dana, S. Steiner, E. B. Kley, S. Gazit, T. Cohen-Hyams, S. Shoham, M. Zibulevsky, I. Yavneh, Y. C. Eldar, O. Cohen, and M. Segev, “Sparsity-based single-shot subwavelength coherent diffractive imaging,” Nat. Mater. 11(5), 455–459 (2012). [CrossRef]   [PubMed]  

43. P. Sidorenko, O. Kfir, Y. Shechtman, A. Fleischer, Y. C. Eldar, M. Segev, and O. Cohen, “Sparsity-based super-resolved coherent diffraction imaging of one-dimensional objects,” Nat. Commun. 6(1), 8209 (2015). [CrossRef]   [PubMed]  

44. Y. Shechtman, Y. C. Eldar, O. Cohen, and M. Segev, “Efficient coherent diffractive imaging for sparsely varying objects,” Opt. Express 21(5), 6327–6338 (2013). [CrossRef]   [PubMed]  

45. H. Jiang, G. Huang, and P. Wilford, “Multi-view in lensless compressive imaging,” APSIPA Trans. Signal Inf. Process 3, e152014.

46. M. L. Moravec, J. K. Romberg, and R. G. Baraniuk, “Compressive phase retrieval,” Proc. SPIE 6701, 670120 (2007).

47. P. Schniter and S. Rangan, “Compressive phase retrieval via generalized approximate message passing,” IEEE Trans. Signal Process. 63(4), 1043–1055 (2015). [CrossRef]  

48. Y. Shechtman, Y. C. Eldar, A. Szameit, and M. Segev, “Sparsity based sub-wavelength imaging with partially incoherent light via quadratic compressed sensing,” Opt. Express 19(16), 14807–14822 (2011). [CrossRef]   [PubMed]  

49. Y. Shechtman, S. Gazit, A. Szameit, Y. C. Eldar, and M. Segev, “Super-resolution and reconstruction of sparse images carried by incoherent light,” Opt. Lett. 35(8), 1148–1150 (2010). [CrossRef]   [PubMed]  

50. Y. Shechtman, E. Small, Y. Lahini, M. Verbin, Y. C. Eldar, Y. Silberberg, and M. Segev, “Sparsity-based super-resolution and phase-retrieval in waveguide arrays,” Opt. Express 21(20), 24015–24024 (2013). [CrossRef]   [PubMed]  

51. E. J. Candès, X. Li, and M. Soltanolkotabi, “Phase retrieval from coded diffraction patterns,” Appl. Comput. Harmon. Anal. 39(2), 277–299 (2015). [CrossRef]  

52. Y. Shechtman, Y. C. Eldar, O. Cohen, H. N. Chapman, J. Miao, and M. Segev, “Phase Retrieval with Application to Optical Imaging: A contemporary overview,” IEEE Signal Process. Mag. 32(3), 87–109 (2015). [CrossRef]  

53. Y. Shechtman, A. Beck, and Y. C. Eldar, “GESPAR: Efficient phase retrieval of sparse signals,” IEEE Trans. Signal Process. 62(4), 928–938 (2014). [CrossRef]  

54. K. Jaganathan, Y. C. Eldar, and B. Hassibi, “STFT Phase Retrieval: Uniqueness Guarantees and Recovery Algorithms,” IEEE J. Sel. Top. Signal Process. 10(4), 770–781 (2016). [CrossRef]  

55. V. Studer, J. Bobin, M. Chahid, H. S. Mousavi, E. Candes, and M. Dahan, “Compressive fluorescence microscopy for biological and hyperspectral imaging,” Proc. Natl. Acad. Sci. U.S.A. 109(26), E1679–E1687 (2012). [CrossRef]   [PubMed]  

56. D. Oren, Y. Shechtman, M. Mutzafi, Y. C. Eldar, and M. Segev, “Sparsity-based recovery of three-photon quantum states from two-fold correlations,” Optica 3(3), 226–232 (2016). [CrossRef]  

57. D. Gross, Y. K. Liu, S. T. Flammia, S. Becker, and J. Eisert, “Quantum state tomography via compressed sensing,” Phys. Rev. Lett. 105(15), 150401 (2010). [CrossRef]   [PubMed]  

58. W. T. Liu, T. Zhang, J. Y. Liu, P. X. Chen, and J. M. Yuan, “Experimental quantum state tomography via compressed sampling,” Phys. Rev. Lett. 108(17), 170403 (2012). [CrossRef]   [PubMed]  

59. D. Cohen and Y. C. Eldar, “Sub-Nyquist Sampling for Power Spectrum Sensing in Cognitive Radios: A Unified Approach,” IEEE Trans. Signal Process. 62(15), 3897–3910 (2014). [CrossRef]  

60. T. Főrster, “10th Spiers Memorial Lecture. Transfer mechanisms of electronic excitation,” Discuss. Faraday Soc. 27(0), 7–17 (1959). [CrossRef]  

61. D. L. Dexter, “A Theory of Sensitized Luminescence in Solids,” J. Chem. Phys. 21(5), 836–850 (1953). [CrossRef]  

62. S. Boyd and L. Vandenberghe, Convex Optimization (Cambridge University, 2004).

63. E. J. Candès, M. B. Wakin, and S. P. Boyd, “Enhancing sparsity by reweighted L1 minimization,” J. Fourier Anal. Appl. 14(5), 877–905 (2008). [CrossRef]  

64. A. Beck and M. Teboulle, “A Fast Iterative Shrinkage-Thresholding Algorithm,” SIAM J. Imaging Sci. 2(1), 183–202 (2009). [CrossRef]  

65. T. Wimalajeewa, Y. C. Eldar, and P. K. Varshney, “Recovery of Sparse Matrices via Matrix Sketching,” arXiv Prepr. arXiv:1311.2448, 2013.

66. D. P. Palomar and Y. C. Eldar, Convex Optimization in Signal Processing and Communications (Cambridge University, 2010).

67. O. Solomon, Y. C. Eldar, M. Mutzafi, and M. Segev, “SPARCOM: Sparsity Based Super-Resolution Correlation Microscopy,” arXiv Prepr. arXiv: 1707.09255, 2017.

68. “Single-Molecule Localization Microscopy.” [Online]. http://bigwww.epfl.ch/smlm/datasets/index.html.

69. H. Kirshner, F. Aguet, D. Sage, and M. Unser, “3-D PSF fitting for fluorescence microscopy: Implementation and localization application,” J. Microsc. 249(1), 13–25 (2013). [CrossRef]   [PubMed]  

70. R. P. J. Nieuwenhuizen, K. A. Lidke, M. Bates, D. L. Puig, D. Grünwald, S. Stallinga, and B. Rieger, “Measuring image resolution in optical nanoscopy,” Nat. Methods 10(6), 557–562 (2013). [CrossRef]   [PubMed]  

71. A. Savitzky and M. J. E. Golay, “Smoothing and differentiation of data by simplified least squares procedures,” Anal. Chem. 36(8), 1627–1639 (1964). [CrossRef]