Sparse deconvolution of high-density super-resolution images

In wide-field super-resolution microscopy, investigating the nanoscale structure of cellular processes, and resolving fast dynamics and morphological changes in cells requires algorithms capable of working with a high-density of emissive fluorophores. Current deconvolution algorithms estimate fluorophore density by using representations of the signal that promote sparsity of the super-resolution images via an L1-norm penalty. This penalty imposes a restriction on the sum of absolute values of the estimates of emitter brightness. By implementing an L0-norm penalty – on the number of fluorophores rather than on their overall brightness – we present a penalized regression approach that can work at high-density and allows fast super-resolution imaging. We validated our approach on simulated images with densities up to 15 emitters per μm-2 and investigated total internal reflection fluorescence (TIRF) data of mitochondria in a HEK293-T cell labeled with DAKAP-Dronpa. We demonstrated super-resolution imaging of the dynamics with a resolution down to 55 nm and a 0.5 s time sampling.


Supplementary File Title Supplementary Notes
Computational strategies for SPIDER

Supplementary Table S1
Comparing the results obtained with the L0-norm penalty and the L1-norm penalty in superresolution imaging Supplementary Table S2 Calculation time of deconvolution methods

Supplementary Table S3
Figures of merit for a set of emitters with a random placement at different densities Supplementary Figure S1 Method for performance evaluation

Supplementary Figure S2
Performance benchmark of CSSTORM, FALCON and SPIDER for a PSF to pixel size ratio of 3.9

Supplementary Figure S3
Insights on the data and results of a series of three juxtaposed circles

Supplementary Figure S4
Fast SPIDER super-resolution for mitochondria in a HEK293-T cell labeled with DAKAP-Dronpa

Supplementary Figure S5
Estimation of the image resolution by Fourier Ring Correlation (FRC) on the HEK293-T cell labeled with DAKAP-Dronpa Supplementary Video S1 Raw movie data for a HEK293-T cell

Supplementary Video S2
SPIDER on an image patch showing a clear recovery of individual emitters Supplementary Software S1 SPIDER algorithm

Theory
The observed image Y can easily be more than 10 5 pixels and thus gives a vector y with more than a hundred thousand elements. The vector x is much larger: with a sub-grid oversampling of 4 by 4, it will have 1.6 million elements and the matrix C will have that many rows and columns. Although it is sparse, it might not fit into memory and even if it did, computation times would be huge.
Fortunately, we do not have to solve the whole system at once: we can split Y into sub-images (patches) and do the regression for each of them. Suppose Y has n by n pixels and the oversampling factor is s. Then X will have m = sn rows and columns and the matrix C will have m 2 = s 2 n 2 . Solving a system of p equations takes time proportional to p 3 , so for the full system the time needed would be proportional to s 6 n 6 . If we divide Y into k by k patches, we need a time proportional to s 6 (n/k) 6 for each of them. There are k 2 patches, so total time would be proportional to s 6 n 6 /k 4 . This is an attractive result: dividing Y in 4 by 4 patches already speeds up the calculations 256 times. This analysis shows the principle of divide-andconquer, but it ignores one important aspect. A proper model for an observed patch of Y should include a border area at all sides with half the width of the PSF, to account for contributions by emitters at the borders. We call this width b. This means that a patch has dimensions n = k + 2b by n = k + 2b and thus at a certain value of k, the contribution of b gets so influential that there will be no gain anymore. We also observed that in any iteration to find the L0 penalized solution, elements of ̃ that are close to zero stay close to zero. So they can be eliminated from the equations beforehand, giving an enormous speed-up. More elements are made zero as a result of the iteration, so we can eliminate more and more of them in the subsequent ones.

Practice
The SPIDER algorithm is, as already said, sensitive to high oversampling of the image patches as the matrices used for calculation can be too large to fit the memory of the computer. Therefore, we have performed a test on the exact times SPIDER needs to perform a calculation (per frame) for different patch sizes (k = 4, 8, 16 and 32 pixels) at different oversampling factors (s = 1, 3, 5 and 7) with 50 repetitions per calculation. Calculations were also performed for CSSTORM. The results are shown in Supplementary  Table S2. The data set used for the calculations are the 50 random simulation maps with a density of 10 µm -2 . One should keep in mind that the density of the data set also plays a role in computation time (more emitters on the map mean more signal and thus a more complicated fit). The overlap (i.e. border area b) between the different patches, to accommodate for emitters on the border of the patches, was 4 pixels on every side of the patch.
The calculations were performed on a computer with an Intel(R) Core(TM) i7-4770 CPU @ 3.40GHz CPU and 16 GB RAM memory. The computer works with a 64-bit Windows 7 professional operating system. The version of Matlab used for these calculations is the R2014a version. Supplementary Table S1. Comparing the results obtained with the L0-norm penalty to the ones obtained with the L1-norm penalty in super-resolution imaging. The results obtained on a simulated map of 20 randomly distributed emitters, which corresponds to a molecule density of 1.25 µm -2 , are reported. All the solutions obtained represent the sparsest solutions with a recall rate of 100 %. The analysis was done with an oversampling factor s increasing from 1 to 8, therefore decreasing the size of the pixels on the superresolution grid from 166 nm to 21 nm. The simulated emitters are located at random positions on the images. One can notice that the L0-norm penalty always provides sparse solutions without false positive emitters (for a lateral tolerance disk with diameter 200 nm). This makes SPIDER an excellent method for quantitative analysis of real data.

SPIDER (in seconds)
Patch size (k) Oversampling factor (s)  ) for SPIDER and CSSTORM on the 50 random simulation maps of density 10 µm -2 . The calculations are done for different patch sizes (k = 4, 8, 16 and 32 pixels) and different oversampling factors (s = 1, 3, 5 and 7). The border area b was taken to be 4 pixels on every side of the patch. MEM: insufficient memory for calculations; *: calculation times > 1800 s. It should be clear that SPIDER is faster than CSSTORM with a factor 10 -100.