1. Introduction

Ever since its birth, localization-based super-resolution fluorescence microscopy (or called localization microscopy) has fascinated a school of scientists due to its attractive capability of providing ultra-high spatial resolution (typically 25-40 nm in xy and ~60 nm in z direction) from relatively simple setup [1, 2]. This technique employs clever strategy to separate close-packed fluorescence emitters into hundreds or even thousands of image frames, so that the locations of individual emitters can be obtained precisely through further molecule localization. These locations are then used to reconstruct a final image whose resolution is not limited by far-field diffraction but by photon-counting statistics [3]. Unfortunately, to obtain such high spatial resolution, localization microscopy greatly scarifies its temporal resolution. This severely limits the versatility and power of the technique in studying dynamic processes inside living samples. Therefore, in recent years there is a strong need to push forward the temporal resolution of localization microscopy while maintaining a desired spatial resolution (thereafter called fast localization microscopy) [4, 5].

To find ways to meet such a demand, it is useful to look into the factors which determine the spatial resolution of localization microscopy. According to Shroff et al [6], these factors include localization precision, which depends on the number of photons detected from individual emitters and is 20~30 nm for typically fluorescence emitters [7], and the Nyquist resolution, which relies on the emitter density (the total number of localized emitters divided by the occupied area of the structures). In a localization microscopy setup where the type of fluorescence emitters is already chosen, the main factor that determines the spatial resolution is thus the total number of localized emitters. Clearly, the key in developing fast localization microscopy is to obtain a sufficient number of localized emitters from a reduced total image acquisition time, which is the product of single frame time and the total number of image frames.

Generally two types of strategies were developed for fast localization microscopy: (1) minimizing single frame time while maintaining the same frame number, and (2) reducing frame number and increasing the total number of localized emitters from a single frame (thereafter called high-density localization microscopy). Using the first type of strategy, Zhuang et al reported a highest achievable imaging speed (0.5 s) in 2D super-resolution imaging with spatial resolution of ~25 nm [4]. For the second type of strategy, several methods have been reported [8–12], but all of them suffered from long data analyzing time. Note that a combination of the two strategies is applicable and will surely provide faster imaging speed. In any case, the goal of developing efficient methods for high-density localization microscopy is very promising, but still beyond reach.

In our previous paper [8], we proposed a method called SSM_BIC which is able to localize multiple emitters with high detection rate from low SBR (signal-background-ratio) images (common in localization microscopy). However, compared to other high-density localization methods, SSM_BIC suffers from tens or even hundreds of times slower in image analyzing which obstructs further uses of this method in localization microscopy. In this paper, we present a more practical method for processing high-density images. This method, which we have named PALMER (PArallel Localization of Multiple Emitters via Bayesian information criterion Recommendation), is the combination of GPU parallel computation, multiple-emitter fitting, and model recommendation via Bayesian Information Criterion (BIC). We verified that the PALMER method exhibits comparable localization accuracy with the SSM_BIC method, while executes >100 times faster. In other words, the PALMER method holds an image analyzing speed comparable to the fastest high-density localization method reported in the literature [9], without compromising the superior localization accuracy from its predecessor, the SSM_BIC method. Moreover, the capability of the PALMER method in analyzing experimental images was also validated.

2. Methods

2.1 General description of the PALMER Method

2.1.1 Strategy of the PALMER method

The PALMER method originates from the SSM_BIC method, and relies mainly on multiple-emitter fitting. Here the fitting model for multiple emitters is the same as that in the SSM_BIC and was described briefly in Section 2.1.2. To take advantage of GPU parallel computation to significantly reduce the entire execution time, the procedures in the SSM_BIC method were modified to have better compatibility with GPU computation. A major modification is that the Structured Sparse Model (SSM) in the SSM_BIC method was replaced with a modified threshold method in the pre-estimation for the initial multiple-emitters model (see detail information in Section 2.2). Other minor modifications for increasing the execution speed were documented in the Results and Discussions. The complete procedures of the PALMER method were described in Section 2.3.

2.1.2 Fitting model for multiple emitters

Point spread function (PSF) is widely used to describe the response of an optical system to individual fluorescence emitters at the focal plane, and can be well approximated by 2D Gaussian function [13, 14]. Assuming that each emitter contributes independently to observed signal at a given pixel, the imaging model for multiple (N) emitters is the convolution of PSF and the function with respect to the positions and the fluorescence intensities of the emitters. Using the 2D Gaussian function to approach PSF, the observed signal can be described as

2.2 Pre-estimation for the imaging model

A modified threshold method was used to pre-estimate the number and positions of candidate emitters in a sub-region after background subtraction. The sum signal from pixel(i,j) and its 4-neighborhood pixels (thereafter called 5-connecting-pixels), L_ij, was compared with a given threshold which is determined by PSF model and emitter brightness.

The method consists of the following steps: Step 1: The distribution of L_ij is calculated for all the pixels in the sub-region. Step 2: If the local maximum L_ij is larger than the given threshold, the center position of this 5-connecting-pixels is recorded into the list of candidate emitters. Step 3: Repeat Step 2 in other 5-connecting-pixels until no L_ij is found to be higher than the threshold or the maximum number of candidate emitters (which was set to be 6 to keep consistent with the SSM_BIC paper) is found. It is worthy of note that the pre-estimation works on the premise that there is no more than one emitter inside one 5-connecting-pixels. Therefore, once a candidate emitter is found, its 4-neighbor pixels will be excluded in the following pre-estimation steps.

The reliability of the pre-estimation depends on the proper setting of threshold value. We notice that (1) for a Gaussian PSF with fixed size, there is a fixed ratio between the sum signal from the 5-connecting-pixels centered at the PSF peak and the total signal from the emitter, and that (2) in a certain experimental condition, the average number of collected photons from the emitter (N_ave, which can be easily obtained by statistics) is also a fixed value. Therefore, in this paper we calculate the threshold value as the product of the ratio and N_ave.

2.3 Complete procedures of the PALMER method

A combined computational framework of CPU and GPU was used in the complete procedures of the PALMER method. As shown in Fig. 1 , five steps were performed successively to get precise emitter positions.

 

Fig. 1 The entire image analysis routine for the PALMER method. Note that multiple images are processed in parallel through this routine.

Download Full Size | PDF

Step 1. De-noising. Initial images were loaded into the CPU memory and then delivered to the GPU memory. The images were convoluted three times with a 3x3 pixels averaging template, and then with a 5x5 pixels annular averaging template [8]. Both convolution templates are separable, thus in this paper we divided them into two consecutive one-dimensional convolution operators, respectively. This treatment reduces the arithmetic complexity and bandwidth usage of the computation for each data point [16].

Step 2. Sub-region extraction. De-noised images were then delivered to the CPU memory. Using a simple threshold method, areas with 9 × 9 pixels centered at local maxima of de-noised images were extracted as sub-regions. The threshold value was usually set to be 5 times the standard deviation of the background noise in the de-noised images. This treatment guarantees that there is at least one active emitter inside a sub-region of 9 x 9 pixels. Once a sub-region was extracted, the signal in the corresponding pixels was removed from the de-noised image to avoid substantial overlapped sub-regions. This step is performed quite efficiently in CPU.

Step 3. Pre-estimation and localization. Sub-regions were sent to the GPU memory and assigned with separate threads. Pre-estimation on the number and positions of candidate emitters in a sub-region was performed in parallel according to Section 2.2. A series of imaging models were generated according to the pre-estimated information in different sub-regions (see Section 2.3 in the SSM_BIC paper for detail), and then Bayesian Information Criterion (BIC) statistics was used to recommend an optimal model. However, since the sub-regions may contain different number of emitters, workload varies in different threads. This significantly decreases the parallel computing performance of GPU. To avoid this, threads were organized logically according to workload, and only those threads with the same number of emitters in PSF model were activated in parallel at the same time point. The number of threads per block is mainly limited by registers per multiprocessor.

Step 4. Molecule clustering and sub-region refinement. The mean shift algorithm [17] for molecule clustering needs to update global positions after each iteration, and thus is not convenient to be implemented with GPU computation. Taking this into consideration, the recommended models obtained in Step 3 were sent back to the CPU memory for further molecule clustering and sub-region refinement. The performance of using CPU in this step is acceptably, evidenced by a small time consumption percentage of this step to the entire routine. It is worthy of note that the refined sub-region is expanded from the cluster center to include more emitters, and is terminated when all the emitters in the cluster are included or the box size reaches 13x13 pixels, which is limited by our GPU hardware platform. The information in the remaining pixels of the cluster will be lost after the refined sub-region reaches its maximum size.

Step 5. Refined sub-region localization. The refined sub-regions and the information for the associated emitters obtained from Step 4 (including the number and positions of emitters) were delivered to the GPU memory. The positions of the emitters can be refined through another round of multiple-emitter fitting and model recommendation, and then returned to the CPU memory for reconstructing a final super-resolution image. The time consumption of this step depends on the thread with largest sub-region size. The refinement in this step always presents better results.

2.4 Simulation data for image analysis

For all simulations, the pixel size at sample plane was set to be 100 nm and all emitted photons from individual fluorescence emitters were assumed to be detected completely by detector. To keep consistent with the SSM_BIC paper, the emitters were generated from a Gaussian PSF model whose standard deviation was set to be 1, so that each emitter occupies 5 × 5 pixels. We compared the effect of different emitter intensity models (same intensity vs Poisson distributed intensity) to the localization performance and found almost no statistics difference. Therefore in this paper we used a simpler model with fixed emitter intensity.

Simulated data sets were generated by placing emitters randomly or within a quarter-scale bitonal Siemens star pattern in an image of 100 × 100 pixels, which corresponds to an area of 100 μm². The emitter densities were controlled by placing different number of emitters within the active area in the image. The images with randomly distributed emitters were used for characterizing the performance of the methods in single molecule detection, while the images with Siemens star pattern were used for quantifying the resolving capability of the methods to fine structures.

2.5 Experimental data for image analysis

The single molecule imaging experiments were performed with a home-built TIRF microscope setup consisting of an inverted microscope (IX 71, Olympus), an oil immersion TIRF objective (100x/NA1.49, UAPON 100XOTIRF, Olympus) and an EMCCD camera (iXon 897, Andor). Fixed BSC-1 cells, where the actin bundles were labeled with fluorescence protein d2EosFP [18], were activated by a 405 nm laser diode and excited by a 561 nm diode pumped solid-state laser (both from CNILaser, China). Data were acquired by the EMCCD camera and the bundled software. The pixel size at sample plane is 160 nm and an Airy disk covers 3 × 3 pixels. Sample drift was not corrected, but was estimated to be less than 20 nm over the acquisition time with independent measurements.

2.6 Hardware platform for image analysis

All programs were written in C language under the Microsoft Visual Studio 2008 environment, using an Intel i5-2310, 2.90 GHz personal computer with 4.0 GB memory. An NVIDIA GeForce GTX 560 graphics card with 2.0 GB memory was used for the GPU-based localization computation. Note that the execution times of the methods depend on the specific hardware platform.

3.Results and discussions

3.1 Pre-estimation for the imaging model: modified threshold method vs SSM

A key step in the SSM_BIC method is the use of the Structured Sparse Model (SSM) to optimally pre-estimate the number and positions of candidate emitters in an extracted sub-region. Considering that active emitters are still sparse distributed in sub-pixel scale in high-density localization microscopy, the SSM_BIC method uses an optimization problem with inequality constraints to describe the sparse distribution characteristic. The inequality constraints correspond to the premise that there is no more than one active emitter in a 5-connencting-pixels. However, this step is the most time-consuming part in SSM_BIC, since solving the optimization problem needs significant computational effort. Therefore in this study, instead of SSM, we adopted a modified threshold method with which we can make a good balance among execution speed, effectiveness in separating signal and noise, and achievable emitter density. We found out that the modified threshold method provided quite acceptable pre-estimation results for further steps, evidenced by the finding that the final localization accuracy obtained from the PALMER method is comparable to that from the SSM_BIC method (Fig. 2 ).

 

Fig. 2 Localization performance of the PALMER, SSM_BIC and DAOSTORM methods in analyzing simulated images with different emitter densities. Note that the background was 100 photons per pixel for all simulations, while the signals were for the total detected photons from individual emitters. Each date point was averaged from five independent measurements. The error bars in (e) and (f) are standard deviations.

Download Full Size | PDF

The effectiveness of the modified threshold method comes from two facts: (1) the premise used in SSM is also effectively taken into account in this threshold method (See Section 2.2), and (2) the pre-estimated information (the number and positions of active emitters) obtained from this modified threshold method are suboptimal solvers of the optimization problem in SSM and the inaccuracy of these information can be compensated from subsequent analysis steps. Moreover, the use of 5-connecting-pixels is appropriate for Gaussian PSF models of 5 × 5 pixels or smaller. For such PSF models, the L_ij at the peak pixel is the largest and occupies the major PSF signal (> 50%), therefore 5-connecting-pixels could be effectively used to separate signal from noise with high robustness.

3.2 Localization accuracy

We quantified the localization accuracy of the PALMER method by analyzing simulated images with randomly distributed emitters. In the simulation, image data sets with two different SBR were generated. Specifically, images with a total detected signal of 350 photons from each emitter were defined as low SBR images to keep consistent with of the SSM_BIC paper, while those images with a total detected signal of 800 photons from each emitter, which is typically obtained in our experiments, were referred to high SBR images. The background was 100 photons per pixel in both cases. The number of emitters in the images varies from 20 to 500, corresponding to image densities from 0.2 to 5 emitters/μm². The localization accuracy was characterized with three popular parameters, including detection rate, false-positive rate and localization precision, according the definitions listed below.

In the previous paper, we verified that the SSM_BIC method exhibits superior localization accuracy than the Sparse Algorithm (SA) in analyzing images with high-density emitters, especially for those images with low SBR. Interestingly, the localization accuracy of the PALMER method is comparable to that of the SSM_BIC method (Fig. 2), especially when the emitter densities are less than 2 emitters/μm². For image densities of 1 molecule/μm² or less, which are popularly regarded as low density [9, 10], the PALMER method exhibits satisfactory localization accuracy, evidenced by detection rate of >85%, false-positive rate of <10% and localized precision (σ) of ~20 nm even in low SBR condition (that is, a Signal/Background of 350/100). For image densities of 3 emitters/μm² or higher, the localization accuracy of the PALMER method deviates slightly from that of the SSM_BIC method, which is due to the size limitation of the refined sub-region (13 × 13 pixels in this study). As mentioned before, the signal beyond the maximum refined sub-region size is not maintained. This treatment happens more frequently with the increase of emitter density, and thus leads to larger deviation in localization accuracy between the PALMER and the SSM_BIC method.

It is worth mentioning that the localization accuracy of the PALMER method depends on the SBR level. For the low SBR case, the detection rate drops to <50% (Fig. 2(a)) and the false-positive rate is ~30% (Fig. 2(c)), when the emitter density increases to 5 emitters/μm². However, at the same emitter density, high SBR images exhibit 20% higher detection rate (Fig. 2(b)) and 5% lower false-positive rate (Fig. 2(d)). This finding suggests that, increasing the image SBR is of great beneficial for high-density localization microscopy, since the permitting emitter density can be much higher with the increase of SBR.

With the same sets of simulated data, we also compared our PALMER method with the previous reported DAOSTORM method [10]. The Python codes provided by the DAOSTORM authors were tested to be executed properly in our computer hardware. We found that the PALMER and SSM_BIC methods show better localization accuracy performance than the DAOSTORM method (Fig. 2). Note that in these comparisons, real positions are defined as localized positions located in less than 50 nm from true positions, while in the original DAOSTORM paper, they are more loosely defined (3 PSF widths, ~300 nm or larger).

3.3 Computational speed

We quantified the computational time of the two methods in the entire image analysis routine (Fig. 3(a) ). To increase GPU computation efficiency, the frame numbers for image stacks with different emitter density were optimized, and fewer frame were used with higher emitter densities. Simulated results show that the PALMER method processes single image frame at two orders of magnitude faster speed than the SSM_BIC method (Fig. 3). Under low emitter densities (1 molecule/μm² or less), the speed gain is >240.

 

Fig. 3 Time consumptions of the PALMER and SSM_BIC methods in analyzing single image frame. Simulated images of 100 × 100 pixels were used, where the total detected signal from individual emitters was 350 photons and the background was 100 photons per pixel. Each data point was calculated from five sets of data.

Download Full Size | PDF

The high speed gain comes mainly from a combination of three factors (listed in order of importance): GPU parallel computation, simplified threshold method in pre-estimation, and basic gradient projection method in emitter localization. First of all, GPU parallel computation was used in all of the time-consuming steps, including de-noising, pre-estimation and localization (Fig. 1). For the remaining steps, sub-region extraction can be performed efficient enough in CPU, while clustering is not convenient to be performed in GPU (see Section 2.3). Then, the simplified threshold method was effectively used to replace the extremely time-consuming fitting-based SSM method. Finally, the total workload of localization was reduced by replacing the generalized projection gradient method [19] with the basic gradient projection method [15], and optimizing the fitting iteration to be 20. Note the SSM_BIC method used in this study also employed the basic gradient projection method, which improves the computational speed by a factor of 2 while maintains the same localization accuracy.

3.4 Resolving capability

We estimated the capability of the PALMER method in resolving fine structures using simulated images with quarter-scale bitonal Siemens star pattern, where 11 line pair structures exist in the pattern and the average line pair separation varies from 18 nm to 350 nm (Fig. 4(a) ). A total number of 1000 image frames with 500 emitters inside each frame were used, corresponding to image density of 10 emitters/μm². Note that the pattern occupies a total area of 50 μm². The SBR of the simulated images was chosen to be compatible with our experimental results, that is, a total detected signal of 800 photons per emitter and a background of 100 photons per pixel (Fig. 4).

 

Fig. 4 Evaluation of the resolving capability of the PALMER method using simulated images with quarter-scale bitonal Siemens star pattern. (a) Image with ground-truth pattern. Wide-field image averaged from 1000 frames (b), a typical image frame (c) and a reconstructed super-resolution image from the PALMER method (d). The normalized intensity profile (green dots, averaged from 11 line pairs) of the green arc in (d) is shown in (e), which is further fitted with rectangular function (purple). (f) The dependence of the contrast on the structure width. The error bars in (f) show the standard errors. Note that the structures within the red dash box in (f) are irresolvable.

Download Full Size | PDF

We analyzed the average intensity profile of line pairs at different widths (Fig. 4(d)) and found that, when the width increases to 48 nm, the corresponding intensity profile exhibits clearly peaks whose positions match with the real positions. To make the resolving capability more convenient to quantify, the intensity profile was normalized and fitted with a rectangular function (Fig. 4(e)). Contrast, which was defined as (I_max - I_min)/(I_max + I_min), can then be used to quantify the resolving capability. For the smallest resolvable structure of 48 nm, the contrast was calculated to be 0.22. Therefore, any sstructures with contrast lower than 0.22 (indicated by the red dash box in Fig. 4(f)) were considered to be irresolvable. The dependence of the contrast on the structure width was calculated and shown in Fig. 4(f). Clearly from Fig. 4(f), the contrast begins to decrease at a width of 250 nm, and reaches the resolving limit at a width of 48 nm (see the red dash box).

3.5 Analysis of experimental data

It would be beneficial to validate the power of the PALMER method in analyzing experimental images. Here the localization performance of the PALMER method was compared with the MaLiang method, an Sparse Algorithm described in detail in reference [20]. TIRF images with high (Fig. 5(a) ) and low (Fig. 5(f)) emitter densities were obtained by stacking a total number of 2000 experimental image frames, respectively. Using the two image analyzing methods, the TIRF images with high emitter densities were processed to obtain super-resolution images (Fig. 5(c) and Fig. 5(d)). Analysis on a representative fine structure in the super-resolution images clearly shows that the feature of the structure can be better resolved with the PALMER method (Fig. 5(e)). For images with low emitter densities, the PALMER method presents comparable or even better localization performance than that from the MaLiang method (Fig. 5(h-j)).

 

Fig. 5 Localization performance of the PALMER and MaLiang methods in analyzing experimental images with high (a-e) and low (f-j) emitter density. TIRF image from a stack of 2000 image frames (a, f) and a single image frame (b, g). Super resolution images were reconstructed by the PALMER (c, h) and the MaLiang (d, i) method, respectively. Cross-sectional distributions of emitter localizations along the direction of the actins in the boxed region are shown (e, j), where the colors in the histograms are the same as those in the corresponding images. Scale bar: 2 μm.

Download Full Size | PDF

On the other hand, it was reported that the Nyquist resolution can be used to present a good quantification on the smallest resolvable structures from a whole reconstructed super-resolution image. According to Shroff et al [6], the Nyquist resolution, Δ_Nyquist = 2/(N_emitter/S)^1/D, depends on the total number of localized emitters (N_emitter), the occupied area of the biological structures (S), and the dimension (D) of the structures. Here, D = 2 for two dimensional images, S was calculated to be 13.7 μm² from Fig. 5(c). The total number of localized emitters (N_emitter) for different images and the corresponding Nyquist resolutions were presented in Table 1 . For the same stack of TIRF images with high emitter densities, the Nyquist resolution obtained from the PALMER method (21 nm) is ~two folds higher than that from the MaLiang method (40 nm). For images with low emitter density, similar findings were obtained (Table 1).

Table 1. Nyquist resolution from reconstructed super-resolution images

View Table

Finally, it is noteworthy that, fast localization microscopy is highly possible from the strategy of combining high-density localization microscopy with practical high-density localization method, evidenced by the finding that within the same amount of acquisition time, the total number of localized emitters from this strategy is ~14 times (119729) that from the strategy of combining low emitter density image acquisition with sparse localization method (8356). In other word, compared to the conventional localization microscopy based on sparse emitter localization, this high-density localization approach allows a speed gain of up to ~14-fold in obtaining a super-resolution image with the same Nyquist resolution.

4. Conclusion

Our PALMER method, based on the combination of GPU parallel computation, multiple-emitter fitting and model recommendation via BIC, is capable of fast analyzing high-density single molecule images. The PALMER method executes more than two orders of magnitude faster than the previously reported SSM_BIC method for analyzing simulation images without obviously compromising the localization accuracy. Through experimental study, we validated that the strategy of combining high-density localization microscopy with practical high-density localization method contributes more than an order of magnitude faster imaging speed for a desired Nyquist resolution than the strategy of combining conventional localization microscopy with sparse emitter localization method. We believe that this PALMER method has great potentials in pushing forward the field of fast localization microscopy.