Single-Particle Tracking with Scanning Non-Linear Microscopy

This study describes the adaptation of non-linear microscopy for single-particle tracking (SPT), a method commonly used in biology with single-photon fluorescence. Imaging moving objects with non-linear microscopy raises difficulties due to the scanning process of the acquisitions. The interest of the study is based on the balance between all the experimental parameters (objective, resolution, frame rate) which need to be optimized to record long trajectories with the best accuracy and frame rate. To evaluate the performance of the setup for SPT, several basic estimation methods are used and adapted to the new detection process. The covariance-based estimator (CVE) seems to be the best way to evaluate the diffusion coefficient from trajectories using the specific factors of motion blur and localization error.


Optical setup
This study was carried out with an inverted scanning microscope. The optical optimization of laser beam is illustrated in Figure S1. The general diagram of the setup is shown in Figure S2. Figure S1. The laser beam is deflected by oscillating mirrors (scanner). Two lenses are used to image the scanning device to the back aperture of the microscope objective. The lens after the scanner is the scan lens with a focal length f1, while the other lens is the tube lens with a focal length f2. The laser is a Ti:Sa femtosecond laser (Tsunami, Spectra-Physics) pumped by a 10 W continuous YAG:Nd laser at 532 nm (Millenia, Spectra-Physics). The maximum gain is obtained for = 750 nm (average emission power 1.4 W). The pulse duration is about 80 fs at the laser output, but after passing through various optical elements, the chromatic dispersion widens each pulse whose effective width is estimated at 100 fs at the sample level. The rate of the laser is 80 MHz.
The laser beam first passes an electromechanical shutter with an automated control, (Uniblitz). An optical fiber is used to take part of the beam scattered on the shutter and to direct it to a compact spectrometer (OceanOptics) in order to visualize the emission spectrum on a dedicated computer, via the OOIBase32 software. The spectral width and the emission wavelength are thus controlled.
The laser beam passes successively through two cells: a Faraday isolator (IO-3-780-HP, Thorlabs) and an ultra-fast Pockels cell (350-80LA, Conoptics). The first cell is designed to prevent any reflection from returning to the laser cavity. The electrically controlled Pockels electro-optical cell is used as an intensity modulator. A half-wave plate λ/2 achromatic (700-1000 nm, Edmund Optics) placed on a motorized rotation stage (PR50, Newport) controls the beam polarization. This stage is connected to one of the eight power controls of a universal controller (XPS, Newport).
The study area is scanned using two galvanometric mirrors (6210H, Cambridge Technology Inc.) controlled by a two-axis driver (MicroMax, Cambridge Technology Inc.). These drivers are voltage-controlled via a signal delivered by an analog card at a maximum rate of 1 MHz and coded on a 12-bit dynamic range (NI PCI-6711, National Instrument). The swept angular space at the output of the mirrors, and by extension the swept area on the sample in the focal plane of the microscope objective, is directly related to the amplitude of this voltage. The different apertures encountered in the microscope strongly limit this amplitude, which must be limited to 2V peak-to-peak in order to obtain a uniform distribution of the laser intensity over the entire scanned surface.
The microscope used is a commercial inverted microscope (IX71, Olympus). A collimated green LED source (Model MTN005303, Thorlabs) is used to make images of the sample in brightfield reflection with a high-sensitivity monochrome CMOS camera (UI-3240ML, Stemmer Imaging). The specimen holder is placed on a motorized positioning stage (XYZ 562-ULTRALIGN and 3 TRA12CC motorized stages, Newport Corp.). The laser beam that is directed onto the microscope is a beam previously deflected by the scanner. An afocal optical mount provides the desired beam broadening factor. The chosen solution is to use a scanning lens placed according to the diagram in Figure S2. In order to reduce geometrical aberrations, a telecentric lens with a focal length of 54 mm (LSM54-850, Thorlabs) was finally implanted.
The reflected beam (SHG or TPEF) is directed onto an ultra-fast photomultiplier (H7422P-40, Hamamatsu) with a rise time of 1 ns and a spectral response ranging from 300 nm to 720 nm, the maximum sensitivity being at 580 nm for a quantum efficiency of 40%. The digitized signal (TTL pulses) is then directed to a photon counting module (C9744, Hamamatsu) with a counting frequency of up to 10 MHz. This module is connected to one of the 8 32-bit counters of a high-speed PCI card (PCI-6602, National Instrument). The timing of the counting module is synchronized with the electronic signal controlling the scan mirrors by a specific high-speed RTSI (Real-Time System Integration, National Instrument) bus through a cable that directly connects the two PCI cards to each other.
The control interface was developed in our laboratory with the Labview design platform (LabVIEW 2010, National Instrument). This interface controls the positioning boards, the laser intensity and polarization, the scanning parameters, it generates the scan signals and the images from the non-linear optical signal obtained by the counting module, performs the synchronization between these two events (scanning and detection) and can save the obtained images as well as all the parameters used for the acquisition.
The scanning settings are as follows: -Image size. The program was developed to generate square images. Therefore, choosing 100 dots for a line generates an image of 100 x 100 = 10,000 pixels.
-Zoom factor: This is the maximum voltage applied to the scan mirrors, which determines the physical size of the scanned area of the sample.
-Exposure time: This is the dwell time, that is the integration time during which the counting takes place. This time corresponds to the laser exposure time of each scanned area of the sample. The minimum value of the exposure time, imposed by the mechanical movement of the mirrors, is 1 µs.

Raster scan methods
Several methods have been tested to scan the area of interest. The scanning imaging particle tracking method requires the optimization of the scanning time to obtain high frame rates. For this purpose, the sampling rate of the electronic counter is set to a value close to 1 MHz, which corresponds to a time-laps close to 1 microsecond, including the dwell time and the time for moving the mirrors to the next position. With a traditional raster scan ( Figure S3a), all lines are scanned in the same direction. This requires the "column" mirror to return very quickly to the beginning of the line when the last position on the previous line has been reached. For fairly large scanned areas, this return time may be longer than the sampling period. In this case, the image obtained is not representative of the scanned area because measurements were taken during the return.
With a "snake" scan ( Figure S3b), travelling back to the beginning of the line is no longer an issue because the lines are scanned alternately in one direction and in the other. This method requires to reconstruct the final image by inverting the even lines, but it has the huge advantage of eliminating dead times after each scanned line. However, there is still a dead time at the end of the scan, which allows the mirrors to return to the initial position for the next scan. Although this dead time is usually not an issue, it is directly linked to the frame rate in the method described here. The issue with this direct return to the initial position is that the time needed to complete the movement is not exactly known. Indeed, for small scan areas, this time will render only a few sampling periods, while for larger scan areas, the number of sampling periods taken up by this displacement increases.
The final solution chosen is represented in Figure S3. It is a snake-like raster scan, but the return to the initial pixel is controlled, by scanning pixel by pixel the first column from the bottom to the top. This dead time is therefore controlled, and its value is perfectly known for each acquisition. Still, it remains a dead time, corresponding to the scanning of an additional column. Because this time is perfectly known and controlled, the exact frame rate is then perfectly known for the future size estimations.
The Figure S3d represents a scanning method which was tested and eventually abandoned. Its purpose is to reduce as much as possible the dead time between each image, by making two successive inverted images with a snake-like scan. With this approach, once the area has been scanned the first time, the mirrors take the reverse path to return to the starting point. A sequence provides the data to create two images, with maximal frame rate and no dead time. This method is the most efficient to obtain quickly images of still objects, but it does not work for the present goal because the time to access a given particle is not the same for even and odd images, making further processing for tracking unnecessarily complicated.

Timing of images
In nonlinear microscopy the illuminating beam should scan an area in order to build up the image pixel by pixel. As a result, each image is not a snapshot of the observed scene, because each pixel is acquired at different times. This is the fundamental difference from the usual particle tracking techniques based on video sequence analysis.
The photons emitted by the fluorescent particles are captured by a photomultiplier tube. Consequently, the information about each pixel is acquired at different times corresponding to the dwelling time per pixel (the time spent by the laser to scan each pixel).
Snake-shaped raster scanning means that no time is wasted by alternating the scanning direction of the lines. A square image of N pixels per side is obtained after a time equal to N 2 × t dwell (green area in Figure S4a). The laser moves then back to the first pixel at the end of the image scanning. This return is done with a controlled time equal to N × t dwell . This time, represented in red in Figure S4a, is a dead time that slightly slows down the frame rate of the video finally obtained but this time is controlled and precisely known whatever the size of the image.
A complete sequence acquisition of 4 successive images is shown in Figure S4b. The time between each image is ∆t. During this time, the particle has moved. With traditional particle tracking methods, the sampling time is fixed and the coordinates (X,Y) of the particle in the two successive images change ( Figure S4c). Here, we access the displacement of the particle between images but also to the variable lag time δt between two acquisitions. This variable, which can be positive or negative, is illustrated in Figure  S4d.

Raw images
The scanning method used to build images limits the signal intensity received on each pixel. The raw images contain very few information because the photomultiplier only counts a few photons for a short dwell time (5µs). The images only use 2-4% of the grey scale dynamic range as represented with an enhanced intensity in Figure S5 in which the maximum raw intensity is 5 over 255.
This image shows the very small amount of information contained in an image of moving objects and the presence of noise coming from the excitation signal. The compromise between image resolution and frame rate mentioned in the main article is well represented here. The low resolution of the image makes difficult to accurately estimate the center of each bright spot. The dynamic range of this enhanced (Raw intensity multiplied by 55) image is shown with the histogram in Figure S6.

TrackMate process
An important step in the detection process is the choice of parameters in the TrackMate plugin. Below is described the procedure of detection from the raw images sequence until the generation of coordinates for every track.
1) In order to record several tracks, a sequence of 1000-2000 images is downloaded in the plugin TrackMate. For the first step of the detection, a type of detector must be chosen, in this study only the LoG (Laplacian of a Gaussian) detector is selected and makes sub-pixel localization possible.
2) To proceed to the detection, an expected blob diameter (purple circle) and a threshold are required to give the initial parameters for the LoG filter. An overestimated threshold will cause several gaps in the track when an underestimated threshold will result in false or erroneous positions detection. This parameter can be adjusted using a detection preview feature. 3) Then all the images of the sequence are searched for blob detection and if two spots are too close, the highest intensity is chosen. 4) When the blobs have been detected and localized, all the successive positions must be tested within a maximum link distance (i.e., search radius) represented by a faded orange circle centered on the blob. If the next image contains a similar blob in this area, then it is considered as the same particle. Every trajectory is build this way until the particle disappears. This maximum search radius is an important parameter depending on the velocity of the particle (i.e., the size of the particle). When particle size is known, it is easy to find this parameter using the maximum diffusion length equation (eq 20 on the main paper): 5) The visualization of all the trajectories is important to detect unusual experimental conditions like important drift, aggregation of several particles, convection movement. But all the short trajectories (less than 100 positions) are not reliable for the estimation of the size. 6) In order to optimize the size estimation, only tracks with more than 100 positions are selected.
7) The last step of the process is to generate a text file with all the positions making up each tracks. This text file is then read by a program on Matlab for the size estimation. Figure S8 and Figure S9 are simple description of the algorithm used to process the raw positions data for the method of constant time lag and variable time lag. The final result gives the estimation on the bead's radius. The drift represents a specific direction of displacement not coming from the brownian motion. Figure S6. Histogram of the raw image Figure S5 represented with a Log scale for the counts. The raw intensity image were multiplied by 55 to be able to see all the bright spots and spread the diagram values.

Principle
As originally proposed by Richards and Wolf in 1959 [1], the excitation focal volume is computed by summing the interfering beams emanating from a single plane wave refracted at multiple angles by the microscope objective.
The interference, that is to say the sum, of all plane waves that converge to point x 2 , y 2 , z 2 of the focal plane can be written : where : • Φ(s x , s y ) take into account eventual aberrations or the focalization through a dielectric medium like the glass coverslip [2]. • s, is the propagating vector of one plane wave. This plane waves are ordered around a cone and can be written in spherical coordinates as : s = (sin θ cos ϕ, sin θ sin ϕ, cos θ) • Ω is the solid angle and dΩ = ds x ds y s z sin θdθdϕ This can be sum up as : (3) where : • A(θ, ϕ) is the normalized amplitude of the optical beam before focalization. As a first approximation, it can be considered as uniform. More usually, it can be described by a gaussian : where f is the focal distance of the objective lens and w is the waist of the laser beam. • B(θ, ϕ) is a term that accounts for energy conservation for an aplanatic optical system : B(θ, ϕ) = √ cos θ [3]. • P(θ, ϕ) describe the polarization state inside the focal volume. It can be written as : with : -T(θ, ϕ) is a transformation matrix that permits to go from the plane right at the exit of the objective to the plane of the focal volume. P 0 (θ, ϕ) is the initial polarization state of the optical beam before focalization. The transformation matrix T(θ, ϕ) can be built from rotation matrices as : where :  Further mathematical simplifications could be done assuming a gaussian profile for the laser beam A(θ, ϕ) and introducing Bessel functions. This would lead to faster calculation at the cost of generality [5].
Proceeding with the calculation of the rotation matrices, one can rewrite the polarization state P(θ, ϕ) inside the focal volume as :

Numerical procedure
Numerical calculation is done via discretization of the variations of θ and ϕ : exp ikn(z 2k cos θ m + sin θ m cos ϕ n x 2i + sin θ m sin ϕ n y 2j ) ∆θ∆ϕ (4) Figure S9. Flow chart of the algorithm using variable time lag to estimate the size of particles.

Unoptimized code
The following unoptimzed python code is loosely based on the matlab one that can be found in Qinggele Li's PhD thesis [6].