Effect of vascular structure on laser speckle contrast imaging.

Laser speckle contrast imaging (LSCI) is a powerful tool for non-invasive, real-time imaging of blood flow in tissue. However, the effect of tissue geometry on the form of the electric field autocorrelation function and speckle contrast values is yet to be investigated. In this paper, we present an ultrafast forward model for simulating a speckle contrast image with the ability to rapidly update the image for a desired illumination pattern and flow perturbation. We demonstrate the first simulated speckle contrast image and compare it against experimental results. We simulate three mouse-specific cerebral cortex decorrelation time images and implement three different schemes for analyzing the effects of homogenization of vascular structure on correlation decay times. Our results indicate that dissolving structure and assuming homogeneous geometry creates up to ∼ 10x shift in the correlation function decay times and alters its form compared with the case for which the exact geometry is simulated. These effects are more pronounced for point illumination and detection imaging schemes, highlighting the significance of accurate modeling of the three-dimensional vascular geometry for accurate blood flow estimates.


Introduction
Laser speckle contrast imaging (LSCI) is a fast, non-invasive full-field optical method for visualizing volume integrated blood flow maps in tissue with high temporal and spatial resolution. LSCI's advantage lies in its ability to rapidly image blood flow over a large field-of-view (FOV) with relatively simple and cost-effective instrumentation. LSCI has been widely adopted in a variety of biomedical applications, such as in skin flap monitoring [1,2], neurovascular blood flow imaging in small animal models [3,4], and more recently, clinical applications in real-time intra-operative cerebral blood flow monitoring during surgery [5][6][7].
LSCI uses spatiotemporal fluctuations in random interference patterns of scattered coherent light, known as speckle, to image blood flow non-invasively at a micron scale. The working principle of LSCI is based on the second moment of speckle intensity fluctuations, which results from dynamic scattering in tissue and is related to the electric field autocorrelation function, g 1 (t) [8]. The decay time of electric field autocorrelation function τ c , has been shown to be quantitatively related to the dynamics of particle motion inside tissue [9] such as red blood cells, although the precise relationship between τ c and red blood cell speed is complex [10]. LSCI aims to measure blood flow changes relative to baseline, by quantifying the amount of spatial blurring in the observed speckle contrast K, defined as the ratio of the standard deviation to the mean of the pixel intensity for a small sliding window (i.e., 7 × 7 pixel size) for a given camera exposure time T [10].
In relating the speckle contrast values to the electric field autocorrelation function and deriving τ c, certain assumptions in terms of the order of the motion (ordered or un-ordered) and scattering (single or multiple scattering) are made. Traditionally, τ c values have been derived assuming ordered, single scattering motion. However, in highly vascular tissue, such as in the brain, where photon transport is in the intermediate regime between single and multiple scattering, these assumptions lead to inaccurate decay time estimates. Davis et al. showed the sensitivity of inverse correlation time (ICT=1/ τ c ) over different regions of the cerebrovascular geometry and found that sampling depth and the form of the autocorrelation function differ for different regions of the geometry (i.e., surface vessel versus parenchyma) [11,12].
A few groups have recently aimed to model the effect of geometry on the electric field autocorrelation function in the diffusion region and relate it to absolute blood flow index [13][14][15]. Sakadzic et al., derived a theoretical model for autocorrelation function in the diffusion region taking both diffusive and correlated advective scatterer motions into account [14]. The results highlight the dependence of the theoretical solution on the detailed knowledge of the vascular morphology and flow profiles, and thereby provided a theoretical basis for modeling vascular morphologies in the diffusion region. However, the effect of tissue geometry remains to be investigated in the sub-diffuse region, where source-detector separations are less than a few photon mean free paths in tissue, as is the case for LSCI.
In this paper, we present a fast and robust numerical method for simulating autocorrelation functions and ultimately speckle contrast images for anatomically correct vascular geometries. To the best of our knowledge, this is the first time a simulation of an entire speckle contrast image has been reported. Using the principles of three dimensional dynamic light scattering Monte Carlo (DLS-MC) model [16], we implemented a scalable platform capable of rapidly simulating speckle contrast images and compared results with experimental data. Our platform, described in section 2, simulates the effects of flow perturbations on speckle contrast within seconds, making it feasible to be iteratively used as a forward model. In section 3, we further examine the effect of changes in the geometry on autocorrelation decay times for homogenized geometries and dissolved vascular structures. Two different schemes for homogenizing the geometry are implemented to assess the sensitivity of decorrelation times to both vessel size and geometry depth. The vasculature are dissolved by randomly redistributing the scattering and velocity vectors for voxels classified as vasculature, keeping volume fractions constant in each axial layer. section 2.6 further illustrates the randomization methodology. Several different geometries and illumination schemes are evaluated to determine the generalizability of the results.

Red blood cell flow motion and electric field autocorrelation
Our simulation platform is built on the theory of DLS-MC [16]. In summary, when a plane wave electric field is incident on the surface of a medium, the resulting backscattered electric field at each detector has a phase shift that is the superposition of the momentum transfer contribution from each detected photon that underwent dynamic scattering. Generally, it is assumed that dynamic scattering arises from the interaction of photons with moving red blood cells (RBC) in tissue.
Starting from the definition of the electric field autocorrelation function [17], and statistically accounting for the ballistic motion of RBCs in blood vessels, the ensemble average for a plane wave electric field reduces to the following [16]: In this formulation, t is the decorrelation lag time, P(Y) is the normalized length-dependent absorption weight for the detected photon, k 0 is the wavenumber, and Y is the dimensionless momentum transfer for each detected photon that underwent dynamic scattering (i.e., scattered inside a vessel at least once on its trajectory). The value of Y can be calculated according to the following equation: wherek f ,n andk i,n are the photon's nth scattering and incident unit vectors respectively and V n is the velocity vector of the corresponding scattering location inside the vessel. The sum is over all scattering locations for a single detected photon. Equation (2) captures the ordered motion of RBCs in vessel, however it does not make any assumptions in terms of number of scattering or degree of correlation. Given that dynamic scattering is assumed to be due to moving red blood cells in vasculature, contribution of extravascular region to momentum transfer Y is negligible and thus V n is set to zero in the extravascular region. Using Eqs. (1) & (2), the electric field autocorrelation function can be calculated if the photon scattering positions as well as velocity vectors at those locations are known.

Blood flow velocity profile
Understanding the RBC dynamics inside a vessel is essential to modeling the blood flow speed accurately and arriving at an accurate simulation model for g 1 (t). The general spatial blood flow profile inside a vessel can generally be assumed to have the following form [18]: where v max is the vessel centerline velocity, R is the radius of the vessel, r is the distance from the centerline, a is a scalar to account for the non-zero velocity near the vessel wall, and m is the measure of the bluntness of the profile. It has been shown that for vessels larger than a single red blood cell, the profile is parabolic. In our simulation, we assumed laminar profile for vessels larger than capillaries (i.e. R>11 µm, m=2), uniform profile (m=0) for smaller vessels, and value of a was set to 1. The centerline velocity distributions were set according to arterial, capillary and venous radius-based velocities presented in Lipowsky et al. [19] and further discussed in Ref. [11].

Image acquisition and vectorization
To accurately model the effect of cerebral RBC dynamics on the electric field autocorrelation function, high resolution images of murine cortex were obtained via two-photon microscopy (2PM) imaging. Figure 1(A) shows a sample maximum intensity, z-axis projection of one of the geometries used in our simulation setup. For this geometry, six 2PM stacks were tiled together to achieve a 3D volume of: [x] 840 µm, [y] 1290 µm, [z] 885 µm [20]. The tile edges appear as shadows on the composite image as visible in Fig. 1(A). These high-resolution images were next vectorized using the Segmentation-Less, Automated, Vascular Vectorization (SLAVV) method recently developed in our lab [21]. Vascular objects are vectorized using three-dimensional, multi-scale, linear filtering of unprocessed image volumes. Vectorized vessel objects contain the volumetric centerline flow field, as well as vessel radii information at each voxel. Vessel vector maps are next rendered according to the spatial flow field equation (Eq. (3)) described in section 2.2.
Our rapid vectorization algorithm allows for extraction of strand objects providing a volumetric vascular connectivity map in addition to statistical information, such as volume fraction and vascular morphology, which are further utilized in the subsequent randomization process described in section 2.6. Figure 1(B) shows a maximum intensity projected rendering of vectorized vascular objects corresponding to the 2PM image, with larger vessels displayed in light gray and capillaries in dark gray.  vessels, where field intensity diminishes towards the vessel wall while it remains uniform across capillaries.

Experimental setup
A dual modality setup including both two-photon microscopy and LSCI was used for the experimental validation. Mice were implanted with 4 mm diameter chronic cranial windows over the somatosensory cortex. Immediately prior to 2PM imaging, LSCI images were captured over the geometry shown in Fig. 1. A 785nm laser diode was incident upon the surface at an oblique angle. The LSCI setup shared the optical instrumentation path with the two-photon microscope where a 4x objective was used for LSCI and a 20x, 1.0 NA objective was used for 2PM. All 2PM images were acquired with intravenous administration of Texas Red that was excited at 1050 nm [22,23]. All simulation parameters were adjusted to match the settings in experimental setup.

Simulation pipeline of g 1 (t)
DLS-MC uses a three dimensional voxel-based Monte Carlo simulation to statistically model the photon trajectories based on the optical properties of the tissue at each voxel and then computes dynamic light scattering quantities for subsets of detected photons based on the dynamic scattering events inside blood vessels [16] . The contribution of this paper, in part, is the implementation of a scalable and highly optimized simulation platform for generating both the autocorrelation decay time as well as speckle contrast images for the various geometries presented.
LSCI images were captured immediately before 2PM. The 2PM intensity images were then vectorized as described in section 2.2. The voxelized geometries of the vessel objects were interpolated and rendered at a 1µm 3 resolution, while three dimensional geometries containing tissue types based on the vessel radii at each voxel were generated as input into the Monte Carlo simulator. Table 1. summarizes the optical properties corresponding to different tissue types set by the radii threshold used in our simulation [11,24]. The extravascular absorption and scattering coefficients (µ a , µ s ) were set based on in vitro measurements reported by Yaroslavsky et al. [25]. In each geometry, all vessels with a diameter less than 11 µm were classified as capillaries.
Additionally, in the vectorization step, three dimensional geometries of flow field images corresponding to the vascular centerline unit vectors were generated to be used in the postprocessing step for the generation of g 1 (t) . The flow fields were adjusted to include laminar or uniform flow profile based on the radii of the vessels as discussed earlier. Parallelized DLS-MC simulations were launched on the Stampede2 Skylake compute nodes on Texas Advanced Computing Center (TACC) using the Message Passing Interface (MPI) protocol to simulate a minimum of 20 × 10 9 photon trajectories for each geometry. A circular collimated wide-field beam was set to illuminate 95% of the smallest axis for each geometry. For photons reflected through the top surface of the geometry, all entry and exit locations as well as the photon trajectories through the volume and photon weights were recorded.
Once photon trajectories and path length dependent absorption weights are calculated, the same model can be perturbed rapidly in a secondary simulation step to calculate g 1 (t) for varying particle flow (i.e., red blood cell speeds) according to Eq. (1). The generation of g 1 (t), decay time, and speckle contrast images were implemented in Python and were parallelized using Python's MPI package (mpi4py) on TACC. The program was optimized to simulate any user-defined illumination and detection pattern with the ability to simulate very large detector grids. A numerical aperture of 0.25 and detector size of 4.2 × 4.2 µm 2 were used in our simulation to reflect the optical configuration of our experimental setup.

Simulated speckle contrast calculation
The underlying physical principle behind LSCI is the relationship between the second moment of speckle contrast (K 2 ) and electric field autocorrelation function g 1 (t) as shown below, which assumes ergodicity and sufficiently long integration time [11,26]: Speckle contrast is usually calculated as the ratio of spatial variance σ s to the mean intensity I of a sliding 7 × 7 window swept across the image captured via a camera, for a given camera exposure time T. In the equation above, β is an instrumentation parameter and accounts for the mismatch between the detector and speckle spot size. A form of g 1 (t) is assumed based on the dynamics of the particle flow (i.e., number of scattering and type of motion) by relating K to the electric field decorrelation times and inferring a blood flow measure. Bandyopadhyay et al. formulated this relationship for different kinds of motion [27]. Assuming single scattering dynamics with diffusive motion, or multiple scattering with ballistic motion, g 1 (t) is shown to have the form exp(−t/τ c ). Substituting g 1 (t) in Eq. (4) results in the following second moment relationship: where x = T/τ c , τ c is the electric field decorrelation time, and T is the camera exposure time. This formulation is the general form traditionally used in LSCI accounting for a single exposure time.
Alternatively, if we assume flow dynamics that resemble multiple scattering and diffusive motion, g 1 (t) has been shown to take the form exp(− t/τ c ), which then yields a different relationship between the second moment of speckle contrast and τ c [27]. Davis et al. previously showed that for cerebral hemodynamics, where a single assumption in terms of number of scattering and type of motion is not valid, the simulated form of g 1 (t) lies between the limits of g 1 (t) shown above [11].
To correctly model the dynamics within our simulated volume, we directly calculated K 2 by integrating the area under the g 1 (t) curve according to Eq. (4). In our simulation, we assumed (β=1) and set the exposure time to 3 ms to replicate the settings used for data acquisition in our experimental setup.

Randomization flow
To quantify the effect of geometry on the autocorrelation decay times, we implemented two different geometry randomization schemes where vascular structures were gradually dissolved. These two schemes involved randomization to dissolve vessels either based on a set radii threshold or axial depth. In the case of the radii-based homogenization, all vessels larger than or equal to a set threshold remained intact, while the rest of the geometry was randomized to maintain the same depth dependent volume fractions.
Our vascular vectorization process allows us to quantify a myriad of volumetric statistics including volume fraction and radii histograms. The radii thresholds were set based on the radii histograms obtained for each geometry and were calculated to dissolve 25% of the vasculature structure in increasing steps from intact to complete homogenization keeping volume fractions constant in each axial slice. Volume fractions were calculated based on the total number of voxels classified as vasculature or capillaries divided by the total number of voxels for a given volume. Volumetric radii maps rendered in the vectorization process were utilized to label the vasculature to be dissolved based on the set threshold. Each voxel in the volumetric geometry contains the velocity vector as well as the optical properties associated with the vascular maps in the geometry. If a vascular object was labeled to be randomized, all voxels corresponding to that segment were randomly redistributed in the axial slice. All velocity vector directions were randomized for the shuffled vascular labeled voxels; however, the velocity scalar values as well as the original optical properties were retained in each updated position to maintain the same average flow across each slice.
Given that the geometries were chosen to have different morphological characteristics, these thresholds varied extensively across the geometries.
In the case of depth-dependent randomization, all vessels below a set axial location were randomized, while maintaining the same volume fraction and average flows laterally. We performed two randomization schemes under wide-field illumination and repeated the depth dependent randomization for the case of spot illumination. Figure 2 illustrates the three Fig. 2. Visualization of the three randomization schemes for a sample geometry. Grid overlay depicts the camera array partitioning (not to scale). (A) Radii based randomization with wide-field illumination (yellow shaded region). All vessels with R<22 µm randomized with larger vessels intact, keeping volume fractions the same in each axial layer. (B) Depth-dependent randomization with wide-field illumination. All vessels at the bottom 200 µm (black arrow) randomized, keeping volume fractions constant in each axial layer. (C) Depth-dependent randomization with spot illumination. All vessels at the bottom 200 µm randomized randomization schemes described above. The grid overlay depicts the camera array partitioning (not to scale) and the highlighted yellow regions illustrate the illumination scheme. In the case of the spot illumination, the beam diameter was set to 40 µm directed in the parenchyma region with optical settings similar to our experimental setup.

Simulated laser speckle contrast image
A simulated LSCI image and comparison with experimental data is presented in Fig. 3. As previously illustrated, this geometry was obtained by stitching six 2PM stacks with dimensions spanning 840 µm x 1290 µm x 885 µm in the x, y, and, z directions respectively. Vectorized vascular objects, tissue type, three-dimensional geometry, radii maps and flow fields were rendered at 1µm 3   The speckle contrast (K) values generally range between 0-1, with smaller K values corresponding to higher flow regions. As shown in Fig. 3, speckle contrast values varied between 0-0.16 for both experimental and simulated results. Speckle contrast values over the parenchyma regions were in good agreement, ranging between 0.08 and 0.16 for both cases. However, the simulated speckle contrast over the large superficial vascular region was slightly lower than the experimental data, suggesting that the radial-based speed profile assumed in this region was higher than the actual flow.

Radii-based randomization under wide-field illumination
Three morphologically distinct murine cortex geometries were chosen to evaluate the effect of randomization on the electric field decay times. Table 2. Summarizes the different geometries used in this study. As shown in Fig. 4, the images represent geometries with varying average vascular radii, ranging from small to large, to test the generalizability of the homogenization process. In each plot, labels R1-R5 represent the geometries that were incrementally randomized, with R1 corresponding to the original intact geometry and R5 the fully homogenized volume. The intermediate radii thresholds were set to incrementally dissolve the geometry, based on radii thresholds shown in Table 2 and procedure described in section 2.6.
The plots in Fig. 4 illustrate the effect of randomization for detectors specified over two regions over the surface, namely the parenchyma (red square) and vessels (purple square) as indicated on each image. As shown for a detector over the parenchyma region (Figs. 4(A), 4(C), 4(E)), the electric field decay times shift left gradually and its form is altered as the geometry is homogenized, predicting faster dynamics as the volume is further dissolved. However, for a detector over both descending and ascending vasculature (Fig. 4(D)) and superficial horizontal vessels (Figs. 4(B), 4(F)), the decay times are relatively unchanged for R1-R4. They only shift right in R5 when the whole geometry is homogenized and vessels directly underneath the detector are dissolved. This suggests that the model predicts slower dynamics in the region if a fully homogenized geometry is assumed.

Depth-dependent randomization
The algorithm for depth-dependent randomization homogenizes the geometry below a certain axial threshold, as described in section 2.6. Figure 5 shows the effect of depth-dependent homogenization on the geometry where vessels are gradually homogenized from bottom up, keeping volume fractions constant at each layer. Two different illumination schemes (wide-field and point source) were implemented for each of the randomized geometries to assess the sensitivity of electric field decay times to randomization at depth. Figure 6 illustrates the effect of depth-dependent randomization of the two illumination schemes on the electric field decay curves. In the case of wide-field illumination, Figs. 6(A), 6(B), D1-D5 correspond to the different simulated geometries, where D1 represents the intact geometry, D5 the fully homogenized volume, and D2-D4 geometries with z >200, 100, and 50 µm randomized, respectively. For spot illumination, Figs. 6(C), 6(D), D1 represents the intact   D). In case of A and B, D1 correspond to the intact geometry and D2-D5 to geometry with z>200, 100, 50, and 0 µm randomized respectively (i.e., geometry was completely randomized in D5). For C and D, the point source entrance position is shown with a yellow dot, where D1 corresponds to intact geometry and D2-D5 for volume where z>200, 100, 50, and 0 µm is randomized respectively.
Over the parenchyma regions g 1 (t) plots gradually shifted left as the geometry was further dissolved, thus predicting faster dynamics in the volume for increased randomization. Similarly, the decay curves over the vessel regions were unchanged until the vessel directly underneath the detector was completely randomized. For the geometry where the vessels were completely randomized (D5), the g 1 (t) curve shifted right, predicting faster flow in the volume. In the case of point source illumination, the beam was directed in the parenchyma region (yellow dot on the photon intensity images). The results show a similar trend over the parenchyma region with a more pronounced shift in the decay times as the geometry is dissolved. Additionally, the form of the g 1 (t) curves appear to alter as the slope changes during each randomization step. Contrary to the wide-field illumination schemes, the detector over the vessel region exhibits a higher level of depth sensitivity apparent by the gradual shift of the decay times for the incremental homogenization of the volume. This is the case even when the vessel underneath the detector is intact.
Since both radii-based and depth-dependent randomization schemes under wide-field illumination exhibited similar behaviors, we have only included analysis for depth-dependent randomization under spot illumination. We expect both radii-based and depth-dependent randomization results to follow similar trends under spot illumination.

Effect of randomization on relative to baseline τ c
Most LSCI studies have been limited to measurement of relative blood flow within a single subject during a single experiment (i.e. acute experiments). Extending such measurements to chronic studies where blood flow values can be compared across time, or allowing comparison of blood flow across subjects requires a form of absolute measurement of flow.
In previous sections we focused on the impact of geometry on absolute τ c errors. In order to analyze the impact of homogeneity assumptions on the accuracy of normalized to baseline measures, we analyzed the sensitivity of τ c to perturbations in individual strands and compared the results for both intact as well as fully homogenized geometries.
In the geometry shown in Fig. 7 (Geometry-3), the large superficial strand (red dotted region) was perturbed to reflect the relative velocity profile shown, and normalized 1/τ c values relative to baseline were calculated for detectors over both vessel (Fig. 7(A)) and parenchyma regions Fig. 7. Sensitivity analysis of Intact (R1) and fully randomized geometry (R5) to strand segment flow perturbations relative to baseline. The large superficial vessel (red dotted strand) was perturbed as shown with relative velocity ranging from 0.5x to 2x the baseline values. (A) Normalized 1/ τ c for detector over the perturbed vessel segment (blue square). The dotted black line shows a sensitivity of 1, drawn for reference. (B) Normalized 1/ τ c for detector over the parenchyma region (red square) ( Fig. 7(B)). As shown for a detector over the perturbed vessel segment, normalized 1/τ c values for both intact and homogenized geometries result in the same sensitivity slope. However for a detector over the parenchyma region, the homogeneity assumption leads to approximately 2x increase in the sensitivity slope when compared with the intact geometry.

Comparison of simulated and experimental laser speckle contrast images
The question of correct statistical model for RBC motion is a topic of numerous studies [13][14][15]. In applications such as Diffuse Correlation Spectroscopy (DCS), where photon directions are fully randomized before exiting the volume, the ensemble average particle motion in the electric field autocorrelation function is characterized by mean squared displacement. In performing LSCI measurements in cerebral tissue, photons sample multiple vasculature with different orientations and flow profiles along their trajectory. However, these velocities are not random and often differ considerably from their extravascular surroundings. This leads to the breakdown of the Gaussian lineshape assumption characterized by multiple scattering and diffusive motion [16].
We computed speckle contrast values by directly integrating the simulated g 1 (t) curves without making any assumptions regarding the form of the autocorrelation function. Figure 3 shows that the simulated LSCI image is in agreement with the experimental data, as discussed above. For both simulated and experimental results, the values of K range between 0 and 0.12 with smaller K corresponding to faster flows. These results are particularly interesting, as DLS-MC theoretically only models the advective motion in volume and assumes that the diffusive motion due to scatterer concentration gradient is negligible during a single imaging session. It has been shown that LSCI is sensitive to both advective and diffusive flux [28]. The strong agreement between the simulated and experimental K values suggests that ordered motion statistics used in deriving the autocorrelation ensemble average adequately models the particle dynamics for LSCI in the sub-diffuse regime.
A source of discrepancy between simulated and experimental speckle contrast data appears around the first branch off from the large vessel in the top left quadrant of Fig. 3(A), where the vessel region seems to be missing in the simulated data. This can be explained by observing the maximum intensity projection of the vectorized geometry shown in Fig. 1(B), which was used as the input geometry to DLS-MC. The light-saturated area around the large diagonal superficial vessel in the 2PM image renders the automated vectorization incapable of detecting the branching point. Thus, the vessel appears as a short floating segment with a smaller vascular velocity and consequently higher K value in the simulated image.
Additionally, the simulated speckle contrast values over superficial vasculature were smaller than the experimental data. Our simulation results indicate that these discrepancies can be mitigated by adjusting the radii-based speeds attributed to the vasculature. This inference is supported by our finding that speckle contrast values over vessel regions are dominated by the photons sampling the vasculature directly underneath the detector. Since the velocities were assigned according to radii-based average values reported in literature for a given hematocrit level, and not through direct measurement, this discrepancy is expected. Given the non-linear relation between flow perturbations inside the vessel and the ensemble average momentum transfer of the detected photon, one would need to reconstruct the velocities by iteratively adjusting the flow speeds in an inverse problem in order to match the experimental and simulated results.

Effect of randomization on g 1 (t) over the parenchyma region
We simulated the effect of randomization under two different illumination schemes for three morphologically distinct murine cortices. We evaluated the effect of geometry by randomizing either based on a set radial threshold or by homogenizing the geometry below a set depth. These two randomization schemes allowed us to compare the sensitivity of the speckle contrast values to vessel sizes across the whole geometry versus vasculature at a given depth.
Our results show that for detectors over the parenchyma region, the autocorrelation function decays faster as the geometry is randomized in all cases. As shown, τ c values reduced up to 10x in some cases as the geometry was dissolved. The extent of variations depended on the range of the vessel dimensions in each geometry and the percentage of the volume randomized in each step. This is evident by comparing the results shown in Figs. 4(B) and 4(C). Geometry-2 ( Fig. 4(B)) consisted mostly of capillaries and smaller vertical vasculature with vessel diameters ranging between 11-35 µm. Fully randomizing geometry-2 resulted in approximately 2x reduction in τ c values. This reduction was approximately 10x for Geometry-3 ( Fig. 4(C)) where the geometry consisted of a mixture of capillaries and large superficial vasculature ranging from 11-100 µm in diameter.
It was shown previously that detectors over parenchyma are sensitive to the top 500 µm of the geometry and 200 µm laterally [11]. Our results indicate that in the case of wide-field illumination, this sensitivity is the highest in the top 150 µm. The sensitivity is further increased by implementing a point source/detector scheme where the photon spatial sampling profile is modulated as a function of source and detector separation. This is apparent by the more uniform shift in the τ c values across all levels of randomization as shown in Fig. 6(C) compared with the results in Fig. 6(A) where the largest shift occurs when the top 150 µm of the geometry was dissolved.
This is in part due to the fact that dissolving vasculature and distributing vessel voxels randomly in tissue will increase the probability of a photon getting dynamically scattered before exiting the volume. This leads to an increased momentum transfer phase shift and ultimately faster decorrelation of g 1 (t). As described earlier, the velocity unit vector directions were randomized as vessel voxels were distributed axially, while maintaining the speeds from the original geometry at each distributed voxel. This would not have a significant impact on the ensemble average if scattering in the original geometry was truly isotropic and all scattering events were independent. However, photons entering any vessel larger than a capillary are likely to undergo multiple consecutive correlated scattering events [13] . Our results indicate that correlated multiple scattering events, on average, account for 50% of all dynamic scattering events for detectors over a vessel and 31% for detectors over a parenchyma region. Thus the calculated q · ì v values differ greatly when photons have sampled multiple vessels with different orientations versus when the voxels have been distributed randomly in the volume. The loss of vascular structure eliminates the probability of multiple scattering in a vessel and thus significantly alters the probability distribution of scattering events as evident by our randomization results.

Effect of randomization on g 1 (t) over the vessel region
It was shown previously that under wide-field illumination photon sampling is heavily weighted towards the top 150 µm axially and 50 µm laterally for surface vessel ROIs [11]. In cases where larger superficial vasculature were embedded lower axially in the geometry, this sensitivity was extended in depth to the location of the vessel. Our results show that under wide-field illumination, the decorrelation times remained unchanged for each level of randomized geometry until the superficial vessels directly underneath the detector were dissolved. As shown in Figs. 4(B), (D) and (F) dissolving the vessel regions resulted in up to 10x increase in the τ c values. This increase was less (5x) in case of Geometry-3 ( Fig. 4(F)) where dissolving the large superficial vessel resulted in many distributed voxels in the top 100 µm of the geometry.
To analyze the impact of randomization on the vessel regions, it is important to note that the photons reaching the detector over a surface vessel ROI all get scattered inside that vessel before exiting the volume. Additionally, as discussed earlier, once a photon enters a large vessel it is likely to undergo multiple correlated scattering events inside that vessel. These two factors together result in a large contribution from the said vessel to the detector's momentum transfer ensemble average. However, once the vasculature is dissolved, the contribution to the ensemble average will be distributed across the geometry instead of being localized in the vessel. This results in a fewer number of photons getting dynamically scattered before reaching the detector. Moreover, the large number of correlated consecutive scattering events will be replaced by fewer isotropic scattering events where ensemble q · ì v values approach a mean of zero. These outcomes together will result in a slower decay of the g 1 curve as shown in the plots.
It is interesting to note that the sensitivity in the vessel ROIs under point source illumination seem to be more distributed across the geometry as shown in Fig. 6(D). This is in part due to the fact that the entrance and exit locations of the photons are restricted in this scheme, thus forcing the photons to sample a certain trajectory along their path. This contrast with the case of wide-field illumination, where the majority of detected photons' entrance positions are localized to areas surrounding the vessel, thus limiting the photon's penetration depth.

Randomization results in presence of optical property variations
We previously reported on the effect of intra and extra vascular µ a and µ s variations on the number of dynamic scattering events [29]. It was shown that photon sampling was relatively insensitive to both vascular and extravascular µ a values. Extravascular scattering coefficients were shown to have negligible effect on the number of dynamic scattering events, thus extravascular µ s variations are not expected to have a significant impact on the overall τ c values. Similarly results showed that doubling the intravascular µ s from 60mm −1 to 120mm −1 did not have a large impact on the relative depth-dependent scattering distributions for both parenchyma and vessel ROIs. While it is reasonable to assume that the uncertainty in µ s of blood would perhaps have a slight effect on the overall τ c values, we expect the shift in τ c values with loss of geometry to follow the same trend even as intravascular µ s values are varied.

Analysis of the relative 1/τ c sensitivity
Comparison of relative 1/τ c sensitivity for detectors over both vessel and Parenchyma regions, (Fig. 7), suggests that relative blood flow measures are more accurate for superficial vascular regions and particularly when the detector is over the perturbed vessel. This is evident given that for both intact and fully homogenized geometries, the change in the normalized to baseline 1/τ c as a function of velocity perturbation follows the same slope. However, over the parenchyma region, randomization assumption has a large impact on the slope and sensitivity of normalized 1/τ c when the surface vessel strand velocity is perturbed. These results are particularly of interest for problems involving inverse volumetric blood flow reconstructions where homogeneity assumptions are made in solving the inverse problem. Our current results indicate that a priori knowledge of larger vascular structures is necessary for accurate volumetric reconstruction of BFI.

Conclusion
We presented a fast model for simulating speckle contrast images using anatomically correct cerebral geometries. Our platform allows for the rapid update of speckle contrast images per vascular flow perturbations, thus making it feasible to be used as a fast numerical forward model. The model is not limited to simulations of LSCI and is equally applicable to DLS measurements since it simulates the full correlation functions for a given geometry and flow distribution within that geometry.
In order to quantify the effect of tissue geometry on g 1 (t) curves, we implemented three different randomization schemes while keeping volume fractions constant throughout the geometry. Our homogenization schemes illustrated the significance of microvascular and vascular structure in determining the shape and decay time of the correlation functions. Traditionally, homogeneity assumptions are made in deriving analytical solutions for calculating g 1 (t) in the diffusion region [30]. Our results suggest that these assumptions may lead to gross (up to 10x) over/under estimation of the actual underlying particle motion, depending on the placement of the detector. The effects were more pronounced for point illumination and detection imaging schemes compared with wide-field illumination paradigms.
In addition to analyzing the effect of randomization on absolute 1/τ c measures, we also examined the sensitivity of relative normalized 1/τ c to velocity perturbations comparing the results for both intact and fully randomized geometry. Our results highlight the need for detailed modeling of cerebral tissue in deriving accurate blood flow estimates in both diffuse and sub-diffuse regimes especially for detectors over the parenchyma region.
These results are particularly of importance in inverse reconstruction problems where homogeneity assumptions are often made in solving the analytical forward mode. We should note that the homogenization schemes often cited in literature, replace the velocity values with a single spatial average that correspond to the original velocity magnitudes weighted by the vascular volumes of the vessels of interest, resulting in a more averaged and smoothed volume [13]. We note that our randomization scheme is a step between the intact and fully homogenized geometry as described above. We expect that assuming spatial averages instead of redistributing the vascular voxels would lead to even larger shifts in the τ c values as the vascular structure is dissolved.
As a future direction, we will investigate utilizing our fast simulation platform to accurately reconstruct three dimensional volumes of cerebral blood flow in the sub-diffuse regime.