Differential pathlength factor in continuous wave functional near-infrared spectroscopy: reducing hemoglobin’s cross talk in high-density recordings

Abstract. Functional near-infrared spectroscopy (fNIRS) estimates the functional oscillations of oxyhemoglobin and deoxyhemoglobin in the cortex through scalp-located multiwavelength recordings. Hemoglobin oscillations are inferred through temporal changes in continuous-wave (CW) light attenuation. However, because of the diffusive multilayered head tissue structures, the photon path is longer than the source–detector separation, complicating hemoglobin evaluation. This aspect is incorporated in the modified Beer–Lambert law where the source–detector distance is multiplied by the differential pathlength factor (DPF). Since DPF estimation requires photons’ time-of-flight information, DPF is assumed a priori in CW-fNIRS. Importantly, errors in the DPF spectrum induce hemoglobin cross talk, which is detrimental for fNIRS. We propose to estimate subject-specific DPF spectral dependence relying on multidistance high-density measurements. The procedure estimates the effective attenuation coefficient (EAC), which is proportional to the geometric mean of absorption and reduced scattering. Since DPF depends on the scattering-to-absorption ratio, EAC limits the spectral dependence assumption to scattering. This approach was compared to a standard frequency-domain multidistance procedure. A good association between the two methods (r2  =  0.69) was obtained. This approach could estimate low-resolution maps of the DPF spectral dependence through large field of view, high-density systems, reducing hemoglobin cross talk, and increasing fNIRS sensitivity and specificity to brain activity without instrumentation modification.


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 (∼650 to 950 nm) light into the scalp and by measuring the diffuse reflectance at different wavelengths (λ), CW-fNIRS allows us to study cortical changes in oxyhemoglobin (O 2 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). 3 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][7][8][9] Standard CW-fNIRS links light modulations to hemoglobin oscillations through a simplified model, defined as the modified Beer-Lambert law (MBLL). 10 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, 11 the DPFs are dependent on the composition and microscopic and macroscopic structures of the tissue investigated; they show variability among subjects and head regions. 12 Since empirical estimation of DPF requires photons' average time-of-flight information, 13 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, 14 or through light phase shifts, acquired through frequency domain multidistance (FDMD) technology. 15 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. 16 This aspect can be detrimental for fNIRS, because of the different, often opposite, behavior of O 2 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 1. it is possible to estimate the EAC by measuring CW light attenuation as a function of distance without the need of photons' time-of-flight information; and 2. reduced scattering coefficient changes as a function of wavelength are generally less variable than absorption changes.
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. 19 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. 20 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. 21 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.

Modified Beer-Lambert Law, Differential Pathlength Factor, and Effective Attenuation Coefficient
In the first approximation, the CW wavelength-dependent diffuse reflectance can be modeled using the MBLL as 10 E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 1 ; 3 2 6 ; 4 1 9 where IðλÞ is the measured wavelength-dependent diffuse reflected light intensity, I o ðλÞ is the injected light, μ a ðλÞ is the absorption coefficient of the probed tissue, d is the geometrical source-detector distance, DPFðλÞ is the differential pathlength factor, and G is a medium-and geometry-dependent constant. Equation (1) can be rewritten as E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 2 ; 3 2 6 ; 3 2 1 where ODðλÞ is defined as the optical density. If multiple chromophores contribute to the total absorption μ a ðλÞ, it is possible to further rewrite Eq. (2) as E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 3 ; 3 2 6 ; 2 4 3 where n is the number of chromophores, ε i ðλÞ is the extinction coefficient at the given wavelength for the i'th substance, and C i is the substance concentration. When the extinction coefficients ε i ðλÞ are known, and the measurements are treated differentially in time, the factor G cancels out and Eq. (3) reduces to E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 4 ; 3 2 6 ; 1 3 9 ΔODðλÞ ¼ X n i¼1 ½ε i ðλÞΔC i dDPFðλÞ; where Δ indicates the time-dependent changes. Standard fNIRS systems generally employ a dual wavelength probing to measure ΔHHb : The above equation needs to be inverted in real scenarios, where concentrations of hemoglobin are derived from empirical measurements of ΔOD: E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 6 ; 6 3 ; 6 0 8 The 2 × 2 matrix inversion yields: E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 7 ; 6 3 ; 5 1 5 Notice that the a priori parameters of the above equation are the extinction coefficients of O 2 Hb and HHb at the two wavelengths ε i ðλÞ and the DPFs at the two wavelengths. Both the geometrical interoptode distance d and the changes in optical densities ΔOD are measured empirically. Equation (7) can be rewritten by extracting from the 2 × 2 inverted matrix one of the two DPFs, as follows: E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 8 ; 6 3 ; 3 1 4 Notice from the above equation that, whereas the absolute value of one of the two DPFs acts by equally modulating the amplitudes of the oscillations of the two hemoglobin forms given ΔODs, the ratio of the two DPFs is involved in the weighted difference of the ΔODs to retrieve the two hemoglobin forms' oscillations. Thus, an error in the ratio of DPFs produces cross talk between the two hemoglobin forms. This dual wavelength DPF ratio concept is generalizable to the multiwavelength spectral dependence of the DPF. Another quantitative parameter when dealing with diffusive media is the EAC often labeled in the literature as μ eff ðλÞ. Here, μ eff ðλÞ expresses the light exponential rate of attenuation caused by the interaction of absorption and diffusion. When d is sufficiently large (above 10 to 15 mm), light attenuation in a diffusive medium is proportional to the factor e −μ eff ðλÞd , where d is the geometrical interoptode distance. With reference to Eq. (1), DPF and μ eff can be linked through the following equation: E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 9 ; 3 2 6 ; 6 9 7 DPFðλÞ ¼ ∂μ eff ðλÞ ∂μ a ðλÞ : In the simplifying assumption of dealing with an infinite homogeneous medium, the EAC can be expressed as E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 1 0 ; 3 2 6 ; 6 3 2 where μ a ðλÞ is the absorption coefficient and μ 0 s ðλÞ is the reduced scattering coefficient. 22 Combining Eqs. (9) and (10), it is possible to obtain E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 1 1 ; 3 2 6 ; 5 6 5 DPFðλÞ The above equation can be rewritten by expressing the DPF as a function of μ eff ðλÞ as E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 1 2 ; 3 2 6 ; 4 9 2 As introduced previously, the variability in DPF as a function of wavelength is particularly important. When considering only two wavelengths, this means assessing the ratio of DPFs. The ratio can then be expressed as E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 1 3 ; 3 2 6 ; 4 0 4 The first of the two equations presented in the above equation is of interest. If, in the first approximation, a linear decay with wavelength in the reduced scattering coefficient is assumed as E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 1 4 ; 3 2 6 ; 2 9 9 where h is a constant representing the intensity of the linear decay (i.e., slope), then the first equation in Eq. (13) can be rewritten as E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 1 5 ; 3 2 6 ; 2 3 5 This means that, if the approximation describing wavelength dependence in reduced scattering is correct, it is possible to estimate the DPF spectral dependence from EAC.

Effective Attenuation Coefficient Computation
EAC values were computed employing the algorithm reported in Ref. 20, where the procedure was described in detail. The algorithm estimates the EAC based on the slope of the log signal as a function of source-detector distance. In fact, using the simplifying assumption that the head can be approximated by a semi-infinite, homogenous medium with zero boundary conditions, 23 the CW light signal IðdÞ at a distance d is linked to the EAC (μ eff ) through the following equation: E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 1 6 ; 6 3 ; 7 3 0 where k is a factor that depends on μ eff but does not depend on distance and is affected by source power, detector efficiency, and optode-to-scalp coupling, and μ eff is the EAC defined as E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 1 7 ; 6 3 ; if it is assumed μ a ≪ μ 0 s (as it normally is in the head tissue in the NIR range).
In real recordings, the SNR of the signal can be linked to light intensity impinging on the detector as E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 1 8 ; 6 3 ; 5 8 8 I ∝ SNR 2 : The above equation holds if noise is mainly due to quantum/shot noise [which is true if the amount of light detected is not extremely low and the SNR is computed using frequencies much higher than the physiological signal spectral range (i.e., f > 10 Hz)]. 20 The SNR of the signal is defined as E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 1 9 ; 6 3 ; 4 9 1 SNR 2 ¼ where var is the variance operator, iðtÞ is the electrical signal recorded by the photodetector in the spectral range of interest (f > 10 Hz), and i avg is the average electrical signal. Thus, Eq. (16) can be rewritten as E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 2 0 ; 6 3 ; 3 9 5 Employing SNR in the above equation dampens the recordings problem of having different sensitivities and operation voltages of different detectors. The variability in the optode-to-scalp coupling, generally accounted for through calibration procedures in quantitative NIRS, is instead dampened on uncalibrated recordings by employing many (few tens) source-detector couples (optical channels) for a single EAC computation and through linear regression employing a least-squares approach. The high numerosity of source-detector couples allows reducing the confidence interval of the estimated EAC (slope of the regression), within acceptable ranges (few point percentage, refer to Ref. 20). Notice that this approach works within the assumption of statistical independence between injected light sources power and distance from detectors.

Participants
A total of 21 healthy adults (6 females) between the ages of 25 and 40 years were recruited for the study that compared the developed procedure to a standard DPF evaluation based on calibrated FDMD. A different sample of six healthy adults (two females), in a similar age range (ages between 25 and 37 years), were further recruited to assess the performances of the algorithm when applied to real fNIRS recordings. All participants were Caucasian. The study was conducted in agreement with the principles described in the Declaration of Helsinki and it was approved by the Research Ethics Board of the local university. Informed consent form was signed by all participants before the experiment and they were able to withdraw from it at any time.

Optical Recordings and Experimental Design
Optical data were acquired with a multichannel frequencydomain 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 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: for 830 nm, and refraction index at both wavelengths η ¼ 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)]. 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., 17;233 − 271 ¼ 16;962, 16;962 − 271 ¼ 16;691, and 16;693 − 271 ¼ 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. 24 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. 25 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. 26 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.

Data Analysis
The effect of errors in DPF was estimated through a simulated dual-wavelength estimation of functional hemoglobin oscillations. Event-related O 2 Hb and HHb temporal changes were simulated by convolving an impulse function (depicting stimulation) with the canonical hemodynamic response function. 27 A maximum positive change in O 2 Hb of 3 μM and a negative change in HHb of 1 μM were assumed. These hemoglobin oscillations are quantitative comparable with those obtained in vivo. 28 Here, Δ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 calib and Ph calib to obtain the known phantom EACs and DPFs. The four rows were independently calibrated through the following equations: 29 E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 2 1 ; 3 2 6 ; 6 4 2 ∂ lnðI calib ðλÞd 2 Þ ∂d ¼ dμ eff_phantom ðλÞ; and E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 2 2 ; 3 2 6 ; 5 7 8 where ω is the modulation frequency of the recording system and v is the speed of light in the medium v ¼ c∕η, where η 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). 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. 10 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 amongsubjects 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 ΔODs, movement corrected, 30 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 Δ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 (DPF 830;690 ¼ 6, DPF ratio ¼ 1), the second procedure estimated DPFs for both wavelengths based on Ref.    Figure 3(b) reports a colormap image where color represents the error in the absolute ratio of O 2 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. Figure 4 shows correlation plots and associated Bland-Altman plots 31 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). 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. 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. 10   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 ∼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. 10 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 O 2 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 O 2 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 O 2 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 O 2 Hb of 10.7% (t ¼ 7.24, df ¼ 5, p ¼ 8 × 10 −4 ) and a negative increase in HHb of 1% (t ¼ 0.58, df ¼ 5, p ¼ N:S:) were obtained. Interestingly, when comparing the EAC-based procedure with the equation model, a further positive increase in O 2 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 ¼ N:S:) were obtained.

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 O 2 Hb and HHb modulations from light recordings (Fig. 3). This cross talk can be detrimental for fNIRS since it is well known that O 2 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. 32 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, 33 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 wellknown 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) 34 and multidistance tissue oximeters that measure tissue oxygenation index 19 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 acrosssubject, 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 acrosssubject 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. 10 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 O 2 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 O 2 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 20 that relied on accounting for the variability in detector sensitivity and optodeto-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 ∼700 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, 33 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. 35 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 O 2 Hb and HHb oscillations in the brain. Correct estimation O 2 Hb and HHb functional oscillations increases fNIRS SNR and reduces false positives and false negatives 36 in the evaluation of brain activity with fNIRS.

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.

Disclosures
The authors declare no competing financial and nonfinancial interests.