1.IntroductionUnderstanding the intricate hemodynamic phenomena occurring in the human brain is crucial for advancing our knowledge of brain functions and related pathologies. The vasomotion, Mayer waves, neurovascular coupling, and influence of cardiac and respiratory activities on brain tissue are among the phenomena that underpin brain functionality.1 Previous studies have explored this field using functional magnetic resonance imaging (fMRI) and functional near-infrared spectroscopy (fNIRS) in the continuous-wave (CW) and frequency-domain (FD) modalities. These studies were conducted both on healthy subjects2,3 and in the presence of various pathologies, such as stroke,4 mild to severe cognitive impairment,5 acute and traumatic brain injuries,6,7 autoregulation dysfunctions,8,9 and autism.10 Multiple fNIRS connectivity studies were also performed.11,12 To the best of our knowledge, only two previous studies attempted to investigate resting-state brain oscillations making use of time-domain (TD) fNIRS. In the first study, by Themelis et al.,13 in-vivo measurements on human and piglet heads were conducted using FD and TD fNIRS systems to detect hemodynamic systemic oscillations within the cortical tissue. The authors presented preliminary results indicating that the modulation of the cortical signal at 830 nm due to the cardiac activity can be recorded with both techniques. However, no power spectral density (PSD) analysis was provided. In the second study, Kacprzak et al.14 employed the “moment” approach15 to extract information on the hemodynamic behavior of different regions within the probed tissue. The authors analyzed the 0th (intensity) and 2nd (variance) centralized moments of the recorded photon distributions of time-of-flight (DTOFs), the former mainly reflecting superficial extra-cerebral variations, and the latter mainly reflecting cerebral compartment variations. Then, they derived the PSD of the intensity and variance time series. This study involved two groups of subjects, both measured during rest: a control group consisting of healthy individuals and a second group comprising patients with severe neurovascular disorders. Notably, statistically significant differences between the two groups were observed when comparing the obtained power spectra. In addition, in both studies, no Fourier-domain analysis was performed on the hemodynamic parameters. Although these investigations have shed light on certain aspects of cerebral hemodynamics, the unique strengths of TD fNIRS have yet to be fully explored. Because of the poor signal-to-noise ratio (SNR), most of the TD fNIRS instruments operate at an acquisition rate of . This is typically enough for monitoring task-related cortical hemodynamic responses, but not for other applications such as studying the resting-state oscillations that occupy wider frequency ranges (up to 5 Hz).16,17 Recently developed TD fNIRS devices address this primary limitation.14,18,19 The enhanced SNR achieved by these devices empowers researchers to delve deeper into the investigation of cerebral hemodynamic oscillations, unlocking new possibilities for understanding brain functionality. In addition, the performance of CW fNIRS in terms of depth-selectivity for this application remains largely unexplored. Simulations found in the literature were mainly designed to characterize the sensitivity of fNIRS to cerebral and extra-cerebral hemodynamics at different values of the source-detector distance (SDD), by evaluating the signal response to instantaneous perturbations imposed in the medium.20,21 In our study, we address these limitations by conducting numerical simulations to evaluate the effect of key operational parameters on the performance of both TD and CW fNIRS in detecting periodical oscillations in the probed tissue hemodynamics, aiming at assisting researchers in defining optimal experimental protocols tailored to their specific study objectives. 2.Methods2.1.Forward SimulationsTo generate the synthetic TD fNIRS signals, namely, a dataset of DTOFs, numerical simulations based on a forward diffusion equation model at different SDDs were implemented using home-made software.22 Two distinct geometries, a homogeneous medium and a bilayer medium, were used to mimic the geometric characteristics of the tissue under investigation while assessing different objectives. For both the homogeneous and the bilayer geometry, the instrument response function15 (IRF) was modeled as an ideal Dirac delta function, and to mimic real experimental conditions, noise following the Poisson distribution was added to each DTOF. 2.1.1.Homogeneous mediumFirst, with the homogeneous medium, we studied the impact on the sensitivity of the technique to periodic perturbations of the optical properties of the medium of the main operational parameters: the SDD, the average total number of photons used to generate the DTOF (), the measurement length (), and the sampling rate (). To limit the complexity of our model, we only considered oxygenated hemoglobin () and deoxygenated hemoglobin (HHb) as the main absorbing chromophores in human tissues in the wavelength range of interest, i.e., 600 to 900 nm. As for their baseline concentrations, we assumed the values of 30 and , respectively.23 Sinusoidal perturbations with a fixed amplitude of 1% of the baseline value and fixed frequency of 1 Hz were imposed to create temporal variations in both hemodynamic parameters. The 1 Hz frequency was selected to mimic the effect of the cardiac muscle contraction on the probed tissue, a phenomenon previously observed in several CW fNIRS studies,24,25 and a relative amplitude of 1% was used to represent the minimum amplitude expected for physiological noise components from the previous literature.24,25 To account for the lack of consensus in the literature regarding the relative phase of oscillations of the same origin between chromophores or across different medium domains, a random phase shift was applied to each perturbation.26,27 Then, by exploiting the known absorption spectra of and HHb, the hemodynamic parameters were translated, by Beer’s law, in the absorption coefficient () at two wavelengths, 690 and 830 nm. The reduced scattering coefficient was set according to the empirical approximation of the Mie theory (), with , , and . Finally, the solution of the diffusion equation for a 5 cm-thick, laterally infinite homogeneous slab with a 1.4 refractive index was used to calculate the dataset of DTOFs28 in the range 0 to 5 ns, with a time resolution of 8 ps. First, three values of were considered: , , and photons (ph), while maintaining the other parameters fixed: , (case H_). Then, three values of were considered: 5, 10, and 15 min, while maintaining , and (case H_). Finally, three values of were considered: 5, 10, and 20 Hz, while maintaining , and (case H_). In these simulations, two SDDs, 1 and 4 cm, were used. An additional simulation was performed using SSD from 1 to 6 cm, with steps of 1 cm, with the other parameters fixed (, , ) to better explore the effect of SSD (case H_SDD). In Table 1, all of the parameters used in each simulation case for the homogeneous medium (H_-H_SDD) are summarized. Table 1Parameters used in each simulation case for the homogeneous medium.
2.1.2.Bilayer mediumThe bilayer medium allowed for the evaluation and comparison of approaches for depth-selectivity enhancement, both in CW and TD fNIRS. This comparison aimed to assess the ability of the techniques to localize perturbations originating at different depths and effectively reject superficial extra-cerebral contributions. The bilayer medium is composed of two stacked cylinders with a radius of with a total thickness of 5 cm and a refractive index of 1.4. The superficial layer (UP) had a thickness of 1 cm, and the deeper layer (DW) had a thickness of 4 cm. Both layers had identical baseline hemoglobin concentrations and optical properties so the two layers are only identified by different and independent perturbations of the hemoglobin concentrations in time. The time duration of measurements and the sampling frequency were fixed based on the results obtained for the homogeneous medium (Sec. 2.1.1). The average total number of photon counts was varied as before by considering three values of , , and and two SDDs of 1 and 4 cm. Preliminary simulations investigated the effect of a single oscillatory component imposed in either the UP layer (case B_UP) or the DW layer (case B_DW) while keeping the concentrations in the other layer constant at baseline values. The same perturbation as in the homogeneous medium was used: an oscillation of both and HHb concentrations with a frequency of 1 Hz and an amplitude of 1% of the baseline value. In addition, we focused on evaluating the effects generated through the superimposition of multiple Fourier contributions simultaneously originating in the UP and DW layers. Four cases were considered:
These simulations were repeated with the applied perturbations in-phase or out-of-phase. The solution of the diffusion equation for the bilayer geometry was used to calculate the dataset of DTOFs28 in the range 0 to 5 ns, with a time resolution of 20 ps to limit the computational time. Table 2 summarizes the parameters used in each simulation case when the bilayer medium was used (B_UP-B_AmpFreq). Table 2Parameters used in each simulation case for the bilayer medium.
2.2.Data Analysis2.2.1.Estimate of optical and hemodynamic parametersFor the simulations carried out using the homogeneous medium, the DTOFs were numerically time integrated to retrieve the measured light intensity, represented as the number of total photons detected within the measurement window, which can be considered an equivalent of the CW fNIRS signal. When the bilayer geometry was considered, a time-windowing approach was used to exploit the depth-selective information encoded in the DTOFs. For our analysis, we used 10 time windows (gates), each with a width of 500 ps, covering the time window between 0 and 5 ns (full width of the largest curve) following the IRF temporal position. Successively, the diffusion equation solution for a bilayer medium28 was used as a model for photon migration in a home-made software22 for non-linear fitting (based on the Levenberg–Marquardt minimization algorithm29). After calculating the baseline optical coefficients and using a homogeneous model fit, the time-dependent mean partial pathlength (TMPP) method30 was used to estimate the variations from the baseline value from each DTOF, to extract the absorption coefficients () in the two layers, or 2, of the medium at every time point. The TMPP method, presented by Zucchelli et al.,30 is based on a time gating of the DTOF, which is used to estimate the time-dependent mean pathlength traveled by photons in each domain of the medium with high accuracy. Once is known, the variations can be retrieved by inverting the following formulation of the time-resolved reflectance curve: . Finally, Beer’s law was used to derive the hemodynamic parameters from the absorption coefficient.28 To enable a direct comparison between TD and CW fNIRS data, the DTOFs were again integrated over the measurement window. The variations of hemoglobin concentrations for the CW fNIRS configuration were then calculated by means of the modified Lambert–Beer law.28,31 Increases or decreases in hemoglobin concentration were expressed relative to a baseline, calculated as the average intensity over the entire measurement. The depth-selectivity in CW fNIRS is generally implemented using multiple acquisition points at different SDDs and exploiting the relationship between the SDD and the photon penetration depth.32 Signals acquired with an SDD of 1 and 4 cm were thus compared in this phase. 2.2.2.Frequency-domain analysisThe PSD of the intensity and hemodynamic signals were calculated using a uniform algorithm across all data, implemented in Matlab (Release 2022a, The MathWorks Inc., Natick, Massachusetts, United States). First, the signals were detrended using a third-order polynomial fit to remove the slow component.10,26 Then, the periodogram estimation of the PSD was performed by employing the Cooley–Tukey fast Fourier transform algorithm.33,34 Because our aim was to evaluate the PSD with the maximum frequency resolution without any variance reduction and all the data series met the stationarity requirements, no windowing was applied. In Fourier-domain analysis, evaluating the significance of spectral peaks is crucial: to discern peaks attributable to underlying oscillatory phenomena from Poisson noise background, an ad-hoc contrast function was introduced. The contrast function was defined to compare the PSD amplitude at frequency to the average amplitude of the PSD in a noise-only interval, In Sec. 3, to quantify the contrast value efficiently, it is expressed in decibels, as needed Due to the characteristics of our synthetic signals, the PSD average noise was typically calculated by averaging the PSD values for frequencies higher than 5 Hz. However, for the case , a modification in the definition of was necessary, as the sampling rate did not permit reconstruction of the PSD for frequencies higher than 2.5 Hz. In this case, was calculated as the average value of the PSD in the range to avoid interference from the second harmonic perturbation at 2 Hz, if present. 3.Results3.1.Homogeneous MediumIn this section, the effect of the main operational parameters on the detected signal is investigated. The results are shown in terms of PSD of the intensity signal, equivalent to the CW fNIRS signal. Therefore, the obtained results are still valid for both TD and CW fNIRS. 3.1.1.Effect of the average total number of photon countsFigure 1 shows the effect of varying on the spectra obtained in case H_. Panels (a) and (b) show the obtained PSD for the SSD of 1 and 4 cm, respectively. Data are shown for a single wavelength (690 nm) as the results of the other wavelengths (830 nm) are similar. The peak at 1 Hz is clearly visible for every value of (different colors). Increasing , higher amplitudes are obtained for the whole spectrum. Figures 1(c)–1(h) show the amplitude of the spectral peak at 1 Hz, the amplitude of the average noise, and the 1 Hz peak contrast as a function of in logarithmic–logarithmic plots, for SSD = 1 cm (4 cm) and both wavelengths. A power law fit is applied to each data series, always returning an value higher than 0.999. It is interesting to note that a second harmonic of the 1 Hz peak (at 2 Hz) is visible for the longer SDD and the higher values of [Fig. 1(b), in green and purple] but remains absent in the short SSD signals [Fig. 1(a)]. Comparing the difference in amplitude between the first and second harmonic (where visible) and assuming this difference to be consistent for all measurements, we confirmed that the second harmonic is covered by noise for all measurements in which it is not visible in the PSD. The results obtained by increasing can be summarized as follows: (i) the 1 Hz peak amplitude increases with ; (ii) the average noise amplitude increases with ; and (iii) the 1 Hz peak contrast value increases with . 3.1.2.Effect of the measurement timeFigure 2 shows the effect of varying the on the spectra obtained in case H_. Panels (a) and (b) show the PSD obtained for the SSD of 1 and 4 cm, respectively. Data are shown for a single wavelength (690 nm) as the results for the other wavelength (830 nm) are similar. The smaller box on the right of each graph represents a zoom on a higher-frequency region, showing the noise level in the three signals. Figures 2(c)–2(h) show the amplitude of the spectral peak at 1 Hz, the amplitude of the average noise, and the 1 Hz peak contrast as a function of , for SSD = 1 cm (4 cm) and both wavelengths. A power law fit was applied to data shown in Figs. 2(c), 2(d), 2(g) and 2(h), resulting in values higher than 0.999. In Figs. 2(e) and 2(f), horizontal lines indicating the 3% variation around the data mean value were added to highlight that the average noise amplitude is invariant with respect to . The results obtained by increasing can be summarized as follows: (i) the 1 Hz peak amplitude increases with ; (ii) the average noise amplitude is invariant with respect to ; and (iii) the 1 Hz peak contrast value increases with . 3.1.3.Effect of the sampling rateFigure 3 shows the effect of varying the on the spectra obtained in case H_. Panels (a) and (b) show the PSD obtained for the SSD of 1 and 4 cm, respectively. Data are shown for a single wavelength (690 nm) as the results for the other wavelength (830 nm) are similar. The smaller box on the right of each graph represents a zoom on a higher-frequency region, showing the noise level in the three signals. Figures 3(c)–3(h) show the amplitude of the spectral peak at 1 Hz, the amplitude of the average noise, and the 1 Hz peak contrast as a function of , for SSD = 1 cm (4 cm) and both wavelengths. In Figs. 3(c) and 3(d), horizontal lines indicating the 3% variation around the data mean value were added to highlight that the 1 Hz peak amplitude is invariant with respect to . A power law fit was applied to the data shown in Figs. 3(e)–3(h), resulting in values higher than 0.99. The results obtained by increasing can be summarized as follows: (i) the 1 Hz peak amplitude is invariant with respect to ; (ii) the average noise amplitude decreases with ; and (iii) the 1 Hz peak contrast value increases with . 3.1.4.Comparison among different SSDsComparing the results obtained for the two SDD values in Secs. 3.1.1–3.1.3, an increase in the peak amplitude and contrast can be observed when increasing the SDD (Figs. 1–3). This effect is mainly due to the increase in the length of the average optical path traveled by the detected photons within the medium.28 To better understand the effect of the SDD on the retrieved PSDs, we performed additional simulations varying SDD (H_SDD). The relative results are shown in Fig. 4. Panels (a)–(c) respectively show the amplitude of the spectral peak at 1 Hz, the amplitude of the average noise, and the 1 Hz peak contrast as a function of SDD, for the two wavelengths. In Fig. 4(b), horizontal lines indicating the 3% variation around the data mean value were added to highlight that the average noise amplitude is invariant with respect to SDD. A power law fit was applied to the data shown in Figs. 4(a) and 4(c), resulting in values higher than 0.999. The results obtained by increasing SDD can be summarized as follows: (i) the 1 Hz peak amplitude increases approximatively with the square of SDD; (ii) the average noise amplitude is invariant with respect to SDD; and (iii) the 1 Hz peak contrast value increases approximatively with the square of SDD. 3.2.Bilayer MediumIn this section, we investigate and compare the depth-selectivity enhancement approaches employed in CW and TD fNIRS. The time duration of measurements and the sampling frequency were selected based on the results presented in Sec. 3.1 ( and , respectively). The average total number of photon counts was varied as in Sec. 3.1.1, although most results presented here are referred to cases , which presents the clearest contrasts, and , the closest to real TD fNIRS in-vivo measurements with an SDD of 4 cm and a 20 Hz sampling rate, performed with a device previously developed at our department.18 3.2.1.Time-windowing of the DTOFsA time-windowing of the DTOFs was applied using 10 gates of 500 ps width, covering a 5 ns time window starting from the IRF temporal position (Fig. 5). Figure 6 provides a comprehensive summary of the results obtained for cases B_UP, B_DW, and B_UPDW, showcasing the contrast of the 1 Hz peak in the PSD of the number of photon counts within each gate, as a function of the gate number, for , both SDDs in Figs. 6(a) and 6(b), respectively, and a single wavelength (690 nm). The results are denoted with green crosses for case B_UP, with yellow triangles for case B_DW, and with purple circles for case B_UPDW. The results obtained for the CW fNIRS signals are also shown on the right. The contrast values are expressed in decibels, and the red horizontal line represents a peak significance threshold of 15 dB, tailored specifically for this analysis, which is the maximum contrast in the noise frequency range (Sec. 2.2.2). Remarkably, a significant peak at 1 Hz was observed in the PSD of each gate that contained a not null portion of the reflectance curve (DTOFs retrieved for the short SDD decayed to zero around 3 ns, leaving the last four gates with only Poisson noise, see Fig. 5) for case B_UP (green crosses) and case B_UPDW (purple circles). This is a notable result because a significant peak is found even in gates covering the very end of the DTOF tail [from around gate #5 in Fig. 6(a) and gate #10 in Fig. 6(b)], in which the average number of photon counts is lower than 20 (Fig. 5). The only exception is represented by the first gate in case B_DW (yellow triangles), for both SDDs. In addition, the curves exhibit differences in peak and barycenter position and slope in the last gates. Another remarkable result is that the PSD calculated for the CW fNIRS signal always shows a substantial peak at the perturbation frequency. Of particular importance is the sensitivity shown by the CW fNIRS signal recorded at 1 cm to the perturbation imposed in depth [Fig. 6(a), case B_DW, yellow triangles], where, even if less pronounced, the peak at 1 Hz remains indeed significant. 3.2.2.Estimated hemodynamic parametersFor TD fNIRS, the TMPP method was applied to retrieve the hemodynamic parameters in the two layers of the tissue. The obtained results are shown in Fig. 7 for case B_UP, in Fig. 8 for case B_DW, and in Fig. 9 for case B_UPDW. The results are presented for and both SDDs. In case B_UP (Fig. 7), oscillations in the hemodynamic parameters of both layers were accurately reconstructed for each SDD as the PSD of both hemoglobin species exhibited a visible peak at the imposed perturbation frequency (1 Hz) in the superficial layer [UP, Figs. 7(a) and 7(b)], but not in the deeper layer [DW, Figs. 7(c) and 7(d)], as expected. In case B_DW (Fig. 8), similar results were observed for the long SDD [4 cm, Figs. 8(b) and 8(d)], with a visible peak at the perturbation frequency in the PSD of the hemoglobin concentrations calculated in the DW layer [Fig. 8(d)], but not in the UP layer [Fig. 8(b)], as expected. However, when the short SDD [1 cm, Figs. 8(a) and 8(c)] was used, different results were obtained. In this case, although the UP layer static behavior was correctly reconstructed [Fig. 8(a)], the perturbation peak in the PSD was no longer visible in the DW layer [Fig. 8(c)], where it should be present. The same observations regarding the lower sensitivity of short SDD measurements to variations in the deeper DW layer can also be done for case B_UPDW [Fig. 9(c)]. When considering SDD = 4 cm, Figs. 9(b) and 9(d) demonstrate a flawless reconstruction of the PSD for both species and both layers. The superimposition of the two contributions had no discernible effects, evident from the perfect overlap between PSDs in Figs. 7(b) and 9(b) [Figs. 8(d) and 9(d)]. Comparing the above PSDs, the amplitude of the 1 Hz peak is indeed equal for both species, except for the variability due to the background noise. Analogous results were obtained by imposing oscillations with different frequencies and/or amplitudes in the two layers of the medium (cases B_Amp, B_Freq, and B_AmpFreq) while considering them either in-phase or out-of-phase (data not shown). For CW fNIRS, the results are shown in Fig. 10 for case B_UP, Fig. 11 for case B_DW, and Fig. 12 for case B_UPDW. The results are shown for to account for the higher intensities usually characterizing CW fNIRS signals compared with TD fNIRS ones. In all of the above cases, a visible peak is present at the perturbation frequency for both species and both SDDs. This is particularly interesting in case B_DW (Fig. 11) when the perturbation is imposed solely in the DW layer. In Fig. 11, the CW fNIRS signals recorded at both SDDs show sensitivity to the perturbation imposed in the DW layer of the medium for both species. In Sec. 3.2.1, the sensitivity shown by the signal recorded at 1 cm [Fig. 11(a)] is a notable and unexpected result. 4.DiscussionIn this work, we conducted a simulation study to assess the viability of employing TD fNIRS for the examination of cerebral hemodynamic oscillations during resting-states in humans. Previous research has delved into this area using fMRI, CW fNIRS, and FD fNIRS techniques, but the adoption of TD fNIRS has been limited by lower SNRs. As highlighted in Sec. 1, to the best of our knowledge, only two attempts of using TD fNIRS in this field of investigation are reported in the literature,13,14 and no studies that systematically evaluate the sensitivity of TD fNIRS to periodical perturbations in the optical properties of the probed medium exist. The aim of our study was to fill an existing literature gap, providing some guidelines for researchers to design optimal experimental protocols when investigating cerebral hemodynamics in humans using fNIRS, contributing to addressing the crucial challenge of standardizing fNIRS experimental procedures. To this purpose, we performed numerical simulations of TD and CW fNIRS acquisitions, considering two different geometries of the probed medium using homogeneous and bilayer models. 4.1.Homogeneous MediumExploiting the homogeneous model, the influence of the number of detected photons per DTOF (), the overall measurement length (), sampling frequency (), and the used SDD on the recorded signal PSD were tested. The sensitivity of fNIRS to sinusoidal perturbations in the hemodynamic parameters of the medium was expressed in terms of an ad-hoc defined spectral peak contrast, calculated from the PSD of the intensity signal. Due to this choice, the reported results are valid for both TD and CW fNIRS. The results show that the contrast value obtained with the two techniques increases proportionally with , , and and approximatively with the square of SDD, when the parameters are varied independently. The same relationships were found to be valid when using different values of perturbation amplitude (0.5% and 5%, data not shown). These findings, together with the successive evaluation of different data analysis pipelines, provide insights into the optimization of fNIRS experimental designs. Depending on the specifications of the used instrument and the characteristics of the investigated phenomena, different combinations of acquisition parameters can be selected to maximize TD and CW fNIRS sensitivity. As an example, when working with an instrument limited in power, for which the achievable is not sufficient to obtain a clear spectrum, this limitation can be compensated for by increasing the value of . Regarding , any increase in this parameter will result in a decrease in the maximum number of photons detected () when the maximum allowable light intensity is used. However, in TD fNIRS, sampling at higher frequencies is particularly useful when using depth-selectivity algorithms based on the time-windowing of the DTOF. The use of corrective methods such as the TMPP reduces the effect due to changes in on the PSD of the derived signal, while the effects of and remain unchanged (data not shown). When the desired effect is to maximize the contrast in the PSD of and HHb of the two layers as calculated using the TMPP method, increasing is therefore an effective method even at the expense of the detected number of photons. Once the minimum threshold for the coefficient of variation35 and the corresponding one for are set,18 the best choice to optimize the quality of the obtained PSD is to use the maximum that guarantees staying above it. 4.2.Bilayer MediumThe bilayer model was chosen for investigating the potentiality of the depth-selectivity enhancement approaches employed in TD fNIRS compared with the multi-distance CW fNIRS approach. First, a time-windowing of the DTOFs was used to evaluate the depth-selectivity of TD fNIRS measurements.15 This approach is based on the evidence that, in TD fNIRS, photons detected with a longer time of flight show a higher probability of having traveled deeper into the medium compared with those detected shortly after the pulse injection.15 Thus, it can be used to isolate superficial and deeper contributions using a single-channel measurement. For this reason, it represents an alternative to the short-SDD regression method used in CW fNIRS, mitigating the challenges associated with the use of multiple SSDs.21 To evaluate the potentiality of the time-windowing approach, the signals obtained in cases B_UP, B_DW, and B_UPDW were analyzed by slicing the DTOF using 10 gates, each with a width of 500 ps, covering the time window between 0 and 5 ns following the IRF temporal position. The PSDs of the number of photons detected within each gate were calculated and compared between the two SDDs and among the three cases. The results indicate that, irrespective of the SDD, photons detected within the first 500 ps are affected by the variation of the optical properties of the medium only when they occur in the superficial layer (Fig. 6, gate 1), confirming that these photons do not penetrate at a depth greater than 1 cm in the tissue (i.e., a typical thickness of the scalp). As a result, the shape of the DTOF in this time region is not influenced by the in-depth behavior of the medium. Consequently, the fNIRS signal collected within the first 500 ps from the pulse delivery contains information solely on the oscillatory pattern in the UP layer. By considering narrower gate widths and non-ideal IRF,18 we also verified that the information content relative to the UP layer remains confined within a time window of to 600 ps from the IRF temporal position (data not shown). Then, the simple gating approach applied to TD fNIRS measurements can effectively identify oscillations originating in the superficial layer. This result opens the way for the implementation of a regression method for the TD fNIRS signal that mirrors the approach used in the CW modality, with the additional advance of removing the requirement for multiple SDD acquisitions.21 Another remarkable result obtained with these simulations is the high sensitivity of the TD fNIRS signal to periodical perturbations in the optical properties of the probed medium. This is demonstrated by the presence of a significant spectral peak in PSDs associated with time windows covering the very end of the DTOF tail, where the detected photons are on the order of some tens or less. After this first analysis phase, two different methods were applied to derive the hemodynamic parameters in the tissue. For TD fNIRS, the TMPP method was used. When the perturbed layer was the superficial one (case B_UP), the PSD of the hemodynamic parameters of both layers was accurately reconstructed for each SDD (1 and 4 cm), as the PSD of both hemoglobin species exhibited a visible peak at the imposed perturbation frequency solely in the UP layer (Fig. 7). Similar results were observed for the long SDD when the perturbed layer was the deeper one (case B_DW), with a visible peak at the perturbation frequency in the PSD of the hemoglobin concentrations in the DW layer, but not in the UP layer, as expected [Figs. 8(b) and 8(d)]. Conversely, different results were obtained for the signal recorded at the short SDD in the same case. Although the UP layer static behavior was correctly reconstructed, the perturbation peak in the PSD was no longer visible in the DW layer, where it should have been present [Figs. 8(a) and 8(c)]. The absence of this peak is attributed to the reduced number of photons that reached the deeper region of the medium, known as late photons. Despite the number of photons used to construct the DTOFs being the same for both SDDs, the distribution shape differs, being narrower for the shorter SDD (Fig. 5). Consequently, although the absolute number of late photons increases at shorter SDD,36 the ratio between the late and early photons, which traveled in the more superficial region of the medium and hence occupy the first portion of the DTOF, is much lower for the short compared with the long SDD.36 This disparity in photon distribution reduces the TD fNIRS signal sensitivity to the deeper layer of the medium relative to the superficial one, limiting its ability to detect the perturbation applied in case B_DW. Additional simulations were performed by imposing a perturbation of reduced amplitude (0.75%, 0.5%, 0.25%, and 0.10% of the baseline values) only in the deeper layer. The obtained results (not shown) were analogous to those discussed in Sec. 3.2.1, even for the lower value of (). The above results demonstrated that, with an adequate SDD, the use of the TMPP method allows us to accurately retrieve the hemoglobin concentration signals from both layers of the medium. Consequently, TD fNIRS has the capability to detect and precisely locate periodical perturbations in the hemodynamic parameters within the probed medium using a single-point measurement approach. A final set of simulations was designed to evaluate the effects generated by the superimposition of multiple Fourier contributions originating in the two layers of the probed medium, having different amplitudes and/or frequencies (B_Amp, B_Freq, and B_AmpFreq). The PSDs obtained for each hemodynamic parameter and layer were compared with the analogous ones obtained when imposing a single oscillatory perturbation in the medium. No significant deviation was found among PSDs obtained for the layer having the same oscillation imposed, independently from the behavior or the other one. As an example, comparing the UP layer in cases B_UP, B_UPDW, and B_Amp, in which the same oscillation is imposed (1% amplitude, 1 Hz frequency), the resulting PSDs are perfectly superimposable (excluding the contribution of noise) at the frequency of interest. The same result is obtained when comparing the DW layer in cases B_DW, B_UPDW, and B_Freq. This result indicates that the TMPP method can be effectively used to identify the distinct oscillatory components in each layer, without any interference or cross-contamination between them. Similar to what was previously done, the above simulations were also repeated, introducing the convolution with a real IRF.18 As before, the resulting PSD showed no significant differences, confirming the robustness of the TMPP method. However, the TMPP approach is limited by the need for an a priori accurate description of the probed medium geometry. Moreover, as the model complexity increases, the computational load escalates rapidly, necessitating simplifications when dealing with measurements on multi-layered structures such as the human brain. Inaccuracies in the medium description (typically the unknown thickness of the UP layer) introduce errors in the calculation of the hemodynamic parameters and restrict the method’s capability to precisely discern the hemodynamic behavior in the different regions of the medium. In this context, exploring the depth-selectivity of the TD fNIRS signal using the time-windowing approach could represent an interesting alternative. 4.3.CW and TD fNIRSTo allow for a direct comparison between TD and CW fNIRS data, DTOFs were numerically integrated over time to obtain the measured light intensity. For CW fNIRS, the hemoglobin concentrations were then recovered from these signals by exploiting the modified Lambert–Beer law. Different from TD, CW fNIRS does not allow for the separation of the information about the oscillations associated with the two layers of tissue using a single SSD. Depth selectivity in this configuration relies on the comparison of the different signals measured simultaneously at different SDDs.37,38 In CW fNIRS studies, one or more short SDD channels are used as a regressor for signals recorded at longer SDDs to correct the latter for the superficial contributions. This approach does not consider the variation in fNIRS signal sensitivity to the different layers of the probed medium when the SDD is varied,21 and the conventional assumption underlying the multi-distance approach is that the penetration depth of photons detected at a given SDD in CW modality is limited to around a half of the SDD value.21,37 However, this assumption is challenged by the results of our study, in which the CW fNIRS signal recorded at 1 cm is affected by perturbations imposed in the DW layer of the medium. Significant peaks attributed to perturbations introduced in the deeper layer are indeed observable in the PSD of both the intensity and hemoglobin concentration signals [Fig. 11(a)]. This effect is probably due to the high sensitivity of fNIRS signals to periodical fluctuations even in the case of very few photons being affected. Owing to this unexpected sensitivity of the CW fNIRS signal at 1 cm to the perturbation in the DW layer, the existing short-distance correction methodologies when a Fourier-domain analysis of the signal is performed are worth re-evaluating. However, the relatively low amplitude of the 1 Hz peak in Fig. 11(a), coupled with the simplifications applied to the imposed perturbation, the measurement system, and the tissue model, suggests that the higher noise levels expected in real-world applications should prevent a significant impact of this effect on actual in-vivo measurements. 4.4.Limitations and Future DevelopmentsThe main limitations of this study lie in the idealized nature of the measurement system and of the perturbations applied and in the simplifications imposed on the geometrical model. In our simulations, the measurement system was considered ideal, generating a delta-shaped IRF and null background noise. Although the broadening of the IRF demonstrated to have no significant effects on the obtained results (as discussed above), the effect of the background noise was not investigated. In real applications, this could lead to lower SNRs in the obtained PSDs, limiting the ability of the technique to detect low-amplitude oscillations. Moreover, the sinusoidal perturbations were modeled with a constant amplitude and frequency, whereas throughout measurements in in-vivo applications, the perturbations of hemoglobin concentrations are expected to exhibit some level of variability due to the natural fluctuations of underlying physiological phenomena. The presence of frequency and phase variations may result in spectral peaks broadening or bifurcation, alongside a concurrent reduction in peak amplitude while maintaining a constant signal energy. These effects introduce complexity into the precise identification of individual spectral peaks and the subsequent classification based on the underlying phenomena. Nonetheless, the limitations arising from the idealized nature of the measurement system and perturbations employed in our study do not invalidate the fundamental findings and conclusions presented thus far. Specifically, our findings on the optimal selection of the main operational parameters and the technique’s efficacy in localizing the detected perturbations within the sample layers remain valid. Finally, we considered the probed medium either as homogeneous or composed of two homogeneous layers, which may not fully represent the complexity of real biological tissues. A multilayer model with a scalp layer, a bone layer, a cerebrospinal fluid layer, a vascular layer on the cerebral cortex, and a brain tissue layer would be necessary to represent more realistically the probed volume. Consequently, the sensitivity of TD fNIRS to fluctuations occurring in vivo could be lower than the one observed in this study. As part of future work, tests will be conducted to assess the performance of TD fNIRS in in-vivo applications. In addition, further simulations will be employed to investigate the impact of inaccuracies in the model used to describe the probed medium when applying the TMPP correction. 5.ConclusionsIn summary, this paper aimed at contributing to the emerging field of the study of brain hemodynamic oscillations through fNIRS. In particular, we focused on the oscillations during resting-state acquisitions with TD and CW fNIRS. Through a series of numerical simulations, we provided valuable guidelines for optimizing four main operational parameters (average photon count rate, measurement duration, sampling frequency, and SSD) to maximize the sensitivity of fNIRS to periodical perturbations of the medium optical properties and for defining the best experimental protocol for detecting resting-state cerebral hemodynamic oscillations. By leveraging the advantages of TD fNIRS, we demonstrated that this technique allows for the detection and depth-localization of periodical fluctuations occurring in the concentrations of and HHb within the probed medium using a single-point acquisition, offering an alternative to multi-distance CW fNIRS setups. Our results demonstrate that the time-windowing of the DTOF is a suitable method for detecting the presence of oscillatory components in the signal, even in the presence of very few detected photons, together with the ability to isolate the behavior of the more superficial layer of the medium. Therefore, the implementation of a method that makes use of the first portion of the DTOF to correct the information encoded in the late photons could allow for obtaining a depth-selective description of the medium hemodynamics in the Fourier domain, even when its geometry is unknown or the photon count rate is insufficient to accurately retrieve the hemodynamic parameters. Moreover, our findings confirm that the use of the TMPP method to retrieve the and HHb concentrations allows us to correctly reconstruct the signals coming from the two layers of the medium with good sensitivity and no cross-talk. Code and Data AvailabilityData underlying the results presented in this paper are not publicly available but may be obtained from the corresponding authors upon reasonable request. AcknowledgmentsThe authors acknowledge Felipe Orihuela-Espina for the fruitful discussions regarding the definition of PSD peak recognition criteria. In addition, we acknowledge the contribution of data previously presented in the proceeding “Assessing TD fNIRS capability to detect hemodynamic oscillations in cerebral cortex,” by L. Contini, R. Re, D. Contini, A. Torricelli, and L. Spinelli, as part of the European Conferences on Biomedical Optics 2023 (DOI: https://doi.org/10.1117/12.2670919). This research was partially funded by the HORIZON EUROPE European Innovation Council (Grant No. 101099093) as part of the project “Preterm Brain-Oxygenation and Metabolic EU-Sensing: Feed the Brain” (Prometeus); by the Ministry of University and Research, “Neuromuscular impairment in aging: a longitudinal study of structural and functional mechanistic bases of age-related alterations (Trajector-AGE),” PRIN 2020 (Grant No. 2020477RW5PRIN); by the NextGenerationEU, National Recovery and Resilience Plan, Age-IT, PE0000015 (Grant No. DD 1557 11.10.2022); and by the NextGenerationEU, National Recovery and Resilience Plan, PRIN 2022 (Grant No. 20227EPKW2). ReferencesM. Li et al.,
“The decoupling between hemodynamic parameters and neural activity implies a complex origin of spontaneous brain oscillations,”
Front. Comput. Neurosci., 17 1214793 https://doi.org/10.3389/fncom.2023.1214793 1662-5188
(2023).
Google Scholar
M. L. Schroeter, O. Schmiedel and D. Y. Von Cramon,
“Spontaneous low-frequency oscillations decline in the aging brain,”
J. Cereb. Blood Flow Metab., 24
(10), 1183
–1191 https://doi.org/10.1097/01.WCB.0000135231.90164.40
(2004).
Google Scholar
Y. Tong, K. P. Lindsey and B. D. Frederick,
“Partitioning of physiological noise signals in the brain with concurrent near-infrared spectroscopy and fMRI,”
J. Cereb. Blood Flow Metab., 31
(12), 2352
–2362 https://doi.org/10.1038/jcbfm.2011.100
(2011).
Google Scholar
S. Muehlschlegel et al.,
“Feasibility of NIRS in the neurointensive care unit: a pilot study in stroke using physiological oscillations,”
Neurocrit. Care, 11
(2), 288
–295 https://doi.org/10.1007/s12028-009-9254-4
(2009).
Google Scholar
D. Yang and K.-S. Hong,
“Quantitative assessment of resting-state for mild cognitive impairment detection: a functional near-infrared spectroscopy and deep learning approach,”
J. Alzheimer’s Dis., 80
(2), 647
–663 https://doi.org/10.3233/JAD-201163
(2021).
Google Scholar
M. H. Othman et al.,
“Resting-state NIRS–EEG in unresponsive patients with acute brain injury: a proof-of-concept study,”
Neurocrit. Care, 34
(1), 31
–44 https://doi.org/10.1007/s12028-020-00971-x
(2021).
Google Scholar
D. Highton et al.,
“Monitoring cerebral autoregulation after brain injury: multimodal assessment of cerebral slow-wave oscillations using near-infrared spectroscopy,”
Anesth. Analg., 121
(1), 198
–205 https://doi.org/10.1213/ANE.0000000000000790
(2015).
Google Scholar
A. V. Andersen et al.,
“Assessing low-frequency oscillations in cerebrovascular diseases and related conditions with near-infrared spectroscopy: a plausible method for evaluating cerebral autoregulation?,”
Neurophotonics, 5
(3), 030901 https://doi.org/10.1117/1.NPh.5.3.030901
(2018).
Google Scholar
J. Bindra et al.,
“Non-invasive monitoring of dynamic cerebrovascular autoregulation using near infrared spectroscopy and the finometer photoplethysmograph,”
Neurocrit. Care, 24
(3), 442
–447 https://doi.org/10.1007/s12028-015-0200-3
(2016).
Google Scholar
H. Cheng et al.,
“Power spectrum of spontaneous cerebral homodynamic oscillation shows a distinct pattern in autism spectrum disorder,”
Biomed. Opt. Express, 10
(3), 1383 https://doi.org/10.1364/BOE.10.001383 BOEICL 2156-7085
(2019).
Google Scholar
S. Sasai et al.,
“Frequency-specific functional connectivity in the brain during resting state revealed by NIRS,”
NeuroImage, 56
(1), 252
–257 https://doi.org/10.1016/j.neuroimage.2010.12.075 NEIMEF 1053-8119
(2011).
Google Scholar
A. Sassaroli et al.,
“Low-frequency spontaneous oscillations of cerebral hemodynamics investigated with near-infrared spectroscopy: a review,”
IEEE J. Sel. Top. Quantum Electron., 18
(4), 1478
–1492 https://doi.org/10.1109/JSTQE.2012.2183581 IJSQEN 1077-260X
(2012).
Google Scholar
G. Themelis et al.,
“Depth of arterial oscillation resolved with NIRS time and frequency domain,”
in Biomed. Top. Meet.,
(2004). https://doi.org/10.1364/BIO.2004.WF2 Google Scholar
M. Kacprzak et al.,
“Frequency analysis of oscillations in cerebral hemodynamics measured by time domain near infrared spectroscopy,”
Biomed. Opt. Express, 10
(2), 761 https://doi.org/10.1364/BOE.10.000761 BOEICL 2156-7085
(2019).
Google Scholar
A. Torricelli et al.,
“Time domain functional NIRS imaging for human brain mapping,”
NeuroImage, 85 28
–50 https://doi.org/10.1016/j.neuroimage.2013.05.106 NEIMEF 1053-8119
(2014).
Google Scholar
S. Cinciute,
“Translating the hemodynamic response: why focused interdisciplinary integration should matter for the future of functional neuroimaging,”
PeerJ, 7
(3), e6621 https://doi.org/10.7717/peerj.6621
(2019).
Google Scholar
Y. Yamada, H. Suzuki and Y. Yamashita,
“Time-domain near-infrared spectroscopy and imaging: a review,”
Appl. Sci., 9
(6), 1127 https://doi.org/10.3390/app9061127
(2019).
Google Scholar
R. Re et al.,
“Reliable fast (20 Hz) acquisition rate by a TD fNIRS device: brain resting-state oscillation studies,”
Sensors, 23
(1), 196 https://doi.org/10.3390/s23010196 SNSRES 0746-9462
(2023).
Google Scholar
H. Y. Ban et al.,
“Kernel flow: a high channel count scalable TD-fNIRS system,”
Proc. SPIE, 11663 116630B https://doi.org/10.1117/12.2582888 PSISDG 0277-786X
(2021).
Google Scholar
J. Selb et al.,
“Sensitivity of near-infrared spectroscopy and diffuse correlation spectroscopy to brain hemodynamics: simulations and experimental findings during hypercapnia,”
Neurophotonics, 1
(1), 015005 https://doi.org/10.1117/1.NPh.1.1.015005
(2014).
Google Scholar
S. Brigadoi and R. J. Cooper,
“How short is short? Optimum source–detector distance for short-separation channels in functional near-infrared spectroscopy,”
Neurophotonics, 2
(2), 025005 https://doi.org/10.1117/1.NPh.2.2.025005
(2015).
Google Scholar
R. Cubeddu et al.,
“Experimental test of theoretical models for time-resolved reflectance,”
Med. Phys., 23
(9), 1625
–1633 https://doi.org/10.1118/1.597739 MPHYA6 0094-2405
(1996).
Google Scholar
L. Gagnon et al.,
“Double-layer estimation of intra- and extracerebral hemoglobin concentration with a time-resolved system,”
J. Biomed. Opt., 13
(5), 054019 https://doi.org/10.1117/1.2982524 JBOPFO 1083-3668
(2008).
Google Scholar
F. Scarpa et al.,
“A reference-channel based methodology to improve estimation of event-related hemodynamic response from fNIRS measurements,”
NeuroImage, 72 106
–119 https://doi.org/10.1016/j.neuroimage.2013.01.021 NEIMEF 1053-8119
(2013).
Google Scholar
V. Bonomini et al.,
“Linear regression models and k-means clustering for statistical analysis of fNIRS data,”
Biomed. Opt. Express, 6
(2), 615 https://doi.org/10.1364/BOE.6.000615 BOEICL 2156-7085
(2015).
Google Scholar
H. Obrig et al.,
“Spontaneous low frequency oscillations of cerebral hemodynamics and metabolism in human adults,”
NeuroImage, 12
(6), 623
–639 https://doi.org/10.1006/nimg.2000.0657 NEIMEF 1053-8119
(2000).
Google Scholar
M. L. Pierro et al.,
“Phase-amplitude investigation of spontaneous low-frequency oscillations of cerebral hemodynamics with near-infrared spectroscopy: a sleep study in human subjects,”
NeuroImage, 63
(3), 1571
–1584 https://doi.org/10.1016/j.neuroimage.2012.07.015 NEIMEF 1053-8119
(2012).
Google Scholar
F. Martelli et al., Light Propagation through Biological Tissue and Other Diffusive Media: Theory, Solutions, and Software, SPIE Press, Bellingham, Washington
(2010). Google Scholar
L. Spinelli et al.,
“Accuracy of the nonlinear fitting procedure for time-resolved measurements on diffusive phantoms at NIR wavelengths,”
Proc. SPIE, 7174 717424 https://doi.org/10.1117/12.808997 PSISDG 0277-786X
(2009).
Google Scholar
L. Zucchelli et al.,
“Method for the discrimination of superficial and deep absorption variations by time domain fNIRS,”
Biomed. Opt. Express, 4
(12), 2893 https://doi.org/10.1364/BOE.4.002893 BOEICL 2156-7085
(2013).
Google Scholar
F. Scholkmann et al.,
“A review on continuous wave functional near-infrared spectroscopy and imaging instrumentation and methodology,”
NeuroImage, 85 6
–27 https://doi.org/10.1016/j.neuroimage.2013.05.004 NEIMEF 1053-8119
(2014).
Google Scholar
G. E. Strangman, Z. Li and Q. Zhang,
“Depth sensitivity and source-detector separations for near infrared spectroscopy based on the Colin27 brain template,”
PLoS One, 8
(8), e66319 https://doi.org/10.1371/journal.pone.0066319 POLNCL 1932-6203
(2013).
Google Scholar
A. V. Oppenheim, R. W. Schafer and J. R. Buck, Discrete-time Signal Processing, 2nd edPrentice-Hall, Inc., Upper Saddle River, New Jersey
(1999). Google Scholar
H. J. Nussbaumer, Fast Fourier Transform and Convolution Algorithms, Springer Berlin Heidelberg, Berlin, Heidelberg
(1982). Google Scholar
A. Pifferi et al.,
“Performance assessment of photon migration instruments: the MEDPHOT protocol,”
Appl. Opt., 44 2104
–2114 https://doi.org/10.1364/AO.44.002104
(2005).
Google Scholar
A. Torricelli et al.,
“Time-resolved reflectance at null source-detector separation: improving contrast and resolution in diffuse optical imaging,”
Phys. Rev. Lett., 95
(7), 078101 https://doi.org/10.1103/PhysRevLett.95.078101 PRLTAO 0031-9007
(2005).
Google Scholar
T. Sato et al.,
“Reduction of global interference of scalp-hemodynamics in functional near-infrared spectroscopy using short distance probes,”
NeuroImage, 141 120
–132 https://doi.org/10.1016/j.neuroimage.2016.06.054 NEIMEF 1053-8119
(2016).
Google Scholar
R. B. Saager and A. J. Berger,
“Direct characterization and removal of interfering absorption trends in two-layer turbid media,”
J. Opt. Soc. Am. A, 22
(9), 1874 https://doi.org/10.1364/JOSAA.22.001874 JOAOD6 0740-3232
(2005).
Google Scholar
BiographyLetizia Contini received her MS degree in biomedical engineering from Politecnico di Milano, Italy, in December 2023. Currently, she is pursuing a PhD in physics at the same institution, and her research delves into time-domain near-infrared spectroscopy (TD-NIRS), focusing on both device advancement and data analysis for applications in clinical and research settings. Caterina Amendola is a dedicated postdoctoral researcher specializing in the field of diffuse optics. She is employed in the Physics Department of Politecnico di Milano (Italy), where she earned her PhD in March 2022. Her research primarily focuses on time-domain near-infrared spectroscopy (TD-NIRS) and diffuse correlation spectroscopy (DCS) techniques. With expertise in developing hybrid TD-NIRS and DCS devices, she applies these advancements in clinical settings while conducting data analysis and physical modeling. Davide Contini received his MS degree in electronic engineering and his PhD in physics from Politecnico di Milano, Italy, in 2004 and 2007, respectively. He is an associate professor in the Department of Physics, Politecnico di Milano. He has authored more than 150 papers in international peer-reviewed journals and conference proceedings. His research activity is focused on time-resolved spectroscopy of highly diffusive media for applications in biology and medicine. Alessandro Torricelli is a full-time professor in the Department of Physics, Politecnico di Milano, Italy. He received his MS degree in electronic engineering from Politecnico di Milano in 1994 and his PhD in physics from Politecnico di Torino, Italy, in 1999. He is the author of more than 200 papers in international peer-reviewed journals. His current research interests include photon migration in diffusive media, functional near-infrared spectroscopy, and noninvasive diffuse spectroscopy with time-domain systems. Lorenzo Spinelli received his MS degrees and PhD in physics from the University of Milan, Italy, in 1994 and 1999, respectively. He devoted his research to the study of structures developing in the section of broad area radiation beams. In 2001, he became a researcher for the Italian Research National Council at the Institute of Photonics and Nanotechnologies. His current research interest is the study of photon migration in turbid media for optical biopsy and imaging. Rebecca Re is an associate professor at the Physics Department of Politecnico di Milano, Italy. Her research activity focuses on light propagation in turbid media with particular interest in biological tissues. In particular, she employs time-domain near-infrared spectroscopy and diffuse correlation spectroscopy to investigate brain hemodynamics and muscle oxidative metabolism and microcirculation in healthy subjects and patients. Her expertise ranges from device development to data analysis and physical models. |
Hemodynamics
Cerebral cortex
Neurophotonics
Brain
Photon counting
Signal detection
Optical properties