Time-Resolved Information-Theoretic and Spectral Analysis of fNIRS Signals from Multi-Channel Prototypal Device

TL;DR

This study uses information-theoretic and spectral methods to analyze fNIRS brain signals during a breath-holding task, revealing how respiratory activity affects scalp hemodynamics.

Contribution

The study introduces a time-resolved information-theoretic framework to analyze fNIRS signals and their modulation by respiratory activity.

Findings

Conditional entropy showed significant modulation during breathing and apnea phases.

Power spectral density revealed distinct vascular and respiratory frequency bands.

High-frequency oscillations were significantly modulated by respiratory activity.

Abstract

Functional near-infrared spectroscopy (fNIRS) is a non-invasive imaging technique that measures brain hemodynamic activity by detecting changes in oxyhemoglobin and deoxyhemoglobin concentrations using light in the near-infrared spectrum. This study aims to provide a comprehensive characterization of fNIRS signals acquired with a prototypal continuous-wave fNIRS device during a breath-holding task, to evaluate the impact of respiratory activity on scalp hemodynamics within the framework of Network Physiology. To this end, information-theoretic and spectral analysis methods were applied to characterize the dynamics of fNIRS signals. In the time domain, time-resolved information-theoretic measures, including entropy, conditional entropy and, information storage, were employed to assess the complexity and predictability of the fNIRS signals. These measures highlighted distinct informational dynamics across the breathing and apnea phases, with conditional entropy showing a significant modulation driven by respiratory activity. In the frequency domain, power spectral density was estimated using a parametric method, allowing the identification of distinct frequency bands related to vascular and respiratory components. The analysis revealed significant modulations in both the amplitude and frequency of oscillations during the task, particularly in the high-frequency band associated with respiratory activity. Our observations demonstrate that the proposed analysis provides novel insights into the characterization of fNIRS signals, enhancing the understanding of the impact of task-induced peripheral cardiovascular responses on NIRS hemodynamics.

Figures

Full text

1. Introduction

Functional near-infrared spectroscopy (fNIRS) is a non-invasive imaging technique which relies on the use of light at specific wavelengths for measuring the hemodynamic activity of brain tissue [1]. Specifically, the near-infrared region (650–950 nm) enables the estimation of the concentration of oxyhemoglobin $[eqn]$ and deoxyhemoglobin $[eqn]$ by exploiting their different absorption spectra at different wavelengths [2,3]. Among the different fNIRS techniques available on the market, the continuous-wave (CW) modality employs constant tissue illumination to quantify light attenuation through the head. Given its reliance on low-cost electronic components, as well as its high portability and ease of use, the CW modality has become the most commonly used technique. Indeed, only recently novel multichannel fNIRS systems were developed to answer the need of a portable system that can monitor hemodynamic activity in ecological and clinical settings [4,5].

In a recent study, an integrated system was developed that combines a novel continuous-wave fNIRS device with a modified commercial electroencephalographic (EEG) system [4]. This hybrid system has been validated for monitoring Event-Related Potentials across three distinct experimental paradigms designed to elicit the following: (i) motor cortex activity via a finger-tapping task [6]; (ii) visual cortex activity using a flickering black-and-white checkerboard at a frequency of 2 Hz [7]; and (iii) frontal cortex activity through a Stroop task [8].

However, it is well established that cognitive and emotional tasks can influence peripheral physiology by altering heart rate, respiration, blood pressure, and skin conductance [9,10]. These physiological changes affect both scalp and brain hemodynamics, thereby contributing to the fNIRS signal in ways that extend beyond neural activity [11]. In particular, systemic physiological variations—especially those driven by autonomic or cardiorespiratory activity—can impact both superficial and intracerebral hemodynamics. Previous research has shown that signals recorded from a given scalp location do not exclusively reflect neurovascular coupling but may also result from hemodynamic changes in extracerebral tissues (e.g., skin, muscles, skull) as well as in deeper brain regions influenced by systemic processes [12]. This issue is especially pronounced in CW fNIRS systems with short source–detector separations (typically < 4 cm), where the sensitivity to extracerebral tissues often exceeds that to cortical areas. Recent studies confirm that the extracerebral compartment—including scalp, skull, and cerebrospinal fluid—absorbs a substantial portion of the incident light, often leading to signals dominated by superficial layers [13,14,15]. This limited depth sensitivity remains a major limitation of CW fNIRS, particularly in adult populations, and must be considered when interpreting signal origin.

As a result, fNIRS measurements are affected by various background physiological components that are often considered “noise” in traditional analyses of brain connectivity, as they can obscure the detection of evoked neural responses. Yet, the influence of these systemic signals on both brain and scalp hemodynamics remains underexplored. Most analysis techniques aim to remove such components from the fNIRS signal [16,17]. Recent developments in the field of network physiology, however, propose a paradigm shift: rather than viewing these fluctuations as artifacts, they can be interpreted as meaningful indicators of physiological interactions [18,19]. According to this perspective, the human body can be seen as an integrated network, where each physiological system, despite having its own regulatory mechanisms, continuously interacts with others to coordinate functions and generate distinct physiological states in both health and disease [20]. In this framework, cardiorespiratory dynamics arise from the coordination of respiratory, cardiovagal, and sympathetic functions, which operate across various levels of the nervous system. These interactions play a critical role in maintaining physiological equilibrium and have only recently been shown to have a direct impact on both scalp and neural signals in specific clinical conditions [21,22]. To investigate these complex interactions, we designed a protocol based on the breath-holding task, in which the presence of respiratory signals is systematically modulated, allowing for a detailed assessment of the influence of respiratory activity on the local concentration of oxy- and deoxyhemoglobin measured over the scalp [23]. While previous studies have demonstrated strong responses in all brain chromophores due to respiratory activity using well-established methods [17], traditional approaches do not fully characterize the regularity, complexity, or transient dynamics of these signals.

Therefore, the aim of this work is to characterize the effects of respiratory activity on scalp hemodynamics as measured by a prototypal fNIRS system. We employ a multi-domain analysis that includes state-of-the-art signal processing techniques, which have not yet been extensively applied in the context of fNIRS signal analysis. These techniques are designed to assess physiological oscillations and their modulation during breath-holding, as well as the transient dynamics during the intermittent presence of respiratory activity. Specifically, to investigate the complexity and predictability of the fNIRS signals dynamics, we leverage information-theoretic measures computed in a time-resolved fashion such as conditional entropy (CE) [24] and information storage (IS) [25]. Moreover, given the presence of oscillatory activity that can be measured over the scalp, we perform a frequency-domain analysis to explore the distribution of power across frequency bands, providing an indirect measure of both the intensity and temporal evolution of neural activity [26].

2. Materials and Methods

2.1. Dataset Description and Pre-Processing

A total of 6 healthy subjects (4 males, 2 females; mean age: 27 ± 2.5 years) were enrolled in this study and provided informed consent prior to participation. The experimental procedure was approved by the ethical committee of the University of Palermo and was conducted in accordance with the ethical standards of the Helsinki Declaration. The subjects were seated in a chair in an isolated room and instructed to breathe normally, avoiding other actions or excessive movements. The experimental protocol consisted of a five-minute recording alternating two different phases: (i) apnea episodes, during which subjects were instructed to hold their breath, starting at predefined time points (8, 50, 110, 170, and 230 s), and (ii) a breathing phase, during which the subject resumed normal breathing before the onset of the subsequent apnea.

The signal acquisitions were performed using a multi-channel CW-fNIRS system. This system includes 16 bicolor light-emitting diodes (LEDs) operating at wavelengths of 735 nm and 850 nm, alongside 16 silicon photomultipliers (SiPMs). Each SiPM is positioned approximately 3 cm away from its corresponding LED source and acquires signals at a sampling frequency ( $[eqn]$ ) of 130 Hz. This is in line with previous studies with a similar setup [22,27], and with a work demonstrating that, under optimal conditions, optode spacings up to 5 cm are usable with NIRS equipment [28]. According to a “banana-shaped” model, which assumes that the penetration is similar for all wavelengths, the resulting wavepath had a maximum penetration depth of about 1.5 cm [29,30]. Thus, these fNIRS measurements are a superposition of signals from multiple tissue layers, with dominant contributions from extracerebral tissues [14,15]. Further details about the prototype electronic components and the realization procedure can be found in [4]. The selected configuration of the optodes is shown in Figure 1, where each optode is labeled with the letter ‘S’ (source) if it is an LED, or ‘D’ (detector) if it is a SiPM. Each optode is also assigned a number to identify its position on the scalp. In particular, the positions of the 6 LEDs are marked in red, while the 6 SiPMs are marked in blue. Together, they provide 16 different channels, which are listed in the table of the same figure.

The pre-processing pipeline employed in this study included the following steps: (i) the raw intensity signals were converted to optical density (OD) signals using the equation $[eqn]$ , where $[eqn]$ is the time-dependent recorded signal intensity and $[eqn]$ is its initial value [4]; (ii) the Temporal Derivative Distribution Repair (TDDR) method was applied to the OD signals as implemented in the fNIRS Brain AnalyzIR Python toolbox [31]. TDDR is a motion correction algorithm designed to eliminate the two most common motion artifacts in fNIRS data: spikes and baseline shifts; (iii) the corrected OD signals were converted to hemoglobin concentration signals ( $[eqn]$ and $[eqn]$ ) using the modified Beer–Lambert law [32], following the procedure described in [4]. The resulting signals were then low-pass-filtered with a cutoff frequency of 0.5 Hz to prevent aliasing and subsequently downsampled to 1 Hz; (iv) finally, a high-pass autoregressive (AR) filter with infinite impulse response (IIR) and zero phase [33], with a cutoff frequency of 0.018 Hz, was applied to emphasize variations in the respiratory frequency band induced by the task. The resulting time series, each consisting of $[eqn]$ samples, were normalized to zero mean and unit variance. Representative trends for one channel and one subject are shown in Figure 2a.1 for $[eqn]$ and Figure 2a.2 for $[eqn]$ , highlighting the variations observed during the breath phases (red windows) and the apnea phases (blue windows).

2.2. Time-Resolved Information Measures

Information-theoretic measures provide powerful approaches to assess the regularity and predictability of the dynamics of a given system. While traditional analyses provide overall measures of complexity and predictability, recent advancements extend these measures to enable time-resolved characterization [24,34,35]. Here, we apply these methods to study the transient evolution of scalp hemodynamics during breath-holding tasks—revealing dynamic responses in extracerebral oxygenation that conventional fNIRS analyses cannot resolve.

The analysis of any given dynamical system can be performed by mapping the system activity with a set of random variables and then studying the statistical dependencies among the observed realizations of the variables collected in the form of time series. In the general field of information theory, starting from the basic concept of entropy introduced by Shannon [36], it is possible to dissect the information processed in a dynamical system into meaningful elements of computation to quantify the complexity and predictability of the aforementioned system.

The $[eqn]$ or $[eqn]$ signal, acquired for each subject and channel, was modeled as a realization of a zero-mean stochastic process X, with $[eqn]$ representing the random variable sampling the process at time n (temporal counter). Then, the Shannon entropy (H) can be used to quantify the amount of information needed to describe $[eqn]$ as follows: $[eqn]$ , where $[eqn]$ refers to a realization of $[eqn]$ and $[eqn]$ is the probability density function of $[eqn]$ measured for the outcome $[eqn]$ , and $[eqn]$ is the expectation operator computing the statistical average over all possible values of $[eqn]$ . Assuming X as a Markov process, its past history can be truncated up to a lag p to obtain the p-dimensional vector $[eqn]$ . To account for the dynamical evolution of the system X, it is possible to quantify the new information carried by the present state $[eqn]$ of the process which cannot be inferred from its past $[eqn]$ , i.e., the Conditional Entropy (CE) [37]: $[eqn]$ where $[eqn]$ refers to a realization of $[eqn]$ and $[eqn]$ is the conditional probability. Specifically, the higher the CE, the more complex the dynamics generated by the system. Entropy and CE serve as measures to introduce the Information Storage (IS), which quantifies the amount of information contained in the present state that can be predicted by the knowledge of its past state [38] and can be defined as

$$[eqn]$$

The IS measures the predictability of the process at time n by quantifying the average level of uncertainty about the current state of the process $[eqn]$ that can be resolved by knowledge of its past states $[eqn]$ .

To describe the time-resolved behavior of the fNIRS signal ( $[eqn]$ , $[eqn]$ ) during breathing and apnea phases, a time-varying analysis of the measures of H, CE, and IS was performed, through the adoption of a time-varying version of an autoregressive (TV-AR) model [24,34]:

$$[eqn]$$

where $[eqn]$ represents the prediction error, while $[eqn]$ denotes the AR coefficient describing the interaction from $[eqn]$ to $[eqn]$ at lag k, relevant to the time instant n. Then, under Gaussian assumption of $[eqn]$ , the time-resolved H can be expressed as follows [24]: $[eqn]$ , where $[eqn]$ represents the variance of $[eqn]$ . Moreover, if $[eqn]$ and $[eqn]$ are jointly Gaussian, the time-resolved CE of $[eqn]$ given $[eqn]$ can be expressed as follows [39]: $[eqn]$ , where $[eqn]$ is the variance of the prediction error $[eqn]$ at the time step n. Thus, the equation given in (1) can be rewritten as

$$[eqn]$$

The identification procedure of the TV-AR model (2) can be performed through the recursive least-squares (RLS) as described in [40]. Briefly, the RLS consists of the following steps: (i) choose a value for the adaption factor $[eqn]$ and an order p of the AR model; (ii) define proper initial conditions for the vector of coefficients at time p, $[eqn]$ and for the correlation matrix of the past state of X stored in $[eqn]$ , $[eqn]$ ; (iii) considering N-p consecutive time steps, for $[eqn]$ to N, repeat the following steps:

$$[eqn]$$

where $[eqn]$ is the so-called gain vector and $[eqn]$ is intended as the a priori estimation error before updating the AR coefficients vector. The parameter (1 − c) controls the memory of the algorithm allowing it to follow possible statistical variations in the property of X in non-stationary conditions. To complete the identification procedure, a recursive estimation of the time-varying innovation variance can be obtained as follows: $[eqn]$ [24,34]. The required recursive estimation of the process variance $[eqn]$ can be directly derived from the structure of the linear TV-AR representation of the process X as described in [24].

2.3. Spectral Analysis of fNIRS Signals

Biomedical signals are often rich with oscillatory content, and therefore naturally lend themselves to spectral representation. Classical approaches integrating the power spectral density profile within the spectral bands of interest attempt to obtain band-specific time-domain powers, while the method of spectral decomposition used here allows for focusing only on the spectral components with frequencies within those bands, thus avoiding spurious contributions due to broadband oscillations [41].

To characterize the oscillatory content of each fNIRS signal ( $[eqn]$ , $[eqn]$ ), a linear model governed by the following difference equation was employed:

$$[eqn]$$

where $[eqn]$ is the prediction error with a variance $[eqn]$ and $[eqn]$ is the AR coefficient describing the interaction from $[eqn]$ to $[eqn]$ at lag k. The linear model (5) can be represented in the Z-domain through its Z-transform yielding $[eqn]$ , where $[eqn]$ is the transfer function relating the Fourier Transform (FT) of U to the FT of the process X. Applying the residue theorem, $[eqn]$ can be expressed as follows [41]:

$$[eqn]$$

where $[eqn]$ , $[eqn]$ , are the p poles of the AR process, while the terms $[eqn]$ are pole-specific factors associated each with a given pole $[eqn]$ , with * indicating the Hermitian transpose. The power spectral density (PSD) of the process X can be expressed in the Z-domain as $[eqn]$ . Then, by using the Heaviside decomposition, it can be expanded into simple fractions corresponding to the poles of the system. These include the poles inside the unit circle, $[eqn]$ , and their reciprocals outside the unit circle, $[eqn]$ , for $[eqn]$ , weighted by the residuals of $[eqn]$ , specifically $[eqn]$ and $[eqn]$ . This results in the following expression [41]:

$$[eqn]$$

where the residuals are given by $[eqn]$ , $[eqn]$ . Evaluating $[eqn]$ on the unit circle of the complex plane, specifically $[eqn]$ , where $[eqn]$ and $[eqn]$ is the sampling frequency (1 Hz), allows for the derivation of the spectral profile, $[eqn]$ , along with its $[eqn]$ component, $[eqn]$ . Each spectral component is characterized by a unique profile, which is determined by the central frequency of the oscillation, derived from the argument of the pole ( $[eqn]$ ), and by the power related to the residual of the pole. For real poles, the power is $[eqn]$ , and for complex conjugate poles, it is $[eqn]$ . Notably, the total variance $[eqn]$ of the process is equal to the sum of the variances of all poles, $[eqn]$ , for $[eqn]$ .

The computation of the PSD of X relies on the identification procedure of the AR model (5) which herein was performed through the ordinary least-squares (OLS) method [42]. By recalling the p-dimensional vector $[eqn]$ containing the past states of X and considering N consecutive time steps, a compact representation of the VAR model (5) can be defined as $[eqn]$ , where $[eqn]$ is the $[eqn]$ vector of unknown coefficients, $[eqn]$ and $[eqn]$ are $[eqn]$ vectors collecting the present states and the residuals, and $[eqn]$ is a $[eqn]$ matrix collecting the regressor terms. Then by using the OLS formula an estimate of the AR coefficient can be obtained as follows: $[eqn]$ . The innovation process can be estimated as the residual time-series $[eqn]$ , whose $[eqn]$ is an estimate of the innovation variance.

2.4. Data Analysis and Statistical Validation

All the analyses were performed by considering separately the concentrations of $[eqn]$ and $[eqn]$ for each channel and for each subject. Regarding their spectral analysis, two distinct identification procedures for the model (5) were performed, separately analyzing the breathing phase and the apnea task. Specifically, given $[eqn]$ as a vector representing the time instants of the apnea task onsets, with the last element indicating the total length of the time series, we defined $[eqn]$ and $[eqn]$ , $[eqn]$ , as the $[eqn]$ -th time window during which the subjects breathed normally (red windows in Figure 2) or performed the apnea task, respectively (blue windows in Figure 2). Let $[eqn]$ denote a generic instant when the subject begins to breathe normally after an apnea task (which may differ for each apnea event and subject). Using this notation, we can define the following time windows: $[eqn]$ , and $[eqn]$ . Then when the breath task was considered, we used $[eqn]$ as the present state, while the matrix of regressors $[eqn]$ was defined by simply switching the content of each window $[eqn]$ of one temporal unit in the past, up to p lags. The same rationale was used for the apnea task where $[eqn]$ . This procedure was implemented to prevent potential discontinuities that could arise from simply concatenating different time windows and to increase the number of observations, thereby reducing the estimation bias [43]. Spectral profiles were computed using (7), and the low-frequency (LF) and high-frequency (HF) components were identified based on poles with central frequencies in the ranges [0.04–0.15] Hz and [0.15–0.4] Hz [44]. The PSD profile relevant to each pole, $[eqn]$ , was then integrated within the whole frequency spectrum to obtain the power content associated with the LF and HF components, respectively. A representative example of the spectral decomposition performed on one representative channel for $[eqn]$ and $[eqn]$ is reported in Figure 2 (panel e: apnea; panel f: breath). The code used for performing the Power Spectral Density Analysis can be found here: https://github.com/YuriAntonacci/LSP_toolbox/tree/main, accessed on 12 December 2024. To test whether the spectral content of fNIRS signals was modulated during the breath-holding task, a statistical comparison between the distributions of central frequency and spectral power obtained across subjects during the breath and apnea phases was performed using the paired non-parametric Wilcoxon signed-rank test $[eqn]$ .

As for the time-resolved analysis of the different information-theoretic measures, the identification procedure of the TV-AR model (2) was performed setting $[eqn]$ and $[eqn]$ . The model order p was selected using the Akaike Information Criterion (AIC) [42] as a guide. Analysis across the entire scalp and for all subjects returned values ranging from four to six. To ensure consistency and avoid duplicate or negative peaks in the subsequent spectral analysis, the model order was fixed at four [44]. As regard for the forgetting factor $[eqn]$ , the selection was performed according with previous studies, which highlighted the interval [0.97, 0.98] as an ideal compromise for the bias–variance trade-off [24,34].

The time-resolved measures of entropy ( $[eqn]$ ), conditional entropy ( $[eqn]$ ), and information storage ( $[eqn]$ ) (panels b–d of Figure 2) were then computed and averaged within each i-th time window corresponding to the analyzed experimental condition (breath and apnea). The obtained values were further averaged across the five available windows for each experimental condition to obtain a representative value for each of the six subjects. Statistically significant differences between the apnea and breath conditions were evaluated for each information-theoretic measure using the Wilcoxon test for paired data $[eqn]$ . The code necessary to compute time-resolved information measures can be found here: https://github.com/YuriAntonacci/Time-VaryingIS, accessed on 12 December 2024.

Therefore, a measure of the effect size was also determined to assess the magnitude of the differences observed among experimental conditions. Specifically, denoting with $[eqn]$ , $[eqn]$ , $[eqn]$ , and $[eqn]$ , the mean and the variance of two distributions $[eqn]$ and $[eqn]$ obtained measuring the mean and the SD of the distribution across subjects of each information theoretic measure, we computed the Cohen’s d measure defined for equally sized groups as follows [45]:

$$[eqn]$$

Typically, a small effect size arises for $[eqn]$ , a medium effect size between 0.2 and 0.8, and large was considered when $[eqn]$  [45]. Finally, we assessed the statistical power of each test to evaluate how the sample size (N = 6) would influence the results. Statistical power ( $[eqn]$ ) represents the probability of correctly rejecting a false null hypothesis (i.e., the likelihood of detecting a true physiological effect if one exists).

3. Results

3.1. Frequency Specific Analysis of fNIRS Signals

Panels e,f of Figure 2 present the spectral decomposition of the PSD for $[eqn]$ (panel 1) and $[eqn]$ (panel 2) signals from one representative subject during the apnea (panel e) and breath (panel f) phases. Regardless of the analyzed time series ( $[eqn]$ or $[eqn]$ ), the PSD trends reveal two distinct peaks centered at approximately 0.05 Hz and 0.2 Hz with the latter showing a modulation when the respiratory activity is considered.

Figure 3 and Figure 4 report the boxplot distributions illustrating the frequency locations of the LF and HF spectral peaks for the apnea (blue) and breath (red) phases, computed separately for $[eqn]$ (panel a) and $[eqn]$ (panel b). Figure 3 displays the distribution of spectral peak frequencies within the [0.04–0.15] Hz frequency band, which reveals negligible modulation between apnea and breathing phases, indicating that respiratory activity has limited influence on low-frequency oscillations.

Figure 4 displays the distribution of spectral peak frequencies within the [0.15–0.4] Hz frequency band, revealing significant task-dependent shifts. Overall, regardless of the time series analyzed (i.e., $[eqn]$ or $[eqn]$ ), there is a consistent trend of increased oscillation frequency observed during the breathing task across all fNIRS channels. This increase becomes statistically significant for $[eqn]$ , occurring in four specific channels (Figure 4b). Remarkably, the corresponding Cohen’s d values consistently exceed 1.2 ( $[eqn]$ 0.6), indicating a strong difference between the two experimental conditions.

Figure 5 and Figure 6 display the boxplot distributions obtained for the spectral power of the fNIRS signal relevant to the LF and HF frequency bands, respectively. The results highlight an increase in both LF and HF power associated with breathing activity. Though this increase is statistically significant only for a few channels, when $[eqn]$ is considered, the corresponding values of Cohen’s d are always greater than 1 ( $[eqn]$ ), highlighting the presence of a strong modulation of the spectral content of fNIRS signals as an effect of respiratory activity.

3.2. Time-Resolved Information-Theoretic Measures

Figure 7, Figure 8 and Figure 9 display the boxplot distributions of the time-resolved information measures of entropy, conditional entropy, and information storage, averaged across five different time windows relevant to the apnea (blue color) and breath (red color) phases, computed separately for $[eqn]$ (panel a) and $[eqn]$ (panel b). The results in Figure 7 suggest a modulation of uncertainty in the oxyhemoglobin and deoxyhemoglobin time series during breathing compared to the apnea condition. However, this difference is not statistically significant, indicating that the amount of information contained in the concentrations of oxyhemoglobin and deoxyhemoglobin remained largely unchanged despite the breath-holding task. Conversely, Figure 8 and Figure 9 reveal statistically significant modulations of the time-resolved conditional entropy and the information storage measures, respectively. In detail, the analysis of time-resolved conditional entropy reveals a statistically significant decrease in the complexity of scalp hemodynamics during breathing periods, regardless of the hemoglobin time series considered, and across nearly all analyzed channels. In contrast, time-resolved information storage shows a statistically significant increase in the predictability of scalp hemodynamics during the same periods. A variability can be observed in these measures across different channels even though Figure 8 and Figure 9 consistently show Cohen’s d values > 0.6 with statistical power $[eqn]$ .

4. Discussion

The main results of this work can be summarized as follows: (a) the presence of two different peaks in the PSD of fNIRS signals centered approximately at 0.05 Hz and 0.2 Hz and mainly related to vascular and respiratory components [46,47,48,49]; (b) the modulation in both the amplitude and frequency of oscillations in the time series of oxyhemoglobin and deoxyhemoglobin concentrations during the breath-holding task, which may be sustained by the strong coupling between blood pressure and the hemodynamic responses measured in extracerebral tissues [22,50]; (c) the presence of a transition in the information processing characteristics of fNIRS dynamics associated with the shift between breathing and apnea phases, which is marked by an increase in signal predictability during respiratory activity—an effect that cannot be captured using simple measures such as Shannon entropy.

The results reported in Figure 3 and Figure 5 highlighted the presence of a vascular component, known as Mayer waves, consisting of rhythmic fluctuations in arterial pressure caused by changes in vasomotor tone, typically occurring within frequencies ≤0.1 Hz, which are reflected in the hemodynamic oscillations of superficial scalp tissues, thereby modulating the optical fNIRS signals [47]. The reduced amplitude of low-frequency components (<0.1 Hz) observed during apnea, even if limited to a few channels, may be attributed to the induction of hypercapnia, i.e., an increase in CO_2_ levels during breath-holding tasks in humans [51]. Indeed, the temporary suspension of breathing activity leads to a rise in carbon dioxide levels, triggering a vasodilation response. This vasodilation increases cerebral blood flow, which directly modulates the amplitude of low-frequency oscillations associated with Mayer waves [52].

The results displayed in Figure 4 and Figure 6 revealed the presence of a spectral component in the 0.15–0.4 Hz range, which can be attributed to respiratory activity. This component reflects respiratory-driven variations in arterial pressure as well as in cerebral and extracerebral blood flow. The observed shift in frequency location, occurring exclusively within the HF band, has been previously reported and is physiologically linked to increased respiratory rates during the breathing phase compared to the apnea phase [53]. Interestingly, the peak relevant to the HF range was also observed during apnea and reflects fundamental respiratory physiology: while voluntary breath-holding suppresses overt chest movements, it cannot interrupt the brainstem’s central respiratory rhythm [54]. We speculate that these HF oscillations can be related to respiratory sinus arrhythmia—a well-documented phenomenon where autonomic nervous system activity causes a person’s heart rate to oscillate in phase with the residual respiratory drive, even during apparent apnea [55]. However, the expected increase in the spectral content within HF components during the breathing phase, as reported in previous studies, was not clearly observed in our analysis—possibly due to the limited duration of the experimental period, which should ideally last at least 28–30 s to reliably capture such changes [56].

Overall, our results revealed dynamic changes in both vascular and respiratory components for oxygenated and deoxygenated hemoglobin, consistent with established mechanisms of respiratory–cardiovascular coupling [54]. Previous studies show that spontaneous hemodynamic oscillations, although not directly elicited by external stimuli or tasks, can still exhibit variations in amplitude and frequency. These fluctuations significantly affect frequency-domain analyses, as they shape the spectral profile of the recorded signals and underscore the influence of underlying physiological rhythms on hemodynamic activity [11,57,58,59].

The results shown in Figure 7, Figure 8 and Figure 9 indicated less-complex fNIRS signal dynamics associated with the presence of respiratory activity. Although not directly comparable with previous studies characterizing brain dynamics and physiological systems using information-theoretic approaches [24,34,38,60], these works consistently show that the emergence of a predominant oscillation—often driven by synchronization phenomena—is associated with increased predictability of brain dynamics across different contexts [61]. These findings are consistent with the PSD analysis of fNIRS signals dynamics herein performed, where the breathing task was associated with an increase in the PSD within both the LF and HF bands. Moreover, since a transition occurs after the onset of apnea, we can link the observed decrease in predictability to previous studies, which have highlighted that such transitions are associated with a reduction of the time-resolved information storage [34]. Lastly, we observed a variability in the information-theoretic measures across fNIRS channels, which may originate from the following: (i) spatial heterogeneity in fNIRS sensitivity (superficial vs. deep tissue contributions) [14]; (ii) regional variations in neurovascular coupling [62]; and (iii) task-specific activation patterns [63]. Despite this channel-wise variability, the measures of time-resolved information storage and conditional entropy consistently showed Cohen’s d values > 0.6 with statistical power $[eqn]$ , indicating robust effects of breath-holding on both time-resolved conditional entropy and information storage.

Final Remarks and Limitations

Our results demonstrated that integrating time-resolved information-theoretic measures (entropy, conditional entropy, information storage) with parametric spectral analysis can yield novel insights into fNIRS signal dynamics. Power spectral density analysis revealed distinct spectral signatures: Mayer waves (0.04–0.15 Hz) exhibited minimal modulation during breath-holding, suggesting autonomic stability, while high-frequency oscillations (0.15–0.4 Hz) persisted during apnea, reflecting the ongoing influence of the central respiratory rhythm and respiratory sinus arrhythmia. Time-resolved measures further showed that the presence of respiratory activity reduces hemodynamic complexity and increases signal predictability compared to apnea, underscoring the role of respiration in shaping scalp-measured hemodynamic activity.

The main limitation of our study is that the presented results cannot be directly attributed to cortical activity, as CW single-distance NIRS is known to be more sensitive to superficial tissues than to brain tissue. The measured signals likely reflect a superposition of responses from multiple depths, with a substantial contribution from scalp hemodynamics [14]. Indeed, both incident and backscattered light traverses highly vascularized extracerebral layers such as the skin, subcutaneous fat, and scalp muscles, making the signals particularly sensitive to hemodynamic fluctuations in these regions [30].

Finally, the generally low statistical power and small effect sizes observed in our analyses suggest that some nonsignificant results may be attributed to insufficient statistical sensitivity, primarily due to the limited sample size. This represents a key limitation of the study and highlights the need for larger datasets to enable more robust statistical assessments of the observed effects.

5. Conclusions

This study highlights the significance of fNIRS monitoring and advanced signal processing techniques in exploring scalp hemodynamics during physiological tasks, specifically the breath-holding task. As a non-invasive and portable modality, fNIRS offers a unique capability to monitor real-time hemodynamic responses, complementing other neuroimaging methods [4]. Our findings support the value of a multi-domain analytical approach for reinterpreting physiological “noise” as evidence of functional interactions within the framework of network physiology, and demonstrate the feasibility of applying advanced signal processing techniques even with prototype CW fNIRS systems. Indeed, combining fNIRS with information-theoretic measures allows for the identification of subtle dynamics in both cerebral and extracerebral functions, thereby deepening the physiological insights that can be gained from fNIRS data.

Our results reveal that novel aspects of respiratory–scalp coupling could enable practical applications in the following: (i) clinical monitoring of patients with impaired respiratory–cerebral coupling, such as those with sleep apnea or cerebrovascular disorders; (ii) passive brain–computer interfaces, where entropy-based measures could provide real-time assessment of cognitive state [64,65]. Future work should focus on the validation of these measures in larger cohorts (N > 30 subjects), comparing performance across different fNIRS systems (including time-domain devices), and developing hybrid fNIRS-EEG protocols to better separate neural and vascular contributions.