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

: A novel approach for time-domain diffuse correlation spectroscopy (TD-DCS) has been recently proposed, which has the unique advantage by simultaneous measurements of optical and dynamical properties in a scattering medium. In this study, analytical models for calculating the time-resolved electric-field autocorrelation function is presented for a multi-layer turbid sample, as well as a semi-infinite medium embedded with a small dynamic heterogeneity. To verify the analytical models, we used Monte Carlo simulations, which demonstrated that the theoretical prediction for the time-resolved autocorrelation function was highly consistent with the Monte Carlo simulation, validating the proposed analytical models. Using these analytical models, we also showed that TD-DCS has a higher sensitivity compared to conventional continuous-wave (CW) DCS for detecting the deeper dynamics. The presented analytical models and simulations can be utilized for quantification of optical and dynamical properties from future TD-DCS experimental data as well as for optimization of the experimental design to achieve maximum contrast for deep tissue dynamics. deeper layer dynamics than CW-DCS.


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][2][3][4][5][6][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][5][6][7][8][9][10][11][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 commonlyused 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 semiinfinite 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.

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 1 (τ,s) for a pathlength s [4- , 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).

TD-DCS in a multi-layer turbid sample
For a pulsed point source located at r s , the time-resolved electric-field autocorrelation function ( , , ) 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 2 6 1 B k D τ  , which generally holds in many practical cases when τ <10 −3 s): Vol. 8, No. 12 | 1 Dec 2017 | BIOMEDICAL OPTICS EXPRESS 5520 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 (n) 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 = 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) . 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, ω).
where ρ lies in the z = constant plane. Then we have 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 1 ), identified by n = 1. The solution of Eq. (3) inside the n'th layer can be written as ( , , , ) exp( ) exp( ) where z 0 ~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 ( , , , ) can be written as a ratio: 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.
The time-resolved electric-field autocorrelation function ( , 0, , ) where J 0 denotes the zero-order Bessel function of the first kind, 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 0 ( , is independent of the polar angle. It might be noted that the more accurate form for the extrapolation length is z 0 = 2AD 1 , 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.

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 1 = (x 1 , y 1 , z 1 = 0), a detector at r 3 = (x 3 , y 3 , z 3 = 0) on the top surface of a semi-infinite medium with a small dynamic heterogeneity at r 2 = (x 2 , y 2 , z 2 ), 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 . 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 The detailed expressions are presented below. We denote the electric-field autocorrelation function in a semi-infinite medium as . When a small inclusion (with volume V) with different dynamics is included, the field autocorrelation function is ( , , ) . According to the perturbation approach, ( , , ) G t τ r can be written as 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 is the distance from the source to the detector on the surface, D = v/(3μ s ' ) is the photon diffusion coefficient.

Verification of the analytical model with the Monte Carlo simulation
To verify the analytical model, we calculated g 1 for a 3-layer turbid sample and a semiinfinite medium embedded with a small spherical dynamic heterogeneity. For each case, the analytical prediction for g 1 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 1 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 1 for the pathlength s was then calculated from Eq. (14), 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).

Simulation with a 3-layer model
When using DCS to study the human brain, the head can be approximately modeled as a 3layer 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). 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 1 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 −7 to 10 −3 s, the prediction of the analytical model for g 1 is highly consistent with the Monte Carlo simulation.
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.

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 −1 , μ s ' = 7.6 cm −1 , and the refractive index n r = 1.4. The scattering particle diffusion coefficient D B = 1 × 10 −8 cm 2 /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 −8 cm 2 /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 1 ). Thus the locations for the source, heterogeneity inclusion and detector were r 1 = (0, 0, 0), r 2 = (0, 0, 1.5 cm) and r 3 = (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.
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 9 launched and 143360025 received photon packages.

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 2 ) was considered, because in real measurements, g 2 is recorded instead of g 1 , and one can derive g 2 from g 1 via the Siegert relation, g 2 = 1 + β|g 1 | 2 . Here β is an intercept on the g 2 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 2 of CW-DCS and TD-DCS (for s = 20 cm) with single-mode fiber detection is shown in Fig. 6. The change in g 2 between the baseline and stimulation (Δg 2 ) is also shown on the right panels in Fig. 6(b) and 6(d). The maximum change in g 2 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. In a TD-DCS measurement, a variety of pathlengths can be selected for g 2 . 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 −8 /cm 2 ps for each incident photon, while at the pathlength s = 30 cm, R would be 5.88 × 10 −9 /cm 2 ps. In actual experiments, lower photon count rates would give rise to lower signal to noise ratio on g 2 .
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 2 , we performed a simulation for g 2 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 2 was 5 s, which were set the same for both the CW-DCS and TD-DCS measurements. The same 3layer 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 −6 /cm 2 ps and 1.528 × 10 −6 /cm 2 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 2 (or g 1 ) is known, e.g., from the 3-layer model), g 2 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 2 ; however, due to the limited photon count rate (28.5 kHz) and integration time (5 s), it hardly differentiates the simulated (or 'measured') g 2 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 2 shows significant difference, in particular, when the lag time τ is within the range from 3 × 10 −6 to 5 × 10 −5 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 2 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)]. 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 1 with longer pathlengths, such as 15 to 20 cm.

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 1 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 2 as the mean squared velocity of the scatterers), the present theory still holds by using (1/6)V 2 τ 2 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. 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 1 ) 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 1 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 −1 to 300 cm −1 ] and ω as [-1.2 × 10 13 Hz to 1.2 × 10 13 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.

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