Sample Drift Correction in 3d Fluorescence Photoactivation Localization Microscopy References and Links

The recent development of diffraction-unlimited far-field fluorescence microscopy has overcome the classical resolution limit of ~250 nm of conventional light microscopy by about a factor of ten. The improved resolution, however, reveals not only biological structures at an unprecedented resolution, but is also susceptible to sample drift on a much finer scale than previously relevant. Without correction, sample drift leads to smeared images with decreased resolution, and in the worst case to misinterpretation of the imaged structures. This poses a problem especially for techniques such as Fluorescence Photoactivation Localization Microscopy (FPALM/PALM) or Stochastic Optical Reconstruction Microscopy (STORM), which often require minutes recording time. Here we discuss an approach that corrects for three-dimensional (3D) drift in images of fixed samples without the requirement for fiduciary markers or instrument modifications. Drift is determined by calculating the spatial cross-correlation function between subsets of localized particles imaged at different times. Correction down to ~5 nm precision is achieved despite the fact that different molecules are imaged in each frame. We demonstrate the performance of our drift correction algorithm with different simulated structures and analyze its dependence on particle density and localization precision. By imaging mitochondria with Biplane FPALM we show our algorithm's feasibility in a practical application.A new wave of cellular imaging," Annu.Ultra-high resolution imaging by fluorescence photoactivation localization microscopy," Biophys.Three-dimensional sub-100 nm resolution fluorescence microscopy of thick samples," Nat.Three-dimensional super-resolution imaging by stochastic optical reconstruction microscopy," Science 319, 810-813 (2008).Three-dimensional localization with nanometer accuracy using a detector-limited double-helix point spread function system," Appl. fluorescent super-resolution microscopy resolves 3D cellular ultrastructure," Proc.Two-color nanoscopy of three-dimensional volumes by 4Pi detection of stochastically switched fluorophores," Nat.Imaging biological structures with fluorescence photoactivation localization microscopy," Nat.Dual-color superresolution imaging of genetically expressed probes within individual adhesion complexes," Proc.Whole-cell 3D STORM reveals interactions between cellular structures with nanometer-scale resolution," Nat.Experimental characterization of 3D localization techniques for particle-tracking and super-resolution microscopy," Opt. Express 17, 8264-8277 (2009).


Introduction
In far-field light microscopy, diffraction limits spatial resolution to about half of the detected wavelength, typically ~250 nm.Super-resolution microscopy has overcome the diffraction limit by taking advantage of fluorescent probe characteristics [1].In particular, techniques such as FPALM, PALM and STORM [2][3][4] exploit stochastic switching of photoconvertible probes to achieve resolution values of ~25 nm in the lateral direction.Their 3D variants can additionally improve the axial resolution from ~600-800 nm to 60-75 nm in the axial direction [5][6]7 ], and down to 10 nm when utilizing two opposing objective lenses [8][9][10].
FPALM and related techniques achieve this resolution improvement by stochastically photoswitching probe molecules between fluorescent states that differ either in emission wavelength or amplitude and are thereby easily distinguishable.Imaging parameters are chosen such that only a small, varying subset of molecules can fluoresce at any time.This leads to a sparse distribution of fluorescent spots that represent the active single molecules in the camera image.The molecule positions are determined by fitting model functions to the intensity distributions and a super-resolution image is compiled from the ensemble of determined molecule positions.Typically, every data set is compiled from a few thousand to more than a million localization events.Since the distribution of molecules in each recorded frame needs to be sparse to localize individual molecules reliably, this requires recording of typically 1,000 to 100,000 camera frames over a time frame of, in the best case, 0.5 s [11], more frequently, however, several minutes [12][13][14].
One of the drawbacks to long time measurements is sample or instrument drift caused by temperature changes or mechanical relaxation effects.Drift can easily be in the range of several hundred nanometers over the course of a few minutes.While this is already problematic in conventional imaging, it becomes unacceptable for super-resolution imaging where drift as low as 10 nm can distort the images.Several approaches for drift correction have been realized in the past.Fiduciary markers, typically gold nanoparticles, quantum dots, or fluorescent beads, can be introduced into the sample.Since they do not bleach significantly over the course of the recording, they can be tracked over the whole time course and their trajectories can be used to correct for drift of the structure of interest [3][4].While this approach is highly reliable, it requires introduction of additional markers into the sample and imaging parameters and instrument need to be adapted to record these markers.A too high marker concentration or markers at the wrong locations can interfere with imaging the structures of interest while a too low marker density leads to failure of the drift correction.In a variant of this approach, antibodies labeled by photoactivatable molecules and non-specifically bound to the cover slip can be imaged by STORM [4].While they blink over the course of imaging, their repeated activation allows for tracking similar to separately introduced fiduciary markers.As an alternative to fiduciary markers, the structure itself can be used to measure and compensate for drift [6,15].Since the shape of the structure does not vary during the measurement, cross-correlation techniques are used to determine shifts between subsequent images.Cross-correlation has been used in energy filtered transmission electron microscopy [16][17], atomic force microscopy [18][19], and scanning tunneling microscopy [18] to correct for sample drift, where subsequent images show the same structure, but differ primarily by their spatial position.
Here, we characterize cross-correlation based drift compensation for 3D super-resolution techniques utilizing stochastic switching.The major hurdle is the fact that different time points over the recording show different subsets of molecules.We show that this fact which intuitively contradicts a cross-correlation approach is no limitation if data processing is performed appropriately.We compare, quantitatively, the performance of our algorithm for two different simulated structures mimicking sub-cellular structures such as bulky tubes of mitochondria or thin cytoskeletal elements.The influences of the number of particles and localization precision of the single molecules are discussed.Experimental one and two-color 3D data of mitochondria recorded with our Biplane FPALM setup [5] confirm the feasibility of our approach for practical applications.

Drift Compensation Algorithm
To determine drift in 3D, we have developed an algorithm utilizing the 3D coordinates as well as the recording time points of the localized particles of a super-resolution data set.No assumptions are made about the microscope and particle localization algorithm.The drift determination algorithm is compatible with, for example, Biplane FPALM [5], astigmatism [6], double-helix [7], or 4Pi-detection based [8][9][10] microscopes.
The assumption underlying our algorithm, based on experimental observations, is that drift is a continuous process but not necessarily linear over the course of the measurement.Drift within a single recorded camera frame is assumed to be negligible.We further assume the general case that drift in the different spatial dimensions is not correlated.
Our algorithm first sorts the particles into T time intervals of equal length.T is chosen large enough so that drift within each time interval can be assumed to be linear, but small enough to include a sufficiently high number of particles in each interval (typically of the order of 1,000).For each time interval t (0 t < T), projections of the 3D data in the x, y, and z-direction are performed and particles binned into pixels.Each pixel value of the resulting three two-dimensional (2D) images therefore represents the number of localized molecules in a certain volume defined by the pixel size (usually 10 nm x 10 nm) and a user-defined depth in direction of the projection.These images are then cross-correlated in 2D with the image of the first time interval t=0.The resulting cross-correlation images which are twice as big as the correlated images are then optionally smoothed with a 2D Gaussian and the maxima of the cross-correlation images are identified.The here described implementation determines maxima positions from the location of the brightest pixel which for the here presented practical applications is sufficiently accurate (±5 nm).Alternatively, the maxima positions can be determined with sub-pixel precision by fitting for example a 2D Gaussian peak function to the cross-correlation images.The two coordinates of each maximum describe the overall drift between time interval 0 and t in two directions.From the three projections, all three drift coordinates can be determined: x and y-drifts are extracted from the z-projection, and z-drift by averaging over the values determined from the x and y-projections.x and y-projections are not used to determine y and x-drift, respectively, since in many practical applications these projections do not contain suitable structures that would allow a reliable maximum localization in the x or y-direction.The drift coordinates are then plotted as a function of t.A cubic spline is optionally fit to the resulting curves to reduce noise.Drift within each time interval is determined by linear interpolation between the drift coordinates obtained for the neighboring time intervals (t-1 and t for particles in the first half of the interval, and t and t+1, for particles in the second half).These drift values are subtracted from the particle coordinates which are then stored as the algorithm's output.

Evaluation of Drift-Compensation Algorithm with Simulated Test Structures
To evaluate the performance of our drift-compensation algorithm, we designed two test structures which feature characteristic properties of mitochondria and the cytoskeleton, respectively, structures that are frequently imaged by super-resolution microscopy.Customwritten LabVIEW (National Instruments, Austin, TX) programs were used to simulate the molecule coordinate lists and time points for these test structures.The first structure consists of 80,000 particles distributed randomly in a 200 nm thick solid tubule arranged as a helix with two full turns around the z-axis (1.6 µm diameter and 1.6 µm height; Figs.1a and b).Figs.1c and d show projections of the second structure resembling a conically deformed spider web.80,000 particles are randomly distributed on circles of 300, 600, 900 and 1,200 nm radii and 12 lines emanating from the tip of the 1.2 µm high cone and ending on the largest circle.
To implement the usually limited localization precision, , all particles were randomly shifted from their ideal position according to a normal distribution with user-selectable full-widths at half-maxima (FWHM).Additionally, each particle was assigned a "recording" time point.For the first simulations (Figs.1-3), 25 nm and 70 nm FWHM were used in x-y and z direction, respectively, which comes close to typically achieved values [5][6].
3D-drift was simulated assuming sinusoidal motion with differing frequencies in each dimension and 500 nm peak-to-peak amplitudes as shown in Fig. 1e.The drift was added to the list of particles represented in Figs.1a-d   In the following, the drift compensation algorithm described above was applied to the data sets.The particles are separated by the algorithm into T=40 equal time intervals and plotted as 2D sum projections, shown in Figs.2a and b for the z-projection of time intervals t=0 and 16 of the mitochondria structure, respectively.The cross-correlation between the initial time interval t=0 and all subsequent intervals is calculated as described above and smoothed with a 2D Gaussian (Fig. 2c; optimal smoothing width depends on structure type).To determine the "zero" drift position in the cross-correlation plots, the peak of the autocorrelation function for the t=0 data is used (shown as cross in Fig. 2c).Figs.2d-f show the same data for the xprojection, Figs.2g-l represent the same procedure for the cytoskeleton structure.
The drift of each time interval is determined from the distance between the corresponding cross-correlation peak and the autocorrelation peak.x and y-drift values are determined from the z-projection cross-correlations and z-drift is determined by averaging the z-values for the x and y-projection cross-correlations.Fig. 3a shows the determined drift values as a function of t (dashed line) overlaid onto the graph of the actual drift (solid line; from Fig. 1e).A cubic spline is fit to each curve to smooth and counter any over or under-corrections.Drift within each time interval is corrected by linear interpolation using the neighboring time intervals as described above.
Fig. 3a shows that the determined drift values closely match the actual drift.To further analyze the quality of the drift compensation, we calculated the difference, d, between the drift values obtained from the algorithm and the original drift values that were fed into the simulation.As can be seen from Figs. 3b and e, drift could be corrected to values only rarely exceeding d=5 nm for the mitochondria simulation (Fig. 3b) and 10 nm for the cytoskeleton data set (Fig. 3e).These values are significantly smaller than the localization precision values, , of 25 nm (x, y) and 70 nm (z) assumed for the simulations.Consequently, applying the obtained drift correction values to the simulated data sets results in practically perfect images as demonstrated in Figs.3c, d, f, and g.
The drift-correction algorithm executes in less than 1 s per time interval (less than 1 min total) on a standard PC.Applying the drift-compensation algorithm iteratively on already corrected data did not lead to any additional improvement.It is to be expected that the drift compensation works best for large particle numbers.Figs.3b-g represent the data shown in Fig. 1 consisting of 80,000 particles (2,000 particles per time interval).For comparison, the same analysis has been performed for only 4,000 (100 per time interval; mitochondria test structure) and 8,000 particles (200 per time interval; cytoskeleton simulation) as shown in Figs.3h-j and Figs.3k-m, respectively.The d-values are significantly larger for these much lower particle numbers, but in most cases still smaller than the localization precision.
To investigate the dependence of our drift compensation performance on particle numbers, n, per time interval and also localization precision values, , more systematically, we determined d for five different n-values and five different -values and both simulated structures.From this data, d , the standard deviation of d over the T=40 time intervals, was calculated for each simulation and plotted as a function of and n (Figs. 4 and 5).Despite the different nature of their particle distributions, similar results were achieved for the bulky mitochondria and the fine microtubule structure.For the first, had almost no influence on the drift correction accuracy (Fig. 4a, c).Only for the highest -values ( z =350 nm vs. 200 nm tubule diameter; Fig. 4c) the z-component of d approximately doubled.For the cytoskeleton simulation, a stronger dependence on can be observed (Figs.5a, c).We conclude that localization precision only plays a major role for drift correction when the observed object itself is dominated by fine structures.Interestingly, drift correction values far better than the assumed tubule diameter can be achieved.
Both, the mitochondria structure and the cytoskeletal structure show a strong dependence of d on the number of particles, n, for n<500.For larger n, d is barely affected (Figs.4b and  d and Figs.5b and d) reaching values of ~5 nm and better for all directions and both structures.From a practical standpoint, where <70 nm is usually obtained, n=500 is therefore sufficient, assuming that the imaged structure resembles the shown simulated structures.

Super-resolution Imaging of Mitochondria
To confirm the feasibility of our drift compensation algorithm in practical applications, we imaged mitochondria networks in human hepatocellular carcinoma (HepG2) cells with our Biplane FPALM [5] instrument and applied the drift compensation algorithm to the data.
For one color imaging of mitochondria networks, mtEos2 was activated with our 405 nm laser and fluorescence read out with the 556 nm laser.To ensure an even distribution of detected molecules over the entire axial range, the objective position was periodically stepped at 200 nm intervals over a range of about 8 µm for the selected cell.Photobleaching of the sample occurs at all axial positions regardless of the focal position on the camera, while localization is limited to a 1 to 2 µm range centered around the two focal planes of the camera [20].The axial scanning procedure was therefore repeated several times with only a fraction of the molecules activated during each cycle.This ensures nearly constant particle yields over the complete axial range.The cell displayed in Fig. 6 has a thickness of ~8 µm in the axial direction.Measurement took place over the course of ~1 hour and resulted in >1.05 million localized particles.To correct for the several 100 nm of drift observed during the recording of the described data set, we ran our drift compensation algorithm.The resulting data is shown in Fig. 6.For display, particles were binned in 50 nm pixels.Fig. 6a represents all particles summed up axially.The color scale in Fig. 6b represents the axial positions of the particles in the projected image.To date, this is to our knowledge the thickest image presented from 3D super-resolution localization-based techniques.
Furthermore, we imaged HepG2 cells co-transfected with mtEos2 and PS-CFP2 fused to mitochondrial transcription factor-A (TFAM; pLenti6.3/V5-DESTvector has been adopted for conjugation of any ORF with the PS-CFP2).TFAM localizes to nucleoids of mitochondrial DNA in the mitochondrial matrix.Apart from this co-transfection step, samples were prepared in the same way as described above.
To image both fluorescent protein populations, we followed a protocol as described by Shroff et al. [13] with the following modification: mtEos2 was successfully activated using 488 nm instead of 405 nm laser illumination.This provides the advantages that first, the population of PS-CFP2-TFAM will not be activated and depleted while imaging mtEos2.Second, activated PS-CFP2-TFAM does not need to be bleached before imaging it since the 488 nm will have bleached any previously active PS-CFP2.A bandpass filter (FF01-594/730-25, Semrock, Rochester, NY), located outside the microscope exit port, blocked the fluorescence from mtEos2 in its green emission state and PS-CFP2-TFAM.mtEos2 was imaged until the population was nearly exhausted.Then the bandpass was removed to image PS-CFP2-TFAM by activating at 405 nm and exciting fluorescence at 488 nm.The same multiline dichroic and bandpass filter cube set was used for both colors.This allowed for an easier alignment of the two-color signal after processing of the data.
Since the two color channels were imaged sequentially, drift-correction is mandatory to draw any conclusions about co-localization.A two-color image before and after driftcorrection is shown in Fig. 7 which represents ~110,000 localized mtEos2 molecules (red label) and ~6,000 localized PS-CFP2-TFAM molecules (green label).Overlaying the two colors was performed using the drift correction algorithm outlined above.Both color-channels were first corrected separately to compensate for drift during the recording of each color channel.In principle, overlaying the red mitochondria and green nucleoids should not work this way because the two probes label two completely different structures.Crosstalk from a small population of mtEos2 molecules still present during the first several hundred frames recorded after switching to PS-CFP2-TFAM imaging (see low level signal in Fig. 7c), however, allowed aligning the two channels by cross-correlating the first time interval of the PS-CFP2-TFAM data to the mtEos2 data recorded before.

Conclusions
Sample drift in super-resolution microscopy has a deleterious effect on a microscope's performance as drift can easily exceed its resolution.The algorithm described here allows compensating for drift in all three dimensions down to a sub-5 nm level for localization-based super-resolution methods.We have demonstrated that drift correction can be performed by applying cross-correlations to different, temporally separated, subsets of localized molecules representing the same, fixed, structure.
Certain types of structures are better suited for drift correction than others.Long filaments oriented in the x-direction, for example, look nearly identical no matter how much the sample drifts in the x-direction; y-drift, in contrast, is easily detectable with the same structure.As a rule of thumb, only molecules that are within a distance a from an object "edge" (a structural feature indicating a strong change in density) contribute to detection of drift of magnitude a.This phenomenon has to be considered when choosing the region of interest (ROI) for the projection included in our algorithm.If the volume chosen for the projection is thick and the object in that volume so dense that the projection does not show dominant structural features anymore, the algorithm fails.It is therefore important that the ROI for the drift correction is chosen appropriately from the recorded data set.For the x and y-projections of the data sets represented in Figs. 6 and 7, ROIs ranged over ~25% of the respective depth while the full depth was used for the z-projection.For the presented simulations, on the other hand, the ROI encompassed the complete data sets for all projections.
Next to these structural influences, we have shown a dependence of the algorithm's performance on localization precision, , and on the number of molecules, n, in each time interval (Figs. 4 and 5).The latter has two effects on drift detection precision: with higher n for a given structure, the nearest neighbor distances (NNDs) decrease.Calculated average NNDs for our simulated mitochondria data ranged from 51.3 nm (n=100) to 11.2 (n=2,000).Comparing the average NND values to the determined drift correction precision confirms a linear relationship between the two measures.Statistical arguments further indicate that doubling n while maintaining the average NND and local structural features (for example by imaging two identical structures next to each other) results in 2 -fold drift compensation improvement.
In practical applications, as demonstrated for diverse samples represented by the simulated structures, the described cross-correlation technique can correct drift of several hundred nanometers to values below 5 nm.This is achieved for typical localization precision values and requires only several hundred to 2,000 localized molecules for each time interval.Our method requires no fiduciary markers and is straightforward to implement.We believe that this technique can easily be applied in a wide variety of super-resolution microscopes.

Fig. 1 .
Fig. 1.Simulated test structures and influence of drift.(a, b) Projections of the "mitochondria" structure (pixel size 25 nm).(c, d) Projections of the "cytoskeleton" structure (pixel size 10 nm).(e) Overview of drift values as a function of the time interval t applied as offsets to the particle positions.(f-i) Projections corresponding to (a-d) after application of drift values.Both structures consist of 80,000 particles.
resulting in the smeared structures shown in Figs.1f-i.

Fig. 2 .
Fig. 2. Examples of projections for different time intervals and the corresponding crosscorrelation functions.(a, b) z-projections of the mitochondria simulation for t=0 and t=16.The particles are plotted in 10 nm pixels.(c) Smoothed cross-correlation function of the data displayed in (a) and (b).The cross represents the peak of the autocorrelation function denoting zero drift.(d-f) Data corresponding to (a-c) for the x-projections.(g-l) Equivalent analysis for the cytoskeletal structure.Scale bars are 500 nm.

Fig. 3 .
Fig. 3. Drift correction for different structures and numbers of particles.(a) Simulated drift (solid lines/circles) and example of drift determined by the drift correction algorithm (dashed lines/hollow circles).(b) Difference between simulated and detected drift for the mitochondria simulation at n=2,000.(c, d) z and x-projections of the data set after correction of drift.(e-g) Same representation as shown in (b-d) for the cytoskeletal structure.(h-m) Analysis equivalent to (b-g) for n=100 (mitochondria) and n=200 (cytoskeleton).Scale bars are 500 nm.

Fig. 4 .
Fig. 4. Standard deviation, d , of d as a function of localization precision and the number of particles n per time interval for the mitochondria simulation.(a, b) Values for x and y-drift.(c, d) Values for z-drift.

Fig. 5 .
Fig. 5. Standard deviation, d , of d as a function of localization precision and the number of particles n per time interval for the cytoskeleton simulation.(a, b) Values for x and y-drift.(c, d) Values for z-drift.

Fig. 6 .
Fig. 6.Mitochondria network in a HepG2 cell labeled by mtEos2 and imaged with a Biplane FPALM instrument after drift correction.Approximately 1 million particles were localized over a 14 µm by 26 µm area and an axial range of ~8 µm.(a) z-projection.The color bar represents the number of detected molecules in each 50 nm pixel.(b) Same data with colorcoded axial positions for each molecule.(c) 3D volume rendering using the Vutara SRX software.Scale bars are 1 µm.

Fig. 7 .
Fig. 7. Biplane FPALM image of mitochondria (red; mtEos2) and nucleoids (green; PS-CFP2-TFAM) in a HepG2 cell prior to and after drift correction, in (a) and (b), respectively.Area denoted by the cyan box in (b) is shown in separate channels in (c, d).The color tables denote numbers of molecules in each 50 nm pixel.The color table is saturated to highlight lower densities.Scale bars are 1 µm.