Voxel-wise information theoretic EEG-fMRI feature integration

NeuroImage 55 (2011) 1270–1286 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 ev i e r. c o m / l o c a t e / y n i m g Voxel-wise information theoretic EEG-fMRI feature integration Dirk Ostwald a,b,c,⁎, Camillo Porcaro a,b,d, Andrew P. Bagshaw a,b a School of Psychology, University of Birmingham, UK Birmingham University Imaging Centre, University of Birmingham, UK c Department of Neurology and Bernstein Centre for Computational Neuroscience, Charite´, Berlin, Germany d Institute of Neuroscience, Newcastle University, Medical School, Newcastle upon Tyne, UK b a r t i c l e i n f o a b s t r a c t We have recently proposed the evaluation of a set of information theoretic quantities (ITQs) for the integration of simultaneously acquired EEG-fMRI data (Ostwald, D., Porcaro, C., Bagshaw, A.P., 2010. An information theoretic approach to EEG-fMRI integration of visually evoked responses. Neuroimage. 49, 498– 516). In our previous experimental evaluation of the information theoretic framework, we defined the data subsets from which to calculate the ITQs using a priori constraints. In the case of EEG, this meant that data were extracted from a single electrode, while for fMRI the analysed data came from voxels contained within a sphere surrounding the most responsive voxel of visual cortex. While this approach was a natural starting point for the evaluation of the framework in the application to combined EEG–fMRI data sets, a more principled approach to data selection is desirable. Here, we propose to combine standard fMRI data preprocessing and low-resolution electromagnetic tomography (LORETA) for the evaluation of ITQs across the entire three-dimensional brain space. We apply the proposed method to a simultaneous EEG–fMRI data set acquired during checkerboard stimulation and assess the topographical informativeness of EEG (time and frequency domain) and fMRI features with respect to the stimulus and each other. The resulting information theoretic effect size maps are supplemented with a statistical evaluation based on Gaussian null model simulations using a false-discovery rate procedure. Given the contamination of EEG recordings by artefacts induced by the MR scanning environment we further assessed the influence of different advanced EEG preprocessing methods (independent component analysis and functional source separation) on the information topography. The results of this analysis provide evidence for the topographically focussed informativeness of both EEG and fMRI features with respect to the stimulus, but for the current feature selection do not detect EEG–fMRI activity dependence. More advanced EEG data pre-processing rendered the feature distributions more stimulus-informative, but did not alter the EEG–fMRI activity and conditional dependencies. © 2010 Elsevier Inc. All rights reserved. Article history: Received 5 July 2010 Revised 7 December 2010 Accepted 8 December 2010 Available online 15 December 2010 Introduction Despite the considerable progress that has been made in acquiring and combining simultaneous electroencephalography (EEG) and functional magnetic resonance imaging (fMRI) data, methods of integrating the data in order to fully exploit their strengths and thereby provide a non-invasive neuroimaging tool with high spatiotemporal resolution are needed. We have recently proposed the evaluation of a set of information theoretic quantities (ITQs) for EEG– fMRI integration (Ostwald et al., 2010). This information theoretic approach provides a unifying framework to address questions of neural coding and neurovascular coupling (Schneidman et al., 2003; Panzeri et al., 2008; Quian-Quiroga and Panzeri, 2009). Information theoretic approaches to functional neuroimaging have a number of advantages. For example, they are inherently single-trial oriented, and ⁎ Corresponding author. School of Psychology, University of Birmingham, UK. E-mail address: dirk.ostwald@bccn-berlin.de (D. Ostwald). 1053-8119/$ – see front matter © 2010 Elsevier Inc. All rights reserved. doi:10.1016/j.neuroimage.2010.12.029 hence make use of all information available in the signal (Debener et al., 2006; Debener et al., 2007). They are model-free (although practically some assumptions about the underlying signals are used); they allow non-linear relationships to be investigated; and do not per se contain assumptions about the Gaussianity of the signals (de Araujo et al., 2003; Nevado et al., 2004; Fuhrmann et al., 2007, 2008; Kriegeskorte et al., 2007; Rolls et al., 2009). In our initial evaluation of the information theoretic framework for simultaneous EEG–fMRI recordings (Ostwald et al., 2010), the data subsets from which to calculate the ITQs were defined using relatively arbitrary a priori constraints. In the case of EEG, this meant that data were extracted from a single electrode showing the maximal visual evoked potential (VEP), while for fMRI the analysed data came from voxels contained within a sphere surrounding the most responsive voxel of visual cortex. While this approach was a natural starting point for the evaluation of the framework in the application to combined EEG–fMRI data sets, a more principled approach to spatial data selection is desirable. In this study, we propose to combine standard fMRI data pre-processing and low-resolution electromagnetic D. Ostwald et al. / NeuroImage 55 (2011) 1270–1286 1271 tomography (LORETA, (Frei et al., 2001; Pascual-Marqui, 2002a,b)) for the evaluation of ITQs across the entire three-dimensional brain space. A full three-dimensional evaluation of the set of ITQs requires single-trial estimates of EEG and fMRI signal features at each location (voxel) of the brain (Fig. 1). For fMRI, these estimates can be obtained using standard data pre-processing including normalization to a standard anatomical space, as the data are inherently threedimensional. For the EEG, however, a data transform from the electrode recordings at the scalp surface to normalized anatomical threedimensional brain space is necessary. Several distributed inverse solutions which can estimate electrical activity features at a predefined three-dimensional grid of anatomical locations exist (see (Michel et al., 2004) for an introduction). LORETA is one such inverse solution, which has been relatively widely used with EEG data in the time domain, especially in its modification as sLORETA (PascualMarqui, 2002b), but also in the frequency domain (Frei et al., 2001; Pascual-Marqui, 2002a). Using this inverse solution, it becomes possible to estimate single-trial electrophysiological features for each location of the three-dimensional brain space. Subsequently, given that both data types are projected into the same standardized anatomical spaces (e.g. MNI space), the estimation of information theoretic measures can proceed based on the data from corresponding locations. The aim of this study is three-fold: First, to demonstrate the proposed approach to three-dimensional EEG–fMRI informationbased feature integration using data collected in a visual stimulation (checkerboard) experiment as reported in (Ostwald et al., 2010). Second, to extend the range of EEG features under consideration in the information theoretic framework to the frequency domain, reflecting a large body of work that exists on the co-variation of EEG frequency domain features and the BOLD signal (de Munck et al., 2007; Feige et al., 2005; Goldman et al., 2002; Goncalves et al., 2006; Laufs et al., 2003, 2006). Third, to undertake a detailed assessment of three EEG pre-processing methods on the estimation of combined EEG–fMRI information, reflecting the question of EEG data quality in simultaneous EEG–fMRI recordings and their influence on inferences about neurophysiological phenomena (Bagshaw and Warbrick, 2007; Debener et al., 2006; Mantini et al., 2007; Porcaro et al., 2010). The single-trial data features under consideration in this study are a) for the EEG time-domain, the amplitude of the visually induced standardized current density power in a time-window of 80–120 ms post-stimulus, corresponding to the P100 amplitude in electrode space, b) for the EEG frequency-domain, the relative change of alpha frequency amplitude between a pre-stimulus and post-stimulus time window of ±1 s, and c) for the fMRI modality, the maximum of the haemodynamic response function in a time-window of 0–16.5 s poststimulus. These features were chosen for the following reasons: first, with respect to the EEG modality, the interplay of visually evoked activity, changes in spontaneous and induced alpha activity, and haemodynamic responses remains largely elusive (Becker et al., 2008; Nierhaus et al., 2009). Second, the amplitude of single-trial haemodynamic responses has previously been shown to be sensitive to and informative about the stimulus contrast (Bagshaw and Warbrick, 2007; Logothetis et al., 2001; Ostwald et al., 2010), and hence is a reasonable choice for the spatial integration attempted here. Clearly many more EEG and fMRI features might potentially be informative about both the stimulus and the complementary modality, but the selected features represent the primary variables which are likely to be most relevant. Finally, for the integration of simultaneously acquired EEG–fMRI data it is vital that the data are of good quality, i.e. that data artefacts caused by the physical interactions of both recording set-ups are attenuated as far as possible. This is especially crucial for the EEG, which is subject to various biological (eye-movements, muscle movements, ballistocardiogram artefact) and non-biological (gradient) artefacts during simultaneous acquisition with fMRI data. Independent component analysis (ICA) is a blind signal processing technique increasingly used for single-trial EEG–fMRI studies and has been successfully incorporated into the information theoretic analysis of EEG–fMRI data previously (Ostwald et al., 2010). Functional source separation (FSS), a semi-blind extension of ICA, has recently been demonstrated to reliably improve single-trial EEG data recorded during concurrent EEG–fMRI (Porcaro et al., 2010). In the current study, the influence of these two pre-processing techniques on the voxel-wise estimation of information theoretic quantities is assessed. All custom written Matlab (The Mathworks, Natick, MA) code used in this report is available from http://www.buic.bham.ac.uk/ downloads/EEG_FMRI_ITQ/EEG_FMRI_ITQ_Analysis.zip and the experimental data are available from openfmri.org. Fig. 1. Conceptual framework of voxelwise information theoretic EEG–fMRI feature integration. A necessary prerequisite for the proposed scheme is a voxel single trial time course for each modality, where voxels from both modalities can be brought into spatial correspondence. In the case of the fMRI modality, the data matches this format upon reconstruction by the MRI scanner console, while for the EEG modality, the sensor readings have to be projected into three dimensional space based on a solution of the EEG inverse problem. Individually for both modalities, time-domain features are extracted from the time course of each trial and voxel. Upon feature extraction, the features from both modalities at corresponding locations enter the information theoretic analysis. The bracketed expressions represent the underlying numerical arrays, where X , Y and Z indicate the number of voxels in x-,y-, and z-dimension of the data, Samples indicates the number of data samples of a single trial of the experiment, and Trials represents the number of individual trials across conditions. ITQ represents the number of information theoretic quantities evaluated. 1272 D. Ostwald et al. / NeuroImage 55 (2011) 1270–1286 Methods Experimental paradigm and data acquisition Both experimental paradigm and data acquisition have been described previously (Ostwald et al., 2010). Briefly, 14 observers were presented with left-hemifield reversing checkerboards of two different contrast levels, high (CMichelson = 1) and low (CMichelson = 0.25). Individual trials of the experiment consisted of a single presentation of the stimulus for 1 s with phase reversal after 500 ms followed by a fixation period of 16.5–21 s. A total of 85 trials per contrast condition were acquired for each observer. To maintain the observer's attention and ensure fixation, observers carried out a simple target detection task. The experiment was conducted at the Birmingham University Imaging Centre using a 3T Philips Achieva Scanner. T2*-weighted functional data were collected with an eight channel phased array SENSE head coil. EPI data were acquired from 20 slices (2.5 × 2.5 × 3 mm resolution, TR 1500 ms, TE 35 ms, SENSE factor 2, flip angle 80°) covering the entire visual cortex. Simultaneously, EEG data were recorded using a 64 channel MR compatible EEG system (BrainAmp MR Plus, Brain Products, Munich, Germany) with a sampling rate of 5000 Hz with an impedance at all recording electrodes of less than 20 kΩ. The EEG data acquisition setup clock was synchronised with the MRI scanner clock using Brain Product's SyncBox. Two of the 14 data sets were excluded due to a technical malfunction of the recording set-up resulting in data sets from n = 12 observers in the analyses reported below. Data pre-processing The EEG data were preprocessed as described previously (Ostwald et al., 2010), which for the basic pre-processing involved gradient artefact removal, pulse artefact correction (Optimal Basis Set method (Allen et al., 1998, 2000; Niazy et al., 2005) as implemented as a plugin to EEGLAB (Delorme and Makeig, 2004)), low-pass filtering at 25 Hz, and down-sampling to 500 Hz. In the following, EEG data preprocessed in this manner will be referred to as ‘Basic EEG’ data. Further, the basic EEG data underwent additional artefact correction using Independent Component Analysis (ICA). The EEG data after this pre-processing step will be referred to as ‘ICA EEG’ data. Finally, the EEG data underwent another pre-processing procedure, namely functional source separation (Porcaro et al., 2010), the resulting data from which will be referred to as ‘FSS EEG’. For the FSS procedure additional information is included to bias the independent component decomposition algorithm towards solutions that satisfy physiological assumptions. Specifically, the aim of FSS is to enhance the separation of relevant signals by exploiting some a priori knowledge without renouncing the advantages of using only information contained in the original signal (Barbati et al., 2006; Porcaro et al., 2008, 2009; Tecchio et al., 2007). The computational underpinnings of FSS are discussed in (Porcaro et al., 2010). The criterion maximized for the independent component solution in the present study was the VEP amplitude in an 80–120 ms post-stimulus time-window. The single-trial VEPs of the basic and retro-projected ICA and FSS EEG data were then subjected to the LORETA algorithms discussed below. Single-trial VEPs were obtained with respect to a 100 ms pre-stimulus baseline for the time-domain features and the entire pre-stimulus period for the frequency domain features, respectively. SPM5 (Friston, 2007) was used for fMRI data pre-processing, namely anatomical realignment, slice scan time correction (reference slice 10), anatomical normalization, and re-interpolation to 7 × 7 × 7 mm voxels to allow MNI space alignment with the EEG inverse solutions (see below). Using the SPM5 general linear model approach, additional high-pass filtering (cut-off 1/128 Hz) and data whitening (AR(1)-model) was implemented by pre-multiplication with the corresponding filtering matrices. As a control, in addition to the in- formation theoretic fMRI feature analyses described below, we also performed a conventional GLM analysis of the data. Two experimental regressors, corresponding to the high and low contrast checkerboard conditions were used. Voxel time-courses were modelled in an eventrelated fashion using regressors obtained by convolving each stimulus onset unit impulse with a canonical haemodynamic response function and its first temporal derivative. Additional nuisance covariates included the realignment parameters to account for residual motion artefacts and session specific means. A mixed-effects analysis was implemented using a summary statistics approach to allow inferences at the population level (Frison and Pocock, 1992; Mumford and Nichols, 2009): upon estimation of the model parameters for each subject, a subject-specific contrast image for each effect of interest was computed. Contrast vectors for the following effects of interest were used: Both stimuli N fixation, high contrast stimulus N low contrast stimulus. The contrast images were then subjected to a onesample t-test on the second level. Results are reported at a statistical threshold of p b 0.05, FDR corrected for multiple comparisons. EEG forward and inverse solutions For the current study a solution to the EEG forward problem, i.e. a lead field matrix, was obtained using Curry 6.0.1 (Compumedics Neuroscan, Charlotte, NC). Specifically, the lead field matrix was based on a standard boundary element model (BEM) created from an MNI space template with a conductivity ratio brain/skull of 1/78 (conductivity parameters 0.33 Sm− 1 for scalp, 0.0042 Sm− 1 for skull, and 0.33 Sm− 1 for brain tissue). Due to the high computational demand of the inverse algorithms used, the resolution of the lead field matrix was set to 7 mm3 isotropic in PAN (L,P,S) space, which was the highest possible spatial resolution for computing the inverse solutions under Matlab 7.0.5 (R2007b) on a 64-bit Linux machine. Since LORETA is inherently a smooth, low resolution inverse solution, it is unlikely that considerable benefits would be gained by using a higher resolution matrix. Two inverse solution algorithms were chosen to obtain estimates of electrical brain activity based on the scalp EEG recordings: first, in the time domain, the standardized current density power for each brain location was estimated using sLORETA (Pascual-Marqui, 2002b). Assuming an instantaneous, distributed, discrete and linear solution to the EEG inverse problem, sLORETA can provide activity estimates for each dipole located at all pre-specified MNI voxel locations over time. The algebraic formulation of sLORETA implemented for the current study is reviewed in Appendix A. Second, in the frequency domain, voxel-wise frequency power estimates for peri-stimulus time-windows were obtained using a LORETA-variant specifically tailored to frequency analyses (Frei et al., 2001; Pascual-Marqui, 2002a). This inverse solution is based on a forward model of the Fourier-transformed impressed source current densities and capitalizes on the evaluation of their cross-spectral matrix. The algebraic formulation of this LORETA-variant implemented for the current study is reviewed in Appendix B. In order to validate the implementation of the LORETA algorithms, a number of analyses were performed on simulated data, as reported in supplementary material S1. For the computation of ITQs based on both imaging modalities at a specific voxel location, it is essential that both data types relate to the same anatomical location. The MNI space as defined in SPM5 was used as reference space (Friston, 2007). Specifically, the fMRI data were reinterpolated to the resolution of the EEG lead field matrix obtained from Curry 6.0.1 (Compumedics Neuroscan, Charlotte, NC), i.e. 7 × 7 × 7 mm voxels ranging in the x-dimension between MNI coordinates −63 and 70 mm, in the y-dimension between −105 and 70 mm, and in the z-dimension between − 42 and 35 mm. A detailed description of the MNI space alignment procedure of SPM5's and Curry's MNI spaces can be found in the supplementary material S2. D. Ostwald et al. / NeuroImage 55 (2011) 1270–1286 Table 1 Feature definitions. For the information theoretic analyses, each feature was extracted from the data corresponding to each voxel of the respective modality. Acronym sCDPP100 Δα HRFmax Definition max {post-stimulus sCDP (80–120 ms)} αprestimulus −αpoststimulus αprestimulus max {post-stimulus percent signal change (0–16.5 s)} 1273 Feature definition In (Ostwald et al., 2010), a large set of time-domain features traditionally associated with the study of VEPs was evaluated, not all of which proved to be informative. The same approach taken here would lead to a combinatorial explosion of the number of computations performed, as well as the number of information theoretic maps to be evaluated. Instead, the discussion focuses on the following features extracted from the peri-stimulus time-courses of each voxel in each modality: a) for the EEG modality in the time-domain, the maximum of the standardized current density time-course in a timewindow of 80–120 ms post-stimulus onset, referred to as sCDPP100, b) for the EEG modality in the frequency domain, the relative difference of 1 s pre-stimulus and 1 s post-stimulus alpha power (average of frequency power of 8–12 Hz) given by αprestimulus −αpoststimulus Δα = αprestimulus ð1Þ required empirical distributions are obtained by marginalization. With respect to the free parameters of the histogram analysis, namely the upper and lower limits of the discretization grid and the total number of bins, the same choices as discussed in (Ostwald et al., 2010) were made: the upper and lower limits of the response grid were set to the maximum and minimum value of each response variable and the number of response bins was chosen to be approximately equal to half the number of trials per stimulus resulting in 36 response bins for the bivariate response variable distributions. The observed probability distributions were then subjected to the respective mutual information equations of Table 2. These so-called plug-in information estimates were subsequently corrected for their limited sampling bias (Panzeri et al., 2007). Specifically, a combination of PT-correction, shuffling correction (IN(S ; R1, R2)) and the subtraction of the remaining biases in the analysis of a linear Gaussian null model with the same parameters as the experiment (number of response variables, number of stimuli, number of trials per stimulus) were employed (Ostwald et al., 2010). Finally, for the current study, all ITQs were evaluated in this manner for each voxel of the reference anatomical space. Statistical evaluation To statistically evaluate the derived information theoretic maps (effect size maps), voxel-wise null hypothesis testing in combination with falsediscovery rate correction for multiple comparison was employed. The expected ITQ distributions cannot be assumed to be Gaussian, because the information theoretic estimates are not symmetrically distributed about zero, since, at least in principle and in the absence of bias, negative information values do not occur. Further, as the computational demands for a permutation based analysis would be very high, ITQ null distributions were derived from a Gaussian null model by sampling 10,000 times and evaluating all ITQs (model M0 from (Ostwald et al., 2010)). This resulted in ITQ distributions reflecting the case of non-informative Gaussian signals. It should be noted that these distributions (illustrated in supplementary Fig. S3) are constant under variation of the variance parameters of the M0 model, as they reflect between-simulation, rather than within-simulation, variance. Based on these observed frequency distributions, the corresponding Kaplan–Meier estimates of the cumulative probability functions for each ITQ were computed, and the critical value for one-tailed hypothesis tests at an alpha error level of 0.05 determined. P-values were then assigned to all MNI space voxels for which data from all observers exists. To correct for multiple comparisons the resulting statistical maps were thresholded using a false-discovery rate (FDR) procedure (Genovese et al., 2002). Controlling the FDR is a well established procedure for false-positive control in the neuroimaging literature (Chumbley and Friston, 2009; Friston, 2007; Genovese et al., 2002), and c) for the fMRI modality, the maximum of the haemodynamic response in a of 0–16.5 s post-stimulus time-window, referred to as HRFmax, where the haemodynamic response was quantified as BOLD percent signal change with respect to the pre-stimulus baseline (average of the three pre-stimulus onset TR values). The feature definitions are summarized in Table 1. It should be noted that given the relatively low-frequency cut-off of the filtered EEG data necessary for MR artefact correction, the sCDPP100 does not contain contributions from frequencies larger than 25 Hz. Information theoretic analysis The information theoretic analysis, i.e. the computation of the mutual information of the joint and marginal distributions of two response signal features R1 and R2 about the stimulus S (IN (S; R1, R2) IN (S; R1) IN (S; R2) and Synergy), and the computation of the response signal features activity dependence IN(R1 ; R2) and conditional dependence IN(R1 ; R2|S)S (see Table 2 for details) from experimental data has been discussed in detail previously (Ostwald et al., 2010). Briefly, upon a histogram-based estimation of the stimulus conditional joint probability distribution of the response features pN(R1, R2|S), the Table 2 Information theoretic quantities (ITQs) considered in this study. ITQ I (S; R1) I (S; R2) I (S; R1, R2) Synergy Definition  ∑ ∑ pðs; r1 Þ log2 s∈S r1 ∈R1 Interpretation pðs; r1 Þ pðsÞpðr1 Þ   Mutual information between stimulus and response signal feature of modality 1 (EEG) ∑ ∑ pðs; r2 Þ log2 s∈S r2 ∈R2 pðs; r2 Þ pðsÞpðr2 Þ   pðs; r1 ; r2 Þ ∑ ∑ ∑ pðs; r1 ; r2 Þ log2 pðsÞpðr1 ; r2 Þ s∈S r1 ∈R1 r2 ∈R2  Mutual information between stimulus and response signal feature of modality 2 (fMRI) Mutual information between stimulus and response signal features of modalities 1 and 2 Normalized difference between the mutual information between stimulus and response variables of modalities 1 and 2 and the sum of the individual mutual information for response signal feature 1 and response signal feature 2 (Information Dependence) Overall mutual information between response signal features of modality 1 (EEG) and modality 2 (fMRI) (Activity Dependence) Average stimulus conditional mutual information between response signal features of modality 1 (EEG) and modality 2 (fMRI) (Conditional Dependence) IðS; R1 ; R2 Þ−ðIðS; R1 Þ + I ðS; R2 ÞÞ IðS; R1 ; R2 Þ  ∑ ∑ pðr1 ; r2 Þ log2 ∑ IðR1 ; R2 j s = si Þ  pðr1 ; r2 Þ pðr1 Þpðr2 Þ I (R1; R2) 〈I (R1; R2|S)〉S r1 ∈R1 r2 ∈R2 1 jSj si ∈S 1274 D. Ostwald et al. / NeuroImage 55 (2011) 1270–1286 and allows control of the upper bound of the expectation of the FDR, which is defined as ratio of the number of false-positives to the total number of statistical tests. Results In the following, the results of the whole-brain information theoretic analysis based on the three features of interest introduced above (i.e. sCDPP100 , Δα and HRFmax) will be discussed. The derived information theoretic maps are the outcome of a series of data conversions, e.g. for the EEG data from electrode space data to MNI space data and from time domain data to frequency domain data, for both modalities from raw to feature data. To assess the influence of these various steps on the information theoretic analyses, the results presented here are supplemented by the results displaying the average feature effect sizes and a comparison between the electrode and brain space EEG data. Further, when appropriate, results are presented for all three EEG data pre-processing techniques (basic, ICA, FSS) to demonstrate the influence of EEG data pre-processing on the final information theoretic analysis outcome. An in-depth quantitative comparison of the effect of ICA and FSS on single-trial EEG data acquired simultaneously with fMRI in a non-information theoretic framework has been presented previously (Porcaro et al., 2010). Average feature effects Evaluating the information contained in brain imaging data features with respect to the stimulus (or another brain imaging data feature) amounts to evaluating the ability of that feature to differentiate between stimulus (or other signal feature) properties for a single experimental trial. It hence provides a statistic for the stimulus' (or other data feature) discrimininability in the presence of noise. Nevertheless, more global signal statistics are often also relevant, e.g. the signal's amplitude on average and its sensitivity to various pre-processing methods. In the following section the average effect sizes for stimulus-induced signal variations for the signal features considered in the information theoretic analysis will be discussed. Since EEG data are more often considered at the electrode than the voxel level, particularly when used for fMRI integration, we will first consider the electrode data before moving on to the voxel space. EEG data in electrode space Fig. 2 depicts the group average EEG time and frequency domain data for electrode PO8 for each EEG pre-processing method. The first column depicts the time-domain signal for a peri-stimulus timewindow of −1 to 1 s. The second column depicts the pre-stimulus induced power spectrum for the time-window −1 to 0 s, while the third column depicts the post-stimulus induced power spectrum for the 0 to 1 s time window. In each panel, the black trace indicates high contrast trials, while the grey trace indicates low contrast trials. For all processing methods, the VEP for the checkerboard onset (0 s) and reversal (0.5 s) are clearly visible. Both ICA and FSS pre-processing result in a slight enhancement of the VEP amplitude contrast for the checkerboard onset P100 component. Across both stimulus conditions ICA and FSS pre-processing result in a slight reduction of P100 amplitude (Table 3). The power spectra in Fig. 2 show the band-pass filtering effects (low frequency cut-off 0.25 Hz, high frequency cut-off 25 Hz) of the basic EEG data pre-processing. Further, the overall frequency power of the basic pre-processed EEG data is strongly reduced by the application of ICA, and even further by FSS (Table 3), probably reflecting a reduction in residual artifacts associated with scanning, and a consequent improvement in data quality, as shown in Fig. 2. EEG data in MNI space While the EEG signal depicted in Fig. 2 represents the data of a single EEG electrode, the features of interest in the information theoretic analysis reported here were extracted from a linear combination of the data from all electrodes, namely their LORETA projections to brain space. For the time domain, the sLORETA conversion results in the non-negative standardized current density power (sCDP). It should be noted that this quantity is qualitatively different from the electrode potential, which varies in both positive and negative directions. The peak of the sCDP around 100 ms poststimulus, which results from the conversion of the P100 peak in the respective electrode potentials, is hence referred to not as P100, but sCDPP100. Fig. 3 depicts the topography of the features, while Fig. 4(A–C) displays the corresponding time and frequency domain representations for a specific voxel, selected as the voxel showing the strongest stimulus-induced effects across EEG pre-processing methods. Specifically, the topography plots depict the average feature amplitude over trials and conditions, normalized by the variability (quantified as the standard deviation) of the time-course at each voxel. All methods result in the highest sCDPP100 feature values being located in right occipital cortex, with a clearer contrast to the remaining voxels apparent for the more advanced pre-processing methods (ICA and FSS). For the Δα feature, the strongest topographical focussed changes in pre- to post-stimulus alpha power are observed for the FSS data, with a clear right occipital focus being apparent (Fig. 3C). While Fig. 2 demonstrated the benefit of more advanced preprocessing methods at the electrode level, Fig. 4 suggests that these effects are even more pronounced when dealing with EEG data in voxel space. Comparing the time domain traces shown in the left hand column of Fig. 4 shows very clearly that effective reduction of noise artefacts from the EEG data is of paramount importance when projecting the data into brain space. With the basic pre-processing the sCDPP100 peak is barely distinguishable above the pre-stimulus noise level, with ICA providing some improvement and FSS additional benefit, in line with previous work (Porcaro et al., 2010). For the frequency domain, comparing Figs. 2 and 4 indicates a much more prominent alpha peak in the power spectrum in the voxel-based data. This benefit is evident for all pre-processing methods. Since the constraint employed in the FSS decomposition was temporal and based on the P100 of the VEP, the FSS data are not optimized for observing changes in the alpha band, which would explain the reduction in alpha power between ICA and FSS data (as discussed above, the ICA pre-processing employed here was quite conservative and involved the removal of artefactual components, rather than the identification of specific components of interest). Despite this, it is clear from Fig. 3 that the FSS source was able to localise changes in alpha power successfully, indicating that it contained the majority of the alpha power. fMRI data in MNI space The fMRI feature amplitude is depicted in Fig. 3D (left column). The fMRI feature was defined as the maximal amplitude in a poststimulus time-window of 0–16.5 s and hence the topography also includes voxels which show high signal amplitudes not affected by the stimulus, as for example the voxels containing the ventricles. This signal variability (quantified as the standard deviation) over a 32 s peri-stimulus time-window is depicted in Fig. 3D (right column) and clearly indicates that right occipital voxels display the strongest stimulus induced signal variation. Extraction of the post-stimulus time-courses for a representative occipital voxel (approximated as displaying the strongest stimulus-induced effects at MNI coordinates [7–100 0]) yields the expected haemodynamic response functions D. Ostwald et al. / NeuroImage 55 (2011) 1270–1286 1275 Fig. 2. Group EEG sensor space data, electrode PO8 for all EEG pre-processing methods (A: basic pre-processing, B: additional ICA pre-processing, C: additional FSS pre-processing). The first column displays the time-domain signals for a peri-stimulus time-window of 1 to 1 s, the second column the frequency power spectrum in the pre-stimulus time-window of − 1 to 0 s, and the third column the frequency spectrum in the post-stimulus time-window of 0 to 1 s. peaking at approximately 4 to 5 s post-stimulus, as depicted in Fig. 4D. Finally, Fig. 3 E displays the results of a mass-univariate GLM analysis of the fMRI data carried out using SPM5. The t-value maps depict the outcome of a random-effects analysis thresholded at p b 0.05 with FDR correction. The left panel of Fig. 3 E depicts the result of the comparison of high vs. low stimulus contrast, while the right panel depicts the result of the comparison of both stimuli vs. fixation Table 3 Quantification of EEG data pre-processing methods on time and frequency domain data in electrode and MNI space. All quantifications were performed on the data displayed in Figs. 2 and 4. For the amplitude features (P100 amplitude in μV, sCDPP100 amplitude in μV/m2) the entries reflect the average signal amplitude in a time-window of 80–100 ms post-stimulus onset across stimuli and subjects ± SEM across subjects For the frequency power features, the entries reflect the average frequency power over the pre- and poststimulus time-windows across all frequencies of interest (1–25 Hz) or the alpha band (8–12 Hz) across stimuli and subjects ± SEM across subjects in power units. Basic Electrode Space P100 Amplitude Power Spectrum Alpha Power MNI Space sCDPP100 Amplitude Power Spectrum Alpha Power 7.7 (± 1.2) 52.1 (± 11.6) 88.2 (± 24.2) ICA 7.5 (± 1.2) 18.5 (± 3.0) 38.1 (± 7.1) FSS 5.6 (± 1.4) 1.4 (± 0.4) 2.0 (± 0.5) 174.2 (± 33.1) 31.2E–03 (± 12.5) 53.2E–03 (± 39.7) 80.9 (± 11.3) 2.0E–03 (± 0.5) 4.9E–03 (± 2.1) 23.6 (± 8.1) 1.5E–05 (± 0.4) 2.0E–05 (± 0.6) baseline. As is evident from the figures, the voxels detected to be modulated by both stimuli (right panel) are found in the right occipital region, comparable to the outcome of the variability assessment of the fMRI signal depicted in Fig. 3D (right column). The comparison of high vs. low contrast stimuli shows comparable results to the topography of fMRI stimulus-information (I(S ; HRFmax)) reported in Fig. 7 below (i.e. a contiguous set of right medial occipital voxels). This result is to be expected, as it can be shown that the likelihood ratio statistic employed within the GLM framework is equivalent to mutual information (Friston, 2007; Ostwald and Bagshaw, under revision). However, it should be noted that the feature extraction approach used in the information theoretic analyses, and the hierarchical general linear model/temporal basis function approach used by SPM5, entail that different aspects of the data enter the statistical evaluation and hence also subtle differences in the results occur (e.g. the maximally informative voxel is located at [14–98 0] (right calcarine gyrus, BA 17), while the most significantly activated voxel for the for the comparison of high vs. low stimulus contrast is located at [14–84 14] (right calcarine gyrus, BA 17)). As is evident from Fig. 3, the strongest stimulus-induced feature effects are topographically slightly displaced for the two imaging modalities. In general, the strongest feature effects for the EEG data are observed slightly more lateral and anterior in comparison to the fMRI feature maxima (Euclidean distance approximately 2 cm, i.e. three voxels at the current resolution). However, this effect is less prominent when comparing the GLM results for the both stimulus vs. 1276 D. Ostwald et al. / NeuroImage 55 (2011) 1270–1286 Fig. 3. Topography of feature amplitudes. Figures A–C depict the respective electrophysiological features, sCDPP100 (left column) and Δα (right column), for the different EEG preprocessing methods (A Basic EEG, B ICA EEG, C FSS EEG). The panels depict the ratio of the respective feature amplitudes (sCDP, frequency power) to their standard deviation, and are hence unit free. Figure D depicts the amplitude of the HRFmax feature (left column) as well as the variability of the fMRI signal for the same voxels (right column) in arbitrary units. All feature amplitudes in panels A–D represent the respective group averages across trials and conditions. Fig. 3E depicts mixed-effects (group level) t-value maps for the comparisons ‘high vs. low stimulus contrast’ (left panel) and ‘both stimuli vs. fixation baseline’ (right panel) as the results of a GLM analysis of the fMRI data thresholded at p b 0.05 corrected for multiple comparisons using the FDR procedure implemented in SPM5. The white inset characters in the first slice indicate the orientation of all slices (A: anterior, P: posterior, L: left, R: right). fixation baseline with the EEG feature topographies. The smooth solution provided by LORETA, however, means that there is considerable overlap between the voxels showing strong feature effect in the two modalities. Information theoretic quantity maps MNI space information theoretic maps derived for all ITQs and features of interest are presented in the following manner (Figs. 5–9): D. Ostwald et al. / NeuroImage 55 (2011) 1270–1286 1277 Fig. 4. Peri-stimulus signal time and frequency domain representations extracted from representative right occipital voxels (MNI coordinates EEG data: [25–80 0]. MNI coordinates fMRI data: [7–100 0], coordinates identified as strongly activated voxels for the respective signal topographies). Panels A–C depict the respective electrophysiological time courses (sCDP) for basic EEG (A), ICA EEG (B) and FFS EEG (C). As in Fig. 2, the first column displays the time-domain signals for a peri-stimulus time-window of − 1 to 1 s, the second column the frequency power spectrum in the pre-stimulus time-window of − 1 to 0 s, and the third column the frequency spectrum in the post-stimulus time-window of 0 to 1 s. Panel D depicts the peri-stimulus hemodynamic response function. the first row of slices of each figure's panel displays the respective voxel-wise information values in bits, i.e. the mutual information effect size. The scaling of the information effect size was determined to emphasize topographical features of the resulting maps. The second row of each panel depicts the statistical evaluation (p-values) of the effect size map based on the Gaussian null model simulations discussed in section 2.6 and Fig. S3. This statistical map is thresholded at p b 0.05 and FDR corrected at q = 0.5 for all information theoretic measures, i.e. the same thresholding procedure is applied to all measures. Finally, in some cases, a more stringent threshold can be chosen which demonstrates the spatial topography more clearly, in which case this was done, and included as a third row of slices. The thresholded statistical maps are overlaid on a high-resolution (2 mm isotropic) group average T1 image (avg152T1.nii as provided by SPM5). For all maps, four consecutive slices, depicting approximately the z-directional centre of the acquired EPI stack, are displayed. All slices are oriented along the posterior–anterior brain axis from left to right and the left to right brain axis from top to bottom. All coordinate labels are MNI coordinates. Table 4 lists the maximally informative voxel coordinates and extracted information values for all feature and method combinations with their anatomical locations. Stimulus-response characteristics: I(S ; R) topography Fig. 5 displays the stimulus related information maps for the sCDPP100 feature and the three types of EEG data pre-processing 1278 D. Ostwald et al. / NeuroImage 55 (2011) 1270–1286 Fig. 5. I (S;P100): Stimulus-related information of the sCDPP100 feature for all pre-processing methods. A) basic EEG pre-processing, B) ICA pre-processing, C) FSS pre-processing. The first row of slices of each panel depicts the information theoretic effect size (mutual information in bits). The second row of slices of each panel depicts the statistical evaluation of the effect size map thresholded at p = 0.05 (FDR corrected at q = 0.5). In the case of C) a more stringent FDR corrected statistical map at q = 10− 4 is included as third row in addition. The white inset characters in the first slice indicate the orientation of all slices (A: anterior, P: posterior, L: left, R: right). (A: Basic, B: ICA, C: FSS). For all cases, a cluster of voxels in the right occipital cortex is informative about the stimulus, and is most pronounced for z-coordinates of 0 to 14 mm. For the basic EEG data, none of these voxels carries information large enough to be statistically significant, while a few voxels survive the statistical correction for the ICA data in the z = 14 plane. The FSS EEG data show the most pronounced cluster of informative voxels and the highest information estimates of all three pre-processing methods. D. Ostwald et al. / NeuroImage 55 (2011) 1270–1286 1279 Fig. 6. I (S; Δα) Stimulus-related information of the Δα feature for all pre-processing methods. A) basic EEG pre-processing, B) ICA pre-processing, C) FSS pre-processing. The first row of slices of each panel depicts the information theoretic effect size (mutual information in bits). The second row of slices of each panel depicts the statistical evaluation of the effect size map thresholded at p = 0.05 (FDR corrected at q = 0.5). The white inset characters in the first slice indicate the orientation of all slices (A: anterior, P: posterior, L: left, R: right). Thresholding this statistical map at q = 0.5 indicated all voxels to be significantly informative. However, using a strongly conservative threshold of q = 10− 4 again implied right occipital voxels to be maximally and most reliably informative. Fig. 6 depicts the stimulus-related information of the Δα frequency power feature for all EEG data pre-processing methods. In comparison to the sCDPP100 feature, the information estimates are lower and less topographically focussed on right occipital cortex. However, as for the time-domain feature, the largest information estimates are observed for the FFS data pre-processing, with a number of voxels displaying statistically significant results when thresholded at q = 0.5. Although in the un-thresholded information maps (Fig. 6C, top row) there is some suggestion of a cluster of informative voxels in the right occipital region, when thresholded for statistical significance the suprathreshold voxels do not appear to be topographically meaningful. It should be noted that the FFS procedure was applied only once and used to optimize the functional signal-to-noise ratio in the time domain. From this single extracted source, both time and frequency domain EEG features were then characterized. A more rigorous approach would be to apply separate time and frequency constraints 1280 D. Ostwald et al. / NeuroImage 55 (2011) 1270–1286 Fig. 7. I (S, HRFmax) Stimulus-related information of the HRFmax feature. The first row of slices of each panel depicts the information theoretic effect size (mutual information in bits). The second row of slices of each panel depicts the statistical evaluation of the effect size map thresholded at p = 0.05 (FDR corrected at q = 0.5). The third row of slices of each panel depicts the statistical evaluation of the effect size map at a more liberal FDR threshold of q = 0.7. The white inset characters in the first slice indicate the orientation of all slices (A: anterior, L: left, R: right). and characterize the respective features on independent sources. From this point of view the results for the alpha band may be suboptimal, and the relationship between sources extracted with time and frequency domain constraints remains to be clarified. Future application of information theoretic concepts to frequency-domain EEG data might shed light on the information gain provided by domain specific FSS data quality optimization. Finally, Fig. 7 displays the stimulus-related information of the HRFmax feature. In comparison to the LORETA derived EEG feature information maps, the informative voxels for the case of fMRI modality are much more focussed and located slightly more ventral in the right occipital cortex (Table 5.3). While for the common threshold of q = 0.5 no voxels show significant informativeness, relaxing the threshold to q = 0.7 yielded a cluster of significantly informative voxel in right occipital cortex (z = 0). Topography of feature combinations across modalities: I(S ; R1, R2) Figs. 8 and 9 depict the stimulus related information estimates derived from the joint distributions of sCDPP100, HRFmax and Δα, HRFmax, respectively (A: Basic EEG, B: ICA EEG, C: FSS EEG). For I(S ; sCDPP100, HRFmax) with the standard threshold of q = 0.5, all maps display large numbers of significantly informative voxels in a topographically unspecific manner (Fig. 8). Using more conservative q-value thresholds as annotated in the Figure, however, implicates mainly right occipital voxels in the representation of stimulus-related information. For I(S ; Δα, HRFmax) very little topographical specificity is observed in any of the maps (Fig. 9). In comparison to their univariate equivalents, the estimated information quantities are larger (Table 3) and more voxels show significant results at lower FDR correction q-values. For example, in the case of the bivariate combination of the basic and ICA preprocessed sCDPP100 feature with the HRFmax feature, a large number of significantly informative voxels can be observed in right occipital cortex at q = 0.08, compared to very few statistically significant informative voxels at q = 0.5 in the case of the univariate EEG features. Comparing the joint distribution information theoretic maps across EEG pre-processing methods indicates that for both timeand frequency-domain bivariate combinations, the information estimates are largest and most topographically focussed for the case of the EEG FSS data (Table 3). Overall, compared to the univariate EEG stimulus-related information feature maps (Figs. 5 and 6), the combination with the fMRI derived feature results in a generally more spatially focussed cluster of most informative voxels, at least in the case of the observed effect sizes (Table 3). In line with our previous evaluation at the electrode level, the evaluation of EEG–fMRI feature Synergy and activity/conditional dependence (I(R1 ; R2)/I(R1 ; R2|S)S) yielded topographical null results which are discussed in the supplementary material (Figs. S4–S6). Discussion In this study, we propose the voxel-wise evaluation of set of ITQs for the integration of EEG and fMRI data following their projection into a standardized anatomical space. The framework was practically applied to a set of three data features, the standardized current density power equivalent of the P100 (sCDPP100), the change in alpha power for a 1 s pre- and post-stimulus time window (Δα) and the maximum of the post-stimulus haemodynamic response function (HRFmax) acquired under low and high contrast checkerboard stimulation. The analyses revealed topographically specific stimulus- Fig. 8. I (S; sCDPP100, HRFmax) Stimulus-related information of the sCDPP100/HRFmax feature joint distribution for all EEG pre-processing methods. A) basic EEG pre-processing, B) ICA pre-processing, C) FSS pre-processing. The first row of slices of each panel depicts the information theoretic effect size (mutual information in bits). The second row of slices of each panel depicts the statistical evaluation of the effect size map thresholded at p = 0.05 (FDR corrected at q = 0.5). Finally, the third row of each panel depicts the statistical map at a more stringent FDR threshold to emphasize topographical features as annotated. The white inset characters in the first slice indicate the orientation of all slices (A: anterior, P: posterior, L: left, R: right). D. Ostwald et al. / NeuroImage 55 (2011) 1270–1286 1281 1282 D. Ostwald et al. / NeuroImage 55 (2011) 1270–1286 D. Ostwald et al. / NeuroImage 55 (2011) 1270–1286 1283 Table 4 Maximal informative voxel coordinates and information values for the different feature and method combinations employed, evaluated over the entire scan volume. The Brodmann Area (BA) and anatomical locations corresponding to the respective MNI coordinates were determined using the SPM Anatomy Toolbox V1.7 (Eickhoff et al., 2005). MNI X I(S; sCDPP100)RAW I(S; sCDPP100)ICA I(S; sCDPP100)FSS I(S; Δα)RAW I(S; Δα)ICA I(S; Δα)FSS I(S; HRFmax) I(S; sCDPP100, HRFmax)RAW I(S; sCDPP100, HRFmax)ICA I(S; sCDPP100, HRFmax)FSS I(S; Δα, HRFmax)RAW I(S; Δα, HRFmax)ICA I(S; Δα, HRFmax)FSS I(sCDPP100; HRFmax)RAW I(sCDPP100; HRFmax)ICA I(sCDPP100; HRFmax)FSS I(Δα; HRFmax)RAW I(Δα; HRFmax)ICA I(Δα; HRFmax)FSS b I(sCDPP100; HRFmax|S) N RAW b I(sCDPP100; HRFmax|S) N ICA b I(sCDPP100; HRFmax|S) N FSS b I(Δα ;HRFmax|S) N RAW b I(Δα ;HRFmax|S) N ICA b I(Δα ;HRFmax|S) N FSS 28 28 0 56 56 − 49 14 14 14 14 14 14 14 28 7 − 28 56 − 21 − 49 42 − 42 49 42 42 42 MNI Y − 77 − 77 − 56 14 0 − 63 − 98 − 98 − 98 − 98 − 98 − 98 − 98 −7 14 35 7 − 77 7 − 42 49 − 14 − 84 − 84 − 84 MNI Z 35 28 35 0 −7 0 0 0 0 0 0 0 0 7 14 −7 14 35 − 21 21 7 − 14 − 14 − 14 − 14 BA – – 7A 44 TE – 17 17 17 17 17 17 17 – – – 44 7A – – – – hOC4v hOC4v hOC4v Anatomical Location Right superior occipital gyrus Right middle occipital gyrus Right precuneus Right inferior frontal gyrus Right superior temporal gyrus Left middle temporal gyrus Right calcarine gyrus Right calcarine gyrus Right calcarine gyrus Right calcarine gyrus Right calcarine gyrus Right calcarine gyrus Right calcarine gyrus Right putamen Right caudate nucleus Left inferior frontal gyrus Right rolandic operculum Left superior occipital gyrus Left temporal pole Right superior temporal gyrus Left middle frontal gyrus Right middle temporal gyrus Right inferior occipital gyrus Right inferior occipital gyrus Right inferior occipital gyrus Information (bits) 0.036 0.034 0.143 0.035 0.028 0.035 0.045 0.116 0.113 0.230 0.132 0.118 0.138 0.043 0.045 0.051 0.036 0.037 0.033 0.144 0.150 0.143 0.161 0.159 0.157 related information for both the features' marginal and joint distributions, for the EEG data predominantly under the use of more advanced pre-processing methods. The inter-modality information analyses remain inconclusive, with no spatial preference of the respective feature combination's activity or conditional dependence. The approach introduced here has the advantage over the previously proposed use of an information theoretic approach to EEG–fMRI integration (Ostwald et al., 2010) that it does not require a priori constraints on the data selection. Previously, data were selected from a single EEG electrode which showed the maximum VEP, and fMRI voxels contained within a sphere centred on the maximum response in a standard GLM analysis. This contains an implicit, strong assumption regarding the spatial correspondence between EEG and fMRI, namely that there is an underlying relationship between the electrical data recorded over the occipital region and the haemodynamic data recorded within visual cortex. While this assumption is at first sight reasonable, it is not clear that it necessarily holds true when examined at the level of single trial data. As well as issues directly related to neurovascular coupling at the macroscopic level, the signal recorded by a scalp electrode contains contributions from a relatively wide area of cortex (Gloor, 1985) not all of which is likely to have a detectable haemodynamic response (Nunez and Silberstein, 2000). By incorporating EEG source localisation into the information theoretic framework, the effect of these issues is reduced, potentially allowing a more principled investigation of the relationship between EEG and fMRI. The results presented in this study demonstrate that such an approach is feasible and that the cortical regions known to be involved in the simple task we employed are highlighted as containing stimulus-relevant information. For the EEG time-domain, the current study analysed the same data as discussed in (Ostwald et al., 2010), however, with important differences in data pre-processing and feature selection (i.e. based on electrode ICA data in (Ostwald et al., 2010) and using FSS and sLORETA in the current study). Perhaps unsurprisingly given the different input data, there are some discrepancies between the studies in terms of magnitude and behaviour of the calculated ITQs. In particular, it was previously found that the fMRI data were more informative than the EEG data regarding the stimulus contrast, whereas the reverse was found in the current study. There are two, probably complementary, explanations for this. On the one hand, the EEG data used in the current study were pre-processed using FSS rather than ICA, which we have shown greatly improves the quality of data acquired in the MRI scanner (Porcaro et al., 2010). In addition, in the previous study we extracted fMRI data from voxels selected on an individual subject basis as having the highest stimulus response in a GLM analysis. In the current analysis, fMRI data were extracted from coordinates defined by the group activation, which is likely to lead to less optimal data regarding stimulus informativeness. The issue here is one of inter-subject spatial heterogeneity in the fMRI response, although the low spatial resolution used in the current study should reduce its effect on the calculation of ITQs. The improvement in EEG data quality using FSS will certainly lead to improved calculation of ITQs, as demonstrated previously for the comparison with raw and ICA-processed data (Porcaro et al., 2010). The results regarding activity dependence (Fig. S5) indicate that, at the level of the false discovery rate threshold used, there were no brain regions in which the P100 amplitude was informative about the HRF amplitude. This is in accordance with our previous work, where the absolute values of activity dependence were low and not significantly different to zero (Ostwald et al., 2010), see also (Whittingstall et al., 2007). This null result highlights the need for a more detailed investigation of the EEG signal space to identify features which carry information about the haemodynamic response. The spatiotemporal framework outlined here allows this issue to be Fig. 9. I (S; Δα, HRFmax) Stimulus-related information of the sCDPP100/HRFmax feature joint distribution for all EEG pre-processing methods. A) basic EEG pre-processing, B) ICA preprocessing, C) FSS pre-processing. The first row of slices of each panel depicts the information theoretic effect size (mutual information in bits). The second row of slices of each panel depicts the statistical evaluation of the effect size map thresholded at p = 0.05 (FDR corrected at q = 0.5). Finally, the third row of each panel depicts the statistical map at a more stringent FDR of q = 0.05. The white inset characters in the first slice indicate the orientation of all slices (A: anterior, P: posterior, L: left, R: right). 1284 D. Ostwald et al. / NeuroImage 55 (2011) 1270–1286 addressed explicitly, to identify when and where the EEG signal is related to the fMRI signal. One feature that is likely to be relevant in this context is EEG activity in the gamma band (~30–100 Hz), which has previously been suggested to be particularly informative about the BOLD response (Muthukumaraswamy et al., 2009; Swettenham et al., 2009). While the current analysis was restricted to frequencies below 25 Hz for technical reasons, it is to be expected that in the future the usable bandwidth of EEG concurrently recorded with fMRI will increase. The framework that has been presented here, exemplified by the analysis of the alpha band activity, is readily generalisable to take account of these technical improvements. We note that the estimation of the conditional dependence (Fig. S6) across the whole-brain volume appears uniformly positively biased, yielding a uniform pattern of statistical significance at the current false-discovery threshold. This result is most likely due to non-optimal information estimation and bias correction. We are currently investigating the application of Gaussian mixture models for the estimation of information from continuous, analogue brain signals and hope to clarify this point in a future communication. Two recent studies have proposed approaches to combined EEG– fMRI data with some similarities to the scheme proposed here (Esposito et al., 2009; Martuzzi et al., 2009). However, while both previous studies propose the projection of EEG electrode data into an anatomical space in alignment with the fMRI data, they differ with respect to how the integration of the two modalities is achieved. Specifically, in (Martuzzi et al., 2009) the authors use a different approach to estimate intracranial source activations, in their case estimated intracranial local field potentials (eLFPs) based on the ELECTRA algorithm (Grave de Peralta et al., 2000; Grave de Peralta et al., 2004). In order to integrate data from both modalities, the authors then compute a similarity index between the statistical activation maps obtained from EEG and fMRI data analysis for (occipital) regions of interest. While this is a sensible approach to compare the outcome of the application of statistical procedures to data from both modalities, the information theoretic framework allows a more direct evaluation of the information provided by the joint observation of the two modalities with respect to the stimulus in the quantities IN(S ; R1, R2) and Synergy. In (Esposito et al., 2009), the authors use another projection of EEG sensor data into anatomical space, namely a cortically constrained distributed inverse model based on the minimum L2-norm approach (Dale et al., 2000). To integrate EEG and fMRI data, the authors then extract single-trial features from the estimated source time-courses and enter these into a mass-univariate GLM analysis. Again, the procedures proposed in the current study are more direct, in that they assume a voxel-byvoxel correspondence of the reconstructed EEG and fMRI signal. However, the general approach of integrating EEG–fMRI data by working in a standard, common space appears to be feasible and have some advantages over more restricted methods. There are a number of ways in which the procedures described here can be improved in future research. First, as discussed in (Ostwald et al., 2010), a bias correction approach tailored to the peculiarities of EEG and fMRI data is desirable. The approach taken in this and the previous study is relatively straightforward to implement, but seems to lead to over-correction of some ITQs. As discussed above, some improvement might be observed by applying entropy estimation procedures for Gaussian mixture models, a topic we are currently investigating. Second, computational constraints limited the analysis to a resolution of 7 mm isotropic, but this is a purely practical issue. Higher resolution analyses will be possible with improved hardware and potentially distributed processing. However, it should be noted that given the relatively low resolution of EEG source localisation, it is debatable whether a considerably higher spatial resolution would be beneficial. Finally, there a considerable number of inverse solutions in addition to LORETA, as well as additional improvements to be gained by the use of individual geometry in the forward model (Fuchs et al., 1998, 2001). Improvements in the spatial accuracy of the EEG source localisation might help to optimize the proposed information theoretic analysis scheme, while inverse solutions more amenable to the localisation of spectral data, such as beamforming (Brookes et al., 2008; Hillebrand and Barnes, 2005) would allow an additional level of inter-modality comparison. It is possible that the low information values observed when calculating the inter-modality comparisons I(R1 ; R2) could have been the result of inaccuracies in the EEG inverse solution leading to a spatial mismatch between EEG and fMRI, an issue that could be improved using more accurate forward models rather than a standard head model. However, as demonstrated in Figs. 3 and 5, there was considerable overlap between the voxels with large effect sizes in fMRI and EEG, primarily as a result of the smoothness of the LORETA solutions, suggesting that this is not the primary explanation for the lack of correspondence between EEG and fMRI. In conclusion, the feasibility of integrating EEG–fMRI data within an information theoretic framework across the entire threedimensional brain space and at the group level was demonstrated. Regions demonstrating significant stimulus-related information were anatomically well localised with both data acquisition modalities. No regions were identified in which the EEG features were activity dependent regarding the amplitude of the HRF, corroborating our previous work at the electrode level (Ostwald et al., 2010). The formalism described here provides a methodology for identifying spatiotemporal or spatiospectral relationships between EEG and fMRI in terms of information content which utilises all of the available EEG data, rather than that selected a priori from a single electrode. This approach may help to uncover a more detailed picture of the regional co-variation between electrical and haemodynamic measures of brain function. Acknowledgments This work was supported by grant number EP/F023057/1 from the UK Engineering and Physical Sciences Research Council (EPSRC). Appendix A. sLORETA Following (Pascual-Marqui, 2002b) we review the time-domain sLORETA inverse solution. Assuming linear superposition of potentials, the electrode potentials at a given time point t  R at the scalp surface are given by V ðt Þ = LJ ðt Þ + c1 ðA:1Þ Denoting the number of electrodes as NE and the number of voxels as NV, Eq. (A.1) describes the electrode potentials V ðt Þ∈RNE ×1 as the matrix product of the lead field matrix L∈RNE ×ð3NV Þ and the current density vector J ðt Þ∈R3NV ×1 plus a constant offset given by c1, where c ∈ R and 1 ∈ RNE is a vector of ones. The experimentally obtained data is represented by V(t), while the data projected into three-dimensional brain space is an estimate for J(t). The lead field matrix L∈RNE ×ð3NV Þ is a time-invariant solution to the EEG forward problem and has the following structure l11 l12 B l21 l22 L=B @… lNE 1 lNE 1 0 1 … l1NV … l2NV C C A … … lN E N V ðA:2Þ D. Ostwald et al. / NeuroImage 55 (2011) 1270–1286 1285 with lij ∈R1×3 for i = 1,..., NE and j = 1,..., NV. The entries lij in the lead field matrix represent the projection from the dipole moments at the jth voxel to the ith electrode. Specifically, lij =  x lij y lij z lij ˆ where J ω ∈ C3NV ×NT and Vω ∈ CNE ×NT . For a single discrete frequency ωi expression (B.2) is equivalent to ˆ J ωi = TVωi 3N ×1 N ×1  ðB:3Þ ðA:3Þ where is the scalp electric potential at the ith electrode, due to a unit strength x-oriented dipole at the jth voxel, and correspondingly for ly ij and lz . ij sLORETA itself is based on the explicit minimization of the function 2 2 lx ij ˆ and and Vωi ∈ C E . As for any discrete Fourier with J ωi ∈ C V transform, although there are NT discrete frequencies, half of them are N ×1 conjugates of the other half, and hence redundant. Let Vωi ∈ C E denote the discrete Fourier transform of a given EEG epoch with NT samples and frequency ωi. From equation (5.15) the cross-spectral matrix of the primary impressed current densities S^ ∈ C3NV ×3NV can Jω i then be derived as follows SJˆ = ωi F ð J ðt ÞÞ = ‖V ðt Þ−LJ ðt Þ−c1‖ + α‖J ðt Þ‖ ðA:4Þ 1 N ˆ ˆÃ ∑ J J φ i = 1 ωi ωi ðB:4Þ where α ≥ 0 ∈ R denotes a regularization parameter and ‖ ‖ denotes an lp-Norm. The explicit solution to this minimization problem is given by the matrix product  h i⊥  ^ T T V ðt Þ J ðt Þ = L H HLL H + αH ðA:5Þ Appendix C. Supplementary data Supplementary data to this article can be found online at doi:10.1016/j.neuroimage.2010.12.029. 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., 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 Mapp. 27, 925–934. Becker, R., Ritter, P., Villringer, A., 2008. Influence of ongoing alpha rhythm on the visual evoked potential. Neuroimage 39, 707–716. Brookes, M.J., Mullinger, K.J., Stevenson, C.M., Morris, P.G., Bowtell, R., 2008. Simultaneous EEG source localisation and artifact rejection during concurrent fMRI by means of spatial filtering. Neuroimage 40, 1090–1104. Chumbley, J.R., Friston, K.J., 2009. False discovery rate revisited: FDR and topological inference using Gaussian random fields. Neuroimage 44, 62–70. Dale, A.M., Liu, A.K., Fischl, B.R., Buckner, R.L., Belliveau, J.W., Lewine, J.D., Halgren, E., 2000. Dynamic statistical parametric mapping: combining fMRI and MEG for highresolution imaging of cortical activity. Neuron 26, 55–67. de Araujo, D.B., Tedeschi, W., Santos, A.C., Elias, J., Neves, U.P.C., Baffa, O., 2003. Shannon entropy applied to the analysis of event-related fMRI time series. Neuroimage 20, 311–317. de Munck, J.C., Goncalves, S.I., Huijboom, L., Kuijer, J.P., Pouwels, P.J., Heethaar, R.M., Lopes da Silva, F.H., 2007. The hemodynamic response of the alpha rhythm: an EEG/ fMRI study. Neuroimage 35, 1142–1151. de Munck, J.C., Goncalves, S.I., Mammoliti, R., Heethaar, R.M., Lopes da Silva, F.H., 2009. Interactions between different EEG frequency bands and their effect on alpha-fMRI correlations. Neuroimage 47, 69–76. 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. Debener, S., Ullsperger, M., Siegel, M., Engel, A., 2007. Towards single-trial analysis in cognitive brain research. Trends Cogn. Sci. 11, 502–503. Delorme, A., Makeig, S., 2004. EEGLAB: an open source toolbox for analysis of singletrial EEG dynamics including independent component analysis. J. Neurosci. Meth. 134, 9–21. Eickhoff, S.B., Stephan, K.E., Mohlberg, H., Grefkes, C., Fink, G.R., Amunts, K., Zilles, K., 2005. A new SPM toolbox for combining probabilistic cytoarchitectonic maps and functional imaging data. Neuroimage 25, 1325–1335. Esposito, F., Mulert, C., Goebel, R., 2009. Combined distributed source and single-trial EEG– fMRI modeling: application to effortful decision making processes. Neuroimage 47, 112–121. Feige, B., Scheffler, K., Esposito, F., Di, S.F., Hennig, J., Seifritz, E., 2005. Cortical and subcortical correlates of electroencephalographic alpha rhythm modulation. J. Neurophysiol. 93, 2864–2872. Frei, E., Gamma, A., Pascual-Marqui, R., Lehmann, D., Hell, D., Vollenweider, F.X., 2001. Localization of MDMA-induced brain activity in healthy volunteers using low resolution brain electromagnetic tomography (LORETA). Hum. Brain Mapp. 14, 152–165. Frison, L., Pocock, S.J., 1992. Repeated measures in clinical trials: analysis using mean summary statistics and its implications for design. Stat. Med. 11, 1685–1704. Friston, K.J., 2007. Statistical parametric mapping. Academic Press, London. Fuchs, M., Wagner, M., Wischmann, H.A., Kohler, T., Theissen, A., Drenckhahn, R., Buchner, H., 1998. Improving source reconstructions by combining bioelectric and biomagnetic data. Electroencephalogr. Clin. Neurophysiol. 107, 93–111. where ⊥ denotes the Moore-Penrose pseudo-inverse, and H = T I− 111 ∈RNE ×NE denotes a centring matrix. 1T Further, the variance of the estimated current density is given by SJˆ = L T  ⊥ T LL + αH L ðA:6Þ sLORETA then corresponds to the following estimate of standardized current density power sCDPj ðt Þ = J ðt Þl ^ ^T h i −1 ^ S^ J ðt Þl J jj ðA:7Þ where J ðt Þl ∈ R3×1 is the current density estimate at the jth voxel and h i SJˆ ∈ R3×3 is the jth diagonal block of matrix SJˆ. jj In essence, given the data of electrical potentials at all electrodes V(t) at time t ∈ R, the standardized current density estimate at all voxels is given by matrix products of V(t) and the lead-field matrix inverse. Appendix B. Frequency-domain LORETA Following (Pascual-Marqui, 2002a), we review the frequencydomain LORETA inverse solution. Let F∈CNT ×NT denote the discrete Fourier transform operator for a time-domain signal of NT time points ˆ and let J ∈R3NV ×NT denote the matrix of all estimated primary impressed current densities at all voxels 1, …, NV and time points 1, …, NT. Further, let V∈RNE ×NT denote the matrix of all sensor potentials at all time points of an EEG epoch of interest. The Fourier transform of the current densities is then given by ˆ J F = TVF ðB:1Þ   Ã⊥  ∈R3NV ×NT (this formulation of the where T = LT H HLLT H + αH pseudo-inverse assumes that the discrete Laplacian matrix B and the lead field normalizing matrix W in (Pascual-Marqui, 2002a) are set to ˆ the respective identity matrices). Equivalently, let J ω and Vω denote the frequency domain representation of the impressed current ˆ densities J and sensor potentials V, respectively, i.e. the columns of ˆ J ω and Vω correspond to all discrete frequencies. Then (5.13) can be equivalently written as ˆ J ω = TVω ðB:2Þ 1286 D. Ostwald et al. / NeuroImage 55 (2011) 1270–1286 Nevado, A., Young, M., Panzeri, S., 2004. Functional imaging and neural information coding. Neuroimage 21, 1083–1095. 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. Nierhaus, T., Schon, T., Becker, R., Ritter, P., Villringer, A., 2009. Background and evoked activity and their interaction in the human brain. Magn. Reson. Imaging 27, 1140–1150. Nunez, P.L., Silberstein, R.B., 2000. On the relationship of synaptic activity to macroscopic measurements: does co-registration of EEG with fMRI make sense? Brain Topogr. 13, 79–96. Ostwald, D., Bagshaw, A.P. Information Theoretic Approaches to Functional Neuroimaging, under revision. Ostwald, D., Porcaro, C., Bagshaw, A.P., 2010. An information theoretic approach to EEG–fMRI integration of visually evoked responses. Neuroimage 49, 498–516. Panzeri, S., Senatore, R., Montemurro, M.A., Petersen, R.S., 2007. Correcting for the sampling bias problem in spike train information measures. J. Neurophysiol. 98, 1064–1072. Panzeri, S., Magri, C., Logothetis, N., 2008. On the use of information theory for the analysis of the relationship between neural and imaging signals. Magn. Reson. Imaging 26, 1015–1025. Pascual-Marqui, R.D., 2002a. Notes on some physics, mathematics, and statistics for time and frequency domain LORETA. Pascual-Marqui, R.D., 2002b. Standardized low-resolution brain electromagnetic tomography (sLORETA): technical details. Meth. Find. Exp. Clin. Pharmacol. 24, 5–12 Suppl D. Porcaro, C., Barbati, G., Zappasodi, F., Rossini, P.M., Tecchio, F., 2008. Hand sensorymotor cortical network assessed by functional source separation. Hum. Brain Mapp. 29, 70–81. Porcaro, C., Coppola, G., Di, L.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., Ostwald, D., Bagshaw, A.P., 2010. Functional source separation improves the quality of single trial visual evoked potentials recorded during concurrent EEG– fMRI. Neuroimage 50, 112–123. Quian-Quiroga, R., Panzeri, S., 2009. Extracting information from neuronal populations: information theory and decoding approaches. Nat. Rev. Neurosci. 10, 173–185. Rolls, E.T., Grabenhorst, F., Franco, L., 2009. Prediction of subjective affective state from brain activations. J. Neurophysiol. 101, 1294–1308. Schneidman, E., Bialek, W., Berry, M., 2003. Synergy, redundancy, and independence in population codes. J. Neurosci. 23, 11539–11553. Swettenham, J.B., Muthukumaraswamy, S.D., Singh, K.D., 2009. Spectral properties of induced and evoked gamma oscillations in human early visual cortex to moving and stationary stimuli. J. Neurophysiol. 102, 1241–1253. Tecchio, F., Porcaro, C., Barbati, G., Zappasodi, F., 2007. Functional source separation and hand cortical representation for a brain-computer interface feature extraction. J. Physiol. 580, 703–721. Whittingstall, K., Stroink, G., Schmidt, M., 2007. Evaluating the spatial relationship of event-related potential and functional MRI sources in the primary visual cortex. Hum. Brain Mapp. 28, 134–142. Fuchs, M., Wagner, M., Kastner, J., 2001. Boundary element method volume conductor models for EEG source reconstruction. Clin. Neurophysiol. 112, 1400–1407. Fuhrmann, A.G., Sun, F., Handwerker, D., D'Esposito, M., Knight, R., 2007. Spatio-temporal information analysis of event-related BOLD responses. Neuroimage 34, 1545–1561. Fuhrmann, A.G., Hein, G., Tsai, N., Naumer, M., Knight, R., 2008. Temporal characteristics of audiovisual information processing. J. Neurosci. 28, 5344–5349. Genovese, C.R., Lazar, N.A., Nichols, T., 2002. Thresholding of statistical maps in functional neuroimaging using the false discovery rate. Neuroimage 15, 870–878. Gloor, P., 1985. Neuronal generators and the problem of localization in electroencephalography: application of volume conductor theoery to electroencephalography. J. Clin. Neurophysiol. 2, 327–354. Goldman, R.I., Stern, J.M., Engel Jr., J., Cohen, M.S., 2002. Simultaneous EEG and fMRI of the alpha rhythm. NeuroReport 13, 2487–2492. Goncalves, S.I., de Munck, J.C., Pouwels, P.J., Schoonhoven, R., Kuijer, J.P., Maurits, N.M., Hoogduin, J.M., Van Someren, E.J., Heethaar, R.M., Lopes da Silva, F.H., 2006. Correlating the alpha rhythm to BOLD using simultaneous EEG/fMRI: inter-subject variability. Neuroimage 30, 203–213. Grave de Peralta, M.R., Gonzalez Andino, S.L., Morand, S., Michel, C.M., Landis, T., 2000. Imaging the electrical activity of the brain: ELECTRA. Hum. Brain Mapp. 9, 1–12. Grave de Peralta, M.R., Murray, M.M., Michel, C.M., Martuzzi, R., Gonzalez Andino, S.L., 2004. Electrical neuroimaging based on biophysical constraints. Neuroimage 21, 527–539. Hillebrand, A., Barnes, G.R., 2005. Beamformer analysis of MEG data. Int. Rev. Neurobiol. 68, 149–171. Kriegeskorte, N., Formisano, E., Sorger, B., Goebel, R., 2007. Individual faces elicit distinct response patterns in human anterior temporal cortex. Proc. Natl Acad. Sci. USA 104, 20600–20605. Laufs, H., Kleinschmidt, A., Beyerle, A., Eger, E., Salek-Haddadi, A., Preibisch, C., Krakow, K., 2003. EEG-correlated fMRI of human alpha activity. Neuroimage 19, 1463–1476. Laufs, H., Holt, J.L., Elfont, R., Krams, M., Paul, J.S., Krakow, K., Kleinschmidt, A., 2006. Where the BOLD signal goes when alpha EEG leaves. Neuroimage 31, 1408–1418. 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. Mantini, D., Perrucci, M.G., Cugini, S., Ferretti, A., Romani, G.L., Del, G.C., 2007. Complete artifact removal for EEG recorded during continuous fMRI using independent component analysis. Neuroimage 34, 598–607. Martuzzi, R., Murray, M.M., Meuli, R.A., Thiran, J.P., Maeder, P.P., Michel, C.M., Grave de Peralta, M.R., Gonzalez Andino, S.L., 2009. Methods for determining frequency- and region-dependent relationships between estimated LFPs and BOLD responses in humans. J. Neurophysiol. 101, 491–502. Michel, C.M., Murray, M.M., Lantz, G., Gonzalez, S., Spinelli, L., Grave de, P.R., 2004. EEG source imaging. Clin. Neurophysiol. 115, 2195–2222. Mumford, J.A., Nichols, T., 2009. Simple group fMRI modeling and inference. Neuroimage 47, 1469–1475. Muthukumaraswamy, S.D., Edden, R.A., Jones, D.K., Swettenham, J.B., Singh, K.D., 2009. Resting GABA concentration predicts peak gamma frequency and fMRI amplitude in response to visual stimulation in humans. Proc. Natl Acad. Sci. USA 106, 8356–8361.