|Home | About | Journals | Submit | Contact Us | Français|
The simultaneous acquisition of electroencephalogram (EEG) and functional MRI (fMRI) signals is potentially advantageous because of the superior resolution that is achieved in both the temporal and spatial domains, respectively. However, ballistocardiographic artifacts along with the ocular artifacts are a major obstacle for the detection of the EEG signatures of interest. Since the sources corresponding to these artifacts are independent from those producing the EEG signatures, we applied the Infomax-based independent component analysis (ICA) technique to separate the EEG signatures from the artifacts. The isolated EEG signatures were further utilized to model the canonical hemodynamic response functions (HRFs). Subsequently, the brain areas from which these EEG signatures originated were identified as locales of activation patterns from the analysis of fMRI data. Upon the identification and subsequent evaluation of brain areas generating interictal epileptic discharge (IED) spikes from an epileptic subject, the presented method was successfully applied to detect the theta- and alpha-rhythms that are sleep onset related EEG signatures along with the subsequent neural circuitries from a sleep deprived volunteer. These results suggest that the ICA technique may be useful for the preprocessing of simultaneous EEG-fMRI acquisitions, especially when a reference paradigm is unavailable.
Functional neuroimaging studies have been conducted to examine the neurophysiological processes that occur during different stages of sleep. EEG has been extensively used to investigate staging of sleep (Flexer et al., 2005) or the detection of abnormal sleep patterns such as narcolepsy (Liu et al., 2008). However, high spatial-resolution, region-specific brain activity cannot be obtained from EEG alone. Recently, positron emission tomography (PET) has been used to detect the changes in regional cerebral blood flow (rCBF) during sleep onset (Kotajima et al. 2005) to detect abnormalities in glucose metabolism in obstructive sleep apnea patients (Antczak et al., 2007).
Functional MRI has gained popularity within the neuroimaging community as a major research tool due to its ability to repeatedly measure and map cortical activity at high-resolution and in a non-invasive manner. Although fMRI can provide information on region-specific brain activity in greater spatial detail, the temporal resolution of the acquisition is rather slow (on the order of seconds) compared to that of EEG and is further compounded by the slow rate of the blood-oxygenation-level-dependent (BOLD) hemodynamic response (on the order of 3–4 sec). In addition, current methods for the characterization of sleep onset and staging are heavily based on the EEG signature. In an effort to enable fMRI to detect sleep-related brain function, a multi-modal approach combining EEG and fMRI is considered (Garreffa et al., 2003; Singh et al., 2003; Diehl et al., 2003; Schabus et al., 2007).
Several MR-compatible EEG systems are commercially available with commonly-shared design features that minimize mutual electronic interference/interaction. The design features include the acquisition of the EEG signals and concurrent conversion (via multiplexing signals) to optical signals after pre-amplification. The optical signals are ported outside of the RF-shield around the MRI room via a fiber optic link and are de-multiplexed for further processing (Ives et al., 1993; Allen et al., 2000; Bonmassar et al., 2002).
EEG acquisition does not typically deteriorate MRI data when proper RF-shielding/filtering is employed (Krakow et al., 2000). However, EEG data is susceptible to contamination by a series of artifacts during MRI scan, which must be removed or compensated for in order to correctly determine the underlying brain activities. The most evident cause of the artifacts is the periodic switching of the rapid magnetic gradients that are necessary during MRI data acquisition and contaminate the acquired EEG signals with signal fluctuation, up to a few millivolts in amplitude (Allen et al., 2000). This type of artifact, due to its predictable and reproducible pattern, can be easily removed through real-time subtraction of time-averaged artifacts (the so called averaged artifact subtraction or AAS; Allen et al., 2000).
More elusive forms of artifacts could be traced from biological origins and are more difficult to correct. They are known as ballistocardiogram artifacts (BCA) and ocular artifacts (Bonmassar et al., 2002; Briselli et al., 2006; Mantini et al., 2007; Debener et al., 2007). BCA, arising from changes in electrical potential associated with subtle scalp motion or cranial/brain perfusion in the middle of a strong magnetic field, appears as a low frequency and periodic noise synchronized with cardiac pulsation (Ives et al., 1993; Allen et al., 1998; Bonmassar et al., 2002). BCA can contaminate the ERP signal or mimic low delta/theta signal bandwidth and is thus the most significant artifact to consider (Mantini et al., 2007). Due to its unpredictable, quasi-periodic nature and its resemblance in some resting slow waves in EEG, BCA remains even after removing the gradient-induced signal fluctuation via AAS. The ocular artifact, which is associated with mechanical eye movement or blinking, is also unpredictable and overpowers the EEG signal, especially from the frontal areas of the brain (He et al., 2004).
Various filtering methods, primarily based on subtraction or adaptive filtering cued by electrocardiogram (ECG) and electrooculogram (EOG) activity, are used to filter out the artifacts. For example, time-locked averaging of an EEG signal during the interval of the ECG signal is used to subtract out the BCA component (Allen et al. 1998). Kim and colleagues (2004) proposed the use of recursive wavelet-based adaptive filtering for the removal of BCA. Bonmassar et al. (2002) utilized an adaptive filtering scheme derived from a Kalman filter in which an additional motion sensor was required to provide a reference signal. All of these methods, however, have been susceptible to deficiencies based on the choice of the modeled filters or reference signal and have lacked consistent reliability for successful artifact removal.
Recently, independent component analysis (ICA) has been applied to detect and successfully remove BCA and ocular artifacts (Srivastava et al., 2005; Briselli et al., 2006; Nakamura et al., 2006; Mantini et al., 2007). ICA is a blind signal-source separation (BSS) technique, which allows the separation of an independent signal from an unknown mixture of signals (Bell & Sejnowski, 1995; Hyvärinen, 2001). The unknown nature of the EEG-artifacts that appears during the fMRI acquisition renders ICA an attractive and effective method for the removal of these artifacts with the potential capability for real-time data processing (Garreffa et al., 2003; Mantini et al., 2007). In comparison to conventional approaches, ICA has consistently shown more promising results in the exploitation of the particularly meaningful signatures of EEG data such as event-related potentials (ERPs), interictal epileptic discharges (IEDs) and frequency-dependent (e.g. delta, theta, alpha, beta, and gamma) rhythms (Srivastava et al., 2005; Briselli et al., 2006; Nakamura et al., 2006; Mantini et al., 2007; Debener et al., 2007).
In this study, we were motivated to examine the feasibility of using ICA to detect the transitions in EEG frequency bands during sleep onset and to characterize the corresponding brain activities during simultaneous EEG-fMRI data acquisition. Since one of the characteristics of sleep onset is an alternating appearance of theta and alpha rhythms in the time domain (De Gennaro et al., 2001; Corsi-Cabrera et al., 2000), these features may be separated from the noisy EEG signals by applying the ICA algorithm. In order to test the robustness of the developed method, we also applied the same software tool to detect interictal epileptic discharges (IEDs) from an epilepsy patient and visualized the IED loci.
The study was conducted according to the rules set forth by the local institutional review committee on the study of human subjects. For the recording of the EEG signals, both an MR-compatible BrainAmp and BrainCap (BrainProducts GmbH, Germany) were used. The BrainCap was an electrode cap containing either 32 or 64 EEG channels including three additional electrodes that could measure two ECG signals and one EOG signal. The reference electrode in the cap was located at the mid-point between Cz and Fz. Each channel was recorded with a sampling rate of 5×103 samples/sec (Hz). Gradient artifacts were eliminated using an AAS method (Allen et al., 2000) implemented in Analyzer 1.05 (BrainProducts, Munich, Germany). The gradient artifact-corrected EEG signals were filtered using a phase-shift-free Butterworth filter with a bandpass range from 0.5 to 55 Hz. The resulting signals were then down-sampled by a factor of 10 (i.e. 500 Hz) to reduce the amount of data, thereby reducing the processing time. Note that the down-sampled multi-channel EEG signals with attenuated gradient artifacts are still contaminated by physiological noise such as ocular artifacts and BCA (Martini et al., 2007; Debener et al., 2007; Briselli et al., 2006; Bonmassar et al., 2002).
A 3.0-Tesla MR system (Magnum, Hitachi Medical, Japan) was used for fMRI data acquisition. For the detection of BOLD signal associated with neural activity, a standard gradient echo planar imaging (EPI) sequence was used to image a volume covering the whole brain (25 axial slices with 4mm thickness per slice with no gap; TR/TE=3000/40ms; field-of-view=22×22 cm2; 64×64 voxels within a plane; flip angle or FA=90°).
We tested the feasibility of ICA in revealing hidden signatures/features from the EEG data obtained during fMRI acquisitions (i.e. simultaneous EEG-fMRI). The hidden signatures/features that were revealed from the EEG data were further utilized to examine the corresponding areas of activation from the fMRI data that were simultaneously obtained. To implement the method, the MATLAB (v7.0, R14, Mathworks, Natick, MA) computational environment was used.
A flow diagram of our analysis scheme is shown in Fig. 1. Firstly, the raw EEG data were processed using AAS and bandpass filtering implemented in Analyzer 1.05 to significantly reduce the gradient-related artifacts and to boost the signals within the frequency range-of-interest (i.e. from the delta to gamma rhythms), respectively. The resulting signals were down-sampled by a factor of 10 (i.e. from 5kHz to 500Hz) to reduce the computational load while maintaining the meaningful characteristics of the EEG signatures. Then, an ICA method was applied in order to blindly separate the source independent components (ICs) containing meaningful features from the noise components including ocular artifacts and BCA. An Infomax algorithm of ICA with an extension of a natural gradient (Bell & Sejnowski, 1995; Amari, 1998) was applied to this processing method by adopting a MATLAB-based toolbox (EEGLAB; http://sccn.ucsd.edu). The default parameters were used and the iterative learning scheme was finished in about one hour using a laptop computer (maximum 500 iterations; Intel Pentium M 2.16GHz; 2GB main memory).
After extracting the ICs, meaningful features of the EEG signals were identified based on their frequency power spectra obtained from a short-time Fourier transform (Oppenheim et al., 1999). These processing stages were implemented as MATLAB-based software modules (1) to examine the time course of the extracted ICs of EEG data to detect events-of-interest, (2) to examine the frequency spectrum of the channel-specific time courses of ICs, and (3) to analyze the detected events within the time-locked interval. Both the onset timings and durations of the events from the identified features were further utilized to model the canonical hemodynamic response functions (HRFs) implemented in SPM2 (www.fil.ion.ucl.ac.uk). The modeled HRFs were employed as regressors of the General Linear Model (GLM) to identify the brain areas (from the fMRI data) that are covariant with the EEG features. Then, a T-contrast map, corresponding to the identified EEG events as compared to the baseline condition, was obtained using SPM2 in which the t-score of each voxel represents the statistical value of activation compared to the null hypothesis (Logan et al., 2007).
Prior to the examination of neural substrates associated with sleep onset, we applied the data analysis scheme to EEG-fMRI data obtained from an epileptic patient (26 year old male). The patient was diagnosed with focal epilepsy originating from within the posterior temporo-occipital areas of both hemispheres (confirmed later by surgical intervention). A total of 10 minutes of EEG data were recorded simultaneously with fMRI data (200 EPI volumes using the same imaging parameters previously described). The first 5 volumes (15 sec) of the fMRI data were excluded from further analysis to allow T1–signal equilibration. The acquired EEG data, along with the preprocessed fMRI data, were then processed following the steps in Fig. 1. After applying the Infomax ICA algorithm, the separated IC showing the events of the IED spikes was identified based on visual inspection of both the temporal and spectral characteristics. Onset timings and durations of the events were also identified.
For the preprocessing of the fMRI data, the timing differences between slices of each EPI volume were compensated for and were subsequently realigned to the first EPI volume for movement correction. The realigned volumes in an individual image space (voxel size=3.4×3.4×4mm3) were normalized to the standardized anatomical MNI (Montreal Neurological Institute) space (voxel size=2×2×2mm3). Finally, an 8mm isotropic full-width-at-half-maximum (FWHM) 3-D Gaussian kernel was applied to the normalized data for spatial smoothing in order to reduce the spatial noise generated during the realignment and normalization processes. These preprocessing steps on the fMRI data were conducted using the SPM2 toolbox implemented based on MATLAB.
A sleep-deprived subject (37 year old male; awake for more than 32 hours) underwent the simultaneous EEG-fMRI experiment. The subject was asked not to consume any caffeine for at least 12 hours prior to the scan. The subject was allowed to sleep during the EEG-fMRI experiment. A total of five minutes of EEG data were recorded simultaneously with the fMRI data. Upon the separation of the independent source signals from the measured EEG data through the adoption of the Infomax ICA algorithm, the ICs containing dominant power within the frequency range of 4~8 Hz (theta-rhythms) and 8~12 Hz (alpha-rhythms) were manually selected. The onset timing and duration corresponding to the events of these rhythms were subsequently identified.
From a total of 100 fMRI volumes, the first 5 volumes were also excluded from further analysis to allow the T1 effect to be stabilized. Again, the acquired fMRI data were preprocessed using SPM2 (in sequential order: slice timing correction, realignment, normalization to the MNI space, and smoothing using an 8mm isotropic FWHM Gaussian kernel). Following the procedure shown in Fig. 1, the activation patterns were generated from the GLM method by employing the canonically modeled HRFs based on the identified EEG events as regressors. Again, a T-contrast map was calculated from the F-contrast maps of the neuronal activations correlated with the events of the identified theta- and alpha-rhythms.
Figure 2 illustrates a 10 sec time frame (which contained one of the IED events) during the acquisition of the EEG data. As shown in Fig. 2a, the AAS algorithm effectively reduced the gradient-related artifacts of EEG data. The resulting separated ICs from the subsequent analysis of this data using the Infomax ICA algorithm are shown in Fig. 2b. Note that the ICs corresponding to ballistocardiogram (red), ocular (green), and residual gradient (blue) artifacts were identified. In this time frame, one of the events of IEDs was detected in the component of IC19 at 387 sec for a duration of 1 sec (rectangular box). Three IED events were also identified at 20, 372, and 414 sec from IC19 in approximately the same duration, which indicates that the Infomax ICA algorithm successfully separated the IED-related source signal from the multi-channel EEG data. Following the reconstruction of the EEG signals from the separated ICs, excluding the artifact-related components (Srivastava et al., 2005; Mantini et al., 2007), the artifact-removed EEG signals were obtained and are shown in Fig. 2c. Note that the IEDs are more dominant in the electrodes located within occipital (O1 & O2), parietal (P3, P4, & P7) and temporal (T7 & T8) areas than those within (pre-) frontal areas (Fp1, Fp2, F3, F4, & Fz).
Four F-contrast maps corresponding to the four identified IED spikes were calculated from the simultaneously obtained fMRI data. Figure 3 shows the activation patterns of the resulting T-contrast map that was subsequently thresholded at p<10−3 and overlaid on the anatomical image. The labels of the substantial activation loci are summarized in the accompanying table. Dominant bilateral activations were clearly identified within the temporal and occipital cortices where the focal epileptic source of the participating patient was confirmed later via corrective surgery.
After confirming the effectiveness of the ICA method on the analysis of IED-related neural activity, the characterization of the EEG data during the simultaneous acquisition of fMRI was obtained and is shown in Fig. 4. Only 11 channels, including one EOG and two ECGs, are shown for simplicity. The significant EPI gradient artifacts within the original EEG data (Fig. 4a) were drastically reduced after the application of AAS and bandpass filtering (Fig. 4b). Further processing, using the Infomax ICA algorithm, successfully separated the ballistocardiogram-related artifacts as well as the hidden EEG signatures (as shown in Fig 4c; IC01~IC04 & IC07 in red: BCA-related; IC05 in green: theta-rhythm related; IC08 in blue: alpha-rhythm related). From the time-frequency plots represented as spectrograms shown in Fig. 4d, we can observe that IC05 and IC08 were dominated by the signal components within the theta- (4~8Hz) and alpha- (8~12Hz) related frequency bands (arrows: exemplary events of the corresponding rhythms). The identified theta and alpha rhythms alternately appear in the time domain and are generally recognized as the EEG features during the lightest stage 1 sleep and relaxed wakefulness or drowsiness, respectively (De Gennaro et al., 2001; Corsi-Cabrera et al., 2000).
Subsequently, the fMRI data were analyzed from the GLM method by using six canonical HRFs as regressors for the theta-rhythm related events and five canonical HRFs as regressors for the alpha-rhythm related events. The T-contrast maps obtained from the resulting six (theta-related) and five (alpha-related) F-contrast maps are shown in Fig. 5 (p<0.01; theta- and alpha-rhythm related activations in green and blue, respectively) along with labels of the activation loci at the bottom. Substantial activations within the putamen, middle frontal, and parietal areas were identified for the corresponding theta-rhythm. Activations in the thalamus, parahippocampus, and precuneus were found to be associated with the alpha-rhythm.
In this study, we presented the feasibility of the application of the ICA method toward the analysis of simultaneously obtained EEG and fMRI data. The ballistocardiogram artifacts in EEG signals due to slight movements of the EEG electrodes inside a strong magnetic field were blindly separated by deploying the Infomax ICA algorithm. Other components of noise such as the ocular and residual gradient artifacts were also segregated in the remaining output components. From the evaluation of the method as applied to the epileptic subject, the onset and duration of IEDs were clearly recognized from the artifact-removed EEG signals compared to the original unprocessed data. In addition, the areas of lesion were successfully identified from the activation patterns of the fMRI data sets that were analyzed from the GLM method by adopting the regressors modeled from the identified EEG signatures. Similarly, from the application to sleep onset detection, the EEG rhythms related to the early stage of sleep (i.e. theta- & alpha-rhythms) were uccessfully separated from the noisy EEG data recorded during the fMRI acquisition.
Based on the activation areas obtained from the simultaneously obtained fMRI data, our ICA-based technique appears to provide reasonable regions-of-interest (ROIs) corresponding to the neural circuitries related with these EEG rhythms. The strong activities corresponding to the theta rhythm in the frontal cortices were consistent with previous studies, whereby the theta rhythm is appreciated as an EEG signature in frontal cortices during sleep following extended wakefulness (De Gennaro et al., 2007; Strijkstra et al., 2003; Cajochen et al., 1999). In addition, the occipital, frontal, parietal, and thalamic activations corresponding to the alpha rhythm were in consonance with previous investigations, whereby the alpha frequencies are known to be strongly related with the occipital ‘visual’ areas or frontal and parietal ‘attentional’ areas (Hughes & Crunelli, 2005; Cantero et al., 2002; Laufs et al., 2006).
Recently in the sleep research area, there are growing interests in utilizing simultaneous EEG-fMRI data for applications such as the analysis of sleep spindles (Schabus et al., 2007; Lund et al., 2007; Laufs, 2008). For example, Schabus and colleagues (2007) have detected sleep spindles from the EEG signals during the non-rapid eye movement (NREM) sleep period. The subsequent neural circuitries involving the generation of the sleep-spindles were then identified from the simultaneously acquired fMRI data. Since the sleep spindles often show quite dominant and distinct EEG signatures, such as the spike bursts within the 11~15 Hz range lasting 0.5~3 sec, relatively simple processing tools such as bandpass filtering were able to detect these signal patterns (Schabus et al., 2007). On the other hand, as shown in Fig. 4b, the transition of standing waves related with sleep onset (i.e. theta- & alpha-rhythms) may not have dominant signal features within the recorded EEG signals. We believe that this study is the first to both demonstrate the transition of these standing waves via the power spectrum analysis from the application of ICA (Fig. 4c & d) and to characterize the neural substrates involved in sleep onset (Fig. 5).
Our present investigation reveals that cortical areas showing activations from the analysis of fMRI data that are in correlation with EEG events. An alternative approach is to apply the artifact-removed EEG signals reconstructed from the separated ICs (Srivastava et al., 2005; Mantini et al., 2007,), and subsequently apply the statistical parametric mapping methods such as low-resolution electromagnetic tomography (LORETA) for the localization of the corresponding cortical areas (Park et al., 2002; Pascual-Marqui et al., 2002). The comparison of the activation patterns based on the fMRI data and EEG signal from the LORETA method will be of great importance and constitute the subjects of our further research.
Another interesting finding from this exploratory study is that a full EEG montage consisting of 32 electrodes may not be necessary to detect a sleep pattern since the EEG rhythms were successfully separated from the artifacts by using the EEG signals from approximately 10 electrodes, excluding those for ECG and EOG signals (data not shown). In this regard, the performance change corresponding to the reduced number of electrodes needs to be investigated so that the computational load for the analysis of the EEG data can benefit from the reduced number of data sets.
The proposed approach may be potentially beneficial in several practical applications deploying a real-time neurofeedback scheme. In neurofeedback applications, individuals are expected to learn either up- or down-regulation of the regional brain activities in the presence of feedback on their performance. Examples of such applications include (1) the treatment of substance abuse and (2) calibration and optimization of brain-computer-interface (BCI) systems. An fMRI-based modality provides superior spatial resolution over an EEG-based modality, and thus region-specific brain activities such as the localized BOLD intensities or the size of activation patterns were fed back to the subject visually as an objective measure of performance (Wisekopft et al., 2003; DeCharms et al., 2005; Yoo et al., 2008). On the other hand, EEG may be more applicable in practical systems compared to fMRI since, for example, EEG system are cost-effective, compact in size, and simple to use (Wolpaw et al., 2002). One of the main problems commonly encountered in EEG-based systems, however, is the difficulty in identifying an appropriate EEG signature corresponding to the specific target brain regions of interest to the task so that these signatures can be fed back to the patient/participant in real-time. In this context, the proposed approach of the blind separation of the hidden EEG features can be gainfully applied, whereby the areas of activations corresponding to the identified EEG features are compared to the ROIs. Thus, the EEG features originating from the target brain area can be further utilized as the EEG signatures during the training sessions of EEG-based real-time neurofeedback.
It is worth to note that the presented study was intended to show the feasibility of the separation of EEG signatures originating from an abnormal brain region (for epilepsy data) and neurophysiological resting state (for sleep onset data) from the noisy EEG data recorded simultaneously with fMRI. The results from multiple subjects’ data sets are thus warranted in order to verify the efficacy of the method, in general. In addition, although the Infomax algorithm appears to provide the capability to successfully separate the unknown EEG signatures from the various artifacts, the result may need to be compared to the results from various learning algorithms for blind source separation such as the fast fixed-point algorithm (Hyvärinen et al., 2001) and the joint approximate diagonalization of eigen-matrices (JADE) algorithm (Cardoso, 1999), so that the performance can be further optimized.
In summary, we presented a method for the analysis of simultaneous EEG-fMRI data by deploying one of the blind source separation (BSS) techniques combined with more traditional GLM-based approach. Based on the results corresponding to the interictal epileptic discharges and sleep onset, the presented method seems to show the capability of identifying EEG signatures along with the corresponding source brain areas. This approach based on the BSS technique is particularly useful when a reference (task-) paradigm is unavailable.
The authors would like to thank Mr. Matthew Marzelli for his editorial work.
Source of Support: This work was partially supported in part by grants from NIH (R01-NS048242 to Yoo, SS and NIH U41RR019703 to Jolesz FA).
*This manuscript has not been published elsewhere and has not been submitted simultaneously for publication elsewhere.