2.1.

Animal Preparation

Three healthy, young-adult female C57BL/6 mice ($\sim 3$ months old; The Jackson Laboratory, Maine, United States) were used in this study. A sealed 3-mm-diameter cranial window was installed over the whisker barrel cortex area (on the left hemisphere; $\sim 2.0\mathrm{mm}$ posterior to bregma and $\sim 3.0\mathrm{mm}$ lateral from the midline). Tracheotomy was performed to enable controlled ventilation and maintain proper physiological status in subjects by providing a mixture of medical air ($\sim 60\mathrm{m}\mathrm{L}/\min$) and isoflurane (1.2% to 1.5%). A femoral artery catheter was installed to facilitate the administration of a contrast agent in 2PM imaging and to monitor physiological parameters, such as blood pressure, heart rate, and blood gases. Mice were kept under isoflurane anesthesia during the subsequent imaging sessions.

2.2.

DLS-OCT Imaging

DLS-OCT was performed using a spectral domain OCT setup⁴² (Thorlabs Inc., New Jersey, United States) [Fig. 1(a)]. The setup utilizes an extended broadband superluminescent diode (SLD) source (LS2000B, Thorlabs Inc., New Jersey, United States) that implements two matched-pair SLDs to offer $>170\mathrm{nm}$ bandwidth at $\sim 1300\mathrm{nm}$ center wavelength. A $10\times$ objective (MPlanApoNIR, $\mathrm{NA}=0.26$, Mitutoyo, Japan) was used in the experiments. An InGaAs line scan camera with 1024 pixels was used to record the interference signal at 46,000 A-scans/s. The system can achieve $\sim 3.5\mu \mathrm{m}$ spatial resolution in all three dimensions within the brain tissue and $\sim 1\mathrm{mm}$ maximum imaging depth. OCT en face images of the pial surface were obtained for selecting a $600\mu \mathrm{m}\times 600\mu \mathrm{m}$ FOV ($400\text{pixels}\times 400\text{pixels}$). In this study, we only imaged down to $\sim 250\mu \mathrm{m}$ beneath the brain surface of each subject considering the $\sim 150\mu \mathrm{m}$ confocal parameter of the system since it provided a sufficient number of arterioles, venules, and capillaries to demonstrate the proposed framework. We first performed OCT angiography and repeated B-scan acquisition protocol was employed. The final 3D angiogram was obtained with 10 times averaging, which took 88 s for data acquisition. For the blood flow speed measurements using DLS-OCT, an $M$-mode data acquisition protocol with 100 A-scan repeats at each transverse location was employed and the data acquisition time for a $600\mu \mathrm{m}\times 600\mu \mathrm{m}\times 250\mu \mathrm{m}$ (voxel size: $1.5\times 1.5\times 2.9\mu {\mathrm{m}}^3$) volume was $\sim 387\mathrm{s}$. Since the measurement of a component of flow velocity that is parallel with the optical axis (Vz) has the largest SNR of all DLS-OCT velocity estimates, we relied on it to calculate the microvascular flow velocity based on the angle between vessel direction and the optical axis of the system.

Fig. 1

The diagrams of the OCT and 2PM imaging setups. (a) Diagram of the spectral domain OCT system. SLDs, superluminescent diodes; FC, fiber coupler; CLs, collimating lenses; DC, dispersion compensation; RM, reference mirror; M, mirror; GM, galvo mirrors; SL, scan lens; TL, tube lens; DM, dichroic mirror; OL, objective lens; DG, diffraction grating; SOL, spectrometer objective lens; LCCD, line-scan CCD; DAQ, data acquisition card; and PC, computer. (b) Diagram of the 2PM system. M1 and M2: mirrors; EOM, electro-optic modulator; SH, shutter; GM, galvanometer mirrors; SL, scan lens; TL, tube lens; DM, dichroic mirror; OL, objective lens; FFs, fluorescence filters; PMTs, photomultiplier tubes; DAQ, data acquisition card; and PC, computer.

2.3.

2PM Angiography

The 2PM angiograms were acquired after DLS-OCT measurements using a custom-built 2PM setup⁴⁵ [Fig. 1(b)]. A mode-locked laser (Insight DeepSee, $\sim 120\mathrm{fs}$ pulse width, 80 MHz pulse repetition rate, Spectra-Physics, California, Unites States) tunable between 680 and 1300 nm was used for the 2P excitation. The output power of the laser was controlled by an electro-optic modulator (EOM) (Model 350-160, ConOptics Inc., Connecticut, United States). Beam scanning in the $XY$ plane was achieved by a set of galvanometer mirrors (6215H, Cambridge Technology Inc., Massachusetts, United States) with laser beam relayed by a scan lens ($f=30\mathrm{mm}$, AC254-030-B, Thorlabs Inc., New Jersey, United States) and a tube lens ($f=180\mathrm{mm}$, Olympus, Japan) to the back aperture of the objective. A water immersion objective lens (XLUMPLFLN20XW, NA = 1.00, Olympus, Japan) was used to focus the beam on the sample. The far-red fluorescent dye Alexa Fluor 680 conjugated to 70 kDa dextran was used to label blood plasma as the contrast agent for 2PM ($400\mu \mathrm{M}$, 0.1 mL). The emission signal of Alexa Fluor 680 was detected using a photomultiplier tube (PMT) (H10770PA-50, Hamamatsu, Japan), a dichroic mirror (FF875-Di01-$25\times 36$, Semrock, New York, United States), and two fluorescence filters (FF01-890/SP-50 and FF01-709/167-25, Semrock, New York, United States). In this study, 2PM survey images of the pial surface were first acquired and, with help of previously acquired OCT en face images, used to define an ROI ($689\mu \mathrm{m}\times 689\mu \mathrm{m}$) for 2PM that contained the ROI of DLS-OCT. The 2PM angiogram was subsequently acquired down to an $\sim 800\mu \mathrm{m}$ imaging depth (voxel size: $1.35\times 1.35\times 2.0\mu {\mathrm{m}}^3$). The use of animals in this study was approved by the Institutional Animal Care and Use Committee at Massachusetts General Hospital. Typical examples of the volumetric datasets acquired from the two imaging sessions are visualized in Figs. 2(a)–2(c).

Fig. 2

Top views of OCT and 2PM volumetric datasets acquired in one mouse. (a) Maximum intensity projection (MIP) of the OCT angiogram stack (FOV: $600\mu \mathrm{m}\times 600\mu \mathrm{m}$). (b) Maximum projection of the DLS-OCT axial blood flow velocity Vz stack (FOV: $600\mu \mathrm{m}\times 600\mu \mathrm{m}$). (c) MIP of the 2PM angiogram stack (FOV: $689\mu \mathrm{m}\times 689\mu \mathrm{m}$). Scale bars: $200\mu \mathrm{m}$.

2.4.

Angiogram Segmentation and Graphing

For quantitative modeling and analysis of cerebral microvascular networks, graph-based representations of the structure are needed.^46,47 The 2PM angiograms were analyzed using Amira software (Thermo Fisher Scientific, Massachusetts, United States) combined with custom software written in C++, running on a graphical processor unit equipped computer, in the following steps. (i) Adaptive contrast enhancement. For each pixel in each image, the distribution of intensities in a surrounding $201\times 201$ window were sampled to create a histogram and the intensities at the 50th and 99th percentiles ($I_{50}$ and $I_{99}$) were calculated. A linear mapping of intensities was performed, with $I_{50}$ mapped to zero and $I_{99}$ mapped to maximum intensity. This technique counteracts the effect of systematic intensity variations through the depth of the stack. (ii) Vesselness filtering. A novel 3D vesselness filter was developed. A set of 28 directions in space was defined, giving approximately uniform coverage of a hemisphere. For each direction, a test function of position was defined as $f=\exp [-1/2(x_i/{\sigma }_x{)}^2-1/2{(y_i/{\sigma }_y)}^2]$, where $x_i$ is the distance from the origin parallel to direction $i$, and $y_i$ is the distance perpendicular to direction $i$, and ${\sigma }_x=3$ and ${\sigma }_y=0.5\text{pixels}$. For a given voxel, this function was convoluted with the image stack for each of the 28 directions, and the maximum result of these was used to construct the processed image at that voxel. This test function has an elongated prolate ellipsoidal distribution, and the procedure preferentially enhances narrow filamentous structures. Because of the small size of the test function, larger structures in the image are not significantly affected by this filter. Larger vessels generally show continuous intensity in the axial direction and do not need enhancement by vesselness filtering. (iii) Fill-in filtering. Some larger vessels show low intensity in their interior, with intensely labeled walls. To fill in such structures, the following method is used. From a given voxel, a set of six lines of 25 voxels extending in each (positive and negative) coordinate direction is defined, and the maximal intensity on each line is identified. If the mean of the six maxima exceeds the intensity at the given voxel, the intensity of the given voxel is increased, by an amount that decreases exponentially with the coefficient of variation of the six maxima. This method ensures that a voxel that is surrounded in multiple directions by voxels of higher intensity receives a maximal boost in intensity. (iv) Segmentation. Using Amira, the stack was thresholded to obtain a 3D solid representing the network. (v) Skeletonization. Using Amira, the solid was skeletonized, to yield a graph consisting of nodes and edges with defined diameters. (vi) Refinement. The skeleton was further processed to combine short, connected edges and to remove tight loops, which occur as artifacts of the skeletonization. The result is a representation of the network as a set of nodes and edges, with several edges forming a segment between vessel branching points. Segment diameters were estimated as the median of edge diameters in the segment. In this study, the graphs that represent the 2PM angiograms of three animal subjects were able to enclose in a connected network more than 99.9%, 99.3%, and 99.5% of the total vessel segments, respectively. A version of the custom software used is available online ( https://github.com/secomb/StackEnhanceV1).

The vessel type (e.g., arteriole, venule, or capillary) was associated with each segment by a combination of manual and automated vessel labeling.^48,49 First, all pial arterioles and venules, as well as the initial segments of diving arterioles and ascending venules close to the pial surface were identified manually based on their morphology and flow direction in their branches that dive in and surface out of the cortex (e.g., penetrating arterioles and surfacing venules, respectively). Then the labeling of the remaining network was performed automatically starting from the labeled pial vascular segments and propagating the same label (e.g., arteriole or venule) down the vascular tree, excluding the microvascular segments that branch to a transverse plane. The remaining microvascular segments were labeled as capillaries. Finally, a color-coded mask of vessel types was generated based on the automatic labeling and overlaid on the angiogram stack for manual inspection and correction [Fig. 3(a)].

Fig. 3

Angiogram segmentation and graphing results. (a) Overlay of vessel types on angiogram. Overlay of vessel branching order on angiogram from the (b) arteriolar side and (c) venular side.

The capillary branching order starting from precapillary arterioles was computed by first assigning branching order zero to all pial and diving arterioles and subsequently assigning branching orders to the rest of the graph network starting from zero-order segments [Fig. 3(b)].⁴⁹ The same procedure was then repeated for assigning the capillary branching order starting from the postcapillary venules [Fig. 3(c)].

2.5.

Coregistration of the 2PM and OCT Angiograms

Angiogram coregistration was performed by transforming the volumetric intensity images obtained by 2PM to the DLS-OCT space. To account for the non-linearity of image spaces caused by distortions, such as field curvature across those two imaging modalities due to the implementation of different optics, an initial global 3D coregistration was combined with multiple regional 2D coregistrations. The global 3D coregistration was implemented over the whole 3D 2PM and OCT stacks with an affine transformation matrix computed based on $>20$ manually selected fiducial marker voxels representing the same locations in two stacks. Since 2PM and OCT image distortions change differently along the axial direction, they are difficult to address using global coregistration only. Therefore, coarse coregistration was followed by dividing the stacks along the depth dimension ($Z$) into a set of $40\text{-}\mu \mathrm{m}$-thick sub-stacks and conducting 2D coregistration of each pair of sub-stacks. The transformation matrices for 2D local coregistrations were computed from additional fiducial marker voxels manually selected in the sub-stacks after global coregistration. The results of angiogram coregistration for two example sub-stacks selected at different depths are shown in Figs. 4(a)–4(c) and Figs. 4(e)–4(g). Clear improvements are observed from applying combined two-step coregistration compared to using global 3D coregistration only [see Figs. 4(d) and 4(h)]. Herein, we applied a quick and intuitive coregistration method that can better account for the non-linear optical distortions of those two imaging modalities than the methods mainly based on affine coregistration. We also conducted a brief comparison of our method with other coregistration methods, such as the ones based on FreeSurfer,⁵⁰ NiftyReg,^51,52 and ANTs^53,54 packages. Our method has shown improved coregistration results compared to the MRI_Robust_Register Function from the FreeSurfer package, which is mainly based on affine coregistration and comparable results to the state-of-the-art methods based on NiftyReg and ANTs, which both utilize non-linear deformable coregistration in addition to the affine coregistration (see the Supplementary Material for detailed results of the comparisons). In addition, our method enables a considerable degree of manual, user control over the coregistration process, which may be helpful to deal with sometimes significant deficiencies in acquired datasets, such as large shadows below the vessels or missing data in large sub-ROIs.

Fig. 4

Angiogram coregistration results. Transformed 2PM pial stack (a) (depth range: 41 to $82\mu \mathrm{m}$) in DLS-OCT space and corresponding coregistration results with DLS-OCT angiogram layers (gray) overlaid by the 2PM angiogram (green) using (b) global 3D coregistration only and (c) combined global 3D coregistration and 2D regional coregistration. (d) Zoomed views of two sub-areas in stack (a) that highlights the improvements with combined 3D and 2D coregistration. Transformed 2PM in-depth stack (e) (depth range: 157 to $197\mu \mathrm{m}$) in DLS-OCT space and corresponding coregistration results with (f) global 3D coregistration only and (g) combined global 3D coregistration and 2D regional coregistration. (h) Zoomed views of two sub-areas in stack (e) that highlights the improvements with combined 3D and 2D coregistration. White arrows highlight the improvements by the regional 2D coregistration. All scale bars are $200\mu \mathrm{m}$ in length.

2.6.

Blood Flow Velocity Mapping

Based on transformation matrices, all nodes from the graph-based network obtained from the 2PM angiogram were co-registered with the DLS-OCT volumes at different depths [Figs. 5(a) and 5(c)]. This enabled mapping of the blood flow velocities obtained by the DLS-OCT to the microvascular segments obtained by the 2PM angiography [Figs. 5(b) and 5(d)]. To assign the blood flow velocities to the microvascular segments represented by the vascular graph, a small cylindrical volume was defined in DLS-OCT space along each graph edge according to the diameter of the corresponding microvascular segment. A mean value of Vz component of the blood flow speed was subsequently extracted from each cylinder. A custom-written software was written in MATLAB (MathWorks, Massachusetts, United States) to enable manual adjustments of cylinder position, orientation, and diameter to correct for the remaining imperfections of coregistration and better match with the DLS-OCT signal from the vessel. In the case of larger vessels ($>8\mu \mathrm{m}$ diameter) that exhibit parabolic-like flow profile, we extracted the Vz component of the maximum axial velocity in the vessel. The maximum axial velocity in larger vessels typically corresponds to the central region of the vessel. Therefore, the region from which velocities were extracted in such vessels was constrained to an $\sim 8\text{-}\mu \mathrm{m}$-diameter cylinder around the center of the vessel. The absolute RBC flow velocities in such vessels were computed using the extracted flow velocities from the cylinder and cosine of the axial tilt angle between the corresponding graph edges and the optical $Z$ axis and defined as mean center-line flow velocities. In the case of capillaries and small arterioles and venules, the extraction region was adjusted to encapsulate most of the vessel cross section since the difference between the RBC velocities from the vessel axis and periphery, if any, could not be resolved in our measurements. Absolute RBC velocities along such vessels were computed similarly using the above method, and the mean center-line flow velocities and mean flow velocities are equivalent for such vessels.

Fig. 5

2PM vascular network graphing and coregistration with the DLS-OCT blood flow velocity measurements. The transformed graph of (a) pial and (c) in-depth 2PM stacks superimposed with top MIPs from the corresponding angiogram layers. Co-registered vascular graph and corresponding DLS-OCT axial flow velocity measurements for (b) pial and (d) in-depth stacks. All scale bars are $200\mu \mathrm{m}$ in length.

DLS-OCT measurements of the Vz component of the RBC velocity exhibit the largest and the smallest SNR when RBCs are travelling parallel and perpendicular to the optical axis ($Z$), respectively. Consequently, in the capillary network, Vz velocities were typically too noisy when vessel center-line formed a large axial tilt angle (approximately larger than 50 deg, see the Supplementary Material for details) and led to a pronounced decrease in the goodness of DLS-OCT flow velocity fitting. Therefore, we chose 35 deg as a conservative threshold for the axial tilt angle of capillary edges in data analysis to maintain a high-accuracy standard in DLS-OCT measurements while still enabling the measurements of velocities in a large number of capillary segments. Only the vascular edges with strong DLS-OCT signals were included in the data processing. However, vascular segments typically contain tens of graph edges and mean segment velocities could be obtained even if some edges within these segments needed to be removed from analysis. Examples of mapping the blood flow velocities onto a graph-based network are shown in Fig. 6.

Fig. 6

Visualization of axial blood flow velocity (Vz) mapping. The top and side projections of a graph-based vascular segment and its corresponding flow velocity distribution from (a) an artery, (b) a capillary, and (c) a vein. Green lines indicate edges along the segment with a white line segment highlighting the current selected edge. The white shadow areas show projections of the blood flow velocity extraction region after manual corrections. All scale bars are $20\mu \mathrm{m}$ in length.

Although most data analysis steps have been performed automatically by the custom-written software, a certain amount of manual corrections and adjustments are still needed in the graphing and coregistration procedures for refined blood flow speed mapping, which is time consuming. In the future, machine learning-based methods may be developed to assist with the manual-correction steps.