Separation of the global and local components in functional near-infrared spectroscopy signals using principal component spatial filtering

Abstract. Global systemic effects not specific to a task can be prominent in functional near-infrared spectroscopy (fNIRS) signals and the separation of task-specific fNIRS signals and global nonspecific effects is challenging due to waveform correlations. We describe a principal component spatial filter algorithm for separation of the global and local effects. The effectiveness of the approach is demonstrated using fNIRS signals acquired during a right finger-thumb tapping task where the response patterns are well established. Both the temporal waveforms and the spatial pattern consistencies between oxyhemoglobin and deoxyhemoglobin signals are significantly improved, consistent with the basic physiological basis of fNIRS signals and the expected pattern of activity associated with the task.

Separation of the global and local components in functional near-infrared spectroscopy signals using principal component spatial filtering Xian Zhang, a, * Jack Adam Noah, a and Joy Hirsch a,b,c,d

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, 1 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][3][4][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][8][9][10][11][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, 13 continuous blood pressure, 5 or the fNIRS signal from skin using a pair of optodes <1 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 14 and temporal filtering methods, 10 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][17][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.

Subjects
Twenty-two healthy subjects participated in a finger-thumb tapping experiment. Among them, six subjects were removed prior to process due to nonphysiological waveforms on more than 30% of the channels, probably resulting from optode movement. Sixteen subjects (11 female/5 male, mean age 24.5 AE 9.0 years) were used for data analysis. Written, informed consent was obtained from all participants in accordance with guidelines approved by the Yale University Institutional Review Board.

Neuroimaging
Signals were acquired using the LABNIRS system (Shimadzu Corp., Kyoto, Japan). Thirty emitter and 29 detector optodes were used, providing a total of 98 acquisition channels distributed across both hemispheres. Each emitter optode connects to laser diodes at three wavelengths (780, 805, and 830 nm) used to discriminate between oxyhemoglobin and deoxyhemoglobin signals. Raw optical density changes were converted into relative chromophore concentration changes through the modified Beer-Lambert law, using standard conversion terms. 19 Signals were sampled every 9.3 ms.

Task
The task used in this study was right finger-thumb tapping. Subjects were instructed to tap fingers and thumb of their right hands repeatedly during the task epochs. A standard 15-s task and 15-s rest block design were used, as previously 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." Zhang, Noah, and Hirsch: Separation of the global and local components in functional. . . described for functional mapping of handmotor areas using fMRI. 20 Each run consisted of six blocks and two runs were performed for a total of 6 min.

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

Functional Near-Infrared Spectroscopy Data Processing
Baseline drift was removed with the wavelet detrending algorithm procedure provided in NIRS_SPM. 21 The first-level general linear model (GLM) analysis 22 was performed using NIRS_SPM. The beta values, i.e., the amplitude of the BOLD signal, for each subject, were transformed into 3-D volume using xjview. 23 Finally, the group analysis on the 3-D data set was performed using SPM12. 22 2.6 Global Component Removal

Decomposition of the signal
We denote by H the n-by-p matrix of fNIRS acquired signals consisting of waveforms for multiple recording channels. Each row of H represents a sample point in time, and each column represents a recording channel. H can be decomposed into PCs using the standard singular value decomposition algorithm: E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 1 ; 3 2 6 ; 6 6 1 Here, U is effectively an n-by-p matrix, where each row corresponds to sample point and each column is the waveform of one PC across time. 24 V is the p-by-p matrix containing spatial information about the PCs. Each column in V corresponds to one PC, and each entry in the column measures the strength or coefficients of that PC at a certain channel. The values in the diagonal matrix Σ represent the relative importance of each PC. Isolating the temporal and spatial patterns of fNIRS data in separate matrices-U and V-allows us to extract the spatial pattern of global components.

Spatial smoothing
To obtain the global component, we removed any localized signal in the spatial pattern of each PC using Gaussian kernel convolution, which is the standard approach for both volume (FMRI) and two-dimensional image processing. Since fNIRS channels are distributed across the curved surface of the skull instead lying of on a flat plane, the distance r between two channels is defined by the arc length of the shortest path joining them along a sphere centered at the anterior commissure. The Gaussian filter is defined as E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 2 ; 3 2 6 ; 3 9 1 G 2D ðrÞ ¼ expð−r 2 ∕2σ 2 Þ; (2) where σ is a parameter controlling the width of the kernel. This width should be sufficiently greater than that of a cortical neuronal activation but sufficiently smaller than that of the global components. Our choice of σ ¼ 50 deg was based on pilot data and, to avoid overfitting, was not optimized for the testing dataset.
Convolution of the spatial pattern with this Gaussian filter eliminates the localized neuronal activity patterns, leaving only the global component intact. The spatial pattern of the global component of each PC can therefore be calculated with E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 3 ; 3 2 6 ; 6 9 7 where v i denotes the i'th column of V, or the spatial pattern of the i'th PC of the data. È is the convolution operator. Note that because the channels do not lie in a flat plane and are discrete, the convolution is done with the discrete convolution. V Ã is the convolved or smoothed spatial patterns of PCs and should only contain the spatial pattern of global components.

Reconstruction of the signal
The waveforms of the global component of the data can be calculated by plugging the smoothed spatial pattern matrix V Ã back into the singular value decomposition formula: E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 4 ; 3 2 6 ; 5 3 7 Subtracting this result from the acquired data yields the localized derived neuronal signal: E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 5 ; 3 2 6 ; 4 7 8

Errors of Waveform and Spatial Pattern Consistency Between Oxyhemoglobin and Deoxyhemoglobin Signals
We can check the validity of the derived H Neuronal based on the assumption that, for a true BOLD response, oxyhemoglobin and deoxyhemoglobin signals exhibit consistent spatial activation pattern and consistent temporal waveform. To estimate the waveform consistency, we use the peak signal within the channels of the motor area (the enclosed area in Fig. 1). We assume that oxyhemoglobin signal is linearly related to the deoxyhemoglobin signal. Therefore, the "ideal" oxyhemoglobin signal can be described using deoxyhemoglobin signal as the model [Eq. (6)]: E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 6 ; 6 3 ; 6 0 2 Oxy waveform ¼ B Ã Dexoy waveform þ ε; (6) where Deoxy waveform and Oxy waveform are vectors representing the waveforms of the event-triggered average with mean value removed. Each number in the array is a sample point of the waveform. B is a scalar parameter that minimizes the residue error ε. The error of waveform consistency between the oxyhemoglobin and deoxyhemoglobin signals is measured as the root-mean square of ε.
We can use a similar model to test the spatial consistency of oxyhemoglobin and deoxyhemoglobin activation. If Oxy spatial and Deoxy spatial are vectors containing the beta values calculated by SPM GLM for each channel within the ROI, then the residual error ε is calculated with E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 7 ; 6 3 ; 4 3 9 Oxy spatial ¼ B Ã Deoxy spatial þ ε: When regions with high oxyhemoglobin activation also exhibit high deoxyhemoglobin levels, as is the case with a true cortical BOLD signal, we expect the root-mean square of ε to be small.

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

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 (H Neuronal ) is a Gaussian patch with σ ¼ 1 cm, which represents a 523-voxel cluster in a NMI standard brain with 2 × 2 × 2 mm resolution. The global component in this simulation (H Global ) has a uniform intensity across the recorded area. The top row of Fig. 5 shows the simulation images or heat map of the H Neuronal , H 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 (σ ¼ 5 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., σ < 3 cm) or the coverage size is not sufficient (e.g., size <9 cm), the size of the derived neuronal component is underestimated and surrounded by artificial negative activity.
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.

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.