Functional source separation improves the quality of single trial visual evoked potentials recorded during concurrent EEG-fMRI

NeuroImage 50 (2010) 112–123 Contents lists available at ScienceDirect NeuroImage j o u r n a l h o m e p a g e : w w w. e l s e v i e r. c o m / l o c a t e / y n i m g Functional source separation improves the quality of single trial visual evoked potentials recorded during concurrent EEG-fMRI Camillo Porcaro a,b,c,d,⁎, Dirk Ostwald a,b, Andrew P. Bagshaw a,b a School of Psychology, University of Birmingham, Birmingham, UK Birmingham University Imaging Centre (BUIC), University of Birmingham, Birmingham, UK ISTC-CNR, Ospedale Fatebenefratelli, Isola Tiberina, 00186 Rome, Italy d ITAB-Institute for Advanced Biomedical Technologies, “G. D'Annunzio” University, Chieti, Italy b c a r t i c l e i n f o a b s t r a c t EEG quality is a crucial issue when acquiring combined EEG-fMRI data, particularly when the focus is on using single trial (ST) variability to integrate the data sets. The most common method for improving EEG data quality following removal of gross MRI artefacts is independent component analysis (ICA), a completely blind source separation technique. In the current study, a different approach is proposed based on the functional source separation (FSS) algorithm. FSS is an extension of ICA that incorporates prior knowledge about the signal of interest into the data decomposition. Since in general the part of the EEG signal that will contain the most relevant information is known beforehand (i.e. evoked potential peaks, spectral bands), FSS separates the signal of interest by exploiting this prior knowledge without renouncing the advantages of using only information contained in the original signal waveforms. A reversing checkerboard stimulus was used to generate visual evoked potentials (VEPs) in healthy control subjects. Gradient and ballistocardiogram artefacts were removed with template subtraction techniques to form the raw data, which were then subjected to ICA denoising and FSS. The resulting EEG data sets were compared using several metrics derived from average and ST data and correlated with fMRI data. In all cases, ICA was an improvement on the raw data, but the most obvious improvement was provided by FSS, which consistently outperformed ICA. The results show the benefit of FSS for the recovery of good quality single trial evoked potentials during concurrent EEG-fMRI recordings. © 2009 Elsevier Inc. All rights reserved. Article history: Received 13 July 2009 Revised 27 November 2009 Accepted 1 December 2009 Available online 16 December 2009 Keywords: Visual evoked potential (VEP) Electroencephalography (EEG) Functional magnetic resonance imaging (fMRI) Independent component analysis (ICA) Functional source separation (FSS) Introduction The simultaneous measurement of EEG and fMRI is an attractive, non-invasive technique for the investigation of human brain function, with the potential to offer a higher spatiotemporal resolution than either method alone. It is increasingly widely used as a tool in cognitive and sensory neuroscience (e.g. Debener et al., 2005; Eichele et al., 2005; Bénar et al., 2007; Mulert et al., 2008; Goldman et al., 2009; Olbrich et al., 2009; Ostwald et al., 2009) and can also shed light on the properties of the underlying neurovascular coupling which, particularly at the macroscopic level at which scalp EEG and whole brain fMRI are measured, are not fully understood (Wan et al., 2006). However, if the potential strengths of EEG-fMRI are to be fully realized, and new methods for data integration developed and exploited, it is vital that good quality EEG and fMRI data are recorded. In particular, EEG data acquired in the MRI scanner are strongly contaminated by artefacts of biological and non-biological origin that ⁎ Corresponding author. School of Psychology, University of Birmingham, Edgbaston, Birmingham B15 2TT, UK. E-mail address: c.porcaro@bham.ac.uk (C. Porcaro). 1053-8119/$ – see front matter © 2009 Elsevier Inc. All rights reserved. doi:10.1016/j.neuroimage.2009.12.002 may prevent the correct determination of the characteristics of the brain signals that are of primary interest. There are several artefacts that contaminate the measurement of neurophysiological EEG and that need to be removed from the recordings before further analysis. Specific to the MRI environment are gradient artefacts (GA) and ballistocardiogram artefacts (BCG), while ocular artefacts (OA) and electrode artefacts (EA) are present in the EEG acquired inside and outside of the scanner. The most widely used techniques to reduce the effects of GA and BCG are variations of template averaging approaches (Allen et al., 1998, 2000), with independent component analysis (ICA) often used as an alternative or secondary step (Debener et al., 2006; Mantini et al., 2007). ICA is a blind signal processing technique that can be used to recover independent sources (or components) that contribute to form the measured signal, upon the assumption that they are linearly mixed (Comon, 1994; Hyvärinen et al., 2001; Hyvärinen, 1999; James and Hesse, 2005). ICA is widely used for characterization of brain activity (Delorme and Makeig, 2004; Makeig et al., 2002;2004a,b; Medaglia et al., 2009) and is increasingly the standard method when performing single trial (ST) EEG-fMRI (Debener et al., 2005; Eichele et al., 2005; Bénar et al., 2007; Mulert et al., 2008). It has also been successfully used for removal of eye blinks, eye movements and C. Porcaro et al. / NeuroImage 50 (2010) 112–123 113 electrode artefacts (Barbati et al., 2004; Mantini et al., 2007; Porcaro et al., 2009). The aim of the current study was to apply a new approach, functional source separation (FSS), to the problem of improving the signal quality of EEG data recorded in the MRI environment. FSS can be seen as a semi-blind extension of ICA that incorporates some prior knowledge about the responses of interest (Barbati et al., 2006; Tecchio et al., 2007). The aim of FSS is to enhance the separation of relevant signals without renouncing the advantages of using only information contained in the original signal waveforms (Tecchio et al., 2007; Barbati et al., 2008; Porcaro et al., 2008; 2009;). Following removal of GA and BCG with standard techniques, data were further processed using ICA and FSS and the results compared both at the level of the average and the single trials. As a final comparison between the EEG preprocessing strategies, the EEG data were correlated with fMRI data extracted from activated voxels to determine whether improvements in EEG data quality also resulted in better correlations with fMRI. Materials and methods Subjects As part of a study to investigate the link between EEG and fMRI at the level of single trials, fourteen subjects were paid for their participation. Of these, ten participants (4 female, 27.6 ± 5.17 years mean age ± SD) had additional EEG data recorded outside of the MRI scanner for the purpose of investigating the EEG signal quality, and are Fig. 1. Experimental setup. (First row) Single-trial experimental design of the visual stimulation experiment. A hemi-field checkerboard of either low or high contrast was presented at time 0 and reversed after 500 ms. After another 500 ms, the checkerboard was removed. On a subset of the trials, the fixation cross changed to an X at the onset of a randomly sampled EPI volume acquisition, and the observer was asked to indicate this change by a button press. (Second row left) Percentage signal change of the fMRI signal for the High (continuous line) and low (dotted line) contrast stimuli, averaged over subjects. (Second row right) fMRI statistical map for a representative subject, high contrast stimuli. (Third row left) VEP for high (continuous line) and low (dotted Line) contrast stimuli, averaged over all subjects. The first vertical dashed line indicates stimulus onset, while the second indicates the reversal of the checkerboard. (Third row right) Localization of the grand average P100 VEP. 114 C. Porcaro et al. / NeuroImage 50 (2010) 112–123 the subjects of the current study. Data from two subjects were subsequently discarded (both male) because of the poor quality of the EEG. The resulting group consisted of four males and four females (28.4 ± 5.24 years mean age ± SD). Written informed consent was obtained from all participants and the protocol was approved by the Research Ethics Board of the Birmingham University Imaging Centre. Stimuli Left hemi-field reversing checkerboards were presented at a spatial frequency of 2 cycles per degree of visual angle at two different contrast levels, high (CMichelson = 1) and low (CMichelson = 0.25). Stimuli were presented with a central fixation cross for 1 s (i.e. one checkerboard reversal, see Fig. 1). Contrasts were randomized, and in total 85 trials of each contrast were presented. For the data recorded outside the scanner, the interstimulus interval (ISI) was randomized between 1.5 and 2.5 s, while for the data recorded inside the scanner the ISI was 16.5–21 s, discretized to 1.5 s (MR repetition time, see below). Design and procedure Individual trials of the experiment consisted of a single presentation of the checkerboard stimulus with phase reversal at a temporal frequency of 2 Hz followed by a fixation period. Since the ISI for the data recorded outside of the scanner was much shorter than that for inside, the data recorded outside were acquired in one block taking approximately 8.5 min. For the data recorded inside the scanner, individual runs consisted of 17 trials per contrast with fixation periods at the beginning and at the end, amounting to a total session length of 441 volumes × 1.5 s, i.e. 11 min. Five of these runs were acquired in each participant, resulting in the 85 trials per contrast. In order to maintain the observer's attention, a simple fixation task was performed: on a random selection of half of the trials of a given session, the fixation cross changed from a plus sign to an X during the fixation period at random time points, discretely (1.5 s) and uniformly sampled from the interval of 4.5–16.5 s after stimulus onset. The observer's task was to report the change in fixation by a button press using the index finger of the right hand (see Fig. 1, first row). Hit rate and number of false alarms were presented to the observer at the end of each session. Stimuli were presented and behavioural data collected using Psychotoolbox3 (Brainard, 1997; Pelli, 1997) for Matlab (The Mathworks, Natick, MA). Stimulus presentation timing was controlled by the MRI scanner volume trigger. MRI data acquisition The experiment was conducted at the Birmingham University Imaging Centre using a 3-T Philips Achieva MRI scanner. An initial T1weighted anatomical (1 mm isotropic voxels) and T2⁎-weighted functional data were collected with an eight channel phased array head coil. EPI data (gradient echo pulse sequence) were acquired from 20 slices (2.5 × 2.5 × 3 mm resolution, TR = 1500 ms, TE = 35 ms, SENSE factor = 2, flip angle = 80°), providing approximately half brain coverage in the dorsal–ventral direction. Slices were oriented parallel to the AC–PC axis of the observer's brain and positioned to cover the entire occipital cortex. EEG data acquisition EEG data were recorded using a 64-channel MR compatible EEG system (BrainAmp MR Plus, Brain Products, Munich, Germany), which incorporates current limiting resistors of 5 kΩ at the amplifier input and in each electrode. The EEG cap consisted of 62 scalp electrodes distributed according to the 10–20 system and two additional electrodes, one of which was attached approximately 2 cm below the left collarbone for recording the ECG, while the other was attached below the left eye (on the lower orbital portion of the orbicularis oculi muscle) for measurement of the electrooculogram. Data were sampled at 5000 Hz. Impedance at all recording electrodes was less than 20 kΩ. The EEG data acquisition clock was synchronized with the MRI scanner clock using Brain Product's SyncBox, resulting in exactly 7500 data points per EPI-TR interval. EEG data preprocessing The GA was removed using the template averaging approach as implemented in Brain Vision Analyzer 1.05 (Brain Products, Munich, Germany). BCG artefact correction was performed using the Optimal Basis Set method (Niazy et al., 2005) as implemented as a plug-in to EEGLAB (Delorme and Makeig, 2004). The data from the five runs were then concatenated, down-sampled to 500 Hz and low-pass filtered (25 Hz). These data from which the gross MRI artefacts had been removed were then further processed using ICA and FSS (Fig. 1). Artefact attenuation in EEG data by independent component analysis (ICA) ICA (Comon, 1994) is a generative ‘latent variable’ model that describes how the observed data are generated by a process of mixing the underlying unknown sources. The sources (independent components [ICs]) are assumed to be statistically independent and nonGaussian. Since the observed mixed signals will tend to have more Gaussian amplitude distributions, ICA strives to find a separation matrix that minimizes the gaussianity of the results, thus optimally separating the signals. The set of EEG signals X is assumed to be obtained as a linear combination (through an unknown mixing matrix A) of statistically independent non-Gaussian sources S (at most, one Gaussian source): X = AS ð1Þ Sources S are estimated (up to arbitrary scaling and permutation) w by independent components S as: w S = WX ð2Þ where the unmixing matrix W is estimated along with the components. In this application, the fastICA algorithm (Hyvarinen et al., 2001; Hyvarinen, 1999) was applied to the EEG data matrix after GA and BCG artefacts rejection. After estimation of ICs, a semiautomatic procedure (Barbati et al., 2004; Porcaro et al., 2009) was applied to identify and eliminate artefactual noncerebral activities, i.e. eye movements, residual BCG, residual GA, etc. In this preprocessing step, no dimensionality reduction was performed in estimating ICs. After the identification of artefactual ICs, the cleaned data at the scalp electrodes were obtained by retroprojecting all the ICs except for those identified as representing artefacts: w EEGrecICA = Ak Sk ð3Þ where Ak is the estimated mixing vector (matrix A of equation 1) for w w the source Sk and EEGrecICA is the resulting Sk retro-projection on the channels space. Functional source separation (FSS) As in the ICA approach, FSS starts from an additive hidden source model of the type in equation 1, where X represents the observed EEG data, S is the underlying unknown sources and A is the source-sensor coupling matrix to be estimated. Additional information to a standard ICA model is used to bias the decomposition algorithm towards solutions that satisfy physiological assumptions. In other words, the aim of FSS is to enhance the separation of relevant signals by C. Porcaro et al. / NeuroImage 50 (2010) 112–123 115 exploiting some a priori knowledge without renouncing the advantages of using only information contained in the original signal waveforms. A modified (with respect to standard ICA) contrast function is defined: F = J + λRFS ð4Þ Table 1 Signal-to-noise ratio across methods. High contrast Raw out S1 S2 S3 S4 S5 S6 S7 S8 Mean SD 51.6 39.0 28.9 31.1 38.5 50.6 38.2 43.7 40.2 8.2 Raw in 22.0 7.9 30.1 11.2 17.0 7.4 34.6 19.9 18.8 10.0 ICA in 19.3 28.1 38.9 15.9 29.2 4.9 29.7 31.8 24.7 10.7 FSS in 62.1 58.2 38.1 46.8 50.1 59.6 44.1 50.6 51.2 8.3 Low contrast Raw out 31.7 16.4 21.4 39.6 12.2 40.4 17.7 27.9 25.9 10.7 Raw in 23.3 19.5 4.2 1.0 1.7 14.2 8.0 23.8 12.0 9.5 ICA in 24.7 15.1 5.7 1.0 10.1 4.7 10.4 29.3 12.6 9.9 FSS in 47.5 33.0 18.4 20.7 16.1 17.9 30.9 53.4 29.7 14.3 where J is the statistical index normally used in ICA, while RFS accounts for the prior information used to extract a single source. According to the weighting parameter, λ, it is possible to adjust the relative weight of these two aspects. In this study, λ was chosen equal to 1000 in all cases, as detailed in Porcaro et al. (2008). Briefly, λ was chosen to both minimize computational time and maximize RFS. Moreover, since prior information about the sources may also be described by a non-differentiable function, the new contrast function F is optimized by means of simulated annealing (Kirkpatrick et al., 1983). This does not require the use of derivatives and performs global optimization, while the gradient-based algorithms usually employed in ICA only guarantee local optimization. In the present work, consistent with the visual evoked potential (VEP) activity under investigation, an ad hoc functional constraint was maximized around the principal peak (P100) of the VEP. The FSS procedure was applied to the same EEG data that was used for the ICA preprocessing step described above. Functional constraints The functional constraint R was defined as: RðFSP100 Þ = tk + Δ2 tk X tk − Δ1 tk The signal-to-noise ratio of the average VEP was calculated for each subject and method. Mean and Standard Deviation (SD) over subjects are also given. also calculated. The signal-to-noise ratio (SNR) of the average VEPs was calculated, with the noise level calculated between −100 and 0 ms and the signals between 1 and 1000 ms. Comparisons were made between the raw data before scanning, raw data in the scanner, data in the scanner after preprocessing by ICA and data in the scanner using FSS (Table 1). Discrepancy The discrepancy was defined as the difference between the EEG data after GA and BCG rejection (i.e. the raw data in the scanner) and the data obtained from equation 3 for ICA data and equation 6 for FSS data: DiscrepancyICA = EEG − EEGrecICA DiscrepancyFSS = EEG − EEGrecFSS ð7Þ jEAðFSP100 ; t Þ j − 0 X − 100 jEAðFSP100 ; t Þj ð5Þ with the evoked activity, EA, computed by averaging signal epochs of the source FSP100, triggered on the visual stimulus (t = 0); tk is the time point with the maximum electric potential around 100 ms after the stimulus onset on the maximal original EEG channel; Δ1tk (Δ2tk) is the time point corresponding to a signal amplitude of 50% of the maximal value before (after) tk. The baseline was computed in the time interval from −100 to 0 ms (Tecchio et al., 2007). The precise value of each latency tk was chosen for each subject, corresponding to the maximum electric potential in the time interval of interest (80–120 ms). The source was then retro-projected to obtain its electric potential distribution at the scalp electrodes: EEGrecFS = AP100 FSP100 ð6Þ To quantify the level of residual response to the visual stimulation after the preprocessing strategies (i.e. removal of artefactual ICs and extraction of the functional source), a mean discrepancy index was calculated. This was obtained as the ratio between the mean R computed on the discrepancy matrix and the mean R across the EEG electrodes after GA and BCG rejection: P DiscrIndex = i RðDiscrepancyMethod Þ2 Á2 PÀ i RðEEGÞ ð8Þ R is the reactivity index defined as: R= R Δ2 tk + Δ1 tk + 1 ð9Þ where AP100 is the estimated mixing vector (matrix A of equation 1) for the functional source and EEGrecFS are the retro-projections on the electrodes of the estimated FSP100. Data evaluation To evaluate the quality of the data following ICA and FSS, five criteria were used: the VEP, discrepancy, normalized cumulative signal-to-noise ratio (ncSNR), localization and single trial variability. To further evaluate the performance of ICA and FSS, these metrics were also applied to the raw data (i.e. only GA and BCG removed using standard techniques) and to that recorded before the MRI scanning session. In order to facilitate comparison, the analyzed data were taken from a single electrode selected on an individual subject basis based on the maximum of the voltage field (equation 3 for ICA method and equation 6 for FSS method) at the latency of P100 peak. VEP The VEPs for high contrast (HC) and low contrast (LC) stimuli were calculated for all subjects, and the grand average across subjects was R is defined in Eq. (5) and tk refers to the latency of each extracted source. The index i runs upon all the channels (see Table 2). Table 2 Discrepancy between ICA and FSS. High Contrast ICA In S1 S2 S3 S4 S5 S6 S7 S8 Mean SD 8.4 2.7 8.9 1.4 4.2 2.1 5.1 7.1 6.0 4.9 FSS In 2.8 1.2 1.5 1.1 1.4 1.1 1.3 2.3 1.6 0.6 Low Contrast ICA In 5.4 4.3 4.8 1.3 2.7 1.6 4.1 5.2 4.2 2.6 FSS In 3.9 1.5 1.0 1.0 1.3 1.1 1.1 2.8 1.7 1.1 The discrepancy index was calculated for each subject and for ICA and FSS. Mean and standard deviation (SD) over subjects are also given. 116 C. Porcaro et al. / NeuroImage 50 (2010) 112–123 Normalized cumulative signal-to-noise ratio (ncSNR) In order to evaluate the performance of the ST data, the SNR for each trial and for all the different data sets was calculated. Again, the performance was evaluated on the single electrode selected on an individual subject basis based on the maximum of the voltage field (Eqs. (3) and (6), respectively, for ICA and FSS methods) on the P100 peak. The SNR was calculated within the same time windows introduced previously for the average VEP. Having calculated the SNR trial by trial, a cumulative sum was computed in order to compare all of the data. This quantity was normalized by the number of trials, in order to give an idea of the magnitude of the ST SNR. Localization and single trial variability The raw, ICA and FSS data were submitted to a source localization algorithm (sLORETA) (Pascual-Marqui 2002) as implemented in CURRY 6 (Neuroscan, Hamburg, Germany, http://www.neuroscan. com/). sLORETA was performed for the P100 peak using a regular grid with a spacing of 4 mm throughout the brain region and a four shell spherical head model. The results were projected onto the template brain of the Montreal Neurological Institute (MNI) within CURRY. To estimate the spatial origin of the ICA and FSS data, the source was first retro-projected to obtain their field distribution at the electrodes. These data were then subjected to sLORETA. Finally, the quality of the single trial responses was visually examined using ERPImage plots as implemented in EEGLAB (Delorme and Makeig, 2004). Correlation between average VEPs and HRFs In order to address the issue of whether improving EEG data quality affects the correlation of the EEG and fMRI data, single trial haemodynamic response functions (HRF) were extracted from spherical volumes of interest (VOI, radius 5 mm) centred on the maximally responsive voxel in the fMRI statistical map. Briefly, fMRI data were preprocessed using standard techniques in SPM5 (Wellcome Trust Centre for Neuroimaging, UCL), including motion and slice timing corrections, spatial normalization and smoothing (5 mm Gaussian kernal). Statistical maps were generated using a general linear model and HRF data were extracted from VOI centred on the voxel with the highest t-value when using the contrast including both stimulus contrasts. Full details of the HRF extraction are given in Ostwald et al., (2010). These ST HRFs were averaged for each subject and compared with the data extracted using the different EEG preprocessing techniques. Two types of correlation were performed: (1) between the amplitudes of the P100 and the HRF, where the precise value of each latency was chosen for each subject, corresponding to the maximum signal in the time interval of interest (80–120 ms) for the VEP and between 2 and 9 s for the HRF, and (2) between the normalized area of the VEP over the same time interval, calculated as the sum and normalized with respect to the window length, and the normalized area of the HRF, also centred on the maximum peak and normalized with respect to the window length. Statistical analysis The effect of preprocessing strategy on the SNR of the average VEPs (Table 1) and the ncSNR (Fig. 4) was assessed with one-way ANOVA, separately for the high and low contrast data. Differences between the discrepancy associated with ICA and FSS (Table 2) were assessed using paired samples t-tests for high and low contrast. For the comparisons between EEG and fMRI data, Pearson correlation coefficients were calculated for the different EEG preprocessing strategies (Fig. 7). This was done for both the amplitudes and areas, as defined above. All statistical tests were performed using SPSS 16.0 (SPSS Inc, Chicago IL, USA) and the significance level was set at p b 0.05. Fig. 2. Visual evoked potential. Time course of the stimulus-averaged VEP in the [− 100 1000] ms time period following high (top row) and low (bottom row) contrast stimulation. Average VEPs are shown for each subject (coloured lines) along with the grand average across subjects (black line). As in Fig. 1, the first vertical dashed line indicates stimulus onset, while the second indicates the reversal of the checkerboard. C. Porcaro et al. / NeuroImage 50 (2010) 112–123 117 Fig. 3. Discrepancy. For a representative subject and for high (top) and low (bottom) contrast stimulation, averaged data are shown in the time window [− 100, 1000] ms (onset at t = 0 ms; reversal at t = 500 ms). Left: raw data. Centre: retro-projected data obtained by ICA (top) or FSS (bottom). Right: raw data minus ICA data (top) and raw data minus FSS data (bottom), i.e. discrepancy. Results Visual evoked potentials (VEPs) Fig. 2 shows the average VEPs for each individual subject (coloured lines) and the grand average over all subjects (thick black). The data for the individual subjects show a considerable amount of variability, which is most evident in the raw data recorded inside the scanner. Comparing the raw data acquired inside and outside of the MRI scanner shows that there is much more variability when recording within the MRI environment. Given that these are the same subjects undergoing the same stimulation paradigm, it is evident that the primary cause of this increased variability is the reduction in EEG data quality caused by the various MRI artefacts. This inter-subject variability is improved with ICA, but the most obvious improvement is between ICA and FSS, with FSS showing similar, or less, variability, than the data recorded outside of the scanner. It is also worth noting that the grand average VEPs are very similar for the different methods, indicating that the grand average VEP is not a good measure of the underlying data quality. In all subjects, the functional source FSP100 was successfully extracted and the data reconstructed with this component shows a clear VEP morphology (see Fig. 2 last column). Moreover, the SNR in Table 1 is higher for FSS than for the other methods (raw and ICA data) and comparable with the data acquired outside the scanner. The statistical analysis confirmed these observations. A one-way ANOVA for the high contrast data indicated a significant difference between the means (F(3,28) = 19.8, p b 0.001), with post hoc tests demonstrating that FSS had a significantly higher SNR than the ICA or the raw data acquired in the scanner (p b 0.001 in both cases). FSS and the raw data recorded outside of the scanner did not differ (p N 0.1). Similarly for the low contrast data, the one-way ANOVA demonstrated a significant difference between the groups (F(3,28) = 5.2, p = 0.005). The FSS data had a significantly higher SNR than the ICA (p = 0.03) and raw data acquired inside the scanner (p = 0.02). Again, there was no significant difference between the SNR of the FSS data and that recorded outside of the scanner (p N 0.1). For neither contrast was there a significant difference between the raw data acquired inside the scanner and the ICA data. 118 C. Porcaro et al. / NeuroImage 50 (2010) 112–123 obvious P100 in the ST data. ICA does not improve this very much but, consistent with the other measures that have been shown, the difference between ICA and FSS is obvious. This is even clearer for the low contrast data (Fig. 6, bottom). The raw and ICA ST data are both extremely noisy, with little evidence of ST VEPs. The FSS ST plot, however, clearly shows a consistent P100 across trials. The low contrast data also demonstrates that improving the ST data quality can affect the localization of the average VEP, since the localization of the FSS data is must more realistic than that of the raw or ICA data. When the data quality is good, as shown for the best subject in Fig. 5, even the raw data show reasonably clear ST VEPs, and this is not improved to any great extent with ICA. Even in this case, however, FSS is able to provide some improvements and additional noise reduction. Correlation between average VEPs and HRFs Fig. 4. Normalized cumulative signal-to-noise ratio. The normalized cumulative signalto-noise ratio (ncSNR) is shown for both high (left) and low (right) contrast, averaged over all subjects. The ncSNR is shown for raw data recorded outside of the scanner (black line: raw data Out), raw data recorded inside of the scanner (green line: raw data in), data recorded in the scanner after preprocessing by ICA (blue line: ICA in) and after FSS (red line: FSS in). Discrepancy Analysis of the discrepancy (Fig. 3, Table 2) demonstrates that FSS source extraction was associated with minimal residual activity, i.e. the extracted source described practically all of the evoked response contained in the original EEG data. Paired sample t-tests indicated that for both contrasts the discrepancy associated with the ICA data was larger than that for the FSS data (t(7) = 3.91, p = 0.006 for high contrast data; t(7) = 4.47, p = 0.003 for the low contrast data), indicating that the FSS source accounted for more of the stimulus-related activity. Normalized cumulative signal-to-noise ratio Fig. 4 shows the ncSNR. This measure summarizes the quality of the data at the level of the individual trials and is of primary importance for the application to ST EEG-fMRI. Examination of Fig. 4 leads to a similar conclusion to that based on the average VEP. The FSS data are of similar quality to that recorded outside of the scanner and considerably better than the ICA data, while the ICA data are an improvement on the raw data. The one-way ANOVA corroborated these observations. For the high contrast data, there was a significant difference between the ncSNR for the different data sets (F(3,280) = 87.9, p b 0.001), with post hoc comparisons indicating that the FSS data had a significantly higher ncSNR than both the raw and ICA data (p b 0.001 in both cases). The difference between the FSS data and that recorded outside of the scanner was not significant (p N 0.1). For the low contrast data, the results were similar, with a significant difference between the data sets (F(3,280) = 47.4, p b 0.001). The ncSNR of the FSS data was significantly larger than all other data sets (p b 0.001 for all paired comparisons). Localization and single trial variability Finally, Figs. 5 and 6 present a final comparison of single trial (ERPImage) plots and source localization of average VEPs. Two individual subjects are shown, selected as those with the best and worst SNR in the raw data VEP. The source localization of the average VEPs appears to be relatively robust, with similar localization for different preprocessing methods. However, the level of single trial variability as assessed by the ERPImage plots is clearly very dependent on the data preprocessing (Figs. 5 and 6 last columns). In the worst case (Fig. 6, top), the raw data are clearly very noisy, with little Fig. 7 shows the comparison between EEG and fMRI data for the different EEG preprocessing strategies. Both the area (Fig. 7a) and the amplitude (Fig. 7b) values of the evoked responses were compared. For visualization and to facilitate comparison between EEG methods, the amplitudes and areas were centred and scaled. No significant correlation was found for the comparison between P100 and HRF amplitude (Fig. 7b). However, the area values were significantly correlated (Fig. 7a) for the FSS data (R = 0.83 and 0.87 for high and low contrast, p = 0.01 and p b 0.01 two-tailed), with a marginally nonsignificant correlation for the ICA data (R = 0.69 and 0.73 for high and low contrast, p = 0.06 and 0.04 two-tailed). There was no correlation for the raw data (Figs. 7a and b). Discussion As it has become technically less demanding to record EEG within the MRI scanner, and commercial systems have become available which allow reasonable quality data to be straightforwardly acquired, the focus of research has shifted toward the development of optimal strategies for data integration. While it is at least conceptually clear that EEG-fMRI can provide improved spatiotemporal resolution compared to existing methods, it is less obvious how the data should be combined in order to achieve this goal. Often, data integration relies upon the use of some features of the EEG data, whether properties of evoked potentials or spectral power variations (Laufs et al., 2003, Debener et al., 2005, 2006, Eichele et al., 2005, Bénar et al., 2007), which are then used to form regressors for a standard general linear model analysis of the fMRI data. The underlying neurophysiology and biophysics relating the macroscopic measures of EEG and fMRI are not currently sufficiently developed to allow a principled identification of the features which best link the two, and this type of approach is complicated by the reduction in data quality caused by simultaneous recording. The development of additional methods to improve the quality of EEG data acquired in the MRI scanner is therefore an ongoing area of research. The most commonly used method for denoising EEG data following gross artefact removal is independent component analysis (ICA). ICA is a completely blind source separation algorithm which incorporates no prior information about the part of the EEG signal that is of interest. Since in many experiments it is clear from the outset whether ST evoked responses, spectral power, etc., are likely to be most informative, a natural extension of ICA is to include this prior information in the data decomposition. This is done by the functional source separation (FSS) algorithm applied in this study, which trades some of the explorative, data-driven strength of ICA in exchange for a more robust decomposition of the signal of primary interest. One expectation might be that this would allow a better separation of sources of interest, while it has the additional advantage over ICA that component selection is not an issue, since only one component is identified. Obviously, this approach assumes that the source of C. Porcaro et al. / NeuroImage 50 (2010) 112–123 119 Fig. 5. Localization and single trial variability (best case). Results for the subject with the highest SNR in the raw data. sLORETA localization and ERPImage plots are shown for the different data sets and preprocessing strategies. sLORETA results are shown superimposed on the MNI brain template. 120 C. Porcaro et al. / NeuroImage 50 (2010) 112–123 Fig. 6. Localization and single trial variability (worst case). Results for the subject with the lowest SNR in the raw data. sLORETA localization and ERPImage plots are shown for the different data sets and preprocessing strategies. sLORETA results are shown superimposed on the MNI brain template. C. Porcaro et al. / NeuroImage 50 (2010) 112–123 121 Fig. 7. Correlation between VEPs and HRFs across subjects. Pearson correlation (two-tailed) between EEG and fMRI data was calculated across subjects both for the amplitudes and the areas of the evoked responses. The data was centred and scaled to facilitate the comparison across the methods. interest in known a priori, which is not always true, in which case the blind decomposition provided by ICA provides a more explorative procedure. However, given the data quality issues associated with recording EEG in the MRI scanner, the benefits provided by FSS may outweigh this loss of a purely data-driven analysis. In the current study, simultaneously acquired EEG-fMRI data were used to compare the performance of FSS with that of ICA over a number of different metrics, with the data recorded outside of the scanner providing a benchmark for the data quality in each subject. These metrics examined the data at the average and ST level, including assessment of source localization provided by sLORETA. A very similar conclusion was reached in all cases: the ICA data were indeed of improved quality with respect to the raw data acquired inside the MRI scanner, but the most pronounced improvement in data quality was consistently provided by FSS. The FSS data approached the quality of the minimally processed data acquired outside of the scanner prior to the scanning session. Particularly notable was the improvement of data that in its raw form was initially of low quality, as demonstrated in Fig. 6. Both the raw and ICA data from this subject were very noisy, and often it is these marginal data sets which may either be excluded from a study or, if included, reduce the strength of any observed effects. The ability to extract good quality ST parameters from a data set such as this is a considerable advantage considering the time, effort and expense required to scan a subject with EEG-fMRI. A natural concern is whether by applying this constraint to the ICA decomposition and extracting a single source, some of the signal of interest has been lost. This was addressed in Fig. 3 and Table 2, which summarize the discrepancy (i.e. the difference between the raw data and the FSS or ICA data). According to this metric, the residual signal not extracted by FSS is comparable to the level of the baseline, indicating that it is not stimulus related, i.e. it is residual noise. In the same figure, it is demonstrated that the discrepancy results from ICA and FSS are comparable, indicating that FSS is not removing stimulusrelated signal that ICA retains. It is of note that although the constraint used in the FSS decomposition was specifically targeted towards P100, it appears from the discrepancy results that the other peaks of the VEP were also extracted into the FSS source. More generally, this may or may not be the case depending whether the different EEG features are generated by the same neuronal pools. 122 C. Porcaro et al. / NeuroImage 50 (2010) 112–123 A related point is whether the reduction in inter- and intrasubject response variability that can be seen in Figs. 2, 5 and 6 would actually be detrimental to studies which rely on this variability, for example if the data were to be used for ST EEG-fMRI integration. Although a full examination of this question is beyond the scope of the current study, Figs. 5 and 6 clearly demonstrate that the ST variability observed in the FSS data is of the same order of magnitude as that seen in the data recorded outside of the scanner. On the other hand, the raw and ICA data recorded inside the scanner have considerably more ST variability. This suggests that the majority of the additional variability in the raw and ICA data is not neuronal in origin, since if it was it would also be seen in the data recorded outside of the scanner, but rather comes from residual artefacts which FSS successfully removes. Fig. 7 addresses the issue of whether the improvement in EEG data quality provided by FSS translates into an improved ability to integrate the EEG and fMRI data. Clearly, this is the primary goal of any combined EEG-fMRI experiment, and Fig. 7a demonstrates a much improved correlation between the area of the EEG and fMRI evoked responses when using the EEG data processed with FSS. Interestingly, no significant correlation between EEG and fMRI amplitudes was found for any method (Fig. 7b), although ST amplitudes are often used in EEG-informed fMRI analyses (Debener et al., 2005, 2006, Eichele et al., 2005, Bénar et al., 2007), highlighting the complexity of the relationship between the two modalities and suggesting that more work is needed to understand which parts of the EEG and fMRI signals are related and how this depends on task. Previous work has demonstrated that the ability to successfully correlate ST variability of the two modalities is improved by additional preprocessing not only of the EEG data but also of the fMRI data (Bagshaw and Warbrick 2007). Similarly to the way in which ICA has been used to denoise EEG data, it can also be used to improve the quality of ST estimates of the haemodynamic response function (Duann et al., 2002). It is also possible to apply the FSS approach to fMRI data, and future work will investigate whether it can lead to improvements in fMRI data quality as for the EEG data reported here. Conclusion The results demonstrate improved performance of FSS with respect to the original data and the data preprocessed by ICA. Continued efforts to improve the quality of EEG data recorded in the MRI scanner are important if EEG-fMRI is to realize its potential and provide a neuroimaging technique with excellent spatiotemporal resolution. Analysis of the type that has been performed in the current study has implications for the future development of EEG-fMRI integration techniques, for example where the emphasis on ST EEG analysis has been motivated by studies which have used ST data as a regressor in the fMRI analysis (Debener et al., 2006). To be fully exploited, this relies upon a thorough understanding of ST variability in both EEG and fMRI data sets (Bagshaw and Warbrick 2007). The results presented here show a consistent and significant improvement in the characterization of ST VEPs using FSS which is clearly the starting point for their use in any type of subsequent analysis. This improvement in the EEG was shadowed by an increased correlation with the fMRI data. The extent to which these results can be generalized to other EEG features of interest and whether improved EEG data quality leads to more successful EEG-fMRI integration are topics for further work. Acknowledgments This work was supported by grant number EP/F023057/1 from the UK Engineering and Physical Sciences Research Council (EPSRC). The authors would like to thank Nina Salman for help with data acquisition. References Allen, P.J., Polizzi, G., Krakow, K., Fish, D.R., Lemieux, L., 1998. Identification of EEG events in the MR scanner: the problem of pulse artifact and a method for its subtraction. Neuroimage 8, 229–239. Allen, P.J., Josephs, O., Turner, R., 2000. A method for removing imaging artifact from continuous EEG recorded during functional MRI. Neuroimage 12, 230–239. Bagshaw, A.P., Warbrick, T., 2007. Single trial variability of EEG and fMRI responses to visual stimuli. Neuroimage 38, 280–292. Barbati, G., Porcaro, C., Zappasodi, F., Rossini, P.M., Tecchio, F., 2004. Optimization of ICA approach for artifact identification and removal in MEG signals. Clin. Neurophysiol. 115, 1220–1232. Barbati, G., Sigismondi, R., Zappasodi, F., Porcaro, C., Graziadio, S., Valente, G., Balsi, M., Rossini, P.M., Tecchio, F., 2006. Functional source separation from magnetoencephalographic signals. Hum. Brain Map. 27, 925–934. Barbati, G., Porcaro, C., Hadjipapas, A., Adjamian, P., Pizzella, V., Romani, G.L., Seri, S., Tecchio, F., Barnes, G.R., 2008. Functional source separation applied to induced visual gamma activity. Hum. Brain Mapp. 29, 131–141. Bénar, C.-G., Schön, D., Grimault, S., Nazarian, B., Burle, B., Roth, M., Badier, J.-M., Marquis, P., Liegeois-Chauvel, C., Anton, J.-L., 2007. Single-trial analysis of oddball event-related potentials in simultaneous EEG-fMRI. Hum. Brain Mapp. 28, 602–613. Brainard, DH., 1997. The psychophysics toolbox. Spatial Vision 10, 433–436. Comon, P, 1994. Independent component analysis: a new concept? Signal Process 36, 287–314. Debener, S., Ullsperger, M., Siegel, M., Fiehler, K., von Cramon, D.Y., Engel, A.K., 2005. Trial-by-trial coupling of concurrent electroencephalogram and functional magnetic resonance imaging identified the dynamics of performance monitoring. J. Neurosci. 25, 11730–11737. Debener, S., Ullsperger, M., Siegel, M., Engel, A.K., 2006. Single-trial EEG-fMRI reveals the dynamics of cognitive function. Trends Cogn. Sci. 10, 558–563. Delorme, A., Makeig, S., 2004. EEGLAB: an open source toolbox for analysis of singletrial EEG dynamics including independent component analysis. J. Neurosci. Methods 134, 9–21. Duann, J.-R., Jung, T.-P., Kuo, W.-J., Yeh, T.-C., Makeig, S., Hsieh, J.-C., Sejnowski, T.J., 2002. Single-trial variability in event-related BOLD signals. Neuroimage 15, 823–835. Eichele, T., Specht, K., Moosmann, M., Jongsma, M.L.A., Quian Quiroga, R., Nordby, H., Hugdahl, K., 2005. Assessing the spatiotemporal evolution of neuronal activation with single-trial event-related potentials and functional MRI. Proc. Natl. Acad Sci. U.S.A. 49, 17798–17803. Goldman, R.I., Wei, C.-Y., Philiastides, M.G., Gerson, A.D., Friedman, D., Brown, T.R., Sajda, P., 2009. Single-trial discrimination for integrating simultaneous EEG and fMRI: identifying cortical areas contributing to trial-to-trial variability in the auditory oddball task. Neuroimage 47, 136–147. Hyvärinen, A., 1999. Fast and robust fixed point algorithms for independent component analysis. IEEE Trans. Neural. Netw. 10, 126–634. Hyvärinen, A., Karhunen, J., Oja, E, 2001. Independent Component Analysis. Wiley, New York. James, C.J., Hesse, C.W., 2005. Independent component analysis for biomedical signals. Physiol. Meas. 26, R15–R39. Kirkpatrick, S., Gelatt Jr., C.D., Vecchi, M.P., 1983. Optimization by simulated annealing. Science 220 (4598), 671–680. Laufs, H., Krakow, K., Sterzer, P., Eger, E., Beyerle, A., Salek-Haddadi, A., Kleinschmidt, A., 2003. Electroencephalographic signatures of attentional and cognitive default modes in spontaneous activity fluctuations at rest. Proc. Natl. Acad. Sci. U.S.A. 100 (19), 11053–11058. Mantini, D., Perrucci, M.G., Cugini, S., Ferretti, A., Romani, G.L., Del Gratta, C., 2007. Complete artifact removal for EEG recorded during continuous fMRI using independent component analysis. Neuroimage 34 (2), 598–607. Makeig, S., Westerfield, M., Jung, T.P., Enghoff, S., Townsend, J., Courchesne, E., Sejnowski, T.J., 2002. Dynamic brain sources of visual evoked responses. Science 295, 690–693. Makeig, S., Debener, S., Onton, J., Delorme, A., 2004a. Mining event-related brain dynamics. Trends Cogn. Sci. 8, 204–210. Makeig, S., Delorme, A., Westerfield, M., Jung, T.P., Townsend, J., Courchesne, E., Sejnowski, T.J., 2004b. Electroencephalographic brain dynamics following manually responded visual targets. PLOS Biol. 6, 0747–0762. Medaglia, M.T., Tecchio, F., Seri, S., Di Lorenzo, G., Rossini, P.M., Porcaro, C., 2009. Contradiction in universal and particular reasoning. Hum. Brain Mapp. 30, 4187–4197. Mulert, C., Seifert, C., Leicht, G., Kirsch, V., Ertl, M., Karch, S., Moosmann, M., Lutz, J., Möller, H.-J., Hegerl, U., Pogarell, O., Jäger, L., 2008. Single-trial coupling of EEG and fMRI reveals the involvement of early anterior cingulated cortex activation in effortful decision making. Neuroimage 42, 158–168. Niazy, R.K., Beckmann, C.F., Iannetti, G.D., Brady, J.M., Smith, S.M., 2005. Removal of FMRI environment artifacts from EEG data using optimal basis sets. Neuroimage 28, 720–737. Olbrich, S., Mulert, C., Karch, S., Trenner, M., Leicht, G., Pogarell, O., Hegerl, U., 2009. EEG-vigilance and BOLD effect during simultaneous EEG/fMRI measurement. Neuroimage 45, 319–332. Ostwald, D., Porcaro, C., Bagshaw, A.P., 2010. An information theoretic approach to EEGfMRI integration of visually evoked responses. Neuroimage 49, 498–516. Pascual-Marqui, R.D., 2002. Standardized low-resolution brain electromagnetic tomography (sLORETA): technical details. Methods Find. Exp. Clin. Pharmacol. 24D, 5–12. Pelli, DG., 1997. The VideoToolbox software for visual psychophysics: transforming numbers into movies. Spatial Vision 10, 437–442. C. Porcaro et al. / NeuroImage 50 (2010) 112–123 Porcaro, C., Coppola, G., Di Lorenzo, G., Zappasodi, F., Siracusano, A., Pierelli, F., Rossini, P.M., Tecchio, F., Seri, S., 2009. Hand somatosensory subcortical and cortical sources assessed by functional source separation: an EEG study. Hum. Brain Mapp 30, 660–674. Porcaro, C., Barbati, G., Zappasodi, F., Rossini, P.M., Tecchio, F., 2008. Hand sensory– motor cortical network assessed by functional source separation. Hum. Brain Mapp. 29, 70–81. 123 Tecchio, F., Porcaro, C., Barbati, G., Zappasodi, F., 2007. Functional source separation and hand cortical representation for brain–computer interface feature extraction. J. Physiol. 580, 703–721. Wan, X., Riera, J., Iwata, K., Takahashi, M., Wakabayashi, T., Kawashima, R., 2006. The neural basis of the hemodynamic response nonlinearity in human primary visual cortex: implications for neurovascular coupling mechanism. Neuroimage 32, 616–625.