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.23.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,22.214.171.124.–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 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 analysis14 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.1617.–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 component3,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.
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 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.
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.
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 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 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.
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.”
|37||7||20||0.61||6||Premotor and supplementary motor cortex|
|38||7||42||0.74||6||Premotor and supplementary motor cortex|
|39||11||59||0.54||6||Premotor and supplementary motor cortex|
|0.34||9||Dorsolateral prefrontal cortex|
|40||10||70||0.81||6||Premotor and supplementary motor cortex|
|0.19||8||Includes frontal eye fields|
|49||55||0.92||6||Premotor and supplementary motor cortex|
|50||67||1||6||Premotor and supplementary motor cortex|
|51||76||1||6||Premotor and supplementary motor cortex|
|60||64||0.62||4||Primary motor cortex|
|0.27||3||Primary somatosensory cortex|
|61||76||0.86||6||Premotor and supplementary motor cortex|
|71||70||0.64||4||Primary motor cortex|
|0.3||3||Primary somatosensory cortex|
|72||80||0.87||4||Primary motor cortex|
|30||18||69||0.66||6||Premotor and supplementary motor cortex|
|0.34||8||Includes frontal eye fields|
|0.15||22||Superior temporal gyrus|
|0.14||6||Premotor and supplementary motor cortex|
|38||1||39||0.61||6||Premotor and supplementary motor cortex|
|0.22||4||Primary motor cortex|
|39||5||59||0.94||6||Premotor and supplementary motor cortex|
|40||4||72||0.98||6||Premotor and supplementary motor cortex|
|0.26||2||Primary somatosensory cortex|
|49||52||0.4||4||Primary motor cortex|
|0.31||6||Premotor and supplementary motor cortex|
|0.28||3||Primary somatosensory cortex|
|50||70||0.8||6||Premotor and supplementary motor cortex|
|51||78||1||6||Premotor and supplementary motor cortex|
|61||75||0.55||4||Primary motor cortex|
|0.45||6||Premotor and supplementary motor cortex|
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) analysis22 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
Global Component Removal
Decomposition of the signal
We denote by the -by- matrix of fNIRS acquired signals consisting of waveforms for multiple recording channels. Each row of represents a sample point in time, and each column represents a recording channel. can be decomposed into PCs using the standard singular value decomposition algorithm:
Here, is effectively an -by- matrix, where each row corresponds to sample point and each column is the waveform of one PC across time.24 is the -by- matrix containing spatial information about the PCs. Each column in 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— and —allows us to extract the spatial pattern of global components.
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 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
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
Reconstruction of the signal
The waveforms of the global component of the data can be calculated by plugging the smoothed spatial pattern matrix back into the singular value decomposition formula:
Subtracting this result from the acquired data yields the localized derived neuronal signal:
Errors of Waveform and Spatial Pattern Consistency Between Oxyhemoglobin and Deoxyhemoglobin Signals
We can check the validity of the derived 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)]:
We can use a similar model to test the spatial consistency of oxyhemoglobin and deoxyhemoglobin activation. If and are vectors containing the beta values calculated by SPM GLM for each channel within the ROI, then the residual error is calculated with
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.
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)].
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 (). 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 results25,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.
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.
|1||50||28||1||46||Dorsolateral prefrontal cortex|
|2||52||0.61||6||Premotor and supplementary motor cortex||4.34|
|3||60||0.25||6||Premotor and supplementary motor cortex||4.46|
|4||70||0.3||3||Primary somatosensory cortex|
|1||12||1||22||Superior temporal gyrus|
|3||34||40||0.72||9||Dorsolateral prefrontal cortex|
|4||16||0.82||22||Superior temporal gyrus|
|5||10||64||0.52||6||Premotor and supplementary motor cortex||5.93|
|6||18||32||56||0.9||8||Includes frontal eye fields||3.83|
|7||80||0.96||4||Primary motor cortex||3.41|
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 () 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 () more consistent for the derived neuronal component than for the acquired data.
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.
|Error for waveform consistency||Error for spatial pattern consistency|
|Subject||Acquired||Derived neuronal||Acquired||Derived neuronal|
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.
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 () is a Gaussian patch with , which represents a 523-voxel cluster in a NMI standard brain with resolution. The global component in this simulation () has a uniform intensity across the recorded area. The top row of Fig. 5 shows the simulation images or heat map of the , , 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 (). 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., ) or the coverage size is not sufficient (e.g., size ), 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.
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.
The research reported in this publication was partially supported by the National Institute of Mental Health of the National Institutes of Health under award number R01MH107513 (PI: Joy Hirsch). The content is solely the responsibility of authors and does not necessarily represent the official views of the National Institutes of Health.
Xian Zhang received his PhD degree in psychology and visual science from Columbia University in New York in 2003. He is an associate research scientist in the Brain Function Laboratory in the Department of Psychiatry, Yale School of Medicine. His research interests include computational neuroscience, signal processing, and neuroimaging technologies, such as EEG, fNIRS, and fMRI and their applications in psychiatry, vision science, social interactions, and decision making.
Jack Adam Noah received his PhD degree in biomedical sciences from the Marshall University School of Medicine in 2007. He is an associate research scientist at Yale School of Medicine in the Department of Psychiatry and the Brain Function Laboratory. His research interests include functional near-infrared spectroscopy and integration of other multimodal and behavioral recording techniques for applications in communication and social interactions, neurofeedback, and cognitive neuroimaging.
Joy Hirsch received her PhD degree in psychophysics and visual science from Columbia University and is now a professor of psychiatry and neurobiology, Yale School of Medicine, and a professor of neuroscience, University College London. She is also the director of the Brain Function Laboratory at Yale University. Her research is focused on investigations of neural circuitry that underlies human social interactions using multimodal neuroimaging techniques including fNIRS, fMRI, EEG, eye-tracking and behavioral measures. Prior to recruitment to Yale, she was a director of the fMRI research center at Columbia University.