Choosing a model for laser speckle contrast imaging

: Laser speckle contrast imaging (LSCI) is a real-time full-field non-invasive technique, which is broadly applied to visualize blood flow in biomedical applications. In its foundation is the link between the speckle contrast and dynamics of light scattering particles–erythrocytes. The mathematical form describing this relationship, which is critical for accurate blood flow estimation, depends on the sample’s light-scattering properties. However, in biological applications, these properties are often unknown, thus requiring assumptions to be made to perform LSCI analysis. Here, we review the most critical assumptions in the LSCI theory and simulate how they affect blood flow estimation accuracy. We show that the most commonly applied model can severely underestimate the flow change, particularly when imaging brain parenchyma or other capillary perfused tissue (e.g. skin) under ischemic conditions. Based on these observations and guided by the recent experimental results, we propose an alternative model that allows measuring blood flow changes with higher accuracy.


Introduction
Laser speckle contrast imaging (LSCI) is a non-invasive optical imaging technique that can provide information on light scattering particles' dynamics. Thanks to the high spatio-temporal resolution and implementation simplicity, LSCI became a key imaging tool to monitor perfusion in the brain [1][2][3], skin [4][5][6], retina [7,8] and other organs. LSCI analyzes speckle pattern arising from random interference of coherent light scattered by the light-scattering particles, such as red blood cells. Movement of the particles results in intensity fluctuations, which occur as blurred speckle pattern when recorded by a camera with a finite exposure time. Pattern blurring is directly related to the particles' dynamics, and, thus, to the blood flow, and can be quantified as speckle contrast K. The form of the contrast-to-blood flow relationship depends on the characteristics of the light scattering process in the sample, such as scattering regime (single or multiple), particle motion type (ordered or unordered) and presence of static scattering. These characteristics are described in the models of the intensity and field temporal autocorrelation functionsg 2 (τ) and g 1 (τ). Under the conventional LSCI theory assumptions, the contrast-to-blood flow relationship becomes 1/K 2 ∝ BFI, where BFI is the blood flow index. This relationship permits real-time data processing because of its numerical simplicity, and thus has been adopted in the majority of LSCI applications [8][9][10][11][12].
However, there is growing evidence that using the relationship as described above results in a discrepancy in the blood flow estimation [13,14]. The discrepancy becomes particularly notable when imaging ischemic stroke -a major LSCI application [1,12,14,15], where conventional LSCI severely underestimates the decrease in blood flow [14]. The reason for the discrepancy is an inaccuracy of assumptions in the contrast-to-blood flow relationship derivation. The key assumptions are related to: (i) the form of the field correlation function g 1 (τ) which is defined by the scattering regime (single or multiple) and the motion type of the particles (ordered or unordered) [16,17]; (ii) the presence of static scattering [18,19]; and (iii) impact of coherence, polarization and speckle averaging effects on the correlation loss which is reflected in the parameter β [20].
The form of the field correlation function was historically assumed to be g 1 (τ) = exp(−τ/τ c ), where τ is the time lag, and τ c is the decorrelation time -a quantitative measure of dynamics. It is the central assumption that was necessary to build the LSCI theory and to study the effects of static scattering [15,18]. However, this form only describes the field correlation for the single scattering unordered motion or multiple scattering ordered motion regimes. It becomes invalid for the cases where single scattering ordered motion or multiple scattering unordered motion regimes are predominant [14,17,21]. Recent research shows that these regimes are unexpectedly abundant in the typical LSCI applications where discrepancies were observed [14,17,22], making it clear that the fundamental assumptions have to be re-visited.
In this study, we derive new LSCI models that correspond to different forms of the field temporal autocorrelation function in order to obtain a more accurate contrast-to-blood flow relationship. We use the derived models to numerically analyze how the key assumptions affect blood flow estimation accuracy. We re-evaluate the effects of static scattering and the impact of the correlations loss. We support our numerical investigation with a comparison between dynamic light scattering imaging (DLSI) [14] and LSCI measurements of the blood flow in the mouse model of ischemic stroke.

Conventional LSCI model
The mathematical relationship between spatial contrast K s and decorrelation time τ c was first derived by Fercher and Briers in 1981 [16] and revised in later studies [23,24]. Briefly, K s is given by where <I> and σ s are the mean and standard deviation of intensity calculated over a pixel's neighbourhood. The spatial variance σ 2 s and, thus, speckle contrast, can then be related to the intensity autocorrelation function g 2 (τ) with the intensity covariance C t (τ): where T is the exposure time, τ is the time lag. There are no assumptions related to the light-scattering process in these equations. However, assumptions are required to relate g 2 (τ) to the field temporal autocorrelation function g 1 (τ) and thus to particles' dynamics. In the original derivation the light scattering medium was assumed to be ergodic, thus allowing the use of the Siegert relation [16]: where β accounts for the correlation loss due to the speckle averaging, polarization and low stability and coherence of the light source. Most importantly, except the most recent studies [14,17], the form of the field correlation function g 1 (τ) is generally defined as: which corresponds to the single light scattering from the unordered motion of scattering particles or to multiple light scattering from the ordered motion of scattering particles [16,17]. Combining Eqs. (1)(2)(3)(4)(5), results in the equation that relates spatial speckle contrast to the decorrelation time: Equation (6) can be solved to find the decorrelation time and get the blood flow estimate in the form of the inverse decorrelation time 1 τ c . However, solving the equation is rarely done in practice as it is computationally expensive and requires pre-calibration of the parameter β. Instead, to allow real-time processing and ease adoption of the technique, the equation was further simplified, and the blood flow index (BFI) was introduced as an estimate of perfusion: Equation (7) relies on the assumption that T ≫ τ c [25], which is true for most applications [14,18,26] and renders the Conveniently, since β depends on the system configuration and does not change within a single experiment, Eq. (8) allows β-independent estimation of the relative blood flow (rBF): where baseline and response are defined either in time (before vs after stimulus) or in space (affected tissue vs non-affected tissue). Calculating rBF with Eq. (9) has become the most common way of analyzing LSCI data, as it provides a system-independent (not affected by β) metric of the relative blood flow change. It is applied in a broad range of biomedical studies, including stroke [9,12,27,28], sensory stimulation [29,30], drug response [31,32], and vasoreactivity research [8,32].

Key assumptions
Assumptions on the form of the field autocorrelation function, ergodicity and the decorrelation time were necessary to derive Eqs. (4-9). Below we discuss each assumption in detail: • Assumption 1: form of the field correlation function is g 1 (τ) = exp(−τ/τ c ). In conventional LSCI, Eq. (5) is used to describe the form of the field autocorrelation function g 1 (τ) and its relation to the decorrelation time. This form is defined by the dynamics scattering regime and corresponds to the single scattering unordered motion, or multiple scattering ordered motion. Recent studies, however, show that these regimes change depending on the imaged tissue and the vessel size [14,17]. Specifically, when imaging capillary perfused tissue (e.g. skin or parenchymal regions of the brain), the observed dynamics are better described with multiple light scattering and unordered motion . In contrast, for vessels larger than 150 µm in diameter (e.g. human retinal vessels), a single scattering and ordered motion is expected. To account for this difference Eq. (6) has to be re-derived using a corresponding form of the g 1 (τ). See Table 1 for a summary of the forms of g 1 (τ).
• Assumption 2: static scattering is absent. The original derivation assumed ergodicity and thus an absence of static scattering in the sample. The question of how static scattering affects the contrast analysis was later addressed with multi-exposure laser speckle imaging (MESI) theory [18]. In MESI, the modified form of the intensity correlation function was introduced: where ρ = (6) is then re-derived as [18]: MESI showed that static scattering is mostly absent when imaging medium or large vessels directly, but must be accounted for when imaging through skull or monitoring perfusion in parenchymal regions of the brain [15]. However, the static scattering contribution was analysed based on the MO/SU n=1 form of the g 1 (τ), and thus requires re-evaluation, particularly in parenchymal regions where MU n=0.5 form is dominating [14].
• Assumption 3: decorrelation time τ c is much shorter than exposure time T. This assumption is critical in the derivation of Eqs. (8) and (9), as it allows discarding the second-order term in Eq. (6). It is true in many of the LSCI applications, as typically used optimal T is on the order of 5ms [26], while τ c is ≤1ms. However, in some applications, such as ischemic stroke [27] or systemic sclerosis [33], flow speed can be reduced by ≈10 times, thus increasing τ c and invalidating the assumption [14]. To keep the assumption valid in such cases, one may consider using a longer exposure time. It, however, will lead to the increased relative contribution of static scattering, making it even more critical to address the assumption 2. Furthermore, increasing T will reduce sensitivity to changes in faster flows [26] and impose stricter frame rate limitations, making it overall impractical in biomedical applications.
• Assumption 4: β has no effect on rBF. This assumption holds in two cases: (i) β ≈ 1 or (ii) static scattering is absent and τ c <<T. The value of β is determined by the experimental setup [20,34,35] and is generally in the range ≈0.3 to 0.6 [14,35], making the first case implausible. The second case, as discussed above, is only true for some of the LSCI applications. For other applications, the deviation of β from 1 should not be ignored, and its contribution to rBF measurements accuracy has to be evaluated. Table 1. Form of the field correlation function g 1 (τ) for different scattering regimes and particle motion types

Newly derived contrast models
To explore how the assumptions affect measurement accuracy, we follow the procedure described above and re-derive Eq. (11) for the different forms of g 1 (τ). This section only presents the final equations, while the detailed derivation can be found in the Supplement 1. Using the multiple scattering with unordered motion form MU n=0.5 of the field correlation function we arrive at: Single scattering with ordered motion SO n=2 results in: where erf (x) is the error function defined by Here and below we refer to the derived Eqs. (12) and (13) as the MU n=0. 5 and SO n=2 models, and the original Eq. 11 and its simplified version Eq. 7 as the MO/SU n=1 and MO/SU n=1,simp models respectively. Figure 1 shows the contrast versus T τ c in the absence of static scattering ( Fig. 1(A)) and with 20% static scattering ( Fig. 1(B)). Even in absence of static scattering, there is a noticeable difference between the MO/SU n=1 and the most commonly used MO/SU n=1,simp models, which becomes visible at T τ c ≈ 10 and grows larger for smaller T τ c , eventually leading to the contrast exceeding the theoretical maximum of 1 for the MO/SU n=1,simp model. For non-ergodic case, the MO/SU n=1,simp model deviates from the detailed models even more, as it still reaches K ≈ 0 at larger T τ c , highlighting the importance of the assumption 2 in its derivation. which is common in parenchyma or skin imaging [14,15] . There is no correlation loss due to pixel averaging or polarization effects (β = 1).

Simulation
To quantify LSCI measurement accuracy, we simulate relative blood flow estimation with different models and calculate the relative error caused by invalid assumptions. The process can be described as follows: • Step 1: Define the baseline decorrelation time τ cb and the relative blood flow rBF true . Following the existing data [14,18], we set τ cb to 50µs, 200µs, and 1ms for the large vessels, medium vessels, and parenchyma, respectively. For the rBF true , we chose a range of values from 0.1 to 3 (10% to 300% of the baseline flow) to cover the flow change from ischemia to hyperperfusion.
• Step 3: Calculate the speckle contrast that corresponds to τ cb and τ cr using the model defined as ground truth. The ground truth model (MU n=0.5 , MO/SU n=1 or SO n=2 ) and its parameters (ρ,β) are defined accordingly to the vessel size and to the tested assumptions.
• Step 4: Convert contrast, calculated in the previous step, to decorrelation time τ cb,test and τ cr,test using the "tested" model. In the tested model, either β,ρ or the form of g 1 (τ) are invalid for the chosen vessel type. •

Animal experiment
To support the simulations, we have performed a series of recordings in the animal model of ischemic stroke and compared LSCI with Dynamic Light Scattering Imaging (DLSI) measurements. All animal procedures were approved by the Boston University Institutional Animal Care and Use Committee and were conducted following the Guide for the Care and Use of Laboratory Animals. N=3 animals (approx 15-week-old male C57BL/6 mice ) were used. Here we provide only a brief description of the animal preparation and imaging procedures, for details see [14,27]. Animals were anaesthetized with isoflurane, and the first craniotomy was made -2mm window was drilled to expose the middle cerebral artery for the further application of ferric-chloride to induce stroke [27]. Until the stroke induction, the craniotomy was covered with aCSF-soaked-gauze for protection. Next, another frontoparietal craniotomy was opened (4 mm in diameter) to visualize the MCA-supplied cortex. The cortex was then covered with 0.7% agarose solution in aCSF, followed by a 5-mm glass coverslip, sealed with dental cement. The animal was then placed under the imaging system. We waited for at least 30 minutes for the cerebral blood flow to stabilize before starting baseline imaging. During the entire length of the surgical procedure, the animal was heated by a homeothermic blanket with rectal-probe feedback to maintain the temperature at 37 C. Arterial oxygen saturation and heart rate were monitored noninvasively with a paw probe (MouseSTAT Jr., Kent Scientific Instruments). Dynamic Light Scattering Imaging system, image acquisition, and data processing are described in detail in [14]. In the present study, we use the same hardware and procedures, thus only provide a brief summary. Unlike LSCI, DLSI utilizes a high-speed camera (in our case Fastec IL5-S, USA) to capture non-blurred speckle intensity at framerates above 20kHz and directly measure intensity autocorrelation function. Measured g 2 (τ) is then fitted with a comprehensive model, based on the Eq. (10), to estimate the decorrelation time, static scattering, correlation loss and dynamics regime (the form of the field correlation function). In the experiment, DLSI data were collected for 2 seconds during the baseline recording (30 minutes after surgery) and after stroke was confirmed (15 minutes after administering ferric-chloride). The same data were used to obtain LSCI images, as described in [14]. Relative blood flow was then calculated for both techniques: rBF DLSI = τ cb,DLSI /τ cr,DLSI for DLSI and rBF LSCI = (K b,LSCI /K r,LSCI ) 2 for LSCI (corresponds to the MO/SU n=1,simp model). Additionally, LSCI images were analyzed using the MU n=0. 5 and MO/SU n=1 models with ρ, β = 1 to test potential accuracy improvements. Further information on the conventional LSCI hardware, image acquisition and contrast calculation can be found in [24,[36][37][38]. Figure 2 shows simulation results, where the ground truth model's g 1 (τ) form is based on the imaged vessel type (see Table 1), and the commonly applied MO/SU n=1 and MO/SU n=1,simp are tested. There is no correlation loss and no static scattering in these simulations (β = 1 and ρ = 1). From Fig. 2(A), it is clear that using the MO/SU n=1 model to interpret flow in the parenchyma leads to a severe underestimation of the flow drop (up to 100% error at the 10% level of the baseline flow). When the most commonly applied MO/SU n=1,simp model is used to interpret the contrast measurements, the error further increases in the parenchymal regions. Increases in flow, on the other hand, results in smaller errors. The reason is that a decrease in flow makes τ c comparable or shorter than the exposure time T, thus breaking assumption 3 and increasing the contribution of high-order terms in Eqs. (11), (12), (13). These results align well with the experimentally observed discrepancy in ischemic stroke imaging, where relative perfusion measured with LSCI is higher than perfusion measured with other techniques [12,14,27]. LSCI typically measures ≈ 25 − 40% of the baseline flow in the ischemic core, while other techniques, such as optical coherence tomography or positron emission tomography, estimate it to be ≈ 5 − 20% of the baseline flow [14,39]. This corresponds with a 100% to 400% relative error and is partially explained with the simulation results presented in Fig. 2.

Static scattering
As was previously shown with MESI, static scattering also contributes to errors in the estimation of blood flow changes if not properly accounted for [15]. In Fig. 3, we simulate effects of ignoring static scattering (ρ test = 1) on the relative flow measurements under the condition that the actual dynamic scattering (ρ true ) contribution varies from 0.1 to 1. Here, we use the same form of g 1 (τ) for both the ground truth and tested models in these simulations, which are chosen according to the imaged vessel size (see Table 1), β = 1, T = 5ms, and only the ρ value is varied. Results show that ignoring the static scattering, similarly to using the wrong form of g 1 (τ), leads to underestimating the flow change. Interestingly, the error is smallest for parenchyma (up to 100%) and increases for larger vessels (400% and 800%). Such errors, however, are unlikely in typical LSCI applications since the amount of static scattering is negligible in larger vessels unless imaging through the skull or other static tissue. Furthermore, even in the parenchyma, the static scattering contribution will range from 1 − ρ = 0.2 in health [14,15] to 1 − ρ = 0.5 in stroke [14], resulting in a maximum error of ≈ 30 − 40%.

Coherence parameter β
When the τ c <<T or ρ = 1 condition is not satisfied (that is assumption 2 or 3 is false), the parameter β will affect not only the absolute value of BFI, but also the accuracy of the rBF estimates. The reason for this is that high-order terms in Eqs. (11),(12),(13) make increased contributions. To study the effect reductions in the coherence parameter (i.e. β<1), we perform simulations similar to that described above, but now ρ is fixed at 1 and only β true is varied. From Fig. 4 we see that using β test = 1 when β true is below 1 also leads to an underestimation of the flow change. The error is largest when imaging parenchyma (≈70-110% for the typical [35] β ≈ 0.4 − 0.6 values) and is reduced when imaging medium or large vessels (down to ≈15% Fig. 2. Relative error (in percent) resulting from using the MO/SU n=1 model (blue) and its simplified version MO/SU n=1,simp (orange). The form of g 1 (τ) and the τ c value for the true models were defined according to the imaged vessel's size as described in the Methods section. Other parameters were set as follows: ρ=1, β=1, T=5 ms. We see an error upwards of 200%, which is caused by using the wrong form of the field correlation function to interpret speckle contrast in the parenchyma. In A-C, the identity line is shown in dashed gray. were used both as ground truth and tested models, according to the imaged vessel size. Other parameters were set as: β = 1, T = 5ms, ρ test = 1. Only ρ true was varied (from 0.1 to 1). All models show an underestimation of the flow change caused by ignoring the static scattering, with the parenchyma being affected less than larger vessels. and ≈2% respectively). These results highlight the importance of calibrating β for accurate rBF measurements in the parenchyma.

Experimental validation
To support the simulation results, we compare blood flow changes during ischemic stroke estimated with DLSI (rBF DLSI ) and different LSCI (rBF LSCI ) models (see Fig. 5). From panels C-F it is clear that rBF LSCI in parenchymal regions deviates from rBF DLSI more than in mediumsized vessels -up to 260% relative error versus 70% relative error. In line with the simulations, applying the MU n=0.5 LSCI model to the parenchymal data even without adjusting for coherence reductions and static scattering (β, ρ = 1) substantially improves the accuracy and reduces the maximum error to 170% (red lines in Fig. 5(E, G)). Using β = 0.85 and ρ = 0.8 evaluated from the baseline DLSI measurements further decreases the maximum error to 70% (green lines in Fig. 5(E, G)). Changing the LSCI model in medium-sized vessels does not have a strong effect on the accuracy, which also agrees well with the simulations. One should notice, however, that unlike the simulation results, the rBF LSCI to rBF DLSI comparison crosses the identity line at rBF ≈ 0.7 instead of rBF = 1 (Fig. 5(E, F)). This behaviour can be explained by one of the critical parameters being different in the baseline and stroke even when rBF = 1. Specifically, it can be caused by: static scattering, which in stroke experiments becomes spatially heterogeneous from ρ = 0.5 − 0.6 in the ischemic core to ρ = 0.8 − 0.9 in healthy tissue [14,15]; dynamics regime that defines the form of the g 1 (τ) which can change in stroke compared to baseline [14]; and exogenous noise (shot noise, dark noise, fixed pattern) which offsets the contrast value [14,18,24].

Conclusion
In the present study, we analyzed how assumptions in the LSCI theory affect the accuracy of blood flow estimation. We have derived speckle contrast models for the different cases of multiple scattering unordered motion and single scattering ordered motion. Using these models, we simulated the effects of dynamics regime, coherence reduction and static scattering on the relative blood flow estimates in physiologically relevant conditions. Our simulations show that the conventional assumptions in LSCI theory generally result in underestimation of the flow change. While the deviation is typically negligible when imaging vessels larger than 30µm in diameter, it becomes substantial for the brain parenchyma or other capillary perfused tissues (e.g. skin). The underestimation increases further in applications where severe blood flow reduction is observed, such as ischemic stroke, systemic sclerosis or tissue damage (i.g. burns) [12,27,33,40]. In such applications, using the wrong form of the field autocorrelation function, ignoring coherence reductions and/or static scattering, each can lead to upward of 100% relative error in the rBF estimation. When combined, the most typical assumptions can lead to an error of 300% and higher (see Supplement 1). Interestingly, while the parenchyma model is more affected by other assumptions, it is less sensitive to the presence of static scattering compared to larger vessels. This finding has to be accounted for when evaluating the static scattering contribution in parenchyma with MESI [15,18] and when imaging larger vessels through the skull.
In practice, the error in the rBF estimation can be reduced by using the MU n=0.5 model instead of MO/SU n=1,simp . This substantially improves the accuracy (up to 100%) in parenchymal measurements, while the error it can introduce in larger vessels is minor (up to 40%, Fig. 3(A) and Fig. S2). It will also reduce errors that occur from the change between MO/SU n=1 and MU n=0.5 regimes in mid-sized vessels associated with severe blood flow reduction [14]. Further accuracy improvement can be achieved by calibrating β using a static scattering phantom, for which ρ = 0, and thus K ≈ √ β [35,41,42]. Such calibration can be used for as long as the hardware configuration, which defines β, remains unchanged. Following β calibration, one can attempt to account for static scattering by calculating contrast at long exposure times ( T τ c >10000), and estimating ρ according to the K ≈ √ β(1 − ρ) for every imaging session. Following these suggestions will allow making rBF measurements with LSCI up to 200-300% more accurate in the relevant applications.