Analytical models for time-domain diffuse correlation spectroscopy for multi-layer and heterogeneous turbid media

Jun Li; Lina Qiu; Chien-Sing Poon; Ulas Sunar

doi:10.1364/BOE.8.005518

1. Introduction

Diffuse correlation spectroscopy (DCS) is an optical technique utilizing the scattering feature of coherent light to probe the dynamical properties of a scattering (or turbid) medium [1–7]. In a conventional DCS, a continuous-wave (CW) laser source is used, which has a long coherent length (e.g., much longer than the typical pathlength of the detected photons), enabling interference of light traveling along different paths [6,7]. The measured temporal autocorrelation function of the detected light intensity is a decay curve with respect to the lag time. The decay curve contains information about the motion of the scatterers; faster decay indicates faster dynamics of the scatterers [4–12]. Since the light detected consists of many photons experiencing various scattering events and traveling along different paths, it is not possible for CW-DCS to differentiate photons with different pathlengths. Therefore, using a CW-DCS system to reveal dynamic heterogeneity in a turbid sample with background dynamics is not straightforward. To deal with this difficulty, various methods including the perturbation model for a semi-infinite medium [4,5,8] and the multi-layer model have been developed [13]. However, due to limited contrast of CW-DCS with respect to deeper dynamics, these approaches are not very robust experimentally when the rate of the detected photons (or photon count rate) is low. In particular, when biological tissue are measured in vivo, the count rate is generally low due to safety standards, which limit the input power of the illuminating light. The higher level of noise on the measured intensity autocorrelation function easily obscures the already limited contrast for deeper dynamics due to the lower count rate.

Although Yodh et al. [14] showed the feasibility of path-length resolved DCS approach previously, the technique used nonlinear optical gating and required high laser powers, which was not suitable for in vivo applications. Very recently, Sutin et al. [15] have reported a novel approach for time-domain (or pathlength-resolved) DCS on phantoms and in a rat brain, which demonstrated the clinical feasibility of this approach. In contrast to a CW-DCS setup, a pulse laser was used as a laser source in TD-DCS. The pulse laser has to have sufficiently short pulse width (tens to hundreds of picoseconds) while maintaining adequate coherence length (several centimeters) [15]. Apart from the laser source, all the other optoelectronic components in this TD-DCS system were very similar to those used in a typical time-resolved spectroscopy (TRS) measurement. Thus, this approach is mostly compatible with commonly-used TRS systems, enabling simultaneous measure of optical parameters of absorption and scattering coefficients and dynamic parameter of scattering particle diffusion coefficient. A combination of such a TD-DCS with a TRS system could provide depth-resolved information on tissue hemodynamics, including blood oxygenation (measured by TRS) and flow (measured by TD-DCS). Thus, it has great potential for biomedical application, such as monitoring cerebral hemodynamics associated with neural activity [15].

To quantify the dynamic parameter from the measured intensity autocorrelation function, various analytical models for CW-DCS have been developed for different measurement geometries [4, 5, 8]. Reflectance geometry is widely adopted in biomedical applications such as brain function monitoring, in which the source and detector are located on the same surface of the sample. For this geometry, the commonly used models include homogeneous semi-infinite and multi-layer models as well as a perturbation model in which a small dynamic heterogeneity embedded in a semi-infinite turbid medium [4,5,8]. However, these currently existing models cannot be directly applied for TD-DCS. Therefore, in this work, we present analytical models for predicting the time-resolved (or pathlength-resolved) electric field autocorrelation function, in particular for a homogeneous multi-layer model and a perturbation model for semi-infinite medium that contains a small dynamic heterogeneity. Monte Carlo simulations were implemented to verify the accuracy of these analytical expressions. By using these analytical and Monte Carlo models we also demonstrate the advantage of TD-DCS over CW-DCS for improved sensitivity of probing deeper dynamics. Thus, we demonstrate that these analytical expressions and Monte Carlo simulations can allow optimizing future experimental protocols for improved depth selectivity.

2. Analytical model for TD-DCS

In a semi-infinite medium characterized by the absorption coefficient μ_a and the reduced scattering coefficient μ_s^’, it is straightforward to obtain the normalized electric-field autocorrelation g₁(τ,s) for a pathlength s [4–10]: $g_{1} (τ, s) = \exp (- \frac{1}{3} k^{2} μ_{s}^{,} s 〈 Δ r^{2} (τ) 〉)$ . k is the wave number in the medium, and <Δr²(τ)> is the mean-squared displacement of the scatterers in time τ. If the scatterers are undergoing diffusive (or Brownian) motion, $〈 Δ r^{2} (τ) 〉 = 6 D_{B} τ$ [4–12], where D_B is the scattering particle diffusion coefficient, then $g_{1} (τ, s) = \exp (- 2 k^{2} μ_{s}^{,} s D_{B} τ)$ , is an exponential function with respect to the lag time τ [15]. For simplicity, in this paper we assume that the scatterers undergo Brownian motion. In the following part of this section, we will focus on implementing the analytical model of TD-DCS for multi-layer turbid sample and a perturbation model for a small heterogeneity embedded within a semi-infinite medium. Here, we only consider reflectance measurement geometry, where the source and detector are located on the same side of the sample surface, since it is commonly used for brain function measurements in both preclinical and clinical studies (see for example the reviews 6-9).

2.1 TD-DCS in a multi-layer turbid sample

For a pulsed point source located at r_s, the time-resolved electric-field autocorrelation function $G (r, t, τ)$ obeys a time-dependent correlation diffusion equation (as long as the diffusion approximation is valid, e.g., photon behavior inside the medium is dominated by scattering, and $6 k^{2} D_{B} τ 〈〈 1$ , which generally holds in many practical cases when τ <10⁻³ s):

[\nabla^{2} - (3 μ_{a} μ {}_{s}^{'}+ 6 μ_{s}^{'}^{2} k^{2} D_{B} τ) - \frac{3 μ_{s}^{'}}{v} \frac{\partial}{\partial t}] G (r, t, τ) = - 3 μ_{s}^{'} δ (r - r {}_{s}) δ (t)

where v is the speed of light in the medium.

Next, we consider a turbid slab consisting of N layers (see Fig. 1) [13], where each layer is characterized by its absorption coefficient μ_a⁽ⁿ⁾ and reduced scattering coefficient μ_s^’(n) (n = 1,2,…,N). The inter-layer boundaries are perpendicular to the z-axis and are located at z = L_n. The front surface of the slab is at z = L₀ = 0 and the back surface is at z = L_N = L. The width of each layer is Δ_n. The source location r_s can be rewritten in polar coordinates as r_s = (ρ = 0,z_s), where z_s = 1/μ_s^’(1) .

Fig. 1 Scheme of the N-layer scattering medium including the position of the source and detector.

Download Full Size | PDF

To solve this equation in a multi-layer medium (along z direction), it is convenient to make a Fourier transform for the real space (ρ, z), as well as the time t, and then solve the equation in the Fourier space (q, z, ω).

\overset{\land}{G} (q, z, ω, τ) = \int_{t} d t \exp (i ω t) \int_{ρ} d^{2} ρ G (ρ, z, t, τ) \exp (i q \cdot ρ)

where ρ lies in the z = constant plane. Then we have

[\frac{\partial}{\partial z^{2}} - 3 μ_{s}^{'} (μ_{a} + 2 μ_{s}^{'} k^{2} D_{B} τ - \frac{i ω}{c}) - q^{2}] \overset{\land}{G} (q, z, ω, τ) = - 3 μ_{s}^{'} δ (z - z_{s})

We divide the top first layer into two sub-layers: layer 0 (0<z<z_s) identified by n = 0, and layer 1 (z_s<z<L₁), identified by n = 1. The solution of Eq. (3) inside the n’th layer can be written as

\overset{\land}{G_{n}} (q, z, ω, τ) = A_{n} \exp (β_{n} z) + B_{n} \exp (- β_{n} z)

where

β_{n}^{2} = q^{2} + 3 μ_{s}^{'}^{(n)} (μ_{a}^{(n)} + 2 μ_{s}^{'}^{(n)} k^{2} D_{B} τ - \frac{i ω}{c})

, A_n and B_n are constant factors for each layer, which can be determined by the boundary conditions of Eq. (3) in the layered medium:

\begin{array}{l} \overset{\land}{G_{0}} (q, z, ω, τ) - z_{0} \frac{\partial}{\partial z} \overset{\land}{G_{0}} (q, z, ω, τ) = 0, z = 0 \\ \overset{\land}{G_{0}} (q, z, ω, τ) = \overset{\land}{G_{1}} (q, z, ω, τ), z = z_{s} \\ \frac{\partial}{\partial z} \overset{\land}{G_{0}} (q, z, ω, τ) = \frac{\partial}{\partial z} \overset{\land}{G_{1}} (q, z, ω, τ) + 3 μ_{s}^{'}^{(1)}, z = z_{s} \\ \overset{\land}{G_{n}} (q, z, ω, τ) = \overset{\land}{G_{n + 1}} (q, z, ω, τ), z = L_{n}, n = 1... N - 1 \\ D_{n} \frac{\partial}{\partial z} \overset{\land}{G_{n}} (q, z, ω, τ) = D_{n + 1} \frac{\partial}{\partial z} \overset{\land}{G_{n + 1}} (q, z, ω, τ), z = L_{n}, n = 1... N - 1 \\ \overset{\land}{G_{N}} (q, z, ω, τ) + z_{N} \frac{\partial}{\partial z} \overset{\land}{G_{N}} (q, z, ω, τ) = 0, z = L_{N} \end{array}

where z₀ ~1/μ_s^’(1) and z_N ~1/μ_s^’(N) are the extrapolation lengths taking into account of internal reflection at the external (z = 0 and z = L, respectively) boundaries of the slab. D_n = v/(3μ_s^’(n)) denotes the photon diffusion coefficient in layer n.

With the above boundary conditions, the constant parameters An and Bn can be determined for each layer, and then $\overset{\land}{G_{n}} (q, z, ω, τ)$ can be obtained. For the top layer (layer 0) $\overset{\land}{G_{0}} (q, z, ω, τ)$ can be written as a ratio:

\overset{\land}{G_{0}} (q, z = 0, ω, τ) = \frac{N u m e r a t o r}{D e n o m i n a t o r}

The detailed expression of the Numerator and Denominator in Eq. (6) is tedious for each case (e.g., depending on the number of layers), thus we give equations below only for N = 1, 2 and 3, with the bottom layer being semi-infinite.N = 1:

\begin{matrix} N u m e r a t o r = 3 μ_{s}^{'}^{(1)} z_{0} \exp (- β_{1} z_{s}) \\ D e n o \min a t o r = 1 + β_{1} z_{0} \end{matrix}

N = 2:

\begin{matrix} N u m e r a t o r = 3 μ_{s}^{'}^{(1)} z_{0} {β_{1} D_{1} \cos h [β_{1} (Δ_{1} - z {}_{s})] + β_{2} D_{2} \sin h [β_{1} (Δ_{1} - z {}_{s})]} \\ D e n o \min a t o r = β_{1} (D {}_{1}+ β_{2} D_{2} z_{0}) \cosh (β_{1} Δ_{1}) + (β_{2} D_{2} + β_{1}^{2} D_{1} z_{0}) \sin h (β {}_{1}Δ_{1}) \end{matrix}

N = 3:

\begin{matrix} N u m e r a t o r = 3 μ_{s}^{'}^{(1)} z_{0} (β_{1} D_{1} \cos h (β_{1} (Δ_{1} - z {}_{s})) (β_{2} D_{2} \cos h (β_{2} Δ_{2}) + β_{3} D_{3} \sin h (β_{2} Δ_{2})) \\ + β_{2} D_{2} (β {}_{3}D_{3} \cos h (β_{2} Δ_{2}) + β_{2} D_{2} \sin h (β_{2} Δ_{2})) \sin h (β_{1} (Δ_{1} - z {}_{s}))) \\ D e n o \min a t o r = β_{2} D_{2} \cos h (β_{2} Δ_{2}) (β_{1} (D {}_{1}+ β_{3} D_{3} z_{0}) \cos h (β_{1} Δ_{1}) + (β_{3} D_{3} + β_{1}^{2} D_{1} z_{0}) \sin h (β_{1} Δ_{1})) \\ + (β {}_{1}(_{β_{3}}^{D_{1}}_{D_{3}}^{+}_{β_{2}}^{^{2}}_{D_{2}}^{^{2}}_{z_{0}}^{)} \cos h (β_{1} Δ_{1}) + (β {}_{2}{^{2}} D_{2}^{2} + β_{1}^{2} β_{3} D_{1} D_{3} z_{0}) \sin h (β_{1} Δ_{1})) \times \sin h (β_{2} Δ_{2}) \end{matrix}

The time-resolved electric-field autocorrelation function $G (ρ, z = 0, t, τ)$ measured on the top surface (z = 0) of the slab is the inverse Fourier transform of $\overset{\land}{G_{0}} (q, z = 0, ω, τ)$

\begin{matrix} G (ρ, z = 0, t, τ) = \frac{1}{2 π} \int_{ω} d ω \frac{1}{{(2 π)}^{2}} \int_{q} d^{2} q \overset{\land}{G_{0}} (q, z = 0, ω, τ) \exp (- i q \cdot ρ) \exp (- i ω t) \\ = \frac{1}{{(2 π)}^{2}} \int_{ω} d ω \int_{q} d q \overset{\land}{G_{0}} (q, z = 0, ω, τ) q J_{0} (ρ q) \exp (- i ω t) \end{matrix}

where J₀ denotes the zero-order Bessel function of the first kind,

J_{0} (ρ q) = \frac{1}{2 π} \int_{0}^{2 π} \exp (- i ρ q \sin θ) d θ

. The last expression in Eq. (10) comes out is because that we can define the location of detector as ρ = (ρ_x, ρ_y) = (0, ρ) due to the axial symmetric geometry, and

\overset{\land}{G_{0}} (q, z = 0, ω, τ)

is independent of the polar angle. It might be noted that the more accurate form for the extrapolation length is z₀ = 2AD₁, A = 2.9486 for n_r = 1.4, A = 2.5153 for n_r = 1.33, respectively, where n_r is the refractive index of the medium.

2.2 TD-DCS analytical model for a semi-infinite medium embedded with a small dynamic heterogeneity

Assume that a laser point source is located at r₁ = (x₁, y₁, z₁ = 0), a detector at r₃ = (x₃, y₃, z₃ = 0) on the top surface of a semi-infinite medium with a small dynamic heterogeneity at r₂ = (x₂, y₂, z₂), as shown in Fig. 2. The medium is characterized by optical parameters including the absorption coefficient μ_a, reduced scattering coefficient μ_s^’, and dynamic parameter D_B. Here we only consider dynamic heterogeneity, which implies that the small inclusion has the same optical properties, but different dynamic parameter D_B + δD_B.

Fig. 2 Geometric scheme for the perturbation model including the position of the source, detector and heterogeneity inclusion.

Download Full Size | PDF

Considering the similarity of the correlation diffusion equation [Eq. (1)] to the photon diffusion equation, the perturbation model developed for photon diffusion in a semi-infinite medium embedded with a small inclusion can be adopted [16], but with a simple replacement $μ_{a} \to μ_{a} + 2 μ_{s}^{'} D {}_{B}k^{2} τ$ [4, 5, 8]. The detailed expressions are presented below.

We denote the electric-field autocorrelation function in a semi-infinite medium as $G^{0} (r, t, τ)$ . When a small inclusion (with volume V) with different dynamics is included, the field autocorrelation function is $G (r, t, τ)$ . According to the perturbation approach, $G (r, t, τ)$ can be written as $G (r, t, τ) = G^{0} (r, t, τ) + δ G (r, t, τ)$ , where the second term is related to the perturbation induced by the inclusion [4,5].

In a semi-infinite medium, the field autocorrelation on the surface (z = 0) is

\begin{matrix} G^{0} (ρ, t, τ) = \frac{v}{{(4 π D t)}^{3 / 2}} \exp (- \frac{ρ^{2}}{4 D t} - (μ_{a} + 2 μ_{s}^{'} D_{B} k^{2} τ) v t) \times \\ {\exp (- \frac{z_{s}^{2}}{4 D t}) - \exp [- \frac{{(z {}_{s}+ 2 z_{0})}^{2}}{4 D t}]} \end{matrix}

_{$ρ = \sqrt{{(x_{3} - x_{1})}^{2} + {(y_{3} - y_{1})}^{2}}$} is the distance from the source to the detector on the surface, D = v/(3μ_s^’) is the photon diffusion coefficient.

\begin{matrix} δ G (r_{1}, r_{2}, r_{3}, t, τ) = - \frac{v^{2}}{{(4 π D)}^{5 / 2} t^{3 / 2}} \exp [- (μ_{a} + 2 μ_{s}^{'} D {}_{B}k^{2} τ) v t] \times \int_{v} d^{3} r_{2} (2 μ_{s}^{'} δ D_{B} k^{2} τ) \times \\ {\begin{cases} (\frac{1}{ρ_{12}^{+}} + \frac{1}{ρ {}_{23}{^{+}}}) \exp [- \frac{{(ρ_{12}^{+} + ρ {}_{23}{^{+}})}^{2}}{4 D t}] - (\frac{1}{ρ_{12}^{+}} + \frac{1}{ρ {}_{23}{^{-}}}) \exp [- \frac{{(ρ_{12}^{+} + ρ {}_{23}{^{-}})}^{2}}{4 D t}] - \\ (\frac{1}{ρ_{12}^{-}} + \frac{1}{ρ {}_{23}{^{+}}}) \exp [- \frac{{(ρ_{12}^{-} + ρ {}_{23}{^{+}})}^{2}}{4 D t}] + (\frac{1}{ρ_{12}^{-}} + \frac{1}{ρ {}_{23}{^{-}}}) \exp [- \frac{{(ρ_{12}^{-} + ρ {}_{23}{^{-}})}^{2}}{4 D t}] \end{cases}} \end{matrix}

where

\begin{matrix} ρ_{12}^{+} = {[{(x_{2} - x_{1})}^{2} + {(y {}_{2}- y_{1})}^{2} + {(z_{2})}^{2}]}^{\frac{1}{2}} \\ ρ_{12}^{-} = {[{(x_{2} - x_{1})}^{2} + {(y {}_{2}- y_{1})}^{2} + {(z_{2} + 2 z_{0})}^{2}]}^{\frac{1}{2}} \\ ρ_{23}^{+} = {[{(x_{3} - x_{21})}^{2} + {(y {}_{3}- y_{2})}^{2} + {(z_{2})}^{2}]}^{\frac{1}{2}} \\ ρ_{23}^{-} = [{(x_{3} - x_{2})}^{2} + {(y {}_{3}- y_{2})}^{2} + {(2 z_{0} + z {}_{2})}^{2}] \end{matrix}

3. Verification of the analytical model with the Monte Carlo simulation

To verify the analytical model, we calculated g₁ for a 3-layer turbid sample and a semi-infinite medium embedded with a small spherical dynamic heterogeneity. For each case, the analytical prediction for g₁ was compared to the Monte Carlo (MC) simulation.

In MC simulation, the survival weight and partial pathlength in each layer (or heterogeneity component) was recorded for each emitted photon package. When calculating g₁ for a certain pathlength s, we chose a pathlength window (e.g., from 0.9∙s to 1.1∙s) and selected all photon packages (e.g., M photon packages in total) whose pathlengths fell within this window. The relative weight (or probability) W(s_i), i = 1:M, was calculated for each photon package among all selected M photon packages. The g₁ for the pathlength s was then calculated from Eq. (14),

g_{1} (τ, s) = \sum_{i = 1}^{M} W (s {}_{i}) \exp (- 2 k^{2} τ \sum_{j}^{N} μ_{s}^{, (j)} s_{i j} D_{B}^{(j)})

The index j denotes each layer (or heterogeneity component) from 1 to N (e.g., for a 3-layer model, N = 3), and s_ij is the partial pathlength of the photon package i in layer j (or component j).

3.1 Simulation with a 3-layer model

When using DCS to study the human brain, the head can be approximately modeled as a 3-layer turbid medium consisting of the scalp, skull and brain [13]. The optical and dynamic parameters used for a head-like 3-layer model are listed in Table 1 [15, 17]. The wavelength of the illuminating light was 785 nm, and the refractive index n_r of each layer was 1.4. In this head-like model, the skull layer was assumed to be static (i.e., D_B = 0).

Table 1. Parameters for the 3-layer model (head-like model)

View Table

Simulations were performed for source-detector separations of 3.0 cm and 0 cm, respectively. Three pathlengths were selected for TD-DCS: 10 cm, 20 cm and 25 cm. Comparison of g₁ between the theory and MC simulation are shown in Fig. 3 for 3 cm and Fig. 4 for 0 cm source-detector separation, respectively. It is clear that for a wide range of lag times, e.g., from 10⁻⁷ to 10⁻³ s, the prediction of the analytical model for g₁ is highly consistent with the Monte Carlo simulation.

Fig. 3 Simulation result of the 3-layer head-like model for the source-detector separation = 3.0 cm: the normalized electric-field autocorrelation function g₁ from the analytical model and Monte Carlo simulation for the three selected pathlengths s = 10 (red), 20 (blue), and 25 cm (black), respectively.

Download Full Size | PDF

Fig. 4 Simulation result of the 3-layer head-like model for the source-detector separation = 0.0 cm: the normalized electric-field autocorrelation function g₁ from the analytical model and Monte Carlo simulation for the three selected pathlengths s = 10 (red), 20 (blue), and 25 cm (black), respectively.

Download Full Size | PDF

The MC simulation for the 3-layer model was performed on a computer [Intel (R) Q9550 (2.83G), Quad-core processor, 4G RAM] with CUDA accelerated Monte Carlo method [18]. This is a fast MC algorithm for the layered model. For the source-detector separations of 0 cm, the running time was about 41.7 s with 44918563 launched and 10021236 received photon packages; while for the separation of 3 cm, the running time was about 1586.3 s with 1351643506 launched and 10002018 received photon packages.

3.2 Simulation for a semi-infinite medium with a small inclusion of dynamic heterogeneity

In this simulation, we assumed the optical parameters for the medium are μ_a = 0.16 cm⁻¹, μ_s^’ = 7.6 cm⁻¹, and the refractive index n_r = 1.4. The scattering particle diffusion coefficient D_B = 1 × 10⁻⁸ cm²/s. A small spherical inclusion was located in the sample right beneath the source, with the same optical parameters, but different dynamic parameter D_BI = 6 × 10⁻⁸ cm²/s. The sphere had a radius of 0.5 cm, with a distance of 1.5 cm from its center to the sample surface. We used the null source-detector separation (i.e., source-detector separation = 0.0 cm) to detect the normalized electric-field autocorrelation function (g₁). Thus the locations for the source, heterogeneity inclusion and detector were r₁ = (0, 0, 0), r₂ = (0, 0, 1.5 cm) and r₃ = (0, 0, 0), respectively. Figure 5 shows the results from the analytical model and Monte Carlo simulations for four selected pathlengths, s = 5, 10, 20, and 40 cm, respectively. The results from the model are again highly consistent with the Monte Carlo simulation.

Fig. 5 Simulation result for a semi-infinite medium embedded with a small spherical dynamic heterogeneity: The source-detector separation = 0 cm. The normalized electric-field autocorrelation function g₁ obtained by the analytical model and Monte Carlo (MC) simulations for the four selected pathlengths s = 5 (red), 10 (pink), 20 (blue) and 40 cm (black), respectively.

Download Full Size | PDF

This simulation was performed on a computer [Intel Xeon E5-2670 (2.3G), 2-core processor, and 64 G memory] with a mesh-based Monte Carlo method (MMCM) [19]. The running time was about 20423.6 s with 10⁹ launched and 143360025 received photon packages.

3.3 Sensitivity to deeper layer dynamics

A remarkable advantage of TD-DCS over CW-DCS is that, with pathlength-resolved measurements, it is possible to differentiate dynamics from different depths in the medium, and at the same time it is also possible to achieve higher contrast to the change in dynamics in deeper layers, such as cerebral blood flow change in the cortex induced by functional activation. To demonstrate this, we used the analytical model presented above for the 3-layer head-like model. The optical and dynamic parameters used in the simulation were the same as those listed in Table 1. We assumed that the cortical (the 3rd layer) dynamics were enhanced by 50% during activation (or stimulation) of the brain. The source-detector separation was set at 3.0 cm. In this simulation, the light intensity autocorrelation function (g₂) was considered, because in real measurements, g₂ is recorded instead of g₁, and one can derive g₂ from g₁ via the Siegert relation, g₂ = 1 + β|g₁|². Here β is an intercept on the g₂ axis that is dependent on the laser coherence length and the detection fiber. If the coherence length is long enough, photon detection with a single mode fiber gives rise to β = 0.5. The simulation of g₂ of CW-DCS and TD-DCS (for s = 20 cm) with single-mode fiber detection is shown in Fig. 6. The change in g₂ between the baseline and stimulation (Δg₂) is also shown on the right panels in Fig. 6(b) and 6(d). The maximum change in g₂ in CW-DCS is about 0.025, while in TD-DCS with s = 20 cm, it is about 0.04, indicating that TD-DCS (with s = 20 cm) has higher sensitivity to deeper layer dynamics than CW-DCS.

Fig. 6 (a) The normalized intensity autocorrelation function g₂-1 of CW-DCS. (b) The change in g₂ between baseline and stimulation for CW-DCS. (c) The normalized intensity autocorrelation function g₂-1 of TD-DCS with a selected pathlength s = 20 cm for the baseline and stimulation. (d) The change in g₂ between baseline and stimulation for TD-DCS.

Download Full Size | PDF

In a TD-DCS measurement, a variety of pathlengths can be selected for g₂. In principal, a longer pathlength can give rise to a higher contrast, as shown in Fig. 7(a). However, with the increase of the photon pathlength, the number of detected diffuse reflected photons dramatically decreases, as seen in Fig. 7(b). For example, at the pathlength s = 20 cm, the diffuse reflectance, R, is 3.35 × 10⁻⁸/cm²ps for each incident photon, while at the pathlength s = 30 cm, R would be 5.88 × 10⁻⁹/cm²ps. In actual experiments, lower photon count rates would give rise to lower signal to noise ratio on g₂.

Fig. 7 (a) The change in g₂ between the baseline and stimulation in the 3-layer model (source-detector separation = 3.0 cm) for different pathlengths. (b) The diffuse reflectance, R (photons/cm²ps) with respect to the photon pathlength.

Download Full Size | PDF

To further demonstrate that TD-DCS provides higher contrast to deeper dynamics compared to CW-DCS, and to illustrate the influence of photon count rate on measured g₂, we performed a simulation for g₂ under real experimental conditions. We assumed the incident power of laser (785 nm) was 50 mW, and the recording time (or integration time) of g₂ was 5 s, which were set the same for both the CW-DCS and TD-DCS measurements. The same 3-layer model (as shown in Table 1) was adopted. During activation, the cortical dynamics were assumed to increase by 50%. In both simulated experiments, a single-mode fiber of 5 micron diameter, numerical aperture (NA) of 0.22 was assumed to detect emitted photons 3 cm away from the source. The quantum efficiency of the photon detector (such as avalanche photodiode, APD) was also considered: 65% for CW-DCS (representing the red-enhanced single-photon counting detector from Excelitas, Quebec, Canada), and 40% (representing the red-enhanced single-photon avalanche diode detector from SPADlab, Politechnico di Milano and MPD Srl) for TD-DCS. The repetition rate of the pulsed laser for the TD-DCS was set at 150 MHz. With these assumptions, we estimated the count rate of detected photons in CW-DCS to be 28.5 kHz. For the TD-DCS measurement, we further assumed the time-window for collecting photons with a selected pathlength as 100 ps. The photons detected for each incident pulse was 8.703 × 10⁻⁶ /cm²ps and 1.528 × 10⁻⁶ /cm²ps for the pathlengths of 20 and 30 cm, respectively.

By generating a photon sequence with a known count rate and correlation between photons (which can be determined if the g₂ (or g₁) is known, e.g., from the 3-layer model), g₂ can be simulated [20], as shown in Fig. 8. In CW-DCS, there is clearly a visible difference between the baseline (blue line) and stimulation (pink line) on the theoretical prediction of g₂; however, due to the limited photon count rate (28.5 kHz) and integration time (5 s), it hardly differentiates the simulated (or ‘measured’) g₂ between the baseline (black line) and stimulation (red line) [Fig. 8(a)]. In TD-DCS, when the selected photon pathlength s = 20 cm [Fig. 8(b)], the two ‘measured’ g₂ shows significant difference, in particular, when the lag time τ is within the range from 3 × 10⁻⁶ to 5 × 10⁻⁵s. This large difference (contrast) enables TD-DCS to more accurately reveal the changes in deep dynamics (such as cerebral blood flow). The difference in the analytical g₂ between the baseline and stimulation is much larger for s = 30 cm than s = 20 cm, because photons with longer paths are likely to reach deeper layer where the changes in dynamics occurs. However due to the lower photon count rates, detected for each incident pulse at s = 30 cm was 17.6% of that when at s = 20 cm, the contrast between the baseline and stimulation at s = 30 cm is worse than s = 20 cm [Fig. 8(c)].

Fig. 8 The autocorrelation function g₂ in (a) CW-DCS and TD-DCS for selected pathlength (b) s = 20 cm and (c) s = 30 cm. In each case, the theoretical prediction from the 3-layer model and simulation on g₂ are presented for the baseline and stimulation.

Download Full Size | PDF

In contrast to CW-DCS, TD-DCS is a pathlength-resolved measurement, which means that even for a fixed source-detector separation; photons from selected pathlengths can be used to generate the autocorrelation function [15]. Therefore, when the influence of the photon count is taken into account, an optimal pathlength may exist for which the best differentiation can be achieved. On the other hand, unlike CW-DCS, TD-DCS may not necessarily require larger source-detector separation to probe for deeper dynamics, as recently shown by Sutin et al. [15]. By selecting a longer path, it is possible to reveal the deeper dynamics even with shorter source-detector separation. To elucidate this, we used the 3-layer head-like model to perform a simulation with zero source-detector separation, as shown in Fig. 9. With zero (null) source-detector separation, as Torricelli et al [21] and Pifferi et all [22] have shown for time resolved spectroscopy, the change in the deeper (3rd layer) dynamics can be readily revealed in g₁ with longer pathlengths, such as 15 to 20 cm.

Fig. 9 The pathlength-resolved normalized electric field autocorrelation function g₁ for the 3-layer head-like model with a source-detector separation of 0 cm for the baseline and stimulation with pathlengths of s = 10, 15, and 20 cm. The parameters used for the stimulation were the same as those used in Figs. 7 and 8, except for the source-detector separation.

Download Full Size | PDF

4. Discussion and conclusion

We have presented analytical expressions for the time-resolved electric-field autocorrelation function for a multi-layer homogeneous turbid media as well as for a semi-infinite medium embedded with a small dynamic heterogeneity. These two models are closely related to the biomedical applications of DCS such as non-invasive measurement of cerebral blood flow associated with brain activation through intact scalp and skull, and the detection of a small tumor in breast tissue. The presented analytical models were validated with the Monte Carlo (MC) simulations, showing that in a wide range of lag times, the theoretical prediction of electric-field autocorrelation functions were highly consistent with the MC simulations.

The present theoretical models are based on the diffusion approximation, which implies the detected photons should undergo a lot of scattering events. In the simulations, the null (zero) source-detector separation [21] was used, while the selected pathlength s for calculating the pathlength-resolved g₁ was at least 5 cm, much longer than the transport mean free path length (e.g., 0.09-0.15 cm). Thus, the diffusion approximation was valid for describing transportation of the detected photons, which was demonstrated by the Monte Carlo simulations. However, in a real experiment, measurement with a null (or very short) source-detector separation may induce saturation of photon detector (such as avalanche photodiode, APD). To deal with this problem, Pifferi et al proposed a gated-detection approach, where an APD was operated in time-gated mode to prevent detection of the early photons for enhanced contribution of late photons [22].

In this study, we only considered Brownian motion as dynamics in the sample. It is because that many studies have demonstrated Brownian motion model can better explain the experimental data measured from living tissues including the human brain [6,8,13]. However, if other dynamics model is adopted such as the random flow (with V² as the mean squared velocity of the scatterers), the present theory still holds by using (1/6)V²τ² to replace D_Bτ in the equations, e.g., Eqs. (1), (11) and (12). It should be noted that in recent years, more accurate models for describing the dynamics in tissue has been investigated [10,11,23]. For example, the hydrodynamic diffusion model which accounts for not only Brownian diffuse motion, but also ballistic motion at short lag times. In this model the mean squared displacement $〈 Δ r^{2} (τ) 〉 = 6 D {}_{H}{τ - τ_{c} [1 - \exp (- τ / τ_{c})]}$ , where D_H is called the effective hydrodynamic diffusion coefficient and τ_c is the time scale required to establish diffusive motion. If one assumes dynamic model for the present theory, one can use the term $D {}_{H}{τ - τ_{c} [1 - \exp (- τ / τ_{c})]}$ instead of D_Bτ in the corresponding equations.

By using these analytical and Monte Carlo models we also demonstrate the advantage of TD-DCS over conventional continuous-wave DCS (CW-DCS) for improved sensitivity of probing deeper dynamics. Since TD-DCS provides pathlength-resolved measurements, it allows discriminating dynamics originating from different layers. We showed that by selecting the source-detector separation and pathlength, one can achieve higher contrast to the changes in dynamics in deeper layers, such as cerebral blood flow change in the cortex induced by functional activation. Furthermore, one can use a fixed source-detector separation and select various pathlengths to generate the autocorrelation function [15]. Thus, by considering optimal photon counts and detector noise levels, one can achieve high sensitivity to deeper dynamics even with short source-detector separations. As Torricelli et al. [21] has shown the concept and advantages of using the null (zero) source-detector separation in time resolved spectroscopy measurements, we used a 3-layer head-like model to perform a simulation with null source-detector separation and have successfully demonstrated sensitivity to the deeper (3rd layer) dynamics.

Although the results obtained by the analytical model for the electric-field autocorrelation function (g₁) matches well with the Monte Carlo simulations, there can be several potential sources of errors. The diffusion approximation used for solving the analytical solution of g₁ might be one reason for this discrepancy. In addition, numerical calculation errors could contribute to the differences between the analytical model and Monte Carlo simulations. For the multi-layer model solution, the integral bound for q theoretically should be from 0 to + ∞, while for ω from -∞ to + ∞ in Eq. (10). However in practice, the numerical integration was performed with limited range for both q and ω. For example, in our calculation we used the integral range for q as [0 cm⁻¹ to 300 cm⁻¹] and ω as [-1.2 × 10¹³ Hz to 1.2 × 10¹³ Hz]. Furthermore, for the small dynamic heterogeneity in a semi-infinite medium analytical solution, the model error could originate from the perturbation approximation, especially when the inclusion volume is larger and closer to the surface..

In conclusion, these analytical solutions provide modeling layered and heterogeneous tissue, closely representing the clinical scenarios like quantification of dynamics in brain and breast tissue. Thus, they can be widely used in future experiments for characterization of dynamic contrasts in a variety of preclinical and clinical settings. In particular, with known optical parameters of tissues such as those obtained from the time-resolved spectroscopy (TRS) measurements, it is possible to quantify accurate blood flow index in deep layers by fitting the measured time-resolved data to the theoretical models. Moreover, the analytical solution based simulations can be utilized for optimization of experimental settings such as choosing optimal source-detector separation, photon pathlength selection, which would allow the optimal contrast for deep tissue. Thus, the presented solutions provide opportunities for quantitative implementation of TD-DCS in clinical and preclinical settings.

Funding

A grant from the Ohio Third Frontier to the Ohio Imaging Research and Innovation Network (OIRAIN, 667750); a grant for Guangdong Provincial Key Laboratory of Optical Information Materials and Technology (Grant No. 2017B030301007), and a grant from National Natural Science Foundation of China (Grant No. 81771876).

Disclosures

The authors declare that there are no conflicts of interest related to this article.

References and links

1. G. Maret and P. E. Wolf, “Multiple light scattering from disordered media: The effect of Brownian motion of scatterers,” Z. Phys. B Condens. Matter 65, 409–413 (1987).

2. D. J. Pine, D. A. Weitz, P. M. Chaikin, and E. Herbolzheimer, “Diffusing wave spectroscopy,” Phys. Rev. Lett. 60(12), 1134–1137 (1988). [PubMed]

3. D. J. Pine, D. A. Weitz, J. X. Zhu, and E. Herbolzheimer, “Diffusing-wave spectroscopy: dynamic light scattering in the multiple scattering limit,” J. Phys. (France) 51, 2101–2127 (1990).

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

5. D. A. Boas and A. G. Yodh, “Spatially varying dynamical properties of turbid media probed with diffusing temporal light correlation,” J. Opt. Soc. Am. A 14, 192–215 (1997).

6. R. C. Mesquita, T. Durduran, G. Yu, E. M. Buckley, M. N. Kim, C. Zhou, R. Choe, U. Sunar, and A. G. Yodh, “Direct measurement of tissue blood flow and metabolism with diffuse optics,” Philos Trans A Math Phys. Eng. Sci. 369, 4390–4406 (2011).

7. D. A. Boas, C. Pitris, and N. Ramanujam, Handbook of Biomedical Optics, (CRC Press, 2011), Chap. 10.

8. 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). [PubMed]

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

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

11. S. A. Carp, G. P. Dai, D. A. Boas, M. A. Franceschini, and Y. R. 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). [PubMed]

12. 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). [PubMed]

13. 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), 44002 (2005). [PubMed]

14. A. G. Yodh, P. D. Kaplan, and D. J. Pine, “Pulsed diffusing-wave spectroscopy: High resolution through nonlinear optical gating,” Phys. Rev. B Condens. Matter 42(7), 4744–4747 (1990). [PubMed]

15. 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). [PubMed]

16. F. Martelli, S. Del Bianco, A. Ismaelli, and G. Zaccanti, Light Propagation through Biological Tissue and Other Diffusive Media: Theory. Solutions, and Software (SPIE Press, Bellingham, 2009).

17. G. Strangman, M. A. Franceschini, and D. A. Boas, “Factors affecting the accuracy of near-infrared spectroscopy concentration calculations for focal changes in oxygenation parameters,” Neuroimage 18(4), 865–879 (2003). [PubMed]

18. E. Alerstam, T. Svensson, and S. Andersson-Engels, “Parallel computing with graphics processing units for high-speed Monte Carlo simulation of photon migration,” J. Biomed. Opt. 13(6), 060504 (2008). [PubMed]

19. Q. Fang, “Mesh-based Monte Carlo method using fast ray-tracing in Plücker coordinates,” Biomed. Opt. Express 1(1), 165–175 (2010). [PubMed]

20. F. Ferri and D. Magatti, “Hardware simulator for photon correlation spectroscopy,” Rev. Sci. Instrum. 74(10), 4273–4279 (2003).

21. A. Torricelli, A. Pifferi, L. Spinelli, R. Cubeddu, F. Martelli, S. Del Bianco, and G. Zaccanti, “Time-resolved reflectance at null source-detector separation: Improving contrast and resolution in diffuse optical imaging,” Phys. Rev. Lett. 95(7), 078101 (2005). [PubMed]

22. A. Pifferi, A. Torricelli, L. Spinelli, D. Contini, R. Cubeddu, F. Martelli, G. Zaccanti, A. Tosi, A. Dalla Mora, F. Zappa, and S. Cova, “Time-resolved diffuse reflectance using small source-detector separation and fast single-photon gating,” Phys. Rev. Lett. 100(13), 138101 (2008). [PubMed]

23. K. Verdecchia, M. Diop, L. B. Morrison, T. Y. Lee, and K. St Lawrence, “Assessment of the best flow model to characterize diffuse correlation spectroscopy data acquired directly on the brain,” Biomed. Opt. Express 6(11), 4288–4301 (2015). [PubMed]

	Thickness (cm)	μ_a (cm⁻¹)	μ_s^’ (cm⁻¹)	D_B (cm²/s)
Layer 1 (scalp)	0.385	0.191	6.60	1 × 10⁻⁸
Layer 2 (skull)	0.769	0.136	8.60	0
Layer 3 (brain)	infinite	0.186	11.10	6 × 10⁻⁸

Analytical models for time-domain diffuse correlation spectroscopy for multi-layer and heterogeneous turbid media

Abstract

1. Introduction

2. Analytical model for TD-DCS

2.1 TD-DCS in a multi-layer turbid sample

2.2 TD-DCS analytical model for a semi-infinite medium embedded with a small dynamic heterogeneity

3. Verification of the analytical model with the Monte Carlo simulation

3.1 Simulation with a 3-layer model

3.2 Simulation for a semi-infinite medium with a small inclusion of dynamic heterogeneity

3.3 Sensitivity to deeper layer dynamics

4. Discussion and conclusion

Funding

Disclosures

References and links

Cited By

Figures (9)

Tables (1)

Equations (14)

Biomedical Optics Express