## 1. Introduction

Epileptologists have now at their disposal a variety of tools for investigating human brain functions. Among the technologies of non-invasive functional imaging that have flowered in the last years, two techniques became particularly popular: the electroencephalogram (EEG) which records electrical voltages from the electrodes placed on the scalp and functional magnetic resonance imaging (fMRI) which records magnetization changes due to variations in blood oxygenation.

Each of these methods have its own advantages and disadvantages and no single method is best suited for all experimental and clinical conditions. EEG is a long-established tool for the non-invasive brain investigation characterized by the high temporal resolution (measured in milliseconds) but very low spatial resolution (measured in square centimetres). In contrast, fMRI provides good spatial resolution (measured in square millimetres) but relatively poor temporal resolution (measured in seconds). Combining EEG and fMRI provides integration of information that results in an enhanced view of the phenomena of interest.

This fusion of information is particularly useful in the context of the study of the epileptic disorders. The EEG was used in the study of epilepsy since it was discovered and it remains nowadays the gold-standard for the diagnosis of epilepsy, the classification of the seizures types and the localization of the generators of the epileptic activity. The EEG measurements recorded on the scalp is visually inspected by the neurophysiologists in order to detect any epileptic pattern such as spikes, spike-wave bursts, seizures, etc. and to diagnose the epilepsy. From a spatial point of view only a topographic localization of the generators of ictal and interictal activity is possible. Because of poor spatial resolution of the EEG technique, in many cases it means just a lateralization of the generators and not their precise localization.

fMRI, first demonstrated in 1990, is a technique that, through the blood oxygen level dependent (BOLD) effect, allows the localization of brain areas in which there is a variation of the level of neuronal activity during an experimental condition compared to a control condition. fMRI is mostly used in the study of sensory, motor and cognitive functions, in which the experimental condition differs from the control condition in a way that is controlled by the experimenter. In the context of epilepsy or spontaneous physiological changes in brain state, one can consider the control condition at the time when the EEG is at baseline and the experimental condition to occur in presence of endogenous electrophysiological phenomena such as an epileptic discharge or a sleep spindle. To define such an experimental condition, it is necessary to combine fMRI with EEG measurements by recording EEG while the subject is in the MR scanner.

In this way, it is possible to determine the region of the brain in which there is a change in the BOLD signal correlated to modified neuronal activity such as an epileptic discharge recorded on the scalp..

In this chapter, we will present the methods that have been developed for integrated processing of EEG and BOLD signals simultaneously acquired focusing our attention on the application of these methods in the context of epilepsy.

More in details, the chapter is organized into three sections dealing with the following issues: analysis of EEG-fMRI data by General Linear Model (GLM) theory and estimation of the haemodynamic response function correlated to the epileptic discharges; analysis of EEG-fMRI data by separation into Independent Component Analysis; connectivity analysis of the brain networks underlying epileptic activity.

## 2. Analysis of the haemodynamic response to the epileptic discharges

One of the most popular approaches to combine EEG and fMRI measurements is to include EEG-derived information as regressors of interest in the GLM definition commonly used for fMRI analysis. This type of analysis is usually known as EEG-informed fMRI analysis.

The most widely used approach to perform EEG-informed fMRI analysis consists in analyzing EEG in the time domain and in using information derived from EEG time course into the GLM analysis.

This type of analysis allows for the spatial localization of brain areas involved in spontaneously arising neural activities in which the “stimuli” are not exogenous (i.e. externally generated and controlled) but are “task-free”, random and endogenous (i.e. internally generated).

This is the case of pathological disorders such as epilepsy in which EEG-informed regressors can be used to identify and localize the epileptogenic and the irritative areas. In this particular application, the temporal sequence of onset times of epileptic disharges is used to construct regressors forming a design matrix that is than fitted to the fMRI image data.

One of the most critical factors limitating the potential of the above mentioned technique and preventing it from becoming a clinical tool in the context of epilepsy is the few knowledge about the haemodynamic response to the epileptic spikes. The commonly used EEG-informed GLM analysis needs, in fact, the definition of an “a priori” model of neurovascular system impulse response (HRF, Hemodynamic Response Function).

During the GLM analysis of fMRI data, a statistical comparison is made between the real time course of the BOLD signal in each voxel of the acquired volume and an expected time course of the BOLD signal in the same voxel. The latter is obtained convolving a train of impulses synchronized with the epileptic events and an a priori modeled HRF. Such statistical comparison results in a fMRI map representing the resemblance between real data and modeled data. The whole procedure is summarised and depicted in Fig.1.

A wide variety of basis functions has been used to describe the hemodynamic response. The simplest one uses a standard HRF, which is the measured response to a brief stimulus, such as an auditory tone (Glover et al., 1999). Clearly, this model assumes that each event is followed by a stereotyped response, not accounting for differences among patients, among brain regions or among sessions within a single patient, although it is known that these effects can be considerable (Aguirre et al., 1998). There are also no evidences to assume that the behaviour of the neuro-vascular system in response to an “endogenous” and spontaneous event such as an epileptic discharge, or any other pathological activity, is the same as the one observed in normal subjects, performing sensory or cognitive tasks, in which exogenous and experimentally controlled stimuli are used.

In order to overcome these limitations, several approaches were proposed in literature. In this section we will describe the existing techniques that were developed to study the HRF related to epileptic discharges and how these can be applied in the context of the study of epileptic disorders.

## 2.1 Parametric vs non-parametric estimation of the haemodynamic response function

The techniques to estimate the temporal dynamics and the spatial variability of the HRF can be classified into two main classes: parametric and non-parametric methods.

The parametric methods fix the shape of the HRF by fitting to fMRI data a particular non-linear function of parameters which models typically the delay and the blurring effect of the HRF.

where *t* is the time, *b*_{1,} *b*_{2,} *d*_{1,} *d*_{2} and *c* are the parameters whose range of variations are derived from physiological information.

A hybrid method with parametric and non parametric characteristics at the same time was more recently introduced (Gossl et al., 2001). More in detail, this approach is parametric on the temporal scale and, on the other hand, is not parametric on the spatial scale since it uses only general prior to model the spatial extension of the signal.

The parametric methods are not able to fully capture the shape variations of HRF within the brain, since they unavoidably introduce a bias in the HRF estimate.

In order to overcome such limitations another class of HRF estimation techniques were recently introduced. This novel class of methods are known as non-parametric methods.

Non-parametric techniques for HRF estimation make no prior hypothesis about the shape of the impulse response of the neurovascular system.

The first type of non-parametric method that was applied is the simple averaging over time of the BOLD response. However, this classical voxelwise analysis is precluded by the low signal-to-noise ratio of fMRI data. For these reasons further non parametric methods were proposed: averaging over regions (Kershaw et al., 2000), selective averaging (Dale et al., 1997), introduction of non-diagonal models for the temporal covariance of the noise (Burock et al., 2000), or introduction of smoothing FIR filters (Goutte et al., 2000).

In such a context, recently, a Bayesian, non-parametric estimation of the HRF has been proposed (Marrelec et al., 2003), (Ciuciu et al., 2003) in which information based on the underlying physiological knowledge was used within a Bayesian framework only to temporally regularize the problem. In this way estimates of the HRF are derived without introducing bias into the estimation, since only very soft regularizing constraints, which are clearly derived from physiological requirements, are imposed. By using a Bayesian approach, the HRF estimate results from a tradeoff between information brought by the data and by our prior knowledge (Marrelec et al., 2003), (Ciuciu et al., 2003).

One of the points of strength of this approach is the possibility to extend it from a voxelwise formulation of the problem of HRF estimation to a regional level (Makni et al., 2008). The development of an estimation method of cerebral haemodynamic response at a regional level allows a joint procedure of HRF estimation and detection of active brain areas. The possibility to use a joint-detection estimation approach, avoiding the classical procedure based on GLM, allows to overcome the necessity of an “a-priori” model of HRF.

For this reasons and since this Bayesian non-parametric estimation of HRF has been shown to be able to provide a quantification of the haemodynamic response to epileptic discharges (Tana et al., 2007), hereinafter we will concentrate our attention on this particular approach describing it in details in the following paragraph.

### 2.2. Bayesian non-parametric estimation of the haemodynamic response function

The non-parametric Bayesian HRF estimation approach proposed by (Ciuciu et al., 2003) was first developed in a voxelwise formulation.

The value of the BOLD signal *y*_{tn} at time *(t*_{n}*=nTR) 1 ≤ n ≤ N* (TR = time of repetition of MR sequence) measured in the voxel *V*_{i} is described by the following convolution model:

which in matrix form becomes:

where *X=*[*x*_{tn,} *x*_{tn-Δt}*,..., x*_{tn-Δt}]^{t} is a binary matrix corresponding to the time of occurrence of the events which, in the particular case of epilepsy, can be both interictal and ictal. *h=*[*h*_{0}*, h*_{Δt}*,..., h*_{KΔt}]^{t} represents the HRF that has to be estimated, *Pl*_{i} is the confounding part (deterministic trends) and *b*_{i} is a zero-mean Gaussian white noise process *b* of unknown variance *r*_{b}.

The likelihood function of this model is given by:

where:

Bayesian formalism is then used to model temporal prior information about the structure of the HRF. Since the underlying process of BOLD response to epileptic events is object of investigation, only basic and soft constraints, that do not contradict current knowledge, can be used (Ciuciu et al., 2003), (Buxton et al., 1997).

More precisely, it is possible to assume that the amplitude of the HRF starts and ends at zero and that its variations are smooth, i.e. that the underlying process evolves rather slowly on the experimental time scale. The first condition is introduced by setting the first and last sample points of the HRF to zero. Quantification of the second condition is achieved by setting a Gaussian prior *θ* stands for the hyperparameter that adjusted the prior weight.

The trade-off between constraints and the information given by data and that is modeled by the hyperparameter *θ*, can be set automatically from the data themselves using a Maximum Likelihood Expectation Conditional Maximization algorithm (ECM) (Ciuciu et al., 2003). Also the coefficients of the functions modelling the low frequencies were automatically tuned by using ECM.

The prior and the likelihood function (3) can be fused by using the Bayes rule into the Gaussian a posteriori distribution whose maximum is the Maximum a Posteriori (MAP) estimate of the HRF *h*^{MAP}:

where

Since the BOLD signal is known, to have some spatial structure, this method of HRF estimation can be also extended, in an appropriate region-based HRF model accounting for the spatial dimension of the data (Ciuciu et al., 2004), (Makni et al., 2004).

### 2.3. HRF estimation applied to the study of interictal epileptic activity

In the context of the study of fMRI response in epilepsy several attempts have been made to study HRF to epileptic activity and to take into account both the interregional HRF variability and HRF variability between healthy and epileptic patients.

The majority of the studies existing in literature are devoted to the investigation of the BOLD response to interictal epileptic discharges (IED).

One of the first works of this type used a linear estimation approach based on a simple averaging over the time of the BOLD response (Bénar et al., 2002). This method consists mainly in performing the following three steps: detecting activated areas by GLM analysis, obtaining the time courses of the BOLD signal for each region of interest as the spatial average over the voxels of the ROI and, at the end, averaging the extracted BOLD signal over the time using IEDs identified on EEG signal as events. (Bénar et al., 2002) found important variations in amplitude and shape between average HRFs across patients (see Fig. 5 from (Bénar et al. 2002)) that probably could reflect in part different pathophysiological mechanisms. More particularly, findings of (Bénar et al., 2002) showed that average HRF presented a wider positive lobe than the Glover model in three patients and a longer undershoot in two patients and allowed to conclude that the HRF for epileptic spikes can be somewhat different from the standard model and is also different from patient to patient. Although only few patients were analyzed in the work of (Bénar et al., 2002), it opened the door to the necessity to improve the standard GLM analysis with more complex analysis in order to take into account the variability of HRF and the fact that it does not seem to be a standard response to all type of epileptic discharges and, therefore, to all types of pathophysiological mechanisms underlying epileptic activity.

Following the way opened by (Bénar et al., 2002), HRF variability was studied in a wider group of subjects (Bagshaw et al., 2004) and the use of patient–specific haemodynamic response following a *parametric* approach was proposed (Kang et al., 2003).

In particularl, in (Bagshaw et al., 2004), a group of 31 patients with focal epilepsy were analysed and variations of the peak time of HRF were taken into account within a GLM framework. Using multiple HRFs (with the following four different peak times: 3, 5, 7, 9, s) composed of a single gamma function, resulted in an increased percentage of data sets with significant fMRI activations, from 45% when using the standard HRF alone, to 62.5%.

The standard HRF was found to be good at detecting positive BOLD responses, but less appropriate for negative BOLD responses, the majority of which were found to be more accurately modelled by an HRF that peaked later than the standard.

(Kang et al., 2003), (Lu et al., 2006), (Lu et al., 2007) proposed to detect fMRI areas of activation using canonical HRF, to subsequently estimate HRF and to use patient-specific (Kang et al., 2003) and voxel-specific HRF (Lu et al., 2006), (Lu et al., 2007) within a GLM framework.

Using patient-specific and voxel-specific HRF, the found active regions are characterized by similar or larger volume extent and higher adjusted coefficient of multiple determination (Razavi et al. 2003) than the regions resulting from GLM analysis with fixed HRF.were found (Kang et al., 2003), (Lu et al., 2006),; and additional activated areas compatible with EEG and anatomical MRI localization of epileptogenic and lesional regions were also found (Kang et al., 2003), (Lu et al., 2006), (Lu et al., 2007).

All the above mentioned work (Bénar et al., 2002), (Kang et al., 2003), (Lu et al., 2006), (Lu et al., 2007) support the hypothesis that the misspecification of the form of the HRF may have an important impact on the probability of detecting significant BOLD responses in epileptic patients.

It is worth mentioning also a recent work of (Lemieux et al., 2008) where *parametric* estimation of haemodynamic response performed on a group of 30 patients revealed that the shape of the haemodynamic response is very different from the canonical HRF only remotely from the suspected focus of epileptic activity whereas it appeared similar to canonical one in the vicinity of the presumable epileptogenic areas.

Recently the *non-parametric* Bayesian approach described in details in the previous theoretical paragraph was proposed to study HRF to interictal spikes (Tana et al., 2007).

In the study of (Tana et al., 2007) *non-parametric* Bayesian estimation of HRF was applied to a two patients with temporal lobe epilepsy and HRF variability among patients and among regions was studied.

(Tana et al., 2007) observed important variations in the time course of the haemodynamic response both between patients and across the different fMRI areas of a same subject. In Fig. 2 an example is shown of HRF estimation obtained with non-parametric Bayesian method in a subject with a clinical diagnosis of focal epilepsy shown in Fig. 3. It can be noted as the HRFs in region of interest are different from canonical haemodynamic response function defined in (Glover, 1999). In the areas congruent with the scalp EEG alterations the haemodynamic response has an initial pattern similar to the one obtained in classical event-related experiments but it is followed by an increase of fMRI signal activity away from the event. This should be related to deep epileptic activity that is not recorded on the scalp (Bénar et al., 2002). A negative fMRI response was also detected in an area far from the localization of the interictal EEG spikes.

Non-parametric Bayesian estimation appears to confirm that the shape of the HRF of the epileptic spikes may differ from the standard model and it is variable across regions. The fact that the haemodynamic response could be also widespread and found in distant cortical regions from the ones related to the scalp EEG findings could be an artefact or can suggest an underlying biological process that extends beyond the area clinically assumed as focus and, therefore, can suggest the possibility of effects of focal EEG spikes on remote but synaptically connected regions (Lemieux et al, 2008), (Tana et al., 2007). It could be possible to further investigate this issue by means of invasive EEG recordings or by means of methodologies of analysis that do not need the detection of epileptic events on scalp EEG (see section 3).

Summarizing all the results described above, it seems possible to conclude that the issue of estimation of the haemodynamic response to epileptic activity and, particularly the potential effect on the detection efficiency remains still an interesting subject of investigation.

This opens the way for new type of analysis of fMRI data that are free not only from whatever type of “a priori” model of HRF but also from the necessity to detect epileptic events on scalp EEG recordings. The following section of this chapter will be devoted to the description of these techniques and to how they are able to significantly contribute to the investigation of haemodynamic correlates of epilepsy.

## 3. Independent component analysis of EEG-fMRI data in epilepsy

In contrast to model-based GLM analyses, data-driven techniques applied in fMRI are not constrained by a fixed hypothesis. The most successful one is the Independent Component Analysis (ICA) (McKeown et al., 1998), which separates the data into a large set of independent components showing brain activation patterns with a common time course. Not imposing a-priori models about the shape of the HRF, ICA analysis might find out more information about the BOLD signal.

Let *y* be the M×N (M=number of scans, N=number of time courses) matrix of the fMRI time series, *x* the L×N matrix whose rows *x*_{i} (i=1,…, L) contain the spatial processes (or sources) that generate *y* (L ≤ N). ICA decomposition of fMRI time series is the estimation of the spatial processes by the following linear statistical model:

*y*=

*Bx*+

*v*

where *v* is the noise contribution and *B* is the unknown mixing matrix. Assuming the hypothesis of the number of sources are minor or equal to the number of observed variables (L ≤ N), and considering the *A* matrix to be full rank (Comon, 1994), the problem is reduced to:

*y*=

*A s*

where *s* is the matrix of sources and *A* is the mixing matrix whose columns *A*_{j} (j=1,…,L) contain the time courses of the L processes.

All the spatial components, with the possible exception of one, are assumed to be independent and non-Gaussian (Hyvärinen & Oja, 2000): the non-Gaussianity is, hence, the criterion for the blind estimation of the original sources. Several measures are proposed for applying non-Gaussianity in ICA estimation, even if there is no published evidence about different performances in fMRI analysis. The estimation of the significant components (i.e., related to brain activation patterns) is not based on a comparison with a BOLD model or different shapes of HRFs, but on the maximizing of a non-Gaussianity contrast function (Comon, 1994), therefore on the statistical feature of fMRI data.

The critical point is the separation of the significant components related to an investigated process (i.e., spike-related BOLD responses) from no-significant ones (MR artefacts, default state mode).

The first proposal to overcome this limitation was presented for the analysis of interictal fMRI in focal epileptic patients (Rodionov et al., 2007). After IC decomposition, components selection was implemented by an automatic IC classification (De Martino et al., 2007, see Fig. 4) resulting in the following set of labels: (1) the ‘BOLD’ class, which included components that are thought to consistently reflect task-related, transiently task-related and brain state-related (e.g. default state) neuronal activity; (2) residual motion artefacts; (3) EPI-susceptibility artefacts; (4) physiological noise; (5) noise at high spatial frequency; and (6) noise at temporal high frequency (Rodionov et al., 2007).

The classification was based on a training dataset from healthy volunteers: it was designed for revealing stereotypical components of normal brain activity, misclassified components related to motion, blood vessels, noise, and components related to interictal epileptiform discharges (IEDs). The method was tested on 63 patients with focal epilepsy, who underwent EEG-fMRI recording (Salek-Haddadi et al., 2006). A mean of 16 over 20 ICs were classified as significant BOLD-related sources. Concordance between the ICA and GLM-derived results was assessed based on spatial and temporal criteria on 8 case studies. The remaining ICs were associated to BOLD patterns of spontaneous brain activity, introducing the possibility of an epileptic activity that was not evident on the scalp EEG.

Using this approach, the potential spike-related components were estimated using the corresponding activation in the GLM analysis and the selection was again depending on the canonical HRF used in the GLM (Moeller et al., 2009). Moreover, the classifier was trained on GLM-activation on a very small number of healthy subjects (De Martino et al., 2007) and the training was implemented on only 7% of the testing dataset, decreasing the statistical significance of the model.

(LeVan and Gotman, 2009) introduced a more independent ICA method using deconvolution for identifying component time courses significantly related to simulated focal spikes without constraining the shape of the HRF. Artificial time courses were obtained by generating spikes at random tims and convolving them with a canonical HRF computed from the difference of two gamma functions (Glover, 1999), and varying the location of the activation, the number of simulated spikes per run, and the HRF amplitude. The robustness of traditional analysis methods based on the GLM when the HRF is mis-specified was evaluated adding an additional dataset based on a non-canonical HRF. Fig. 5 shows a schematic representation of the data processing and the obtained results related to the percentage of simulations with constant and variable HRF amplitudes, and the related mean number of falsely activated voxels in concordant components.

Components matching the simulated activation regions were found in 84.4% of simulations, while components at discordant locations were found in 12.2% of simulations; large artefacts occurring simultaneously with spikes are the majority of the false activations. This method mainly depends on the simulation parameters because, when the number of spikes was low, concordant components could only be identified when HRF amplitudes were large.

An application of this method was in detecting dynamic ictal BOLD responses in focal seizures (LeVan et al., 2009), for investigating HRFs with clear peaks - but varying latency – and differentiating the ictal focus from propagated activity. Components related to seizures of 15 patients (suffering of focal ictal discharges or generalised spike and slow waves and sharp and slow wave discharges) were identified by fitting an HRF to the component time courses at the time of the ictal EEG events. HRFs with a clear peak were used to derive maps of significant BOLD responses and their associated peak delay (LeVan et al., 2009). The prominence of the HRF peak was defined considering as baseline the data more than 5 s before or after the peak, and computing the ratio of the peak amplitude to the standard deviation of the baseline. The so-obtained ICA maps were significantly correlated with the GLM maps for each patient (Spearman's test, p<0.05), for a further confirm. The ictal BOLD responses identified by ICA consisted of the presumed epileptogenic zone, but in more widespread area (about 20.3% in addition of the average value). The introduction of an ICs classification method based on the peak delay showed that BOLD response clusters corresponded to early HRF peaks were concordant with the suspected epileptogenic focus, while late HRF peaks to ictal propagation (Fig. 6).

A last attempt of ICA analysis in epilepsy was presented by Moeller et al. (Moeller et al., 2011), where patients with idiopathic generalized epilepsy (IGE) and generalized spike wave discharges (GSW) were studied for the particularity of having hemodynamic behaviour not easily identified with a standard HRF (Moeller et al., 2008). Moreover, GSW are often related with more robust results in regular areas, thus allowing a good comparison between GLM and ICA methods.

After fMRI preprocessing and ICA decomposition (Moeller et al., 2011), ICA time courses were modelled as blocks with the same timings and durations as the GSWs, and convolved with a Fourier basis set (Josephs et al., 1997), assuming that BOLD response to the GSWs could be contained within an interval from 10 s before to 20 s after the marked events, and accommodating the HRF variability. From each temporal IC, a deconvolved HRF was produced by fitting the Fourier basis set. Components significantly related to the GSW were then identified by an F-test (P < 0.05, corrected for the number of components), with the motion parameters used as confounds as in the following GLM analysis. HRFs fitted to the GSW components were then investigated to determine the sign of the HRF peaks, as in (LeVan & Gotman, 2009).

In 12 epileptic patients, comparison of GLM maps and ICA maps showed significant correlation and revealed BOLD responses in the thalamus, caudate nucleus, and default mode areas. Few areas of BOLD signal changes that were only detected by ICA in 8 patients (Fig. 7) and showed variable shapes, different from the canonical HRF in most cases, while one component that was only detected by ICA showed an HRF resembling the canonical HRF.

The particular result is that, in patients with a low rate of discharges per minute, GLM maps detected BOLD signal changes within the thalamus and the caudate nucleus, not present in the ICA ones. Even if these results demonstrated that the BOLD response largely resembles the standard HRF and confirmed the adequacy of GLM analysis, it is worth noting that the number of investigated subjects and their variability in epileptic diseases might compromise the final outcome.

## 4. Brain connectivity and propagation of epileptic activity

One of the challenges of the neurologists in the study of epileptic disorders is the understanding of the propagation of epileptic abnormal activity inside the brain.

The propagation of epileptic seizure and interictal activity is a key concept in epilepsy which indicates the observation of similar patterns or of signals with different patterns but all suspected of reflecting a common underlying phenomenon, on an increasing number of EEG recording channels. (Lemieux et al., 2011).

The techniques used to extract information about interacting brain areas involved in the phenomenon of propagation, can be classified according to two broad typologies of approaches: functional connectivity and effective connectivity. Functional connectivity is defined as the “temporal correlation between spatially remote neurophysiological events” (Friston et al., 1993) and effective connectivity is defined as "the causal influence that a system exerts over one other" (Friston et al., 1993) and it reveals the strength and the direction of the flow of information between fMRI areas.

Functional connectivity can only partly account for the wide variety of the interaction patterns that can be expressed by the effective connectivity, (Friston et al. 1993). The full understanding of the network interaction structure need of information about the directionality of flows provided only by effective connectivity.

The issue of effective connectivity can be approached by two main typologies of analysis techniques: model-based methods (e.g. Dynamical Causal Modeling or DCM) and data-driven methods (e.g. Granger Causality Analysis or GCA). Both approaches try to estimate directed casual influences between cerebral structures by extracting useful information from the temporal dynamics in the EEG and fMRI signals. Both DCM and GCA approaches have advantages and disadvantages and, in the case of the fMRI signal, it is a current open issue, strongly debated in literature, to establish which is the most suitable method for the investigation of connectivity (David 2009), (Roebroeck et al., 2009).

In this section we will review and describe both methods and their current applications on the investigation of the propagation of epileptic activity.

### 4.1. Model-based effective connectivity: Dynamical Causal Modeling

The central idea of DCM is to treat the brain as a deterministic non-linear dynamic system that is subject to inputs and produces outputs. DCM relies a dynamic neuronal model of interacting brain regions, whereby neuronal activity in a given brain region causes changes in neuronal activity in other regions according to a graphical model. A further forward model (the so-called balloon model (Friston et al., 2003)) of the relationships between haemodynamic response and neural activity supplements the above mentioned neural model. The use of balloon model allows to include in the analysis the effect of the haemodynamic convolution and to quantify the strength of the interactions within brain networks directly at neural level. The parameters of the model are then inferred using a Bayesian inference scheme, (Stephan et al., 2009).

The point of strength of DCM is that it is able to model the effect of experimental, external, modulatory inputs on network dynamics and but one of its critical features is the effective ability of the proposed neuro-vascular coupling model to capture the real relationships between blood flow changes and oxygen metabolism changes during activation (Marrelec et al. 2006), (Aubert et al., 2002).

Since DCM takes dynamics and modulations into account in the model, its mathematical framework, which is also able to capture nonlinearities and temporal correlations into account, is very complex and, as a consequence, DCM is computationally limited by the number of regions that can be included in the analysis (maximum of eight according to Penny et al. (2004a); three in Mechelli et al. (2003), Penny et al. (2004b), and Ethofer et al. (2006); three and five in Lee et al. (2006)).

In order to overcome this problem, (Penny et al. (2004b) proposed an extension of the DCM framework to perform model comparison within a set of graphs given a priori. How this approach can be generalized to allow for blind model selection from the whole set of structural models (i.e., with no structural model required a priori) remains a central, yet complex, issue (Marrelec et al. 2006).

### 4.2. Data-driven effective connectivity: Granger Causality Analysis

The most important disadvantage of model-based effective connectivity is that it requires an “a priori” specification of a structural model in the form of a directed graph and that it allows to test hypothesis about connectivity only within a small set of models assumed to be applicable. To overcome this problem GCA was introduced on the basis of the Granger causality concept according to which the activity of a region of interest (ROI)-1 “causes” the activity of a ROI-2 if the knowledge of past values of the ROI-1 time series improves the prediction of the current value of the ROI-2 time-series

GCA does not require any pre-specification or a priori knowledge about the connectivity structure and was successfully applied to fMRI data measuring BOLD response both in bivariate (Roebroeck et al., 2005); (Abler et al. 2006) and multivariate GC models (Deshpande et al., 2008); (Sato et al., 2009); (Deshpande et al. 2009). (Sato et al., 2010), (Havlicek et al., 2010).

Although the choice between GCA and model-based method like DCM is currently debated in literature (David, 2009), (Roebroeck et al., 2009), at the state of art no theoretical reasons exist to exclude the effectiveness of GCA analysis to infer connectivity on BOLD signal (Roebroeck et al., 2009). GCA approach was indeed successfully applied to fMRI data for studying brain connectivity during cognitive (Roebroeck et al., 2005), (Demirci et al., 2009), (Sato et al., 2010), (Ide et al., 2011), (Shippers et al., 2011), (Seger et al., 2011), sensory (Deshpande et al., 2008), (Stilla et al., 2007), (Stilla et al., 2008), (Havlicek et al., 2010), and motor tasks (Abler et al., 2006), (Sato et al, 2006), (Chen et al., 2009) and for investigating resting-state networks, (Liao et al., 2011), (Jiao et al., 2011) and in pathological conditions like epilepsy (Tana et al., submitted).

One of the methodology more widely used to calculate Granger causality are is based on multivariate autoregressive models which are fitted to the signals (EEG or BOLD) of interest.

The multivariate autoregressive (MVAR) recorded from a set of k signals can be expressed as:

where *Y(t)=[Y*_{1}*(t), Y*_{2}*(t),... Y*_{k}*(t)]* is a vector wherein *Y*_{k} is the time series corresponding the *k*-channel, *A(i)* is the matrix of the model parameters, *E(t)* is the vector containing white noise processes and *p* is the model order which can be determined using information theory criteria, (Akaike, 1974).By transforming (7) to the frequency domain, we can obtain the following formulation:

where *H(f)* is the transfer matrix of the model whose element *H*_{ij} represents the connection from the *j*-th to the *i*-th channel.

Normalizing the transfer matrix of the model in (8) with respect to the inflows into channel *i*, we can obtain the Directed Transfer Function (DTF):

where *H*_{ij}*(f*) is the minor produced by removing *i*-th row and the *j*-th column from transfer matrix *H(f)*.

Values of DTF equal to 1 between a pair of channels show maximum direct causal relationships and values of DTF near to 1 indicate that most of the signal in channel *i* consists of signal from channel *j*. DTF values of 0 or close to 0 means an absence of causal interrelations between the two channels considered.

DTF gives information only about the causal relationships underlying the networks of the investigated system and it is not able to discriminate direct transmission for a given pair of regions from indirect propagation of information mediated by other regions. To solve this problem, the estimator dDTF can be introduced, (Korzeniewska et al., 2003). It is given by the following equation:

where *η*_{ij}*(f)* is the DTF normalized over the full frequency band also known as full frequency direct Directed Transfer Function (ffDTF), (Korzeniewska et al., 2003), and *χ*_{ij}*(f)* is the partial coherence defined as in (Korzeniewska et al., 2003). The partial coherence is a measure of how much the interaction between two brain regions is mediated by other areas.

Multiplying DTF by partial coherence, we can combine the information about the directionality of the interactions within the networks with the information provided by partial coherence function and can clearly identify causal direct connections and distinguish them from causal indirect connections.

Another method to identify causal direct connections is the calculation of the Partial Direct Coherence (PDC) introduced by (Baccalà et al., 2001). PDC is a frequency domain representation of the key concept of Granger causality and is defined as:

where *i*, *j*th element and *j*th column of the matrix

PDC_{ij}*(f)* represents the frequency domain GC from the *j*-th time series to the *i*-th time series at frequency *f*.

Significance of both DTF and PDC values can be assessed employing surrogate data (Deshpande et al., 2008), (Korzeniewzka et al., 2003). To this aim, surrogate data can be generated by transforming the data to the frequency domain, randomizing their phases and transforming back to the time domain, (Theiler et al., 1992). The resulting time series have the same power spectrum but random phases with respect to the original signals.

A null distribution can then be obtained by generating a set of sufficient number (for examples 1000 (Deshpande et al., 2008) of surrogate data and calculating the DTF or PDC from these datasets. The DTF or PDC value obtained for each connection from the original time series have then to be compared with the null distribution for a test of significance (e.g. a two-tailed test (Deshpande et al. 2008).

### 4.3. Application to epilepsy

Both interictal and ictal activities can be propagated themselves inside the brain and the progressive recruitment of brain areas far from the origin of the epileptiform activity appears as a spread of pathological EEG pattern.

The causes of the progressive spread of the epileptic activity depends on both the connections existing locally and at a wider range scale and also on the capability of the brain region to be recruited by the abnormal neural activity characteristic of an epileptic discharge (Lemieux et al., 2011).

Sometimes, the pattern of propagation can be studied on the basis of general knowledge of anatomical information and correlated with the temporal evolution of clinical symptoms of the ictal episodes. The majority of the times these general clinical considerations are insufficient to understand the real pattern of propagation which can also vary from event to event (Lemieux et al., 2011) and the question of how an initially spatially localized epileptic focus can spread to involve a larger portion of the cortex remains an open issue.

The identification of the pathways for the spread of partial seizures is, therefore, one of the most relevant and interesting field on which the above described brain connectivity technique can be applied.

The majority of studies existing in literature regarding the connectivity in epilepsy are mainly related to the EEG signal as mean to measure and evaluate neural activity.

Regarding functional connectivity, as reported in (Wendling et al., 2010), the first attempts in the field of epilepsy, were done in the middle of the twentieth century (Barlow and Brazier, 1954), just after the introduction of the Fast Fourier Transform algorithm.

In the 70s the propagation of interictal events was studied by calculating cross-correlation in time domain or equivalently the coherence in the frequency domain using first invasive EEG (Brazier et al., 1972) and later scalp EEG recordings (Lopes da Silva et al., 1977).

Coherence function was then used also to study ictal events, and in particular to study the evolution of the partial seizures and their propagation between the two hemispheres (Gotman et al., 1982). The concept of coherence was extended in 90s in time-varying context in order to study the evolution of the degree of synchronization of interictal and ictal activity (Haykin et al., 1996) and (Franaszcuk et al., 1999).

The methods for studying functional connectivity based on cross-correlation function or on its equivalent form in the frequency domain, that is the coherence function, are based on the assumption that the interaction between EEG signals are linear and, in order to capture also non-linear interactions further techniques were developed. Among these, we can mentioned: mutual information (Mars et al., 1983), non linear regression analysis (Wendling et al, 2001), phase synchronization methods (Rosemblum et al., 2004), generalized synchronization methods (Stam et al., 2003). A recent work of (Wendling et al., 2009) have compared the performance of ten of these methods (which can be grouped into three main families: non linear regression, phase synchronization and generalized synchronization) using simulated data in which the degree of coupling can be controlled. The results obtained by applying the various types of methods to simulated data showed that there is not a “universal” method that performs better than the other ones whatever the considered situation (Wendling et al., 2009), (Wendling et al., 2010).

Regarding effective connectivity, GC has been used to study temporal lobe epileptic seizures in the form of DTF applied to electrocortigram (ECoG) recorded with subdural grids and stereotactic EEG (SEEG) recorded with deep electrodes (Franaszczuk et al. 1998).

(Franaszczuk et al. 1998) extends, in particular, the application of the DTF method to compare patterns of flow of seizures with different sites of origin. Analysis of a seizure originating from mesial temporal structures is compared with a seizure originating from lateral temporal neocortex. The DTF method has the potential to determine patterns of flow of activity, including periods when visual analysis of the intracranial ictal EEG may not allow for definitive source localization.

More recently GC was also applied in the form of PDC to scalp EEG signal in order to study temporal lobe epilepsy (Baccalà et al., 2004).

In particular PDC was applied to mesial temporal epileptic seizures and it has been shown that it is able to localize the epileptogenic focus via the simultaneous analysis of multiple EEG channels thanks to the determination of the direction of information flow among signals and thanks to its representation by means of directed graphs, where focal electrodes are associated with high observed rates of pertinence to strongly connected subgraphs.

It is possible to affirm that both DTF and PDC are able to study effective connectivity in epilepsy, to infer epileptic seizure propagation and to identify the focus of epileptic activity (Cadotte et al., 2009).

Differently from EEG, nowadays there are still few application of fMRI effective connectivity studies in the field of epilepsy.

Effective connectivity studies using DCM have been performed in (David et al., 2008) and (Vaudano et al., 2009).

DCM and Baysian model comparison was used in (Vaudano et al., 2009) to investigate the role of thalamus, prefrontal cortex and precuneus in seizure generation. EEG-fMRI data recorded in a group of seven patients with idiopathic generalized epilepsy (IGE) with frequent generalized SW discharges (GSWD) and significant GSWD-correlated haemodynamic signal changes in the thalamus, the prefrontal cortex and the precuneus.

In order to perform Bayesian model selection three dynamic causal models were constructed: GSWD was modelled as autonomous input to the thalamus (model A), ventromedial prefrontal cortex (model B), and precuneus (model C). Bayesian model comparison revealed that, although model C (GSWD as autonomous input to precuneus) is the best in five patients while model A (GSWD as autonomous input to thalamus) prevailed in two cases, at the group level model C dominated.

The findings lead the authors to hypothesize a role for the precuneus as a form of modulator of generalized SW activity, and by extension, of the occurrence of absence seizures, linking spontaneous fluctuations in brain state as reflected by the so-called Default-Mode Network of brain activity to the occurrence of epileptic discharges (Vaudano et al., 2009), (Lemieux et al., 2011).

In a rat model of absence epilepsy (David et al., 2008) performed simultaneous EEG and fMRI measurements, and subsequent intracerebral EEG (iEEG) recordings in regions strongly activated in fMRI (first somatosensory cortex, thalamus and striatum). (David et al., 2008) showed that using DCM, instead of GCA, it is possible to spatially localize the origin of spontaneous spike-and-wave (SW) discharges in the first somatosensory cortex.

More recently, a study of epileptic seizure propagation using GCA was performed in (Tana et al, submitted). In this study, GCA analysis was applied to networks of brain areas showing fMRI activation during epileptic seizures. EEG-fMRI recording was performed on a group of four case studies related to patients with different epileptic pathologies and GCA analysis was applied to networks obtained using different parcellation strategies for the definition of the nodes. In Figure 8, the connectivity pattern of a patient with a diagnosis of occipital lobe epilepsy is shown. The patient shows an area of almost continuous spiking under left occipital lobe revealed by ECoG measurements (not shown in figure) and fMRI activation mainly localized in the left occipital and right parietal lobe. GCA results showed that the source of the connectivity network is the left occipital fusiform gyrus (LOFG) confirming intracranial EEG findings.

The role of LOFG as possible starting point of seizure propagation and, therefore, the ability of GCA to recognize LOFG as epileptogenic focus is also confirmed by the fact that the area identified as GC network source is located within the brain volume removed in successful surgery for epilepsy (Fig. 8b).

According to GC results, the seizure starting in the left occipital lobe propagates to the contralateral occipital, parietal and temporal regions. The sink of the network (i.e. the area in which the propagation of the seizure terminates) is localized in the contralateral temporo-parietal area.

## 5. Conclusion

This chapter summarizes various techniques for integrated analysis of EEG and fMRI signals for the investigation of epilepsy. In particular, we describe how the information coming from the EEG can be used for triggering the analysis of the BOLD signal, in order to establish a direct connection between EEG events and haemodynamic activation, and for a precise localization of the brain areas involved in the specific events. This type of analysis is limited by the current little knowledge about the HRF that is necessary to define a model of expected BOLD signal to construct GLM regressors. In order to improve the understanding of the impulse response function of neurovascular system, several techniques have been developed for the estimation of the HRF temporal and spatial characteristics.

To overcome the limitations of classical GLM approach that are also due to the fact that EEG describes only cortical activity and is not always able to detect deep brain events, other methodologies have been proposed in literature, such ICA, calculated on the BOLD signal only.

After the localization of the activated cortical areas during epileptic seizure, it is of extreme interest to study the temporal interactions inside the network formed by activated brain regions and in order to identify the pattern of seizure propagation. Applying connectivity techniques like GCA and DCM to both EEG and, more recently, BOLD signal it is possible to obtain useful information about the localization of the epileptogenic focus, the propagation of the seizure and sometimes to localize also how the seizure terminates.

Even if some problems are still open and several features are under investigation, all the methodologies described in this chapter are promising tools that will allow a deeper investigation of the mechanisms involved in epileptic activity and can found large application in the clinical field as well as in research.