Comparative Study of Different Pulse Artefact Correction Techniques during Concurrent EEG-FMRI using FMRIB

Simultaneous electroencephalography (EEG) and functional magnetic resonance imaging (fMRI) has become a widely used technique for studying brain function [1,2], capitalizing on the high spatial resolution of the former and excellent temporal resolution of the later. Over the last decade the applications of this technique have grown rapidly as new methods for improving data quality have been developed. Current applications include the study of resting state brain networks and the correlation of natural variations in externally stimulated neuronal responses measured using EEG and fMRI [3-6]. To date, the most widely explored clinical application is the identification of epileptic foci non-invasively [7,8]. Recently, the multi-modal EEGfMRI technique has also been used to investigate sleep and has been shown to have potential uses in the study of sleep disorders [9-11].


Introduction
Simultaneous electroencephalography (EEG) and functional magnetic resonance imaging (fMRI) has become a widely used technique for studying brain function [1,2], capitalizing on the high spatial resolution of the former and excellent temporal resolution of the later. Over the last decade the applications of this technique have grown rapidly as new methods for improving data quality have been developed. Current applications include the study of resting state brain networks and the correlation of natural variations in externally stimulated neuronal responses measured using EEG and fMRI [3][4][5][6]. To date, the most widely explored clinical application is the identification of epileptic foci non-invasively [7,8]. Recently, the multi-modal EEG-fMRI technique has also been used to investigate sleep and has been shown to have potential uses in the study of sleep disorders [9][10][11].
EEG measurements made in an MR scanner during concurrent magnetic resonance imaging are, however, plagued by large artefacts from two different sources: the pulse artefact (PA) caused by the pulsatile motion of blood, and gradient artefacts (GA) produced by temporally varying magnetic fields associated with the switched gradient waveforms used in MRI. The process by which this artefact is generated is well understood and is predictable because of its periodic nature. The temporal forms of these artefacts have been well studied, and their periodic nature has allowed the development of artefact correction methods based on subtraction of an average artefact template for the recording from each channel [12]. The average artefact subtraction (AAS) approach has proved extremely powerful, particularly when applied to the correction of the gradient artefact and is employed in the analysis of almost all EEG data acquired in combined EEG-fMRI experiments. The efficacy of this approach has been demonstrated in the literature, though a number of quality and practicality issues still remain [7,13]. Firstly, some residual artefacts remain on some channels. Allen proposed the use of adaptive noise cancellation (ANC) to remove these residuals; however, this approach does not remove all residual artefacts. Secondly, in order to minimize the residuals, a high sampling frequency is needed. Some unsatisfactory results are often obtained from commercial implementation of this algorithm even at sampling rates of 10 kHz. However, even if better quality data were to be achieved at such high sampling rates, the amount of generated data (especially in high electrode density experiments) limits the length of the experiments and causes practical problems when the data need to be analysed using third party software such as Matlab. In fMRI artefact slice template removal (FASTR), a unique artefact template for each slice artefact in each EEG channel is constructed and then subtracted. Each slice template is constructed as the local moving average plus a linear combination of the basic functions that describe the variation of residuals [14]. The basic functions are derived by performing temporal principal component analysis (PCA) on the artefact residuals and selecting the dominant components to serve as a basis set. This technique is demonstrated to be superior to AAS and applicable at a sampling rate as low as 2048 Hz.
The pulse artefact, which is linked to the cardiac cycle, is less well understood and significantly less predictable in nature than the GA.
In particular, the PA shows considerable differences when compared across subjects and can also vary in form across cardiac cycles in an individual subject [15][16][17]. The periodic nature of the PA means that it is also amenable to correction using AAS, but variation of the artefact across cardiac cycles reduces the efficacy that can be achieved

Abstract
Simultaneous acquisition of electroencephalography (EEG) and functional magnetic resonance imaging (fMRI) aims to combine the strength of both methods by investigating brain activity with the high temporal resolution of the former and the good spatial precision of the later. However, simultaneous EEG-fMRI measurements are technically challenging because both methods interfere with the other modality. The main confounding factors are the pulse artefact (PA) caused by pulsatile motion linked to the cardiac cycle and the gradient artefact (GA) produced by the temporally varying magnetic fields required for MR imaging. Both of these artefacts are orders of magnitude larger than the neuronal activity of interest, but their inherent periodicity facilitates reasonable artefact correction by post-processing techniques such as average artefact subtraction (AAS), Optimal Basis Set (OBS). While the different post processing methods for PA artefact are being used to-date, but their comparison were not studied yet. In this work, different methods -OBS, simple mean (AAS) and Gaussian-weighted mean (GWM) -implemented in open source FMRIB tool-box, were used for PA correction and to compare their performance in artefact correction while retaining the neuronal information. It has been found that, of the three different PA removal methods, OBS is better in preserving bio-signal while removing PA successfully.

Experiment
Experiment was carried out on a healthy human subject with the consent of local ethic committee to test the performance of artefact removal tool FASTR in removing GA and PA from the EEG data collected concurrently with fMRI. An EEG cap of 32 electrodes was placed over the subject's head and according to 10/20 system two electrodes (O1 and O2) were appeared to be useful in detecting EEG signal from the visual cortex. All channels were referenced to a common electrode at the FCz location. Thirty electrodes were being used for recording EEG signal whilst remaining two channels were EoG and ECG channels. EEG amplifier and MR scanner were synchronized so that artefact template can be reproducible which is necessary for improved artefact reduction.
EEG data were recorded inside the MR scanner while the subject was asked to open and close his eyes in periods of 10 seconds during the execution of the multi-slice EPI sequence. Slice-timing triggers ('slice' marker) from the MRI machine were recorded in the EEG data. The data were sampled at 2048 Hz, more than twice the maximum gradient artefact frequency where a typical EPI sequence is about 700-800 Hz. Forty fMRI volumes were collected with total scanning period of 2 minutes. In the EEG, each slice was marked by 'slice' marker which was, then, used to segment the EEG data to produce artefact template. The ECG trace was also recorded during EEG recording using a single bipolar channel (where pair of electrodes was attached to the chest -one at the mid-line and the other to the left of the heart) in conjunction with a SD32 MRI amplifier and using System98 software. The ECG trace was used for the R-peak detection that is required for PA correction. The EEG data were first cleaned for GA correction using FASTR and then residual artefact based on optimal basis set using PCA. Finally channels associated with visual sensitivity were used for validation purpose. This was done by a comparing the PA correction through evaluating the alpha power of one visual cortex channel (O1) after PA correction using the three correction methods-simple mean, Gaussian-weighted average and OBS methods.

Analysis
EEG data have been analysed in EEGLAB v.11.0.4.3b, a Matlab based open source software for biomedical data analysis, along with its plug-in 'FMRIB' [14,21]. In the EEG data recorded inside MR scanner, gradient artefacts were visible while scrolling the raw data in EEGLAB. The raw EEG data were exported to do further analysis in EEGLAB and Matlab. All the analysis were done on 30 channels excluding EoG and ECG channel which were not carrying brain signals.

Gradient artefact (GA) removal
FASTR offers a feature for removal of residual artefacts based on an optimal basis set using principal component analysis [14]. It is based on constructing a unique template for each artefact segment, in each channel, generated during the acquisition of a single fMRI slice. For artefact removal, it requires a slice timing event (or a section/volume trigger) is present in the data i.e. an event for each fMRI slice acquired. Low pass filtering of cut off, 100 Hz was used to remove high frequency gradient artefact. As required by algorithm for artefact removal correctly, data have been interpolated to 10 folds (20.48 kHz). Increasing in averaging window size provides the better approximation of the true artefact but good adaptation in any changes in the shape of the artefact (e.g., due to head movement) can be achieved with smaller window size with a degradation of biological information. Thus a trade off in window size of 30 has been selected. Adaptive noise cancellation (ANC) was used as suggested by Allen. The in PA correction via AAS, compared with the performance that can be achieved in correcting the GA. The variation of the PA waveform over time means that the artefact template used in AAS is generally formed by averaging over a small number of cardiac cycles. However, if the averaging is done over too few cycles then neuronal signals of interest may also be attenuated in the correction procedure; as a compromise, a sliding window template based on the average of around ten repetitions of the cardiac cycle is therefore, typically used for PA correction via AAS [12]. Given the limitations of using AAS for PA correction, it is not surprising that significant effort has been dedicated to devising improved techniques for PA correction. Much of this effort has been focused on blind source separation methods, such as independent component analysis (ICA) and optimal basis sets (OBS) analysis [14,[17][18][19]. Although several groups have reported success when using ICA for PA correction, other studies have shown less positive results using ICA [13,[17][18][19][20]. A possible reason for this discrepancy is the different field strengths of the scanners used in these studies. In summary, the spatial filtering approaches for PA correction developed till to-date such as ICA and PCA may not be as efficient as template methods. Moreover, due to the wide variety of post processing methods for PA artefact corrections are being used, the comparison of techniques were not studied yet. In this work, different PA correction methods (OBS, simple mean (AAS) and Gaussian-weighted mean), implemented in FMRIB tool-box available in open source EEGLAB software, were compared to evaluate their performance in artefact correction while retaining the neuronal information [14].

Data acquisition
EEG and electrocardiograph (ECG) data were recorded using the System PLUS EEG system and an SD32 MRI amplifier. The system is capable of recording from 30 common reference EEG channels and two bipolar channels to be used for electromyogram (EMG), ECG or electro-oculogram (EOG) recordings. All channels had 10 kV current limiting resistors and 0.15 Hz, 40 dB/decade high-pass filters to avoid dc coupling whereas 600 Hz, 20 dB/decade low-pass filters to protect against radio frequency (RF) artefact. A Sigma Delta analog-to-digital converter (ADC) with anti-aliasing filtering was used. The system had a sampling frequency of 2048 Hz, a dynamic range of 25.6 mV and a resolution of 12.2 nV. The head box containing the amplifiers, filters and ADC hardware is placed in the scanner room ( Figure 1) and the signal is transmitted via optical fibres to the acquisition workstation in the console room. FMRI was collected using a 3-T Varian Inova scanner (Palo Alto, CA). A standard, axial multi-slice EPI sequence was implemented with TR = 3 s and 21 axial slices per volume.

Pulse artefact (PA) removal
Gradient artefact corrected EEG data were then PA corrected using the methods available in FMRIB tool. This tool identifies QRS complexes locations of ECG data and then recorded the QRS events as 'qrs' marker on the data set which can be used for making artefact template. PA artefact was removed separately by the previously mentioned three methods (Simple mean, Gaussian weighted average and OBS). In simple mean and Gaussian weighted mean methods, no user-input is required other than the 'qrs' marker; however, the number of PC has to find out in OBS method. In simple mean method, averages of successive pulse artefacts around a contaminated data segment were taken and then subtract the result from the data which is consequently implementing average artefact subtraction (AAS). However, in Gaussian-weighted mean the artefacts were averaged after multiplying by a Gaussian window weights to emphasise the current artefact shape and reduce the effect of artefacts further. In OBS, PCA on a matrix of all the pulse artefacts in a channel was calculated, then take the first N PCs to form an optimal basis set describing the variations in the artefact. The OBS was then fitted and subtracted from each artefact. In this work, the default value of PC (4) has been used.
Fast fourier transforms (FFTs) of the one visual cortex representative channel O1 data were carried out for raw, GA and PA corrected data to verify the findings in frequency domain. The performance of pulse artefact correction was evaluated by comparing the time domain and frequency domain representation of PA corrected data for both eyes open and eyes close data epochs. Finally for further quantification, the mean root mean square (RMS) values over time for all channels were calculated using Matlab and artefact attenuation after different correction techniques were also calculated. Figure 2 shows a segment of the raw EEG data and GA corrected EEG data using FASTR method. It can be visualize clearly that the signal quality of raw data when there is no artefact correction has been performed is completely obscured by different artefacts and it is not possible to extract any neuronal activity from any of the channels (Figure 2). The magnitude of the artefacts is several millivolts whereas the magnitude of actual brain-EEG signal is around 200 µV. However same figure shows that after GA correction the magnitude of the EEG signal greatly reduced for each channel making it possible to be used for further processing (Figure 2). The PAs were revealed if checked carefully in some of the channels once the high frequency gradient artefact has been removed using the post-processing method in FASTR. It is also obvious that the magnitude of PA is small compared to GA at lower frequency ensuring the linked of this artefact to the cardiac cycle. Although alpha oscillation can be found on the above mentioned channels but in crucial case, still any decision on neuronal signals could be error prone.

Results and Discussion
In Figure 3, a comparison of the EEG data quality that can be achieved after gradient and pulse artefact correction using FMRIB toolbox (Figure 3). PA was corrected using OBS, simple mean and Gaussian weighted average respectively after detecting the QRS complex of the cardiac waveform from the ECG trace. It is clear that the amplitude of the remaining signals is much smaller and therefore neuronal signals are no longer obscured. Figure 3 show that alpha oscillation starts at 40 th second on different channels (e.g., O1, O2, P3, Pz, P4, PO3 and PO4 etc.) (Figure 3).
In Figure 4, the frequency spectrum of the raw, GA corrected (FASTR) and PA corrected (OBS) EEG data were shown (Figure 4). It clearly shows that the GA in the raw data occurs at distinct frequencies which were harmonics of the frequency of slice acquisition in the fMRI sequence, spanning the entire frequency range of the recording whereas the pulse artefact has a lower frequency than the gradient artefact (mainly below 10 Hz) as it is linked to the cardiac activity. It is also apparent from the Figure 4 that the FMRIB tool effectively removes artefact from the raw data to make it easier to extract biological information (Figure 4). It is also noticeable from Figure 4 that the low frequency variability in the EEG data has been greatly reduced after GA correction (Figure 4).  To investigate the performance difference between different PA techniques, a time series of 10 second segment of EEG data containing eyes open and close for channel O1 and O2 using three methods were plotted in Figure 5. It can be visualized from the figure that all three methods are doing well in preserving alpha signal while the subject has closed his eyes however, OBS preserves the signal shape better than other two methods; i.e., preserves the brain signal more accurately than other two methods. Figure 6A clearly reveals the fact of preserving brain signal in frequency domain (8)(9)(10)(11)(12)(13) Hz alpha band) while reducing artefacts in channel O1 and Figure 6B shows no trace of alpha while the eyes were opened ( Figure 6A and 6B). However, in both the case OBS was shown to attenuate residual artefacts more than other two methods. To test the effectiveness of the different PA correction methods, a quantitative measurement was also done using Matlab for the raw, GA corrected and PA corrected data using different methods. The mean RMS values over 30 channels were calculated using MATLAB and the corresponding attenuation for different correction method with respect to raw RMS were also calculated. Mean attenuation for GA corrected data and simple mean, Gaussian-weighted mean and OBS corrected data were 29.56, 29.91, 29.90 and 29.96 dB respectively. This clearly reveals that OBS outperforms over simple mean and Gaussian-weighted mean in correcting pulse artefact.

Conclusion
Simultaneous EEG-fMRI is a powerful tool for studying brain function, as the high temporal resolution of EEG can be combined with the high spatial resolution of fMRI. To date, a number of studies have used this multi-modal approach to gain a better understanding of brain function. In clinical studies the main use of the technique has been to investigate the foci of inter-ictal epileptic discharges which are inherently difficult to localize non-invasively. These examples show the power of this multimodal imaging tool. However, EEG data acquired during simultaneous fMRI are affected by several artifacts, including the gradient artefact (due to the changing magnetic field gradients required for fMRI), the pulse artefact (linked to the cardiac cycle). While numerous post-processing methods are being used or under developed for correcting the gradient and pulse artifacts but the comparison of different types of PA correction methods (OBS, AAS and GWM) were not studied yet, which was studied in this work. It has been found that, of the three methods of removing PA, implemented in FMRIB tool box, OBS is recommended for PA removal, as it is better in preserving bio-signal while removing PA successfully. Moreover, a significant amount of research work is being carried out in the field of combined EEG-fMRI to improve the data acquisition and analysis methods. The outcomes mentioned here will help to explicate important queries in neuroscience, combination of the high spatial resolution of fMRI and excellent temporal resolution of EEG. Comparison of OBS with the other types of PA correction techniques such as OBS-ICA will be studied in future work.