1.

Introduction

Cerebral blood flow (CBF) is an important indicator of brain function and health.^1–6 Diffuse correlation spectroscopy (DCS) serves as a major optical technique that is capable of measuring blood flow non-invasively and continuously.^7–11 In DCS, coherent light is incident on the surface of the scattering medium, and the re-emitted scattered light is collected by a detector at a certain distance away from the source position, typically in the range of 10 to 30 mm. Dynamics in the tissue, mainly arising from blood flow induces phase variations of the scattered light that alters the interference pattern of the partial waves of the re-emitted light, namely, the speckle pattern, thus causing light intensity fluctuations with time at the detector. It is common to quantify these temporal fluctuations using a temporal intensity autocorrelation function. The decay rate of this autocorrelation function provides a measure of the blood flow index. Numerical and theoretical tools have been developed to guide interpretations of experimental measurements. For traditional continuous-wave DCS (CW-DCS), the analytical expressions of the autocorrelation functions have been well established in diffusion theory,^7,12,13 and Monte Carlo simulations have been developed to numerically obtain the field autocorrelation function for inhomogeneous media.^8,14 Together with the noise model obtained analytically,¹⁵ the performance of a CW-DCS system for a particular measurement geometry can be theoretically predicted.

Recently, a time-domain variant of DCS (TD-DCS) has been developed, which utilizes a pulsed laser and a gating strategy to select photons detected at different arrival times.^16–20 This technology has the potential to differentiate blood flow information from different layers, such as the skull and the brain tissue using a single source–detector separation. The theoretical basis for TD-DCS has been established taking into account the effects of the instrument response function (IRF), which includes the incident pulse shape and the broadening of the pulse arising from the instrument, and the coherence length of the light source on the shape of the intensity autocorrelation function.¹⁸ Ideally, a narrower width of the IRF will result in better specificity to different photon pathlengths. However, the width of IRF is fundamentally limited by the coherence length, and a narrow pulse in time will result in a short coherence length, which will degrade signal to noise ratio (SNR) in TD-DCS measurements.¹⁸ A pioneering study on the performance of TD-DCS has been conducted assuming the same noise model¹⁵ that is developed for CW-DCS.²¹ A proper noise model that quantifies the performance of TD-DCS measurements in terms of SNR and contrast to noise ratio (CNR) for any IRF, coherence length, measurement geometry, and brain activation pattern is still lacking.

We have constructed a numerical model which is capable of simulating both CW- and TD-DCS measurements from first principles. We have integrated Monte Carlo simulations with wave calculations to obtain the intensity fluctuations induced by tissue dynamics. The speckle intensity is converted to photon counts within each bin time mimicking the experimental measurements. The numerical results are compared with existing theoretical expressions of the average intensity autocorrelation function for CW-DCS¹⁴ and TD-DCS.¹⁸ The simulated noise is also validated with the existing noise model for CW-DCS.¹⁵ We have demonstrated the value of using our model to predict the performance for any TD-DCS systems and measurement geometry in terms of CNR due to brain activation. We have found that for a given set of parameters, brain geometry, and the IRF profile utilized in this paper, CW-DCS at a large source–detector separation (30 mm) out-performs TD-DCS. We can increase the detected photon flux per bin time of TD-DCS at a small source–detector separation (10 mm) by increasing the gate width at the cost of reduced specifity to a particular photon pathlength or depth to achieve a performance comparable to CW-DCS at a medium source–detector separation (20 mm). We could not find realistic measurement conditions under which TD-DCS outperformed CW-DCS. However, since the performance of TD-DCS is highly dependent on the IRF profile²¹ and the other measurement parameters, the quantitative results predicted in this paper can be different for researchers using different input parameters in our model. The performance for a given TD-DCS system can be easily calculated by changing the parameters of the numerical model, which is made publicly available.²² TD-DCS is also capable of differentiating blood flow at different depths using a single source–detector separation. This model provides guidance for the experimental development of novel TD-DCS systems.

2.

Method

2.1.

Theory

DCS quantifies the dynamics of the sample by calculating the autocorrelation function of the speckle intensity fluctuations of light reemitted from tissue. For a single scattering event shown in Fig. 1(a), the electric field measured at the detector is

Eq. (1)

$$\begin{gathered} \overrightarrow{E}\propto {\overrightarrow{E}}_0\sum _{n=1}^N\exp (i{\overrightarrow{k}}_{\text{in}}·{\overrightarrow{r}}_n)\exp (i{\overrightarrow{k}}_{\text{out}}·({\overrightarrow{R}}_d-{\overrightarrow{r}}_n)) \\ \propto {\overrightarrow{E}}_0\exp (i{\overrightarrow{k}}_{\text{out}}·{\overrightarrow{R}}_d)\sum _{n=1}^N\exp (-i{\overrightarrow{q}}_n·{\overrightarrow{r}}_n). \end{gathered}$$

Fig. 1

(a) Illustration of the single scattering event and the parameters used in Eq. (1). (b) The geometry of the Monte Carlo simulation used in this manuscript and the illustration of the gating strategy for TD-DCS. The size of the detector is 2 mm, the thickness of the first layer or first tissue type is 15 mm, the reduced scattering coefficient is ${\mu }_{s^{\prime }}=1{\mathrm{mm}}^{-1}$ and absorption coefficient ${\mu }_a=0.01{\mathrm{mm}}^{-1}$. The optical properties ${\mu }_{s^{\prime }}$ and ${\mu }_a$ are set to be the same for the two layers. The dynamics in the second layer is increased when brain activation is introduced. The source detector separation is set to be $\rho =20\mathrm{mm}$ unless otherwise stated. Measurements at a later gate will be able to select photons with longer pathlengths, thus it is more sensitive to dynamics at deeper depths. (c) Illustration of the calculation of the photon arrival time for the $n$th photon as the sum of the IRF time and the transit time within the sample $t_{na}=t_{n\mathrm{RIF}}+t_{nL}$. (d) The IRF profile we have used in this manuscript obtained experimentally. The full width half maximum of the IRF is 0.31 ns. The dashed lines indicate the levels of the intensity of the IRF that corresponds to 0.1, 0.01, and 0.001 of the peak intensity, and the widths that correspond to these levels are 0.56, 0.82, and 1.18 ns, respectively.

Here, $\overrightarrow{E}$ is the electric field reaching the detector located at ${\overrightarrow{R}}_d$; ${\overrightarrow{r}}_n$ is the position of the $n$’th scatterer; $N$ is the total number of the scatterers; ${\overrightarrow{E}}_0$ is the incident field; and ${\overrightarrow{q}}_n={\overrightarrow{k}}_{\text{out}}-{\overrightarrow{k}}_{\text{in}}$ is the momentum transfer, where ${\overrightarrow{k}}_{\text{out}}$ and ${\overrightarrow{k}}_{\text{in}}$ are the output and incident wave vectors, respectively. We only consider elastic scattering such that $|{\overrightarrow{k}}_{\text{out}}|=|{\overrightarrow{k}}_{\text{in}}|=k$.

For the case of multiple light scattering, Eq. (1) is revised as

Eq. (2)

$$\overrightarrow{E}\propto \sum _{n=1}^N\sum _{s=1}^{N_{sn}}\exp (-i{\overrightarrow{q}}_{ns}·{\overrightarrow{r}}_{ns}).$$

Here, $s$ denotes the $s$’th scattering event, and $N_{sn}$ is the total number of scattering events for the $n$’th component of the wave, or the $n$’th photon in our Monte Carlo simulations as will be described in Sec. 2.2. The field autocorrelation function is calculated from $g_1(\tau )=⟨E^{*}(t)E(t+\tau )⟩/⟨{|E(t)|}^2⟩$. For a semi-infinite medium, the analytical form of $g_1$ for CW-DCS system has been calculated as^7,9

Eq. (3)

$$g_1(\rho ,\tau )=\frac{3C{\mu }_{s^{\prime }}}{4\pi }[\frac{\exp (K{r}_1)}{r_1}-\frac{\exp (K{r}_2)}{r_2}].$$

Here, $K^2=3{\mu }_a{\mu }_{s^{\prime }}+6{\mu }_{s^{\prime }}^2{k}_0^2{D}_B\tau$, $\rho$ is the source–detector distance, $r_1={({\rho }^2+z_0^2)}^{1/2}$, $r_2={({\rho }^2+{(z_0+2z_b)}^2)}^{1/2}$, $z_0=1/{\mu }_{s^{\prime }}$, $z_b=(5/3){\mu }_{s^{\prime }}$, $k=2\pi /\lambda$, $\lambda$ is the wavelength, $k_0$ is the central wave vector, $C$ is the normalization constant, $D_B$ is the diffusion coefficient, and ${\mu }_a$ and ${\mu }_{s^{\prime }}$ are the absorption and reduced scattering coefficients, respectively. In real time measurements, a simpler expression is preferred for fitting and Eq. (3) can be simplified to a single exponential decay at small $\tau$ limit^15,23

Eq. (4)

$$g_1(\tau )=\exp (-\tau /{\tau }_c),$$

where ${\tau }_c$ is the decay time constant that is utilized to quantify the tissue dynamics. Experimentally, instead of $g_1(\tau )$, the intensity correlation function $g_2(\tau )=⟨I(t)I(t+\tau )⟩/{⟨|I(t)|⟩}^2$ is measured, and $g_2(\tau )$ is related to $g_1(\tau )$ via the Siegert relation²⁴

Eq. (5)

$$g_2(\tau )=1+\beta g_1{(\tau )}^2,$$

with $\beta =1$ indicating complete coherence of the detected photons, and $\beta <1$ accounting for loss of temporal coherence and/or detection of multiple modes such as the two polarization states of the electromagnetic field. The effect of the coherence length $l_c$ has been estimated as²⁵

Eq. (6)

$$g_2(\tau )=1+{\int }_0^{\infty }\mathrm{d}L{\int }_0^{\infty }\mathrm{d}{L}^{\prime }P(L)P(L^{\prime })g_1(L,\tau )g_1(L^{\prime },\tau )e^{-16{[(L-L^{\prime })/(\pi l_c)]}^2}.$$

Here, $P(L)$ is the photon pathlength distribution. It is noteworthy that this relation is only valid when the incident light spectrum is a Gaussian function. For a semi-infinite highly scattering medium relevant for human brain measurements, $P(L)$ can be calculated as^26,27

Eq. (7)

$$P(L)=\frac{vS}{{(4\pi DL/v)}^{3/2}}\exp (-{\mu }_{a}L)\left[\exp \right(-\frac{z_{-}^2}{4DL})-\exp (-\frac{z_{+}^2}{4DL}\left)\right].$$

Here, $v$ is the speed of light in the sample, $z_{-}^2=z_0^2+{\rho }^2$, $z_{+}^2={(z_0+2z_b)}^2+{\rho }^2$, $z_0=1/{\mu }_{s^{\prime }}$, $z_b=2D(1+R_{\mathrm{eff}})/(1-R_{\mathrm{eff}})$, $R_{\mathrm{eff}}$ the boundary reflectivity, $D=1/[3({\mu }_a+{\mu }_{s^{\prime }})]$, ${\mu }_a$ the absorption coefficient, ${\mu }_{s^{\prime }}$ the reduced scattering coefficient, and $S$ the normalization factor such that ${\int }_0^{\infty }P(L)\mathrm{d}L=1$. Here, $g_1(L,\tau )$ is the pathlength-dependent normalized field temporal autocorrelation function

Eq. (8)

$$g_1(L,\tau )=\exp (-2{\mu }_{s^{\prime }}{D}_B{k}_0^{2}L\tau ),$$

At $\tau =0$, $g_1(L,0)=1$, which provides the relation between $\beta$ and $l_c$ in Eq. (6).

The noise model for CW-DCS systems has been established analytically.¹⁵ The standard deviation of $(g_2(\tau )-1))$, $\sigma (\tau )$, at each time delay $\tau$ is

Eq. (9)

$$\begin{gathered} \sigma (\tau )=\sqrt{T_b/T}[{\beta }^2\frac{(1+e^{-2\mathrm{\Gamma }{T}_b})(1+e^{-2\mathrm{\Gamma }\tau })+2m(1-e^{-2\mathrm{\Gamma }{T}_b})e^{-2\mathrm{\Gamma }\tau }}{1-e^{-2\mathrm{\Gamma }{T}_b}} \\ +2{⟨n⟩}^{-1}\beta (1+e^{-2\mathrm{\Gamma }\tau })+{⟨n⟩}^{-2}(1+\beta e^{-\mathrm{\Gamma }\tau }){]}^{1/2}. \end{gathered}$$

Here, $T_b$ is the bin time, $m$ is the bin index, $T$ is the measurement time window, $⟨n⟩$ is the average number of photons in a bin time $T_b$, $g_2(\tau )-1=\beta e^{-2\mathrm{\Gamma }\tau }$, where $\mathrm{\Gamma }=1/{\tau }_c$.

For TD-DCS, a gating strategy is utilized to select a portion of the photons that arrive at the detector. For a gate delay time $t_s$ and gate width $t_w$, only the photons that arrive within a measurement window are selected to calculate the electric field in Eq. (2), as shown in Fig. 1(b). Ideally, if both the widths of the IRF and the gate width are infinitely small, the gating strategy will only select photons with a single pathlength $L=t_s*v$, and $g_1(\tau )$ is reduced to $g_1(L=t_s*v,\tau )$ in Eq. (8), where $v$ is the speed of the light in the tissue. However, both the IRF and the gate width will increase the width of the distribution of the photon pathlengths detected within a measurement window. The expression of $g_2(\tau )$ for TD-DCS has been analytically obtained and the effect of the coherence length $l_c$ and IRF are taken into account.¹⁸ At $t_s$, the gate delay time dependent intensity autocorrelation function $g_2(t_s,\tau )$ is

Eq. (10)

$$\begin{gathered} g_2(t_s,\tau )=1+{\int }_0^{\infty }\mathrm{d}L{\int }_0^{\infty }\mathrm{d}{L}^{\prime }P(L)I_p(t_s-L/v)P(L^{\prime }) \\ {I}_p(t_s-L^{\prime }/v)g_{1L}(L,\tau )g_{1L}(L^{\prime },\tau )e^{-16{[(L-L^{\prime })/(\pi l_c)]}^2}. \end{gathered}$$

Here the effective pathlength distribution $P(t_s,L)=P(L)I_p(t_s-L/v)$ is the pathlength distribution of photon trajectories measured at $t_s$, where the normalized IRF profile is denoted as $I_p(t)$. The definition of coherence length is shown in Sec. 2.2, Step 4. It is noteworthy that the definition of the coherence length can be different in other studies.^18,25

2.2.

Monte Carlo Wave Model

We have established a numerical recipe to simulate photon counts received at the detector as a function of time for both CW-DCS and TD-DCS measurements from first principles. The steps needed to construct the model are explained below.

• Step 1: Simulate photon migration within a scattering medium using Monte-Carlo simulations. The code for the Monte Carlo simulation is a derivation of that used in previous publications.^14,23,28 In the Monte Carlo simulation, we specify the source–detector separation and the optical properties of the sample including the scattering coefficient ${\mu }_s$ and anisotropy factor $g$ in different tissue types. We have used a two-layer sample geometry denoted as two tissue types in the Monte Carlo code as shown in Fig. 1(b). Here we have set $g=0$ and ${\mu }_s={\mu }_{s^{\prime }}=1{\mathrm{mm}}^{-1}$ for both tissue types in this study, but these can be easily configured to be different by researchers for future work. In biological tissue, the value of $g$ obtained experimentally is typically above 0.9.²⁹ But in the diffusive regime as we are interested, the DCS results are only governed by the reduced scattering coefficient as indicated in Eq. (3). Thus, we can implement isotropic scattering with $g=0$ and use a smaller ${\mu }_s$ by assigning ${\mu }_s={\mu }_{s^{\prime }}$ in the Monte Carlo simulation to improve computational efficiency, which is equivalent to $g=0.9$, ${\mu }_{s^{\prime }}={\mu }_s(1-g)$. The effect of absorption will be taken into account at a later stage in Step 8. We will assign different blood flow dynamics to the two layers to demonstrate the DCS signal variation induced by blood flow change due to neural activation in the deeper layer. The number of scattering events $N_{sni}$ and pathlength $L_{nti}$ for the $n$’th photon in the $i$ $ti$’th tissue type are recorded. The values of other parameters we have used in the model are specified in the caption of Fig. 1(b).

• Step 2: Calculate the arrival time $t_{na}$ of the $n$’th photon. We set the time at the peak of the IRF to be $t=0$. The transit time that a photon spent within the medium is $t_{nL}=L_n/v$, where $v=c/n$ is the speed of light within the medium, $c=300\mathrm{mm}/\mathrm{ns}$ is the speed of the light in vacuum and the refractive index of the tissue is set to be $n=1.33$. The IRF time $t_{n\mathrm{IRF}}$ is randomly drawn from the probability density distribution determined by the IRF profile, by obtaining the time at which a random number drawn within the interval [0 1] matches the cumulative probability. This takes into account the effect of the IRF on the distribution of the photon pathlengths detected within a specific measurement window. The arrival time of the photon is $t_{na}=t_{nL}+t_{n\mathrm{IRF}}$. In this paper, we have used an experimentally measured IRF as shown in Fig. 1(d). In the experiment, a PicoQuant VisIR laser at 765 nm was used as the source with a PPG512 to trim the pulse, and a PicoQuant PMA42 photomultiplier was used as the detector to reduce the long tail. Any IRF profile obtained experimentally can be utilized in this model to predict the performance for a particular TD-DCS system.

• Step 3: Select photons that fall within the specified measurement gate. For a gate delay time $t_s$ with a gate width $t_w$, only the photons with $t_{na}$ that fall within the time span $[t_s-t_w/2,t_s+t_w/2]$ are selected for the following calculations. It is noteworthy that if both the width of the IRF and the gate width $t_w$ are infinitely small, the photons with a single pathlength $L_n=t_s*v$ are selected, which is the ideal case described by Eq. (8), but is generally not achievable experimentally.

• Step 4: For each photon, assign the spectrum of the input light that corresponds to a particular coherence length $l_c$. Here we have utilized a Gaussian spectrum, and the relation between the Gaussian spectrum in real and $k$ space is related via

  Eq. (11)

  $$S(k)=e^{-{(k-k_0)}^2/(2k_c^2)},S(L)=e^{-L^2/(2l_c^2)},k_c^2={\pi }^2/(4l_c^2).$$

Here, we have defined the coherence length $l_c$ as the width of the Gaussian. It is noteworthy that the definition of the coherence length can be different in other studies. We have combined the photon description of light with the wave description. In principle, we cannot assign a wavelength to a photon. However, each photon can be considered as a light trajectory for wave propagation, which allows us to calculate the accumulated phase along a light path for different wavelengths. In our model, any spectrum shape $S(k)$ measured experimentally can be incorporated, which is one advantage of our numerical model over existing theoretical models,^18,25 where a Gaussian spectrum profile is always assumed.

• Step 5: Assign an initial phase randomly drawn within $[02\pi ]$ for each photon trajectory ${\varphi }_n$.

• Step 6: Obtain the displacements of the scatterers as functions of time. Here we denote the displacement of the scatterer in the $s$’th scattering event of the $n$’th photon trajectory $\mathrm{\Delta }{\overrightarrow{r}}_{ns}(t)$. The displacement is obtained from a second Monte Carlo simulation of random particle motion with a diffusion coefficient $D_B$. The displacement of a scatterer $\mathrm{\Delta }{\overrightarrow{r}}_{ns}(t)=(\overrightarrow{r}(t)-\overrightarrow{r}(0){)}_{ns}$ can be calculated from the evolution of $\overrightarrow{r}(t)=(x(t),y(t),z(t))$

  Eq. (12)

  $$\begin{gathered} x(t+\mathrm{\Delta }t)=x(t)+N(\mathrm{0,2}{D}_B\mathrm{\Delta }t), \\ y(t+\mathrm{\Delta }t)=y(t)+N(\mathrm{0,2}{D}_B\mathrm{\Delta }t), \\ z(t+\mathrm{\Delta }t)=z(t)+N(\mathrm{0,2}{D}_B\mathrm{\Delta }t). \end{gathered}$$

Here, $N(\mathrm{0,2}{D}_B\mathrm{\Delta }t)$ is a normal distribution with mean 0 and variance $2D_B\mathrm{\Delta }t$. It is noteworthy that the motion of the scatterers or red blood cells is in general a combination of convective flow and shear induced diffusion. It has been experimentally observed and numerically demonstrated that the DCS signal is dominated by the diffusive behavior of RBCs red blood cells (RBCs).^14,30,31 Thus, we only model the diffusive motion of RBCs as in Eq. (12). The diffusion coefficient at the baseline, i.e. without brain activation, is set to be $D_B=1*{10}^{-6}{\mathrm{mm}}^2/\mathrm{s}$ in the model, which provides a decay time constant close to experimental measurements.^9,23

• Step 7: Calculate the momentum transfer ${\overrightarrow{q}}_{\mathrm{nis}}=({\overrightarrow{k}}_{\text{out}}-{\overrightarrow{k}}_{\mathrm{in}}{)}_{\mathrm{nis}}$ for each scattering event. Here $i$ denotes the $i$’th wave vector $k$ for the $s$’th scattering event of the $n$’th photon. The angle of a scattering event can be sampled from the phase function. However, since we are assuming $g=0$ where the scattering is isotropic as described in Step 1, we can simply generate unit vectors ${\overrightarrow{n}}_{\text{out}}$ with a random direction and calculate ${\overrightarrow{n}}_{\text{out}}-{\overrightarrow{n}}_{\text{in}}$ using a coordinate space where ${\overrightarrow{n}}_{\text{in}}=(\mathrm{0,0},1)$. The momentum transfer is then ${\overrightarrow{q}}_{\mathrm{nis}}=k_i*({\overrightarrow{n}}_{\mathrm{out}}-{\overrightarrow{n}}_{\mathrm{in}}{)}_{\mathrm{ns}}$.

• Step 8: At every time point $t$, obtain the contribution of the $n$’th photon and $i$’th $k$ vector to the electric field as

  Eq. (13)

  $$E_{n,i}(t)=\sqrt{S(k)}{e}^{i({\varphi }_n+k_i{L}_n)+\sum _{s}i{\overrightarrow{q}}_{ins}*\overrightarrow{\mathrm{\Delta }}{r}_{ns}(t)}{\sqrt{e^{-{\mu }_a{L}_n}}}^{i({\varphi }_n+k_i{L}_n)-\sum _{s}i{\overrightarrow{q}}_{ins}·\overrightarrow{\mathrm{\Delta }}{r}_{ns}(t)}\sqrt{e^{-{\mu }_a{L}_n}}.$$

It is noteworthy that if ${\mu }_a$ is different in different tissue types, ${\mu }_a{L}_n=\sum _{ti}{\mu }_{ati}{L}_{nti}$, where ${\mu }_{ati}$ and $L_{nti}$ are the absorption coefficient and photon pathlength in the $ti$’th tissue type

• Step 9: Calculate the intensity fluctuations $I(t)$. This is calculated by summing over the contribution from all the photon trajectories and wavelengths

  Eq. (14)

  $$I(t)=\sum _i{|\sum _n{E}_{n,i}(t)|}^2.$$

Here we have ignored the term $e^{i\omega t}$, where $\omega$ is the angular frequency of the light, and we have performed an incoherent sum over the $k$ vectors since typically the measurement time is long enough such that light waves with different angular frequencies or wavelengths do not interfere. An example of the simulated intensity temporal fluctuations $I(t)$ is shown in Fig. 2(a).

• Step 10: Convert intensity $I(t)$ to photon counts $N(t)$. If the average photon flux within a bin time $\mathrm{\Delta }t$ is $N_{\text{bin}}$, we then normalize $I(t)$ such that the average value of ${⟨I(t_n)⟩}_n=N_{\text{bin}}$, where $I(t_n)$ is the intensity within the $n$’th bin time $t_n$. Then we perform a Poisson draw for each bin time as, $N(t_n)=Pois(I(t_n))$. This process converts the continuous function $I(t)$ to discrete values as well as incorporates shot noise into the numerically generated photon counts $N(t)$. Other sources of noise can also be added to $N(t)$, including dark noise, but this is out of the scope of this paper as dark noise is usually negligible in experimental systems. An example of the conversion from $I(t)$ in Fig. 2(a) to $N(t)$ is shown in Fig. 2(b). The Poisson draw here is not unique due to the nature of shot noise.

• Step 11: Calculate the intensity autocorrelation function $g_2(\tau )=⟨N(t)N(t+\tau )⟩/{⟨N(t)⟩}^2$. It is noteworthy that if we skip Step 10, $I(t)$ can also be utilized to calculate the intrinsic noise in $g_2(\tau )=⟨I(t)I(t+\tau )⟩/{⟨I(t)⟩}^2$ that arises from speckle statistics due to a finite measurement time,³² which does not include the effect of shot noise.

For CW-DCS simulations, we skip Steps 2 and 3, and use all the recorded photons to calculate the speckle fluctuations utilizing Eqs. (13) and (14). In Sec. 3, we compare our modeling results with existing theoretical predictions described in Sec. 2.1.

Fig. 2

Example of (a) $I(t)$ and (b) $N(t)$ generated from the Monte Carlo wave model assuming the average photon flux is $100,000\text{photons}/\mathrm{s}$. Here (b) is generated from (a) following the procedure in Step 10.

3.

Results

We first simulate CW-DCS measurements with the measurement configuration described in Fig. 1(b). We compare our results with the theoretical predictions of both $g_2(\tau )$ in Eq. (4) and the noise model in $g_2(\tau )$ and $\sigma (\tau )$ in Eq. (9) as shown in Figs. 3(a) and 3(b). We see that our numerical results are in good agreements with the theoretical predictions. These results are obtained for a single wavelength at 800 nm under the assumption that the coherence length is infinitely long.

Fig. 3

(a) $g_2$ obtained from the Monte Carlo wave model averaged over 500 independent simulations compared with the theoretical expression of $g_2(\tau )=1+(\exp (-2\tau /{\tau }_c))$. (b) $\sigma (\tau )=\mathrm{std}(g_2(\tau )-1)$ obtained from Monte Carlo wave model with 500 independent simulations compared with the theoretical prediction given in Eq. (9). The parameters used are $T_b=1\mu \mathrm{s}$, $T=10\mathrm{ms}$, $\beta =1$, $⟨n⟩=0.1$, and ${\tau }_c\approx 32\mu \mathrm{s}$ obtained from fitting in (a). The total number of photons used in the Monte Carlo simulation is ${10}^8$.

Next, we simulate CW-DCS measurements using multiple wavelengths with central wavelength 800 nm to compare with the theoretical predictions of the reduction of $\beta$ due to loss of coherence described by Eq. (6), by varying the width of the source spectrum (i.e., the coherence length). The pathlength distribution $P(L)$ of all the photons detected at the detector agrees with the theoretical prediction in Eq. (7), as shown in Fig. 4(a). With the same $P(L)$, we have calculated $\beta =g_2(0)-1$ as a function of the coherence length $l_c$ using our numerical model and compared with the results calculated from the theoretical prediction in Eq. (6). In these calculations, we have used $I(t)$ instead of $N(t)$ to obtain $g_2(\tau )$ in Step 11 of Sec. 2.2, since the effect of shot noise is not included in Eq. (6).²⁵

Fig. 4

(a) Pathlength distribution of all the photons that arrive at the detector obtained from the Monte Carlo wave model and the theoretical prediction from Eq. (7). The total number of detected photons is $N_p=34759$. (b) Comparison of the Monte Carlo wave model prediction for the dependence of $\beta$ on coherence length with the analytical model in Eq. (9).

After validating our modeling results with existing theoretical predictions for CW-DCS systems, we now simulate TD-DCS measurements with and without the IRF profile in Fig. 1(d) taken into account, using a coherence length of $l_c=90\mathrm{mm}$ for the input light source. For the case without the effect of the IRF, the IRF profile is assumed to be infinitely narrow, i.e., $t_{n\mathrm{IRF}}=0$, $t_{na}=t_{nL}=L/v$ in Step 2 of Sec. 2.2. The gate delay time dependent pathlength distributions $P(t_s,L)$ for $t_s=0.5$, 1.5, and 2.5 ns and $g_2(\tau )$ with and without the IRF are shown in Figs. 5(a) and 5(b). We see that if the effect of the IRF is ignored, the pathlength distribution $P(t_s,L)$ is solely determined by the gates including gate delay $t_s$ and gate width $t_w$. Compared to Fig. 5(a), the pathlength distributions in Fig. 5(b) are broadened and their shapes are determined by the IRF profile in Fig. 1(d). The corresponding results of $g_2(\tau )$ are shown in Figs. 5(c) and 5(d). Compared to Fig. 5(c), the values of $\beta =g_2(0)-1$ in Fig. 5(d) are smaller due to a broader distribution of pathlengths, as we have expected. To validate the modeling results, we compare the $g_2(\tau )$ curves obtained from the Monte Carlo wave model (solid lines) with those obtained from the analytical expression in Eq. (10) (dashed lines) as shown in Figs. 5(c) and 5(d). We see that the modeling and analytical results are in good agreement. It is noteworthy that for Figs. 5(a) and 5(c), we have assumed a finite coherence length with an infinitely small width of the IRF, which is generally not physically realizable. This numerical analysis serves as a sanity check against established concepts.

Fig. 5

The histogram for the gate delay time dependent pathlength distribution $P(t_s,L)$ for the case when (a) the effect of IRF is not taken into account and (b) the effect of IRF is taken into account. (c) $g_2(\tau )$ obtained from the Monte Carlo wave model (solid lines) and analytical predictions (dashed lines) in Eq. (10) for the case when (c) the effect of IRF is not taken into account and (d) the effect of IRF is taken into account. We have calculated results for $t_s=0.5$, 1.5, 2.5 ns, coherence length $l_c=90\mathrm{mm}$, gate width $t_w=0.2\mathrm{ns}$, source–detector separation $\rho =20\mathrm{mm}$. Other parameters are the same as in Fig. 3.

Next, we perform a noise analysis when only the speckle noise is considered, i.e., in the high photon flux limit. We first obtain $g_2(\tau )=⟨I(t)I(t+\tau )⟩/{⟨I(t)⟩}^2$ and then obtain ${\tau }_c$ from fitting using Eq. (4). The signal is obtained as the average of the ${\tau }_c$ values, which is denoted as ${\tau }_c$ for simplicity. Noise is obtained from the $\mathrm{std}({\tau }_c)$ by running 100 instances of the Monte Carlo wave model simulations using the baseline sample properties. The SNR is defined as $\mathrm{SNR}={\tau }_c/\mathrm{std}({\tau }_c)$. To calculate the CNR for estimating changes in blood flow, we produce an activated state by changing the diffusion coefficient from the baseline value $D_B=1\times {10}^{-6}{\mathrm{mm}}^2/\mathrm{s}$ to $D_B=2.25\times {10}^{-6}{\mathrm{mm}}^2/\mathrm{s}$ in the second layer (Tissue 2) shown in Fig. 1(b). The value of ${\tau }_{c^{\prime }}$, which is the average ${\tau }_c$ at the activated states, is obtained from the average of 10 instances of the wave model. The contrast is $\mathrm{\Delta }{\tau }_c={\tau }_c-{\tau }_{c^{\prime }}$, and the CNR is $\mathrm{CNR}=\mathrm{\Delta }{\tau }_c/\mathrm{std}({\tau }_c)$. Here, we have used a large change of the diffusion coefficient to minimize errors in the contrast induced by noise to save computational time, since this paper is mainly focused on demonstrations of the model instead of calculating the realistic contrast due to brain activation. A smaller change of $D_B$ can be applied, while more configurations used for averaging are required to obtain both ${\tau }_c$ and ${\tau }_{c^{\prime }}$. The comparison of the performance of TD-DCS with and without the effect of IRF, and CW-DCS in the speckle noise limit are shown in Fig. 6. The value ${\tau }_c$ decreases with increasing $t_s$ in Fig. 6(a) as expected, since photons with longer pathlengths that undergo more scattering events are selected at larger $t_s$. The SNR increases with $t_s$ as shown in Fig. 6(b). This is because for a fixed measurement time ($T=10\mathrm{ms}$ in this case), more speckle evolutions are averaged at longer $t_s$ because of the shorter ${\tau }_c$. To quantify the specificity to brain activation in the second layer in Fig. 1(b), we obtain the contrast with activation only in the second layer, i.e., the diffusion coefficient is increased from baseline value $D_B=1\times {10}^{-6}{\mathrm{mm}}^2/\mathrm{s}$ to $D_B=2.25\times {10}^{-6}{\mathrm{mm}}^2/\mathrm{s}$ in the second layer (Tissue 2), and compared with the case when brain activation is applied to the full medium, i.e., the diffusion coefficient is increased from baseline value $D_B=1\times {10}^{-6}{\mathrm{mm}}^2/\mathrm{s}$ to $D_B=2.25\times {10}^{-6}{\mathrm{mm}}^2/\mathrm{s}$ in the whole simulation region with its contrast denoted as $\mathrm{\Delta }{\tau }_{c,h}$. We see in Fig. 6(c) that the specificity to brain activation in the second layer, i.e., $\mathrm{\Delta }{\tau }_c/\mathrm{\Delta }{\tau }_{c,h}$, increases with increasing $t_s$, and that at large $t_s$, $\mathrm{\Delta }{\tau }_c/\mathrm{\Delta }{\tau }_{c,h}$ is larger for TD-DCS than for CW-DCS. This demonstrates the ability of TD-DCS measurements at larger $t_s$ to achieve better sensitivity to tissue dynamics at deeper layers. Thus, TD-DCS can tune the depth specificity by changing the measurement gate delays with a fixed source–detector geometry, which is convenient for technologies such as high-density DCS measurements. Also, CNR is higher for TD-DCS compared to CW-DCS in this high-photon flux, speckle noise limit as shown in Fig. 6(d).

Fig. 6

Comparison of the performance of TD-DCS and CW-DCS in the speckle noise limit. Since CW-DCS results do not vary with $t_s$, they are denoted as dashed lines with constant y values. (a) ${\tau }_c$, (b) SNR, and (c) the ratio of the contrast for activation in the second layer and activation in the full regime $\mathrm{\Delta }{\tau }_c/\mathrm{\Delta }{\tau }_{c,h}$, and (d) CNR for CW-DCS, and TD-DCS with and without the effect of the IRF taken into account. The parameters used here are the same as in Fig. 5.

Now, we consider the more realistic case for human brain measurements where the photon flux is very low and the DCS signal is shot noise limited instead of speckle noise limited. The photon flux difference between TD-DCS and CW-DCS measurements needs to be taken into account. We assume the incident photon flux to be $100,000\text{photons}/\mathrm{s}$ for CW-DCS at $\rho =20\mathrm{mm}$. We then scale the photon flux accordingly by the number of photons detected in a measurement gate in Step 3 for TD-DCS and the total number of the detected photons in Step 1 of the Monte Carlo simulation for CW-DCS, as shown in Table 1. For the low photon fluxes typical of human measurements, long measurement times $T$ are needed to get estimates of $g_2(\tau )$ with sufficient quality that can be fit to obtain an estimate of ${\tau }_c$. Running the wave model for long $T$ is computationally expensive. Further, we want to keep $T$ fixed to be 10 ms to facilitate comparison with prior examples. We thus utilize the following trick to improve computational efficiency that gives us an estimate of the shot noise contribution to $\mathrm{std}({\tau }_c)$ (i.e., $\mathrm{std}{({\tau }_c)}_{\mathrm{shot}}$) that is distinct from the finite speckle sampling noise contribution to $\mathrm{std}({\tau }_c)$ (i.e., $\mathrm{std}{({\tau }_c)}_{\text{speckle}}$). We can then get the total noise from the sum of the variances. Our trick is to sample the wave model for $T=10\mathrm{ms}$ to get $I(t)$. We then obtain $N_{\mathrm{avg}}=1000$ independent Poisson draws of this $I(t)$ get 1000 instances of $N(t)$. We obtain $g_2(\tau )$ for each instance of $N(t)$ and average the 1000 correlation functions to get a more smooth $g_{2,\mathrm{avg}}(\tau )$. We then fit to get ${\tau }_{c,\mathrm{avg}}$ from $g_{2,\mathrm{avg}}(\tau )$. We repeat this process 100 times to obtain $\mathrm{std}({\tau }_{c,\mathrm{avg}})$. Finally, we obtain $\mathrm{std}{({\tau }_c)}_{\text{shot}}=\sqrt{N_{\mathrm{avg}}}*\mathrm{std}({\tau }_{c,\mathrm{avg}})$, since shot noise contribution is proportional to $1/\sqrt{N_{\mathrm{avg}}}$.

With these photon flux, we then calculate the CNR for TD-DCS at $\rho =\mathrm{10,20}\mathrm{mm}$ with the effect of the IRF and coherence length taken into account, and CW-DCS at $\rho =\mathrm{10,20},30\mathrm{mm}$ as shown in Fig. 7(a). We see that for TD-DCS, the maximum CNR for the cases we have calculated occurs at $\rho =10\mathrm{mm}$, $t_s=1.5\mathrm{ns}$, when $t_w$ is fixed to be 0.2 ns, but its CNR is still lower than CW-DCS at $\rho =20\mathrm{mm}$. The best performance for the cases we have simulated is CW-DCS at the largest source–detector separation of $\rho =30\mathrm{mm}$. Thus, when photon flux is taken into account, the relative performance of TD-DCS as compared to its CW-DCS counterpart has degraded. One way to increase the photon flux in TD-DCS is to increase the gate width $t_w$. We see that by increasing $t_w$, the CNR of TD-DCS at $t_w=0.5\mathrm{ns}$ has surpassed the performance of CW-DCS at $\rho =20\mathrm{mm}$. But at large $t_w$, the specificity to a particular photon pathlength $L$ and thus the ability to differentiate tissue depths decreases.

Fig. 7

CNR comparisons of TD-DCS and CW-DCS in the shot noise limit. (a) CNR of TD-DCS as a function of $t_s$ at $t_w=0.2\mathrm{ns}$, $\rho =\mathrm{10,20}\mathrm{mm}$, compared with CW-DCS at $\rho =\mathrm{10,20},30\mathrm{mm}$. (b) CNR of TD-DCS as a function of $t_w$ at $t_s=1.5\mathrm{ns}$, $\rho =10\mathrm{mm}$, compared with CW-DCS at $\rho =20\mathrm{mm}$.

Table 1

The photon flux used for all the simulation cases in Fig. 7.

4.

Discussion

We have developed a Monte Carlo wave model that is capable of simulating both TD-DCS and CW-DCS measurements from first principles. This model can be utilized to calculate the noise in TD-DCS and compare with CW-DCS measurements. We have validated our model with existing analytical expressions of $g_2(\tau )$ for CW-DCS, TD-DCS, and the noise model of CW-DCS. Our numerical model is capable of simulating any spectral profile, while existing theories always assume a Gaussian spectrum.^18,25 Our system is capable of dealing with different brain activation geometries other than the two layer system we have specified in this manuscript. From the comparison of TD-DCS and CW-DCS, we have seen that the performance of TD-DCS is not necessarily better compared to CW-DCS, with the metric given by the CNR induced by tissue dynamics change in the deeper layer of a two layer sample, mainly due to the relatively low photon flux in TD-DCS since only a portion of the detected photons are selected. Methods that can increase the photon flux for a measurement gate, such as temporal focusing, could significantly improve the performance of TD-DCS, as we have seen that in the high photon flux limit, TD-DCS out-performs CW-DCS. The performance of TD-DCS is highly sensitive to the measurement parameters including gate delay time $t_s$, gate width $t_w$, and source–detector separation $\rho$. It also depends on the hardware properties including the wavelength and coherence length of the input laser light, the size of the detector, and the IRF profile. An analysis of the full parameter space is out of the scope of this manuscript. For a give TD-DCS measurement and brain activation geometry, the optimized parameters can be determined using the code of our Monte Carlo wave model. Our model will guide the future experimental design of novel TD-DCS systems.