Abstract
Brain-computer interface (BCI) is one technology that allows a user to communicate with external devices through detecting brain activity. As a promising noninvasive technique, functional near-infrared spectroscopy (fNIRS) has recently earned increasing attention in BCI studies. However, in practice fNIRS measurements can suffer from significant physiological interference, for example, arising from cardiac contraction, breathing, and blood pressure fluctuations, thereby severely limiting the utility of the method. Here, we apply the multidistance fNIRS method, with short-distance and long-distance optode pairs, and we propose the combination of independent component analysis (ICA) and least squares (LS) with the fNIRS recordings to reduce the interference. The short-distance fNIRS measurement is treated as the virtual channel and the long-distance fNIRS measurement is treated as the measurement channel. Least squares is used to optimize the reconstruction value for brain activity signal. Monte Carlo simulations of photon propagation through a five-layered slab model of a human adult head were implemented to evaluate our methodology. The results demonstrate that the ICA method can separate the brain signal and interference; the further application of least squares can significantly recover haemodynamic signals contaminated by physiological interference from the fNIRS-evoked brain activity data.
1. Introduction
Near-infrared spectroscopy (NIRS), employing continuous wave (CW) instrumentation with several wavelengths and a variety of source-detector configurations, has attracted growing interest recently [1]. Changes in the concentration of oxyhaemoglobin and deoxyhaemoglobin (HHb), evoked by an appropriate stimulus, are calculated using the modified Lambert-Beer law (MLBL) [2, 3]. In real situations, the physiological events related to the cardiac cycle, breathing, and blood pressure pulsations often have low correlation with the functional response evoked by the task execution [4]. The measured signals are inevitably affected by such interference, which can arise from the superficial layers and the brain itself, including scalp, skull, cerebrospinal fluid (CSF), gray matter and white matter.
Near-infrared, spectroscopy uses stimuli to evoke physiological responses, and it has developed into a commonly used method called functional near-infrared spectroscopy (fNIRS). The concentration of HbO2 and HHb changes in brain activities can be measured simultaneously and this allows a user to interact with the outside world through the measurement of correlates of neural activity associated with mental processes. Brain-computer interfaces (BCIs) can be characterized in a number of ways, such as functional magnetic resonance imaging (fMRI), magnetoencephalography (MEG), positron emission tomography (PET), and electroencephalography (EEG) [5, 6]. Comparing with these methods, the advantages of fNIRS for the investigation of brain activity include portability, fewer physical restrictions, good temporal resolution, safety, and inexpensive instrumentation [7]. These characteristics make fNIRS being a promising method to measure the concentration changes of HbO2 and HHb, the two primary absorbing chromophores in the brain tissue, and accurately describes the cerebral haemodynamics. But the presence of physiological interference is one of the main problems in using fNIRS. These physiological interferences come from capillaries and vessels in superficial layer and internal layers in brain, served as “global interference” or “systemic physiological interference.”
In an attempt to decorrelate the physiological fluctuations from the evoked haemodynamic response, various groups have made extensive research on fNIRS data. Bandpass filtering, low-pass filtering, moving averaging, and Wiener filtering [8] have been developed to eliminate high-frequency instrument noise and low-frequency drift. The interference caused by cardiac oscillations can be effectively removed by these frequency-based algorithms. However, the specific physiological noise signals such as respiratory and blood pressure variation are left over since these fluctuations are difficult to be distinguished from the haemodynamic response to brain activity by frequency characteristics alone. A growing number of algorithms are still being developed for fNIRS noise reduction and signal improvement. Saager and Berger [9] utilized a dual-detector system and least squares method to remove top-layer-only fluctuations and validated the effect of the methodology by performing Monte Carlo simulations based on a two-layer turbid media model. Zhang et al. [6] proposed eigenvector-based methods to separate an activity-evoked response from systemic physiological interference in diffuse optical imaging data. They usually assumed that signal due to noise is broadly spatially distributed compared to signal due to neural activity and it is difficult to apply to real-time processing. Morren et al. [10] and Zhang et al. [11] adopted adaptive filtering to remove global interference. The difference between the two approaches is that the former used signals from an additional hardware as a reference and the latter used signals from a channel with a very short distance between the emitter and detector as a reference. A Kalman filtering model has also been used to analyze interference components [12, 13]. Abdelnour and Huppert proposed an adaptive general linear model based on a Kalman filtering algorithm for real-time assessment of brain function [14]. Despite the fact that the general linear model in their papers exhibited the potential for removal of physiological signals from cardiac, respiratory, and Mayer wave fluctuations, the authors expressly indicated that prior distributions must be assumed for both the process and observed noise in the Kalman filtering model.
Recently, independent component analysis (ICA) has been shown to be able to separate and identify interference arising from some cardiac events and breathing and has been widely used in the EEG and fMRI research community [15]. One important requirement in ICA is that all of sources are assumed to be mutually temporally independent and spatially stationary. For fNIRS recordings with a reference signal, the independence between the reference signal and sources contributing to the deep brain tissue should at least be approximately true based on the previous study. The requirement of spatial stationarity of sources is a limitation of ICA, but for fNIRS is not generally a problem if one can avoid the pulsatile movement of the brain and probe. To the best of our knowledge, ICA has not been applied to common multidistance fNIRS measurement recordings except for in a recent paper dealing with spatial eigenfiltering algorithm.
In our study presented here, we propose a new method of isolating useful information about hemodynamics in the cerebral layer. We used a multidistance measurement method and a theoretical analysis of global interference reduction based on ICA and the least squares (LS) criterion. The short-distance fNIRS measurement is treated as the reference channel comprising superficial haemodynamic changes induced by physiological fluctuations and the long-distance fNIRS measurement is treated as the measurement channel containing both the functional haemodynamic response and global interference. We aim to remove global interference that is correlated with superficial haemodynamic fluctuations, evoked by cardiac contractions, breathing, blood pressure, and so forth. ICA is a powerful method to separate the mixed multichannel recordings and recover them into their constituent latent sources. By decomposing the long-distance measurement with the ICA algorithm, we separated the signal into different components based on the virtual channel. The least squares criterion was then used to adjust the corresponding weighting coefficients and the direction to estimate the real signal with the obtained component. Monte Carlo simulations of a five-layer human model were used to investigate the performance of the ICA-LS for removing global interference in brain activity measurement.
2. Theory
2.1. Multidistance Probe Configuration and Modified Lambert-Beer Law
The use of near-infrared (NIR) light (700~1000 nm) for functional measurement in biological tissues depends, firstly, on the fact that biological tissues are relatively transparent at these wavelengths and, secondly, on the oxygen dependence of haemoglobin and cytochrome aa3 optical absorption properties [16]. A multilayer tissue model, that consists of scalp, skull, cerebrospinal fluid (CSF), gray matter, and white matter, was used in our study to represent the human head. A multidistance NIRS probe arrangement together with a schematic of five-layer head model is illustrated in Figure 1. One light source with two wavelengths was located on the surface of the head and two detectors were positioned at different distances from the source to collect light emerging from the head after diffuse reflectance. The use of particular source-detector separations allows us to distinguish between processes occurring at different tissue depths, as indicated by the banana-shaped regions that encompass statistically determined photon paths [17]. We use the short distance with source-detector 1 to probe the superficial tissue layer and long distance with source-detector 2 to probe the deeper tissue layer.
When we obtained the change of optical density with the source and the detectors, the concentration changes of HbO2 and HHb at a given time could be estimated using the modified Lambert-Beer law (MLBL) as follows:
2.2. The Standard ICA Model with Virtual Channels
ICA is a multivariate statistical technique which is used to estimate a set of signals of which only a mixture is available. The sources and the mixing process are both unknown, and the sources are estimated on the assumption that they are statistically independent from one another [18].
For a one-dimensional signal with additive noise, the mathematical model can be represented by the following equation: wherein is the observed signal, is the source, and is the mixture of multiple noise components. In the case of possessing two or more separate noise components, can be expressed as where is the number of the noise and is the weight of the th noise. Introducing noise components as the virtual noise channels, then we have the following expression: According to (3) and (4), (2) can be expressed by the following standard ICA model: From (3), the desired signal, , can be obtained by introducing the virtual channels. In our paper, we use the matrix below to apply the ICA method with virtual channels in which is a one-dimensional vector standing for the observed signal and is a one-dimensional vector stand for the mixture of multiple noise components: But problems still exist since the independent component analysis needs to do pretreatment on the data to make sure that they are statistically independent from one another. The range of matrix has changed after pre-treatment, such as centralization and whitening, so we cannot get the accurate result directly by ICA method. Here, we applied least squares method to do the fitting of desired signal and the observed signal so that the range of the data can be confirmed. These are the elements in and in which the number stands for the length of the vector: The functional relationship between and can be described in the form below: In order to obtain the optimal values of and , we use the mean-square error [19]: And is the sum of : If we want to get the minimum of , the partial derivation on and must be zero: By consolidating (9) and (10), we get the functional relationship between and : The coefficients and can be calculated from the above two equations. By using least squares method, the range of the signals can be recovered and physiological interference can be reduced.
3. Methods
3.1. Monte Carlo Simulations
A layered model of the adult human head was used to investigate the light propagation and attenuation. The model was composed of five layers including scalp, skull, cerebrospinal fluid (CSF), gray matter, and white matter [20]. Our study used two source-detector pairs and two wavelengths of light (see Figure 1). The generation of the NIRS time series was conducted with Monte Carlo simulations [21] based on our previous analysis of functional hemodynamic changes [16]. By using two source-detector pairs and two wavelengths, we can get the mixed signals and the virtual channels. Since the mixed signal is a one-dimensional matrix, the virtual channels are necessary for the use of independent component analysis.
The Monte Carlo code used here is an extension of the general multilayer, three-dimensional, weighted photon Monte Carlo codes developed by Wang et al. [21]. Scattering anisotropies were assumed to be 0.9 and Fresnel reflection at the tissue-air boundary was also considered. We assumed the same refractive index of for all layers [22]. The standard parameters used in the simulation are given in Table 1. Thickness, transport scattering coefficient , and HbO2 and HHb baseline concentrations were taken from published data [7, 20, 23]. The baseline concentrations of HbO2 and HHb assumed that the oxygen saturation in the head was 70%. Absorption coefficients were calculated with the HbO2 and HHb baseline concentrations and the molar extinction coefficients. The molar extinction coefficients at 750 nm and 830 nm were obtained from the literature [24].
3.2. Haemodynamic Response Function and Physiological Interference Simulation
To further quantify the utility of ICA-LS in the removal of physiological noise, and therefore the improved recovery of a functional response, we introduced a series of simulated haemodynamic functional responses function. The functional haemodynamic responses in gray matter were defined as the convolution of the stimulation , for rest period, and 1 for stimulation and a prototypical haemodynamic impulse [11, 20, 25]: The parameters and in the gamma variate function were set at 8.6 and 0.56, respectively, which corresponds to recent findings [20, 25]. The evoked haemodynamic response was the convolution of and : where is a simple scaling factor in our model.
In order to make the simulation as realistic as possible, the haemodynamic changes were simulated as a combination of the functional haemodynamic responses and the physiological interference. The physiological interference is generated by a combination of cardiac fluctuation , respiratory fluctuations , low-frequency oscillation , very-low-frequency oscillation , and independent fluctuation induced by the temperature changes and the sweat on the skin. The haemodynamic changes in each layer are denoted as follows: where , , , and represent the concentrations of HbO2 and HHb in each layer as a function of time, with the superscripts 1 to 5 indicating the layer index for scalp, skull, CSF, gray and white matters, respectively. and HHbbase represent average or baseline concentrations. The cardiac fluctuation , respiratory fluctuation , low-frequency oscillation , and the very-low-frequency oscillation were all simulated as the summation of sinusoidal wave and white Gaussian noise. The coefficients , , , , and with layer index as superscript and HHb or HbO2 as subscript are the haemodynamic variation amplitude control parameters, which can be found in Table 2. The independent interference is generated by biased and low-pass filtered Gaussian white noise. The independent interference in the scalp layer and the different weight on cardiac, respiratory, and low-frequency and very-low-frequency fluctuations in each layer were to simulate a certain amount of uncorrelated changes in the superficial layers compared with deep layers. The parameters used are based on Scholkmann et al. [26] and Zhang et al’s. work [11].
The simulated haemodynamic changes were used to calculate the optical measurement by Monte Carlo method. We derived the simulated optical measurements by launching 108 photon packets and the running time of Monte Carlo simulation is about 10 hours on the desktop (Intel Core i5-2320 CPU). The sampling rate was set to 10 Hz and the whole time series for the changes of optical density were acquired under the assumption that the scattering properties of the head do not vary with time. The experiment is designed as a 5-epoch block and each individual epoch consisted of a series of 400 points, 200 points of rest, and 200 points of stimulation. The short source-detector spacing was set to 5 mm, which only considers penetration of light into the extra cerebral tissue. The long source-detector spacing was set to 45 mm, which is long enough for penetration into the cerebral cortex. Therefore, the short-distance optode pair was used as the virtual channel and the long-distance optode pair was used as the measurement channel in the ICA model. Using (3) and the derived formulae, the evoked brain signals could be separated.
4. Results
The simulated optical measurements here were obtained by employing S-D1 with a separation of 5 mm and S-D2 with a separation of 45 mm. The concentration changes of oxyhaemoglobin, , and deoxyhaemoglobin, were then derived with the MLBL. The concentration changes of oxyhaemoglobin, , and the concentration changes of deoxyhaemoglobin, , were separated using the ICA algorithm and the results are presented in Figures 2 and 3, respectively.
(a)
(b)
(c)
(d)
(a)
(b)
(c)
(d)
ICA is a signal discrimination method that extracts independent components from multiple signals without knowledge of the obtained signal by utilizing the statistical independence of the source components [27].
In Figure 2, where the result is presented, the virtual channel dataset is the calculated from S-D1 with 5 mm source-detector separation (Figure 2(a)), the measurement channel dataset is the calculated from S-D2 with 45 mm source-detector separation (Figure 2(b)). In theory, the measurement channel dataset should contain the evoked hemodynamic response; however, the raw time series present significant interference. In fact, its similarity in the raw time series to the reference dataset Figures 2(a) and 2(b) indicates that it is dominated by global interferences. Using the data obtained with S-D1 and S-D2 source-detector pairs, ICA was used to extract independent components at each S-D2 distance. Figures 2(c) and 2(d) plot the time course of independent components extracted with ICA.
In Figure 3, we show the equivalent results of Figure 2 but for [HHb]. During the implementation of independent component analysis for [HHb], we chose the same algorithm as that used for . The difference between Figures 2 and 3 is that the residual interference for [HHb] is less obvious than that for . This result could be explained on the basis that the interference component in the [HHb] signal is relatively smaller than that in the signal, and this is in line with known facts. Both and [HHb] results shown in Figures 2 and 3 demonstrated that the ICA algorithm exhibited a good performance in physiological interference separation.
The results of hemodynamic changes and the associated processed results are shown in Figures 2 and 3. We found that the hemodynamic changes (Figures 2(a), 2(b), 3(a), and 3(b)) are contaminated by the physiological interference. Using the ICA method, physiological interference could be separated from the original data, as shown in Figures 2(c) and 3(c). The brain signals extracted with ICA are shown in Figure 2(d) and Figure 3(d). It is obvious that the physiological interference has been suppressed.
However, the magnitudes of and [HHb] are underestimated because the activated volume is smaller than the sampling volume. This is generally referred to as the partial volume effect (PVE) [28]. By means of Monte Carlo simulations, the ratio of the optical path length in the activated volume to the optical path length in the sampling volume can be achieved to compensate the PVE effect. Here, we compensate the PVE so as to compare the filter output quantitatively with the real evoked haemodynamic response.
Concentration changes of HbO2 during the simulated measurement are shown in Figure 4. Thereinto, Figure 4(a) shows the time series results extracted with ICA from a 45 mm source-detector distance that were calculated with the MLBL and Figure 4(b) shows the block average result. Similarly, Figures 4(c) and 4(d) show the [HHb] time series extracted with ICA from a 45 mm source-detector distance and block average results, again calculated with the MLBL. After signal processing with the ICA-LS, as seen in Figures 4(a) and 4(c), the evoked haemodynamic changes are clear and most of global interference has been removed. After we have compensated for the PVE effect in the recovered results, we can compare the recovered results with the real evoked haemodynamic changes in the gray matter. This comparison can be seen in Figures 4(a) and 4(c). In these two figures, the solid line denotes the recovered results with PVE compensation and the dashed line denotes the true evoked haemodynamic changes used in the simulation. Although some fluctuations still remain, the recovered processed with the ICA-LS algorithm provides an obvious evoked response. By calculating ensemble average for the whole time series, the results demonstrate that the proposed methodology could remove approximately 95% of global interference.
(a)
(b)
(c)
(d)
5. Discussion
Several methods have been explored in the literature to attempt to remove interference from fNIRS recordings. If a reference signal for physiological interference is available, such as the pulse oximeter, electrocardiogram (ECG), it can be subtracted from the fNIRS data after scaling it by an appropriate factor determined by regression either in the time domain or frequency domain. However, the recorded reference signal tends to be contaminated as well; regressing it out would thus attenuate the fNIRS value, which is undesirable. Moreover, reference signals are not available for other interferences such as the low-frequency oscillation. Another approach for interference elimination is to use digital filters; however, the frequency spectrum of most interferences overlaps with that of the brain activity signal; hence, the interferences can only be partially removed by this method.
More recently, principal component analysis (PCA) [6] has been used to separate the fNIRS signal into uncorrelated components. Physiological interference can be identified and removed, especially if their amplitude is high. However, PCA restricts the extracted components to have orthogonal spatial topographies, which is a difficult assumption. ICA uses a stronger assumption of statistical independence, which is more appropriate for fNIRS recordings. However, in the case of fNIRS heavily contaminated by artifacts, some components might be a mixture of multiple strong sources. The ICA is thus still able to remove several interference components from the fNIRS recordings.
We have used a Monte Carlo method of a five-layered model of the human adult head to simulate the brain activity experiments. Monte Carlo modeling has been extensively utilized in the biomedical fields. Different approaches have been used to develop the Monte Carlo methods such as white Monte Carlo [29], hybrid Monte Carlo-diffusion method [30], and the graphics processing units (GPU) [31] technique. A further matter of importance is that the GPU method has been widely adopted to accelerate Monte Carlo simulation because GPU can be used for parallel computing. Recently, many GPU-based Monte Carlo computational tools such as online-object-oriented and peer-to-peer Monte Carlo simulations have been deeply developed by Doronin and Meglinski [32, 33]. As yet we have not carried out GPU-based Monte Carlo study. This is the subject of our ongoing research and will be used for ICA-LS algorithm in due course.
The present paper introduces a new approach to physiological interference cancellation from fNIRS recordings. The inclusion of the reference channel has succeeded in ensuring that one component is primarily determined by physiological interference. We have used Monte Carlo simulations to assess the use of independent component analysis in global interference removal from fNIRS brain activity data. Least squares algorithm was further employed to determine the amplitude of the independent components. Our results have shown that the combination of independent component analysis and least squares can reduce global interference within the measurements of and [HHb].
6. Conclusions
Studies of independent component analysis of fNIRS data with a virtual channel have demonstrated that the physiological interference in the fNIRS signals can be substantially suppressed. These results were evaluated with multidistance fNIRS measurement using Monte Carlo simulation. Least squares has been applied to deal with the separating components calculated with ICA and the results show that brain activity response can be separated in fNIRS signals since the method is effective to remove global interference induced not only by heartbeat and respiration but also by low-frequency oscillation and very-low-frequency oscillation, and other correlated interferences between superficial and deep layer. The advantage of ICA-LS in multidistance fNIRS measurement compared with other possible methods also arises from its convenient implementation; it neither requires an auxiliary measurement instrument nor the dependence on a priori knowledge of the global interference frequency. The experiment also indicated that, by introducing the virtual noise channel properly, the brain activity response can be extracted even in the strong noise background. Thus, the methodology has clear potential for use in fNIRS brain-computer interface measurement.
Acknowledgments
The authors are grateful for the support from the National Science Foundation of China (Grant no. 61201017), the Natural Science Foundation of Heilongjiang Province (Grant no. QC2011C097), China Postdoctoral Science Foundation (Grant no. 2013M531027), Heilongjiang Postdoctoral Fund (Grant no. LBH-Z12093), the Fundamental Research Funds for the Central Universities (Grant no. HIT.NSRIF.2013010), and OBHL, UK. They also thank colleagues in China, Italy, UK, and Japan for their helpful comments on this work.