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 D_B, 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, g₁(τ_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 g₁(τ_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, g_dyn, and fewer dynamic scattering events [38]. Since RBCs have g_dyn∼0.98 and we need to analyze g₁(τ_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, q²(θ) is the momentum transfer as a function of dynamic scattering angle θ, <Δr²(τ_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 g₁(τ_s,τ_d) even with few dynamic scattering events and high g_dyn [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[〈Δr²(τ_d)〉⁴]. Incorporating hybrid motion [40,41], 〈Δr²(τ_d)〉=6D_Bτ_d+v²τ_d², where v² is the second moment of the Gaussian velocity distribution and arranging different powers of τ_d (see full details in Appendix A), we have

O(τ_d⁴) 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}$v_cτ_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 V_pp 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 V_pp 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.

 

Fig. 1. Dual wavelength iNIRS setup. All fibers are single mode with angle polished connectors. I/C controller: current controller; DFB: distributed feedback laser; L,L’: lenses; APP: anamorphic prism pair; WDM: wavelength division multiplexer; DBD: dual balanced detector. Bold lines denote optical paths, while dotted lines denote electrical paths. Path 1 (red) is for 855 nm, while path 2 (blue) is for 773 nm.

Download Full Size | PDF

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,t_d) as described in Ref. [30]. Here t_d 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 g_dyn∼0.6, while paths in biological tissue involve scattering from static (blue) and dynamic (red) scatterers [Fig. 2(b)], the latter with g_dyn∼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 v_c, wavelength-dependence of ξ(τ_s) should arise from k²= (2πn/λ)² 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 λ₁ from ${\xi ^{{\lambda _2}}}({{\tau_s}} )$ at wavelength λ₂ via

In Eq. (11), $\mu_{s, d y n}^{\prime}$(λ₁) and $\mu_{s, d y n}^{\prime}$(λ₂) are reduced dynamic scattering coefficients at wavelengths λ₁ and λ₂, 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),

 

Fig. 2. (a) Light paths in an emulsion such as Intralipid include only dynamic scattering events. (b) Light paths in vascularized biological tissues include both static and dynamic scattering events. Therefore, for an emulsion, $\mu_{s, d y n}^{\prime}$=$\mu_{s}^{\prime}$ and α=1, whereas in biological tissue with blood flow, $\mu_{s, d y n}^{\prime}$≠$\mu_{s}^{\prime}$ and α≠1, in general.

Download Full Size | PDF

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}$:

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, τ_d² and τ_d³ are represented as ξ₁(τ_s), ξ₂(τ_s) and ξ₃(τ_s), respectively. To provide a simple framework to interpret the coefficients in Eq. (13), we take v = 0 in Eq. (8) and neglect O(τ_d⁴):

Following the constant $\mu_{s, d y n}^{\prime}$-scaling approach and using Eqs. (5), (6), (7) with Eq. (14), the coefficients ξ₁(τ_s), ξ₂(τ_s) and ξ₃(τ_s) must be multiplied by λ², λ⁴ and λ⁶ respectively to remove their wavelength-dependency, as the ratios of the quantities 〈sin²(θ/2)〉, 〈sin⁴(θ/2)〉 and 〈sin⁶(θ/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 λ₁ and λ₂ respectively, $\xi _2^{{\lambda _1}}({{\tau_s}} )$ and $\xi _2^{{\lambda _2}}({{\tau_s}} )$ correspond to the coefficients of τ_d² at wavelengths λ₁ and λ₂ 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 λ₁ and λ₂ 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.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 (R²) [43]. R² quantifies the goodness of agreement between the predicted and measured $g_1^{\lambda_2}$(τ_s,τ_d). We used 1-R² to compare scaling approaches, where a lower value signifies better agreement between predicted and measured quantities.

3. Results

We begin by fitting G₁(τ_s,τ_d) from Intralipid (g_dyn∼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., g₁(τ_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].

 

Fig. 3. Measured G₁ and single exponential fit for Intralipid mixture at different τ_s (855 nm) (EARLY τ_s: 60 ps, MID τ_s: 330 ps and LATE τ_s: 500 ps). The IRF and corresponding TPSF are shown in the inset of Fig. 3(a).

Download Full Size | PDF

 

Fig. 4. Same as above but for 773 nm. The IRF and corresponding TPSF are shown in the inset of Fig. 4(a).

Download Full Size | PDF

Unlike the Intralipid phantom, G₁(τ_s,τ_d) from in-vivo measurements (g_dyn∼0.98) do not decay exponentially. Thus, G₁(τ_s,τ_d) from both wavelengths were fitted to the unnormalized form of Eq. (13) [Fig. 5, Fig. 6]. For in-vivo data, we fit g₁(τ_s,τ_d) ranging from 1 to 0.1. Since a single exponential model does not fit in-vivo G₁(τ_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.

 

Fig. 5. Measured G₁, single exponential fit, and Eq. (13) fit for human forearm at different τ_s (855 nm) (EARLY τ_s: 60 ps, MID τ_s: 330 and LATE τ_s: 500 ps). The IRF and corresponding TPSF are shown in the inset of Fig. 5(a).

Download Full Size | PDF

 

Fig. 6. Same as Fig. 5 but for 773 nm. The IRF and corresponding TPSF are shown in the inset of Fig. 6(a).

Download Full Size | PDF

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, G₁(τ_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.

 

Fig. 7. Comparison of Intralipid $g_1^{855}\left(\tau_s, \tau_d\right)$ and $g_1^{773}\left(\tau_s, \tau_d\right)$ at early τ_s of 60 ps with (a) original τ_d -axes, (b) scaled τ_d -axis for $g_1^{855}\left(\tau_s, \tau_d\right)$ [Eq. (12)] assuming constant α and original τ_d -axis for $g_1^{773}\left(\tau_s, \tau_d\right)$, (c) scaled τ_d -axis for $g_1^{855}\left(\tau_s, \tau_d\right)$ [Eq. (12)] assuming constant $\mu_{s, d y n}^{\prime}$ and original τ_d -axis for $g_1^{773}\left(\tau_s, \tau_d\right)$, (d) Bland-Altman plots showing the differences between $g_1^{855}\left(\tau_s, \tau_d\right)$ and $g_1^{773}\left(\tau_s, \tau_d\right)$ for a-c.

Download Full Size | PDF

 

Fig. 8. Same as Fig. 7 but for MID τ_s of 330 ps.

Download Full Size | PDF

 

Fig. 9. Same as Fig. 7 but for LATE τ_s of 500 ps

Download Full Size | PDF

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.

 

Fig. 10. Same as Fig. 7 but for in-vivo data (EARLY τ_s: 60 ps)

Download Full Size | PDF

 

Fig. 11. Same as above but for MID τ_s of 330 ps

Download Full Size | PDF

 

Fig. 12. Same as above but for LATE τ_s of 500 ps. Differences in scaling approaches are less apparent due to low SNR.

Download Full Size | PDF

Finally, to comprehensively compare the scaling techniques, we calculated 1-R² (where R² 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.

 

Fig. 13. Ability to predict $g_1^{773}\left(\tau_s, \tau_d\right)$ from $g_1^{855}\left(\tau_s, \tau_d\right)$ by scaling the τ_d -axis for (a) Intralipid phantom and (b) human forearm. Note that 1-R² is shown, where R² is the coefficient of determination. The constant α-scaling approach performs best for the purely dynamic Intralipid phantom whereas the constant $\mu_{s, d y n}^{\prime}$-scaling approach performs best in-vivo, with mixed static and dynamic scattering. At late τ_s, differences between scaling approaches are less apparent due to low SNR.

Download Full Size | PDF

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 G₁ at all TOFs, and that the constant α-scaling approach works best for Intralipid.

 

Fig. 14. (a) Recovered ξ(τ_s) with 95% confidence interval (CI) from single exponential fit of G₁(τ_s,τ_d) from Intralipid phantom for both wavelengths. (b) Same as (a) but ξ(τ_s) at 855 nm was scaled using Eq. (11).

Download Full Size | PDF

Next, we proceeded to consider scaling of terms in the τ_d-cumulant expansion for in vivo measurements. The in-vivo G₁(τ_s,τ_d) were fitted to Eq. (13) for the coefficients ξ₁(τ_s), ξ₂(τ_s) and ξ₃(τ_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.

 

Fig. 15. (a) Recovered ξ₁(τ_s) from Eq. (13) fitting of G₁(τ_s,τ_d) from human forearm for both the wavelengths. (b) Same as (a) but ξ₁(τ_s) at 855 nm was scaled using Eq. (15).

Download Full Size | PDF

 

Fig. 16. a) Recovered ξ₂(τ_s) from Eq. (13) fitting of G₁(τ_s,τ_d) from human forearm for both the wavelengths. (b) Same as (a) but ξ₂(τ_s) of 855 nm wavelength was scaled using Eq. (16).

Download Full Size | PDF

 

Fig. 17. a) Recovered ξ₃(τ_s) from Eq. (13) fitting of G₁(τ_s,τ_d) from human forearm for both the wavelengths. (b) Same as (a) but ξ₃(τ_s) at 855 nm was scaled using Eq. (17).

Download Full Size | PDF

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)², referred here as λ²-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].

 

Fig. 18. (a) Fitting $g_1^{773}$ at early τ_s of 60 ps with Eq. (13); (b) Predicting $g_1^{773}$ at early τ_s by scaling the recovered 855 nm ξ₁(τ_s), ξ₂(τ_s) and ξ₃(τ_s) based on theory [Eqs. (15), (16) and (17)]; (c) Same as (b) but scaling each constant with (855/773)²; and (d) Bland-Altman plot showing difference between fitting and measurement of (a), between theory-based prediction and measurement of (b), and between λ² scaling-based prediction and measurement of (c). (EARLY τ_s: 60 ps, MID τ_s: 330 ps and LATE τ_s: 500 ps)

Download Full Size | PDF

 

Fig. 19. Same as Fig. 18 but for mid τ_s of 330 ps.

Download Full Size | PDF

 

Fig. 20. Same as above but for late τ_s of 500 ps. Differences in scaling approaches are less apparent due to low SNR.

Download Full Size | PDF

Finally, we summarized the error (1-R²) 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-R² 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 τ_dⁿ coefficient is scaled by λ²ⁿ outperforms the approach where all coefficients are scaled by λ². This suggests that higher-order τ_d terms scale more strongly with wavelength than the first-order (exponential decay) term.

 

Fig. 21. Comparison of the prediction error (1-R²) of higher cumulant-scaling (based on theory, and λ²-scaling). The error of a direct $g_1^{773}\left(\tau_s, \tau_d\right)$ fit to Eq. (13) serves as a benchmark for comparison. At late τ_s, differences between scaling approaches are less apparent due to low SNR.

Download Full Size | PDF

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=αD_B is also wavelength-dependent (as clearly, D_B is wavelength-independent).

 

Fig. 22. (a) Spectral dependence of the scattering coefficient (μ_s) of whole blood from Ref. [34]. (b) Derived $\mu_{s}^{\prime}$ = μ_s(1-g) of whole blood using anisotropy factor (g) and μ_s from Ref. [34]. (c) Spectral dependence of tissue $\mu_{s}^{\prime}$ with three different scattering powers (b) [46]. (d) Derived probability of dynamic scattering, α, from (b)-(c), assuming a 5% blood volume fraction.

Download Full Size | PDF

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⁻⁷ (±3.14 × 10⁻⁸) mm²/sec for 855 nm and 1.57 × 10⁻⁷ (±4.72 × 10⁻⁸) mm²/sec for 773 nm, across three subjects who participated in the study.

Table 1. BFI of Intralipid and in-vivo forearm with their ratio across wavelengths

View Table

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}$D_B, 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 (g_dyn∼0.6 [49–51]), but inadequately captures the slowly decaying TOF-resolved field autocorrelation tails in vivo (g_dyn∼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 g₂= 1+|g₁|² (intensity autocorrelation) as compared to g₁ (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 [ξ₁(τ_s), ξ₂(τ_s), and ξ₃(τ_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.