Non-neuronal evoked and spontaneous hemodynamic changes in the anterior temporal region of the human head may lead to misinterpretations of functional near-infrared spectroscopy signals

Abstract. Several functional near-infrared spectroscopy (fNIRS) studies report their findings based on changes of a single chromophore, usually concentration changes of oxygenated hemoglobin ([O2Hb]) or deoxygenated hemoglobin (HHb). However, influence of physiological actions may differ depending on which element is considered and the assumption that the chosen measure correlates with the neural response of interest might not hold. By assessing the correlation between [O2Hb] and [HHb] in task-evoked activity as well as resting-state data, we identified a spatial dependency of non-neuronal hemodynamic changes in the anterior temporal region of the human head. Our findings support the importance of reporting and discussing fNIRS outcomes obtained with both chromophores ([O2Hb] and [HHb]), in particular, for studies concerning the anterior temporal region of the human head. This practice should help to achieve a physiologically correct interpretation of the results when no measurements with short-distance channels are available while employing continuous-wave fNIRS systems.

Non-neuronal evoked and spontaneous hemodynamic changes in the anterior temporal region of the human head may lead to misinterpretations of functional near-infrared spectroscopy signals Guilherne Augusto Zimeo Morais, a, * Felix Scholkmann, b Joana Bisol Balardin, c,d Rogério Akira Furucho, c Renan Costa Vieira de Paula, c Claudinei Eduardo Biazoli Jr., c and João Ricardo Sato c 1 Introduction 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 (O 2 Hb) and deoxygenated (HHb) forms of hemoglobin. 1,2 Similar to functional magnetic resonance imaging (fMRI) and its blood-oxygen-leveldependent (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 O 2 Hb, 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 O 2 Hb 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 O 2 Hb and HHb, two wavelengths are general enough to yield meaningful measures when properly chosen. [6][7][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 O 2 Hb or HHb, even though studies have shown that these measures are individually influenced by the underlying physiology. [9][10][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 O 2 Hb ([O 2 Hb]) 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 [O 2 Hb] 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 Scholkmann 15 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 (PaCO 2 ), 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 nonneuronally 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 [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 available 24 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 [O 2 Hb], [HHb], and total hemoglobin (½tHb ¼ ½O 2 Hb þ ½HHb) concentration changes via the modified Beer-Lambert law. 5 For the measurement, 17 light sources and 16 detectors have been arranged in a 3 × 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 25 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 → substraction to be performed: 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 [mM ¼ mmol∕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%).

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 25 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 toolbox 26 to obtain the channels distribution, based on which we generated the montage representation that is shown in Fig. 1(b).
Resting-state time series were collected from 45 participants (mean age: 24.2 years, standard deviation: 3.8 years, males: 39) during ∼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 CVð%Þ ¼ 100 × standard deviation ðdataÞ∕mean ðdataÞ.
The CV has been used as a metric for signal quality control. All channels with a 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, 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 ([O 2 Hb] and [HHb]) changes based on the modified Beer-Lambert law. 5 28 was calculated based on the general equation proposed by Scholkmann and Wolf 29 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).

Data Analysis
As previously mentioned, both datasets underwent analysis with the goal to assess the correlation between [O 2 Hb] 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.

[O 2 Hb] to [HHb] correlation
To quantify the correlation between [O 2 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 grouplevel 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 Table 1 Defined regions of interest for resting-state data and corresponding channels. Individual channel positions are shown in Fig. 1(a). 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 [−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.

Region of interest Channels
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.

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  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.  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.

Dataset 1: Task-Evoked Activity
The block averages also illustrate the phenomenon observed with the boxplots distribution, i.e., how the task-related responses gradually shift from the expected [O 2 Hb] increase and [HHb] decrease to a concomitant increase of both chromophores over the anterior temporal area. Fig. 2 Task-evoked results for four illustrative channels with color correspondence to the spatial distribution in Fig. 1(a) 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 (R 2 ). Although the mean scatter plot presents a relatively high absolute ρ (0.42), the R 2 is rather low due to the greater distance between fitted line and spread points.
Interestingly, the green points (channels 28 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, Fig. 3 Scatter plots of different metrics of the BVPA (y -axis) with respect to the [O 2 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 (ρ) 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).   Fig. 5(b)].
Particularly for [O 2 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 presented the highest coefficient of determination (R 2 ) for the quadratic curve fit (0.88 and 0.90, respectively).

Discussion and Conclusions
In this study, we observed consistent positive correlation between [O 2 Hb] 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 [O 2 Hb] 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 [O 2 Hb] 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 optimal 8 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 vessels 39 had expected a different spatial localization for the positive [O 2 Hb] 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 [O 2 Hb] 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 demonstrated 19,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 desaturation 45 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 [O 2 Hb] and [HHb] increase in channels positioned over the anterior temporal area. 21 Their report of positive [O 2 Hb] 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 (1.32 swallows∕ min). 48 Additionally, a tonic activation of the temporalis muscle during resting state might be involved in maintaining the mandibular rest position to compensate the   [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 [O 2 Hb] 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 O 2 Hb 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).

Disclosures
GAZM was a full-time employee at NIRx Medizintechnik GmbH during work development and manuscript preparation. All other authors declare no conflicts of interest.