VivoFollow 2: Distortion-Free Multiphoton Intravital Imaging

Intravital multiphoton microscopy has become one of the central tools used in the investigation of dynamic cellular activity and function in living animals under nearly physiological conditions and is particularly important for studying the dynamic immune system. During intravital imaging in mice, periodic motion of tissue caused by respiration, induces significant shifts of the imaged region. In slow laser scanning imaging modalities, such as multiphoton microscopy, this movement can lead to considerable distortion and discontinuity in three dimensions of the acquired images. Here, we introduce VivoFollow 2, a toolkit that concurrent with image acquisition performs a precise measurement of the respective image distortion, enabling subsequent automated correction of the imaging data. Recovery of one three-dimensional image stack, corresponding to the tomographic tissue sectioning by the optical plane from each single raw image stack, preserves the time continuity within each image stack. Implementation of VivoFollow 2 thus enables a minimization of motion artifacts in tissues exposed to periodic movements and allows for long-term time-lapse imaging and subsequent precise image analysis of the dynamics of cellular and humoral factors in vivo.

Intravital multiphoton microscopy has become one of the central tools used in the investigation of dynamic cellular activity and function in living animals under nearly physiological conditions and is particularly important for studying the dynamic immune system. During intravital imaging in mice, periodic motion of tissue caused by respiration, induces significant shifts of the imaged region. In slow laser scanning imaging modalities, such as multiphoton microscopy, this movement can lead to considerable distortion and discontinuity in three dimensions of the acquired images. Here, we introduce VivoFollow 2, a toolkit that concurrent with image acquisition performs a precise measurement of the respective image distortion, enabling subsequent automated correction of the imaging data. Recovery of one three-dimensional image stack, corresponding to the tomographic tissue sectioning by the optical plane from each single raw image stack, preserves the time continuity within each image stack. Implementation of VivoFollow 2 thus enables a minimization of motion artifacts in tissues exposed to periodic movements and allows for long-term time-lapse imaging and subsequent precise image analysis of the dynamics of cellular and humoral factors in vivo.

INTRODUCTION
In the last two decades, intravital multiphoton microscopy has become one of the central tools in basic and preclinical biological research. It allows for the investigation of dynamic cellular activity and function in living animals under nearly physiological conditions (reviewed in [1] and [2]). Due to the high tissue penetration depth, low phototoxicity, and the possibility of imaging cellular behavior over an extended period of time in living animals, 2-photon intravital microscopy (2P-IVM) has become a sophisticated technique in the field of neurosciences and immunology [3]. For instance, 2P-IVM has been used to visualize a wide variety of biological processes, such as the homing and interaction of dendritic cells with T cells in lymph nodes [4,5], the migration of T cells across the blood-brain barrier into the central nervous system (CNS) [6][7][8], the processes of tumor angiogenesis and metastasis (reviewed in [9] and [10]), and the post-arrest behavior of immune cells on endothelial cells lining blood vessels during their extravasation. A common hurdle is to get optical access to the region of interest, for example, when aiming to image immune cell migration in the CNS, which is shielded by bone, the dura mater, and the leptomeningeal layers [11]. Nevertheless, different surgical preparations have been developed allowing for imaging of the rodent brain and spinal cord, such as cranial windows, thinned skull preparations, and lumbar and thoracic spinal cord windows [12][13][14][15][16][17]. Recently, we have developed a novel window preparation that allows visualization of immune cell migration across microvessels in the cervical spinal cord of mice using 2P-IVM [18].
A major challenge of 2P-IVM is to minimize motion artifacts in the region of interest that are present during the imaging. In the conventional intravital imaging setting, drift of the tissue and movement caused by the vital functions of the experimental animal can compromise the precise analysis of, for example, motility parameters of individual cells, or even make the acquired data completely unanalyzable. In recent years, several approaches have been developed to overcome this problem. When tissue motion is minor, offline processing tools such as StackReg [19] can be used following image acquisition to correct for the shift caused by the tissue drift. Yet, when the tissue shift over the course of time-lapse acquisition becomes comparable or greater than the size of the imaged volume, these tools are not applicable.
Previously, we have developed the VivoFollow feedback microscopy system that performs tracking of the sample drift in real time and adjusts the imaging position directly during the acquisition [20]. VivoFollow has thus solved the issue of severe tissue drift. While originally developed for 2P-IVM, this system can be applied to a wide range of bioimaging settings.
However, the periodic tissue motion caused by respiration (usual frequency 2-4 Hz) poses a very different problem. This issue is especially present while imaging the cervical spinal cord region in a living animal due to its anatomical position. In this case, the shift of the region of interest becomes significant already within a single image stack of a three-dimensional (3-D) tomographic time-lapse imaging stack for the fast imaging modalities (when image stack acquisition speed is >20 frames/s). In the commonly used 2P-IVM employing the galvanometric scanning mirrors, the frame rate is 0.5-2 frames per second, which corresponds to 1-8 periods of tissue motion within a single image frame. Effectively, this causes a significant distortion in each of the acquired images of a 3-D stack. Moreover, uncoherent distortion of each Z frame leads to the discontinuity of the acquired 3-D volume.
For some studies, tissue motion can be reduced by physical restraints by employing coverslips, suction holders, or other stabilization holders [21], but their application presents other significant issues. They can cause obstruction of the physiological blood flow in the tissue or induce tissue damage and inflammation and are thus a trade-off between minimizing the impact on physiological tissue condition and stabilization efficacy. On the other hand, a number of motion stabilization approaches have been developed from the image analysis side. Several Image-J [22] plugins have been introduced for the correction of motion artifacts in time-lapse videos, for the separation of cellular movement from slow tissue movement [23], or for removing distorted and out-of-focus images caused by tissue motion [24]. In another study [25], the authors proposed a method of recovering an undistorted image by composing fragments of the original dataset. Another technique employing time gating of the image acquisition was proposed [26] that enables acquiring image data only during the time when the sample is stationary. Unfortunately, these tools are not easily applicable to 3-D time-lapse imaging of dynamic processes, such as cell migration in vivo.
As it will be shown, a simple synchronization of image acquisition with the respiration of the animal generally solves the issue of 3-D discontinuity, yet the acquired image stacks remain significantly distorted, compromising precise measurements of, for example, cell motility parameters. Moreover, in this case, one cannot directly estimate the value of distortion and corresponding systematic error of the position measurements.
Here, we present a system named VivoFollow 2, which allows precise distortion measurement and subsequent distortion correction by employing a real-time processing pipeline and a novel approach for image acquisition. This approach relies on the synchronization of the image acquisition with the periodic movement of the tissue of anesthetized mice. For the particular examples reported in the present study, the movement was caused by the mouse respiration, and the synchronization was achieved by monitoring the periodic motion of the inhalation anesthesia ventilation unit used for the endotracheal intubation. VivoFollow 2 enables imaging at a close-to-nominal frame rate, requires minimal changes to the 2P-IVM setup, and recovers the images corresponding to the tomographic tissue sectioning by the optical plane, enabling high-precision 3-D time-lapse imaging.

VivoFollow 2 System
When a static sample is imaged, the acquired images correspond to the tomographic tissue sectioning by the optical plane ( Figure 1A). But if the sample is subjected to periodic movement with a frequency comparable to or higher than image acquisition frame rate, each image does not correspond to a plane, but instead to a complicated sheet sourced from a range of depths in the sample. Furthermore, when acquisition of sequential images in an image stack is not synchronized, these surfaces interleave, leading to the discontinuity of the acquired 3-D volume, as shown in Figure 1B.
The VivoFollow 2 system recovers the 3-D image stack corresponding to the optical plane tomographic tissue sectioning in three steps. First, image acquisition is synchronized with the respiration of the mouse, ensuring recovery of the 3-D continuity in the acquired volume ( Figure 1C). To this end, we have developed a specialized Triggering VivoFollow device, named "TrigViFo, " which enables synchronization of individual image frame acquisition with respiration by performing a readout from the piston position of the mouse ventilator, that provides anesthetic gas mixture to the animal, employing a slotted optical switch (see Hardware Synchronization for details). This device allows for a remote automated control of the phase of the respiration at which the image acquisition is triggered.
Second, to recover the optical sectioning of the tissue, we perform a short calibration measurement of the periodic tissue movement at the beginning of the image acquisition. This calibration is done by acquiring 32 "calibration" image stacks at eight different phases uniformly distributed over the respiration period, employing the TrigViFo that can trigger acquisition at any phase of the periodical movement. The tissue motion profile is then built based on the anatomical landmark position in the sequence of the calibration image stacks. This motion profile and the corresponding image distortion map are computed in real time during the image acquisition using these calibration stacks. Finally, by correcting the image stacks for distortion, the tomographic tissue sectioning by the optical plane is reconstructed ( Figure 1D; see Distortion Mapping and Correction).
After the calibration measurement is completed, all of the subsequently acquired image stacks are corrected for the distortion in real time. Building upon the original VivoFollow [20], the present system employs the distortion corrected images to additionally perform tissue drift correction during the course of image acquisition.
Following the experiment, the distortion correction is iteratively refined, and the final corrected image stack sequence is generated. This reduces the displacement of the tissue on the images from typical value of 15 µm to 0.4-0.8 µm, which is comparable with the error of measuring the coordinates of individual cells. The accuracy achieved by VivoFollow 2 allows also the first 32 calibration image stacks to be used for further biological analysis, for example, analysis of the immune cell motility.

Hardware Synchronization
To perform triggering at different phases of the respiration cycle, a dedicated TrigViFo device, allowing for synchronization signal readout, period measurement, triggering pulse generation at adjustable phase, and remote control, was designed. Our setup relies on the use of the Raspberry Pi 3 B+.
Schematic representation of the circuit is given in Figure 2A.
The IR LED-PIN photodiode pair is used as a slotted optical sensor ( Figure 2B) to read the position of the piston of the MiniVent (Harvard Apparatus) ventilator that supplies the isoflurane anesthesia mixed with the oxygen to the anesthetized mouse. The sensor is attached to the transparent tube of the ventilator pump and fixed with a plastic screw ( Figure 2C). The analog signal from the photodiode is amplified and then digitized using the Adafruit ADS115 Analog-to-Digital Converter (ADC). Pulses on one of the microcontroller's GPIO pins are generated to trigger the image acquisition on the microscope ( Figure 2D).
The custom-designed software running on the TrigViFo performs a continuous readout of the signal from the ADC, so that signal sampling corresponds to the ADC frequency, 860 Hz. This signal is then filtered for readout duplicates and smoothened by averaging with a sliding window of size 5. Afterward, the rising edges are detected at the level of a manually selected threshold. An example of the filtered waveform is shown in Supplementary Figure 1. We estimate the time resolution of the device as σ t = σ V dt dV re ≈ 0.2 ms, where σ V is the signal voltage standard deviation at the linear part of the rising edge and dV re dt is the voltage slope at the rising edge. The breathing period τ is measured using a sliding window of the last n = 10 periods, corresponding to the period measurement resolution of σ τ = σ t 2 n ≈ 90 µs. For usual acquisition conditions, according to the motion profiles reconstructed with our system, the maximal tissue motion speed due to animal respiration can be estimated as 200 µm/s. Consequently, the period measurement resolution corresponds to an error of 0.02 µm in tissue position estimation, which is below the resolution of 2-photon microscopy imaging modality.
The rising edge detection threshold can be manually adjusted to the approximate level of the steepest part of the waveform using the graphical user interface (GUI) on the touch screen. The GUI allows the operator to see the current period and the waveform profile to facilitate sensor installation.
The triggering pulse is then generated at time t t = t re + f τ + t where t re is the time point of the rising edge, the factor f defines the phase of the breathing cycle at which the triggering occurs, and delay t = 20 ms is introduced to ensure stable operation. In practice, it can be reduced for applications where the motion period τ is significantly smaller than in our case.
TrigViFo can be controlled remotely using the Remote Procedure Call (RPC) interface operating over the TCP network layer. This allows setting the value of the phase factor f as well as switching the triggering on and off. Acquisition of each image frame is triggered by the pulse from this device and thus is synchronized to the respiration cycle. Diagram of the algorithm running on the device is shown in Supplementary Figure 2. It is implemented in C++ and is running on a standard Raspbian build.

Image Acquisition
Image acquisition was performed using a TrimScope-II 2-Photon microscope (2PM) system equipped with an Olympus BX51WI fluorescence microscope (LaVision BioTec) using a Nikon CFI Apo LWD 25× water immersion objective. Images were acquired using excitation wavelength of 920 and 1,045 nm. The laser beam scanner was running at 800 Hz, and acquisition was performed in a bidirectional mode.
The computational setup for the real-time correction consists of two PC and the TrigViFo device interconnected via 1 Gbit local network. One PC runs the commercial 2PM acquisition program ImSpector (LaVision BioTec), which controls the microscope hardware and image acquisition. This program allows for access to the acquired images, and customization of the acquisition pipeline by a Python script. We have designed a script defining an image acquisition pipeline that controls the triggering phase and the scanning position. The acquisition consists of two sequential measurements. First, we acquired a time-lapse image stack sequence for mapping the distortion and then continued the normal image acquisition with the synchronization as well as with the tissue drift correction. A flowchart of the pipeline is shown in Supplementary Figure 3.

Distortion Mapping and Correction
To perform the distortion mapping, we acquired a sequence of the calibration 3-D image stacks at n ph different phases along the respiration cycle. This is achieved by setting the triggering phase factor f for each of the image stacks sequentially to The total number of acquired image stacks is n st = n ph · n r , where n r is the number of repetitions.
In our experiments, we set n ph = 8 and n r = 4. These numbers were chosen empirically to balance robust operation and keep the time of the calibration stacks acquisition short. Acquired image stacks I t are then filtered for the out-of-range noise, smoothed, and downscaled by a factor of 2 in XY plane. Next, they are split into overlapping image crops I c, t of width 32 pixels with step 8 pixels along the slow scanning axis ( Figure 2E).
For each crop c, we then find the relative position of all n st image stacks. This is performed in two steps. We denote the offset in 3-D between crops I c, i and I c, j at time points i, j as r c, i→j . At first, we find relative offsets between crops at all time points pairwise ( Figure 2F) using fine pattern matching in 3-D as previously described [20]. For this, we rely on the image channel containing the anatomical landmarks. The measurement error is defined by the standard deviation of the peak center measurement on the Pearson correlation map σ p; c, i→j and the maximum correlation value C peak as where C min = 0.6 is the minimum accepted correlation value and ǫ = 10 −4 . This approach allows us to keep track of the alignment error when the distortion causes the appearance of the two image crops to differ significantly.
Next, we perform global solving to find relative positionsr c, i of the crops c at all time points i > 0 with respect to the position of the corresponding crop at the first time pointr c, 0 . This is performed by the weighted least square minimization of L = i,j; i > j The tissue can be subjected to the slow drift during the calibration stacks acquisition often caused by subtle changes in tissue tone. The offset during each repetition k ∈ {0, . . . , n r − 1} due to the slow drift can be estimated as r k = E c (r c, (k+1)·nph − r c, k·n ph ), where E c denotes weighted average with weight 1/σ c,ij , σ 2 c,ij = σ 2 c,i + σ 2 c,j to account for measurement error of each crop's absolute position. Such weighting is mandatory since errors σ c,i can significantly vary when the anatomical landmarks cover the imaging region non-uniformly, which is often the case in practice. The corresponding drift speed is then v k = r k n ph · t , where t is the time step size between consecutive time frames. We then correct the crop positions to account for the linear motion component occurring due to the slow tissue drift: r c,k·n ph +i → r c,k·n ph +i − v k · t · i, ∀k ∈ [0, n r ) , ∀i ∈ [0, n ph ). Afterward, we shift all crop positions to make average positions across all time points equal to zero: r c,i → r c,i − j r c,j n st , so that we will obtain r c,i oscillating around zero ( Figure 2G).
Finally, all crops c are pulled together to obtain the sample motion profile (Figure 2H). For this, we fit the motion profile with a periodic third-order spline with n s nodes uniformly distributed along the periodr m (φ| r s ) =r m y, k ph |r s , where φ = 2π k ph n ph + 2π y λ y is the motion phase, k ph -phase index of the image stack, y is the coordinate on the image along the slow scanning axis (axis Y in our experiments), λ y is the effective "wavelength" of the motion along the slow scanning axis, and r s is the set of all the spline nodes positionsr i , i ∈ {0, . . . , n s − 1}. Please note that while the wavelength λ y is in fact defined by the breathing period τ and line scanning speed v scan as λ y = τ v scan , at this point, it is more convenient to treat it as an independent variable. To find the wavelength λ y , we first fit the Z motion profile with a simple sine function: z y, k ph = R sin(2π k ph n ph + 2π y λ y + φ 0 ). Distortion correction is then performed by resampling the image taken at phase k ph by shifting each pixel at pixel coordinatesr pixel to the new coordinatesr pixel −r m (y p , k ph ). This image resampling is performed employing trilinear interpolation.
The most resource-intensive operation is finding the relative offsets between pairs of image crops. To enable performing distortion correction in real time, in addition to relying on the GPU processing, we have developed a parallel processing scheme. We perform pairwise offset evaluation before all image stacks are available, namely, the offsets r c, i→j , i ≤ j and corresponding errors are calculated right after the time frame j is acquired and during the acquisition of the next time frame j + 1. This allowed us to perform the global solving and fitting of the motion profile right after the acquisition of the last calibration image stack I n st −1 , so that the parameters for distortion correction are ready before the acquisition of the first stack of normal data taking I n st is completed.
Finding the distortion map and performing the distortion correction in real time allow us to apply the tissue drift correction procedure during the acquisition on the images corrected for the distortion.
Finally, we perform additional post-processing after the image acquisition to achieve image quality suitable for further biological data analysis. For this, we iteratively repeat the distortion mapping and correction to find the residual distortion until convergence. During the first iteration, we correct only the Z offset because significant shifts in Z can impair the accuracy of the distortion mapping in XY. Iteration continues until the change in the distortion map remains above 0.1 µm but not more than five iterations to prevent potential divergence in poor landmark signal or landmark coverage conditions. This iterative procedure improves the distortion measurement by 10% (Supplementary Figure 4) and enables performing the distortion correction with subpixel resolution.
The image processing server and the post-processing tool are implemented in C++ and can be run on the Windows or Linux platforms with CUDA-enabled GPU.

Evaluation on Synthetic Data
To assess performance of the distortion measurement and correction procedure, we have first evaluated it on a synthetic dataset. For this, we have generated three image stacks of resolution 512 × 512 × 22 pixels with randomly distributed Gaussian spots (width 10 pixels in XY, 2 pixels in Z) of density 240, 30, and 10 spots/field of view (FoV), corresponding to good, average, and poor landmark coverage in the sample ( Figure 3A). We then distorted them by applying image resampling according to a simulated motion profile r(φ) = A i = 1..4 sin(i φ) i 2 , where φ is the phase of periodic motion. In case of laser scanning imaging, the phase depends on the position on the image: φ = 2π k ph n ph + 2π y λ y , where k ph is the phase index of the image stack, n ph is the number of used phases, y is the coordinate on the image along the slow scanning Y axis, and λ y is the effective wavelength of the motion along the slow scanning axis (see Distortion Mapping and Correction for details). We have chosen an amplitude along all of the three axes of A = 4 pixels that approximately corresponds to the usual experimental conditions. We generated 32 image stacks for each density sample corresponding to four repetitions of eight image stacks at different phases (n ph = 8) and added random Poisson pixel noise to the images. The periodic motion wavelength λ y was chosen to be 343.3 pixels, so that the image height is 1.5 λ y .
The synthetic dataset was then processed, and the image distortion was iteratively reconstructed as described in the Methods section. In Figure 3B, the 8-node spline fit of the reconstructed periodic motion profiles for the three datasets and the ground truth (GT) used to distort the images are shown. Comparison between different landmark coverage for spline fit with 16, 8, and 4 nodes is shown in Supplementary Figure 5.
One can see that the developed method can reconstruct the distortion pattern with accuracy up to 0.1 pixel for a wide range of experimental conditions. In Figure 3C

Evaluation on the 2P-IVM Data
We demonstrate the performance of the system in 2P-IVM imaging of the cervical spinal cord vasculature of mice [18]. To this end, we used CX3CR1-GFP//CCR2-RFP reporter mice systemically injected with a non-function-blocking anti-PECAM-1-Alexa Fluor 555 antibody. This allowed to visualize tissue resident myeloid cells via the GFP signal and infiltrating monocytes by the RFP signal. Systemic injection of the anti-PECAM-1 antibody readily labeled circulating immune cells.
In Figure 4A, we show the MIPs of an image stack of size 395µm × 395µm × 84µm acquired without image acquisition synchronization. One can clearly see the duplicates of the fine GFP + cellular processes of the CNS myeloid cells in XY projection caused by the motion of the sample, as shown in Figure 4C. Due to significant motion in Z direction, the 3-D structural information is lost in this case, as seen on the XZ and YZ projections in Figures 4A,C. In Figures 4B,D, we show the MIPs of the image stack in the same region acquired with triggering by TrigViFo, synchronizing the start of each Z frame acquisition with the respiration cycle. One can see that the consecutive Z planes are aligned, and the 3-D structure is recovered. Since the period of the motion (τ ≈ 0.54s) is, however, comparable with the acquisition time of a single Z frame (0.69 s), the acquired volume is effectively distorted. To measure the distortion pattern, we performed acquisition of several Z stacks, with triggering the start of each Z frame acquisition at different phases of the respiration cycle, as described in the Methods section. In Figure 4E, the examples of the YZ MIPs acquired with triggering phases 0 8 τ , 2 8 τ , 4 8 τ , and 6 8 τ are presented. The difference in the position of the tissue on the respective images is clearly visible. The peak-to-peak amplitude amounts in this case to about 15 µm. We performed the first iteration of the distortion measurement using the myeloid cell CX3CR1-GFP signal as landmark in real time during the acquisition. The corresponding YZ MIPs after the final iterative offline distortion correction of the same dataset are presented in Figure 4F. The tissue position after the distortion correction is the same, confirming that the correction quality is sufficient to include the original calibration image stacks, used for distortion measurement, into the biological analysis. Hence, the first 16 min of acquisition are not lost and remain accessible for subsequent image analysis. The XY, XZ, and YZ corrected MIPs are shown in Figure 5A. We compared the distortion measurement performance using 4-, 8-, and 16-node spline fit in Figure 5B. The four-node fit only slightly (<1 µm) affects the precision of the distortion measurement. Performing the distortion correction in real time enables subsequent accurate real-time tracking of the region of interest in the tissue. This is done by performing continuous offset measurement using fine pattern matching of the acquired image stacks in 3-D and adjusting the stage position where the next time frame is to be scanned in real time during the course of the imaging [20]. The tissue drift profile for a 60min acquisition following the distortion measurement (16 min) is shown in Figure 5C. valid volume after distortion correction for the two datasets, as well as time-lapse of the corresponding MIPs are presented in Supplementary Videos 3-6.
We then used a dataset with three types of anatomical landmarks (Figure 5D), namely, the endogenous VE-Cadherin-GFP reporter, allowing to see adherens cellular junctions of endothelial cells lining the vascular walls, meningeal collagen layers (second harmonic generation, SHG), and an intravascular dextran tracer (TRITC) to assess the distortion mapping performance. These landmarks have different quality and coverage in the imaging volume. For example, we deliberately used a low concentration of the TRITC tracer to assess the performance on poor landmarks. In Figure 5E, the periodic motion profiles along the X, Y, and Z axes of the tissue motion reconstructed using the VE-Cadherin-GFP, SHG, and TRITC signals are shown. We used the four-node spline fit for this dataset. Using an eight-node fit led to divergence of the fit curve because not enough spatial information is available. This is consistent with the synthetic dataset results (Supplementary Figure 4). Flexibility in the number of the spline fit nodes allows us to apply this method to a wide variety of imaging settings. The motion pattern spline fit allows the method to perform robustly also in cases when anatomical landmark coverage is sparse. For example, the SHG signal covers <50% of the FoV. One can also see that the reconstructed Z motion profile, which has the highest impact on the final motility parameter estimation, is reconstructed well in all cases, even for the low signal-to-noise ratio TRITC signal (Figure 5E, red line). The reconstructed distortion corrected image stacks enabled subsequent tracking of the motility of immune cells in the bloodstream, as visualized by tdTomato + CD8 T cells (Figure 5D and Supplementary Videos 7, 8).

DISCUSSION AND OUTLOOK
In this study, we have developed the VivoFollow 2 system, which allows for the measurement and correction of distortion caused by periodic motion, such as animal breathing activity during relatively slow (i.e., acquisition frame rate is comparable with breathing rate) 2P-IVM image acquisition. This system enables further precise localization of objects of interest over time, which is, for example, a prerequisite for the precise measurement of cell motility parameters, movement of humoral factors, or sites of cellular extravasation from blood vessels. While previously described methods discard distorted frames [24] or recover undistorted images by combining tiles of multiple image stacks using pattern matching [25] or gating schemes [26], our method recovers one image stack corresponding to the optical plane tomographic tissue sectioning from each single raw image stack. Preserving the time continuity within each 3-D image stack allows this method to be used for time-lapse imaging of dynamic processes in vivo, such as cell interactions and cell migration. Adopting this method requires minimal changes to the usual 2P-IVM setup, namely, triggering the single frames acquisition and access to the pump of a suitable animal anesthesia ventilation system to perform readout of its motion. In addition, by relying mainly on image processing, the method does not require a piezo actuator for image acquisition. The image acquisition speed is slightly reduced by the increased delay between image frames due to the necessity of synchronization and an increased imaging depth. However, since the method does not involve a gating scheme, the effect on the overall image acquisition speed remains minimal (less than one respiration cycle per frame).
We have demonstrated the applicability of the developed have successfully used for distortion mapping and correction. Most importantly, the implemented spline fit of the motion pattern allows the system to be applicable even in cases when the anatomical landmark coverage in the region of interest is sparse. Finally, by measuring and applying the distortion correction in real time, we were able to perform the tissue drift correction in regions subjected to significant periodic motion.
While in this work we have focused on a single imaging region, the presented method is directly applicable to any imaging of a periodically moving sample with the motion frequency comparable to the image acquisition frame rate. For example, abdominal organs or even the lung are subjected to periodic motion caused by the animal breathing, and the acquired images exhibit the same motion aberrations. Furthermore, this technique can be applied to faster imaging modalities, such as spinning disc confocal microscopy or 2P-IVM employing the resonant laser scanning, when motion is negligible within a single Z frame but remains significant in a Z stack. In this case, one could perform synchronization on the whole Z stack acquisition instead of synchronizing each single frame. By applying the same procedure, the distortion map would be found as a function of the depth Z. Alternatively, fast imaging modalities can benefit from this method when imaging tissue regions subjected to fast periodic motion, such as the heartbeat. The only difference in this case would be the synchronization to the electrocardiogram of the animal, instead of respiration.
One current limitation of the method is the requirement of inhalation anesthesia with intubation, allowing for external breathing control, and consecutively stable periodic tissue motion. While many experiments employ isoflurane anesthesia with intubation, a significant amount of experiments relies on injection anesthesia and recently also perform imaging of awake mice. Demanding the use of intubation would put additional hurdles on the often already complex experimental setups. We plan to address this issue in our further work by performing synchronization for experimental setups, which do not employ ventilators. Also, as previously mentioned, the employed triggering scheme can affect acquisition speed. The increase of acquisition time depends on the imaging settings, ranging between zero and one respiration period per image frame. Researchers should bear this effect in mind when choosing imaging parameters. Finally, since the excitation laser power is adjusted with depth for each image frame during deep tissue imaging, the fluorescent signal brightness will exhibit some degree of inconsistency due to sample movement, the same way as in usual triggered acquisition. While currently we did not address the compensation of the fluorescent signal intensity, with our approach, this becomes feasible due to knowledge of the absolute position of the tissue at each location. This will open the possibility to perform quantitative analysis of the fluorescent signals and, for example, study the T-cell activation dynamics employing fluorescent reporter mouse lines in a moving tissue.

MATERIALS AND EQUIPMENT
Mice C57BL/6J mice were obtained from Janvier Labs (Genest Saint Isle, France). VE-cadherin-GFP knock-in mice express a C-terminal GFP fusion protein of VE-cadherin in the endogenous VE-cadherin locus [27,28]. VE-cadherin-GFP knock-in mice were used for visualizing and imaging of the endothelial adherens junctions [18]. CX3CR1 GFP/+ ::CCR2 RFP/+ mice were kindly provided by Dr. Israel Charo and Dr. Richard Ransohoff [29] and were generated by crossbreeding CX3CR1 GFP/+ mice (B6.129P-Cx3cr1 tm1Litt /J) with CCR2 RFP/+ mice [B6.129(Cg)-Ccr2 tm2.1Ifc /J] [30]. OT-I tdTomato mice were kindly provided by Dr. Jens Stein (University of Fribourg, Switzerland) and were generated by crossing OT-I mice [C57BL/6-Tg(TcraTcrb)1100Mjb/J] [31] with Ai14 mice [B6.Cg-Gt (ROSA)26Sor tm14(CAG−tdTomato)Hze /J]. Before crossing to OT-I mice, the stop cassette was removed from the Ai14 mice by breeding to Tg(ZP3-cre)93Knw. Prior to the experiment, the Ai14.1 allele was bred to homozygosity. All transgenic mice were backcrossed to a C57BL/6J background for at least eight generations. Mice were housed in individually ventilated cages under specific pathogen-free conditions. Adult female and male animals, 8-10 weeks old were used in the experiments. All animal experiments were performed in accordance with the Swiss legislation on the protection of animals and were approved by the veterinary office of the Canton of Bern (permission BE31/17).

Cervical Spinal Cord Window Preparation
Preparation of the cervical spinal cord window was performed as described previously [18]. Briefly, a carotid catheter pointing toward the aortic arch was placed in the right carotid artery of an anesthetized with fentanyl/midazolam/medetomidine mouse for systemic injection of exogenously labeled cells, fluorochromelabeled plasma markers, or antibodies via the arterial blood supply prior to the imaging. Tracheotomy was done to place the tracheal cannula and connect it via a tube to a mouse ventilator (MiniVent, Harvard Apparatus) providing the humidified gas mixture to the mouse (1% isoflurane in a mixture of 40% O 2 and 60% air, tidal volume of 175 µl, inspiration rate of 120 strokes/min). Laminectomy was then performed to expose the mouse cervical spinal cord at the level of C2-C6. The dura mater was kept intact to protect the spinal cord tissue. The physiological parameters (ECG and temperature) of the anesthetized mouse were monitored during the preparative surgery and imaging.

Fluorescent Labeling of Endogenous Blood Cells
To image the interaction of fluorescently labeled endogenous blood cells with the spinal cord microvessels, a non-functionblocking rat anti-mouse PECAM-1 antibody (mAb 390) [32] conjugated to Alexa Fluor 555 was systemically injected via the carotid artery catheter into the surgically prepared mouse (30 µg/mouse). mAb 390 was kindly provided by Dr. Steven  Albelda (Philadelphia, USA). The systemic injection of the PECAM-1 antibody allows for direct labeling of endogenous blood cells, including platelets, monocytes, neutrophils, and subsets of T cells that are circulating in the cervical spinal cord vasculature.
Adoptive Transfer of CD8 + T Cells Expressing a Fluorescent Protein (tdTomato Reporter) The interaction of adoptively transferred tdTomato CD8 + T cells with the postcapillary venules of the mouse spinal cord was imaged. To provide the in vitro activated CD8 + T cells for the adoptive transfer, the single cell suspension was prepared from the spleen and lymph nodes of T-cell receptor (TCR) transgenic OT-I tdTomato mice. The purification and in vitro activation of CD8 + T cells were performed as previously described [31]. Prior to the imaging, 4-5 × 10 6 in vitro activated OT-I CD8 + T cells were suspended in pre-warmed (37 • C) 0.9% isotonic NaCl and then injected via the carotid catheter into the surgically prepared VEcadherin GFP mice. In addition to visible endothelial adherens junctions of the blood vessels in VE-cadherin GFP mice, 1% TRITC-conjugated dextran (MW = 155,000, Sigma-Aldrich, Switzerland) was applied as a plasma marker via intra-carotid artery injection [18].

Image Processing
Image and video preparation for this publication was performed with ImageJ2 and Imaris 9.2 software packages.

DATA AVAILABILITY STATEMENT
The raw data supporting the conclusions of this manuscript will be made available by the authors, without undue reservation, to any qualified researcher.

ETHICS STATEMENT
The animal study was reviewed and approved by Veterinary office of the Canton of Bern (permission BE31/17).

AUTHOR CONTRIBUTIONS
MV implemented the hardware and software for the described system, helped with experimental design and imaging, and wrote the manuscript. NH performed experiments and contributed to manuscript writing. EK performed the experiments. AA and BE provided guidance, coordinated this study, and wrote the manuscript.