Evaluation of fNIRS signal components elicited by cognitive and hypercapnic stimuli

Functional near infrared spectroscopy (fNIRS) measurements are confounded by signal components originating from multiple physiological causes, whose activities may vary temporally and spatially (across tissue layers, and regions of the cortex). Furthermore, the stimuli can induce evoked effects, which may lead to over or underestimation of the actual effect of interest. Here, we conducted a temporal, spectral, and spatial analysis of fNIRS signals collected during cognitive and hypercapnic stimuli to characterize effects of functional versus systemic responses. We utilized wavelet analysis to discriminate physiological causes and employed long and short source-detector separation (SDS) channels to differentiate tissue layers. Multi-channel measures were analyzed further to distinguish hemispheric differences. The results highlight cardiac, respiratory, myogenic, and very low frequency (VLF) activities within fNIRS signals. Regardless of stimuli, activity within the VLF band had the largest contribution to the overall signal. The systemic activities dominated the measurements from the short SDS channels during cognitive stimulus, but not hypercapnic stimulus. Importantly, results indicate that characteristics of fNIRS signals vary with type of the stimuli administered as cognitive stimulus elicited variable responses between hemispheres in VLF band and task-evoked temporal effect in VLF, myogenic and respiratory bands, while hypercapnic stimulus induced a global response across both hemispheres.

Over the past two decades, functional near infrared spectroscopy (fNIRS) has emerged as an effective and viable modality to study brain activity during various cognitive and motor tasks or at resting states in healthy and diseased populations with acquired pathologies (i.e. Traumatic Brain Injury, Stroke) or induced stimuli (i.e. Carbon Dioxide -CO 2 , sevoflurane) [1][2][3][4][5][6] . Increases in metabolic demand within the brain is coupled with an increase in oxygen consumption and cerebral blood flow (CBF), which leads to highly correlated changes in oxygenated hemoglobin (HbO) and deoxygenated hemoglobin (HbR) that is readily detectable by fNIRS 7,8 . Although many studies have validated the sensitivity of fNIRS in relation to functional magnetic resonance imaging (fMRI) results, some studies have noted false positives and false negatives in fNIRS activation maps 9,10 . These results have been attributed to variation in head anatomy, physiology, and functioning and the near infrared (NIR) light propagation within this complex structure.
In cognitive studies, task-related vasodilation stemming from neuronal origins within the cerebral cortex as a result of neurovascular coupling (NVC) is the signal of interest 11 . Alternatively, in cerebral autoregulation (CA) or cerebrovascular reactivity (CVR) studies, changes originating from systemic factors, such as blood pressure or partial arterial CO 2 (PaCO 2 ) are the signals of interest 12,13 . However, regardless of stimuli, fNIRS measurements are composed of signal components originating from different: (i) physiological causes (i.e. neuronal vs systemic), (ii) spatial locations (axially for tissue layers i.e. cerebral vs extracerebral or laterally i.e. left vs right hemispheres), and (iii) experimental or task related elements (i.e. evoked vs resting state or with vs. without stressors) [14][15][16][17][18] . Therefore, extraction of the signal of interest requires application of appropriate signal processing techniques, which first entails thorough characterization and analysis of signal components arising because of the task being studied.
Systemic factors associated with cardiac, respiration, and myogenic (oscillations of the sympathetic vasomotor tone and blood pressure regulation mechanism by autonomic nervous system-ANS and baroreceptors) contribute at least 35% to fNIRS signals 19 . Cardiac, respiratory, and myogenic activities have been shown to appear within unique frequency bands of 0.6 to 1.6 Hz, 0.2 to 0.4 Hz, and 0.06 to 0.1 Hz, respectively 20,21 . Systemic activity associated with PaCO 2 shares the same band of 0.01 to 0.06 Hz as functional/neuronal activity, and has been projected to contribute at least 1.5% to fNIRS signals [22][23][24] . The simplest and the most utilized class of methods to remove frequencies associated with cardiac, respiratory and myogenic activities, consist of applying a band pass (cut-off frequencies at 0.01 and 0.09 Hz) or low pass (cut-off frequency at 0.09 Hz) Finite Impulse Response, or Butterworth filters 25 . Other filters such as Gaussian, Kalman, Wiener and hemodynamic response filters have also been commonly used to remove aforementioned systemic noises in real time or off-line signal processing applications 26 . However, such filters could suppress some of the task-related neuronal effects due to non-ideal filter characteristics and assumptions used in filter designs. Methods utilizing additional systemic physiological sensors such as electrocardiography, laser doppler flowmetry and others can be used for similar purposes [27][28][29][30] . However, addition of such sensors requires time synchronization between different sensors and may be cumbersome on the participant. Furthermore, these methods are unable to remove PaCO 2 , or extracerebral contributions 31,32 .
Due to the physics of photon migration within the head, sensitivity of NIR light is affected by various extracerebral (skin, scalp, skull, and cerebrospinal fluid) and cerebral layers (gray and white matter) of different optical characteristics and thickness. A comprehensive examination on the percentage of light absorbed by different brain layers in healthy and clinical conditions revealed that at best, only ~ 18% of light was absorbed by cerebral layers 33 . Without additional hardware requirements, blind source separation techniques such as Independent Component Analysis (ICA) or Principle Component Analysis (PCA) can be used to suppress extracerebral signals using multi-channel fNIRS measurements [34][35][36] . These techniques require assumptions regarding uncorrelatedness (as in PCA), orthogonality and statistical independence (as in ICA) between components to be guaranteed. However, such assumptions are typically violated when skin perfusion changes are also modulated by the task 28 . Alternatively, methods utilizing the addition of short source-detector separation (SDS) channels (0.5-1.0 cm separation) to the fNIRS sensor alongside the long SDS channels (2-4 cm separation) allow for direct separation of signals that have penetrated the cerebral layer and those that have not. The most direct techniques that utilize short SDS channels are the least squares adaptive filters and regression modelling [37][38][39][40][41] . Although these techniques improve contrast-to-noise ratio by 70%, they assume that extracerebral contributions are static and spatially homogeneous, which is not always accurate 28,35,42 . Independent of task, characteristics of extracerebral signals may vary spatially across different regions of the head due to the differences in terms of geometry and composition of tissue layers 34,43 . Kalman filter and state-space model-based methods help account for the nonstationary nature of hemodynamics 44,45 . Addition of multiple short SDS channels placed at varying length from a single detector allow for accountancy of spatial heterogeneity 42,46 .
The magnitude and frequencies of all the aforementioned fNIRS signal components may be influenced by the specific experimental tasks being used to investigate cognition, psychosocial stress, emotional state, or induced physiological stressors 28,29,[47][48][49] . Furthermore, it is proposed that neuronal and systemic factors arising from cerebral and extracerebral layers do not take place as isolated processes but come from an interactive network of interlinked processes 50,51 . For example, studies have demonstrated significant changes in myogenic activity during motor and cognitive tasks in cerebral and extracerebral layers 52,53 . Others have shown that PaCO 2 induced CBF changes can lead to an underestimation of NVC resulting CBF changes during out loud and inner speech-oriented cognitive tasks 54 . In a cognitively challenging n-back study, not removing extracerebral contribution from fNIRS measurements led to insignificant changes in functional activation with task difficulty 55 . However, removal of extracerebral contribution led to a decrease in HbO measures due to increases in skin conductance and heart rate, and an increase in HbO in the cerebral compartment due to NVC. While in some experimental protocols it is imperative to remove these systemic or extracerebral effects to obtain task-evoked cerebral neuronal signals, in other cases they can be treated as additional signals of interest as they may carry complementary information regarding the state of the subject 56,57 . For example, one study added time and frequency domain features of cardiac oscillations extracted from fNIRS signals to a machine learning architecture and reported increases in classification accuracy of stress by at least 10% when compared against standard models containing only the neuronal activity related information 58 . Another study reported that including end-tidal CO 2 (EtCO 2 -surrogate for PaCO 2 measures) in the interpretation of fNIRS signals improved ability to monitor pain 59 . These findings indicate that to comprehensively evaluate cerebral hemodynamic changes due to a task/stimulus, it is important to not just remove the non-neuronal effects, but to analyze contribution of each of the signal components to the fNIRS measurements, separately and in relation to each other.
Time frequency analysis serves as a good analytical tool to investigate fNIRS signal components within the same plane. The wavelet transform (WT) provides windows of adjustable length thereby ensuring good resolution for both high and low frequency content 60 . WT has been previously used in fNIRS signal processing to remove motion artifacts, and global drifts 61,62 . In fact, comparisons of the wavelet filter against the band pass filter, indicated the wavelet filter improved temporal consistency and spatial specificity over conventional methods 63 . WT analysis of fNIRS signals collected from prefrontal cortex (PFC) of healthy individuals have reported presence of local maxima in cardiac, respiratory, myogenic, and very low frequency (VLF -PaCO 2 or neuronal related activity) bands 29,[64][65][66] . Alterations in amplitude and frequency of these factors as indicated by WT analysis have been reported in studies investigating healthy aging, stroke, and cerebral infraction 23,[67][68][69] . However, existing studies have not used WT analysis on measurements taken from cerebral and extracerebral layers or applied the technique to compare evoked versus non-evoked activity or evaluated differences across regions of the PFC. www.nature.com/scientificreports/ The main aim in this work was to quantify the contribution of signal components from different origins including physiological causes, tissue layers and hemispheres on fNIRS recordings collected under different experimental paradigms. To achieve this goal, we used WT, long and short SDS channels, and multi-channel measurements. The data used here was acquired from two different studies. One study gathered data from a complex cognitive task, while the other one collected data while inducing hypercapnia 70,71 . The cognitive study employed realistic search and surveillance task scenarios on novice unmanned aircraft system (UAS) sensor operators. This experiment was chosen because preliminary results have shown strong NVC activity in left and right PFC regions associated with working memory and sustained attention 72 . In the second study, hypercapnia was the chosen stimulus, since it is known to result in a homogenous effect across all regions of PFC, does not require subject engagement and does not induce evoked NVC 71 . More importantly, the hypercapnic task allows for quantification of PaCO 2 's effect.

Methods
Study I: cognitive stimuli. Thirteen participants between the ages of 19 and 40 (3 female) were recruited to take part in this study. All participants provided written informed consent approved by the Institutional Review Board (IRB) of Drexel University. The research was performed in accordance with relevant regulations. Inclusion criteria required no prior UAS piloting simulator experience, normal or corrected to normal vision, and right handedness verified via the Edinburgh Handedness assessment 73 . A high-fidelity ground station simulator, the Simlat's C-STAR (Simlat Inc., Miamisburg, Ohio), was used to implement realistic scenarios of a sensor operator tasks. Five unique scenarios were created in the simulator, with the primary difference being the difficulty level (3 easy and 2 hard scenarios). Within each difficulty level, the scenarios were randomized. The exact administration of the task is provided in Fig. 1A. During each scenario (referred to as task condition from now on), participants were instructed to complete two synchronous tasks: (1) operate the UAS sensor camera to scan along the designated flight path, and (2) identify and track red civilian bus. Further details related to the simulation and experimental design can be found in 70 . Study II: hypercapnic stimuli. Five healthy controls between the ages of 24 and 41 (2 female) voluntarily consented to participate in two University of Pennsylvania IRB approved studies. The research was performed in accordance with relevant regulations. All participants had no history of traumatic brain injury (TBI), pregnancy, and pre-existing disabling neurological or psychiatric disorders. Hypercapnic and normocapnic conditions were administered to each subject via the fixed inspired method using a Douglas Bag 74 . The fixed inspired method delivers room air to the subject at the baseline/normocapnic condition and CO 2 enriched air during hypercapnic condition, which is kept at 5% CO 2 . The stimulus was provided in a block design manner, where normocapnia and hypercapnia were induced in an alternating manner. Each block was maintained for 60 s, and the experiment was conducted for a total of 7 min (see Fig. 1B). Detailed description of instrumentation set up and how the protocol was administered can be found in 71 . fNIRS instrumentation. During both studies, hemodynamic changes within the PFC of the participants were monitored using the fNIR Imager 1200 (fNIR Devices LLC, Potomac, MD). This continuous wave fNIRS system is composed of three parts: a flexible sensor pad that is placed over the forehead to monitor the underlying dorsal and inferior frontal cortical areas; a control box for hardware management; and a computer that runs the data acquisition and visualization software -COBI Studio.
The sensor consists of four surface mount light emitting diodes (LED), and twelve silicone photodiodes with integrated trans-impedance preamp housed on a flexible printed circuit board covered with a soft foam for comfort, insulation, durability, flexibility, and safety. The LEDs illuminate light at peak wavelengths of 730 nm and (a) In cognitive task, each participant started with a brief instructional session regarding UAS controls and task objectives, followed by five 12 min task conditions of varying difficulties. The first three task conditions were of easier difficulty, while the last two were of harder difficulty. A 5 and 15-min breaks were given between instructional and first easy condition, and last easy and first hard condition. (b) In hypercapnic task, each participant began with a 30 s baseline period (Room Air -0.04% CO 2 ), followed by three hypercapnic (5% CO 2 ) and normocapnic (Room Air) conditions of 60 s each. www.nature.com/scientificreports/ 850 nm, specifically selected to allow spectroscopic conversion without crosstalk 75 . Raw intensity measurements at 730 and 850 nm wavelengths and at dark current condition (when no light is shone) are collected by the system serially from each channel at every 100 ms (10 Hz sampling rate). Ten of the twelve detectors are placed 2.5 cm away from each of the four sources, and the remaining two detectors are placed 1 cm away from the middle two sources (see Fig. 2A). With the given source-detector configuration and sensor placement on the forehead, 16 long SDS and 2 short SDS measurements can be collected from the prefrontal areas of the head (see Fig. 2B).

Preprocessing of data for artifact removal and hemodynamic conversion. We used Matlab
(MathWorks, R2019b) to preprocess and prepare data for statistical analysis. Channels whose raw intensity measures were saturated (> 4500), had high dark current values (> 200), or had high correlation with dark current measurements for each channel (r > 0.7), suggesting poor coupling with the skin, were removed from further analysis 76,77 . Following channel rejection, wavelet-based motion artifact removal (the tuning parameter was set to 0.1) was applied to remove abrupt spikes 61 . Then, the modified Beer-Lambert law (MBLL) was used to calculate the changes in the concentrations of HbO and HbR during the task relative to a baseline ( �c HbO and �c HbR ) from the intensity measurements ( I ) at two wavelengths, using Eqs. (1) and (2).
OD is change in light attenuation or optical density at a given wavelength ( ) (730 or 850 nm) during a task condition relative to baseline, ε λ X (ε 730 HbO = 0.390, ε 730 HbR = 1.102, ε 850 HbO = 1.058, and ε 850 HbR = 0.691) is the wavelength dependent molar extinction coefficient of chromophore X (HbO or HbR), d is the distance between source and detector, DPF λ is the wavelength dependent differential pathlength factor that relates the distance that light travels from the source to the detectors with scattering and absorption effects 78,79 . DPF λ was corrected for age 16 . Spline interpolation method was used on HbO, HbR and total hemoglobin (HbT = HbO + HbR) to correct abrupt signal shifts caused by movement of the sensor or the head 80 . Measurements from channels 3 to 6 and 11 to 14 were averaged and labelled as left and right long SDS regions, respectively. Alternatively channels 17 and 18 were used to represent left and right short SDS regions, respectively. Lastly, to ensure reliable extraction of wavelet coefficients, fNIRS measures were detrended by subtracting a third order polynomial fit and bandpass filtered with cut-off frequencies at 0.009 Hz, and 2 Hz as described in 81 . Spectral analysis by wavelet transform. Continuous wavelet transform (CWT) allows for representation of signals in the time-frequency plane. CWT of a signal x(t) can be obtained by using Eq.
is the scaled and translated mother wavelet function.
Similar to other studies that have investigated cardiovascular dynamics using CWT, the present study used Morlet wavelet as the mother wavelet 29,60 . Morlet wavelet is a Gaussian function, modulated by a sine wave, whose time and frequency domain representation are shown in Eqs. (4) and (5), respectively, where f 0 is the frequency and time resolution parameter.
(1) OD = log I baseline I task �c HbR �c HbO

Statistics.
In order to statistically evaluate wavelet coefficients across bands, tissue layers, hemispheres, and stimuli, we used quantitative measures established by 60 . The first quantitative measure is the time-averaged energy density ( ε ) within a given frequency band, which was calculated using Eq. (6) where f 1 , f 2 are minimum and maximum frequencies associated with each band (Ex. Cardiac: f 1 = 0.4 Hz , f 2 = 2 Hz).
Overall magnitude of ε is sensitive to subject-dependent factors, such as activation in hemisphere with task condition. These effects may lead to an underestimation of band and tissue layer effects. Therefore, to mitigate these effects, we calculated relative energy density ( relε ) using Eq. (7) where ε total average is the average energy density across the entire frequency spectrum (0.009 Hz to 2 Hz).
An example of changes in time averaged wavelet coefficients across 0.009 to 2 Hz of long-and short-SDS measurements taken from one subject is shown in Fig. 3A. Strong activity can be seen in the VLF band, followed by myogenic, cardiac, and respiratory bands. The magnitude in the VLF band was greater in long-SDS measurements, while all other bands had greater magnitudes in short-SDS measurements. Apart from the VLF band, local maxima appear to occur at the same frequencies in both long-and short-SDS measurements. Specifically, oscillations appear at 1.14 Hz, 0.19 Hz, and 0.08 Hz of the cardiac, respiration and myogenic bands, respectively. In the VLF band a large peak emerged at 0.02 Hz in the long-SDS measurements, but not in short-SDS measurements. The changing periodicity of the local maxima is visualized in the plotted scalogram (see Fig. 3B).
The results indicate that cardiac and VLF peaks remained stationary for the most part, while respiratory and myogenic peaks varied with time.
We used R (R Core Team, 2019), lme4, lmerTest and emmeans to perform linear mixed effects (LME) analyses [83][84][85] . Model development investigated the main and interaction effects of band (VLF, myogenic,   (H2)) was added as another factor to investigate task-evoked systemic and neuronal effects. Model terms were built on a-priori knowledge and selected based on backward elimination of non-significant predictors from a saturated model 84 . The final model chosen for cognitive and hypercapnic stimuli are coded using lme4 syntax from R and are shown in Eqs. (8) and (9), respectively. Significance of fixed effect terms were evaluated using likelihood ratio tests, where the full effects model was compared against a model without the effect in question (E.g., 1 + Band + (1|Participant) vs 1 + (1|Participant) when evaluating significant effect of Band). Post hoc analysis was conducted to evaluate differences between levels per each model term. Seven planned comparisons were performed when evaluating both stimuli, where three were for "Band" (VLF vs myogenic, respiratory, or cardiac) and four were for "Band: SDS" (e.g., long vs short from VLF). For cognitive stimulus, an additional 56 planned comparisons were performed, where eight were for "Band: SDS : Hemisphere" (e.g. left vs right from VLF of short), and 48 were for "Band : SDS : Condition" (i.e., E1 vs E2, E1 vs E3, H1 vs H2, E1 vs H1, E2 vs H2, and E3 vs H1 from VLF of short) terms. Example of code and resulting output of the likelihood ratio and post hoc tests used are shown in Supplemental Table 1.
Homogeneity of variance and normality tests were conducted using visual inspection of residuals and Q-Q plots for all models. If model predictions showed heteroscedasticity or a non-normal distribution, then log10 transformations were performed on the response variables. However, model estimates ( β ) and confidence intervals ( CI ) were back transformed to original units. Satterthwaite approximation of degrees of freedom was used in post hoc analyses 86 . For all statistical analyses, the level of significance was set at α = 0.05. Adjustments using false discovery rate (FDR) were made on p-values to account for Type I error inflation. Lastly, Cohen's d was calculated per post hoc comparison to examine effect sizes 87 . Cohen's d of 0.2, 0.5 and > 0.8 reflect small, medium and large effects, respectively. Coding used to calculate these effect sizes are shown in Supplemental Table 1.

Results
Study I: cognitive stimulus. Figure 4 shows the comparison of relε from short-and long-SDS channels across all subjects for each biomarker. Over the frequency range of 0.009 to 2 Hz, both SDSs presented similar average wavelet coefficients with noticeable local maxima at ~ 1 to 1.5 Hz, ~ 0.25 Hz, ~ 0.08 and ~ 0.01 to 0.02 Hz in cardiac, respiratory, myogenic and VLF bands, respectively. These local maxima were visible in HbO and HbT, but were not visible in HbR, especially those associated with myogenic and respiration activity.
Statistical evaluation revealed significant effects of band per biomarker for relε (HbO: χ 2 (3) = 1178.0, p < 0.001; HbR: χ 2 (3) = 1178.1, p < 0.001; HbT: χ 2 (3) = 1133.1, p < 0.001) and ε (HbO: χ 2 (3) = 979.8, p < 0.001; HbR: χ 2 (3) = 1170.0, p < 0.001; HbT: χ 2 (3) = 890.7, p < 0.001). Regardless of the biomarker, power was contained primarily in VLF band, followed by myogenic, cardiac, and respiratory bands, as shown by decreasing β values per biomarker (see Estimate Band column in Table 1). Post hoc comparisons between bands revealed significant differences between VLF and other bands (see Post-hoc Band column of   Interaction between band, SDS and condition was significant across all biomarkers for relε (HbO: χ 2 (32) = 69.9, p < 0.001; HbR: χ 2 (32) = 86.1, p < 0.001; HbT:χ 2 (32) = 67.9, p < 0.001) and ε (HbO: χ 2 (32) = 117.2, p < 0.001; HbR: χ 2 (32) = 153.3, p < 0.001; HbT:χ 2 (32) = 116.6, p < 0.001). ε increased across conditions of same difficulty level in VLF, myogenic and respiratory bands for all biomarkers (see Fig. 5 and Supplemental Table 3). Alternatively, ε decreased in cardiac bands of HbO and HbT biomarkers and increased in HbR. Differences between recordings of the third and last easy (E3) and the first hard task condition (H1) was significant across all SDS measurements and biomarkers of VLF, myogenic and respiratory bands, with E3 having higher values. No significant differences were found between conditions of different difficulty (i.e., E1 vs H1 or E2 vs H2). Differences within the same difficulty conditions varied with band and SDS measurement. In the VLF band,  Table 1. Differences in relative energy density ( relε ) values between bands and between SDS measurements for each band per biomarker of cognitive stimulus. β Estimate, CI confidence interval;adj.p FDR adjusted p values with bolded values representing significance at α < 0.05, d Cohen's d. In post-hoc band column, the metrics represent differences between very low frequency (VLF) and other bands (e.g., Myogenic -VLF), where positive d means that the effect was more pronounced in VLF band. Alternatively, in post-hoc SDS column the metrics represent differences between short-and long-SDS measurements, where negative d values mean the effect was more pronounced in long-SDS measurements. Study II: hypercapnic stimulus. Figure 6 shows average relε of long-and short-SDS measurements per biomarker, when subjects were administered with hypercapnic stimulus. Noticeable local maxima at ~ 1.2 Hz, ~ 0.2 Hz, and ~ 0.1 can be seen in cardiac, respiratory and VLF bands of HbO and HbT, respectively. Alternatively, multiple local maxima of small amplitude can be seen in myogenic band.   Table 2), relε was concentrated in VLF band, followed by myogenic, respiratory, and cardiac band for HbR and HbT biomarker. In HbO, relε was concentrated in the VLF band, followed by myogenic, cardiac and respiratory bands. Post hoc comparisons between bands revealed significant differences between VLF and other bands, except for VLF and myogenic band (see Post-hoc Band column of Table 2 Table 2). The respiratory activity was dominant in long-SDS measurements, while cardiac activity was leading in short-SDS measurements.

Discussion
In this study, we utilized wavelet transform analysis and implemented multiple long-and short-source-detector separation (SDS) channels to evaluate and quantify the contribution of fNIRS signal components with respect to physiological origins, tissue layers and hemispheres during cognitive and hypercapnic stimuli. Overall, our results indicate that fNIRS signal components have unique temporal, spatial and frequency responses to different types of stimuli. The following discussion of these results begins by identifying physiological factors affecting fNIRS signals, followed by how much these factors contribute to various fNIRS markers, and how their contributions vary axially, laterally, and temporally with different stimuli.
Characteristic frequencies were observed on average at around 0.02, 0.09, 0.20 and 1.20 Hz (see Figs. 4 and 6). These detected peaks have appeared at similar frequencies in other fNIRS, fMRI and laser Doppler studies, and have been reported to be due to cardiac, respiratory, myogenic, neuronal and partial arterial carbon dioxide (PaCO 2 ) activities 21,88,89 . Specifically, peaks at 0.02 Hz have been attributed to nitric oxide (NO)-related endothelial activity. During increased metabolic demand, release of NO leads to vasodilation of cerebral vessels to allow for constant and optimal nutrient and oxygen supply as well as effective removal of waste products such as CO 2 . Although vasodilation does not occur solely due to NO, studies have reported that NO is responsible for at least 61% of CBF increase 90 . During hypercapnic stimuli, neuronal activity is minimal, therefore the activity within this band can be attributed directly to PaCO 2 -related changes. This change is commonly referred to as cerebrovascular reactivity (CVR) 13,91 . Alternatively, during a cognitive task, since the dilation is being induced by increased neuronal activity, the phenomenon is termed neurovascular coupling (NVC) 92 . However, a cognitive task may elicit both PaCO 2 and neuronal activity related vasodilation 24 . In fact, PaCO 2 -related effects may Table 2. Differences in relative energy density ( relε ) values between bands and between SDS measurements for each band per biomarker of hypercapnic stimulus. β Estimate, CI confidence interval,adj.p FDR adjusted p values with bolded values representing significance at α < 0.05; d Cohen's d. In post-hoc band column, the metrics represent differences between very low frequency (VLF) and other bands (e.g., Myogenic -VLF), where positive d means that the effect was more pronounced in VLF band. Alternatively, in post-hoc SDS column the metrics represent differences between short-and long-SDS measurements, where negative d values mean the effect was more pronounced in long-SDS measurements. www.nature.com/scientificreports/ contribute to at least 1.5% increase in fNIRS biomarkers when cognitive tasks involving out-loud speech or those requiring mental thought are performed 54 . Since the cognitive task we employed does not require internal or external conversation, we assume that the effects of PaCO 2 were minimal and that changes reflected are solely due to NVC. Local maxima at ~ 0.09 Hz have been reported to occur due to a combination of Mayer waves and vasomotion. Mayer waves are associated with arterial pressure and sympathetic response, while vasomotion is defined as the oscillation in the tone of blood vessels that is independent of respiration, pulsation and neuronal activity 21,88,93 . Due to high correlation and coherence between vasomotion and Mayer waves, the exact physiological phenomenon behind the observed peaks is difficult to isolate 19 . Peak frequencies at ~ 0.20 Hz and ~ 1.20 Hz correspond to respiration rate and heart rate, respectively 21,29,88,89 . Contribution of each of the aforementioned factors to the overall fNIRS signal varied with stimuli. Specifically, during cognitive stimulus the VLF band had the highest activity followed by myogenic, cardiac, and respiratory bands (see Table 1). Alternatively, during hypercapnic stimulus, activity was highest in the VLF band followed by cardiac, respiratory, and myogenic bands (see Table 2). Observed contributions of VLF band during both stimuli are in agreement with other spectral analysis studies 29,64,[66][67][68][69] , and form the physiological foundation behind fNIRS usage to study NVC and CVR. However, the effect size associated with this band was equal in hypercapnic (HbO: d = 0.58; HbR: d = 0.97; HbT: d = 0.55) and cognitive tasks (HbO: d = 0.59; HbR: d = 1.27; HbT: d = 0.55). This finding is in conflict with a transcranial doppler study, which reported that CO 2 increased cerebral blood flow velocity (CBFV) by 25% while NVC only increased it by 5% 94 . On the other hand, the study utilizing BOLD-fMRI and CBF-ASL, reported comparable increases in measures due to a cognitive task (complex Stroop) and hypercapnic tasks in a young cohort similar to the findings of our study 95 . Differences in myogenic band contributions to fNIRS signals based on stimuli administered are likely due to Mayer waves being attenuated during hypercapnia, but not during cognitive tasks 22 . In fact, in cognitive studies, contributions due to myogenic activity have been shown to reduce the accuracy of extracted neuronal related hemodynamic effects by at least threefold 82 . The respiratory peak had the weakest presence in cognitive task and a noticeable presence in the hypercapnic task. Such changes in respiration activity during cognitive and hypercapnic tasks on fNIRS signals have been reported previously 22 . During the hypercapnic stimulus, increased PaCO 2 levels stimulate central chemoreceptors to increase the rate and depth of breathing in order to return CO 2 concentrations to normal resting levels 96 . Alternatively, during functional activation as evoked by cognitive tasks, changes in PaCO 2 are typically not as drastic, therefore not stimulating chemoreceptors and resulting in weaker respiration activity 22 . Cardiac peaks had the second largest contribution to fNIRS signal during hypercapnic stimuli, but less so in cognitive task. This may be due to the cardio-respiratory coupling being stronger during hypercapnia stimuli than during cognitive stimuli, however further research needs to be conducted to investigate the reasoning behind this finding.
Systemic activity dominated short-SDS measurements, while neuronal-related activity dominated long-SDS measurements in the cognitive task (see Post-hoc SDS columns of Table 1). Such differences were not significant in the hypercapnic task (see Post-hoc SDS columns of Table 2), however trends observed in Fig. 6, indicate that respiratory activity was more dominant in short-SDS measurements, while PaCO 2 and cardiac activity were stronger in long-SDS measurements. Magnitude of the activity within the VLF band from short-SDS measurements was not negligible for either stimulus. In fact, the band power was greater than that of other systemic bands in both stimuli (see Estimate short-SDS columns of Tables 1 and 2). Furthermore, the VLF band power from long-SDS measurements were greater than that of short-SDS measurements (see Post-hoc SDS columns of Table 1) in the cognitive task but did not differ in the hypercapnic task (see Post-hoc SDS columns of Table 2). Studies have indicated that such activity is likely due to scalp blood flow (SBF) [34][35][36][37][38][39][40][41] . Based on this, the results in this study indicate that measurements taken from long SDS channels reflect a combination of hemodynamic changes associated with NVC and SBF in the cognitive task and strictly SBF in the hypercapnic task. The results from the cognitive task are in agreement with other studies, where activity within VLF band of short-SDS channels have been reported to not only share similar spectral characteristics as that of long-SDS channels, but also contribute at least 80% of the signal detected by at the long-SDS channels 22,33 . However, presence of SBF during the hypercapnic task is surprising, since studies have reported that PaCO 2 does not induce SBF changes 23,97 . A possible explanation underlying such findings may be that SBF changes in extracerebral layers are temporally delayed, rather than minimal or not present 98,99 . In cognitive tasks, systemic activities such as myogenic, respiratory and cardiac oscillations were found to have higher contributions in the short-SDS measurements than the long-SDS measurements (see Post-hoc SDS columns of Table 1). These findings have been reported by other studies, where higher coherences were observed between physiological signals collected from peripheries (blood pressure, respiratory or heart rate) and short-SDS measurements than long-SDS measurements 29,38 . In hypercapnic tasks, respiratory activity was found to have a significantly higher presence in short-SDS measurements, while myogenic activity had no difference across SDS measurements and cardiac activity having higher presence in long-SDS measurements.
Activity varied laterally across hemispheres during the cognitive task, but not during the hypercapnic task. In fact, hemisphere was not a significant predictor for the hypercapnic task. This is supported by CVR literature, which suggests that PaCO 2 is a global stimulus that induces similar magnitude of change in CBF across all regions of PFC to a fixed input 71 . For the cognitive task, significant differences in activity were observed only in VLF band of long-SDS measurements. Firstly, no significant difference observed between hemispheres for other bands of long-SDS measurements and for all bands of short-SDS measurements is in accordance with numerous cognitive and motor studies, which have reported homogeneity of systemic activities and extracerebral fNIRS signals between right and left middle frontal areas 29,42,43,47 . Secondly, this finding verifies that the activity observed and reported in our prior publications was mainly due to task-related neuronal activity and not due to task-related extracerebral activity 100 www.nature.com/scientificreports/ A task-moderated effect was seen within the VLF, myogenic and respiratory bands in both long-and short-SDS measurements during the cognitive task (see Fig. 5 and Supplemental Table 3). Specifically, regardless of SDS channel, VLF band activity significantly increased across tasks of the same difficulty (E1 vs E3 or H1 vs H2) but not across tasks of the same order and different difficulties (E1 vs H1 or E2 vs H2). Differences across harder tasks (H1 vs H2) had lower effect sizes than those associated with easier conditions (E1 vs E3). Such within difficulty differences in long-SDS measurements can be explained by focus, while the across difficulty can be explained by failure to meet demands posed by the more difficult task condition 102 . Specifically, if the task is too taxing, then the likelihood of the subject to make errors and more specifically disengage with the task is high 103 . Based on this, we can speculate that some subjects stopped engaging with the harder task conditions, therefore resulting in an overall lower average change across participants. Preliminary evidence of this disengagement has been shown in a previous study utilizing the same cognitive task,where behavioral outputs became worse with the harder task and were positively correlated with decreased hemodynamic changes 104 . Alternatively, similar results from short-SDS measurements may be attributed to stress relief 28,29,105,106 . These effect sizes of differences between conditions were greater in magnitude in short-SDS measurements than long-SDS measurements, suggesting a large task-evoked systemic contribution. Additionally, more task condition pairs were found to be significant in short SDS measurements, indicating the need to remove not just global scalp activity, but also task-specific systemic activity before deriving subject-disengagement related interpretations. Similar task-evoked activity as seen in the VLF band, was seen in the myogenic and respiratory bands. One interpretation of the results is that increased neuronal activity may have collectively elicited myogenic responses and the intrinsic load elicited by lack of a break given between conditions of the same difficulty may have led to over breathing, elevated consumption of oxygen and release of CO 2 . Interestingly, a significant decrease was observed from E3 to H1 in both long and short-SDS measurements across all bands. This decrease in activity may be explained by the presence of a fifteen-minute break between the two difficulty conditions. The presence of the break may have allowed for re-establishment of baseline/resting conditions prior to beginning of hard task block.
In comparison to HbO and HbT, HbR had the largest magnitude within the VLF band for both tasks. HbR has been demonstrated to be a more reliable biomarker for neuronal activation than HbO 28,107 . This statement is further supported by the fact that during cognitive tasks HbR has higher correlation with BOLD-fMRI than HbO or HbT 108 . Higher correlations between HbR and BOLD-fMRI during CO 2 tasks have also been reported 71,99 . However, due to its low dynamic range, HbR is more prone to motion artifacts and therefore results in an underestimation of hemodynamic activity associated with NVC or CVR 30 . The data utilized in this study had relatively high signal to noise ratios, therefore we predict that the high magnitude of peaks associated with HbR are reflective of neuronal activity and not noise. HbO and HbT indicated stronger presence of systemic factors than HbR. This is because oscillations in cardiac output induce changes in blood volume of the arteries, but not in the veins 34 . Thus, the HbR component of the fNIRS signal is generally less affected by physiological activities other than the cerebral function 28,64,107 . However, HbR is not devoid of systemic effects (see Figs. 4 and 6). Therefore, HbR should not be left uncorrected 63 . Systemic activity related to SBF was present in all biomarkers. However, unlike many studies, which have reported that SBF effects are minimal in HbR 28,41,106 , we did not observe such trend. In fact, we observed smaller differences between short and long-SDS channels, suggesting a larger extracerebral contribution to HbR signal. This finding is relatively new but has been supported by recent papers 63,109 . However, application of short SDS did not improve extraction of NVC related HbR signals 41,44,106 . Studies have reported that this may be due to a contribution that is not related to SBF 41,63 . Therefore, further investigation regarding systemic confounders on activity within the VLF band of HbR needs to be done.
Despite the promising methodology, this study results are subject to a few limitations. Due to the short duration of the hypercapnia experiments, wavelet coefficients below 0.02 Hz could not be extracted with high accuracy. Therefore, findings within VLF bands need to be further investigated with a longer protocol. The short protocol length also prevented evaluation of effects of hypercapnic vs normocapnic conditions on physiological causes, since splitting would lead to a reduction in data length to 3 min, which would prevent extraction of wavelet coefficients below 0.06 Hz. As previously mentioned a protocol of at least 15 min is needed to extract wavelet coefficients as low as 0.009 Hz 81 . Furthermore, only data from five participants were analyzed for the hypercapnic task, therefore the findings reported here are preliminary in nature. SBF and myogenic activity have been reported to be heterogeneous within extracerebral layers of frontal areas 29,42,47 . Specifically, medial frontal regions consist of scalp draining veins which lead to increased SBF effects than lateral regions. The instrumentation that was used in our study did not allow for investigation of lateral extracerebral regions. The presented CO 2 and cognitive studies were employed in two different groups of subjects; therefore, we are not able to relate individual contributions of PaCO 2 and VLF-induced hemodynamic changes. Hence a future study investigating both NVC and CVR within the same subject population needs to be performed. We did not collect data during the resting period, which is needed if we are to quantify non-evoked versus evoked activity within each band. We realize that the contributions observed may be a result of the long-term correlated 1/f fluctuations and therefore the physiological origin of the energy density at a certain band may be more complicated than the generally accepted and implemented separate bands in fNIRS applications as in this study 110 . To further refine the understanding of the effect of the physiological causes on fNIRS signals, in the future we plan on subtracting the spectral component of 1/f fluctuations from the frequency spectra and re-assessing with the proposed methodology. We realize that the long-SDS measurements are a combination of extracerebral and cerebral activity. Therefore, to truly assess contributions arising from each tissue layer requires separation of signals from cerebral and extracerebral layers by implementing one of the conventionally used short SDS regression methods on long-SDS measurements which we plan to perform as a future work. Lastly, to further investigate task-evoked effects, verify frequency bands in different systemic bands and identify subject's state (in terms of stress, perceived workload, etc.), we plan to conduct similar analysis between activated and non-activated channels, coherence analysis between fNIRS and www.nature.com/scientificreports/ other sensors data (e.g., electrocardiography, laser doppler flowmetry) and correlational analysis between fNIRS and survey results (e.g., state trait anxiety, NASA-TLX).
In conclusion, this study indicated that spectral (physiological causes), spatial (axially-tissue layers and laterally-regions of middle frontal areas in left and right hemispheres), and temporal (task conditions) characteristics of fNIRS signals vary with the type of stimuli administered. Furthermore, we utilized continuous wavelet transform and employed multiple long-and short-source SDS channels that allows for identification of the source in signal variability, in turn, enables determination of the appropriate signal processing techniques needed to extract the signal of interest. Lastly, this study also indicates that task-evoked systemic and extracerebral effects may serve as supplementary signals of interest regarding the state of the participant while engaging with a real-world task.

Data availability
The data that supports the findings of this study are available from the corresponding authors, P.R and K.I, upon request.