Expand this Topic clickable element to expand a topic
Skip to content
Optica Publishing Group

Wavelet-based method for removing global physiological noise in functional near-infrared spectroscopy

Open Access Open Access

Abstract

Functional near-infrared spectroscopy (fNIRS) is a fast-developing non-invasive functional brain imaging technology widely used in cognitive neuroscience, clinical research and neural engineering. However, it is a challenge to effectively remove the global physiological noise in the fNIRS signal. The global physiological noise in fNIRS arises from multiple physiological origins in both superficial tissues and the brain. It has complex temporal, spatial and frequency characteristics, casting significant influence on the results. In the present study, we developed a novel wavelet-based method for fNIRS global physiological noise removal. The method is data-driven and does not rely on any additional hardware or subjective noise component selection procedure. It consists of two steps. Firstly, we use wavelet transform coherence to automatically detect the time-frequency points contaminated by the global physiological noise. Secondly, we decompose the fNIRS signal by using the wavelet transform, and then suppress the wavelet energy of the contaminated time-frequency points. Finally, we transform the signal back to a time series. We validated the method by using simulation and real data at both task- and resting-state. The results showed that our method can effectively remove the global physiological noise from the fNIRS signal and improve the spatial specificity of the task activation and the resting-state functional connectivity pattern.

© 2018 Optical Society of America under the terms of the OSA Open Access Publishing Agreement

1. Introduction

Because of low-cost, portability and high ecological validity, functional near-infrared spectroscopy (fNIRS) becomes a fast-developing non-invasive functional brain imaging technology. fNIRS has been widely used in various fields, such as cognitive neuroscience [1], clinical research [2] and neural engineering [3]. fNIRS measures the hemodynamic activity in the cerebral cortex. A typical fNIRS measurement channel consists of a near-infrared light source and a light detector placed on the scalp, forming a “banana-shaped” light path in the head due to the biological tissue’s diffusional scattering property. The distance between the source and the detector (S-D distance) is generally around 3 cm, which ensures the “banana-shape” reach adequate depth from the scalp and pass through the cerebral cortex [4]. Because the oxygenated hemoglobin (HbO) and the deoxygenated hemoglobin (HbR) have distinct absorption spectrum to the near-infrared light, so the concentration change parameter of the HbO and the HbR in the cerebral blood can be calculated from the intensity change of dual wavelength near-infrared lights using the modified Beer-Lambert law (MBLL) [5].

However, in practice it is a major challenge that the fNIRS measurement is often contaminated by the global physiological noise. This globally distributed physiological noise usually casts significant influence on various kinds of experiment. For example, in the task-based study, it may cover up the evoked neural activity from being detected [6–8], whereas it may lead to spurious global connectivity pattern in the resting-state functional connectivity study [9, 10]. The global physiological noise in fNIRS mainly arises from the physiological fluctuation of the blood supply system in the superficial tissue layers including the scalp and the skull, and it also arises from the physiological vascular oscillation inside the brain [11–13]. It is usually the dominant component of the fNIRS signal because the photon density of the “banana-shaped” light path decays severely with the depth increases, so the fNIRS signal is much more sensitive to the hemoglobin concentration change in the superficial layers than that in the cortical layers (about 10—20 times [12]). Therefore, essentially it is a difficult weak signal detection problem to recover the neural activity from the fNIRS signal contaminated by the global physiological noise.

Several methods have been developed to remove the fNIRS global physiological noise (see [14] for a review). One frequently-used method is the band-pass filtering which removes the noise by suppressing the power in the noise frequency band [15]. The problem of this approach is that it is necessary to know the exact frequency characteristics of all physiological noise which may vary over time, space and across individuals. Another issue is if the frequency band of the noise overlaps with that of the neural activity, the filter may destroy the experimental effects of interest.

The short source-detector distance channel (sSD channel) based correction is a kind of effective method for removing the fNIRS global physiological noise [8, 13, 16–20]. By using relatively short S-D distance (about 0.5 to 1.5 cm [12, 21]), the “banana-shaped” light path can be constrained within the superficial layers rather than reaching the cerebral cortex. Therefore, the signal from the sSD channel is usually assumed as a good measure of the superficial physiological noise and is used to remove the superficial physiological noise in the normal S-D distance channel. Although this method is very effective, it has several practical shortcomings. Firstly, not all fNIRS apparatus support sSD channel. Secondly, configuring sSD channel increase the setting up time. Thirdly, the superficial physiological noise is not homogeneously distributed on the head [11], but the number of fNIRS channels is usually limited. It is impossible to allocate each normal channel with an additional sSD channel to record the noise. It is still a difficult problem to optimize the arrangement of the limited sSD channels [18, 21].

The blind source separation (BSS) algorithms such as principal component analysis (PCA) [22–25] and independent component analysis (ICA) [26, 27] are another option for global physiological noise removal. These methods firstly decompose the fNIRS signal into a series of source components with either the orthogonal (PCA) or the statistical independent (ICA) assumption. Then the noise components are identified by means such like energy sorting or spatial uniformity checking, based on different assumptions about the noise. Finally, the signal is separated from the noise and reconstructed. The BSS methods are good alternatives for global physiological noise removal, especially when sSD method is unavailable. However, the performance of these methods highly depends on the noise component identification procedure, which is usually under human’s supervision or even needs manual selection which requires rich experience from the user. Another potential problem of the BSS methods is that they cannot handle the time lags among the physiological noise of different channels.

The wavelet transform (WT) is a powerful mathematical tool for non-stationary signals’ time-frequency analysis. The wavelet transform can expand a time series into time-frequency space by using both time and frequency localized basis functions (wavelet basis). It uses low time resolution at low frequencies and high time resolution at high frequencies and thus is time-frequency adaptive. Based on WT, the wavelet transform coherence (WTC) can further give the coherence between two signals in the time-frequency space. In previous studies, WT and WTC were widely used in fNIRS data preprocessing and analysis. For example, WT was used in fNIRS signal detrending [28] and motion-induced artifact removal [29–32]. While WTC was used to calculate the time-frequency coherence between the fNIRS signal and the peripheral physiology recording signal to identify and quantify the physiological components in fNIRS signal [33, 34], or between the fNIRS data and the predicted brain response to detect task activation [35]. Moreover, in recent years, WTC was also frequently used in fNIRS hyperscanning studies to investigate the inter-brain synchronization during social interaction between two or more brains (e.g, [36–39]).

However, to our knowledge, the wavelet-based fNIRS global physiological noise removal remains vacant. In the present study, we proposed a novel wavelet-based method to remove the fNIRS global physiological noise. Our method consists of two steps. In the first step, we use WTC to automatically detect the time-frequency points contaminated by the global physiological noise. In the second step, the fNIRS signal is decomposed with WT, the wavelet energy of the contaminated time-frequency points was then suppressed, and the signal was finally reconstructed. We conducted simulation and real experiments of both task- and resting-state to validate the performance of our method.

2. Method

2.1 Wavelet transform and wavelet transform coherence

Wavelet transform is a powerful time-frequency analysis technique. A wavelet is a function with zero mean and is both frequency and time localized. For example, the Morlet wavelet, employed in the present study and other fNIRS studies [33–35] for its good performance in time-frequency analysis, is defined as

ψ0(η)=π1/4eiω0ηe12η2,
where ω0 is dimensionless frequency (ω0 > 5) and η is dimensionless time.

ψ0(η) is usually known as a “mother” wavelet. By adjusting the scale (s) of the mother wavelet (η=st) and normalizing it to have unit energy, a series of wavelet basis can be obtained. The wavelet transform can be viewed as applying a bandpass filter to the time series, in which the time series (xn, n=1,2,,N) is convolved with the scaled and normalized wavelets (usually in Fourier space for fast speed [40]), writing as

Wxn(n,s)=δtsn'=1Nxn'ψ0[(n'n)δts],
where δt is uniform time steps and |Wxn(n,s)|2 is defined as the wavelet power [41]. The above wavelet transform also has its corresponding inverse wavelet transform (IWT) which can reconstruct the signal using time- and frequency- localized components (for mathematical details, see [42] and [43]).

Based on the wavelet transform, the wavelet transform coherence can give the coherence between two time series at each point of the wavelet-divided time-scale plane, and thus can be used to analyze the co-varying characteristics between two time series with both temporal and frequency resolution. The wavelet transform coherence between two time series xn and yn is defined as

R2(n,s)=|S[s1Wxn,yn(n,s)]|2S[s1|Wxn(n,s)|2]S[s1|Wyn(n,s)|2],
where Wxn,yn=WxnWyn* and * denotes complex conjugation. |Wxn,yn|2 is defined as the cross wavelet power. S is a smoothing operator consisted of a temporal smoother and a scale smoother [44]. The computational implements of WTC and WT/IWT have been provided in the MATLAB software (R2017a, MathWorks Inc., Massachusetts, USA).

2.2 The WT-based method for global physiological noise removal

Generally, the basic concept of the wavelet denoising is to decompose the signal into time-frequency space, and then identify those noise-related time-frequency points by using appropriate criterion and set their wavelet coefficients to zero, and finally transform it back to a time series. Our method uses the same concept. The key point is that we use WTC to identify those globally co-varying time-frequency points as contaminated by the global physiological noise. In our method, the global physiological noise is removed per channel. For each channel, the noise is removed by using a two-stage procedure. In stage one, the time-frequency points contaminated by the global physiological noise are detected by using WTC. Specifically, firstly we calculate the time-frequency distributed WTC map (also called scalogram [35]) between the current denoising channel’s signal and every other channel’s unfiltered signal. These WTC maps are further binarized according to the significance of the WTC value at each time-frequency pixel. Then all these WTC maps are averaged, forming a globally co-varying time-frequency map. The value of a given pixel in this averaged map reflects how globally the current channel co-varies with other channels at this time-frequency point. Finally, a denoising mask for the current channel is generated by applying a threshold k to the globally co-varying time-frequency map. In stage two, the signal of the current denoising channel is decomposed into the time-frequency space by using WT. Then the wavelet energy at those noise-contaminated time-frequency points is suppressed by applying the mask obtained in stage one to the wavelet coefficients. Finally, the signal is reconstructed by using the inverse wavelet transform. The above two stages are repeated per channel to complete the global physiological noise removal. (Fig. 1).

 figure: Fig. 1

Fig. 1 The flow chart of the method. The chart illustrates the two-stage procedure for removing a single channel’s global physiological noise. The black section of the mask is the blocked portions. The procedure is repeated per channel to complete the entire noise removal processing.

Download Full Size | PDF

This method utilizes the global distribution property of the physiological noise. Using this method, the time-frequency components with high globally co-varying feature are labeled as contaminated and removed, while the interested neural components are reserved because they are relatively local. Also note that because the WTC can track the co-varying feature of two signals even if there is a time delay between them, so our method can deal with the possible time lags between the physiological noise of different channels.

The selection of the threshold k is important. Indeed, k is a quantification definition of the globality of the global physiological noise. Ideally, k should be equal to 100% if assuming the global physiological noise is significantly coherent among all channels. However, the global physiological noise is not completely spatial homogeneous. Therefore, too strict definition of the globality (high k) can decrease the noise removal performance of the algorithm. On the other hand, too loose definition of the globality (low k) may lead loss of the interested neural activity because the neural components within a same functional area may be also coherent. Therefore, theoretically, if assuming the neural components are coherent within the entire interested brain regions, the threshold k should satisfy

k>NROINT,
where NROI indicates the number of the interested channels (regions of interests, ROI) and NT means the number of all the channels. In the present study, we used a simulated experiment with set/known number of NROI and NT to validate this lower bound of k. We also conducted a real experiment to test the influence of different k value selection on the denoising performance.

3. Experiments

3.1 Simulation experiment

In the simulation experiment, we simulated a 66-channel fNIRS measurement (6 rows by 11 columns) on the head surface. The measurement contained a 14-channel ROI (region of interest). All the channels were filled with simulated global physiological noise consisted of multiple sinusoidal signals of different frequency, including the heart beat (1 Hz), the respiration (0.3 Hz), the Mayer waves (around 0.1 Hz) and the very low-frequency oscillations (below 0.01 Hz) [35, 45, 46]. Moreover, to test whether our method can deal with the possible time delay of the global physiological noise spreading among different measurement channels on the head surface, a maximum of 6 s time delay was simulated, spreading from an origin point (the bottom-middle channel) outward to the other channels. We set the time lag of a channel according to its Chebyshev distance (d1,2 = max(|x1 - x2|, |y1 - y2|)) from the origin point (Fig. 2(A)).

 figure: Fig. 2

Fig. 2 The simulation experiment. (A) The simulated 66-channel fNIRS measurement (6 rows by 11 columns). All the channels were filled with simulated global physiological noise. The red frame indicates the ROI in which the simulated neural activity was added. A maximum of 6 s time-lag of the physiological noise was simulated, spreading from the bottom-middle channel to the outer channels as indicated by grey scale. (B) The simulated task-evoked neural signal (left) and the spontaneous resting-state neural signal (right) as the ground truth. (C) The noise contaminated signal (the blue line). The ground truth signal was also plotted (the black lines). Note the extremely low signal-to-noise ratio of the mixture signal.

Download Full Size | PDF

Both task- and resting-state data set were simulated to validate the proposed method. The task-evoked neural activity signal was simulated as a boxcar function convolved with the canonical hemodynamic response function (HRF). The simulated resting-state spontaneous neural activity signal was a combine of a series of sinusoidal functions with random phase. The frequency range of these sinusoidal functions was from 0.01 Hz to 0.1 Hz, stepped by 0.001 Hz, with amplitude distributed in 1/f manner [45] (Fig. 2(B)). The simulated task- and resting-state signals were set only in the ROI channels and were mixed with the simulated global physiological noise. The amplitude ratio of the neural signal and the physiological noise was set to from 1: 8 to 1: 10, randomly across channels. A 10 dB random Gauss white noise was also added to each channel’s simulated signal. The length of the simulated time series was 300 s, with a sampling rate of 10 Hz (Fig. 2(C)).

Both the task- and resting-state data sets were denoised by using the proposed method with different globality threshold k to remove the global physiological noise. The recovered signals were then compared with the ground truth signals. To examine the lower bound of the globality threshold parameter k, the recovered signal in ROI derived from each k was then regressed by the ground truth signal to see whether the recovered signal can well reproduce the ground truth signal.

3.2 Real experiment

We then conducted a finger-tapping motor task experiment and a resting-state experiment to further validate the proposed method. Forty healthy right-handed young adult participants (21.7 ± 2.5 years of age, 22 males and 18 female) were recruited from Shenzhen University. All participants underwent a 10-minute resting-state session followed by a 5.6-minute task session of fNIRS recording. In the resting-state session, the participants were instructed to keep still with their eyes closed and relax their mind. In the task session, participants performed a sequential bilateral finger-tapping task containing seven task blocks with a pseudo-randomized block length of 20–30 s [47]. Informed consent was obtained from all subjects. The study protocol was approved by the Institutional Review Board at Shenzhen Key Laboratory of Affective and Social Neuroscience, Shenzhen University.

The fNIRS measurement was conducted by using a NIRScout continuous wave fNIRS system (NIRx Medical Technologies, USA). The absorptions of the near-infrared lights at two wavelengths (785 nm and 830 nm) were measured with a sampling rate of 7.8125 Hz. Two 4 × 4 probe sets were used in this study with each probe set consisting of eight laser sources and eight detectors, forming 24 measurement channels (48 channels in total). The source-detector distance was 30 mm. The two probe sets were respectively placed on the head with their center at C3 and C4 of the international 10–20 system [48] to cover the bilateral sensorimotor areas. The cortex localization of the channels was obtained by using a 3-dimentional digitizer and the NIRS-SPM software [49, 50]. The oxygenated (HbO) and the deoxygenated (HbR) signals were calculated with the modified Beer–Lambert law [5], with differential pathlength factor (DPF) of 7.25 and 6.38 for 785 nm and 830 nm, respectively [51].

All the fNIRS signals were firstly detrended by fitting the first and second order polynomial functions to remove the linear and bilinear trends [27]. Then the global physiological noise was removed by using the proposed method. For the task activation calculation, the general linear model (GLM) approach was used. The regressor was made by convolving the block design with the canonical HRF. The group-level activation t-map was derived by conducting a one sample t-test on all individual β values [47, 52]. For the resting-state functional connectivity (RSFC) calculation, the spontaneous neural activity was extracted by using a 0.01—0.08 Hz band-pass filter [53]. The seed channels (the sensorimotor area, left channel 5 and 9) were determined according to the task activation results. The average time course of the two seed channels was used as the seed time course to calculate the Pearson’s correlation coefficients with the time course of each channel. The Fisher’s z-transform was applied to the Pearson’s correlation coefficients to increase the Gaussianility [54]. Then the group-level RSFC t-map was derived by using the one sample t-test [27]. For comparison, the task activation map and the RSFC map were also calculated without denoising. The receiver operating characteristic (ROC) curve and the area under the ROC curve (AUC) [55] were used to evaluate the performance of our method in detecting the task activation and the RSFC of the sensorimotor area. The sensorimotor template derived from the cortex localization was used as the ground truth to draw the ROC curve [56]. By using varied k values, we also explored the influence of the k value selection on the denoising performance in the real fNIRS data set.

4. Results

4.1 Simulation experiment

Figure 3 illustrates the recovered neural signal (the red line) and the corresponding ground truth signal (the black line) of a randomly selected ROI channel for the task- and resting-state data set, respectively, using k of 50%. The Pearson’s correlation coefficient between the reconstructed signal and the ground truth signal was 0.60 for the task-state and 0.43 for the resting-state, respectively. The results showed that the recovered signal well reproduced the ground truth signal (Fig. 3). The high-frequency oscillation in the recovered signal was the residual Gauss white noise which can be easily removed by low-pass filtering.

 figure: Fig. 3

Fig. 3 The recovered neural signal. The red line indicates the recovered neural signal and the black line indicates the ground truth signal. Only for display purpose, a 5-point moving average of the recovered signal was applied to control the amplitude of the Gauss white noise.

Download Full Size | PDF

The denoising results showed that when k was greater than or equal to the ratio of the number of ROI channels to the number of total channels (20% in our simulation experiment), the recovered signal can well reproduce the ground truth signal (Fig. 4). However, when k is less than this ratio, the recovered signal cannot reproduce the ground truth signal. These results suggested that the globality threshold k should at least be the ratio of the area of the interest to the area of the total measurements to avoid damaging the interested neural activity.

 figure: Fig. 4

Fig. 4 The denoising results of different globality threshold k in the simulation experiment. The vertical coordinate was the p-value of the regression. The horizontal coordinate was the globality parameter k, represented by a ratio (the ratio means for the currently denoising channel, with how many channels of the other 65 channels a time-frequency component co-varies should it be labeled as contaminated in the present study). Note that when k was less than 20% (13/65, the number of the ROI channels divided by the number of all channels, except the currently denoising ROI channel), the recovered signal cannot significantly reproduce the ground truth signal (the dashed line indicates the significance level of 0.05).

Download Full Size | PDF

4.2 Real experiment

The activation maps (Fig. 5, left panel) and the RSFC maps (Fig. 5, right panel) derived from the HbO signal with using our method to remove the global physiological noise (Fig. 5, second row) showed more consistency with the sensorimotor template (the first row) compared to those maps without using our method (Fig. 5, third row). ROC curve analysis showed that higher AUC indices in both the task activation map and the RSFC map derived from the corrected data using our method than those from the uncorrected data (0.81 vs. 0.78 for task activation, and 0.82 vs. 0.71 for RSFC; Fig. 6).

 figure: Fig. 5

Fig. 5 The task activation maps and the RSFC maps. The left column and the right column respectively show the group-level task activation results and the RSFC results, derived from the corrected HbO data (the second row) and the uncorrected HbO data (the third row), with the globality threshold k = 50%. The first row shows the sensorimotor template.

Download Full Size | PDF

 figure: Fig. 6

Fig. 6 The receiver operating characteristic (ROC) curve. The left and the right panel show the ROC curves of the task activation map and the RSFC map, respectively (globality threshold k = 50%). The red line indicates the ROC curve of the map derived from the corrected HbO data, and the blue line indicates the ROC curve of the map derived from the uncorrected HbO data.

Download Full Size | PDF

The denoising analyses with varied threshold k showed that the best performance of our WT-based method was at k = 50% in the present data set (Fig. 7). When k < 50%, performance of the WT-based method decreased as a function of the k value. This result was consistent with that of the simulation experiment. As we predicted in the theoretical analysis, when k > 50%, the performance of the algorithm decreased as well. This may be because the physiological noise is spatially heterogeneous, therefore an over-estimated globality threshold may be too strict for noise detection. Similar result could also be found on the HbR data (see Appendix).

 figure: Fig. 7

Fig. 7 The denoising performance of different globality threshold k on the real experiment. The vertical coordinate is the AUC index of the task activation map (the blue line) and the RSFC map (the red line) of the HbO data. The horizontal coordinate is the globality parameter k from 20% to 100%, stepped by 10%. The dashed lines show the AUC values derived from the uncorrected data.

Download Full Size | PDF

5. Discussion and conclusion

The current study developed a wavelet transform based method for fNIRS global physiological noise removal. The method consisted of two parts. The first part was the temporal-frequency global noise detection based on the wavelet transform coherence, and the second part was the signal decomposition and noise removal based on the wavelet transform and reconstruction. The results from the simulation experiment and the real experiment showed that the method can effectively remove the global physiological noise in fNIRS signal and improve the spatial specificity of the task activation and the resting-state functional connectivity.

Our method has several advantages. Firstly, compared with the traditional frequency-based filtering, the exact frequency distribution of the noise is unnecessary to know for our method because it adequately utilizes the global distribution property of the noise. The globally co-varying time-frequency components in the signal can be detected by using the wavelet transform coherence. Moreover, our method can also deal with the potential time variance of the noise based on the time-frequency analysis ability of the wavelet transform. Secondly, our method can deal with the time-lag of the physiological noise among different brain regions. The global physiological noise shows a significant spatially spreading pattern. Tong and Frederick found that the whole circulatory process of the low-frequency physiological signal averagely takes about 6 s spreading over the head [57]. And a recent study made by Zhang et al. [21] also showed that the physiological interference in superficial layers is highly consistent in a local brain region, but the consistency is significantly lower between distant brain regions. These findings suggested that the physiological noise is not time-aligned among different fNIRS channels, which may bring difficulty to the blind source separation methods such as ICA and PCA when extracting a common physiological noise component from all channels. In our method, the wavelet transform coherence can capture the co-varying characteristics between two signals with time delay between them and thus can solve this problem. Thirdly, our method is completely data-driven. It is independent of any extra hardware to record the physiological noise or manually picking up the noise component. It is simple to use and can be applied to most of the multi-channel fNIRS apparatus. These advantages make our method a promising preprocessing tool for fNIRS studies.

Our method also has some limitations. For example, the current method is based on utilizing the global property of the physiological noise. The globality threshold k plays a key role in the method. But the selection of k based on a-priori knowledge of the expected size of the ROI also relies on rich user experience. On one hand, selecting a too small k will reduce the specificity of the noise detection, might even lose the signal of the neural activity. According to the theoretical analysis and the experimental results, k should not be less than the ratio of the area of the interested functional region to the area of the total measurement. Therefore, the area of the ROI should be estimated as accurately as possible. This is relatively easy when the ROI is clearly determined. However, when the ROI is not explicit, for example, in some exploratory studies, the k value should be selected with caution. In this situation, a larger k may reduce the risk of removing the interested neural activity. On the other hand, a too large k will reduce the sensitivity of the noise detection. In the present study, the results from the real data suggested that a k of 50% can give the best performances for both task activation and RSFC detections. The results also indicated that higher k values (from 60% to 90%) could decrease the performance. However, one can still benefit from the method even if higher k values are used when compared to nonuse. In addition, the task activation AUC is relatively independent of k above 50% when compared to the RSFC AUC. This may be because the RSFC derived from HbO signal is more sensitive to the global physiological noise. Further studies are still needed to explore the optimal k value selection for different experiment types. Another limitation is that the current method can only be applied to off-line data. In the future, we plan to develop the on-line version of this method for real-time applications such as the brain-computer interface and the neurofeedback.

In conclusion, the present study proposed a promising method for fNIRS global physiological noise removal which has unique advantages. Our method can be applied to both task-based and resting-state fNIRS data, and its effectiveness was validated by using both simulation experiment and real experiment. Our method can be complementary to the sSD and the BSS methods to remove the global physiological noise, which represents a necessary evolution in fNIRS signal processing.

Appendix

Appendix listed the corresponding results of the HbR data (Fig. 8, 9 and 10).

 figure: Fig. 8

Fig. 8 The task activation maps and the RSFC maps of HbR data. The left column and the right column respectively show the group-level task activation results and the RSFC results, derived from the corrected HbR data (the second row) and the uncorrected HbR data (the third row), with the globality threshold k = 50%. The first row shows the sensorimotor template.

Download Full Size | PDF

 figure: Fig. 9

Fig. 9 The receiver operating characteristic (ROC) curves of HbR data. The left and the right panel show the ROC curve of the task activation map and the RSFC map, respectively (globality threshold k = 50%). The red line indicates the ROC curve of the map derived from the corrected HbR data, and the blue line indicates the ROC curve of the map derived from the uncorrected HbR data.

Download Full Size | PDF

 figure: Fig. 10

Fig. 10 The denoising performance of different globality threshold k on the real experiment (HbR data). The vertical coordinate is the AUC index of the task activation map (the blue line) and the RSFC map (the red line) of the HbR data. The horizontal coordinate is the globality parameter k from 20% to 100%, stepped by 10%. The dashed lines show the AUC values derived from the uncorrected data.

Download Full Size | PDF

Funding

National Natural Science Foundation of China (61503030, 31500920, 31530031, 81471376); the National Key Basic Research Program of China (973 Program, 2014CB744600); Natural Science Foundation of Shenzhen University (2017075); Shenzhen Peacock Program (827-000235); Key Program of the Natural Science Foundation of Tianjin (18JCZDJC36300); the Science and Technology Plan of Tianjin (14RCGFGX00847).

Acknowledgments

The authors are very grateful to the anonymous reviewers for their significant and constructive critiques and suggestions which greatly improved the manuscript. The authors also thank Zhonghua Tang and Xue Ma for their assistance in performing the experiment.

Disclosures

The authors declare that there are no conflicts of interest related to this article.

References and links

1. S. Cutini, S. B. Moro, and S. Bisconti, “Review: Functional near infrared optical imaging in cognitive neuroscience: an introductory review,” J. Near Infrared Spectrosc. 20(1), 75–92 (2012). [CrossRef]  

2. Z. Qian and W. Li, Functional Near-Infrared Spectroscopy (fNIRs) Clinical Applications Review (Life Science Instruments, 2013).

3. N. Naseer and K. S. Hong, “fNIRS-based brain-computer interfaces: a review,” Front. Hum. Neurosci. 9, 3 (2015). [PubMed]  

4. M. Ferrari and V. Quaresima, “A brief review on the history of human functional near-infrared spectroscopy (fNIRS) development and fields of application,” Neuroimage 63(2), 921–935 (2012). [CrossRef]   [PubMed]  

5. M. Cope and D. T. Delpy, “System for long-term measurement of cerebral blood and tissue oxygenation on newborn infants by near infra-red transillumination,” Med. Biol. Eng. Comput. 26(3), 289–294 (1988). [CrossRef]   [PubMed]  

6. S. B. Erdoğan, M. A. Yücel, and A. Akın, “Analysis of task-evoked systemic interference in fNIRS measurements: insights from fMRI,” Neuroimage 87, 490–504 (2014). [CrossRef]   [PubMed]  

7. E. Kirilina, A. Jelzow, A. Heine, M. Niessing, H. Wabnitz, R. Brühl, B. Ittermann, A. M. Jacobs, and I. Tachtsidis, “The physiological origin of task-evoked systemic artefacts in functional near infrared spectroscopy,” Neuroimage 61(1), 70–81 (2012). [CrossRef]   [PubMed]  

8. Q. Zhang, G. E. Strangman, and G. Ganis, “Adaptive filtering to reduce global interference in non-invasive NIRS measures of brain activation: How well and when does it work?” Neuroimage 45(3), 788–794 (2009). [CrossRef]   [PubMed]  

9. E. Sakakibara, F. Homae, S. Kawasaki, Y. Nishimura, R. Takizawa, S. Koike, A. Kinoshita, H. Sakurada, M. Yamagishi, F. Nishimura, A. Yoshikawa, A. Inai, M. Nishioka, Y. Eriguchi, J. Matsuoka, Y. Satomura, N. Okada, C. Kakiuchi, T. Araki, C. Kan, M. Umeda, A. Shimazu, M. Uga, I. Dan, H. Hashimoto, N. Kawakami, and K. Kasai, “Detection of resting state functional connectivity using partial correlation analysis: A study using multi-distance and whole-head probe near-infrared spectroscopy,” Neuroimage 142, 590–601 (2016). [CrossRef]   [PubMed]  

10. B. R. White, A. Z. Snyder, A. L. Cohen, S. E. Petersen, M. E. Raichle, B. L. Schlaggar, and J. P. Culver, “Resting-state functional connectivity in the human brain revealed with diffuse optical tomography,” Neuroimage 47(1), 148–156 (2009). [CrossRef]   [PubMed]  

11. L. Gagnon, R. J. Cooper, M. A. Yücel, K. L. Perdue, D. N. Greve, and D. A. Boas, “Short separation channel location impacts the performance of short channel regression in NIRS,” Neuroimage 59(3), 2518–2528 (2012). [CrossRef]   [PubMed]  

12. T. Takahashi, Y. Takikawa, R. Kawagoe, S. Shibuya, T. Iwano, and S. Kitazawa, “Influence of skin blood flow on near-infrared spectroscopy signals measured on the forehead during a verbal fluency task,” Neuroimage 57(3), 991–1002 (2011). [CrossRef]   [PubMed]  

13. Q. Zhang, E. N. Brown, and G. E. Strangman, “Adaptive filtering for global interference cancellation and real-time recovery of evoked brain activity: a Monte Carlo simulation study,” J. Biomed. Opt. 12(4), 044014 (2007). [CrossRef]   [PubMed]  

14. S. Tak and J. C. Ye, “Statistical analysis of fNIRS data: a comprehensive review,” Neuroimage 85(Pt 1), 72–91 (2014). [CrossRef]   [PubMed]  

15. T. J. Huppert, S. G. Diamond, M. A. Franceschini, and D. A. Boas, “HomER: a review of time-series analysis methods for near-infrared spectroscopy of the brain,” Appl. Opt. 48(10), D280–D298 (2009). [CrossRef]   [PubMed]  

16. R. B. Saager and A. J. Berger, “Direct characterization and removal of interfering absorption trends in two-layer turbid media,” J. Opt. Soc. Am. A 22(9), 1874–1882 (2005). [CrossRef]   [PubMed]  

17. L. Gagnon, K. Perdue, D. N. Greve, D. Goldenholz, G. Kaskhedikar, and D. A. Boas, “Improved recovery of the hemodynamic response in diffuse optical imaging using short optode separations and state-space modeling,” Neuroimage 56(3), 1362–1371 (2011). [CrossRef]   [PubMed]  

18. L. Gagnon, M. A. Yücel, D. A. Boas, and R. J. Cooper, “Further improvement in reducing superficial contamination in NIRS using double short separation measurements,” Neuroimage 85(Pt 1), 127–135 (2014). [CrossRef]   [PubMed]  

19. R. Saager and A. Berger, “Measurement of layer-like hemodynamic trends in scalp and cortex: implications for physiological baseline suppression in functional near-infrared spectroscopy,” J. Biomed. Opt. 13(3), 034017 (2008). [CrossRef]   [PubMed]  

20. R. B. Saager, N. L. Telleri, and A. J. Berger, “Two-detector Corrected Near Infrared Spectroscopy (C-NIRS) detects hemodynamic activation responses more robustly than single-detector NIRS,” Neuroimage 55(4), 1679–1685 (2011). [CrossRef]   [PubMed]  

21. Y. Zhang, F. Tan, X. Xu, L. Duan, H. Liu, F. Tian, and C. Z. Zhu, “Multiregional functional near-infrared spectroscopy reveals globally symmetrical and frequency-specific patterns of superficial interference,” Biomed. Opt. Express 6(8), 2786–2802 (2015). [CrossRef]   [PubMed]  

22. T. Wilcox, H. Bortfeld, R. Woods, E. Wruck, and D. A. Boas, “Hemodynamic response to featural changes in the occipital and inferior temporal cortex in infants: a preliminary methodological exploration,” Dev. Sci. 11(3), 361–370 (2008). [CrossRef]   [PubMed]  

23. T. Wilcox, J. A. Haslup, and D. A. Boas, “Dissociation of processing of featural and spatiotemporal information in the infant cortex,” Neuroimage 53(4), 1256–1263 (2010). [CrossRef]   [PubMed]  

24. Y. Zhang, D. H. Brooks, M. A. Franceschini, and D. A. Boas, “Eigenvector-based spatial filtering for reduction of physiological interference in diffuse optical imaging,” J. Biomed. Opt. 10(1), 011014 (2005). [CrossRef]   [PubMed]  

25. J. Virtanen, T. Noponen, and P. Meriläinen, “Comparison of principal and independent component analysis in removing extracerebral interference from near-infrared spectroscopy signals,” J. Biomed. Opt. 14(5), 054032 (2009). [CrossRef]   [PubMed]  

26. S. Kohno, I. Miyai, A. Seiyama, I. Oda, A. Ishikawa, S. Tsuneishi, T. Amita, and K. Shimizu, “Removal of the skin blood flow artifact in functional near-infrared spectroscopic imaging data through independent component analysis,” J. Biomed. Opt. 12(6), 062111 (2007). [CrossRef]   [PubMed]  

27. H. Zhang, Y. J. Zhang, C. M. Lu, S. Y. Ma, Y. F. Zang, and C. Z. Zhu, “Functional connectivity as revealed by independent component analysis of resting-state fNIRS measurements,” Neuroimage 51(3), 1150–1161 (2010). [CrossRef]   [PubMed]  

28. K. E. Jang, S. Tak, J. Jung, J. Jang, Y. Jeong, and J. C. Ye, “Wavelet minimum description length detrending for near-infrared spectroscopy,” J. Biomed. Opt. 14(3), 034004 (2009). [CrossRef]   [PubMed]  

29. B. Molavi and G. A. Dumont, “Wavelet based motion artifact removal for Functional Near Infrared Spectroscopy,” in International Conference of the IEEE Engineering in Medicine & Biology, 2010), 5. [CrossRef]  

30. A. M. Chiarelli, E. L. Maclin, M. Fabiani, and G. Gratton, “A kurtosis-based wavelet algorithm for motion artifact correction of fNIRS data,” Neuroimage 112, 128–137 (2015). [CrossRef]   [PubMed]  

31. R. J. Cooper, J. Selb, L. Gagnon, D. Phillip, H. W. Schytz, H. K. Iversen, M. Ashina, and D. A. Boas, “A systematic comparison of motion artifact correction techniques for functional near-infrared spectroscopy,” Front. Neurosci. 6, 147 (2012). [CrossRef]   [PubMed]  

32. S. Brigadoi, L. Ceccherini, S. Cutini, F. Scarpa, P. Scatturin, J. Selb, L. Gagnon, D. A. Boas, and R. J. Cooper, “Motion artifacts in functional near-infrared spectroscopy: a comparison of motion correction techniques applied to real cognitive data,” Neuroimage 85(Pt 1), 181–191 (2014). [CrossRef]   [PubMed]  

33. E. Kirilina, N. Yu, A. Jelzow, H. Wabnitz, A. M. Jacobs, and I. Tachtsidis, “Identifying and quantifying main components of physiological noise in functional near infrared spectroscopy on the prefrontal cortex,” Front. Hum. Neurosci. 7, 864 (2013). [PubMed]  

34. P. Pinti, D. Cardone, and A. Merla, “Simultaneous fNIRS and thermal infrared imaging during cognitive task reveal autonomic correlates of prefrontal cortex activity,” Sci. Rep. 5(1), 17471 (2015). [CrossRef]   [PubMed]  

35. X. Zhang, J. Yu, R. Zhao, W. Xu, H. Niu, Y. Zhang, N. Zuo, and T. Jiang, “Activation detection in functional near-infrared spectroscopy by wavelet coherence,” J. Biomed. Opt. 20(1), 016004 (2015). [CrossRef]   [PubMed]  

36. H. Xue, K. Lu, and N. Hao, “Cooperation makes two less-creative individuals turn into a highly-creative pair,” Neuroimage 172, 527–537 (2018). [CrossRef]   [PubMed]  

37. T. Nozawa, Y. Sasaki, K. Sakaki, R. Yokoyama, and R. Kawashima, “Interpersonal frontopolar neural synchronization in group communication: An exploration toward fNIRS hyperscanning of natural interactions,” Neuroimage 133, 484–497 (2016). [CrossRef]   [PubMed]  

38. J. Jiang, C. Chen, B. Dai, G. Shi, G. Ding, L. Liu, and C. Lu, “Leader emergence through interpersonal neural synchronization,” Proc. Natl. Acad. Sci. U.S.A. 112(14), 4274–4279 (2015). [CrossRef]   [PubMed]  

39. Y. Pan, X. Cheng, Z. Zhang, X. Li, and Y. Hu, “Cooperation in lovers: An fNIRS-based hyperscanning study,” Hum. Brain Mapp. 38(2), 831–841 (2017). [CrossRef]   [PubMed]  

40. C. Torrence and G. P. Compo, “A practical guide to wavelet analysis,” Bull. Am. Meteorol. Soc. 79(1), 61–78 (1998). [CrossRef]  

41. A. Grinsted, J. C. Moore, and S. Jevrejeva, “Application of the cross wavelet transform and wavelet coherence to geophysical time series,” Nonlinear Process. Geophys. 11(5/6), 561–566 (2004). [CrossRef]  

42. E. A. Lebedeva and E. B. Postnikov, “On alternative wavelet reconstruction formula: a case study of approximate wavelets,” R. Soc. Open Sci. 1(2), 140124 (2014). [CrossRef]   [PubMed]  

43. E. B. Postnikov, E. A. Lebedeva, and A. I. Lavrova, “Computational implementation of the inverse continuous wavelet transform without a requirement of the admissibility condition,” Appl. Math. Comput. 282, 128–136 (2016). [CrossRef]  

44. C. Torrence and P. J. Webster, “Interdecadal changes in the ENSO–Monsoon system,” J. Clim. 12(8), 2679–2690 (1999). [CrossRef]  

45. Y. J. Zhang, L. Duan, H. Zhang, B. B. Biswal, C. M. Lu, and C. Z. Zhu, “Determination of dominant frequency of resting-state brain interaction within one functional system,” PLoS One 7(12), e51584 (2012). [CrossRef]   [PubMed]  

46. Y. Zhao, R. N. Dai, X. Xiao, Z. Zhang, L. Duan, Z. Li, and C. Z. Zhu, “Independent component analysis-based source-level hyperlink analysis for two-person neuroscience studies,” J. Biomed. Opt. 22(2), 027004 (2017). [CrossRef]   [PubMed]  

47. C. M. Lu, Y. J. Zhang, B. B. Biswal, Y. F. Zang, D. L. Peng, and C. Z. Zhu, “Use of fNIRS to assess resting state functional connectivity,” J. Neurosci. Methods 186(2), 242–249 (2010). [CrossRef]   [PubMed]  

48. H. H. Jasper, “The 10–20 Electrode System of the International Federation,” Electroencephalogr. Clin. Neurophysiol. 10, 371–375 (1958).

49. A. K. Singh, M. Okamoto, H. Dan, V. Jurcak, and I. Dan, “Spatial registration of multichannel multi-subject fNIRS data to MNI space without MRI,” Neuroimage 27(4), 842–851 (2005). [CrossRef]   [PubMed]  

50. J. C. Ye, S. Tak, K. E. Jang, J. Jung, and J. Jang, “NIRS-SPM: statistical parametric mapping for near-infrared spectroscopy,” Neuroimage 44(2), 428–447 (2009). [CrossRef]   [PubMed]  

51. M. Hiraoka, M. Firbank, M. Essenpreis, M. Cope, S. R. Arridge, P. van der Zee, and D. T. Delpy, “A Monte Carlo investigation of optical pathlength in inhomogeneous tissue and its application to near-infrared spectroscopy,” Phys. Med. Biol. 38(12), 1859–1876 (1993). [CrossRef]   [PubMed]  

52. A. P. Holmes and K. J. Friston, “Generalisability, random effects and population inference,” Neuroimage 7, S754 (1998).

53. B. Biswal, F. Z. Yetkin, V. M. Haughton, and J. S. Hyde, “Functional connectivity in the motor cortex of resting human brain using echo-planar MRI,” Magn. Reson. Med. 34(4), 537–541 (1995). [CrossRef]   [PubMed]  

54. V. D. Calhoun, T. Adali, G. D. Pearlson, and J. J. Pekar, “A method for making group inferences from functional MRI data using independent component analysis,” Hum. Brain Mapp. 14(3), 140–151 (2001). [CrossRef]   [PubMed]  

55. S. M. LaConte, S. C. Ngan, and X. Hu, “Wavelet transform-based Wiener filtering of event-related fMRI data,” Magn. Reson. Med. 44(5), 746–757 (2000). [CrossRef]   [PubMed]  

56. L. Duan, Y. J. Zhang, and C. Z. Zhu, “Quantitative comparison of resting-state functional connectivity derived from fNIRS and fMRI: a simultaneous recording study,” Neuroimage 60(4), 2008–2018 (2012). [CrossRef]   [PubMed]  

57. Y. Tong and B. D. Frederick, “Time lag dependent multimodal processing of concurrent fMRI and near-infrared spectroscopy (NIRS) data suggests a global circulatory origin for low-frequency oscillation signals in human brain,” Neuroimage 53(2), 553–564 (2010). [CrossRef]   [PubMed]  

Cited By

Optica participates in Crossref's Cited-By Linking service. Citing articles from Optica Publishing Group journals and other participating publishers are listed here.

Alert me when this article is cited.


Figures (10)

Fig. 1
Fig. 1 The flow chart of the method. The chart illustrates the two-stage procedure for removing a single channel’s global physiological noise. The black section of the mask is the blocked portions. The procedure is repeated per channel to complete the entire noise removal processing.
Fig. 2
Fig. 2 The simulation experiment. (A) The simulated 66-channel fNIRS measurement (6 rows by 11 columns). All the channels were filled with simulated global physiological noise. The red frame indicates the ROI in which the simulated neural activity was added. A maximum of 6 s time-lag of the physiological noise was simulated, spreading from the bottom-middle channel to the outer channels as indicated by grey scale. (B) The simulated task-evoked neural signal (left) and the spontaneous resting-state neural signal (right) as the ground truth. (C) The noise contaminated signal (the blue line). The ground truth signal was also plotted (the black lines). Note the extremely low signal-to-noise ratio of the mixture signal.
Fig. 3
Fig. 3 The recovered neural signal. The red line indicates the recovered neural signal and the black line indicates the ground truth signal. Only for display purpose, a 5-point moving average of the recovered signal was applied to control the amplitude of the Gauss white noise.
Fig. 4
Fig. 4 The denoising results of different globality threshold k in the simulation experiment. The vertical coordinate was the p-value of the regression. The horizontal coordinate was the globality parameter k, represented by a ratio (the ratio means for the currently denoising channel, with how many channels of the other 65 channels a time-frequency component co-varies should it be labeled as contaminated in the present study). Note that when k was less than 20% (13/65, the number of the ROI channels divided by the number of all channels, except the currently denoising ROI channel), the recovered signal cannot significantly reproduce the ground truth signal (the dashed line indicates the significance level of 0.05).
Fig. 5
Fig. 5 The task activation maps and the RSFC maps. The left column and the right column respectively show the group-level task activation results and the RSFC results, derived from the corrected HbO data (the second row) and the uncorrected HbO data (the third row), with the globality threshold k = 50%. The first row shows the sensorimotor template.
Fig. 6
Fig. 6 The receiver operating characteristic (ROC) curve. The left and the right panel show the ROC curves of the task activation map and the RSFC map, respectively (globality threshold k = 50%). The red line indicates the ROC curve of the map derived from the corrected HbO data, and the blue line indicates the ROC curve of the map derived from the uncorrected HbO data.
Fig. 7
Fig. 7 The denoising performance of different globality threshold k on the real experiment. The vertical coordinate is the AUC index of the task activation map (the blue line) and the RSFC map (the red line) of the HbO data. The horizontal coordinate is the globality parameter k from 20% to 100%, stepped by 10%. The dashed lines show the AUC values derived from the uncorrected data.
Fig. 8
Fig. 8 The task activation maps and the RSFC maps of HbR data. The left column and the right column respectively show the group-level task activation results and the RSFC results, derived from the corrected HbR data (the second row) and the uncorrected HbR data (the third row), with the globality threshold k = 50%. The first row shows the sensorimotor template.
Fig. 9
Fig. 9 The receiver operating characteristic (ROC) curves of HbR data. The left and the right panel show the ROC curve of the task activation map and the RSFC map, respectively (globality threshold k = 50%). The red line indicates the ROC curve of the map derived from the corrected HbR data, and the blue line indicates the ROC curve of the map derived from the uncorrected HbR data.
Fig. 10
Fig. 10 The denoising performance of different globality threshold k on the real experiment (HbR data). The vertical coordinate is the AUC index of the task activation map (the blue line) and the RSFC map (the red line) of the HbR data. The horizontal coordinate is the globality parameter k from 20% to 100%, stepped by 10%. The dashed lines show the AUC values derived from the uncorrected data.

Equations (4)

Equations on this page are rendered with MathJax. Learn more.

ψ 0 ( η ) = π 1 / 4 e i ω 0 η e 1 2 η 2 ,
W x n ( n , s ) = δ t s n ' = 1 N x n ' ψ 0 [ ( n ' n ) δ t s ] ,
R 2 ( n , s ) = | S [ s 1 W x n , y n ( n , s ) ] | 2 S [ s 1 | W x n ( n , s ) | 2 ] S [ s 1 | W y n ( n , s ) | 2 ] ,
k > N R O I N T ,
Select as filters


Select Topics Cancel
© Copyright 2024 | Optica Publishing Group. All rights reserved, including rights for text and data mining and training of artificial technologies or similar technologies.