Frequency‐domain analysis of fNIRS fluctuations induced by rhythmic mental arithmetic

Abstract Functional near‐infrared spectroscopy (fNIRS) is an increasingly used technology for imaging neural correlates of cognitive processes. However, fNIRS signals are commonly impaired by task‐evoked and spontaneous hemodynamic oscillations of non‐cerebral origin, a major challenge in fNIRS research. In an attempt to isolate the task‐evoked cortical response, we investigated the coupling between hemodynamic changes arising from superficial and deep layers during mental effort. For this aim, we applied a rhythmic mental arithmetic task to induce cyclic hemodynamic fluctuations suitable for effective frequency‐resolved measurements. Twenty university students aged 18–25 years (eight males) underwent the task while hemodynamic changes were monitored in the forehead using a newly developed NIRS device, capable of multi‐channel and multi‐distance recordings. We found significant task‐related fluctuations for oxy‐ and deoxy‐hemoglobin, highly coherent across shallow and deep tissue layers, corroborating the strong influence of surface hemodynamics on deep fNIRS signals. Importantly, after removing such surface contamination by linear regression, we show that the frontopolar cortex response to a mental math task follows an unusual inverse oxygenation pattern. We confirm this finding by applying for the first time an alternative method to estimate the neural signal, based on transfer function analysis and phasor algebra. Altogether, our results demonstrate the feasibility of using a rhythmic mental task to impose an oscillatory state useful to separate true brain functional responses from those of non‐cerebral origin. This separation appears to be essential for a better understanding of fNIRS data and to assess more precisely the dynamics of the neuro‐visceral link.


| INTRODUCTION
Based on the neurovascular coupling principle, functional near-infrared spectroscopy (fNIRS) aims at detecting the hemodynamic changes evoked by neuronal oxygen consumption. NIRS is a non-invasive optical imaging technology and it has been widely used to measure brain, mostly cortical, activity through relative concentration changes in oxygenated (HbO) and deoxygenated (HbR) hemoglobin (for reviews see: Obrig & Villringer, 2003;Pinti et al., 2020).
A major challenge in fNIRS research is to reliably disentangle the hemodynamic response due to neurovascular coupling from other confounding components . fNIRS changes caused by brain activity are naturally low in amplitude and unfortunately also overlap with other fluctuations that do not originate in the cerebral cortex, mainly: (i) systemic hemodynamic activity detectable in both cerebral and extracerebral regions (Bauernfeind et al., 2014;Minati et al., 2011;Tachtsidis et al., 2009), (ii) local blood flow changes in superficial tissue layers across the head (Kirilina et al., 2012), and (iii) instrumental noise and other artifacts. The first two, far from being simply spontaneously generated, can also be evoked by cognitive, emotional or physical tasks. If the modulation of these non-cortical task-related components of the signal mimics the dynamics of the brain activation of interest they could become an important source of interference and noise (Nambu et al., 2017;Näsi et al., 2013;Zimeo Morais et al., 2017). So much so that Takahashi et al. (2011) showed that the task-related skin blood flow changes could explain over 90% of the NIRS signal on a verbal fluency experiment, while Minati et al. (2011) further demonstrated the strong confounding effect of arterial blood pressure fluctuations.
To better infer the presence of a functional response, experimental protocols attempt to increase statistical power by repeating the stimuli a sufficient number of times, interspersed with contrast conditions in which a different response (or none) is expected. To this end, fNIRS experiments often used blocked or event-related designs, depending on whether one wants to analyze sustained or transient responses, respectively (Pinti et al., 2020). Event-designs use short-duration stimuli, normally randomized in order and separated by a constant or jittered inter-stimulus interval. Block-designs attempt to maintain mental engagement by presenting stimuli within a condition for a long enough time interval, followed by a different condition or a resting inter-stimulus interval (Amaro & Barker, 2006). To investigate the interaction between "sustained" and "transient" responses, mixed designs can also be used (Petersen & Dubis, 2012).
Depending on the stimulus presentation strategy, different analysis methods have been developed to make inference about the functional hemodynamic response and isolate it from confounding interferences (for a review see Tak & Ye, 2014). Although classic averaging strategies provide robust results, the usual averaged-based statistical tests, such as t test or ANOVA, do not allow estimating the shape or time-course of fNIRS signals, so they have been progressively replaced by more powerful methods. These include the general linear model (GLM) framework (Friston et al., 2007;Schroeter, Bücheler, et al., 2004), data-driven approaches as principal component analysis and independent component analysis (Kohno et al., 2007;Zhang et al., 2005) and dynamic state-space modeling Kolehmainen et al., 2003). GLM is one of the most widely adopted statistical framework to quantify how well the measured fNIRS signals fit a hemodynamic model that reflects the expected neural response. It exploits the good temporal resolution of fNIRS and allows to include different covariates within the regression model (e.g., physiological signals). In its most basic form, the model is obtained by convolving a hemodynamic response function (HRF) with a stimuli function that encodes the hypothesized time course of the neuronal response (Koh et al., 2007;Tak & Ye, 2014). Therefore, GLM is a hypothesis-driven approach that requires the combination of a specific HRF (often taken from fMRI studies) and other nuisance regressors to construct the linear model, which might not be obvious depending on the task type, brain region and participant's idiosyncrasy. Moreover, GLM demands special caution when applied to fNIRS signals due to some statistical issues (Huppert, 2016;Huppert et al., 2009;Koh et al., 2007). In contrast, "principal" and "independent" component analysis methods only rely on general statistical assumptions as orthogonality and independence, respectively. Although useful for separating the mixed components that make up the fNIRS signals, they require additional processing to elucidate which of them are task-related and which are not, particularly difficult when extracerebral and cerebral responses are correlated (Zhou et al., 2020). State-space models, mainly based on the Kalman filter, allow building complex hemodynamic models to describe the time varying characteristics of the fNIRS signals and estimate the HRF. Although dynamic analysis appears to provide better estimates of the HRF and better account for non-stationary signals, it still requires improvements in model specifications and statespace estimators.
Regardless of the strengths and weaknesses of each experimental method, all benefit from the inclusion of short-distance recordings to obtain a reference of the contribution of superficial layers to fNIRS signals (for reviews see: Fantini et al., 2018;Tak & Ye, 2014). Multi-distance measurements are considered particularly effective in isolating the actual cerebral response. Nevertheless, some open questions remain regarding, for example, the ideal range for source-detector distances, the optimal number of short-channels and their arrangement relative to long-channels. Ideally, and because there is growing evidence of the heterogeneous nature of surface hemodynamics (Wyser et al., 2020), each long-channel should be paired with at least one nearby short-channel. Unfortunately, such a precise spatial configuration of the paired measurements is not currently possible with the most commonly used NIRS devices today.
Our interest is focused on investigating the coupling between hemodynamic changes arising from the superficial and deep layers during mental effort, and isolating the task-evoked cortical response from other confounders. The rationale behind is the assumption that extracerebral and cerebral responses are the result of a coordinated effect product of different interrelated processes (Caldwell et al., 2016;. Rather than viewing surface fluctuations only as annoying confounders that must be removed, we regard them as carriers of valuable information. Information, that might prove essential not only to gain a better understanding of the fNIRS data but also, and perhaps as important, to more accurately assess the full dynamics of the neuro-visceral link. This goal, however, is hampered by the inherent difficulty in differentiating spontaneous from task-locked fluctuations. One possible solution is to deliberately induce periodic oscillations, so that they can be easily located and analyzed by using well-established frequency-domain methods. In previous studies, a number of protocols have been used to generate hemodynamic oscillations at particular frequencies to investigate cerebrovascular regulation, including cuff inflation-deflation (Aaslid et al., 2007;Kainerstorfer et al., 2014), head-up tilting (Hughson et al., 2001), squat-stand exercises (Claassen et al., 2009), paced breathing (Pierro et al., 2014), and even in visual and motor studies (Obrig et al., 2000;Schroeter, Bücheler, et al., 2004;Wolf et al., 2002).
Here, we hypothesized that performing a cyclic cognitive task would also induce periodic hemodynamic fluctuations measurable in fNIRS recordings. To test this hypothesis, the current study was designed to generate an oscillatory state suitable for effective analysis in the frequency-domain. To this end, we employed a mental arithmetic task within a cyclic block-design at a specific frequency, while performing dense multi-distance recordings on the forehead using a newly developed multichannel NIRS device. Concurrently, we recorded cardiac activity by continuously monitoring heart rate.
The analysis of the magnitude and phase relationships between signals in the frequency-domain would allow to: (i) identify common task-related oscillatory activity, (ii) estimate the contribution of shallow and deep tissue layers to fNIRS signals, (iii) separate task-related surface hemodynamics from the putative cortical response, and (iv) measure the relative timing between HbO and HbR changes to better interpret the underlying physiological processes. Furthermore, we used the empirical transfer function as an alternative method to estimate the functional brain activity and assess the timed-coordination between extracerebral and cerebral responses. To our knowledge, such an approach has never been previously tested in fNIRS studies.

| METHOD
All data processing was done off-line with MATLAB (Version R2019a, Mathworks, Natick, MA, USA), using native functions, self-made scripts and open source packages.

| Participants
All procedures performed in this study were approved by the Ethics Committee of the University Miguel Hernandez, in accordance with the declarations of Helsinki. Participants did not receive any remuneration. A total of twenty-four healthy young adults volunteered, 10 males and 14 females (mean age: 22.3; SD: 4.2), were recruited for this study. All of them were instructed beforehand about the purpose of the experiment and provided informed written consent prior to study enrollment. After that, the participants practiced the task for 10-15 min to make sure they understood and got used to it (to minimize stress responses). They were seated in a comfortable position while performing the task.

| Mental task
In this work, the task was based on a block protocol designed as a cyclical pattern of mental effort, alternating phases of mental math with phases of pause of the same duration, that is, regular repetitions of activation-rest. The idea behind this was to induce periodic hemodynamic changes in the form of cycles of some kind of response followed by a return to basal levels. In this way, such an oscillatory pattern may be analyzed by conventional spectral methods. As illustrated in Figure 1, the experimental session was organized into three consecutive uninterrupted recordings: (i) 300 s of baseline in resting condition, (ii) 300 s of task, and (iii) 300 s of recovery in a relaxed state. Participants were asked to keep their eyes on the screen throughout the experiment. During specific time intervals, an image was presented on the screen to instruct the participants to relax mentally. This image was a soft-colored paint depicting an inanimate scene (it was rated as neutral by the participants in a parallel study with affective pictures; results not shown here). The NIRS probe remained in the same position throughout the session. The task consisted of 10 consecutive 30-s trials. Each trial began with 15 s of mental calculation, followed by a 15-s pause of relaxation. To perform the mental math participants were asked to iteratively subtract a small number (between 5 and 9) from a three-digit number (between 100 and 199), as fast as possible. Both numbers, chosen randomly in each trial, were presented on a 21.5″ display monitor, 80 cm. away from the participants' eyes. Afterwards, the pause started by presenting the question "Result?" for 5 s, which prompted the participants to inform verbally of the final result of their mental calculations (to allow scoring the performance and ensuring that the participants were paying attention), followed by the "relax" image. Two seconds before presenting the subtraction operands, a fixation cross was displayed in the middle of the computer screen to announce the beginning of the mental calculation. Note that the 30-s period of the trials corresponds to a frequency of 0.033 Hz, which we will refer to as f t (task frequency) throughout this manuscript. This frequency was chosen so that it did not overlap with well-known spontaneous fluctuations such as blood pressure (0.08-0.12 Hz) (Huppert et al., 2009;Julien, 2006), or very slow endothelial activity (0.01-0.02 Hz) (Stefanovska et al., 1999). Furthermore, the 15-s duration of mental effort accommodates that of a typical hemodynamic response (Buxton et al., 2004;Friston et al., 1998), while the next 15-s pause allows a return to baseline levels, being an optimal inter-event interval to minimize overlaps between consecutive hemodynamic responses (Aarabi et al., 2017;Dale, 1999). Here, we tried to minimize the generation of distress by using a subtraction arithmetic task that was cognitively challenging, but did not exceed the participants' mental abilities. Further, we constantly emphasized the importance of the mental effort, and not the amount or accuracy of the operations performed.

| Heart rate measurement
ECG was registered using a BIOPAC MP36 physiological monitoring system and the AcqKnowledge software 4.1 (Biopac Systems, Inc., Goleta, CA, USA) at a sampling rate of 500 Hz. ECG was recorded in lead II configuration with disposable electrodes. The MP36R device was digitally synchronized through the I/O port with the PC running the stimuli presentation, so event markers were recorded as well. At the end of the experiment, the raw data were post-processed with AcqKnowledge to: (i) create an R-R tachogram from the ECG signal, using the implemented Pan-Tompkins algorithm (Pan & Tompkins, 1985) for R wave detection, and (ii) extract the instantaneous heart rate from the reciprocal of the tachogram. Finally, the data was exported to MATLAB and then resampled to 10 Hz using cubic spline interpolation.

| Heart rate exclusion criteria
Many studies have investigated how cognitive performance correlates with the stress level induced by mental workload, in most cases through different physiological measures, heart rate being one of common use (Charles & Nixon, 2019;Hakimi, 2018;Mandrick et al., 2016;Tao et al., 2019). However, there is not a clear threshold to differentiate between heart rate changes due to a pure mental effort and those due to a stressful situation. As the present work focuses on the first, it was necessary to estimate a maximum increase in heart rate, beyond which the influence of stress was considered disproportionate. We took as a reference the results of other studies that applied the Trier Social Stress Test, a standard protocol for stress induction in healthy people (Kudielka et al., 2007). In their review work, Kudielka et al. (Kudielka et al., 2007) reported that the mean heart rate increases to the test are about 15-25 bpm. Kirschabum et al. (1993), found increases about 26 bpm during the test. On the basis of the aforementioned literature, we decided to apply a limit well below the reported values. Compared to the mean heart rate at baseline, the threshold was set to a maximum F I G U R E 1 Schematic representation of the experimental procedure increase of 12 beats per minute (bpm) during the task. Two participants were excluded for exceeding this limit, leaving a sample of N = 22. Thus, only those participants showing reasonably stable heart rate throughout the baseline and task periods were further considered. The Wilcoxon-signed rank test was used to resolve for differences between the maximum heart rate reached during baseline and task.

| fNIRS recordings
In this study we used a multichannel, wireless, continuouswave NIRS device (Brainspy28, Newmanbrain, S.L., Elche, Spain), which employs four sources and ten detectors forming a rectangular grid of 80x20 mm. Each source housed two light-emitting-diodes (LED) at wavelengths 740 nm and 850 nm. Through a precise switching cycle, the device combines pairs of optodes at different separation distances, providing 16 short-channels and 12 longchannels that corresponds to a source-detector distance of 14 and 32 mm, respectively ( Figure 2a). Moreover, the device measures and corrects the ambient light contribution.
Also it incorporates a 3-axis accelerometer to account for head motion. It transfers data wirelessly (via Bluetooth) at a sample rate of 10 Hz. The NIRS probe was placed onto the forehead, centered on AFpz according to the international 10-5 system, mainly covering the frontopolar area of the prefrontal cortex (PFC) (Figure 2b). The optodes contact the skin through an intermediate convex lens pressing the skin when the probe is held firmly, in order to reduce cutaneous blood flow and, therefore, its hemodynamic interference (Takahashi et al., 2011).

| Signal quality check-Channels and participants exclusion criteria
To ensure that only clean signals pass to further analysis, we conducted some preliminary quality tests. To account for instrumental noise (Huppert et al., 2009;Orihuela-Espina et al., 2010) we evaluated first the raw optical data to identify channels exhibiting extreme values (bellow 5% or above 95% of the device dynamic range), or an excessive coefficient of variation >7.5% (calculated as the percentage ratio between the standard deviation and the mean) F I G U R E 2 Probe geometry and placement. (a) Optode arrangement (yellow squares: Sources; white circles: Detectors). Sixteen short-channels (black lines with numbers) and twelve longchannels (green lines with letters). (b) Probe position on the forehead. Green shaded areas roughly mark the regions explored by the long-channels. Black rectangles outline the long-channels averaged within each ROI (right, medial and left) (Zimeo Morais et al., 2017). Thus, we discarded the recordings suffering from poor signal to noise ratio, saturation or unphysiological noise contamination. By visual inspection, we rejected recordings affected by motion artifacts by identifying sharp changes in the fNIRS signals aligned with abrupt shifts in the accelerometer data. Two participants were excluded after the quality check. Therefore, the final sample consisted of N = 20 participants (8 males and 12 females; mean age: 21.7; SD: 3.8).

| fNIRS data preprocessing
Preprocessing was carried out using the Homer2 NIRS package (Huppert et al., 2009) based in MATLAB. The raw optical data were converted to optical density, and then into oxy-([HbO]) and deoxyhemoglobin ([HbR]) relative concentration changes via the modified Beer-Lambert law (Delpy et al., 1988;Kocsis et al., 2006). We used a differential pathlength factor calculated in accordance with the general equation described in Scholkmann and Wolf (2013), which takes into account the participant's age and the wavelength. No partial volume correction was used.
The fNIRS data were digitally low-pass filtered by using a zero-phase, 5th-order Butterworth filter, cut-off 2.5 Hz (MATLAB Signal Processing Toolbox); no high-pass filtering was applied. Therefore, we remove only highfrequency noise while preserving cardiac, respiratory, blood pressure and vasomotor components. After preprocessing, for each chromophore we obtained 16 time-series from the short-channels plus 12 from the long-channels that we refer to as "shallow-" and "deep-signals", respectively, to indicate how much the light penetrated during their corresponding records. Thus, we had a set of multidistance recordings to effectively address the problem of extracerebral contamination (Pfeifer et al., 2018;Saager & Berger, 2005;Scholkmann et al., 2014;Yücel et al., 2015), by assuming that short-separation recordings are sensitive only to extra-cerebral changes, while long-separation recordings are sensitive to both extracerebral and cerebral activity (Brigadoi & Cooper, 2015;Saager & Berger, 2005;Scarpa et al., 2013;Yücel et al., 2017).
When using multi-distance recordings regression can be performed assuming that the physiological noise has comparable time courses in both shallow-and deepsignals, and that the focal, task-evoked, cerebral hemodynamics is independent, that is, uncorrelated (Saager et al., 2011;Saager & Berger, 2008). However, the heterogeneous nature of the local superficial fluctuations cannot be dismissed, which raises the need to collect the shallowsignals at recording sites as close as possible to the deepsignal to be decontaminated from surface hemodynamics (Gagnon et al., 2012). Here, we used a NIRS device that allows each deep-signal to have three shallow-signals that meet the proximity requirements: two obtained close to the long-channel's detector and source, respectively, and one close to its center (see Figure 2a). Adopting the "double short separation measurements" approach recommended in (Gagnon et al., 2014), we used the combination (the sum) of the two shallow-signals recorded near the detector and source to estimate the extracerebral component to be suppressed from the corresponding deep-signal. Thus, for example, the signal from long-channel "A" was regressed on the sum of the signals from the short-channels 1 and 5 ( Figure 2). Over the entire time courses, for all the deep-signals we computed: where 0 and, 1 are the regression coefficients, S deep is the deep-signal, S shallow is the combined shallow-signal and S clean the desired "clean signal" (in fact, the raw residuals). We solved linear regression by applying the MATLAB function "robustfit", which uses an iteratively reweighted least squares algorithm and is less sensitive to outliers than ordinary least-squares (Holland & Welsch, 1977). After regression, each S deep has its associated pair S clean and S shallow , making a total of 12 signal triplets for each chromophore. Noteworthy, although S clean probably represents the neural component, at this point we prefer the term "clean" without making assumptions about its actual origin. This precaution is based on the fact that, because our mental task attempts to induce periodic oscillations the requirement of non-correlation between extracerebral and cerebral hemodynamics may not be met, compromising the performance of regression Saager & Berger, 2008). Later, in this work, we will apply additional analysis to verify the nature of these regression-estimated signals.
Finally, as the channel positions are not fully consistent across individuals due to the variability in head shape and size (Tak et al., 2016), the signal-to-noise ratio and signal reliability can be improved by spatial clustering (Plichta et al., 2006;Schecklmann et al., 2008). To this end and to avoid interpreting isolated channels, for every single participant we averaged across the signals belonging to three regions of interest (ROI), left, medial, and right ( Figure 2b). As such, the signal triplets associated with the four leftmost long-channels (A, B, C, D) yield the corresponding averages of S shallow , S deep and S clean . The same procedure was applied for the medial (E, F, G, H) and the rightmost (I, J, K, L) four long-channels. Thus, each of the three ROIs now reduces to just three average signals that, from now on, we denote by SS (shallow), DS (deep), and CS (clean). All further processing was done on these signals, which display comparable signal-to-noise ratio across the three ROIs because they have been obtained by averaging the same number of neighboring signals in each frontopolar region. For each ROI, we computed the averaged time courses of HbO and HbR for the three signals across all the participants. Furthermore, we obtained the grand average of the trials epoched in the interval −15 to 30 s relative to the onset (10 trials × 20 participants = 200 trials). Standard error of the mean (SEM) was calculated as well.

| Identification of task-induced frequencies
Firstly, we evaluated if the task successfully induced oscillations by estimating the power spectral density (PSD) of each participant's NIRS data during the task and the baseline condition. The PSDs were computed using Welch's averaged periodogram method (Welch, 1967) with a hamming window of length 2000 samples and 50% of overlap to achieve a good frequency resolution (0.005 Hz), spectral smoothness, and reduced noise variance (Ilvedson, 1998). To allow comparisons, the PSDs were normalized to relative percentage values by calculating the power ratio of each frequency bin to the total power of the entire spectrum (Aarabi & Huppert, 2016). For further analysis, we focused on the frequency range 0.005 to 0.08 Hz under the assumption that the expected task-induced oscillations would fall within that range.
We assessed the presence of significant task-induced oscillations by contrasting the PSDs during task and baseline. Along frequency bins we performed one-tail paired t-tests to check whether within-subject PSD values were statistically higher for the task compared with baseline. The observed t-statistics were corrected for multiple comparisons by following the cluster-based nonparametric approach given in (Maris & Oostenveld, 2007). We computed Monte Carlo cluster tests over 1000 permutations of the same t-test where the condition, that is, task and baseline, was randomly shuffled within subject. Then, we estimated the so-called permutation p-value from the proportion of random realizations that have larger cluster-statistic than the observed one. We set a critical alpha-level = 0.01 to identify the frequency bins that showed significantly higher PSD values during the task. Finally, the PSDs were averaged across participants (N = 20) to obtain the average normalized PSD of HbO and HbR for each signal type and ROI. In addition, the 95% confidence interval (CI) for the mean at each frequency bin was calculated by bootstrapping over 1000 resamples.

| SS and DS relationships
Here, SS data were obtained from fNIRS channels exploring the frontopolar region with a source-detector distance of 14 mm. The scalp-cortex distance is known to be increased in this region (15 mm to 17 mm) compared to more lateral frontal areas (Cui et al., 2011;Haeussinger et al., 2011) which, together with the presence of the frontal sinuses, decreases the cerebral fNIRS sensitivity (Haeussinger et al., 2014). Moreover, as pointed out by Zhang et al. (2015), shortchannels with a source-detector distance in the range of 14 to 16 mm have a sensitivity to the brain of only about 0.47%. Although it is unlikely that SS picked up cortical signals, it was still important to examine the commonalities and differences between SS and DS to reasonably ensure that SS data are primarily dominated by shallow hemodynamics, while DS also contain deeper components that, likely, stem from the cortical layer. The rationale for this analysis was to demonstrate that the link SS-DS would be altered at the task frequency if another deep oscillatory process (e.g., neurovascular response) appears in the DS signal. To this end, for the time-series pairs HbO ss /HbO ds and HbR ss /HbR ds we performed the following within-subject analysis: (i) crossspectrum to identify shared fluctuations, (ii) frequencydomain correlation to identify significant covariation, and (iii) transfer function to evaluate the relationship in magnitude and phase. Following the recommendations given in (Claassen et al., 2015) for transfer function analysis, neither detrending nor high-pass filtering was used, and a triangular smoothing window (coefficients ¼, ½, ¼) was applied to both auto-and cross-spectra.
To identify synchrony (or shared fluctuations) in the frequency-domain between shallow and deep signals, a cross-spectrum analysis was performed to compare the two signals. We computed the cross-power spectral density (CPSD) of the bivariate time-series using the MATLAB function "cpsd", based on the Welch's averaged periodogram method (Welch, 1967); as for PSD we set a hamming window of length 2000 samples and 50% of overlap. From the complex-valued result, we obtained the magnitude and phase values to find the shared power and phase shift between both signals at particular frequencies.
The magnitude peaks that showed significantly higher values during the task were located by following the same approach as for PSD in Section 2.5; the values were normalized and averaged in the same way.
To ensure the reliability of the cross-spectral estimates, we evaluated whether the signals showed significant stability in their relative amplitude and phase at particular frequencies. To this end, we estimated their frequencydomain correlation by computing the magnitude-squared coherence (MSC) as a function of the PSDs and the CPSD (Zhang et al., 1998): Coherence values vary between 0 (no correlation) and 1 (perfect correlation), describing the linearity of the relationship between both signals in the frequency domain. The coherence values were also averaged across participants for each signal pair and ROI. Only if significant coherence exists, the cross-spectral estimates have a useful meaning (Claassen et al., 2015). Values of about 0.5 have been commonly considered as a threshold for significance  but, as the threshold may depend on the frequency, a better alternative for assessing significance is the surrogate data approach (Faes et al., 2004;Paluš, 1997). In this work, following the nomenclature given in (Paluš, 1997), we chose the FT1 surrogates method as it conserves the power spectrum of the original signals (Faes et al., 2004). FT1 surrogates were constructed by substituting the phase of the Fourier transform of the original signals with random values in the range [−π, π] while the modulus is preserved, and then returning to time-domain by inverting the Fourier transform. We generate 1000 surrogate series and for each one the averaged MSC was computed as for the original signals, realizing the null hypothesis that the averaged MSC stem from pairs of signals that fluctuate asynchronously at the same frequencies of the experimental signals. Afterwards, the mean MSC values obtained from the actual signals were compared with the FT1 distribution to estimate the MSC threshold levels for each frequency bin at α = .05. A transfer function describes the dynamic relationship between the output signal of a system and the input signal. Under the assumption of linearity, the transfer function can be estimated from the frequency-domain representation of the experimental input-output signals. Transfer function models have become a popular approach to investigate the dynamic of cerebrovascular autoregulation (Claassen et al., 2015;Van Beek et al., 2008), and they have also been used to remove systemic physiological noise from fNIRS signals (Bauernfeind et al., 2013;Florian & Pfurtscheller, 1997). In this work, assuming that SS has energy in the frequency range of interest and contain quasi-periodic oscillations, the transfer function H (ƒ) was approximated from the experimental fNIRS data as (Zhang et al., 1998): For each time-series pair, shallow and deep signals were, respectively, the input and output data used to obtain an approximation of the transfer function at particular frequencies (e.g., input HbO ss , output HbO ds ). From the complex-valued result, we obtained the magnitude (gain), which represents the relative change in μM between input and output, and the phase that carries their temporal coupling (phase difference or time-lag), while coherence values indicate the reliability of these measures. For reporting, the gain data were converted into percentage values. Then, at the group level and for each signal pair and ROI, we computed the averages of gain and circular phase angle; phase statistics were managed by means of the Matlab Toolbox CircStat (Berens, 2009). In addition, we calculated the 95% confidence interval (CI) band around the mean by bootstrapping over 2000 resamples. Finally, to assess whether the task induced a consistent phase coupling across participants, we applied a bootstrapped Rayleigh test (Oden, 1983) on the phase values at the task frequency. Consistent coupling should be reflected as a narrow distribution of phase values around a preferred angle. On the contrary, a poor inter-subject synchronization should display a more uniform distribution throughout the 360° circle.

| Estimating deep component from transfer function
We previously expressed some concerns relating to the reliability of the regression-estimated CS signals. To address them, we applied an alternative novel approach to estimate CS from the transfer function. Supposing that DS solely contains the same fluctuations seen in SS, the relation between both signals should be highly linear and coherent, only altered by the small differences in the volume sampled by our short-and long-channels. Therefore, the transfer function should yield fairly constant gain values across frequencies (depending of the fractional part of DS power that is produced by SS) and a phase shift close to zero (in-phase). Furthermore, coherence values should be close to one. However, when another oscillatory process is added to DS we can expect some level of disturbance in gain and/or phase at certain frequencies. At a specific frequency, cyclic hemodynamic oscillations can be approximated as sinusoids, completely defined by their values of amplitude, phase and frequency. Adding two sinusoids of common frequency results in a sinusoid with same frequency but with amplitude and/or phase altered. In our case, if the task elicits independent cyclic oscillations in both shallow and deep layers, the observed deep sinusoid X ds (t) at the task-frequency would result from the sum of the two contributing sinusoids X ss (t) and X uk (t), that is: where X ss (t) is the observed shallow sinusoid and X uk (t) is the unknown deep component. All of these sinusoids are characterized by their values of amplitude A, angular frequency and phase , thus, in sinusoidal form: Because they represent harmonic oscillations of common frequency (i.e., ds = ss = uk = task frequency), they do not depend on or t, but only on A and , which makes it possible to convert them into phasors (or complex numbers). In exponential form would be: �� ⃗ X ss and �� ⃗ X uk are the phasor representation of X ds (t), X ss (t) and X uk (t) respectively. According to phasor algebra: hence, we can estimate the unknown component by computing the subtraction of the two complex numbers corresponding to the known phasors, that is,: To this end, we firstly obtained the phasors' values for phase and amplitude as follows: (1) We designated �� ⃗ X ss as the "reference phasor" and, hence, ss = 0. (2) From the transfer function we obtained the phase value between SS and DS at the task frequency f t , so ds = arg H f t , which indicates the phase-shift of �� ⃗ X ds with respect to �� ⃗ X ss . (3) A ds was obtained from the PSD by calculating the RMS amplitude at f t and converting it to the peak amplitude of a sinusoid that is, where FR is the frequency resolution (0.005 Hz in our case). Due to spectral leakage, better amplitude estimates are obtained by summing the PSD values within the frequency interval f t − FR, f t + FR . (4) A ss was estimated in the same way and then scaled to the theoretical value that it should reach in DS by itself. This is a crucial step in the procedure. The scaling factor was the gain value of the transfer function at f t but during "baseline", which represents the fraction of SS magnitude present in DS when no significant deep component contributes. Thus, A ss was multiplied by the baseline gain to yield its scaled amplitude. Then, we performed the phasor subtraction (Equation 1) to obtain the amplitude A uk and the phase uk of the phasor �� ⃗ X uk . Finally, �� ⃗ X uk was converted to a sinusoid, which represents the alternative CS estimated by transfer function. This procedure was performed for every single-participant's signals and averages were computed for each ROI and chromophore.

| Simulations
To assess the feasibility of CS estimation by phasors, we performed a fairly realistic simulation by using the actual data during task rather than artificial or baseline signals. SS was used as is. DS was obtained by scaling SS to the amplitude that it should reach in the deep-recording, that is, multiplying by the transfer function gain during baseline (See Section 2.7, step 4). We generated a synthetic neural signal as a sinusoidal wave oscillating at 0.033 Hz, with amplitude and phase obtained from the averaged CSs estimated by regression (amplitude = 0.05 μM and 0.04 μM, time-lag = 13 and −9.5 s for HbO and HbR, respectively, and relative to CS) (See Results 3.1 and Figure 4). Then, the sinusoid was added to DS to build the simulated deepsignal. This procedure was applied to the medial-ROI data of every participant and then subjected to both regression and phasor estimation of the neural component. Because phasor-estimated result represents the average over the task, the regression-estimated signal was averaged across trials to allow intra-subject comparisons. The quality of the neural estimates was quantified using the root mean square error (RMSE) between the true synthetic signal and the recovered one. A paired t test was applied to resolve for statistical differences.

| Measuring HbO/HbR coupling
To investigate significant frequency-domain correlation between HbO and HbR, we computed their MSC and CPSD during the baseline and task conditions. In this case, we were only interested in the SS and CS signals because we assume that DS is nothing more than the linear combination of the first two. MSC and CPSD data, as well as their averages and significance thresholds, were obtained by the same procedure previously detailed. In addition, from the complex-valued CPSD we extracted the phase data (Müller et al., 2003;Reinhard et al., 2006) at the task-frequency to assess the temporal relation (timeshift) between both chromophores. As usual, this procedure was applied for every single participant. Next, we computed the circular mean of phase angles at the group level and then, to assess whether it was significantly oriented in a preferred direction, we applied a bootstrapped Rayleigh test (Oden, 1983) through 2000 resamples. The 95% CI was also computed by bootstrapping.
2.10 | Measuring HbO/heart rate coupling As in Section 2.9, we investigated the relationships between HbO and heart rate by MSC and CPSD. We only used the HbO data of the SS signals because the coupling between heart rate and the rest of signals can be inferred from the results obtained in other previously performed tests.
A ds cos ds t + ds = A ss cos ss t + ss + A uk cos uk t + uk .
�� ⃗ X ds = A ds e j ds , �� ⃗ X ss = A ss e j ss and �� ⃗ X uk = A uk e j uk ,

| RESULTS
We analyzed how fNIRS signals fluctuate in response to a cyclic mental arithmetic task, which can induce physiological stress responses (Charles & Nixon, 2019;Kudielka et al., 2007). To reduce any putative stress response that could have influenced our results we have taken several cautions. First, we have used a non-strenuous task duration, only 30 s, including a calculation and a recovery phase. Second, volunteers practiced the task for 10-15 min to make sure they understood and got used to it and the experimental setting. Third, the participants included in this study had a stable mental performance and reasonably low heart-rate fluctuations during the entire experimental session. In fact, they showed no significant differences between their heart rate values across the two conditions (baseline and task) of the experiment, (Wilcoxon's signed rank test, p = .091). On average, during task, the mean heart rate increased only slightly, 7 bpm (9%), with respect to the baseline (Table 1).

| Oscillations alignment to the task frequency
fNIRS signals showed clear oscillations in accordance with the task frequency and such oscillations were consistently observed, over multiple ROIs, on both shallow and deep layers. Figure 3 shows the results of the PSD analysis at the group level. The normalized PSDs were averaged across participants for each of the three ROIs, for each signal type and for each condition (baseline and task). The DC component (frequency = 0) was set to zero and the upper value of the displayed frequency range was limited to 0.08 Hz. A ubiquitous peak can be seen at ~0.01 Hz both at rest (Figure 3, black traces) and during mental task (Figure 3, colored traces), although the cluster-permutation test found no significant differences between both conditions. However, a clear peak at f t was detected exclusively during the task and marked as significant by the cluster-test (p < .01) (Figure 3, shaded vertical rectangles). This peak corresponds closely to the task frequency and was observed in all ROIs and signal types. In some cases, a secondary peak around 0.066 Hz was also found, which most likely represents a harmonic of the fundamental task-frequency. It seems obvious that the task successfully induced cyclic fluctuations of the HbO and HbR, which were present in the shallow-signals ( Figure 3SS), deep-signals ( Figure 3DS) and in cleansignals ( Figure 3CS). Figure 4 shows the group-level averaged temporal traces of HbO and HbR after band-pass filtering around the task frequency, using a filter width of 0.015 Hz. It can be seen how SSs ( Figure 4SS) start to oscillate in sync with the task-trials, although showing an evident time-shift between both chromophores. Figure 4 also depicts the grand averages of the trials across participants (smaller plots next to temporal traces), which for SS display a common response consisting of: (i) shortly before the trial onset, the HbO strongly increases reaching a maximum at ~11.2 s, (ii) then HbO gently returns to previous levels during the subsequent rest, and (iii) HbR changes were less pronounced and lead HbO by ~4 s. Regarding the DSs, a similar periodic response can be observed ( Figure 4DS), but HbO peaks are slightly anticipated to SS (~ 0.3 s) while HbR lags SS by ~2 s. Finally, CSs also show cyclic fluctuations aligned to the task ( Figure 5CS), but in this case HbO oscillates in counter-phase with respect to SS. Furthermore, HbO and HbR show similar amplitudes (~0.05 μM) and evolve anti-correlated with each other, HbO showing a valley at ~12.6 s and HbR a peak at ~18.5 s. At this point, it seems likely that the observed time-shift in DS, relative to SS, is due to the summation of the CS component. Interestingly, the pattern seen in CS would correspond to an inverted HbO/HbR response, that is, decrease in HbO together with an increase in HbR. Noteworthy, the inverted pattern does not appear immediately but progressively reaches stability during the first few trials. This observation might be of interest for a tentative physiological interpretation (see Discussion).
Because fNIRS signals can be highly individual-specific, we have provided additional information about hemodynamic response and spectral data for each participant in the Supplementary Material, with the aim of illustrating the individual differences and the variable contribution of surface tissues to the fNIRS signals in different regions of the forehead.

| SS/DS relationships
To answer the question of to what extent shallow fluctuations contribute to deep-signals, we performed coherence, cross-spectra and a transfer function estimate between SS and DS data ( Figure 5). The analysis revealed a significant cross-spectral peak during the task centered at f t (p < .01), and that was present in all ROIs for HbO and HbR (Figure 5a, left axis). These peaks indicate that SS and DS oscillate at the task frequency with a remarkable shared power. Moreover, a higher peak was located at ~0.01 Hz but showing no differences during task compared to baseline. Coherence values were above the significance threshold along the entire range of explored frequencies in both conditions (Figure 5a, right axis), describing a consistent linear relationship between SS and DS in the frequency domain. Coherence levels were particularly high for HbO (on average > 0.8), indicating a stronger correlation and, likely, that HbO is more influenced by shallow hemodynamics than HbR. It should be noted that coherence also peak during the task around f t (quite evident for HbR), meaning that the task-induced oscillations were more coherent than the spontaneous ones.
During the task, transfer function analysis revealed a decrease in magnitude around f t for HbO, especially pronounced in the medial ROI (Figure 5b). Phase values also showed certain perturbation at the same frequency (Figure 5c). HbR exhibited a slight increase in magnitude, mainly observed in the medial and left ROIs (Figure 5b), in parallel with strong disturbances in the phase values (Figure 5c). Again, the resting condition did not show such changes. A t test revealed no differences at the group-level between gain values of HbO and HbR during rest (p = 0.2). However, during task the gain values were significantly lower for HbO (p < .01) in medial and left ROIs, which agrees with the decreasing gain values shown in Figure 5b. The paired t test showed no differences for either HbO or HbR when comparing task with rest, indicating that on average the change during task is very subtle.

F I G U R E 3
Grand average (20 participants) of the normalized PSDs for each ROI (right, medial and left), chromophore and signal type in the frequency range 0.005 to 0.08 Hz. Solid black curves refer to baseline and colored curves to task. Dashed lines represent the 95% CI of the mean. Shaded rectangles delimit the frequency ranges that show significantly higher power during task. (SS) results for shallow-signals. (DS) deep-signals. (CS) clean-signals F I G U R E 4 Averaged time courses of HbO (red traces) and HbR (blue traces) across participants during task (plus a portion of baseline and recovery to the left and right respectively) for each ROI and signal type. The small plots next to time courses show the grand average of trials. SEMs are depicted by thin lines. Gray boxes mark the 15-sec of mental math of each trial. Each row show the results for shallow (SS), deep (DS) and clean-signals (CS) for the three ROIs (right, medial and left). Circular histograms show that during rest the phase angles concentrate around 0° at the task frequency in all the cases (Figure 5g), being the Rayleigh test highly significant (p < .01) and the 95% CI quite narrow. This implies that SS and DS oscillate almost in-phase at that frequency (no time-lags). However, during the task the phases of HbR were shifted clearly to positive values (DS lags SS) in the medial and left ROIs by 22° and 30°, respectively (Rayleigh test, p < .01), which correspond to time-lags of 1.8 and 2.5 s at f t (Figure 5d). Phase changes in HbO were less apparent than those of HbR, a slight shift to negative values (DS leads SS) was observed in the same ROIs (−8° and −4°, time-lags 0.6 and 0.3 s). Noteworthy, the time-lags calculated from the phase values coincide with those obtained from the averaged time-series (see Figure 4).
As expected, these results corroborate that shallowand deep-signals are highly correlated, underlining the strong influence of surface hemodynamics on deep recordings, which is particularly true for HbO (also observable within a single individual, compare Figures S1 and S2 of supplementary material). Fortunately, the analysis also disclosed changes in magnitude and phase related to the task, pointing to the contribution of a deep component, uniquely present in DS, as responsible for the observed disturbances. Therefore, these findings confirm that our deep recordings captured other oscillatory processes that are different from the superficial ones.

| Neural signal estimation by transfer function
We applied an alternative method to estimate the putative neural signal by using phasor representations of the magnitude and phase data obtained from transfer function. Under the assumption that signals were quasi-stationary over the sort period of time determined by the task, we modeled them as sinusoids oscillating at 0.033 Hz, that is, the task frequency. Phasor algebra was used to extract the hidden deep component that explained the disturbances observed in DS. As a representative example, Figure 6 compares the CS obtained by regression with that estimated from transfer function data in the medial ROI. As reported in Section 3.2, the transfer function of HbO showed a clear decrease in magnitude around f t in parallel with a phaseshift of about −8°, which is illustrated in the middle-plot of Figure 6a After performing the phasor subtraction, the interfering deep component emerged as the sinusoid depicted by the red trace in the plot to the right of Figure 6c. In the case of the HbR, a small increase in magnitude and a phase-shift of about 22° were measured at ~f t . Similarly, the right-hand plot of Figure 6b illustrates the differences between the predicted (solid line) and observed (dashed line) DS. Phasor algebra pointed to the blue sinusoid drawn to the right of Figure 6c as responsible for the disturbances. The left-hand plot of Figure 6c shows the CSs obtained by means of regression. It can be seen the good match with the deep components estimated by phasors. For the three ROIs (right, medial and left), the estimated time lags between SS and CS were on average 11.5, 13 and 10.7 s for HbO, and −9.9, −9.5 and −9.5 s for HbR ( Figure  S6). Phasor analysis corroborates that the task induced deep hemodynamic fluctuations in form of an inverted HbO/HbR response, which from now on we will consider as the neural signal we were looking for. Despite individual differences, accounting for the inter-subject variability, we found that 90% of the participants (n = 18) showed this type of response in at least one ROI (see Figure S3). Figure 7a shows a segment of the bandpass filtered SS, simulated DS and synthetic neural response for each chromophore of a representative participant. Simple visual inspection reveals time courses comparable to the groupaveraged experimental signals shown in Figure 4. The recovered signals averaged across all participants are shown in Figure 7b for each estimation method. In the case of HbO, it can be seen that the shapes of the recovered signals are very similar for both methods, but the phasors fit better to the true synthetic signal showing significantly lower RMSE values (paired t test, p < .01). HbR time courses were also comparable, but again phasors performed better (p < .01).

| Simulation results
These findings suggest that, compared to phasors, regression underestimates amplitude, mainly for HbO. Here, we found higher SS-DS coherence for HbO than for HbR, indicating a strong correlation that could affect more regression performance. Overall, phasors seem to work as well as or even better than regression, and could help to independently verify the results in any case.

| HbO/HbR coupling
During the task, CPSDs showed a significant peak of shared power at f t in both SS ( Figure 8SS, left axis, red trace) and CS ( Figure 8CS, left axis, red trace) over all ROIs. MSC measures revealed significant frequency-domain correlation levels at several frequencies, as indicated by the asterisks along the MSC curves drawn in the right-axis of the spectral plots. It can also be seen that MSC levels were higher during the task (red trace) than at rest (black trace), suggesting that the task increased coupling not only in its frequency but also beyond. Since Interpreting coherence when spectral power is very low could be risky and results outside the task specific frequency were beyond the scope F I G U R E 7 Construction of synthetic data and simulation results. (a) Simulated DS time courses for HbO (left) and HbR (right) data of a representative participant. DSs (dashed traces) were constructed by adding a synthetic neural signal (magenta traces) to the theoretical scaled version of the observed SS (solid traces in red and blue). (b) Averaged recovered signals across all participants for HbO (left) and HbR (right) by regression (green) and phasors (magenta), compared to the true synthetic signal (black) F I G U R E 8 Averaged CPSDs, MSCs and phases between HbO and HbR for each ROI, comparing baseline and task. (SS) results for shallow-signals. CPSDs are drawn on the left axis for baseline (black curves) and task (red curves), dashed lines depict the 95% CI of the mean. Gray boxes indicate the frequencies that showed significantly higher cross-power during task. MSCs are shown on the right axis; asterisks indicate frequencies with significant coherence. Circular histograms show phase differences at f t during baseline (top) and task (bottom); significant phase concentration statistical values are labeled in magenta text. (CS) results for regression-estimated clean-signals of the present work, here we focus solely on the peak of interest.
MSC also peaked around f t , reaching the highest levels in the medial and left ROIs (MSC > 0.9, p < .01), which suggest strong correlation between HbO and HbR. Concerning temporal coupling for SS, circular histograms show how phase differences at f t significantly concentrate during task around −39°, −39° and −47° in the right, medial and left ROIs, respectively (Rayleigh test, p < .01), with HbR leading HbO by 3.2 to 3.9 s ( Figure 8SS, circular plots at bottom). For baseline data, phase analysis was not significant for any ROI ( Figure 8SS, circular plots at top). Notably, CS revealed that, on average, HbO and HbR oscillate almost in phase-opposition in the medial and left ROIs (−162° and 176°, respectively, p < .01), showing a time lag around 15-s but with an unclear precedence. In the right ROI the significance level was not reached in any case. These results imply that the task induced coherent fluctuations between HbO and HbR in both shallow and deep layers, but with some differences. For SS, the temporal coupling was consistent across participants over all the ROIs, showing an outof-phase relationship. However, CS phases were consistent at the group level only in the middle and left regions, moreover showing counter-phase fluctuations, in line with an inverted hemodynamic response.

| HbO/heart rate coupling
CPSD and MSC measures were also performed to investigate HbO and heart rate couplings. We focused solely on significant peaks in shallow signals. As illustrated in Figure 9, a peak of shared power was located in the three ROIs at f t . In parallel, MSC also peaked at the same frequency, indicating strong correlation at the task frequency (p < .01). Notably, the phase mean did not reach significant levels in any case (p > .05). These results indicate that at the single-subject level HbO and heart rate oscillate well coupled, resulting in significant coherence at the averaged group level. However, individual phase values were different enough to disperse the angular mean, reflecting an evident inter-subject variability in the temporal coupling of the two signals (see Figure S1). Figure 9 also shows the averaged fluctuations of the band-pass filtered HbO (bottom plots, red trace) and heart rate (bottom plots, black trace). The grand average of trials across participants (small plots next to time courses) also revealed that both signals start to increase few seconds before the trial onset.

| DISCUSSION
The main purpose of this work was to assess the feasibility of using rhythmic cognitive tasks to induce periodic hemodynamic fluctuations suitable for effective frequencyresolved measurements. First, we investigated whether power spectral analysis can distinguish task engagement from rest. Next, we measured the phase-amplitude coupling between different signal pairs to estimate their linear relationship, aiming to differentiate the functional brain response from extra-cerebral confounders and to infer the nature of the underlying processes.
F I G U R E 9 Averaged time courses, CPSDs, MSCs and phases between superficial HbO and heart rate (HR) for each ROI, comparing baseline and task. Top row plots show the CPSDs (left axis), MSCs (right axis) and circular histograms of phase differences at f t . Baseline values are drawn in black and task values in red. Gray boxes delimit the frequencies with significant higher cross-power during task. Bottom row plots show HbO (red) and heart rate (black) time courses along the entire task; gray boxes mark the 15-s of mental math of each trial. The small plots next to time courses show the grand average of trials. Thin lines depict SEMs We found that mental arithmetic successfully evoked cyclic changes in NIRS signals measured in the frontopolar region of the forehead, in the form of highly characteristic spectral peaks centered on the task frequency (i.e., 0.033 Hz). These peaks are clearly discernible from spontaneous activity in resting-state and apparently independent of a stress response. As mental arithmetic tasks have been used frequently as stressors (Al-Shargie et al., 2016;Hakimi, 2018;Takamoto et al., 2013), and mental overload increases stress (Mandrick et al., 2016;Tao et al., 2019), the experimental protocol used here was tailored to minimize such an effect.
Predictably, we corroborated the strong influence of surface hemodynamics on deep fNIRS signals, as reflected by the high coherence levels found between signals obtained from multi-distance recordings. In this study, rather than simply considering the surface contribution as unwanted noise, we took advantage of the shared oscillatory state imposed by the task to extract amplitude and phase data, and use it to separate the true deep signal from that originating in extracerebral tissues. Notably, our results revealed that the deep signals extracted follow an unusual pattern of HbO decrease accompanied by HbR increase, which is reversed (i.e., symmetrically opposite) with respect to the canonical brain activation response.
Previous works have reported that fNIRS signals oscillate close to the experimental stimulation frequency. Franceschini et al. (2000), for instance, using a motor task comprising sequences of 10 s of tapping and 17 s of rest, found spectral peaks at the frequency of the 27-s task period (0.037 Hz). Schroeter, Schmiedel, et al. (2004), exploring the visual cortex, identified spectral peaks for HbO and HbR at 0.023 Hz, close to the 0.028 Hz corresponding to the visual stimulation cycle of 35 s (18-s of stimulation plus 17-s of rest). Likewise, Zhang et al. (2007) also found a peak in visual cortex related to the stimulation frequency of 0.033 Hz. Regarding mental tasks, Kirilina et al. (2012) observed periodic changes in fNIRS, fMRI and other physiological signals, coupled to the single block period of an n-back (30-s) and a semantic task (34-s). Also using a sematic categorization task, Kirilina et al. (2013) identified coherent oscillations between skin blood flow and HbO corresponding to the 34-s period of stimulation. Nonetheless, with the exception of Zhang et al. (2007), none of the cited studies used task/rest intervals of the same duration. We argued that a block-design consisting of exactly regular cycles would better induce an oscillatory state, stationary enough throughout the duration of the task to allow reliable measurements. Although no consensus exists as to the most appropriate stimulus interval in block-design experiments, fNIRS studies often fall in the range of very low frequencies (0.02 to 0.08 Hz), referred as "activation-band" by Kirilina et al. (2013).
However, for frequency-analysis purposes, large task cycles lead to extremely low frequencies, which are difficult to identify and could overlap with spontaneous very slow waves (Stefanovska et al., 1999). For example, Vermeij et al. (2014), using a verbal n-back working-memory task of 180 s (0.005 Hz), only reported fluctuations in the range 0.02 to 0.07 Hz, correctly concluding that they could not be attributed to the task cycle. Likewise, Obrig et al. (2000) used similarly large periods of 120 s, finding peaks at 0.1-0.04 Hz, far from the 0.008 Hz predicted by the task frequency. On the contrary, short periods could fall within the range of spontaneous blood pressure waves (i.e., Mayer waves) that could obfuscate (or override) the functional response (Yücel et al., 2016). Furthermore, it might also be desirable to choose periods shorter enough to accommodate a single response (i.e., not several successively overlapped). Here, we used repetitive 15-s cycles of mental math plus 15-s pause, highlighting the oscillatory activity of 30-s period. We successfully found significant task-locked oscillations, separable from spontaneous activity, and showing time courses compatible with isolated, single responses. Further research is needed to test other stimulation periods that might be even more appropriate.
In almost all the aforementioned studies, the authors discussed the contribution of non-neural components to the observed fNIRS changes, stressing the importance of separating brain activation from these potential confounders. This concern should be particularly addressed when exploring the frontopolar region due to the influence of task-related skin blood flow changes, mainly on HbO signals (Haeussinger et al., 2014;Sato et al., 2013;Takahashi et al., 2011). It seems that, in most functional experiments, the superposition of extracerebral and cerebral hemodynamic responses could not be avoided, as they are not independent but inter-related processes (Caldwell et al., 2016;. Our results validate these considerations by showing that, on the frontopolar region, deep-recordings are strongly influenced by superficial activity. In fact, by solely analyzing deep-signals without applying proper corrections, an activation response evoked by the task can be erroneously deduced. Such a misinterpretation is more likely when only HbO is taken into account, for example as seen in Figure 4DS where HbO shows a clear increase/decrease pattern locked to the task. However, the HbR time course is less conclusive, not supporting a typical activation response and underscoring the need to assess both chromophores to convey more realistic interpretations . Our results also provide complementary evidence that surface hemodynamics influences the different hemoglobin species in a differential way, with HbO being more affected than HbR. This is in agreement with previous studies reporting a stronger influence of confounding factors on HbO (Bauernfeind et al., 2014;Gagnon et al., 2011;Heinzel et al., 2013;Kirilina et al., 2012).
To address the contamination issue, we employed a multi-distance approach in which each deep-recording was cleaned from surface influence by using two specific short-recordings as reference, thus controlling for inhomogeneous hemodynamic activity in the scanned surface area. As commented in Section 2.4.2, such a strategy is commonly accepted as very effective to remove noise, and even considered as a gold standard (Zhou et al., 2020). After performing conventional linear regression, we obtained a clean deep-signal for each ROI that also oscillates at the task frequency, but apparently showing deactivation instead of an expected activation ( Figure 4CS). In view of this unusual pattern, and despite the thickness of the tissues that cover and protect the brain (see Section 2.6), it was necessary to verify whether our short 14 mm recordings were also picking up cortical signals, leading to poorly estimated regression residuals. By exploiting the induced oscillatory state, we employed the empirical transfer function and phasor representation of HbO and HbR fluctuations ( Figure 6) to explain the amplitude and phase disturbances observed in deep-signals. The results confirmed the presence of deep HbO and HbR components that fluctuate coupled to the task, following an inverted pattern with both chromophores almost in counter-phase. The convergence of two independent analysis, regression and transfer function, strongly supports the finding of a putative cortical response in the form of deoxygenation/ oxygenation cycles; noteworthy, it was found in 90% of our sample, a considerable higher occurrence than those reported in, for example, motor imagery fNIRS studies (Abdalmalak et al., 2020;Holper et al., 2011). Therefore, we demonstrated that our short-and deep-recordings captured different hemodynamic components.
Moreover, simulated data revealed that combining transfer functions and phasors might provide better estimates of the amplitude of functional cortical responses. It is worth noting that if shallow and neural signals had oscillated highly correlated (positive or negatively), the regression could have failed because the scaling factor (beta) would have caused the subtraction to flatten the residuals, leading to accidental removal of the cerebral signal. This drawback could be avoided using the method proposed here. Phasors were first proposed by Zheng et al. (2010) to explain the phase relationships between hemoglobin species but, to our knowledge, this is the first time they have been used to analyze multi-distance fNIRS recordings. Further improvements on the method are currently under way.
Inverse oxygenation responses have been previously reported in fNIRS studies using different task modalities, such as motor imagery (Abdalmalak et al., 2020;Holper et al., 2011), visual stimulation (Maggioni et al., 2015), working memory n-back (Haeussinger et al., 2014;Kirilina et al., 2012), emotional stimulation , and mental arithmetic (Bauernfeind et al., 2008;Pfurtscheller et al., 2010). This phenomenon has only been partially explained so far, being of considerable interest to understand the underlying neuro-vascular mechanisms and gain insight into the negative BOLD response observed in fMRI studies (see Holper et al. [2011] and Maggioni et al. [2015] for in-depth discussions of possible explanations). Regarding fNIRS studies probing the PFC during working memory tasks, some researchers found a decrease in HbO in frontopolar region using fixed inter-optode distances (i.e., no multi-distance correction), while they found no significant HbR changes (Haeussinger et al., 2014;Kirilina et al., 2012). They suggested that task-evoked sympathetic vasoconstriction drives skin blood flow changes, which in turn impair fNIRS long-recordings, leading to an apparent decrease in oxygenation response. Noteworthy, in both studies, fMRI data showed deactivation in the medial region of the PFC. In contrast, Takahashi et al. (2011) found a positive pattern (i.e., increased oxygenation) using a verbal-fluency task, which they attributed to skin vessels dilatation by comparing short-(5 mm) and long-distance (30 mm) signals and showing that the effect disappears when pressure is applied to the skin. In the present work, we used two short-channels to clean each long recording, which is a very effective method to remove the components (of local or systemic origin) that are common to the shallow and deep-signals . However, if a systemic component appears at different times in the surface and deep layers (e.g., blood flow delay), the regression itself may render a false neuronal response. Time differences in vascular reactivity driven by the task (e.g., delayed autonomic mediation) may also cause delays between layers. Wyser et al. (2020) reported an average time lag of ~0.51 s between shallow and deep signals for Mayer pressure waves, similar to the values found by Kirilina et al. (2013). Tong and Frederick (2010) estimated a time of ~6 s for a pressure wave to pass the whole brain. Using nose tip temperature as a proxy to assess autonomic activity during a mental math task, Pinti et al. (2015) found mean time lags of less than 4 s, albeit showing significant individual variability, between changes of cutaneous blood flow and prefrontal fNIRS signals. During a cognitive task, a mean time lag of ~6 s between skin blood flow and prefrontal fNIRS signals was found in (Kirilina et al., 2013). The phasor method do not solve the problem of different timing either, but they provide independent information about the delays between signals and avoid regression constraints. We found clearly longer delays than those we have mentioned (see Section 3.3 and Figure S6), and we believe that they are unlikely to be due to a delayed physiological response, otherwise they would reflect a considerable lack of coordination between autonomic and prefrontal activity. Assuming that the observed response is an artifact due to time delays, it can be expected that for a specific individual the delay will be the same in the three ROIs of the narrow frontopolar region explored. However, we found that in some participants showing consistent surface signals across ROIs, the estimated neural signal was clearly different depending on the region (e.g., participants 1, 2, 3 and 4, Figures S1 and S3). Therefore, we conclude that such an inverted pattern actually represents the functional cortical response. However, the existence of some delay cannot be fully ignored and follow-up work is required.
Overall, our results are more in line with those reported by Pfurtscheller et al. (2010). Using a similar task, they also found an HbO decrease in the medial area of the PFC, but no significant changes for HbR. They corrected for extracerebral contamination by employing a common average reference spatial filter, which relies on subtracting the averaged signals from each fNIRS channel. After correction, we found significant fluctuations for both HbO and HbR, which is of great interest for interpretation purposes and implies that our methodological approach is better suited to discriminate the functional brain response. Other fNIRS studies have reported reduced activation in medial areas of the frontal cortex related to arithmetic subtraction in adolescents (Artemenko et al., 2018) or to task difficulty (Verner et al., 2013). However, since they did not use multi-distance corrections, comparing results could be risky.
Currently, fNIRS inverse oxygenation response remains an open topic, as the underlying neurovascular mechanisms are only partially understood. Some studies have been devoted to relate it with the fMRI negative BOLD response (NBR) in visual (Maggioni et al., 2015) and motor cortices (Abdalmalak et al., 2020), investigating its potential cause. Maggioni et al. (2015) found a consistent spatial correlation between NBR and inverse fNIRS response, while Abdalmalak et al. (2020) attribute the phenomenon to motion artifacts. The first discuss indepth on the origin of NBR, favoring the idea of neural deactivation (Mullinger et al., 2014) against the "blood stealing" explanation promoted by other authors (Shmuel et al., 2002). Pfurtscheller et al. (2010) also explained the inverse response in terms of a "focal activation/surround deactivation" pattern. As we did not find any delay in the hemodynamic response, our results also speak against a sequestration of blood from neighboring areas to active areas. Our NIRS probe covered a relatively small area of the PFC, so our data are insufficient to support or reject the idea of concurrently activated and deactivated areas. However, based on the delay need to reach consistent counter-phase changes in HbO/HbR and the strong rhythmicity imposed by the task, we suggest a simpler explanation. It may be plausible that such an inverted response actually expresses a cyclical brain activation state, which after reaching stability appears as a period of oxygen consumption (during mental effort) in a previously well oxygenated brain area, followed by a subsequent reoxygenation (during the pause). Thus, in the steady-state, the mental effort period starts under a condition of O 2 excess (or compensated) that leads to use the currently available O 2 until a new supply of fresh blood is needed. In this line, Wylie et al. (2009) proposed that different HbO/HbR (and total Hb) combinations might be present in activated visual cortex areas. Nevertheless, more research is needed to fully elucidate the complex coupling between O 2 consumption and blood flow/volume changes in the PFC.
We also found that only the medial and left ROIs showed a consistent inverse response at the group level, which is supported by the significant concentration of HbO/HbR phase angles around 180° ( Figure 8CS) and the stronger disturbances in gain and phase detected by transfer function during the task ( Figure 5). This finding suggests that, on average, frontopolar activity was slightly lateralized to the left. A meta-analyses conducted by Arsalidou and Taylor (2011) indicates that, among others, frontopolar area seems to be generically engaged in mental arithmetic, sustaining working memory functions that are necessary to achieve good mathematical performance. They proposed that, for calculation tasks, this area manages the successive executive steps that can lead to the final calculated result. Noteworthy, they also reported that activity in the left part of frontopolar cortex was concordant among studies involving calculation tasks. The characteristics of our arithmetic task fit with that functional specialization as it is necessary to coordinate each of the iterative subtractions (steps), while holding the result of the previous operation in the working memory, and then combining both to obtain the final result (main goal). Nonetheless, as we did not control for task difficulty or used calculation modalities other than subtraction, we acknowledge that complementary research is needed for more rigorous comparisons and interpretations. In any case, an important point to be also considered is whether the observed response can be attributed to the "task-negative" activity of the default-mode network (DMN) (Raichle et al., 2001;Raichle & Snyder, 2007). As the brain region interrogated by our NIRS probe overlaps with the medial prefrontal cortex (DMN-associated region), our findings could reflect decreased neuronal activity in this key region due to task engagement. Since DMN also include deeper structures, unreachable for fNIRS, exploring this possibility will require further research using fMRI imaging methods.
Concerning shallow hemodynamic, we found a consistent task-locked pattern of monotonic HbO increase followed by a decrease during the task pause, which was shortly preceded by similar HbR changes of lesser amplitude. As previously discussed, these fluctuations greatly impair the deep-recordings and potentially lead to their misinterpretation as brain activation, especially when only accounting for HbO. Therefore, we suggest special caution should be taken when interpreting HbO changes reported by fNIRS studies employing mental arithmetic tasks and using only long-recordings, as in (Çiftçi et al., 2008;Tanida et al., 2004;Verner et al., 2013;Yang et al., 2009).
Aiming to elucidate the origin and influence of the extracranial confounds present in fNIRS signals, some studies monitored concurrently skin blood flow, heart rate and arterial blood pressure during cognitive tasks (Haeussinger et al., 2014;Kirilina et al., 2012;Takahashi et al., 2011). Although with variable results, they suggested that taskinduced sympathetic outflow (leading to increased cardiac output and arterial blood pressure), together with skin blood flow/volume changes due to local vasomotion (constriction and/or dilatation), are the mechanisms responsible for these extracranial hemodynamic changes. In our work here, we consider that skin blood flow contribution has been small or negligible due to: (i) the optodeskin interface of our NIRS device reduces the skin blood flow (see Section 2.4), and (ii) our short-channels are long enough (14 mm) to allow light penetrate deeper into the subcutaneous and muscle tissues of the forehead (perhaps even the skull), with little contribution of the skin compared with the volume illuminated. Therefore, albeit we cannot exclude local vasodynamics effects, we point to a systemic drive as the major cause. Our findings revealed that the task induced highly coherent heart/HbO fluctuations at the single-subject level, but with considerable individual variability in their temporal coupling, which leads to inconsistent phase values at the group level ( Figure 9). Many studies have emphasized the link between fNIRS signals, heart rate and blood pressure during rest and under functional stimulation (Franceschini et al., 2000;Franceschini et al., 2006;Kirilina et al., 2012;Minati et al., 2011;Tachtsidis et al., 2008Tachtsidis et al., , 2009Takahashi et al., 2011). In our work, we did not monitor blood pressure, thus preventing the possibility of linking it with heart rate and HbO. However, based on the aforementioned studies, we reasonably assume that the task used here may also have induced blood pressure oscillations contributing to HbO fluctuations.
The observed HbO/HbR surface pattern seems compatible with an oxygenation effect due to an increased arterial inflow in the micro-vascular bed, which in turn leads to parallel HbR changes. Thus, the phase difference between HbO and HbR (−39° to 47°) would reflect the complex contribution of capillary transit time, blood flow and blood volume changes in shallow layers (Elting et al., 2020;Zheng et al., 2010). However, mainly for HbR, other mechanisms as forehead venous volume changes (Kirilina et al., 2012) or even superficial O 2 consumption could overlap. Another interesting possibility is that surface signals also overlap changes in cerebral blood flow arising from increasing metabolic demands and cerebral autoregulation. This is plausible as we placed the NIRS probe over a forehead region mainly supplied by the supraorbital and supratrochlear arteries that ultimately connect (via the ophthalmic artery) to the frontopolar branch of the anterior cerebral artery, which plays a key role in blood supply to the frontal lobes. Prior studies have suggested that blood pressure measured from the supraorbital artery may reflect cerebral perfusion pressure (Lee & Westenskow, 1998;Narus et al., 1995), or that the clamping of the internal carotid affects both supraorbital blood flow and frontal lobe oxygenation (Hove et al., 2006). Jenkins and Brown (2014) also postulated the relationship between frontal activity asymmetry and forehead blood flow in a study using EEG and infrared thermography.
Another interesting finding was that the HbO response precedes the trial onset by some seconds. Previous studies reported that the PFC increases oxygenation a few seconds prior to the onset of voluntary exercise, independently of its actual intensity Ishii et al., 2018;Matsukawa et al., 2015). It has been hypothesized that a feedforward mechanism (termed "central command"), involving higher brain centres, sends descending signals that adjust physiological systems, as the cardiovascular one, to the upcoming effort (Goodwin et al., 1972;Williamson, 2010). Furthermore, there is growing evidence that preparing for a mental challenge induces activity in certain cerebral areas as the Anterior Cingulate Cortex and PFC (Sohn et al., 2007;Vassena et al., 2014Vassena et al., , 2019. We suggest that, here, the anticipatory effect may have been enhanced by the use of such a rhythmic and predictable task. It is tempting to speculate that the aforementioned studies and our results point to the same task-related arousal mechanism that brings fresh arterial blood, full of oxygen to the cortex in preparation for upcoming cognitive demands.
Task-related arousal mechanisms requires a close interaction between cognitive function and autonomic control (Forte, De Pascalis, et al., 2019;Nicolini et al., 2014;Thayer & Lane, 2009;Wang et al., 2016). Thus, the autonomic control appears to be associated with activity levels in executive brain regions, which allows an adaptive response to environmental demands. Conversely, autonomic dysfunctions may be related to the deterioration of certain cognitive functions, specifically of executive functions (Forte, De Pascalis, et al., 2019;. This close coordination of extracerebral and cerebral responses with the task period may have great functional value. The correct coupling between physiological resources may be the sign of proper cognitive and/or cardiovascular function and, its disruption, a potential early marker of cognitive decline and/or cardiovascular disease.