# Drift correction for single-molecule imaging by molecular constraint field, a distance minimum metric

- Renmin Han
^{1, 2}, - Liansan Wang
^{3}, - Fan Xu
^{1, 2}, - Yongdeng Zhang
^{4}, - Mingshu Zhang
^{4}, - Zhiyong Liu
^{5}, - Fei Ren
^{1}Email author and - Fa Zhang
^{1}Email author

**8**:1

**DOI: **10.1186/s13628-014-0015-1

© Han et al.; licensee BioMed Central. 2015

**Received: **22 October 2014

**Accepted: **26 November 2014

**Published: **13 January 2015

## Abstract

### Background

The recent developments of far-field optical microscopy (single molecule imaging techniques) have overcome the diffraction barrier of light and improve image resolution by a factor of ten compared with conventional light microscopy. These techniques utilize the stochastic switching of probe molecules to overcome the diffraction limit and determine the precise localizations of molecules, which often requires a long image acquisition time. However, long acquisition times increase the risk of sample drift. In the case of high resolution microscopy, sample drift would decrease the image resolution.

### Results

In this paper, we propose a novel metric based on the distance between molecules to solve the drift correction. The proposed metric directly uses the position information of molecules to estimate the frame drift. We also designed an algorithm to implement the metric for the general application of drift correction. There are two advantages of our method: First, because our method does not require space binning of positions of molecules but directly operates on the positions, it is more natural for single molecule imaging techniques. Second, our method can estimate drift with a small number of positions in each temporal bin, which may extend its potential application.

### Conclusions

The effectiveness of our method has been demonstrated by both simulated data and experiments on single molecular images.

### Keywords

Image reconstruction techniques Resolution Fluorescence microscopy Superresolution## Background

Far-field optical microscopy has played an important role in the biological sciences. However, the imaging capacity of a conventional light microscope is fundamentally limited by the wavelength of the light. Recent developments in optical microscopy, such as fluorescence photoactivation localization microscopy(FPALM/PALM) [1] and stochastic optical reconstruction microscopy(STORM) [2,3], have overcome the diffraction limit [4] and achieved resolution values down to 20 nm, a factor of ten lower than the resolution of conventional light microscopy. These techniques exploit the stochastic photoswitching of fluorescent probe molecules and carefully choose imaging parameters to guarantee that no molecules closer together than the point spread function (PSF) width would emit signal at the same time. Because these techniques should readout a sequence of camera frames in a diffraction limited situation for the reconstruction of a super-resolution image, it is common that the image acquisition will take a long time. Typically, 1,000-100,000 individual camera frames [1,3,5,6] are needed, with exposure times of 1-100 ms per frame; this translates to several minutes to hours in total.

Long acquisition times increase the risk of sample drift during imaging [7]. Sample drift on the nanometer scale is hard to avoid as it may be caused by a variety of reasons, such as mechanical instability of instrument and vibration. Without correction, the lateral drift will smear the image and degenerate the resolution. The drift correction problem is a special case of image registration that has been discussed for decades in medical image processing and leads into a variety of solutions. Recently, several methods have been proposed to solve the drift correction problem in the field of single molecule imaging. Fiducial markers, typically gold nanoparticles or fluorescent beads, can be incorporated into the sample. Because the brightness of fiducial markers has a high contrast to that of probe molecules [1,3] and is stable during measurement, the positions of markers can be easily tracked and their trajectory can be used to correct the drift, which makes the method quite reliable. However, the fiducial marker-based method needs beads introduced into the sample and limits our observation area. Fiducial markers also may interfere with the sample and degenerate the localization precision of the nearby probe molecules [8]. To overcome the limitation of this marker-based method, cross-correlation is an alternative option. Cross-correlation utilizes the biological structure itself, which is stationary during the measurement, to estimate the drift. Though cross-correlation avoids the problems of introducing fiducial markers, there are also several shortcomings. Generally, the positions of molecules should be spatially binned as a 2D histogram to carry out the cross-correlation analysis [9,10]. The behavior of cross-correlation is close relative to the localized fluorescent position density of the 2D histogram. There are two decisive factors: the parameters of spatial binning and temporal binning. First, we should carefully choose the spatial binning parameter. If the spatial binning is too coarse, i.e., if the pixel width of the 2D histogram is too wide, smaller features of the molecular distribution will be lost. If the density of molecules in a pixel is too low, the histogram will not represent the global distribution of molecules, which may lead to the failure of the cross-correlation analysis. Second, the choice of the temporal binning parameter also determines the number of localized positions in a temporal bin. Typically, the number of positions in a temporal bin should be sufficient to give a meaningful description of molecule distribution. However, for applications where there are a small number of localized florescent events in each frame, the temporal bins will be too large to ignore the drift within one interval, which leads to an inaccurate estimation of the overall drift.

In this paper, we propose a new metric based on the distance between molecules to solve the drift correction. The proposed metric directly uses the position information of molecules to estimate the frame drift. Thus, the introduction of fiducial markers and the problems of binning the frames into 2D histogram are avoided. Since the metric directly operates on the positions, it is more natural and conceptually straight-forward for the single molecule imaging techniques. Though the choice of parameters of temporal binning is still a problem in our metric, by directly using the distance metric, just a small amount of information on the stable molecules can uncover the relationship between different frames. Therefore, the proposed method could estimate the drift with a small number of positions in each temporal bin, which may extend the potential application scenarios of our metric. We also designed an incremental algorithm to implement the metric for the general application of drift correction. Different experiments based on simulated data and real data have been carried out and the results demonstrate the effectiveness of our method.

## Methods

The stochastic switching of probe molecules is the fundamental mechanism of single molecule imaging. The switching on and switching off phenomenon of a single molecule obeys a Poisson process in the temporal domain [10]. Because the switching events of probes are asynchronous, the observation (frame) at different time points becomes a random sampling of the molecule distribution in the structure. In addition, considering the difference between the readout rate of microscopy and switching rate of probe, a florescent event observed in a frame may change or not be recorded in the next frame. Though the molecule position can be determined by various methods [11], the essence of these localization methods is fitting the distribution of the measured photon spot to that of the ideal point spread function (PSF) [12], which will unavoidably introduce localization error. Therefore, it is difficult to trace an identical molecule in different frames and we will not make any assumption about the traceability of probe molecules in our model.

The super-resolution image is reconstructed from the molecule positions. The localization precision of probe molecules has a great effect on the reconstruction quality. There are three kinds of uncertainties in the localization of probe molecules. The first uncertainty comes from the localization technique itself. The fitting methods are based on a probability model and certain assumption about the ideal PSF, but the fundamental assumption may not be completely consistent with the emitter properties in each application. In addition, the observed photon distribution may be affected by background noise or mutual interference of molecules, all of which affect the localization accuracy. The second uncertainty comes from the thermodynamic instability of probe molecules. The long acquisition of single molecule imaging gives a molecule enough chance to undergo random motion. The effect of random motion is intrinsic and inevitable in the period of measurement. The third uncertainty is the sample drift. These drifts will degenerate a reconstruction image to a great degree. Fortunately, they can be compensated by the drift correction method.

### Imaging model

*D*is the sample drift. In general, the calculation error obeys a Gaussian distribution, without considering the model bias. The probable positions of a probe molecule in random motion are also assumed to follow a Gaussian distribution, according to the most recent study and experiment by Nieuwenhuizen [13]. For convenience of discussion, we ignore their difference in mechanism and define them as a combination \(\vec {\phi }\), i.e. \(\vec {\phi }=\vec {\phi }_{\textit {fit}}+\vec {\phi }_{\textit {ran}}\) (\(\vec {\phi }\) also obeys a Gaussian distribution). That is to say, if we already know the ideal position \(\vec {p}_{\textit {ora}}\) of a florescent event and the drift

*D*, the probability that a measured position

*p*corresponds to the ideal position is:

where *σ* is the standard deviation of \(\vec {\phi }\).

where *δ*(·) is the Dirac delta function. The Dirac delta function can depict the signal distribution exactly, but there comes a problem when processing molecule position sets in which the sample drift has not been corrected. The Dirac delta function can not be used to judge the similarity between two position sets, because it does not formulate the uncertainty in the measurement. Here, we propose a novel metric to solve the problem. Our metric is based on the distance between molecules and can directly measure the similarity of signals by molecule positions.

- 1.
For different position sets that depict the identical biological structure, the distribution of their superset should be as sharp as possible.

- 2.
For a given structure (distribution) that is depicted by molecule positions, the more positions there are, the better the signal strength output will be.

*p*, the function is given by:

*l*2−norm and \(\vec {p}_{j}\) represents the point located in the region. The certain form of the kernel function FUN(·) has not been given out here. FUN(·) should be a function symmetric about the

*y*-axis and depends on the real application. The summation of kernel function maps the distance between two points into a value describing the sparseness. According to the discussion of uncertainty of the image model, Gaussian function is chosen as the kernel function. Therefore, the generalized PMCF of position

*p*is defined as follows:

*σ*restricts the effective width of kernel function. It is not hard to extend the idea from a molecule position to a position set. Assuming there are two position sets Φ

_{ A }and Φ

_{ B }, the Molecular Constraint Field (MCF) of the given two position sets is defined as:

where \(\vec {p}_{i}\in \varPhi _{A}, i=1,2,\ldots N\) and \(\vec {p}_{j}\in \varPhi _{B}, j=1,2,\ldots M\). MCF is a metric based on the conception of position sets and directly outputs the similarity of two signal distributions represented by the Dirac delta function. Note that the properties of MCF are coincident with the above basic rules. The conception of MCF is similar to cross-correlation analysis and is connected with the concept of the Parzen-Window Density Estimation [14,15]. However, MCF directly operates on the molecule positions to judge the similarity of two signal distributions, which is more conceptually straight-forward and easy to operate compared with the cross-correlation analysis of a 2D histogram.

- 1.
Reciprocity: If \(\vec {p}_{i}\) contributes value to the PMCF of \(\vec {p}_{j}\), \(\vec {p}_{j}\) will contribute equal value to the PMCF of \(\vec {p}_{i}\).

- 2.
Additivity: Every molecule position contributes to the PMCF of a certain region. If there are several difference regions, the molecule will have contribution to each.

- 3.
Regionally restraint: The output of the kernel function of MCF declines to zero when the value of the distance exceeds a threshold.

### Drift correction

The super-resolution image is reconstructed from a series of frames. For convenience in our later discussion, we define the positions of all the frames, i.e., the positions composing the super-resolution image, by Φ_{
all
}, and define the positions of the *i*th frame by Φ_{
f
r(i)}. Therefore, Φ_{
f
r(i)} is a subset of Φ_{
all
}, i.e., ∀*i*,Φ _{
f
r(i)}⊆Φ _{
all
}. The position sets of several successive frames can be combined (temporally binned) into a union. If the frames are binned at an interval of *t*, we defines the union set of the *i*th interval by Φ_{
T(i)}, i.e., \(\varPhi _{T(i)}=\varPhi _{fr(i\cdot t-t+1)}\bigcup \varPhi _{fr(i\cdot t-t+2)}\bigcup \ldots \bigcup \varPhi _{fr(i\cdot t)}\).

_{ A }and Φ

_{ B }(\(\vec {p}_{i}\in \varPhi _{A}, i=1,2,\ldots N\) and \(\vec {p}_{j}\in \varPhi _{B}, j=1,2,\ldots M\)), the cost function

*M*

*C*

*F*(

*L*) of drift correction is given by:

where *L*(·) is a transformation function, Φ_{
B
} is movable and Φ_{
A
} is fixed. Our aim is to find an *L*(·) that maximizes the MCF. The choice of the shape of *L*(·) depends on the method of image acquisition. In most applications, the drift is assumed to be linear [9,10], so we will make that assumption (though affine projection may occur in some cases). Therefore, in the following contents, the transformation is defined as linear, i.e., \(L(\vec {p})=\vec {p}+\vec {d}\), where \(\vec {d}=(\alpha,\beta)^{T}\), *α* is the drift compensation of the *x*-axis and *β* is the drift compensation of the *y*-axis.

*t*. Our algorithm takes {Φ

_{ T(i)}} and an empty reference set Φ

_{ ref }as input. The reference position set Φ

_{ ref }is initialized by Φ

_{ T(1)}and the counter

*i*is set to two. Then position set Φ

_{ T(i)}is selected to compare with Φ

_{ ref }and the sample drift is estimated. The drift of Φ

_{ T(i)}will be compensated and the positions of Φ

_{ T(i)}will be merged into Φ

_{ ref }. Repeat the selection and drift compensation process until all the subset are corrected. Finally, the reference position set Φ

_{ ref }will be output as the data for the super-resolution image reconstruction. Many techniques could be utilized to solve the maximum problem, for example, Grid Searching, Conjugate Gradient. We have implemented the algorithm in C++, based on GNU Scientific Library (GSL); it is freely available from the authors upon request.

In the process of drift estimation between Φ_{
T(i)} and Φ_{
ref
}, all the information from the previously corrected position sets are used. Therefore, according to the study of Geisler [10], the deviation of drift estimation error can be controlled by this incremental strategy of our algorithm. Two parameters, the *σ* of MCF and the temporal binning interval *t*, should be initialized. According to the Nyquist sampling theorem, the value of *σ* should be more than double the standard deviation of the molecular localization uncertainty. The choice of temporal binning interval *t* is less restricted than that in the method of cross-correction analysis. In the following experiments, we prove that our method can achieve high drift correction accuracy with very few molecule positions in each Φ_{
T(i)}, assuming the sampling has no bias.

## Results

### Simulated data

*T*=40 equal step length bins (i.e., there are 40 groups of molecules and in each group there are 1000 molecules) with drift imposed along the x-axis. The linear drift in time

*t*

_{ j },

*j*=1,2,…

*T*is illustrated in Figure 3(b) and the structure after the drift is shown in Figure 3(a). The second structure, named ‘Grid’, consists of 60,000 molecules distributed on a grid structure, the thickness of the lines of the grid is 120 nm. The molecules are divided into

*T*=40 equal step length bins (i.e., there 40 groups of molecules and in each group there are 1500 molecules) with drift imposed along x-axis and y-axis. The drift of x-axis and y-axis in time

*t*

_{ j },

*j*=1,2,…

*T*is illustrated in Figure 3(d) and the structure after the drift is shown in Figure 3(c). The third structure, named ‘Radio’, consists of 60,000 molecules labeled on a structure composed of two ring structures and several lines radiating from the center, where the thickness of the ring structure is 120 nm and the thickness of the line is 80 nm. The molecules are divided into

*T*=40 equal step length bins (i.e., there are 40 groups of molecules and in each group there are 1500 molecules) with drift imposed along x-axis and y-axis. The drift of x-axis and y-axis in time

*t*

_{ j },

*j*=1,2,…

*T*is illustrated in Figure 3(f) and the structure after drift is shown in Figure 3(e).

*n*) and localization precision (denoted by

*ξ*) of data set Grid and Radio were carried out to illustrate the impact. The mean

*Δ*

_{ d }and the standard deviation

*σ*

_{ d }of residual over the

*T*=40 temporal intervals were calculated for each simulation. Figures 7 and 8 illustrate the curves of

*Δ*

_{ d }and

*σ*

_{ d }as a function of

*n*and

*ξ*. For both Figures 7 and 8, sub-fig (a) (c) represent the same

*Δ*

_{ d }and sub-fig (b) (d) represent the same

*σ*

_{ d }, in two different views. The tendencies of the mean residual and the standard deviation of the residual under changes in

*n*and

*ξ*are consistent. Judging from Figure 7(a) and (b), we can find that our method behaves well in the situations in which the localization precision

*ξ*≤40 for data set Grid. Under the condition that localization precision

*ξ*≤80 and molecule numbers per bin

*n*≥500, the accuracy of drift compensation by our method can reach to 2∼10 nm. Comparing Figure 8(a) (b) to Figure 7(a) (b), we find that the performance of our method in data set Radio is not as good as that in Grid when

*ξ*=80. The reason why there is such a difference is that when the localization precision reaches

*ξ*=80, the localization uncertainty is so large compared to the width of the line structure in Radio (80 nm) that it degenerates the sampling of the molecule distribution. Our method performs not good when

*ξ*=120 nm (the width of line in Grid) and the molecule number per bin

*n*is very small. The performance of our method is impacted by bias sampling or mal-sampling, because our method is based on the direct recovery of position distributions. Nevertheless, according to recent studies [1,13], localization precisions are usually controlled to 10∼35 nm in 2D, which makes the low localization precision not a problem for our method. Judging from sub-fig (c) and (d), we can find that

*n*=500 is sufficient to support a good performance of our method. Additionally, if the localization precision

*ξ*≤20 nm, judging from Figures 7 and 8, we can conclude that our method can compensate the sample drift to a residual of 20∼40 nm with 40 molecule positions per temporal bin and to a residual of about 15 nm with 100 molecule positions per temporal bin. Readers who are interest in the details of the performance of our method in the cases where

*n*=40 and

*ξ*=0 nm can refer to the Appendix.

### Super-resolution experimental data

To confirm the effectiveness of our drift correction method for practical applications, real super-resolution experimental data is utilized.

COS-7 cells were cultured in DMEM complete medium (Gibco) supplemented with 10*%* fetal bovine serum and maintained at 37°C and 5*%* CO2 in a humidified incubator (Thermo). The cells were fixed with 3*%* (w/v) paraformaldehyde and 0.5*%* glutaraldehyde in PBS for 15 to 40 min at 37°C and washed 3 to 5 times with filtered PBS. After that, the fixed COS-7 cells were permeabilized for 10 min with 0.1*%* Triton X-100 and then blocked in 5*%* bovine serum albumin (BSA, AMRESCO) diluted in PBS for 60 min. The mouse anti- *β* tubule monoclonal antibody (Sungenebiotech) was diluted 1:200 in PBS containing 2.5*%* BSA and added to the COS-7 cells in 37°C incubator for 60 min. After three times rinsing with PBS, COS-7 cells were then incubated for another 60 min with the Alexa Fluor 647 labeled rabbit anti-mouse secondary antibody (Invitrogen), which was also diluted 1:200 in PBS. At last, the cells were rinsed four times with PBS and kept in a dark place. The STORM imaging buffer contained imaging buffer base (10*%* glucose (m/v), 50 mM Tris (pH 8.0) and 10 mM NaCl), an oxygen scavenger system (0.5 mg ml ^{−1} glucose oxidase (Sigma-Aldrich), 40 *μ*g ml ^{−1} catalase (Sigma-Aldrich)) and 10 mM MEA. Two fluorescent beads were embedded on the surface of the sample.

STORM imaging of microtubules was performed as described in [16]. We used an Olympus IX71 inverted microscope equipped with a 150×1.45 numerical aperture (NA) oil objective (Olympus PLAN APO). Two lasers (405 nm and 647 nm (OBIS, Coherent)) were controlled by an acousto-optic tunable filter (AA Optoelectronic). For excitation, the power of the 647 nm laser was 29.71 mW, measured near the rear pupil of the objective. The intensity of the 405 nm laser, typically 10-30 *μ*W, was adjusted so that a low density of molecules was activated at each frame. A *λ*/4 plate was used to produce circular polarization excitation light. The fluorescence signals were acquired using an electron-multiplying charge-coupled device (EMCCD) camera (Andor iXon DU-897). The images were acquired at a frame rate of 50 Hz and the EM gain of EMCCD was set to 300. Fluorescent beads were embedded into the sample for our convenience in comparing them. The super-resolution image reconstruction was performed as described in [17]. In order to minimize the influence of the background, TIRF (total internal reflection) illumination was used in this study.

## Discussion and conclusion

In this paper, we proposed a new metric based on distance minimums to cope with the drift correction in single-molecule imaging. Our method is based on the restoration of the distribution of molecular position. There are two advantages of our method. First, because our method does not require space binning of positions of molecules but directly operates on the positions, it is more natural to single molecule imaging techniques. Second, our method estimates the drift with a small number of positions in each temporal bin, which may extend the potential application of our metric. The result of simulated test data proved the effectiveness of our method. Our method can be easily extended into 3D super-resolution photography, and the Gaussian kernel can also be replaced in some special situations.

## Appendix

## Declarations

### Acknowledgements

This work was supported by grants National Natural Science Foundation of China (61232001, 61202210, 61472397, 31270910 and 31170818); Beijing Natural Science Foundation (7131011).

## Authors’ Affiliations

## References

- Betzig E, Patterson GH, Sougrat R, Lindwasser OW, Olenych S, Bonifacino JS, Davidson MW, Lippincott-Schwartz J, Hess HF. Imaging intracellular fluorescent proteins at nanometer resolution. Science. 2006; 313(5793):1642–1645.ADSView ArticleGoogle Scholar
- Hess ST, Girirajan TPK, Mason MD. Ultra-high resolution imaging by fluorescence photoactivation localization microscopy. Biophys J. 2006; 91(11):4258–4272.View ArticleGoogle Scholar
- Rust MJ, Bates M, Zhuang X. Sub-diffraction-limit imaging by stochastic optical reconstruction microscopy (storm). Nat Methods. 2006; 3(10):793–796.View ArticleGoogle Scholar
- Hell SW. Far-field optical nanoscopy. Science. 2007; 316(5828):1153–1158.ADSView ArticleGoogle Scholar
- Egner A, Geisler C, Von Middendorff C, Bock H, Wenzel D, Medda R, et al. Fluorescence nanoscopy in whole cells by asynchronous localization of photoswitching emitters. Biophys J. 2007; 93(9):3285–3290.View ArticleGoogle Scholar
- Bates M, Huang B, Zhuang X. Super-resolution microscopy by nanoscale localization of photo-switchable fluorescent probes. Curr Opin Chem Biol. 2008; 12(5):505–514.View ArticleGoogle Scholar
- Gould TJ, Verkhusha VV, Hess ST. Imaging biological structures with fluorescence photoactivation localization microscopy. Nat Protoc. 2009; 4(3):291–308.ADSView ArticleGoogle Scholar
- Zessin P, Kruger C, Malkusch S, Endesfelder U, Heilemann M. A hydrophilic gel matrix for single-molecule super-resolution microscopy. Opt Nanoscopy. 2013; 2(1):4.View ArticleGoogle Scholar
- Mlodzianoski MJ, Schreiner JM, Callahan SP, Smolková K, Dlasková A, Šantorová J, Ježek P, Bewersdorf J. Sample drift correction in 3d fluorescence photoactivation localization microscopy. Opt Express. 2011; 19(16):15009–15019.ADSView ArticleGoogle Scholar
- Geisler C, Hotz T, Schönle A, Hell SW, Munk A, Egner A. Drift estimation for single marker switching based imaging schemes. Opt Express. 2012; 20(7):7274–7289.ADSView ArticleGoogle Scholar
- Deschout H, Zanacchi FC, Mlodzianoski M, Diaspro A, Bewersdorf J, Hess ST, Braeckmans K. Precisely and accurately localizing single emitters in fluorescence microscopy. Nat Methods. 2014; 11(3):253–266.View ArticleGoogle Scholar
- Thompson RE, Larson DR, Webb WW. Precise nanometer localization analysis for individual fluorescent probes. Biophys J. 2002; 82(5):2775–2783.View ArticleGoogle Scholar
- Nieuwenhuizen RP, Lidke KA, Bates M, Puig DL, Grünwald D, Stallinga S, Rieger B. Measuring image resolution in optical nanoscopy. Nat Methods. 2013; 10(6):557–562.View ArticleGoogle Scholar
- Rosenblatt M. Remarks on some nonparametric estimates of a density function. Ann Math Stat. 1956:832–837.
- Parzen E. On estimation of a probability density function and mode. Ann Math Stat. 1962; 33(3):1065–1076.View ArticleMATHMathSciNetGoogle Scholar
- Dempsey GT, Vaughan JC, Chen KH, Bates M, Zhuang X. Evaluation of fluorophores for optimal performance in localization-based super-resolution imaging. Nat methods. 2011; 8(12):1027–1036.View ArticleGoogle Scholar
- Zhang M, Chang H, Zhang Y, Yu J, Wu L, Ji W, et al. Rational design of true monomeric and bright photoactivatable fluorescent proteins. Nat. 2012; 9(7):727–729.Google Scholar

## Copyright

This is an Open Access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/2.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly credited.