1.

Introduction

Functional near-infrared spectroscopy (fNIRS) measures relative changes by radiating near-infrared lights of two or more wavelengths and records the reflected lights. Emitted infrared lights are differentially absorbed by oxyhemoglobin and deoxyhemoglobin and relative concentrations of blood oxyhemoglobin and deoxyhemoglobin are measured,¹ providing a noninvasive technique for measuring blood oxygen level dependent (BOLD) responses to neuronal activation in the cerebral cortex.^2,3 Since the relative concentrations of oxyhemoglobin and deoxyhemoglobin change in opposite directions, their signals would be expected to have similar waveforms with opposite polarity, and the amplitudes of oxyhemoglobin and deoxyhemoglobin would have similar spatial patterns provided that the fNIRS signal originates only from brain regions associated with the task. However, fNIRS data are often not consistent with this expectation and raise concerns for interpretation.

Systemic effects, such as blood pressure, respiratory and blood flow variation, alter relative blood oxyhemoglobin and deoxyhemoglobin concentrations. It is well established that fNIRS signals contain both global (systemic) and cortical BOLD effects.^2–5 In this paper, we use the term global components to denote the combined systemic components not associated with task-specific activity in fNIRS data. The global component may be more pronounced in fNIRS than in functional magnetic resonance imaging (fMRI) acquisitions because the near-infrared lights pass through superficial layers of blood vessels in skin and yield significant signal changes primarily due to systemic components. The removal of the global component continues to be the central focus of fNIRS research.^4,6,7–12 Several approaches have been previously proposed.

Most approaches involve using either the waveforms or the temporal frequency characteristic of the global components. The waveform of the global components can be measured with laser-Doppler flowmetry,¹³ continuous blood pressure,⁵ or the fNIRS signal from skin using a pair of optodes $<1\mathrm{cm}$ away.^13,14 A general linear regression algorithm has been used to remove these global components from fNIRS data,^5,15 and conventional signal processing methods, such as independent component analysis¹⁴ and temporal filtering methods,¹⁰ have also been applied. These methods assume that the waveform of the task-related neuronal signal is not correlated with the waveform of the global components. However, in some cases, the waveform of the global component is found to be highly correlated with task-related fNIRS signals.^16–18 For example, blood pressure and blood flow may show a waveform similar to the expected cortical BOLD signal responding to a task with a block design.^5,13 To overcome the challenge that the global component can be highly correlated with the task-induced BOLD signal, an alternative method using the spatial distribution for the removal of global component^3,11 has been proposed. Since the only common components between resting state data and the experimental data are the global components, such as respiratory and the blood pressure variation, which have been shown to be highly correlated with the first and second principal components (PCs) of the resting state signal, in theory, we can obtained the spatial pattern of the first and second PCs from the resting state data and remove them from the experimental dataset using a linear regression method.^3,11 The mathematical concept of this idea is elegant. However, the global components are relatively weak during the resting state. Therefore, to obtain high-quality global components for subsequent data analysis, prolonged resting state data may be needed. In addition, more than three PCs are often needed to fully describe global components. However, according to the authors, removing more PCs may cause reduction of neural signals related to the tasks.

Here, we present an alternative spatial filtering approach for isolating the global components in fNIRS signals to address some of these limitations. Neurologically originated fNIRS signals have a much narrower spatial distribution than that of global components, based on the organizing principal that specific functions engage specific neuronal patterns of activity. Therefore, the global systematic components can be separated from the cortical BOLD signal by Gaussian spatial filtering. There are three principle differences between our approach and the previous spatial filtering approach. Resting state data are not required; the spatial locations of each channel, in terms of Montreal Neurological Institute (MNI) coordinates, are utilized; and the global components are defined as both temporal waveforms and spatial patterns and are extracted from the entire dataset. The technique, however, requires optode coverage over an area much larger than the expected pattern of activity. We use a smoothing kernel 10 times the size of the standard spatial smoothing kernel for fMRI (8 mm). Such a large-sized smoothing kernel not only ensures the algorithm to be robust to random noise but also ensures that the smoothed spatial pattern contains little task-related neuronal components and thus enables the cortical BOLD signal to be derived from acquired data. We test the approach using fNIRS responses to a finger-thumb tapping task, which yields a well-established activity pattern, and by comparison of the oxyhemoglobin and deoxyhemoglobin signals.

2.

Methods

2.4.

Localizing the Optodes and the Definition of Region of Interest

The locations of both emitters and receivers, along with land marks such as inion, nasion, Cz and left and right ears, were determined using a three-dimensional (3-D) digitizer (Polhemus Tech, Vermont). The MNI coordinates for the recording channels were obtained using the NIRS_statistical parametric mapping (SPM) program,²¹ and the corresponding anatomical locations of each channel were determined with the atlas provided in NIRS_SPM. Figure 1(a) shows normalized channel locations superimposed on a brain surface. The region of interest (ROI), left hemisphere motor area, was established a priori based on all channels with $>25\%$ chance of being part of the left motor or premotor region, determined separately for each individual subject (enclosed areas in Fig. 1). Table 1 lists the anatomical information for channels within the ROI for two typical subjects. Figure 1(b) shows the schematic channel layout used for subsequent event-triggered average plots.

Fig. 1

(a) Locations and number identifications of recording channels superimposed on standard MNI brain surface for two typical subjects, shown as red dots. The dark red dots indicate the channels identified as primary or premotor cortex on the left hemisphere. (b) The schematic channel layout used for the event triggered averages (see Fig. 2). The circled ROIs include the anatomical regions (channels shown as dark red dots above, see Table 1). Note the individual difference in anatomical locations for their channels.

Table 1

The channel list for left motor areas for two typical subjects (see Fig. 1). The channel identification number is indicated on the left column and the X, Y, and Z columns indicate the MNI coordinates of the channel. “Prob.” indicates the probability that a channel is associated with the anatomical region, listed on the right column. Brodmann area is indicated by column “BA.”

S1,Ch#	X	Y	Z	Prob.	BA	Anatomy
37	$-64$	7	20	0.61	6	Premotor and supplementary motor cortex
				0.27	43	Subcentral area
38	$-56$	7	42	0.74	6	Premotor and supplementary motor cortex
39	$-42$	11	59	0.54	6	Premotor and supplementary motor cortex
				0.34	9	Dorsolateral prefrontal cortex
40	$-22$	10	70	0.81	6	Premotor and supplementary motor cortex
				0.19	8	Includes frontal eye fields
49	$-52$	$-3$	55	0.92	6	Premotor and supplementary motor cortex
50	$-32$	$-2$	67	1	6	Premotor and supplementary motor cortex
51	$-13$	$-4$	76	1	6	Premotor and supplementary motor cortex
60	$-46$	$-16$	64	0.62	4	Primary motor cortex
				0.27	3	Primary somatosensory cortex
61	$-24$	$-16$	76	0.86	6	Premotor and supplementary motor cortex
71	$-39$	$-28$	70	0.64	4	Primary motor cortex
				0.3	3	Primary somatosensory cortex
72	$-13$	$-31$	80	0.87	4	Primary motor cortex
S2,Ch#	X	Y	Z	Prob.	BA	Anatomy
30	$-12$	18	69	0.66	6	Premotor and supplementary motor cortex
				0.34	8	Includes frontal eye fields
37	$-66$	$-2$	17	0.54	43	Subcentral area
				0.16	48	Retrosubicular area
				0.15	22	Superior temporal gyrus
				0.14	6	Premotor and supplementary motor cortex
38	$-60$	1	39	0.61	6	Premotor and supplementary motor cortex
				0.22	4	Primary motor cortex
39	$-44$	5	59	0.94	6	Premotor and supplementary motor cortex
40	$-22$	4	72	0.98	6	Premotor and supplementary motor cortex
48	$-67$	$-13$	30	0.55	43	Subcentral area
				0.26	2	Primary somatosensory cortex
49	$-56$	$-12$	52	0.4	4	Primary motor cortex
				0.31	6	Premotor and supplementary motor cortex
				0.28	3	Primary somatosensory cortex
50	$-36$	$-11$	70	0.8	6	Premotor and supplementary motor cortex
51	$-12$	$-8$	78	1	6	Premotor and supplementary motor cortex
61	$-26$	$-23$	75	0.55	4	Primary motor cortex
				0.45	6	Premotor and supplementary motor cortex

3.

Results

Our task is expected to elicit activity in left motor cortex areas, as well as related premotor and somatosensory areas. The expected locations of motor cortex areas are circled based on the predetermined anatomical locations of each channel for two typical subjects in Fig. 2 showing event triggered average data. The difference between the acquired signals [Figs. 2(a) and 2(d)] and the extracted global components [Figs. 2(b) and 2(e)] represent the derived neuronal signals [Figs. 2(c) and 2(f)].

Fig. 2

The event-triggered average fNIRS data responding to the finger thumb tapping task are shown for two typical subjects and for both oxyhemoglobin (red) and deoxyhemoglobin (blue) signals. The time span for each trace is 30 s, where the task period is from 0 to 15 s and the rest period is from 15 to 30 s. (a) and (d) The acquired fNIRS signal [H in Eq. (5)] is shown. (b) and (e) The global component [${\mathrm{H}}_{\text{Global}}$ in Eq. (5)] is shown. (c) and (f) The derived neuronal component [${\mathrm{H}}_{\text{Neuronal}}$ in Eq. (5)] is shown.

To estimate the effect of the global component across all subjects, we first performed SPM GLM analysis using the canonical BOLD model. The results of GLM estimation, or the beta values, for each individual subject were projected onto the standard MNI brain and, the group analysis was performed using SPM 12. Figure 3 shows the results of SPM group analysis for both the acquired signals [Figs. 3(a) and 3(c)] and the derived neuronal component [Figs. 3(b) and 3(d)] for both oxyhemoglobin and deoxyhemoglobin signals ($p<0.005$). The acquired oxyhemoglobin fNIRS data [Fig. 3(a)] display a wide-spread, nonspecific positive activation pattern, consistent with individual event-triggered averages [red traces in Figs. 2(a) and 2(e)]. Interestingly, the global component of the deoxyhemoglobin data appears to be nearly absent [Fig. 3(c)]. By contrast, the derived neuronal components show expected localized neural activity in the group analysis data [Figs. 3(b) and 3(d)]. Information about the active clusters is listed in Table 2. The derived neuronal oxyhemoglobin signal shows significant negative activity in clusters identified in Table 2. Although previous fMRI results^25,26 have suggested negative BOLD signals can occur in the contralateral hemisphere for a finger-thumb tapping task, we cannot interpret the nature of these negative signals. A highpass spatial filter may result in artificial negative activity or exaggerate true negative activity.

Fig. 3

The SPM group analysis results for right finger-thumb tapping task for (a and c) both the acquired signal and (b and d) the derived neuronal component. (a and b) The oxyhemoglobin signal; (c and d) the deoxyhemoglobin signal.

Table 2

Cluster report for the derived neuronal components correlated with the right finger-thumb tapping task [Figs. 3(b) and 3(d)]. The columns X, Y, and Z represent the MNI coordinate for the peak of the cluster.

To further demonstrate the effectiveness of our approach, we measure the consistencies of the waveform and spatial pattern between the oxyhemoglobin and deoxyhemoglobin signals. For the waveform consistency measure, for each subject, we calculate the waveform consistency using the channel with the highest overall activation within the ROI, since we suppose this channel is the most representative of the task-related response. To avoid bias, this selection is made independently for the acquired and derived neuronal components. The left columns in Table 3 show the error in waveform consistency between oxyhemoglobin and deoxyhemoglobin. The oxyhemoglobin and deoxyhemoglobin waveforms are significantly ($P<0.000003$) more consistent in the derived neuronal component than in the acquired data. The right columns in Table 2 show the error in spatial pattern consistency between oxyhemoglobin and deoxyhemoglobin within the motor areas. The beta values of BOLD estimation for oxyhemoglobin and deoxyhemoglobin are significantly ($P<0.003$) more consistent for the derived neuronal component than for the acquired data.

Table 3

Waveform and spatial pattern consistency between oxyhemoglobin and deoxyhemoglobin data. The errors for both temporal and spatial consistency represent the discrepancy between the oxyhemoglobin and deoxyhemoglobin signals. The smaller error represents an improvement in consistency between the two signals.

Figure 4 presents the mean and standard error for the waveforms entered into the waveform consistency analysis. Note that Fig. 4(a) shows typical acquired fNIRS waveforms, where the oxyhemoglobin is much stronger than deoxyhemoglobin signal. Figure 4(b) shows that the global component has a waveform very close to the canonical BOLD response in a block design, consistent with previous findings.^5,13 Figure 4(b) also illustrates the challenge for removing the global mean due to its high correlation with the expected cortical BOLD signal. Figure 4(c) shows the derived neuronal components as represented by both the oxyhemoglobin and deoxyhemoglobin responses as expected based on known physiology for a task-related effect.

Fig. 4

The peak activation channel was determined for waveform consistency measure for each of the 16 subjects. The mean event-triggered average fNIRS signal (lines) and standard error range (shaded color) for both oxyhemoglobin (red) and deoxyhemoglobin (blue) signal are presented ($N=16$). (a) The acquired fNIRS signal [H in Eq. (5)] is shown. (b) The global component [${\mathrm{H}}_{\text{Global}}$ in Eq. (5)] is shown. (c) The derived neuronal component [${\mathrm{H}}_{\text{Neuronal}}$ in Eq. (5)] is shown.

4.

Discussion

Our algorithm to separate global from local neuronal effects is based on a novel spatial filtering approach. We present two lines of validation. First, what has been removed from the original data, i.e., the waveforms of the global component shown in Figs. 2(b) and 2(e), are uniform across all channels, as determined by the large sized spatial filter. Thus, there is little chance that such a component reflects task-related neural responses. Second, the results shown in Fig. 3 confirm that the spatial distributions of the oxyhemoglobin and deoxyhemoglobin are more similar for the derived cortical BOLD signal compared to that of the acquired data, consistent with the basic concept of task-related effects.

Our approach has several limitations because it relies on a highpass spatial filter with a large-sized kernel. First, any localized artifacts such as those introduced by large blood vessels cannot be removed. Second, a sufficiently large-sized recording coverage is required because the size of the spatial filter kernel is 50 deg. In other words, the area of activation must be a small portion of the recording coverage. In addition, any highpass spatial filter approach may result in artificial negative activity. To quantify those potential pitfalls, we performed a simulation where the true neuronal component (${\mathrm{H}}_{\text{Neuronal}}$) is a Gaussian patch with $\sigma =1\mathrm{cm}$, which represents a 523-voxel cluster in a NMI standard brain with $2\times 2\times 2\mathrm{mm}$ resolution. The global component in this simulation (${\mathrm{H}}_{\text{Global}}$) has a uniform intensity across the recorded area. The top row of Fig. 5 shows the simulation images or heat map of the ${\mathrm{H}}_{\text{Neuronal}}$, ${\mathrm{H}}_{\text{Global}}$, and acquired data (H). The middle row shows the heat map of derived neuronal components using various sizes of spatial filters. The bottom row shows the heat map of derived neuronal component using various coverage sizes (e.g., the range of optode placement). We show the result for each coverage size with a large filter size ($\sigma =5\mathrm{cm}$). Mathematically, the effective filter size is limited to be smaller than the coverage size. Figure 5 shows that when either the size of the spatial filter kernel is not sufficient (e.g., $\sigma <3\mathrm{cm}$) or the coverage size is not sufficient (e.g., size $<9\mathrm{cm}$), the size of the derived neuronal component is underestimated and surrounded by artificial negative activity.

Fig. 5

Simulated fNIRS activity assuming the true neuronal activity ${\mathrm{H}}_{\text{Neuronal}}$ is a single positive Gaussian patch: $H_{\text{Neuronal}}=e^{-(x^2+y^2)/(2{\sigma }^2)}$, where $\sigma =1\mathrm{cm}$. The global component ${\mathrm{H}}_{\text{Global}}$ is assumed to be constant. The top row shows the simulated spatial pattern of ${\mathrm{H}}_{\text{Neuronal}}$, ${\mathrm{H}}_{\text{Global}}$, and acquired data H. The middle row shows the derived neuronal component ${\mathrm{H}}_{\text{Neuronal}}$ using various spatial filter sizes. The bottom row shows the derived neuronal component ${\mathrm{H}}_{\text{Neuronal}}$ using various coverage sizes with the largest possible filter size. For all panels, each pixel or color block represents a 1-cm square.

In summary, when the size of the smoothing kernel is sufficient, there will be no localized artificial negative activity in the derived neuronal component. Therefore, both the size of the smoothing kernel and the coverage size are determined by the spatial extent of the expected activity, which can often be estimated a priori.

5.

Conclusion

We have developed a PC spatial filtering technology to address the challenge of separating the neuronal from global components. Large vessel-induced local artifacts remain a challenge. Our approach may provide a promising tool to enhance the neuroscience application for fNIRS, when the optode coverage is widely distributed and the expected effects are specifically localized.