1. Introduction

Dynamic light scattering (DLS) and laser speckle contrast imaging (LSCI) are two complementary label-free optical techniques that can be utilized to study cerebral blood flow (CBF) of small animal brains such as mice [1–8], as well as blood flow in other organs including skin [9–11], and retina [12–16]. Both techniques exploit the speckle pattern statistics arising from the interference of coherently back-scattered light since the evolution of the speckle patterns is directly related to the dynamics of the scattering samples [8,17–23], as illustrated in as shown in Fig. 1, and they have significantly impacted, for instance, the studies of cortical spreading depression [24–26] and ischemic stroke [27–30]. In DLS, sometimes referred to as DLS imaging (DLSI) [31] for animal brain imaging, or diffuse correlation spectroscopy (DCS) [32–35] for human brain measurements, the temporal intensity auto-correlation function $g_2^t(\tau )$ is calculated, and the dynamics of the scattering sample or blood flow is quantified as the decorrelation time constant $\tau _c$ of the auto-correlation functions [31,36–39]. To temporally resolve $g_2^t(\tau )$, high measurement rate detectors, with the time between adjacent frames much less than $\tau _c$, are needed. Meanwhile, in LSCI, the spatial contrast $K^s$ or the temporal contrast $K^t$ is calculated, and the reduction of $K^s$ or $K^t$ within a certain camera exposure time $T$ reflects the degree of speckle pattern blurring induced by tissue dynamics which can be utilized to quantify CBF variations [1,29,40–47]. Unlike DLS, LSCI does not require high frame rate detectors since $T$ needs to be larger than $\tau _c$ for the blurring of speckle patterns to occur. Multi-exposure LSCI that measures $K^s$ versus $T$ [29,48–50], and a recent model that utilizes a generalized spatial contrast $K(\tau )$ [51] have been developed to quantify CBF using a conventional CMOS camera. However, the subtle distinction between the spatial and temporal statistics of speckle patterns is often ignored. For instance, the spatial contrast $K^s(T)$ measured in LSCI is often described to be directly related to the temporally averaged auto-correlation function $g_2^t(\tau )$ measured in DLS instead of the spatially averaged auto-correlation function $g_2^s(\tau )$. The static component in $g_2^t(\tau )$ and $K^t$ is spatially varying but not emphasized, which can result in inaccurate measurements of blood flow if not taken into account [31]. Additionally, formulas derived for diffusive and ergodic systems are sometimes applied to animal brain measurements where these assumptions may not be valid. A proper guide is desired to bridge the gap between fundamental studies of speckle pattern evolution to the biomedical applications of measuring CBF and/or other tissue dynamics from speckle fluctuations.

Our goal is to provide a comprehensive guide for choosing the model for DLS and LSCI measurements, taking into account the effects of non-ergodicity, the spatial and temporal statistics of speckles, scatterer motion types, and laser coherence properties. Nowadays, the spatial intensity auto-correlation function $g_2^s(\tau )$ as well as the temporal contrast $K^t$ can also be conveniently measured, thanks to the rapid development of CMOS detector technologies that enable imaging at high frame rates while maintaining a large number of pixels. We summarize the theoretical models and the derivations for both spatial and temporal DLS and LSCI measurements, i.e. analytical expressions for $g_2^t(\tau )$, $g_2^s(\tau )$, $K^s(T)$, $K^t(T)$, from first-principles. We validate the analytical results using our recently developed dynamic scattering model (DSM) that simulates the evolution of speckle patterns [52]. We then provide examples of temporal and spatial DLS and LSCI measurements in in-vivo mouse brain at a frame rate of 30 kHz and compared them with the analytical expressions.

2. DLS theory

In DLS, the speckle field auto-correlation function is defined as $g_1(\tau ) = \frac {\langle E^*(t)E(t + \tau )\rangle }{\langle |E(t)|^2\rangle }$ and the intensity auto-correlation function is $g_2(\tau ) = \frac {\langle I(t)I(t + \tau )\rangle }{\langle I(t)\rangle ^2}$. Here $E$ is the speckle electric field and $I=|E(t)|^2$ is the measured speckle intensity. For ergodic systems, i.e. systems with spatially and temporally sampled speckle intensities following the same statistics, the relation between $g_1(\tau )$ and $g_2(\tau )$ is described by the Siegert relation [53–57]

Here $\beta$ is the coherence parameter taking into account the loss of coherence, with $\beta <1$ accounting for the loss of coherence and/or detection of multiple modes of light. The effect of the speckle-to-pixel size ratio, i.e., when the pixel size does not meet the Nyquist criterion to resolve the speckles, can also be incorporated into $\beta$ [52,58–60]. Here $\langle \cdots \rangle$ is used to denote averaging, but the distinction between spatial and temporal averaging for non-ergodic systems is not emphasized, and temporal averaging is typically used in DLS measurements by default. Throughout this manuscript, we use $\langle \cdots \rangle _s$, $\langle \cdots \rangle _t$ to denote spatial and temporal averaging respectively.

For many biomedical applications such as in-vivo mouse brain imaging, there is often a fraction of photons that have not experienced scattering from red blood cells, which is referred to as a static component that breaks down ergodicity [61–66]. When a static component is present, the speckle electric field can be described as $E(t)=E_f (t)+E_c$, where the subscripts $f$ and $c$ represent fluctuating and static/constant terms respectively, and thus, the temporal and spatial $g^{t,s}_1 (\tau )$ can be calculated from its definition,

The fluctuating component is ergodic since the spatial and temporal averaging are the same, i.e., $\langle I_f\rangle _s=\langle I_f\rangle _t$. Therefore the subscripts $t$ and $s$ are omitted for all the $f$ components. Here $g_{1f}(\tau )=\frac {\langle E_f^*(t)E_f(t + \tau )\rangle }{\langle |E_f(t)|^2\rangle }$ is the fluctuating component of $g_1(\tau )$ that contains information about CBF. The difference between $g_{1}^s$ and $g_{1}^t$ is that $\langle I_c\rangle _{t}$ is constant over time but is spatially varying, since $I_c$ is a static speckle pattern that contains bright and dark intensity regions due to constructive and destructive interference, while $\langle I_c\rangle _{s}$ provides the average intensity of the static component which is nominally constant over time and space, assuming the window size over which the spatial averaging occurs is large enough.

The measurement of $g_1(\tau )$ requires a phase-resolved technique as has been done in DLS-OCT [67–69]. For most techniques that measure light intensity, $g_2(\tau )$ can be calculated to quantify CBF. To derive $g_2(\tau )$, we first consider the case for fully coherent waves, where the intensity with a static component present is

Here $\theta$ (i.e. $\theta (t)$ with $t$ omited for simplicity) is the phase difference between $E_f(t)$ and $E_c$. $I_f$ is a simplified representation of $I_f(t)$. The numerator of $g_{2}^{t,s}$ or the un-normalized intensity auto-correlation function $G_{2}^{t,s}(\tau )$ is defined as

Again, we have omitted the $t,s$ subscripts for the fluctuating component due to ergodicity, and $\theta (\tau )$ is the phase difference between $E_f(t+\tau )$ and $E_c$. The difference between temporal and spatial averaging lies in the first term $\langle I_c^2\rangle _{t,s}$ where

Here $\rho _{t,s}=\frac {\langle I_f\rangle }{\langle I\rangle _{t,s}}$, and the Siegert relation in Eq. (1) with $\beta =1$ for the fluctuating term is utilized. At large $\tau$, $g_2^t$ approximately goes to $1$ and $g_2^s$ goes to $1+(1-\rho _s)^2$. Note that $real[g_{1,f}]$ is kept here but in many cases it is simplified to $g_1$, since for un-correlated phase shifts $real[g_1]=g_1$, which is true for multiple scattered light or un-ordered motion of scatterers encountered in most DLS measurements. The subscript $f$ is usually ignored in the literature, but we keep it here to distinguish it from $g_1(\tau )$ in Eq. (2).

When there is loss of coherence, Eq. (5) can be rewritten as [18]

For the last term in Eq. (4), the reduction of coherence can be modeled as phase fluctuations [71,72] where

Thus, Eq. (8) can be written as

Here we consider the loss of coherence to be the same for both the static and dynamic components. These assumptions might not be true in some experimental cases and need to be tested. For instance, if a large fraction of $E_c$ comes from the reflection of the optical components in the system, the speckle statistics can be impacted since the specular reflections do not generate speckles and need to be corrected. Also, in some cases, $\sqrt {\beta }$ instead of $\beta$ is applied in the last term of $g_2^t(\tau )$ considering $I_c$ not fluctuating over time, i.e., no $exp[i\Phi _j(t)]$ is applied to $E_c$ components. However, we consider the case that the loss of coherence comes from the temporal fluctuations of the laser light [73,74] that apply to both $I_c$ and $I_f$. We will not explore the details of the differences between the spatial and temporal coherence of a light source. But if there exists any difference, the analytical expressions become more complex and we recommend using numerical models to simulate the laser coherence properties by adding spatially or temporally correlated phase fluctuations in the light source [52,75].

Since the introduction of the diffusing wave spectroscopy, the field autocorrelation function $g_{1,f}(\tau )$ has been modeled as an integral transform of the probability density function of the interfering waves in the scattering medium [54,76]. Assuming the number of scattering events and the types of motions, the analytical expression for $g_{1,f}(\tau )$ can be calculated and written as [12,31,33,43,77,78]

The temporal averaging of $g^s_2(\tau )$, i.e., $\langle g^s_2(\tau ) \rangle _t$, can be used to improve the signal-to-noise ratio. However, as described in Supplement 1, $\langle g^t_2(\tau ) \rangle _s$ will change the functional form of $g^t_2(\tau )$, mainly due to the non-linear averaging of the spatially varying $\rho _t$.

3. LSCI theory

With the definition of the temporal and spatial $g_{2}^{t,s}(\tau )$, it is straightforward to construct the temporal and spatial contrast from the integration [1,12,36]

And

Here $I_T$ is the integrated intensity with a certain camera exposure time $T$ that can be calculated from $I_T=\int _0^T I(t)dt$. If we consider the $n=1$ model of $g_{1,f}(\tau )=e^{-\tau /\tau _c}$, integrating the above equations give

Similar to $\langle g_2^{t,s}\rangle _{s,t}$ described above, $\langle K^s\rangle _t$ can be used to improve the signal-to-noise ratio at the expense of reduced measurement rate, but $\langle (K^t)^2\rangle _s$ changes the functional form of Eq. (21) due to the spatially varying $\rho _t$. In particular, to improve the signal-to-noise ratio in the estimates of $\tau _c$ using multiple $K^t$ measurements, one can estimate $\tau _c$ for each $K^t$ measurements and then average over the $\tau _c$s. We have included a detailed discussion in Supplement 1.

4. Slow dynamics

DLS and LSCI can also be utilized to measure slower dynamics other than blood flow, such as intracellular motility with much slower speed than red blood cells with a decorrelation time on the order of 100 ms to a few seconds [81–86]. Considering the slow dynamics component as well as the static component present in the speckle measurement, the speckle electric field can be described as $E(t) = E_{f1}(t) + E_{f2}(t) + E_c$, where $E_{f1}(t)$ and $E_{f2}(t)$ represent fast and slow varying light fields respectively. Following similar steps of derivation as indicated in Sections 2 and 3, the spatial and temporal $g_2(\tau )$ function and $K(T)$ model with the slow dynamics component can be obtained (see detailed derivations of $g_2(\tau )$ in Supplement 1):

Here $\tau _{c2}$ is a second decorrelation time constant describing the slow dynamics and $\rho _2$ is the associated fraction of light that has been scattered from the slow-varying scatterers; $x_1=\frac {T}{\tau _{c1}}$, $x_2=\frac {T}{\tau _{c2}}$, $x_{12}=T(\frac {1}{\tau _{c1}}+\frac {1}{\tau _{c2}})$. In the diffusive regime, both the fast decay function $g_{1,f1}(\tau )$ and slow decay function $g_{1,f2}(\tau )$ can be modeled as exponential decay function, i.e., $g_{1,f1}(\tau )= e^{-(\frac {\tau }{\tau _{c1}})^n}$, $g_{1,f2}(\tau )= e^{-(\frac {\tau }{\tau _{c2}})^n}$, with n = 0.5, 1 or 2 depending on the motion types of the scatters as discussed in Section 2. More dynamic components can be incorporated similarly into the models if needed. However, experimentally we have only observed two dynamic components so far [86].

5. Numerical validations

We utilized numerical simulations to validate the equations for $g_2^{t,s}(\tau )$ and $K^{t,s}(T)$ in a non-ergodic system derived from first principles, as shown in Fig. 2. We use our recently developed dynamic speckle model to simulate the temporal evolution of the speckle pattern with details described in [23,52]. In short, we use $N_p = 5000$ particles acting as point sources emitting spherical waves that are moving in time. The particles are randomly distributed in a volume of size [1000 $\mu m$ 1000 $\mu m$ 1000 $\mu m$] acting as the scattering medium. The decorrelation time constant is set to be 10 ms and the wavelength of light $\lambda$ is set to be 0.785 $\mu m$. The diffusion coefficient in the model is therefore derived from $D = 1/\sqrt {\tau _c k^2}$, where k is the wavenumber of light: $2\pi /\lambda$. We use 1 ms as the simulation time step and the particle positions are updated at each time step with the standard Monte-Carlo method of Gaussian random walks. The camera is located at 3000 $\mu m$ away from the center of the volume with a sensor size of 17.4 $\mu m$ by 17.4 $\mu m$ and a pixel size of 0.435 $\mu m$. Four data points are sampled at each speckle to ensure enough sampling of speckles from the Nyquist theorem. The speckle pattern calculated from the superposition of the spherical waves as a function of time is recorded at the sensor.

 

Fig. 1. Illustration of DLS and LSCI measurements. (a) DLS theory: (a1) Examples of speckle patterns captured at high frame rate; spatial correlation within a window (red rectangle) is calculated between frames with a time lag $\tau$ to obtain $g^s_2(\tau )$; temporal correlation is calculated on a time series of intensity fluctuations for individual pixels (blue rectangle) to obtain $g^t_2(\tau )$. (a2) Examples of $g^t_2(\tau )$ and $g^s_2(\tau )$ calculated from the same non-ergodic system with $\beta = 1, \rho = 0.5, \tau _c = 30\,\mu s$. (b) Multi-exposure LSCI theory: (b1) speckle patterns are captured with different exposure time; spatial contrast is calculated on a spatial window (red rectangle) and temporal contrast is calculated on a time series of intensity fluctuations for a single pixel (blue rectangle). (b2) Examples of $K^t(T)$ and $K^s(T)$ calculated from the same non-ergodic system with $\beta = 1, \rho = 0.5, \tau _c = 30 \,\mu s$. Examples of time-evolving speckle patterns are generated from dynamics speckle modeling described in 5.

Download Full Size | PDF

 

Fig. 2. Numerical validation: (a) $g_2(\tau )$ and $K(T)$ for the case of $\beta = 1, \rho = 0.5, \tau _c = 10$ ms. $\beta$ is fixed at 1. Fitting results: with $g^t_2(\tau )$ model, $\rho _t = 0.59, \tau _c = 10$ ms; with $g^s_2(\tau )$ model, $\rho _s = 0.49, \tau _c = 10$ ms; with $K^t(T)$ model, $\rho _t = 0.60, \tau _c = 10$ ms; with $K^s(T)$ model, $\rho _s = 0.50, \tau _c = 10$ ms; (b) $g_2(\tau )$ and $K(T)$ for the case of $\beta = 0.5, \rho = 0.5, \tau _c = 10$ ms. $\beta$ is fixed at 0.5. Fitting results: with $g^t_2(\tau )$ model, $\rho _t = 0.78, \tau _c = 9.8$ ms; with $g^s_2(\tau )$ model, $\rho _s = 0.48, \tau _c = 9.8$ ms; with $K^t(T)$ model, $\rho _t = 0.80, \tau _c = 10$ ms; with $K^s(T)$ model, $\rho _s = 0.49, \tau _c = 9.8$ ms. In both (a) and (b), blue curves represent temporal $g^t_2(\tau )$ or $K^t(T)$; red curves represent spatial $g^s_2(\tau )$ or $K^s(T)$; dot curves represent simulation results obtained from dynamic speckle model; solid curves represent fitted results with the models obtained from Eqs. (15), (16), (20) and (21).

Download Full Size | PDF

To simulate a non-ergodic system with $\rho = 0.5$, we first generate a static speckle pattern and then add it to a time-evolving speckle pattern as the dynamic scattering component described above with the same average intensity as the static speckle pattern. To simulate $\beta =0.5$, we add the electric fields from two independently generated speckle patterns both with $\rho =0.5$. $g^s_2(\tau )$ and $K^s(T)$ are computed from 40 by 40 pixel speckle patterns over a duration of 2 s and $g^t_2(\tau )$ and $K^t(T)$ are calculated from speckle patterns observed for 200 s at a single pixel position to ensure enough averaging of spatial and temporal statistics of speckles. We fit $\rho$ and $\tau _c$ with $\beta$ fixed at 1 and 0.5 respectively. In practice, $\beta$ is usually calibrated by the spatial contrast measurement of a static phantom since crosstalks often exist between $\beta$ and $\tau _c$ when fitting [74,87]. We see that the simulated results are in good agreement with the analytical expressions described in Eqs. (15), (16), (20) and (21) as shown in Fig. 2. It should be noticed that there is a difference between $\rho ^t$ and $\rho ^s$ since $\rho ^t$ is spatially varying and depends on the intensity distribution of the static speckle pattern as discussed in Section 3.

To simulate the slow dynamics, another dynamic speckle pattern with $\tau _2 = 1 \, s$ is added to the existing speckle pattern described in Fig. 2(a) with the same average intensity as the fast-dynamics speckle pattern as well as the static speckle pattern. The simulation and fitting results for the spatial and temporal $g_2(\tau )$ and $K(T)$ with the slow dynamics component are illustrated in Fig. 3. To reduce the crosstalk among the fitting parameters, we fix both $\beta$ and $\tau _{c1}$ and fit the resting parameters with models described in Eqs. (22)–(25). Experimentally this can be done by obtaining $\beta$ and $\tau _{c1}$ at short time scales, then obtain $\tau _{c2}$ at longer time scales [86]. Similarly, $\rho ^t_1$ and $\rho ^t_2$ are not the same as $\rho ^s_1$ and $\rho ^s_2$ but the fitting results of $\rho ^t_1$ and $\rho ^t_2$ are consistent with the results obtained from $g^t(\tau )$ and $K^t(T)$ models since they are evaluated at the same spatial position.

 

Fig. 3. Numerical simulation for slow dynamics: $g_2(\tau )$ and $K(T)$ for the case of $\beta = 1, \rho _1 = 0.33, \rho _2 = 0.33, \tau _{c1} = 10$ ms, $\tau _{c2} = 1$ s. In both (a) and (b), blue curves represent temporal $g^t_2(\tau )$ or $K^t(T)$; red curves represent spatial $g^s_2(\tau )$ or $K^s(T)$; dotted curves represent simulation results obtained from dynamic speckle model; solid curves represent fitted results with the models obtained from Eqs. (22)–(25). For all the fittings, $\beta$ is fixed at 1 and $\tau _{c1}$ is fixed at 10 ms. Fitting results: with $g^t_2(\tau )$ model, $\rho _{1t} = 0.41, \rho _{2t} = 0.41, \tau _{c2} = 1.0$ s; with $g^s_2(\tau )$ model, $\rho _{1s} = 0.34, \rho _{2s} = 0.36, \tau _{c2} = 1.1$ s; with $K^t(T)$ model, $\rho _{1t} = 0.41, \rho _{2t} = 0.43, \tau _{c2} = 1.1$ s; with $K^s(T)$ model, $\rho _{1s} = 0.35, \rho _{2s} = 0.31, \tau _{c2} = 0.95$ s;

Download Full Size | PDF

6. Experimental results

We have performed fast speckle pattern measurements on the mouse brain cortex for in-vivo experimental demonstration. A high-speed camera (Photron Nova S6, 20x20 $\mu m^2$ pixels) is used to record the backscattered light, through a 6.3x objective (Leica UVI 6.3x, NA = 0.13), tube lens (TTL200-B) and a long-pass filter (cut on at 750 nm). The recordings are performed at a framerate of 30 kHz at an exposure time of 33 $\mu s$. The number of active pixels is reduced to 512x384 to achieve this imaging speed, leading to $\approx$3$\times$2.25 mm$^2$ field of view. A fiber-coupled volume holographic grating stabilized laser diode (785 nm, Thorlabs FPV785P) is mounted on the arm perpendicular to the imaging axis, operated at the recommended settings. The coherent light from the laser diode is delivered on the sample co-axially through the polarising beamsplitter positioned behind the objective. The backscattered light collected by the objective passes through the polarizing cube, which acts as a cross-polarizer to increase the contrast.

To obtain the temporal $g^t_2(\tau )$ function, the intensity auto-correlation function is calculated from a time series of data of 7.5 ms (containing 225 time points) for each pixel. The spatial or ensemble averaging is conducted for each frame on a window of 15 $\times$ 15 pixels to obtain the spatial $g^s_2(\tau )$ function. The time lag $\tau$ ranges from 0 to 0.8 ms. The same dataset used in the $g_2(\tau )$ measurement is also utilized to obtain the spatial and temporal contrast. Similarly, the temporal contrast $K^t(T)$ is calculated from 225 time points for each pixel and the spatial contrast $K^s(T)$ is calculated from a 15 $\times$ 15 window for each frame. By applying frame binning, we mimic the camera intensity measurements at multiple exposure times with $T$ varying from 0.033 ms to 3 ms. To match the spatial resolutions of temporal measurements and spatial measurements, 2D mean filters of 15 $\times$ 15 are applied to both $g^t_2(\tau )$ and $K^t(T)$ for visualization purposes. But as described in Sections 2, 3, the spatially smoothed data cannot be used for fitting. The mean values of $g^s_2(\tau )$ and $K^s(T)$ are calculated from 225 time points to enhance the signal-to-noise ratio of the spatial measurements. Fig. 4(a1) shows the mouse cortex spatial maps of $g^t_2(\tau )$ and $g^s_2(\tau )$ at time lags of $\tau$ = 0.1 and 0.8 ms respectively. Fig. 4(b1) shows the spatial maps of $K^t(T)$ and $K^s(T)$ at $T$ = 1 and 3 ms respectively. Higher residual static components are observed in $g^s_2(\tau )$ and $K^s(T)$ than that in $g^t_2(\tau )$ and $K^t(T)$. The time courses of $g^s_2(\tau )$ and $K^s(T)$ averaged from a 15$\times$15 window (orange box in Fig. 4(a1) and (b1)) in the tissue region and the time courses of $g^t_2(\tau )$ and $K^t(T)$ at a single pixel within the indicated window are demonstrated in Fig. 4(a2) and (b2). We use Eqs. (15), (16), (21) and (20) to fit the spatial and temporal $g_2(\tau )$ and $K(T)$ respectively as shown in Fig. 4(a2) and (b2). Here $\tau _c$ is measured to be $0.24 \pm 0.015$ ms, indicating consistent fast dynamics obtained by fitting with different models. Experimental measurement of slow dynamics requires optical systems to enhance the sensitivity to the static and slow components compared to wide field measurements. We have demonstrated a short source-detector separation speckle contrast optical spectroscopy system to measure the slow tissue dynamics in our previous work [86] at a spatial point, and will construct multiple source-detector pairs to obtain images of the slow dynamics in the future.

 

Fig. 4. In-vivo mouse cortex measurements: (a) $g^{t,s}_2(\tau )$ measurement: (a1) spatial maps of $g^t_2(\tau )$ and $g^s_2(\tau )$ at the specific time lags of $\tau = 0.1$ and 0.8 ms; (a2) Examples of $g_{2}(\tau )$ curves for a specific region (orange box in a1). (b) $K^{t,s}(T)$ measurement: (b1) spatial maps of $K^t(T)$ and $K^s(T)$ at the specific exposure times of $T = 1$ and 3 ms; (b2) Examples of $K^{t,s}(T)$ curves for a specific region (orange box in b1). In both (a2) and (b2), blue curves represent temporal $g^t_2(\tau )$ or $K^t(T)$; red curves represent spatial $g^s_2(\tau )$ or $K^s(T)$; dot curves represent in-vivo measurements. solid curves represent fitting results obtained from Eqs. (15), (16), (20) and (21); the white line represents 200 $\mu m$. For all the fittings, $\beta$ is fixed at 0.24. Fitting results: with $g^t_2(\tau )$ model, $\rho _t = 0.59, \tau _c = 0.21$ ms; with $g^s_2(\tau )$ model, $\rho _s = 0.61, \tau _c = 0.24$ ms; with $K^t(T)$ model, $\rho _t = 0.67, \tau _c = 0.23$ ms, with $K^s(T)$ model, $\rho _s = 0.64, \tau _c = 0.29$ ms.

Download Full Size | PDF

7. Conclusion

We have demonstrated the derivation from first-principles of the models for DLS and LSCI measurements and shown experimentally the calculations of $g^{t,s}_2(\tau )$ and $K^{t,s}(T)$ from high-speed speckle imaging measurements. To the best of our knowledge, our work is the first experimental demonstration of wide-field imaging of $g^s_2(\tau )$ in vivo. Note that there are other factors in the experiments that we have not fully explored, for instance, the noise level. The uncorrelated noise could reduce the value of the first point of $g_2^{t,s}(\tau )$, as well as add a bias term to $K^{t,s}(T)$. For mouse brain measurements, the photon flux is relatively high thus the noise contributions are ignored. But they can be modeled as needed [52,88]. Also, for temporal and spatial averaging, we are assuming enough samples, i.e., an infinite number of speckles to be used for averaging. Experimentally, this might not be true. For instance, when measuring the slow dynamics as described in Section 4, if the measurement time $t$ is not much larger than the slow dynamics ($\sim s$), there will not be enough averaging to obtain $\tau _{c2}$ and slow dynamics component appears to be the same as the static component, thus the single $\tau _c$ model described in Sections 2, 3 should be used (See Supplementary Fig. S2 for an example). The quantification of noise due to finite spatial samples and its corrections have been analyzed in [89–91]. We found that $g^s_2(\tau )$ is more robust to model fitting than $g^t_2(\tau )$ due to the higher value at larger time lags $\tau$ that contains information of the static and slow dynamics components. The spatial $g^s_2(\tau )$ is capable to resolve the slow dynamics even with insufficient measurement time as shown in Supplementary Fig. S3, suggesting the potential advantage of analyzing the spatial statistics ($g^s_2(\tau )$ or $K^s(T)$) over the temporal statistics ($g^t_2(\tau )$ or $K^t(T)$) to measure the slow dynamics. Furthermore, this paper presents the speckle contrast model $K(T)$ based on the $MO/SU_{n=1}$ functional form of $g_{1,f}(\tau )$. It is worth noting that speckle contrast models, derived for different scattering regimes and motion types such as $MU_{n=0.5}$ and $SO_{n=2}$ as detailed in [87], can become highly intricate, especially when slow dynamics is taken into account. Therefore we suggest implementing the numerical integral of $g_{2}(\tau )$ as shown in Eqs. (18) and (19) to obtain the corresponding $K^{t,s}(T)$ rather than explicitly deriving the analytical expressions of $K^{t,s}(T)$ for more complex motion types. We have made our model fitting code publicly available with all the motion types included [92]. In theory, extracting information from $g_2(\tau )$ and $K(T)$ is considered equivalent, given that the $K(T)$ model is an integral of $g_2(\tau )$. In practice, analyzing $g_2(\tau )$ can be more straightforward because of the less complex expression compared to the $K(T)$ models but requires high temporal sampling rate to resolve $\tau _c$.

In summary, in this paper we have reviewed the derivations of the models for DLS and LSCI measurements from first-principles. Our work elucidates the one-to-one relationship between $g^{t,s}_2(\tau )$ and $K^{t,s}(T)$ and emphasizes the distinctions in temporal and spatial speckle statistics, particularly in the presence of non-ergodicity. Our work offers theoretical support for the selection of accurate analytical formulas for in vivo DLS and LSCI measurements.