Localisation microscopy with quantum dots using non-negative matrix factorisation

We propose non-negative matrix factorisation with iterative restarts (iNMF) to model a noisy dataset of highly overlapping fluorophores with intermittent intensities. We can recover high-resolution images of individual sources from the optimised model, despite their high mutual overlap in the original data. Each source can have an arbitrary, unknown shape of the PSF and blinking behaviour. This allows us to use quantum dots as bright and stable fluorophores for localisation microscopy. We compare the iNMF results to CSSTORM, 3B and bSOFI. iNMF shows superior performance in the challenging task of super-resolution imaging using quantum dots. We can also retrieve axial localisation of the sources from the shape of the recovered PSF. © 2014 Optical Society of America OCIS codes: (180.2520) Fluorescence microscopy; (100.6640) Superresolution. References and links 1. A. Small and S. Stahlheber, “Fluorophore localization algorithms for super-resolution microscopy,” Nat. Methods 11(3), 267–279 (2014). 2. A. R. Small, “Theoretical limits on errors and acquisition rates in localizing switchable fluorophores,” Biophys. J. 96(2), 16–8 (2009). 3. U. Resch-Genger, M. Grabolle, S. Cavaliere-Jaricot, R. Nitschke, and T. Nann, “Quantum dots versus organic dyes as fluorescent labels,” Nat. Methods 5(9), 763–775 (2008). 4. F. D. Stefani, J. P. Hoogenboom, and E. Barkai, “Beyond quantum jumps: blinking nanoscale light emitters,” Phys. Today 62(2), 34 (2009). 5. F. Huang, S. L. Schwartz, J. M. Byars, K. A. Lidke, and F. Huang, “Simultaneous multiple-emitter fitting for single molecule super-resolution imaging,” Biomed. Opt. Express 2(5), 1377–1393 (2011). 6. S. J. Holden, S. Uphoff, and A. N. Kapanidis, “DAOSTORM: an algorithm for high-density super-resolution microscopy,” Nat. Methods 8(4), 279–280, (2011). 7. L. Zhu, W. Zhang, D. Elnatan, and B. Huang, “Faster STORM using compressed sensing,” Nat. Methods 9(7), 721–723 (2012). 8. K. A. Lidke, B. Rieger, T. M. Jovin, and R. Heintzmann, “Superresolution by localization of quantum dots using blinking statistics,” Opt. Express 13(18), 7052 (2005). 9. A. Barsic and R. Piestun, “Super-resolution of dense nanoscale emitters beyond the diffraction limit using spatial and temporal information,” Appl. Phys. Lett. 102(23), 231103 (2013). 10. 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, (2011). #217119 $15.00 USD Received 16 Jul 2014; revised 11 Sep 2014; accepted 11 Sep 2014; published 30 Sep 2014 (C) 2014 OSA 6 October 2014 | Vol. 22, No. 20 | DOI:10.1364/OE.22.024594 | OPTICS EXPRESS 24594 11. S. Geissbuehler, N. L. Bocchio, C. Dellagiacoma, C. Berclaz, M. Leutenegger, and T. Lasser, “Mapping molecular statistics with balanced super-resolution optical fluctuation imaging (bSOFI),” Opt. Nanoscopy 1(1), 4 (2012). 12. D. D. Lee and H. S. Seung, “Algorithms for non-negative matrix factorization,” Adv. Neural Inf. Process. Syst. 13, 556–562 (2001). 13. O. Mandula, “Super-resolution methods for fluorescence microscopy,” PhD thesis, University of Edinburgh (2012). 14. J. Kim and H. Park, “Toward Faster Nonnegative Matrix Factorization: A New Algorithm and Comparisons,” 2008 Eighth IEEE International Conference on Data Mining (2008), pp. 353–362. 15. P. O. Hoyer, “Non-negative matrix factorization with sparseness constraints,” J. Mach. Learn. Res. 5, 1457–1469 (2004). 16. MATLAB. (2010). version 7.10.0 (R2010a). Natick, Massachusetts: The MathWorks Inc. 17. https://www.dropbox.com/s/731u3u4htpdndgr/inmf.zip 18. https://github.com/aludnam/inmf 19. A. G. York, A. Ghitani, A. Vaziri, M. W. Davidson, and H. Shroff, “Confined activation and subdiffractive localization enables whole-cell PALM with genetically expressed probes,” Nat. Methods 8(4), 327–333 (2011). 20. T. Dertinger, R. Colyera, G. Iyera, S. Weissa, and J. Enderleind, “Fast, background-free, 3D super-resolution optical fluctuation imaging (SOFI). Supporting Information,” PNAS 106(52), 22287–22292 (2009). 21. http://bigwww.epfl.ch/smlm/datasets/index.html?p=real-hd 22. N. Olivier, D. Keller, P. Gonczy, and S. Manley, “Resolution doubling in 3D-STORM imaging through improved buffers,” PLoS One 8(7), (2013). 23. M. Everingham, L. Gool, C. K. I. Williams, J. Winn, and A. Zisserman, “The Pascal Visual Object Classes (VOC) Challenge,” Int. J. Comput. Vis. 88(2), 303–338 (2009). 24. B. Zhang, J. Zerubia, and J.-C. Olivo-Marin, “Gaussian approximations of fluorescence microscope point-spread function models” Appl. Opt. 46(10), 1819–1829 (2007). 25. A. Hyvärinen and E. Oja, “Independent component analysis: algorithms and applications,” Neural Netw. 13(4-5), 411–430 (2000). 26. D. Baddeley, M. B. Cannell, and C. Soeller, “Three-dimensional sub-100 nm super-resolution imaging of biological samples using a phase ramp in the objective pupil,” Nano Res. 4(6), 589–598 (2011).


Introduction
Localisation microscopy (LM) is a conceptually simple and accessible technique for superresolution imaging of fluorescent samples. LM analyses a time-lapse sequence of images of a specimen labelled with fluorophores stochastically switching between bright (ON) and dark (OFF) states, and estimates the positions of the individual sources with sub-pixel precision, resulting in a super-resolution image of the specimen.
Conventional LM methods (PALM, fPALM, STORM, dSTORM) [1] actively drive a vast majority of the fluorophores into the OFF state ensuring a low density of ON sources in each frame. This effectively avoids the overlaps between the individual point spread functions (PS-Fs). The localisation of the well-separated PSFs is simple but only few sources are localised in each frame. Due to this extremely low throughput several tens of thousands frames have to be acquired to reconstruct an image.
Higher density of the ON sources in each frame can increase the throughput of the data and hence decrease the acquisition time [2], but one has to use more sophisticated algorithms to localise the overlapping PSFs. Accounting for overlapping PSFs allows us to consider quantum dots (QDs) as fluorescent markers for LM despite that they stay in the ON state for a significant fraction of time. In addition to an order of magnitude greater photostablity, QDs have an order of magnitude higher brightness compared to organic dyes or fluorescent proteins [3] which can further decrease the required acquisition time per frame. All these properties make QDs very attractive for fluorescence microscopy. Under continuous excitation QDs exhibit stochastic blinking [4] which simplifies the acquisition procedure as no active ON-OFF switching (fPALM, STORM), nor presence of specialised buffers (dSTORM) is required. However, the blinking cannot be controlled and the QD-labelled data typically contain highly overlapping PSFs. Super-resolution imaging with QDs remains a challenging task.
Several methods dealing with overlapping sources have been proposed recently. Among methods analysing each frame of the LM dataset individually [5,6] CSSTORM [7] has shown a superior performance. Another group of algorithms models the whole LM dataset as a collection of blinking emitters. These techniques take the reappearance of the fluorophores in different frames into account, and allow the localisation of higher densities of ON sources. Among other techniques [8,9] a Bayesian treatment of the sources' position and intensity estimation (3B [10]) has received attention recently. Finally, bSOFI [11] recovers the distribution of fluorophores with high resolution from the analysis of the intensity fluctuation higher order statistics.

Non-negative matrix factorisation
We propose non-negative matrix factorisation (NMF) as a model for data with overlapping blinking sources. NMF finds an approximate factorisation of the dataset d(x,t) into a set of K individual non-negative PSFs w k (x) with corresponding non-negative intensity time traces h k (t): where x and t correspond to the spatial (pixel index) and time (frame index) coordinates, respectively. The matrix form of Eq.(1) restricted to matrices with non-negative entries, shows the matrix factorisation explicitly. In Eq.
(2) T frames of the dataset d(x,t) were reshaped into a N × T matrix D D D and the individual sources w k (x) and intensities h k (t) are reshaped into columns and rows of matrices W W W (N × K) and H H H (K × T ), respectively. N denotes the total number of pixels in each frame. The non-negativity is a strong natural constraint for both the PSFs and the intensity time profiles. An important assumption is the reappearance of the sources in different time frames throughout the dataset. This creates the redundancy, which together with the non-negativity constraints allows us to recover both the shape and the intensity sequence for each individual source. These can be different and unique and are inferred directly from the data. Note the difference of NMF to blind deconvolution, where the goal is to estimate one unknown PSF shared by all sources. NMF is illustrated in Fig.1(a).
We used multiplicative updates [12] for W W W and H H H to solve Eq. (2) (the symbol " " denotes the element-wise division of matrices): Replace first j columns of Winit with first j columns of sorted W . (Correspondingly for rows of Hinit.) which makes them particularly suitable for data with low photon counts due to short acquisition time.
The NMF model is fitted to data iteratively using the updates Eq. (3) sequentially: where j denotes the iteration of the update. There is a scaling indeterminacy between W W W and H H H in the NMF model. We fix this by forcing the L 1 norm of each column of W W W sum to 1: ∑ j w jk = 1. The background fluorescence in the images is modelled as one flat component 1/N with corresponding intensity h h h K . The spatial part of the background component w w w K is not updated during the optimisation, while the temporal part h h h K is updated to account for changes in background levels during the data acquisition due to bleaching or fluctuation of the excitation light, for example.

NMF with iterative restarts (iNMF)
NMF becomes challenging when applied to large datasets. Besides computational time and memory requirements, local minima further complicate the NMF optimisation [14]. We address this partly by dividing the data into partially overlapping patches (typically 25×25 pixels with 5 pixels overlap), so that NMF is applied to each patch individually and then stitching the patches back together, blending the border pixels with weights linearly decreasing towards the border. Note that the individual patches are evaluated independently which permits a simple parallelization of the evaluation. We also address this problem via the iterative NMF algorithm described in the next paragraph.
Although the Lee and Seung updates are convex with respect to W W W and H H H separately, they are non-convex in both simultaneously [12]. Multiple restarts can be used to address the problem of local optima, but we have not found good solutions with this approach. Instead, we exploit some prior knowledge about the problem, namely that the PSFs are likely to have a fairly compact structure. We therefore rank the columns w w w k of the matrix W W W according to the L 2 norm (L q = (∑ j w q jk ) 1/q ). The sources with large L 2 norm tend to have "sparser" structure according to Hoyer's definition of sparseness [15]. Note that w w w k are normalised to have the L 1 norm equal to one. This leads to an iterative NMF (iNMF) algorithm summarised in Fig. 1(b), which we empirically found to perform well on the LM data. The iNMF algorithm runs as follows: on iteration ( j + 1) the first j L 2 -sorted sources {w w w k k k } j k=1 are used as initial values for the first j columns of W W W . The corresponding rows {h h h k k k } j k=1 are reused for the first j rows of initial values of H H H. The remaining K − j components are re-initialised from the uniform random distribution. The initial values of W W W for the ( j + 1)th iteration are therefore composed of the j sparsest components of the previous iteration and (K − j) randomly initialised components. The NMF algorithm is then run with this initialisation, updating all K components. The iterative procedure runs until j = K.
iNMF represents a "soft" sparsity enhancement, which is in contrast with Hoyer's sparse NMF algorithm [15], where the desired sparsity of the components is explicitly enforced. Note that iNMF allows the recovery of the sources each having different sparseness.
NMF requires knowledge of the number of sources K in the dataset. We used a crude estimation of the upper bound on K from the sorted principal values of the data matrix D D D because, as shown on simulated data, iNMF progressively recovers the optimal number of PSFs if K is overestimated [13]. The redundant components are used for approximation of background noise. Our iNMF algorithm written in MATLAB [16] can be downloaded from [17] or can be cloned from github [18].

Localisation and visualisation
The recovered individual sources w w w k k k can be localised in 3D with, for example, a technique based on cross-correlation with measured or simulated PSFs [19]. All the estimated locations can be plotted to create a super-resolution image.
For the samples with mostly in-focus fluorophores we propose a simple alternative visualisation, which does not require fitting of the individual emitters. By taking the pixel-wise power p > 1 of the estimated sources w w w p we achieve "shrinking" of the individual in-focus PSFs while keeping some characteristics of each estimated source's shape (elongation, multiple sources in one w w w k ). If we normalise the L 1 norm of w w w p to one (∑ x w w w p (x) = 1), we can reconstruct a "superresolution" image by summing all w w w p k , weighted by the corresponding mean intensity mean(h h h k k k ). Note that up-sampling of w w w k s is needed before taking the powers p. We used bi-cubic interpolation for the four times up-sampled w w w k s images. A false colour-map can be used to enhance the dim features of the reconstructed image.
This visualisation method is only one, rather convenient, way of displaying the iNMF results. However, it is not applicable to 3D samples, where the PSFs of each source can significantly differ. These results have to be visualised in a conventional way by plotting the estimated (3D) locations. Note that for some applications the main interest is in the explicit locations of the individual sources rather then in an image of the structure.

Simulations
An emission wavelength of λ = 625 nm and a 1.3 NA oil objective was considered for the simulation of the PSF. We used the maximum emitted intensity of 5 · 10 3 photons/source/frame (perhaps slightly strong estimate for QD data), a uniform background 100 photons/pixel and an ON/OFF time ratio of 1.0 (on average, 50% of sources were ON in each frame which results in high overlap of the PSFs). We generated a movie of blinking fluorophores from the simulated PSFs and the intensity profiles. We add a uniform background and each pixel has been corrupted with Poisson noise.
The distance between the parallel lines in the artificial structure in Fig. 4(a) was set to approximately half the Airy disk's radius (150 nm for an emission wavelength of λ = 625 n-m and 1.3 NA). The linear density of the QDs attached to the collinear structure was set to µ = 15 µm −1 corresponding to approximately 67 nm spacing between the adjacent sources.
The blinking of the QDs was simulated as a down-sampled telegraph process. The telegraph process is a Markov process, which takes only two values, in our case ON and OFF fluorescence states. To avoid an unrealistic binary blinking, we simulated the intensity on a q = 100 times up-sampled time axis as a telegraph process with a rate ofγ = 0.5/q. This up-sampled signal was binned over q bins. Note that our simulated ON and OFF time probability density follows an exponential decay in contrast to the QDs where the ON and OFF time are known to follow an inverse power law [4].

Experimental data
HEp-2 cells (ATCC CCL-23) (Fig. 4(b), and 5(a)) were grown on high precision coverslips No 1.5 (Zeiss) in DMEM medium. The cells has been fixed and embedded in 1% w/v Mowiol (Sigma Aldrich) following the protocol from [20]. Imaging was with a 63× 1.4 NA oil objective lens (image pixel size 80 nm) with laser excitation light at 405 nm.
Randomly scattered QDs (Fig. 6) were obtained by depositing a droplet of diluted QD565 Invitrogen on a cover slip. We used a 100× 1.4 NA oil objective lens (83 nm pixel-size) to acquired 10 3 frames with 100 ms/frame acquisition.
Alexa 647 labeled data ( Fig. 5(b)) have been downloaded from [21] and used with kind permission of Suliana Manley. Data (500 frames) have been taken using an unoptimized buffer containing COT [22] with 1.3NA oil objective with 40 ms/frame acquisition time.
A neurone with neurotransmitter receptor subunits labeled with QD605 ( Fig. 7(a)), was imaged with 100×, 1.4NA oil objective. Data were kindly provided by Anja Huss.

Quantitative comparison with other methods
In this section we compare iNMF with CSSTORM, bSOFI and 3B on their ability to detect and correctly localise sources. The average precision (AP) measure [23] allows a quantitative performance comparison for different algorithms, summarising both localisation precision and the ability to recover the individual sources. The individual estimated sources w w w k were localised by ML fitting of a Gaussian approximation of the PSF. We used a greedy algorithm to assign the estimated locations (E) to their nearest true positions (T ). Only one estimated position c" d" was assigned to each true position. If the distance between the estimated and the true position was smaller than a threshold r, then the source was considered as a true positive (TP). Each true position with no estimated source within a disk of radius r was counted as false negative (FN), whereas an estimated position further than r from any true position was considered as false positive (FP), see Fig. 2(a). M estimated sources in the proximity of one true source were counted as 1TP and (M − 1)FP. One estimated source in proximity of M true sources gives 1TP and (M − 1)FN, see Fig. 2(b). We set the threshold r = σ /2, where σ = √ 2 2π λ em NA is the standard deviation of the in-focus PSF Gaussian approximation [24]. For the parameters used in our simulations the threshold corresponds to r = 56 nm (0.7 pixels). The estimated positions E were ranked according to the source's mean brightness retrieved from the matrix H H H as a mean along the rows. The interval [l min , l max ] between the dimmest l min = min k (b k ) and the brightest l max = max k (b k ) source intensity was divided into a number of intervals (confidence levels) l i defined by the steps in the sorted intensities of all sources. For each confidence level l i only the sources with b k above l i were considered. True positives ( TP i ), false negatives ( FN i ) and false positives ( FP i ) were computed for each confidence level l i . The precision P and recall R were computed from TP(l i ), FP(l i ) and FN(l i ) for each confidence level l i : An example of precision P(l i ) and recall R(l i ) curves for different confidence levels is shown in Fig. 2(c). Following Everingham et al. [23], the precision/recall (PR) curve P(R) was interpolated for 11 equally spaced recall levelsR i ∈ [0 : 0.1 : 1] by taking the maximum precision for which the corresponding recall exceedsR i (Fig. 2(d)): The precision/recall (PR) curve was interpolated in order to reduce the impact of "wiggles" in the PR curve ( Fig. 2(d)). Note that to obtain a high AP, the method must have precision at all levels of recall. This penalises the methods that can accurately estimate only few very bright sources. Average precision (AP) is then defined as the mean of interpolated precision: 6. Results

Evaluations of the data with other methods
We compared iNMF with CSSTORM [7], 3B [10] and bSOFI [11] on both simulated and experimental data. The same data sets were used for all methods. In Fig. 3(a)   the ability to recover the individual emitters; higher values of AP are better. For simulated data (10 3 frames) the true background value of 100 photons/pixel/frame was subtracted (clipping any negative values to zero) before CSSTORM and 3B evaluation. The true PSF was provided to the CSSTORM, bSOFI and 3B algorithms. The 3B algorithm was also provided with the true number of sources. The estimated background value and PSF corresponding to the parameters of the imaging setup were used for evaluation of the experimental data. If not specified otherwise, we left the parameters set to the published code default values.
CSSTORM processes each input frame individually, independent of the rest of the dataset. This method tries to recover a sparse distribution of the active (ON) fluorophores in each frame considering a known PSF (shared by all sources). The output for each frame is an image showing the possible positions of the sources on a sub-pixel grid (8 times oversampled). For AP estimation we processed the sum of all CSSTORM output frames, which summarises all the estimated sources. The summed image was filtered with a Gaussian kernel (σ = 0.8 pixel, which corresponds to σ = 8 nm). In the evaluation of the simulated data of randomly scattered sources we identified the local maxima stronger than 5% of the global maximum (green crosses in Figs. 3(b)-3(c)). We chose the threshold of 5% because for this value the number of local maxima roughly corresponds to the true number of sources K true . The positions of the local maxima were used as the estimated sources' positions for computation of the AP.
For the 3B analysis the prior parameters for the size of the PSF were adjusted to the true values. Also the true number of sources K true was used as an initial number of spots in the model. The 3B algorithm was run for at least 30 iterations (we typically stopped the algorithm after 12 hours of evaluation). Following [10], the output coordinates of the 3B analysis were placed on a 100 times oversampled grid (0.8 nm pixel-size) and convolved with a Gaussian (σ = 10 pixels, which corresponds to σ = 8 nm). Similar to the analysis with the CSSTORM, we identified local maxima in the reconstructed images. Only the maxima above a threshold were considered for evaluation (green crosses in Figs. 3(b)-3(c)). The threshold was set individually for each image, such that the number of local maxima roughly corresponds to the number of sources considered by the 3B analysis after the last iteration. The bSOFI evaluation was performed by Marcel Leutenegger and Stefan Geissbuehler using the SOFI software package. Following the bSOFI authors standard procedure, the data sequence for bSOFI evaluation in Fig. 4 were divided into two shorter subsequences each having half number of frames. We show the average of the evaluation of the two subsequences of the fourth order bSOFI image. For comparison shown in Fig. 3 the positions of the local maxima (stronger then 5% of the global maximum) were used as the estimates of the sources locations (green crosses in Figs. 3(b)-3(c)).

Comparison of iNMF with other methods
iNMF shows superior performance especially at higher density levels (blue marks in Fig. 3(a)). However, very densely packed emitters can sometimes be approximated by a lower number of individual PSFs (blue arrows in Fig. 3(c)), which can create discontinuities in the visualised structure. Different runs of iNMF can recover different subsets of fluorophores due to random initialisation of matrices W W W and H H H. We can interpret this as different random samples from the distribution of the fluorophores, and below we show an average of several iNMF reconstructions which results in smoother representation of the specimen structure. For AP comparison in Fig. 3(a) we always used only one iNMF evaluation to allow a fair comparison. Fig. 4(a) shows a comparison of the four methods on simulated data of an artificial structure (10 3 frames). 3B completely fails to recover the double parallel lines (red arrows in Fig. 4(a)), replacing them with one intensity crest in between the lines (blue arrow). CSSTORM shows the double line structure at the periphery of the specimen (red arrows) but the double lines join into a single line close to the centre of the cross (blue arrow). The hole in the middle of the specimen (green arrow) is completely unresolved and is replaced by an intensity maximum in the CSSTORM image. The bSOFI image suggests a double-line structure at some parts of the cross (red arrows) but the structure is discontinuous, especially at the lower part of the cross. iNMF reveals the double line structure in the full length (red arrows in Fig. 4(a)) and the hole in the middle of the structure is clearly resolved (green arrow).
We further tested the performance of the methods on a real biological specimen -tubulin fibres of a of HEp-2 cell imuno-labelled with QDs (QD625, Invitrogen) with emission peak  at λ = 625 nm. We analysed 10 3 frames with 50 ms/frame. The total acquisition took approximately 1 min. A comparison of the methods on a small patch is shown in Fig. 4(b). The wide-field image and iNMF suggests the crossing of two fibres. The structure recovered by CSSTORM is blurry near the presumed crossing and the bSOFI structure is discontinuous, similar to the simulation results above. The 3B analysis replaces the close-by fibres by a single intensity crest similar to Fig. 4(a). The iNMF results in Fig. 4 can be reproduced with the provided iNMF software containing corresponding LM datasets. The iNMF reconstruction of a larger area of the tubulin fibres is shown in Fig. 5(a) with highlighted regions showing high-resolution features. For our computer (Intel ® Core™2 Duo @ 2GHz processor with 3GB of RAM) the computational time for the simulated dataset (10 3 frames with 21 × 21 pixels 2 ) was: • iNMF: ∼ 20 mins for K = 50 sources and K restarts.
Note that the iNMF results in Fig. 4 are averages of 10 iNMF evaluation. The computation time is therefore comparable (∼ 200 mins) to CSSTORM.
To demonstrate that iNMF can also be used with standard fluorophores we analysed images of tubulin structures labelled with Alexa 647 (500 frames with 40ms/frame acquisition time). The reconstructed image is shown in Fig. 5(b) with enlarged areas showing resolution improvement in the reconstructed structure.
An interesting feature of iNMF is the possibility to recover sources with different unknown PSFs. In Fig. 6(b) we demonstrate the recovery of several out-of-focus PSFs from blinking QDs randomly dispersed on a cover-slip Fig. 6(a). The individual out-of-focus PSFs have been recovered without any prior knowledge about their shape. Independent component analysis (ICA) proposed by Lidke et al. [8], allows recovery of different individual PSFs. However, ICA's performance is poor when applied to noisy data due to the negative entries and the lack of the noise statistic consideration [13]. ICA evaluation of the out-of-focus QDs is shown in Fig. 6(c). Prior to the ICA evaluation (FastICA algorithm [25]) the background was subtracted (clipping any negative values to zero). The number of sources was set to K = 22. We used "tanh" as the nonlinearity option in the fixed-point algorithm. The blue pixels indicate the negative values in the estimated sources. Most of the estimated components clearly combine several overlapping PSFs together. The overall quality is inferior to the iNMF results which can be reproduced with the iNMF software containing corresponding LM data.
The shape of each PSF can be used for axial localisation of the source, using a technique based on cross-correlation with simulated PSFs [19] (see upper right inset of Fig. 7(a)). Fig. 7(a) shows a reconstruction of a neurone with neurotransmitter subunits labelled with QDs acquired on a standard wide-field microscope. The estimated axial position is colour-coded. Only the sources where the correlation index with the simulated PSF exceeds 0.7 were considered for visualisation. The area in the center of the cell contains too high density of out-of-focus sources which prevents iNMF from recovering the individual sources.
The PSF varies very slowly close to the focal point making the axial localisation challenging. An interesting alternative is to use a specially designed PSF such as PRILM [26] where the PSF changes its shape dramatically with slight defocus. Fig. 7(b) shows reconstruction of simulated data using a PRILM tailored PSF.
Note that analysis of data with overlapping different individual PSFs is beyond ability of 3B, CSSTORM and bSOFI. The 3B analysis adjusts the width of each PSF but the same PSF shape is shared across all fluorophores. CSSTORM and bSOFI use one fixed PSF for all sources.

Conclusion
iNMF enlarges the family of localisation microscopy techniques, and enables using quantum dots as fluorescent labels for localisation microscopy. iNMF shows superior performance for simulated data of blinking QDs with highly overlapping PSFs when compared to CSSTORM, bSOFI and 3B. The access to the shape of each individual PSF allows us to localise each source  in 3D, making iNMF a promising technique for super-resolution images of three-dimensional samples.