Abstract
Blood flow index (BFI) is an optically accessible parameter, with unit distance-squared-over-time, that is widely used as a proxy for tissue perfusion. BFI is defined as the dynamic scattering probability (i.e. the ratio of dynamic to overall reduced scattering coefficients) times an effective Brownian diffusion coefficient that describes red blood cell (RBC) motion. Here, using a wavelength division multiplexed, time-of-flight- (TOF) - resolved iNIRS system, we obtain TOF-resolved field autocorrelations at 773 nm and 855 nm via the same source and collector. We measure the human forearm, comprising biological tissues with mixed static and dynamic scattering, as well as a purely dynamic scattering phantom. Our primary finding is that forearm BFI increases from 773 nm to 855 nm, though the magnitude of this increase varies across subjects (23% ± 19% for N = 3). However, BFI is wavelength-independent in the purely dynamic scattering phantom. From these data, we infer that the wavelength-dependence of BFI arises from the wavelength-dependence of the dynamic scattering probability. This inference is further supported by RBC scattering literature. Our secondary finding is that the higher-order cumulant terms of the mean squared displacement (MSD) of RBCs are significant, but decrease with wavelength. Thus, laser speckle and related modalities should exercise caution when interpreting field autocorrelations.
© 2024 Optica Publishing Group under the terms of the Optica Open Access Publishing Agreement
1. Introduction
Diffuse Correlation Spectroscopy (DCS) [1–3] is an established optical technique for non-invasive estimation of biological tissue blood flow index (BFI) [1]. BFI is defined as the product of α, the ratio of the dynamic reduced scattering coefficient ($\mu_{s, d y n}^{\prime}$) to the total tissue reduced scattering coefficient ($\mu_{s}^{\prime}$), and DB, an effective red blood cell (RBC) Brownian diffusion coefficient. BFI is a proxy for nutritive blood flow, and correlates with conventional perfusion metrics [4–13]. BFI is typically measured in the wavelength range from 767-855 nm [14–24], though there is a growing push in the DCS field to measure around 1064 nm, a wavelength range which confers several benefits [25–28]. To the extent that BFI is a valid proxy for tissue perfusion, one expects that it should not depend on the wavelength of light used to probe tissue. Indeed, most studies in the field have treated BFI as wavelength-independent, and at least one study explicitly states this assumption [23]. However, a wavelength-dependent α could introduce inconsistencies in BFI measurements [23]. Directly testing wavelength dependence of α experimentally is challenging as there are numerous confounding factors. For instance, in multi-wavelength DCS, different sources or collectors for the two wavelengths may introduce variability in the tissue volume probed. Also, DCS autocorrelations depend not only on BFI, but also the time-of-flight distribution (DTOF), which is also wavelength-dependent. Thus, differences in BFI can be masked by other experimental factors in multi-wavelength DCS [23].
Interferometric near-infrared spectroscopy (iNIRS) was recently introduced to provide autocorrelations from biological tissue with TOF resolution [21,29–33]. Here, to address the above questions, we further advance the iNIRS approach with wavelength-division multiplexing [33]. The two system wavelengths, 773 and 855 nm, span the most frequently-used DCS wavelengths. The dual-wavelength iNIRS setup was meticulously designed to encompass key pivotal features: shared source and collector for both wavelengths, simultaneous autocorrelation measurements at both wavelengths without time multiplexing, and TOF-resolved autocorrelations with matched TOF resolution at both wavelengths. These elements contribute to a robust framework for recovering BFI experimentally, avoiding many issues that may confound BFI comparisons in DCS.
Our primary finding is a wavelength-dependence of α in biological tissue. This finding, while unexpected for the DCS community, is in fact already implied by the RBC light scattering literature [34–36]. This result should inform interpretation of DCS in its numerous applications. Our secondary finding; that that higher-order cumulant terms are significant but decrease with wavelength, should inform interpretation of laser speckle and related modalities.
2. Methods
2.1 Theory
iNIRS recovers a TOF-resolved field autocorrelation, g1(τs,τd), where τs is TOF and τd is time lag [21,29,31,32]. Allowing for a finite TOF resolution, the iNIRS experimental data closely relates to DWS theory [37], which is based on simple fundamental principles. According to DWS, the field autocorrelation, $g_1^{DWS}$(τs,τd), is given by
Though well-known and widely used, Eq. (2) includes a number of assumptions: approximation of a distribution of the number dynamic scattering events by a single mean value, neglect of higher order MSD cumulants, and neglect of random flow. Therefore, Eq. (2) may not describe g1(τs,τd) well at early τs and late τd [38], as is required for this work. On the other hand, Bonner and Nossal’s seminal dynamic light scattering theory [39], if extended to include TOF [38], incorporates a Poisson-distributed number of dynamic scattering events and implicitly includes higher MSD cumulants by integrating over the dynamic phase function [38]. Their theory can be further modified to allow both Brownian motion and random flow, a feature missing from Eq. (2), to drive the mean squared displacement (MSD) of dynamic scatterers [38]. The error arising from exclusion of the second cumulant of MSD increases with higher dynamic scattering anisotropy, gdyn, and fewer dynamic scattering events [38]. Since RBCs have gdyn∼0.98 and we need to analyze g1(τs,τd) at early τs, we use an expression from [38] which extends Bonner and Nossal’s theory and incorporates fewer assumptions than Eq. (2).
In Eq. (3), $\bar{m}$(τs) is the average number of dynamic scattering events at τs, q2(θ) is the momentum transfer as a function of dynamic scattering angle θ, <Δr2(τd)〉 is the MSD of the RBCs, and angle brackets $\langle\cdot\rangle_{\theta}$ denotes an angular weighted average over θ. In contrast to Eq. (2), the more general Eq. (3) can describe g1(τs,τd) even with few dynamic scattering events and high gdyn [38]. Hereafter, we will drop the subscript θ in the angular average for sake of readability.
In simulation, Eq. (3) was shown to better describe TOF-resolved autocorrelations, particularly in the slowly-decaying “tails” where DWS theory of Eq. (2) fails [38]. Taylor series expansion of the exponential term in Eq. (3) up to three MSD cumulants yields
Note that the fourth and the higher order terms are contained in O[〈Δr2(τd)〉4]. Incorporating hybrid motion [40,41], 〈Δr2(τd)〉=6DBτd+v2τd2, where v2 is the second moment of the Gaussian velocity distribution and arranging different powers of τd (see full details in Appendix A), we have
O(τd4) contains fourth and higher orders of τd. It is also worthwhile noting that MSD-cumulants [Eq. (4)] do not necessarily correspond with τd-cumulants [Eq. (8)] when hybrid motion is assumed. Importantly, in Eq. (8), different order terms for Brownian motion and random flow contribute to both the second- and third-order τd terms, with opposite signs. Moreover, setting v to 0 and neglecting second- and higher orders of τd in Eq. (8) yields conventional DWS [Eq. (2)], given $\bar{m}$(τs)=$\mu_{s, d y n}$vcτs, where $\mu_{s, d y n}$ is the dynamic scattering coefficient.
2.2 Experimental setup
Dual-wavelength iNIRS [Fig. 1] comprises two distributed feedback (DFB) lasers with center wavelengths of 773 nm (Eagleyard EYP-DFB-0773-00075-1500-TOC03-0002) and 855 nm (Eagleyard EYP-DFB-0855-00150-1500-TOC03-0000). We achieved synchronous wavelength tuning by sinusoidally modulating the currents via two independent current controllers (SRS LDC-501), with modulation input voltages from a common function generator (TTi TGF4042). Wavelength tuning with 2.5 Vpp from the function generator provided a TOF resolution of 75 ps FWHM for 773 nm and 67 ps FWHM for 855 nm. To equalize the TOF resolution, we introduced a series resistor to create a voltage divider at the input of the current controller driving the 855 nm laser. This enabled us to achieve a TOF resolution of 58 ps for both wavelengths with 3.3 Vpp from the function generator. Both configurations were used in this study. The sinusoidal drive frequency was 100 kHz, yielding 5 μs field autocorrelation resolution [30].
The lasers emitted light with elliptical spatial beam profiles. To couple this light into a single-mode (SM) fiber, for each wavelength, we collimated each beam using a convex aspheric lens and adjusted the beam shape to be circular with an anamorphic prism pair. Subsequently, for each wavelength, the collimated beam passed through a free space isolator before being focused into a SM fiber-based 90:10 splitter, with 90% of the light directed to the sample arm and the remaining 10% to the reference arm.
The two 90% ports were connected to a custom wavelength-division multiplexer (WDM) characterized by 98% transmission and ∼20 dB isolation. The WDM combined both wavelengths before illuminating the sample. The remitted light was collected through a similar WDM, this time to separate the two wavelengths. This approach ensured that both wavelengths illuminated and collected light from the same spot on the tissue, in accordance with our experimental requirements. Slight differences in the probed tissue regions due to variations in optical properties at the two wavelengths were unavoidable, however.
Subsequently, the light from each arm of the splitter WDM was recombined with the appropriate reference arm light of the same wavelength via a 50:50 combiner. Finally, both arms of both 50:50 combiners were connected to New Focus 1807 dual-balanced detectors (DBD) which provided the interference signal. At this stage, any cross-talk in the WDMs was further suppressed by interferometric gating whereby only the sample light of the correct wavelength was amplified with reference light. The DBD output voltages were sampled by an Alazar Technologies 9416 fast data acquisition card (DAQ) at 100 MHz.
The DAQ card derived its clock from the reference clock output of the function generator. All measurements presented in this paper were continuously acquired following a trigger signal, generated from the function generator, and sent to the DAQ to commence data acquisition. The trigger signal to the DAQ is a square wave, synchronous to the sinusoidal signal used to tune the DFB lasers. This mode of triggered and synchronized acquisition, coupled with low timing jitter of the function generator, allowed us to omit recording of the trigger signal, thereby simplifying our data processing and conserving hard drive space.
2.3 Experimental measurements
With the dual wavelength iNIRS system, measurements were performed on an Intralipid mixture, prepared by mixing 21 mL of 20% Intralipid in 500 mL water. No absorbing agent was added. Measurements were performed at room temperature. The source-detector separation was 11 mm for the Intralipid measurements.
Measurements were also performed on the forearms of three healthy males at rest. The procedures and protocols were approved by the NYU Langone and UC Davis Institutional Review Boards. The source-detector separation was 6 mm for the human measurements. The probe was placed on the dominant forearm right above the brachioradialis muscle. The combined skin and adipose tissue thickness over the brachioradialis muscle was measured by a skinfold caliper for each subject. Measurements at both wavelengths were acquired simultaneously for 40 seconds. Though major findings were ultimately verified across subjects (age 31.67 yrs. ± 4.04 yrs. for N = 3), we mainly present measurements and analysis from one subject (age: 28 yrs).
2.4 Data acquisition and processing
iNIRS provides time-of-flight (TOF) resolution through wavelength tuning of the source [21,29–33]. To ensure accurate tracking of instantaneous wavelength variations, a calibration measurement, hereafter referred to as the instrument response function (IRF) measurement, was conducted. During the IRF measurement, the sample in Fig. 1 is replaced with a single-mode (SM) patch cord of known length. This process generates a fringe pattern at the output of the dual-balanced detectors (DBD) for each complete wavelength sweep.
We initiate the IRF processing by identifying stationary points, which correspond to time instances where the instantaneous wavelength reaches either the maximum or minimum of the sweep. The phase information derived from the Hilbert transform of the recorded fringes serves as a proxy for the optical frequency. A stationary point occurs where the time derivative of the phase becomes zero. It is important to note that for accurate wavelength tracking using phase information, the time derivative of the phase must change sign at every alternate sweep. We designate a sweep as a “forward sweep” when the optical frequency increases from minimum to maximum and a “backward sweep” when it decreases from maximum to minimum. Since our DFB lasers exhibited non-identical forward and backward sweeps, separate phase-versus-time calibration curves were generated for each sweep type by averaging the phase over many sweeps. The fringe envelope for each sweep type was also determined from the IRF measurement, and later used for Gaussian shaping of the measurement data.
With the envelope and phase-versus-time calibration curves for the sweep types in hand, we are prepared for measurements using the setup shown in Fig. 1. As previously mentioned, triggered acquisition ensures that the first stationary point always occurs a fixed time after acquisition starts, and always for the same sweep transition (e.g. forward-to-backward sweep). A low phase noise function generator ensures subsequent stationary points occur after each sweep cycle. Once stationary points are identified, the time axis of each sweep is transformed into phase (or by proxy, optical frequency) using the corresponding phase-versus-time curve from the IRF measurement.
Thereafter, we combined the measured data from the two sweep types, followed by resampling, Gaussian shaping and finally, Discrete Fourier Transform, to obtain the complex mutual coherence function Γrs(τs,td) as described in Ref. [30]. Here td refers to delay time. Finally, the un-normalized iNIRS field autocorrelation function, $G_1\left(\tau_s, \tau_d\right)$, is obtained from Eq. (9) [29,30]:
Background measurements are acquired in the same way as sample interference measurements, except while blocking light transmission through the sample arm. Data processing steps for the background measurement mirror those of sample interference measurements. Subtracting the background autocorrelation $G_1^{B G}\left(\tau_s, \tau_d=0\right)$ from the measurement $G_1\left(\tau_s, \tau_d=0\right)$ yields the desired temporal point spread function (TPSF), used to recover optical properties of the probed region.
2.5 Scaling approaches for wavelength-independent DCS measurements
This work investigates both Intralipid and biological tissue. The key difference is that paths in Intralipid involve only dynamic scattering [Fig. 2(a)] with gdyn∼0.6, while paths in biological tissue involve scattering from static (blue) and dynamic (red) scatterers [Fig. 2(b)], the latter with gdyn∼0.98. For both media, the dual wavelength TOF-resolved unnormalized field autocorrelation functions $G_1\left(\tau_s, \tau_d\right)$ were first fitted to the single exponential DWS model of Eq. (2) to obtain the TOF-resolved decay constants, ξ(τs),
Assuming negligible wavelength-dependence of vc, wavelength-dependence of ξ(τs) should arise from k2 = (2πn/λ)2 and potentially, $\mu_{s, d y n}^{\prime}$. We assume that n, the phase refractive index, does not vary significantly with wavelength. Therefore, we may predict ${\xi ^{{\lambda _1}}}({{\tau_s}} )$ at wavelength λ1 from ${\xi ^{{\lambda _2}}}({{\tau_s}} )$ at wavelength λ2 via
In Eq. (11), $\mu_{s, d y n}^{\prime}$(λ1) and $\mu_{s, d y n}^{\prime}$(λ2) are reduced dynamic scattering coefficients at wavelengths λ1 and λ2, respectively. Similarly, to predict $g_1^{\lambda_2}$(τs,τd) from $g_1^{\lambda_1}$(τs,τd), the τd-axis of $g_1^{\lambda_1}$(τs,τd) is scaled according to Eq. (12),
However, to the best of our knowledge, only $\mu_{s}^{\prime}$, and not $\mu_{s, d y n}^{\prime}$, can be directly measured in-vivo. To address this, we adapt two mutually exclusive scaling approaches to infer wavelength-dependence of $\mu_{s, d y n}^{\prime}$:
- 1. constant α: This scaling approach assumes $\mu_{s, d y n}^{\prime}$(λ)=α$\mu_{s}^{\prime}$(λ). This assumption is valid per force in purely dynamic phantoms [Fig. 2(a)] such as Intralipid (where $\mu_{s, d y n}^{\prime}$=$\mu_{s}^{\prime}$, and α=1 independent of wavelength). This is the standard assumption for DCS of biological tissues, even though α≠1.
- 2. constant $\mu_{s, d y n}^{\prime}$: This scaling approach assumes $\mu_{s, d y n}^{\prime}$ is constant across wavelengths. A corollary of this assumption is that α(λ)= $\mu_{s, d y n}^{\prime}$ /$\mu_{s}^{\prime}$(λ) is wavelength-dependent in general. To our knowledge, assuming a wavelength-dependent α is very rare.
DWS-based scaling approaches consider early τd. To account for the tails of the autocorrelation at late τd, we also fitted measurements to a more accurate expression, inspired by the higher-order τd terms seen in Eq. (8):
Coefficients of τd, τd2 and τd3 are represented as ξ1(τs), ξ2(τs) and ξ3(τs), respectively. To provide a simple framework to interpret the coefficients in Eq. (13), we take v = 0 in Eq. (8) and neglect O(τd4):
Following the constant $\mu_{s, d y n}^{\prime}$-scaling approach and using Eqs. (5), (6), (7) with Eq. (14), the coefficients ξ1(τs), ξ2(τs) and ξ3(τs) must be multiplied by λ2, λ4 and λ6 respectively to remove their wavelength-dependency, as the ratios of the quantities 〈sin2(θ/2)〉, 〈sin4(θ/2)〉 and 〈sin6(θ/2)〉 at the two wavelengths are found to be nearly the same for a Henyey-Greenstein phase function. Thus, for inter-wavelength transformations, we used the following relations:
Here, $\xi _1^{{\lambda _1}}({{\tau_s}} )$ and $\xi _1^{{\lambda _2}}({{\tau_s}} )$ correspond to the coefficients of τd at wavelengths λ1 and λ2 respectively, $\xi _2^{{\lambda _1}}({{\tau_s}} )$ and $\xi _2^{{\lambda _2}}({{\tau_s}} )$ correspond to the coefficients of τd2 at wavelengths λ1 and λ2 respectively, and $\xi _3^{{\lambda _1}}({{\tau_s}} )$ and $\xi _3^{{\lambda _2}}({{\tau_s}} )$ correspond to the coefficients of $\tau _d^{3}$ at wavelengths λ1 and λ2 respectively. Scaling of τd-cumulants described by Eqs. (15), (16) and (17) was compared against a simpler approach where all cumulants were scaled by the ratio of λ2.
2.6 Data analysis
In order to compare the scaling approaches described in the previous section, we need to quantify the degree of agreement between the predicted $g_1^{\lambda_2}$(τs,τd), obtained from scaling the the τd -axis of $g_1^{\lambda_1}$(τs,τd), and the measured $g_1^{\lambda_2}$(τs,τd), obtained from our dual wavelength iNIRS setup. To facilitate this comparison, we utilize a Bland-Altman plot [42] to visually evaluate agreement between predicted and measured values. This method involves plotting the mean of the predicted and measured $g_1^{\lambda_2}$(τs,τd) values on the X-axis and their difference on the Y-axis for all τd at a fixed τs. In this analysis, a smaller difference in the Bland-Altman plot signifies a better agreement between the predicted and measured values. Here, the Bland-Altman plots offer insights into the degree of agreement for every τd. Note that the difference at τd = 0 will invariably be zero due to normalization.
The Bland-Altman plots offer detailed analysis at individual τs. To comprehensively assess agreement across all TOFs (τs), we employed the coefficient of determination (R2) [43]. R2 quantifies the goodness of agreement between the predicted and measured $g_1^{\lambda_2}$(τs,τd). We used 1-R2 to compare scaling approaches, where a lower value signifies better agreement between predicted and measured quantities.
3. Results
We begin by fitting G1(τs,τd) from Intralipid (gdyn∼0.6) to the simple DWS exponential model [unnormalized form of Eq. (2)], routinely used in the DCS field. For Intralipid data, we fit the entire autocorrelation curve i.e., g1(τs,τd) ranging from 1 to 0. TOF-resolved autocorrelations decay exponentially for both wavelengths [Fig. 3, Fig. 4], consistent with DWS theory. Here the absolute medium “EARLY”, “MID” and “LATE” τs are 60, 330 and 500 ps, respectively. This definition of time gates is consistent with previous DCS literature [28,44].
Unlike the Intralipid phantom, G1(τs,τd) from in-vivo measurements (gdyn∼0.98) do not decay exponentially. Thus, G1(τs,τd) from both wavelengths were fitted to the unnormalized form of Eq. (13) [Fig. 5, Fig. 6]. For in-vivo data, we fit g1(τs,τd) ranging from 1 to 0.1. Since a single exponential model does not fit in-vivo G1(τs,τd), we will not use this model for fitting in-vivo measurements further. Instead we will use the autocorrelation functions directly to test both proposed scaling approaches (Section 2.5) on the Intralipid and in-vivo data.
For Intralipid, we find that constant α-scaling performed best for predicting $g_1^{773}\left(\tau_s, \tau_d\right)$ from $g_1^{855}\left(\tau_s, \tau_d\right)$ [Figs. 7–9]. To scale the τd-axis [Eq. (12)], $\mu_{s}^{\prime}$ were obtained from fitting the measured TPSF, G1(τs,τd = 0), to a diffusion model, yielding 0.84 mm-1 and 0.76 mm-1 for 773 and 855 nm, respectively. Note that where specified the τd-axis of 855 nm is scaled.
For the in-vivo forearm measurements, constant $\mu_{s, d y n}^{\prime}$-scaling performed best in predicting $g_1^{773}\left(\tau_s, \tau_d\right)$ from $g_1^{855}\left(\tau_s, \tau_d\right)$ [Figs. 10–12]. The τd -axis of 855 nm was scaled using Eq. (12) with recovered $\mu_{s}^{\prime}$ of 1.23 and 1.44 mm-1 for 855 and 773 nm respectively.
Finally, to comprehensively compare the scaling techniques, we calculated 1-R2 (where R2 is the coefficient of determination) across different τs to show how well they predicted $g_1^{773}\left(\tau_s, \tau_d\right)$ from $g_1^{855}\left(\tau_s, \tau_d\right)$. For the Intralipid mixture, we find that the constant α-scaling approach is better than constant $\mu_{s, d y n}^{\prime}$ -scaling for every τs [Fig. 13(a)]. In the human forearm in-vivo [Fig. 13(b)], we find that constant $\mu_{s, d y n}^{\prime}$ -scaling is better.
The preceding analysis scaled the τd-axis of the autocorrelations at one wavelength to predict the other wavelength. However, theory [Eqs. (14)-(17)] predicts a more complex relationship where higher-order cumulants scale more strongly with wavelength. So, we proceeded to consider scaling of individual terms in the τd-cumulant expansion. As the first τd-cumulant is sufficient for Intralipid, the results are very straightforward. Namely, the constant α-scaling approach applied to Eq. (11) successfully predicted ξ(τs) at 773 nm from ξ(τs) at 855 nm [Fig. 14]. These results on the exponential fit parameter are unsurprising given the results already presented, which showed that an exponential describes Intralipid G1 at all TOFs, and that the constant α-scaling approach works best for Intralipid.
Next, we proceeded to consider scaling of terms in the τd-cumulant expansion for in vivo measurements. The in-vivo G1(τs,τd) were fitted to Eq. (13) for the coefficients ξ1(τs), ξ2(τs) and ξ3(τs) [Fig. 15(a), Fig. 16(a) and Fig. 17(a), respectively], for both wavelengths. The prediction of these three coefficients at 773 nm from their counterparts at 855 nm was accomplished through scaling based on Eqs. (15), (16) and (17), respectively. Predicted coefficients [Fig. 15(b), Fig. 16(b) and Fig. 17(b)] approximate experimental coefficients reasonably well.
Next, we investigated the ability to predict 773 nm autocorrelations from the τd-cumulant expansion order at 855 nm. Constant $\mu_{s, d y n}^{\prime}$ -scaling [Eqs. (15)-(17), also referred here as theory-based scaling] and a simpler approach where all the constants were scaled by the ratio of the wavelengths squared, here (855/773)2, referred here as λ2-scaling, were investigated. For both approaches, collectively referred to as higher cumulant-scaling, scaled coefficients were plugged into Eq. (13) to generate predictions. A direct Eq. (13) fit of $g_1^{773}\left(\tau_s, \tau_d\right)$ serves as a benchmark. We found that theory-based scaling of the three τd-cumulant coefficients at 855 nm predicted $g_1^{773}\left(\tau_s, \tau_d\right)$ best [Fig. 18, Fig. 19, and Fig. 20].
Finally, we summarized the error (1-R2) of the $g_1^{773}\left(\tau_s, \tau_d\right)$ prediction from $g_1^{855}\left(\tau_s, \tau_d\right)$ using the higher cumulant-scaling methods across τs [Fig. 21]. The 1-R2 values obtained from Eq. (13) fit of $g_1^{773}\left(\tau_s, \tau_d\right)$ serves as a benchmark. We see that the approach where the τdn coefficient is scaled by λ2n outperforms the approach where all coefficients are scaled by λ2. This suggests that higher-order τd terms scale more strongly with wavelength than the first-order (exponential decay) term.
4. Discussion
This study employs a new dual-wavelength iNIRS setup to compare TOF-resolved autocorrelations across wavelengths for the first time. Our investigation uncovers a previously-overlooked wavelength-dependence of α, the probability of dynamic scattering. The immediate consequence of this primary finding is that BFI, the optical index of blood flow, is wavelength-dependent. Secondary findings highlight the importance of incorporating higher MSD cumulants in modeling the later time lags of field autocorrelations (“tails”), particularly at early and even moderate TOFs.
4.1 Wavelength dependence of α
We applied two scaling methods to estimate 773 nm autocorrelations from 855 nm autocorrelations: the constant α-scaling approach and the constant $\mu_{s, d y n}^{\prime}$-scaling approach.
We observed that the constant α-scaling approach accurately predicted 773 nm autocorrelations from 855 nm autocorrelations in an Intralipid mixture [Fig. 7, Fig. 8, Fig. 9, and Fig. 13(a)], consistent with α=1. As expected, the constant $\mu_{s, d y n}^{\prime}$approach performed worse in this purely dynamic scattering phantom, as Intralipid reduced scattering does vary with wavelength [45]. Yet in-vivo, the constant $\mu_{s, d y n}^{\prime}$ approach better predicted 773 nm autocorrelations from 855 nm autocorrelations [Fig. 10, Fig. 11, Fig. 12, and Fig. 13(b)]. Though surprising, this result in fact aligns with RBC measurements in the literature [34–36], where reduced scattering is similar at 773 nm and 855 nm. To further underscore this point, α was derived [Fig. 22] from prior RBC scattering literature [34], under various assumptions for wavelength-dependent tissue scattering [46]. Though derived purely from the literature, Fig. 22 recapitulates the major experimental result of our work: the wavelength-dependence of α. This result implies that BFI=αDB is also wavelength-dependent (as clearly, DB is wavelength-independent).
BFI was derived from the autocorrelation decay rate at τs = 400 ps [31]. Consistent with the above discussion, in-vivo forearm BFI ratios between wavelengths differed from 1, and also varied between subjects [Table 1]. This variation may occur because different skin types yield different wavelength dependencies of forearm reduced scattering [47,48]; yet dynamic (RBC) scattering should not depend on skin type. The mean and standard deviation of BFI for the in-vivo measurements were 1.86 × 10−7 (±3.14 × 10−8) mm2/sec for 855 nm and 1.57 × 10−7 (±4.72 × 10−8) mm2/sec for 773 nm, across three subjects who participated in the study.
A wavelength-dependent α implies that static and dynamic scattering need not co-vary with wavelength. Fortuitously, the constant $\mu_{s, d y n}^{\prime}$ approach worked well for the two chosen wavelengths in our in-vivo forearm experiments. Yet, from a physical standpoint, as refractive index variations leading to scattering are different for static and dynamic tissue, we can envision situations (such as different tissues or wavelengths) where neither α nor $\mu_{s, d y n}^{\prime}$ are invariant with wavelength. To enable meaningful comparisons, one solution is to evaluate $\mu_{s}^{\prime}$BFI=$\mu_{s, d y n}^{\prime}$DB, as a substitute for BFI alone. Such an approach mitigates concerns arising from differential behavior of $\mu_{s, d y n}^{\prime}$ and $\mu_{s}^{\prime}$ with wavelengths, but is still subject to variability in $\mu_{s, d y n}^{\prime}$ with wavelength, if any.
4.2 Better autocorrelation modelling by including higher MSD cumulants
Building on theoretical predictions [38], in this work we experimentally demonstrated the importance of considering higher-order MSD cumulants in autocorrelation modeling. We observed that the common single exponential model of DWS works well for Intralipid (gdyn∼0.6 [49–51]), but inadequately captures the slowly decaying TOF-resolved field autocorrelation tails in vivo (gdyn∼0.98 for RBCs [52–54]). This limitation is particularly pronounced at early TOF and is visible even at moderate TOFs of 330 ps [Fig. 5 and Fig. 6]. Our findings underscore the need for considering higher MSD cumulants to describe autocorrelation tails in laser speckle, when determining camera exposure times, and DCS, especially with short source-collector separations, when selecting fitting ranges for autocorrelation functions. Interestingly, higher-order MSD cumulants are almost never considered in the field. Admittedly, their impact would be lessened in g2 = 1+|g1|2 (intensity autocorrelation) as compared to g1 (field autocorrelation). Addressing these considerations is crucial for achieving more accurate and comprehensive interpretations of tissue dynamics.
The multiple wavelengths in this study shed further light on the nature of deviations from the single exponential DWS model [Eq. (2)]. As anticipated by theory [Eq. (8)], we found more prominent autocorrelation tails for the shorter wavelengths (773 nm) versus longer wavelengths (855 nm). Interestingly, we observed a fortuitous alignment of the in-vivo autocorrelation tails (late τd) between the two wavelengths in Fig. 10(a), Fig. 11(a) and Fig. 12(a). This alignment appears to result from two cancelling errors: first, the faster initial decay slope at the shorter wavelength and second, the more prominent tail at the shorter wavelength.
Last, to further test the importance of higher-order cumulants, we scaled the 855 nm τd-cumulant coefficients to predict autocorrelation functions at 773 nm (Fig. 18, Fig. 19, Fig. 20 and Fig. 21). Results reinforce the theoretical models [e.g. Eqs. (4), (8), and (14)] that predict stronger scaling of higher cumulants with wavelength. It is noteworthy that the scale factors were derived from Eq. (14) while neglecting random flow; however, random flow may still impact the tails of the autocorrelation function. Likewise, while the signs of the τd-cumulant coefficients [ξ1(τs), ξ2(τs), and ξ3(τs)] align with expectations from a pure Brownian diffusion model [Fig. 15, Fig. 16 and Fig. 17], we cannot conclusively exclude the impact of random flow on the autocorrelations [38]. We also cannot exclude the impact of multiple layers (skin, fat, muscle) with different dynamics which may also affect autocorrelation tails [31].
Finally, it is important to note that the results on the effects of α are based mainly on the first cumulant (DWS theory) which describes the initial decay at early τd of the autocorrelation function. This is distinguished from the higher-order cumulants which affect the autocorrelation “tails” at late τd. Hence the two results are discussed separately and independently.
5. Conclusions
Our experimental findings substantiate an increase in the dynamic scattering probability, α, by 23 ± 19% from 773 nm to 855 nm in the forearms of 3 human subjects. This implies that BFI is wavelength-dependent, and depending on scattering, may be subject dependent as well. Therefore, we advocate for a meticulous approach to interpreting DCS-derived BFI that accounts for a possibly wavelength-dependent α, among numerous other factors that can confound comparisons between DCS systems. Additionally, our experiments fortify prior theoretical predictions underscoring the importance of higher MSD cumulants in accurate autocorrelation modeling at early to intermediate TOFs.
Appendix A
Here we show how we arrive at Eq. (8) from Eq. (4). We begin by taking natural logarithm on both sides of Eq. (4) and multiply them with -1/$\; \bar{m}$(τs) to obtain:
Assuming the hybrid motion, i.e., 〈Δr2(τd)〉 = 6DBτd +v2τd2, we obtain:
Now, rearranging the right-hand side (R.H.S) of Eq. (A.2) in terms of different orders (or powers) of τd, we have
Taking
Now multiplying both the sides with -$\bar{m}$(τs) and then taking exponential on both the sides of Eq. (A.4), we obtain Eq. (8).
Funding
Research to Prevent Blindness (unrestricted grant to Ophthalmology Department); National Institutes of Health (EB029747, EB032840, EY031469).
Acknowledgments
We specially thank Ruoyu Meng and Dr. Santosh Aparanji for aiding in experiments.
Disclosures
The authors declare no conflicts of interest.
Data availability
Data underlying the results presented in this paper are not publicly available at this time but may be obtained from the authors upon reasonable request.
References
1. T. Durduran and A. G. Yodh, “Diffuse correlation spectroscopy for non-invasive, micro-vascular cerebral blood flow measurement,” NeuroImage 85, 51–63 (2014). [CrossRef]
2. 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]
3. J. Selb, K.-C. Wu, J. Sutin, et al., “Prolonged monitoring of cerebral blood flow and autoregulation with diffuse correlation spectroscopy in neurocritical care patients,” Neurophotonics 5(4), 1 (2018). [CrossRef]
4. S. Carp, G. P. Dai, D. A. Boas, et al., “Validation of diffuse correlation spectroscopy measurements of rodent cerebral blood flow with simultaneous arterial spin labeling MRI; towards MRI-optical continuous cerebral metabolic monitoring,” Biomed. Opt. Express 1(2), 553–565 (2010). [CrossRef]
5. T. Durduran, C. Zhou, E. M. Buckley, et al., “Optical measurement of cerebral hemodynamics and oxygen metabolism in neonates with congenital heart defects,” J. Biomed. Opt. 15(5), 037004 (2010). [CrossRef]
6. G. Yu, T. F. Floyd, T. Durduran, et al., “Validation of diffuse correlation spectroscopy for muscle blood flow with concurrent arterial spin labeled perfusion MRI,” Opt. Express 15(3), 1064–1075 (2007). [CrossRef]
7. C. Zhou, S. A. Eucker, T. Durduran, et al., “Diffuse optical monitoring of hemodynamic changes in piglet brain with closed head injury,” J. Biomed. Opt. 14(3), 034015 (2009). [CrossRef]
8. E. M. Buckley, N. M. Cook, T. Durduran, et al., “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]
9. N. Roche-Labarbe, S. A. Carp, A. Surova, et al., “Noninvasive optical measures of CBV, StO2, CBF index, and rCMRO2 in human premature neonates’ brains in the first six weeks of life,” Hum. Brain Mapp. 31(3), 341–352 (2010). [CrossRef]
10. M. N. Kim, T. Durduran, S. Frangos, et al., “Noninvasive measurement of cerebral blood flow and blood oxygenation using near-infrared and diffuse correlation spectroscopies in critically brain-injured adults,” Neurocrit. Care 12(2), 173–180 (2010). [CrossRef]
11. M. Diop, K. Verdecchia, T.-Y. Lee, et al., “Calibration of diffuse correlation spectroscopy with a time-resolved near-infrared technique to yield absolute cerebral blood flow measurements,” Biomed. Opt. Express 2(7), 2068–2081 (2011). [CrossRef]
12. V. Jain, E. M. Buckley, D. J. Licht, et al., “Cerebral oxygen metabolism in neonates with congenital heart disease quantified by MRI and optics,” J. Cereb. Blood Flow Metab. 34(3), 380–388 (2014). [CrossRef]
13. M. Giovannella, B. Andresen, J. B. Andersen, et al., “Validation of diffuse correlation spectroscopy against 15O-water PET for regional cerebral blood flow measurement in neonatal piglets,” J. Cereb. Blood Flow Metab. 40(10), 2055–2065 (2020). [CrossRef]
14. D. R. Busch, R. Balu, W. Guo, et al., “Detection of Brain Hypoxia Based on Noninvasive Optical Monitoring of Cerebral Blood Flow with Diffuse Correlation Spectroscopy,” Neurocrit. Care 30(1), 72–80 (2019). [CrossRef]
15. T. Durduran, C. Zhou, B. L. Edlow, et al., “Transcranial optical monitoring of cerebrovascular hemodynamics in acute stroke patients,” Opt. Express 17(5), 3884–3902 (2009). [CrossRef]
16. D. Wang, A. B. Parthasarathy, W. B. Baker, et al., “Fast blood flow monitoring in deep tissues with real-time software correlators,” Biomed. Opt. Express 7(3), 776–797 (2016). [CrossRef]
17. A. I. Zavriyev, K. Kaya, P. Farzam, et al., “The role of diffuse correlation spectroscopy and frequency-domain near-infrared spectroscopy in monitoring cerebral hemodynamics during hypothermic circulatory arrests,” JTCVS techniques 7, 161–177 (2021). [CrossRef]
18. S. Y. Lee, K. R. Cowdrick, B. Sanders, et al., “Noninvasive optical assessment of resting-state cerebral blood flow in children with sickle cell disease,” Neurophotonics 6(3), 1 (2019). [CrossRef]
19. H. S. Yazdi, T. D. O’Sullivan, A. Leproux, et al., “Mapping breast cancer blood flow index, composition, and metabolism in a human subject using combined diffuse optical spectroscopic imaging and diffuse correlation spectroscopy,” J. Biomed. Opt. 22(4), 045003 (2017). [CrossRef]
20. L. Dong, M. Kudrimoti, D. Irwin, et al., “Diffuse optical measurements of head and neck tumor hemodynamics for early prediction of chemoradiation therapy outcomes,” J. Biomed. Opt. 21(8), 085004 (2016). [CrossRef]
21. D. Borycki, O. Kholiqov, V. J. Srinivasan, et al., “Reflectance-mode interferometric near-infrared spectroscopy quantifies brain absorption, scattering, and blood flow index in vivo,” Opt. Lett. 42(3), 591–594 (2017). [CrossRef]
22. W. Zhou, O. Kholiqov, S. P. Chong, et al., “Highly parallel, interferometric diffusing wave spectroscopy for monitoring cerebral blood flow dynamics,” Optica 5(5), 518–527 (2018). [CrossRef]
23. D. Tamborini, P. Farzam, B. Zimmermann, et al., “Development and characterization of a multidistance and multiwavelength diffuse correlation spectroscopy system,” Neurophotonics 5(1), 1 (2017). [CrossRef]
24. M. Giovannella, D. Contini, M. Pagliazzi, et al., “BabyLux device: a diffuse optical system integrating diffuse correlation spectroscopy and time-resolved near-infrared spectroscopy for the neuromonitoring of the premature newborn brain,” Neurophotonics 6(2), 025007 (2019). [CrossRef]
25. S. Carp, D. Tamborini, D. Mazumder, et al., “Diffuse correlation spectroscopy measurements of blood flow using 1064 nm light,” J. Biomed. Opt. 25(9), 097003 (2020). [CrossRef]
26. C.-S. Poon, D. S. Langri, B. Rinehart, et al., “First-in-clinical application of a time-gated diffuse correlation spectroscopy system at 1064 nm using superconducting nanowire single photon detectors in a neuro intensive care unit,” Biomed. Opt. Express 13(3), 1344–1356 (2022). [CrossRef]
27. N. Ozana, A. I. Zavriyev, D. Mazumder, et al., “Superconducting nanowire single-photon sensing of cerebral blood flow,” Neurophotonics 8, 035006 (2021). [CrossRef]
28. N. Ozana, M. Lue, M. Renna, et al., “Functional time domain diffuse correlation spectroscopy,” Front. Neurosci. 1(16), 932119 (2022). [CrossRef]
29. D. Borycki, O. Kholiqov, S. P. Chong, et al., “Interferometric Near-Infrared Spectroscopy (iNIRS) for determination of optical and dynamical properties of turbid media,” Opt. Express 24(1), 329–354 (2016). [CrossRef]
30. O. Kholiqov, D. Borycki, V. J. Srinivasan, et al., “Interferometric near-infrared spectroscopy (iNIRS): performance tradeoffs and optimization,” Opt. Express 25(23), 28567–28589 (2017). [CrossRef]
31. O. Kholiqov, W. Zhou, T. Zhang, et al., “Time-of-flight resolved light field fluctuations reveal deep human tissue physiology,” Nat. Commun. 11(1), 391 (2020). [CrossRef]
32. O. Kholiqov, W. Zhou, T. Zhang, et al., “Scanning interferometric near-infrared spectroscopy,” Opt. Lett. 47(1), 110–113 (2022). [CrossRef]
33. D. Mazumder, O. Kholiqov, V. J. Srinivasan, et al., “Dual-wavelength interferometric near-infrared spectroscopy (iNIRS),” in Bio-Optics: Design and Application, (Optica Publishing Group, 2023), DW4A. 7.
34. N. Bosschaart, G. J. Edelman, M. C. G. Aalders, et al., “A literature review and novel theoretical approach on the optical properties of whole blood,” Lasers Med. Sci. 29(2), 453–479 (2014). [CrossRef]
35. J. Laufer, Clare Elwell, Dave Delpy, et al., “In vitro measurements of absolute blood oxygen saturation using pulsed near-infrared photoacoustic spectroscopy: accuracy and resolution,” Phys. Med. Biol. 50(18), 4409–4428 (2005). [CrossRef]
36. A. Roggan, M. Friebel, K. Do Rschel, et al., “Optical properties of circulating human blood,” in Optical Diagnostics of Biological Fluids III, (SPIE, 1998), 70–82.
37. D. J. Pine, D.A. Weitz, J. X. Zhu, et al., “Diffusing-wave spectroscopy: dynamic light scattering in the multiple scattering limit,” J. Phys. 51(18), 2101–2127 (1990). [CrossRef]
38. V. Du Le and V. J. Srinivasan, “Beyond diffuse correlations: deciphering random flow in time-of-flight resolved light dynamics,” Opt. Express 28(8), 11191–11214 (2020). [CrossRef]
39. R. Bonner and R. Nossal, “Model for laser Doppler measurements of blood flow in tissue,” Appl. Opt. 20(12), 2097–2107 (1981). [CrossRef]
40. H. M. Varma, C. P. Valdes, A. K. Kristoffersen, et al., “Speckle contrast optical tomography: A new method for deep tissue three-dimensional tomography of blood flow,” Biomed. Opt. Express 5(4), 1275–1289 (2014). [CrossRef]
41. T. Binzoni and F. Martelli, “Assessing the reliability of diffuse correlation spectroscopy models on noise-free analytical Monte Carlo data,” Appl. Opt. 54(17), 5320–5326 (2015). [CrossRef]
42. J. M. Bland and D. Altman, “Statistical methods for assessing agreement between two methods of clinical measurement,” The Lancet 327(8476), 307–310 (1986). [CrossRef]
43. D. C. Montgomery, et al., Introduction to Linear Regression Analysis (John Wiley & Sons, 2021).
44. D. Mazumder, M. M. Wu, Nisan Ozana, et al., “Optimization of time domain diffuse correlation spectroscopy parameters for measuring brain blood flow,” Neurophotonics 8(3), 035005 (2021). [CrossRef]
45. B. Aernouts, R. Van Beers, R. Watté, et al., “Dependent scattering in Intralipid® phantoms in the 600-1850nm range,” Opt. Express 22(5), 6086–6098 (2014). [CrossRef]
46. S. L. Jacques, “Optical properties of biological tissues: a review,” Phys. Med. Biol. 58(11), R37–R61 (2013). [CrossRef]
47. R. B. Saager, M. Balu, V. Crosignani, et al., “In vivo measurements of cutaneous melanin across spatial scales: using multiphoton microscopy and spatial frequency domain spectroscopy,” J. Biomed. Opt. 20(6), 066005 (2015). [CrossRef]
48. T. Kono and J. Yamada, “In vivo measurement of optical properties of human skin for 450–800 nm and 950–1600 nm wavelengths,” Int. J. Thermophys. 40, 1–14 (2019).
49. 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]
50. H. J. Van Staveren, C. J. M. Moes, J. van Marie, et al., “Light scattering in lntralipid-10% in the wavelength range of 400–1100 nm,” Appl. Opt. 30(31), 4507–4514 (1991). [CrossRef]
51. S. T. Flock, S. L. Jacques, B. C. Wilson, et al., “Optical properties of Intralipid: a phantom medium for light propagation studies,” Lasers Surg. Med. 12(5), 510–519 (1992). [CrossRef]
52. M. Friebel, J. Helfmann, U. Netz, et al., “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(2), 034001 (2009). [CrossRef]
53. A. Kienle, M. S. Patterson, L. Ott, et al., “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]
54. M. Friebel, A. Roggan, G. Müller, et al., “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 functions,” J. Biomed. Opt. 11(3), 34021 (2006). [CrossRef]