Establishing the quantitative relationship between diffuse speckle contrast analysis signals with absolute blood flow.

Diffuse speckle contrast analysis (DSCA) measures blood flow in deep tissues by taking advantage of the sensitivity of the speckle contrast signal to red blood cells (RBCs) motions. However, there has yet to be presented a clearly defined relationship between the absolute blood flow BFabs and the measured speckle contrast signal. Here, we derive an expression of linear approximation function for speckle contrast, taking into account both shear-induced diffusive and correlated advective RBCs motions in the vessels. We provide a linear relationship between the slope k slope of this linear function and BFabs. The feasibility of this relationship is validated by Monte Carlo simulations of heterogeneous tissue with varying vessel radii. Furthermore, based on this quantitative relationship, we can determine the relative contributions of diffusive RBCs motion on the reduction of speckle contrast, considering different vascular morphology and flow profiles.

from the dynamics, and, therefore, extracts a quantitative estimate of blood flow in vivo [22]. Besides, the fast cameras with a high frame rate [23,24] can rapidly improve acquisition speed of speckle images.
In fact, DSCA and DCS are examining different aspects of the same entity [25], i.e., the temporal electric autocorrelation function of the diffusely reflected speckles caused by the movement of red blood cells (RBCs). Currently, the most common model for describing the RBCs motion is a Brownian diffusion-like process, which provides a better characterization of the autocorrelation function signals in a range of subjects and tissue types than the expected random flow model which describes the RBCs motion is ballistic. Based on Brownian motion model, we have developed and tested a thorough speckle contrast model [26][27][28] in DSCA for quantitative measurement of the RBCs motion. Although the blood flow index (BFI) in Brownian motion model has been shown to be approximately proportional to absolute blood flow BF abs , the interpretation of the BFI is complicated due to its units of cm 2 /s and the value of this proportionality between BFI and BF abs is not clear in DSCA. Therefore, establishing the quantitative relationship between BFI in DSCA and BF abs is urgently necessary and significant.
For this quantification, firstly it is necessary to clearly understand intravascular RBCs dynamics in the tissues. It has been found [29,30] that RBCs undergo shear-induced displacement in blood flow, and the calculated electric autocorrelation function largely depends on this shear-induced diffusive RBCs motions, while the effect of correlated advective RBCs motion is relatively smaller [31][32][33]. The challenge of understanding the effect of diffusive and advective RBCs motions on DSCA signals stems from the fact that the speckle contrast signal is obtained from dynamic scattering inside the vessels and the vessels comprise only a small fraction of the tissue volume. Therefore, it is necessary to separate the intravascular scattering events from the extravascular scattering, and explicitly record the location of each dynamic scattering event in the vessels. Monte Carlo simulation [5,31] can accurately model dynamic scattering event in the realistic and complex tissue.
In this paper, the reduction in speckle contrast due to the relative contribution of diffusive and correlated advective RBCs motions is studied by theory and Monte Carlo simulations. We derive the linear approximation expression for speckle contrast in the heterogeneous tissue with varying vessel radii and blood flow profiles, and establish the quantitative relationship between the slope k slope of this linear function with the absolute blood flow BF abs . Meanwhile, we employ three-dimensional Monte Carlo simulations of photons propagation through tissues to support these results. This linear approximation expression can be also used to model the speckle contrast signals in the realistic tissue with the heterogeneous blood flow profiles and a complex vascular network.

Intravascular red blood cell motion
The question of the proper statistical model for describing intravascular RBCs motion and ultimately relating this model to the measured electric autocorrelation function or speckle contrast have been the subject of numerous studies [30,[34][35][36]. Usually, RBCs displacement has been modeled as random flow with <Δr 2 (τ)> = V 2 τ 2 , where <Δr 2 (τ)> is the mean-squared displacement (MSD) of RBCs in time τ and V 2 is the second moment of the velocity of RBCs, or as Brownian motion with <Δr 2 (τ)> = 6D B τ, where D B is the effective Brownian diffusion coefficient. Unexpectedly, Brownian motion model leads to better description of the experimental results in validation studies compared with random flow model.
It should be noted that both models mentioned above have generally assumed that the dynamic scattering events with RBCs in vessels are uncorrelated. This assumption is valid if the photons scatter only once inside a vessel and then encounter another vessel before their direction becomes randomized by the scattering from extravascular tissue. In fact, the scattering length of the photons inside the vessel at the common laser wavelengths (790-810nm) used in DCS or DSCA is ~12μm [37], which is smaller than some vessel diameters in the tissue. This leads to the breakdown of this uncorrelated scattering assumption because the photons in a large vessel will most likely undergo multiple scattering events which are in fact highly correlated. When multiple scattering in a single vessel is present, the contribution of each scattering event to dynamic information depends on the velocity difference between subsequent correlated scattering events [38].
To account for the effect of correlated RBCs motions on the speckle contrast signal, it is necessary to understand RBCs dynamic in the vessels. Usually the spatial flow speed profile of the vessels has the form as [39] where v max is the maximum speed at the center line of the vessel, R is the radius of the vessel, r is the distance from the vessel center line, a is a scale factor that determine the nonzero velocity at the vessel wall, and m is the bluntness of flow profile. Equation (1) describes a general form of advective flow in vessels and is experimentally verified by the method of dynamic light scattering -optical coherence tomography (DLS-OCT) [33]. Usually, the flow profile of bigger vessels (R>20μm) [33] can be considered to be parabolic (i.e., m = 2) and there is more bluntness in smaller vessels. Meanwhile, recent studies [31,32] found that in addition to advective motions RBCs also undergo shear-induced diffusive motions in the vessels. The shear-induced diffusion coefficient D(r) is proportional to the shear rate α ss [40], i.e., where α ss has a magnitude of ~0.1 to 0.5 × 10 −6 mm 2 [33]. Based on Eq. (1) and (2), we can determine the velocity of the location of each scattering event in the vessel taking into account both diffusive and advective RBCs motions. In the following section, we will show the effect of correlated RBC motion on electric autocorrelation function and speckle contrast.

Diffuse speckle contrast analysis (DSCA)
A speckle pattern is a random interference pattern generated when photons from coherent light are scattered in a random medium. If the scatterers in the medium are mobile, as in the case of RBCs, the time integrated speckle pattern recorded by the camera becomes blurred. DSCA qualifies this degree of blurring by the ratio of the measured standard variance ( I σ ) to the mean intensity ( I μ ) in spatial or temporal domain, as , , where ρ is the source-detector (SD) separation, T is the exposure time of the camera, and K is the speckle contrast. In fact, the speckle contrast is essentially related to the normalized electric autocorrelation function g 1 (r,τ) [15,16] as where β is a constant determined by the experimental setup [41], g 1 (ρ,τ) = G 1 (ρ,τ)/G 1 (ρ,0), G 1 (ρ,τ) = <E(ρ,t)E*(ρ,t + τ)>, τ is the delay time and E(ρ,t) is the light electric field.
Formally, the transport of G 1 (ρ,τ) in multiple scattering media is well modeled by the correlation diffusion equation (CDE) [7,42]. Sakadžić et al. [32] recently considered the effect of the correlated RBCs motion with diffusive and advective dynamics on the CDE, and derived the expression of G 1 (ρ,τ) in a realistic vascular network with a heterogeneous distribution of the vessels with different radii and flow velocities. For simplicity, we consider only one vessel type with the same radius R and blood flow in the tissue and the analytic solution of G 1 (ρ,τ) for the semi-infinite media is given by [32] ( ) Here, 2 , R eff is the effective reflection coefficient accounting for the index mismatch between the tissue and surrounding medium [43], a μ is the average absorption coefficient, s μ′ is the average reduced scattering coefficient, k 0 is the wavenumber of light in the media, l tr is the mean free path in the vessel, δ ves accounts for the probability of the dynamic scattering from RBCs. Note that the average optical properties ( a μ and s μ′ ) consider the heterogeneity of the tissue and can be measured by multi-distance DCS [44,45] or frequency-domain spatially resolved spectroscopy [46]. Besides, the parameter δ ves is currently difficult to determine experimentally and usually is estimated by the volume fraction of the vessels in the tissue. D av and v av are the average value of diffusion coefficient D(r) and advective motion v(r), respectively. This allows us to write [32] ( ) ( ) Based on Eq. (4)-(9), we can determine the theoretical dependence of speckle contrast on the SD separation or exposure time under the condition of the combination of diffusive and convective RBCs motion. Meanwhile, we note that it is difficult to obtain the analytical expression for the speckle contrast function due to the complex computation process of integral equation, and the necessary steps to be taken to overcome it.

Linear approximation for diffuse speckle contrast analysis
Our previous works [27] have demonstrated that the relation between 1/K 2 and exposure time can be described by a linear approximation equation and the slope k slope is equal to the inverse of the correlation time (τ c ) of autocorrelation function g 1 (ρ,τ). In this section, we derive the expressions for the correlation time τ c and this linear approximation model considering diffusive and advective motion.
To obtain correlation time in Eq. (5), we firstly note that the function G 1 (ρ,τ) is derived from the radiative transfer equation using diffusion approximation. Since diffusion approximation is usually valid when the SD separation ρ is much larger than the transport mean-free path z 0 , r 2 in Eq. (5) can be approximately given by Meanwhile, we can introduce another approximation ( ) ( ) which is usually satisfied at the larger SD separation. This procedure yields the expression for G 1 (ρ,τ c ) So g 1 (ρ,τ c ) can be expressed as Equation (13) can be subsequently written as Thus, Eq. (14) can be further represented as Equation (17) is a general result that determines the correlation time in DCS considering the combination of diffusive and advective RBCs motion. Meanwhile, we have demonstrated that the slope k slope of the linear approximation model for DSCA is equal to the inverse of the correlation time [27], i.e., k slope = 1/τ c . Thus, the corresponding speckle contrast 1/K 2 can be approximatively expressed as where b in is the intercept and can be obtained by the linear fitting.
In order to clearly understand the relationship between the blood flow index k slope in Eq. (19) and absolute blood flow BF abs , it is important to have a clear definition of BF abs used in DCS and DSCA. Usually BF abs is defined as the volume of blood transiting through a given cross-sectional area per second, which is given by cross-sectional area times average RBCs speed, i.e., 2 abs av Therefore, a clear linear relationship between k slope and BF abs can be established by From Eq. (20), it is apparent that the relative change in BF abs can be rapidly obtained by calculating the corresponding change in k slope . Meanwhile, knowledge of this proportionality in Eq. (20) is required to determine BF abs form the measured k slope . However, this proportionality is complicated and related to multiple parameters such as the vessel radius, RBCs rheology (a, m, and α ss ), optical properties of the tissue, δ ves and SD separation, etc. Some parameters in tissue are now not accessible by experimental method and further studies are necessary to discuss this problem. Besides, an important assumption for deriving Eq. (20) is that the first order expansion ( ) ( ) is a good approximation.
In fact, this assumption is valid for the large SD separation, i.e., 1 s ρ μ′  . Analytical and numerical computations of the influence of this error due to this approximation on k slope are given in the following section.

Monte Carlo simulation
For studying photon traversal inside the heterogeneous tissue and considering the effect of correlated RBCs motion on the speckle contrast, Monte Carlo simulation can be used to track the consecutive scattering events in the tissue. Our Monte Carlo code is based on a three dimensional voxelized model developed in Ref [47]. Each voxel of the Monte Carlo geometry is assigned the optical properties that determine the distribution of photon scattering lengths and the degree of the absorption. The Henyey-Greenstein phase function with corresponding anisotropy g [48,49] within the voxel containing the scattering event is used to calculate the scattering angle.
The  1) and (2). Based on these recorded values, the correlation function G 1 (ρ,τ) can be given by summing over all detected photons [31] ( ) , ,, , , , where N P is the number of photons detected at the specific SD separation, N s,n is the number of scattering event in N v,n vessels for the n'th detected photon, N s,n,l is the number of consecutive scattering events of the n'th detected photon in the l'th vessel, N tis is the number of tissue type, , is performed in the same vessel to consider the degree of correlation between N s,n,l consecutive scattering events. This equation considers the influence of RBCs displacements due to the correlated advective and diffusive motion on the decay of G 1 (ρ,τ). Using Eq. (4) and the information about RBCs speed and shear-induced diffusion, we can further calculate the expected speckle contrast values by MC simulation. To demonstrate the flexibility and validity of linear approximation model for DSCA in the heterogeneous media, a voxelized tissue-mimicking geometry in Fig. 1 is used in conjunction with Monte Carlo simulation to generate the speckle contrast measurements. Figure 1 shows the (x, z) cross section of the geometry with a size of 6 × 3cm. For simplicity, the blood vessels have the same radius R and are all oriented along the y-axis with an equal spacing h = 200μm [31] in x-axis and z-axis. The values of the radius R and the vessel spacing h determine the volume fraction of the vessels. If needed, other values for h can be further used to consider different volume fraction of the vessels. Therefore, the geometry has an infinite extent in y-axis due to the translation symmetry. In our simulations, the blood flow profiles in Eq. (1) are assumed to be constant in time and approximately parabolic flow for the vessels with a larger radius (R>20μm), i.e., a = 1 and m = 2. In general, the vessels with a smaller radius tend to have a uniform flow profile and in this case RBCs motion is dominated by diffusive motion. Increasing the vessel radius x z and considering parabolic flow profiles can find the influence of correlated RBCs motion on speckle contrast signals. If needed, other values for a and m can be further used to consider different blood flow profiles in the future. Besides, the value of shear rate α ss is set to 0.24 × 10 −6 mm 2 from Ref [33]. measured in vivo. For each simulation, 10 8 photons at a point are used to illuminate the surface of the geometry along z-axis and a 0.5mm × 0.5mm square region of interest (ROI) at different distances ρ from the source point is used as the detector. All of the information for calculating Eq. (21) and (4) are recorded by our MC simulation. Besides, the intravascular and extravascular optical properties at 800nm for a blood hematocrit 40% [37] are shown in Table  1.

Corrected k slope expression
Firstly, we test the accuracy of the expression for correlation time τ c using theoretical data as shown in Fig. 2. An important assumption in this expression is that there are smaller errors arising from truncating SD separation terms in the Taylor Series expansion This assumption is valid when using large SD separations ρ, i.e., 1 s ρ μ′  . Therefore, it is necessary to quantify the range of SD separation for which this expression can be accurately employed. The relative error in τ c due to this assumption is defined as ( ) 1 Error in ln , 1. In order to correct this error due to this assumption, the modified expression for correlation time τ cr is introduced by Note that knowledge of tissue optical properties ( a μ and s μ′ ) and SD separation can determine

Validation with MC simulation data
In this section, we validate the linear approximation model by MC simulation results. We vary RBCs speed (v max from 1mm/s to 6mm/s with a step of 1mm/s) and vessel radius (R = 20, 25, 30, 35, 40μm) with a fixed vessel spacing (h = 200μm) in MC simulations to test the dependence of k slope on BF abs . For all simulations, the intravascular and extravascular optical properties at 800nm for a blood hematocrit 40% are shown in Table 1 Here, we use R = 30μm vessels and SD separation is 2cm. Figure 3(a) and 3(b) show K 2 and 1/K 2 results from MC simulations and theoretical calculations, respectively. In this example, the vessel radius is set to 30μm and SD separation is 2cm. Note that the theoretical results are numerically obtained by Eq. (4)- (9). These results show that the theoretical speckle contrast model fits MC simulation results well and the linear relation between 1/K 2 and exposure time is obvious. Figure 3(c) shows the fitted values of k slope versus BF abs , compared with those calculated by the corrected relation using Eq. (24). These results in Fig. 3 demonstrate the validity of linear approximation model for DSCA and confirm that k slope increases linearly with BF abs . Figure 4 shows the calculated absolute blood flow BF abs from the linear fitting k slope results at different SD separations. Note that the changes in the absolute blood flow BF abs are accomplished by varying the vessel radius (from 20 to 40μm) with a fixed flow speed v max = 3mm/s and vessel spacing h = 200μm. The solid lines in Fig. 4 represent the expected BF abs results for SD separations of 1.0, 1.5 and 2.0cm. It is clearly found that this slope k slope of the linear approximation model predict well the absolute blood flow BF abs . The change in BF abs can be simply obtained by calculating the change in k slope .

Discussion
Diffuse speckle contrast analysis (DSCA) has been employed extensively in blood flow measurement, in large part because of its simplicity. Usually the calculated speckle contrast fluctuations are driven by the RBCs motion and thus by blood flow. The question as to which model best characterizes RBCs motion within vessels has been widely studied. In this paper, we derive linear approximation expression for 1/K 2 taking account both shear-induced diffusive and advective RBCs motions. To validate this theoretical model, we have performed Monte Carlo simulations in a simple tissue-mimicking geometry with parabolic flow profile and varying vessel radius. These theoretical and simulation results allow us to write a linear relation between the measured k slope and absolute blood flow BF abs . Meanwhile, in order to connect k slope calculation to BF abs some effects on this model need to be further discussed in actual applications.
It should be noted that our previous works [27] have deduced an analytical expression k slope-DB for diffusive RBCs motions, i.e., 1 It is worthwhile to determine the relative contribution of diffusive motion on k slope in Eq. (24). We simply quantify this contribution by calculating the ratio between k slope-DB and k slope as shown in Fig. 5. Here we have used a fixed vessel volume fraction δ ves = 0.02, parabolic flow profile, a shear rate α ss = 0.24 × 10 −6 mm 2 and the same optical properties in Table 1. Figure 5 shows that increasing the vessel radius and decreasing SD separation both result in reduced importance of diffusive RBCs motion. Especially for the large vessel radius at the short SD separation [50], the advective RBCs motion may dominate k slope . For the small vessel radius, k slope mainly depends on diffusive RBCs motion. Furthermore, this contribution of diffusive motion also depends on other parameters of the vessels, such as shear rate α ss and blood flow profile (a and m). So far the measured values of α ss have large difference from in vivo and ex vivo [33,40]. One can note that the higher the shear rate α ss , the larger contribution of RBCs motion. The influence of blood flow profile on this contribution can be 30  calculated by Eq. (24) and k slope-DB . Equation (24) suggests an approach to determine the relative contributions of diffusive and advective RBCs motions on the reduction in speckle contrast. The source-detector (SD) separations used in MC simulation results are such that SD direction is perpendicular to the blood flow direction as shown in Fig. 1. In fact, the actual vessel directions in the tissue are random and complex. Figure 6(a) shows the effect of specific SD direction, which is perpendicular or parallel to vessel direction, and all used SD directions on the obtained speckle contrast 1/K 2 results from the same MC simulations. It is obvious that the results 1/K 2 tend to be the same and the corresponding absolute percentage changes in 1/K 2 compared with the perpendicular case are shown in Fig. 6(b). The relative error in 1/K 2 is no more than 2%. Meanwhile, Ref [31] have verified that the preferred direction of the vessels only introduce a small bias for the correlation function g 1 (ρ,τ) compared with a truly random direction of the vessels. It is apparent that the calculated 1/K 2 at a certain SD separation doesn't depend on the relative direction of SD separation with respect to vessel direction. Fig. 6. Effect of SD direction with respect to vessel direction on the obtained speckle contrast 1/K 2 results from MC simulations. The red and black lines in (a) indicate 1/K 2 from SD separations that is perpendicular and parallel to vessel direction, respectively. The blue lines in (a) indicate 1/K 2 from SD separations in all directions. The corresponding percentage change in 1/K 2 at different SD separations compared with perpendicular case is shown in (b). Here, we use R = 30μm vessels with a fixed flow speed v max = 3mm/s and a fixed vessel spacing h = 200μm.
For simplicity, we neglect the variation in radius and in flow speed among different vessels, and consider only one type of the vessel in this paper. In fact, when performing DSCA or DCS measurements in the tissue, the vessels encountered are different and  complicated. For instance, capillaries with a small radius (<5μm) tend to be have similar velocities, but veins, venules, arterioles, and arteries all have different radii and flow speeds. Usually the photons come across multiple vessels along their path, each with different speeds and radiuses, and often interact with extravascular compartments. These contributions taken together cause the decay of speckle contrast function. Note that the analytic solutions for G 1 (ρ,τ) [32] in Eq. (5)-(7) also can be used to predict DSCA signals in a realistic vascular network with a heterogeneous distribution of vessels with different radii and blood flow velocities. In this case, F(τ) in Eq. (7) can be expressed as [32] ( ) ( where δ ves [R,v ves,av (R)] and D ves,av [R, v ves,av (R)] are the volume fractions and average diffusion coefficient of the vessel with radius R and average velocity v ves,av (R), respectively. Therefore, the average diffusion coefficient D av and advective speed v av in the linear approximation model have a complex nonlinear dependence on the distribution of vessel radius and RBCs speed distribution. It is important to point out that the differences of vessel radius and flow speed don't affect the validity of linear approximation model. This linear approximation expression can be used to model the speckle contrast signals in the realistic tissue with the heterogeneous blood flow profiles and a complex vessel geometry. Further measurements and numerical modeling of a realistic vascular network are needed to study the relation between the DSCA measurements and blood flow in the future.
Our results show that we can linearly relate the measured k slope to absolute blood flow BF abs . Note that the proportionality is related to the vascular morphology, RBCs rheology, optical properties of the tissue, and SD separation. In order to accomplish the quantitative measurement of BF abs from the measured k slope , we need to have a detailed knowledge of these parameters. However, some parameters now are not readily accessible. But the progress in the experimental methods may make these parameters available in the future. Further research is needed to study the quantitative measurement of BF abs .

Conclusion
In summary, we have presented the theoretical derivations for linear approximation model in DSCA that consider both diffusive and correlated advective RBCs motion, and demonstrated that the slope k slope of this linear model has a linear relation with absolute blood flow BF abs , which are in agreement with our Monte Carlo simulation results. Through this expression for k slope , we can determine the relative contribution of diffusive motion on the reduction of speckle contrast. We expect further studies to use this relation to obtain a direct measurement of BF abs in the more complex and realistic configurations of the tissue.