Expand this Topic clickable element to expand a topic
Skip to content
Optica Publishing Group

Numerical approach to quantify depth-dependent blood flow changes in real-time using the diffusion equation with continuous-wave and time-domain diffuse correlation spectroscopy

Open Access Open Access

Abstract

Diffuse correlation spectroscopy (DCS) is a non-invasive optical technique that can measure brain perfusion by quantifying temporal intensity fluctuations of multiply scattered light. A primary limitation for accurate quantitation of cerebral blood flow (CBF) is the fact that experimental measurements contain information about both extracerebral scalp blood flow (SBF) as well as CBF. Separating CBF from SBF is typically achieved using multiple source-detector channels when using continuous-wave (CW) light sources, or more recently with use of time-domain (TD) techniques. Analysis methods that account for these partial volume effects are often employed to increase CBF contrast. However, a robust, real-time analysis procedure that can separate and quantify SBF and CBF with both traditional CW and TD-DCS measurements is still needed. Here, we validate a data analysis procedure based on the diffusion equation in layered media capable of quantifying both extra- and cerebral blood flow in the CW and TD. We find that the model can quantify SBF and CBF coefficients with less than 5% error compared to Monte Carlo simulations using a 3-layered brain model in both the CW and TD. The model can accurately fit data at a rate of <10 ms for CW data and <250 ms for TD data when using a least-squares optimizer.

© 2022 Optica Publishing Group under the terms of the Optica Open Access Publishing Agreement

1. Introduction

Adequate cerebral blood flow (CBF) is required for delivery of oxygen and other nutrients to brain tissue, allowing CBF to serve as a biomarker for brain function and health [1]. The quantitation of CBF is therefore important for managing and diagnosing brain injury or other diseases associated with ischemia or inadequate vascular autoregulation [2]. Diffuse correlation spectroscopy (DCS) is an optical technology that noninvasively measures blood flow by quantifying temporal speckle fluctuations of a coherent light field illuminating the tissue surface, usually in the near-infrared spectrum [35]. Detected photons in such applications must propagate through the extra-cerebral scalp and skull layers to reach brain tissue [6,7]. The backscattered light intensity temporally fluctuates and is quantified via its temporal autocorrelation function [2,8]. Analytical solutions of the correlation equation are then obtained using diffusion theory and are used to fit measured autocorrelation functions to extract a quantitative blood flow index (BFi, $\mathrm {cm}^2/$s), which has been shown to be a good surrogate for CBF measured in vivo [2,3,911].

Traditionally, DCS systems have used continuous-wave (CW) light sources which require long source-detectors separations to increase sensitivity to CBF [12,13] as the CW-DCS signal is significantly impacted by shallower extracerebral scalp blood flow (SBF) changes [2,6]. More recently, using picosecond pulsed lasers as sources has allowed DCS to operate in the time-domain (TD), which can increase sensitivity to deeper tissues by gating the photons time of flight to measure the autocorrelation at different time windows [1417]. Several studies have demonstrated the potential advantage of TD-DCS compared to CW-DCS for measuring deeper flow changes at shorter source-detector separations [14,1618]. However, although TD-DCS signals can select for longer pathlengths, these photons still traverse the superficial layers at least twice which must be accounted for in data analysis [19].

A major challenge in using DCS for neuromonitoring is that the detected light intensity interacts with all of the tissue layers in the propagation path (e.g., scalp, skull, cerebrospinal fluid (CSF), brain) [7,9,12,20]. Homogeneous diffusion models are typically used for their simplicity and speed, but extracted BFi metrics have been shown to non-linearly depend on the flow properties of each layer [12]. Several analysis methods have been developed to account for these partial volume effects by using differential path-length factors [2123], selecting longer photon pathlengths by weighting fits to shorter correlation times [6,22], or by using pressure modulation procedures [24,25]. Although useful for enhancing sensitivity to CBF, these methods do not eliminate top layer flow contributions as they still assume some level of tissue homogeneity [23]. Alternatively, using heterogeneous brain models simulated with Monte Carlo (MC) methods [6,7] or analytical solutions to the multi-layered correlation diffusion equation [9,20,22,23,26,27] have shown promise for more robustly separating SBF from CBF [7,9,20].

The increase in available computational resources has led to the MC approach becoming more popular due to its higher accuracy at short source-detector separations and ability to accommodate complex tissue geometries [7,19,28]. Additionally, the MC method may be more robust than layered analytical models which can be numerically unstable [7]. However, for DCS to be employed to monitor patients in clinical settings, data analysis tools must be able to continuously and accurately quantify CBF changes for rapid treatment in the clinic [7]. A major limitation of MC is its high computational cost which mostly precludes its use in real-time data analysis [7,29,30].

Layered analytical models can provide a good compromise between reconstructive accuracy [31] and computational speed, yet they still remain significantly slower to compute than homogeneous models and are not typically used in inverse problems [7]. Forward models for TD-DCS in layered media have also been derived [32], but to the best of our knowledge have not yet been applied as inverse solvers, as TD solutions require orders of magnitude higher computational effort than CW approaches [33]. Currently, most applications employ homogeneous tissue models with diffusion theory for both CW- and TD-DCS [14], or use heuristic exponential based fitting for quantifying the BFi [19].

In this report, we develop and apply diffusion theory based analytical solutions of DCS signals in layered media to overcome existing challenges and validate its: (1) physical accuracy compared to MC simulations, (2) numerical stability and computational speed, and (3) use in inverse solutions to recover BFi from signals measured in layered head-tissue models. We present the theoretical framework for both CW- and TD-DCS based on the diffusion equation for an N-layered turbid medium [34,35]. The model’s performance is verified by comparisons to MC simulations across a range of source-detector separations for both the CW and TD. We examine using the forward model as an inverse solver and benchmark quantitative recovery of both SBF and CBF from a 3-layered brain model. We conclude by discussing the advantages of CW and TD systems from the results, establishing the sensitivities of CW- and TD-DCS in recovery of SBF and CBF, and investigating the impact of source-detector separations and delay times on these sensitivities.

2. Methods

2.1 Diffusion theory in layered media

The motion of scattering particles is associated with the unnormalized electric field temporal autocorrelation function $G_1$ which in highly scattering turbid media such as biological tissue can be approximated with the correlation diffusion equation [3]. Solutions to the correlation diffusion equation can be obtained by first solving the regular photon diffusion equation

$$D \nabla^2 \Phi(\vec{r}) - \mu_a \Phi(\vec{r}) ={-}S(\vec{r})$$
and then replacing the absorption coefficient with a dynamic absorption term $\mu _a^d \rightarrow \mu _a + 2\mu _s'D_B k_0^2 \tau$ [8]. $\Phi$, $D$ , $S(\vec {r})$ and $\mu _a$ denote the fluence rate, the photon diffusion coefficient, the source configuration, and the absorption coefficient, respectively, where $D=(3\mu _s')^{-1}$ and $\mu _s'$ the reduced scattering coefficient. $D_B$ is the Brownian flow-diffusion coefficient, $k_0 = 2 \pi n / \lambda$ is the wavenumber and $\lambda$ is the wavelength of the light source [36]. It is assumed that the mean-square displacement of the scattering particles over time $\tau$ undergoes Brownian motion (i.e., $\langle \Delta r^2 (\tau ) \rangle = 6 D_B \tau$).

Because solutions to the photon diffusion equation can be used to calculate $G_1$ by a simple replacement, existing methods to solve the diffusion equation for specific geometries and boundary conditions can be extended for calculating autocorrelation using diffusion theory [8,14]. Here, we extend our computationally efficient method [35] to solve the diffusion equation using extrapolated boundary conditions for fast calculation of the autocorrelation function in an N-layer tissue model made up of axially layered cylinders, each with independent $\mu _a$, $\mu _s'$, and thickness $l$ [34]. Assuming an incident beam onto the center of the topmost cylindrical layer can be approximated by an isotropic point source located at a distance of $z_0 = 1/\mu _s'$ from the incident surface, Eq. (1) can be solved in real space [35] by

$$\Phi(\rho, z = 0) \approx \Phi^{SI}(\rho) + \frac{1}{\pi a'^2} \sum_{n=1}^{\infty} \hat{G}(s_n, z)J_0(s_n \rho)J_{1}^{{-}2}(a' s_n)$$
where $J_m$ is the Bessel function of first kind and order $m$ and $s_n$ is determined from the roots of $J_m$ for $J_m(a' s_n)=0$ where $n=1,2,\ldots$ is the $n^{th}$ root of $J_m$. $z$ is the location of the detector within the medium in cylindrical coordinates and $\Phi ^{SI}(\rho )$ is the steady-state fluence for a semi-infinite medium [37]. An extrapolated boundary is determined using $a' = a + z_b$ where $a$ is the radius of the cylinder and $z_b = 2AD$ where $A$ is proportional to the fraction of photons that are internally reflected at the boundary [38]. For clarity on nomenclature, $G_1$ represents the unnormalized temporal autocorrelation function whereas $\hat {G}$ represents the Green’s function in the top layer due to an isotropic source located in the top layer.

For the solutions presented here, we utilize the approximate forms in Eq. (2) given previously [35] which are accurate when $z \approx z_0$. We have shown that given $\mu _{s1}' > 2 \, \mathrm {cm}^{-1}$ and $a >> \rho$ at $z=0$ these forms are accurate to at least 16 digits which is exact in double precision arithmetic [35]. The expression for the Green’s function $\hat {G}$ in the top layer in Eq. (2) is

$$\begin{aligned} \hat{G}(s_n, z) ={} & \frac{e^{\alpha_1(z+z_0-2l_1)} (1 - e^{{-}2\alpha_1 (z_0 + z_{b1})})(1 - e^{{-}2\alpha_1(z+z_{b1})})}{2D_1 \alpha_1} \\ & \times\frac{D_1 \alpha_1 n_1 ^2 \beta_3 - D_2 \alpha_2 n_2 ^2 \gamma_3}{D_1 \alpha_1 n_1 ^2 \beta_3 (1 + e^{{-}2\alpha_1 (l_1 +z_{b1})}) + D_2 \alpha_2 n_2 ^2 \gamma_3 (1 - e^{{-}2\alpha_1 (l_1 +z_{b1})})} \end{aligned}$$
where $\alpha _k = \sqrt {(\mu _{ak} + 2\mu _{sk}'D_{Bk} k_0^2)/D_k + s_n^2}$. The subscript $k$ signifies the $k^{th}$ layer of the cylinder with distinct absorption $\mu _{ak}$ and scattering $\mu _{sk}'$ in each layer. In general, the quantities $\beta _3$ and $\gamma _3$ are obtained by downward recurrence relations with start values
$$\begin{aligned} \beta_N ={} & D_{N-1} \alpha_{N-1} n_{n-1}^2 (1 + e^{{-}2\alpha_{N-1}l_{N-1}})(1 - e^{{-}2\alpha_{N}(l_{N} + z_{b_2})}) \\ & + D_{N} \alpha_{N} n_{n}^2 (1 - e^{{-}2\alpha_{N-1}l_{N-1}})(1 + e^{{-}2\alpha_{N}(l_{N} + z_{b_2})}) \\ \gamma_N ={} & D_{N-1} \alpha_{N-1} n_{n-1}^2 (1 - e^{{-}2\alpha_{N-1}l_{N-1}})(1 - e^{{-}2\alpha_{N}(l_{N} + z_{b_2})}) \\ & + D_{N} \alpha_{N} n_{n}^2 (1 + e^{{-}2\alpha_{N-1}l_{N-1}})(1 + e^{{-}2\alpha_{N}(l_{N} + z_{b_2})}) \end{aligned}$$
with the downward recurrence given by
$$\begin{aligned} \beta_{k - 1} ={} & D_{k - 2} \alpha_{k - 2} n_{k - 2}^2 (1 + e^{{-}2\alpha_{k - 2}l_{k - 2}}) \beta_k + D_{k - 1} \alpha_{k - 1} n_{k - 1}^2 (1 - e^{{-}2\alpha_{k - 2}l_{k - 2}}) \gamma_k \\ \gamma_{k - 1} ={} & D_{k - 2} \alpha_{k - 2} n_{k - 2}^2 (1 - e^{{-}2\alpha_{k - 2}l_{k - 2}}) \beta_k + D_{k - 1} \alpha_{k - 1} n_{k - 1}^2 (1 + e^{{-}2\alpha_{k - 2}l_{k - 2}}) \gamma_k \end{aligned}$$
We note that we seek just the terms $\beta _3$ and $\gamma _3$ which must be determined recursively if the total number of layers $N$ is larger than 3. In that case, Eq. (4) is used to generate starting values in the recurrence relation, then Eq. (5) is used recursively until $\beta _3$ and $\gamma _3$ are obtained. If $N=2$, $\beta _3 ={} 1 - e^{-2 \alpha _2 (l_2 + z_{b2})}$ and $\gamma _3 ={} 1 + e^{-2 \alpha _2 (l_2 + z_{b2})}$. For $N=3$, only Eq. (4) is needed to calculate $\beta _3$ and $\gamma _3$. Presenting the coefficients in terms of exponentially decaying functions prevents numerical overflow in the inverse transforms and increases numerical stability [35].

For solutions in the time-domain, an inverse Laplace transform is performed using the substitution $\mu _a \rightarrow \mu _a + \bar {s}/c$ with $c$ being the speed of light in the medium. The time-dependent fluence $\Phi (\rho, t)$ can be expressed by

$$\Phi(\rho, t, z=0) \approx \Phi^{SI}(\rho, t) + \frac{1}{\pi a'^2} \sum_{n=1}^{\infty} \frac{1}{2\pi i} \left[\int_{B} e^{\bar{s}t}\hat{G}(s_n, \bar{s}) \,d\bar{s}\right] \times J_0(s_n \rho)J_{1}^{{-}2}(a' s_n)$$
where $\Phi ^{SI}(\rho, t)$ is the time-domain fluence in a semi-infinite medium [37]. $B$ denotes the Bromwich path where $\bar {s}$ is a complex number along the contour. We use a hyperbola contour as detailed previously [35,39]. Such a contour forces rapid decay in the integrand leading to more accurate and faster computation of the inverse time transform [35].

At this point, we need to relate the unnormalized electric field temporal autocorrelation function $G_1$ to the fluence rates previously derived. In this study, as well as most experimental settings, we are interested in the reflected light intensity at the top boundary ($z=0$). The diffuse reflectance is typically calculated as the current across the boundary or the flux [37]

$$R_f(\rho) ={-}D\nabla \Phi(\rho, z) \cdot (-\textbf{z})|_{z=0}$$
Alternatively, the reflectance can be written as a combination of both the fluence and flux terms [37] and for a medium of refractive index 1.4, the combined reflectance can be written as
$$R(\rho) = 0.118\Phi(\rho) + 0.306R_f(\rho)$$
The unnormalized temporal autocorrelation of the electric field function $G_1$ is simply calculated from the reflectance by replacing $\mu _a \rightarrow \mu _a^d$
$$G_1(\rho, \tau, D_b) = R(\rho, \mu_a^d)$$
The normalized electric field temporal autocorrelation function $g_1$ is then calculated by
$$g_1(\rho, \tau, D_b) = G_1(\rho, \tau) / G_1(\rho, \tau=0)$$
For solutions in the time-domain, the TD fluence given in Eq. (6) is used to calculate the reflectance in either Eq. (7) or Eq. (8), which are then translated to $G_1$ using Eq. (9).

2.2 Baseline flow tissue models using Monte Carlo simulations

To assess the accuracy and robustness of the described diffusion theory model to extract depth-dependent flow parameters, we consider a physiologically relevant 3-layer brain model consisting of the scalp, skull, and brain layer, as reported previously [9]. In this work, we only consider variations in BFi coefficients for scalp blood flow (SBF) and cerebral blood flow (CBF). All other parameters including the optical properties and layer thicknesses were held constant as shown in Fig. 1. The refractive index of each layer was also kept constant at 1.4 and the external medium’s refractive index equal to 1.0 to simulate air. Blood flow in the skull layer was assumed to be negligible (BFi = 0).

 figure: Fig. 1.

Fig. 1. Schematic of the optical properties and layer thicknesses of the 3-layer brain model used in the Monte Carlo simulations. The scalp and cerebral blood flow coefficients are varied while keeping all labeled properties constant.

Download Full Size | PDF

Using the optical properties and layer thicknesses shown in Fig. 1, the MC method was used to simulate a set of CW and TD normalized autocorrelation functions $g_1$. A previously developed MC code for photon propagation in laterally infinite layered media [40] was used to simulate the unnormalized autocorrelation function $G_1(\tau,\rho, t)$. The bin widths for $\Delta t$ and $\Delta \rho$ were $\pm$ 5 ps and 0.5 mm, respectively. 500 million photons were used for each simulation. For CW-DCS, $G_1(\tau, \rho )$ was obtained by summing $G_1(\tau, t)$ across the total time-of-flight simulated (from 0 ns to 7.5 ns) [14]. Following previous reports [9], we used 10 evenly spaced values for CBF $\in [20, 90]\, \mathrm {cm}^2/\mathrm {Gs}$ and 10 values of SBF $\in [2.5, 30]\, \mathrm {cm}^2/\mathrm {Gs}$ such that SBF $\in [1/8, 1/3]\times$ CBF. For the case when SBF $=30\, \mathrm {cm}^2/\mathrm {Gs}$, only 6 values of CBF from $20 - 58 \, \mathrm {cm}^2/\mathrm {Gs}$ were used. This resulted in 96 total combinations of tissue flow models that were stored for source-detector separations $\rho$ $\in [0.1, 3.05]$ in step sizes of 0.05 cm. The normalized autocorrelation functions $g_1(\tau )$ was obtained by dividing $G_1(\tau, t)$ by $G_1(0, t)$.

Simulated MC $g_1(\tau )$ were then fit to recover the SBF and CBF of the top-most and bottom-most layer using Eq. (10) with both reflectance models $R_f(\rho )$ given by Eq. (7) and $R(\rho )$ with Eq. (8). These fits were performed for varying $\rho$ in the CW domain and for different $t$ values in the TD. Each fit set the range of fitting for $\tau$ where $0.02 < g_1(\tau ) < 0.98$ as the fit window. Estimates of SBF and CBF were retrieved using a non-linear fitting tool in the Julia programming language [41] based on the Levenberg-Marquardt algorithm [42]. Forward solutions were computed using the numerical model (Eq. (7) and Eq. (8)). The optical properties and layer thicknesses of all the layers were considered known and were required as inputs to the inverse solver.

2.3 Noise model and sensitivity analysis

We briefly consider the sensitivity of $g_1(\tau )$ to a change in the analytical model’s respective BFi (e.g., SBF and CBF). The effect of noise on the sensitivity metric is considered using a previously accepted noise model [43,44]. We use a simplified version where the standard deviation of the noise in $g_1(\tau )$ is proportional to

$$\sigma(\tau) \approx \mathrm{exp}({\rho}) \cdot \sqrt{\frac{(1 + \mathrm{exp}(-\gamma \tau))}{T_{int}}}$$
where $T_{\mathrm {int}}$ is a combination of the correlator bin time interval and integration time (i.e., measurement duration), $\rho$ is the source-detector separation, and $\gamma$ is the decay rate of $g_1(\tau )$ [44]. The noise was added to $g_1(\tau )$ after multiplying the standard deviation with a random number between [-1, 1]. In Fig. 2, we show examples of the added noise on $g_1(\tau )$ for two source-detector separations, $\rho = 1$ and $3$ cm, with constant optical properties, $\mu _s' = 10 \, \mathrm {cm}^{-1}$ and $\mu _a = 0.1 \, \mathrm {cm}^{-1}$. We consider homogeneous flow values of $1$ and $20$ $\mathrm {cm}^2/\mathrm {Gs}$ in $g_1(\tau )$ and $g^*_1(\tau )$, respectively, for better visualization. Two noise models, $\sigma _1(\tau )$ and $\sigma _2(\tau )$, were calculated using Eq. (11) with $T_{\mathrm {int}}$ values of $1$ and $5$. $\gamma$ was calculated from the decay rate of $g_1(\tau )$.

 figure: Fig. 2.

Fig. 2. Example flow models at two source detector separations, $g_1(\tau, \rho =1$ cm) and $g^*_1(\tau, \rho =3$ cm), with two levels of added noise, $\sigma _1(\tau )$ and $\sigma _2(\tau )$, computed with Eq. (11) using $T_{\mathrm {int}} = 1$ and $5$, respectively.

Download Full Size | PDF

The sensitivity of $g_1(\tau, BFi)$ to a flow change, $\Delta$BFi, was estimated by calculating the relative change in $g_1(\tau )$, $(1 - g_1(\tau, BFi + \Delta BFi) / g_1(\tau, BFi))$, due to the perturbation. This allows for sensitivity estimates of CW- and TD-DCS to separate changes in SBF and CBF with noise. We note that the sensitivity can also be affected by $\tau, \rho, \mu _a,$ and $\mu _s'$. Here, we compare the sensitivity for $\rho \in (0.1, 3.5)$ cm with fixed optical properties, $\mu _s' = 10 \, \mathrm {cm}^{-1}$ and $\mu _a = 0.1 \, \mathrm {cm}^{-1}$. The dependence on $\tau$ was accounted for by averaging the relative change across a logarithm scaled array of $\tau$ values corresponding to $0.05<g_1(\tau )<0.95$.

3. Results

In Fig. 3, we compare the CW normalized autocorrelation function $g_1(\rho, \tau )$ simulated with the MC method and those calculated with diffusion using Eq. (7) and Eq. (8). Comparisons are shown for two different flow models at two source-detector separations of $0.4$ and $3$ cm. For Fig. 3(a), flow parameters of $6$ and $20$ $\mathrm {cm}^2/\mathrm {Gs}$ were used for the SBF and CBF, respectively whereas a SBF and CBF of $27$ and $20$ $\mathrm {cm}^2/\mathrm {Gs}$, respectively are shown in Fig. 3(b). Optical properties and layer thicknesses are as shown in Fig. 1.

 figure: Fig. 3.

Fig. 3. Comparison of the continuous-wave normalized autocorrelation function $g_1(\tau, \rho )$ simulated with the (symbols) Monte Carlo method and diffusion theory using the (blue solid line) $R(\rho )$ or (orange dash line) $R_f(\rho )$ reflectance models at two source-detector separations and at flow rates of (a) SBF = $6$ $\mathrm {cm}^2/\mathrm {Gs}$; CBF = $20$ $\mathrm {cm}^2/\mathrm {Gs}$ and (b) SBF = $27$ $\mathrm {cm}^2/\mathrm {Gs}$; CBF = $20$ $\mathrm {cm}^2/\mathrm {Gs}$. The model using $R(\rho )$ in Eq. (8) showed better agreement relative to $R_f(\rho )$ in Eq. (7) for the shorter source-detector separations and larger $\tau$.

Download Full Size | PDF

The largest difference between the two reflectance models can be observed at the shorter source-detector separations. Excellent agreement between MC and diffusion theory is observed when using Eq. (8), even at short distances of $\rho =4$ mm. A larger difference is observed when using just the flux term described by Eq. (7) at short distances, however this expression becomes more accurate farther away from the source. At the longer distance $\rho =3$ cm, both expressions show good agreement to the MC method. However, there is a small deviation for small $g_1$ values (<0.2) when using $R_f$ whereas $R$ is able to match the MC simulations over the full $\tau$ window.

In Fig. 4, we compare the solutions for the TD-DCS $g_1(\rho, \tau, t)$ when simulated with the MC method and calculated from diffusion theory using both reflectance models $R$ and $R_f$. Data in each figure are shown for a tissue model that has a SBF and CBF of $18$ and $43$ $\mathrm {cm}^2/\mathrm {Gs}$, respectively for two time values of 0.79 and 2.49 ns. The figures show these data at two source-detector separations, 0.2 and 2 cm (as noted in axis labels) in Fig. 4(a) and Fig. 4(b), respectively. In contrast to the CW domain, the choice in using Eq. (7) or Eq. (8) showed little difference between each other to calculate $g_1(\rho, \tau, t)$. Either calculation had very good agreement to the MC method showing only small deviations from each other when $g_1(\rho, \tau, t) < 0.1$. Larger differences between Monte Carlo and diffusion theory were observed for shorter source-detector distance and longer time gates (Fig. 4(a), $t = 2.49$ ns). However, it is also useful to note that MC simulations exhibit higher noise for longer times which could be causing these differences.

 figure: Fig. 4.

Fig. 4. Comparison of the time-domain normalized autocorrelation function $g_1(\tau, \rho, t)$ simulated with the (symbols) Monte Carlo method and diffusion theory using the (blue solid line) $R(\rho )$ or (orange dash line) $R_f(\rho )$ reflectance models at two time values and at source-detector separations of (a) $\rho = 0.2$ cm and (b) $\rho = 2$ cm. A single flow model is considered where SBF = $18$ $\mathrm {cm}^2/\mathrm {Gs}$ and CBF = $43$ $\mathrm {cm}^2/\mathrm {Gs}$. Both reflectance models showed similar accuracy.

Download Full Size | PDF

In Fig. 5, the results of the inverse fitting procedure for the CW domain are shown. Here, the diffusion model is used as a forward model to fit for both SBF and CBF using the MC method as baseline simulations. In Fig. 5(a) and Fig. 5(b), the recovered SBF and CBF are shown, respectively for all 96 of the flow tissue models recovered at two source-detector separations of 1 and 2.5 cm. For these figures we show just the fits using the $R_f$ reflectance model as only small differences between the two reflectance models are observed at these two source-detector separations. There was a larger difference from true flow coefficients as determined by the MC simulations (black dashed lines) and recovered coefficients using the diffusion model when SBF was larger, however this did not appear to affect the accuracy of recovered CBF. A larger $\rho = 2.5$ cm was able to better recover CBF, though, both $\rho = 1$ and $2.5$ cm showed sensitivity to SBF and CBF.

 figure: Fig. 5.

Fig. 5. Results of inverse fits using CW-DCS diffusion theory fitted to MC simulated $g_1(\tau, \rho )$ curves. The recovered (a) SBF and (b) CBF coefficients are shown for two values of $\rho$ where the black dashed lines represent the true flow values used in the MC simulations, for each tissue model, when recovered with the $R_f$ reflectance model. The percent error between recovered absolute values of (c) SBF and (d) CBF is averaged over all the tissue models at four source-detector separations using both reflectance expressions (Eq. (7) and Eq. (8)). The average error in relative flow changes of (e) SBF and (f) CBF are shown when recovered values of the first tissue model (as shown in (a) and (b)) are normalized to one.

Download Full Size | PDF

In Fig. 5(c) and Fig. 5(d), the average percent error $100\times ((\mathrm {BFi}_{MC} - \mathrm {BFi}_{DT})/\mathrm {BFi}_{MC})$ between the recovered SBF and CBF, is shown averaged across all 96 tissue models for four source-detector separations when using the two reflectance models. As shown in Fig. 3, at short distances ($\rho = 0.5$ cm), $R(\rho )$ produces significantly more accurate results. Interestingly, the reflectance model using just the flux $R_f(\rho )$ produced slightly more accurate results at the longer distances ($\rho > 1$ cm). The error in recovered SBF also slightly increases at longer distances as the sensitivity to the top layer decreases as we move farther away from the source. On the other hand, the recovered CBF significantly improves as we move to larger separations (Fig. 5(d)).

The percent error in recovery of relative SBF and CBF are shown in Fig. 5(e) and Fig. 5(f), respectively for the two reflectance models at five source-detector separations. The relative error is calculated by normalizing the values of the true and recovered BFi of the first tissue model shown in Fig. 5(a) and Fig. 5(b) (SBF=2.5 $\mathrm {cm}^2/\mathrm {Gs}$; CBF = 20 $\mathrm {cm}^2/\mathrm {Gs}$), and averaging the percent error across all the tissue models. The choice to normalize to the first tissue model is rather arbitrary, however the best absolute accuracy was observed for low SBF so this resulted in shifting absolute values the least. Good recovery of SBF is observed even at very short distances ($\rho = 2.5$ mm) while error increases after 1 cm due to the decrease in sensitivity to SBF. In this case, using $R(\rho )$ produced more accurate estimates of SBF. For CBF, large errors are observed at short distances ($\rho = 5$ mm) due to the very limited sensitivity to deeper layers at short separations. This error rapidly improves to less than 5% for $\rho > 1$ cm for both reflectance models.

In Fig. 6, the results of the inverse fitting procedure for the TD are shown. Similar to the CW domain, the inverse problem fits for two parameters, SBF and CBF. However, we show results using only the flux reflectance model, $R_f$, as the two reflectance models produced errors within 3% of each other for all simulations. Therefore, we focus on the recovery of the flow coefficients at time values of 0.5, 1.5, and 2.5 ns. In Fig. 6(a) and Fig. 6(b), we report results for all 96 tissue models at a single source-detector separation of $\rho =1$ cm for the three time values to recover SBF and CBF, respectively. Similar to the CW results, the recovery of SBF depended on the absolute SBF value as larger SBF values produced larger errors particularly for late time values. However, recovery of CBF was not affected by higher errors in recovery of SBF. As expected, shorter times were able to recover SBF better, while later times were better able to recover CBF. Rather surprisingly, the short time ($t=0.5$ ns) also showed large sensitivity to CBF.

 figure: Fig. 6.

Fig. 6. Results of the inverse fitting process when fitting the time-domain diffusion theory model to simulated $g_1(\tau, \rho, t)$ with the MC method. The recovered (a) SBF and (b) CBF coefficients are shown for three values of $t$ when $\rho =1$ cm where the black dashed lines represent the true flow values used in the Monte Carlo simulations. The percent error between recovered absolute values of (c) SBF and (d) CBF is averaged over all the tissue models at three source-detector separations for three $t$ values. We show just the error using the $R_f(\rho, t)$ reflectance model as it was similar to $R(\rho, t)$. Averaged relative error normalized to the first tissue model is shown for values of (e) SBF and (f) CBF.

Download Full Size | PDF

The resulting percent error shown in Fig. 6(a) and Fig. 6(b) between the true MC values and recovered values using diffusion theory were then averaged across all tissue models for source-detector separation of $\rho = 0.5, 1$ and $2$ cm. The average error in recovered SBF and CBF for the three time values are shown in Fig. 6(c) and Fig. 6(d), respectively. The late time value, $t=2.5$ ns, showed the highest and most volatile error in recovered SBF. This can be explained by the inconsistent recovery at higher absolute values of SBF as shown in Fig. 6(a). The error in recovered SBF using $t=0.5$ and 1.5 ns was less than 10% in all the source-detector separations shown.

The recovered CBF was less than 4% at all time values when $\rho =0.5$ cm, however for $t=0.5$ ns the error increased at larger source-detector separations to over 10% when $\rho =2$ cm. In Fig. 6(e) and Fig. 6(f) we show the percent error in relative changes in SBF and CBF, respectively when normalizing the starting values to the first tissue model. In contrast to the CW results, the relative errors were actually worse than absolute recovery due to the much noisier recovery of SBF at faster flow rates. Though, recovery of CBF using $t=1.5$ and 2.5 ns was less than 2%.

To better understand the results presented in Figs. 5 and 6, in Fig. 7 we show an example flow model with baseline SBF and CBF of $20$ and $60$ $\mathrm {cm}/ \mathrm {Gs}$, respectively in the (Fig. 7(a) and Fig. 7(b)) CW domain at two source-detector separations ($\rho =0.5$ and $4.5$ cm) and in the (Fig. 7(c) and Fig. 7(d)) TD at two delay times ($t=0.5$ and $2.5$ ns) at a single source-detector separation ($\rho =1$ cm). The baseline optical properties are as shown in Fig. 1. We next increase either SBF or CBF by three times to show the sensitivity of $g_1(\tau )$ on the two flow parameters. In the CW domain, increasing SBF by 3 times has an effect at both short ($\rho =0.5$ cm) and long ($\rho =4.5$ cm) source-detector separations. However, the effect is more uniform at shorter distances while longer distances are more affected at later delay times (Fig. 7(a)). On the other hand, the short distances show little effect to changes in CBF (Fig. 7(b)) except at shorter delay times. This effect is more pronounced at the longer source-detector distance where more sensitivity to CBF is observed. In the TD, the shorter delay time ($t=0.5$ ns) shows a large sensitivity to SBF changes while showing less sensitivity to CBF. The later time ($t=2.5$ ns) shows a much larger sensitivity to CBF while only showing a minimal though non-negligible sensitivity to SBF.

 figure: Fig. 7.

Fig. 7. The sensitivity of the CW $g_1(\tau )$ to an increase of either (a) SBF or (b) CBF by three times at two source-detector separations. A similar example is shown in (c) and (d) for the TD $g_1(\tau )$ at a single source-detector separation ($\rho =1$ cm) and two time values for SBF and CBF, respectively.

Download Full Size | PDF

In Fig. 8, we quantitatively compared the sensitivity in CW- and TD-DCS from a relative change in $g_1(\tau )$ due to a doubling in the respective flow coefficient of interest (e.g., SBF and CBF) with other variables constant as a function of source-detector separation $\rho \in [0.1, 3.5]$ cm. The baseline optical properties and flow parameters are as described in Fig. 1 and Fig. 7. To account for variable sensitivities to SBF and CBF as a function of $\tau$, the relative change was averaged in the window $0.05<g_1(\tau )<0.95$. In Fig. 8(a) and Fig. 8(b) we show the sensitivity of CW-DCS to an increase in SBF and CBF, respectively. In addition to the direct calculation of the forward model, we compute the sensitivity metric after adding noise using Eq. (11) with two measurement collection times, $T_{\mathrm {int}}=1$ and 5 labeled as $\sigma _1(\tau )$ and $\sigma _2(\tau )$, respectively. The error bars correspond to the standard deviation after averaging over 50 samples of random noise generation. In Fig. 8(c) and Fig. 8(d) the same is shown considering the sensitivity of TD-DCS to an increase in SBF and CBF at two photon delay times of $t= 0.5$ and $2.5$ ns as labeled directly for each group of noise models. We note the black lines represent the sensitivity directly from the analytical model without any added noise and do not have error bars.

 figure: Fig. 8.

Fig. 8. The average relative change of $g_1(\tau )$ from increasing the (left column) SBF and (right column) CBF by two times in the (top row) CW and (bottom row) TD at two photon arrival times ($t = 0.5$ and 2.5 ns) as a function of source-detector separation. Two noise models, $\sigma _1(\tau )$ and $\sigma _2(\tau )$, were considered using a collection time, $T_{\mathrm {int}}$, of 1 and 5 seconds, respectively, as computed with Eq. (11). The CW domain shows a varying sensitivity to both SBF and CBF as a function of source-detector separation whereas the TD shows a relatively constant sensitivity.

Download Full Size | PDF

For the SBF sensitivity shown in Fig. 8(a), we can see that the CW model shows the highest sensitivity to SBF at shorter distances while steadily decreasing at longer source-detector separations. However, the TD shows a relatively constant sensitivity to SBF over the whole range of $\rho$. In fact, the sensitivity slightly increases as we increase the source-detector distance. This can be explained by the fixed photon path length within the medium and that to maintain the same path length in the medium at longer distances from the source the average photon has to travel more in the shallower layers. On the other hand, the CW model shows an increasing sensitivity to CBF as we increase the source-detector separation (Fig. 8(d)). However, the sensitivity appears to have only logarithmic growth. As in the SBF case, the TD CBF sensitivity is relatively independent to the source-detector separation and only decreases slightly at larger distances. It was also worthwhile to observe that the TD sensitivity seems to saturate quickly (for $t>2.5$ ns), whereby increasing the time gates only marginally increased CBF sensitivity. At larger source-detector separations the added noise model significantly affects the expected sensitivity to both SBF and CBF.

4. Discussion

Layered analytical models have been employed with success in CW-DCS applications to recover flow coefficients [9,20,23,27] but are not typically used for TD-DCS analysis due to their increased complexity and computational cost [32]. We have validated a layered diffusion model based on our previous work [35] that allows for solving the inverse problem in both CW- and TD-DCS. The main advantages of this approach is its better numerical stability, ability to provide error tolerances, computational speed, and ability to provide fast gradient information [35]. This allows for forward simulation of the CW and TD autocorrelation function in $\approx 20$ and $\approx 150$ microseconds, respectively. Typically, the inverse problem for the fits shown here took $<10$ and $<250$ milliseconds to fit for both SBF and CBF in the CW and TD, respectively. Although this is orders of magnitude faster than both MC simulations and previously reported layered models [7,33,34], using homogeneous models that only fit a single parameter are about 10x and 60x faster in the CW and TD, respectively. These reported times were all done using a single core on a personal laptop. The computational time is also linearly dependent on the number of $\tau$ values considered in the inverse problem and should be kept as sparse as possible for faster data processing. Additionally, we only considered fitting two parameters (SBF and CBF) while other coefficients were considered known. The effect of optical property errors have been previously investigated [9,20,36]. However, the forward model simultaneously computes both the TD reflectance along with $g_1(\tau )$, allowing it be used with TD near-infrared spectroscopy measurements to fit for optical properties along with DCS measurements to fit for flow coefficients [45,46] with similar forward computational times.

A disadvantage of the presented model are the limitations given by the validity of diffusion theory [7]. The Monte Carlo method has become more popular as computing resources are more readily available to overcome some of the concerns of diffusion theory at short source-detector separations [7,28]. However, we have found that the diffusion model can accurately recover absolute flow coefficients at very short separations ($\rho = 4$ mm) in the CW domain while being able to recover relative errors at even shorter separations, which is inline with previous research [40]. We investigated two equations that seek to model the reflectance as just the flux or a combination of the fluence and flux [37]. We have found that for CW measurements, using both the fluence and flux (Eq. (8)) provides more accurate results at shorter source-detector separations and smaller values of $g_1(\rho, \tau )$. This translated to better recovered error (Fig. 5 c and d) at shorter separations, however this did not provide better error for source-detector separations longer than 1 cm. On the other hand, we did not find any significant differences between the two reflectance models when computing the TD autocorrelation. At short source-detector separations ($\approx 2$ mm), we found the time-domain model could be used more freely to model $g_1(\tau, \rho, t)$ when selecting for later arriving photons $t>0.5$ ns (Fig. 4).

One of the primary reasons to use a layered model is the ability to recover and separate both extracerebral and cerebral flow changes. In the CW domain, the diffusion theory model was able to recover SBF within 10% of the Monte Carlo simulations for source-detector separations between 0.4 and 3 cm (Fig. 5(c)). The error slightly increases as we increase in source-detector separation which can be explained by the decrease in sensitivity to SBF at longer distances (Fig. 8(a)). On the other hand, increasing the source-detector separation rapidly improved the estimates of CBF which was recovered within 10% for $\rho >2$ cm. Perhaps most interesting is the rate at which the error improves (Fig. 5(d)) in recovered CBF as we increase $\rho$, which could be explained by the rate at which the sensitivity to CBF increases. In Fig. 8(b), we see an almost logarithmic growth in sensitivity in the CW domain where a rapid rise in sensitivity is seen as we increase $\rho$ up to 2 cm with just a smaller marginal increase after. We briefly considered a simplified noise model applied to both the CW and TD in Fig. 8. Although we showed steady improvements in recovered CBF at longer distances with no noise, adding noise can significantly affect sensitivities as we increase source-detector separation and decrease the measurement duration. Therefore, the marginal increase in recovered accuracy at longer source-detector separations and the decrease in expected sensitivity due to noise suggests that there is an optimal distance for measurement. In other words, the optimal distance will be after a large gain in sensitivity but short enough to maintain adequate signal to noise.

TD-DCS provides potential for better sensitivity to deeper tissues, and therefore CBF, by selecting for later arriving photons [14]. When considering just $g_1(\tau, \rho, t)$, we also demonstrated the larger sensitivity to CBF when compared to CW $g_1(\tau, \rho )$. Figure 6(d) shows that the TD is able to recover CBF within less than 5% error between source-detector separations of 0.5 and 2 cm. A primary advantage is even using relatively short time windows of 0.5 and 1.5 ns contrast to CBF is still high. As shown in Fig. 8(b), using very late time windows doesn’t substantially increase the CBF sensitivity. This is experimentally important as using shorter time windows improves signal to noise contrast, thus using moderate time gates (1-2 ns) could provide similar sensitivities to very late gates while providing more contrast. Additionally, we found that the sensitivity and CBF recovery accuracy is mostly independent of the source-detector separation. We do note that the data and sensitivities reported here are specific to the tissue geometry considered. Therefore, having a larger extracerebral layer or higher scattering coefficients could influence these results. Finally, making substantial claims on the ability of TD measurements to better quantify CBF compared to CW is difficult [47]. In Fig. 8(b), we show a theoretical comparison based on $g_1(\tau )$ which does show that TD measurements have a much larger sensitivity to CBF when $\rho <2.5$ cm. However, experimentally we measure $g_2(\tau )$ which can be related to $g_1(\tau )$ by the Siegert relation [3]. The effect of $\beta$ in the Siegert relation can dramatically affect the sensitivities of each domain to CBF that we calculate in Fig. 8(b) and Fig. 8(d) [47,48]. Therefore, comparing the CW and TD for experimental measurements must consider the effects of the instrument response function, coherence length, time gate, and noise. The simplified noise models we discuss in Fig. 8 must also consider the effect of photon arrival time and gate width which we did not consider. Future studies will focus on more rigorous noise models and their effect within the inverse problem as well as incorporating instrument response functions and finite coherence to properly compare absolute sensitivity of CW and TD measurements.

We also found that the diffusion model in both the CW and TD showed a slight offset in recovered SBF as the absolute value of SBF approached that of CBF. However, we did not find that this had any effect on recovered values of CBF. This could be explained by a more difficult inverse problem when the two flow coefficients have similar values, however we did not observe any crosstalk in the two coefficients as CBF was largely unaffected by these increased errors in SBF. We also note the importance of absolute SBF values in sensitivity calculations. As SBF increases, the sensitivity of both CW and TD $g_1(\tau )$ become more sensitive to changes in SBF. We show a typical physiological flow rate of SBF and CBF in Fig. 7 where SBF is a third of the speed of CBF. In Fig. 7 we can see that every $g_1(\tau )$ has some sensitivity to SBF even in the CW domain at long source-detector separations ($\rho = 4.5$ cm) and in the time-domain at long delay times ($t=2.5$ ns). This effect is more pronounced for higher SBF values, though if SBF was much less than a third of CBF we would see less of a sensitivity to SBF. Although this fact is less important when using layered models as they can more robustly model these changes, using homogeneous models will show variable sensitivities to both SBF and CBF depending on their absolute flow rates. These findings suggest that more complicated models might be necessary even when selecting for photons of longer pathlengths and larger source-detector separations.

5. Conclusion

We developed and verified a diffusion based layered analytical model that can be used to recover estimates of both scalp and cerebral blood flow in real-time when using either continuous-wave or time-domain diffusion correlation spectroscopy. Utilizing a reflectance model that uses both the fluence and flux allows for absolute recovery of flow parameters at source-detector separations as short as 4 mm whereas relative recovery was demonstrated down to separations of 2 mm for continuous-wave measurements. The verified model was able to simultaneously recover SBF and CBF within 5% of Monte Carlo simulations in under 10 ms and 250 ms in the continuous-wave and time-domains, respectively allowing its use for real-time data analysis. All analysis tools are openly available online [49].

Acknowledgments

The authors are greatly appreciative of Dr. André Liemert for helpful discussions on developing approximate forms of the theoretical solutions.

Disclosures

The authors declare no conflicts of interest.

Data Availability

No data were generated or analyzed in the presented research. The layered analytical model presented here is freely available online for use at [49].

References

1. T. Durduran and A. G. Yodh, “Diffuse correlation spectroscopy for non-invasive, micro-vascular cerebral blood flow measurement,” NeuroImage 85, 51–63 (2014). [CrossRef]  

2. E. M. Buckley, A. B. Parthasarathy, P. E. Grant, A. G. Yodh, and M. A. Franceschini, “Diffuse correlation spectroscopy for measurement of cerebral blood flow: future prospects,” Neurophotonics 1(1), 011009 (2014). [CrossRef]  

3. T. Durduran, R. Choe, W. B. Baker, and A. G. Yodh, “Diffuse optics for tissue monitoring and tomography,” Rep. Prog. Phys. 73(7), 076701 (2010). [CrossRef]  

4. S. Y. Lee, J. M. Pakela, M. C. Helton, K. Vishwanath, Y. G. Chung, N. J. Kolodziejski, C. J. Stapels, D. R. McAdams, D. E. Fernandez, J. F. Christian, J. O’Reilly, D. Farkas, B. B. Ward, S. E. Feinberg, and M.-A. Mycek, “Compact dual-mode diffuse optical system for blood perfusion monitoring in a porcine model of microvascular tissue flaps,” J. Biomed. Opt. 22(12), 1 (2017). [CrossRef]  

5. R. H. Wilson, K. Vishwanath, and M.-A. Mycek, “Optical methods for quantitative and label-free sensing in living human tissues: principles, techniques, and applications,” Adv. Phys. 1(4), 523–543 (2016). [CrossRef]  

6. J. J. Selb, D. A. Boas, S.-T. Chan, K. C. Evans, E. M. Buckley, and S. A. Carp, “Sensitivity of near-infrared spectroscopy and diffuse correlation spectroscopy to brain hemodynamics: simulations and experimental findings during hypercapnia,” Neurophotonics 1(1), 015005 (2014). [CrossRef]  

7. M. M. Wu, S.-T. Chan, D. Mazumder, D. Tamborini, K. A. Stephens, B. Deng, P. Farzam, J. Y. Chu, M. A. Franceschini, J. Z. Qu, and S. A. Carp, “Improved accuracy of cerebral blood flow quantification in the presence of systemic physiology cross-talk using multi-layer monte carlo modeling,” Neurophotonics 8(01), 015001 (2021). [CrossRef]  

8. D. A. Boas, L. Campbell, and A. G. Yodh, “Scattering and imaging with diffusing temporal field correlations,” Phys. Rev. Lett. 75(9), 1855–1858 (1995). [CrossRef]  

9. H. Zhao, E. Sathialingam, and E. M. Buckley, “Accuracy of diffuse correlation spectroscopy measurements of cerebral blood flow when using a three-layer analytical model,” Biomed. Opt. Express 12(11), 7149–7161 (2021). [CrossRef]  

10. D. A. Boas, S. Sakadžić, J. J. Selb, P. Farzam, M. A. Franceschini, and S. A. Carp, “Establishing the diffuse correlation spectroscopy signal relationship with blood flow,” Neurophotonics 3(3), 031412 (2016). [CrossRef]  

11. S. Carp, G. Dai, D. A. Boas, M. A. Franceschini, and Y. Kim, “Validation of diffuse correlation spectroscopy measurements of rodent cerebral blood flow with simultaneous arterial spin labeling mri; towards mri-optical continuous cerebral metabolic monitoring,” Biomed. Opt. Express 1(2), 553–565 (2010). [CrossRef]  

12. R. B. Saager and A. J. Berger, “Measurement of layer-like hemodynamic trends in scalp and cortex: implications for physiological baseline suppression in functional near-infrared spectroscopy,” J. Biomed. Opt. 13(3), 034017 (2008). [CrossRef]  

13. J. Selb, K.-C. Wu, P.-Y. Ivy Lin, J. Sutin, P. Farzam, S. Bechek, A. Shenoy, A. B. Patel, D. A. Boas, M. A. Franceschini, and E. S. Rosenthal, “Prolonged monitoring of cerebral blood flow and autoregulation with diffuse correlation spectroscopy in neurocritical care patients,” Neurophotonics 5(04), 1 (2018). [CrossRef]  

14. J. Sutin, B. Zimmerman, D. Tyulmankov, D. Tamborini, K. C. Wu, J. Selb, A. Gulinatti, I. Rech, A. Tosi, D. A. Boas, and M. A. Franceschini, “Time-domain diffuse correlation spectroscopy,” Optica 3(9), 1006–1013 (2016). [CrossRef]  

15. M. Pagliazzi, S. K. V. Sekar, L. Colombo, E. Martinenghi, J. Minnema, R. Erdmann, D. Contini, A. Dalla Mora, A. Torricelli, A. Pifferi, and T. Durduran, “Time domain diffuse correlation spectroscopy with a high coherence pulsed source: in vivo and phantom results,” Biomed. Opt. Express 8(11), 5311–5325 (2017). [CrossRef]  

16. L. Colombo, M. Pagliazzi, S. K. V. Sekar, D. Contini, A. Dalla Mora, L. Spinelli, A. Torricelli, T. Durduran, and A. Pifferi, “Effects of the instrument response function and the gate width in time-domain diffuse correlation spectroscopy: model and validations,” Neurophotonics 6(03), 1 (2019). [CrossRef]  

17. D. Mazumder, M. M. Wu, N. Ozana, D. Tamborini, M. A. Franceschini, and S. A. Carp, “Optimization of time domain diffuse correlation spectroscopy parameters for measuring brain blood flow,” Neurophotonics 8(03), 035005 (2021). [CrossRef]  

18. D. Tamborini, K. A. Stephens, M. M. Wu, P. Farzam, A. M. Siegel, O. Shatrovoy, M. Blackwell, D. A. Boas, S. A. Carp, and M. A. Franceschini, “Portable system for time-domain diffuse correlation spectroscopy,” IEEE Trans. Biomed. Eng. 66(11), 3014–3025 (2019). [CrossRef]  

19. S. Samaei, P. Sawosz, M. Kacprzak, Ż. Pastuszak, D. Borycki, and A. Liebert, “Time-domain diffuse correlation spectroscopy (td-dcs) for noninvasive, depth-dependent blood flow quantification in human tissue in vivo,” Sci. Rep. 11(1), 1817 (2021). [CrossRef]  

20. L. Gagnon, M. Desjardins, J. Jehanne-Lacasse, L. Bherer, and F. Lesage, “Investigation of diffuse correlation spectroscopy in multi-layered media including the human head,” Opt. Express 16(20), 15514–15530 (2008). [CrossRef]  

21. A. Duncan, J. H. Meek, M. Clemence, C. E. Elwell, L. Tyszczuk, M. Cope, and D. Delpy, “Optical pathlength measurements on adult head, calf and forearm and the head of the newborn infant using phase resolved optical spectroscopy,” Phys. Med. Biol. 40(2), 295–304 (1995). [CrossRef]  

22. T. Durduran, G. Yu, M. G. Burnett, J. A. Detre, J. H. Greenberg, J. Wang, C. Zhou, and A. G. Yodh, “Diffuse optical measurement of blood flow, blood oxygenation, and metabolism in a human brain during sensorimotor cortex activation,” Opt. Lett. 29(15), 1766–1768 (2004). [CrossRef]  

23. K. Verdecchia, M. Diop, A. Lee, L. B. Morrison, T.-Y. Lee, and K. S. Lawrence, “Assessment of a multi-layered diffuse correlation spectroscopy method for monitoring cerebral blood flow in adults,” Biomed. Opt. Express 7(9), 3659–3674 (2016). [CrossRef]  

24. W. B. Baker, A. B. Parthasarathy, T. S. Ko, D. R. Busch, K. Abramson, S.-Y. Tzeng, R. C. Mesquita, T. Durduran, J. H. Greenberg, D. K. Kung, and A. G. Yodh, “Pressure modulation algorithm to separate cerebral hemodynamic signals from extracerebral artifacts,” Neurophotonics 2(3), 035004 (2015). [CrossRef]  

25. R. C. Mesquita, S. S. Schenkel, D. L. Minkoff, X. Lu, C. G. Favilla, P. M. Vora, D. R. Busch, M. Chandra, J. H. Greenberg, J. A. Detre, and A. G. Yodh, “Influence of probe pressure on the diffuse correlation spectroscopy blood flow signal: extra-cerebral contributions,” Biomed. Opt. Express 4(7), 978–994 (2013). [CrossRef]  

26. F. Jaillon, S. E. Skipetrov, J. Li, G. Dietsche, G. Maret, and T. Gisler, “Diffusing-wave spectroscopy from head-like tissue phantoms: influence of a non-scattering layer,” Opt. Express 14(22), 10181–10194 (2006). [CrossRef]  

27. J. Li, G. Dietsche, D. Iftime, S. E. Skipetrov, G. Maret, T. Elbert, B. Rockstroh, and T. Gisler, “Noninvasive detection of functional brain activity with near-infrared diffusing-wave spectroscopy,” J. Biomed. Opt. 10(4), 044002 (2005). [CrossRef]  

28. M. M. Wu, K. Perdue, S.-T. Chan, K. A. Stephens, B. Deng, M. A. Franceschini, and S. A. Carp, “Complete head cerebral sensitivity mapping for diffuse correlation spectroscopy using subject-specific magnetic resonance imaging models,” Biomed. Opt. Express 13(3), 1131–1151 (2022). [CrossRef]  

29. C. Zhu and Q. Liu, “Review of monte carlo modeling of light transport in tissues,” J. Biomed. Opt. 18(5), 050902 (2013). [CrossRef]  

30. Q. Fang and D. A. Boas, “Monte carlo simulation of photon migration in 3d turbid media accelerated by graphics processing units,” Opt. Express 17(22), 20178–20190 (2009). [CrossRef]  

31. J. J. Selb, T. M. Ogden, J. Dubb, Q. Fang, and D. A. Boas, “Comparison of a layered slab and an atlas head model for monte carlo fitting of time-domain near-infrared spectroscopy data of the adult head,” J. Biomed. Opt. 19(1), 016010 (2014). [CrossRef]  

32. J. Li, L. Qiu, C.-S. Poon, and U. Sunar, “Analytical models for time-domain diffuse correlation spectroscopy for multi-layer and heterogeneous turbid media,” Biomed. Opt. Express 8(12), 5518–5532 (2017). [CrossRef]  

33. S. Geiger, D. Reitzle, A. Liemert, and A. Kienle, “Determination of the optical properties of three-layered turbid media in the time domain using the p 3 approximation,” OSA Continuum 2(6), 1889–1899 (2019). [CrossRef]  

34. A. Liemert and A. Kienle, “Light diffusion in a turbid cylinder. ii. layered case,” Opt. Express 18(9), 9266–9279 (2010). [CrossRef]  

35. M. Helton, S. Zerafa, K. Vishwanath, and M.-A. Mycek, “Efficient computation of the steady-state and time-domain solutions of the photon diffusion equation in layered turbid media,” Sci. Rep. 12(1), 18979 (2022). [CrossRef]  

36. D. Irwin, L. Dong, Y. Shang, R. Cheng, M. Kudrimoti, S. D. Stevens, and G. Yu, “Influences of tissue absorption and scattering on diffuse correlation spectroscopy blood flow measurements,” Biomed. Opt. Express 2(7), 1969–1985 (2011). [CrossRef]  

37. A. Kienle and M. S. Patterson, “Improved solutions of the steady-state and the time-resolved diffusion equations for reflectance from a semi-infinite turbid medium,” J. Opt. Soc. Am. A 14(1), 246–254 (1997). [CrossRef]  

38. D. Contini, F. Martelli, and G. Zaccanti, “Photon migration through a turbid slab described by a model based on diffusion approximation. i. theory,” Appl. Opt. 36(19), 4587–4599 (1997). [CrossRef]  

39. A. Liemert and A. Kienle, “Application of the laplace transform in time-domain optical spectroscopy and imaging,” J. Biomed. Opt. 20(11), 110502 (2015). [CrossRef]  

40. K. Vishwanath and S. Zanfardino, “Diffuse correlation spectroscopy at short source-detector separations: Simulations, experiments and theoretical modeling,” Appl. Sci. 9(15), 3047 (2019). [CrossRef]  

41. J. Bezanson, A. Edelman, S. Karpinski, and V. B. Shah, “Julia: A fresh approach to numerical computing,” SIAM Rev. 59(1), 65–98 (2017). [CrossRef]  

42. P. K. Mogensen and A. N. Riseth, “Optim: A mathematical optimization package for Julia,” J. Open Source Softw. 3(24), 615 (2018). [CrossRef]  

43. C. Zhou, G. Yu, D. Furuya, J. H. Greenberg, A. G. Yodh, and T. Durduran, “Diffuse optical correlation tomography of cerebral blood flow during cortical spreading depression in rat brain,” Opt. Express 14(3), 1125–1144 (2006). [CrossRef]  

44. L. Cortese, G. L. Presti, M. Pagliazzi, et al., “Recipes for diffuse correlation spectroscopy instrument design using commonly utilized hardware based on targets for signal-to-noise ratio and precision,” Biomed. Opt. Express 12(6), 3265–3281 (2021). [CrossRef]  

45. M. Diop, K. Verdecchia, T.-Y. Lee, and K. St Lawrence, “Calibration of diffuse correlation spectroscopy with a time-resolved near-infrared technique to yield absolute cerebral blood flow measurements,” Biomed. Opt. Express 2(7), 2068–2081 (2011). [CrossRef]  

46. W. B. Baker, R. Balu, L. He, V. C. Kavuri, D. R. Busch, O. Amendolia, F. Quattrone, S. Frangos, E. Maloney-Wilensky, K. Abramson, E. M. Gabrielli, A. G. Yodh, and W. A. Kofke, “Continuous non-invasive optical monitoring of cerebral blood flow and oxidative metabolism after acute brain injury,” J. Cereb. Blood Flow Metab. 39(8), 1469–1485 (2019). [CrossRef]  

47. X. Cheng, H. Chen, E. J. Sie, F. Marsili, and D. A. Boas, “Development of a monte carlo-wave model to simulate time domain diffuse correlation spectroscopy measurements from first principles,” J. Biomed. Opt. 27(08), 083009 (2022). [CrossRef]  

48. X. Cheng, D. Tamborini, S. A. Carp, O. Shatrovoy, B. Zimmerman, D. Tyulmankov, A. Siegel, M. Blackwell, M. A. Franceschini, and D. A. Boas, “Time domain diffuse correlation spectroscopy: modeling the effects of laser coherence length and instrument response function,” Opt. Lett. 43(12), 2756–2759 (2018). [CrossRef]  

49. M. Helton, “Lightpropagation.jl,” Github, 2021, https://github.com/heltonmc/LightPropagation.jl.

Data Availability

No data were generated or analyzed in the presented research. The layered analytical model presented here is freely available online for use at [49].

49. M. Helton, “Lightpropagation.jl,” Github, 2021, https://github.com/heltonmc/LightPropagation.jl.

Cited By

Optica participates in Crossref's Cited-By Linking service. Citing articles from Optica Publishing Group journals and other participating publishers are listed here.

Alert me when this article is cited.


Figures (8)

Fig. 1.
Fig. 1. Schematic of the optical properties and layer thicknesses of the 3-layer brain model used in the Monte Carlo simulations. The scalp and cerebral blood flow coefficients are varied while keeping all labeled properties constant.
Fig. 2.
Fig. 2. Example flow models at two source detector separations, $g_1(\tau, \rho =1$ cm) and $g^*_1(\tau, \rho =3$ cm), with two levels of added noise, $\sigma _1(\tau )$ and $\sigma _2(\tau )$, computed with Eq. (11) using $T_{\mathrm {int}} = 1$ and $5$, respectively.
Fig. 3.
Fig. 3. Comparison of the continuous-wave normalized autocorrelation function $g_1(\tau, \rho )$ simulated with the (symbols) Monte Carlo method and diffusion theory using the (blue solid line) $R(\rho )$ or (orange dash line) $R_f(\rho )$ reflectance models at two source-detector separations and at flow rates of (a) SBF = $6$ $\mathrm {cm}^2/\mathrm {Gs}$; CBF = $20$ $\mathrm {cm}^2/\mathrm {Gs}$ and (b) SBF = $27$ $\mathrm {cm}^2/\mathrm {Gs}$; CBF = $20$ $\mathrm {cm}^2/\mathrm {Gs}$. The model using $R(\rho )$ in Eq. (8) showed better agreement relative to $R_f(\rho )$ in Eq. (7) for the shorter source-detector separations and larger $\tau$.
Fig. 4.
Fig. 4. Comparison of the time-domain normalized autocorrelation function $g_1(\tau, \rho, t)$ simulated with the (symbols) Monte Carlo method and diffusion theory using the (blue solid line) $R(\rho )$ or (orange dash line) $R_f(\rho )$ reflectance models at two time values and at source-detector separations of (a) $\rho = 0.2$ cm and (b) $\rho = 2$ cm. A single flow model is considered where SBF = $18$ $\mathrm {cm}^2/\mathrm {Gs}$ and CBF = $43$ $\mathrm {cm}^2/\mathrm {Gs}$. Both reflectance models showed similar accuracy.
Fig. 5.
Fig. 5. Results of inverse fits using CW-DCS diffusion theory fitted to MC simulated $g_1(\tau, \rho )$ curves. The recovered (a) SBF and (b) CBF coefficients are shown for two values of $\rho$ where the black dashed lines represent the true flow values used in the MC simulations, for each tissue model, when recovered with the $R_f$ reflectance model. The percent error between recovered absolute values of (c) SBF and (d) CBF is averaged over all the tissue models at four source-detector separations using both reflectance expressions (Eq. (7) and Eq. (8)). The average error in relative flow changes of (e) SBF and (f) CBF are shown when recovered values of the first tissue model (as shown in (a) and (b)) are normalized to one.
Fig. 6.
Fig. 6. Results of the inverse fitting process when fitting the time-domain diffusion theory model to simulated $g_1(\tau, \rho, t)$ with the MC method. The recovered (a) SBF and (b) CBF coefficients are shown for three values of $t$ when $\rho =1$ cm where the black dashed lines represent the true flow values used in the Monte Carlo simulations. The percent error between recovered absolute values of (c) SBF and (d) CBF is averaged over all the tissue models at three source-detector separations for three $t$ values. We show just the error using the $R_f(\rho, t)$ reflectance model as it was similar to $R(\rho, t)$. Averaged relative error normalized to the first tissue model is shown for values of (e) SBF and (f) CBF.
Fig. 7.
Fig. 7. The sensitivity of the CW $g_1(\tau )$ to an increase of either (a) SBF or (b) CBF by three times at two source-detector separations. A similar example is shown in (c) and (d) for the TD $g_1(\tau )$ at a single source-detector separation ($\rho =1$ cm) and two time values for SBF and CBF, respectively.
Fig. 8.
Fig. 8. The average relative change of $g_1(\tau )$ from increasing the (left column) SBF and (right column) CBF by two times in the (top row) CW and (bottom row) TD at two photon arrival times ($t = 0.5$ and 2.5 ns) as a function of source-detector separation. Two noise models, $\sigma _1(\tau )$ and $\sigma _2(\tau )$, were considered using a collection time, $T_{\mathrm {int}}$, of 1 and 5 seconds, respectively, as computed with Eq. (11). The CW domain shows a varying sensitivity to both SBF and CBF as a function of source-detector separation whereas the TD shows a relatively constant sensitivity.

Equations (11)

Equations on this page are rendered with MathJax. Learn more.

D 2 Φ ( r ) μ a Φ ( r ) = S ( r )
Φ ( ρ , z = 0 ) Φ S I ( ρ ) + 1 π a 2 n = 1 G ^ ( s n , z ) J 0 ( s n ρ ) J 1 2 ( a s n )
G ^ ( s n , z ) = e α 1 ( z + z 0 2 l 1 ) ( 1 e 2 α 1 ( z 0 + z b 1 ) ) ( 1 e 2 α 1 ( z + z b 1 ) ) 2 D 1 α 1 × D 1 α 1 n 1 2 β 3 D 2 α 2 n 2 2 γ 3 D 1 α 1 n 1 2 β 3 ( 1 + e 2 α 1 ( l 1 + z b 1 ) ) + D 2 α 2 n 2 2 γ 3 ( 1 e 2 α 1 ( l 1 + z b 1 ) )
β N = D N 1 α N 1 n n 1 2 ( 1 + e 2 α N 1 l N 1 ) ( 1 e 2 α N ( l N + z b 2 ) ) + D N α N n n 2 ( 1 e 2 α N 1 l N 1 ) ( 1 + e 2 α N ( l N + z b 2 ) ) γ N = D N 1 α N 1 n n 1 2 ( 1 e 2 α N 1 l N 1 ) ( 1 e 2 α N ( l N + z b 2 ) ) + D N α N n n 2 ( 1 + e 2 α N 1 l N 1 ) ( 1 + e 2 α N ( l N + z b 2 ) )
β k 1 = D k 2 α k 2 n k 2 2 ( 1 + e 2 α k 2 l k 2 ) β k + D k 1 α k 1 n k 1 2 ( 1 e 2 α k 2 l k 2 ) γ k γ k 1 = D k 2 α k 2 n k 2 2 ( 1 e 2 α k 2 l k 2 ) β k + D k 1 α k 1 n k 1 2 ( 1 + e 2 α k 2 l k 2 ) γ k
Φ ( ρ , t , z = 0 ) Φ S I ( ρ , t ) + 1 π a 2 n = 1 1 2 π i [ B e s ¯ t G ^ ( s n , s ¯ ) d s ¯ ] × J 0 ( s n ρ ) J 1 2 ( a s n )
R f ( ρ ) = D Φ ( ρ , z ) ( z ) | z = 0
R ( ρ ) = 0.118 Φ ( ρ ) + 0.306 R f ( ρ )
G 1 ( ρ , τ , D b ) = R ( ρ , μ a d )
g 1 ( ρ , τ , D b ) = G 1 ( ρ , τ ) / G 1 ( ρ , τ = 0 )
σ ( τ ) e x p ( ρ ) ( 1 + e x p ( γ τ ) ) T i n t
Select as filters


Select Topics Cancel
© Copyright 2024 | Optica Publishing Group. All rights reserved, including rights for text and data mining and training of artificial technologies or similar technologies.