Easy-DHPSF 2.0: open-source software for three-dimensional localization and two-color registration of single molecules with nanoscale accuracy

Automated processing of double-helix (DH) microscope images of uorescent single molecules (SMs) streamlines the protocol required to obtain super-resolved three-dimensional (3D) reconstructions of ultrastructures in biological samples by single-molecule active control microscopy (SMACM). Here, we present a suite of MATLAB subroutines, bundled with an easy-to-use graphical user interface (GUI), that facilitates 3D localization of single emitters (e.g. SMs, uorescent beads, or quantum dots) with typical precisions of tens of nanometers in multi-frame movies acquired using a wide-eld DH epiuorescence microscope. The algorithmic approach is based upon template matching for SM recognition and least ‐ squares tting for 3D position measurement, both of which are computationally expedient and precise. Overlapping images of SMs are ignored, and the precision of least-squares tting is not as high as maximum likelihood-based methods. Version 2.0 of Easy-DHPSF augments the approach of Version 1.0 by incorporating code that facilitates combining localization data from two spectral channels using a locally-weighted quadratic 3D registration function.


Introduction
Three-dimensional localization of single emitters for super-resolution imaging and single-particle tracking continues to be a vibrant eld of study 1 . The double-helix microscope 2 provides 3D information from 2D images using PSF engineering to generate images of single emitters which have two lobes, where the zposition information is contained in the angle between the two lobes. The Easy-DHPSF software provides an easy-to-use set of MATLAB programs for analysis of the "Double-Helix PSF" data generated from this kind of microscope, and with minor changes, can be adapted to other PSFs which use lobe arrangement to encode position 3 .
Overview of version 2.0 update Easy-DHPSF version 2.0 utilizes the same algorithm for identi cation of SMs in a DH microscope image as version 1.0: template matching in each region of interest via phase correlation, then double-Gaussian tting of the DHPSF 4 (here called DH for brevity) images through non-linear least squares minimization ( Figure 1). This algorithm represents a good compromise between tting accuracy and computational complexity. For a more detailed description of this algorithm, the reader should consult the documentation for Easy-DHPSF v1.0 5 as the basic details will not be repeated here. Version 2.0 introduces functions for two-color registration and for de ning the calibration curves from many points, such as in a regularly-spaced nanohole array (NHA) 6 .
Registration of two color channels with a locally-weighted quadratic transform One key aspect in most two-color wide-eld imaging experiments is the procedure used to quantitatively map one channel onto the other to provide one consistent coordinate system. Simple implementations have used a ne transforms or other global mapping approaches. Although global mapping functions have the advantage that they can be estimated from only a few corresponding control point pairs, they may not be exible enough to accurately register high-resolution SM localization data. An accurate mapping function de ned over the entire 3D eld of view is essential for wide-eld high-resolution microscopy.
The mapping function implemented in Easy-DHPSF v2.0 is a 3D locally-weighted quadratic transformation. The mapping function t transforms points from one coordinate system into another 7 , t(x)=x'. For coordinate transformations in 3D, t(x) has 3 components, t x' , t y' , and t z' , that separately generate the three transformed coordinates x', y', and z' (see Figure 15). The component functions t x' are given as the weighted mean of polynomial transformations L i,x' using the local weights W i associated with control point i with coordinates x i (see Figure 16), where the sum is over the full set of N control points in the image domain.
The polynomial transformation functions L i,x' are de ned exactly for each pair of control points: i.e., L i,x' (x i ) = x i ' without error. The x' component of a 3D quadratic transformation is given by the equation below, and y' and z' components are similarly de ned (see Figure 17).
The functions L i,x' are su cient to de ne the map between the control points, i.e. the points x i and x i ', while the weighting coe cients W i are used to interpolate between these points. The weight associated with control point i is given as a 3D Gaussian sphere centered at x i . The width of the Gaussian sphere is de ned by a global width parameter s, and a local width parameter σ i that varies from point to point to accommodate the local density of control points (see Figure 18).
For the automated computation of the complete transformation function, the user needs to specify three parameters: the number of neighboring control points n used to calculate the coe cients of the polynomial transformation by least squares estimation, the global width parameter s, and the generating function used to de ne σ i . These three parameters are systematically determined by the optimization program for each set of control point pairs to produce the lowest overall registration errors as described below in Step 8. The design of the registration algorithm is discussed further in Ref. 7.
Below, we demonstrate the use of open-source Easy-DHPSF v2.0 by analyzing an SM dataset of microtubules immunolabeled with blinking Alexa647 (Life Technologies) and CF568 (Biotium) in BSC-1 cells. This code has been validated in several different imaging conditions (Refs. 6,8,9) and is provided free-of-charge under the 3-clause New BSD License at https://opensource.org/licenses/BSD-3-Clause.
Speci cally, this code is provided "as is" without guarantees of any kind.

Reagents Equipment
Calibration and control points -beads vs. NHA Calibration and control-point scans can be performed with either uorescent beads immobilized on a coverslip 8 or with a metal nanohole array (NHA) lled with a solution containing the uorophore labels of interest 6 . Whereas use of beads on a coverslip allows for faster tting and generation of the transformation function, the NHA provides ner sampling of the eld of view. This allows both the DH calibration and two-color registration function to account for eld-dependent errors arising from aberrations in the imaging system. Moreover, by scanning the NHA laterally, we can sample the eld of view in a ner and more controllable way than with beads randomly dispersed onto a coverslip.
Below, to generate the set of control points used to calculate the registration function, we imaged bright uorescent beads (200 nm Tetraspeck, Invitrogen) dispersed on a coverslip with emission simultaneously visible in both detection channels. We repeatedly scanned the positions of beads throughout a 3D eld of view to assemble a large set of control point pairs to nely calibrate the spatial differences in DH response at many locations in the two color channels. A similar procedure of axial and lateral scanning of a uorophore-lled NHA can also be performed.

Software
The Easy-DHPSF package consists of MATLAB functions and subroutines and requires the image processing, optimization, statistics, and wavelet toolboxes. Easy-DHPSF has been tested on workstations running MATLAB 2017b. There are no known issues with earlier versions of MATLAB, but some built-in functions (e.g. fft2, imread) may have reduced performance in earlier versions.
The rst image stacks are the raw data le(s) of the two color channels with SMs to be localized. If desired, the eld of view may contain one or more ducials to correct for sample drift, ideally located near the edge of the eld of view and not overlapping with the region containing single molecules to be localized. The two-color imaging may be performed on a standard two-color SM microscope with each color channel imaged on separate cameras, or the two color channels may be imaged onto different regions of the same camera.
The second image stack is a calibration scan. This should comprise a series of frames at different zpositions of either beads on a coverslip or a NHA that are detected in both color channels simultaneously. As an example: The calibration scan is taken by translating the piezoelectric z stage in 50 nm steps from -1 micron to +1 micron relative to the focus (with ~20 acquired frames at each step to improve precision of the measurement), ensuring that there are enough frames between each step to allow for relaxation of the stage or oil and the sample. A sequence log le should be generated with the sequence of steps used to generate this calibration scan. An example sequence log le for this calibration step is provided (sequenceLog_calibrationScan.dat).
The third image stack is a control point scan. The control point scan is taken by translating the z stage in known step sizes across the range of the DHPSF, i.e., 200 nm steps from -1 micron to +1 micron, with multiple frames (~20) at each step to improve precision. After each complete axial scan, the stage is then translated laterally and the axial scan is repeated to ensure ne sampling of the eld of view. The user should ensure that the relative brightness of the DHPSFs are similar in both color channels. The user should also ensure that this scan is taken with the same lters as those used for SM data. A sequence log le should be generated with the sequence of steps used to generate this calibration scan. An example sequence log for this control point scan step is provided (sequenceLog_controlPointScan.dat).
The last image stacks are 'dark offsets' for the calibration and the data les. The frames in the dark offset .tif stacks are averaged and subtracted before data processing to account for the offset and dark counts from the EMCCD. Acquire a series of a few hundred frames with the camera shutters closed, and use the same acquisition parameters (e.g. integration time, EM gain, shift speed) as those used for the calibration and data acquisition. If these parameters are the same between the calibration and raw data images, only one dark offset le is necessary.

Procedure
Step 0. Image acquisition and program setup Acquire the raw data, calibration scan, control point scan, and dark offset images described in Equipment. Ensure that you have saved a sequence log le for the calibration and control point scans.
Step 1. Two-color DHPSF calibrations This section will outline the steps to perform two-color calibrations of the DHPSF either with beads or with a NHA.

1.1
Open an instance of the Easy-DHPSF GUI by entering 'easy_dhpsf' in the MATLAB command line.

1.2
Before running any of the modules, set the conversion gain and pixel size of the detector using the elds under 'Setup.' Note that changing these parameters after processing the data does not retroactively change the calculated results ( Figure 2).

1.3
Under Setup > Channel, select a channel 'G' for green or 'R' for red.

1.5
Follow the prompts if you would like to automatically save the partly processed dataset after this module.

1.6
Open the calibration images le, the dark offset image le matching the calibration, and the sequence log le (either .mat or .dat) for the calibration. The program will prompt you for these les via dialog boxes. Select a region of interest (ROI) containing at least one bright uorescent object by adjusting the box and double-clicking inside it. Even though the calibration bead may be isolated, it is helpful to select an entire wide-eld excitation region so that the background uorescence can be estimated accurately. Select uorescent objects by clicking the center of each DHPSF ( Figure 3).

1.8.2
Select an ROI by adjusting the box and double-clicking inside it. Select 4 holes that form a large square, then select ducials to outline the edge of your eld of view ( Figure 4). Make sure to try to click the center of each DHPSF, halfway between the two lobes. Hit enter when nished DHPSFs.

1.9
Each frame of raw data, as well as the accompanying double Gaussian t, is displayed as it is processed.

1.10
After processing, an output of calibration information will be generated in the 'calibration [DateAndTime]' folder. Inspect this for each bead to ensure that the angle vs. z position curve is roughly linear, and that the variation of angle measurements in each step is not unreasonably high (error bars in upper-left plot) (see Figure 5 below). Additional information from the calibration, such as the observed shift in x and y position with axial position or wobble, may be useful when aligning the DH phase mask.

1.11
Select the most reliable calibration bead in the Easy-DHPSF main window using the drop-down menu. The calibration information and templates generated from this bead will be used in all later processing steps.

1.12
Go back to step 1.2 and repeat the calibration process for the other color channel. When done, click the save icon (located on the top left of the window) to save your calibrations. You can open this calibration le for use in the subsequent steps.
Step 2. Single-molecule detection calibration for SMACM data In this module, the templates generated from the calibration subroutine are used to generate a large array of matches to the raw SMACM image data. This step serves to assess what value of phase correlation is typical for a good match to an SM image. After this step, the user will de ne a threshold for each template such that only DHPSF images of reasonable quality will be analyzed with the double-Gaussian tting algorithm.

2.1
Under Setup > Channel, select a channel 'G' for green or 'R' for red.

2.2
Click the 'Run' button in the 'Calibrate SM identi cation' panel.

2.3
Follow the prompts if you would like to automatically save the GUI after this module.

2.4
Set the EM gain used, load in the SM data .tif le(s), select what frames to perform the SM detection calibration on, and select the templates to be used. By default, templates are chosen by the program such that the angle of the line connecting the two lobes is given by {-60°, -30°, 0°, 30°, 60°, 90°}, and this selection should generally be appropriate ( Figure 6).

2.5
Load in the dark offset .tif le with the same parameters as the SM data.

2.6
The program will then prompt the user to choose a background processing algorithm such as median ltering, wavelet ltering, no background subtraction, or select if the data is an NHA, Figure 7.
Median ltering parameters should be lled in if choosing this option.

2.7
Select a region containing the SMs of interest. This exact ROI will be used for the singlemolecule tting module and ideally should not contain uorescent objects that are much (≥10x) brighter than the SMs you wish to analyze. Template matches are indicated as circles drawn over the raw data. Stronger matches are drawn as larger circles, and each color corresponds to a different template match.

2.9
Select appropriate thresholds for each template, and enter these into the 'threshold' eld in the GUI. These thresholds should be chosen such that the two lobes of the DHPSF are faintly visible and that there is a low rate of false matches for higher correlation values. You should not be concerned if there are a few incorrect matches, say to one lobe of the DHPSF from a very bright molecule, because these will be rejected by the subsequent double-Gaussian t module.
Step 3. Fiducial tracking (optional) This module tracks the movement of one or more uorescent beads or other stationary markers. Stage drift during a SM experiment can be thus be removed by subtracting the movement of the ducial marker from the SM localizations.

3.1
Click the 'Run' button in the 'Track duciaries' panel.

3.2
Follow the prompts if you would like to automatically save the GUI after this module.

3.3
Select a region of interest containing at least one bright ducial by adjusting the box and double-clicking inside it.

3.4
Select one or more ducials by clicking in the center of each DHPSF.
Step 4. Single-molecule localization for SMACM data Using the results of the previous processing steps, this module identi es single molecules that score above the template matching threshold identi ed in step 2.9, then ts them using nonlinear least squares minimization to a double-Gaussian function.

4.1
Click the 'Run' button in the 'Localize DHPSF SMs' panel.

4.2
Follow the prompts if you would like to automatically save the GUI after this module.

4.3
Choose the frames from the raw data to process, as well as whether you would like to estimate the laser pro le and nd true lobe widths (sigmas). By default, the entire image stack is processed.

4.4
Whenever starting analysis of a new experiment, monitor this process as it proceeds to ensure that the threshold used includes the single molecules visible in the raw data, and to ensure that the ts are performed successfully. Only matches that score above the thresholds are t, and good ts are plotted in the frame-by-frame reconstruction. Here the de nition of a good double-Gaussian t is a t that meets several straightforward criteria, hard-coded into the module. The tests performed to de ne a good t are listed below under Troubleshooting (Step 4).

4.5
Save this Easy-DHPSF le using the save icon at the top of the GUI, ensuring that the Channel selected is the correct one.

4.6
Go back to Step 2, this time selecting the other color channel.

4.7
Ensure that you have two saved easy-dhpsf les, one for each color.
Step 5. Single-molecule detection calibration for control points 5.1 Load in the Easy-DHPSF GUI le from Step 1, saved after calibrating both the green and red channels.

5.2
Under Setup > Channel, select a channel 'G' for green or 'R' for red.

5.4
Follow the prompts if you would like to automatically save the GUI after this module.

5.5
Set the EM Gain used, load in the control point scan .tif le(s), select what frames to perform the control point detection calibration on, and select the templates to be used. A subset of the frames may be used, but ensure that the frames you selected are ones that include at least 1-2 full axial scans so that different z positions are included in the calibration.

5.6
Load in the sequence log le (.dat or .mat) with the sequence of steps used for the control point scan.

5.7
Load in dark offset .tif with the same parameters as the control point scan.

5.8
The program will then prompt the user to choose background median ltering, wavelet ltering, no background subtraction, or select if the data is a NHA array. Median ltering parameters should be lled in if choosing this option. If an NHA is used, select NHA so that the program will t localizations found in a grid -this will increase the number of PSFs detected. If NHA is selected, median background subtraction is performed.

5.9
A frame of the scan will pop up, and the user will be prompted to select an ROI of the correct color channel using a square ROI. The user will then be able to select an ROI with a polygon.

5.10
After the calibration, follow Steps 2.4-2.5 described above for selecting thresholds.
Step 6. Single-molecule localization for control points This will outline the steps to t control points in the two channels.
6.1 Click the 'Run' button in the 'Localize DHPSF SMs' panel.

6.2
Choose the frames you would like to t, as well as whether you would like to estimate the laser pro le and nd true lobe sigmas.

6.3
As this step is for tting the control points, enter how many stationary frames there are for each axial position and click 'OK.' 6.4 If using an NHA, select 4 ducials that form a box that outlines the eld of view (Figure 8). The program will nd the rest of the DHPSFs located on a grid for each axial scan. Otherwise, this step does not show up for bead tting (Figure 9).

6.5
Once the program is nished tting, lter and export your results. Under 'Output DHPSF SM localization', click 'Filter Output.'

6.6
If using a NHA, choose 'Yes' when the program asks whether you would like to use spatially dependent calibrations. Doing this will reduce eld-dependent aberrations by applying nearest-neighbor calibrations to your data with the NHA.

6.7
A different calibration may be used during this step. The program will prompt the user with this option.

6.8
Filter your control points as prompted by the program (for instance, lter by frame, z position, lobe distance, photons, localization precision, etc.).

6.9
Choose the region to be plotted in the scatter plot, inspect the localizations, then save your data when you are happy with the results. The ltered localizations will be in the le 'Output.mat', for use in the next step.

6.10
Repeat steps 5 and 6 for the other color channel. Take note of where you saved the ltered output les, as they will be needed for the next step.
Step 7. Identi cation of control points This module will identify candidate control points between two channels. Open the 'Output.mat' data le with the molecule ts in transmitted/RED channel.

7.4
Open the .mat or .dat sequence log le for the control point scan.

7.5
If available, open a previous transform to use as an initial guess. Otherwise, click 'Cancel' 7.6 Input parameters of how many stationary frames were taken at each z position, how many frames to wait before starting the averaging, and if using an old transform, how closely the initial transform should match old transform points. The program will prompt you for these values.

7.7
The program will rst nd control point (CP) candidates, merging localizations at given z positions.

7.8
A gure will pop up with tting statistics for the molecules ts you loaded. Enter ltering parameters for σ x , σ y , σ z , and number of measurements.

7.9
Inspect averaged bead positions in both channels.

7.10
If an initial transform was not given in step 3.5, the program will ask you to hand click control points. Aim to hand click at least ~15 control points. Follow the prompts in the gure (Figure 10, Figure   11). Both the color and size of each localization encodes the relative brightness (number of photons detected) of that molecule. If the program is still prompting you to select control points but you have already selected >15 and would like to stop, get a new frame then hit 'esc' to end.

7.11
When nished hand clicking control points, the program will run and generate a preliminary transformation le that you will use for the next step. This le will be located in the same directory where your green channel 'Output.mat' is saved.
Step 8. Generating and evaluating nal 3D transform from control points This module will evaluate the initial transform calculated in Step 7 above and generate a nal transform you can use in the next step.

8.3
Enter the range of parameters you would like to evaluate. The program will evaluate these parameter ranges then return the calculated target registration errors (TRE) and ducial registration errors (FRE) 9 . An example result is in Figure 12. Inspect the results by clicking on the gure. Ideally, the user should aim for errors ≤ 10 nm.

8.4
Select the value of nearest neighbors, sigma, and kth nearest neighbor, and input these parameters into the dialog box that appears. The program will then calculate a transformation using the parameters you provided. For the ~1000 control point pairs in the example here, the mapping function took ~1 hour to generate using an Intel Xeon E5-1620 processor.

8.5
The program will return a histogram of TREs and FREs, as well as the spatial distribution of these metrics. The nal transformation is '3D_Transform_lwquadratic.mat' and will be located in the same folder as 'Identify_ControlPoints_3D_output.mat.' Step 9. Transforming SMACM data This module will smooth the ducial tracks, perform the 3D two-color registration, correct for sample drift, apply an index mismatch correction, and calculate localization precision. Indicate whether the green data or red data was acquired rst, or if the two-color acquisition was interleaved.

9.3
Load in the green Easy-DHPSF le. The program will perform smoothing via wavelet ltering of the ducial tracks and apply this correction to your data. Follow the prompts that appear.

9.4
Load in the red Easy-DHPSF le, and follow prompts that appear to smooth ducial track data.

9.5
Load in the '3D_Transform_lwquadratic.mat' le generated in the previous step.

9.6
The program will prompt whether or not you would like to crop each dataset before the transformation.

9.7
Select the location and enter a le name for where to save your registered dataset.

9.8
The registered data is stored in the MATLAB structure 'dataSets,' with 'dataSets(1)' containing the green data and 'dataSets(2)' containing the red data. The elds xFid, yFid, and zFid_n contain the ducial-corrected and index-corrected localizations. The rest of the elds are described in Figure 19.

9.9
Select whether you would like to export the localizations into a csv le. If yes, choose the location and le name of each color le.

Anticipated Results
Below is a 2D scatterplot of Alexa647-and CF568-immunolabeled tubulin in a BSC-1 cell as measured by SMACM blinking of the labels (Figure 13). Both colors are labeling the same structure which contains microtubules in a thin region of the cell. Note that this 2D scatterplot representation of 3D data is quick to render for inspection purposes but does not show depth (z) or density, overemphasizing sparse regions of the sample and often appearing noisier than a density plot. To avoid misinterpreting results, we strongly recommend exporting your two-color localizations to a different program for better visualization with a more reasonable reconstruction, such as that shown in Figure 14. Yellow boxes indicate steps that involve both color channels. Dashed outlines indicate that user input is needed during these modules. The nal result is the red/green box, "Registered 2-color 3D data." The Easy-DHPSF GUI before any of the modules are run. The user can switch between colors using the dropdown menu under Setup > Channel. DHPSF calibration with beads. An example of the GUI that prompts user to select DHPSFs during the calibration step if beads are used for calibration. The program will prompt the user to select at least 1 bead.

Figure 4
DHPSF calibration with NHA. An example of the GUI that prompts user to select DHPSF during the calibration step if an NHA is used for calibration. The program will prompt the user to select four ducials that form a square (DHPSFs 1-4 in the image above), then to select ducials that out the edge of the eld of view.  Typical DHPSF templates. Each template is chosen from the aforementioned calibration scan, where the rotation of the double helix between templates is evenly spaced at ~30°.   Example still image of beads on a coverslip during a control point scan. Still image of a eld of view of uorescent beads in the green (left) and red (right) channels during a control point scan.

Figure 10
Example of the GUI that prompts user to select control point pairs when beads are used. There is a re ection between the re ected (green) and transmitted (red) channels in this example due to the optics needed to place both colors on the same camera. The light blue y = x line is added as an additional visual reference.

Figure 12
Example output from evaluating parameter ranges for transformation.
Page 25/28     Table of elds in output