Beyond diffuse correlations: deciphering random flow in time-of-flight resolved light dynamics

V. N. Du Le; Vivek J. Srinivasan; Vivek J. Srinivasan

doi:10.1364/OE.385202

1. Introduction

Diffusing Wave Spectroscopy (DWS) was formulated in the late 1980s to probe dynamics of scatterers in soft matter based on intensity fluctuations of multiply scattered light. DWS models the field autocorrelation as an integral over photon paths, assuming uncorrelated motion at different scattering events [1–3], in the light diffusion approximation. In biological tissues, DWS and its differential formulation of diffuse correlation spectroscopy (DCS) can measure red blood cell (RBC) dynamics, which correlate empirically with blood flow [4–6]. The assumptions underlying DWS theory are powerful and widely used, even in distant fields such as laser speckle flowmetry [7,8]. However, given recent experimental advances in measuring time-of-flight- (TOF-) resolved light dynamics [9–11] in detail for the first time, some assumptions in DWS/DCS of blood flow merit further scrutiny. Here we focus on two important and inter-related assumptions; 1) the nature of RBC motion, and 2) the first cumulant approximation of the average over dynamic scattering angles.

Particle motion can be ordered, as in advective, randomly oriented, flow, or disordered, as in Brownian motion, or a combination of both [7]. DWS/DCS measurements of particle motion in biological tissues such as the rat brain [12,13], piglet brain [14], human brain [6,15–18], breast [19,20], head and neck [21], and skeletal muscle [22] have shown that experimental data is best fitted by ignoring advection of RBCs, and only considering Brownian motion. Thus, the blood flow index (BFI) in DWS/DCS has units of a Brownian diffusion coefficient.

The starting point of dynamic light scattering theories [23,24] is the normalized field autocorrelation of m^th order paths: $\langle \textrm{exp}\left[ {i\mathop \sum \nolimits_{j = 1}^m {{\boldsymbol q}_j} \cdot \Delta {{\boldsymbol r}_j}} \right] \rangle $, where 〈〉 denotes an average over particle displacement vectors $\Delta {{\boldsymbol r}_j}$ and dynamic momentum transfer vectors ${{\boldsymbol q}_j}$, which are assumed to be independent, and $i = \sqrt { - 1} $. In the multiple scattering regime, assuming uncorrelated dynamic scattering events while neglecting both random flow and higher order terms in the cumulant expansion [23], DWS expresses the normalized TOF-resolved field autocorrelation function, g₁, as:

(1)$$g_1^{DWS}({\tau _s},{\tau _d}) = \exp \left[ { - \overline m ({\tau_s}){{\left\langle {{q^2}(\theta )} \right\rangle }_\theta }{D_B}{\tau_d}} \right]$$

In Eq. (1), $\bar{m}$ is the average number of dynamic scattering events at time-of-flight (TOF), τ_s, 〈q²(θ)〉_θ is the angle-averaged momentum transfer for dynamic scattering, D_B is the Brownian diffusion coefficient of dynamic particles, and τ_d is autocorrelation time lag.

It is notable that most studies of blood flow in biological tissues utilize Eq. (1), although this expression was first derived for uniform colloidal suspensions, foams, and gels with unordered (Brownian) motion. A major biological function of RBCs is the transport of gases in the circulation; thus while hydrodynamic diffusion of RBCs is well-documented [25,26], RBCs should also possess characteristics of ballistic motion [24]. Therefore, it seems incongruous that blood flow, in this conventional sense, is typically assumed not be observable through light scattering techniques.

In this work, we present a plausible explanation to help resolve these discrepancies. We simulate field autocorrelations with a custom extended Monte Carlo method that explicitly allows for a mixture of static and dynamic particles with different scattering phase functions. We attempt to accurately model both the high scattering anisotropy and low volume fraction of RBCs in tissue. We assess the results to critically examine the first cumulant approximation, or first cumulant expansion, which is common in DWS. We show that under certain circumstances, errors due to the cumulant approximation, on one hand, and neglecting random flow, on the other, can cancel, enabling the conventional theory in Eq. (1) to describe autocorrelations well, even while ignoring advection. We then propose a more accurate approach that avoids these assumptions and might enable direct quantification of RBC advection in dynamic light scattering measurements.

2. Field autocorrelation models

In this section, we develop candidate models for TOF-resolved field autocorrelations, highlighting key assumptions in each.

2.1 Nature of scatterer motion

In DWS/DCS, the mean-squared displacement at time lag τ_d, 〈Δr²(τ_d)〉, of RBCs is modeled as either random flow, given by 〈Δr²(τ_d)〉 = v²τ_d², where v² is the second moment of the Gaussian velocity distribution, or as Brownian motion, given by 〈Δr²(τ_d)〉 = 6D_Bτ_d, where D_B is the Brownian diffusion coefficient. In this work, we allow for RBCs to exhibit a “hybrid” of both random flow and Brownian motion, given by 〈Δr²(τ_d)〉 = 6D_Bτ_d + v²τ_d², where the two types of motion are assumed to be independent.

2.2 Adapting Bonner and Nossal’s theory

To assess the validity of the cumulant approximation (or first cumulant expansion), we start from Bonner and Nossal’s theory [24], which implicitly includes all cumulant terms by directly integrating over the dynamic scattering phase function. This theory envisions tissue as a mixture of static (tissue matrix) and dynamic scatterers (i.e. RBCs and other scattering blood components). We extend this theory as follows: 1) the lower limit of number of dynamic scattering events, m, is set to zero to include purely static scattering paths, 2) the Henyey-Greenstein (and later, the more general Gegenbauer kernel) scattering phase function for dynamic scattering is assumed in deriving the normalized autocorrelation for a single dynamic scattering event (g₁^ss), 3) the average number of dynamic scattering events, $\bar{m}$(τ_s), is assumed to be $\bar{m}$(τ_s) = p_dynµ_sτ_sc/n, where p_dyn is the probability of dynamic scattering, µ_s is the static scattering coefficient, c is speed of light in vacuum, and n is the index of refraction of medium, and 4) as described above, we explicitly allow for RBC Brownian motion, in addition to random flow considered Bonner and Nossal’s original theory.

It is appropriate to remark on p_dyn, the probability that a scattering event is dynamic. Note that p_dyn is given by p_dyn = V_Bµ_s,blood/µ_s, where V_B is blood volume fraction and µ_s,blood is blood scattering coefficient, and does not depend on anisotropy factor. Given µ_s = 100 cm⁻¹, at a hematocrit of 33%, µ_s,blood/µ_s ≈ 7.8 at 850 nm [27], and p_dyn = 10% corresponds to V_B = 1.3%. On the other hand, the DCS probability of dynamic scattering, α [14,17], is given by α = V_Bµ_s’_,blood/(V_Bµ_s’,_blood+µ_s’) = p_dyn (1-g_dyn)/[p_dyn(1-g_dyn)+(1-g_stat)], where g_dyn and g_stat are anisotropy factors of dynamic and static scatterers, respectively. If p_dyn(1-g_dyn) << (1-g_stat), α ≅ p_dyn(1-g_dyn)/(1-g_stat). For example, α ≅ 2% when p_dyn = 10%, g_dyn = 0.98 (typical for RBCs [27]) and g_stat =0.9.

Bonner and Nossal [24] assumed that multiple static scattering in tissue is sufficient to achieve random illumination of dynamic particles and uncorrelated dynamic scattering events. This enabled them to construct the field autocorrelation in the multiple dynamic scattering regime from the single dynamic scattering field autocorrelation, g₁^ss. Our extension of their expression is:

(2)$$g_1^{BN}({\tau _s},{\tau _d}) = {\sum\limits_{m = 0}^\infty {{P_m}({\tau _s})[{g_1^{ss}({\tau_d})} ]} ^m}.$$

In Eq. (2), P_m is the probability that a photon will experience m scattering events with dynamic particles before leaving the medium, g₁^ss is the autocorrelation for a single scattering event from a moving particle with random illumination, as a function of time lag, τ_d:

(3)$$\begin{aligned} g_1^{ss}({\tau _d}) &= {\left\langle {\exp \left[ { - \frac{{{q^2}(\theta )}}{6}\left\langle {\triangle {r^2}({\tau_d})} \right\rangle } \right]} \right\rangle _\theta }\\ & \textrm{ } = \frac{{\int\limits_0^\pi {{p_{HG}}(\theta )\exp \left[ { - \frac{{{q^2}(\theta )}}{6}\left\langle {\triangle {r^2}({\tau_d})} \right\rangle } \right]\sin (\theta )d\theta } }}{{\int\limits_0^\pi {{p_{HG}}(\theta )\sin (\theta )d\theta } }}. \end{aligned}$$

In Eq. (3), the angle brackets 〈〉_θ denote an angular weighted average over dynamic scattering angle θ, and p_HG(θ) is the Henyey-Greenstein scattering phase function for dynamic scatterers with an anisotropy of g_dyn [28–31]:

(4)$${p_{HG}}(\theta ) = {{({1 - g_{dyn}^2} )} \mathord{\left/ {\vphantom {{({1 - g_{dyn}^2} )} {[{4\pi {{({1 + g_{dyn}^2 - 2{g_{dyn}}\cos \theta } )}^{3/2}}} ]}}} \right.} {[{4\pi {{({1 + g_{dyn}^2 - 2{g_{dyn}}\cos \theta } )}^{3/2}}} ]}},$$

and

(5)$$\int\limits_0^\pi {{p_{HG}}(\theta )\sin (\theta )d\theta } = {\raise0.7ex\hbox{$1$} \!\mathord{\left/ {\vphantom {1 {2\pi }}} \right.}\!\lower0.7ex\hbox{${2\pi }$}}$$

Figure 1(a) compares the Henyey-Greenstein (HG) phase function for g_dyn = 0.6 (typical for Intralipid [32]) to g_dyn = 0.98 (typical for RBCs [27]), with polar plots in the inset. Forward scattering dominates for g_dyn = 0.98. Assuming that the number of dynamic scattering events P_m has a Poisson distribution [24], as detailed in Appendix A, Eq. (2) becomes:

(6)$$g_1^{BN}({\tau _s},{\tau _d}) = \exp [{\overline m ({\tau_s})\{{g_1^{ss}({\tau_d}) - 1} \}} ].$$

Fig. 1. The cumulant approximation error [Eq. (10)] increases with anisotropy of dynamic scattering, g_dyn. (a) Exemplary Henyey-Greenstein phase functions with normalized polar plots in inset. (b) Cumulant approximation error (R) versus g_dyn, for different $\bar{m}$ values, p_dyn = 10%, and μ_s = 100 cm⁻¹.

Download Full Size | PDF

2.3 Critical examination of the cumulant approximation

The cumulant approximation conveniently brings the angular average of Eq. (3) into the exponent, providing a very simple and intuitive expression for the autocorrelation function. In this section, we revisit the cumulant expansion, paying special attention to errors incurred by the high scattering anisotropy of RBCs. Commonly, in the derivation of DWS [23], the normalized field autocorrelation is succinctly expressed as:

(7)$$g_1^{}({\tau _s},{\tau _d}) = \left\langle {\exp \left[ {({ - 1/6} ){q^2}(\theta )\left\langle {\triangle {r^2}({\tau_d})} \right\rangle } \right]} \right\rangle _\theta ^{\overline m ({\tau _s})} = {[{g_1^{ss}({\tau_d})} ]^{\overline m ({\tau _s})}}$$

In contrast to the Bonner and Nossal summation [Eq. (2)], which posits an integer distribution of dynamic scattering events, Eq. (7) assumes that the number of dynamic scattering events is equal to the (possibly non-integral) average value. Note that while both Eqs. (6) and (7) yield the same first cumulant term, higher order cumulant terms are different. Here, we examine the cumulant expansion of the extended Bonner and Nossal expression in Eq. (6), which incorporates fewer assumptions. In this case, Taylor expansion can be applied to the natural logarithm of Eq. (6), as detailed in Appendix B, to yield the first and second cumulant terms:

(8)$$\begin{array}{l} g_1^{BN}({\tau _s},{\tau _d}) = \exp \left\{ { - \overline m ({\tau_s})\textrm{ }\left[ {\frac{{{{\left\langle {{q^2}(\theta )} \right\rangle }_\theta }}}{6}\left\langle {\triangle {r^2}({\tau_d})} \right\rangle - \frac{1}{{2!}}\frac{{{{\left\langle {{q^4}(\theta )} \right\rangle }_\theta }}}{{{6^2}}}{{\left\langle {\triangle {r^2}({\tau_d})} \right\rangle }^2} - O\left( {{{\left\langle {\triangle {r^2}({\tau_d})} \right\rangle }^3}} \right)} \right]} \right\}\\ \textrm{ } \end{array}$$

Note that third and higher cumulant terms are contained in O(〈Δr²(τ_d)〉³), and the subscript θ in angular average has been dropped for readability. We examine the second cumulant to determine the error in the first cumulant approximation. Then, Eq. (8) can be simplified as:

(9)$$\begin{array}{c} g_1^{BN}({\tau _s},{\tau _d}) = \exp \left[ { - \textrm{ }\frac{1}{6}\overline m ({\tau_s})\left\langle {{q^2}(\theta )} \right\rangle \left\langle {\triangle {r^2}({\tau_d})} \right\rangle \left( {1 - \frac{1}{6}\overline m ({\tau_s})R\left\langle {{q^2}(\theta )} \right\rangle \left\langle {\triangle {r^2}({\tau_d})} \right\rangle } \right)} \right]\textrm{ }\\ \textrm{ } \times \textrm{ }\exp \left[ {\overline m ({\tau_s})\textrm{ }O\left( {{{\left\langle {\triangle {r^2}({\tau_d})} \right\rangle }^3}} \right)} \right], \end{array}$$

where R is the magnitude ratio of the second cumulant to first cumulant-squared, given by:

(10)$$R = {{\left\langle {{q^4}(\theta )} \right\rangle } \mathord{\left/ {\vphantom {{\left\langle {{q^4}(\theta )} \right\rangle } {\left[ {2\overline m {{\left\langle {{q^2}(\theta )} \right\rangle }^2}} \right]}}} \right.} {\left[ {2\overline m {{\left\langle {{q^2}(\theta )} \right\rangle }^2}} \right]}},$$

where

(11)$$\left\langle {{q^2}(\theta )} \right\rangle = \left\langle {{2^2}{k^2}{{\sin }^2}(\theta /2)} \right\rangle = 2{k^2}({1 - {g_{dyn}}} )$$

and

(12)$$\left\langle {{q^4}(\theta )} \right\rangle = \left\langle {{2^4}{k^4}{{\sin }^4}(\theta /2)} \right\rangle .$$

As shown in Fig. 1(b), R slowly increases with g_dyn when g_dyn < 0.9, and rapidly increases with g_dyn when g_dyn > 0.9. For example, when $\bar{m}$ = 214.3 (τ_s = 1000 ps and μ_s = 100 cm⁻¹), R increases 3.2 times from g_dyn = 0.6 to g_dyn = 0.9 and increases 4.8 times from g_dyn = 0.9 to g_dyn = 0.98, or overall 15.4 times from g_dyn = 0.6 to g_dyn = 0.98. Furthermore, the second cumulant term vanishes as R approaches zero [Eq. (9)], after many dynamic scattering events. This happens more rapidly for smaller g_dyn (Fig. 1(b)). Thus, the error in the first cumulant approximation merits particular attention for large g_dyn.

Ultimately, the functional form, or “shape” of the autocorrelation decay with respect to time lag, τ_d, depends on which polynomial orders of τ_d contribute in the exponent. With hybrid motion, different cumulant orders can contribute to the same polynomial order of τ_d. Thus, another way to simplify Eq. (8) is to group terms according to the order of τ_d:

(13)$$\begin{array}{ll} g_1^{BN}({\tau _s},{\tau _d}) &= \exp \left[ { - \frac{1}{6}\overline m ({\tau_s})\left\langle {{q^2}(\theta )} \right\rangle 6{D_B}{\tau_d}} \right]\\ &\textrm{ } \times \textrm{ }\exp \left[ { - \overline m ({\tau_s})\tau_d^2\left( {\frac{{\left\langle {{q^2}(\theta )} \right\rangle }}{6}{v^2} - \frac{{\left\langle {{q^4}(\theta )} \right\rangle }}{{72}}{{({6{D_B}} )}^2}} \right) + O({\tau_d^3} )} \right] \end{array}$$

The term O(τ_d³) in Eq. (13) includes contributions from higher order τ_d terms. Although written in different ways, Eqs. (9) and (13) both represent expressions for the cumulant expansion. Importantly, Eq. (13) clearly shows that both Brownian motion and random flow contribute to the τ_d² term, representing the lowest order deviation from the pure exponential decay of DWS.

2.3.1 Special case 1: Brownian motion dominates random flow (A1)

As Brownian motion is unordered on time scales larger than the collision time, the mean squared displacement increases as τ_d, whereas for ordered motion, the mean squared displacement increases as τ_d². When τ_d << 6D_B/v², Brownian diffusion dominates random flow (A1), i.e. 〈Δr²(τ_d)〉 ≈ 6D_Bτ_d. This approximation is implicit in Eq. (1). Many studies have empirically found that the Brownian motion model leads to a better fit to experimental data than the random flow model [12–22], but did not explicitly verify that τ_d << 6D_B/v². For example, for v = 1.25 mm/s and D_B = 6 × 10⁻¹¹ m².s⁻¹, this condition corresponds to τ_d << 220 μs. Thus, given A1, we modify Bonner and Nossal’s theory, which implicitly includes all cumulants [Eq. (6)], to describe Brownian diffusion.

2.3.2 Special case 2: first cumulant dominates second cumulant (A2)

If the first cumulant term dominates the second cumulant term, Eqs. (9) and (13) reduce to Eq. (14). This occurs when the mean-squared displacement of the moving particles satisfies 〈Δr²(τ_d)〉 << 6 / [$\bar{m}$ R〈q²(θ)〉]. For example, with g_dyn = 0.98 (i.e. RBCs), v = 1.25 mm/s and D_B = 6 × 10⁻¹¹ m².s⁻¹, the condition is satisfied when τ_d << 141 μs. Thus, given A2, we modify conventional DWS theory, which neglects higher order cumulants, to describe hybrid motion.

(14)$$g_1^{BN}({\tau _s},{\tau _d}) = \exp \left[ { - \frac{1}{6}\overline m ({\tau_s})\left\langle {{q^2}(\theta )} \right\rangle ({{v^2}\tau_d^2 + 6{D_B}{\tau_d}} )} \right]$$

Note that the lag time range for neglecting higher order cumulants (τ_d << 141 μs) coincides with the lag time range for neglecting random flow (τ_d << 220 μs), thus this special case has little practical significance, given our assumed parameters. Nevertheless, it is included for completeness to isolate each individual assumption.

2.3.3 Special case 3: first order τ_d term dominates (A1&2)

While special case 1 compares magnitudes of ordered and unordered motion, in special case 3, we examine τ_d terms which determine the ultimate shape of the autocorrelation. From Eq. (13), the general condition for the first order τ_d term to dominate the τ_d² term is given by:

(15)$${\tau _d} < < {{{D_B}} \mathord{\left/ {\vphantom {{{D_B}} {\left|{\frac{{{v^2}}}{6} - \overline m ({\tau_s})R\left\langle {{q^2}(\theta )} \right\rangle D_B^2} \right|}}} \right.} {\left|{\frac{{{v^2}}}{6} - \overline m ({\tau_s})R\left\langle {{q^2}(\theta )} \right\rangle D_B^2} \right|}}$$

Note that Eq. (15) can be more or less stringent than the condition in special case 1, depending on the relative magnitudes of ordered and unordered motion in the denominator. When Eq. (15) holds, Eq. (13) can be simplified, up to O(τ_d³), to

(16)$$g_1^{BN}({\tau _s},{\tau _d}) \approx \exp \left[ { - \overline m ({\tau_s})\left\langle {{q^2}(\theta )} \right\rangle {D_B}{\tau_d}} \right] = g_1^{DWS}({\tau _s},{\tau _d})$$

Importantly, Eq. (16) is equivalent to conventional DWS [Eq. (1)], which is derived by neglecting both higher order cumulant terms and random flow (A1&2). Table 1 summarizes all three special cases.

Table 1. Summary of special cases in the second cumulant approximation

View Table | View all tables in this article

2.3.4 First cumulant for random flow cancels second cumulant for Brownian diffusion

From Eq. (13), the coefficient of τ_d² includes contributions from the first cumulant term for random flow and second cumulant term for Brownian motion. The velocity distribution may be selected such that these two contributions cancel, minimizing deviations from an exponential decay at small τ_d. The distribution standard deviation, v*, is then given by:

(17)$$v\ast{=} {D_B}\sqrt {{{3\left\langle {{q^4}(\theta )} \right\rangle } \mathord{\left/ {\vphantom {{3\left\langle {{q^4}(\theta )} \right\rangle } {\left\langle {{q^2}(\theta )} \right\rangle }}} \right.} {\left\langle {{q^2}(\theta )} \right\rangle }}} = {D_B}\sqrt {6\overline m ({\tau _s})R\left\langle {{q^2}(\theta )} \right\rangle }$$

Note that this value of v* ensures that Eq. (16) is valid across lags. For an HG phase function with g_dyn =0.98 and D_B = 6 × 10⁻¹¹ m².s⁻¹ [26,33,34], v* has a value of 1.25 mm/s. This is a reasonable value of RBC speed [35,36]. Thus, while our velocity distribution is chosen to eliminate the second cumulant, it remains well within physiological norms.

2.4 Monte Carlo simulations

Because experimental systems with Brownian motion and random flow are challenging to contrive, in this study, Monte Carlo is taken as the gold standard. A previous Monte Carlo algorithm simulating the time resolved reflectance of a semi-infinite medium [37] was modified to calculate the normalized field autocorrelation. In the modified version, both dynamic scattering (moving particles, with a customized anisotropy g_dyn) and static scattering (background tissue scattering, with a fixed anisotropy, g_stat = 0.9) are explicitly simulated, with a probability of dynamic scattering p_dyn. Only dynamic scattering events contribute to the recorded total accumulated momentum transfer Y_dyn = Σ_j(1-cosθ_j) [34,37,38], where θ_j is the dynamic scattering angle at dynamic scattering site j (Fig. 2). By explicitly distinguishing dynamic and static scattering events, this approach obviates assumptions, including the cumulant approximation and similarity relation that are often invoked in standard Monte Carlo approaches [39]. Even though their momentum transfer is not stored, static scattering events are important in determining the Y_dyn distribution. Unless otherwise stated, the Henyey-Greenstein phase function [Eq. (4)] is assumed for dynamic scattering.

Fig. 2. Tissue model and Monte Carlo simulation: (a) Scattering from the static tissue matrix (blue) changes light paths without imparting dynamic phase shifts, while scattering from dynamic particles, such as RBCs (red), causes dynamic phase shifts. (b) Our Monte Carlo records accumulated momentum transfer, just from dynamic scattering events, versus TOF. A weight histogram with p_dyn = 10% and g_dyn= 0.98 is shown.

Download Full Size | PDF

Monte Carlo simulations stochastically estimate the probability density P(Y_dyn,τ_s). When a photon reaches the detector, Y_dyn associated with that photon is scored in a P(Y_dyn,τ_s) histogram. Then, the overall normalized field correlation is determined as the weighted average of field correlations across all dimensionless momentum transfers experienced by detected photon paths [3,37,40]:

(18)$$g_1^{MC}({\tau _s},{\tau _d}) = \frac{{\int\limits_0^\infty {P({Y_{dyn}},} {\tau _s})\exp \left( { - \frac{1}{3}{Y_{dyn}}({\tau_s}){k^2}\left\langle {\triangle {r^2}({\tau_d})} \right\rangle } \right)d{Y_{dyn}}}}{{\int\limits_0^\infty {P({Y_{dyn}},} {\tau _s})d{Y_{dyn}}}}$$

Equation (18) assumes uncorrelated motion at different scattering events, as the dynamic momentum transfer is summed in the exponent, which is equivalent to multiplying autocorrelations for individual dynamic scattering events. Note that the displacements at all scattering events are assumed to be Gaussian, independent, and identically distributed with variance 〈Δr²(τ_d)〉. Table 2 summarizes the simulation parameters, selected from previous studies of human brain [41–44] at a wavelength of 850 nm. To study the dependency of DWS/DCS validity on anisotropy of dynamic scatterers (g_dyn), a range of g_dyn values (between 0.6 and 0.98) and corresponding speeds (v*), selected to eliminate the τ_d² term in Eq. (17), were simulated. The range of values for g_dyn encompasses particles of interest such as Intralipid (g_dyn ≈ 0.6) [32,45,46], 1 µm diameter polystyrene microspheres (g_dyn ≈ 0.9) [47–51], and red blood cells (g_dyn ≈ 0.98) [27,52,53]. Unless mentioned, other parameters are kept constant as shown in Table 2.

Table 2. Monte Carlo simulation parameters

View Table | View all tables in this article

The propagation of a photon in the medium continues until the photon escapes from the medium, or until TOF exceeds a pre-defined 1 ns threshold. Photons that escape the medium within 5 cm from the source are collected. Absorption in the medium is included by multiplying the photon weight by the albedo at every scattering event [37,54]. The TOF bin width for the simulation is 5 ps, and the autocorrelation is assigned to the center point of each time bin. The normalized TOF-resolved field autocorrelation should not depend on absorption for our narrow TOF bins. Each simulation launched 300 million photons and took an average of 35 hours to run on an Intel i9 Model 9900K 3.6-GHz processor.

3. Results

3.1 First cumulant approximation (A2)

First the consequence of ignoring higher order cumulants on the field autocorrelation is examined by comparing g₁ derived from Monte Carlo simulations (MC hybrid) to those derived from the first cumulant approximation (A2) [Eq. (14)], with hybrid motion, at two different TOF values, τ_s = 47.5 ps (Fig. 3(a)) and 197.5 ps (Fig. 3(b)), and various g_dyn values. These TOF values correspond to photon path lengths of 1.01 cm and 4.23 cm in the medium, respectively. In general, both methods agree well at small τ_d whereas the long time lag tail of A2 is lower than that of MC hybrid, with more disagreement at early TOF (Fig. 3(a)) and higher g_dyn (Fig. 3(b)-(c)). While inaccuracy of the cumulant approximation at early TOF (few dynamic scattering events) is well-known [23], its dependency on anisotropy of dynamic scatterers has not been rigorously investigated. Surprisingly, discrepancies in the tails of the autocorrelation extend beyond τ_s ∼ 500 ps for g_dyn = 0.98 (Fig. 3(c)). A zoomed-in linear scale of Fig. 3(c) shows that, for high g_dyn, neglecting higher order cumulants significantly underestimates the tails even at large TOFs (Fig. 3(d)).

Fig. 3. The cumulant approximation (A2) attenuates long autocorrelation tails. The normalized field autocorrelation function, g₁, using the first cumulant approximation for hybrid motion (A2) and Monte Carlo with hybrid motion (MC hybrid) for five g_dyn values at τ_s =47.5 ps (a) and τ_s= 197.5 ps (b), and for g_dyn = 0.98 at τ_s =497.5 ps (c-d). A linear time lag scale is shown in (d) with an equivalent dashed box in (c).

Download Full Size | PDF

The role of g_dyn in lowest order error term in the cumulant approximation is described by R, the ratio of the second cumulant term to the square of the first cumulant term [Eq. (10) and Fig. 1(b)]. Larger R values at early TOF and larger g_dyn (Fig. 1(b)) predict larger errors in the first cumulant approximation.

3.2 Neglecting random flow (A1)

Next the consequence of ignoring random flow on the field autocorrelation is examined by comparing Bonner and Nossal’s theory for Brownian motion alone (A1) to MC hybrid (Fig. 4). As described earlier, A1 accounts for all cumulant terms, but ignores random flow. A1 agrees well with MC hybrid at small τ_d, but decorrelates slower than MC hybrid at large τ_d. In other words, ignoring random flow increases the autocorrelation tails. This phenomenon is more obvious at shorter τ_s (Fig. 4(a)) and at higher g_dyn (Fig. 4(b)-(d)) because for fixed particle dynamics, the autocorrelation decays slower for high g_dyn (e.g. g_dyn = 0.98) than for low g_dyn (e.g. g_dyn = 0.6).

Fig. 4. Random flow attenuates long autocorrelation tails: g₁ using Bonner and Nossal’s theory for Brownian diffusion (A1), and Monte Carlo with hybrid motion (MC hybrid) for five selected g_dyn values (a-b) at τ_s =47.5 ps (a) and τ_s= 197.5 ps (b), and for g_dyn = 0.98 at τ_s =497.5 ps and 997.5 ps (c-d). A linear time lag scale is shown in (d) with an equivalent dashed box in (c).

Download Full Size | PDF

3.3 Conventional DWS (A1&2)

As mentioned earlier, conventional DWS ignores both random flow and higher cumulant terms (A1&2). Even though A1 and A2 are individually incorrect, Eq. (16) is still approximately satisfied as the velocity distribution width was chosen [Eq. (17)] so that it is valid up to O(τ_d³). Thus A1&2 surprisingly agrees well with the gold standard MC hybrid, especially when g_dyn < 0.8 (Fig. 5(a)), or when g_dyn = 0.98 and τ_s > 200 ps (Fig. 5(c)-(d)). The fact that A1&2 performs better than A1 and A2 is surprising because A1 and A2 include only one incorrect assumption each, whereas A1&2 includes two incorrect assumptions. This result can be explained by the fact that the autocorrelation tails can be reduced either by neglecting higher order cumulant terms (observed in A2 in Fig. 3) or by including random flow (suggested in A1 in Fig. 4). Thus, for A1&2, ignoring flow increases the tail whereas ignoring higher order cumulants reduces the tail. Therefore, in this case, two incorrect assumptions (Fig. 5) provide better predictions than just one (Fig. 3 and Fig. 4). Residual minor errors in Fig. 5 are due to O(τ_d³) terms in Eq. (13).

Fig. 5. Errors due to two incorrect assumptions can cancel: g₁ using DWS theory with the first cumulant approximation for Brownian diffusion (A1&2), and Monte Carlo with hybrid motion (MC hybrid) for five selected g_dyn values (a-b) at τ_s =47.5 ps (a) and τ_s = 197.5 ps (b), and for g_dyn =0.98 at τ_s =497.5 ps and 997.5 ps (c-d). A linear time lag scale is shown in (d) with an equivalent dashed box in (c).

Download Full Size | PDF

To illustrate the role of the first cumulant approximation in DWS theory for pure Brownian motion or diffusion, in Fig. 6, Bonner and Nossal’s theory (A1) and Monte Carlo simulations (MC diff.) are compared to DWS theory (A1&2) In this case, all cumulant terms are accounted for in MC diff. and A1, while only the first order cumulant term is accounted for in A1&2. As expected, A1 and MC diff. agree well while the tail of A1&2 is much lower, even at τ_s = 997.5 ps (Fig. 6(c)-(d)). This further confirms that neglecting higher order cumulants reduces the autocorrelation tails.

Fig. 6. For pure Brownian motion, the cumulant approximation attenuates long autocorrelation tails: g₁ using DWS theory with first cumulant approximation for Brownian diffusion (A1&2), Bonner and Nossal’s theory for Brownian diffusion (A1), and Monte Carlo for Brownian diffusion (MC diff.) for five selected g_dyn values (a-b) at τ_s =47.5 ps (a) and τ_s= 197.5 ps (b), and for g_dyn =0.98 at τ_s =497.5 ps and 997.5 ps (c). A linear time lag scale is used in (d) with an equivalent dashed box in (c).

Download Full Size | PDF

3.4 Full cumulant expansion (Bonner and Nossal) with hybrid motion

To avoid assumptions in DWS/DCS theory, Bonner and Nossal’s approach of directly integrating over the scattering phase function can be used. Thus, a B&N hybrid model accounts for all cumulant terms, and for both Brownian diffusion and random flow of dynamic scatterers. As shown in Fig. 7, B&N hybrid is an improvement over A1&2, agreeing well with MC hybrid regardless of anisotropy of dynamic scatterers (Fig. 7(a)-b), and the agreement is observed at short (Fig. 7(a)) and long τ_s (Fig. 7(c)-(d)). Thus, eliminating incorrect assumptions altogether achieves excellent agreement with the gold standard.

Fig. 7. Elimination of assumptions yields excellent agreement with MC: g₁ using DWS theory with first cumulant approximation for Brownian diffusion (A1&2), as well as Bonner and Nossal’s theory and Monte Carlo for hybrid motion (B&N hybrid and MC hybrid) for five selected g_dyn values (a-b) at τ_s = 47.5 ps (a) and τ_s= 197.5 ps (b), and for g_dyn =0.98 at τ_s =497.5 ps and 997.5 ps (c-d). A linear time lag scale is used in (d) with an equivalent dashed box in (c).

Download Full Size | PDF

4. Discussion

Using a customized Monte Carlo of light scattering in realistic tissue models as a gold standard, we highlight that errors due to two common, incorrect assumptions, one regarding the nature of RBC motion and the other regarding the form of the field autocorrelation, can cancel, yielding much better predictions of TOF-resolved field autocorrelations than just one incorrect assumption. This suggest that goodness-of-fit alone is a flawed measure by which to assess assumptions underlying a theory. By carefully examining each assumption, we arrive at a more accurate approach to assess random flow and further improve goodness-of-fit.

4.1 Accuracy of the cumulant approximation

The cumulant approximation, or first cumulant expansion, is widely used in light scattering theory, because it yields simple analytic expressions [1,7,8,23,55]. Here we argue that errors in the cumulant approximation require special attention, particularly at late time lags for highly anisotropic dynamic particles such as RBCs and, as discussed further below, in tissues with low probabilities of dynamic scattering. As such, we encourage the view that the cumulant approximation, which is central in the diffusion approximation for correlation transport, be viewed and evaluated essentially independently of the diffusion approximation for radiative transport. While the cumulant approximation can be improved by restricting fitting to early time lags (τ_d), our results show significant disagreements at late time lags when τ_s ≈ 200 ps (Fig. 3(b) and Fig. 6(b)). Importantly, as described in the next section, assessing RBC random flow requires analyzing longer time lags where the accuracy of the cumulant approximation becomes questionable.

4.2 Revealing random flow

The appropriate description of RBC motion (ordered vs. unordered) has been the topic of continued debate in the field of laser speckle [8,55]. Certainly, as discussed in Section 2.3.1, at sufficiently short time lags (i.e. τ_d << 220 μs), Brownian motion dominates random flow. Meanwhile, random flow just impacts the tails of TOF-resolved field autocorrelations (for instance, where A1 is higher than MC hybrid in Fig. 4). Therefore, assessing random flow is thus only feasible at later time lags.

As discussed in Section 2.3.2, the cumulant approximation is invalid if τ_d ≈ 141 μs or larger, precisely the lag time range needed to assess random flow. This revelation adds a new wrinkle to the debate on the nature of RBC dynamics, as the cumulant approximation is assumed in standard models for ordered and unordered flow in the literature [7]. Given our assumptions (g_dyn =0.98, v = 1.25 mm/s and D_B = 6 × 10⁻¹¹ m².s⁻¹) we hypothesize that the first cumulant approximation cannot accurately assess random flow.

To test this hypothesis, we treated the MC hybrid model as experimental data and fit either the A2 or B&N hybrid model to this data, with D_B and v as free parameters. Figure 8 compares the fitted D_B and v values to the expected (actual) values. Two different fitting ranges of the MC field autocorrelation are examined: g₁^MC > 0.01 (Fig. 8(a)-(c)), and g₁^MC > 0.5 (Fig. 8(d)-(f)). B&N hybrid accurately recovers D_B and v while A2 underestimates D_B and v. As A2 ignores higher order cumulants, the fitting procedure reduces D_B to approximate the long time lag tail in g₁^MC. With a restricted fitting range (to exclude the tail of g₁^MC), A2 more accurately estimates D_B (Fig. 8(d)). However, notably, A2 always yields v∼0 mm/s, suggesting that A2 is a poor model for quantifying random flow.

Fig. 8. The commonly-used cumulant approximation (A2) cannot accurately recover dynamics from Monte Carlo simulations (Section 2.4) with random flow and diffusion (MC hybrid). Fitted D_B (a, d). and v (b, e) values are underestimated by A2 (blue line). Eliminating the cumulant approximation with a B&N hybrid model enables accurate estimates (black line), in agreement with expected values used in the simulations (red line).

Download Full Size | PDF

Thus, particular caution is warranted in applying the first cumulant approximation for small numbers of dynamic scattering events. Thus, for laser speckle flowmetry in parenchymal regions with sparse blood vessels, the range of autocorrelation lags (or the range of exposure times) must be chosen judiciously to avoid model errors.

4.3 Dynamic scattering probability dependence

In the above sections, we show that one confound in applying DWS/DCS theory to tissue measurements is the exceptionally high g_dyn of RBCs. Here, we highlight another important confound; the low dynamic scattering probability (p_dyn) of tissues. As mentioned in Section 2.4, a smaller V_B implies a smaller p_dyn. To illustrate the effect of p_dyn on DWS/DCS accuracy, p_dyn = 30% (V_B = 3.9%, Fig. 9(b)) was compared to p_dyn = 10% (V_B = 1.3%, Fig. 9(a)), used heretofore in this work. Obviously, smaller p_dyn values result in a slower autocorrelation decay due to the smaller number of dynamic scattering events. For both p_dyn, A1&2, representing conventional DWS/DCS, underestimates the autocorrelation tails of MC hybrid, which are well-represented by B&N hybrid. Importantly, disagreement between A1&2 and MC hybrid persists to later TOFs for smaller p_dyn, due to the smaller number of dynamic scattering events. For example, A1&2 agrees well with MC hybrid when p_dyn = 30% at 497.5 ps (Fig. 9(b)), but still decays faster when p_dyn = 10% (Fig. 9(a)).

Fig. 9. Accuracy of DWS/DCS assumptions depend on the probability of dynamic scattering, p_dyn: g₁ from less accurate models (the first cumulant approximation for Brownian motion (A1&2) and B&N diffusion (A1)), disagrees with more accurate models (B&N hybrid and MC hybrid, the gold standard) up to longer τ_s for p_dyn = 10% (a) than for p_dyn = 30% (b) when g_dyn₌0.98.

Download Full Size | PDF

To quantify agreement between various models as p_dyn changes, the coefficient of determination (r-squared) versus τ_s is shown (Fig. 10). The MC hybrid model is taken as the gold standard, and hence has an r-squared of unity. Different ranges of the field autocorrelation function are examined, corresponding to g₁^MC > 0.01, 0.2, and 0.5. As expected, B&N hybrid agrees best, with an r-squared approaching unity in all cases, while A1&2 (conventional DWS/DCS) is the second best, followed by A1 (B&N diffusion). Limiting the fitting range greatly improves the r-squared for A1, especially at early TOF. For example, the lowest r-squared value for A1 is well below 0.4 (0.9) if the tail of autocorrelation function is included (Fig. 10(a),(d)), and above 0.4 (0.9) if the tail is excluded (Fig. 10(c),(f)) when p_dyn = 10% (30%). These results are expected since the first cumulant term for Brownian motion dominates early time lags.

Fig. 10. Quantification of results in Fig. 9. Coefficient of determination (r-squared) based on comparing the gold standard MC hybrid model to different models: A1, B&N hybrid, and A1&2 when p_dyn = 10% (a-c), and p_dyn = 30% (d-f). Note that by definition, the MC hybrid model has an r-squared of unity across TOFs.

Download Full Size | PDF

Increasing p_dyn (the number of dynamic scattering events) improves the r-squared for A1 and A1&2. For instance, A1 yields r-squared below 0.6 at early TOFs when p_dyn = 10% (Fig. 10(a)), and well above 0.6 at early TOFs when p_dyn = 30% (Fig. 10(d)). The improvement in A1 is expected since with more dynamic scattering events, the first cumulant term for Brownian motion dominates more of the autocorrelation decay. The improvement in A1&2 is also expected since higher order cumulants are less important with more dynamic scattering events.

4.4 Proposed improved expression

A common approach (COM) to determine field autocorrelations is a weighted average over total dimensionless momentum transfer Y_tot [4,40,56]:

(19)$$g_1^{\textrm{COM}}({\tau _s},{\tau _d}) = \frac{{\int\limits_0^\infty {P({Y_{tot}},} {\tau _s})\exp \left( { - \frac{1}{3}\alpha {Y_{tot}}({\tau_s}){k^2}\left\langle {\triangle {r^2}({\tau_d})} \right\rangle } \right)d{Y_{tot}}}}{{\int\limits_0^\infty {P({Y_{tot}},} {\tau _s})d{Y_{tot}}}}$$

In homogenous media, note that this approach is equivalent to DWS, if Y_tot = μ_s(1-g_stat)s, where s=τ_sc/n is the photon path length in the medium. In Eq. (19), Y_tot is the total dimensionless momentum accumulated from all scattering events and can be obtained from a standard Monte Carlo method that does not explicitly simulate dynamic scattering [38]. To account for dynamic scattering events, the momentum from simulations is multiplied by α (note that α ≈ 0.2p_dyn when g_dyn = 0.98 and g_stat =0.9). This approach implicitly assumes the cumulant approximation and a simple similarity relationship between dynamic and total momentum transfer: Y_dyn = αY_tot.

Here, we argue that the most accurate approach is to use a Monte Carlo simulation that includes both static and dynamic scattering events (Section 2.4), as this approach correctly performs angular averaging over the scattering phase function, implicitly including all cumulant orders. To achieve accurate results with a standard Monte Carlo code that does not explicitly include separate dynamic and static scattering events, we propose an improved expression (IMP):

(20)$$g_1^{\textrm{IMP}}({\tau _s},{\tau _d}) = \frac{{\int\limits_0^\infty {P({Y_{tot}},} {\tau _s}){{[{g_1^{ss}({\tau_d})} ]}^{\overline m ({{Y_{tot}}} )}}d{Y_{tot}}}}{{\int\limits_0^\infty {P({Y_{tot}},} {\tau _s})d{Y_{tot}}}}$$

Rewriting Eqs. (19) and (20), respectively, to take a weighted average of field autocorrelations over all detected photon paths, each with weight w_i [37,57], we arrive at

(21)$$g_1^{\textrm{COM}}({\tau _s},{\tau _d}) = \frac{{\sum\limits_{ - \delta {\tau _s} \le {\tau _{si}} - {\tau _s} < \delta {\tau _s}}^{} {{w_i}\exp \left( { - \frac{1}{3}{k^2}\left\langle {\triangle {r^2}({\tau_d})} \right\rangle \alpha {Y_{tot,i}}} \right)} }}{{\sum\limits_{ - \delta {\tau _s} \le {\tau _{si}} - {\tau _s} < \delta {\tau _s}}^{} {{w_i}} }}$$

and

(22)$$g_1^{\textrm{IMP}}({\tau _s},{\tau _d}) = \frac{{\sum\limits_{ - \delta {\tau _s} \le {\tau _{si}} - {\tau _s} < \delta {\tau _s}}^{} {{w_i}({\tau _s}){{[{g_1^{ss}({\tau_d})} ]}^{{{\overline m }_i}}}} }}{{\sum\limits_{ - \delta {\tau _s} \le {\tau _{si}} - {\tau _s} < \delta {\tau _s}}^{} {{w_i}} }}.$$

In Eqs. (21) and (22), w_i is the weight and τ_si is the TOF for the detected photon i, g₁^ss is the normalized single scattering autocorrelation function [Eq. (2)], $\bar{m}$_i = αY_tot_,i /(1-g_dyn) is the average number of dynamic scattering events (Alternatively, similar results are obtained with $\bar{m}$_i = p_dynµ_ss_i). To compare the common method [Eq. (21)] to our more accurate customized Monte Carlo method (MC hybrid), and to validate the proposed improvement [Eq. (22)], we record Y_tot using the GPU-based code Monte Carlo eXtreme (MCX) developed by Fang et al. [38]. The MCX simulation, which accumulates Y_tot from all scattering events [39], is used to determine the field autocorrelation for both the common approach [Eq. (21)] and the improved approach [Eq. (22)]. The simulation launches 10⁸ photons at the source, which took approximately 5 minutes on an NVIDIA Dual GeForce RTX 2080 GPU. The simulation returns a list of all detected photons and their individual momentum transfers and partial paths. Absorption is included by applying the Beer-Lambert law to determine the photon weight, w_i. To approximate the 5 cm detection radius in the gold standard MC simulations (Section 2.4), in the MCX simulations, 11 detectors with radii of 0.5 mm were evenly spaced from 0 to 5 cm from the source. All detected photons were included when applying in Eqs. (21) and (22).

As shown in Fig. 11, the proposed method (IMP), with a more accurate autocorrelation (adapted from Bonner and Nossal), closely approximates our customized Monte Carlo method (MC hybrid) which explicitly simulates dynamic and static scattering, resulting in the best agreement with near unity r-squared both when p_dyn = 10% (α = 2%) [Fig. 11(e)-(f)], and when p_dyn = 30% (α = 6%) [Fig. 11(g)-(h)]. On the other hand, the common approach (COM hybrid), which assumes the cumulant approximation, underestimates the autocorrelation tails when τ_s is less than 300 ps. While agreement improves for τ_s > 300 ps, these discrepancies are significant because the cumulant approximation is widely assumed to be valid in DCS/DWS across a range of source-detector separations.

Fig. 11. Performance of proposed method to eliminate DWS/DCS assumptions and improve accuracy. (a-d) g₁ for hybrid motion with our gold-standard Monte Carlo simulations (MC hybrid), the common Monte Carlo method (COM hybrid), and our proposed method (IMP hybrid) for g_dyn₌0.98 when (a,b) p_dyn = 10% (α = 2%) and (c,d) p_dyn = 30% (α = 6%). (e-h) Corresponding r-squared values, determined by comparing COM hybrid and IMP hybrid to MC hybrid, for two autocorrelation ranges when p_dyn = 10% (e,f) and p_dyn = 30% (g,h). A linear time lag scale is used in (b,d) with equivalent dashed boxes in (a,c).

Download Full Size | PDF

In summary, our proposed method [Eq. (22)]: 1) improves the underestimation of autocorrelation tails observed in the common method (and DWS/DCS theory), and 2) enables the use of highly optimized existing Monte Carlo tools. Another potential advantage of the proposed method is that it can be extended to multi-layered tissues, by tracking momentum transfer or partial path length in each tissue type, as is possible in MCX.

4.5 Gegenbauer phase function

Bonner and Nossal’s theory and the MC method incorporate all cumulants and depend on the entire scattering phase function, unlike DWS theory (A1&2), which just depends on 〈cosθ〉. In this section, we examine the effect of phase function choice on the autocorrelation. Early goniophotometric measurement of RBC suspensions showed that light scattering by RBCs is well-described by a Gegenbauer (GE) phase function [58–61]. The GE phase function which is given by:

(23)$${p_{GE}}(\theta ) = \frac{{a{g_{GE}}}}{\pi }\frac{{{{({1 - g_{GE}^2} )}^{2a}}}}{{[{{{({1 + {g_{GE}}} )}^{2a}} - {{({1 - {g_{GE}}} )}^{2a}}} ]{{[{1 + g_{GE}^2 - 2{g_{GE}}\cos \theta } ]}^{1 + a}}}}$$

Equation (23) requires |g_GE| ≤ 1 and a >- 0.5. Note that g_GE ≠ 〈cosθ〉 in Eq. (23), and that only if a = 0.5, Eq. (23) reduces to Eq. (4), with g_GE = g_dyn = 〈cosθ〉. Thus, the GE phase function is a generalization of the HG phase function. If 〈cosθ〉 = 0.98 and a = 1.2 [58], we obtain that g_GE = 0.893. As shown in Fig. 12(a), HG and GE phase functions are approximately the same when 〈cosθ〉 = 0.6 but are significantly different when 〈cosθ 〉 = 0.98, especially at large scattering angles. The cumulant approximation error, R, is 2.35 times lower at 〈cosθ 〉 = 0.98 for GE (Fig. 12(b)).

Fig. 12. Cumulant approximation error decreases when the Gegenbauer phase function (GE, a = 1.2) is used. (a) Henyey-Greenstein (HG) and Gegenbauer phase functions for 〈cosθ 〉 = 0.6 and 0.98, (b) R versus 〈cosθ 〉 for different $\bar{m}$ values, p_dyn = 10%, and μ_s = 100 cm⁻¹. The inset in (b) shows the cumulant approximation error ratio between HG and GE, which is independent of $\bar{m}$. The upper x-axis in (b) shows the g_GE value corresponding to 〈cosθ 〉 on the lower x-axis.

Download Full Size | PDF

As above, the velocity distribution v* is selected so that the first cumulant term for random flow cancels second cumulant term for Brownian diffusion. Like R, v* depends on phase function [Eq. (17)], and v* = 0.817 mm/s when a = 1.2, g_GE = 0.893 and D_B = 6 × 10⁻¹¹ m².s⁻¹ [26,33,34]. Figure 13 shows that while the autocorrelation tails are greater for the HG phase function (B&N hybrid, HG) than the GE phase function (MC hybrid, GE and B&N hybrid, GE), all have longer tails than DWS (A2), which ignores higher order cumulants. The smaller cumulant approximation error for the GE phase function is explained by smaller R values in Fig. 12(b).

Fig. 13. Gegenbauer (GE) phase function (a = 1.2, g_GE = 0.893) reduces autocorrelation tails compared to the Henyey-Greenstein (HG) phase function even though 〈cosθ 〉 = 0.98 for both: g₁ from DWS theory with first cumulant approximation for hybrid motion (A2) decays faster than accurate models (B&N hybrid and MC hybrid), but disagrees less with GE autocorrelations, as expected based on Fig. 12.

Download Full Size | PDF

Figure 14 also confirms that ignoring random flow (A1, GE) increases the autocorrelation tail. Thus assuming a GE rather than HG phase function does not change the essential conclusions regarding autocorrelation tail; namely, that neglecting higher order cumulants reduces the tail, neglecting random flow increases the tails, and that both assumptions in tandem lead to better agreement with the gold standard MC simulation. Moreover, Fig. 14 shows that conventional DWS/DCS (A1&2) with two incorrect assumptions agrees well with the more accurate models that use the GE phase function for both p_dyn = 10% and p_dyn = 30%.

Fig. 14. A replicate of Fig. 9 with the GE phase function (a = 1.2, g_GE = 0.893): once again g₁ from less accurate models (A1&2 and A1), disagrees with more accurate models (B&N hybrid and MC hybrid, the gold standard) up to longer τ_s for p_dyn = 10% (a) than for p_dyn = 30% (b) when 〈cosθ 〉 = 0.98

Download Full Size | PDF

Figure 15 shows r-squared versus τ_s for different ranges of the field autocorrelation. As before, the MC model is taken as the gold standard. Like Fig. 10, Fig. 15 shows that MC hybrid, GE agrees well with A1&2 (two incorrect assumptions) at all TOFs but disagrees with A1, GE (one incorrect assumption) at earlier TOFs, and that r-squared for A1, GE improves in tissues with higher p_dyn, and when the range is limited to early time lags.

Fig. 15. Quantification of results of Fig. 14. Coefficient of determination (r-squared) based on comparing the gold standard MC hybrid, GE model to different models: A1, GE; B&N hybrid, GE; and A1&2 when p_dyn = 10% (a-c), and p_dyn = 30% (d-f). Note that by definition, the MC hybrid, GE model has an r-squared of unity across TOFs. As in Fig. 10, A1&2 (2 errors) agrees better with the gold standard at both p_dyn values than A1, GE (1 error), and here, notably, overlaps with B&N hybrid, GE.

Download Full Size | PDF

4.6 Simplifying assumptions

To arrive at our proposed improved expression for field autocorrelations in biological tissue [Eq. (22)], several assumptions were adapted from the original Bonner and Nossal theory [24]. Specifically, it was assumed that RBCs move independently, that averaging over photon paths and scattering events is equivalent to ensemble averaging, that the distances between static tissue source points and moving cells is large compared to displacements for decorrelation, and that static scattering randomizes the direction of incident light between dynamic scattering events. In incorporating a hybrid model of RBC motion, we have also assumed that advective and Brownian components of the motion are independent.

Importantly, vascular geometry was not explicitly simulated in the current study. RBCs were assumed to be uniformly distributed throughout the tissue volume, parameterized by the probability of dynamic scattering (p_dyn). In a previous study with a proscribed parallel vessel geometry, p_dyn was empirically found to be proportional to vessel radius-squared and inversely proportional to vessel spacing-squared [34]. While our simplified tissue model enabled us to investigate the cumulant approximation without assuming a specific geometry, future studies should investigate realistic vascular geometries, accounting for the issues highlighted in the present work.

5. Conclusion

We have shown that field autocorrelations in realistic biological tissue models with hybrid motion are accurately predicted by our adaptation of Bonner and Nossal’s original theory for RBC scattering. Importantly, neglecting higher order cumulants reduces the field autocorrelation tails whereas neglecting random flow increases them. Thus, for particular choices of random flow distribution and Brownian diffusion coefficient, DWS theory with pure Brownian motion might agree with experimental data, leading to the incorrect conclusion that random flow is absent. Furthermore, if experimental data are fit using the DWS cumulant approximation, allowing for hybrid motion, dynamics are underestimated. Here, we argue that dynamics can be accurately assessed by analyzing field autocorrelations with a model that includes higher order cumulants. Altogether, this work supports that random flow may be assessable in TOF-resolved field autocorrelations, but higher order cumulants should be incorporated into models to ensure its accurate assessment.

Appendix A

Here, we derive Eq. (6). In Bonner and Nossal’s original theory [24], the normalized field autocorrelation function can be expressed as

(A.1)$$g_1^{BN,orig}({\tau _s},{\tau _d}) = \sum\limits_{m = 1}^\infty {{{{P_m}({\tau _s}){I_m}({\tau _d})} \mathord{\left/ {\vphantom {{{P_m}({\tau_s}){I_m}({\tau_d})} {({1 - {P_0}} )}}} \right.} {({1 - {P_0}} )}}} ,$$

where I_m is the autocorrelation function for a path, and P_m is the probability of such a path. When the lower limit in the summation is set to zero to include purely static scattering paths, we re-write Eq.(A.1):

(A.2)$$g_1^{BN}({\tau _s},{\tau _d}) = \sum\limits_{m = 0}^\infty {{P_m}({\tau _s}){I_m}({\tau _d})} .$$

If we incorporate our symbol for the autocorrelation for a single dynamic scattering event, g₁^ss, then

(A.3)$${I_m}({\tau _d}) = {[{g_1^{ss}({\tau_d})} ]^m}.$$

We then rewrite Eq. (A.2) as:

(A.4)$$g_1^{BN}({\tau _s},{\tau _d}) = \sum\limits_{m = 0}^\infty {{P_m}({\tau _s}){{[{g_1^{ss}({\tau_d})} ]}^m}} .$$

Assuming that P_m follows a Poisson distribution, we rewrite Eq. (A.4):

(A.5)$$g_1^{BN}({\tau _s},{\tau _d}) = \sum\limits_{m = 0}^\infty {{{\frac{{{e^{ - \mathop m\limits^\_ ({\tau _s})}}\left[ {\mathop m\limits^\_ ({\tau_s})g_1^{ss}({\tau_d})} \right]}}{{m!}}}^m}} .$$

Finally, applying a Taylor series, Eq. (A.5) simplifies to Eq. (6).

Appendix B

Here, we show how Eqs. (1) and (9) are derived by applying cumulant expansion to Bonner and Nossal’s theory [Eq.(6)]. First, we take logarithm both sides of Eq.(6):

(B.1)$$\ln ({g_1^{BN}({\tau_s},{\tau_d})} )= [{\overline m ({\tau_s})\{{g_1^{ss}({\tau_d}) - 1} \}} ]= \left[ {\overline m ({\tau_s})\left\{ {{{\left\langle {\exp \left[ { - \frac{{{q^2}(\theta )}}{6}\left\langle {\triangle {r^2}({\tau_d})} \right\rangle } \right]} \right\rangle }_\theta } - 1} \right\}} \right]$$

Next, Taylor series expansion is applied to Eq. (B.1) up to second order to give:

(B.2)$$\begin{array}{c} \ln ({g_1^{BN}({\tau_s},{\tau_d})} )= \overline m ({\tau _s})\left[ {{{\left\langle {1 - \frac{{{q^2}(\theta )}}{6}\left\langle {\triangle {r^2}({\tau_d})} \right\rangle + \frac{1}{2}\frac{{{q^4}(\theta )}}{{{6^2}}}{{\left\langle {\triangle {r^2}({\tau_d})} \right\rangle }^2} + O\left( {{{\left\langle {\triangle {r^2}({\tau_d})} \right\rangle }^3}} \right)} \right\rangle }_\theta } - 1} \right],\\ \textrm{ } ={-} \overline m ({\tau _s})\frac{{{{\left\langle {{q^2}(\theta )} \right\rangle }_\theta }}}{6}\left\langle {\triangle {r^2}({\tau_d})} \right\rangle \left[ {1 - \frac{1}{6}\frac{{{{\left\langle {{q^4}(\theta )} \right\rangle }_\theta }}}{{2{{\left\langle {{q^2}(\theta )} \right\rangle }_\theta }}}\left\langle {\triangle {r^2}({\tau_d})} \right\rangle } \right] + \overline m ({\tau _s})O\left( {{{\left\langle {\triangle {r^2}({\tau_d})} \right\rangle }^3}} \right)\\ \end{array}, $$

or

(B.3)$$\textrm{ }\ln ({g_1^{BN}({\tau_s},{\tau_d})} )={-} \overline m ({\tau _s})\frac{{{{\left\langle {{q^2}(\theta )} \right\rangle }_\theta }}}{6}\left\langle {\triangle {r^2}({\tau_d})} \right\rangle \left[ {1 - \frac{1}{6}\overline m ({\tau_s})R{{\left\langle {{q^2}(\theta )} \right\rangle }_\theta }\left\langle {\triangle {r^2}({\tau_d})} \right\rangle } \right] + \overline m ({\tau _s})O\left( {{{\left\langle {\triangle {r^2}({\tau_d})} \right\rangle }^3}} \right),$$

where $R = {{{{\left\langle {{q^4}(\theta )} \right\rangle }_\theta }} \mathord{\left/ {\vphantom {{{{\left\langle {{q^4}(\theta )} \right\rangle }_\theta }} {\left[ {2\overline m ({\tau_s}){{\left\langle {{q^2}(\theta )} \right\rangle }_\theta }^2} \right]}}} \right.} {\left[ {2\overline m ({\tau_s}){{\left\langle {{q^2}(\theta )} \right\rangle }_\theta }^2} \right]}}$. Note that third and higher cumulant terms are contained in O(〈Δr²(τ_d)〉³). Taking the exponential of both sides of Eq. (B.3), and retaining only the first order cumulant term yields:

(B.4)$$g_1^{BN}({\tau _s},{\tau _d}) \approx \exp \left[ { - \overline m ({\tau_s})\frac{{{{\left\langle {{q^2}(\theta )} \right\rangle }_\theta }}}{6}\left\langle {\triangle {r^2}({\tau_d})} \right\rangle } \right]$$

Equation (B.4) is equivalent to Eq. (1) in the case of Brownian motion. Meanwhile, if both first and second cumulants are retained, the exponential of Eq. (B.3), after dropping θ subscripts, yields Eq. (9).

Funding

National Institutes of Health (R01NS094681, R03EB023591, R21NS105043).

Acknowledgments

We thank Oybek Kholiqov and Wenjun Zhou for valuable discussions.

Disclosures

The authors declare no conflicts of interest.

References

1. D. J. Pine, D. A. Weitz, P. M. Chaikin, and E. Herbolzheimer, “Diffusing-wave spectroscopy,” Phys. Rev. Lett. 60(12), 1134–1137 (1988). [CrossRef]

2. G. Maret and P. E. Wolf, “Multiple light scattering from disordered media. The effect of Brownian motion of scatterers,” Z. Phys. B: Condens. Matter 65(4), 409–413 (1987). [CrossRef]

3. 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(1), 192–215 (1997). [CrossRef]

4. T. Durduran, R. Choe, W. B. Baker, and A. G. Yodh, “Diffuse optics for tissue monitoring and tomography,” Rep. Prog. Phys. 73(7), 076701 (2010). [CrossRef]

5. J. Li, G. Dietsche, D. Iftime, S. E. Skipetrov, G. Maret, T. Elbert, B. Rockstroh, and T. Gisler, “Noninvasive detection of functional brain activity with near-infrared diffusing-wave spectroscopy,” J. Biomed. Opt. 10(4), 044002 (2005). [CrossRef]

6. T. Durduran, G. Yu, M. G. Burnett, J. A. Detre, J. H. Greenberg, J. Wang, C. Zhou, and A. G. Yodh, “Diffuse optical measurements of blood flow, blood oxygenation and metabolism in human brain during sensorimotor cortex activation,” Opt. Lett. 29(15), 1766–1768 (2004). [CrossRef]

7. R. Bandyopadhyay, A. S. Gittings, S. S. Suh, P. K. Dixon, and D. J. Durian, “Speckle-visibility spectroscopy: A tool to study time-varying dynamics. Review of scientific instruments,” Rev. Sci. Instrum. 76(9), 093110 (2005). [CrossRef]

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

9. J. Sutin, B. Zimmerman, D. Tyulmankov, D. Tamborini, K. C. Wu, J. Selb, A. Gulinatti, I. Rech, A. Tosi, D. A. Boas, and M. A. Franceschini, “Time-domain diffuse correlation spectroscopy,” Optica 3(9), 1006–1013 (2016). [CrossRef]

10. D. Borycki, O. Kholiqov, and V. J. Srinivasan, “Reflectance-mode interferometric near-infrared spectroscopy quantifies brain absorption, scattering, and blood flow index in vivo,” Opt. Lett. 42(3), 591–594 (2017). [CrossRef]

11. O. Kholiqov, W. Zhou, T. Zhang, V. N. D. Le, and V. J. Srinivasan, “Time-of-flight resolved light field fluctuations reveal deep human tissue physiology,” Nat. Commun. 11(1), 391 (2020). [CrossRef]

12. C. Cheung, J. P. Culver, K. Takahashi, J. H. Greenberg, and A. G. Yodh, “In vivo cerebrovascular measurement combining diffuse near-infrared absorption and correlation spectroscopies,” Phys. Med. Biol. 46(8), 2053–2065 (2001). [CrossRef]

13. C. Zhou, G. Yu, D. Furuya, J. H. Greenberg, A. G. Yodh, and T. Durduran, “Diffuse optical correlation tomography of cerebral blood flow during cortical spreading depression in rat brain,” Opt. Express 14(3), 1125–1144 (2006). [CrossRef]

14. C. Zhou, S. A. Eucker, T. Durduran, G. Yu, J. Ralston, S. H. Friess, R. N. Ichord, S. S. Margulies, and A. G. Yodh, “Diffuse optical monitoring of hemodynamic changes in piglet brain with closed head injury,” J. Biomed. Opt. 14(3), 034015 (2009). [CrossRef]

15. E. M. Buckley, N. M. Cook, T. Durduran, M. N. Kim, C. Zhou, R. Choe, G. Yu, S. Shultz, C. M. Sehgal, D. J. Licht, and P. H. Arger, “Cerebral hemodynamics in preterm infants during positional intervention measured with diffuse correlation spectroscopy and transcranial doppler ultrasound,” Opt. Express 17(15), 12571–12581 (2009). [CrossRef]

16. G. Yu, T. Floyd, T. Durduran, C. Zhou, J. J. Wang, J. A. Detre, and A. G. Yodh, “Validation of diffuse correlation spectroscopy for muscle blood flow with concurrent arterial-spin-labeling perfusion,” Opt. Express 15(3), 1064–1075 (2007). [CrossRef]

17. T. Durduran, C. Zhou, B. L. Edlow, G. Yu, R. Choe, M. N. Kim, B. L. Cucchiara, M. E. Putt, Q. Shah, S. E. Kasner, and J. H. Greenberg, “Transcranial optical monitoring of cerebrovascular hemodynamics in acute stroke patients,” Opt. Express 17(5), 3884–3902 (2009). [CrossRef]

18. J. Li, M. Ninck, L. Koban, T. Elbert, J. Kissler, and T. Gisler, “Transient functional blood flow change in the human brain measured noninvasively by diffusing-wave spectroscopy,” Opt. Lett. 33(19), 2233–2235 (2008). [CrossRef]

19. C. Zhou, R. Choe, N. S. Shah, T. Durduran, G. Yu, A. Durkin, D. Hsiang, R. Mehta, J. A. Butler, A. E. Cerussi, and B. J. Tromberg, “Diffuse optical monitoring of blood flow and oxygenation in human breast cancer during early stages of neoadjuvant chemotherapy,” J. Biomed. Opt. 12(5), 051903 (2007). [CrossRef]

20. T. Durduran, R. Choe, G. Yu, C. Zhou, J. C. Tchou, B. J. Czerniecki, and A. G. Yodh, “Diffuse optical measurement of blood flow in breast tumors,” Opt. Lett. 30(21), 2915–2917 (2005). [CrossRef]

21. U. Sunar, H. Quon, T. Durduran, J. Zhang, J. Du, C. Zhou, G. Yu, R. Choe, A. Kilger, R. A. Lustig, and L. A. Loevner, “Non-invasive diffuse optical measurement of blood flow and blood oxygenation for monitoring radiation therapy in patients with head and neck tumors: a pilot study,” J. Biomed. Opt. 11(6), 064021 (2006). [CrossRef]

22. G. Yu, T. Durduran, G. Lech, C. Zhou, B. Chance, E. R. Mohler, and A. G. Yodh, “Time-dependent blood flow and oxygenation in human skeletal muscle measured with noninvasive near-infrared diffuse optical spectroscopies,” J. Biomed. Opt. 10(2), 024027 (2005). [CrossRef]

23. D. J. Pine, D. A. Weitz, J. X. Zhu, and E. Herbolzheimer, “Diffusing-wave spectroscopy: dynamic light scattering in the multiple scattering limit,” J. Phys. 51(18), 2101–2127 (1990). [CrossRef]

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

25. J. M. Higgins, D. T. Eddington, S. N. Bhatia, and L. Mahadevan, “Statistical dynamics of flowing red blood cells by morphological image processing,” PLoS Comput. Biol. 5(2), e1000288 (2009). [CrossRef]

26. J. J. Bishop, A. S. Popel, M. Intaglietta, and P. C. Johnson, “Effect of aggregation and shear rate on the dispersion of red blood cells flowing in venules,” Am. J. Physiol.-Heart C. 283(5), H1985–H1996 (2002). [CrossRef]

27. M. Friebel, J. Helfmann, U. J. Netz, and M. C. Meinke, “Influence of oxygen saturation on the optical scattering properties of human red blood cells in the spectral range 250 to 2000nm,” J. Biomed. Opt. 14(3), 034001 (2009). [CrossRef]

28. S. L. Jacques and S. A. Prahl, “Henyey-Greenstein scattering function,” https://omlc.org/classroom/ece532/class3/hg.html

29. A. J. Welch, M. J. C. van Gemert, and W. M. Star, “Definitions and Overview of Tissue Optics,” in Optical-Thermal Response of Laser-Irradiated Tissue, A. J. Welch and M. J. C. van Gemert, eds. (Springer Science, 2011)

30. S. A. Prahl, “A Monte Carlo model of light propagation in tissue,” Dosimetry of laser radiation in medicine and biology. Vol. 10305. International Society for Optics and Photonics (1989).

31. T. Binzoni, T. S. Leung, A. H. Gandjbakhche, D. Ruefenacht, and D. T. Delpy, “The use of the Henyey-Greenstein phase function in Monte Carlo simulations in biomedical optics.”,” Phys. Med. Biol. 51(17), N313–N322 (2006). [CrossRef]

32. B. W. Pogue and M. S. Patterson, “Review of tissue simulating phantoms for optical spectroscopy, imaging and dosimetry,” J. Biomed. Opt. 11(4), 041102 (2006). [CrossRef]

33. H. L. Goldsmith and J. C. Marlow, “Flow behavior of erythrocytes. II. Particle motions in concentrated suspensions of ghost cells,” J. Colloid Interface Sci. 71(2), 383–407 (1979). [CrossRef]

34. D. Boas, S. Sakadžić, J. J. Selb, P. Farzam, M. A. Franceschini, and S. A. Carp, “Establishing the diffuse correlation spectroscopy signal relationship with blood flow,” Neurophotonics 3(3), 031412 (2016). [CrossRef]

35. J. Autio, H. Kawaguchi, S. Saito, I. Aoki, T. Obata, K. Masamoto, and I. Kanno, “Spatial frequency-based analysis of mean red blood cell speed in single microvessels: investigation of microvascular perfusion in rat cerebral cortex,” PLoS One 6(8), e24056 (2011). [CrossRef]

36. I. Gurov, M. Volkov, N. Margaryants, A. Pimenov, and A. Potemkin, “High-speed video capillaroscopy method for imaging and evaluation of moving red blood cells,” Opt. Laser. Eng. 104, 244–251 (2018). [CrossRef]

37. S. Prahl, “ Monte Carlo Light Scattering Programs,” https://omlc.org/software/mc/

38. Q. Fang and D. A. Boas, “Monte Carlo simulation of photon migration in 3D turbid media accelerated by graphics processing units,” Opt. Express 17(22), 20178–20190 (2009). [CrossRef]

39. E. Sathialingam, S. Y. Lee, B. Sanders, J. Park, C. E. McCracken, L. Bryan, and E. M. Buckley, “Small separation diffuse correlation spectroscopy for measurement of cerebral blood flow in rodents,” Biomed. Opt. Express 9(11), 5719–5734 (2018). [CrossRef]

40. D. J. Durian, “Accuracy of diffusing-wave spectroscopy theories,” Phys. Rev. E 51(4), 3350–3358 (1995). [CrossRef]

41. F. Bevilacqua, D. Piguet, P. Marquet, J. D. Gross, B. J. Tromberg, and C. Depeursinge, “In vivo local determination of tissue optical properties: applications to human brain,” Appl. Opt. 38(22), 4939–4950 (1999). [CrossRef]

42. J. J. Selb, D. A. Boas, S. Chan, K. C. Evans, E. M. Buckley, and S. A. Carp, “Sensitivity of near-infrared spectroscopy and diffuse correlation spectroscopy to brain hemodynamics: simulations and experimental findings during hypercapnia,” Neurophotonics 1(1), 015005 (2014). [CrossRef]

43. K. L. Leenders, D. Perani, A. A. Lammertsma, J. D. Heather, P. Buckingham, T. Jones, M. J. R. Healy, J. M. Gibbs, R. J. S. Wise, J. Hatazawa, and S. Herold, “Cerebral blood flow, blood volume and oxygen utilization: normal values and effect of age,” Brain 113(1), 27–47 (1990). [CrossRef]

44. S. L. Jacques, “Optical properties of biological tissues: a review,” Phys. Med. Biol. 58(11), R37–R61 (2013). [CrossRef]

45. H. J. van Staveren, C. J. Moes, J. van Marie, S. A. Prahl, and M. J. van Gemert, “Light scattering in Intralipid- 10% in the wavelength range of 400-1100 nm,” Appl. Opt. 30(31), 4507–4514 (1991). [CrossRef]

46. S. T. Flock, S. L. Jacques, B. C. Wilson, W. M. Star, and M. J. van Gemert, “Optical properties of Intralipid: a phantom medium for light propagation studies,” Lasers Surg. Med. 12(5), 510–519 (1992). [CrossRef]

47. N. Rajaram, T. H. Nguyen, and J. W. Tunnell, “Lookup table–based inverse model for determining optical properties of turbid media,” J. Biomed. Opt. 13(5), 050501 (2008). [CrossRef]

48. Q. Wang, K. Shastri, and T. J. Pfefer, “Experimental and theoretical evaluation of a fiber-optic approach for optical property measurement in layered epithelial tissue,” Appl. Opt. 49(28), 5309–5320 (2010). [CrossRef]

49. N. Rajaram, T. J. Aramil, K. Lee, J. S. Reichenberg, T. H. Nguyen, and J. W. Tunnell, “Design and validation of a clinical instrument for spectral diagnosis of cutaneous malignancy,” Appl. Opt. 49(2), 142–152 (2010). [CrossRef]

50. V. N. Du Le, Z. Nie, J. E. Hayward, T. J. Farrell, and Q. Fang, “Measurements of extrinsic fluorescence in intralipid and polystyrene microspheres,” Biomed. Opt. Express 5(8), 2726–2735 (2014). [CrossRef]

51. K. Vishwanath and S. Zanfardino, “Diffuse Correlation Spectroscopy at Short Source-Detector Separations: Simulations, Experiments and Theoretical Modeling,” Appl. Sci. 9(15), 3047 (2019). [CrossRef]

52. A. Kienle, M. S. Patterson, L. Ott, and R. Steiner, “Determination of the scattering coefficient and the anisotropy factor from laser Doppler spectra of liquids including blood,” Appl. Opt. 35(19), 3404–3412 (1996). [CrossRef]

53. M. Friebel, A. Roggan, G. Müller, and M. Meinke, “Determination of optical properties of human blood in the spectral range 250 to 1100 nm using Monte Carlo simulations with hematocrit-dependent effective scattering phase function,” J. Biomed. Opt. 11(3), 034021 (2006). [CrossRef]

54. S. L. Jacques, “Monte Carlo modeling of light transport in tissue (steady state and time of flight),” in Optical-thermal response of laser-irradiated tissue, A. J. Welch and M. J. C. van Gemert, eds. (Springer, 2011).

55. D. D. Duncan and S. J. Kirkpatrick, “Can laser speckle flowmetry be made a quantitative tool?” J. Opt. Soc. Am. A 25(8), 2088–2094 (2008). [CrossRef]

56. A. Middleton and D. Fisher, “Discrete scatterers and autocorrelations of multiply scattered light,” Phys. Rev. B 43(7), 5934–5938 (1991). [CrossRef]

57. J. Li, L. Qiu, C. S. Poon, and U. Sunar, “Analytical models for time-domain diffuse correlation spectroscopy for multi-layer and heterogeneous turbid media,” Biomed. Opt. Express 8(12), 5518–5532 (2017). [CrossRef]

58. L. O. Reynolds and N. J. McCormick, “Approximate two parameter phase function for light Scattering,” J. Opt. Soc. Am. 70(10), 1206–1212 (1980). [CrossRef]

59. M. Hammer, D. Schweitzer, B. Michel, E. Thamm, and A. Kolb, “Single scattering by red blood cells,” Appl. Opt. 37(31), 7410–7418 (1998). [CrossRef]

60. A. N. Yaroslavsky, I. V. Yaroslavsky, T. Goldbach, and H. Schwarzmaier, “Influence of the scattering phase function approximation on the optical properties of blood determined from the integrating sphere measurements,” J. Biomed. Opt. 4(1), 47–54 (1999). [CrossRef]

61. M. Hammer, A. N. Yaroslavsky, and D. Schweitzer, “A scattering phase function for blood with physiological haematocrit,” Phys. Med. Biol. 46(3), N65–N69 (2001). [CrossRef]

Case number	Description	Condition	Expression
1	Brownian motion dominates random flow (A1)	$τ_{d} ≪ 6 D_{B} / v^{2}$	$g_{1}^{B N} (τ_{s}, τ_{d}) = \exp [\bar{m} (τ_{s}) {g_{1}^{s s} (τ_{d}) - 1}],$ where $g_{1}^{s s} (τ_{d}) = {⟨ \exp [- D_{B} τ_{d} q^{2} (θ)] ⟩}_{θ}$
2	first cumulant term dominates second cumulant term (A2)	$⟨ △ r^{2} (τ_{d}) ⟩ ≪ \frac{6}{\bar{m} (τ_{s}) R ⟨ q^{2} (θ) ⟩}$	$g_{1}^{B N} (τ_{s}, τ_{d}) = \exp [\begin{array}{l} - \frac{1}{6} \bar{m} (τ_{s}) ⟨ q^{2} (θ) ⟩ \\ \times (v^{2} τ_{d}^{2} + 6 D_{B} τ_{d}) \end{array}]$
3	first order τ_d term dominates second order τ_d term (A1&2)	$τ_{d} ≪ \frac{D_{B}}{\| \frac{v^{2}}{6} - \bar{m} (τ_{s}) R ⟨ q^{2} (θ) ⟩ D_{B}^{2} \|}$	$\begin{array}{ll} g_{1}^{B N} (τ_{s}, τ_{d}) & = \exp [- \bar{m} (τ_{s}) ⟨ q^{2} (θ) ⟩ D_{B} τ_{d}] \\ = g_{1}^{D W S} (τ_{s}, τ_{d}) \end{array}$

Constants	Description	Value
D_B (m².s⁻¹)	Brownian diffusion coefficient	6 × 10⁻¹¹
p_dyn	probability of dynamic scattering	10%
λ (nm)	wavelength in vacuum	850
n	refractive index of medium	1.4
k (m⁻¹)	wave number in medium	1.035 × 10⁷
v* (mm.s⁻¹)	selected root-mean-squared speed	1.47	1.36	1.30	1.28	1.25
g_dyn	anisotropy factor of dynamic scatterers	0.60	0.80	0.90	0.94	0.98
µ_s (cm⁻¹)	scattering coefficient of medium	100
µ_a (cm⁻¹)	absorption coefficient of medium	0.01
g_stat	anisotropy factor of static medium	0.9
2δτ_s (ps)	time-of-flight (TOF) bin width	5
Y_dyn bin	total dimensionless momentum transfer bin	0.01
N	number of photons launched	300,000,000

Case number	Description	Condition	Expression
1	Brownian motion dominates random flow (A1)	$τ_{d} ≪ 6 D_{B} / v^{2}$	$g_{1}^{B N} (τ_{s}, τ_{d}) = \exp [\bar{m} (τ_{s}) {g_{1}^{s s} (τ_{d}) - 1}],$ where $g_{1}^{s s} (τ_{d}) = {⟨ \exp [- D_{B} τ_{d} q^{2} (θ)] ⟩}_{θ}$
2	first cumulant term dominates second cumulant term (A2)	$⟨ △ r^{2} (τ_{d}) ⟩ ≪ \frac{6}{\bar{m} (τ_{s}) R ⟨ q^{2} (θ) ⟩}$	$g_{1}^{B N} (τ_{s}, τ_{d}) = \exp [\begin{array}{l} - \frac{1}{6} \bar{m} (τ_{s}) ⟨ q^{2} (θ) ⟩ \\ \times (v^{2} τ_{d}^{2} + 6 D_{B} τ_{d}) \end{array}]$
3	first order τ_d term dominates second order τ_d term (A1&2)	$τ_{d} ≪ \frac{D_{B}}{\| \frac{v^{2}}{6} - \bar{m} (τ_{s}) R ⟨ q^{2} (θ) ⟩ D_{B}^{2} \|}$	$\begin{array}{ll} g_{1}^{B N} (τ_{s}, τ_{d}) & = \exp [- \bar{m} (τ_{s}) ⟨ q^{2} (θ) ⟩ D_{B} τ_{d}] \\ = g_{1}^{D W S} (τ_{s}, τ_{d}) \end{array}$

Constants	Description	Value
D_B (m².s⁻¹)	Brownian diffusion coefficient	6 × 10⁻¹¹
p_dyn	probability of dynamic scattering	10%
λ (nm)	wavelength in vacuum	850
n	refractive index of medium	1.4
k (m⁻¹)	wave number in medium	1.035 × 10⁷
v* (mm.s⁻¹)	selected root-mean-squared speed	1.47	1.36	1.30	1.28	1.25
g_dyn	anisotropy factor of dynamic scatterers	0.60	0.80	0.90	0.94	0.98
µ_s (cm⁻¹)	scattering coefficient of medium	100
µ_a (cm⁻¹)	absorption coefficient of medium	0.01
g_stat	anisotropy factor of static medium	0.9
2δτ_s (ps)	time-of-flight (TOF) bin width	5
Y_dyn bin	total dimensionless momentum transfer bin	0.01
N	number of photons launched	300,000,000

Beyond diffuse correlations: deciphering random flow in time-of-flight resolved light dynamics

Abstract

1. Introduction

2. Field autocorrelation models

2.1 Nature of scatterer motion

2.2 Adapting Bonner and Nossal’s theory

2.3 Critical examination of the cumulant approximation

2.3.1 Special case 1: Brownian motion dominates random flow (A1)

2.3.2 Special case 2: first cumulant dominates second cumulant (A2)

2.3.3 Special case 3: first order τ_d term dominates (A1&2)

2.3.4 First cumulant for random flow cancels second cumulant for Brownian diffusion

2.4 Monte Carlo simulations

3. Results

3.1 First cumulant approximation (A2)

3.2 Neglecting random flow (A1)

3.3 Conventional DWS (A1&2)

3.4 Full cumulant expansion (Bonner and Nossal) with hybrid motion

4. Discussion

4.1 Accuracy of the cumulant approximation

4.2 Revealing random flow

4.3 Dynamic scattering probability dependence

4.4 Proposed improved expression

4.5 Gegenbauer phase function

4.6 Simplifying assumptions

5. Conclusion

Appendix A

Appendix B

Funding

Acknowledgments

Disclosures

References

Cited By

Figures (15)

Tables (2)

Equations (32)

Optics Express

Abstract

1. Introduction

2. Field autocorrelation models

2.1 Nature of scatterer motion

2.2 Adapting Bonner and Nossal’s theory

2.3 Critical examination of the cumulant approximation

2.3.1 Special case 1: Brownian motion dominates random flow (A1)

2.3.2 Special case 2: first cumulant dominates second cumulant (A2)

2.3.3 Special case 3: first order τd term dominates (A1&2)

2.3.4 First cumulant for random flow cancels second cumulant for Brownian diffusion

2.4 Monte Carlo simulations

3. Results

3.1 First cumulant approximation (A2)

3.2 Neglecting random flow (A1)

3.3 Conventional DWS (A1&2)

3.4 Full cumulant expansion (Bonner and Nossal) with hybrid motion

4. Discussion

4.1 Accuracy of the cumulant approximation

4.2 Revealing random flow

4.3 Dynamic scattering probability dependence

4.4 Proposed improved expression

4.5 Gegenbauer phase function

4.6 Simplifying assumptions

5. Conclusion

Appendix A

Appendix B

Funding

Acknowledgments

Disclosures

References

Cited By

Figures (15)

Tables (2)

Equations (32)

Optics Express

2.3.3 Special case 3: first order τ_d term dominates (A1&2)