Single trial variability of EEG and fMRI responses to visual stimuli

www.elsevier.com/locate/ynimg NeuroImage 38 (2007) 280 – 292 Single trial variability of EEG and fMRI responses to visual stimuli Andrew P. Bagshaw⁎ and Tracy Warbrick School of Psychology and Birmingham University Imaging Centre (BUIC), University of Birmingham, Edgbaston, Birmingham B15 2TT, UK Received 24 January 2007; revised 26 June 2007; accepted 16 July 2007 Available online 11 August 2007 Recent EEG–fMRI studies have suggested a novel method of data fusion which uses single trial (ST) estimates of event-related potentials in the fMRI analysis. This is potentially very powerful, but rests on the assumption that the ST variability observed in EEG is reflected in the fMRI signal. The current study investigated this assumption and compared two different data processing strategies for each modality. Five subjects underwent separate EEG and fMRI sessions with checkerboard stimuli at two contrasts. EEG data were preprocessed using wavelet denoising and independent component analysis (ICA), whilst the general linear model and ICA were used for fMRI. Amplitudes and latencies of the P1 and N2 components of the visual evoked potential (VEP) were calculated for each trial. For fMRI, the amplitudes and latencies of the ST haemodynamic responses (HR) were calculated. Within modality, the results for the two processing methods were significantly correlated in the majority of data sets. Across modality, the average amplitudes of the VEPs and HRs were also significantly correlated. Examination of ST variability demonstrated that the amplitudes of the mean VEPs and HRs are both influenced by the latency variability of the ST responses to a greater extent than the amplitude variability. For high contrast stimuli the latency variability in EEG and fMRI was significantly correlated, with a similar trend seen for the low contrast stimuli. The results confirm the validity of examining both the EEG and fMRI signals on an ST basis and suggest an underlying neuronal origin in both modalities. © 2007 Elsevier Inc. All rights reserved. Introduction The combination of electroencephalography (EEG) and functional magnetic resonance imaging (fMRI) holds the promise of millisecond temporal resolution and millimetre spatial resolution and thus a more comprehensive view of cortical processing than can be provided by either method alone. However, the question of how best to fuse EEG and fMRI data sets has not been fully resolved. Kilner et al. (2005) described three methods of data fusion: through prediction (using a property of the EEG signal to ⁎ Corresponding author. Fax: +44 121 414 4897. E-mail address: a.p.bagshaw@bham.ac.uk (A.P. Bagshaw). Available online on ScienceDirect (www.sciencedirect.com). 1053-8119/$ - see front matter © 2007 Elsevier Inc. All rights reserved. doi:10.1016/j.neuroimage.2007.07.042 predict changes in the BOLD response); through constraints (constraining EEG source localisation based on fMRI results); and through common forward models (a comprehensive model or set of models which describes the mechanisms generating the two signals and the links between them). The successful development of any of these data fusion techniques requires a much more detailed knowledge of their correspondence than is presently available. Recent work has suggested a novel method of fusion through prediction. Unlike previous studies which examined the EEG signal averaged over a number of trials (e.g. Arthurs et al., 2000; Bledowski et al., 2006; Janz et al., 2001; Mulert et al., 2005; Singh et al., 2003), Debener et al. (2005) and Eichele et al. (2005) used single trial (ST) estimates of the amplitude of a given event-related potential (ERP) peak as regressors in the fMRI analysis. This is potentially a very powerful method of data fusion, but rests on the implicit assumption that the ST EEG variability is reflected in the BOLD signal. This assumption seems valid given the ability of the ST EEG regressors to generate meaningful fMRI statistical maps. However, further investigation is required, especially since the widespread use of averaging, both of EEG and fMRI data, is based on the opposite assumption: namely, that variability from trial to trial is a result of processes which are not relevant to the task. Several EEG and MEG studies have questioned this view and demonstrated that a considerable amount of information is lost by averaging (Liu and Ioannides, 1996; Jung et al., 2001; Laskaris et al., 2003; Makeig et al., 2002, 2004). A recent fMRI study has also highlighted the benefit of examining the haemodynamic response (HR) on a trial by trial basis to understand the coupling between the stimulus evoked response and ongoing brain processes (Fox et al., 2006). Examination of ST data requires preprocessing. Two main methods have been used for EEG: wavelet denoising (WD, Quian Quiroga and Garcia, 2003; Eichele et al., 2005) and independent component analysis (ICA, Makeig et al., 2002, 2004; Debener et al., 2005). Similarly, fMRI data have been widely analysed using the general linear model (GLM) and ICA. In terms of ST analysis, the different methods provide time courses from which ST parameters can be extracted (amplitudes, latencies, etc.). The current study acquired EEG and fMRI data in response to the same visual stimulus in separate sessions. Two analysis methods were applied to each imaging modality (WD/ICA for EEG data, GLM/ A.P. Bagshaw, T. Warbrick / NeuroImage 38 (2007) 280–292 281 ICA for fMRI data) in order to examine the concordance between the ST parameters thereby generated. For each modality, one of the methods was selected for further analysis and a comparison performed between EEG and fMRI data. The purpose of the study was two fold: (1) to compare different analysis strategies for ST EEG and fMRI analysis and (2) to investigate the extent to which the trial by trial variability observed in EEG is reflected in fMRI. Methods Experimental protocol Five subjects (1 female, ages 26–33) were paid for their participation. Written informed consent was obtained and the protocol was approved by the Research Ethics Board of the Birmingham University Imaging Centre. The experimental protocol consisted of separate sessions of EEG and fMRI recording with an identical full field checkerboard stimulus reversing at 2 Hz (i.e. the checkerboard was presented for 500 ms before reversing). This frequency allowed individual, non-overlapping visual evoked potentials (VEPs) to be recorded following each reversal. The timing of stimulus presentation was optimised for each modality (see below). Presentation (Neurobehavioral Systems, Inc., CA, USA) was used to display the stimuli. The EEG session was conducted first in all subjects using a 128 channel BioSemi Active Two EEG system (BioSemi, Amsterdam, Netherlands), with electrodes placed in a nylon cap according to the 10–5 system (Oostenveld and Praamstra, 2001). Following application of the electrodes, subjects were seated in a comfortable position approximately 80 cm from a computer monitor. Checkerboards reversing at 2 Hz were presented in blocks of 5 s with a ten-second gap between successive blocks. Participants were instructed to fixate on a central red dot which was present between checkerboard stimuli and to minimise blinking whilst the checkerboard was presented. Each run consisted of twenty blocks, and took approximately 5 min. Alternate runs consisted of high contrast (black and white) and low contrast (two shades of grey) checkerboards, with three runs of each acquired. In total, 600 VEPs were collected for each contrast (10 VEPS per block, 20 blocks per run, 3 runs per contrast). EEG data were sampled at 2048 Hz with a linked mastoid reference. The acquisition of two stimulus contrasts allowed the differentiation of the analysis methods based on the comparison of data for which there was an a priori expectation of a relationship. Both EP and HR magnitudes are dependent on contrast (Logothetis et al., 2001; Odom et al., 2004), and it would therefore be expected that measures derived from the two contrasts would be correlated (Fig. 1). For example, selection of the analysis method to be compared across modality was based on the extent to which the data for high and low contrast stimuli were correlated. Following these six blocks of EEG, the electrodes but not the electrode cap were removed. The BioSemi active electrodes are not MR compatible since they contain preamplifier electronics. The participant was taken immediately to the MRI scanner (3 T Achieva, Philips Medical Systems, Best, Netherlands), where a standard T1-weighted anatomical scan was performed with the electrode cap still in place. Following this, the participant was taken out of the scanner, the electrode cap removed and the participant replaced in the scanner. By acquiring an anatomical scan with the electrode cap in place accurate electrode positions could be recovered for subsequent EEG source localisation (data on the spatial aspect of the comparison between EEG and fMRI are not reported here). The cap was removed before fMRI scanning for comfort. fMRI data were acquired using an eight channel phased array head coil. Subjects viewed a screen in the rear of the magnet bore via a mirror placed on the head coil (eye to screen distance approximately 65 cm). High and low contrast checkerboards were presented in alternate runs. To allow detection of HRs to individual stimuli without overlap, the 2 Hz checkerboard was presented for 1 s with 30 s between stimuli. Two to four runs of each stimulus contrast were acquired per participant, depending on the willingness of the subject to remain in the scanner (subject A had 2 runs for each stimulus, subject B three runs, and subjects C–E 4 runs). Ten stimuli were presented per run, which lasted just over 5 min. To encourage fixation and maintain attention, subjects were instructed to press a button whenever the central fixation dot changed color. fMRI data were acquired from 8 slices oriented parallel to and centred on the calcarine sulcus with a TR of 500 ms (voxel size 2.5 × 2.5 × 3 mm, TE = 35 ms, flip angle 50°, SENSE factor of 2 (Pruessman et al., 1999)). The short TR was chosen to allow good characterisation of the HR. For registration purposes a whole head EPI image with the same matrix size and slice orientation was collected. In addition, one run was collected with a 30-second block design consisting of four repeats of the high contrast stimulus. The whole experiment, from the start of the application of EEG electrodes to the subject being removed from the scanner after all fMRI runs, lasted around 4 h per subject. EEG data analysis Epochs were extracted from the continuous EEG data (50 ms pre- to 400 ms post-stimulus) using BrainVision Analyzer (BrainProducts, Munich, Germany). Data were band pass filtered from 0.1 to 70 Hz and baseline corrected using the first 50 ms of the epoch, and epochs were rejected in which an absolute deviation of 100 μV was observed. Two components of the VEP were examined: a positive deviation peaking at approximately 100 ms (denoted P1) and a later negative deviation peaking at around 140 ms (denoted N2). Some studies (e.g. Kashikura et al., 2001; Di Russo et al., 2001) have labelled the negative of these two peaks N1, but in this study it is labelled N2 to distinguish it from an earlier negative peak at around 70 ms (Bonmasser et al., 2001; Laskaris et al., 2003, labelled C1 in Di Russo et al., 2001). The peak at 70 ms was not reliably observed in the current data, probably because a full field stimulus was used (Di Russo et al., 2001; Ohla et al., 2005) and is not considered further. Epochs were further processed for ST analysis in two ways: with wavelet denoising (WD) and independent component analysis (ICA). For WD, data were taken from a single electrode which was chosen individually based on the topography of the evoked potential (EP) response at the peak of the P1. Considerable variability was noted in the EP topography between subjects, suggesting that the use of a single standard electrode, for example O1 or O2, would not be optimal. The selected electrodes for subjects A–E were Oz, PPO6h, PO3, O1 and 282 A.P. Bagshaw, T. Warbrick / NeuroImage 38 (2007) 280–292 Fig. 1. Several comparisons between EEG data preprocessed with WD and ICA: (a) mean P1–N2 amplitude for high contrast vs. low contrast stimuli, (b) mean P1–N2 amplitude for WD vs. ICA, (c) Trial amplitudes for high contrast vs. low contrast stimuli and (d) ST amplitude variability for high contrast vs. low contrast stimuli. POO1 for the high contrast stimulus, and PPO2h, POO2, POO1, O1 and POO2 for the low contrast stimuli respectively. The nomenclature is according to the 10–5 electrode system (Oostenveld and Praamstra, 2001). WD was accomplished for each of the runs separately using EP_den_v2 (http://www.vis.caltech.edu/~rodri/EP_den/EP_den_ home.htm). Following the method of Quian Quiroga and Garcia (2003), the data were decomposed into five wavelet scales, and wavelet coefficients chosen to optimise the average signal over the first 250 ms post-stimulus. The denoising scheme optimised for the average EP was then applied to the ST data and the resulting denoised trials saved for further analysis. For ICA, data from the full 128 electrodes were exported from BrainVision Analyzer into EEGLAB (http://sccn.ucsd.edu/eeglab/, Delorme and Makeig, 2004). As for WD, ICA was applied to each run of the EEG data separately. Data were rereferenced to an average reference and subjected to the ICA algorithm runICA. The results of this were saved for further analysis. Analysis of ICA data can proceed via two basic routes: selection of task-related components or rejection of non-taskrelated components. In the current study, task-related components were extracted in order to minimise the variance contribution from non-neuronal sources. The danger is that not all components which have a neuronal origin will be selected, leading to an underestimation of the amount of inter-trial neuronal variability. Components were selected which survived a dual criterion of spatial and temporal correlation with the average EP for each run (Lee et al., 2003). The correlation threshold was set to 0.6, that is, components were selected whose time course had a correlation of greater than 0.6 with the average EP and whose spatial map had a correlation greater than 0.6 with the spatial map of the average EP. This resulted in 3–13 components per run. These components were then projected to the same channel as selected for WD analysis, resulting in the data that were analysed further (that is, to recover the data with correct polarity and in units of microvolts, the product of the ICA inverse weight matrix and the A.P. Bagshaw, T. Warbrick / NeuroImage 38 (2007) 280–292 283 component activation matrix were calculated, with the rows of the activation matrix for non-selected components set to zero; see Jung et al., 2001 for further details). These two preprocessing methods resulted in two sets of 600 trials for each contrast which subsequently underwent the same analysis in Matlab (MathWorks, USA). For each trial, the amplitudes of P1 and N2 were taken as the maximum and minimum values within a time window of 50 ms centred on the time of the respective peaks in the average EP. The ST latencies were taken as the times of these responses. Two measures of the EP amplitude were calculated for both average and ST data. The first was the peak-to-peak amplitude of the P1 and N2 components (‘P1–N2 amplitude’). This relative measure of two peaks is widely used for EEG and MEG data and provides a robust measure of the evoked response (e.g. Arthurs et al., 2000; Debener et al., 2005). In addition, the root mean squared (RMS) value of each trial within a time window from 50 ms to 250 ms post-stimulus was calculated as a measure of the overall stimulus-related activity (expressed as ‘Trial amplitude’, similar to the measure used by Ploner et al., 2006 and Koch et al., 2006). As a measure of the performance of the preprocessing methods in allowing visualisation of ST EP components, the number of trials which demonstrated a response was determined (i.e. the number of trials which had an ST P1 or N2 above the level of the background noise). Following the method of Sander et al. (2005), a histogram was created of the ST P1 and N2 amplitudes, fitted with a Gaussian profile and thresholded at the level of the average noise. The noise level was calculated as the RMS value from 0 to 25 ms, calculated on a ST basis and averaged over all trials. fMRI data analysis fMRI data were analysed with the GLM as implemented in FSL (FMRIB, Oxford, UK). Data were smoothed with a 3 mm kernel, motion and slice-timing corrected and analysed using a gamma HRF plus temporal derivative. Statistical maps were thresholded using clusters determined by Z N 2.3 and a (corrected) cluster significance threshold of p = 0.05 (Worsley et al., 1992). For each run separately, the time course of all activated voxels was extracted and split into epochs of 30 s. Each epoch therefore contained the HR to a single 1 s stimulus. In addition to the GLM analysis, fMRI data underwent ICA using MELODIC, also part of the FSL package. The same preprocessing steps were applied as for the GLM analysis. For all components, voxels with an estimated posterior probability of activation greater than a threshold of 0.5 were designated as activated (Beckmann and Smith, 2004). As for the EEG data, task-related components were selected. In this case, they were selected based on the Percentage Overlap Ratio (POR, Duann et al., 2002) of the spatial map with that of the main task-related component from the 30 s block design. The POR measures the number of commonly activated voxels as a percentage of the total number of activated voxels. The main taskrelated component of the block design data was identified as the one with the maximum correlation with the expected BOLD response given the task design. This two stage process allowed components to be selected independently of their time course. In the majority of cases, a single task-related component was found. In some runs two were identified, in which case the component with the highest POR was selected. As for the GLM analysis, a time course was extracted and split into 30 s epochs. The amplitude of the ST haemodynamic response was defined as the maximum in the first 12 s post-stimulus. Since standard HRFs peak after 5–6 s (Glover, 1999) this allows for a considerable amount of variation whilst ensuring that the primary response is identified. The latency of the ST response was taken as the time point of the maximum response. fMRI data were not fitted or interpolated, meaning that ST latencies were available with a temporal resolution of 500 ms. Intra- and inter-modality comparison As discussed above, EEG data were preprocessed with WD and ICA, whilst the GLM and ICA were used for fMRI data. These analyses were compared within modality (that is, for EEG, comparisons were made between WD and ICA, whilst for fMRI, the results from GLM vs. ICA were contrasted). Based on the intramodality results and in order to reduce the number of inter-modality comparisons, data from one of the preprocessing methods were selected and comparisons were made across modalities. Ratios (high/ low contrast) were calculated for the EEG data to standardise the individual measurements and allow comparisons across subjects. While this will reduce effects due to individual anatomy, geometry, source location etc., it has the potential to mask some effects and to be sensitive to noise since the responses to the low contrast stimuli were generally of lower amplitude. Both EEG and fMRI analyses followed standard practices, meaning that EEG WD amplitudes were expressed in microvolts whilst fMRI ICA data were variance normalised and expressed in normalised units (Beckmann and Smith, 2004). This means that fMRI data have undergone a normalisation procedure to take account of inter-subject differences and are not expressed as a ratio between the responses to the two contrasts. The results are presented in four steps. Firstly, within modality, several comparisons are made of the responses generated using the two preprocessing methods. Data from the most effective processing method are then selected (WD for EEG, ICA for fMRI), and more detailed properties of the single trial data are extracted, again within modality. To reinforce the validity of the data processing, average EEG and fMRI signals are then presented across modality to confirm their correlation, as expected from previous studies (e.g. Arthurs et al., 2000; Mulert et al., 2005). Finally, the ST properties (amplitude and latency variability) are compared across modality to investigate the extent to which ST variations in EEG data are replicated in fMRI data. Results EEG: comparison of WD and ICA analysis Fig. 1 plots a number of comparisons between average EEG data preprocessed with WD and ICA. In this and other similar Table 1 Number of trials for which P1 and N2 amplitude was above the noise level (out of a total of 600 for each contrast) Subject WD P1 High A B C D E 514 403 423 459 545 Low 376 459 496 503 515 N2 High 431 364 332 415 476 Low 317 312 411 402 442 ICA P1 High 457 412 462 452 547 Low 396 426 520 451 487 N2 High 323 393 422 347 463 Low 357 414 300 353 442 284 A.P. Bagshaw, T. Warbrick / NeuroImage 38 (2007) 280–292 demonstrates a significant correlation between high and low contrast Trial amplitudes for the WD (R2 = 0.86, p b 0.05, two tailed) but not the ICA data, in agreement with the P1–N2 data in Fig. 1a. The P1–N2 amplitude was also assessed on an ST basis, allowing its variation to be calculated (variances for P1 and N2 amplitudes over all trials were added to calculate the variance of the P1–N2 amplitude, which was then square rooted to get the pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi standard deviation i.e. rP1ÀN2 ¼ r2 þ r2 ). This is plotted in P1 N2 Fig. 1d, which shows a significant correlation (R2 = 0.97, p b 0.01, two tailed) for the WD data, but not for the ICA data. Table 1 contains the number of trials which showed a P1 and N2 amplitude above the noise threshold using the method of Sander et al. (2005). There was no significant difference between the numbers of trials identified for WD and ICA data (two tailed t test, p N 0.4 in all cases). An example of the results for the two preprocessing methods is shown in Fig. 2. For subject E, the 600 trials acquired with the high contrast stimulus are plotted with the amplitude colour coded (note the different amplitude scales for WD and ICA data: data are scaled symmetrically to the maximum response). No inter-trial smoothing has been applied. Both P1 and N2 can be observed, although they appear to be clearer in the WD data, consistent with the analysis discussed above. In Fig. 2c the average EPs generated by the two methods are plotted. Although there is a difference in the amplitude of the WD and ICA data, as demonstrated in Figs. 2a and b, the overall shape of the mean EP is very similar, particularly for the P1–N2 response. Table 2 contains Pearson correlation coefficients for the comparison between ST P1–N2 amplitudes calculated with WD and ICA. Correlations were highly significant for all subjects and contrasts (p b 0.001, two tailed, except for subject C, low contrast stimulus, p b 0.05). fMRI: comparison of GLM and ICA analysis In subjects C and D, two runs and one run respectively of the low contrast stimulus did not lead to significantly activated voxels in the GLM analysis. In addition, ICA analysis of two low contrast runs from subject C did not identify any components with Table 2 Pearson correlation coefficients for EEG (P1–N2 amplitude, ICA and WD) and fMRI (HR amplitude, GLM and ICA) data Subject EEG (N = 600) High Fig. 2. For subject E, the ST EEG data from the 600 trials acquired with the high contrast stimulus are plotted with the amplitude colour coded following preprocessing with (a) WD and (b) ICA. In panel c the average EPs are plotted with normalised amplitude. Both WD and ICA allow visualisation of ST responses, whilst the average data are very similar. A B C D E 0.24 b 0.20 b 0.23 b 0.15 b 0.44 b Low 0.37 b 0.17 b 0.10 d 0.41 b 0.17 b fMRI (N a) High 0.63 c 0.53 c 0.14 e 0.68 c 0.87 c (20) (30) (40) (40) (40) Low 0.78 c (20) 0.64 c (30) 0.77 c (10) 0.67 c (30) 0.25 e (40) figures, each point represents the data from a single subject. In Fig. 1a, P1–N2 amplitude is plotted for high contrast vs. low contrast stimuli, for both WD and ICA data. There was a significant correlation for the WD data (R2 = 0.85, p b 0.05, two tailed), but not the ICA data (R2 = 0.23). Fig. 1b demonstrates that the WD and ICA amplitudes were correlated for the high (R2 = 0.93, p b 0.01, two tailed) but not the low contrast data (R2 = 0.68), whilst Fig. 1c a N is given in brackets for each subject and contrast. It varies between subjects because of the number of runs of data acquired, plus the effect of some runs not resulting in activated voxels. Data were only compared when estimates of the ST HR were available for both preprocessing methods, which required that for a given run both methods resulted in significantly activated voxels. b Significant at p b 0.001. c Significant at p b 0.01. d Significant at p b 0.05. e Non-significant (p N 0.1 in both cases). A.P. Bagshaw, T. Warbrick / NeuroImage 38 (2007) 280–292 Fig. 3. For subject B, the ST fMRI data are plotted for high and low contrast stimuli and preprocessing via the GLM and ICA. The right hand panels show the average HRs for the two methods. For both ST and average HRs the main features of the data are replicated for the two preprocessing methods. Positive values are rendered on a red colour scale, whilst negative changes are on a blue colour scale. 285 286 A.P. Bagshaw, T. Warbrick / NeuroImage 38 (2007) 280–292 appreciable POR (POR b 30%). In these cases, the data from the remaining runs are reported. ST and average HRs for the two processing strands are compared for subject B in Fig. 3. The ST data are plotted in the same way as in Fig. 2, with the amplitude colour coded. In this case, data are normalised to the maximum positive response (i.e. the maximum response has a value of unity in all cases). Generally, and as seen in this subject, the GLM and ICA results were qualitatively similar, with the overall pattern of the responses replicated in both analyses (see, for example, the top panel which shows increased latency in response to the high contrast stimulus in trials 15–20 for both GLM and ICA data). This was reflected in the average HRs, which were also extremely similar (Fig. 3, right hand panels). Fig. 4a demonstrates the tendency in both data sets for the subjects with the highest amplitude high contrast stimuli HRs to be those with the highest amplitude low contrast stimuli HRs, although certainly for the ICA data this relationship is nonlinear. As can be seen from Fig. 4b, a significant correlation was found between GLM and ICA mean amplitudes (R2 = 0.82, p b 0.05, two tailed) for high but not low contrast (R2 = 0.71) stimuli. Fig. 4c demonstrates a significant relationship (R2 = 0.89, p b 0.05, two tailed) between the latency variability in response to high and low contrast stimuli for ICA data, but not when using the GLM. Latency variability from ICA and GLM data was also correlated (R2 = 0.99, p b 0.001, two tailed) for high contrast stimuli but not for low contrast stimuli (Fig. 4d). Data from the low contrast stimuli were generally lower in amplitude, and there was less agreement between the GLM and ICA analyses as seen in Figs. 4b and d. Table 2 contains Pearson correlation coefficients for the comparison between ST HR amplitudes calculated with the GLM and ICA. Correlations were significant (p b 0.01, two tailed) for all but two of the fMRI data sets (subjects C and E, high and low contrast stimuli respectively, correlations non-significant (p N 0.1)). EEG and fMRI: ST analyses within modality Fig. 5 demonstrates that the amplitude of the P1 component of the mean EP is proportional to the latency variability (R2 = 0.79, Fig. 4. Several comparisons between fMRI data preprocessed with WD and ICA: (a) mean HR amplitude for high contrast vs. low contrast stimuli, (b) mean HR amplitude for GLM vs. ICA, (c) ST latency variability for high contrast vs. low contrast stimuli and (d) ST latency variability for GLM vs. ICA. A.P. Bagshaw, T. Warbrick / NeuroImage 38 (2007) 280–292 287 p b 0.05, two tailed) but not the amplitude variability of the ST data. The data are plotted for P1, rather than the combined P1– N2 values that were presented in Fig. 1, because of the ambiguity in defining latency for the combined peaks. Fig. 5c plots similar data for the Trial amplitude (RMS value from 50 to 250 ms poststimulus). A significant positive correlation (R2 = 0.87, p b 0.05, two tailed) was observed between the mean value of the Trial amplitude and the standard deviation on an ST basis. Similar plots for the fMRI data are presented in Fig. 6. Fig. 6b demonstrates that the amplitude of the mean HR was inversely proportional to the latency variability for both stimuli (Fig. 6b, R2 = 0.79 and 0.80 for high and low contrasts, p b 0.05, two tailed in both cases; for clarity the linear fits have not been plotted). In addition, there is a tendency for subjects with higher amplitude HRs to have less ST amplitude variability (Fig. 6a). EEG–fMRI: comparison of average signals Figs. 5 and 6, as well as Figs. 1 and 4, lend support to the view that the ST parameters being extracted from the data were meaningful and not dominated by noise since if the extracted parameters were sampled from noise significant correlations between variables within modality would not be expected. The fact that such correlations could be demonstrated was reassuring and encouraged the view that performing inter-modality comparisons was meaningful. Fig. 7 compares the average values across modality. Fig. 7a demonstrates that the P1–N2 amplitude ratio is strongly correlated (p b 0.01, two tailed) with the amplitude of the mean HR for both high and low contrast stimuli. Figs. 7b and c plot similar data, but include the estimate of mean EP P1–N2 amplitude provided by the ST data. The difference between this measure (labelled ‘ST data’) and the amplitude derived from the mean EP is that it includes signal strength that is not precisely phase-locked from trial to trial. These two measures of P1–N2 amplitude are significantly correlated (data not shown, R2 = 0.84, p b 0.05, two tailed), as would be expected. The P1–N2 amplitude derived from the mean EP was more strongly correlated with the mean HR amplitude for both high and low contrast stimuli (high contrast: R2 = 0.94 and 0.85 for mean EP and ST data (significant at p b 0.01 and p b 0.05); low contrast: R 2 = 0.93 and 0.73 (significant at p b 0.01 and non-significant p N 0.05)). Linear fits are not shown in Figs. 7b and c for clarity. Fig. 7d demonstrates that the correlation between P1–N2 and HR amplitudes was specific: no such relationship is observed between the Trial amplitude and the HR. EEG–fMRI: comparison of ST signals Fig. 7a demonstrates that the average P1–N2 amplitude is strongly correlated with the amplitude of the mean HR. However, Fig. 8a shows that the variability of the P1–N2 amplitude on an ST basis is not significantly correlated with the variability of the HR amplitude. On the other hand, there is a relationship between the latency variability of P1 and that of the HR (Fig. 8b). This relationship is significant for the high contrast data (linear fit not shown for clarity, R2 = 0.81, p b 0.05, two tailed). Although it is not significant for the low contrast data (R2 = 0.54), the same tendency is seen for subjects with Fig. 5. The amplitude of the P1 component of the mean EP (a) is not related to the amplitude variability of the ST data but (b) is related to the latency variability of the ST data. (c) The mean Trial amplitude is proportional to the ST amplitude variability. 288 A.P. Bagshaw, T. Warbrick / NeuroImage 38 (2007) 280–292 Fig. 6. The amplitude of the mean HR (a) is not related to the amplitude variability but (b) is related to the latency variability of the ST data. more variability in the P1 latency to have more variability in the timing of the HR. Discussion This study investigated the extent to which useful information can be extracted from trial by trial variability using EEG and fMRI responses to a robust visual task. This is an important question in light of the expanding use of EEG–fMRI in neuroscience, and particularly given the potential of recent analysis methods which make use of ST ERP amplitudes as regressors in the fMRI analysis (Debener et al., 2005, 2006; Eichele et al., 2005). The study had two main purposes. The first was to compare commonly applied analysis strategies for EEG and fMRI data since, in order to extract useful information from inherently noisy ST data, some preprocessing is required. For EEG data, two methods have been used quite widely: wavelet denoising (WD) and independent component analysis (ICA). The current study applied both methods and demonstrated highly correlated results (Table 2) and that both can lead to improved visualisation of the ST waveforms (Fig. 2). In general, WD and ICA gave interpretable results with WD giving improved correlations between data for high and low contrast stimuli (Fig. 1). However, although the WD and ICA amplitudes were significantly correlated, the correlation coefficients were relatively low (Table 2), and the extracted values were not identical, as can readily be appreciated from the data presented in Fig. 2. The effect this would have on ST EEG–fMRI analysis of the type proposed by Debener et al. (2005) and Eichele et al. (2005) requires further investigation in order to optimise and generalise the technique. One explanation for the relatively low correlations between WD and ICA preprocessing of the EEG data, and between GLM and ICA preprocessing of the fMRI data, stems from the way in which independent components were selected. In this study, task-related components were selected, rather than the alternative which is to remove non-task-related components. This is quite a strong constraint and is likely to have a much greater effect on the data than WD or GLM analysis. This approach was followed in order to focus on variability of neuronal origin. For example, fMRI data sets contain variability related to respiratory artefacts (Thomason et al., 2007), whilst in EEG artefacts related to blinking and muscle activity can be observed (Jung et al., 2001). However, this component selection means that the ICA data analysed on a ST basis conform to quite a narrow definition of what is task-related. This has the advantage that other sources of artefacts or variability that were not directly related to the task (e.g. background fluctuations) were excluded, but it potentially loses information and underestimates the true level of the variance. It is also responsible for the differences in the magnitudes of the ST responses that were obtained with different preprocessing methods (see, for example, Fig. 2 for the EEG data). For example, the EEG component selection will reduce the contribution from non-phaselocked activity because the average time course of the component had to correlate with that of the mean EP. (Note however that this does not affect the following discussion concerning the relative contribution of evoked and induced responses to the BOLD signal since that was predicated on results from the WD data.) A more conservative approach, whereby obviously artefactual components are removed, may be more appropriate, particularly for EEG data (Jung et al., 2001). Alternatively, Iyer and Zouridakis (2007) suggest component selection via an iterative ICA method which seems to outperform WD. Despite the complications of component selection when using ICA for noise reduction, an issue that is not unique to the current study, it is clear that the main characteristics of both the mean EPs and mean HRs were maintained following ICA (Figs. 2 and 3) and that visualisation of ST responses was improved, particularly for fMRI. For fMRI data, the GLM and ICA were compared. As for EEG data, results were highly correlated (Table 2) as well as being qualitatively similar. ICA had some advantages particularly for low contrast stimuli, apparently conserving physiological variability whilst reducing noise. More runs had to be excluded from the GLM analysis because of a lack of significantly activated voxels, indicating an advantage to the use of ICA when the signal to noise ratio is low. In addition, for examination of time course data, ICA has the advantage that it places no constraints on the form of the haemodynamic response (HR). However, with both methods analysis of the HR on a trial by trial basis was feasible, even A.P. Bagshaw, T. Warbrick / NeuroImage 38 (2007) 280–292 289 Fig. 7. Cross modality comparisons of amplitude and latency data: (a) the P1–N2 amplitude ratio is strongly correlated with the amplitude of the mean HR for both high and low contrast stimuli, panels b and c show comparisons between the P1 and N2 amplitude derived from the mean EP and the ST data and the mean HR amplitude for high and low contrast stimuli, (d) no correlation was observed between the Trial amplitude and the HR. with the short stimulus that was used. In general, the methods applied to both EEG and fMRI data allowed ST parameters to be extracted which resulted in meaningful results within modality (Figs. 5 and 6). The second purpose of this study was to investigate the extent to which the ST variability observed in EEG data is replicated in fMRI data. Analysis within modality demonstrated the validity of the parameters extracted from the ST data, and comparison of the average EEG and fMRI responses confirmed that there was a strong relationship between the data sets (Fig. 7), in agreement with a number of previous studies (Arthurs et al., 2000; Di Russo et al., 2001; Janz et al., 2001; Singh et al., 2003; Henning et al., 2005). Analysis at the level of the ST revealed that the variability in the latency of P1 was related to the variability in the latency of the HR, although this was only significant for the high contrast stimulus. No significant correlations were observed between the amplitude variability observed in EEG and fMRI data. The practical consequences of this are unclear, and efforts are currently underway to examine this issue in simultaneously acquired EEG– fMRI. This will allow a more direct comparison of the trial by trial variability, providing electrical and haemodynamic responses to each stimulus and thus a more rigorous examination of the relationship between them. A number of interesting results have arisen from the current work. Figs. 7b and c suggest that the amplitude of the BOLD response is related to the phase-locked component of the EEG signal (i.e. the evoked rather than induced response). An improved correlation is found with the amplitude of the mean EP compared to the mean amplitude of the ST responses: these measures differ in that the latter includes ST responses within a 50 ms window centred on the overall peak. Since previous studies have not performed an ST analysis, they have by default focussed exclusively on the phase-locked component of the signal. There is also little theoretical work into the relative contributions of phase and non-phase-locked EEG responses to the BOLD signal. Robinson et al. (2006) predicted that the magnitude of the BOLD 290 A.P. Bagshaw, T. Warbrick / NeuroImage 38 (2007) 280–292 Fig. 8. (a) The variability of the P1–N2 amplitude is not significantly correlated with the variability of the HR amplitude but (b) the latency variability of P1 is correlated with that of the HR. response should be proportional to the time integrated neuronal activity, rather than its peak value. This may suggest a dependence of the BOLD response on the non-phase-locked component of the ERP response, which is not supported by Figs. 7b and c, nor by the fact that the Trial amplitude (RMS from 50 to 250 ms poststimulus) was not significantly correlated with the BOLD response (Fig. 7d). Application of the Bayesian framework described in Friston et al. (2006), which allows estimation of the evoked and induced components of the response, may help to disentangle their relative contributions to the BOLD signal. Further empirical and theoretical work is clearly needed to characterise the interdependencies of the EEG and fMRI signals, and analysis on an ST basis can provide a useful contribution to this (Debener et al., 2006). In principle, the mean amplitudes of the EPs and HRs can be affected both by amplitude and latency variability on an ST basis. The data presented here suggest that for both modalities latency variability is much more correlated with the mean amplitude than is the amplitude variability (Figs. 5 and 6). This observation has previously been made for EEG and MEG data, and latency corrected single trials are sometimes used to calculate mean EPs (Quian Quiroga and Garcia, 2003; Jin et al., 1997; Jung et al., 2001; Mayhew et al., 2006; Quian Quiroga et al., in press). Such issues have received less attention in fMRI, where HRs to single events are rarely examined and indeed are often difficult to examine at field strengths below 3 T. However, given the similarity between the effect of latency variability on EP and HR mean amplitudes observed in this study, it may prove fruitful to address this issue further, particularly with simultaneous EEG–fMRI measurements. For example, it would be interesting to determine whether the increased latency observed between trials 15 and 20 in Fig. 3 (top left panel) has a correlate in EEG data. This type of ST effect may have implications for the analysis of EEG and fMRI data, but may also shed light on the neurophysiological underpinnings of the two signals. Analysis of the amplitude distribution of the ST EEG responses, following the method of Sander et al. (2005), did not show any significant difference between WD and ICA in terms of the percentage of trials with a response above background (Table 1). When applied by Sander et al. to early somatosensory evoked field MEG data, up to 97% of trials were found to yield a measurable response. For the high contrast WD data, an average of 78% of trials (range 67–91%, first column of Table 1) demonstrated a P1 response, which was the most robust of the peaks generated by the visual stimulus used in this study (Table 1 and Fig. 2). Whether this discrepancy is a result of methodological differences, is related to differences between EEG and MEG or points to underlying regional variations between the visual and somatosensory systems remains to be seen. To date there has been little systematic examination of ST EEG data in response to different sensory stimuli, but the methodology employed in the current study would allow investigation of such region and modality specific effects. This work was motivated by recent studies which have used ST EEG data recorded during simultaneous EEG–fMRI in the fMRI analysis. This technique is at an early stage, but has considerable potential for integrating EEG and fMRI data and exploiting the complementary strengths of the two techniques. The data presented in this study demonstrate highly significant correlations between the EEG ST component amplitudes extracted using WD and ICA, supporting the use of both in simultaneous EEG–fMRI experiments. Useful information was also gained from fMRI trial by trial variations, suggesting an underlying neuronal origin for this variability. This point is rarely addressed in fMRI analysis, where differences between the average HR and a canonical function are dealt with, for example, by including the temporal derivative in the regression (Friston et al., 1999) without attempting to understand its scope and basis. Regarding the EEG data, in addition to the amplitude and latency information that was examined in this and other studies, EEG data are rich in spectral information, which can also be extracted and used on an ST basis (Lee et al., 2003). ST EEG–fMRI provides a methodology for exploiting the wealth of spatiotemporal information that the two modalities provide, but further work is A.P. Bagshaw, T. Warbrick / NeuroImage 38 (2007) 280–292 291 needed to understand the extent and limitations of the relationship between them. References Arthurs, O.J., Williams, E.J., Carpenter, T.A., Pickard, J.D., Boniface, S.J., 2000. Linear coupling between functional magnetic resonance imaging and evoked potential amplitude in human somatosensory cortex. Neuroscience 101 (4), 803–806. Beckmann, C.F., Smith, S.M., 2004. Probabilistic independent component analysis for functional magnetic resonance imaging. IEEE Trans. Med. Imag. 23 (2), 137–151. Bledowski, C., Kadosh, K.C., Wibral, M., Rahm, B., Bittner, R.A., Hoechstetter, K., Scherg, M., Maurer, K., Goebel, R., Linden, D.E.J., 2006. Mental chronometry of working memory retrieval: a combined functional magnetic resonance imaging and event-related potentials approach. J. Neurosci. 26 (3), 821–829. Bonmasser, G., Schwartz, D.P., Liu, A.K., Kwong, K.K., Dale, A.M., Belliveau, J.W., 2001. Spatiotemporal brain imaging of visual-evoked activity using interleaved EEG and fMRI recordings. NeuroImage 13, 1035–1043. 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 single-trial EEG dynamics including independent component analysis. J. Neurosci. Methods 134, 9–21. di Russo, F., Martinez, A., Sereno, M.I., Pitzalis, S., Hillyard, S.A., 2001. Cortical sources of early components of the visual evoked potential. Hum. Brain Mapp. 15, 95–111. 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.LA., 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. 102 (49), 17798–17803. Friston, K.J., Zarahn, E., Josephs, O., Henson, R.N.A., Dale, A.M., 1999. Stochastic designs in event-related fMRI. NeuroImage 10, 607–619. Friston, K., Henson, R., Phillips, C., Mattout, J., 2006. Bayesian estimation of evoked and induced responses. Hum. Brain Mapp. 27, 722–735. Fox, M.D., Snyder, A.Z., Zacks, J.M., Raichle, M.E., 2006. Coherent spontaneous activity accounts for trial-to-trial variability in human evoked brain responses. Nat. Neurosci. 9 (1), 23–25. Glover, G.H., 1999. Deconvolution of impulse response in event-related BOLD fMRI. NeuroImage 9, 416–429. Henning, S., Merboldt, K.-D., Frahm, J., 2005. Simultaneous recordings of visual evoked potentials and BOLD MRI activations in response to visual motion processing. NMR Biomed. 18, 543–552. Iyer, D., Zouridakis, G., 2007. Single-trial evoked potential estimation: comparison between independent component analysis and wavelet denoising. Clin. Neurophysiol. 118, 495–504. Janz, C., Heinrich, S.P., Kornmayer, H., Bach, M., Hennig, J., 2001. Coupling of neural activity and BOLD fMRI response: new insights by combination of fMRI and VEP experiments in transition from single events to continuous stimulation. Magn. Reson. Med. 46, 482–486. Jin, Y., Potkin, S.G., Patterson, J.V., Sanman, C.A., Hetrick, W.P., Bunney, W.E., 1997. Effects of P50 temporal variability on sensory gating in schizophrenia. Psychiat. Res. 70, 71–81. Jung, T.-P., Makeig, S., Westerfield, M., Townsend, J., Courchesne, E., Sejnowski, T.J., 2001. Analysis and visualisation of single-trial eventrelated potentials. Hum. Brain Mapp. 14, 166–185. Kashikura, K., Kershaw, J., Yamamoto, S., Zhang, X., Matsuura, T., Kanno, I., 2001. Temporal characteristics of event-related BOLD response and visual-evoked potentials from checkerboard stimulation of human V1: a comparison between different control features. Magn. Reson. Med. 45, 212–216. Kilner, J.M., Mattout, J., Henson, R., Friston, K.J., 2005. Hemodynamic correlates of EEG: a heuristic. NeuroImage 28, 280–286. Koch, S.P., Steinbrink, J., Vollringer, A., Obrig, H., 2006. Synchronisation between background activity and visually evoked potential is not mirrored by focal hyperoxygenation: implications for the interpretation of vascular brain imaging. J. Neurosci. 26 (18), 1948–4940. Laskaris, N.A., Liu, L.C., Ioannides, A.A., 2003. Single-trial variability in early visual neuromagnetic responses: an explorative study based on the regional activation contributing to the N70m peak. NeuroImage 20, 765–783. Lee, P.-L., Wu, Y.-T., Chen, L.-F., Chen, Y.-S., Cheng, C.-M., Yeh, T.-C., Ho, L.-T., Chang, M.-S., Hsieh, J.-C., 2003. ICA-based spatiotemporal approach for single-trial analysis of postmovement MEG beta synchronisation. NeuroImage 20, 2010–2030. Liu, L., Ioannides, A.A., 1996. A correlation study of averaged and single trial MEG signals: the average describes multiple histories each in a different set of single trials. Brain Topogr. 8 (4), 385–396. Logothetis, N.K., Pauls, J., Augath, M., Trinath, T., Oeltermann, A., 2001. Neurophysiological investigation of the basis of the fMRI signal. Nature 412, 150–157. 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–694. Makeig, S., Delorme, A., Westerfield, M., Jung, T.-P., Townsend, J., Courchesne, E., Sejnowski, T.J., 2004. Electroencephalographic brain dynamics following manually responded visual targets. PLOS Biol. 2 (6), 0747–0762. Mayhew, S.D., Iannetti, G.D., Woolrich, M.W., Wise, R.G., 2006. Automated single-trial measurement of amplitude and latency of laserevoked potentials (LEPs) using multiple linear regression. Clin. Neurophysiol. 117, 1331–1344. Mulert, C., Jäger, L., Propp, S., Karch, S., Störmann, S., Pogarell, O., Möller, H.-J., Juckel, G., Hegerl, U., 2005. Sound level dependence of the primary auditory cortex: simultaneous measurement with 61-channel EEG and fMRI. NeuroImage 28, 49–58. Odom, J.V., Bach, M., Barber, C., Brigell, M., Marmor, M.F., Tormene, A.P., Holder, G.E., Vaegen, 2004. Visual evoked potentials standard (2004). Documenta Opthalmologica 108, 115–123. Ohla, K., Busch, N.A., Dahlem, M.A., Herrmann, C.S., 2005. Circles are different: the perception of glass patterns modulates early event-related potentials. Vis. Res. 45, 2668–2676. Oostenveld, R., Praamstra, P., 2001. The five percent electrode system for high resolution EEG and ERP measurements. Clin. Neurophysiol. 112, 713–719. Ploner, M., Gross, J., Timmermann, L., Pollok, B., Schnitzler, A., 2006. Oscillatory activity reflects excitability of the human somatosensory system. NeuroImage 32, 1231–1236. Pruessman, K.P., Weiger, M., Scheidegger, M.B., Boesiger, P., 1999. SENSE: sensitivity encoding for vast MRI. Magn. Reson. Med. 42, 952–962. Quian Quiroga, R., Garcia, H., 2003. Single-trial event-related potentials with wavelet denoising. Clin. Neurophysiol. 114, 376–390. Quian Quiroga, R., Atienza, M., Cantero, J.L., Jongsma, M.L.A., in press. What can we learn from single-trial event-related potentials? Chaos Complexity Lett. Robinson, P.A., Dysdale, P.M., van der Merwe, H., Kyriakou, E., Rigozzi, M.K., Germanoska, B., Rennie, C.J., 2006. BOLD responses to stimuli: dependence on frequency, stimulus, form, amplitude and repetition rate. NeuroImage 31, 585–599. 292 A.P. Bagshaw, T. Warbrick / NeuroImage 38 (2007) 280–292 Thomason, M.E., Foland, L.C., Glover, G.H., 2007. Calibration of BOLD fMRI using breath holding reduces group variance during a cognitive task. Hum. Brain Mapp. 28, 59–68. Worsley, K.J., Evans, A.C., Marrett, S., Neelin, P., 1992. A threedimensional statistical analysis for CBF activation studies in human brain. J. Cereb. Blood Flow Metab. 12, 900–918. Sander, T.H., Burghoff, M., Curio, G., Trahms, L., 2005. Single evoked somatosensory MEG responses extracted by time delayed decorrelation. IEEE Trans. Signal Process. 53 (9), 3384–3392. Singh, M., Sungheon, S., Kim, T.-S., 2003. Correlation between BOLDfMRI and EEG signal changes in response to visual stimulus frequency in humans. Magn. Reson. Med. 49, 108–114.