Expand this Topic clickable element to expand a topic
Skip to content
Optica Publishing Group

Localisation microscopy with quantum dots using non-negative matrix factorisation

Open Access Open Access

Abstract

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

1. Introduction

Localisation microscopy (LM) is a conceptually simple and accessible technique for super-resolution 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 (PSFs). 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.

2. Method description

2.1. 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 wk(x) with corresponding non-negative intensity time traces hk(t):

d(x,t)k=1Kwk(x)hk(t),
where x and t correspond to the spatial (pixel index) and time (frame index) coordinates, respectively. The matrix form of Eq. (1)
D×WH,
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 and the individual sources wk(x) and intensities hk(t) are reshaped into columns and rows of matrices W (N × K) and 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).

 figure: Fig. 1:

Fig. 1: (a) Illustration of NMF. Each frame in the dataset d(x, t) is described as the sum of K individual PSFs wk(x) with appropriate intensity hk(t). (b) iNMF algorithm applied to a single patch.

Download Full Size | PDF

We used multiplicative updates [12] for W and H to solve Eq. (2) (the symbol “∅” denotes the element-wise division of matrices):

wxk=xxkt=1Thkt[(DWH)H]xkhkt=hktx=1Nwxk[W(DWH)]kt.
minimising the (generalised) Kullback-Leibler divergence between data matrix D and its factorised approximation WH. It can be shown [13] that these updates seek to maximise the log-likelihood of the model WH for data D corrupted with Poisson noise
logp(D|W,H)=xt(dxtlogk=1Kwxkhktk=1Kwxkhkt)+const.,
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:

Wj+1(Wj,Hj,D)Hj+1(Wj+1,Hj,D),
where j denotes the iteration of the update.

There is a scaling indeterminacy between W and H in the NMF model. We fix this by forcing the L1 norm of each column of W sum to 1: ∑j wjk = 1. The background fluorescence in the images is modelled as one flat component 1/N with corresponding intensity hK. The spatial part of the background component wK is not updated during the optimisation, while the temporal part hK is updated to account for changes in background levels during the data acquisition due to bleaching or fluctuation of the excitation light, for example.

2.2. 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 and 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 wk of the matrix W according to the L2 norm ( Lq=(jwjkq)1/q). The sources with large L2 norm tend to have “sparser” structure according to Hoyer’s definition of sparseness [15]. Note that wk are normalised to have the L1 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 L2-sorted sources {wk}k=1j are used as initial values for the first j columns of W. The corresponding rows {hk}k=1j are reused for the first j rows of initial values of H. The remaining Kj components are re-initialised from the uniform random distribution. The initial values of W for the (j + 1)th iteration are therefore composed of the j sparsest components of the previous iteration and (Kj) 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 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].

3. Localisation and visualisation

The recovered individual sources wk 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 wp we achieve “shrinking” of the individual in-focus PSFs while keeping some characteristics of each estimated source’s shape (elongation, multiple sources in one wk). If we normalise the L1 norm of wp to one (∑x wp(x) = 1), we can reconstruct a “super-resolution” image by summing all wkp, weighted by the corresponding mean intensity mean(hk). Note that up-sampling of wks is needed before taking the powers p. We used bi-cubic interpolation for the four times up-sampled wks 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.

4. Simulated and experimental data

4.1. 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 · 103 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 nm 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.

 figure: Fig. 2:

Fig. 2: (a–b) Classification of true positives (TP), false positives (FP) and false negatives (FN) from the true (red dot) and estimated (green cross) locations. (c) Example of the precision P(li) (blue) and recall R(li) (green) curve against the confidence levels li for one iNMF evaluation. (d) The precision/recall curve P(R) (blue) with interpolated precision Pinterp() (red) reducing the wiggles. The average precision (green dashed line) was estimated from averaging of the red curve.

Download Full Size | PDF

 figure: Fig. 3:

Fig. 3: Simulated data of randomly scattered emitters. (a) Average precision which summarises both localisation precision and the ability to recover the individual sources (higher values are better). Error bars show standard deviation for repeated simulations. (b–c) Red dots in the wide-field (WF) image show the true positions of the emitters. Green crosses are the identified locations for iNMF (visualisation of one run of iNMF), CSSTORM, bSOFI and 3B reconstructions. Red circles indicate the threshold distance r from the true locations. A green cross within a red circle was counted as a true positive (TP). A green cross outside a red circle was counted as false positive (FP). Red circle without any green cross inside was counted as a false negative (FN). Red arrows point to close unresolved sources. Green arrow point at sources resolved only in the iNMF reconstruction. Blue arrows show where multiple emitters were approximated by a single source in the iNMF image.

Download Full Size | PDF

 figure: Fig. 4:

Fig. 4: Comparison of wide field (WF), iNMF, CSSTORM, bSOFI and 3B techniques. (a) Reconstruction of a simulated structure. The true sources’ positions are indicated as red dots in the WF image. The red arrows point at sub-resolution double line structure recovered with iNMF, CSSTORM and bSOFI. The green arrows point to a sub-resolution hole visible only in the iNMF image. The blue arrows show the area where 3B and CSSTORM failed to resolve the double line of the artificial specimen. (b) Tubulin structure labelled with QDs. Scale bars 400nm.

Download Full Size | PDF

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

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

 figure: Fig. 5:

Fig. 5: iNMF reconstructed images. (a) Tubulin labelled with QDs. Red dashed box shows the region used for comparison of the methods in Fig. 4(b). (b) Tubulin labelled with Alexa 647. Scale bars 500nm.

Download Full Size | PDF

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 103 frames with 100 ms/frame acquisition.

 figure: Fig. 6:

Fig. 6: Recovery of the different PSFs. (a) Movie of blinking QDs randomly scattered on an out-of-focus plane. 103 frames in total. (b) Individual PSFs wk recovered with iNMF (K = 22). The bars under each frame show the mean intensity of the source. (c) wk recovered with ICA (K = 22). Blue colour indicate the pixels with negative values. Scale bars 1μm.

Download Full Size | PDF

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.

 figure: Fig. 7:

Fig. 7: (a) A neurone with neurotransmitter receptor subunits labeled with QD605. Inset in the upper right corner shows the PSF obtained from measurement of the axial profile of a fluorescence bead (WF), iNMF reconstructed region of low QD density taken at different axial positions (iNMF) and simulated PSF (Sim) used for localisation. The axial localisation is colour-coded (colorbar with corresponding PSF slices) with transparency proportional to the mean brightness of each estimated source. Data were kindly provided by Anja Huss. (b) Simulated data with PRILM [26] tailored PSF. Experimental PSF (shown on the left) has been used for data simulations. Data kindly provided by David Baddeley. The PSFs are shown on the same scale as the large images.

Download Full Size | PDF

5. 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 wk 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 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 σ=22πλemNA 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

bk=meant(hkt),
retrieved from the matrix H as a mean along the rows. The interval [lmin, lmax] between the dimmest lmin = mink(bk) and the brightest lmax = maxk(bk) source intensity was divided into a number of intervals (confidence levels) li defined by the steps in the sorted intensities of all sources. For each confidence level li only the sources with bk above li were considered. True positives (TPi), false negatives (FNi) and false positives (FPi) were computed for each confidence level li.

The precision P and recall R were computed from TP(li), FP(li) and FN(li) for each confidence level li:

P(li)=TP(li)TP(li)+FP(li)
R(li)=TP(li)TP(li)+FN(li).
An example of precision P(li) and recall R(li) 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 levels i ∈ [0 : 0.1 : 1] by taking the maximum precision for which the corresponding recall exceeds i (Fig. 2(d)):
Pinterp(R˜)=maxR;RR˜P(R).
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:

AP=111R˜Pinterp(R˜).

6. Results

6.1. 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) we show the average precision (AP) Eq. (10) as a quantitative measure of the algorithm performance on simulated data of randomly scattered sources. AP summarises both localisation precision and the ability to recover the individual emitters; higher values of AP are better.

For simulated data (103 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 Ktrue. 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 Ktrue 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)).

6.2. 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 and 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 (103 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 imunolabelled with QDs (QD625, Invitrogen) with emission peak at λ = 625 nm. We analysed 103 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 (103 frames with 21 × 21 pixels2) was:

  • iNMF: ∼ 20 mins for K = 50 sources and K restarts.
  • CSSTORM: ∼ 260 mins.
  • bSOFI: ∼ 1 min.
  • 3B: ∼ 12 hours for 30 iterations.
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.

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

Acknowledgments

We thank Stefan Geissbuehler and Marcel Leutenegger for providing the bSOFI algorithm and help with the bSOFI evaluation and David Baddeley and Anja Huss for providing us with data of three dimensional samples. We thank to Susan Cox, Martin Kielhorn and Kai Wicker for interesting discussions. This work was supported in part by the grants EP/F500385/1 and BB/F529254/1 to the University of Edinburgh School of Informatics Doctoral Training Centre in Neuroinformatics and Computational Neuroscience (www.anc.ac.uk/dtc) from the UK Engineering and Physical Sciences Research Council (EPSRC), UK Biotechnology and Biological Sciences Research Council (BBSRC), and the UK Medical Research Council (MRC).

References and links

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

2. A. R. Small, “Theoretical limits on errors and acquisition rates in localizing switchable fluorophores,” Biophys. J. 96(2), 16–18 (2009). [CrossRef]  

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

4. F. D. Stefani, J. P. Hoogenboom, and E. Barkai, “Beyond quantum jumps: blinking nanoscale light emitters,” Phys. Today 62(2), 34 (2009). [CrossRef]  

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

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

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

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

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

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

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

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). The MathWorks IncNatick, Massachusetts.

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

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

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), 69004 (2013). [CrossRef]  

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

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

25. A. Hyvärinen and E. Oja, “Independent component analysis: algorithms and applications,” Neural Netw. 13(4–5), 411–430 (2000). [CrossRef]   [PubMed]  

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

Cited By

Optica participates in Crossref's Cited-By Linking service. Citing articles from Optica Publishing Group journals and other participating publishers are listed here.

Alert me when this article is cited.


Figures (7)

Fig. 1:
Fig. 1: (a) Illustration of NMF. Each frame in the dataset d(x, t) is described as the sum of K individual PSFs wk(x) with appropriate intensity hk(t). (b) iNMF algorithm applied to a single patch.
Fig. 2:
Fig. 2: (a–b) Classification of true positives (TP), false positives (FP) and false negatives (FN) from the true (red dot) and estimated (green cross) locations. (c) Example of the precision P(li) (blue) and recall R(li) (green) curve against the confidence levels li for one iNMF evaluation. (d) The precision/recall curve P(R) (blue) with interpolated precision Pinterp() (red) reducing the wiggles. The average precision (green dashed line) was estimated from averaging of the red curve.
Fig. 3:
Fig. 3: Simulated data of randomly scattered emitters. (a) Average precision which summarises both localisation precision and the ability to recover the individual sources (higher values are better). Error bars show standard deviation for repeated simulations. (b–c) Red dots in the wide-field (WF) image show the true positions of the emitters. Green crosses are the identified locations for iNMF (visualisation of one run of iNMF), CSSTORM, bSOFI and 3B reconstructions. Red circles indicate the threshold distance r from the true locations. A green cross within a red circle was counted as a true positive (TP). A green cross outside a red circle was counted as false positive (FP). Red circle without any green cross inside was counted as a false negative (FN). Red arrows point to close unresolved sources. Green arrow point at sources resolved only in the iNMF reconstruction. Blue arrows show where multiple emitters were approximated by a single source in the iNMF image.
Fig. 4:
Fig. 4: Comparison of wide field (WF), iNMF, CSSTORM, bSOFI and 3B techniques. (a) Reconstruction of a simulated structure. The true sources’ positions are indicated as red dots in the WF image. The red arrows point at sub-resolution double line structure recovered with iNMF, CSSTORM and bSOFI. The green arrows point to a sub-resolution hole visible only in the iNMF image. The blue arrows show the area where 3B and CSSTORM failed to resolve the double line of the artificial specimen. (b) Tubulin structure labelled with QDs. Scale bars 400nm.
Fig. 5:
Fig. 5: iNMF reconstructed images. (a) Tubulin labelled with QDs. Red dashed box shows the region used for comparison of the methods in Fig. 4(b). (b) Tubulin labelled with Alexa 647. Scale bars 500nm.
Fig. 6:
Fig. 6: Recovery of the different PSFs. (a) Movie of blinking QDs randomly scattered on an out-of-focus plane. 103 frames in total. (b) Individual PSFs wk recovered with iNMF (K = 22). The bars under each frame show the mean intensity of the source. (c) wk recovered with ICA (K = 22). Blue colour indicate the pixels with negative values. Scale bars 1μm.
Fig. 7:
Fig. 7: (a) A neurone with neurotransmitter receptor subunits labeled with QD605. Inset in the upper right corner shows the PSF obtained from measurement of the axial profile of a fluorescence bead (WF), iNMF reconstructed region of low QD density taken at different axial positions (iNMF) and simulated PSF (Sim) used for localisation. The axial localisation is colour-coded (colorbar with corresponding PSF slices) with transparency proportional to the mean brightness of each estimated source. Data were kindly provided by Anja Huss. (b) Simulated data with PRILM [26] tailored PSF. Experimental PSF (shown on the left) has been used for data simulations. Data kindly provided by David Baddeley. The PSFs are shown on the same scale as the large images.

Equations (10)

Equations on this page are rendered with MathJax. Learn more.

d ( x , t ) k = 1 K w k ( x ) h k ( t ) ,
D × WH ,
w x k = x x k t = 1 T h k t [ ( D WH ) H ] x k h k t = h k t x = 1 N w x k [ W ( D WH ) ] k t .
log p ( D | W , H ) = x t ( d x t log k = 1 K w x k h k t k = 1 K w x k h k t ) + const . ,
W j + 1 ( W j , H j , D ) H j + 1 ( W j + 1 , H j , D ) ,
b k = mean t ( h k t ) ,
P ( l i ) = TP ( l i ) TP ( l i ) + FP ( l i )
R ( l i ) = TP ( l i ) TP ( l i ) + FN ( l i ) .
P interp ( R ˜ ) = max R ; R R ˜ P ( R ) .
AP = 1 11 R ˜ P interp ( R ˜ ) .
Select as filters


Select Topics Cancel
© Copyright 2024 | Optica Publishing Group. All rights reserved, including rights for text and data mining and training of artificial technologies or similar technologies.