Characterization of continuous wave ultrasound for acousto-optic modulated diffuse correlation spectroscopy (AOM-DCS)

: Intra and post-operative blood ﬂow monitoring of tissue has been shown to be eﬀective in the improvement of patient outcomes. Diﬀuse correlation spectroscopy (DCS) has been shown to be eﬀective in measuring blood ﬂow at the bedside, and is a useful technique in measuring cerebral blood ﬂow (CBF) in many clinical settings. However, DCS suﬀers from reduced sensitivity to blood ﬂow changes at larger tissue depths, making measurements of CBF in adults diﬃcult. This issue can be addressed with acousto-optic modulated diﬀuse correlation spectroscopy (AOM-DCS), which is a hybrid technique that combines the sensitivity of DCS to blood ﬂow with ultrasound resolution to allow for improved spatial resolution of the optical signal based on knowledge of the area which is insoniﬁed by ultrasound. We present a quantitative model for perfusion estimation based on AOM-DCS in the presence of continuous wave ultrasound, supported by theoretical derivations, Monte Carlo simulations, and phantom and human subject experiments. Quantiﬁcation of the inﬂuence of individual mechanisms that contribute to the temporal ﬂuctuations of the optical intensity due to ultrasound is shown to agree with previously derived results. By using this model, the recovery of blood-ﬂow induced scatterer dynamics based on ultrasound-modulated light is shown to deviate by less than one percent from the standard DCS measurement of scatterer dynamics over a range of optical scattering values and scatterer motion conditions. This work provides an important step towards future implementation of AOM-DCS setups with more complex spatio-temporal distributions of ultrasound pressure, which are needed to enhance the DCS spatial resolution.


Introduction
Hemodynamic monitoring of patients in the intra and post-operative periods has been shown to be effective in guiding treatment and reducing negative outcomes such as organ failure [1]. Guiding treatment based on measurements of lactate concentration, central venous oxygen saturation, and extraction of oxygen were shown to improve patient outcomes. A tool that can allow for the quantification of these parameters continuously at the bedside could be extremely helpful in guiding therapeutic interventions. Diffuse optical techniques have been shown to allow for non-invasive monitoring of tissue at the bedside for these relevant parameters [2][3][4]. One such technique is diffuse correlation spectroscopy (DCS). Developed in the 1990's, DCS is an optical technique that quantifies the blood flow in tissue through the analysis of the temporal evolution of speckle intensity [5]. Laser light with a long temporal coherence length is launched into the tissue, and a speckle pattern is generated at the surface of the skin. The temporal evolution of the speckle intensity allows the interrogation of dynamic motion in the tissue through the analysis of the intensity temporal autocorrelation function, known in the literature as g 2 (τ). This function has a characteristic decay that is due to the motion of the scattering particles in the tissue, dominated by red blood cells (RBC), and allows for the quantification of perfusion in terms of a blood flow index (BF i ) [6]. Hence DCS can quantify the motion of RBCs in the microvasculature, but as with near-infrared spectroscopy (NIRS) and other optical methods, it loses sensitivity to RBC motion with increased depth [7]. This limitation is particularly relevant when accurate measures of cerebral blood flow are desired, where the scalp and skull can have a thickness of 1 cm or more and changes in brain blood flow are largely masked by concurrent changes in the superficial tissue blood flow in the scalp, making it difficult to extract meaningful changes relevant to cerebral physiology. Here, we utilize a multi-modal approach through the combination of ultrasound tagging and diffuse correlation spectroscopy, known as acousto-optic modulated DCS (AOM-DCS), in an attempt to overcome this limitation and improve depth sensitivity. Past work has shown that the combined use of ultrasound and light is useful in the quantification of cerebral blood flow using a cross-correlation technique of input ultrasound pressure and the modulation of the measured speckle intensity [8]. This technique is effective in quantifying relative flow differences at different depths in tissue, though the use of the cross-correlation allows for only a single correlation parameter to be calculated at each depth, and doesn't allow for the extraction of the entire autocorrelation function, g 2 (τ). Other work combining DCS and ultrasound has shown quantification of fluid flow in a capillary tube embedded in a scattering phantom in the transmission geometry [9], as well as the effects of acoustic radiation force on the diffuse correlation spectroscopy signal [10]. In this work, we extract the entire autocorrelation function of the tagged light from the modulation present in the intensity autocorrelation function and apply DCS theory to the analysis of the tagged light signal to quantitatively determine the scatterer dynamics along the paths of tagged photons. As a first step, we investigate the interaction between continuous wave ultrasound and the speckle fluctuations in the DCS signals, through the use of theoretical derivations, Monte Carlo simulations, experiments with tissue mimicking phantoms, and in vivo pilot experiments. Comparisons to previous derivations of the influence of different mechanisms of ultrasound modulation of light in a scattering media show good agreement in the simulation results, and BF i extracted from the ultrasound-modulated light signal matches the BF i measured by traditional DCS. Though the full benefits of increased spatial localization given by ultrasound tagging of light are not realized in this work, as continuous wave ultrasound is used, this work acts as a basis for the future development of the AOM-DCS setups with more complex ultrasound pressure distributions, like those seen in focused or pulsed ultrasound.

DCS theory without an ultrasound field
DCS is a technique that is intrinsically sensitive to alterations in the interference of light, and it is useful in sensing the motion of optical scatterers, mainly red blood cells, in biological tissue [11,12]. In DCS, the blood flow is quantified by estimating the blood flow index (BF i ) from the measured temporal autocorrelation of the optical intensity [5]. Equation (1) is used to estimate the autocorrelation of the electric field (e.g., g 1 (τ) = E(t)E * (t + τ) ) under the weak-scattering approximation that the optical mean-free path is much greater than the optical wavelength, in which case the transfer of light can be described by ladder diagrams [13]. By following the diffusing wave spectroscopy (DWS) approach [14], the following expression for g 1 (τ) was obtained under the assumption that only scatterer motion contributes to the temporal variation of the electric field [15]: In Eq. (1), P(s) is the probability density of a light pathlength s, µ s ' is the reduced scattering coefficient, k is the wavenumber of the light in the scattering media, and ∆r 2 (τ) is the mean squared displacement during the time τ. For an ergodic sample, from the measured intensity, we can relate the normalized intensity autocorrelation, g 2 (τ) = I(t)I(t + τ) , to the normalized electric field autocorrelation, g 1 (τ), using the Siegert relation [16] The term of interest, which generates the temporal decay of the autocorrelation, is the mean squared displacement, ∆r 2 (τ) . Throughout the development of DCS, the most suitable form of this parameter in tissue has been found to be one resulting from a diffusive motion in a multiple scattering regime, given by ∆r 2 (τ) = 6αD eff τ, where D eff is the effective scatterer diffusion coefficient and α is the fraction of scattering events that occur from moving scatterers. [17,18] The measured g 2 (τ) is then fit to obtain the blood flow index (BF i ) parameter, defined as BF i = αD eff . The BF i parameter has been found to correlate well with the blood flow measured with ASL-MRI, Xenon-CT, fluorescent microspheres, PET, and transcranial doppler ultrasound [19][20][21][22][23].

Ultrasound-induced changes of the optical phase along the photon paths
When a DCS measurement is performed in a scattering media in the presence of ultrasound, the temporal autocorrelation function of the optical intensity, g 2 (τ), includes a component that depends on the ultrasound induced motion of the scatters [24]. Furthermore, as light travels through the tissue, ultrasound-induced periodic modulation of the index of refraction causes temporal changes of the optical phase [25]. The combination of these two mechanisms of ultrasound modulation of the optical waves yields a g 2 (τ) that carries modulation at the ultrasound frequency. If the optical scatterers are moving only due to presence of the continuous ultrasound, the amplitude of the g 2 (τ) modulation is not expected to change with time. However, if the optical scatterers in addition to ultrasound-induced oscillation also exhibit diffusive and/or advective motions, which is expected in the perfused biological tissue, then additional temporal decorrelation will occur, and both g 2 (τ) and the amplitude of the g 2 (τ) modulation by the ultrasound will decay with time. Due to similar decorrelation mechanisms involved in the decay of both g 2 (τ) and its modulation amplitude, it should be possible to use either one of them to quantify the same stochastic properties of the scatterers motion and make inferences about the blood flow. However, the quantification of the decay of the g 2 (τ) modulation may have an advantage over the standard g 2 (τ) decay as it can be used to measure the flow based on the 'tagged photons', which are spatially and/or temporally localized to the interaction region between the light and the ultrasound. Here, our goal is to establish a first step towards developing a DCS setup that is effective in making quantitative measurements of blood flow enabled by spatially and temporally heterogenous ultrasound pressure distributions. In order to characterize the effects of ultrasound on the DCS signal, here we start by first considering the continuous wave ultrasound, and propose a theory to describe the interaction between ultrasound and light in the context of DCS measurements. For brevity, an abbreviated derivation will be presented below with a more complete summary given in the appendix.
The contributions of tagged photons to the modulation present in the autocorrelation function are reliant on the degree of interaction that they have with the ultrasound pressure. For a continuous plane wave ultrasound, the ultrasound pressure at a particular position and time is given by P(r, t) = P 0 cos(ω u t − k u · r + φ), where P 0 is the pressure amplitude of the ultrasound, ω u is the angular frequency of the ultrasound, k u is the ultrasound wavevector with magnitude k u , and φ is the relative phase shift of the ultrasound wave. For simplicity, the ultrasound pressure attenuation term has been omitted in this derivation, though it is included in the Monte Carlo simulations. The ultrasound pressure manipulates the position of scattering particles as well as modulates the index of refraction of the media, which results in phase shifts of the electric field along the photon paths, given below, respectively for displacement modulation (Eq. (3)) and index of refraction modulation (Eq. (4)), Where ∆φ d and ∆φ n are the modulations in phase due to displacement and index of refraction modulation due to the modulation in position, ∆r i , and index of refraction, ∆n(r, t), respectively, q i is the momentum transfer at a particular scattering center, and k 0 is the wavenumber of the light in vacuum. Finally, if we assume that blood flow-induced motion of the RBCs can be modeled as a random walk, the total optical phase increment along the photon path, ∆φ T = ∆φ US + ∆φ B , will include the components due to ultrasound, ∆φ US = ∆φ d + ∆φ n , and the Brownian motion [24,[26][27][28]. A full summary expanding the terms given above can be found in the appendix.

Electric field and optical intensity temporal autocorrelation functions in the presence of both blood flow and ultrasound field
When the phase variation ∆φ T (t) along the path is much less than unity, we can simplify the calculation of g 1 (τ, s) along the path of length s as [14] where F s (τ) = ∆φ 2 T (t, τ) s , ∆φ T (t, τ) = ∆φ T (t + τ) − ∆φ T (t), and averaging s is performed over time t and all configurations of the pathlength with length s [28,29]. Here, we provide the simplified solution that assumes isotropic scattering, to which the anisotropic case reduces over pathlength distances larger than the transport mean free path (l tr ) via the similarity relation as [30]: where F B,s (τ) = 4n 2 0 k 2 0 D B s l tr τ is the term due to the blood flow and F US,s (τ) = W US,s sin 2 ω u τ 2 is the term due to ultrasound modulation. where In Eq. (7), G = arctan(k u l tr ) k u l tr , and we assumed that k u l tr and sl −1 tr are much greater than one, such that higher order terms, not proportional with s l tr , can be neglected. For a small ultrasound modulation along the path, which is a reasonable assumption in a low-pressure CW ultrasound field or when using pulsed and focused ultrasound, we can approximate exp(−F US,s (τ)/2) as 1 − . This allows us to express g 1 (τ, s) as g 1 (τ, s) = is the electric field autocorrelation function for the pathlength s, without ultrasound. The term g 1,US (τ, s) = − F US,s (τ)g 1,0 (τ,s) 2 is the perturbation of g 1,US (τ, s) due to ultrasound, obtained by scaling g 1,0 (τ, s) with the oscillating, ultrasound-induced term , which causes further decorrelation of g 1 (τ, s). We show that this approximation is valid for the low-pressure CW ultrasound that is explored in this work in Fig. 6 (appendix).
By adding the contributions from all pathlengths s, g 1 (τ) in the presence of ultrasound can also be expressed as a sum of two components: g 1 (τ) = g 1,0 (τ) + g 1,US (τ), where g 1,0 (τ) is given by Eq. (1) and it represents the electric field autocorrelation function due to blood flow-induced scatterer motion only. The g 1,US (τ) is perturbation due to ultrasound, given by Since F US,s (τ) = W US,s sin 2 ω u τ 2 , we can further simplify the expression for g 1 (τ) if we write This allows us to better group g 1 (τ) terms into oscillating and non-oscillating components as where for the small ultrasound modulation we assumed that g 1,0 (τ) g 1,US (τ) and g 1,0 (τ) − g 1,US (τ) ≈ g 1,0 (τ). Please note that because ultrasound modulation W US,s is larger along the longer pathlengths (Eq. (9)), the weight of longer pathlengths in expression for g 1,US (τ) is greater and g 1,US (τ) decorrelates faster with τ than g 1,0 (τ).
In order to obtain the expression for the intensity autocorrelation function g 2 (τ) in the presence of both blood flow and ultrasound, we apply the Siegert relation to the g 1 (τ) given by Eq. (9) and obtain where g 2,0 (τ) = 1 + β|g 1,0 (τ)| 2 is the intensity autocorrelation function due to blood flowinduced scatterer motion only. The ultrasound modulation of g 2 (τ) is given by the term M(τ) = M 0 (τ) cos(ω u τ), where modulation amplitude (sometimes also referred to as modulation Please note that in order to obtain Eqs. (10) and (11), we neglected small terms that are both non-oscillating and oscillating at two times the ultrasound frequency. A full treatment of the modulation can be found in the appendix, though fitting with only the dominant term is sufficient to accurately resolve blood flow. It is also important to note that modulation amplitude M 0 (τ) can be obtained experimentally, and it carries the information about the blood perfusion within the space region localized by the ultrasound. However, M 0 (τ) decays with τ faster than g 2,0 (τ) because g 1,US (τ) decorrelates faster with τ than g 1,0 (τ), which may need to be taken into account when fitting M 0 (τ) for BF i . In the theoretical derivation, we assumed that all optical scatterers along the photon path are exhibiting motion due to both the ultrasound and blood flow, which is the case when the blood perfusion and the plane ultrasound wave are present in the entire media. In typical DCS experiments, only a fraction of scattering events is from the RBCs and the ultrasound is confined to a fraction of the tissue volume. This means that only some portion of the photon paths may be associated with the ultrasound and/or Brownian motion-induced optical phase increments. Fortunately, dealing with such complexity of interactions may be relatively straightforward by using Monte Carlo simulations.

Monte Carlo simulation
Monte Carlo simulations were performed using the open source photon simulation software MCX. [31] A semi-infinite, homogenous media with optical scattering properties similar to what might be measured in vivo, µ s = 3.0 − 9.0 cm −1 , µ a = 0.05 cm −1 , n = 1.34, was simulated. Each simulation consisted of 10 8 launched photons. Ultrasound pressure attenuation of the media was given as 0.8 dB cm −1 MHz −1 for all simulated optical scattering properties [32]. Photon packets were collected at a source-detector separation of 1.8 cm, and the positions of scattering events along the paths of the detected photons were saved. Ultrasound induced phase modulation along each photon path was computed using Eqs. (3) and (15).
In Monte Carlo simulations, we can simulate the optical phase modulation along each path, and build up the electric field autocorrelation function as a sum of the individual autocorrelation functions along different paths as where w i is the photon weight for a path i through the tissue and g 1,i (τ) is the electric field autocorrelation function for the same photon path. Similar to expression in Eq. (5), and averaging is performed over time t. F i (τ) can be further expressed as a sum of two components: one due to the Brownian motion, F B,i (τ), and the other one due to the ultrasound modulation, . For a small ultrasound modulation along the path, g 1,0 (τ) and g 1,US (τ) can be simulated as g 1,0 , similar to expressions for the same parameters in Eqs. (8) and (9). In order to better understand the effect of ultrasound-modulation on g 1,US (τ), before summation over all simulated photon paths i, we regrouped the data based on the photon pathlength s to estimate P(s), W US,s , and F B,s (τ). This enabled us to obtain the values of g 1,US (τ), g 1,US (τ), and g 1,0 (τ) based on Eqs. (1) and (8), and to analyze how the difference between the decorrelation rates of g 1,0 (τ) and g 1,US (τ) is affected by the W US,s dependence on the pathlength. We subsequently computed the values of g 2,0 (τ), M 0 (τ), and g 2 (τ) based on Eqs. (10) and (11).
To compare the relative effects of the two mechanisms of ultrasound modulation (e.g., due to motion of the scatterers and due to changes in the index of refraction), autocorrelations were simulated for the following cases: i) no ultrasound modulation, ii) individual mechanisms of ultrasound modulation, and iii) combination of both mechanisms. The spatial distribution of ultrasound was also varied in MC simulations. When comparing between just simulation results, a plane wave ultrasound geometry filling the half-space was used. When comparing MC simulations to experimental results, the insonified region was modified to resemble the distribution that would be generated using the real probe, as seen below in Fig. 1(a), where only the area beneath the piezoelectric ring would experience modulation.
The spatial sensitivity of standard DCS to changes in flow is characterized by the difference in the measured change in flow from g 2,0 (τ) for a known flow perturbation at a given location. Due to the differences in the weighting of the pathlengths as predicted by Eq. (7), it might be expected that AOM-DCS has a different spatial sensitivity profile than standard DCS. To compare the spatial sensitivities of g 2,0 (τ) and M 0 (τ), spatial sensitivity maps were generated for the plane directly beneath the source and detector. In the MC simulations, 1 mm 3 cubic voxels were given flow perturbations individually, and g 2,0 (τ) and M 0 (τ) curves were generated. To quantify spatial sensitivity, sensitivity for a particular voxel was defined as S(r) = ∆BF i /∆BF i,r , where ∆BF i is the measured change in BF i and ∆BF i,r is the simulated change in BF i in the voxel at position r. Spatial sensitivity is compared between the full planar maps for g 2,0 (τ) and M 0 (τ), as well as the sensitivity with respect to depth. In order to validate the BF i measurement based on tagged photons, we compared BF i values estimated from M 0 (τ) to the BF i values obtained by a standard DCS without ultrasound. In addition, to evaluate the significance of applying a proper theoretical model when fitting M 0 (τ) for BF i , we fit the modulation depth for BF i , in two different ways: i) by using a proper expression for M 0 (τ) (Eq. (11); please also see Eq. (24)), and ii) 'naively', as though M 0 (τ) was the standard g 2,0 (τ) intensity autocorrelation function without ultrasound.

Phantom experiments
Experimental comparisons to the Monte Carlo simulations were made with a custom-built AOM-DCS system. Though these experiments are usually done with liquid phantoms, here, to reduce the ultrasound induced drive flow, semi-solid phantoms were chosen. Semi-solid, gelatin phantoms were mixed with 20% intralipid to obtain reduced scattering coefficients ranging from 3.0 to 9.0 cm −1 . Distilled water was mixed with gelatin powder (Knox) and heated to 80°C for 1 hour with continuous stirring. The mixture was allowed to cool to room temperature, and 20% intralipid was added to reach the desired scattering properties. Following mixing, the phantoms were refrigerated for 12 hours before measurements were taken. Light from a 785 nm long coherence length laser [DL785-100-S, Crystalaser] was delivered to the phantoms using a 62.5 µm GRIN multimode fiber and collected with four 5 µm single mode fibers at a source detector separation of 1.8 cm. Light from the detection fibers was sent to four single photon counting detectors [SPCM-AQRH14, Excelitas], and arrival times of the photons were collected and transmitted for later analysis. The source fiber was placed in the center of a piezoelectric ring transducer [SMR28D9T1111, STEMINC] and the detector fibers were placed at the edge, as seen in Fig. 1(a). The piezoelectric transducer had an inner diameter of 9 mm, an outer diameter of 29 mm, and a fundamental frequency of 2.08 MHz. The piezoelectric transducer was connected to a power amplifier [325LA, ENI], which amplified a continuous wave, sinusoidal waveform at 2 MHz from a function generator [SDG 5162,Siglent]. For all in vitro and in vivo experiments, the pressure generated by the ultrasound was measured to be 40 kPa by a hydrophone [HGL0085, Onda] in a water bath. Autocorrelation functions were calculated from the collected photon arrival timestamps, and the modulation depth was computed as the envelope of the modulation of the autocorrelation function.

In vivo experiments
For in vivo measurements, institutional ethical approval was obtained and three human subjects were measured using the AOM-DCS system. To ensure the safety of the subjects, the laser power was limited to 28 mW (ANSI limit for a laser sport size > 3 mm 2 at 785 nm). The ultrasound pressure generated by the probe was verified to be well below the mechanical index (MI) threshold for safety (MI = 1.9). As measured by a hydrophone in a water bath, the pressure was 40 kPa, giving an MI value of 0.0283 at 2 MHz. The AOM-DCS probe was placed on the subject's forearm, and a blood pressure cuff was applied to the upper arm. To measure a range of BF i values, the cuff was inflated in a step wise manner while measurements were being taken. Pressure remained constant for 1 minute at each level, and the ultrasound was turned ON and OFF every 10 seconds, giving three periods of 10 seconds of ultrasound ON and OFF at each pressure value. The results were quantified in the same manner as given above, and BF i was fit from both the autocorrelation curve and the modulation depth.

Results
3.1. Comparing simulated g 2,0 (τ) and M 0 (τ) for different ultrasound modulation mechanisms The examples of the g 2 (τ) curves generated by the MC simulations using an ultrasound pressure of 25 kPa at a frequency of 2 MHz for a reduced scattering coefficient of 6 cm −1 , BF i of 1.48 *10 −9 cm 2 /s, and a source-detector separation of 1.8 cm are presented in Fig. 1(b). The extracted modulation depths from simulated g 2 (τ) curves for different combinations of Brownian motion and individual ultrasound modulation mechanisms are presented in Fig. 1(c). In addition, the min-max normalized g 2,0 (τ) and M 0 (τ) curves are presented in Fig. 1(d), showing the differences in their decay rates. These differences in decay rates are quantified by fitting g 2,0 (τ) and M 0 (τ) for BF i using the correlation diffusion equation based on Eqs. (1) and (2), i.e., ignoring the effects of ultrasound modulation. A comparison of the simulated results as a function of reduced scattering coefficient is presented in Fig. 2(a). The BF i estimated by naively fitting M 0 (τ) exceeds the BF i from the g 2,0 (τ) by more than 18% in all simulated conditions. As derived in Eqs. (8)- (11), the degree to which the ultrasound interacts with the photons of different pathlengths is driving the difference between the decay rates of M 0 (τ) and g 2,0 (τ). In addition, contribution of different ultrasound modulation mechanisms to M 0 (τ) decay rate is also different for each pathlength s. The pathlength distributions and the mean squared optical phase accumulation W US,s for individual and combined ultrasound modulation mechanisms are compared in Figs. 2(b) and 2(c), respectively. The increased rate of decorrelation of M 0 (τ) in comparison to g 2,0 (τ) as well as the difference between the decorrelation rates of M 0 (τ) when including the influence from the ultrasound induced displacement of scatterers or index of refraction modulation or both, can be simply explained by considering the pathlength distributions P(s) weighted by the ultrasound modulation term W US,s . In both mechanisms, generally, longer pathlengths have a greater phase modulation due to ultrasound (Eq. (7)), and so are weighted more strongly in the adjusted pathlength distribution P(s)W US,s . The increased contribution of the longer paths to g 1,US (τ) causes the faster M 0 (τ) decay with τ (Eqs. (8)-(11)). When fitting M 0 (τ) for BFi, if weighting the photon pathlength distribution in g 1,US (τ) by W US,s is ignored and assumed to be a pathlength distribution predicted only by light propagation (i.e., if g 1,US (τ) is considered equal to g 1,0 (τ)), the fitted BF i values will be higher, reflecting the faster than expected decay. This is important to note, as the naively fitted BF i from M 0 (τ) is ∼20% greater than expected, which is a large error. Though when ultrasound modulation contribution to g 1,US (τ) is properly accounted for, accurate BFi estimates can be extracted from the M 0 (τ) (Fig. 2(a); MC fit).
The increased ultrasound modulation due to the index of refraction changes for shorter pathlengths (Fig. 2(c)) can be attributed to photons that travel largely superficially, parallel with the tissue surface and perpendicular to the ultrasound propagation direction. Because these photon paths are confined to a single or a few spatial periods of ultrasound, the modulation of optical phase along those paths is large. This effect is not captured by theoretical derivations (Eq. (7)), which assume diffusive light propagation, and it illustrates the importance of using MC simulations to properly model complex experimental geometries. The relative amplitudes of the two ultrasound modulation mechanisms are also important to consider when fitting for the BF i . Theoretical derivation of the relative contributions of the two mechanisms were given in Sakadzic et al. [26] as a function of the ultrasound wavenumber multiplied by the reduced mean free path. The ratio of the peak modulation due to each of the modulation mechanisms for a range of simulated µ s shows a good agreement with the theoretical derivation [23] (Fig. 2(d)). From these results and previous derivations, it can be concluded that the modulation due to the index of refraction changes represents the dominant effect for this implementation of AOM-DCS with relatively low ultrasound frequencies (1)(2)(3)(4)(5) and for scattering coefficients representative of soft tissue at NIR wavelengths.

Comparing the spatial sensitivity of DCS and AOM-DCS derived from Monte Carlo simulation
As seen in the previous section, the use of ultrasound gives an increased sensitivity to photons with longer pathlengths. To quantify the effects of the changes in the weighting of the pathlengths on spatial sensitivity, maps were derived from the simulated standard autocorrelation function, g 2,0 (τ), and modulation, M 0 (τ), fitted with their respective proper theories, and can be seen in Fig. 3(a) and 3(b). Both maps are presented in log scale, and isolines are drawn at labeled values of the sensitivity for comparison. To quantify the depth sensitivity, the maps were averaged in the direction parallel to the source and detector, and the percent difference of the depth sensitivity was calculated, and can be seen in Fig. 3(c). From both presentations of the data, two regimes of differences can be identified, 1) the reduction in sensitivity of AOM-DCS at shallower depths relative to standard DCS and 2) the monotonic increase in AOM-DCS sensitivity relative to standard DCS as the depth of the change increases. This comparison connects the observed increase in sensitivity to photons with longer pathlengths to the possible benefits of utilizing ultrasound modulation for increased sensitivity to cerebral blood flow in vivo, even in the case of CW ultrasound modulation.

BF i measurements in the gelatin phantoms having different scattering coefficient and temperature
To validate the M 0 (τ) fitting for BF i in vitro, the same range of µ s values explored in MC simulations was used in the semi-solid gelatin phantoms. The measurements were performed using the experimental setup shown in Fig. 1(a). The examples of measured g 2 (τ), g 2,0 (τ), and For an MC simulation with µ a = 0.05 cm −1 , µ s ' = 6.00 cm −1 , and a source-detector separation of 1.8 cm, (a) the corresponding spatial sensitivity maps for standard DCS and AOM-DCS. Averaging the percent difference in the X direction between the two maps gives the results seen in (b), which indicates that AOM-DCS is less sensitive to changes in the more superficial layers and more sensitive to changes in the deeper layers than is standard DCS. The errorbars represent the standard deviation of the percent difference along the X direction.
M 0 (τ) are shown in Fig. 4(a). Signal to noise ratio (SNR) of g 2,0 (τ) and M 0 (τ) were verified to be sufficiently high in each measurement to ensure accurate fitting. An example of the SNR of the curves is presented in Fig. 4(d). The SNR of g 2,0 (τ) and M 0 (τ) were quantified as SNR g2,0 (τ) = (g 2,0 (τ) − 1)/σ g2,0 (τ) and SNR M0 (τ) = M 0 (τ)/σ M0 (τ), respectively, where σ(τ) is the standard deviation of the curve at a lag time τ. The SNR of M 0 (τ) is seen to be highly sensitive to the pressure used, even for small changes in the pressure magnitude. This is beneficial, as the SNR can be increased through slight increases in the pressure, while maintaining the assumption that ultrasound induced phase modulation is small. As predicted by both theoretical derivations and MC simulations, the faster decay rate of M 0 (τ) in comparison to g 2,0 (τ) was observed for each µ s (data not shown). In order to accurately extract the BF i , the fitting of the M 0 (τ) depth was performed using Eqs. (8)- (11), which properly accounts for the ultrasound tagging of light. The contributions of individual photon pathlengths were quantified by the Monte Carlo simulations, taking into account both the optical properties of the phantoms and the experimental geometry. A good agreement between BF i values based on g 2,0 (τ) and M 0 (τ) was observed for all µ s (Fig. 4(b)). Similarly accurate measurements of BF i was obtained when gelatin phantom temperature was changed from 4 to 16°C, which increased scatterer diffusivity while keeping the µ s constant at 6 cm −1 (Fig. 4(c)). Although comparison to the Stokes-Einstein equation would have been ideal over the temperature range measured, this comparison is complicated by the use of polydisperse scatterer size, and the gel-like nature of the phantom, so a linear trend line is used to estimate the temperature dependence. Note that the variation in the actual BF i of the phantoms at different scattering values (Fig. 4(b)) is simply due to variations in the phantom fabrication procedure and does not relate to the change in scattering.

In vivo demonstration of AOM-DCS
To validate the proposed AOM-DCS theory and experimental setup in vivo, the BF i was measured in the forearm of 3 human subjects over a range of perfusion values created by inflating a blood pressure cuff with different pressures. The BF i estimated from the g 2,0 (τ) was compared to the BF i estimated from the M 0 (τ) using both the naïve and the proper fitting, shown in Figs. 5(a) and 5(b), respectively. These results are demonstrating that BF i can be measured accurately based on the ultrasound-tagged light. In addition, the use of the proper theoretical model results in a more accurate BF i measurements. Though the difference between naïve and proper fitting may be slight in this case, in a more optically heterogeneous area in the body, such as the head, with a more complex spatio-temporal ultrasound distribution and higher blood flow, the differences in the fits could become more significant, and the use of proper theoretical method may allow for more accurate quantification of BF i .

Discussion
In the MC simulations, by separating the ultrasound modulation mechanisms, the influence of each mechanism on the modulation depth was assessed. We confirmed previous observations that the index of refraction modulation is a dominant mechanism under simulated conditions, selected to match conditions that would be seen in vivo. Dependence of the modulation depth was also assessed as a function of optical scattering to determine the influence of different scattering conditions, with dependencies being found in the relative intensities of the two mechanisms, as was predicted by previous theory [26], as well as the rate at which each method drove decorrelation. While the index of refraction modulated signal was found to come to a steady state of decorrelation rate as a function of scattering, at low scattering, the displacement modulation was found to decorrelate faster, indicating that a certain number of scattering events are most likely required to reach a steady state of measured decorrelation rate. Fortunately, in the lower scattering case, the relative weighting of the displacement modulation component is low compared to the index of refraction modulation, and the changes to the overall decorrelation rate are minimal. This finding is important in understanding the relative contributions of the different mechanisms and their effects of the measured modulation. The observation that the modulation depth decorrelates faster than g 2,0 (τ) is also an important one, indicating that accurate quantification of perfusion from ultrasound modulated DCS requires considering the difference in the mechanisms driving modulation depth decay vs. intensity auto-correlation decay in the absence of ultrasound. The derivations presented provide a framework with which to analyze the modulated signals, and the concordant BF i measurement results of standard DCS and AOM-DCS techniques provide evidence that, for the particular conditions these measurements were taken under, namely a small degree of ultrasound modulation of the optical phase, the assumptions presented to allow for separation of the modulated and unmodulated components of g 2 (τ) are valid. Analysis of the relative weighting of pathlengths provides a reasonable explanation, that longer pathlengths are tagged preferentially relative to shorter ones, which was confirmed in simulation and phantom experiments. The relatively increased weighting of longer pathlengths could potentially provide a mechanism to probe dynamics deeper in tissue, which can be helpful for sensing the blood flow in the brain, even in the case of CW ultrasound tagging.

Conclusion
Here we have developed an analytical model to use the ultrasound modulated light for estimation of BF i in DCS. The model was validated using Monte Carlo simulation for different values of the reduced scattering coefficient, and experimentally in solid gelatin phantoms and three human subjects using a novel AOM-DCS system. This work provides a basis by which to move forward in the development of a spatially resolved, acousto-optic modulated diffuse correlation spectroscopy system to measure blood flow in tissue.

A. Appendix
Ultrasound-induced changes of the optical phase along the photon paths Here we detail the full derivation of the ultrasound induced, optical phase modulation. The displacement of a scattering particle, ∆r s (t) = r s (t) − r s , in the media due to ultrasound, where r s is the unperturbed location of the scatterer, can be expressed as where S u and ϕ u are the amplitude and phase deviations of the particle from the media moving around it,Ω u is the unit vector in the ultrasound propagation direction, ρ is the density of the media, and ν u is the ultrasound speed in the media [26,27]. The optical phase increment that is due to the scatterer i motion induced by the ultrasound is given by the dot product of the momentum transfer, q i = k i − k i+1 , with the scatterer displacement, ∆r i (t), induced by the ultrasound. The phase increment generated along the path due to ultrasound-induced displacements of the scatterers is then given by the sum of the phase increments over all scattering events along the path, given by Eq.
Where η is the elasto-optic coefficient of the media, equal to η = ρ ∂n ∂ρ , which for water is approximately equal to 0.32 [26]. The ultrasound-induced optical phase increment along the free path between two scattering particles r i and r i+1 is given as ∆φ i (t) = ∫ r i+1 r i k 0 ∆n(r, t)dr, where ∆n(r, t) = n 0 η ρν 2 u P(r, t). For an entire photon pathlength through tissue, the phase modulation due to index of refraction modulation can be given as an expanded form of Eq. (4), Combining and expanding the expressions given in Eqs. (3) and (15) for the individual mechanisms of ultrasound modulation gives the entire phase modulation along the path by ultrasound, referred to as ∆φ US (t), and is given by cos(φ P (t, r))dr where φ P (t, r), the phase term for the ultrasound pressure, is given by φ P (t, r) = ω u t − k u · r + φ.
For the case of CW ultrasound insonification in the full half space, the solution to the integral for each free path is solved and the terms deriving from the motion of the particles are expanded to properly include geometry of the path. If we assume that amplitude and phase of scatterers oscillation due to ultrasound are equal to those of the surrounding media, the ∆φ US (t) can be expressed as: Further expanding the momentum term gives the change in phase accumulation just in terms of scatterer position, and simplifies the phase calculation for the Monte Carlo (MC) simulation: Combining Eq. (18) with the term describing phase accumulation due to Brownian motion gives the total change in optical phase along a path as a function of a lag time τ, ∆φ This expression can be used for each of the optical paths to generate the change in phase at any lag time, τ. In the main text, we assume that ∆φ T (t, τ) is much less than unity for the regime of dynamic scattering that we are measuring, allowing for the approximation given in Eq. (7) to be valid. Here we show a comparison of the electric field temporal autocorrelation functions calculated with the approximation and calculated directly as ) for photon paths selected from the MC simulation, shown below in Fig. 6 to demonstrate the similarity of the autocorrelation curves and to validate the use of the approximation. Fig. 6. Electric field temporal autocorrelation function computed directly from the phase of the electric field (solid lines) for three photon pathlengths compared to the electric field temporal autocorrelation function computed as given in Eq. (5) in the main text (circles), with the inset showing the initial portion of the g 1 (τ) decay. The correspondence between the two over the large majority of g 1 (τ) indicates that this approximation should allow for accurate characterization of the ultrasound induced phase using this model.

Deriving the expression for the modulation depth (M(τ)) from the summation of pathlengths
An abbreviated derivation of the modulation depth from the intensity autocorrelation function (g 2 (τ)) is given in the main text, and will be elaborated on here. As given in Eq. (1), the electric field autocorrelation function (g 1 (τ)) can be described as an integral over the pathlength distribution of autocorrelations for individual pathlengths, recapitulated here in Eq. (20).
Applying the approximation for ∆φ T (t) much less than unity, we can express Eq. (20) in terms of the mean squared phase accumulation from the Brownian motion and ultrasound modulation, as is expressed in Eq. (8). For small ultrasound phase modulation, we can approximate the electric field autocorrelation as This integral can be seen to be a sum of the unperturbed g 1 (τ), referred to in the main text as g 1,0 (τ), and a tagged g 1 (τ), which is weighted by the influence of the ultrasound for a pathlength s. Separating the integrals and expressing F US,s as W US,s sin 2 ω u τ 2 followed by the separation of oscillating and non-oscillating terms gives The second term in Eq. (23) is what is called g 1,US (τ) in the main text. To express the experimentally measured quantity, g 2 (τ), we apply the Siegert relation and substitute in both g 1,0 (τ) and g 1,US (τ) to give g 2 (τ) = 1 + β(g 1,0 (τ) − g 1,US (τ)) 2 + 2β cos(ω u τ)g 1,US (τ)(g 1,0 (τ) − g 1,US (τ)) +βcos 2 (ωτ)(g 1,US (τ)) 2 (24) As was done in the main text, we assume that g 1,0 (τ) g 1,US (τ), so g 1,0 (τ) − g 1,US (τ) ≈ g 1,0 (τ) and the component oscillating at two times the ultrasound frequency is found to be small compared to the component oscillating at the ultrasound frequency. Simplifying Eq. (24) with these two steps gives Eqs. (10) and (11), which are used to fit the ultrasonically tagged autocorrelation function.
Computing the modulation of phase in Monte Carlo and using the results to fit the measured modulation of the autocorrelation function (M 0 (τ)) In the MC simulations, ultrasound phase accumulations at each scattering center and along the free paths between scattering centers are computed based on Eqs. (3) and (15). Contributions from each mechanism were computed, and mean squared phase increment differences were calculated for each photon path, i, and given as ∆φ 2 US,d,i (t, τ) , ∆φ 2 US,n,i (t, τ) , and ∆φ 2 US,d+n,i (t, τ) .
To compute the results presented in Figs. 2(b) and 2(c), the mean squared phase increment differences were binned along the photon pathlengths s, into bins such that where ∆φ 2 US (s k ≤ s<s k+1 ) is the averaged mean squared phase modulation by ultrasound between two pathlength bin edges s k and s k+1 . Multiplying each mean squared phase term, ∆φ 2 US (s) , by the pathlength distribution, P(s), results in the modified pathlength distributions given in Fig. 2(b). The modulation of the autocorrelation function was fit using these weighted pathlength distributions with the expression given in Eq. (11) in the main text. Estimates of the ultrasound tagged electric field autocorrelation, g 1,US (τ), were calculated using an expression similar to Eq. (12) in the main text, where the weighting term is given by w(s) = P(s) ∆φ 2 US (s), and the individual path electric field autocorrelation function is given as g 1,s (τ) = exp − 1 3 µ s k 2 ∆r 2 (τ)s . The estimate of the unperturbed, electric field autocorrelation, g 1,0 (τ), was calculated using Eq. (1) in the main text, and the product of the two was used to generate the estimate of the modulation depth, M 0 (τ) = A * g 1,US (τ) * g 1,0 (τ), for a given Brownian, mean squared displacement, ∆r 2 (τ) = 6D b τ, and modulation amplitude, A. The estimated modulation, M 0 (τ), was compared to the measured modulation, and the sum squared error, (M 0 (τ) − M 0 (τ)) 2 , was minimized by changing the estimated diffusion coefficient, D b , and the modulation amplitude, A, to quantify the flow measured from the modulation decay.