Bayesian multi-modal model comparison: a case study on the generators of the spike and the wave in Generalised Spike-Wave complexes.

We present a novel approach to assess the networks involved in the generation of spontaneous pathological brain activity based on multi-modal imaging data. We propose to use probabilistic fMRI-constrained EEG source reconstruction as a complement to EEG-correlated fMRI analysis to disambiguate between networks that co-occur at the fMRI time resolution. The method is based on Bayesian model comparison, where the different models correspond to different combinations of fMRI-activated (or deactivated) cortical clusters. By computing the model evidence (or marginal likelihood) of each and every candidate source space partition, we can infer the most probable set of fMRI regions that has generated a given EEG scalp data window. We illustrate the method using EEG-correlated fMRI data acquired in a patient with ictal generalized spike-wave (GSW) discharges, to examine whether different networks are involved in the generation of the spike and the wave components, respectively. To this effect, we compared a family of 128 EEG source models, based on the combinations of seven regions haemodynamically involved (deactivated) during a prolonged ictal GSW discharge, namely: bilateral precuneus, bilateral medial frontal gyrus, bilateral middle temporal gyrus, and right cuneus. Bayesian model comparison has revealed the most likely model associated with the spike component to consist of a prefrontal region and bilateral temporal-parietal regions and the most likely model associated with the wave component to comprise the same temporal-parietal regions only. The result supports the hypothesis of different neurophysiological mechanisms underlying the generation of the spike versus wave components of GSW discharges.


Introduction
During the past decade, EEG-correlated fMRI has been used to map haemodynamic changes correlated with the occurrence of epileptiform discharges in focal and generalised epilepsies (Salek-Haddadi et al., 2003;Salek-Haddadi et al., 2006;Laufs 2007;Gotman et al., 2006).In focal epilepsy, the localising information thus obtained has been found to be broadly concordant with the location of the sources inferred from other electro-clinical data and EEG source reconstruction (Lemieux et al., 2001;Benar et al., 2006).In contrast, generalised epileptiform discharges, such as generalised spike-wave (GSW) complexes seen on the EEG traces, are commonly associated with widespread haemodynamic changes in the neocortex and sub-cortical regions (Laufs et al., 2006;Hamandi et al., 2006).
Given the temporal resolution of fMRI, of the order of a few seconds, one must assume that these maps are representative of haemodynamic changes taking place over entire discharges.As a consequence, they are not informative with respect to the electrophysiological processes that underlie the different components of such GSW1 complexes.On the other hand, the problem of localizing these processes from the EEG remains a difficult challenge due to the widespread nature of the pattern which is thought to reflect rapidly propagating neural activity over a large part of the cortex.This is in addition to a fundamental difficulty of the EEG inverse problem namely that the underlying current sources cannot be estimated uniquely from EEG scalp measurements without invoking priors or constraints.EEG/MEG inverse solutions differ in the nature of their priors, which should ideally be specific to the neuronal phenomenon at hand.For instance, assuming that the underlying active network consists of a few focal sources has been used to justify equivalent current dipoles (ECD) localisation methods (Sherg and Ebersol, 1993).ECD solutions applied to the analysis of focal epileptic spikes have proven relatively concordant with both fMRI statistical maps (see e.g.Korjenova 2001) and intracranial electrode recordings (Merlet and Gotman, 1999).However, this class of inverse solutions can be misleading in presence of spatially extended sources (see e.g.Kobayashi et al., 2005).In contradistinction, distributed linear (DL) methods aim at estimating the amplitude of a predefined highly dense ensemble of dipoles, typically spread over the cortical sheet (Dale and Sereno, 1993).Usually, additional spatial and/or temporal constraints are used to finesse the underdetermination of DL inverse solutions (see e.g.Daunizeau et al., 2007b).In this context, fMRI-derived spatial priors are thought to significantly improve the spatial resolution of DL inverse solutions, particularly if the inverse technique is able to account for the potential mismatch between EEG and fMRI sources (see e.g.Daunizeau et al., 2008 for a comprehensive review of EEG/fMRI information fusion).
Recently, such probabilistic DL source reconstruction methods using spatial priors derived from fMRI statistical parametric maps (SPMs) have been developed (see e.g.Daunizeau et al., 2005;Sato et al., 2005;Daunizeau et al., 2007a).Most of these techniques fall into a Bayesian framework, which provides an estimation of the current sources (posterior probability maps or PPMs, see e.g.Friston et al., 2007) and allows for generic model comparison, through the computation of the model evidence/marginal likelihood (see e.g.Trujillo-Barreto et al., 2005or Mattout et al., 2006).
Source reconstruction using DL and ECD models from high-density EEG have already been applied to GSW, suggesting a focal frontal origin for the spike and broader frontal generator for the wave (Holmes et al., 2004;Tucker et al., 2007), broadly in line with earlier work (Lemieux and Blume 1983;Niedermeyer and Lopes da Silva, 2004).This encouraging preliminary result, along with previous work on both EEG-correlated fMRI data analysis and fMRI-constrained EEG source reconstruction, has led us to develop a principled probabilistic technique dedicated to identifying the respective networks involved in the respective generation of the spike and wave components of GSW complexes from EEG and fMRI data.
In Daunizeau et al., 2005, we have proposed a Bayesian model comparison scheme for assessing the relevance of fMRI-derived spatial priors in probabilistic EEG source reconstruction.In Grova et al., 2008, we have applied this method in the context of focal interictal spike localization, using EEG-correlated fMRI statistical parametric maps (SPMs).
In this work, we systematize and extend this method to cope with multiregions EEG-correlated fMRI SPMs, by means of the Expectation-Maximization (EM) algorithm (see e.g.Friston et al., 2008).The method is designed to take advantage of the spatial resolution of fMRI and the temporal resolution of EEG, and consists of the three following steps: Deriving the static EEG-correlated fMRI SPM.Due to the low temporal resolution of fMRI compared to EEG, the resulting map of activations can be taken to reflect multiple aspects of the EEG events of interest such as the generators of the spike and the wave for the specific case of GSW.
(ii) Building the cortical source space partitions that are associated with every combination of the activated regions.
(iii) Using Bayesian model comparison to identify the most likely source space partition with respect to the spike and the wave component of the EEG scalp measurement of GSW complexes, respectively.
We demonstrate the potential of this method in an analysis of multi-modal EEG-fMRI data acquired in a patient affected by idiopathic generalized epilepsy (IGE) with frequent absence seizures.

Patient clinical history
We studied a right-handed 23 year-old man (written informed consent and ethics committee approval obtained) affected by juvenile absence epilepsy (JAE; onset age 12y) with frequent absence seizures (2-3 episodes per week) and rare (fewer that one per year) generalized tonic clonic seizures.He was born by a caesarean section three weeks preterm; he was well at birth and his developmental milestones were within normal limits.He has no history of febrile convulsions, brain injury or other risk factors for the development of epilepsy.His father's mother was diagnosed with temporal lobe epilepsy at the age of 58.
Neurological examination and morphological MRI scans (at 1.5T and 3T) were normal.
Previous EEG recordings showed a normal background interrupted by 2.5-3 per second generalized spike-and-wave activity with anterior predominance, facilitated by hyperventilation.There was no response to photic stimulation.
The patient was treated with AED, in mono-therapy or association (Ethosuximide, Lamotrigine, Levetiracetam and Valproate) without achieving complete control of the absences.
At the time of our investigations, the patient was taking Levetiracetam 2500mg/day and Valproate 1000mg/day.The EEG showed frequent spontaneous and hyperventilation-related generalized spike-and-wave (GSW) discharges lasting between <1 and 20 s; clinically, the longest discharges (more than 15 seconds) were accompanied by psychomotor arrest and eyelid blinking.
The patient was asked to rest with his eyes closed and to keep still.Two 20minute fMRI sessions were acquired.A high-resolution T1-weighted scan was also acquired.

EEG data pre-processing
The BrainVision software package (http://www.brainproducts.com/)was used to correct the EEG traces from both gradient-and pulse-related artefacts.
Both the MRI gradient and the cardiac pulse artefact removal algorithms are based on artefact template subtraction (the latter using the ECG signal for modelling the artefact template).
For the purpose of the fMRI analysis, the onset and offset times of GSW discharges were marked and recorded on the artefact-corrected EEG traces.
For the purpose of the EEG source reconstruction, we analysed those events in an "event-related potential" (ERP) fashion, as follows: The GSW typical discharges lasted for a few seconds (see the GSW discharges lasting for 18 seconds on Figure 1) and contained a series of spike-wave events.We identified two different event types, namely (i) the spike and (ii) the slow wave.We manually selected a representative spike and slow wave, using the SPM8 EEG data review functionality, as templates for detection of spikes and waves throughout the recording.A new list of events of each type was then obtained for all data window in the EEG traces whose correlation coefficient with each of the templates exceeded 0.95.Finally, we averaged the detected events, obtaining typical ERP responses.

fMRI analysis
The FMRI data were processed and analysed using the SPM5 software package (www.fil.ion.ucl.ac.uk/spm/).After discarding the first four image volumes, the EPI time series were realigned, and spatially smoothed with a cubic Gaussian Kernel of 8 mm full width at half maximum and normalised to MNI space.
A general linear model was used to assess the presence of regional GSWrelated BOLD changes.The marked GSW events were represented as variable-duration blocks from GSW onset to cessation (block design).
Motion-related effects were included in the GLM in the form of 24 regressors representing the 6 scan realignment parameters and a Volterra expansion of these, plus Heaviside functions for large motion effects (Friston et al., 1996;Salek-Haddadi et al., 2006;Lemieux et al., 2007).An additional set of confound regressors was included to account for pulse-related signal changes (Liston et al., 2006).The model is based on an over-complete basis set expressing a linear relationship between cardiac-related MR signal and the phase of the cardiac cycle and has been validated anatomically and its effect on efficiency of the estimation of the effects of interest.It is generally considered good practice to model as many confounding effects as possible, leading to increased confidence in the results [Lund et al, 2006], which may be particularly important in studies of individual patients.
The GSW event blocks were convolved with the canonical hemodynamic response function (HRF), its temporal derivative (TD) and dispersion derivatives (DD), to form regressors testing for GSW-related BOLD changes.
Significant BOLD signal changes correlated with GSW were assessed using an F-contrast across the three regressors of interest.The resulting SPM was thresholded at p< 0.05, corrected for multiple comparisons (Friston et al., 1991).

Building empirical fMRI-derived priors for the EEG inverse problem
We performed fMRI-constrained source reconstruction for the averaged spike and wave, respectively, on the cortical surface, taken to be the canonical cortical mesh provided by the SPM software package.
According to the DL framework, each EEG dataset y is assumed to be generated from a linear mixture of d dipoles of unknown amplitudeθ , whose positions and orientations are those of the vertices of the SPM canonical cortical surface.Within that framework, the locations of the underlying extended sources are defined by those connected set of vertices (spatial components) which have a significant activity.This means that prior knowledge about the position of the underlying extended sources is translated into higher prior activity power for the corresponding spatial components.As fully detailed in (Daunizeau et al., 2005), we can use an increased prior variance over these spatial components to cast fMRI-derived source location knowledge within a Bayesian treatment of the DL framework.In this work though, we slightly depart from this perspective and associate an unknown variance hyperparameter on the prior variance of each of these spatial components.The Bayesian probabilistic generative model m is then fully specified by the number and composition of spatial components of the source space.It is these spatial components and the ensuing model space we want to explore.This means we want to identify the combination of spatial components that is the more plausible with respect to the measured scalp EEG data.This is important, because we a priori do not know which subset of fMRI clusters have generated the EEG data.From a Bayesian perspective however, this is simply a matter of model comparison: we can use the model evidence to identify the source space partition that is the most likely to have generated the EEG data.
First, we interpolated the 3D volumic thresholded fMRI SPM on the canonical cortical surface of the SPM software.The interpolation kernels were based on Voronoï cells centred on each cortical mesh vertex (see e.g.Grova et al., 2006), and constrained to lie within the limits of a 3D-volumic grey matter mask.Then, we identified the connected components of the thresholded fMRI activation map on the cortical manifold.To do so, we applied standard mathematical morphology (closing and erosion on the cortical manifold, see e.g.Soille 1999) to obtain K anatomically connected clusters around each local maximum of the interpolated SPM.This furnished a set of K spatial components with compact support, each components corresponding to an active cluster extracted from the fMRI SPM.
Having obtained K cortical patches, we then build the 2 K generative models m (c) ( c = 1,..., 2 K ) that correspond to each and every combination of these clusters.These models contain from one to K spatial components, each of which is associated with a diagonal covariance component, having non-zero elements only for dipoles belonging to the corresponding cluster (see Appendix I).In addition to the fMRI-derived spatial components, each generative model m (c) includes two "whole-brain" prior covariance components, i.e. (i) the identity matrix (yielding the standard "minimum norm" DL source reconstruction algorithm) and (ii) the discrete Laplacian operator (yielding the "maximum smoothness" inverse solution).
This parameterization of the models m (c) assumes that the structure of cortical activity is composed of a sum of independent spatial processes, i.e.
both a smooth and a rough active field (which are spread over the whole cortical surface) and a set of patchy sources (with bounded spatial support) whose power profile is given by the fMRI activation score2 (Daunizeau et al., 2005).
We refer the interested reader to the Appendix 1 for more details about the construction of the Bayesian probabilistic generative model that associates weighted prior covariance components to the above fMRI-derived and "wholebrain" spatial components.
The contribution of each of these processes to the underlying cortical activity structure is unknown a priori, and is estimated from the data.We have done this using the ReML (restricted Maximum Likelihood) algorithm of the SPM software package, which is a stand-alone MATLAB code that aims at inverting this class of generative models (see e.g.Mattout et al., 2006, Friston et al., 2007).In brief, we use ReML to estimate covariance hyperparameters at both the sensor and source levels, yielding both an estimate of the cortical sources and an approximation to the model evidence ( ) The latter is then used for comparing the different source partitions m (c) , and to derive the best model source subset for the spike and the slow wave of GSW, respectively.This means we invert (using ReML) the 2 K generative models m (c) ( c = 1,..., 2 K ) that correspond to the different combinations of spatial components, and identify the more plausible source partition ( *) c m in terms of its model evidence ( ) p y m , as approximated by ReML for both datasets (spike and slow wave).
We refer the reader to the Appendix 2 for a simulated experiment highlighting the expected properties of the proposed probabilistic approach.

EEG-correlated fMRI results:
Good quality EEG was obtained after off-line artefact subtraction.During the first EEG-fMRI session one typical 2.5-3 Hz GSW discharge, 18 seconds long, was recorded (see Figure 1).
The rest of the EEG showed a normal awake background with continuous 9-Hz posterior alpha rhythm.Following the experiment, the patient said that he had a "seizure" during the scanning session.

fMRI-constrained EEG source reconstruction:
K = 7 anatomically connected BOLD clusters were obtained (see Figure 2), which were used for fMRI-constrained EEG source reconstruction.We then calculated the posterior probabilities of each of the 2 K = 128 source space partitions (consisting of all different combinations of the 7 BOLD clusters), conditional upon the spike and the slow wave, respectively.The results are shown in Figure 3, and a summary table ( First, we note that the non fMRI-constrained EEG source reconstructions prove significantly less probable than the best fMRI-constrained source reconstructions.However, the non fMRI-constrained inverse solution was more likely than most of the other fMRI-constrained inverse solutions ( F = −170.8for the spike and F = −168.4for the slow wave).This was the case for both the spike and the slow wave components of the GSW discharge.This is important, because the plausibility of the best fMRIconstrained inverse solution cannot be explained only by the reduction of the effective degrees of freedom (due to the prior spatial constraint).This means that the best fMRI-derived source space partition is likely to have generated the EEG scalp measurement (see Daunizeau et al., 2005).Second, among the 128 models tested, the posterior probability distribution over source space partitions significantly identifies one best model for both the spike and the slow wave components of the GSW complex.For both the spike and the wave components, the best model contains the left and right middle temporal gyri.
However, the spike component most likely source space partition also contains the left medial frontal gyrus (the second best model adds the right frontal medial gyrus, see Figure 3).In other terms, the frontal activity present during the spike seems to be inhibited during the slow wave component of the ictal GSW discharges.
For completeness, Figure 4 shows the non fMRI-constrained source reconstructions for both the spike and the slow wave.In summary, the level of matching of the non fMRI constrained source reconstruction with the most plausible sets of fMRI regions is questionable, but there seems to be a similar Figure 3).

Discussion
In this work we presented a novel approach to identify the generators of brain activity captured using simultaneous EEG-fMRI based on probabilistic Bayesian EEG source model comparison using fMRI-derived regional priors.
Application of the method to the spike and wave components during an ictal GSW discharges demonstrated a different origin of the two components.That is, both the spike and the wave components were generated by bilateral temporal-parietal cortex activity, but the left medial frontal gyrus source (indentified during the spike) disappeared during the following slow wave of the GSW complex.This is consistent with early involvement of the ventromedial frontal cortex during the spike discharges of absence seizure in line with previous studies (see Tucker et al 2007, Holmes et al 2004).
To the best of our knowledge, this work demonstrates the first application of Bayesian multimodal EEG-fMRI modelling to the fine spatio-temporal characterization of the neural correlates of generalized epilepsy ictal activity.
The method combines asymmetrical EEG-correlated fMRI statistical analysis and fMRI-constrained probabilistic EEG source reconstruction.Note that these results could not have been obtained using standard (non fMRIconstrained) source reconstruction, which was significantly less plausible than the best fMRI-constrained inverse solutions.This is very likely to be due to the under-determination of the EEG inverse problem, which causes the source estimates to be highly uncertain (which is taken into account by Bayesian model likelihood measures).In comparison, fMRI-constrained inverse solutions show a lower degree of uncertainty, at the cost of constraining the solution to resembling the fMRI profile of activation.Bayesian model comparison works because Bayesian model likelihood quantifies the potential conflict between the prior and the likelihood.In other terms, the best fMRIderived source partition confirms the spatial information that can be extracted from the EEG data, which makes it more plausible than the non fMRIconstrained inverse solution.
In general, the precision of these results, in terms of effective spatio-temporal resolution, can then be thought of as a consequence of a well-balanced combination of the respective spatio-temporal resolutions of EEG and fMRI.
However, one cannot expect the same spatial resolution to hold for both fMRI analysis and EEG source reconstruction (even constrained by fMRI spatial information).This is because there are numerous confounds that intrinsically limit the spatial resolution of any source reconstruction technique.First of all, the definition of the spatial model (the so-called gain matrix, see Appendix) relies on a number of well-known approximations, e.g.imperfect spatial realignment of the electrodes, crude geometrical model of the tissue conductivity, potential misspecification of the position and orientation of the distributed dipoles on the cortical sheet.These approximations contribute to the loss of spatial resolution that could theoretically be achieved by any source reconstruction technique.Secondly, the inverse technique itself is limited by its sensitivity to the underlying prior assumptions and by the level of measurement noise corrupting the data.In the general case, it is very difficult, if not impossible, to quantify the expected spatial resolution of source reconstruction (see e.g.Baillet, Riera et al. 2001or Darvas 2004).When using fMRI spatial information, we argue that the minimum requirements of good practice are (i) to use fMRI clusters that match the expected spatial resolution, (ii) to construct a test statistics that accounts for the spatial uncertainty.A weaker form of (i) is simply met by homogenizing the size of the clusters to the average fMRI cluster size, which we did.The latter requirement is more difficult, since we cannot quantify the expected spatial precision.
However, the Bayesian marginal likelihood accounts at least for the uncertainty arising from the inverse problem difficulties.Since the added spatial uncertainty arising from the forward problem is identical for all compared models, it should only lead to an overconfident (as opposed to biased) model comparison.This means that the test statistics are correct in average, but artificially inflate the evidence in favour of the more plausible model.
It is also worth mentioning the potential difficulties related to the practical implementation of this method.Among them, we found very difficult to derive a fully automated pre-processing step, for both the EEG scalp data and the interpolation and cortical manipulation (dilatation/erosion) of fMRI clusters.
Concerning the effect of quality degradation of the EEG recorded during scanning on the inverse solution, we have recently demonstrated the validity of source estimation based on EEG data recorded during fMRI using the same artefact correction methodology as employed in this study 3 [Vulliemoz   et al, 2009].Furthermore, the fMRI activation map does not have the same topology when considered in 3D-volumic (native) space or in 2D-surfacic (cortical) space.So far, the default mathematical morphology applied on the cortical surface teased apart the different parts of the interpolated fMRI clusters that were covering opposite sides of a sulcus.These highly nonlinear operations can nonetheless be sensitive to the topology of the fMRI activation map in its 3D-volumic (native) space4 .These effects also interact with the actual definition of the spatial covariance components of the generative models (see Appendix 1).On the whole, we believe that improvements in preprocessing should lead to increased robustness of the proposed methodology.
Also, we did not discuss so far the potential influence of "missed" sources, i.e.
of source that could underlie the EEG data but are not part of the set of activated fMRI clusters.To our knowledge, the probabilistic framework we propose seems to be partly robust to potential missed sources, in the sense that we expect it to favour the "whole-brain background" (no cluster) model whenever the missed source is clearly expressed in the data (see simulated experiment in Appendix 2).This reproduces the results highlighted in Daunizeau et al., 2005 andin Phillips et al., 2005, that both use similar hierarchical Bayesian techniques.However, it should be noted that this is a matter of sensitivity, in the sense that the fMRI-missed source could still potentially be missed if its contribution to the EEG data is weak.Conducting a comprehensive analysis of the sensitivity of this method is beyond the scope of the present article, but we expect the method to be further extended and assessed in future publications.
Although limited to a single case, the result supports the hypothesis of different neurophysiological mechanisms underlying the generation of spike versus wave components of GSW discharges.Earlier electrophysiological studies in generalized penicillin epilepsy of the cat indicated that the spike corresponds to short periods of increased cortical excitation whereas the wave component comes from longer-lasting periods of intense cortical inhibition (see Gloor, 1977).
Our results suggest that left medial frontal region contains hyperexcitable neurons.This cortical hyperexcitation was confirmed by a transcranial magnetic stimulation study (see Gianelli et al. 1994)  Density during 25 absences and demonstrated independent fields for the spike and slow wave, the former anterior (prefrontal and fronto-polar) and the latter posterior (parieto-occipital).Lemieux and Blume (1983) reported that the "slow wave were more diffuse, more symmetrical distribution and more posteriorly centred than either spikes and troughs." The neurophysiological mechanisms to explain the findings presented in this study remain speculative.It is of interest that the identified network only includes two major sources (frontal and temporal-parietal).This seems to support, as previously reported (Holmes et al., 2004;Tucker et al., 2007), that absences are not truly generalized but they involved selective cortical networks.Moreover, our results are consistent with abundant evidence, both in animals and humans (see Bancaud J et al., 1974;Lüders H et al., 1984;Amor et al., 2005) suggesting a primary role of ventromedial frontal cortex in the generation of GSW, although our analysis does not address causality directly.Lastly, our findings are coherent with the putative role of the prefrontal cortex during absences5 (Pavone and Niedermayer, 1995).
In conclusion, our novel approach to EEG-fMRI data fusion has proved useful in identifying fMRI-derived EEG source partitions more plausible than the non fMRI-constrained inverse solutions.Model comparison within the Bayesian framework has allowed us to identify specific and distinct generator partitions for the spike and wave components of the 3Hz spike-wave complex.
Appendix I: Bayesian modelling of the EEG inverse problem.
In this section, we describe the probabilistic (generative) model through which we introduce fMRI-derived spatial priors in the EEG inverse problem, and its inversion.
We start with a 2-levels hierarchical linear model of EEG data Y ∈°n × s over n channels and s samples: where is a known gain or lead-field matrix and unknown source dynamics at d dipoles.In the following, the gain matrix has been computed according to a three-sphere conductor head model (De Munck 1988), given the known electrode positions on the scalp and the canonical cortical mesh of SPM.The electrode positions were realigned to the canonical space using standard anatomical landmarks (nasion, left hear and right hear fiducials) whose 3D positions with respect to the electrodes were also given.The subsequent realignment accuracy is of the order of the centimeter.This is not critical, since this is bellow the expected spatial resolution of the probabilistic source reconstruction method.The relevance of these spatial precision concerns to the proposed EEG-fMRI fusion approach is discussed further in the discussion section.

The terms
) Note that the derivatives required to actually implement this EM scheme simplify greatly in our case, because the data are a linear function of the source intensities θ (see e.g.Mattout et al., 2006).
In the context of fMRI-constrained EEG source reconstruction, we associate each of the interpolated clusters k C ( 1,..., k K = ) with a diagonal covariance component, having non-zero elements only for dipoles belonging to the corresponding cluster, i.e.: ( ) ( ) We also added two standard covariance components, namely the identity matrix and the discrete Laplacian operator ∆ , to further regularize the inverse solution (see e.g.Mattout et al., 2006).
All these covariance components then enter the source prior in the generative model required by the EM scheme above.Any combination of fMRI clusters (i.e.any source space partition) will then be associated with a posterior probability (on model space), given each and every dataset (spike and wave): where i j F is the free energy of the i th source space partition, conditional on the j th dataset (j=1: GSW spike, j=2: GSW slow wave).
Both free-energies and posterior probabilities are given in Figure 3.

Appendix II: simulated experiment
In this section, we describe a simulated experiment we designed to inspect the properties of the proposed Bayesian model comparison method.This simulated experiment is a simplification of the evaluation strategy that was developed in Daunizeau et al. 2005, in which we assessed the response of the statistical framework to various perturbations of the fMRI-derived location prior.
We randomly picked two spatially extended sources over the cortical surface, and generated scalp EEG data (SNR = 10) corresponding to a mixture of these two sources (see Equation 1).The gain matrix was similar to that used for the real EEG data whose analysis is described in the main body of this manuscript.The two clusters and the corresponding EEG scalp data can be seen on Figure 5.
We then constructed four models, namely: cluster 1, cluster 2, clusters 1+2, and no cluster.We then applied the proposed method to two distinct data time windows (see Figure 5): The first data window was expressing the only contribution of cluster 1, whereas the second data window only contained signal originating in cluster 2.
Figure 6 shows the results of the Bayesian model comparison for the first data window.As expected, the more probable model was the model containing only the first cluster.Moreover, the model containing only the second cluster was less plausible than the model containing no cluster.This means that the method preferred the whole-brain background component alone to a model containing a spatial component that was not found in the scalp data (cluster 2).This is important, since the no fMRI prior (whole-brain background model) would have been favoured to the wrong fMRI-derived prior (cluster 2 + wholebrain background).Note also that the estimated cortical activity under the whole-brain background model spreads over the cluster 1 and extends further apart from it, due to the nature of the regularization (see Figure 6).In contradistinction, the estimated cortical activity under the model containing cluster 1 (and the whole-brain background) is sharp and focused on cluster 1, again due to the nature of the regularization (note also the difference in terms of maximum power of the two solutions).This means that the cortical activity is much better estimated under the latter model.
Finally, Table 2 summarizes the results for both datasets (1 and 2): It can be seen that for whatever the dataset, the model containing clusters 1 and 2 is always the second most plausible model, and that the "whole-brain background" model (no cluster) is always second least plausible.In addition, the true model (cluster 1 for dataset 1 and cluster 2 for dataset 2) is always favoured.Also, the wrong model (cluster 2 for dataset 1 and cluster 1 for dataset 2) is always the least probable of all models.This form of double dissociation basically means that: if the spatial components include the true source, the true source (and only the true source) is favoured; -if the spatial components do not include the true source, the "wholebrain background" model is favoured.
These results both reproduce and extend the findings reported in Daunizeau et al. 2005, since the method reported in the above reference boils down to the relative comparison of the full (clusters 1+2, with identical associated hyperparameters) model to the "whole-brain background" (no cluster) model.
However, we are now in a position to assess the specificity of the method with respect to each of the individual fMRI clusters, which could not have been done before.Bottom-left: the SPM canonical cortical mesh, onto which the 3D-volumic SPM was interpolated.
Top-right: the SPM field interpolated on the cortical surface, shown on the inflated SPM canonical cortical mesh.
Bottom-right: The 7 fMRI clusters that were used for fMRI-constained EEG source reconstruction, after morphological operations (thresholding, closing and erosion of the thresholded field).Top-right: the two best source space partitions for the spike component of the GSW complex.These two solutions only differ by the laterality of the active medial frontal gyrus.
Bottom-right: the best source space partitions for the slow wave component of the GSW complex.This set of bilateral temporal-parietal clusters is common to both he spike and the wave components.Top-right: the estimated cortical activity power under the most plausible source space partition (cluster 1 only).
Bottom-right: the estimated cortical activity power under the "whole-brain background" (no cluster) model.This is basically similar to a standard minimum norm estimate (MNE) of cortical activity.
trend in terms of the specific activation of the orbitofrontal cortices (only during the spike).First, the activity spreads over almost the whole cortical sheet, which is due to the nature of the regularization.Second, both the spike and the slow wave reconstruction exhibit patterns of activity in the left and right anterior temporal lobes and around the left and right prefrontal cortices.These patterns might be partially explained by the two (left and right) parieto-occipital fMRI sources.Third, the spike source reconstruction shows a pattern of activity on the left and right orbitofrontal cortices, which corresponds to the two frontal fMRI clusters (most likely partition and 2 d most likely partition, in patients with IGE, where motor evoked potentials were recorded simultaneously with the spike and the wave component of the GSW complex.While a decrease stimulation threshold has been found when the stimulus was time-locked to the spike, a threshold increase has been documented during the slow wave.This data confirms the experimental results reporting a wave component corresponding to long periods of inhibition.Besides the different electrophysiological genesis, several GSW source analyses have shown distinctive brain regions implied in the generation of the component spike versus the component slow wave.Studying the poly-spike and wave complex of patients affected by juvenile myoclonic epilepsy (JME), Rodriguez and colleagues (2002) identified a bilateral current source in the medial frontal gyrus corresponding to the spike and multiple sources in different cortical regions corresponding to the wave.Interestingly, as in our case, temporal lobe sources were observed when the slow wave were analysed.Advanced methods of EEG source analysis have been applied to identify the brain regions involved in generation of spike versus slow wave within GSW discharges.Recently,Tucker and colleagues   (2007)  applied advanced methods of electrical source analysis to dense array EEG (256-channel) recordings of GSW in 5 patients with absence spells.They demonstrated a highly stereotyped localization of the spike component within the midline frontal cortex in all cases.Similar results were obtained in 5 patients with IGE (seeHolmes et al., 1994)  using a different source analysis method.Both studies also demonstrated a wide and symmetrical frontotemporal network engagement during the slow wave component of the GSW complex.We found involvement of the posterior temporal regions in relation to the slow wave.The conclusions of other authors on the electrical mapping of GSW are in line with our result.Ernst Rodin (1999) mapped Current Source square root of the interpolated SPM F-score at the i th vertex of the SPM canonical cortical mesh.Using the local fMRI significance score allows accounting for modulations in the spatial profile of the spatial components.However, alternative definitions of the spatial components could be proposed, including flat spatial profile (no spatial weighting) or local smoothness (using a finite support Laplacian operator).To our knowledge, this choice does not significantly impact on the final model comparison results.

Figure 1 .
Figure 1.Spike and slow wave: scalp EEG data.Left: 32-channels EEG recorded inside the scanner after MR artefact subtraction, showing the prolonged (20 sec) generalized spike and wave discharge (3 Hz).Top right: coregistrated GSW complexes (channel P7) and subsequent evoked (average) response.Bottom right: respective scalp topologies of the spike (left) and slow wave (right) evoked responses at their respective peak.Note the overall (reversed) similarity of the components, probably indicating a partial common set of generators.

Figure 3 .
Figure 3. Spike and slow wave: results.Top-left: Free-energies (EM approximation to the log-Bayesian model likelihood, see Appendix I) of the 127 different source space partitions derived from each fMRI cluster combination (blue: spike, red: slow wave).The stars indicate the significantly most plausible fMRI-derived source partition.For the

Figure 4 .
Figure 4. Spike and slow wave: non fMRI-constrained estimates Top row: top, bottom, front and back views of the "whole-brain background" (non fMRI-constrained) source reconstruction of the spike component.Bottom row: top, bottom, front and back views of the "whole-brain background" (non fMRI-constrained) source reconstruction of the slow wave component.The estimated cortical activity power under the "whole-brain background" (no cluster) model is basically similar to a standard minimum norm estimate (MNE) of cortical activity.

Figure 6 .
Figure 6.Simulated experiment: results for data window 1 Top-left: Free-energies (EM approximation to the log-Bayesian model likelihood, see Appendix 1) of the 3 different source space partitions derived from each cluster combination (cluster 1, 2 and 1+2).The stars indicate the significantly most plausible source partition.The green dashed lines show the free-energy of non fMRI-contrained inverse solution.Note that, despite being significantly less likely than the most plausible fMRI-derived source partitions, the non fMRI-constrained inverse solutions still win over the worst source partition (cluster 2 only).Bottom-left: Posterior probabilities of the 3 different source space partitions derived from each cluster combination.These are basically normalized model likelihoods (see Appendix 1).The stars indicate the significantly most plausible source partition.
Figure 1 Click here to download high resolution image

Table 2 :
log-model evidences of each pair of model and dataset.