Measuring neuronal activity with diffuse correlation spectroscopy: a theoretical investigation

Abstract. Significance: Diffuse correlation spectroscopy (DCS) measures cerebral blood flow non-invasively. Variations in blood flow can be used to detect neuronal activities, but its peak has a latency of a few seconds, which is slow for real-time monitoring. Neuronal cells also deform during activation, which, in principle, can be utilized to detect neuronal activity on fast timescales (within 100 ms) using DCS. Aims: We aim to characterize DCS signal variation quantified as the change of the decay time of the speckle intensity autocorrelation function during neuronal activation on both fast (within 100 ms) and slow (100 ms to seconds) timescales. Approach: We extensively modeled the variations in the DCS signal that are expected to arise from neuronal activation using Monte Carlo simulations, including the impacts of neuronal cell motion, vessel wall dilation, and blood flow changes. Results: We found that neuronal cell motion induces a DCS signal variation of ∼10−5. We also estimated the contrast and number of channels required to detect hemodynamic signals at different time delays. Conclusions: From this extensive analysis, we do not expect to detect neuronal cell motion using DCS in the near future based on current technology trends. However, multi-channel DCS will be able to detect hemodynamic response with sub-second latency, which is interesting for brain–computer interfaces.


Introduction
Diffuse correlation spectroscopy (DCS) is an optical imaging method that measures blood flow non-invasively and continuously. It quantifies a blood flow index by measuring the temporal autocorrelation function of the speckle intensity fluctuations of diffusive light remitted from tissue. [1][2][3][4] A change in tissue dynamics results in a change of the decay time of the temporal autocorrelation function. Thus, DCS can be utilized to detect tissue dynamics arising from neuronal activities. The variation of the decay time is often only attributed to a change in cerebral blood flow (CBF). 5,6 The peak of the CBF often occurs at a time delay of a few seconds with respect to the onset of neuronal activation, which is slow and not feasible for real-time monitoring of brain activation in applications such as brain-computer interfaces. *Address all correspondence to Xiaojun Cheng,xcheng17@bu.edu Other mechanisms that cause tissue dynamics due to neuronal activation can also contribute to a change in the DCS signal. Studies have demonstrated that the optical properties of brain tissue vary due to neuronal activation. 7,8 The neural mechanisms that may contribute to the optical signal change may be related to the firing of action potentials. During an action potential, ions are exchanged through the membrane of the neuron. This may cause a change in the neuron shape (e.g., swelling as ions enter the cell), and one hypothesis is that these changes contribute to phase changes of light as it reflects or passes through the shifting cells. [9][10][11][12] This fast optical signal associated with brain activation has been measured in vivo, [13][14][15] but negative results have also been reported. 15,16 Reported delay time of the cell dynamics with respect to the onset of neuronal activation ranges from within 1 to 100 ms, 12,15 and our recent in vivo mouse brain measurements have demonstrated that most of the values of the delay time are within 100 ms, 17 as illustrated in Fig. 1. Apart from changes to the neurons themselves, cells such as pericytes and glial cells may also reshape to support neuronal activity. In addition to direct cellular signals, Fig. 1 Illustration of neurovascular coupling and the onsets of neuronal activation, neuronal cell dynamics, and vessel wall dynamics. Upon brain activation, the neuron generates an action potential (a) that leads to a deformation of the neuron cell body that potentially relates to the exchange of ions through the cell membrane. (b) The onset of this neuronal cell dynamics occurs on the order of 100 ms timescale. Meanwhile, the neuron sends the signal mediated by astrocytes and pericytes to trigger blood vessel wall dilation. (c) The blood vessel wall movement leads to hemodynamic changes in the CBV and CBF. Such neurovascular coupling dynamics contribute to motions that, in principle, could be picked up by the DCS signal. neurovascular coupling causes blood vessels to dilate during neuronal activation, which delivers more oxygen to the excited region. This hemodynamic response causes a change in the phase of light scattered from the vessels, as well as cerebral blood volume (CBV) and CBF variations. The onset of hemodynamic changes has a typical delay time of 450 ms with respect to the onset of neuronal activation, 18 which is a slower timescale compared with fast cellular signals. All of these mechanisms can result in changes in the DCS signal; as such, DCS can potentially provide important measures of real-time neuronal activity non-invasively. A comprehensive analysis of the DCS signal variation induced by neuronal activity will benefit the instrument development of highly sensitive DCS systems to study brain function.
We extensively analyzed the DCS signal variation due to neuronal activation by taking into account all of the above-mentioned mechanisms. Specifically, our model considers the impact of neuronal cell movement, blood vessel wall dilation, and blood flow and volume changes related to neuronal activities. We have found that the DCS signal change induced by neuronal cell dynamics is beyond the sensitivity of currently available DCS systems. However, sub-second detection of neuronal activation utilizing the early behavior of hemodynamicsinduced DCS signal variations is technological feasible with SPAD cameras currently under development. 19 This analysis of the mechanisms that underlie the DCS signal variation is important to the development of highly sensitive DCS systems for various applications, including studies of brain functions, monitoring of brain states at the bedside, and brain-computer interfaces.

Methods
In this section, we describe the use of Monte Carlo simulations to calculate the DCS signal variations induced by neuronal cell dynamics and hemodynamics. Other mechanisms that could contribute to DCS signal changes are also discussed. The calculation of the noise is also presented to obtain the contrast-to-noise ratio (CNR) to quantify the performance of a particular measurement.

Monte Carlo Simulations and the Calculation of DCS Signals
We used Monte Carlo simulations to model photon migration through a semi-infinite 3D dynamic scattering medium. The Monte Carlo code is a derivation of that utilized in previous publications. [20][21][22][23][24] In the simulation, a total number of 10 8 photons were launched on the sample surface and normal to it in the z direction, as illustrated in Fig. 2(a). For a photon that reaches the detector, the pathlength L n and the accumulated dimensionless momentum transfer Y n were recorded. A local region (shaded with orange color) is specified at z ¼ 15 mm, x ¼ 15 mm with respect to the position of the source, with size in x and y of 10 mm and size in z = of 4 mm. A detected photon's total momentum transfer Y n2 and pathlength L n2 within this local region are also recorded. (b) The representative baseline DCS signal obtained from Monte Carlo simulations as compared with the prediction of diffusion theory obtained from Eq. (4). Here the source-detector separation ρ ¼ 30 mm, detector radius r ¼ 2 mm, μ 0 s ¼ 1 mm −1 μ a ¼ 0.01 mm −1 , αD ¼ 10 −6 mm 2 ∕s, and λ ¼ 800 nm.
The remitted photons were collected at a detector placed ρ ¼ 30 mm away from the source and with a radius of r ¼ 2 mm. The reduced scattering coefficient of the medium was set to μ 0 s ¼ 1 mm −1 . A non-zero μ a ¼ 0.01 mm −1 value was incorporated when calculating the temporal field autocorrelation function g 1 ðτÞ ¼ hE Ã ðtÞEðt þ τÞi∕hjEðtÞj 2 i as shown below, where τ is the correlation time and EðtÞ is the measured electric field at the detector. For the n'th photon arriving at the detector, we recorded its total pathlength L n and the accumulated dimensionless momentum transfer Y n . The momentum transfer is the change of the wavevector during a single scattering eventq ¼k out −k in , wherek in andk out are the incident and scattered wavevectors, respectively. Here, we considered only elastic scattering, thus jk in j ¼ jk out j ¼ k 0 . The accumulated dimensionless momentum transfer of a detected photon is the normalized sum of the square of the momentum transfer of all of the scattering events Y n ¼ P s q 2 s ∕2k 2 0 , where s denotes the s'th scattering event of the n'th photon. With Y n and L n obtained from the Monte Carlo simulations, the temporal field autocorrelation function js calculated as E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 1 ; 1 1 6 ; 5 7 1 where C is a normalization factor such that g 1 ð0Þ ¼ 1, N p is the total number of the detected photons, and α is the probability that a scattering event happens at a particular type of scatterer of interest, such as RBCs or neurons. The motion of the scatterers was assumed to be uncorrelated. The functional form of the mean square displacement hΔr 2 i depends on the nature of the dynamics of the scattering particles. For scatterers exhibiting ballistic motion, where v is the speed of the particles, whereas for diffusive motion hΔr 2 ðτÞi ¼ 6Dτ, where D is the diffusion coefficient. The motion of the RBCs is a combination of ballistic flow and shearinduced diffusion. 25 It has been experimentally observed and numerically demonstrated that DCS signals are dominated by the diffusive behavior of RBCs. 21,[25][26][27] Thus, we only modeled the contribution from the diffusive behavior of the moving RBCs in this paper. As described later in Secs. 2.2 and 2.3, the dynamics of the neuronal cell motion and vessel wall dilation are modeled as ballistic motions, for we have not seen existing literature that suggests random diffusive behavior of these movements. The baseline DCS signal in the absence of neuronal activation is expressed as E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 2 ; 1 1 6 ; 3 4 9 expð−2αY n k 2 0 DτÞ expð−μ a L n Þ: Here, we used αD ¼ 10 −6 mm 2 ∕s, which provides a decay time of g 2 close to experimental observations, 3 wavelength λ ¼ 800 nm, and k 0 ¼ 2π∕λ. The representative baseline DCS signal before neuronal activation that was applied is shown in Fig. 2(b). Experimentally, the intensity temporal autocorrelation function g 2 ðτÞ ¼ hIðtÞIðt þ τÞi∕hjIðtÞji 2 is measured instead of g 1 , and g 2 ðτÞ is related to the theoretically modeled g 1 ðτÞ via the Siegert relation: 28 E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 3 ; 1 1 6 ; 2 3 7 with β ¼ 1 indicating complete coherence of the detected photons and β < 1 accounting for loss of coherence and detection of multiple modes of the electromagnetic field. The representative DCS signal computed with Monte Carlo simulations is compared with analytical results obtained from the correlation diffusion equation for a semi-infinite medium: 1,3 E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 4 ; 1 1 6 ; 1 5 7 Here, s . The DCS signals obtained from the Monte Carlo simulation utilizing Eq. (2) and from the theoretical prediction obtained using the correlation diffusion equation [Eq. (4)] are in good agreement as shown in Fig. 2(b).
Neuronal activation can occur in different regions of the brain. Here, we refer to the case of neuronal activation being applied to the full simulation region as global activation. To account for localized neuronal activity within a small region in our model, we specified a second tissue type in the Monte Carlo simulation for a typical local neuronal activation measurement, 29 with the size in x − y of 10 mm and in z of 4 mm and being 15 mm away from the source in the x direction and 15 mm beneath the sample surface as shown in Fig. 2(a) (orange color), which is utilized to calculate the local neuronal activation that we refer to in this paper. The accumulated dimensionless momentum transfer Y n2 and the total pathlength L n2 within this local region were also recorded to calculate the DCS signal change induced by neuronal activation only within this local region. The autocorrelation function is then expressed as the contributions from these two tissue types E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 5 ; 1 1 6 ; 6 0 4 with i ¼ 1 and i ¼ 2 representing the scattering events and photon trajectories outside and within this local region, respectively, and To quantify the DCS signal variation, we need to identify a parameter that characterizes the decay rate of the g 1 ðτÞ curves. Note that, since the medium is semi-infinite, g 1 ðτÞ is no longer a single exponential decay function, as is the case for an infinite medium. The full expression of Eq. (4) can be used to obtain a blood flow index D. However, a simpler function used for fitting is highly preferred for real-time measurements. The exact solution in Eq. (4) for Brownian motion is recast into For the parameters that we are using, τ s ∕τ c ¼ 2.6, and τ c ¼ 46 μs is estimated from fitting. We used the first 70 μs of the g 1 ðτÞ curve to fit the functional form of g 1 ðτÞ ¼ expð−τ∕τ c Þ, which ensures that τ < τ s , to obtain the decay time τ c . Finally, we use the single parameter τ c to characterize the dynamics of the brain tissue.

DCS Signal Change Induced by Neuronal Cell Dynamics
Some studies have demonstrated that action potential propagation induces phase changes of light passing through or reflected from neuronal tissue, corresponding to a membrane displacement on the order of a few nanometers, 12,31-35 which normally happens within 100 ms with respect to the onset of neuronal activation as illustrated in Fig. 1. This phase change is potentially caused by the change of the cell size or the refractive index within the cell. For the purpose of modeling this effect on the DCS signal, we only need to know the effective phase change per unit time, which we discuss in terms of movement of the cell membrane only. We use the average speed of the cell membrane movement of v ¼ 1 nm∕ms, which is consistent with the literature of ex vivo studies of various types of cells 12,31-35 and our recent in vivo measurements of the fast optical signals in the mouse brain using optical coherence tomography. 17 To account for the effect of this neuronal cell motion on the DCS signal, we revised the total temporal field autocorrelation function as E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 6 ; 1 1 6 ; 2 3 7 g 1n;blood ðτÞg 1n;neuronal ðτÞ expð−μ a L n Þ; Here, C is the normalization factor such that g 1 ð0Þ ¼ 1, and g 1n denotes the contribution from a single photon. We considered the neuronal scattering probability α neuronal to be 1, which provides the best-case scenario prediction and g 1n;neuronal ðτÞ ¼ expð− 1 3 Y n;neuronal k 2 0 v 2 τ 2 Þ. In reality, this probability can vary with brain regions where the densities of the neurons differ. [36][37][38] However, as discussed in Sec. 3, the induced DCS signal change for this best-case scenario is beyond the sensitivity of currently available DCS systems. Therefore, we do not delve into the details of the variations of the neuronal cell densities in different brain regions in this paper.

DCS Signal Change Induced by Vessel Wall Dynamics
In addition to neuronal cell movement, blood vessels also undergo dynamics due to neurovascular coupling on a slow timescale, as illustrated in Fig. 1. The increase of vessel diameter increases blood volume and decreases vascular resistance, thus resulting in an increase in blood flow. Each of these effects impacts the DCS signal in the following three ways. First, with the blood volume increase, the probability of scattering from the moving RBCs increases and thus α in Eq. (2) increases. The baseline RBC scattering probability is set to be α ¼ 2%, which is similar to the volume fraction of the vessels, and we let the scattering probability change in direct proportion to the change in CBV during brain activation. Second, the blood flow speed increases, which causes a proportional increase in the diffusion coefficient D in Eq. (2). 21 Third, the vessel wall movement results in a phase change of light scattered from the vessel wall. The average speed of the vessel wall movement v vessel can be obtained from the time course of the vessel diameter change as demonstrated below.
We incorporated the phase changes induced by vessel wall movement (the third effect) into the formalism the same way as the phase changes induced by neuronal cell motion by replacing g 1n;neuronal with g 1n;vessel ðτÞ ¼ expð− 1 3 αY n k 2 0 v 2 vessel τ 2 Þ in Eq. (6). We found that the magnitude of this (third) effect is at least 5 orders of magnitude smaller compared with that of the blood volume (first) and the flow speed changes (second), and the value of τ c remains the same for the precision used; thus it is ignored for the rest of the work presented here. We only consider the variation in g 1n;blood ðτÞ ¼ expð− 1 3 αY n;blood k 2 0 Ã 6DτÞ, and g 1 ðτÞ ¼ C N p P N p n¼1 g 1n;blood ðτÞ expð−μ a L n Þ. To obtain the changes of CBV and CBF, we utilized a single compartment vascular model that considers the full vascular network to be one compartment to estimate the time courses of rCBVðtÞ and rCBFðtÞ due to neuronal activation. The letter r denotes the relative value normalized to the baseline value. For example, rCBFðtÞ ¼ CBFðtÞ∕CBFð0Þ. Here, t ¼ 0 denotes the onset of vessel wall dilation. Note that there is a time delay between the onset of vessel wall dilation and the onset of neuronal activation, which is estimated to be about 450 ms, 39 whereas the time delay for neuronal cell dynamics is typically within 100 ms, 17 as illustrated in Fig. 1. Thus, when the hemodynamics-induced DCS signal change is considered, the definition of the fast timescale is within 100 ms with respect to the onset of vessel wall dilation or within 450 to 550 ms with respect to the onset of neuronal activation. The functional form of the relative change of the vessel diameter d is given as an input to the model 39 as follows: E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 7 ; 1 1 6 ; 2 8 3 The constant parameters are set to be Δd ¼ 0.07 and σ d ¼ 1.83 s, following the recommended values in Refs. 39 and 40. The time course of the rCBVðtÞ is related to rdðtÞ via rCBVðtÞ ¼ rdðtÞ 2 , due to the cylindrical geometry of the vessels and given that the total length of the vessel is invariant. When the total pressure across the compartment is fixed, rCBFðtÞ ¼ rdðtÞ 4 from Poiseuille's law. The illustration of the model and the time courses of rd, rCBV, and rCBF are shown in Fig. 3. This is a simpler approximation of the blood flow dynamics following neuronal activation as compared with the multi-compartment model 40 more commonly used to model fMRI 41,42 and fNIRS 18,43 measurements. Depending on the hemodynamics model used, the peak value and time of the rCBFðtÞ and rCBVðtÞ can vary due to the extra parameter of the delay time between the peak of CBV and CBF in the multi-compartment models. However, for the early time behavior that does not depend on the delay time we are mostly interested in, this simplified model is sufficient, and the magnitude of the results remain the same as compared with using a more sophisticated model. We explored the DCS signal change due to the variation of α and D within the measurement time T when we start the measurement at time t d after the onset of vessel wall dilation. Note that CBF is the volumetric flow that includes both the effects of the (1) volume and (2) flow speed changes. For the effect of the variation of α and D within T, we consider that, when we start the measurement at a time delay t d , the dynamics can be characterized by a time varying diffusion coefficient Dðt d ; τÞ and scattering probability αðt d ; τÞ with αðt d ; τÞDðt d ; τÞ ¼ αðt d ÞDðt d Þð1 þ R rCBF τÞ, where the rate of change is defined as R rCBF ¼ ΔrCBF∕T, T is the measurement time window, and ΔrCBF is the change of rCBF within T. Strictly speaking, the exact formalism for (1) needs to be calculated from integration for a D value that varies with time. Here, we used a linear approximation as above since the change is small, which can be well described by the first-order term in the Taylor expansion. The product αD at the time delay t d is αðt d ÞDðt d Þ ¼ αð0ÞDð0Þ Ã rCBFðt d Þ. The contribution from a single photon is then g 1n;blood ðτÞ ¼ expð− 1 3 αðt d ; τÞY n;blood k 2 0 Ã 6Dðt d ; τÞτÞ. We see that the effect of the variation of D and α within the measurement time T to the DCS signal change is a few orders of magnitude smaller compared with the effect of the delay time t d for the range of t d values that we explored, i.e., t d from 50 ms to 2 s. This range of t d is chosen since it provides a DCS signal variation that can be potentially detected using multi-channel DCS systems currently under development. Thus, we have assumed that αðt d ; τÞDðt d ; τÞ ¼ αðt d ÞDðt d Þ for a non-zero t d value used in this paper, leading to g 1n;blood ðτ; t d Þ ¼ expð− 1 3 αð0ÞY n;blood k 2 0 Ã 6Dð0Þ Ã rCBFðt d ÞτÞ. Local vessel wall dynamics induced by local neuronal activation were also considered, where the increasing values of α and D and the speed of vessel wall dilation were only applied to the i ¼ 2 component in Eq. (5).

Noise Model for DCS Measurements
To determine whether a DCS signal change can be detected experimentally, we compared the DCS signal change induced by neuronal cell movement or hemodynamics with the noise level of the measurement. The noise in g 2 ðτÞ for DCS has been analytically calculated. 30 The standard deviation of ðg 2 ðτÞ − 1ÞÞ, σðτÞ, at each correlation time τ is E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 8 ; 1 1 6 ; 6 6 8 Here, T b is the bin time, m is the bin index, T is the measurement time window, and  Fig. 4. One way to reduce the noise level is to use multiple channels or multiple instances for averaging. When N c channels are used, the noise level is reduced to σ N ðτÞ ¼ σðτÞ∕ ffiffiffiffiffiffi N c p from the central limit theorem. Comparison between the noise level and the neuronal activity-induced DCS signal variation determines the CNR of a measurement. However, as discussed, instead of the change of g 2 ðτÞ, which is a function of the correlation time τ, a single variable τ c and its variation are used to characterize the dynamics of the tissue. Thus, the relation between the theoretically calculated σðτÞ and the noise-induced variation of τ c , stdðτ c Þ, is desired. To calculate stdðτ c Þ, we numerically generated noisy g 2 ðτÞ curves by adding random fluctuations drawn from a Gaussian distribution with mean zero and standard deviation σðτÞ at each τ value. An example of the g 2 ðτÞ curves before and after adding noise is shown in Fig. 5(a), using the parameters that give σðτÞ as in Fig. 4. The distribution of the values of τ c is obtained from the fitting of the noisy g 2 ðτÞ curves, which provides the estimation of stdðτ c Þ∕τ c induced by noise, as shown in Fig. 5(b). For this particular example stdðτ c Þ∕τ c ¼ 13 μs∕46 μs ¼ 0.28. Thus, this system can only resolve activities that induce a signal change of Δτ c ∕τ c > 0.28. Increasing N c reduces the noise level and thus decreases stdðτ c Þ. We define the CNR of the measurement as CNR ¼ Δτ c ∕stdðτ c Þ ¼ C Ã SNR, where the relative contrast C ¼ Δτ c ∕τ c and the signal-to-noise ratio, SNR ¼ τ c ∕ stdðτ c Þ. To estimate the N c required to achieve a CNR ∼1 for neuronal cell movement, vessel wall dilation, and hemodynamics as discussed in Sec. 3, we need to obtain the relation between stdðτ c Þ and N c . Using σ N ðτÞ to generate noisy g 2 ðτÞ curves, we calculate stdðτ c Þ as a function of N c as shown in Fig. 6. We see that 1∕stdðτ c Þ is proportional to ffiffiffiffiffiffi N c p . Thus, CNR ¼ Δτ c ∕stdðτ c Þ . This relation determines the required N c to detect an activation, with CNR ¼ 1.

Results
In this section, we demonstrate the results of the DCS signal variations induced by neuronal cell motion and hemodynamic changes arising from neuronal activations. We first show the DCS signal change that arises on the slow timescale of a few seconds due to the hemodynamic changes, which is what has been typically measured. As an example, we obtained g 2 ðτÞ at baseline and at the peak of the rCBFðtÞ response shown in Fig. 3. The results of the induced DCS signal change arising from global and local changes in blood flow are shown in Fig. 7. The decay times τ c obtained from fitting are indicated in the legends. The induced fractional change of τ c , i.e., Δτ c ∕τ c , is 0.25 for global and 0.01 for local flow and volume changes. For a given measurement time window T ¼ 10 ms, the number of instances/channels required for averaging is N c ¼ 2 and N c ¼ 3805 for global and local activation, respectively, to achieve a CNR ¼ 1 for the particular set of measurement parameters as used in the noise model in Fig. 4. This can feasibly be detected in real time with the state-of-the-art SPAD cameras. 19 We can also increase  On the fast timescale within 100 ms with respect to neuronal activation, the changes in the DCS signal arise from neuronal cell motion. The results using the speed of the cell membrane movement of 1 nm∕ms are shown in Fig. 8. The induced fractional change Δτ c ∕τ c is on the order of 10 −5 for global and 10 −7 for local activation. Note that we have assumed that the best-case scenario is all of the scattering events globally or locally happening at a neuronal cell, i.e., we assume α ¼ 1. In reality, α ≤ 1, so in reality the induced DCS signal change could be smaller than what we predicted here. Unlike hemodynamics-induced changes that last for a few seconds, the signal induced by neuronal cell motion only lasts for <100 ms. 17 Thus, increasing the measurement time is not feasible for the detection of neuronal cell motion in real-time measurements. The number of instances/channels required for averaging to achieve CNR ∼1 is ∼10 9 and ∼10 13 for global and local neuronal activation, respectively, for the particular set of measurement parameters as used in the noise model in Fig. 4. This does not seem to be achievable in the near future based on current detector technology trends. 19 The DCS signal changes induced by fast scale hemodynamics for t d ¼ 0 (450 ms with respect to the onset of neuronal activation, see Sec. 2.3) were also obtained. Compared with the results of slow signals shown in Fig. 7 in which the measurement took place at the peak of the rCBFðtÞ curve, we now discuss the results if the measurement took place right after the for t d ¼ 0 is roughly 1 to 2 orders of magnitude smaller than the effect of neuronal cell on fast timescales, and it cannot be detected by any foreseeable DCS system. We also calculated the DCS signal variation on slow timescales induced by hemodynamics to estimate the earliest time delay t d with respect to vessel wall dilation that a neuronal activation is measurable for a particular DCS system. As opposed to Fig. 7 where the DCS signal variation is obtained at the peak of the hemodynamic changes, we obtain Δτ c ∕τ c measured at different time delays t d with respect to the onset of the vessel wall dilation, i.e., with a measurement time window spanning between t ¼ t d and t ¼ t d þ T, as shown in Fig. 10(a). For example, at a time delay of 100 ms, Δτ c ∕τ c reaches 10 −2 to 10 −3 for global activation and 10 −3 to 10 −4 for local activation. We ignored the DCS signal change induced by changes of D and α within the measurement T, which is at least 4 orders of magnitude smaller. We also estimated the number of channels N c required to reach CNR ∼1 as shown in Fig. 10(b) for the noise level in Fig. 4.
For the current technology that utilizes a kilopixel SPAD camera, for example, it is possible to reach an SNR ¼ τ c ∕stdðτ c Þ of the order of 300. 19 To detect the hemodynamic signal, we need Δτ c ∕τ c > stdðτ c Þ∕τ c ¼ 1∕SNR. From Fig. 10(a), the earliest t d that an activation can be detected is at t d ¼ 100 ms for global activation or t d ¼ 700 ms for local activation. Adding the 450 ms, which is the delay time of the onset of vessel wall dilation with respect to neuronal activation ( Fig. 1), the earliest time that a neuronal activation is detectable is 550 ms (global) and 1.15 s (local) with respect to the onset of the activation. An even earlier detection of hemodynamic signal can be potentially reached by adopting a megapixel SPAD camera for DCS detection in the future. 44

Discussion
We have extensively analyzed the mechanisms that contribute to DCS signals during neuronal activation, including neuronal cell movement and vessel wall movement, though traditionally the variation of the DCS signal is often attributed only to the change of the blood flow, with the diffusion coefficient calculated using Eq. (4). In addition to neurons, other cells such as pericytes and glial cells also reshape due to neuronal activity. Pericytes are contractile cells that act to control the size of the capillaries. 45,46 This effect was already incorporated in the hemodynamics as described in Sec. 2.3 since the capillary size change is essentially covered by our modeling of the vascular diameter changes. The glial cells also swell during neuronal activation. However, as the literature suggests, the timescale of the glial cell motion due to the opening of aquaporins is on the order of a few minutes. 9,[47][48][49] This is much slower compared with the effects of neuronal cell motion in response to membrane potential or vessel wall dynamics. We therefore did not model the contributions from pericytes and glial cells dynamics to the DCS signal. We also obtained the DCS signal variation on both fast and slow timescales. On the fast timescale of 10 to 100 ms after neuronal activation, the largest DCS signal change Δτ c ∕τ c induced by neuronal cell motion is on the order of or smaller than 10 −5 . This is beyond the sensitivity of any foreseeable DCS system based on current detector technologies due to the large number of channels required to achieve a super high SNR.
A point worth mentioning is that Δτ c ∕τ c is smaller than we had expected. As we can see in Fig. 11, the decay time induced solely by neuronal cell motion is 38 ms, which gives a decay rate 10 −3 that of the decay rate of the baseline DCS signal, but the overall induced Δτ c ∕τ c is 10 −5 . This is due to the difference in the functional forms of hΔr 2 ðτÞi in Eq. (1), which is hΔr 2 ðτÞi ¼ v 2 τ 2 for the ballistic motion of neuronal cell motion and hΔr 2 ðτÞi ¼ 6Dτ for the diffusive behavior of blood flow. We tested if imposing a diffusive behavior of the neuronal cell motion with the same decay time as in Fig. 11 would induce Δτ c ∕τ c ∼ 10 −3 as expected. This indicates that, if there exists diffusive-like motion associated with neuronal cell motion similar to that of blood flow, the DCS signal change will be 2 to 3 orders of magnitude larger than what was predicted in Fig. 8. To the best of our knowledge, diffusive-like motion has not been reported for neuronal cell motion during activation. The reported time courses of the phase change of light passing through the cells due to activation does not seem to suggest a diffusive behavior. 12 The change of the DCS signal arising from hemodynamic responses due to neurovascular coupling is another mechanism that can be utilized to detect neuronal activation at a reasonably short latency. We need to add another ∼450 ms to t d in Fig. 10, if the time delay with respect to the onset of the neuronal activation is of interest. This is a slower process compared with cellular motion, but it provides an opportunity to detect sub-second neuronal activation in real time with multi-channel DCS systems.
We analyzed DCS signal changes induced by global and local neuronal activation. The global activation that we refer to here does not have to be an activation that occurs for the whole brain. We define global activation for an activation region large enough to cover the span of most photon trajectories for the particular source-detector separation. For the local activation, we specified a particular small activation region and detection geometry. Apparently, the signal change is different for various geometries, which we did not explore here. We expect the order of magnitude to be similar to what is reported in this work. However, for further exploration, it is straightforward to adjust the parameters in the Monte Carlo simulations to mimic a particular experimental condition. The CNR analysis also depends on the parameters used in the noise model, which can be adjusted accordingly.
We utilized the decay time τ c as a single parameter to characterize brain dynamics in real time. One limitation for this technique is that it requires the single-photon detector to have high time resolution that is capable of resolving the decay time. In addition to τ c , other parameters such as β, μ s , and μ a 40 in Eq. (4) and the spatial contrast as measured in laser speckle contrast imaging 50 may also be utilized to detect neuronal activation. We will compare the performance of using DCS to detect brain activities with other techniques in the future.
In summary, we have extensively analyzed the components of the DCS signal variation during neuronal activation using Monte Carlo simulations. This study has enhanced our fundamental understanding of the underlying mechanisms of brain tissue dynamics that contribute to the DCS signal variations. Our results also provide guidance for the instrument development of DCS systems to detect brain activation with a given latency, which may be relevant for applications such as understanding brain functions, therapeutic monitoring of blood flow, and non-invasive brain-computer interfaces.

Disclosures
Dr. Boas consulted with Facebook Inc. on topics related to diffuse correlation spectroscopy that ultimately led to the ideas presented in this paper. Dr. Boas's interests were reviewed and are managed by Boston University in accordance with their conflict of interest policies. of the journal Neurophotonics published by SPIE. He was awarded the Britton Chance Award in Biomedical Optics in 2016.
Francesco Marsili is a research scientist manager at FRL developing diffuse optical techniques for the measurement of biosignals. Previously, he led research in single photon detectors, optical communication, quantum optics, superconductivity, and nanofabrication at the Jet Propulsion Laboratory, the National Institute of Standards and Technology, and MIT. He received his PhD in physics from the École Polytechnique Fédérale de Lausanne, Switzerland, in 2009.