Sensitivity of laser speckle contrast imaging to flow perturbations in the cortex

Mitchell A. Davis; Louis Gagnon; David A. Boas; Andrew K. Dunn

doi:10.1364/BOE.7.000759

1. Introduction

Laser speckle contrast imaging (LSCI) is a powerful, low cost method of imaging dynamic motion with high spatial and temporal resolution. LSCI has been heavily adopted for neuroscience applications such as blood flow imaging of neurovascular pathologies [1, 2], functional activation [3], and even human cortical blood flow imaging during neurosurgery [4–6].

The dynamic imaging capability of LSCI arises from interactions between coherent photons and tissue. When these photons interact with moving particles, the variance of the detected intensity fluctuations of the speckle pattern changes, which cause spatial blurring when averaging over a fixed exposure duration. This blurring is directly linked to the change in the intensity autocorrelation function g₂(t), which in turn is related to sample motion through the field autocorrelation function g₁(t). The goal of most speckle contrast models is to relate changes in speckle contrast to changes in the autocorrelation decay time, τ_c. The autocorrelation time has been said to be inversely related to the speed of the scatterers [7], with multiple scattering theories including a weighting term for each additional dynamic scattering event [8–10]. Sequential intravascular scattering has been studied in the context of diffuse correlation spectroscopy, but remains to be examined for LSCI [11].

Attempts to develop a more thorough, neurovascular-specific understanding of speckle imaging rely on revisiting assumptions used in dynamic light scattering models and including realistic information about the complex spatial structure of the vascular network. Recent work has shown that the degree of intravascular multiple scattering during speckle imaging is lower than the photon diffusion limit and is significantly different when imaging surface vessels and parenchyma [12]. Furthermore, speckle visibility models that incorporate multiple scattering as a function of vessel caliber have been demonstrated to reliably predict the relationship between $\frac{1}{τ_{c}}$ and red blood cell (RBC) velocity changes in surface vessels in vivo [10]. None of these methods have examined the sensitivity of the LSCI signal to flow changes in specific regions of the vascular bed.

The speckle contrast signal results from an ensemble average of all dynamic scattering events experienced by detected photons. Determining the sensitivity of speckle contrast imaging to flow changes therefore requires a description of intravascular scattering not only directly below the detector, but also in every vessel that a detected photon might have traveled through. Using our recently developed method for calculating the autocorrelation function for ordered flow in 3D geometries [13], the sensitivity of LSCI to changes in underlying velocity is examined using Monte Carlo simulations of light scattering in cortical vasculature and a corresponding description of blood velocity in each vessel. Sensitivity is determined by virtually reducing blood flow velocities in the perturbed volume and calculating the resulting change in $\frac{1}{τ_{c}}$ . Differences in sensitivity when imaging over surface vessels and parenchyma are examined. Several different geometries are evaluated to determine the generalizability of the results.

2. Theory

To numerically simulate the sensitivity of LSCI to flow perturbations, the effect of the dynamic scattering of light from moving particles must be considered. Traditionally, the spatial contrast calculation for LSCI has been calculated from the raw speckle image as:

K_{s} = \frac{σ_{s}}{〈 I 〉}

where σ_s is the spatial standard deviation of the intensity and 〈I〉 is the mean intensity of a N × N pixel sliding window [14, 15]. The speckle contrast is a function of both the camera exposure time, T, and the autocovariance, C_t(t), of the fluctuations in intensity at a single pixel. Provided the integration time is sufficiently long, this relationship is demonstrated through the second moment of the speckle contrast value:

K_{s}^{2} = \frac{σ_{s}^{2}}{{〈 I 〉}^{2}} = \frac{1}{T {〈 I 〉}^{2}} \int_{0}^{T} (1 - \frac{t}{T}) C_{t} (t) d t .

The autocovariance, C_t(t), is related to the intensity autocorrelation function through the following relation:

g_{2} (t) = 1 + \frac{C_{t} (t)}{{〈 I 〉}^{2}} .

The quantity g₂(t) is calculated because of its relation to the electric field autocorrelation function, g₁(t), which is related to g₂(t) through the Siegert relation:

g_{2} (t) = 1 + β {| g_{1} (t) |}^{2}

where β is an instrumentation factor that accounts for loss of correlation related to the ratio of the detector size to speckle size, as well as the polarization of the source and detector [16]. It then follows that speckle contrast may be directly related to g₁(t) by:

K_{s}^{2} = \frac{1}{T} \int_{0}^{T} β {| g_{1} (t) |}^{2} (1 - \frac{t}{T}) d t .

The autocorrelation function g₁(t) depends on the ensemble average of the velocities probed by photons traveling through the tissue. Relating speckle contrast to sample velocities therefore critically depends on an accurate model of the particle motion and the amount of photon-particle interactions.

2.1. Velocity distribution

Typically, the form of g₁(t) has been calculated based on an assumed number of scattering events and the type of particle motion. Table 1 shows four forms that are commonly used for dynamic light scattering methods. Note that while the correlation time in each case is defined as τ_c, the interpretation of that value is dependent on the scattering regime and particle motion [23]. For particle motion, unordered (Brownian) motion is characterized by diffusion coefficient D and ordered (linear) motion is characterized by mean squared velocity v_rms. Since Fercher and Briers, the majority of LSCI studies have assumed a Lorentzian line shape, which corresponds to the exponential form of g₁(t) and is appropriate for single dynamic scattering and unordered motion [14].

Table 1. Forms of g₁(τ)

View Table | View all tables in this article

Duncan and others have recently pointed out that the motion probed in biological applications arises primarily from blood flow, which is not unordered [17–19]. They assert that the correct line shape, and therefore the correct form of g₁(t), falls between the limiting behaviors of the Lorentzian (unordered motion) and Gaussian (ordered motion) line shapes.

Much of the discussion on the correction of the form of g₁(t) is based on a single dynamic scattering assumption. However, in our previous work, we have shown that single dynamic scattering is never an accurate assumption of the scattering characteristics in the cerebral cortex [12]. This is likely also true of imaging in other highly vascularized tissues. Therefore, the true form of g₁(t) lies not only between the two single dynamic scattering forms of g₁(t), but also between the single and multiple scattering forms as well.

To add an additional complication, the motion models assume that dynamic scattering events along a single photon path are uncorrelated. This assumption has been justified by asserting that photons scatter only once inside vessels, and then travel outside vessels for a sufficient amount of time for their direction to become randomized before encountering another vessel [14]. However, the photon mean free path length inside blood vessels is 8–16 μm, which is much smaller than many arteries, veins, arterioles and venules [30]. This leads to the breakdown of the un-correlated dynamic scattering assumption because with the high scattering anisotropy in blood vessels, successive dynamic scattering events would in fact be highly correlated. To reformulate a model in which each scattering event along a photon path can be considered independently, we must return to the basic principles behind the formulation of the autocorrelation function g₁(t).

2.2. Simulation of g₁(t)

The dynamic light scattering Monte Carlo (DLS-MC) model was used, which allows g₁(t) to be directly calculated if both the photon paths through a sample are known in addition to the speed and direction of flow at every intravascular scattering location [13]. The form of the calculation is:

g_{1} (t) = 〈 E (0) E^{*} (t) 〉 = \int_{- \infty}^{\infty} P (Y) exp (- j k_{0} Y t) d Y,

where k₀ is the wave number in the sample and Y is the total amount of phase shift experienced by photons undergoing dynamic scattering events. The value of Y depends on the incident and scattering wavevectors, k_i and k_f, and velocity, V, at each scattering event n along a photon’s path:

Y = \sum_{n} (({\hat{k}}_{f, n} - {\hat{k}}_{i, n}) \cdot V_{n}) .

Monte Carlo simulations can be used to determine the photon path through the sample if a geometry that is derived from the sample is used along with appropriate optical properties. Methods such as the recently developed vascular anatomical network (VAN) model, or dynamic light scattering optical coherence tomography (DLS OCT) may be used to determine the vector direction of flow at each point in the sample [20, 21]. With this information, g₁(t) can be directly calculated for an imaging scenario by using the appropriate illumination and detection parameters in the DLS Monte Carlo simulation. The simulated autocorrelation time constant, τ_c, is then determined by finding the time when g₁(t) crosses $\frac{1}{e}$ , which is consistent with the time constant definitions for each of the models shown in Table 1.

2.3. Sensitivity of autocorrelation time to velocity changes

The relationship between changes in velocity and the measured correlation time has not been well established. Using DLS-MC the relationship between the inverse of the autocorrelation time constant $\frac{1}{τ_{c}}$ and velocity may be directly examined. The sensitivity function, s(r), described as

s (r) = \frac{\frac{\frac{\partial 1}{τ_{c}}}{\frac{1}{τ_{c}}}}{\frac{\partial v (r)}{v (r)}} = \frac{v (r)}{\frac{1}{τ_{c}}} \frac{\frac{\partial 1}{τ_{c}}}{\partial v (r)},

relates the change in the width of the correlation function to changes in particle motion at r, and can be determined from DLS-MC by changing the velocity in a given volume about r and then calculating the resulting change in τ_c. The degree of velocity change can be varied to determine the relationship between

\frac{\frac{\partial 1}{τ_{c}}}{\frac{1}{τ_{c}}}

and

\frac{\partial v (r)}{v (r)}

.

2.4. Sensitivity of speckle contrast to $\frac{1}{τ_{c}}$ changes

Experimentally, once the speckle contrast value is obtained, the autocorrelation time is a result of fitting according to the assumed form of g₁(t). Using an assumption of unordered/Brownian motion and single dynamic scattering, $g_{1} (t) = exp (- \frac{t}{τ_{c}})$ , Bandyopadhyay et al. derived the following form of the speckle contrast value:

K^{2} (τ_{c}, T) = β \frac{e^{- 2 x} - 1 + 2 x}{2 x^{2}}

where

x = \frac{T}{τ_{c}}

and T is the exposure duration of the imaging camera [23]. Parthasarathy et al. demonstrated the ability to extract more accurate τ_c values by incorporating the fraction of statically scattered light (1 − ρ), adding terms describing noise (v_n) and non-ergodic variance(v_ne), and imaging over several exposure durations [22]:

K^{2} (τ_{c}, T) = β ρ^{2} \frac{e^{- 2 x} - 1 + 2 x}{2 x^{2}} + 4 β ρ (1 - ρ) \frac{e^{- x} - 1 + x}{x^{2}} + v_{n e} + v_{n} .

Using DLS-MC, the fraction of statically scattered light and the experimental noise terms can be approximated as zero. In this case, the form in Eq. (10) reduces to Eq. (9).

Using the ordered motion, single scattering assumption; or $g_{1} (t) = exp (- {(\frac{t}{τ_{c}})}^{2})$ , results in the following speckle contrast relation [23]:

K^{2} (τ_{c}, T) = β \frac{e^{- 2 x^{2}} - 1 + \sqrt{2 π} x erf (\sqrt{2} x)}{2 x^{2}} .

While the sensitivity of the autocorrelation time to velocity changes is of the most interest, the validity of that relationship with respect to speckle contrast relies on the assumption that the the assumed forms of g₁(t) implicit in the speckle contrast models in Eq. (9) and Eq. (11) do not negatively affect the ability to determine changes in $\frac{1}{τ_{c}}$ . This relationship may be determined with DLS-MC by fitting the simulated g₁(t) values to the assumed forms underlying the common speckle models. While the assumed form of g₁(t) may render the model incapable of accurately capturing the baseline $\frac{1}{τ_{c}}$ value, the change in the fitted value of $\frac{1}{τ_{c}}$ with velocity may be determined and compared to the change in the DLS-MC derived $\frac{1}{τ_{c}}$ .

3. Methods

3.1. Image acquisition and processing

Two photon fluorescence microscopy stacks were taken of FITC-labeled blood plasma in mice (N=5) [24,25]. The image stacks were enhanced using a 3×3×3 median filter and then vectorized using a well-established method developed by Tsai et al. and implemented in the VIDA suite software package [26]. Figure 1(a) shows an example three dimensional rendering of a two-photon stack. A rough outline of the vectorization process follows.

Fig. 1 (a) Three dimensional rendering of mouse microvasculature acquired with two-photon fluorescence microscopy. (b) Segmentation of microvascular map into center-lines and radii, followed by reconstruction.

Download Full Size | PDF

The stacks were first intensity normalized, and then split into separate 100 × 100 × 100 μm sub-blocks for processing as shown in Fig. 1(b). Each sub-block was match filtered using rods to increase the intensity of pixels inside contiguous vessels. The radii of the enhanced gray scale blocks were then thresholded and the radii of each vessel were recorded. The binary vasculature was then shaved down to its centerlines to produce the vector outline of the vasculature. The vector description of the vasculature was recorded by pairing the previously recorded radii with each respective segment of the vasculature; which was defined as the contiguous centerline pixels between two branches. The sub-blocks were then recombined to form a vectorized map of the entire two-photon image stack. This method is quite robust, and has been used in several applications including producing a full 3D network of an entire mouse brain and assisting with automated histology [27, 28]. Extensive manual correction was then performed to ensure all vessel segments were interconnected. Capillary segments which were cut by the imaging field of view were removed. The results of this analysis were the center-line coordinates and radii of each vessel segment that was fully contained within the field of view, which were then used to reconstruct a binary vascular map as shown in the last image of Fig. 1(b).

3.2. Optical properties

The intravascular absorption coefficients were generated based on the extinction coefficients of hemoglobin. The concentration of hemoglobin in the vasculature was assumed to be 2.3 mM [29]. The range of intravascular scattering coefficients were interpolated from hematocrit dependent measurements done by Meinke et al [30]. Extravascular absorption and scattering coefficients were based on the in vitro measurements by Yaroslavsky et al [31]. It was necessary to use in vitro measurements of the extravascular tissue because blood was assumed to only be present in intravascular space. The in vitro measurements were taken in the absence of blood. Table 2 lists the optical properties used in these simulations. Vessels with a diameter less than 11 μm were assumed to be capillaries with a hematocrit (Hct) of 15%, while larger vessels were assumed to have a Hct of 45% [11]. Figure 2(a) shows an example geometry with capillaries and non-capillaries separated by color.

Fig. 2 (a) Example microvascular geometry (725 μm×725 μm laterally, 670 μm deep). Colors represent different optical properties. (b) Three dimensional velocity data generated using a vascular anatomical network model.

Download Full Size | PDF

Table 2. Optical properties for microvasculature geometry

View Table | View all tables in this article

3.3. Blood flow velocity

Using a closed network model of vasculature has been shown by Lorthois et al. to result in accurate flow distributions [32]. The resistence of each vascular segment was calculated using Poiseuilles law corrected for hematocrit as described by Pries et al [33]. The boundary conditions required were the flow speeds of pial arteries and the blood pressure of the pial arteries and veins. Inflowing pial arterial velocities were calculated using the arterial diameter and an assumed perfusion of 100 mL/100g/min. The blood pressure of pial vessels was set using values from Lipowsky et al [34]. Finally, the blood flow distribution was computed using matrix equations from Boas et al [20].

Though the blood flow profile demonstrates variation over the course of each cardiac cycle [35], here the spatial blood flow profiles were assumed to be constant in time and approximately parabolic for vessels larger than a single blood cell:

V (r) = V_{\max} (1 - {(\frac{r}{R})}^{2}) .

where R is the radius of the vessel and r is the distance from the vessel center line. Vessels the size of a single blood cell or smaller were assumed to have a uniform flow profile. Figure 2(b) shows a maximum intensity projection of a 3D parabolized velocity map of one of the five geometries used in this study.

3.4. DLS-MC simulations

For each of the five microvasculature geometries, Monte Carlo simulations were performed using 4 × 10¹⁰ launched photons. A circular collimated beam with a diameter equal to 95% of the geometry width was used to illuminate the surface of each geometry. Several (5–10) 30 μm×30 μm regions of interest (ROI) were chosen in areas with and without surface vessels in each geometry. The coordinates of each scattering event for each detected photon were saved in order to calculate the autocorrelation function as in Eq. (6). The baseline autocorrelation time, τ_cb, was calculated using the DLS-MC method as described in our recent paper [13].

3.4.1. Three dimensional sensitivity

A volumetric study of the sensitivity function $\frac{\frac{\partial 1}{τ_{c}}}{\partial v (r)}$ was performed in one representative surface vessel and parenchyma ROI. To determine the sensitivity of g₁(t) to flow changes in a given volume of the geometry, g₁(t) was recalculated with the speed of all flow within the perturbed volume reduced by 5%. The change in τ_c from the baseline was then calculated as:

Δ \frac{1}{τ_{c}} = 1 - \frac{τ_{cb}}{τ_{cp}}

where τ_cb is the baseline autocorrelation time and τ_cp is the autocorrelation time of g₁(t) calculated when the velocity of the volume of interest is perturbed. The change in τ_c relative to a small change in velocity for a given volume was necessary to determine the sensitivity as defined in Eq. (8).

A sliding 30 μm cube was used as the perturbation volume. The 30 μm cube was scanned over the geometry in three dimensions to determine S(r) for each 30 μm cubic volume in the geometry. Cubic interpolation was used to attain an S(r) value at each 1 μm volume in the geometry and the 3-dimensional S(r) map was masked using the binary geometry to properly visualize the regions which are responsible for the change in $\frac{1}{τ_{c}}$ .

3.4.2. Relationship between velocity change and the resulting $\frac{1}{τ_{c}}$ change

A separate study to determine relationship between the change in velocity and the change in $\frac{1}{τ_{c}}$ was examined using two perturbation volumes for both the parenchyma and surface vessel ROIs. For the parenchyma ROI, the first volume was a cylinder with 50 μm radius centered at the ROI and extending through the entire simulation geometry. This volume was chosen to represent what is generally thought to be the probed volume in parenchyma regions: the depth integrated capillary volume directly below and slightly to the sides of the image ROI. For the surface vessel ROI, the first volume was the entire surface vessel. For both ROIs, the second perturbation volume was the top 200 μm of the geometry, which results in an examination of the sensitivity to surface vasculature over the entire image.

The velocity change in these perturbed regions was varied from 1%–200% of baseline and the resulting change in $\frac{1}{τ_{c}}$ was calculated.

3.4.3. Axial and lateral sensitivity

The sensitivity of g₁(t) to flow changes in both the axial (z) and lateral (x, y) dimensions were calculated using the volume perturbation method. Since each sensitivity calculation is highly dependent on the geometry and the surrounding vasculature, ROI from both surface vessels (N=14) and parenchyma (N=18) were chosen from among the five geometries to allow for generalization of the results. The flow velocity was reduced to 5% of the baseline velocity for each perturbed region.

The axial sensitivity to flow changes was calculated using perturbed volumes which spanned the entire geometry in x and y and were 25 μm deep. The lateral sensitivity was determined using cylindrical volumes that spanned the entire depth of the geometry. The smallest perturbed cylinder had a radius of 25 μm. Each successively larger cylindrical volume had a 25 μm larger radius, where the perturbed volume was the outer 25 μm of the cylinder (i.e. the difference between the inner and outer cylinder).

3.4.4. Accuracy of commonly used speckle models

The ability of the speckle models in Eq. (9) and Eq. (11) to describe changes in $\frac{1}{τ_{c}}$ were examined by calculating the speckle contrast value as in Eq. (5) using 15 exposure durations ranging from 0.05–80 ms and fitting to the models to generate an approximate experimental τ_c value [36]. The changes in the fitted value of $\frac{1}{τ_{c}}$ were compared to the changes in the simulated $\frac{1}{τ_{c}}$ for velocity changes ranging from 1%–200% of baseline as above.

4. Results

The base Monte Carlo simulations of light scattering in each geometry required 200 processors and 7 hours of wall-clock time to run on the TACC Lonestar machine. Each Monte Carlo simulation generated a database of photon scattering positions that could be combined with velocity information in postprocessing for rapid computation of g₁(t) for a given 3D velocity distribution [13].

4.1. Form of g₁(t)

Results showing the form of g₁(t) and the multiple scattering characteristics for both ROIs can be seen in Fig. 3. The simulated autocorrelation function and degree of multiple scattering corresponding to the surface vessel and parenchyma ROIs are shown in black in Figs. 3(a) and 3(b) and Figs. 3(c) and 3(d), respectively. Three forms of g₁(t) are plotted for comparison. The Gaussian line shape shown in green represents a unordered motion and single dynamic scattering assumption. The exponential form in red represents both ordered motion and single scattering, as well as unordered motion in the multiple scattering limit. Finally, the additional form shown in blue represents ordered motion in the multiple scattering limit. Note that a log time scale was used in order to more effectively visualize the relative shape of each curve, although it does exaggerate the error of the fits as opposed to plotting on a linear time scale.

Fig. 3 Surface vessel ROI: (a) Simulated autocorrelation function compared with common forms of g₁(t). (b) Histogram of the number of dynamic scattering events for each collected photon. (c) and (d) show the same data for a parenchyma ROI.

Download Full Size | PDF

The g₁(t) curve from the surface vessel ROI shows two distinct behaviors. The first part of the curve shares a similar shape to the unordered motion, single scattering (exp(−(t/τ)²)) form of g₁(t). The second part of the curve changes shape and lies between the unordered motion/multiple scattering (exp(−t/τ)) and ordered motion/multiple scattering ( $exp (- \sqrt{t / τ})$ ) forms.

The parenchyma ROI also changes form, but less drastically than the surface vessel ROI. The first half of the parenchyma curve decays more rapidly than the exponential case, but only slightly. The longer time behavior of the curve behaves like the surface vessel curve: bounded by the exponential curve and the ordered motion/multiple scattering curve.

The histograms of multiple scattering between the surface vessel and parenchyma ROI are shown in Fig. 3(b) and Fig. 3(d) respectively. The scattering arising from photons traversing through the surface vessel can be clearly seen in the first 5–10 scattering events in the histogram.

4.2. Three dimensional sensitivity of τ_c to changes in particle velocity

The three dimensional sensitivity functions calculated using 30 μm perturbation volumes are shown in Fig. 4. The sensitivity was calculated using a 5% velocity perturbation using Eq. (8). The surface vessel ROI in Figs. 4(a) and 4(b) demonstrates a significant localization of sensitivity to the surface vessel, while the parenchyma ROI in Figs. 4(c) and 4(d) is much more de-localized. From the XZ view in 4(b), the surface vessel also contributes substantially to the depth-dependent sensitivity as it descends into the geometry.

Fig. 4 Three dimensional sensitivity function of surface vessel ROI (a) XY view (b) XZ view, and the parenchyma ROI (c) XY view (d) XZ view. (log scale)

Download Full Size | PDF

While there is a maximum sensitivity in the parenchyma ROI as seen in the capillary just below the surface in 4(d), it is not a sharp maximum and there are several regions both laterally and in depth that display similar sensitivities. Additionally, the sensitivity to lateral surface vasculature in the parenchyma ROI is higher in all surface vessels except the vessel where the surface vessel ROI is placed.

4.3. Sensitivity of $\frac{1}{τ_{c}}$ as a function of velocity change

Figure 5 shows the change in $\frac{1}{τ_{c}}$ for a change in velocity for the surface vessel and parenchyma ROI. The sensitivity as described in Eq. (8) is represented as the slope of the line, with a constant slope over a range of velocities indicating a linear relationship. The targeted regions in Fig. 5(a) were set to represent the region that is being investigated: the entire surface vessel was perturbed for the surface vessel ROI, and a 25 μm cylinder centered at the ROI was perturbed for the parenchma ROI. Figure 5(a) shows that for targeted perturbation regions, both the surface vessel and parenchyma ROI demonstrate nonlinear behavior as the velocity decreases from baseline to zero. However, for small decreases from baseline or for increases from baseline, the relationship between velocity and $\frac{1}{τ_{c}}$ is approximately linear.

Fig. 5 The relationship between $\frac{1}{τ_{c}}$ and particle velocity in a representative surface vessel and parenchyma ROIs. The slope represents the sensitivity function S(r) and the dotted line represents a sensitivity of 1. The plot in (a) shows targeted perturbation regions: the velocity in the entire surface vessel was adjusted for the surface vessel ROI, and a 25 μm radius cylinder was perturbed for the parenchyma ROI. In (b) the top 200 μm of the geometry was perturbed for both ROI.

Download Full Size | PDF

In the perturbed 200 μm layer shown in Fig. 5(b), the surface vessel region displays only slight nonlinearity as the velocity decreases from baseline to zero, but again shows approximately linear behavior for small decreases from baseline and for increases from baseline. The parenchyma ROI demonstrates the same trend, but the nonlinearity is more pronounced, likely due to a more significant amount of signal originating from below the 200 μm layer

4.4. Axial and lateral sensitivity

The results from the axial sensitivity study are shown in Fig. 6(a) and 6(b). Sensitivities for these figures were calculated based on a 95% perturbation of baseline velocities (i.e. 5% of baseline) using Eq. (8). The surface vessel average sensitivity and standard error were calculated using 14 distinct surface vessel ROI from among the five microvasculature geometries. The parenchyma calculations used 18 different parenchymal ROI. The standard error of the top 150 μm of the surface vessel ROI in Fig. 6(a) is significantly larger than the error attributed to any other depth of the tissue regardless of ROI type. This is the result of highly variable surface vessel radii and depths. Due to the higher number of scattering events, surface vessels with large radii produced larger S(r) values than those with smaller radii. Many surface vessels had little sensitivity in the top 50 μm because they were completely contained in the layer 100 μm into the geometry.

Fig. 6 Average sensitivity and standard error of (a) surface vessel ROI (N=14) and (b) parenchyma ROI (N=18) to changes in particle velocity in each 50 μm layer of the geometry. The bottom row shows the lateral sensitivity of (c) surface vessel ROI (N=14) and (d) parenchyma ROI (N=18) to changes in particle velocity in concentric cylinders with radii increasing by 25 μm.

Download Full Size | PDF

The average lateral sensitivity is shown in Fig. 6(c) and 6(d). The sensitivities were calculated using 95% perturbations and concentric cylinders, the smallest of which had a radius of 25 μm. Each successively larger cylinder had a 25 μm larger radius. The surface vessel sensitivities in Fig. 6(c) are significantly more localized around the ROI than the parenchyma ROI in Fig. 6(d). The sensitivity of the inner 25 μm radius cylinder is smaller than the successive concentric cylinders partly because of the differences in the perturbed volumes.

4.5. Sensitivity of K to changes in τ_c

The charts in Fig. 7 show the accuracy of the assumed forms of g₁(t). Figure 7(a) shows the same graph of inverse autocorrelation time change versus velocity shown in Fig. 5, but includes the change in $\frac{1}{τ_{c}}$ as predicted by the three assumed forms of g₁(t) shown in Table 2. Figure 7(b) shows the relative baseline values for 32 ROIs across the five geometries. The baseline $\frac{1}{τ_{c}}$ values derived from fitting the three assumed forms of g₁(t) were divided by the $\frac{1}{τ_{c}}$ values from the simulated g₁(t).

Fig. 7 The accuracy of typically assumed forms of g₁(t). (a) is the same graph of sensitivity to flow perturbations shown in 5(a), but with matching $\frac{1}{τ_{c}}$ values extracted using the models shown in Table 2. (b) shows the accuracy of the assumed forms of g₁(t) at determining the true $\frac{1}{τ_{c}}$ value in 32 ROIs.

Download Full Size | PDF

Figure 7(a) demonstrates that all three assumed forms of g₁(t) are nearly equal in their ability to extract the change in $\frac{1}{τ_{c}}$ of the true underlying autocorrelation function. The results show that using the assumed forms does slightly under-predict the change in $\frac{1}{τ_{c}}$ when the changes are large, such as in the surface vessel ROI or in the parenchyma ROI when the velocity is close to zero.

The comparison of baseline $\frac{1}{τ_{c}}$ values in Fig. 7(b) demonstrate that the choice of g₁(t) form significantly affects the baseline values of the autocorrelation time. Using the form g₁(t) = exp(−(t/τ_c)²) which represents single scattering and ordered motion, appeared to strongly under-predict the baseline $\frac{1}{τ_{c}}$ value. In contrast, the form g₁(t) = exp(−(t/τ_c)^1/2), which represents multiple scattering and diffusive motion, over-predicted the $\frac{1}{τ_{c}}$ by a similar margin. The exponential form, while still under-predicting the baseline inverse autocorrelation time value, provided the closest match to the simulated $\frac{1}{τ_{c}}$ baseline value.

5. Discussion

5.1. Correct form of g₁(t)

The assumed form of g₁(t) in tissues has historically depended on the two factors in Table 1: particle motion and the degree of scattering. These two factors are sufficient for establishing the limits of the expected behavior, but none of the forms in themselves describe the simulated forms of g₁(t) shown in Fig. 3.

The simulated g₁(t) curve in both the surface vessel and parenchyma ROI displayed similar behavior. In the shorter time scales, they were bound by the g₁(t) form expected from singly scattered/ordered motion and multiply scattered/ordered motion forms. This is expected primarily because the DLS-MC method only simulates ordered motion, so despite g₁(t) associated with experimental measurements having sampled diffusive motion, here the only contribution is from photons interacting with particles moving at constant velocity. This begs an explanation from the longer time behavior. For both ROI, the g₁(t) curves are bound by the exponential function and the form associated with multiple scattering and unordered motion.

This behavior is the result of a further breakdown of the assumptions implicit in the forms of g₁(t). The forms which describe ordered motion and single scattering are based on a Gaussian line shape. If this assumption is accurate, then the scattering vectors and velocity vectors that together determine the electric field phase shift for a dynamic scattering event are decoupled from each other and the distribution has a zero mean value. This drastically simplifies the analytical form and results in the commonly used second-moment based g₁(t) forms, which are characterized by the mean squared displacement.

However, when performing LSCI measurements in the cortex, the velocities encountered are varied but not random. Capillaries tend to have similar velocities, but surface arteries, veins, arterioles, and venules all have highly varied flow speeds. Photons sample multiple vessels along their path, each with a different velocity and orientation, and often interact with vessels whose velocities are outliers with respect to surrounding vasculature. These effects taken together cause a breakdown in the Gaussian line shape assumption and render the associated forms incapable of accurately describing the particle dynamics.

5.2. Non-linearity of sensitivity with velocity

The results in Fig. 5 show non-linearity as the velocity is reduced from baseline to zero for both the surface vessel and parenchyma ROIs. However, given how g₁(t) arises, this should come as no surprise. The contribution to g₁(t) from each dynamic scattering event depends on q · v as it travels through the geometry.

Since the measured autocorrelation function is an ensemble average of all photons, the contribution from one particular region in the tissue contributes in two ways. First, all photons which have traveled through the region and then been detected will have different phases. Any change resulting from these photons will therefore only impact the fraction of the ensemble average that they represent. Secondly, of the photons that do travel through the perturbed region, the fraction of the total phase does not decline linearly since the phase shifts from unperturbed regions are still included. Since both these effects are nonlinear, the dependence of τ_c on flow velocity in any given vessel or region is inherently nonlinear when imaging in the cortex.

It is important to point out that the underlying autocorrelation function from which the speckle contrast value is calculated is never dependent on a single vessel. A decrease in flow in a surface vessel, for example, should not be expected to affect the autocorrelation time in a linear manner. Instead, the process of the vessel velocity changing from baseline to zero in a single vessel is reflected by a shift in the ensemble average of the dynamic phase shifts to a new ensemble average representing all the other sampled vessels. The remaining vessels’ weights in the ensemble average are higher and therefore their sensitivity is higher than it was when the surface vessel was present. This is why the axial sensitivity functions of the surface vessel and parenchyma demonstrate such a stark contrast.

5.3. Sensitivity to flow perturbations

The most immediately obvious conclusion from the 3D sensitivity and the depth-based axial sensitivity studies is that surface vessels and parenchyma ROIs have very different sensitivity. The sensitivity function in the surface vessel ROIs are very spatially confined to the vessel it represents. The log scale 3D sensitivity map in Fig. 4(b) and 4(c) demonstrates this clearly: the sensitivity to the vessel is over two orders of magnitude stronger than sensitivity to any other vasculature. The parenchyma ROIs on the other hand display a more uniform sensitivity function, with τ_c being nearly equally sensitive to capillaries and descending vessels beneath the surface, as well as to surface vessels lateral to the ROI. The parenchyma ROI results show that for any particular capillary region (25 μm cube), the sensitivity of $\frac{1}{τ_{c}}$ to changes in velocity is very small: between 0.1–1% change for a 100% change in velocity. This property effectively precludes LSCI from being usable to detect changes on an individual capillary level.

The differences between surface vessel and parechyma ROIs are also reflected in the average axial and lateral sensitivities in Fig. 6. The strong localization of surface vessel ROIs result in large but highly variable depth dependent sensitivity in the first 150 μm, followed by very little sensitivity to deeper vasculature. The more uniform parenchyma ROIs also display higher sensitivities to surface vessels, but the variation is smaller because nearly all surface vasculature lateral to the ROI contributes to the sensitivity. The parenchyma ROIs also display higher sensitivity to lower depths relative to the surface vessel ROI.

As alluded to in the previous section, these effects are ultimately the result of the ensemble average that generates the g₁(t) curve. In surface vessels, most detected photons must travel through the vessel to reach the detector. This results in a large contributor to the ensemble average that minimizes the effect of other vasculature on the resulting behavior of g₁(t). The photons reaching a parenchyma ROIs are not strongly weighted towards any specific vessel, which allows them greater sensitivity to changes deeper into the tissue as well as in laterally positioned surface vessels.

5.4. Speckle contrast sensitivity

The sensitivity of speckle contrast calculations are ultimately dependent on whether the assumed form of g₁(t) implicit in the model can accurately represent changes in the autocorrelation time of the true underlying g₁(t) curve. Fortunately, as shown in Fig. 7(a) it appears that using any of the three most commonly assumed forms results in almost the same sensitivity to underlying flow as the DLS-MC calculated g₁(t) function. This implies that the results and discussion on the sensitivity of g₁(t) to flow changes in the cortex are equally applicable to speckle contrast calculations.

Though the sensitivity to flow changes might be the same, in order to connect speckle contrast calculations to absolute flow values in the tissue, the models must be able to accurately extract the true baseline $\frac{1}{τ_{c}}$ value. This was a primary criticism advanced by Duncan et al. regarding attempts to make LSCI more quantitative [17]. In this case, the results in Fig. 7(b) show that none of the assumed forms are capable of determining the correct baseline $\frac{1}{τ_{c}}$ . However, the assumption of an exponential form of g₁(t) comes the closest. Among the 32 ROIs examined, the error in the baseline $\frac{1}{τ_{c}}$ for the exponential form was on average 12% and at most slightly over 20%.

6. Conclusion

Using dynamic light scattering Monte Carlo simulations, the sensitivity of the autocorrelation function and speckle contrast to flow changes in the cortex was extensively examined. It was found that the commonly used forms of g₁(t), while they happen to provide limiting behaviors for the shape of the curve, rely on assumptions regarding the type and amount of scattering that are not accurate.

Furthermore, the sensitivity studies demonstrate that the sensitivity of τ_c to changes in a surface vessel region of interest (ROI) is strongly localized to the single vessel, while parenchymal ROI have a larger sensitivity to flow changes deeper into the tissue. The axial sensitivity trends hold in the lateral sensitivity case, with surface vessel sensitivity being most strongly localized to the surrounding 25 μm, but parenchymal ROI demonstrating significant sensitivity to vasculature up to 200 μm or more lateral to the ROI.

Additionally, it was shown that utilizing the commonly used speckle contrast models resulted in nearly the same sensitivity to underlying flow. This demonstrates that the sensitivity results shown in this paper are applicable to speckle contrast images processed using the most commonly assumed forms of g₁(t). However, none of the speckle contrast models were able to accurately extract the baseline $\frac{1}{τ_{c}}$ . While this is currently a major limitation preventing the accurate quantification of absolute flow, it is shown that assuming an exponential form of g₁(t) results in the smallest error of the three examined forms.

Acknowledgments

We gratefully acknowledge support from the Consortium Research Fellows Program, NIH (EB-011556, NS-078791, NS-082518), American Heart Association (14EIA8970041), Coulter Foundation, and the ORISE Research Participation Program. The authors also acknowledge the Texas Advanced Computing Center (TACC) at the University of Texas at Austin for providing HPC resources that have contributed to the research results reported within this paper. URL: http://www.tacc.utexas.edu.

References and links

1. G. A. Armitage, K. G. Todd, A. Shuaib, and I. R. Winship, “Laser speckle contrast imaging of collateral blood flow during acute ischemic stroke,” J. Cereb. Blood Flow Metab. 30, 1432–1436 (2010). [CrossRef] [PubMed]

2. S. M. S. Kazmi, L. M. Richards, C. J. Schrandt, M. A. Davis, and A. K. Dunn, “Expanding applications, accuracy, and interpretation of laser speckle contrast imaging of cerebral blood flow,” J. Cereb. Blood Flow Metab. in press (2015). [CrossRef] [PubMed]

3. A. K. Dunn, A. Devor, H. Bolay, M. L. Andermann, M. A. Moskowitz, A. M. Dale, and D. A. Boas, “Simultaneous imaging of total cerebral hemoglobin concentration, oxygenation, and blood flow during functional activation,” Opt. Lett. 28, 28–30 (2003). [CrossRef] [PubMed]

4. N. Hecht, J. Woitzik, J. P. Dreier, and P. Vajkoczy, “Intraoperative monitoring of cerebral blood flow by laser speckle contrast analysis,” Neurosurg. Focus 27, E11 (2009). [CrossRef] [PubMed]

5. A. B. Parthasarathy, E. L. Weber, L. M. Richards, D. J. Fox, and A. K. Dunn, “Laser speckle contrast imaging of cerebral blood flow in humans during neurosurgery: a pilot clinical study,” J. Biomed. Opt. 15, 066030 (2010). [CrossRef]

6. L. M. Richards, E. L. Towle, D. J. Fox, and A. K. Dunn, “Intraoperative laser speckle contrast imaging with retrospective motion correction for quantitative assessment of cerebral blood flow,” Neurophotonics 1, 015006 (2014). [CrossRef]

7. R. Bonner and R. Nossal, “Model of laser Doppler measurements of blood flow in tissue,” Appl. Opt. 20, 2097–2107 (1981). [CrossRef] [PubMed]

8. D. A. Boas and A. G. Yodh, “Spatially varying dynamical properties of turbid media probed with diffusing temporal light correlation,” J. Opt. Soc. Am. A 14, 192–215 (1997). [CrossRef]

9. D. Pine, D. Weitz, P. Chaikin, and E. Herbolzheimer, “Diffusing wave spectroscopy,” Phys. Rev. Lett. 60, 1134–1137 (1988). [CrossRef] [PubMed]

10. S. M. S. Kazmi, E. Faraji, M. A. Davis, Y.-Y. Huang, X. J. Zhang, and A. K. Dunn, “Flux or speed? examining speckle contrast imaging of vascular flows,” Biomed. Opt. Express 6, 2588–2608 (2015). [CrossRef] [PubMed]

11. S. A. Carp, N. Roche-Labarbe, M.-A. Franceschini, V. J. Srinivasan, S. Sakadžić, and D. A. Boas, “Due to intravascular multiple sequential scattering, Diffuse Correlation Spectroscopy of tissue primarily measures relative red blood cell motion within vessels,” Biomed. Opt. Express 2, 2047–2054 (2011). [CrossRef] [PubMed]

12. M. A. Davis, S. M. S. Kazmi, and A. K. Dunn, “Imaging depth and multiple scattering in laser speckle contrast imaging,” J. Biomed. Opt. 19, 086001 (2014). [CrossRef] [PubMed]

13. M. A. Davis and A. K. Dunn, “Dynamic light scattering Monte Carlo: a method for simulating time-varying dynamics for ordered motion in heterogeneous media,” Opt. Express 23, 17145–17155 (2015). [CrossRef] [PubMed]

14. A. Fercher and J. Briers, “Flow visualization by means of single-exposure speckle photography,” Opt. Commun. 37, 326–330 (1981). [CrossRef]

15. J. Briers and S. Webster, “Laser speckle contrast analysis (LASCA): a nonscanning, full-field technique for monitoring capillary blood flow,” J. Biomed. Opt. 1, 174–179 (1996). [CrossRef] [PubMed]

16. P. Lemieux and D. Durian, “Investigating non-Gaussian scattering processes by using n th-order intensity correlation functions,” J. Opt. Soc. Am. A 16, 1651–1664 (1999). [CrossRef]

17. D. D. Duncan, S. J. Kirkpatrick, and R. K. Wang, “Statistics of local speckle contrast,” J. Opt. Soc. Am. A 25, 9 (2008). [CrossRef]

18. J. Ramirez-San-Juan, “Impact of velocity distribution assumption on simplified laser speckle imaging equation,” Opt. Express 16, 3197–3203 (2008). [CrossRef] [PubMed]

19. D. Briers, D. D. Duncan, E. Hirst, S. J. Kirkpatrick, M. Larsson, W. Steenbergen, T. Stromberg, and O. B. Thompson, “Laser speckle contrast imaging: theoretical and practical limitations,” J. Biomed. Opt. 18, 066018 (2013). [CrossRef] [PubMed]

20. D. Boas, S. Jones, A. Devor, T. Huppert, and A. Dale, “A vascular anatomical network model of the spatio-temporal response to brain activation,” Neuroimage 40, 1116–1129 (2008). [CrossRef] [PubMed]

21. J.-Y. Lee, M.-P. Lu, K.-Y. Lin, and S.-H. Huang, “Measurement of in-plane displacement by wavelength-modulated heterodyne speckle interferometry,” Appl. Opt. 51, 1095–1100 (2012). [CrossRef] [PubMed]

22. A. B. Parthasarathy, W. J. Tom, A. Gopal, X. Zhang, and A. K. Dunn, “Robust flow measurement with multi-exposure speckle imaging,” Opt. Express 16, 1975–1989 (2008). [CrossRef] [PubMed]

23. R. Bandyopadhyay and A. Gittings, “Speckle-visibility spectroscopy: A tool to study time-varying dynamics,” Rev. Sci. Instrum. 76, 093110 (2005). [CrossRef]

24. L. Gagnon, S. Sakadžić, F. Lesage, E. T. Mandeville, Q. Fang, M. A. Yaseen, and D. A. Boas, “Multimodal reconstruction of microvascular-flow distributions using combined two-photon microscopy and doppler optical coherence tomography,” Neurophotonics 2, 015008 (2015). [CrossRef] [PubMed]

25. L. Gagnon, S. Sakadžić, F. Lesage, J. J. Musacchia, J. Lefebvre, Q. Fang, M. A. Yücel, K. C. Evans, E. T. Mandeville, J. Cohen-Adad, et al., “Quantifying the microvascular origin of bold-fmri from first principles with two-photon microscopy and an oxygen-sensitive nanoprobe,” J. Neurosci. 35, 3663–3675 (2015). [CrossRef] [PubMed]

26. P. S. Tsai, J. P. Kaufhold, P. Blinder, B. Friedman, P. J. Drew, H. J. Karten, P. D. Lyden, and D. Kleinfeld, “Correlations of neuronal and microvascular densities in murine cortex revealed by direct counting and colocalization of nuclei and vessels,” J. Neurosci. 29, 14553–14570 (2009). [CrossRef] [PubMed]

27. D. Kleinfeld, A. Bharioke, P. Blinder, D. D. Bock, K. L. Briggman, D. B. Chklovskii, W. Denk, M. Helmstaedter, J. P. Kaufhold, W.-C. A. Lee, H. S. Meyer, K. D. Micheva, M. Oberlaender, S. Prohaska, R. C. Reid, S. J. Smith, S. Takemura, P. S. Tsai, and B. Sakmann, “Large-scale automated histology in the pursuit of connectomes,” J. Neurosci. 31, 16125–16138 (2011). [CrossRef] [PubMed]

28. P. Blinder, P. S. Tsai, J. P. Kaufhold, P. M. Knutsen, H. Suhl, and D. Kleinfeld, “The cortical angiome: an interconnected vascular network with noncolumnar patterns of blood flow,” Nat. Neurosci. 16889–897 (2013).

29. F. E. Robles, S. Chowdhury, and A. Wax, “Assessing hemoglobin concentration using spectroscopic optical coherence tomography for feasibility of tissue diagnostics,” Biomed. Opt. Express 1, 310–317 (2010). [CrossRef]

30. M. Meinke, G. Müller, J. Helfmann, and M. Friebel, “Empirical model functions to calculate hematocrit-dependent optical properties of human blood,” Appl. Opt. 46, 1742–1753 (2007). [CrossRef] [PubMed]

31. A. N. Yaroslavsky, P. C. Schulze, I. V. Yaroslavsky, R. Schober, F. Ulrich, and H. J. Schwarzmaier, “Optical properties of selected native and coagulated human brain tissues in vitro in the visible and near infrared spectral range,” Phys. Med. Biol. 47, 2059–2073 (2002). [CrossRef] [PubMed]

32. S. Lorthois, F. Cassot, and F. Lauwers, “Simulation study of brain blood flow regulation by intra-cortical arterioles in an anatomically accurate large human vascular network. Part II: flow variations induced by global or localized modifications of arteriolar diameters,” Neuroimage 54, 2840–2853 (2011). [CrossRef]

33. A. R. Pries, T. W. Secomb, P. Gaehtgens, and J. F. Gross, “Blood flow in microvascular networks. Experiments and simulation,” Circ. Res. 67, 826–834 (1990). [CrossRef] [PubMed]

34. H. H. Lipowsky, “Microvascular rheology and hemodynamics,” Microcirculation 12, 5–15 (2005). [CrossRef] [PubMed]

35. T. P. Santisakultarm, N. R. Cornelius, N. Nishimura, A. I. Schafer, R. T. Silver, P. C. Doerschuk, W. L. Olbricht, and C. B. Schaffer, “In vivo two-photon excited fluorescence microscopy reveals cardiac- and respiration-dependent pulsatile blood flow in cortical blood vessels in mice,” Am. J. Physiol. Heart Circ. Physiol. 302, H1367–H1377 (2012). [CrossRef] [PubMed]

36. S. M. S. Kazmi, S. Balial, and A. K. Dunn, “Optimization of camera exposure durations for multi-exposure speckle imaging of the microcirculation,” Biomed. Opt. Express 5, 2157–2171 (2014). [CrossRef] [PubMed]

g₁(τ) =	Light scattering	Particle motion
exp( $- \frac{t}{τ_{c}}$ )	Single	Unordered
exp( $- {(\frac{t}{τ_{c}})}^{2}$ )	Single	Ordered
exp( $- \frac{t}{τ_{c}}$ )	Multiple	Ordered
exp( $- \sqrt{\frac{t}{τ_{c}}}$ )	Multiple	Unordered

	μ_a (mm⁻¹)	μ_s (mm⁻¹)	g
Capillaries	0.20	60	0.98
Non-capillaries	0.20	95	0.98
Extravascular	0.02	10	0.90

g₁(τ) =	Light scattering	Particle motion
exp( $- \frac{t}{τ_{c}}$ )	Single	Unordered
exp( $- {(\frac{t}{τ_{c}})}^{2}$ )	Single	Ordered
exp( $- \frac{t}{τ_{c}}$ )	Multiple	Ordered
exp( $- \sqrt{\frac{t}{τ_{c}}}$ )	Multiple	Unordered

	μ_a (mm⁻¹)	μ_s (mm⁻¹)	g
Capillaries	0.20	60	0.98
Non-capillaries	0.20	95	0.98
Extravascular	0.02	10	0.90

Sensitivity of laser speckle contrast imaging to flow perturbations in the cortex

Abstract

1. Introduction