3.1.

Numerical Simulations

We investigated the improvement obtained by adding DOCT measurements to the recovery procedure using a synthetic vessel network. The synthetic network contained 175 nodes and 240 vessel segments, including 1 inflowing arteriole and 3 outflowing venules, as shown in Fig. 2(a). The capillary nodes connected three capillaries, which is the typical configuration measured experimentally.²⁰ The diameter of the arteriole and the venules was set to $20\mu \mathrm{m}$, whereas the diameter of the capillaries was fixed to $5\mu \mathrm{m}$. The length of all vessels was set to $50\mu \mathrm{m}$. Pressure boundary condition was set to 105 mmHg in the inflowing arteriole and to 10 mmHg in the outflowing venules. The pressure at the end node capillaries was set between 30 and 80 mmHg depending on the position in the network. Flow was first computed given this set of boundary conditions. The noise was then introduced by randomly varying the values for the pressure boundary conditions and the resistance of the vessels. Noise level was set to 25% (max/mean) for pressure boundary conditions which is of the order of the variability measured experimentally.¹⁸ Noise in the arterial and venous resistance was set to 25%, which corresponds to an experimental error in the diameter estimation of 5% to 6% (15 to $20\mu \mathrm{m}$ $\pm 1\mu \mathrm{m}$). For the capillary resistances, the error in the experimental diameter estimation is also $1\mu \mathrm{m}$ but since the mean diameter of capillaries is about $5\mu \mathrm{m}$, this corresponds to an error of 20% which transforms into an error of 100% in the resistance calculation. We, therefore, used 100% for the noise of the capillary resistance in our simulations. Flow was then recovered using the erroneous pressure boundary conditions and resistance values without DOCT priors using Eq. (11). A fraction of the vessels were then randomly selected and DOCT priors were set for these vessels. Flow was then reconstructed using these priors from Eq. (16). The fraction of vessels with DOCT priors ranged from 2 to 10%. The value for the DOCT measurements was set to the true (simulated) flow value with noise added. Different levels of noise for the DOCT values were tested ranging from 5% to 25%. The performance was quantified by computing the mean squared error (MSE) between the simulated and the recovered flow in both cases (i.e., with and without DOCT priors). To insure reproducible results, the simulations were run 500 times with different noise instances for boundary conditions and resistance values. MSE values were averaged across all 500 simulations.

Fig. 2

Simulations showing the improvement in accuracy obtained by integrating DOCT measurements in the computation of flow. The fraction of vessels with DOCT was 10% and noise in the DOCT measurement was 10%. (a) Segments containing DOCT priors. (b) Error map in percent for flow recovered with and without priors. (c) Comparison of the MSE obtained for different configurations simulated. Each bar represents the average MSE obtained for 500 simulations with different DOCT prior configurations and different noise instances. Simulations were run with 10% DOCT priors located only in the capillaries, with 10% DOCT priors located anywhere but including one, two and three outflowing venules, and with 10% DOCT priors located anywhere but including the inflowing arteriole.

3.2.

Animal Preparation

Animal protocols were approved by the Massachusetts General Hospital Subcommittee on Research Animal Care. We anesthetized a C57BL/6 mouse (male, 25 g) with isoflurane (1-2% in a mixture of ${\mathrm{O}}_2$ and air) under constant temperature (37°C). We then opened a 2.5-mm-wide cranial window over the primary somatosensory cortex with the center of the window in the parietal bone 2 mm posterior from the bregma and 2 mm lateral from midline, removed the dura, and sealed the window with a $150\text{-}\mu \mathrm{m}\text{-}\text{thick}$ microscope coverslip. We used a catheter in the femoral artery to administer the dyes, to continuously monitor the blood pressure and heart rate, and to sample systemic blood gases (${\mathrm{PCO}}_2$ and ${\mathrm{PO}}_2$) and pH. During the measurement period, the isoflurane anesthesia was reduced to 0.7% to 1.2% and the following physiological parameters were maintained: arterial ${\mathrm{PO}}_2=130\pm 20\mathrm{mm}\mathrm{Hg}$, arterial ${\mathrm{PCO}}_2=37\pm 3\mathrm{mm}\mathrm{Hg}$, and mean arterial $\text{blood pressure}=80–100\mathrm{mm}\mathrm{Hg}$.

3.3.

Multimodal Microscope

Our multimodal microscope, described in Sakadžić et al.,²¹ has independent TPLSM and OCT scanning arms that share the same imaging objective. This allows for easy switching between sequential TPLSM and OCT measurements.

The TPLSM optical beam was scanned in the $x\text{-}y$ plane by galvanometer scanners (6215H, Cambridge Technology, Inc.). An electro-optic modulator (ConOptics, Inc.; extinction ratio 500) served to gate the output of a Ti:Sapphire oscillator (840 nm, 80 MHz, 110 fs, Mai-Tai, Spectra-Physics). The pulse duration was 350 fs (assuming a ${\mathrm{sech}}^2$ pulse shape), as measured at the sample. The emission was reflected by a dichroic mirror (LP 735 nm, Semrock) and detected with a large collection efficiency ($\sim 15{\mathrm{mm}}^2\mathrm{sr}$) by a detector array, consisting of four independent photomultiplier tubes.

DOCT is a path length-resolving optical imaging technique that can measure velocity axial projections in well-defined microscopic volumes. These axial projections can be used to perform absolute measurements of blood flow. DOCT has been applied to brain microvasculature¹⁵ and validated against hydrogen clearance methods.¹⁶ DOCT can achieve a spatial resolution of $7\mu \mathrm{m}$ in the transverse plane and 3 to $4\mu \mathrm{m}$ in the axial direction. The penetration in cortical tissue can reach beyond 1 mm. The temporal resolution for CBF measurement in a single vessel is on the order of 0.1 s and the time required to image an entire volume of $1{\mathrm{mm}}^3$ is $\sim 10\mathrm{s}$. A superluminescent diode (model EX8705-2411, Exalos Inc.) with a center wavelength of 856 nm, bandwidth of 54 nm, and a power of 5 mW was used as the light source. The axial (depth) resolution was $\sim 6.4\mu \mathrm{m}$ in air ($4.8\mu \mathrm{m}$ in tissue) after spectral shaping. The power on the sample was 1.5 mW, enabling a maximum sensitivity of 99 dB. A home-built spectrometer using a 2048 pixel, 12-bit line scan camera (Aviiva SM2 camera, e2v Semiconductors) enabled an imaging speed of 22,000 axial scans per second and an axial imaging range of 2.5 mm in air. The transverse resolution was $11\mu \mathrm{m}$.

Two pairs of galvanometer mirrors (Cambridge Technology, Inc.) were used to independently scan the TPLSM and OCT beams. The galvanometer mirrors were relay imaged onto the back focal plane of a microscope objective (Olympus XLUMPLFL20XW/IR, 0.95 NA, water immersion, 2 mm working distance). The TPLSM and OCT optical beams were coupled into the microscope objective by a movable dichroic mirror. A motorized stage controlled the focal position by moving the objective along the vertical axis ($Z$). Both TPLSM and OCT systems were controlled by custom-designed software written in LabView (National Instruments).

3.4.

Experimental Measurements

We first performed a DOCT scan protocol for quantitative flow measurements. The details of the scanning protocol can be found in Srinivasan et al.¹⁵ The DOCT scan generated a 3-D map of the axial ($z$) projection of velocity. The $X-Y$ resolution of our OCT was 10 to $15\mu \mathrm{m}$, preventing velocity estimates in capillaries because of partial volume effects. In addition, the system could only measure velocity projections in the range $\pm 3.6\mathrm{mm}/\mathrm{s}$ which can potentially create phase wrapping in large arteries. Because of these two limitations, DOCT measurements were only considered in larger ascending venules in this work. When the $X\text{-}Y$ resolution of the DOCT scan is smaller than the diameter of the vessels, the flow at specific locations in a vessel is obtained by summing the velocity axial projection over the vessel area in the en face (also known as transverse or $xy$) plane.¹⁵ In this work, an alternative method was used to avoid partial volume effects. The flow in a given vessel was computed from the maximal velocity measured in this vessel and the diameter of the vessel measured from the TPLSM angiogram, assuming a parabolic flow profile. The maximal velocity was obtained from the DOCT velocity axial projection by correcting for the angle between the longitudinal direction of the vessel and the axial direction of the OCT beam. Since ascending venules are almost parallel to the OCT beam, these angular corrections were very small.

Blood plasma was then labeled with fluorescein isothiocyanate-conjugated dextran (FITC) at 500 nM concentration in blood. We then acquired a $600\times 600\times 660\mu {\mathrm{m}}^3$ stack of the FITC-labeled vasculature using TPLSM ($1.17\times 1.17\times 2.0\mu {\mathrm{m}}^3$ voxel size) as described in Sakadžić et al.²¹ Although OCT can also be used to generate angiograms, the spatial resolution and contrast-to-noise ratio is generally higher with FITC TPLSM, which makes it easier to graph (i.e., vectorize) the vascular stack in the subsequent steps. This data was previously published in Sakadžić et al.²¹ and Gagnon et al.²²

3.5.

Data Processing

The TPLSM angiogram was graphed with a custom software using a combination of automatic segmentation procedures and manual corrections.^5,23 Vessel diameter was estimated at each graph node by thresholding the image at a low value of $\sim 2\%$ of the maximum image intensity, considering lines through the node point oriented every 3 deg in the local plane perpendicular to the vessel axis, and finding the minimum distance from vessel edge to vessel edge. To reconstruct flow in the experimental vessel network, we set the blood pressure boundary condition for the arterioles and venules crossing the sides of the imaged volume using tabulated values from Ref. 18 given the vessel diameter and the vessel type (arteriole or venule). Similarly to our simulations, pressure boundary condition was set to 105 mmHg in the inflowing arteriole and to 10 mmHg in the outflowing venules. This set of boundary conditions has been shown to result in accurate and realistic flow distributions in previous studies.^5,12

4.1.

Simulations Results

The results of our simulations are shown in Fig. 2. In the case of this specific output, 10% of the vessel segments contained DOCT priors. A random noise was added to each of those DOCT priors and the maximum amplitude of the noise was set to 10% of the flow value for each segment. Figure 2(a) shows the specific segments that contained a DOCT prior for this specific simulation. Figure 2(b) shows the $\mathrm{\mu }\mathrm{CBF}$ error map in percent with and withtout DOCT priors. We first observe that the residual error in both cases is not located around particular vessels. It is rather spread across the entire network. We also observe a reduction in the error of the reconstructed flow using DOCT priors. The error reduction is spread across the network and is not limited around vessels containing a DOCT prior or restricted to inflowing or outflowing vessels.

To account for the variability introduced by the noise, simulations were run 500 times with different noise instances and different configurations for the DOCT priors. The average residual error in the reconstructed flow was 300% without DOCT measurements and 150% when 10% of the segments contained DOCT measurements, even if those measurements contained 25% noise. The high level of residual error in the reconstructed flow is the result of a high noise in the resistance of the capillaries which was set to 100%.

Figure 2(c) shows the comparison of the MSE between the simulated $\mathrm{\mu }\mathrm{CBF}$ and the recovered $\mathrm{\mu }\mathrm{CBF}$ with and without the DOCT priors. Each bar represents the average MSE obtained from 500 simulations with different noise instances and different DOCT prior configurations. Simulations were run with 10% DOCT priors located only in the capillaries, with 10% DOCT priors located in any vessels including 1, 2 and 3 outflowing venules, respectively, and with 10% DOCT priors located in any vessels including the inflowing arteriole. Note that vessels with DOCT priors did not contribute to the MSE calculation and thus did not bias the results. We observe that constraining the flow reconstruction with DOCT priors in 10% of the vessels reduced the MSE to 50% of the MSE obtained when no DOCT priors were used. When the vessels with priors included 1, 2 and 3 outflowing venules, the MSE was further reduced to 45, 42, and 39% of the MSE obtained with no DOCT prior, respectively. The same reduction to 39% of the initial MSE was obtained when the vessels with priors included the inflowing arteriole. Constraining the flow value in the inflowing arteriole or an outflowing venule will null the impact of the pressure boundary condition for this specific vessel and several downstream vessels (or upstream vessels for the case of a venule) since flow must be conserved along the vascular tree. We note that constraining the flow in a single outflowing venule still reduces the MSE compared to the case where the DOCT priors were located only in the capillaries. Moreover, constraining the flow in the inflowing arteriole or all the outflowing venules will fix the regional blood flow (or the global perfusion) of the network and results in global lower MSE. These results indicate that getting DOCT measurements in large inflowing arterioles or large outflowing venules will have the most impact on the flow computation.

Simulations were run with different noise levels and different fractions of segments containing DOCT priors. Averaged results are presented in Fig. 3. Each data point corresponds to an average over 500 simulations. Panel (a) corresponds to simulations where the vessels with DOCT priors were located in the capillaries only, whereas panel (b) corresponds to simulations where the vessels with DOCT priors included either the inflowing arteriole or all three outflowing venules. Note that the top line corresponds to DOCT priors only in the inflowing artery or in all the outflowing veins (no priors in capillaries). Greater reduction in MSE was obtained by increasing the number of vessel segments with DOCT priors compared with reducing the DOCT noise. We also note that the difference obtained by including the inflowing arteriole or the outflowing venules is more meaningful if the number of DOCT priors is lower. Specifically, the improvement doubled when 2% of the segments contained DOCT priors (from 12% to 24%) compared to a less meaningful difference (from 50% to 60%) when 10% of the segments contained DOCT measurements.

Fig. 3

Averaged results of the simulations. The improvement in MSE with DOCT priors is shown for different levels of noise and different fraction of segments containing DOCT priors. (a) Vessels with DOCT priors were located only the capillaries. (b) Vessels with DOCT priors included the inflowing arteriole or all three outflowing venules.

4.2.

Reconstruction of μCBF from Experimental Measurements on a Mouse

The FITC-labeled vasculature measured with TPLSM is shown in Fig. 4(a) where a maximum intensity projection is displayed over the first $300\mu \mathrm{m}$ of the mouse cortex. This vascular stack and the DOCT measurements were previously published in Sakadžić et al.²¹ and Gagnon et al.²² A high contrast-to-noise ratio can be appreciated and this feature allowed the vascular network to be graphed (see Sec. 3). A slice of the corresponding DOCT velocity projection volume is shown in Fig. 4(b). Segments colored in blue correspond to blood flowing down into the brain (usually penetrating arterioles) while segments colored in red correspond to blood flowing out of the brain (generally ascending venules), although in this example we do see undulations of large surface vessels. Nine outflowing venous segments contained DOCT flow measurements. From our simulation results [the top line of Fig. 3(b)], we estimate the MSE in the reconstructed blood flow to be 80% of the MSE that we would have obtained without DOCT measurements. This corresponds to an improvement of 20%.

Fig. 4

Reconstruction of flow in the mouse cortex. (a) Two-photon microscopy FITC angiogram. (b) Doppler OCT velocity image. (c) $\mathrm{\mu }\mathrm{CBF}$ reconstructed from the model with constraints from the DOCT measurements.

The resulting $\mathrm{\mu }\mathrm{CBF}$ distribution reconstructed from the multimodal approach (i.e., using the DOCT constraints in ascending venules) is shown in Fig. 4(c). Values of 15 to $25\mathrm{nL}/\min$ were reconstructed for the absolute $\mathrm{\mu }\mathrm{CBF}$ in the ascending venules which is in agreement with values previously reported for rats anesthesized with isofluorane (20 to $150\mathrm{nL}/\min$ and 25 to $130\mathrm{nL}/\min$).^15,16 The reconstructed $\mathrm{\mu }\mathrm{CBF}$ was larger in penetrating arterioles compared to ascending venules since the number of penetrating arterioles was lower than the number of ascending venules—a feature of the mouse cortical microvasculature previously reported in the literature.²⁰ We note that the reconstructed $\mathrm{\mu }\mathrm{CBF}$ in the capillary bed was very heterogeneous, which is consistent with recent measurements of capillary flux across the rat cortex.²⁴

Our multimodal approach enables us to estimate blood pressure at every branching point of the vascular network and blood velocity in all segments. The reconstructed blood pressure drop from arterioles to venules is shown in Fig. 5(a) while the velocity profile is plotted in Fig. 5(b). The values reconstructed are superimposed with experimental values of blood pressure measured on cats¹⁸ and blood velocity measured on mice.²⁵ As shown in Fig. 5, a good agreement was obtained between the reconstructed values and the values previously reported. Our velocities for capillaries are also in good agreement with values reported in the rat cortex (between 0.5 and $1.0\mathrm{mm}/\mathrm{s}$).^25–28 Sources of error for capillary velocities include the estimation of the diameter of the capillaries. The mean diameter of the capillaries was $5\mu \mathrm{m}$ in our experimental vascular network.

Fig. 5

(a) Blood pressure and (b) blood velocity reconstructed in the mouse cortex from the multimodal data set as a function of vessel type and diameter. Error bars represent standard error on the mean. Values from Lipowsky¹⁸ and Santisakultarm et al.²⁵ are also displayed for comparison.

The reconstructed blood pressure dropped sharply from the arteriole to the first branches of the precapillary arterioles indicating that arterioles are a major source of resistance in this microvascular network. The velocity profile decreased from the arterioles to the capillary bed and increased from the capillaries to the venules. This feature is explained by the fact that the total cross-section of the vasculature is larger in the capillary bed compared to arterioles and venules. This behavior was also obtained in previous works simulating blood flow in two-dimensional synthetic vascular networks.⁶ The velocity in the venules stayed lower compared to the arterioles because of the higher cross-sectional area of venules compared to arterioles in a mouse.

4.3.

Calculation of Regional Blood Flow

The reconstructed $\mathrm{\mu }\mathrm{CBF}$ presented here in Fig. 4 can be used to calculate regional blood flow of the entire vascular stack. This can be achieved by summing either the flow of all inflowing vessels or the flow of all outflowing vessels and assuming or measuring the cortical thickness. This can be viewed as a forward or direct problem and is very straightforward to perform. Assuming a cortical thickness of 1.2 mm,²⁹ the result for the stack of Fig. 4 was $88.7\mathrm{mL}/100\mathrm{g}\mathrm{tissue}/\min$ which is good agreement with values reported in the literature.³⁰ Back-calculating $\mathrm{\mu }\mathrm{CBF}$ given a regional blood flow measurement represents the corresponding inverse problem and is more difficult to perform. One thing to note is that back-calculating flow in individual inflowing or outflowing vessels given a regional blood flow measurement does not require the entire reconstruction of $\mathrm{\mu }\mathrm{CBF}$ as presented in this work. In the case where a single inflowing arteriole feeds the vascular stack, one can directly compute the flow in the single inflowing arteriole from its diameter (assuming or given a measurement of the cortical thickness). The same thing happens if a single outflowing venule drains the entire stack. However, as the number of inflowing arteries or outflowing venules become $n>1$, there exist multiple solutions to the inverse problem. In this case, flow measurements in $n-1$ individual inflowing or outflowing vessels are needed as well as diameter measurements in the remaining inflowing or outflowing vessel.

4.4.

Potential Applications

The change in CBF during healthy aging is an active area of research.^27,31 The numerical method presented in this work has phe otential to investigate age-related differences in CBF at the microvascular level based on measurements of relevant CBF parameters.²⁷ Specifically, our method will allow us to reconstruct $\mathrm{\mu }\mathrm{CBF}$ with better accuracy in the aged rodent cortex which might reveal differences in microvascular circulation between young and aged mice. This will potentially allow us to make predictions about macroscopic flow parameters that can be measured noninvasively in humans studies.

Although our multimodal reconstruction method reduces the error on flow reconstruction by 50%, our simulations indicate that the residual error is still on the order of 150%. Most of this error comes from the capillary bed and is dominated by uncertainty in capillary diameter estimation. This high level of noise in the reconstructed $\mathrm{\mu }\mathrm{CBF}$ limits the application of our method to study phenomena that only modify capillary blood flow distribution.³ However, the residual error in the regional blood flow is lower and is of the order of the noise in the experimental DOCT measurements in the inflowing or outflowing vessels. Our method can, therefore, be used to study phenomena which modify regional blood flow such as aging.^27,31

4.5.

Limitations and Future Studies

As mentioned before, the $X\text{-}Y$ resolution of our OCT system was 10 to $15\mu \mathrm{m}$, which did not permit measurement of flow in the capillary bed. New methods were recently developed^24,32–34 to measure $\mathrm{\mu }\mathrm{CBF}$ in the capillary network. These measurements could be included in our computational framework in future studies. Another limitation of our system is the low readout rate of our OCT camera limiting the measured velocity projections to $\pm 3.6\mathrm{mm}/\mathrm{s}$. This can potentially create phase wrapping in large arterioles, therefore, we limited our DOCT measurements to venules. A faster OCT system like the one used in Srinivasan et al.¹⁶ could resolve this issue and allow us to constrain flow reconstruction in both arterioles and venules.

It is worth mentioning that the simulations performed in the first part of this work did not assume any specific modality for the measurement of CBF in individual vessels. Rather, the simulations focused on the improvement that can potentially be obtained by constraining the flow reconstruction with any experimental flow measurement. In our case, we chose to use DOCT for the flow measurements since our multimodal microscope allowed for such experiments. Nevertheless, one could have chosen to perform the flow measurement with TPLSM line scan measurements²⁶ instead. We will mention here that DOCT has some advantages over TPLSM such as faster acquisitions.¹⁵

New experimental methods permit capillary flow speed measurements on several vessels simultaneously.^24,34,35 These methods will allow us to validate our computational method in future studies and potentially improve it by providing more measurements that can be used to constrain the flow reconstruction. Moreover, these measurements will allow us to measure temporal variation in $\mathrm{\mu }\mathrm{CBF}$ while our computational method assumes a steady-state flow.