1.

Introduction

The most common application of near-infrared (NIR) light for studying the human brain is functional near-infrared spectroscopy (fNIRS).^1,2 By shining constant [continuous-wave (CW)] NIR ($\sim 650$ to 950 nm) light into the scalp and by measuring the diffuse reflectance at different wavelengths ($\lambda$), CW-fNIRS allows us to study cortical changes in oxyhemoglobin (${\mathrm{O}}_2\mathrm{Hb}$) and deoxyhemoglobin (HHb) concentrations. These changes are caused by a variety of physiological events (e.g., pulse wave propagation within the head structure, periodic changes in vein pressure induced by respiration, and changes in arterial pressure induced by oscillations in baroreceptors and chemoreceptor reflex control systems, i.e., Mayer waves) including local vasodilatation and hyperemia induced by neural activity (neurovascular coupling phenomenon).³ This means that, similar to functional magnetic resonance imaging (MRI), which measures blood oxygen-level-dependent effects, fNIRS allows estimation and inference of brain activity.^4,5 To retrieve hemoglobin modulations through scalp-located measurements, models of light propagation through head structures are required. In fact, due to the diffusive nature of biological structures in such spectral range, NIR light propagation through tissues is a complex process.^6–9

Standard CW-fNIRS links light modulations to hemoglobin oscillations through a simplified model, defined as the modified Beer–Lambert law (MBLL).¹⁰ The MBLL is an exponential law that links the measured light to hemoglobin through the extinction coefficients (substance absorption per unit concentration) at the wavelengths employed and the effective distance traveled by the photons to reach the detection point. Importantly, the effective distance traveled by the photons in diffusive media is not equal to the source–detector geometrical distance. In fact, because of the scattering of photons in the tissue, the effective distance is several times the interoptode distance (between four and seven times longer) and, importantly, it changes as a function of light wavelength. This characteristic is incorporated in the MBLL by estimating the effective distance through simple multiplication between the geometrical distance and a scalar that is defined as the differential pathlength factor (DPF). Whereas the extinction coefficients are specific for the given substance and are tabulated,¹¹ the DPFs are dependent on the composition and microscopic and macroscopic structures of the tissue investigated; they show variability among subjects and head regions.¹²

Since empirical estimation of DPF requires photons’ average time-of-flight information,¹³ DPF is assumed a priori in CW recordings, often not accounting for intersubject and interhead region variability, or accounting only for a limited portion of variables affecting DPF (e.g., subject age). In fact, to provide quantitative estimations of DPF, photons’ average time-of-flight information is required. This information can be assessed either through photons’ distribution of time of flight (or temporal point spread function), acquired through time-domain technology,¹⁴ or through light phase shifts, acquired through frequency domain multidistance (FDMD) technology.¹⁵ However, these technologies are complex and expensive, and they are far less common than CW technology, particularly for functional neuroimaging applications.

Importantly, whereas an error in the average value of DPFs among different wavelengths influences the reconstructed oscillation amplitudes of hemoglobin forms, an error in the spectral dependence of DPF produces hemoglobin forms cross talk when applying MBLL.¹⁶ This aspect can be detrimental for fNIRS, because of the different, often opposite, behavior of ${\mathrm{O}}_2\mathrm{Hb}$ and HHb during local neural activity.

Another quantitative optical parameter linked to DPF is the effective attenuation coefficient (EAC). EAC quantifies the light exponential rate of attenuation and it is approximately a function of the product of absorption and reduced scattering (i.e., it is proportional to the geometric mean of absorption and reduced scattering).^17,18 Since DPF is approximately a function of the ratio of reduced scattering and absorption, it is possible to express the DPF as a ratio of reduced scattering and EAC. This means that, by measuring the EAC, the effect of absorption on DPF can be empirically accounted for, reducing the a priori assumption of DPF to the reduced scattering coefficient. Thus, assuming the spectral dependence of DPF, the a priori assumptions are reduced to the spectral dependence of reduced scattering. This is interesting because

In fact, whereas absorption spectra depend on the chromophores, and which chromophores show peculiar spectral features and which concentrations can vary among tissues and subjects, once scattering is determined at a given wavelength, several equations have been developed to describe its spectral characteristics. For example, a simple linear decay at increasing wavelength is described in the literature and is commonly employed in commercial systems to retrieve tissue saturation from multiwavelength analysis of the EAC.¹⁹ Thus, if the approximation describing wavelength dependence of scattering is correct, it is possible to estimate the DPF wavelength dependence from EAC estimates, ultimately dampening the problem of hemoglobin cross talk in CW-fNIRS introduced by an erroneous DPF spectrum assumption.

Chiarelli et al.²⁰ recently proposed a methodology for measuring EAC which relies on high-density CW-fNIRS recordings. The procedure measures the log light decay with distance employing common uncalibrated fNIRS systems. The algorithm is based on dampening the effects of variability in source intensity and optode-to-scalp coupling by measuring signals at several distances with many source–detector couples (channels, few tens). Specifically, light intensity is first estimated from high-frequency signal-to-noise ratio (SNR) in each channel, and, after, EAC is estimated as the slope of the log intensity of light reaching the detector as a function of source–detector distance through linear regression, providing an accurate and unbiased EAC estimate. Indeed, the algorithm requires a high-density optical array to work; however, based on the rapid development and increased usage of high-density CW-fNIRS systems, this approach could be useful in fNIRS to provide reliable estimates of hemoglobin oscillations through CW recordings.

In this paper, a validation of this approach is presented. Specifically, a measurement on the forehead of several subjects was performed employing a high-density, multidistance optical array and frequency-domain technology. Employing the algorithm described in Ref. 20, EAC was estimated from the uncalibrated CW component of the frequency-domain signal. The estimated EAC and DPF spectral dependencies were then compared to the results obtained employing standard FDMD measurements, where EAC was estimated from the phantom calibrated log light decay with distance, whereas DPF was estimated from the calibrated light phase shift with distance.²¹ Finally, the performances of the proposed procedure were assessed in the in-vivo fNIRS recordings acquired on the forehead from an independent sample of subjects. The subjects performed a mathematical task that supposedly induced functional activity in the frontal cortex.

2.

Methods

2.4.

Optical Recordings and Experimental Design

Optical data were acquired with a multichannel frequency-domain oximeter (ISS Imagent™, Champaign, Illinois) equipped with 32 laser diodes (16 emitting light at 690 nm and 16 emitting light at 830 nm) and 4 photomultiplier tubes (PMTs). Time-division multiplexing was employed so that each detector picked up light from 32 different sources at different times within a multiplexing cycle. The sampling rate was 39.0625 Hz. NIR light was carried to the scalp using single optic fibers (0.4-mm core) and from the scalp back to the PMTs using fiber bundles (3-mm diameter). The fibers were held in place using custom-built optical probe. The probe had an optical array geometry constituted of four rows (at a distance of 5 mm) of one detector and four dual-wavelength sources at a distance from the detector of 20, 25, 30, and 35 mm. The total area covered by the optical probe was $700{\mathrm{mm}}^2$ [Fig. 1(a)]. All source–detector couples were measured and employed for further analysis for a total of 64 dual-wavelength channels with an interoptode distance ranging from 20 to 46 mm. For comparison with the FDMD procedure, the subjects were requested to rest on an examination table. For FDMD calibration, signals were first acquired on the optical phantom for 1 min. The phantom had optical parameters equal to: ${\mu }_a=0.0107{\mathrm{mm}}^{-1}$ ${\mu }_s^{\prime }=0.97{\mathrm{mm}}^{-1}$ for 690 nm, ${\mu }_a=0.0103{\mathrm{mm}}^{-1}$ ${\mu }_s^{\prime }=0.90{\mathrm{mm}}^{-1}$ for 830 nm, and refraction index at both wavelengths $\eta =1.4$. Soon after, data were acquired on the subjects’ forehead for another 1 min. The probe was held in place by an operator [Fig. 1(b)].

Fig. 1

(a) Schematic representations of the optical probe employed in the study. The probe was constituted of an optical array of four rows (at a distance of 5 mm) of one-detector, four dual-wavelengths sources at a distance from the detector of 20, 25, 30, and 35 mm. All source–detector couples were recorded for a total of 64 channels with an interoptode distance ranging from 20 to 46 mm. (b) Example of an in-vivo measurement. The probe was held in place during the recording by an operator. (c) Simulation of the average light sensitivity of the 64 channels of the optical probe when located on the forehead. The simulation employed a FEM simulation based on a structural MRI of one participant in the study. The channels’ average sensitivity is displayed up to an attenuation of 60 dB (1000 times) from its maximum value. The probe investigated a localized portion of the head.

The fNIRS data were acquired with the same optical probe located on the same head location but on different subjects. The experimental paradigm included eight blocks. Each block was constituted of a 20-s prestimulus period (rest); a task phase, lasting a maximum of 32 s, and a poststimulus recovery phase, lasting 20 s. The task consisted of responding to three consecutive arithmetic subtractions (five digits number minus three digits number, e.g., $\mathrm{17,233}-271=\mathrm{16,962}$, $\mathrm{16,962}-271=\mathrm{16,691}$, and $\mathrm{16,693}-271=\mathrm{16,420}$). Each subtraction was presented for 5 s, and the participants were asked to choose the proper answer from a list of four possible results. A maximum of 7 s was allowed to respond and guess the correct subtraction value. The stimulation was developed using Psychtoolbox.²⁴ The execution of the mathematical task supposedly induced functional brain activity in frontal cortex, right under the area where the high-density optical probe was located.²⁵

Figure 1(c) displays a simulation assessing the average light sensitivity of the 64 channels of the optical probe when located on the forehead. The simulation employed a finite element method (FEM) approach.²⁶ FEM simulation was performed based on a structural MRI of one participant in the study. The channels’ average sensitivity is displayed up to an attenuation of 60 dB (1000 times) from its maximum value. The probe investigated a localized portion of the head.

2.5.

Data Analysis

The effect of errors in DPF was estimated through a simulated dual-wavelength estimation of functional hemoglobin oscillations. Event-related ${\mathrm{O}}_2\mathrm{Hb}$ and HHb temporal changes were simulated by convolving an impulse function (depicting stimulation) with the canonical hemodynamic response function.²⁷ A maximum positive change in ${\mathrm{O}}_2\mathrm{Hb}$ of $3\mu \mathrm{M}$ and a negative change in HHb of $1\mu \mathrm{M}$ were assumed. These hemoglobin oscillations are quantitative comparable with those obtained in vivo.²⁸ Here, $\mathrm{\Delta }\mathrm{ODs}$ were calculated through the MBLL assuming two wavelength recordings at 690 and 830 nm, with an interoptode distance of 30 mm and DPFs of 6 and 5, respectively. The effect of errors in the DPFs was estimated by inverting the MBLL employing wrong DPFs values, with a maximum percentage error in the DPFs of 20%.

Calibrated estimations of the EACs and DPFs were acquired through combining phantom and in-vivo recordings. Calibration coefficients were evaluated on the optical phantom by measuring the average CWs and phase shifts as a function of interoptode distance $d$ and by adjusting their value to $I_{\mathrm{calib}}$ and ${\mathrm{Ph}}_{\mathrm{calib}}$ to obtain the known phantom EACs and DPFs. The four rows were independently calibrated through the following equations:²⁹

Eq. (21)

$$\frac{\partial \ln (I_{\mathrm{calib}}(\lambda )d^2)}{\partial d}=d{\mu }_{\mathrm{eff}\_\text{phantom}}(\lambda ),$$

and

Eq. (22)

$$\frac{\partial \ln ({\mathrm{Ph}}_{\mathrm{calib}}(\lambda ))}{\partial d}=\frac{{\mathrm{DPF}}_{\text{phantom}}(\lambda )\omega }{v},$$

where $\omega$ is the modulation frequency of the recording system and $v$ is the speed of light in the medium $v=c/\eta$, where $\eta$ is the refraction index. Calibrated DPFs and EACs were estimated in vivo by applying the correction coefficients on the measured CWs and phase shifts. Each subject DPFs and EACs were estimated based on the average value across the different optodes’ rows. An example of calibrated CWs and phase shift signals for a particular row of the optical probe as a function of interoptode distance are reported in Fig. 2(a).

Fig. 2

(a) Example of calibrated CWs and phase shift signals for a particular row of the optical probe as a function of interoptode distance. Least-square regression line is also reported. (b) Example of the uncalibrated CW high-frequency SNR as a function of interoptode distance for all the 64 channels employed. Least-square regression line is also reported. The absolute slope of the regression represents the estimated EAC [Eq. (20)].

Uncalibrated estimates of EAC were obtained employing Eq. (21). The optical CW intensities were normalized and high-pass filtered above 10 Hz (fifth order, zero-lag, and Butterworth digital filter). The SNR was estimated over the in-vivo recording period of 1 min. The approach reported in Ref. 20 computes the EAC for each channel (source–detector pair) using subsets of channels near (neighbors) to the channel being estimated (within the neighborhood radius distance) in a multidistance configuration.

Because of the intrinsic variability introduced by source power, detector efficiency, and optode-to-scalp coupling, at a fixed source–detector distance range, the more channels employed in a single EAC computation from uncalibrated recordings the smaller the EAC confidence interval. Based on extensive analysis reported in Ref. 20, a numerosity of channels (source–detector couples) of few tens is sufficient to reduce the confidence interval of EAC to few points’ percentage.

In this study, a large neighborhood radius was employed to include all of the 64 source–detector couples in the EAC computation. An example of computation of EAC employing the uncalibrated procedure is reported in Fig. 2(b). The EAC values for all the 64 channels can be considered the EAC value of the area investigated by the optode layout (central forehead region) for the given wavelength [Fig. 1(c)].

Interestingly, because EAC is estimated from experimental evaluation of the SNR, the EAC confidence interval tends to be a monotonic decreasing function not only of the number of channels employed in EAC computation but also of the signal integration time (at a given sampling frequency). Indeed, signal variability and hence EAC can be estimated at a shorter integration time compared to the 1-min time window employed for further analysis in this study. Hence, to better highlight and investigate the EAC dependence on integration window, the variability of the EAC estimate as a function of integration time was also evaluated.

An among-subjects comparison between calibrated FDMD and uncalibrated EAC was performed to assess the developed algorithm performances; Further, the differences between calibrated and uncalibrated ratio of DPFs (spectral dependence of DPFs when considering two wavelengths) were assessed and compared to the expected values extracted from an equation model that estimated DPF based on light wavelength and subject’s age.¹⁰

Moreover, an analysis reporting the effect of assuming different constants $C$ [Eq. (15), i.e., different spectral dependencies of reduced scattering] on the unexplained amount of the among-subjects variance of the DPF spectral dependence was reported to address the validity of the assumptions of the work and the stability of the proposed approach.

Finally, fNIRS data were analyzed to assess the procedure performance in real fNIRS recordings. Raw CW signals were converted into $\mathrm{\Delta }\mathrm{ODs}$, movement corrected,³⁰ downsampled to 4 Hz, and bandpass filtered between 0.01 and 0.5 Hz to highlight functional activity (fifth order, zero-lag, and Butterworth digital filter). Hemoglobin oscillations were estimated from $\mathrm{\Delta }\mathrm{ODs}$ at the two wavelengths through the MBLL. Importantly, hemoglobin oscillations were evaluated based on the three different procedures for estimating DPFs. The first procedure assumed equal values of DPF at the two wavelengths employed (${\mathrm{DPF}}_{\mathrm{830,690}}=6$, $\mathrm{DPF}\text{ratio}=1$), the second procedure estimated DPFs for both wavelengths based on Ref. 10, and the third procedure assumed the absolute value of one of the two DPF from Ref. 10, whereas the other was extracted based on the DPF ratio retrieved by the algorithm developed. Block-average oscillations in ${\mathrm{O}}_2\mathrm{Hb}$ and HHb induced by the mathematical task were assessed from 5 s prior to 40 s after the task onset. The performance of procedures in estimating the DPF was assessed by comparing the reconstructed amplitudes of task-induced functional responses of ${\mathrm{O}}_2\mathrm{Hb}$ and HHb among the different algorithms. Channels showing the strongest functional activation were selected for the analysis.

3.

Results

Figure 3(a) reports an example of simulated hemodynamic response to an event-like stimulation. The classical temporal characteristics of functional activity are visible. Original ${\mathrm{O}}_2\mathrm{Hb}$ and HHb oscillations are reported with continuous lines, whereas the reconstructed ${\mathrm{O}}_2\mathrm{Hb}$ and HHb oscillations are reported with dashed lines. The reconstructed signals (i.e., dashed lines) were evaluated by inverting the MBLL with the introduction of an error of 10% in the dual-wavelength spectral dependence of DPFs. Figure 3(b) reports a colormap image where color represents the error in the absolute ratio of ${\mathrm{O}}_2\mathrm{Hb}$ and HHb as a function of the error in DPFs for the two wavelengths employed in the simulation. The figure clearly shows that hemoglobin cross talk is introduced only when there is a modification in the wavelength–DPF ratio, i.e., the spectral dependence of DPF is altered. Notice that the errors in the DPFs introduced in the simulation are compatible with errors that can be easily performed in-vivo when a priori values of DPFs are assumed without considering intersubject and interhead region variability. These errors clearly show a possible hemoglobin cross talk of the same order of magnitude of the error in the DPF ratio.

Fig. 3

(a) Example of the simulated functional brain activation induced by an event stimulation as measured by fNIRS. The continuous line reports the original ${\mathrm{O}}_2\mathrm{Hb}$ and HHb oscillations, whereas the dashed line depicts the estimated ${\mathrm{O}}_2\mathrm{Hb}$ and HHb oscillations when inverting the MBLL with errors in DPFs. (b) Colormap image showing the error in the absolute ratio of ${\mathrm{O}}_2\mathrm{Hb}$ and HHb as a function of the error in DPFs at the two wavelengths employed in the simulation.

Figure 4 shows correlation plots and associated Bland–Altman plots³¹ reporting the uncalibrated versus calibrated EAC values. Figure 4(a) reports the 690-nm wavelength, whereas Fig. 4(b) reports the 830-nm wavelength. Significant correlations between the calibrated and uncalibrated EACs for both wavelengths were obtained ($r=0.94$, $p=1.4·{10}^{-10}$ for 690 nm and $r=0.85$, $p=9.0·{10}^{-7}$ for 830 nm). For both wavelengths, the slope was not significantly different from 1 ($p=0.40$ for 690 nm and $p=0.41$ for 830 nm). The calibrated versus uncalibrated EAC differences were always within the acceptance criteria of Bland–Altman plot, except for one subject. The average among subjects’ EAC difference was not significantly different from 0 ($t=1.77$, $df=20$, $p=0.09$ for 690 nm and $t=-0.99$, $df=20$, $p=0.33$ for 830 nm).

Fig. 4

Correlation plots and associated Bland–Altman plots reporting the uncalibrated versus the calibrated EAC values for (a) 690-nm wavelength and (b) 830-nm wavelength.

Figure 5 reports among-subject average (and related standard error) EAC estimation variability [standard deviation (STD)] with respect to 1-min integration time (%) as a function of integration time window. A clear monotonic decrease with integration time occurs, with an average variability ranging from about 4% at a 1 s to about 1% at 32-s integration time. This means that EAC can indeed be estimated at a fast pace at the expense of a loss of confidence in its actual value (reduced number of points in SNR evaluation). However, decreasing the integration time window does not drastically increase the confidence interval of EAC estimate, in theory allowing fast evaluation of EAC.

Fig. 5

Among-subject average (and related standard error) EAC estimation variability (STD) with respect to 1-min integration time (%) as a function of integration time window.

Figure 6 reports the correlation plot and associated Bland–Altman plot between the DPF ratio at the two wavelengths of interest computed from absolute DPF measurements from calibrated FDMD phase analysis and uncalibrated EAC-based computation. The EAC-based DPF ratio estimation was obtained with an optimized scattering-related constant $C$ [Eq. (15)] of 0.76 (refer to Fig. 7). Bland–Altman plot shows the consistency of the two evaluation procedures with no between measurement difference on average ($t=0.91$, $df=20$, $p=0.38$). The correlation plot also reports the expected DPF ratio estimated employing an equation model commonly used in fNIRS analysis that estimates the subjects’ DPFs based on light wavelength and subjects’ age.¹⁰

Fig. 6

Correlation plot and associated Bland–Altman plot considering the ratio of the DPFs at the two wavelengths of interest computed from calibrated absolute DPFs and uncalibrated EAC measurements. The correlation plot also reports the expected ratios between DPFs estimated, employing an equation model commonly used in CW-fNIRS analysis that estimates the DPFs based on light wavelength and subject age.

Fig. 7

Effect of assuming different $C$ constant [Eq. (15)] in the proposed algorithm on the RMSE of estimated DPF ratio with respect to the FDMD within the experimental data sample.

Figure 7 reports the analysis of the effect of assuming different $C$ constants [Eq. (15)] in the proposed algorithm. In fact, since the correct ratio of DPFs (measured through calibrated light phase delay) was known, it was possible to estimate the residual variance [variance in the across-subject DPF ratio or root mean square error (RMSE)] of the model, and to compare that to the original variance (or STD) of the real DPF ratio. The experimentally measured across-subject variability (STD) of the real DPF ratio was 0.15 (15%). The model reduced the unexplained variance of the ratio (RMSE between real and predicted values) at a minimum of $\sim 0.04$ (4%) at a constant $C$ of 0.76 (value used for Fig. 6 analysis). Below and above this value, the algorithm’s performance, as assessed by the RMSE metric, decreased.

Figure 8 reports the results from real fNIRS experiments in which the performance of the developed algorithm was compared to procedures where a constant wavelength dependence of DPF was assumed or a priori values of DPFs from a published equation model were selected.¹⁰ Figure 8(a) displays the classical functional responses to the task in selected highly activated channels of the optical probe for each of the six subjects who underwent the experiment. Interestingly, for all the subjects examined, an increase in the positive amplitude of ${\mathrm{O}}_2\mathrm{Hb}$ and a negative amplitude of HHb is visible, when sequentially applying equal DPF values between wavelengths, the a priori equation model, and the developed procedure. To statistically quantify the improved performances of the developed algorithm, the maximum amplitudes of the absolute oscillation of ${\mathrm{O}}_2\mathrm{Hb}$ and HHb were extracted from each subject and procedure. Figure 8(b) shows the boxplots depicting the percentage increase in absolute maximum response in hemoglobin with respect to assuming equal DPFs between wavelengths when the developed algorithm and the a priori equation model were employed. With respect to assuming constant DPF, when applying the developed EAC-based algorithm, an average positive increase in maximum ${\mathrm{O}}_2\mathrm{Hb}$ oscillations during task of 26.9% ($t=5$, $df=5$, $p=0.004$) and a negative increase in maximum HHb oscillations of 7.4% ($t=1.44$, $df=5$, N.S.) were obtained. When applying the a priori equation model, an average positive increase in ${\mathrm{O}}_2\mathrm{Hb}$ of 10.7% ($t=7.24$, $df=5$, $p=8\times {10}^{-4}$) and a negative increase in HHb of 1% ($t=0.58$, $df=5$, $p=\mathrm{N.S.}$) were obtained. Interestingly, when comparing the EAC-based procedure with the equation model, a further positive increase in ${\mathrm{O}}_2\mathrm{Hb}$ of 16.2% ($t=3.58$, $df=5$, $p=0.015$) and a negative increase in HHb of 6.4% ($t=1.67$, $df=5$, $p=\mathrm{N.S.}$) were obtained.

Fig. 8

Results from the real fNIRS experiments in which the performance of the developed algorithm was compared to procedures where a constant wavelength dependence of DPF was assumed or a priori values of DPFs from published equation model were selected.¹⁰ (a) Classical functional responses to the task in selected highly activated channels of the optical probe for each of the six subjects that underwent the experiment. (b) Boxplots depicting the percentage increase in absolute maximum response in hemoglobin with respect to assuming equal DPFs between wavelengths when the developed algorithm and the a priori equation model were employed.

4.

Discussion

Correct estimation of DPF spectral dependence is crucial in CW-fNIRS, since errors in the evaluation of its spectrum creates cross talk in the retrieval of ${\mathrm{O}}_2\mathrm{Hb}$ and HHb modulations from light recordings (Fig. 3). This cross talk can be detrimental for fNIRS since it is well known that ${\mathrm{O}}_2\mathrm{Hb}$ and HHb have very different (often opposite) evoked functional responses induced by neural activity. In fact, the sensitivity to functional brain activity in fNIRS is often assumed when there is an anticorrelation between the two hemoglobin forms.³²

However, DPF has to generally be assumed a priori in CW-fNIRS. Although several works that attempted to provide reliable tabulated values or equation models of DPF to be employed for CW-fNIRS measurements,^10,21 the intersubjective variability and DPF dependence on a large number of factors, such as age, sex, head region,³³ and disease, make it particularly challenging to provide reliable estimation. Although it is not possible to quantitatively evaluate DPF from standard CW-fNIRS recordings because of the lack of photons’ time-of-flight information, in this paper a procedure to estimate DPF spectral dependence was presented. Multidistance measurements of the EAC (which is a function of the product of absorption and reduced scattering) from high-density CW-fNIRS recordings was employed to empirically estimate DPF (which is a function of the ratio of reduced scattering and absorption) spectrum by limiting the a priori assumption to the spectral dependence of reduced scattering, which is generally more predictable than absorption spectrum. This aspect requires further commenting. Indeed, the well-known decrease in wavelength of the DPF is strictly associated with the (also well-known) decrease in reduced scattering coefficient with wavelength; if, for example, reduced scattering is assumed constant with wavelength, an opposite trend of DPF spectrum is expected (i.e., DPF increases with wavelength). In fact, biological tissues are generally more absorbing at shorter wavelength (because of HHb, melanin, etc.); this aspect, combined with the effect of absorption of reducing the average time-of-flight of photons in tissues [absorption coefficient as denominator in Eq. (11)], causes an increase in DPF with wavelength when a constant reduced scattering is assumed. In reality, the decrease with wavelength of reduced scattering [nominator in Eq. (11)] is stronger than the increase in absorption [denominator of Eq. (11)] causing the experimentally found decrease of the DPF with wavelength. However, the average effect of reduced scattering and absorption on the expected value of DPF must be conceptually uncoupled from the across-subject, across-head region variabilities that they induce. In fact, although some of the variance of the across-subject, across-head region DPF spectrum indeed cannot be explained without a direct assessment of reduced scattering, a good portion of the variance can be accounted for with measures of EAC. The basic idea of the work is that the behavior of reduced scattering as a function of wavelength is more predictable than absorption. The assumption of fixing the spectral dependence of reduced scattering is in fact commonly employed in pulse oximeters (which are generally calibrated in vivo)³⁴ and multidistance tissue oximeters that measure tissue oxygenation index¹⁹ through spatially resolved spectroscopy (i.e., that estimate tissue saturation from EAC measurements).

The DPF spectral dependence was investigated on a homogeneous population (with age between 25 and 40 years) of the same ethnicity (Caucasian), on the same head region (forehead), to assess the intrinsic intersubjective variability which, in general, cannot be summarized by a single equation accounting for light wavelength and subject age. By using a commonly employed equation for estimating scattering dependence [Eq. (14)], it has been possible to obtain a reliable estimate of the DPF change with wavelength (Fig. 6), from EAC measurements [$r=0.83$, Eq. (15)]. Notice that this result corroborated the assumptions of the procedure developed. In fact, the 69% of the variance was explained employing the assumption of a linear association between variables [well-defined constant $C$ in Eq. (15)]. This means that, in the population examined, the across-subject variability in the constant $C$ is not introducing a large effect on the total across-subject variance of DPF ratio.

However, there are indeed limitations to the proposed approach. In fact, reduced scattering and absorption contribute to the DPF in a multiplicative and not in a summation manner. This means that, the across-subject, across-head region variance explained by the two coefficients cannot be uncoupled in a simple summation of effects. This results in the fact that, if the across-subject, across-head region constant $C$ is highly incorrect, the predicted ratio of DPF can be highly incorrect, introducing an error which is higher than the original across-subject, across-head region variability. An analysis reporting the effect of assuming different $C$ constants in the algorithm is reported in Fig. 7. In fact, since the correct ratio of DPFs (measured through calibrated light phase delay) was known, it was possible to estimate the residual variance (variance in the across-subject DPF ratio or RMSE) of the model, and to compare that to the original variance (or STD) of the real DPF ratio. The measured across-subject variability (STD) of the real DPF ratio was 0.15 (15%), which could induce a same order of magnitude hemoglobin cross talk. The developed procedure reduced the unexplained variability (RMSE) of the ratio at a minimum of 0.04 (4%) at a constant $C$ of 0.76. Notice that for a range of $C$ values between 0.66 and 0.86, the model still reduced the unexplained variance with respect to the original variance of the data. This means that, for this dataset, the model helped to improve the DPF ratio estimate with respect to chance. A $C$ constant between 0.66 and 0.86 encompasses a large portion of variability in the drop of reduced scattering. Notice, however, that below and above this value, the model worsens the result compared to chance. This effect is introduced by the multiplication of the coefficient involved in the DPF ratio evaluation. This means that, although the experimentally established $C$ coefficient should not vary substantially as a function of subjects and head regions, it would be helpful in future work to systematically map across-regions and across-subject variability of the $C$ coefficient to help and reduce the across-subjects’ across-head regions’ unexplained variance in the DPF ratio of the procedure to a minimum. This would result in introducing a hybrid approach where a combination of a priori assumption on $C$ constant and experimentally evaluated EAC might improve the prediction accuracy of DPF spectral dependence with respect to a completely a priori assumption on DPFs.

Importantly, the subject-based spectral dependence of DPF was not predicted employing a general equation that estimates the DPF based on wavelength and subject age as the only parameters.¹⁰ However, notice that the results presented in this paper are not in disagreement with the procedure and results reported in Ref. 10. In fact, the general equation for predicting DPF relies on expressing the DPF, at a certain wavelength, as a function of age only, not accounting for possible other sources of variance. Moreover, whereas this equation might be accurate in predicting the DPF value, the error of prediction of the ratio of DPF as a function of wavelength (i.e., the spectral dependence) can be amplified. In this paper, because of the small age range examined, the possible residual variance not accounted for by the equation model was investigated. This means that an orthogonal dimension of the DPF spectral dependence variability was examined. Interestingly, the almost constant ratio of DPF within the population retrieved by the equation model was in fact exactly equal to the across-subject average value found employing the proposed approach, further corroborating this hypothesis. For example, it would be possible to fix one DPF employing the general equation of Ref. 10 (to provide quantitative estimation of hemoglobin concentrations oscillations), and to estimate the other DPF through the described approach. Indeed, the two approaches should be defined as complementary. In fact, this rationale was also applied for the real fNIRS validation of the approach.

Importantly, the fNIRS results confirmed the utility of the developed algorithm. In fact, among subjects significant increases in the functional task-induced hemoglobin oscillations were obtained employing the developed procedure with respect to both assuming constant wavelength dependence of DPF and employing a priori values of DPFs. The increase in the retrieved functional activity ranged from 16.2% to 26.9% for ${\mathrm{O}}_2\mathrm{Hb}$ and was highly statistically significant. The increase in functional activity of HHb ranged from 6.4% to 7.4% and, although not being statically significant, its trend was clear. The lower effect of the procedure on HHb could be associated with a lower SNR of HHb with respect to ${\mathrm{O}}_2\mathrm{Hb}$. Notice that, these increases in retrieved functional oscillations were not derived from a simple scaling of the functional response, but they were actually associated with an increase in the SNR of the fNIRS measurement.

Interestingly, the EAC estimation used for DPF spectral dependence computation was extracted through uncalibrated CW-fNIRS recordings based on a method²⁰ that relied on accounting for the variability in detector sensitivity and optode-to-scalp coupling based on large number of multidistance high-density optical channels and statistical procedures. In fact, the EAC estimates reported here were evaluated employing a high-density configuration featuring overlapping channels. The optical probe covered a small region of the forehead; thus, it provided a local estimate of the metric of interest [forehead, Fig. 1(c)] and was not assumed to provide a whole-head estimate. The purpose of this work was to demonstrate the capabilities and performances of the developed algorithm in providing local estimate of the metric of interest without assuming the retrieved value to be valid for the whole head. An accurate mapping of the EAC across the head is beyond the scope of the reported work. By employing large field-of-view high-density optical array, this algorithm could be employed in a straightforward manner to provide whole head maps by subdividing the optical probe in multiple head regions and by iterating the procedure within different portions of the optical probe (in similarity with Ref. 20).

This algorithm opens the possibility of estimating EAC (and DPF spectral dependence), on standard CW-fNIRS recordings with the only requirements of having a large number of optodes in a high-density configuration. Notice that, whereas the high number of optodes is required for proper statistical analysis and a high confidence in the retrieved EAC values, the high density is required to provide sufficient localization to the EAC and DPF ratio estimations limiting possible partial volume effects. For example, in the current study, 64 channels over an area of $\sim 700{\mathrm{mm}}^2$ were employed for a single EAC computation; thus, the underlying assumption of the procedure is that the EAC spatial changes within the investigated area are minor. Although it is known that DPF changes as a function of brain region,³³ these changes should be sufficiently smooth to employ such a multichannel algorithm relying on sufficient optode density.

Another interesting application of the developed approach would be to study the time-varying properties of DPF spectral dependence as a function of functional brain activity. In fact, brain activity indeed causes changes to the baseline optical parameters of several point percentages that influence DPF and are typically overlooked.³⁵ These changes might by assessed by the developed procedure through a sufficiently short integration time window of the algorithm. For this reason, the stability of the estimated EAC as a function of integration time was evaluated. As displayed in Fig. 7, EAC estimation variability (STD) decreases with integration time from about 4% at a 1-s integration to about 1% at 32-s integration time. This means that EAC can be estimated at a fast pace at the expense of a loss of confidence in its actual value. However, decreasing the integration time window does not drastically modify the confidence in its value, in theory allowing fast evaluation of EAC.

Because of the rapid development of high-density systems for diffuse optical tomography analysis of functional brain activity through CW-fNIRS, the procedure opens the possibility to provide quantitative evaluation of baseline optical properties that can be useful for correct estimation of ${\mathrm{O}}_2\mathrm{Hb}$ and HHb oscillations in the brain. Correct estimation ${\mathrm{O}}_2\mathrm{Hb}$ and HHb functional oscillations increases fNIRS SNR and reduces false positives and false negatives³⁶ in the evaluation of brain activity with fNIRS.

5.

Conclusion

In this paper, a procedure that estimates the spectral dependence of DPFs in CW-fNIRS was validated. Although it is not possible to quantitatively estimate DPF from CW light measurements, by employing uncalibrated high-density recordings, it is possible to evaluate the rate of decay of light as a function of distance (quantity defined EAC) and to limit the a priori assumption of DPF to reduced scattering. Because reduced scattering spectral dependence is more predictable than absorption, this procedure can be useful to dampen hemoglobin cross talk in fNIRS induced by erroneous DPF spectral dependence, ultimately increasing the reliability and usefulness of fNIRS in investigation of brain functional activity. Importantly, if a high-density optode configuration is provided, this procedure does not require alteration of the CW-fNIRS recordings.