2.1.

Correlation Diffusion Equation

DCS measurements of blood flow are generally analyzed using the CDE.^2,3 The CDE is derived from the correlation transfer equation (CTE)^3,16 under the assumption that the probability of light scattering is much greater than the probability of light absorption, i.e., reduced scattering coefficient (${\mu }_s^{\prime }$) is much greater than the absorption coefficient (${\mu }_a$). The CTE¹⁷ permits calculation of the electric field temporal autocorrelation function $G_1(\tau )$ under more general conditions of photon migration. The CTE is similar to the traditional radiative transfer equation (RTE), which describes the propagation of light intensity through scattering media. The difference between the CTE and the RTE is that the CTE describes the propagation of the time-varying specific intensity, which represents an angular spectrum of the mutual coherence function. The temporal autocorrelation function is obtained as an integral of the time-varying specific intensity over all solid angles.

Generally, an analytic solution of the CDE for a semi-infinite medium is used to fit the experimental measurements, i.e.,

It is customarily expected that the scatterer motion measured by DCS arises almost exclusively from RBCs. This is justified by experimental studies in animals, which showed that cessation of blood flow results in a nearly 100-fold reduction in the observed dynamics.⁴ In this paper, we will assume that RBCs are the only source of dynamic scattering events, and we will use $P_{\mathrm{RBC}}$ instead of $\alpha$ to denote the probability of scattering from an RBC.

The mean square displacement (MSD) $\langle \mathrm{\Delta }{r}^2(\tau )\rangle$ of the RBCs is generally given by

2.2.

Monte Carlo Simulations

Investigations with DCS have generally assumed that the motion of tissue scatterers is uncorrelated. Specifically, the assumption is that the phase shifts accumulated by the motions of scattering particles is uncorrelated, which depends on the motion of the scattering particle and the scattering angle. This assumption is generally valid when photons are scattering from particles undergoing Brownian motion. In the case of moving RBCs, this assumption is valid if the photon direction is randomized between scattering events from RBCs. This is true if the scattering is isotropic (i.e., the average cosine of the scattering angle is $g=0$) or if the photons scatter no more than once inside a vessel and the direction is randomized before scattering from another RBC in another vessel. The optical properties of blood at 800 nm indicate,²¹ however, that for a typical hematocrit of 40%, $g\approx 0.98$ and the scattering length $l_s$ is $\sim 12\mu \mathrm{m}$, whereas the absorption length $l_a$ is over 250 times larger, at around 3.3 mm. Thus, photons entering any blood vessel larger than a capillary will most likely undergo multiple consecutive scattering events within the vessel with correlated linear motions of the RBCs. This breakdown of the uncorrelated motion assumption raises the possibility that the typically used solutions of the CTE, and thus the CDE, are not good models for analyzing DCS measurements of blood flow.

To properly account for the effect of correlated motions on the decay of $G_1(\tau )$, we use Monte Carlo simulations of photon migration through highly scattering dynamic media and record details of the scattering angles for consecutive scattering events within vessels. Our Monte Carlo code is derived from what Boas originally utilized in Ref. 16 based on earlier work by Middleton and Fisher,²² Durian,²³ and Koelink et al.²⁴ The original Monte Carlo simulation propagates photons according to the scattering length and scattering anisotropy of the tissue, which we describe in more detail in Ref. 25. Photons are launched from a source position and recorded at discrete detector locations. For each detected photon, the total path length $L$ is recorded as well as the total dimensionless momentum transfer accumulated over all scattering events $Y$. The momentum transfer for a photon scattering event is given by $\overrightarrow{q}={\overrightarrow{k}}_{\text{out}}-{\overrightarrow{k}}_{\text{in}}$, where ${\overrightarrow{k}}_{\text{out}}$ and ${\overrightarrow{k}}_{\text{in}}$ are the photon wavevectors scattered from and incident on the scattering center, respectively, and the total dimensionless momentum transfer is given by $Y=\sum _{}^{}{q}^2/2k_0^2$, where the sum is over scattering events. The temporal field autocorrelation function is then calculated by

The principal advance in the Monte Carlo simulations used in this paper is the introduction of blood vessels with different optical and dynamic properties than the surrounding tissue and with specified internal flow profiles, as shown in Figs. 1 and 2. Photon migration is simulated in a semi-infinite medium with the air tissue interface in the ($x,y$) plane at $z=0$. The blood vessels share a common radius $\alpha$ and are all oriented along the $y$-axis with an even spacing in $x$ and $z$ of $h_{\text{space}}$. The resulting translational symmetry allows for significant savings in the amount of memory needed to represent the tissue structure for the Monte Carlo simulation. We have verified that this preferred orientation of the vessels only introduces a small bias in our results compared to a truly random orientation of vessels, as expected, given that we typically simulate photons that have propagated a net-distance of several centimeters and have scattered hundreds of times.

Fig. 1

Simulation geometry used in this study. (a) Three-dimensional view and (b) ($x,z$) plane cross section. Vessels of radius $R$ are aligned with the $y$-axis and are uniformly spaced in the ($x,z$) plane every $h_{\text{space}}$ $\mu \mathrm{m}$. Photons are injected at a location within the ($x,y$) plane and detected a distance $\rho$ away.

Fig. 2

Photon propagation tracking. For each scattering event, we record the photon location and the magnitude of the momentum transfer $|\overrightarrow{q}|$. If the scattering event occurs inside a vessel, we also record the radial position $r$ and the component of $\overrightarrow{q}$ along the vessel direction. In conjunction with a specified flow profile, the radial position $r$ determines the RBC convection velocity $v_{\mathrm{RBC}}(r)$ and the shear rate experienced by the RBC $s(r)$.

Photons are launched at a point on the surface along the $z$-axis and are detected through a 2-mm-diameter circular aperture, a distance $\rho$ from the source point. Figure 2 shows a segment of a photon path and the parameters we record during the photon propagation. For each scattering step during the simulation, we explicitly test for the photon crossing the boundary of a cylindrical vessel. If a crossing occurs, the remaining propagation length for the scattering step is renormalized based on the new scattering coefficient as we describe in Ref. 25. For each detected photon, we record for every scattering event inside a vessel the index of the vessel, the radial position of the scattering event, the magnitude of the momentum transfer $|\overrightarrow{q}|$, and the component of $\overrightarrow{q}$ along the vessel direction ($y$-axis). Recording the vessel index enables us to determine consecutive scattering events within a vessel. Knowing the radial position of the scattering event permits us to consider a radial distribution for the convective and diffusive dynamics of the RBCs. We record the $y$ component of $\overrightarrow{q}$ as it is needed to consider the correlated convective motion of the RBCs along the $y$-axis. All Monte Carlo simulations in this paper use a source–detector separation of 2 cm, an ($x,z$) cross section of $6\times 3\mathrm{cm}$, and make use of translational symmetry to represent an infinite extent in $y$. Each simulation run launches ${10}^8$ photons.

To calculate the temporal field autocorrelation function, we start with the definition of $G_1(\tau )$ and expand it to sum over all detected photons

We are considering the convective and shear-induced diffusive movement of RBCs such that the MSD of the RBCs is given by Eq. (3). Although the strength of the diffusive motion correlates with the magnitude of the local shear, the actual diffusive displacements are spatially uncorrelated with the convective displacements. Therefore, note that ${\overrightarrow{q}}_{n,m}·\overrightarrow{\mathrm{\Delta }}{r}_{n,m}(\tau )$ is generally much smaller than 1, and Eq. (5) can be rewritten as

2.3.

Red Blood Cell Laminar Flow Profile and Shear-Induced Diffusion

Our Monte Carlo simulations permit us to incorporate radial distributions for the flow speed and diffusion coefficient. A general form of the flow speed profile in vessels has been proposed to be²⁶

RBC dynamics, while flowing within a vessel, have been studied in numerous rheological studies including video microscopy to track the motion of RBCs ex vivo¹⁵ and whole blood in microchannels,²⁷ as well as in vivo in rat venules.²⁸ These studies have revealed that shear induces a diffusive behavior in the movement of the RBCs with the diffusion coefficient proportional to the shear rate, i.e.,

Of note, a prior study by Wu et al.²⁹ has shown that the gradient in the speed of scattering particles flowing within a channel can be the leading factor in the decay of the autocorrelation function. This is an informative study, but the parameters used are quite different from those encountered in biological tissue, leading to negligible levels of shear-induced diffusion. This is due to the use of a single flow channel with a much smaller shear rate than that experienced by RBCs in blood vessels, as well as a much smaller particle concentration and scattering coefficient. Furthermore, this study did not consider multiple flow channels in an otherwise static scattering medium.

2.4.

Relating the Diffuse Correlation Spectroscopy Blood Flow Index to Absolute Blood Flow

We are now in a position to perform Monte Carlo simulations to calculate $g_{1,\mathrm{MC}}(\tau )$ for RBCs flowing at varying speeds in vessels of varying diameters and spacings. We can then fit the result with $g_{1,\mathrm{CDE}}(\tau )$ obtained from the CDE to determine if the CDE remains a good model and to establish a quantitative relationship between the CDE blood flow index ${\mathrm{BF}}_i$ and the absolute blood flow ${\mathrm{BF}}_{\text{abs}}$. For a semi-infinite medium, we use $g_{1,\mathrm{CDE}}(\tau )$ given in Eq. (1) and estimate $K$ from Eq. (2). The optical properties and wavenumber of light are generally known in Eq. (2), leaving us to estimate $P_{\mathrm{RBC}}\langle \mathrm{\Delta }{r}^2(\tau )\rangle$. We define the blood flow indices ${\mathrm{BF}}_i^{\mathrm{conv}}=\sqrt{P_{\mathrm{RBC}}\langle v^2\rangle }$ and ${\mathrm{BF}}_i^{\mathrm{dif}}=P_{\mathrm{RBC}}\langle D\rangle$ to be estimated by fitting either the convective or diffusive MSD model to $g_{1,\mathrm{MC}}(\tau )$, where either ${\mathrm{BF}}_i$ is expected to be linearly proportional to flow speed. The angle brackets for $⟨v⟩$ and $⟨D⟩$ are an explicit recognition that the ${\mathrm{BF}}_i$ measures the average speed $v$ or diffusion coefficient $D$. Note that the probability of scattering from an RBC is not generally known a priori. In the context of photon diffusion theory, we would expect that $P_{\mathrm{RBC}}$ is given by the product of the volume fraction of blood, $f_{\text{blood}}=\pi R^2/h_{\text{space}}^2$ and the ratio of the blood to average tissue reduced scattering coefficient, ${\mu }_{s,\text{blood}}^{\prime }/{\mu }_{s,\mathrm{avg}}^{\prime }$, i.e.,

We present results below for fitting $g_{1,\mathrm{MC}}(\tau )$. To establish the quantitative relationship between ${\mathrm{BF}}_i$ and ${\mathrm{BF}}_{\mathrm{abs}}$ for convective and diffusive RBC displacements independently, we obtain fitting results alternately setting $D=0$ and $v_{\mathrm{RBC}}=0$, respectively. We will distinguish these temporal field autocorrelation functions as $g_1^{\mathrm{conv}}(\tau )$ and $g_1^{\mathrm{dif}}(\tau )$ for convective and diffusive RBC displacements, respectively.

In the case of convective motion only, we expect $⟨v⟩$ to be the spatially weighted average of RBC speed within the blood vessel. As the absorption length and photon transport length inside the vessel are generally larger than the vessel diameter, we expect uniform spatial sampling of the vessel and thus expect $\langle v\rangle =v_{\text{max}}/2$ for a parabolic laminar flow profile. We thus arrive at an equation for ${\mathrm{BF}}_i^{\mathrm{conv}}$

In the case of diffusive motion only, we expect $⟨D⟩$ to be the spatially weighted average of the RBC diffusion coefficient within the blood vessel, i.e., $\langle D\rangle =4{\alpha }_{\text{shear}}{v}_{\text{max}}/3R$ given a parabolic laminar profile and uniform photon sampling of the blood vessel. We thus arrive at an equation for ${\mathrm{BF}}_i^{\mathrm{dif}}$

We define absolute blood flow as the volume of blood transiting through a unit volume of tissue per second, which is given by the flow in a vessel (cross-sectional area times average RBC speed) times the number of vessels per unit volume of tissue, i.e.,

Finally, we present results below considering the combination of convective and diffusive RBC displacements to establish which dominates the decay of $g_1^{\text{total}}(\tau )$. Given the overwhelming indication from experimental studies,⁴ we expect to observe that the diffusive motion dominates the decay.

3.1.

Impact of Scattering Length Within a Vessel

Table 1 summarizes representative optical properties of blood and extravascular tissue at 800 nm for a blood hematocrit of 40%.²¹ We use $g$ and the scattering coefficient ${\mu }_s=1/l_s$ in our Monte Carlo simulations. For the analytical solution provided by the CDE, we use the reduced scattering coefficient ${\mu }_s^{\prime }={\mu }_s(1-g)$.

Table 1

Optical properties for blood and extravascular tissue at 800 nm.

As a preamble, we explicitly test the prediction offered by the CDE that $g_1(\tau )$ is invariant to changes in ${\mu }_s$ and $g$ provided that ${\mu }_s^{\prime }$ is constant. To this end, we ran Monte Carlo simulations for three different combinations of ${\mu }_{s,\text{blood}}$ and $g_{\text{blood}}$ maintaining a constant ${\mu }_{s,\text{blood}}^{\prime }=1.9{\mathrm{mm}}^{-1}$. Specifically, we use ${\mu }_{s,\text{blood}}=1.9$, 19, and $190{\mathrm{mm}}^{-1}$, and $g_{\text{blood}}=0,0.9$, and 0.99, respectively. The rest of the optical parameters match the values in Table 1, except that we use ${\mu }_{s,\mathrm{tis}}^{\prime }=1.0{\mathrm{mm}}^{-1}$ and $g_{\text{tis}}=0$, which are typical values we use in our Monte Carlo simulations for the human head to accelerate the computation.²⁵ We show results for $80\text{-}\mu \mathrm{m}$-diameter vessels, spaced every $300\mu \mathrm{m}$, with $v_{\mathrm{RBC},\max }=2\mathrm{mm}/\mathrm{s}$.

The results for $g_1^{\mathrm{conv}}(\tau )$ and $g_1^{\mathrm{dif}}(\tau )$ are shown in Fig. 3. These results confirm the prediction from the CDE that, at least for early decay times, $g_1(\tau )$ is invariant as long as ${\mu }_s(1-g)$ remains constant. The behavior at long decay times does exhibit some dependence on ${\mu }_s$ of the blood [see enlargement in panel b)]. Likely, as ${\mu }_s$ of the blood is decreased, the probability increases that a photon traverses a vessel without scattering from a moving RBC. Importantly, we observe that $g_1^{\mathrm{dif}}(\tau )$ decays more quickly than $g_1^{\mathrm{conv}}(\tau )$, indicating that the shear-induced diffusion of the RBCs produces more dynamic fluctuations in the electromagnetic field than the flow speed convective displacements. Based on these results, for the remainder of the paper, we use $g_{\text{blood}}=0.9$ since this gives us effectively the same results we would obtain with $g_{\text{blood}}=0.977$ and leads to a more efficient execution of the Monte Carlo simulations.

Fig. 3

(a) The decay of the temporal field autocorrelation function at a 2-cm source–detector separation is shown for ${\mu }_{s,\text{blood}}^{\prime }=1.9{\mathrm{mm}}^{-1}$ with $g=0,0.9$, and 0.99. We used $80\text{-}\mu \mathrm{m}$-diameter vessels, spaced every $300\mu \mathrm{m}$, with $v_{\mathrm{RBC},\max }=2\mathrm{mm}/\mathrm{s}$. The dotted, dashed, and solid lines indicate $g_1^{\mathrm{conv}}(\tau )$, $g_1^{\mathrm{dif}}(\tau )$, and $g_1^{\mathrm{tot}}(\tau )$, respectively. (b) An enlargement of the ${10}^{-5}<\tau <2\times {10}^{-4}$ range in (a).

3.2.

Blood Flow Indices Estimated from Fitting Monte Carlo Data Using the Correlation Diffusion Equation

In this section, we present results to validate the relations for the blood flow indices ${\mathrm{BF}}_i^{\mathrm{conv}}$ and ${\mathrm{BF}}_i^{\mathrm{dif}}$ given by Eqs. (10) and (11), respectively. We do so by fitting the solution of the CDE, given by Eq. (1), to $g_{1,\mathrm{MC}}(\tau )$, obtained from the Monte Carlo simulations. As detailed in Table 2, we varied RBC speed, vessel radius, vessel spacing, and hematocrit in these simulations to fully test the dependencies of ${\mathrm{BF}}_i^{\mathrm{conv}}$ and ${\mathrm{BF}}_i^{\mathrm{dif}}$ on each parameter. Note that we fit $g_1^{\mathrm{conv}}(\tau )$ to estimate ${\mathrm{BF}}_i^{\mathrm{conv}}$ and we fit $g_1^{\mathrm{tot}}(\tau )$ to estimate ${\mathrm{BF}}_i^{\mathrm{dif}}$. We only fit $g_1^{\mathrm{tot}}(\tau )$ with the diffusion term because our simulations indicate that the ${\tau }^2$ dependence arising from the convective displacement of the RBCs has negligible impact on $g_1^{\mathrm{tot}}(\tau )$ for our parameters of interest. For all cases presented in this section, the blood optical properties were ${\mu }_s^{\prime }=1.9{\mathrm{mm}}^{-1}$, $g=0.9$, and ${\mu }_a=0.30{\mathrm{mm}}^{-1}$. For the tissue, we set ${\mu }_s^{\prime }=1.0{\mathrm{mm}}^{-1}$, $g=0.0$, and ${\mu }_a=0.002{\mathrm{mm}}^{-1}$.

Table 2

Vascular parameter ranges explored in this study. Each parameter was varied separately while the others were held at a fixed value.

We first present results increasing the speed $v_{\text{max}}$ of the RBCs. Figures 4(a) and 4(b) show $g_{1,\mathrm{MC}}(\tau )$ and the CDE fit for $v_{\text{max}}=\mathrm{1,2}$, and $4\mathrm{mm}/\mathrm{s}$ for $g_1^{\mathrm{tot}}(\tau )$ and $g_1^{\mathrm{conv}}(\tau )$, respectively. In this example, the vessel diameter is set to $80\mu \mathrm{m}$, the vessels are spaced every 300 μm and the hematocrit is set to 40%. These results indicate that the CDE fits $g_1^{\mathrm{tot}}(\tau )$ and $g_1^{\mathrm{conv}}(\tau )$ well. Figures 4(c) and 4(d) show the fitted values of ${\mathrm{BF}}_i^{\mathrm{dif}}$ and ${\mathrm{BF}}_i^{\mathrm{conv}}$ versus speed, compared with those predicted using Eqs. (10) and (11). These results confirm that the ${\mathrm{BF}}_{i^{}}$s increase linearly with RBC speed and confirm the accuracy of Eqs. (10) and (11).

Fig. 4

(a) and (b) The CDE fit to the Monte Carlo results for $g_1^{\mathrm{tot}}(\tau )$ and $g_1^{\mathrm{conv}}(\tau )$. The Monte Carlo results for different RBC speeds are indicated by the symbols and the diffuse equation fits are given by the colored solid lines. The symbols in (c) and (d) show the blood flow index ${\mathrm{BF}}_i$ determined from fitting $g_1^{\mathrm{tot}}(\tau )$ with a diffusive motion model and from fitting $g_1^{\mathrm{conv}}(\tau )$ with a convective motion model, respectively, versus the true RBC speed. The solid lines in (c) and (d) are the theoretical predictions for the ${\mathrm{BF}}_{i^{}}$s based on Eqs. (10) and (11). Here, we used 80 μm vessels spaced every 300 μm and set the hematocrit to 40%.

Figure 5 shows the ${\mathrm{BF}}_i$ fitting results for variations in the absolute volumetric blood flow achieved by varying the vessel diameter at a constant maximum RBC speed $v_{\text{max}}$. The vessel diameter was varied from 30 to $100\mu \mathrm{m}$ and results are shown for vessel spacings of 200 and 300 μm, respectively. As with the previous results, we fit $g_1^{\mathrm{tot}}(\tau )$ and $g_1^{\mathrm{conv}}(\tau )$ to estimate ${\mathrm{BF}}_i^{\mathrm{dif}}$ and ${\mathrm{BF}}_i^{\mathrm{conv}}$, respectively. We then used Eqs. (14) and (15) to determine the expected ${\mathrm{BF}}_i$ values and find that this predicts well the measured ${\mathrm{BF}}_i$.

Fig. 5

The symbols in (a) and (b) show the diffusive and convective blood flow indices, respectively, estimated using the CDE versus the true blood flow. The blood flow variation was accomplished by increasing the vessel diameter from 30 to $100\mu \mathrm{m}$ with a fixed RBC flow speed of $1\mathrm{mm}/\mathrm{s}$. The solid lines indicate the expected blood flow index values predicted by Eqs. (15) and (14), respectively. These results were obtained for a hematocrit of 40%

In Fig. 6, we plot fitted ${\mathrm{BF}}_i$ values for changing the blood hematocrit by $\pm 20\%$ from our baseline value of 40%. In this range, the reduced scattering coefficient of blood is changing linearly with hematocrit.²¹ We also show the expected ${\mathrm{BF}}_i$ values and find excellent agreement.

Fig. 6

The symbols in (a) and (b) show the diffusive and convective blood flow indices, respectively, estimated using the CDE for vessel scattering coefficients corresponding to hematocrit values ranging from $-20\%$ to $+20\%$ of our 40% baseline for a vessel diameter of $80\mu \mathrm{m}$, with $300\text{-}\mu \mathrm{m}$ vessel spacing, and an RBC flow speed of $2\mathrm{mm}/\mathrm{s}$. The solid lines indicate the expected blood flow index values predicted by Eqs. (15) and (14), respectively.