Functional near-infrared spectroscopy (fNIRS) is an optical neuroimaging technique enabling to quantify changes in hemodynamics and oxygenation within the intracerebral and extracerebral tissue layers of the human head based on the absorption of near-infrared light by the oxygenated () and deoxygenated (HHb) forms of hemoglobin.1,2 Similar to functional magnetic resonance imaging (fMRI) and its blood-oxygen-level-dependent (BOLD) signal, the inference on neural activation relies on the principle of neurovascular coupling and thus a close relation with changes in cerebral blood flow.3 However, while traditional fMRI inference is based solely on the BOLD signal, fNIRS has expanded the sensitivity to intrinsic chromophores, including , HHb, and the oxidized form of the cytochrome c-oxidase (oxCCO),4 in addition to further dependent quantities, including total hemoglobin, tissue oxygen saturation, and exchange efficiency.
Based on the modified Beer–Lambert law,5 the number of different intrinsic chromophores that can be measured by fNIRS systems is mathematically limited by the number of wavelengths employed for the measurement. Different studies have shown that for a proper measure of and HHb, two wavelengths are general enough to yield meaningful measures when properly chosen.67.–8 Although the great majority of commercially available fNIRS systems employs at least two wavelengths for the measurement,1 many reports are still based on a single chromophore, usually or HHb, even though studies have shown that these measures are individually influenced by the underlying physiology.910.–11 Activation inference based on neurovascular coupling foresees an increase of oxygen metabolism followed by an outweighed increase of cerebral regional blood flow, thus resulting in an increase in the concentration of () expected to be accompanied by a decrease in concentration of HHb ([HHb]).12 Thus, the interpretation of the results should be based on both chromophores.
Report findings based on a single chromophore assume that the chosen measure correlates with the neural response of interest. However, the observation of an increase of  alone may not be enough to make proper inferences about the underlying neural activity, as an accompanying increase of [HHb] could be an indicator of hemodynamic changes not evoked by neural activity, happening in the intracerebral and/or extracerebral tissues.13,14 In fact, Tachtsidis and Scholkmann15 recently summarized different physiological sources that may lead to false positives and false negatives in fNIRS measures and the reality is that the fNIRS signal assumed to always be neuronally evoked may be present in many other physiological components. These include changes in partial pressure of carbon dioxide in the arterial blood (),16,17 blood pressure,18,19 and activity of the sympathetic nervous system.19,20 As shown in these studies, the physiological confounds may lead to non-neuronally evoked responses, therefore, corresponding result interpretations should be done carefully because the hemodynamic responses measured are not due to neurovascular coupling.
In the present study, we report observations that identify a spatial dependency of non-neuronal hemodynamic changes to be of significance for fNIRS studies: hemodynamics on the anterior temporal region, which might be caused by “head muscles activity,” as also recently reported by Volkening et al.21 Based on the expected negative correlation between  and [HHb] for task-related responses, we observed that the anterior temporal region may be especially susceptible to this behavior. In addition, we demonstrate that the negative  to [HHb] correlation is sustained for resting-state spontaneous fluctuations for all cortical regions on the left hemisphere but not for the anterior temporal area.
Materials and Methods
In this study, we have assessed the correlation between  and [HHb] in two different datasets of experiments involving the measurement of (i) task-related and (ii) resting-state brain activities. The materials and methods corresponding to each dataset are described in detail in the following sections. The spatial distribution of channels used to obtain each set of data is shown in Fig. 1.
Dataset 1: task-evoked activity
The first dataset used for the present study is publicly available24 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.22
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 , [HHb], and total hemoglobin () concentration changes via the modified Beer–Lambert law.5
For the measurement, 17 light sources and 16 detectors have been arranged in a 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 system25 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: to be performed: , , , 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  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%).
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 system25 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.
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.25 We imported them into the SPM for fNIRS toolbox26 to obtain the channels distribution, based on which we generated the montage representation that is shown in Fig. 1(b).
Defined regions of interest for resting-state data and corresponding channels. Individual channel positions are shown in Fig. 1(a).
|Region of interest||Channels|
|Frontal anterior||1, 2, 3, 4, 6, 8|
|Frontal superior||5, 9, 10, 15, 16|
|Frontal inferior||7, 11, 12, 13, 19, 20|
|Temporal anterior||14, 21, 28, 29, 30|
|Temporal posterior||36, 42, 43, 44|
|Premotor||17, 18, 22, 23|
|Motor||24, 25, 26, 31|
|Parietal inferior||34, 35, 37, 38, 40|
|Parietal superior||27, 32, 33, 39, 41|
|Occipital||45, 46, 47, 48, 49|
Resting-state time series were collected from 45 participants (mean age: 24.2 years, standard deviation: 3.8 years, males: 39) during . 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 .
The CV has been used as a metric for signal quality control. All channels with a 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,27 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 ( and [HHb]) changes based on the modified Beer–Lambert law.5 The extinction coefficients used for 760 nm were 1.4866 () and 3.8437 ([HHb]), whereas for 850 nm they were 2.5264 () and 1.7986 ([HHb]) in units of . The different pathlength factor (DPF)28 was calculated based on the general equation proposed by Scholkmann and Wolf29 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).
As previously mentioned, both datasets underwent analysis with the goal to assess the correlation between  and [HHb]. In addition, we employed further metrics to better quantify and understand the relation of these results with other behaviors of interest. For the dataset with task-evoked data, we also evaluated the grand average of the blood volume pulse amplitude (BVPA); and for the resting-state, channel-wise correlation for each chromophore was calculated for all channels as seeds.
[O2Hb] to [HHb] correlation
To quantify the correlation between  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.
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.32
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 [, 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.
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 analysis35,36).
Similar to  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, (, [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 entry in the matrix, therefore, represents the median correlation between channel and channel 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.
Dataset 1: Task-Evoked Activity
At the anterior temporal region, the [O2Hb] 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  and [HHb] on channels over the anterior temporal region in both hemispheres. The results for all channels are available in the supplementary material.27
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  increase and [HHb] decrease to a concomitant increase of both chromophores over the anterior temporal area.
At the anterior temporal region, the [O2Hb] 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  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 () values as well as adjusted coefficient of determination (). Although the mean scatter plot presents a relatively high absolute (0.42), the is rather low due to the greater distance between fitted line and spread points.
Interestingly, the green points (channels 28 and 29) both presented low  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)].
Dataset 2: Resting-State Activity
At the anterior temporal region, the [O2Hb] 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,27 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 instead, whereas their third quartile is considerably below the null correlation, too.
At the anterior temporal region, the [O2Hb] to [HHb] correlation strength depends on the channel-wise correlation strength of [O2Hb], [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  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  to [HHb] correlation results obtained in Sec. 3.2.1 [Fig. 5(b)].
Particularly for  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  to [HHb] correlation also present a lower influence of global components. While global components seem to be most prominent for  and [tHb], all chromophores presented a strong negative correlation (lower than ) between the global component metric (-axis) and the  to [HHb] correlation. Also,  and [tHb] presented the highest coefficient of determination () for the quadratic curve fit (0.88 and 0.90, respectively).
Discussion and Conclusions
In this study, we observed consistent positive correlation between  and [HHb] both for task-evoked responses and resting-state data in channels placed over the anterior temporal area of the human head. The positive  to [HHb] correlation phenomenon has been previously suggested to occur due to (i) extracerebral and/or systemic signals,13 (ii) motion artifacts and bad coupling,37 (iii) cross talk between chromophores,8,38 (iv) predominance of capillaries contribution to the measured hemodynamic signal,39 or (v) hemodynamic and oxygenation in muscles.21 In the following paragraphs, we will critically discuss the possible contributions of each one of these factors for the interpretation of the findings we observed.
Based on the methods we applied to analyze the data, we believe that our results are likely not due to motion artifacts. The Spearman’s correlation is more robust in the presence of outliers (spike artifacts) in the data and, in addition to that, we have calculated the median correlation over blocks to also increase the robustness to outliers. Also, our results suggest that bad coupling could not explain our findings either, as the illustrative block average signals for the task-related data appear rather clean and there also appears to be a smooth transition of responses patterns while approaching the anterior temporal area. For the resting-state data, we have controlled the coefficient of variation to exclude any channels with potential bad coupling prior to the analysis.
Cross talk between  and [HHb] chromophores has been suggested to occur either when the wavelengths pairs that are used to measure the intensity data and estimate the concentration changes are not optimal8 or when the DPF values used for the modified Beer–Lambert law are not precise.38 This concern, however, does not apply as the pairs that are employed for NIRx NIRScout (760 and 850 nm) and Hitachi ETG-4000 (695 and 830 nm) are pairs that yield optimal measures.8,40 Also, while for the task-related dataset the concentration changes were based on unknown pathlength factor, for the resting-state data we used a general equation 29 to estimate the optimal DPF based on both the wavelength as well as the mean subjects’ age from this study, therefore, we do not expect this to yield cross talk either. Although a dependence of the DPF on anatomical variability is expected,41 the large number of subjects examined makes it unlikely this would be the cause of a significant bias.
Finally, the study suggesting a theoretical model to explain the increase of [HHb] as a task-evoked response in regions with a greater density of capillaries in comparison to large vessels39 had expected a different spatial localization for the positive  to [HHb] correlation, namely the Broca area. We believe this spatial localization difference may be due to individual anatomical differences, as their study only measured a single subject and considering that the Broca area (inferior frontal gyrus) is very close to the anterior temporal region. However, as their theory does not expect the positive correlation to be present over the temporal region, in our opinion, it could not explain our findings either.
In face of the considerations presented and discussed above, we believe that the spatially localized positive  to [HHb] correlation observed in our study might be better explained by extracerebral changes and/or muscle oxygenation. Figure 6 shows potential anatomical peculiarities concerning the anterior temporal region: a major presence of extracranial vessels and head muscles.
The impact of the extracranial vessels on the measured hemodynamic signals has been previously demonstrated19,43 and this could support the origins of the observed response in our study based on the higher density of superficial vessels as shown in Fig. 6(a). Nevertheless, we consider that the temporalis muscle might yield the most prominent influence on the signals measured over the anterior temporal region because of an intrinsic muscle activity and corresponding hemodynamic change, as it has been demonstrated, e.g., measuring the anterior tibial artery flow with fMRI,44 and the vastus lateralis muscle desaturation45 as well as resting and exercising skeletal muscles metabolism with NIRS.46
It was also recently shown that teeth clenching, which is expected to trigger the activation of temporalis muscle and related hemodynamic changes, yields a pattern of concurrent  and [HHb] increase in channels positioned over the anterior temporal area.21 Their report of positive  to [HHb] correlation following particular head muscles activity might be further evidence of the potential source of our task-related findings over the anterior temporal region. Also, in resting-state data, the activation of the temporalis muscle can occur during spontaneous saliva swallowing,47 although its frequency is relatively low ().48 Additionally, a tonic activation of the temporalis muscle during resting state might be involved in maintaining the mandibular rest position to compensate the gravitational forces pulling the mandible down, as observed by Yilmaz et al.49
In addition to the positive  to [HHb] correlation, we explored other metrics in the present study that provide supporting evidence to particularities in the anterior temporal region. With the task-evoked responses, we have observed a close correlation of the observed positive  to [HHb] correlation with the BVPA. This suggests that the simultaneous increase of both  and [HHb] in that area presents non-neuronal origins, as it has been shown that the scalp blood flow exhibits a strong blood volume pulsation corresponding to the heart rate.50 Also, we observed that the spontaneous resting-state global fluctuations of both  and [tHb] are not as prominent in the anterior temporal region channels as in all other cortical areas over the left hemisphere, which further supports the idea that there should be a localized source of extracortical signals impacting the measures over the particular anterior temporal area.
Nevertheless, our study presents a few limitations that should be addressed. The first is related to the mean population age assessed with both datasets: 26.0 years for the task-related and 24.2 years for the resting-state data, both with low standard deviation. Therefore, generalization of present findings to other populations, e.g., infants or elderly people, should be done carefully, as an age dependency of the skeletal muscles’ BOLD signal was previously reported.51 The second limitation concerns the absence of direct measures of muscle activity (e.g., concurrent electromyography) to be correlated with the hemodynamic changes measured in both mental arithmetic task and resting state. The last limitation is related to the absence of short-distance measurement channels to separate the extracortical information from the cortical signals to assess our theory for the observed findings concerning the temporal extracortical layers.
We conclude that the examination of individual components of hemoglobin signal alone can lead to misinterpretation of fNIRS findings because influence of physiological actions may differ depending on which element is considered. The report of results based on both and HHb is, therefore, important for proper interpretation of anterior temporal-related outcomes. Inclusion of short-distance measurement channels into this region to separate neuronal signals from extracortical components may further improve results interpretation. These guidelines seem to be important not only for studies focused on task-evoked responses, but also for resting-state studies (e.g., connectivity analysis).
GAZM was a full-time employee at NIRx Medizintechnik GmbH during work development and manuscript preparation. All other authors declare no conflicts of interest.
The authors express their appreciation to Dr. Randall Barbour for his review and many helpful comments. The authors would also like to thank Dr. Günther Bauernfeind for his prompt and detailed clarifications on the task-related data. The authors are also grateful for the support provided by Mr. Jackson Cionek (Brain Support Corporation).