Dynamic light scattering and laser speckle contrast imaging of the brain: theory of the spatial and temporal statistics of speckle pattern evolution

Dynamic light scattering (DLS) and laser speckle contrast imaging (LSCI) are closely related techniques that exploit the statistics of speckle patterns, which can be utilized to measure cerebral blood flow (CBF). Conventionally, the temporal speckle intensity auto-correlation function g2t(τ) is calculated in DLS, while the spatial speckle contrast K s is calculated in LSCI measurements. Due to the rapid development of CMOS detection technology with increased camera frame rates while still maintaining a large number of pixels, the ensemble or spatial average of g2s(τ) as well as the temporal contrast K t can be easily calculated and utilized to quantify CBF. Although many models have been established, a proper summary is still lacking to fully characterize DLS and LSCI measurements for spatial and temporal statistics, laser coherence properties, various motion types, etc. As a result, there are many instances where theoretical models are misused. For instance, mathematical formulas derived in the diffusive regime or for ergodic systems are sometimes applied to small animal brain measurements, e.g., mice brains, where the assumptions are not valid. Therefore, we aim to provide a review of the speckle theory for both DLS and LSCI measurements with detailed derivations from first principles, taking into account non-ergodicity, spatial and temporal statistics of speckles, scatterer motion types, and laser coherence properties. From these calculations, we elaborate on the differences between spatial and temporal averaging for DLS and LSCI measurements that are typically ignored but can result in inaccurate measurements of blood flow, particularly the spatially varying nature of the static component in g2t(τ) and K t . We also obtained g2s(τ) maps in in vivo mouse brain measurements using high frame rate CMOS cameras which have not been demonstrated before, and compared with g2t(τ) and Ks,t. This work provides a useful guide for choosing the correct model to analyze spatial and temporal speckle statistics in in-vivo DLS and LSCI measurements.


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][2][3][4][5][6][7][8], as well as blood flow in other organs including skin [9][10][11], and retina [12][13][14][15][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][18][19][20][21][22][23], as illustrated in as shown in Fig. 1, and they have significantly impacted, for instance, the studies of cortical spreading depression [24][25][26] and ischemic stroke [27][28][29][30].In DLS, sometimes referred to as DLS imaging (DLSI) [31] for animal brain imaging, or diffuse correlation spectroscopy (DCS) [32][33][34][35] for human brain measurements, the temporal intensity auto-correlation function g t 2 (τ) is calculated, and the dynamics of the scattering sample or blood flow is quantified as the decorrelation time constant τ c of the auto-correlation functions [31,[36][37][38][39].To temporally resolve g t 2 (τ), high measurement rate detectors, with the time between adjacent frames much less than τ 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][41][42][43][44][45][46][47].Unlike DLS, LSCI does not require high frame rate detectors since T needs to be larger than τ c for the blurring of speckle patterns to occur.Multi-exposure LSCI that measures K s versus T [29,[48][49][50], and a recent model that utilizes a generalized spatial contrast K(τ) [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 t 2 (τ) measured in DLS instead of the spatially averaged auto-correlation function g s 2 (τ).The static component in g t 2 (τ) 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 s 2 (τ) 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 t 2 (τ), g s 2 (τ), 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.

DLS theory
In DLS, the speckle field auto-correlation function is defined as 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 (τ) and g 2 (τ) is described by the Siegert relation [53][54][55][56][57] Here β is the coherence parameter taking into account the loss of coherence, with β<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 β [52,[58][59][60].Here ⟨• • • ⟩ 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.
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][62][63][64][65][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 (τ) can be calculated from its definition, The fluctuating component is ergodic since the spatial and temporal averaging are the same, i.e., ⟨I f ⟩ s = ⟨I f ⟩ t .Therefore the subscripts t and s are omitted for all the f components.Here is the fluctuating component of g 1 (τ) that contains information about CBF.
The difference between g s 1 and g t 1 is that ⟨I c ⟩ 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 ⟨I c ⟩ 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 (τ) requires a phase-resolved technique as has been done in DLS-OCT [67][68][69].For most techniques that measure light intensity, g 2 (τ) can be calculated to quantify CBF.To derive g 2 (τ), we first consider the case for fully coherent waves, where the intensity with a static component present is Here θ (i.e.θ(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 t,s 2 or the un-normalized intensity auto-correlation function G t,s 2 (τ) is defined as Again, we have omitted the t, s subscripts for the fluctuating component due to ergodicity, and θ(τ) is the phase difference between E f (t + τ) and E c .The difference between temporal and spatial averaging lies in the first term ⟨I 2 c ⟩ t,s where since I c is constant over time, while spatially it varies and follows the statistics of speckle patterns with a probability distribution P(I c ) ∼ exp(−I c /⟨I c ⟩ s ), with which ⟨I 2 c ⟩ s can be calculated from ∫ P(I c )I 2 c dI c [18,70].The second term 2⟨I c I f ⟩ t,s can be written as with the assumption that the static and dynamic components are not correlated.The third term by definition is where g 2,f (τ) is the intensity auto-correlation function for the fluctuating component that contains information of CBF.The last term can be written as where we have used cos the phase of the speckle pattern is randomized and thus the electric field averages to zero; , and plug Eqs. ( 5)-( 8) into Eq.( 4), Here ρ t,s = ⟨I f ⟩ ⟨I ⟩ t,s , and the Siegert relation in Eq. ( 1) with β = 1 for the fluctuating term is utilized.At large τ, g t 2 approximately goes to 1 and g s 2 goes to 1 ] 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 (τ) 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 here E coh,j is a coherent electric field for a particular light path j with pathlength denoted L j , and Φ j (t) accounts for the loss of coherence for each light path.It can be calculated that [72] Thus, Eq. ( 8) can be written as where E c and E f both have phase fluctuations as described by Eq. ( 12).Plug Eqs. ( 1), ( 6), ( 7), (11), (14) into Eq.( 4), 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, √ β instead of β is applied in the last term of g t 2 (τ) considering I c not fluctuating over time, i.e., no exp[iΦ 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 (τ) 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 (τ) can be calculated and written as [12,31,33,43,77,78] g where τ c is the decorrelation time constant, and n depends on the motion type of the scatterers.For multiple scattering from unordered motion (MU), n = 0.5; for multiple scattering from ordered motion (MO) or single scattering from unordered motion (SU), n = 1; for single scattering from ordered motion (SO), n = 2.It should be noticed that the interpretation of τ c also depends on the scattering regime and particle motion.Unordered (Brownian) motion is characterized by the diffusing coefficient D and ordered motion is characterized by the mean square velocity v rms .In reality, g 1,f (τ) may also be a linear or nonlinear combination of all these components depending on the scattering regime.More discussions can be found in [12,22,31,43,79].For simplicity, we focus on the case of n = 1 for the rest of this manuscript.All the other cases can be derived similarly.Theoretical calculations also assume that the camera exposure time in DLS measurements is much smaller than the decorrelation time constant τ c .However, it might not be true experimentally and the impact of the exposure time in DLS measurements can be incorporated as an extra reduction of β [31].The measurements of g s 2 (τ) have only been done in phantom samples using single photon correlators [61,80].Nowadays the advancement of CMOS detectors enables measurements of the spatio-temporal profiles of speckle patterns in vivo.
The temporal averaging of g s 2 (τ), i.e., ⟨g s 2 (τ)⟩ t , can be used to improve the signal-to-noise ratio.However, as described in Supplement 1, ⟨g t 2 (τ)⟩ s will change the functional form of g t 2 (τ), mainly due to the non-linear averaging of the spatially varying ρ t .

LSCI theory
With the definition of the temporal and spatial g t,s 2 (τ), 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 = ∫ T 0 I(t)dt.If we consider the n = 1 model of g 1,f (τ) = e −τ/τ c , integrating the above equations give where x = T/τ c .The analytical expressions of K s and K t have been described in the literature, but the subtle differences between ρ t and ρ s are often ignored.For instance, the difference between (K s ) 2 and (K t ) 2 has been used to obtain the static component (1 − ρ s ) 2 , but it is only true when a pixel size is much larger than the speckle size so enough speckle averaging is performed at a single pixel level such that ρ t =ρ s .But in most LSCI experiments, this is not true and The spatially varying ρ t needs to be taken into account for fitting to obtain an accurate estimation of τ c when K t is measured.Similar to ⟨g t,s 2 ⟩ s,t described above, ⟨K s ⟩ t can be used to improve the signal-to-noise ratio at the expense of reduced measurement rate, but ⟨(K t ) 2 ⟩ s changes the functional form of Eq. ( 21) due to the spatially varying ρ t .In particular, to improve the signal-to-noise ratio in the estimates of τ c using multiple K t measurements, one can estimate τ c for each K t measurements and then average over the τ c s.We have included a detailed discussion in Supplement 1.

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][82][83][84][85][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 , where E f 1 (t) and E f 2 (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 (τ) function and K(T) model with the slow dynamics component can be obtained (see detailed derivations of g 2 (τ) in Supplement 1): Here τ c2 is a second decorrelation time constant describing the slow dynamics and ρ 2 is the associated fraction of light that has been scattered from the slow-varying scatterers; In the diffusive regime, both the fast decay function g 1,f 1 (τ) and slow decay function g 1,f 2 (τ) can be modeled as exponential decay function, i.e., g ) 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].

Numerical validations
We utilized numerical simulations to validate the equations for g t,s 2 (τ) 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 µm 1000 µm 1000 µm] acting as the scattering medium.The decorrelation time constant is set to be 10 ms and the wavelength of light λ is set to be 0.785 µm.The diffusion coefficient in the model is therefore derived from where k is the wavenumber of light: 2π/λ.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 µm away from the center of the volume with a sensor size of 17.4 µm by 17.4 µm and a pixel size of 0.435 µ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.
To simulate a non-ergodic system with ρ = 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 β = 0.5, we add the electric fields from two independently generated speckle patterns both with ρ = 0.5.g s 2 (τ) and K s (T) are computed from 40 by 40 pixel speckle patterns over a duration of 2 s and g t 2 (τ) 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 ρ and τ c with β fixed at 1 and 0.5 respectively.In practice, β is usually calibrated by the spatial contrast measurement of a static phantom since crosstalks often exist between β and τ 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 ρ t and ρ s since ρ 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 τ 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 (τ) and K(T) with the slow dynamics component are illustrated in Fig. 3. To reduce the crosstalk among the fitting parameters, we fix both β and τ c1 and fit the resting parameters with models described in Eqs. ( 22)- (25).Experimentally this  15), ( 16), (20) and (21).can be done by obtaining β and τ c1 at short time scales, then obtain τ c2 at longer time scales [86].Similarly, ρ t 1 and ρ t 2 are not the same as ρ s 1 and ρ s 2 but the fitting results of ρ t 1 and ρ t 2 are consistent with the results obtained from g t (τ) and K t (T) models since they are evaluated at the same spatial position.

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 µ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 µs.The number of active pixels is reduced to 512x384 to achieve this imaging speed, leading to ≈3×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 (τ) 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 × 15 pixels to obtain the spatial g s 2 (τ) function.The time lag τ ranges from 0 to 0.8 ms.The same dataset used in the g 2 (τ) 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 × 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 × 15 are applied to both g t 2 (τ) 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 (τ) 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 (τ) and g s 2 (τ) at time lags of τ = 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 (τ) and K s (T) than that in g t 2 (τ) and K t (T).The time courses of g s 2 (τ) and K s (T) averaged from a 15×15 window (orange box in Fig. 4(a1) and (b1)) in the tissue region and the time courses of g t 2 (τ) 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 (τ) and K(T) respectively as shown in Fig. 4(a2) and (b2).Here τ c is measured to be 0.24 ± 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.

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 (τ) 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 (τ) 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 t,s 2 (τ), 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 (∼ s), there will not be enough averaging to obtain τ c2 and slow dynamics component appears to be the same as the static component, thus the single τ c model described in Sections 2, 3 should

Fig. 1 .
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 τ to obtain g s 2 (τ); temporal correlation is calculated on a time series of intensity fluctuations for individual pixels (blue rectangle) to obtain g t 2 (τ).(a2) Examples of g t 2 (τ) and g s 2 (τ) calculated from the same non-ergodic system with β = 1, ρ = 0.5, τ c = 30 µ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 β = 1, ρ = 0.5, τ c = 30 µs.Examples of time-evolving speckle patterns are generated from dynamics speckle modeling described in 5.