2.1.1.

Dataset 1: task-evoked activity

The first dataset used for the present study is publicly available²⁴ and was originally collected to evaluate the task-evoked responses related to a simple mental arithmetic paradigm. Details on the original experimental design and results obtained can be found in the corresponding report.²²

The dataset has been acquired with a continuous-wave NIRS system (ETG-4000, Hitachi Medical Corporation, Osaka, Japan) that employs two different wavelengths (695 and 830 nm) for intensity data acquisition and calculation of [${\mathrm{O}}_2\mathrm{Hb}$], [HHb], and total hemoglobin ($[\mathrm{tHb}]=[{\mathrm{O}}_2\mathrm{Hb}]+[\mathrm{HHb}]$) concentration changes via the modified Beer–Lambert law.⁵

For the measurement, 17 light sources and 16 detectors have been arranged in a $3\times 11$ grid covering frontal and temporal regions of the cortex, as shown in Fig. 1(a), yielding 52 channels in total. The sampling frequency was set to 10 Hz and the source–detector separation was fixed to 30 mm. The lowest line of sources and detectors was placed along the FP1–FP2 line of the 10 to 20 international system²⁵ with channel 48 exactly at FP1 position.

The experimental design consisted of six blocks (per run) of a mental arithmetic task performed during 12 s followed by a resting period of 28 s. Participants were instructed to mentally subtract a one-digit number from a two-digit number as quickly as possible within the block period. The first subtraction was presented as a visual cue on the screen and the following subtractions should be based on the two-digit result and the initial one-digit number presented (e.g., cue: $97-4\to \text{substraction}$ to be performed: $97-4=93$, $93-4=89$, $89-4=85$, etc.).

The dataset downloaded consists of unfiltered concentration changes obtained from eight subjects (mean age: 26.0 years, standard deviation: 2.8 years, males: 3) separated in different runs. Subjects 01 to 03 contain three runs, whereas the datasets for subjects 05 to 08 have four runs. Data are provided with the scale unit of [mM mm] representing the molar concentration [$\mathrm{mM}=\mathrm{mmol}/\mathrm{L}$] multiplied by the unknown pathlength of each wavelength [mm]. Information on the trigger onsets of each trial of task and rest are also provided. Before proceeding to data analysis based on correlation (as described in Sec. 2.2.1), we applied a bandpass filter of [0.01 to 0.20] Hz (filter type: linear-phase FIR, roll-off width: 15%).

2.1.2.

Dataset 2: resting-state activity

The second dataset was collected at the Universidade Federal do ABC following approval by the local Ethics Committee and upon written informed consent provided by all participants.

A continuous-wave NIRS system (NIRScout, NIRx Medical Technologies, Glen Head, New York) that uses two different wavelengths (760 and 850 nm) was employed. We used 16 light sources and 16 detectors to cover most of the left hemisphere of the participants, as shown in Fig. 2(a). This layout resulted in 49 channels measured with a sampling frequency of 3.91 Hz. Sources and detectors positions were based on the 10-5 international system²⁵ and the source–detector separations ranged from 20 (channels 17 and 24) to 30 mm (all other channels). The upper separation limit was achieved with the aid of stabilizing links provided by the manufacturer and placed for all channels between their forming source and detector positions.

Fig. 2

Task-evoked results for four illustrative channels with color correspondence to the spatial distribution in Fig. 1(a). Boxplots with correlation results over all subjects are shown on the top. Block average results (over runs) of a representative subject (#08) are depicted on bottom. [${\mathrm{O}}_2\mathrm{Hb}$] time series are plotted as red and [HHb] as blue signals. The shaded area represents the standard deviation.

Ten different regions of interest were visually defined in respect to underlying anatomical regions (Table 1). To retrieve the channel positions, we considered the MNI coordinates of the 10-5 international system positions as provided by Jurcak et al.²⁵ We imported them into the SPM for fNIRS toolbox²⁶ to obtain the channels distribution, based on which we generated the montage representation that is shown in Fig. 1(b).

Table 1

Defined regions of interest for resting-state data and corresponding channels. Individual channel positions are shown in Fig. 1(a).

Resting-state time series were collected from 45 participants (mean age: 24.2 years, standard deviation: 3.8 years, males: 39) during $\sim 5\min$. For consistency purposes, data from all participants were truncated during offline processing from the first data point to the 1172nd, therefore, resulting in exactly 5 min long time series. Prior to data acquisition, subjects were instructed to keep their eyes opened and not to perform structured thought processes during the measurement.

The coefficients of variation (CVs) of the truncated intensity data were then calculated for each channel and from each subject according to $\mathrm{CV}(\%)=100\times \text{standard deviation}(\text{data})/\text{mean}(\text{data})$.

The CV has been used as a metric for signal quality control. All channels with a $\mathrm{CV}>7.5\%$ have been excluded from the data analysis because they contained unphysiological noise. Data collected from 38 to 45 subjects were considered per each channel. In the supplementary material,²⁷ we have included a table with the number of subjects considered per channel of interest.

After performing a signal quality check, intensity data have been bandpass filtered within the range (0.01 to 0.20) Hz (filter type: linear-phase FIR, roll-off width: 15%). The lower cutoff frequency was set to remove very low frequency oscillations and the higher to eliminate both respiratory and cardiac components of the signal.

The filtered data were then converted to optical density data and then to ([${\mathrm{O}}_2\mathrm{Hb}$] and [HHb]) changes based on the modified Beer–Lambert law.⁵ The extinction coefficients used for 760 nm were 1.4866 ([${\mathrm{O}}_2\mathrm{Hb}$]) and 3.8437 ([HHb]), whereas for 850 nm they were 2.5264 ([${\mathrm{O}}_2\mathrm{Hb}$]) and 1.7986 ([HHb]) in units of $[(1/\mathrm{cm})/(\mathrm{mmol}/\mathrm{L})]$. The different pathlength factor (DPF)²⁸ was calculated based on the general equation proposed by Scholkmann and Wolf²⁹ considering a mean age of 24 years for all subjects and yielding 6.12 for 760 nm and 5.06 for 850 nm. Finally, the baseline was defined as the mean of the whole time series (5 min).

2.2.1.

[O₂Hb] to [HHb] correlation

To quantify the correlation between [${\mathrm{O}}_2\mathrm{Hb}$] and [HHb], the Spearman’s correlation coefficient was chosen because this parameter does not expect that the signals have a linear relationship with each other, and it is also more robust against outliers (e.g., motion artifacts within the window of interest).

For the task-related data, we calculated the Spearman’s correlation for each mental arithmetic trial within a 30-s window following the onset time to accommodate the complete hemodynamic response until the signal returned to its baseline. We computed the median of the correlation results over trials to obtain the correlation per run. Similarly, the correlation for a given subject was obtained with the median over runs. This resulted in a matrix with eight rows (subjects) and 52 columns (channels).

For the resting-state dataset, as there is no defined task-window as a reference to compute the correlation, we rather calculated it based on a 20 s sliding window. The length of the window has been chosen arbitrarily. The correlation for a given subject was based on the median over all windows within the time series. Similarly, this resulted in a matrix with 45 rows (subjects) and 49 columns (channels). However, in this case, entries corresponding to the excluded channels were replaced with not-a-number (NaN) values to be ignored in the group-level correlation and distribution results. The procedure of replacing results obtained from excluded channels with missing data was done to avoid bias of nonsatisfactory signal quality in the results and achieved using MATLAB functions “nanmedian” and “boxplot,” which ignore NaN values.

The matrix of correlation results for both task-related and resting-state datasets was then used to generate boxplots for each channel to better illustrate the distribution over subjects. In addition, we computed the median correlation over subjects.

2.2.2.

Blood volume pulse amplitude

Concerning the hemodynamics related to the task-related dataset, we hypothesized that a further metric to evaluate the influence of systemic changes in each channel would be the BVPA. The blood volume pulsation due to the periodic heart activity is generally present in fNIRS signals.^30,31

First, we considered the unfiltered concentration changes and applied a different bandpass filter of [0.50 to 2.50] Hz (filter type: linear-phase FIR and roll-off width: 15%) to extract the blood volume pulsation from the signals. Then, we computed the BVPA as the difference between the upper and lower envelopes of the heartbeat oscillation from [tHb]. Upper and lower envelopes have been computed based on the peak-detection algorithm described by Scholkmann et al.³²

As our hypothesis is based on the BVPA response as evoked by the task, we computed block averages of the BVPA responses within the window [$-5$, 30] s, in which 0 s corresponds to the mental arithmetic trial onset. The window upper limit corresponds to that used for the correlation computation (see Sec. 2.2.1) and the lower limit is to allow for 5 s prior to the onset to be counted as a baseline. The baseline mean was subtracted from each data point within the following 30 s of response block to properly normalize the blocks’ prior averaging.

Similar as for the correlation calculation, we computed the mean of block averages results over runs and then over subjects. We preferred to use the “mean” function (instead of “median”) to preserve the continuity of the signal within the response window for the grand average.

Finally, we computed further metrics of interest based on the grand average results, now within the window [0, 30] s, as follows: (a) peak-to-peak, (b) mean, (c) variance, and (d) median absolute deviation. These were calculated to have a single representative value for the systemic influences in each channel to be related to the correlation results obtained for the same channel, as shown in Sec. 3.1.2.

2.2.3.

Channel-wise correlation

For the resting-state dataset, we computed the channel-wise correlation as a metric to assess the global influence of expected spontaneous vascular and metabolic fluctuations that are present in both task and resting states.^33,34 As the primary analysis goal was to evaluate the presence of global fluctuations, we did not apply any further correction (e.g., principal component analysis^35,36).

Similar to [${\mathrm{O}}_2\mathrm{Hb}$] to [HHb] correlation analysis (Sec. 2.2.1), we used a sliding window of 20-s length here to calculate the Spearman’s correlation between all channels for each chromosphere, ([${\mathrm{O}}_2\mathrm{Hb}$], [HHb], and [tHb]), for each subject. After that, we computed the median correlation value over all subjects. This procedure resulted in a correlation matrix of 49 rows and 49 columns (49 is the number of measured channels) for each chromophore. An $(X,Y)$ entry in the matrix, therefore, represents the median correlation between channel $X$ and channel $Y$ for a given chromophore over all subjects.

The global influence of spontaneous signals’ fluctuations was assessed with the median correlation value over all channels, yielding a vector of length 49 for each chromophore. An entry is a single representative value on the correlation between a given channel and all other measured channels.

3.1.1.

At the anterior temporal region, the [O₂Hb] response is positively correlated with the [HHb] response

According to the correlation analysis described in Sec. 2.2.1, we obtained results that demonstrate a consistent positive correlation between [${\mathrm{O}}_2\mathrm{Hb}$] and [HHb] on channels over the anterior temporal region in both hemispheres. The results for all channels are available in the supplementary material.²⁷

Figure 2 shows the correlation distribution of four particular channels (28, 29, 40, and 51) that are highlighted in Fig. 1 as green, yellow, and red circles. Boxplots and their central marks gradually shift toward positive correlation values while the measurement channels get closer to the anterior temporal area. In particular, channel 51 (red) presents not only the median value, but also shows the first quartile in the positive area.

The block averages also illustrate the phenomenon observed with the boxplots distribution, i.e., how the task-related responses gradually shift from the expected [${\mathrm{O}}_2\mathrm{Hb}$] increase and [HHb] decrease to a concomitant increase of both chromophores over the anterior temporal area.

3.1.2.

At the anterior temporal region, the [O₂Hb] to [HHb] correlation strength is correlated with the blood volume pulse amplitude

As outlined in Sec. 2.2.2, we have hypothesized the task-evoked response of the BVPA to be an indicator related to the influence of systemic changes in that channel. Therefore, we calculated the grand average results of BVPA for each channel available and considered four different metrics to explore their relation to the previously obtained [${\mathrm{O}}_2\mathrm{Hb}$] to [HHb] correlation over all subjects.

From all four metrics explored in Fig. 3, the quadratic fit for peak-to-peak, variance, and median absolute deviation has resulted in the greatest correlation ($\rho$) values as well as adjusted coefficient of determination ($R^2$). Although the mean scatter plot presents a relatively high absolute $\rho$ (0.42), the $R^2$ is rather low due to the greater distance between fitted line and spread points.

Fig. 3

Scatter plots of different metrics of the BVPA ($y$-axis) with respect to the [${\mathrm{O}}_2\mathrm{Hb}$] to [HHb] correlation ($x$-axis) for task-evoked data. From left to right and then top to bottom, the BVPA metrics are: (a) peak-to-peak, (b) mean, (c) variance, and (d) median absolute deviation (mad). For each metric, a polynomial fit has been applied to the data (red), whereas the orange curves represent the 95% confidence interval. A linear fit was empirically found to better fit the mean, whereas the quadratic fit for all the other metrics. Each graph presents the correlation ($\rho$) between the variables as well as the adjusted coefficient of determination ($R^2$). Colored points (green, yellow, and red) relate to color-coded channels in Fig. 1(a). Channels highlighted in Fig. 2 are marked with a white number (28, 29, 40, and 51).

Interestingly, the green points (channels 28 and 29) both presented low [${\mathrm{O}}_2\mathrm{Hb}$] to [HHb] correlation and low BVPA results [metrics (a), (c), and (d)], whereas the yellow (channel 40) lies on a transition portion and the red (channel 51) presents relatively high values for both correlation and BVPA results [(a), (c), and (d)].

3.2.1.

At the anterior temporal region, the [O₂Hb] oscillations are positively correlated with the [HHb] ones

As described in Sec. (2.2.1), for the resting-state dataset, we have obtained a correlation matrix with 45 rows (subjects) and 49 columns (channels). To avoid plotting all 49 channels individually, and as these results can be found in the supplementary material,²⁷ we preferred to group them in 10 different regions of interest, which consisted of groups from four to six channels (Table 1). Therefore, before generating the boxplots results with the distribution over all subjects, we first calculated the mean results within each region of interest.

From the results presented in Fig. 4, the most important remark is related to the central mark of each boxplot (median value) for each region of interest. While the median result for the anterior temporal region represents a positive correlation, the central mark of all other regions of interest is around $-0.5$ instead, whereas their third quartile is considerably below the null correlation, too.

Fig. 4

Boxplots illustrating the [${\mathrm{O}}_2\mathrm{Hb}$] to [HHb] resting-state correlation results over all subjects for each region of interest according to the underlying anatomical regions. The red box highlights the anterior temporal region (channels 14, 21, 28, 29, and 30).

3.2.2.

At the anterior temporal region, the [O₂Hb] to [HHb] correlation strength depends on the channel-wise correlation strength of [O₂Hb], [HHb], and [tHb]

According to the methods described in Sec. 2.2.3, for the resting-state data, we have also assessed the influence of global components on each channel and explored its relation to the previously obtained [${\mathrm{O}}_2\mathrm{Hb}$] to [HHb] correlation.

We calculated the channel-wise median correlation for each chromophore over all subjects [Fig. 5(a)] and then computed the median value over channels as a metric of global influence. The latter has been plotted against the [${\mathrm{O}}_2\mathrm{Hb}$] to [HHb] correlation results obtained in Sec. 3.2.1 [Fig. 5(b)].

Fig. 5

(a) The median channel-wise correlation over all subjects (resting-state data) is depicted with scaled colors for different chromophores. Orange arrows have been added on top of each image to indicate the channels related to the anterior temporal region. (b) The scatter plots of the median channel-wise correlation for each chromophore in respect to the [${\mathrm{O}}_2\mathrm{Hb}$] to [HHb] correlation. Red curves indicate the quadratic fit and the orange lines the 95% confidence interval. On the left superior corner, corresponding $\rho$ and $R^2$ values for each plot are shown.

Particularly for [${\mathrm{O}}_2\mathrm{Hb}$] and [tHb], the channels on the anterior temporal region presented a relatively low channel-wise correlation in comparison to the results obtained for all other channels. This can be evidenced from the contrast of colors in the depicted images, which can represent a lower influence of global spontaneous fluctuations on these channels’ signals.

The latter observation is confirmed with the scatter plots showing that channels with stronger [${\mathrm{O}}_2\mathrm{Hb}$] to [HHb] correlation also present a lower influence of global components. While global components seem to be most prominent for [${\mathrm{O}}_2\mathrm{Hb}$] and [tHb], all chromophores presented a strong negative correlation (lower than $-0.6$) between the global component metric ($y$-axis) and the [${\mathrm{O}}_2\mathrm{Hb}$] to [HHb] correlation. Also, [${\mathrm{O}}_2\mathrm{Hb}$] and [tHb] presented the highest coefficient of determination ($R^2$) for the quadratic curve fit (0.88 and 0.90, respectively).