A mean shift algorithm for drift correction in localization microscopy

Single-molecule localization microscopy techniques transcend the diffraction limit of visible light by localizing isolated emitters sampled stochastically. This time-lapse imaging necessitates long acquisition times, over which sample drift can become large relative to the localization precision. Here, we present an efficient and robust method for estimating drift, using a simple peak-finding algorithm based on mean shifts that is effective for single-molecule localization microscopy in two or three dimensions.


INTRODUCTION
Stochastic super-resolution microscopy techniques, such as STORM (1,2) and PALM (3,4), exploit photoswitching of fluorescent probes to enable imaging of densely labeled samples with resolutions an order of magnitude smaller than the diffraction limit of visible light. Sparsely distributed point spread functions (PSFs) of single emitters are identified in individual image frames, and their centroids are determined according to an appropriate fitting algorithm. The axial position of molecules can be encoded in their PSFs through engineering measures utilizing astigmatism (5,6), multifocal plane imaging (7), or a double helix PSF (8). The final reconstruction is typically a twodimensional (2D) or three-dimensional (3D) histogram of these single-molecule positions.

RESULTS
A graphical illustration of the MS algorithm as applied to sample 2D localizations is presented in Fig. 1. The localizations all lie in one of two data sets which sample the same uniformly distributed emitters, but with a constant relative shift r shift in space. The first step of the algorithm is to extract pairwise displacements between all localizations across the two data sets. When individual displacements are plotted as points ( Fig. 1  b), displacements arising from the same labeled objects (magenta points) cluster around r shift , whereas displacements arising from different objects (green points) distribute randomly over space. The MS algorithm determines the center of the peak of the distribution through iteration (27)(28)(29). At each iteration, all pairs within the radius of consideration are extracted, and the updated shift estimate is the centroid of these pairs. The uniformly distributed background will tend to bias the centroid toward the center of the observation window, whereas the peak moves the mean toward r shift . The observation window is then redrawn around the new mean and the process is repeated until the peak is centered in the observation window. Three iterations of the algorithm are visualized in Fig. 1 c.
Although the emitters of Fig. 1 are distributed uniformly in space, leading to the uniform distribution of the pairs from different emitters, the MS method does not depend on this assumption. In samples in which emitters are organized into structures or randomly clustered, the pairs arising from different emitters are also more likely to be at shorter distances, so that the distribution of green points in Fig. 1 b will also be peaked at r shift . However, in our experience, pairs of localizations from the same emitter are more important for the MS and other drift estimates. We also note that our analysis assumes that emitters that are localized in one data set remain within the field of view in the second data set, and vice versa. This may not always be the case, and could in principle lead to bias in shift estimates, but in practice this is typically a negligible effect. Roughly, the contribution of fluorophores near the edge of the field of view may be biased by up to about the localization precision, and the fraction of fluorophores that are affected is restricted to those that lie within about a localization precision from the edge of the field of view, in the direction of the drift. So, for example, in a 100 mm field of view with localization precision of 15 nm, we would expect this bias to be on the order of a picometer.
To benchmark this MS approach, we evaluated the ability of the algorithm to detect known shifts of simulated data sets of a circular test cell, as summarized in Fig. 2. Shifts were estimated by both the MS algorithm and by NLLS fitting of a Gaussian to the spatial cross-correlation function of the two data sets, as implemented in the supporting software provided with (21). The performance of each algorithm was similar for easy cases that produce a welldefined peak at r shift . An extremely easy case is depicted at the top of Fig. 2 a, in which single molecules are well spaced (surface density ¼ 5/mm 2 ) and their positions are well sampled in both frames (twice per molecule on average). In this case, the shift can be clearly identified by eye, and both algorithms reliably and accurately identify the displacement between frames. The simulation depicted at the bottom of Fig. 2 a represents a much harder case, in which molecules are present at higher surface density (20/mm 2 ) and only 1 in 20 molecules are imaged on average in a given data set. In this case, MS modestly outperforms NLLS fitting, both by locating the peak with improved precision and by more reliably finding the peak overall. These trends hold over simulations conducted over a broad range of molecular densities and localizations per molecule (Fig. 2 b). We also estimated shifts from the overall center of mass of each data set, which yielded precisions more than an order of magnitude worse than both the MS and NLLS methods. Moreover, MS is more computationally efficient than NLLS, largely because fast Fourier transforms are not computed in the MS approach. This improvement in speed is enabled through the use of a particularly efficient algorithm from the R package spatstat (30) to extract pairwise displacements between nearby points (see Supporting methods and materials).
For large displacements, both the MS and NLLS algorithms applied in Fig. 2 require an initial step to identify an approximate starting point for the higher accuracy calculation. Fig. S1 shows the failure rate of each algorithm as a function of the distance of the start point from the true shift. MS robustly identifies the main peak over a broad range of simulation conditions as long as it resides within the initial observation window, so large shifts can be identified simply using a large window in the first iteration. This window is typically 100 nm for experimental localizations and 150 nm for the simulations of Fig. 2. NLLS robustly identifies the main peak when the starting point for the computation falls within the localization precision of the peak of the cross-correlation function. In many practical cases, the peak is much farther from the origin than the localization precision, so a separate method is needed to identify a suitable starting point. Here, this is accomplished using a particularly effective algorithm that identifies the global maximum in a smoothed cross-correlation function, as described in the supporting material of (21). The robustness of the NLLS fitting approach is dependent on the ability of this algorithm to identify a suitable starting point over a broad range of simulation conditions. Note that data sets for which the emitter distribution is highly structured or clustered typically lead to improved performance of the start point identification routine by introducing a broad peak in the cross-correlation function in addition to the sharp peak that represents repeat localizations of the same fluorophore.
This MS approach is applied to SMLM localizations that experience continuous drift by distributing localizations into nonoverlapping temporal bins with equal numbers of frames, and displacement estimates are tabulated between all possible pairs of bins. The number of frames in each temporal bin is an important parameter; short temporal bins have few localizations per molecule, so individual displacements may be estimated imprecisely. Long temporal bins have more localizations per molecule and more precise drift estimates but reduce the time resolution of the drift estimate. A linear least-squares fitting algorithm is then used to generate a trajectory that passes through control points positioned, at times, centered on each temporal bin, as described previously (21), taking advantage of the high redundancy to improve precision of the control points. We have slightly modified this past approach by including weights in the linear least-squares fitting, where weights are determined directly from data using a relation that approximates error in the mean displacement (described in Supporting methods and materials), as demonstrated in simulated data sets (Fig. 2 c). Briefly, errors are reduced when there are more pairs originating from the same molecules (magenta points in Fig. 1) and errors increase when more pairs originate from different molecules within the observation window (green points in Fig. 1). Estimated errors can also act as a proxy for overall reliability of the algorithm. Fig. 2 c also shows that MS reliably finds the desired peak when the estimated error remains smaller than one-quarter of the localization precision. This observation can act as a guide when selecting the number of frames included in temporal bins.
It is tempting to distribute frames into overlapping temporal bins, which in principle could improve time resolution while retaining a sufficient number of localizations to accurately determine displacements. However, we find that drift estimates from overlapping time bins are subject to substantial bias, underestimating the actual displacements accrued over time (Fig. S2 a). This occurs because the same localizations are present in adjacent bins, biasing the result toward r shift ¼ 0. Similar bias can arise even in the absence of overlapping time bins because SMLM data frequently contain time-correlated localizations arising from the finite off rates of fluorescent blinking (PALM/dSTORM) or binding (PAINT). These factors mean that pairs of localizations from the same fluorophore are mostly from time separations that are shorter than the time difference between the bin centers and therefore underestimate the average drift between the bins.
The MS approach is applied to a 2D experimental data set of nuclear pore complexes (NPCs) in Fig. 3. NPC assemblies are labeled with primary and secondary antibodies against Nup210 within the nuclear envelope of intact primary mouse neurons, and Fig. 3, a-d shows reconstructed images at various magnifications. Fig. 3 b is a reconstruction produced without drift correction, in which localizations from single NPCs are smeared over a large area, highlighting the importance of drift correction.
The performance of the MS algorithm was tested on this data set by generating multiple drift trajectories through binning with different temporal resolutions. These trajectories were each applied to the full SMLM data set, and Fourier ring correlation (FRC) (31,32) was used to quantify image resolution (Fig. 3 e). For comparison, we conducted drift corrections using the redundant cross-correlation NLLS approach, as described previously (21). In this case, MS modestly outperforms NLLS fitting, allowing for accurate drift correction with smaller temporal bins and modestly improving the resolution of the reconstructed image. We used this data set to explore possible bias introduced because of temporal correlations of single-fluorophore blinking by running the linear least-squares algorithm and including or excluding adjacent pairs of bins on the data of Fig. 2. We found no significant difference between the two cases ( Fig. S2 b), indicating that the impact of this bias is negligible within experimental errors. Additional diagnostics for the MS and NLLS approaches are shown in Fig. S3.
Drift trajectories are shown in Fig. 3 f for temporal bin widths that produce accurate FRC metrics for the MS and NLLS approaches. For MS, a temporal bin slightly larger than the minimum from the FRC curve is used because this produces smaller errors on individual control points. As expected, the drift trajectories follow the same general shape, but the trajectory generated from MS has improved time resolution. In parts of the trajectory, the errors of the control points are smaller than the distance between the trajectories. In these regions, higher time resolution yields improved spatial resolution in the final reconstructed image. Although the differences in the trajectories are significant, their impact is not apparent when viewing reconstructed images of entire nuclei or collections of NPCs, as in Fig. 3, a and c. Differences become more apparent in images of individual pores, in which displacements of several nanometers shift the relative positions of labeled subunits (Fig. 3 d).
The MS algorithm is easily extended to localizations acquired in 3D, in which performance improvements are more evident compared with the established NLLS approach. Because the MS algorithm uses points instead of reconstructed images and fast Fourier transforms, it can be extended into 3D without needing expanded memory resources that limit the practical application of NLLS in 3D. Instead, the 3D application of NLLS drift correction is typically accomplished by generating 2D projections that contain less information than the 3D localizations from which they are produced (21). To see why, consider a pair of emitters that are close together in x-y but far apart in z. Pairs of localizations from this pair of fluorophores will be included in 2D MS drift estimation when using data projected into the x-y plane but excluded from the full 3D drift estimation method by virtue of their large separation in z. We compare the precision and robustness of 3D MS and NLLS on simulated localizations spread over a cylindrical volume in Fig. S4, in which the NLLS correction is performed on projections into the x-y, x-z, and y-z planes, as described in (21). We also directly compare the x-y performance of the full 3D MS method with the 2D MS method performed on data projected into the x-y plane (Fig. S5). Fig. 4 applies the MS approach to an experimental SMLM data set of labeled B cell receptors on the ventral membrane of B cells imaged using a phase mask in the emission path to localize fluorophores in 3D (8). As was the case for simulated data sets, the differences in the performance of the MS and NLLS fitting methods are more pronounced than in the 2D data set of Fig. 3. Additional diagnostics for the 3D case are shown in Fig. S6.
In summary, a mathematically simple MS algorithm modestly outperforms cross-correlation-based estimates of drift correction in 2D and more significantly improves the time resolution of drift corrections in 3D. The approach is computationally efficient, is robust without sophisticated methods to estimate starting points, and does not require image reconstruction with memory and pixelation limitations. The metric provided to estimate error and predict robustness directly from data provides users with a means to evaluate the quality of a drift correction within an SMLM analysis pipeline. For the example data sets explored, modest improvements in resolution lead to adjustments of localized molecule positions relevant for evaluating the structure of protein complexes in cells.

AUTHOR CONTRIBUTIONS
The method was devised by T.R.S. in consultation with F.J.F. and S.L.V. F.J.F. wrote the majority of the code and performed the majority of analyses in consultation with T.R.S. and S.L.V. F.J.F., T.R.S., and S.L.V. wrote the text. S.K. prepared and imaged nuclear pore complex samples for Fig. 3. R.A.B. prepared and imaged B cell samples for Fig. 4.

DECLARATION OF INTERESTS
The authors declare no competing interests.