Thèse Dirigée Par

Christian Jutten

Hamid Soltanian-Zadeh

M Fabrice Wendling

Rapporteur M Mohammad

Bagher Shamsollahi

Examinateur Rapporteur M Habib Benali

M Gholam

Examinateur Ali Hossein-Zadeh

M Christian

Co-Directeur De Thèse

M Hamid Soltanian-Zadeh

Intégration de signaux EEG et IRM pour l'identification des régions cérébrales actives fondée sur la décomposition de sources parcimonieuses

L'électroencéhalographie (EEG) est une technique d'imagerie cérébrale non invasive importante, capable d'enregistrer l'activité neuronale avec une grande résolution temporelle (ms), mais avec une résolution spatiale faible. Le problème inverse en EEG est un problème difficile, fortement sous-déterminé : des contraintes ou des a priori sont nécessaires pour aboutir à une solution unique. Récemment, l'intégration de signaux EEG et d'imagerie par résonance magnétique fonctionnelle (fMRI) a été largement considérée. Les données EEG et fMRI relatives à une tâche donnée, reflétent les activités neuronales des mêmes régions. Nous pouvons donc supposer qu'il existe des cartes spatiales communes entre données EEG et fMRI. En conséquence, résoudre le problème inverse en EEG afin de trouver les cartes spatiales des sources EEG congruentes avec celles obtenues par l'analyse de signaux fMRI semble être une démarche réaliste. Le grand défi reste la relation entre l'activité neuronale électrique (EEG) et l'activité hémodynamique (fMRI), qui n'est pas parfaitement connue á ce jour. La plupart des études actuelles reposent sur un modèle neurovasculaire simpliste par rapport á la réalité. Dans ce travail, nous utilisons des a priori et des faits simples et généraux, qui ne dépendent pas des données ou de l'expérience et sont toujours valides, comme contraintes pour résoudre le problème inverse en EEG. Ainsi, nous résolvons le problème inverse en EEG en estimant les sources spatiales parcimonieuses, qui présentent la plus forte corrélation avec les cartes spatiales obtenues par fMRI sur la même tâche. Pour trouver la représentation parcimonieuse du signal EEG, relative à une tâche donnée, on utilise une méthode (semi-aveugle) de séparation de sources avec référence (RSS), qui extrait les sources dont la puissance est la plus corrélée à la tâche. Cette méthode a été validée sur des simulations réalistes et sur des données réelles d'EEG intracrânienne (iEEG) de patients épileptiques. Cette représentation du signal EEG dans l'espace des sources liées à la tâche est parcimonieuse. En recherchant les fonctions d'activation de fMRI similaires à ces sources, on déduit les cartes spatiales de fMRI très précises de la tâche. Ces cartes fournissent une matrice de poids, qui impose que les voxels activés en fMRI doivent être plus importants que les autres voxels dans la résolution du problème inverse en EEG. Nous avons d'abord validé cette méthode sur des données simulées, puis sur des données réelles relatives à une expérience de reconnaissance de visages. Les résultats montrent en particulier que cette méthode est très robuste par rapport au bruit et à la variabilité inter-sujets.

iii

Dedicated to Hermin and Helia

The important thing in science is not so much to obtain new facts as to discover new ways of thinking about them. Sir William Lawrence Bragg (1890 -1971) widely considered.

The EEG and fMRI data of a specific task, eventually reflect the neurological events of the same activation regions. Therefore, we expect that there exist common spatial patterns in the EEG and the fMRI data. Therefore, solving the EEG inverse problem to find the spatial pattern of the EEG sources which is congruent with the fMRI result seems to be close to the reality. The great challenge is the relationship between neural activity (EEG) and hemodynamic changes (fMRI), which is not discovered by now. Most of the previous studies have used simple neurovascular model because using the realistic model is very complicated. Here, we use general and simple facts as constraints to solve the EEG inverse problem which do not rely on the experiment or data and are true for all cases. Therefore, we solve the EEG inverse problem to estimate sparse connected spatial sources with the highest correlation with the fMRI spatial map of the same task. For this purpose, we have used sparse decomposition method. For finding sparse representation of the EEG signal, we have projected the data on the uncorrelated temporal sources of the activity. We have proposed a semiblind source separation method which is called reference-based source separation (R-SS) and extracts discriminative sources between the activity and the background. R-SS method has been verified on a realistic simulation data and the intracranial EEG (iEEG) signal of five epileptic patients.

We show that the representation of EEG signal in its task related source space is sparse and then a weighted sparse decomposition method is proposed and used to find the spatial map of the activity. In the weighted sparse decomposition method we put fMRI spatial map in the weighting matrix, such that the activated voxels in fMRI are considered more important than the other voxels in the EEG inverse problem.

We validated the proposed method on the simulation data and also we applied the method on the real data of the face perception experiment. The results show that the proposed method is stable against the noise and subject variability.

Introduction

Integrating information across multiple neuroimaging modalities during the same task has been considered tremendously during the past decade. Each modality has its advantages as well as limitations. The dream of integrating multiple modalities is to exploit the common as well as unique information from complementary modalities to achieve better estimation of the localization and dynamics of neurological activity. In practice, a trade-off exists between the estimation accuracy of the common information estimation and the estimation of all activities involved in any of the modalities. In this work, electroencephalography (EEG) and functional magnetic resonance imaging (fMRI) integration is studied.

The recording of electrical potentials from electrodes applied to the human scalp surface is known as electroencephalography (EEG). The human electroencephalogram was first recorded by Berger in the 1920's and since then has been widely used in both research and clinical applications. The main sources of scalp EEG signals are post-synaptic cortical currents associated with synchronously activated pyramidal cortical neurons, which are oriented perpendicular to the cortical surface [1]. EEG detects directly the rapid electrical change of neurophysiological processes, which allows studies of the dynamics of neural networks that occur on the order of tens of milliseconds. Although, scalp EEG recording is noninvasive and cheap, its localization is very complicated [2]. Its complexity originates from its weak and blurred signal and limited number of observations. The scalp EEG recordings are blurred because of the different electrical conductivity of the brain, skull and scalp. To overcome this problem, realistic head models obtained from magnetic resonance imaging (MRI) are used to model the propagation of the potential from the cortex to the scalp. Also, high resolution EEG recording technologies are introduced to record EEG with high number of electrodes (a cap sample is shown in Figure 1.1). Despite the efforts on increasing the accuracy of EEG modeling and recording, it suffers from low spatial solution due to its limited number of spatial measurements and ambiguity of electromagnetic inverse problems [3].

Functional magnetic resonance imaging or functional MRI (fMRI) is a functional neuroimaging procedure using MRI technology that measures neural activity indirectly through oxygena- tion blood flow and metabolism changes. This technique relies on the fact that cerebral blood flow and neuronal activation are coupled. Since the 1890s it has been known that changes in blood flow and blood oxygenation in the brain (collectively known as hemodynamics) are closely linked to neural activity. When an area of the brain is in use, blood flow to that region also increases. The primary form of fMRI uses the Blood-oxygen-level dependent (BOLD) contrast, discovered by Seiji Ogawa. fMRI data provide a measure of brain function on a millimeter spatial scale and a subsecond (and delayed) temporal scale. EEG and fMRI integration is studied intensively and more than other possibilities, because combining the complementary temporal and spatial resolutions of EEG and fMRI may reveal high resolution reconstruction of the spatiotemporal structures of the neural activity. Comparison of spatiotemporal resolution of different imaging method is shown in Figure 1.2. Another reason is that simultaneous acquisition of EEG and fMRI is possible. The main question here is: As EEG and fMRI detect different aspects of neural activity, what is the relationship between these two observation sets? Reports in the literature do not provide a clear picture of the link between EEG and BOLD signals. Figure 1.3 shows a schematic of EEG-fMRI relation. In particular, contradictory results have been presented regarding the dependency of BOLD changes on the EEG power and spectral profiles. The co-existence of visible and invisible sources in a same behavioral task for the E/MEG and fMRI techniques are discussed in [4,5,6,7].

Different models of the transfer function between EEG and BOLD signals are compared in [8], in the prediction of the fMRI data, in a visual stimulation experiment with human healthy subjects. The models explored included the EEG total power (TP; [9]); linear combinations of the power from different frequency bands [10]; and several variations of a heuristic model proposed by [11] in which BOLD changes are assumed to be proportional to the root mean square frequency (RMSF) of the EEG spectrum. The results obtained showed a clear superiority of the RMSF metrics in predicting the BOLD signal, when compared to power-weighted metrics.

Their results show that changes in BOLD are indeed associated with changes in the spectral profile of neural activity and they do not arise from one specific spectral band.

Although a great amount of both experimental and theoretical work has been dedicated to the clarification of the relationship between neural activity and associated hemodynamic changes, neurovascular coupling mechanisms remain an active area of research [8]. The principal limitations on multimodal EEG-fMRI integration are imposed by: � Physiology -One reason why EEG and fMRI sources might be dislocated is the distance between the neuronal population whose electrical activity is generating the EEG signal and the vascular tree, which provides the blood supply to these neurons [12].

-Some source activity may be located or oriented such that there is little electromagnetic field outside the head. Examples of this are radially oriented sources in MEG and deep 'closed field' sources in EEG, for which the activity patterns are such that the total macroscopic current cancels out. All of these could generate significant fMRI but they are not visible in EEG or MEG.

-The fMRI activations can be detected where there are no neuronal activities because the fMRI signal is sensitive to parameters reflecting energy consumption. In practice, the brain consumes energy for many more processes which are not directly linked to the neuronal activities, e.g. neurotransmitter release and uptake, vesicular recycling, maintenance of membrane potentials and so on [13,14]. These kinds of sources are usually referred to as 'fMRI extra sources' or 'EEG invisible sources'.

-Some E/MEG sources cannot be detected in fMRI, which have been usually referred to as 'EEG extra sources' or 'fMRI invisible sources'. Some neuronal sources which are active only for a short time period may be detected in EEG or MEG, but do not appear in fMRI results since fMRI integrates brain activity over time [15].

� Experimental limitations

Another important potential source of bias in EEG-fMRI integration is experimental variability. In some situations, it might be necessary to acquire the EEG and fMRI data in separate sessions. In this case, habituation effects, variations in the stimulation paradigm, or any other difference between the sessions might lead to the differential activity of neural networks [5,16].

Simultaneous scalp EEG and fMRI

The desire of epileptologists to better spatial localization of electrical epileptic sources in patients undergoing presurgical evaluation motivated the idea of EEG-fMRI integration [17]. While, the first EEG monitoring during echo planar MRI were accomplished in 1993 [18] and now, it is generally used in presurgical evaluation of patients with epilepsy.

Developing methods for the integrative analysis of simultaneously recorded EEG and fMRI data is difficult because of the artifacts caused by each method on the other one. (i) As Faraday's law says: a time varying magnetic field in a wire loop induces an electromotive force (EMF) proportional in strength to the area of the wire loop and to the rate of change of the magnetic field component orthogonal to the area. Strong ambient magnetic field of fMRI recording instrument induces a significant EMF on the EEG electrodes. This artifact which is caused by the switching of the magnetic field gradients used in the image acquisition is called "scanning artifact". This is usually the largest in amplitude (in the order of mV) but the most stable over time [START_REF] Allen | A Method for Removing Imaging Artifact from Continuous {EEG} Recorded During Functional {MRI}[END_REF]. (ii) Even a slight motion of the EEG electrodes within the strong static field of the magnet can induce a significant EMF [START_REF] Hill | E{EG} during {MR} imaging: differentiation of movement artifact from paroxysmal cortical activity[END_REF][START_REF] Kruggel | Hemodynamic and electroencephalographic responses to illusory figures: recording of the evoked potentials during functional {MRI}[END_REF]. For instance, native pulsatile motion related to a heart beat yields a ballistocardiographic artifact in the EEG that can be roughly the same magnitude as the EEG signals themselves [18,[START_REF] Goldman | Acquiring simultaneous {EEG} and functional {MRI}[END_REF]. This artifact is called "pulse artifact".

These artifacts are a real concern for concurrent EEG-fMRI acquisition. However, hardware and algorithmic solutions may be able to unmask the EEG signal from MR disturbances.

� Scanning artifact can be removed efficiently by subtracting an average artifact template [START_REF] Goncalves | Artifact removal in co-registered {EEG}/f{MRI} by selective average subtraction[END_REF][START_REF] De Munck | Novel artefact removal algorithms for co-registered EEG/fMRI based on selective averaging and subtraction[END_REF]. The effectiveness of this method is related to the time variations of the MR artifact waveform which can reduce the success of this method [START_REF] Cohen | Simultaneous {EEG} and f{MRI} made easy[END_REF][START_REF] Cohen | Method and apparatus for reducing contamination of an electrical signal[END_REF]. The problem can be resolved through hardware modification that increases the precision of the synchronization of MR and EEG systems [START_REF] Anami | Stepping stone sampling for retrieving artifactfree electroencephalogram during functional magnetic resonance imaging[END_REF] or during post-processing by using precise timings of the MR pulses during EEG waveform averaging [START_REF] Cohen | Simultaneous {EEG} and f{MRI} made easy[END_REF].

� Spectral domain filtering, spatial laplacian filtering, PCA, and ICA (see [START_REF] Sijbers | Restoration of {MR}-induced artifacts in simultaneously recorded {MR}/{EEG} data[END_REF][START_REF] Bonmassar | Motion and ballistocardiogram artifact removal for interleaved recording of {EEG} and {EP}s during {MRI}[END_REF][START_REF] Garreffa | Real-time {MR} artifacts filtering during continuous {EEG}/f{MRI} acquisition[END_REF][START_REF] Negishi | Removal of time-varying gradient artifacts from {EEG} data acquired during continuous f{MRI}[END_REF][START_REF] Srivastava | I{CA}-based procedures for removing ballistocardiogram artifacts from {EEG} data acquired in the {MRI} scanner[END_REF][START_REF] Mantini | Complete Artifact Removal for EEG Recorded During Continuous fMRI Using Independent Component Analysis[END_REF][START_REF] Amini | MR Artifact Reduction in the Simultaneous Acquisition of EEG and fMRI of Epileptic Patients[END_REF]) can be used to remove the artifacts.

� Interleaved EEG/fMRI: the experiment protocol consists of time blocks and only a single modality is acquired during each time-block [START_REF] Bonmassar | Spatiotemporal brain imaging of visual-evoked activity using interleaved {EEG} and f{MRI} recordings[END_REF][START_REF] Makiranta | B{OLD}-contrast functional {MRI} signal changes related to intermittent rhythmic delta activity in {EEG} during voluntary hyperventilation-simultaneous {EEG} and f{MRI} study[END_REF]. This means that every stimulus has to be presented at least once per modality. To analyze ERP and fMRI activations, the triggered fMRI protocol can be used with every stimulus presentation so that EEG and MR are sequentially acquired in order to capture a clean E/MEG signal followed by the delayed HR [START_REF] Sommer | Combined measurement of event-related potentials ({ERP}s) and f{MRI}[END_REF].

� Electrode leads are also important equipments. Creating electrode leads which can have direct contact with the subject's scalp to prevent the development of nuisance currents is important. Also, EEG equipments should not disturb the homogeneity of the magnetic field and distort the resulting MR images [18,[START_REF] Lazeyras | Functional {MRI} with simultaneous {EEG} recording: feasibility and application to motor and visual activation[END_REF]. These artifacts can be effectively avoided by using specially designed EEG equipment [START_REF] Goldman | Acquiring simultaneous {EEG} and functional {MRI}[END_REF]: specialized geometries, and "MR-safe" materials (carbon fiber, plastic) for the leads.

Why do we still work on separated EEG and fMRI acquisition?

Despite advances in simultaneous EEG-fMRI hardware and software, the signal to noise ratio (SNR) of these signals is still significantly lower than the corresponding unimodal paradigms.

Although artifact removal algorithms have been somehow successful [7], some EEG information may be removed during artifact removal step. Also, each modality still imposes limitations on the other, i.e. not any EEG paradigm can be run with fMRI or inside an MR scanner, and not every fMRI paradigm is suitable for EEG analysis. The conduction of separate EEG and fMRI sessions should hence always be considered. If an experiment can be repeated more than once with a high degree of reliability of the data, separate E/MEG and fMRI acquisition may be appropriate.

Co-registration

In EEG-fMRI integration, a common geometrical framework has to be derived in order to locate appropriately the voxels whose EEG responses is high and voxels whose hemodynamic response is increased/decreased during the task performance. The issue of deriving a common geometrical framework for the data obtained by different imaging modalities is called the "co-registration" problem. Several techniques have been used to produce an optimal match between the realistic head reconstruction obtained in the high resolution EEG/MEG by the MRIs of the experimental subject and the fMRI image coordinates. The first group of techniques is based on the presence of landmarks on the both images used for the coregistration. Corresponding landmarks have to be determined in both modalities [START_REF] Hill | E{EG} during {MR} imaging: differentiation of movement artifact from paroxysmal cortical activity[END_REF]. A second group of techniques is based instead on the matching of surfaces belong to the same head structure, as obtained by the different image modalities. In these techniques a prerequisite is the segmentation of the structures whose surfaces have to be matched [16]. With the volume-based registration technique no additional information as landmarks or surface detection is necessary [START_REF] Wells | Multi-modal volume registration by maximization of mutual information[END_REF].

In the case in which the multimodal EEG and fMRI is performed simultaneously, the setup of a common geometrical framework becomes simpler. In this case registration can be performed based on a scanner coordinate system.

Method

The main goal of EEG and fMRI integration is to fuse the observed data to find out the common data which are hidden in both. However, the most appropriate way in which to combine the data in order to achieve this goal is not clear. The majority of previous EEG-fMRI studies have used aspects of the EEG time-series to inform the fMRI data analysis within the framework of the generalized linear model (GLM) [START_REF] Debener | Trial-by-trial coupling of concurrent electroencephalogram and functional magnetic resonance imaging identifies the dynamics of performance monitoring[END_REF][START_REF] Eichele | Assessing the spatiotemporal evolu-tion of neuronal activation with singletrial event-related potentials and functional MRI[END_REF][START_REF] Goldman | Single-trial discrimination for integrating simultaneous EEG and fMRI: identifying cortical areas contributing to trial-to-trial variability in the auditory oddbal ltask[END_REF], called EEG-informed fMRI method, while others have used the opposite strategy of informing EEG source analysis using spatial constraints extracted from fMRI data [START_REF] Dale | Spatiotemporal mapping of brain activity by integration of multiple imaging modalities[END_REF][START_REF] Esposito | Distributed analysis of simultaneous EEG-fMRI time-series: modeling and interpretation issues[END_REF]. These methods are called fMRI-informed EEG methods. Some recent developments have been directed towards the generation of common forward models for both modalities to be applied to the bimodal space-time series [START_REF] Brookings T Ortigue | Using ICA and realistic BOLD models to obtain joint EEG/fMRI solutions to the problem of source localization[END_REF][START_REF] Daunizeau | Symmetrical eventrelated EEG-fMRI information fusion in a variational Bayesian framework[END_REF]11,[START_REF] Sotero | Biophysical model for integrating neuronal activity, EEG, fMRI and metabolism[END_REF][START_REF] Valdés-Sosa | Model driven EEG/fMRI fusion of brain oscillations[END_REF]. These methods are called symmetric methods or EEG-fMRI fusion.

The main ideas of different groups of EEG-fMRI integration, their advantages and their limitations are described in this section. [START_REF] Huster | Methods for Simultaneous EEG-fMRI: An Introductory Review[END_REF][START_REF] Jorge | EEG-fMRI integration for the study of human brain function[END_REF] review EEG-fMRI integration methods. All EEG-fMRI integration methods are shown in a schematic in figure 1.4.

EEG-informed fMRI Methods

In this method, temporal features extracted from the EEG signal is used to form one or more regressors in the generalized linear model (GLM) of fMRI (Figure 1.5). The obtained E/MEG features first get convolved with a hypothetical hemdynamic response function (HRF) and are then subsampled to fit the temporal resolution of fMRI.

An open question in this context is which EEG features are ideally suited for this type of analysis [START_REF] Bagshaw | Single trial variability of EEG and fMRI responses to visual stimuli[END_REF][START_REF] Ostwald | An information theoretic approach to EEG-fMRI integration of visually evoked responses[END_REF]. Several features of the EEG are suitable, such as ERP amplitudes [START_REF] Debener | Trial-by-trial coupling of concurrent electroencephalogram and functional magnetic resonance imaging identifies the dynamics of performance monitoring[END_REF], ERP latencies [START_REF] Bénar | Singletrial analysis of oddball event-related potentials in simultaneous EEG-fMRI[END_REF], EEG synchronization and phase coherence [START_REF] Mizuhara | Long-range EEG phase synchronization during an arithmetic task indexes a coherent cortical network simul-taneously measured by fMRI[END_REF][START_REF] Jann | {BOLD} correlates of {EEG} alpha phase-locking and the fMRI default mode network[END_REF], or the power within specific EEG frequency bands [START_REF] Scheeringa | Trial-by-trial coupling between EEG and BOLD identifies networks related to alpha and theta EEG power increases during working memory maintenance[END_REF].

Event-related potentials (ERPs) are computed from scalp recordings by computer averaging epochs (recording periods) of EEG time locked to repeated occurrences of sensory, cognitive, or motor events. The spontaneous background EEG fluctuations, which are typically random relative to when the stimuli occurred, are averaged out, leaving the event-related brain potentials.

The ERP thus reflects, with high temporal resolution, the patterns of neuronal activity evoked by a stimulus.

EEG power in a specific frequency is also used as an EEG feature. In [START_REF] Goldman | Simultaneous EEG and fMRI of the alpha rhythm[END_REF] the correlation between the BOLD signal and the variations of the alpha wave is investigated in a rest data.

In this study the alpha wave power is extracted from occipital electrodes and are convolved with HRF, then the correlation between the result and the BOLD signal is computed and high correlation is reported for some regions.

The stability of these methods are investigated in different studies against various factors. [START_REF] Laufs | Electroencephalographic signatures of attentional and cognitive default modes in spontaneous brain activity fluctuations at rest[END_REF], various subjects and experiment sessions [START_REF] Gonçalves | Correlating the alpha rhythm to BOLD using simultaneous EEG/fMRI: inter-subject variability[END_REF] and physiological factors [START_REF] De Munck | The Relation between Alpha Band Power, Heart Rate and Fmri[END_REF] are studied. Independent component analysis (ICA) is also used to extract single (EEG) trial variability as predictors for fMRI signal changes [START_REF] Debener | Trial-by-trial coupling of concurrent electroencephalogram and functional magnetic resonance imaging identifies the dynamics of performance monitoring[END_REF][START_REF] Eichele | Assessing the spatiotemporal evolu-tion of neuronal activation with singletrial event-related potentials and functional MRI[END_REF]. [START_REF] Ostwald | An information theoretic approach to EEG-fMRI integration of visually evoked responses[END_REF] explicitly states quantitative degrees of dependency.

Stability against artifacts

In special cases, special approaches are considered. For example, to study epilepsy, after artifact correction, the epileptiform activity is identified by an expert on the EEG traces. These events are then convolved with a hemodynamic response function (HRF), and used as a regressor in the standard GLM [START_REF] Bagshaw | analysis of the EEG-fMRI response to prolonged bursts of interictal epileptiform activity[END_REF]. To study resting-state network spontaneous fluctuations of power in specific frequency bands are quantified in the EEG signal. Time-dependent power in each of these frequency bands is used as the feature. Some multivariate extensions of this type of approach have also been proposed. Essentially, these techniques try to find a linear decomposition of fMRI data which is correlated with a time-frequency decomposition of the EEG [START_REF] Martínez-Montes | Concurrent EEG/fMRI analysis by multiway Partial Least Squares[END_REF][START_REF] Moosmann | Joint independent component analysis for simultaneous EEG-fMRI: principle and simulation[END_REF]. In [START_REF] Martino | Predicting EEG Single Trial Responses with Simultaneous fMRI and Relevance Vector Machine Regression[END_REF] a multivariate machine learning method based on the Relevance Vector Machine (RVM) regression is proposed to consider the multivariate nature of the fMRI data. In [START_REF] Leite | Transfer function between EEG and BOLD signals of epileptic activity[END_REF] a transfer function between EEG and BOLD signal is used.

fMRI-informed EEG Methods

The most straightforward way to impose the fMRI constraint upon the EEG source localization is to restrict the source spaces to locations exceeding a threshold predetermined for fMRI statistical analyses [START_REF] George | Mapping function in the human brain with magnetoencephalography, anatomical magnetic-resonance-imaging, and functional magnetic-resonance-imaging[END_REF]. The main goal of these methods is to use fMRI-derived spatial priors in the EEG source reconstruction problem to achieve high spatial resolution and at the same time utilize the fast dynamics of neuronal activity measured by EEG. This approach can be divided into two classes, associated with the EEG source model employed: (i) the equivalent current dipole (ECD) model [START_REF] Kiebel | variational inversion of the equivalent current dipole model in {EEG/MEG}[END_REF] and (ii) the distributed source model (known as imaging methods) [START_REF] Friston | Multiple sparse priors for the m/eeg inverse problem[END_REF].

In the ECD model, few current equivalent dipoles are assumed to generate the EEG signal.

The spatial information of the fMRI data is used to set the position of the current dipoles at the fMRI hotspots, allowing to the localization procedure to rotate freely only the direction of the dipoles to fit the EEG or MEG data [4,16,[START_REF] Bledowski | Mental chronometry of working memory re-trieval: a combined functional magnetic reso-nance imaging and event-related potentials approach[END_REF].

Because of invisible fMRI sources, it may not be possible to simply restrict the source model solutions to areas where fMRI shows activation. To overcome the displacement between locations obtained from fMRI analysis and ECD modeling, [5], in the first step, compares the detected activations across the two modalities to increase the reliability of dipole localization alone. In [16], additional weighting by the distance from the ECD to the corresponding fMRI activation foci is introduced to guide ECD optimization. Further, fMRI invisible activations can be accommodated by introducing free dipoles without the constraint on dipole location. In addition, many ECDs are required to model spatially extended regions correctly [START_REF] Shiraishi | Application of Magnetoencephalography in Epilepsy Patients with Widespread Spike or Slow-Wave Activity[END_REF].

In imaging techniques (distributed source model), weighted regularization techniques are used to enter the fMRI information in the EEG inverse problem. In these techniques, the model uses the fMRI activation as a prior on the spatial profile of cortically-distributed sources. As a consequence, the estimation penalizes sources whose fMRI-derived activation probability is low. This approach has been shown to estimate the position and extent of underlying sources robustly; whenever the fMRI-derived constraints are congruent with the reality [START_REF] Mantini | Multimodal Integration of fMRI and EEG Data for High Spatial and Temporal Resolution Analysis of Brain Networks[END_REF]. Hence, although approaches for fMRI-informed EEG are used on simultaneously acquired data, their strengths are most apparent in the context of separate recordings. This approach is very sensitive to some generators of EEG or MEG signals that are not detected by fMRI, which have been usually referred to as 'fMRI invisible sources'. In this case, the EEG source reconstruction is strongly biased, which is why many variants of the fMRI penalty term have been proposed [START_REF] Liu | Spatiotemporal imaging of human brain activity using functional MRI constrained magnetoencephalography data: Monte Carlo simulations[END_REF][START_REF] Babiloni | Multimodal integration of EEG, MEG and fMRI data for the solution of the neuroimage puzzle[END_REF][START_REF] Babiloni | Multimodal Imaging from Neuroelectromagnetic and Functional Magnetic Resonance Recordings[END_REF][START_REF] Ou | Multimodal functional imaging using fMRI-informed regional EEG/MEG source estimation[END_REF]. Most of the previous works assumed BOLD signal as the power of EEG, but [START_REF] Brookings T Ortigue | Using ICA and realistic BOLD models to obtain joint EEG/fMRI solutions to the problem of source localization[END_REF] interpret the BOLD signal as the � 1 norm of the EEG source signal. [START_REF] Liu | Spatiotemporal imaging of human brain activity using functional MRI constrained magnetoencephalography data: Monte Carlo simulations[END_REF] revealed that the distortion by the fMRI invisible sources could be reduced considerably by just giving a constant weighting factor to the diagonal terms of source covariance matrix in a linear inverse operator. W f M RI = w ii , where w ii = 1 for fMRI activated voxels and w ii = w 0 ∈ [0, 1] for voxels which are not revealed by fMRI analysis. A Monte Carlo simulation showed that w 0 = 0.1 (which corresponds to the 90% relative fMRI weighting) leads to a good compromise with the ability to find activation in the areas which are not found active by fMRI analysis and to detect active fMRI spots (even superficially). [3] used non-thresholded fMRI activation maps with other factors instead of one in the diagonal elements.Therefore, the

W f M RI is reformulated to (W f M RI ) ii = w 0 + (1 -w 0 ) δ i
δmax , where δ i corresponds to the relative change of the fMRI signal in the i-th voxel, and δmax is the maximal detected change. This way the relative E/MEG-fMRI scheme is preserved and locations of stronger fMRI activations have higher prior variance. [START_REF] Babiloni | Multimodal integration of high-resolution EEG and functional magnetic resonance imaging data: a simulation study[END_REF] proposed to use non-diagonal matrix for W f M RI . They proposed two different weighting matrices as follows:

w -1 ii = g(δ i ) 2 w -1 ij = g(δ i )g(δ j ).corr ij and w -1 ii = g(δ i ) 2 �G .i � -2 w -1 ij = g(δ i )g(δ j )�G .i � -1 �G .j � -1 .corr ij
where g(δ i ) is a function of the statistically significant percentage increase of the fMRI signal during the task, �G .i � is the � 2 norm of the ith column of the lead field matrix, and cor ij is the correlation of the fMRI signal of the ith and jth sources. The simulations in this study show that using the norm of the columns of the lead field matrix leads to better performance.

In [START_REF] Babiloni | Multimodal integration of EEG, MEG and fMRI data for the solution of the neuroimage puzzle[END_REF] a balancing parameter is defined to balance the effect of fMRI information. These techniques require the tuning of the balancing parameter which regulates the weight of the penalty term, in relation to the accuracy or model fit. [START_REF] Im | A technique to consider mismatches between fMRI and EEG/MEG sources for fMRI-constrained EEG/MEG source imaging: a preliminary simulation study[END_REF] proposed a method to define the balancing parameter adaptively according to the mismatches between EEG and fMRI using partial power of the EEG data in different time windows. [START_REF] Liu | Effects of fMRI-EEG mismatches in cortical current density estimation integrating fMRI and EEG: a simulation study[END_REF] proposed two level algorithm. First, they estimate the dipole magnitudes from the EEG inverse problem hardly biased with fMRI results. In the second step, they used the distance between the dipole magnitudes and the estimated dipole magnitudes in the first step as the regularization term in the EEG inverse problem. Simulations of this study show that this type of regularization leads to smaller localization error, more stability against fMRI invisible sources, fMRI extra source, and fMRI displacement in compare with the previous ones which have used the weighted � 2 norm as the regularization term.

Some authors [START_REF] Mattout | MEG source localization under multiple constraints : an extended Bayesian framework[END_REF][START_REF] Daunizeau | assessing the relevance of fMRI-based prior in the EEG inverse problem : a Bayesian model comparison approach[END_REF] have proposed principled Bayesian techniques to optimize standard regularization procedures, by estimating the balancing parameter from the EEG data. In [START_REF] Daunizeau | assessing the relevance of fMRI-based prior in the EEG inverse problem : a Bayesian model comparison approach[END_REF],

a Bayesian model comparison method is proposed to decide whether one should use the fMRI constraint or not. This approach has been applied successfully to clinical epilepsy data [START_REF] Grova | Assessing the Concordance Between Distributed EEG Source Localization and Simultaneous EEG-fMRI Studies of Epileptic Spikes[END_REF].

Symmetric Methods

Unlike data integration methods, which tend to use information from one modality to improve the other, symmetric methods, known as data fusion techniques, incorporate both modalities in a combined analysis, thus allowing for true interaction between the different data types.

Symmetric approaches can be divided into two groups; the first group is based on Datadriven analysis and second one is model-based methods. Model-based approaches investigate the goodness-of-fit of the data to the prior knowledge about the experimental paradigm and the properties of the data. In practice, very few fusion approaches have relied on realistic neurophysiological models [START_REF] Riera | Fusing EEG and fMRI based on a bottom-up model: Inferring activation and effective connectivity in neural masses[END_REF]. This is because the complexity of real metabolic-hemodynamic cascades renders the estimation of their parameters a difficult problem. In contrast, other researchers have relied on simplified variants of the neurovascular coupling model by restricting its parameters to model some common properties exhibited by "active" areas contributing to both event-related EEG and fMRI measurements [11,[START_REF] Laufs | Where does the BOLD signal goes whenalpha EEG leaves[END_REF][START_REF] Daunizeau | Symmetrical eventrelated EEG-fMRI information fusion in a variational Bayesian framework[END_REF].

However, any model-based EEG-fMRI fusion procedure will suffer from the usual limitation of modeling which is arisen from the fact that whether the assumptions of the model are satisfied

or not in a given experimental context.

Data-driven analysis methods, such as blind source separation (BSS) based on independent component analysis (ICA), are very useful in the study of a brain function. Generally, datadriven methods are very useful when the dynamics are hard to model and underlying assumptions about the data have to be minimized. Some Data-driven methods have been proposed for EEG-fMRI fusion. We will review them below.

Trilinear Partial Least Squares (tri-PLS2)

In [START_REF] Martínez-Montes | Concurrent EEG/fMRI analysis by multiway Partial Least Squares[END_REF] Partial Least-Squares (PLS) regression was generalized into the tri-PLS2 model, which

represents the E/MEG spectrum as a linear composition of trilinear components. Each component is the product of spatial (among E/MEG sensors), spectral and temporal factors, where the temporal factors have to be maximally correlated with the corresponding temporal component of the similar fMRI signal decomposition into bilinear components: products of the spatial and temporal factors. Analysis the data from [START_REF] Goldman | Simultaneous EEG and fMRI of the alpha rhythm[END_REF] using tri-PLS2 found a decomposition into 3 components corresponding to alpha, theta and gamma bands of the EEG signal. The fMRI components had a strong correlation only in alpha band component (Pearson correlation 0.83 (p = 0.005)), while the theta component also showed a linear correlation of 0.56 (p = 0.070). It is interesting to note, that spectral profiles of the trilinear EEG atoms received with and without fMRI influence were almost identical, which can be explained either by the non-influential role of fMRI in tri-PLS2 decomposition of EEG, or just by a good agreement between the two. On the other hand, EEG definitely guided fMRI decomposition, so that the alpha rhythm spatial fMRI component agreed very well with the previous findings [START_REF] Goldman | Simultaneous EEG and fMRI of the alpha rhythm[END_REF].

Joint ICA

Joint independent component analysis (jICA) has been successfully used for the fusion of data from the two modalities such as fMRI and EEG [START_REF] Calhoun | Neuronal chronometry of target detection: Fusion of hemodynamic and event-related potential data[END_REF].

Given two feature data sets X 1 and X 2 , the jICA approach involves concatenating the data sets alongside each other and then performing ICA on the concatenated data set as in [X

1 X 2 ] = A[S 1 S 2 ]
. Joint-ICA assumes that the sources have a common modulation profile (A) across subjects or modalities, which is a strong constraint considering that the data come from two different modalities.

Parallel ICA

Parallel ICA has been proposed in [START_REF] Liu | Parallel independent component analysis for multimodal analysis: an application to fMRI and EEG data[END_REF][START_REF] Liu | Mapping the bilateral visual integration by EEG and fMRI[END_REF]. It performs separate analysis on each modality and enhances the inter-subject co-variations by constraining the correlations between modulation profiles. It is flexible as it does not assume identical inter-subject variations.

In [START_REF] Wu | Parallel independent component analysis using an optimized neurovascular coupling for concurrent EEG-fMRI sources[END_REF], a multivariate parallel ICA decomposition incorporating dynamic neurovascular coupling for concurrent EEG-fMRI recordings is proposed. This method simultaneously optimizes both the maximum independence of EEG and fMRI sources in addition to the maximum neurovascular coupling which is derived from the correlation between the HRF-convolved EEG and the fMRI time courses. The main difference between mCCA and ICA-based methods (jICA and Parallel ICA) is that mCCA searches for uncorrelated sources while the other methods search for independent sources.

Relaxing the independence criterion may enable multimodal-CCA to explore the connectivity between different networks of functional areas associated with both hemodynamic response and electric potential.

mCCA and jICA are compared in [START_REF] Correa | Canonical correlation analysis for feature-based fusion of biomedical imaging modalities and its application to detection of associative networks in schizophrenia[END_REF] and the results show that they have similar performance when two modalities ahave strong connections and mCCA outperfoms when the connections are weak. Connections between the two modalities are simulated by generating correlation between profile pairs, formed across modalities. See figure 1.7.

STEFF

[90] proposed a hybrid approach, in which group ICA is used for generating temporal and spatial priors from the EEG data recorded simultaneously with the fMRI data, respectively. Then, the temporal and spatial priors are introduced to a simple model-driven approach. The study compares MNE (Minimum Norm Estimation), LORETA, MSP (Multiple Sparse Priors [START_REF] Friston | Multiple sparse priors for the m/eeg inverse problem[END_REF]), fMRI-weighted MNE [START_REF] Liu | Spatiotemporal imaging of human brain activity using functional MRI constrained magnetoencephalography data: Monte Carlo simulations[END_REF] and STEFF (Spatial-Temporal EEG/fMRI Fusion). The comparison shows that the method with the best spatial resolution is STEFF and the others are sorted as MSP, fMRI-weighted MNE, LORETA and MNE after STEFF.

Applications

Epilepsy

Simultaneous EEG-fMRI recordings are used greatly in the evaluation of epilepsy, namely in the characterization of brain networks related to epileptic activity. In EEG-fMRI studies, epileptic events are usually described as boxcar signals based on the timing information retrieved from the EEG, and subsequently convolved with a hemodynamic response function to model the associated Blood Oxygen Level Dependent (BOLD) changes [START_REF] Grova | Assessing the Concordance Between Distributed EEG Source Localization and Simultaneous EEG-fMRI Studies of Epileptic Spikes[END_REF][START_REF] Carney | Siblings with refractory occipital epilepsy showing localized network activity on EEG-fMRI[END_REF][START_REF] Walker | EEG-fMRI in Adults with Focal Epilepsy[END_REF].

Cognitive Study

Different integration methods of EEG-fMRI data are used in cognitive tasks. Some examples of these studies are the processing of target detection (oddball paradigm) [START_REF] Eichele | Assessing the spatiotemporal evolu-tion of neuronal activation with singletrial event-related potentials and functional MRI[END_REF][START_REF] Goldman | Single-trial discrimination for integrating simultaneous EEG and fMRI: identifying cortical areas contributing to trial-to-trial variability in the auditory oddbal ltask[END_REF], performance monitoring [START_REF] Debener | Trial-by-trial coupling of concurrent electroencephalogram and functional magnetic resonance imaging identifies the dynamics of performance monitoring[END_REF], decision making [START_REF] Mulert | Single-trial coupling of EEG and fMRI reveals the involvement of early anterior cingulate cortex activation in effortful decision making[END_REF], working memory maintenance [START_REF] Scheeringa | Trial-by-trial coupling between EEG and BOLD identifies networks related to alpha and theta EEG power increases during working memory maintenance[END_REF][START_REF] Hanslmayr | The Relationship between Brain Oscillations and BOLD Signal during Memory Formation: A Combined EEG-fMRI Study[END_REF] and simple auditory responses [START_REF] Mayhew | EEG signatures of auditory activity correlate with simultaneously recorded fMRI responses in humans[END_REF].

Resting-State Brain Activity

The study of [START_REF] Mantini | Electrophysiological signatures of resting state networks in the human brain[END_REF] used co-registered EEG-fMRI to detect the electrophyiological signatures of resting state networks. In [START_REF] Gonçalves | A Data and Model-Driven Approach to Explore Inter-Subject Variability of Resting-State Brain Activity Using EEG-fMRI[END_REF] the generalized linear model is used to determine which of the identified patterns correlates significantly to the spontaneous variations of the alpha rhythm. [START_REF] Biazoli | Application of Partial Directed Coherence to the analysis of resting-state EEG-fMRI data[END_REF] extracts partial directed coherence (PDC) measure from the EEG data which is used to find the nodes of resting-state networks in fMRI.

Chapter 2

Motivation

We are interested in studying the brain function to find out how the human brain works. It would help the physicians to find a treatment for mental diseases and it would also work as an inspiration for inventors to design intelligent systems and robots.

Several imaging modalities have been created to study human brain functionality. To use the advantage of different modalities and recover the disadvantage of one modality, multimodality analysis has been a focus of several researches in the recent years. One possible and more general modality fusion is EEG-fMRI integration. Integration of EEG and fMRI is of particular interest because of two main reasons. First, the feasibility of their simultaneous acquisition is helpful to have synchronized data. Second, because of their complementary characteristic, high spatiotemporal source localization becomes achievable. Overall, the purpose of the EEG-fMRI integration is to achieve high spatial resolution of fMRI and high temporal resolution of EEG simultaneously.

EEG-fMRI integration approaches can be divided in two groups. The first group uses one modality as an auxiliary data for analyzing the other modality data [START_REF] Debener | Trial-by-trial coupling of concurrent electroencephalogram and functional magnetic resonance imaging identifies the dynamics of performance monitoring[END_REF][START_REF] Im | A technique to consider mismatches between fMRI and EEG/MEG sources for fMRI-constrained EEG/MEG source imaging: a preliminary simulation study[END_REF][START_REF] Ou | Multimodal functional imaging using fMRI-informed regional EEG/MEG source estimation[END_REF]. The second group, which is called symmetric or fusion, uses both EEG and fMRI data jointly to find spatiotemporal sources [START_REF] Moosmann | Joint independent component analysis for simultaneous EEG-fMRI: principle and simulation[END_REF][START_REF] Calhoun | Parallel independent component analysis for multimodal analysis: an application to fMRI and EEG data[END_REF][START_REF] Liu | Parallel independent component analysis for multimodal analysis: an application to fMRI and EEG data[END_REF][START_REF] Wu | Parallel independent component analysis using an optimized neurovascular coupling for concurrent EEG-fMRI sources[END_REF][START_REF] Daunizeau | Symmetrical eventrelated EEG-fMRI information fusion in a variational Bayesian framework[END_REF]. Each group also can be divided into different groups. In this work, the spatial map of the fMRI data is used as prior information to solve the EEG inverse problem such that the most common and uncommon sources can be detected. Therefore, the proposed method settles down in the fMRI-informed EEG analysis group. The early methods in this group add fMRI results in the weighting matrix in the regularization terms in the EEG inverse problem [START_REF] Babiloni | Integration of high resolution EEG and functional magnetic resonance in the study of human movement-related potentials[END_REF][START_REF] Babiloni | Cortical source estimate of combined high resolution EEG and fMRI data related to voluntary movements[END_REF]. These methods interpreted fMRI spatial map as the power of EEG data and used weighted � 2 norm regularization approach. Therefore, their result is very smooth and sensitive to the mismatch of the EEG-fMRI data [START_REF] Liu | Effects of fMRI-EEG mismatches in cortical current density estimation integrating fMRI and EEG: a simulation study[END_REF]. Here, we use the fMRI results as a spatial hint for the EEG inverse problem, which does not need to find the real physiological relation between EEG and fMRI data. Although we also add fMRI data in the weighting matrix in the regularization term, the regularization term is constructed from the spatial map and it is not based on the dipole currents directly. Therefore the interpretation of the way we use fMRI prior information is directly the spatial prior information.

The aim of our method is applying fMRI result or any other prior information in the inverse problem of EEG to improve the spatial resolution of the source localization. This method would be useful when we do not have the fMRI data for each subject but only a group fMRI analysis or a reliable fMRI analysis result which can be used as prior information is available. The proposed method in this survey is fast with low computations in comparison to previous ones and it is more stable to invisible sources in each modality.

First, we propose a reference based source separation method which has been evaluated on the intracranial EEG signals of 5 epileptic patients who were seizure free after an operation. We present a fast method to extract the sources related to the interictal epileptiform state. The method is based on general eigenvalue decomposition using two correlation matrices during: 1) periods, including a reference activation which can be interictal epileptiform discharges (IED), and 2) background activity periods which can be periods excluding IEDs or abnormal physiological signals. After extracting the most similar sources to the reference or IED state, the related sensors are estimated by using a multi-objective optimization method. The method is evaluated using both realistic simulated data and actual intracerebral electroencephalography (iEEG) recordings of patients suffering from focal epilepsy. These patients are seizure-free after the resective surgery. Quantitative comparisons of the proposed IED regions with the visually inspected ictal onset zones by the epileptologist and another method of identification of IED regions reveal good performance.

Secondly, we introduce an EEG-fMRI integration method. R-SS method is used to extract the sources from the scalp EEG signal. The method is based on the spatial sparsity of neural activity using weighted sparse decomposition. An optimization method is met which is also evaluated on a simple synthetic data. Here, we analyze fMRI data using generalized linear model (GLM) and the results are used in solving EEG inverse problem as a weighting matrix in the weighted sparse decomposition. We also use the � 0 norm as a measure of sparsity instead of the � 1 norm. The performance of the proposed method is evaluated quantitatively along with multiple sparse priors (MSP) method with no prior information, and MSP uses fMRI results as prior information. Source distribution index (SDI) and localization bias (LB) are used to measure the performance of different localization approaches with or without a variety of fMRI-EEG mismatches. The method is also applied to the real data of face perception of 16 subjects and its stability against subject variations is evaluated with the number of false positive and true positive regions. Simulation results show that the proposed method has the lowest localization bias and it is more stable to noise. Although the existence of an extra region in the fMRI data enlarges the localization bias, the proposed method outperforms the other methods. Conversely, a missed region in the fMRI data does not affect the localization bias of the common sources in the EEG-fMRI data. Results on the real data are congruent with the previous studies and produce two clusters in the fusiform and occipital face areas (FFA and OFA, respectively).

Moreover, it shows more stable localization against variations in different subjects.

The main contributions of the current work are: � Reference-based source separation method estimates a set of relevant sources of the state of interest.

� � 2,0 regularization finds the sparse correlated sources.

� EEG-fMRI integration method which estimates the sparse decomposition of the EEG signal in its relevant source space.

� Using the spatial map estimated from the fMRI data as a weighting matrix in the sparse decomposition.

� Introduce Pareto optimization for the final decision of the brain active regions.

The rest of the thesis is organized as follows:

1. Part II: an algorithm for � 2,0 regularization is proposed in this part, which would be used in Chapter 6.

2. Part III: Reference-based source separation method basics and its evaluation on the simulation and actual data are explained in this part including two chapters: (a) Chapter 4: Reference-based source separation method is introduced in this chapter.

(b) Chapter 5: in this chapter the R-SS method is applied to the realistic simulation data and also it is applied to the intracranial EEG data of 5 epileptic patients and it is also evaluated quantitatively.

Introduction

The synchronized neuron activations which produce currents within the brain can be modeled with some dipoles. The effect of these dipoles on voltages measured at the scalp electrodes using Maxwell's equations is estimated by the linear forward model.

X = GJ + E n (3.1)
where the lead field G (N ×M ) would be calculated, using structural MRI and boundary element method (BEM) X is a N × T matrix that contains scalp readings, N is the number of channels and T is the number of time points, J M ×T contains dipole magnitudes at different time instants, The inverse problem consists in finding an estimate of the dipole temporal courses, which are the rows of matrix J, given the electrode positions, scalp readings X and the gain matrix G. The EEG inverse problem is under-determined (since N � M ) and lacks a unique solution, which is due to the fact that there are fewer observations than the unknown variables.

Two main approaches to the inverse problem are non-parametric (imaging) and parametric (scanning) methods [START_REF] Baillet | Electromagnetic Brain Mapping[END_REF]. A brief review of these methods can be found in [START_REF] Bolstad | Space-time event sparse penalization for magneto-electroencephalography[END_REF]. Parametric or scanning methods use a small number of dipoles [START_REF] Michel | EEG source imaging[END_REF][START_REF] Bai | On the estimation of the number of dipole sources in EEG source localization[END_REF][START_REF] Shun | Efficient Dipole Parameter Estimation in EEG Systems With Near-ML Performance[END_REF], multi-poles [START_REF] Mosher | EEG and MEG: forward solutions for inverse methods[END_REF][START_REF] Bandyopadhyay | A numerical approach for automatic detection of the multipoles responsible for ill conditioning in generalized multipole technique[END_REF], or cortical patches [START_REF] Limpiti | Cortical patch basis model for spatially extended neural activity[END_REF][START_REF] Wagner | Human sleep EEG under the influence of pulsed radio frequency electromagnetic fields. Results from polysomnographies using submaximal high power flux densities[END_REF] and scan over locations to find the best set of sources to represent the data. Examples of scanning methods include MUSIC [START_REF] Mosher | EEG and MEG source localization using recursively applied (RAP) MUSIC[END_REF], beamforming [START_REF] Van Veen | Minimum variance beamforming with soft response constraints[END_REF][START_REF] Fuchs | Anatomically constrained minimum variance beamforming applied to EEG[END_REF], and maximum likelihood estimation [START_REF] �nan | AR spectral analysis of EEG signals by using maximum likelihood estimation[END_REF]. Non-parametric methods are also referred to as Imaging methods. In these models several dipoles with fixed locations and possibly fixed orientations are distributed in the whole brain volume or cortical surface. Here, we focus on imaging methods as our proposed method is an imaging method.

Obtaining a physiologically feasible solution can be done through reasonable constraints. The brain neural activities are sparsely localized [START_REF] Simoncelli | Natural images statistics and neural reperesentation[END_REF][START_REF] Op De Beeck | Interpreting fMRI data: Maps, modules, and dimensions[END_REF], thus a relevant constraint is to use the sparse characteristics of the spatial distribution of the sources [START_REF] Gramfort | Mixed-norm estimates for the M/EEG inverse problem using accelerated gradient methods[END_REF][START_REF] Ding | Sparse cortical current density imaging in motor potentials induced by finger movement[END_REF][START_REF] Ding | Reconstructing Cortical Current Density by Exploring Sparseness in the Transform Domain[END_REF][START_REF] Friston | Multiple sparse priors for the m/eeg inverse problem[END_REF]. Mostly used penalty terms are � 2 and � 1 norm of J, which are respectively referred as Tikhonov regularization [START_REF] Tikhonov | Solutions of ill-posed problems[END_REF] and LASSO [START_REF] Tibshirani | Regression shrinkage and selection via the lasso[END_REF]. These penalty terms and different methods are reviewed in [START_REF] Bolstad | Space-time event sparse penalization for magneto-electroencephalography[END_REF].

The � 2 norm is the most common regularization method, which is known as Minimum Norm Estimation (MNE). For the problem:

X = GJ + E n E n ∼ N (0, C n )
The optimization problem leads to:

Ĵ = arg min J � 1 2 �E -GJ� 2 2 + λ�WJ� 2 � . (3.2)
Minimum Norm Estimation has an analytic solution which leads to a simple matrix multiplication [START_REF] Tikhonov | Solutions of ill-posed problems[END_REF]. This makes the estimation extremely fast.

J = (W � W) -1 G � [G(W � W) -1 G � + λC n ] -1 X (3.3)
This corresponds to a penalized maximum likelihood estimate assuming the sources are Gaussian and normally distributed, with a diagonal covariance matrix [START_REF] Wipf | A unified Bayesian framework for MEG/EEG source imaging[END_REF]. Now, consider the following 2-level hierarchical linear model:

X = GJ + E n E n ∼ N (0, C n ) J = 0 + E j E j ∼ N (0, C j )
where the notation N specifies a Gaussian distribution. Then the likelihood p(Y|J) = N (GJ, C n ) with prior as p(J) = N (0, C j ) leads to maximum a posteriori (MAP) estimate as follows:

p(J|X) ∝ p(X|J)p(J)

J M AP = C j G � [GC j G � ] -1 X (3.4)
Comparing Eq. 3.4 with Eq. 3.3 shows that they are equivalent if

C j = (W � W) -1 .
The hypothesis in the � 2 norm regularization method is the smoothness of the solution.

Therefore, the estimated sources with this method have too smooth spatial maps and the extent of active regions is often over-estimated. Therefore, this method fails to recover high spatial frequencies. To address this limitation many alternatives to MNE have been proposed using different weighting matrices.

Minimum current estimate(MCE) uses the � 1 regularization, which has been proposed to find a focal spatial map and avoid smoothness [START_REF] Uutela | Visualization of Magnetoencephalographic Data Using Minimum Current Estimates[END_REF]. This corresponds to a penalized maximum likelihood estimate assuming the sources are exponentially distributed (E j ∼ exp(-|E j |)). One reason is that, under certain conditions, it has been proved that sparsity could enable the perfect resolution of ill-posed problems [START_REF] Donoho | Compressed sensing[END_REF][START_REF] Candès | The dantzig selector: statistical estimation when p is much larger than n[END_REF]. The � 1 norm is a strong hypothesis which promotes sparse solutions, which means that the solution should only have a small number of non-zero coefficients. In the � 1 regularization context, many algorithms have been proposed. As an example, FOCUSS algorithm has been proposed which uses Iterative Re-weighted Least Squares (IRLS) [START_REF] Gorodnitsky | Neuromagnetic source imaging with focuss: a recursive weighted minimum norm algorithm[END_REF]. The drawback of the � 1 regularization compared with the � 2 regularization, is that the optimization problem does not have an analytic solution. In the last few years, the machine learning and signal processing communities have devoted a lot of efforts into the improvement of the optimization methods that help to solve non-differentiable problems arising when considering sparse priors. Among the list of algorithms that have been proposed, IRLS methods, similar to the FOCUSS algorithm consist in iteratively computing weighted MN solutions with weights updated after each iteration [START_REF] Daubechies | Iteratively re-weighted least squares minimization: Proof of faster than linear rate for sparse recovery[END_REF]. The LARS-LASSO algorithm [START_REF] Tibshirani | Regression shrinkage and selection via the lasso[END_REF][START_REF] Efron | Least Angle Regression[END_REF][START_REF] Tibshirani | Sparsity and Smoothness via the Fused Lasso[END_REF] is an extremely powerful method for solving the � 1 problem. Simple coordinate descent methods [START_REF] Friedman | Regularization Paths for Generalized Linear Models via Coordinate Descent[END_REF] or blockwise coordinate descent, also called Block Coordinate Relaxation (BCR) [START_REF] Bruce | Block coordinate relaxation methods for nonparamatric signal denoising[END_REF], are also possible strategies. Alternatively, methods based on projected gradients have been proposed in [START_REF] Daubechies | Iteratively re-weighted least squares minimization: Proof of faster than linear rate for sparse recovery[END_REF][START_REF] Beck | A fast iterative shrinkage-thresholding algorithm for linear inverse problems[END_REF].

Most of these methods belong to the general Bayesian framework proposed in [START_REF] Wipf | A unified Bayesian framework for MEG/EEG source imaging[END_REF], which is explained in 3.2.

Bayesian Modeling with General Gaussian Scale Mixtures and Arbitrary Covariance Components

In this section we present a general purpose Bayesian framework for source localization. We begin with the noise model:

p(X|J) ∝ exp(- 1 2 �X -GJ� 2 Σ -1 � ) (3.5)
where �X� Σ -1

�

denotes the weighted matrix norm

� trace[X � Σ -1
� X] and Σ � is the noise covariance. Next the following source prior is assumed:

p(J|γ) ∝ exp(- 1 2 trace � J � Σ -1 s J � ), Σ s = dγ � i=1 γ i C i (3.6) where γ = [γ 1 , • • • , γ dγ ]
� is a vector of d γ nonnegative hyperparameters that control the relative contribution of each covariance basis matrix C i . While the hyperparameters are unknown, the set of components

C = {C i : i = 1, • • • , C dγ } (3.7)
is assumed to be fixed and known. In the simplest case, Σ s = γ 1 C 1 = γ 1 I , where I is an identity matrix, leads to a weighted minimum � 2 norm solution. If we assume the underlying source currents are formed from a collection of dipolar point sources located at each vertex, then we

may choose C = {e i e � i : i = 1, • • • , C dγ }
, where e i is a standard index vector of zeros with a '1' for the i-th element. Focal Underdetermined System Solution (FOCCUSS) [START_REF] Cotter | Sparse solutions to linear inverse problems with multiple measurement vectors[END_REF][START_REF] Gorodnitsky | Neuromagnetic source imaging with focuss: a recursive weighted minimum norm algorithm[END_REF], MCE, sparse Bayesian learning (SBL) [START_REF] Tipping | Sparse bayesian learning and the relevance vector machine[END_REF] and Multiple Sparse Priors (MSP) [START_REF] Friston | Multiple sparse priors for the m/eeg inverse problem[END_REF] have used the same assumptions. [START_REF] Wipf | A unified Bayesian framework for MEG/EEG source imaging[END_REF] analyzes and extends several broad categories of Bayesian inference directly applicable to source localization.

Mixed Norm Regularization

The � 1 norm regularization leads to spread sparse sources. However, the activation sites of the brain activities are localized and each localization contains correlated dipoles, but correlated sources are not estimated by the � 1 norm regularization. In order to go beyond these limitations, there has been a growing interest for the mixed norm constraints, so that the advantages of these norms are integrated into the solution [START_REF] Valdés-Sosa | EEG source imaging with spatio-temporal tomographic nonnegative independent component analysis[END_REF][START_REF] Vega-Hernández | Penalized least squares methods for solving the EEG inverse problem[END_REF]. The � 2 norm regularization has the advantage of recovering correlated sources, but its spatial resolution is low as the result is very smooth. The � 1 norm regularization has the advantage of providing sparse results, but cannot find correlated sources. Integrating both of them in the penalty term could lead to sparse correlated sources.

The mixed norm constraints can be divided into two groups which are reviewed in sections 3.3.1 and 3.3.2, respectively. The first group, which is called Mixed-Norm Estimates (MxNE), has been proposed in [START_REF] Gramfort | Mixed-norm estimates for the M/EEG inverse problem using accelerated gradient methods[END_REF][START_REF] Gramfort | Time-frequency mixednorm estimates: Sparse M/EEG imaging with non-stationary source activation[END_REF]. The second group, which is called Weighted Elastic Net (WEN), has been proposed in [START_REF] Zou | Regularization and variable selection via the elastic net[END_REF].

We propose a new mixed norm constraint which is more accurate and more stable against the noise and it is simpler than previous ones. The method is introduced in section 3.5.

Mixed-Norm Estimates (MxNE)

In this method, the � 1 norm of the � 2 norm of the time courses of the current dipoles are used as a penalty term as follows:

f (J) = � i � � k j ik 2 (3.8)
where j ik is the entry of the matrix J located at the i-th row and k-th column. � k j ik 2 represents the power of the i-th dipole temporal course. The idea behind this method is to estimate J with a minimum power of temporal courses in the sense of the � 1 norm. The � 1 norm leads to the sparsity of the dipoles power which is actually true, because the power of most of the current dipoles are too small because few of them are located in the active regions. In addition, in the time domain, the � 2 norm will apply the smoothness constraint. In this method, the time series of the sources and their spatial map will be estimated simultaneously. The related optimization problem can be solved by SOCP (Second Order Cone Program) [START_REF] Gramfort | Mixed-norm estimates for the M/EEG inverse problem using accelerated gradient methods[END_REF] but with a high computational cost.

Weighted Elastic Net

A newly developed variable selection method [START_REF] Zou | Regularization and variable selection via the elastic net[END_REF], called elastic net (EN), can produce a sparse model with good prediction accuracy, while encouraging a grouping effect. Empirical results and simulations have demonstrated superiority of the elastic net over LASSO [START_REF] Zou | Regularization and variable selection via the elastic net[END_REF].

Weighted elastic net (WEN) has also been developed in [START_REF] Hong | Weighted Elastic Net Model for Mass Spectrometry Imaging Processing[END_REF][START_REF] Jun-Tao | An Improved Elastic Net for Cancer Classi ¯cation and Gene Selection[END_REF]. Weighted elastic net for A fast algorithm has been proposed in [START_REF] Babaie-Zadeh | Weighted Sparse Signal Decomposition[END_REF], which is called weighted smooth � 0 (WSL0), to solve (3.10). The main idea of smooth � 0 (SL0) is to use a smooth approximation of the � 0 norm. Here, we briefly review SL0 and WSL0 as we need to know their ideas in the following.

Assume that f σ () is a continuous function satisfying lim σ→0 f σ (b) = 1 -|b| 0 , for example f σ (b) = exp ( -b 2 2σ 2 )
. Then |b| 0 can be approximated as

|b| 0 ≈ 1 -f σ (b)
, where σ determines the accuracy of the approximation: the smaller σ, the better approximation, and the larger σ, the smoother approximation. So, the optimization problem is as follows:

maximize b F σ = � f σ (b i ) = � exp (-b 2 i /2σ 2 ) subject to a = Gb (3.11)
[142] used graduated non-convexity (GNC) approach for maximizing the above optimization problem. GNC starts with a very large σ, for which F σ is nearly concave and its maximization is easy, and then gradually decreases σ, and for each σ it starts the search for the maximizer of F σ from the maximizer for the previous (larger) σ. Using such an annealing process, it is expected (but not mathematically guaranteed) to escape from being trapped into the local maximum. [START_REF] Babaie-Zadeh | Weighted Sparse Signal Decomposition[END_REF] showed that GNC-based SL0 solver is directly applicable to the weighted minimization. In this

case, |s| 0,w ≈ � m i=1 w i - � m i=1 w i f σ (b i )
, and so, the following term should be maximized:

maximize b F w σ = � w i f σ (b i ) = � w i exp (-b 2 i /2σ 2 )
subject to a = Gb (3.12)

The final algorithm would be as shown in Figure 3.1. There are a few more points to be mentioned about this algorithm: For smaller σ's, F w σ is more fluctuating and we should use a smaller step-size of the gradient ascent loop for its maximization. For reasons detailed in [START_REF] Eftekhari | Robust-{SL0} for stable sparse representation in noisy settings[END_REF], a good choice is to decrease the step-size proportional to σ 2 as µσ 2 . Also a good initialization is to use the maximizer of F w σ (b) for σ → ∞ which is the minimum weighted � 2 norm [START_REF] Eftekhari | Robust-{SL0} for stable sparse representation in noisy settings[END_REF].

WSL0:

1. Initialization:

� Set b 0 equal the minimum weighted � 2 norm solution of a = Gb.
� Choose a suitable decreasing sequence for σ :

[σ 1 • • • σ J ]. 2. For j = 1, • • • , J: (a) Let σ = σ j .
(b) Solve (3.12) using L Gradient-Projection iterations:

� Initialization: b = bj-1 . � For � = 1, 2, ..., L -Let s ← s + (µσ 2 )∇F w σ (b). -Project s back onto the feasible set b|Gb = a: s ← b -G � (Gb -a). � Set bj = b.
3. Final answer is b = bJ .

� 2,0 Regularization

In this section, we propose a new method of optimization with sparsity constraint, where we replace a l1-norm term by a l0-norm, inspired by [START_REF] Eftekhari | Robust-{SL0} for stable sparse representation in noisy settings[END_REF]. A fast stable method has been proposed for the � 0 norm regularization in [START_REF] Eftekhari | Robust-{SL0} for stable sparse representation in noisy settings[END_REF] and since then it can be used in many sparse decomposition problems. The � 0 norm regularization leads to a very sparse response. The idea of combining the � 2 norm and the � 0 norm is raised for using the advantages of both to have correlated sparse regions in the response. The idea is the same as elastic net but, we replace the � 1 norm with the � 0 norm in the weighted elastic net (3.9).

b = arg min

b �a -Gb� + λ 2 �W 2 b� + λ 1 �W 1 b� 0 (3.13) or minimize b �W 1 b� 0 subject to �a -Gb� + λ 2 �W 2 b� < δ
It can be simplified using the following change of variables [START_REF] Hong | Weighted Elastic Net Model for Mass Spectrometry Imaging Processing[END_REF]:

G * = (1 + λ 2 ) -1 2 � G √ λ 2 W 2 � a * = � a 0 � b * = � 1 + λ 2 b
then, the optimization problem can be written as:

minimize b �W 1 b * � 0 subject to a * = G * b * (3.14)
Now, WSL0 can be used for this problem directly. The variables in (3.14) are greater in size than the initial variables as they are concatenated versions of the initial variables. For example

G is a N × M matrix while G * is a (N + M ) × M matrix. In underetermined problems N � M
and this would result in high computation load. To recover the high computation load we have used a multi-objective optimization approach.

� 2,0 :

1. Initialization:

� Set b 0 equal the minimum weighted � 2 norm solution of a = Gb. � Choose a suitable decreasing sequence for σ : [σ 1 • • • σ J ]. 2. Forj = 1, • • • , J: (a) Let σ = σ j .
(b) Solve ( 8) using L Gradient-Projection iterations:

� Initialization: b = bj-1 . � For � = 1, 2, ..., L -Let s ← s + (µσ 2 )∇F w σ (b). -Project s back onto the feasible set {b|Gb = a}: s ← b-G � (Gb-a). � For � = 1, 2, ..., K -Let s ← s -µ k s. -Project s back onto the feasible set {b|Gb = a}: s ← b -G � (Gb - a). -if F w σ (b) > (1 -δ 100 )(F w σ (b �-1 )) i. b = b �-1 ii. µ k = ρµ k � Set bj = b.
3. Final answer is b = bJ . Multi-objective optimization consists of optimizing a set of objectives. A p-objective optimization problem can be formally described as:

minimize f(x) = f 1 (x), f 2 (x), ..., f p (x) subject to g j (x) ≤ 0 , i = 1, • • • , m (3.15)
where x = [x 1 , x 2 , ..., x n ] � is the vector of n decision variables and g i (x) represent a set of m inequality constraints.

One way to tackle multiple objectives is lexicographic method, a technique that requires the decision-maker to establish a priority for each objective [START_REF] Collette | [END_REF]. With the lexicographic method, the objective functions are arranged in order of importance. Then the following optimization problems are solved one at a time:

minimize x∈X f i (x) , i > 1, i = 1, • • • , p subject to f j (x) ≤ f j (x * j ) j = 1, • • • , i -1
Here, i represents a function's position in the preferred sequence, and f j (x * j ) represents the optimum of the j-th objective function, found in the j-th iteration. The lexicographic approach is usually useful when dealing few objectives (two or three). It should also be reminded that sometimes its performance is tightly subject to the ordering of objectives.

A constraint relaxation is represented by [START_REF] Osyczka | Multicriterion Optimization in Engineering with FORTRAN Programs[END_REF]. It is induced by increasing the right hand side of the constraint by a percentage of f j (x * j ) as follows:

f j (x) ≤ (1 + δ j 100 )f j (x * j ) (3.16)
where δ j ranges between 0 and 100.

The � 2,0 optimization problem can also be considered as a multi-objective optimization problem as follow:

minimize f(b) = [�W 1 b� 0 , �W 1 b�] subject to a = Gb
where, the � 0 norm is the first objective and the � 2 norm is the second one. The WSL0 method minimizes � 0 iteratively by decreasing σ value. For each σ value, after minimizing the � 0 norm, we try to minimize the � 2 norm as shown in Figure 3.2. The δ value is equal to 5.

Data and Results

We conducted a simulation to compare weighted WSL0 and � 2,0 . In this simulation, we randomly

G = [G 1 , G 2 , a -G 2 b � 2 ], b 1 = (b � 1 , 0 1×50 ) � , and b 2 = (0 1×50 , b � 2 , 1) � .
Finally, the columns of G were normalized to have unit � 2 -norm, and the entries of b 1 and b 2 were scaled accordingly. The output ( b) was then compared to b 1 and b 2 using their correlation coefficient (ρ 1 and ρ 2 , respectively). For ρ 1 > 0.9, it is expected that b 1 is estimated by the algorithm, while for ρ 2 > 0.9 it is expected that the algorithm result would be close to b 2 ; otherwise, it was declared that the algorithm has failed to estimate either one of the solutions. This experiment was repeated 1,000 times with different randomly generated systems.

Conclusion

In this chapter, we have introduced previous regularization methods which have been used in EEG source localization. We have also proposed a new regularization term which is the combination of the � 0 norm and the � 2 norm. It is called � 2,0 regularization. Theoretically, It is more accurate than the similar methods proposed before, because of using the � 0 norm. In a simple simulation, we have shown that the � 2,0 regularization tends to estimate the sparse connected regions and it is more accurate than weighted elastic net. As the optimization algorithm is also simple, it seems that the � 2,0 regularization would be a good choice for our problem and it would be used in chapter IV. Some experiments can be considered as future works listed as follows:

� The � 2,0 regularization algorithm should be verified on different simulations which simulate different situations in the real data. For example, considering different types of noise.

� Choosing the best value for the repetition number of the � 2 norm minimization (K) in the

� 2,0 regularization (Figure 3.2). Part III Reference Based Source Separation (R-SS) Chapter 4 R-SS Method 4.

Introduction

In this chapter, reference-based source separation (R-SS) is introduced. Its solution is computed by generalized eigenvalue decomposition (GEVD) which is fast since it admits an exact analytic solution. A similar idea has already been applied for extracting fetal ECG from maternal ECG [START_REF] Sameni | Multichannel Electrocardiogram Decomposition Using Periodic Component Analysis[END_REF][START_REF] Lathauwer | SVD-based methodologies for fetal electrocardiogram extraction[END_REF]. Spatial filters computed by GEVD to maximize the power of the temporal sources during the reference state using GEVD, are used for enhancing spatial resolution of recordings. As a result, GEVD provides temporal sources sorted according to their decreasing similarity with the reference. Selecting the right number of temporal sources is a tricky problem.

In this chapter, as a second processing step, we propose an original method, based on Bayes rule, is then applied for selecting the number of sources similar enough to the reference.

GEVD principles are reviewed in Section 4.2. Data structure and the heart of the R-SS method are presented in Sections 4.3 and 4.4, respectively. Source classification method which is the last step in R-SS is explained in Section 4.7.

GEVD Principles

Generalized eigenvalue decomposition (GEVD) of one pair of symmetric and positive definite matrices (R 1 ,R 2 ) can be stated as follows:

R 1 W = R 2 WΛ (4.1)
where Λ is the diagonal matrix of generalized eigenvalues. The generalized eigenvectors build the columns w of matrix W. The ratio

λ(w) = w � R 1 w w � R 2 w (4.2)
where (•) � denotes transpose, is known as the Rayleigh quotient, and its maximization is usually performed by GEVD. In fact, the first derivative of λ whose roots are the extreme points of λ is:

∇λ(w) = 2 w � R 2 w (R 1 w -λR 2 w). (4.3)
Eq. ( 4.3) shows that the greatest eigenvalue in (B.1) will be the maximum value of Rayleigh quotient.

Defining various pairs of matrices (R 1 ,R 2 ) and using GEVD for maximizing the Rayleigh quotient can be used in different applications, associated with different criteria [START_REF] Sameni | Multichannel Electrocardiogram Decomposition Using Periodic Component Analysis[END_REF]. GEVD has some useful properties which are defined in the following two theorems.

GEVD Properties

Theorem 4.1. For symmetric positive definite matrices (r 1 and r 2 ), we have:

w � i R 2 w j = � 0 i � = j β i > 0 i = j (4.4)
and

w � i R 1 w j = � 0 i � = j λ i β i > 0 i = j (4.5)
Proof : For i-th eigenvector we have

R 1 w i = λ i R 2 w i (4.6)
The scalar product with another eigenvector gives:

w � j R 1 w i = λ i w � j R 2 w i (4.7)
and of course also:

w � i R 1 w j = λ j w � i R 2 w j (4.8)
Since R 1 and R 2 are symmetric:

λ i w � i R 2 w j = λ j w � i R 2 w j (4.9)
and hence:

(λ i -λ j )w � i R 2 w j = 0 (4.10)
If λ iλ j � = 0 then:

w � i R 2 w j = 0 (4.11)
and if i = j, as R 2 is positive definite we have:

w � i R 2 w j = β i > 0 (4.12)
In the same way we have:

w � i R 1 w j = � 0 i � = j λ i β i > 0 i = j (4.13)
Theorem 4.2. The eigenvectors are linearly independent.

Proof : Suppose w i are not linearly independent. This would mean that we could write an eigenvector w k as

w k = � j� =ii γ j w j (4.14)
This means that for j � = k,

λ i w � j R 2 w k = λ j w � j R 2 w j � = 0 (4.15)
which violates equation 4.4, the eigenvectors are linearly independent.

Data

Let us denote the observation data as X ∈ R N ×T , and its rows as

x � i = [x i [1], • • • , x i [T ]], i = 1, • • • , N is a zero-mean 1 × T matrix
corresponding to the T samples recorded in electrode lead i, and N is the number of channels (electrode leads).

Assuming a linear model, X = AS, let's denote the columns of A as

a j = [a 1 [j] . . . a N [j]] �
and the rows of S as

s � i = [s i [1], • • • , s i [T ]]
where they represent the j-th spatial pattern and i-th temporal sources, respectively. Here, a i [j] shows the contribution of the ith temporal source, s � i , to the jth electrode lead, x � j .

Heart of Reference-Based Source Separation Method

Our objective is to extract the sources related to a reference activation model. We consider two states, denoted C 1 and C 2 , which correspond to the reference and non-reference activations, respectively. Denoting T � , � = 1, 2, the set of time samples related to each state, we can build the corresponding segment matrix, X � ∈ R N ×M � .The correlation matrix of data for each state can be estimated as:

� R � = 1 M � X � X � � . (4.16)
The spatial filters, W (whose columns are generically denoted w), for which the temporal sources, S = W � X have maximum similarity with the reference activation state, i.e., maximum variance in the reference state compared to the other state, is computed as:

max W w � � R 1 w w � � R 2 w (4.17)
Using (4.3) for solving (4.17

) leads to GEVD of ( � R 1 , � R 2 ): � R 1 W = � R 2 WΛ (4.18)
Using W, the spatial patterns, A = (W � ) -1 , and the temporal sources, S = W � X, are extracted.

As explained before, the maximum eigenvalue in (4.18) is related to the maximum power ratio in (4.17). We rank the eigenvalues in decreasing order. This implies ranking of the estimated temporal sources, according to their resemblance to the reference activation state.

Sources are uncorrelated

In R-SS the observations are segmented into two parts. In EEG analysis for reducing the noise and the background the average of the epochs is used. Therefore, we naturally have two parts in the observation. The equivalent observation matrix in R-SS is produced as:

X = [X 1 X 2 ] (4.19)
The source matrix and the transform matrix W is etimated by R-SS.

S = WX = W[X 1 X 2 ] (4.20)
Here, we will show that the estimated sources are uncorrelated.

SS � = WXX � W � = W[X 1 X 2 ] � X 1 � X 2 � � W � = W(R 1 + R 2 )W � (4.21)
Using theorem 4.1 it is proven that the correlation matrix of the sources is diagonal and without loss of generality the temporal sources are normalized, which means that SS � = I, where I is an N × N identity matrix.

Comparison with previous applications of GEVD

GEVD has been used before in different applications like common spatial patter (CSP), nonstationary source separation and periodic independent component analysis. Here we want to review these methods to clarify the idea behind them. This classification is taken from [START_REF] Sameni | A Deflation Procedure for Subspace Decomposition[END_REF].

� SNR Maximization: Suppose that the observation x = x � i , the ith row of X in Eq. 3.1, has an additive model as follows:

x = x s + x n
where x s is the desired part of the observation, and x n is the undesired part, considered as noise and/or interference. The goal is to find a linear transform y = w � x to obtain a transform with the highest SNR value.

SN R(w)

= y � s y s y � n y n = w � R xs w w � R xn w where R xs = 1 T x � s x s and R xn = 1 T x � n
x n are the correlation matrix of desired and undesired parts of the observation. The maximum value of the SNR is achieved by GEVD of the matrix pair (R xs , R xn ).

� Nonstationarity Maximization: Suppose that the observations are the nonstationary mixtures of the desired and undesired signals. For example, in Brain Computer Interface (BCI), the goal is to find the brain patterns related to two different tasks like left and right hand movements. In this application, Common Spatial Patterns (CSP) are estimated using GEVD of the correlation matrices of different epochs of the EEG signal [START_REF] Blankertz | Optimizing Spatial filters for Robust EEG Single-Trial Analysis[END_REF]. In our work the reference and non-reference states are two different time intervals with different statistical measures. Extracting artifacts from EEG which is recorded simultaneously with fMRI also can be considered as nonstationary signals [START_REF] Amini | MR Artifact Reduction in the Simultaneous Acquisition of EEG and fMRI of Epileptic Patients[END_REF].

In all these applications the spatial filter is derived by GEVD of the correlation matrices of two time intervals. The choice of time intervals is dependent on the application.

� Spectral Contrast Maximization: Sometimes undesired part of the observation is separable in other space like frequency domain. In such a case we use GEVD of the cross-spectrum of the observations in a specific frequency bandwidth and the covariance of the whole observations considering all frequencies.

Eq. 4.6 can be transformed as follows:

Y (f ) = F {y(t)} = w � F {x(t)} = w � X(f )
and GEVD is written as:

σ(w) = E ν {|Y (ν)| 2 } E f {|Y (f )| 2 } = w � S x w w � C x w
where S x = E ν {X ν X H ν } is the cross-spectrum of the observation vector averaged over the bandwidth of interest. See [START_REF] Sameni | A Deflation Procedure for Subspace Decomposition[END_REF] for more details.

� Periodicity Maximization: When the desired or undesired parts of the observations have a periodic or quasi-periodic structure, one may seek for a linear transform that maximizes a measure of periodicity, while keeping the signal energy bounded:

�(w) = E t {y(t)y(t + τ t )} E t {y(t) 2 } = w � P x w w � C x w
The MUSE algorithm assumes a constant value for τ t [START_REF] Tong | Indeterminacy and Identifiability of Blind Identification[END_REF] and [START_REF] Sameni | Multichannel Electrocardiogram Decomposition Using Periodic Component Analysis[END_REF] assumes a time varying value for τ t .

Source Classification

Blind source separation (BSS) methods are powerful tools for separating sources in many applications. Estimating the number of sources provides useful information for BSS methods.

Applying Principle Component Analysis (PCA) and choosing the number of dominant eigenvalues of the spatial correlation matrix as the number of sources is a well known method. However, it is difficult to distinguish dominant eigenvalues. Various methods based on Information Theoretic Criteria (ITC) have been proposed to estimate the principle dimension of multivariate data.

Several methods based on ITC have been demonstrated to be attractive for model order selection in signal processing, including Akaike's information criterion(AIC) [?], Kullback-Leibler information criterion(KIC) [START_REF] Cavanaugh | A large-sample model selection criterion based on Kullback's symmetric divergence[END_REF], the minimum description length (MDL) criterion [START_REF] Rissanen | Modelling by the shortest data description[END_REF], Bayesian information criterion (BIC) [START_REF] Schwarz | Estimating the Dimension of a Model[END_REF] and a Laplace approximation to Bayesian Criterion based on model evidence [START_REF] Minka | Automatic choice of dimensionality for PCA[END_REF]. All of these methods have high computations. These criteria are investigated in [START_REF] Hui | An Empirical Comparison of Information-Theoretic Criteria in Estimating the Number of Independent Components of fMRI Data[END_REF].

Here, we propose a simple method to define the number of sources. Remember that the simplicity of this method is due to the usage of experiment information in the source separation step and it is achieved by a special interpretation of the eigenvalues which will be explained below.

After obtaining the discriminative sources (s i ) between the reference and non-reference states ranked according to their similarity to the reference state. Now, we need to select the sources, which are similar enough to the reference, for being considered as belonging to the reference class. To this end, we propose the following procedure.

The probability of the reference class (ω 1 ) membership is calculated as follows:

p(s � i ∈ ω 1 ) = λ i max(λ j=1,...,N ) (4.22)
where λ i , i = 1, . . . , N indicate the eigenvalues or the diagonal elements of Λ in (B.1). For classification of the sources, the error probability using Bayes rule is defined as:

p error = N � j=1 � p(s � j ∈ ω 2 |ω 1 )p(ω 1 ) � + N � j=1 � p(s � j ∈ ω 1 |ω 2 )p(ω 2 ) � where p(s j ∈ ω 2 |ω 1 ) = 1-p(s j ∈ ω 1 ) and p(s j ∈ ω 1 |ω 2 ) = 1-p(s j ∈ ω 2 )
as ω 1 and ω 2 constitute a partition. p(ω 1 ) and p(ω 2 ) are the prior probabilities of the reference and non-reference classes, respectively. We remind that GEVD sorts the separated sources in decreasing order of similarity with respect to the reference. Therefore, if we assume that only the first i sources belong to the reference class (and consequently the Ni others belong to the non-reference class), then the total error probability (false positive plus false negative errors) can be written as follows:

p error (i) = � i j=1 p(s � j ∈ ω 2 |ω 1 )p(ω 1 ) + � N j=i+1 p(s � j ∈ ω 1 |ω 2 )p(ω 2 ) (4.23)
where p(ω 1 ) = i N and p(ω 2 ) = N -i N . Thus, the minimum of p error (i) provides the number i * of

X 1 m M 1 m N m = 1 . . . L 1 R 1 m N N m = 1 . . . L 1 R 1 N N X 2 m M 2 m N m = 1 . . . L 2 R 2 m N N m = 1 . . . L 2 R 2 N N GEVD S T N s i Figure 4.1: R-SS method. R 1 m and R 2 m , m = 1, .
. . , L � denote the correlation matrices of different IED and non-IED time intervals, respectively. The average of these matrices denoted as R 1 and R 2 are the inputs of GEVD. Using GEVD, the discriminative temporal sources, s i , i = 1, . . . , N , are estimated. These sources are ranked according to their similarity to the IED class.

the sources in the reference class, i.e., i * = arg min i p error (i).

Conclusion

We propose a new reference-based source separation method using generalized eigenvalue decomposition for identification of brain regions involved in a reference brain state (like interictal events). Using GEVD, we estimate the spatial filter that maximizes the power of reference over non-reference state, and provides temporal sources ranked according to their similarity to reference state. R-SS method can be considered as a semi-blind source separation as it uses the states in which the observations are recorded. R-SS is based on GEVD algorithm. Using GEVD has two advantages: First, its solution can be computed analytically. Therefore, it is fast. Second, GEVD in a manner that is used in R-SS method gives extra information about the sources. It sorts the sources in increasing order of their similarity to the reference state. R-SS uses this extra information in source classification step to find the sources of interest. Here, we use epileptic EEG to evaluate the R-SS method to find out how much it would succeed to localize the epileptic sources. The drug-resistant focal epileptic patients are recommended for resective surgery. The goal of this surgery is to remove the brain regions responsible for the epileptic seizures [START_REF] Rosenow | Presurgical Evaluation Of Epilepsy[END_REF]. The classic method used by the epileptologists for the presurgery evaluations is studying seizure onset zones (SOZ), i.e. the regions where the first electrophysiological changes are detected at ictal onset. However, if the seizures do not occur, or occur rarely, the presurgical evaluation cannot be performed or is prolonged. Moreover, as seizures are rare events, their study may not lead to statistically robust results [START_REF] Alarcon | Origin And Propagation Of Interictal Discharges In The Acute Electrocorticogram. Implications For Pathophysiology And Surgical Treatment Of Temporal Lobe Epilepsy[END_REF]. For these reasons, studying the interictal epileptiform discharges (IED) is very valuable as IEDs are frequent events during EEG recording. The relationship between SOZs and IED regions is an open question and there are several studies wondering whether identification of IED regions can be used to guide the resection surgical decision [START_REF] Alarcon | Electrophysiological Aspects Of Interictal And Ictal Activity In Human Partial Epilepsy[END_REF][START_REF] Hufnagel | Clinical Relevance Of Quantified Intracranial Interictal Spike Activity In Presurgical Evaluation Of Epilepsy[END_REF][START_REF] Bourien | A Method To Identify Reproducible Subsets Of Co-Activated Structures During Interictal Spikes. Application To Intracerebral EEG In Temporal Lobe Epilepsy[END_REF][START_REF] Marsh | Interictal EEG Spikes Identify The Region Of Electrographic Seizure Onset In Some, But Not All, Pediatric Epilepsy Patients[END_REF].

The complexity of IED identification lies within reducing the effect of background activity using a robust method. Toward solving this problem, a method based on functional connectivity was proposed in [START_REF] Amini | Directed Differential Connectivity Graph of Interictal Epileptiform Discharges[END_REF]. Contrary to previous methods, both IED and non-IED time intervals were used in [START_REF] Amini | Directed Differential Connectivity Graph of Interictal Epileptiform Discharges[END_REF] to develop a differential connectivity graph (DCG). To extract the statistically robust connections, DCG uses a permutation method applied to a large number of IED and non-IED time intervals which is computationally heavy, while reference-based source separation (R-SS) method efficiently reduces the computational load. R-SS method has been presented in Chapter 4.

Here, for estimating the IED regions which are not necessarily unique, we applied a multiobjective optimization algorithm. Evaluation of our method is based on both realistic simulated data [START_REF] Cosandier-Rimélé | A Physiologically Plausible Spatio-Temporal Model For EEG Signals Recorded With Intracerebral Electrodes In Human Partial Epilepsy[END_REF] and on actual data. For actual data, the IED regions estimated based on our method are compared with visually inspected SOZ and leading IED regions estimated by DCG-based method introduced in [161] using the same iEEG recordings.

We introduce the parameters of the R-SS method and also the optimization method in Section 5.2. The proposed method is evaluated using two sets of data: simulated and actual data in Section 5.3 and 5.4, respectively.

Method

The steps of the method are as follows: 1) Preprocessing: manual labeling, and band-pass filtering; 2) R-SS: source separation, and source classification; 3) Feature Extraction; 4) Optimization.

In the following, each step of the method is explained in details.

Preprocessing

In the following, we explain the preprocessing steps. More details about these steps can be found in [START_REF] Amini | Directed Differential Connectivity Graph of Interictal Epileptiform Discharges[END_REF].

Manual labeling

The IED and non-IED time intervals are manually identified by the epileptologist for each patient considering all the iEEG channels. An IED period is a time interval, including a single IED or burst of IEDs. A non-IED period is a time interval without any IED or abnormal events. IED periods studied for each patient are homogeneous in terms of their characteristics. However, the IED periods may include single IED or burst of IEDs.

Band-pass filtering

The estimated power spectrum density of IED intervals is large in the frequency range of 4-64 Hz. Therefore, a band-pass filter with a passband from 4 to 64 Hz is applied to the temporal signals. According to the manual labeling of previous step, the filtered signal is segmented into L 1 IED and L 2 non-IED segments.

R-SS

In this step, we use the R-SS method to extract the most similar sources to the reference model.

Data

The observation data is denoted as

X = [x 1 • • • x N ] � ∈ R N ×T , where x j = [x j [1], • • • , x j [T ]] � , j = 1, • • • , N is a zero-mean T × 1 matrix
corresponding to the T samples recorded on electrode lead j, and N is the number of channels (electrode leads).

Here, our objective is to extract the sources related to a reference activation model or IED event. We consider two states, denoted C 1 and C 2 , which correspond to the reference and nonreference activations like IED and non-IED, respectively. Denoting T � , � = 1, 2, the set of time samples related to each state, we can define each column of the corresponding segment matrix, X � ∈ R N ×M � , as:

χ i = [x 1 [i], • • • , x N [i]] � , i ∈ T � , M � = card(T � ) (5.1)
where card(.) indicates the cardinality of a set. The correlation matrix of data for each state can be estimated as:

� R � = 1 M � X � X � � . (5.2)
The spatial filters, W, for which the temporal sources, S = W � X have maximum similarity with IEDs and the number of the IED related sources (i * ), are computed by the R-SS method.

The scheme depicted in Fig. 4.1 demonstrates the block diagram of the reference-based source separation method. The L 1 reference (or IED) and L 2 non-IED segment matrices are denoted as

X � m ∈ R N ×M � , where m = 1, . . . , L � is the index of each segment. Each column of X � m is χ i , i ∈ T � m , M � m = card(T � m ).
T � m indicates the set of time samples of the mth segment related to each state. � is equal to 1 and 2 for indicating IED and non-IED, respectively. The correlation matrix for each of these IED and non-IED segment matrices are calculated using (B.2) by substituting X � with X � m and M � with M � m . These correlation matrices are denoted as R � m ∈ R N ×N . The average of L 1 IED correlation matrices and L 2 non-IED correlation matrices denoted as R 1 and R 2 are the inputs of GEVD.

The output of GEVD is the ranked spatial patterns and temporal sources according to their similarity to reference or IED state. The sources s i , i = 1, . . . , i * are identified as IED class members. After obtaining i * , as we are not interested in the non-IED class members, the probability of the IED class membership (4.22) can be rewritten as:

p(s i ∈ ω 1 ) = � λ i λ 1 i = 1, • • • , i * 0 i = i * + 1, • • • , N.
(5.3)

Feature extraction

In the previous step, the source members of IED class are identified. Now, for the identification of IED electrode leads, we need to transfer this information from the source space to the observation space.

According to the linear relationship between observations and sources and to the constraint �w� = 1 considered in (4.2), the probability of activation of each source, s i , in each observation,

x j (associated with electrode lead j), corresponds to its relative power contribution to x j :

x j = N � i=1 a i [j]s i ⇒ p(x j |s i ) = a i [j] 2 � N i=1 a i [j] 2
(5.4) Using (5.4) and (5.3), we define the membership probability of each observation for the IED class, p(x j ∈ C 1 ), as follows:

p(x j ∈ C 1 ) = i * � i=1 p(x j |s i )p(s i ∈ ω 1 ) = i * � i=1 p i j (5.5)
One can estimate the IED regions by maximizing (5.5) as a single objective function, i.e. e = arg max

j p(x j ∈ C 1 ) (5.6)
where e indicates the electrode lead number involved in the IED event. Estimating the IED electrode leads through (5.6) has the following two problems. First, (5.6) is a single objective optimization problem and provides one single IED lead, while most often there exist a network of structurally and functionally connected brain regions [START_REF] Bragin | Chronic Epileptogenesis Requires Development of a Network of Pathologically Interconnected Neuron Clusters: A Hypothesis[END_REF], thus we expect a set of IED leads close to these regions. However, to obtain a set of IED leads, one can select the leads whose IED class membership probability (p(x j ∈ C 1 )) is greater than a given threshold, which provides threshold dependent results. Second, in (5.5) the probability of activation of different sources of IED class (s i , i = 1, . . . , i * ) in each electrode lead are averaged, while we are interested in optimizing the activation of each source in each electrode. These problems of single objective optimization methods are addressed in [START_REF] Deb | Multi-Objective Evolutionary Algorithms: Introducing Bias Among Pareto-Optimal Solutions[END_REF][START_REF] Branke | Multiobjective Optimization, Interactive And Evolutionary Approaches[END_REF]. Historically, Pareto introduced multi-objective optimization methods in which all the objective functions are considered simultaneously.

Optimization

The most important reason for rejecting the single objective optimization solution is that, we are interested in optimizing the contribution of each source to each electrode lead. Let us denote the membership probability of each observation or node, j, for the i * sources of IED class as

p j = � p 1 j , p 2 j , . . . , p i * j � �
, p i j = p(x j |s i )p(s i ∈ ω 1 ), j = 1, . . . , N . In an ideal situation, the objective function values, p i j should achieve the maximum value for each individual dimension [START_REF] Deb | Multi-Objective Evolutionary Algorithms: Introducing Bias Among Pareto-Optimal Solutions[END_REF], i.e.

∃z ∈ {1, . . . , N },

p i z = max j (p i j ), ∀i = 1, . . . , i * p z = [p 1 z , . . . , p i * z ] � .
(5.7)

However, practically the ideal point z does not exist in our search space. Instead, more often there exists a set of solutions (leads), where each solution receives a greater contribution from at least one of the IED sources. This set of solutions is called a set of non-dominated solutions (or layer) in Pareto optimization. This concept of Pareto optimality [START_REF] Deb | Multi-Objective Evolutionary Algorithms: Introducing Bias Among Pareto-Optimal Solutions[END_REF][START_REF] Branke | Multiobjective Optimization, Interactive And Evolutionary Approaches[END_REF], has been used in many applications as in economics, management science, mechanical engineering, etc.

The multi-objective optimization problem, in Pareto sense, gets the following form: maximize

� p 1 j , p 2 j , . . . , p i * j � subject to p j ∈ P (5.8)
consisting of i * objective functions that are aimed to be maximized simultaneously. The decision/variable vector p j ∈ R i * belongs to the search space P ⊂ R i * . We classify the search space P according to the Pareto concept of non-domination [START_REF] Deb | Multi-Objective Evolutionary Algorithms: Introducing Bias Among Pareto-Optimal Solutions[END_REF]: a node is a member of nondominated layer if either it dominates the others, or there is no other node dominating it. Node j dominates node j 0 , if ∀i, p i j ≥ p i j 0 , and ∃i 0 , p i 0 j > p i 0 j 0

1 . The first non-dominated layer, D(P ) is obtained from the nodes of the entire search space P . In the following, we explain how to estimate D(P ) using Pareto optimization algorithm [START_REF] Deb | Multi-Objective Evolutionary Algorithms: Introducing Bias Among Pareto-Optimal Solutions[END_REF].

Let us consider N i * -dimensional decision vectors, p j , as N nodes in the search space P .

1. Initialize D(P ) with the first node (j = 1) with value of p 1 . The initialization can be done with any node of search space.

2. Choose a new node (j = j + 1):

(a) If any node in D(P ) dominates node j go to step 3.

(b) Else add node j to D(P ) and remove any nodes of D(P ) that node j dominates.

3. If j is not equal to N go to step 2.

The members of D(P ) build the first non-dominated layer. We can then compute a second nondominated layer, etc. each one corresponding to a less probable set of solutions. For the kth non-dominated layer, we remove the nodes of first (k -1)th non-dominated layers from search space P , and we use the same algorithm on the remaining nodes, P k .

In our problem, we consider the estimated results as a hint for the epileptologists for presurgical evaluations. In our analysis, by classifying the search space into different non-dominated layers [START_REF] Deb | Multi-Objective Evolutionary Algorithms: Introducing Bias Among Pareto-Optimal Solutions[END_REF], we can estimate different sets of close electrode leads to at least one of the epileptic sources according to their probability. From the first to the last layer, the probability of being close to at least one epileptic source decreases. Of course, the first layer includes the most probable solutions in comparison with other layers. Therefore, we propose the first layer, and to enlarge the set of solutions, we suggest the second layer conditioning that it is spatially close enough to the first one. If the latter condition is satisfied, the union of the first and the second layer nodes are defined as the "Pareto optimal set" or the set of estimated IED electrode leads.

The distance between first and second layers is computed using Hausdorff distance (d max ), which is the supremum of minimum Euclidean distances between the first and second non-dominated layers' nodes. By comparing d max with a relative margin (threshold), we test if first and second layers are close enough. In the following, we explain how the relative margin is calculated.

Practically, ideal point z (5.7) might or might not be in search space P . However, for each search space P , we can calculate its related �p z � that is dependent on the maximum objective function values: max j (p i j ), and eventually on the first layer nodes. A relative margin: � of �p z � (� is a small positive number) is considered to measure the closeness of first and second layers.

Therefore, if d max ≤ � �p z �, i.e. the second layer is within the relative margin of first layer (100�% of �p z �), the set of estimated IED leads are suggested as the union of the first and second layer nodes. Otherwise, only the first layer nodes are considered.

The Pareto optimal solution provides a set of electrode leads, while epileptologists would like to get the brain regions. In our data, we assume that intracerebral electrodes are inserted into the clinically suspected brain area. This assumption is validated based on the expertise of the epileptologists. In this work, the IED leads are associated with the brain regions (where the leads are located) by the epileptologists using the implantation scheme and other clinical information. For more accurate association, a source localization method could be applied to the intracerebral data of the selected electrode leads to solve the associated inverse problem as it has been done in [START_REF] Chang | Dipole Localization Using Simulated Intracerebral EEG[END_REF]. This would be more robust but it is out of the scope of this paper.

Simulated Data

The efficiency of the proposed method is evaluated using computer simulations. As illustrated in Fig. 5 ). Model parameters (related to neuronal excitability) are tuned such that epileptic source dipoles are assigned an epileptic time-course (spiking activity), while non-epileptic source dipoles are assigned a 'normal' time-course (background activity). For epileptic activity, it is assumed that spikes are originated in source e 1 and propagated to source e 2 with a delay of 30-50 ms. The EEG signals, produced by brain sources at all depth-electrode leads, are simulated by solving the 'EEG forward problem' (see [START_REF] Cosandier-Rimélé | A Physiologically Plausible Spatio-Temporal Model For EEG Signals Recorded With Intracerebral Electrodes In Human Partial Epilepsy[END_REF] for details). Briefly, we assume an infinite homogeneous medium with conductivity σ = 33 × 10 -5 S/mm. At each electrode lead E , the electric potential V produced at time t by a current dipole is V E (t) = (m(t).u r )/(4πσr 2 ), where m(t) = m(t)d is the dipole moment (with magnitude m(t) and orientation d), u r is a unit vector oriented from the source to the electrode lead, r is the distance between them, and the potentials produced by individual sources add up linearly.

The duration of the simulations is 600 s (sampling frequency: f s = 512 Hz). In practice, for each simulation, 100 IED time intervals (length: 300 samples each), centered on the peaks of IEDs, are extracted. Each IED time interval includes a single IED (here, a spike). For each of the three simulations, 100 non-IED time intervals (same length as the IED intervals) are also extracted.

The method is applied to the three simulated data. The results are shown in Fig. 5.1, using thick and thin frames for the first and second Pareto layers, respectively. It can be seen that the union of first and second Pareto layers includes the closest electrode leads to at least one of the epileptic sources for the three simulated data. There is no major difference between the results related to the three data of different oriented epileptic sources. Electrode lead A 2 is not selected as there are other non-epileptic sources close to this electrode, thus the contribution of epileptic source e 1 to A 2 is less than its contribution to A 0 and A 1 . More discussion on the results in terms of Pareto layers can be found in Section 5.

The proposed method has been evaluated using more sophisticated simulations in 3-D and a number of electrodes (not reported in this paper due to lack of place), which provided congruent results.

To study the effect of signal to interference ratio2 (SIR) on the proposed method, the simulations are repeated for different SIR values. Different SIRs are generated by changing the contribution of the non-epileptic sources (generating background activity) to the simulated iEEG.

For SIR ∈ [-2, +20] dB, we get congruent results in all directions, although the IEDs are not visible in the iEEG signals with the lowest SIR value. More specifically, the identified electrode leads are identical for all SIR values with some changes in the assignment to the first or the second Pareto layer. Indeed, for example, in the simulated data with source orientation D0, for the highest SIR value, all the identified electrode leads (Fig. 2) are assigned to the first layer.

Decreasing SIR causes A2 to be re-assigned from the first layer to the second layer, and for the lowest SIR value, B0 is also assigned to the second layer. Table 5.1 shows the detail of the results for different SNR values.

We have also simulated 3D electrode implantation as in figure 5.2 for the same sources in the 2D simulation and SIR=5dB. The electrode leads which are detected as epileptic ones are shown in table 5.2.

Using simulated data presented in this paper, we qualitatively demonstrated that the proposed algorithm works properly. However, a complete study, including the analysis of different parameters and their related quantitative evaluations will be done in our future work.

Actual Data

The iEEG recordings were obtained from five patients suffering from focal epilepsy. The types of focal epilepsy of the patients are brought in Table 5.3. The patients underwent presurgery evaluations with iEEG recordings. They are seizure free after resective surgery. Eleven to fifteen semirigid multi-lead intracerebral electrodes with 0.8 mm diameter were bilaterally implanted in suspected seizure origins based on clinical considerations. The multi-lead electrodes (Dixi, Besançon, France) include 5, 10, 15, or 18 leads. Each lead has 2 mm length and is evenly spaced with interspace of 1.5 mm. The iEEG was recorded with an audio-video-EEG monitoring system (Micromed, Treviso, Italy) with a maximum of 128 channels and sampled at 512 Hz.

Electrode lead number

Low SIR D 0 1st Layer A 0 B 0 C 9 2nd Layer A 1 D 1 1st Layer A 0 2nd Layer B 0 C 9 D 2 1st Layer A 0 C 9 2nd Layer A 1 B 0 B 1 Medium SIR D 0 1st Layer A 0 B 0 C 9 2nd Layer A 1 D 1 1st Layer A 0 C 9 2nd Layer B 0 C 8 D 2 1st Layer A 0 C 9 2nd Layer A 1 B 0 B 1 High SIR D 0 1st Layer A 0 B 0 C 9 2nd Layer A 1 C 8 D 1 1st Layer A 0 2nd Layer B 0 C 9 D 2 1st Layer A 0 C 9 2nd Layer A 1 B 0 B 1 C 8

Comparison with other methods

We compare the IED regions estimated by R-SS method with two methods: visually inspected SOZ (vSOZ) by the epileptologist and leading IED regions (�IED) [START_REF] Amini | Directed Differential Connectivity Graph of Interictal Epileptiform Discharges[END_REF].

Comparison with vSOZ

In Table 5.5, the IED regions estimated by R-SS method and Pareto optimization method are reported. Here, we experimentally choose � = 0.3 from our data, i.e. the second Pareto layer can be considered as the neighborhood of the first Pareto layer if d max is smaller than or equal to 30% of �p z � (5.7). The values of d max / �p z � for our five patients are as follows: 0.16, 0.09, 0.01, 0.17, 0.6. In our experiments, d max ≤ � �p z � is satisfied for patients 1 to 4.

Therefore, the Pareto optimal solutions reported in Table 5.5 for these patients are the union of the first and second non-dominated layer members while, for patient 5, the regions are associated with the first non-dominated layer members. Patient 5 is a particular patient with very focused vSOZ, which may explain why the Pareto optimal solution is limited to only the first non- dominated layer and eventually the related IED regions are very focused.

In the following we present the results in terms of regions where the electrode leads are located since it is more significant. The abbreviations used in Table 5.5 are as follows: amygdala (amyg); anterior/posterior/internal/superior (ant/post/int/sup); entorhinal cortex (entC); hippocampus (Hc); parahippocampal gyrus (pHcG); temporal (T); temporal pole (TP); mesial (m); gyrus (G); middle short gyrus of insula (midInsG); patient i (P i ).

The comparison reported in Table 5.5 shows the congruency between vSOZ and the estimated IED regions by our method except for the patient P 2 . For P 2 , there are vSOZs which are not estimated by R-SS. It is important to mention that the suggested vSOZs, which we considered as ground truth, are based on EEG and extra clinical information such as semiology. Moreover, since all patients are seizure-free after resective surgery, we deduce that the removed regions included the necessary regions for creating the seizures, but the removed regions might be more However, the estimation of IED regions, including the vSOZ (zero false positive, or precision = 100%, see Table 5.6) is valuable since vSOZs are always included in the removed brain regions.

Precision and sensitivity are defined later in this Section.

Comparison with �IED

In Table 5.5, the IED regions estimated by our method are compared with the �IED regions estimated in [START_REF] Amini | Directed Differential Connectivity Graph of Interictal Epileptiform Discharges[END_REF] using directed DCG (dDCG), where the causal relationships were considered.

The �IED regions are the estimates of leading or source regions involved in IED events where the IED signals are assumed to be originated. The results of our method and �IED are congruent (Table 5.5) except for a few regions: mesial temporal pole in P 1 , parahippocampal gyrus in P 3 , and Amygdala in P 4 . One may interpret that the two former regions that are not included in �IED regions could be transit or sink regions which are involved in the IED event, and not necessarily originating IED signals [START_REF] Alarcon | Origin And Propagation Of Interictal Discharges In The Acute Electrocorticogram. Implications For Pathophysiology And Surgical Treatment Of Temporal Lobe Epilepsy[END_REF][START_REF] Amini | Development of Differential Connectivity Graph for Characterization of Brain Regions Involved in Epilepsy[END_REF]. Amygdala is not selected by our method in P 4 .

However, mesial temporal pole in P 1 , parahippocampal gyrus in P 3 , and amygdala in P 4 are included in vSOZ and were removed during surgery.

Our method provides a greater congruency with vSOZ in comparison with �IED regions.

This result is shown quantitatively in our method and �IED regions are compared with vSOZ by assuming that vSOZs are the ground truth. For this purpose, the following measures are used.

1. P recision = T P T P +F P (in %), where T P and F P indicate true positive and false positive in terms of brain regions. T P is the number of common brain regions between vSOZ and our estimated IED regions, while F P is the number of uncommon regions which were falsely detected by our method.

Sensitivity = T P

T P +F N (in %), where F N indicates the false negative in terms of brain regions. F N is the number of regions missed by our method. Conversely, ovp = 0% shows that we could not obtain any IED electrode lead in proximity of any of vSOZ electrode leads. 5. ovp2: the average percentage of the number of vSOZ electrode leads that are in the neighborhood (< 1.5cm) of at least one of the IED electrode leads. ovp2 is the same as ovp except that the set of vSOZ electrode leads is replaced with the set of IED electrode leads.

Evaluation of the proposed method in terms of robustness

The robustness of the method is evaluated according to the influence of the number of IED time intervals and errors in identification of IED time intervals.

1) Number of IED time intervals:

The proposed GEVD-based method is fast due to its exact analytic solution, but the correct estimation of correlation matrices is important for the reliability of the results. Large enough number of data samples in each IED or non-IED time interval is needed for a proper estimation of correlations. The length of each IED time interval depends on the length of single or bursts of IEDs. Here, the mean of minimum and maximum length of IED time intervals over patients are equal to 236 and 3.3 × 10 3 samples, respectively which provides statistically reliable estimation of correlation matrices. The length of the required recorded data for processing is dependent on the number of IED and non-IED time intervals that exist in this length of data. It means that if in a given set of data there does not exist enough number of IED and non-IED time intervals, we need longer set of data to obtain the sufficient number of IED and non-IED time intervals. In our recordings, the mean of the length of selected data is about 1 hour. To reduce the number of IED or non-IED time intervals (not the length of each time interval), and eventually reduce the required time length of data for recording, we test the reliability of the method in terms of the number of IED time intervals. For this purpose, we use the jackknife method as follows. First, a percentage (L) of available IED and non-IED time intervals is selected randomly. Next, the IED electrode leads related to these time intervals (IED # ) are estimated. Finally, the related false positive and false negative (F P e and F N e ) values are calculated in terms of electrode leads (for more details, see the Appendix). By repeating this procedure 100 times, the mean and standard deviation of F P e and F N e are calculated (left column of Fig. 5.4). We assume that the set of IED electrode leads using L = 100% of available IED and non-IED time intervals is the ground truth. L is set to 30%, and 70%. Therefore, we consider that the results of the method are reliable for L < 100%, if there is no significant change between the results related to L < 100% and L = 100%. Results of Fig. 5.4 (left column) can show the reliability of the method for a small number of IED time intervals (about 50) for most of the patients (P 1 -P 4 ). For P 5 , the standard deviation is higher in comparison with other patients, which reveals its sensitivity for small number of IEDs.

× × × × �IED × × × × R-SS × × × × × P 2 antHC postHC pHcG amyg vSOZ × × × × �IED × R-SS × P 3 antHC postHC pHcG vSOZ × × × �IED × × R-SS × × × P 4 antHC postHC amyg mTP entC vSOZ × × × × × �IED × × × × R-SS × × × P 5 midInsG vSOZ × �IED × R-SS × Table 5
2) Errors in identification of IED time intervals: To test the effect of possible errors in identification of IED time intervals (due to expert variability, or in the future the error in automatic labeling), we exchange the labels of a randomly selected percentage (E) of available IED and non-IED time intervals, and estimate the related IED electrode leads (IED # ). The same jackknife-based method explained above is applied. Here, we assume the ground truth is the set of IED electrode leads estimating by using original IED and non-IED labels, i.e. E = 0%.

E is set equal to 10%, 20%, 30%, 40%, and 50%. Therefore, if there is no decrease of performance for each of the former error percentages in comparison with E = 0%, we can conclude that the method can tolerate this percentage of error in IED labelling. Fig. 5.4 (right column) shows the mean and standard deviation of F P e and F N e for all the patients. The result shows the importance of the discriminative information between the two IED and non-IED states and how it is crucial for obtaining the proper results. According to Fig. 5.4 (right column), one can conclude that 10% error in labeling again provides correct results for most of patients: only the results of P 1 are very sensitive to the error of labeling.

Discussion

Simulation Evaluation

The is a simplified example of real data, where we assume the presence of two epileptic sources in the volume around the electrodes, supposed to be inserted in a clinically suspected brain area.

We are not looking for a unique solution in our method, instead we are interested in extracting a set of solutions, which the number of solutions is practically unknown. Using Pareto optimality concept and eventual classification of the search space into different non-dominated layers, we aim to estimate the set of electrode leads that are the closest to at least one of the epileptic sources. As it can be seen in Fig. 5.1, the first non-dominated layer includes mostly the electrodes that are the closest to epileptic sources in the three simulations with different orientations of epileptic sources. These leads are not necessarily close to a single source, instead these are the leads which receive a greater contribution from at least one of the epileptic sources, and eventually they have a higher membership probability to IED class. Since d max ≤ � �p z � (see section 5.2.4), and thus the second layer is close enough to the first layer, we enlarge the set of results by considering the former one. Therefore, we obtain the extra electrode leads in the neighborhood, such as A 1 , C 8 , and B 0 . By considering the second layer results, we can see if there is any solution in the neighborhood (even with less membership probability in comparison with first layer) which might be useful for the epileptologists. In the current work, we do not identify the location and orientation of the dipole sources. However, these extra information is very helpful for future works for localizing the dipole sources in respect to the location of neighborhood electrode leads.

For the three simulations of different orientations, most of the estimated IED leads are common. Some of the differences between the set of results are electrode leads A 1 and A 2 . A 1 is identified for the three orientations except D 1 . This is the specific condition in which the angle between the source dipole e 1 and position vector of A 1 are perpendicular. A 2 is less specific since there are other non-epileptic sources close to this lead.

The method performs well over a large range of SIR values (SIR ∈ [-2, +20] dB), as the same set of electrode leads are always identified. Although decreasing SIR decreases the number of electrodes in the first layer, all the significant electrodes (A 0 and C 9 ) are identified for even the lowest SIR value.

The results of 3D implantation are congruent with the 2D simulation. It shows that although when the electrical field of the sources are perpendicular to any electrode lead, it would not record any signal, but the positions of the electrode leads can be considered as a good estimation of the places of the sources. Using different electrodes with different orientations will arise the accuracy of the estimation accuracy.

Advantages of using GEVD

In most usual BSS methods, identification of sources related to the desired state is done in two steps. The sources are first estimated, and secondly the sources of interest are selected usually by correlation with the reference. On the contrary, in GEVD the two tasks are done in a unique step, since the estimated sources are ranked according to their similarity to the reference model. Therefore, we only need to identify the number of related sources. In addition, GEVD does not require the assumption of either independence or non-Gaussianity as it is done in independent component analysis (ICA) methods. However, it would obviously lead to decorrelated sources.

Since in our application, the sources are not necessarily independent, GEVD is more flexible.

Finally, solving GEVD is very fast and only based on second order statistics, contrary to ICA which requires higher order statistics.

Comparison between dDCG-based and R-SS methods

dDCG-based identification method considers On the contrary, R-SS is simple and fast as its solution is based on an analytical mathematical problem, i.e. GEVD. This property is valuable for an online process of iEEG recordings. The processing time of R-SS method for about 40 minutes of recording which provides 195 IED and non-IED time intervals and for 111 channels is about 4 minutes, while it is about 100 times longer for dDCG (using a shared 3 GHz, 4 core Xeon 64 bit processor). Another advantage of R-SS over dDCG is that R-SS can provide reliable results for less number of IED and non-IED time intervals (about 50 time intervals). For construction of dDCG that is based on multiple test using permutations, to obtain statistically reliable results, about 150 time intervals are needed for false alarm rate equal to 0.05, 100 connections, and 10,000 permutations. R-SS just requires enough number of samples in each time interval for a statistically reliable estimation of correlation matrices.

One current limitation of both dDCG and R-SS methods is using the manual labeling of IED and non-IED time intervals. This information is crucial for estimating the IED regions for both methods correctly. However, R-SS can tolerate 10% of errors in the labeling.

Conclusion

Using the R-SS method leads to the temporal sources and their corresponding spatial patterns, each electrode lead is represented as an i * -dimensional feature vector. The number of sources of IED class, i * is estimated automatically based on Bayes probability error. By applying Pareto optimization method on the feature vector values of all the electrode leads, we obtain the electrode leads with significantly large IED activity that build the estimated IED regions.

The method is applied to the iEEG recordings of five patients suffering from epilepsy. All patients are seizure-free after resective surgery. The IED regions estimated by our method are congruent with SOZ visually inspected by an epileptologist and other automatic method [START_REF] Amini | Directed Differential Connectivity Graph of Interictal Epileptiform Discharges[END_REF] for these five patients. The method requires two sets of intervals, i.e. IED (or reference) and non-IED. The correct labeling of these intervals is crucial for obtaining the congruent results, but 10% of errors in labeling does not make a significant change in the results. The proposed method is also reliable for small number (about 50) of IED and non-IED time intervals. However, each time interval should include enough number of samples for statistically reliable estimation of correlations.

The efficiency of the method is also evaluated qualitatively using simple simulated data. A complete study on more number of simulated data for the analysis of the effect of different parameters, and their related quantitative evaluations are our first perspective. One limitation of the method is the manual labeling of IED and non-IED time intervals: automatic detection is thus a second perspective. As a third perspective, the proposed method must be applied to a larger number of simulations and patients for a more complete validation. A forth perspective could be using R-SS method as a preprocessing step for directed DCG construction for decreasing its computation time. For this purpose, we can apply R-SS to find IED related electrode leads and then calculate dDCG for these electrode leads to estimate the robust differential network between them. Therefore, for a limited number of electrode leads, the computations for dDCG would be much faster. Finally, localization of the epileptic sources by solving the inverse problem can also be considered as future works to obtain better localization for epileptic sources. Labeling the epileptic intracranial EEG automatically as IED and Non-IED intervals can be considered as an extension to this work. Two main approaches for solving the inverse problem are non-parametric (imaging) and parametric (scanning) methods [START_REF] Baillet | Electromagnetic Brain Mapping[END_REF]. A brief review of these methods can be found in [START_REF] Bolstad | Space-time event sparse penalization for magneto-electroencephalography[END_REF]. In this paper, we develop an imaging method.

Part IV EEG-fMRI Integration

The forward problem of EEG is modeled as:

X = GJ + E n (6.1)
where X is a N × T matrix that contains scalp readings, N is the number of channels and T is the number of time points, G N ×M is the gain matrix, J M ×T contains dipole magnitudes at different time instants, E n represents the noise, and M is the number of mesh vertices and, practically, N � M . Under this notation, the inverse problem consists in finding an estimate of the dipole temporal courses, which are the rows of matrix J, given the electrode positions and scalp readings X and using the gain matrix G, which can be calculated with boundary element method (BEM).

Volumetric structures derived from the tessellation procedure are used to create a realistic geometry of the head, which is crucial in the forward modeling of EEG fields. Spatial information is especially important for EEG forward modeling due to the fact that it is more strongly affected by the conductivities of the skull and the scalp than the MEG forward model. There are four numerical methods available to solve the EEG modeling problem, and the Boundary Elements Method (BEM) [START_REF] Hämäläinen | Realistic conductivity geometry model of the human head for interpretation of neuromagnetic data[END_REF] is the most commonly used when isotropy (direction independence) of the matters is assumed, so that only boundary meshes obtained by the tessellation process are required. It was shown, however, that anisotropy of the skull [START_REF] Marin | Influence of skull anisotropy for the forward and inverse problem in EEG: simulation studies using FEM on realistic head models[END_REF] and white-matter [START_REF] Wolters | Influence of tissue conductivity anisotropy on EEG/MEG field and return current computation in a realistic head model: a simulation and visualization study using high-resolution finite element modeling[END_REF] can bias EEG and MEG forward models. To solve the forward problem in the case of an anisotropic medium, the head volume is presented by a large assembly of small homogeneous tetrahedrons, and a Finite Elements Method (FEM) [START_REF] Miller | Finite element analysis of bioelectric phenomena[END_REF] is used to approximate the solution.

Another possible way is to use the Finite Difference Method (FDM) on a regular computational mesh [START_REF] Saleheen | New finite difference formulations for general inhomogeneous anisotropic bioelectric problems[END_REF]. Forward modeling of EEG signal relies on the knowledge of matter conductivities.

Common values of conductivities for different tissues can be found in the literature [START_REF] Geddes | The specific resistance of biological material-a compendium of data for the biomedical engineer and physiologist[END_REF], or can be estimated on a per-subject basis using Electrical Impedance Tomography (EIT) [START_REF] Gonçalves | In vivo measurement of the brain and skull resistivities using an EIT-based method and the combined analysis of SEF/SEP data[END_REF] or Diffusion Tensor (DT) [START_REF] Tuch | Conductivity tensor mapping of the human brain using diffusion tensor MRI[END_REF] MRI. Here, we use BEM with default values in SPM8 for simplicity.

The other methods can be considered in future works.

The EEG inverse problem is underdetermined (since N � M ) and lacks a unique solution, which is due to the fact that there are fewer observations than the unknown variables. Obtaining a physiologically feasible solution can be done through reasonable constraints. Here, we use the spatial sparsity of the brain activity as a constraint. As the brain activities are localized, this assumption is congruent with reality.

EEG is not potentially sparse. Therefore, for using sparsity, the data should be projected onto another space or a sparse function of the data should be used in the penalty term. Some basis has been used to project EEG data onto a more sparse space, like wavelet [START_REF] Strohmeier | MEG/EEG source reconstruction based on Gabor thresholding in the source space[END_REF] and principle components (uncorrelated sources with high power contribution in the power of the observations which are estimated by principal component analysis (PCA)) [START_REF] Ou | A Distributed Spatio-temporal EEG/MEG Inverse Solver[END_REF]. Experiment or task information is not used in any of these methods, so an additional statistical method is essential to find the task related regions. In this work, we use the information about the experiment to extract the sources relevant to the task of interest as the basis. Then, to achieve more sparsity in the EEG inverse problem, we project the EEG data onto its source space. Source space is found by applying the reference-based source separation (R-SS) method introduced in chapter 4. R-SS is a semi-blind source separation method which extracts the discriminative sources of one state in comparison with another. Therefore, the source space here is the task related sources. In other words, we reformulate the EEG inverse problem to achieve a sparse unknown variable which is the spatial map of sources. Therefore, sparse decomposition methods can be used to estimate the spatial map of sources. For this purpose, we use the � 2,0 regularization method which is proposed in Chapter 3. In our work, for the same task, we have both EEG and fMRI recordings, so we can use the results of fMRI as a weighting matrix in the penalty term.

Another point is that we extract a few sources from the source separation method so we estimate different spatial maps related to different sources. We should judge these spatial maps to find the final results. One way is to find the power in each point and threshold it, but here, we use Pareto optimization [START_REF] Deb | Multi-Objective Evolutionary Algorithms: Introducing Bias Among Pareto-Optimal Solutions[END_REF] which is a multicriteria optimization method able to find the dominant 1 points and which does not need any thresholding process.

The main contributions of the current work are: 1) project the forward problem of EEG onto the task related source space to have a sparse space; 2) replace � 1 norm by � 0 norm in the elastic net optimization problem; 3) use the fMRI spatial map as weight in the penalty term to solve the inverse problem of EEG; 4) exploit Pareto optimization to find the final results. 1 The concept of dominance will be defined at the end of Section 6.3.

SPM 1 h() Linear Model R-SS × Pareto WEN0 fMRI Data Spatial Map SM 1 h(SM ) X = GJ X s = A s S J s = B s S S =    s � 1 . . . s � i *    A s = G s B s B s B optimum Projection Figure 6.1: Diagram of the proposed method (WEN0). h(SM ) = 1 (1-α)+αSM
The architecture of the proposed method is shown in Figure 1. and the rest of this chapter is organized as follows. Reference-based source separation method is reviewed in Section 6.2. The proposed method is explained in Section 6.3.

Source Space Estimation

Here, we use R-SS method introduced in chapter 4 and [START_REF] Samadi | Reference-Based Source Separation Method For Identification of Brain Regions Involved in a Reference State From Intracerebral EEG[END_REF] to extract the uncorrelated temporal courses which have high power in the task state. For this purpose, we assume that the EEG observations are a linear combination of uncorrelated sources as:

X = AS (6.2)
where, A is a N × N mixing matrix, a j is its j-th column, S is a N × T matrix, containing uncorrelated sources, and s � i is its i-th row. As depicted in section 4.5, the sources are uncorrelated and we have SS � = I.

The variables of the R-SS method are as follows: for the simulated data X 1 is the mean of task trials and X 2 is the mean of control trials; for the clinical data X 1 is the mean of face trials and X 2 is the mean of scrambled trials.

Denoting i * the number of relevant estimated sources (see Section 6.2), equation (6.2) can be rewritten as:

X = [a 1 . . . a i * ] [s 1 . . . s i * ] � + [a i * +1 . . . a N ] [s i * +1 . . . s N ] �
We rename the relevant part (i.e. the part related to the i * relevant sources) of A as A s = [a 1 . . . a i * ] and the rest is renamed as A s = [a i * +1 . . . a N ]. So (6.2) can be rewritten as:

X = X s + X s = A s [s 1 . . . s i * ] � + A s [s i * +1 . . . s N ] � (6.3)
where the dimensions of A s and A s are N × i * and N × (Ni * ), respectively.

As X and J have a linear relation according to (6.1), J can also be estimated from S as:

J = BS (6.4)
where B is an M × i * unknown mixing matrix which predicts J from S.

We can partition J in two parts estimated by relevant sources and irrelevant sources, as done for X.

J = J s + J s = B s [s 1 . . . s i * ] � + B s [s i * +1 . . . s N ] � (6.5)
where the dimensions of B s and B s are M × i * and M × (Ni * ), respectively.

We put the relevant uncorrelated temporal courses, estimated with R-SS in the rows of the matrix Ŝ = [ŝ 1 • • • ŝi * ] � and call it from now on the source space.

Spatial Localization

Source Space Projection

In this step, we project the EEG forward problem onto the source space by multiplying the forward model with the matrix Ŝ. Therefore, the projected version of the forward problem onto the source space would be as:

X Ŝ� = GJ Ŝ� + n Ŝ� (6.6)
Using ( 13) and ( 15) and the fact that sources are uncorrelated and normalized (S Ŝ� ∼ I) and sources are uncorrelated with noise (n Ŝ� ∼ 0), we would have:

A s I + A ŝ0 = GB s I + GB ŝ0 + 0 (6.7)
where I is an i * × i * identity matrix and 0 is a zero matrix. Therefore, A s and B s are related by the following equation:

A s = GB s (6.8)
The i-th column of B s contains the projection of the current dipole signal of all vertices onto the i-th source. Therefore, the active vertices related to the i-th source have high value in the i-th column of B s . In other words, each column of B s represents the corresponding spatial map of the sources. Equation (6.8) shows A s , the projection of the observations onto the source space, is a linear combination of the spatial maps. As the brain activity is localized, most of the vertices are not active and have a small value (near zero) in B s , so B s is a sparse matrix.

� 2,0 Regularization

The � 2,0 regularization algorithm has been developed in Section 3.5. � 2,0 regularization is considered as a multi-objective optimization problem as follows:

minimize f(b) = [�W 2 b� 0 , �W 1 b�]
subject to a = Gb (6.9)

We use lexicography method to solve the above multi-objective optimization to optimize the � 0 regularization using the fast algorithm has been proposed in [START_REF] Babaie-Zadeh | Weighted Sparse Signal Decomposition[END_REF] and in each iteration close to its � 0 response, we find the minimum � 2 norm point.

We have used this routine to estimate b s 1 , • • • , b s i * which are the columns of B s . In our problem, a s k (an N × 1 vector) is the k-th column of A s which is the projection onto the k-th source, and the weighted matrices W 2 and W 1 are M × M diagonal matrices initialized by the inverse of the spatial map extracted from the fMRI data. For this purpose, the spatial map of fMRI is normalized by its maximum value and it is registered on the mesh vertices. It is defined as SP M f M RI , which is a M × 1 vector. The diagonal elements of W 1 and W 2 are initialized as:

diag(W 1 ) i = diag(W 2 ) i = 1 (1 -α) + αSP M f M RI i (6.10)
where α can be estimated experimentally. α = 0 removes the effect of the spatial map of the fMRI data and α = 1 constrains the results strictly to the spatial map of the fMRI. SP M f M RI provides maximum values while the EEG localization is a minimization problem (6.9) so the inverse of the fMRI spatial map is used in the weighting matrix. In other words, a low value for a given vertex implies that it is not crucial in the sparseness of the variable and the algorithm will focus on decreasing the values of the other vertices.

For i * times (number of relevant sources), we solve (6.9) for b = b s k and a = a s k , k = 1 • • • i * , which are the columns of B s and a s , respectively. For finding the activation site or the place of large coefficients in B s , we use Pareto optimization.

Pareto Optimization

B s is a M × i * matrix whose columns present the contributions of the related sources in each mesh vertex. At the first glance, we can compute the power of each vertex in the source space and use thresholding to find active vertices. In this method, the threshold affects the results.

To avoid this limitation, we use multi-objective optimization method, Pareto [START_REF] Deb | Multi-Objective Evolutionary Algorithms: Introducing Bias Among Pareto-Optimal Solutions[END_REF].

A multi-objective optimization problem, in Pareto sense, has the following form:

maximize (b s i ) � subject to (b s i ) � ∈ P ⊂ � i * (6.11)
consisting of i * objective functions that are aimed to be maximized simultaneously. (b s i ) � is the i-th row of the matrix B s which shows the contributions of the sources at the ith vertex. From the geometrical point of view, each (b s i ) � can be considered as a node in a i * -dimensional space. We classify the search space P according to the Pareto concept of non-domination [START_REF] Deb | Multi-Objective Evolutionary Algorithms: Introducing Bias Among Pareto-Optimal Solutions[END_REF]: a node is a member of non-dominated layer if either it dominates the others, or there is no other node

dominating it. Node (b s i ) � dominates node (b s k ) � , if ∀l, b s i (l) ≥ b s k (l), and ∃l * , b s i (l) > b s k (l * ) 2 .
The set of all non-dominated nodes is called non-dominated layer and defined as D(P ) ⊂ P . In the following, we explain how to estimate D(P ) using Pareto optimization algorithm [START_REF] Deb | Multi-Objective Evolutionary Algorithms: Introducing Bias Among Pareto-Optimal Solutions[END_REF]. Let us consider N i * -dimensional decision vectors, (b s i ) � , as N nodes in the search space P . The non-dominated layer is obtained using the following Pareto optimization algorithm [START_REF] Deb | Multi-Objective Evolutionary Algorithms: Introducing Bias Among Pareto-Optimal Solutions[END_REF]:

1. Initialize D(P ) with the first node (i = 1) with value of (b 1 ) � . This can be any node. 3. If i is not equal to N go to step 2.

Conclusion

In this chapter, we proposed a new method to integrate EEG and fMRI for identification of active regions of the brain with high spatial resolution from scalp EEG recordings. Using the R-SS method, the temporal sources, which are responsible for the activity of interest were estimated.

Using these temporal sources, the EEG signal was projected onto the source space. The inverse problem of EEG in the source space achieves conditions of stability and uniqueness of sparse component analysis. For sparse decomposition, � 0 was used because of its advantages to � 1 . Also, to maintain the group effects, we have used the elastic net idea by combining � 2 term with � 0 term which is called � 2,0 regularization. fMRI data is used to guide the EEG spatial pattern by using a weighting matrix in the wighted � 2,0 regularization. The inverse of the fMRI spatial map is used as the weighting matrix.

In this section, we compare the proposed method with other similar methods theoretically and in the next chapter the proposed method is investigated practically.

A similar idea has been used in [START_REF] Daunizeau | Bayesian spatio-temporal approach for EEG source reconstruction: conciliating ECD and distributed models[END_REF][START_REF] Bolstad | Space-time event sparse penalization for magneto-electroencephalography[END_REF] and evaluated on simulated data. They estimate B in equation 6.4 using auxiliary variables by Bayesian statistical inference framework and Expectation Maximization (EM) method, respectively. The regularization terms of these studies are the � 2 and the � 1 norms, respectively. [START_REF] Bolstad | Space-time event sparse penalization for magneto-electroencephalography[END_REF] seems to have better spatial resolution, while both methods are sensitive to the SNR value of the observations. Here, we use reliable sources which are related to the desired task with high probability and we also use the � 0 norm which is more accurate than the � 1 and the � 2 norms. Using WSL0 leads to stability against the noise.

In EEG source localization group MUSIC (Multiple Signal Classification) and FOCUSS (Focal Underdetermined System Solution) are somehow similar to the proposed method. MUSIC belongs to the ECD group and FOCUSS belongs to the imaging group. The imaging methods can estimate the spread sources which is their advantage in comparison with ECD methods.

MUSIC [START_REF] Mosher | EEG and MEG source localization using recursively applied (RAP) MUSIC[END_REF] is searching the best dipole localization by projecting the data onto the source space derived from SVD of the observations. In MUSIC, SVD is used to extract the uncorrelated temporal sources, which are considered as the activation time courses. As such, the correlated sources cannot be estimated. In the proposed method the uncorrelated temporal sources are considered as the basis of the activation space. Therefore, the activation time courses can be correlated. Also the number of dipoles in the ECD methods is an open problem. While in the proposed method by using R-SS method the number of temporal sources is estimated automatically.

FOCUSS [START_REF] Grave De Peralta | Distributed Source Models: Standard Solutions and New Developments[END_REF] is an iterative method which in each iteration the weight is computed using the solutions estimated in the previous step, leading to a non-linear solution. In this method a physical constraint is applied such that the currents are bounded to the brain volume and thus the radial components should go to zero when approaching the surface of the brain. This assumption leads to overcome the problem of the surface-restricted minimum norm (MN) algorithm. This weighting algorithm is based on purely mathematical operations without any physiological basis while the proposed method is based on the physiological basis and it used the experimental information to extract the solution. Another drawback of FOCUSS is its sensitivity to initialization and noise [START_REF] Wipf | A unified Bayesian framework for MEG/EEG source imaging[END_REF].

In the integration methods MSP (Multiple Sparse Priors) [START_REF] Friston | Multiple sparse priors for the m/eeg inverse problem[END_REF] and SBL (Sparse Bayesian Learning) [START_REF] Tipping | Sparse bayesian learning and the relevance vector machine[END_REF] are the most similar methods to the proposed one, as they use the spatial sparsity of the activated regions. These methods use the Bayesian framework described in section 3.2.

In these methods the active sources can be correlated and any extra information which gives information about the covariance matrix of the sources can be used e.g. the fMRI results. In these methods the information about the experiment is not used, while we have used it in the proposed method. The behavior of these methods will be compared with ours in the future chapter.

The advantages of the proposed method in comparison to the previous ones can be summarized as follows.

� Using � 0 norm leads to a more stable and accurate solution.

� The projection basis used in the proposed method directly results in the spatial map of interest.

� Using Pareto optimization helps us avoid additional statistical methods and thresholding. noisy measurements SNR= 5dB. In fMRI, SNR is defined as

SN R = 20 log( S σ N ) (7.1)
where S represents the average magnitude of the signal, and σ N stands for the standard deviation of the noise [START_REF] Krüger | Physiological Noise in Oxygenation-Sensitive Magnetic Resonance Imaging[END_REF]. The noise is a mixture of different random signals: 1) Rician system noise;

2) temporal noise of order 1 ; 3) low-frequency drift; 4) physiological noise; and 5) task-related noise [START_REF] Krüger | Physiological Noise in Oxygenation-Sensitive Magnetic Resonance Imaging[END_REF]. The background is set to the mean image of the clinical data.

From each fMRI data two sources are estimated: one is related to the task state and the other to the control state (see Figure 2). The center of the spatial sources related to the task state is chosen randomly and the sources are extended with a radius of 10 mm. The control source is located at (-10,-29,55) mm in the Montreal Neurological Institute (MNI) coordinates and is extended as a sphere with a radius of 7 mm and a fading ratio of 0.01. The task sources are located at least 6 mm far from the control source.

EEG

For each fMRI data, the related EEG data is generated using the linear forward model (6.1).

The gain matrix is calculated using the structural MRI of a real subject and the position of � Activation site. We assume 50 random sites to put a dipole related to the activation of interest (task). These random sites are the same for EEG and fMRI. We also assume a fixed site for a dipole related to the control state (like in fMRI it is located at (-10,-29,55) mm in the MNI coordinates).

� Activation time courses. For simulating the activation time courses, we extract EEG source time courses from the clinical data. To extract the task and control time courses, the R-SS method is used on the clinical data with the following variables. These time sources are shown in Figure 2.

-Task: X 1 =task trial and X 2 =control trial.

-Control: X 1 =control trial and X 2 =task trial.

Background time course also is extracted from clinical data. It can be found by extracting the common sources between the task and control trials. Common source extraction is done by canonical correlation analysis (CCA). CCA and its relation with GEVD is explained in appendix B. We use the first common source to estimate the background EEG signal. By back reconstruction, we compute the EEG background in the electrode space.

EEG Preprocessing

The steps of preprocessing are as follows:

� Epoching from 500 to +1000 ms poststimulus onset.

� Down-sample to 250 Hz (using an antialiasing low-pass filter with a cutoff frequency of approximately 100 Hz).

� Reject the epochs with the EOG exceeded 150 µV.

� Re-reference the EEG data to the average over non-bad channels.

MRI + fMRI Preprocessing

The fMRI preprocessing steps are as follows:

� Realignment. Motion correction: Adjust for movement between slices Analysis of all MRI data were performed with SPM8 (http://www.fil.ion.ucl.ac.uk/ spm). Statistical analysis was performed using the usual summary statistic approach to mixed effects modeling. General linear model (GLM) was used to model BOLD response and maximum likelihood estimation was used to estimate the parameters. A statistical parametric map (SPM) of the F-statistic was thresholded for p < 0.05 (family-wise error corrected across the wholebrain) and regions of at least 10 contiguous voxels were survived as active regions. These steps are shown in Figure 7.3. This analysis produced clusters in the fusiform and occipital face areas (FFA and OFA, respectively). The result is shown in Figure 7.4.

Results

The proposed method is evaluated using simulated and clinical data. The parameter α is equal to 0.9 unless it is specified.

Simulated Data

The efficiency of the proposed method is evaluated using computer simulations, and compared with two other methods. The first method, implemented in SPM8 (http://www.fil.ion.ucl.

ac.uk/spm), is the multiple sparse priors (MSP) [START_REF] Friston | Multiple sparse priors for the m/eeg inverse problem[END_REF]. We also use MSP with the fMRI spatial map as prior information to see how fMRI information will help solving the EEG inverse problem.

MSP belongs to the Bayesian framework and it models the source covariance as a sum of unknown basis covariances. It uses Restricted Maximum Likelihood (ReML) for optimization step. In MSP, fMRI information would be used to construct a C i as a member of C in equation 3.7.

Another method is Sparse Bayesian Learning (SBL) [START_REF] Tipping | Sparse bayesian learning and the relevance vector machine[END_REF] which belongs to Bayesian framework.

It is known as the best method in source localization [START_REF] Wipf | A unified Bayesian framework for MEG/EEG source imaging[END_REF]. It estimates more sparse and more accurate sources. Here, we use fMRI information as γ to raise the contribution of the active fMRI regions in the source covariance computed in equation 3.6.

Spatial Accuracy

To compare the performance of the methods, two performance measures are defined. Localization bias (LB) is defined as the Euclidean distance between mean site of localized activity ( μ) and the real site of simulated activity, µ.

μ = � p �j p �r p � p �j p � (7.2)
where, r p is the MNI coordinate of a vertex.

LB = �µ -μ� (7.3)
Source distribution index (SDI) is defined as the localization spread: Mean values and standard deviations of performance measures are calculated (see Figures 5 and6).

SDI = � p �j p ��r p -μ� � p �j p � (7.4)

Extra Source in EEG

For evaluating the effect of extra EEG sources which cannot be detected in the fMRI data, we have added another active region in the EEG data with the weak signal power which is not detected in fMRI data. Repeating the Monte Carlo simulation for 50 random common source localization and fixed extra source in EEG shows that the results would not change. It means that our method will miss the regions which are in EEG, but not in fMRI. For SN R EEG = 0 dB and α = 0.9 the distance between the estimated active site and real activation site would be 13 mm, and no extra region has been found in 45 experiments. α can balance the role of fMRI information in the EEG inverse problem (see [START_REF] Hill | E{EG} during {MR} imaging: differentiation of movement artifact from paroxysmal cortical activity[END_REF]. To see the effect of α, we have repeated the experiment for α = 0.8 and α = 0.7. Figure 7.8 shows LB and SDI for three different values of α for the common source. About the EEG extra source, Table 7.1 shows the percentage of its detection in the 50 EEG-fMRI sets of data. Detection happens when the proposed method finds a region less than 30 mm far from the EEG extra source. We do the same experiment for SBL, and it seems that SBL is very sensitive to the mismatches between two modalities. The result of SBL is independent of the α value. It cannot estimate the common source, but it estimates the extra source with around 27mm average localization bias.

MSP is also sensitive to the uncommom sources, but it is more stable than SBL. The result of MSP in the presence of extra EEG source is identical to the result of MSP when we don not use fMRI prior. In MSP there is no variable to adjust the effect of fMRI prior.

Extra Source in fMRI

For evaluating the effect of extra fMRI sources which cannotbe detected in the EEG data, we have added another active region in the fMRI data, but not in the EEG data. Repeating the Monte Carlo simulation shows that the results would not change. It means that our method will miss the regions which are in fMRI, but not in EEG. For SN R EEG = 0 dB the distance between the estimated active site and the actual activation site would be 18 mm, and no extra region has been found. Using the same experiment on SBL and MSP show the sensitivity of the method once more as it could not detect any of the EEG sources common and uncommon).

Clinical Data

Spatial Sources

Spatial sources for group analysis of 16 subjects extracted by MSP and the proposed method is shown in Figure 7.9 and 9, respectively. For group analysis the EEG data are concatenated in a large matrix and the proposed method is applied to it like an individual subject. The results of the group analysis of their fMRI data also is used as prior information.

Stability

To evaluate the stability of the proposed method, we compute the ROC curve using the idea of [START_REF] Strother | Evaluating fMRI preprocessing pipelines[END_REF]. For this purpose, we partition the EEG data into two distinct groups each with 8 subjects. For each group, the spatial map is extracted using the proposed method and also MSP. One group is assumed as the ground truth and the other as the test one. To evaluate true positives and false positives, the localization bias of the regions in the test image are compared with the localization bias of the regions in the ground truth one. If any region in the test data is at least 20 mm far from all regions in the ground truth, it would be false positive. Otherwise, it would be true positive. 20 mm is chosen because the best resolution of MSP in the simulation was around 20 mm and it means that the regions that are less than 20 mm far from each other are acceptable as close regions. This test is repeated 20 times for 20 different partitions. The results for the two methods are shown in Figure 7.12. As the values of the false positive ratio (FPR) are not large and also in many partitioning tests the true positive ratio (TPR) is equal to one, the ROC curve plot is not informative enough. Therefore, we show the TPR and FPR in different figures.

Concatenating the data and using one lead field matrix related to the first subject would ignore the variability of the lead field matrix for different subjects while in fact they are not the same as the electrode positions in the MNI coordination are not the same for all subjects. Therefore, we have changed the strategy of our group analysis. We have analyzed each subject's data and SPM is used to do group analysis. One-sample t-test is used and the t-statistic is thresholded for p < 0.05 and regions of at least 10 contiguous voxels survive as active regions for each group.

For stability check of the new group analysis method, we partition the subjects randomly into two non-overlapped groups each with 8 subjects. Each group is analyzed using t-statistics and the results of the two groups are compared as before. This test is repeated 20 times for 20 different partitioning. The results show that TP is one and FP is zero in all 20 partitioning. It means that the new method estimates approximately the same regions in all participants.

Discussion

Simulated Data

Figure 5 shows the mean localization bias of different methods in different SNR values. Comparing the results of MSP and SBL method when fMRI information is used or not shows that using fMRI information results in better resolution. For MSP it is observable especially in lower SNR values. The figure also shows that the proposed method finds regions closer to the activation site compared with MSP, but SBL method with fMRI information has the best localization results. [START_REF] Wipf | Solving Sparse Linear Inverse Problems: Analysis of Reweighted � 1 and � 2 Methods[END_REF] proves that SBL can be considered as the reweighted � 1 that approximates the � 0 norm (Details can be found in [START_REF]A New View of Automatic Relevance Determination[END_REF]). In the reweighted � 1 methods, highly sparse solution is produced at every iteration and so early stopping is always feasible if desired [START_REF]A New View of Automatic Relevance Determination[END_REF].

In addition, the proposed method is more robust to noise. This is because the referencebased source separation method removes the common sources between the two states and thus, it removes the background effect as much as possible. Although, SBL has better spatial resolution in complete EEG-fMRI concordance, but it is too sensitive to mismatches. This fact is caused by using the fMRI information as a hard constraint on EEG source localization. Because we have used the fMRI results as γ in equation 3.6 and it means that the algorithm should find sources with predefined covariance matrix.

The proposed algorithm is also examined for an activation site which is invisible in fMRI.

The activation regions which are missing in fMRI data may be related to the neural activity which have the short activation time or involve a few numbers of neurons that cannot produce a detectable increase of cerebral blood flow [START_REF] Babiloni | Multimodal integration of EEG, MEG and fMRI data for the solution of the neuroimage puzzle[END_REF]. Missed regions can also be due to analytical method. The result shows that the effect of missed regions in fMRI is related to the balancing variable α in [START_REF] Hill | E{EG} during {MR} imaging: differentiation of movement artifact from paroxysmal cortical activity[END_REF]. With α = 0.9 fMRI invisible sources are also missed in the EEG-fMRI result. 7.1 show that by decreasing the α value, the probability of the extra source detection would increase while the mean of LB would not change much. It means that the proposed method finds the common active regions in the EEG-fMRI data set with better resolution and by decreasing α value extra regions in the EEG data can be found. Some brain cells (stellate cells) and cortical regions (thalamus) produce closed electromagnetic field. This field cannot be detected on the scalp electrodes, although their activation requires high blood flow. Therefore, these neuronal populations present high metabolism requirements that can be detected by the fMRI technique, while at the same time they are invisible in the EEG modality [START_REF] Babiloni | Multimodal Imaging from Neuroelectromagnetic and Functional Magnetic Resonance Recordings[END_REF]. We have evaluated the proposed method under the extra source in fMRI that is invisible in EEG. More accurate EEG localization is the goal of the fMRI-informed EEG analysis. Obviously, fMRI can just help in common regions. The best method should be such that the extra sources do not affect the common sources in the fusion method. In this case, the proposed method increases the mean of the localization bias around 8 mm (from 10 mm to 18 mm), but again, it is better than other methods, which in the best case, MSP using fMRI information achieves a mean localization bias of 26 mm.

Clinical Data

Spatial sources are shown in Figure 9. This activation area is a subset of the EEG analysis individually and they are located in FFA and OFA. These clusters are in general agreement with previous studies reporting a similar contrast of the faces versus the scrambled faces [START_REF] Henson | MEG and EEG data fusion: Simultaneous localisation of face-evoked responses[END_REF][START_REF] Henson | A Parametric Empirical Bayesian framework for the EEG/MEG inverse problem: generative models for multisubject and multimodal integration[END_REF].

Study of face perception on prosopagnosia patients [START_REF] Rossion | A network of occipito-temporal face-sensitive areas besides the right middle fusiform gyrus is necessary for normal face processing[END_REF] and TMS studies on healthy participants suggest that OFA is a crucial region for face perception and represents faces prior to subsequent processing in FFA [START_REF] Nguyen | Fusing concurrent EEG-fMRI with dynamic causal modeling: Application to effective connectivity during face perception[END_REF]. The FFA, on the other hand, is known to show the greatest specificity for face perception [195,[START_REF] Grill-Spector | The fusiform face area subserves face perception, not generic within-category identification[END_REF]. It is suggested that, while OFA represents faces more by their physical features [START_REF] Liu | Perception of face parts and face configurations: an FMRI study[END_REF], the FFA is involved in more integrated analysis to represent faces more holistically [START_REF] Kanwisher | The fusiform face area: a cortical region specialized for the perception of faces[END_REF]. Furthermore, studies of structural white-matter connectivity show that functionality defined face-selective regions of the OFA and FFA are strongly interconnected anatomically [START_REF] Gschwind | White-Matter Connectivity between Face-Responsive Regions in the Human Brain[END_REF]. The right hemisphere is consistently shown to be more dominant than the left for face processing [195,[START_REF] Nguyen | Fusing concurrent EEG-fMRI with dynamic causal modeling: Application to effective connectivity during face perception[END_REF].

The most important factor here is that this activation area should be similar in all subjects.

In other words, the activation site should have little spatial variation in different subjects. To examine this feature, we have estimated the false positive ratio and true positive ratio using the algorithm proposed in [START_REF] Strother | Evaluating fMRI preprocessing pipelines[END_REF]. Figure 7.12 compares FPR and TPR of the proposed method with the MSP method. It shows that, for the proposed method, the TPR is higher than MSP and its FPR is lower than MSP. Therefore, the proposed method leads to more stable localization for different subject groups. In other words, the comparison shows the low spatial variation of the method over different subject groups. This is a very valuable factor as the proposed method also detect more focal sources in comparison with MSP.

Conclusion

The integration method is applied to simulated data and on the actual EEG-fMRI data of a face perception experiment. The results of the simulated data and its comparison with two methods: 1) MSP without fMRI, and 2) MSP with fMRI, shows that the spatial map estimated by the proposed method is much closer to the real activation sites and the results are more stable in comparison to others.

The estimated activation sites on the clinical data were located in FFA and OFA regions which are congruent with the previous studies. Also, the results are stable for different subjects, examining the ROC curve using different partitioning of the subjects.

In the simulation results, it has been shown that the missed regions in fMRI would also be missed in the proposed method. For overcoming this problem, we could investigate the following idea: after applying the proposed method to localize the common activation sites of EEG and fMRI, we could use the forward model to compute the scalp signal of common EEG and fMRI regions in the electrodes and delete it from observations. We also can remove the mesh vertices related to the common EEG and fMRI regions from the forward model and use EEG localization methods to find extra activation sites in EEG in the residual data after removing common regions and their signals. The activation sites discovered by this method would lead to the regions which may be related to the task or background, or it would be a false positive. The idea should be studied in depth to find reliable activation sites.

Chapter 8

Conclusion and Future Work

Conclusion

The human brain is a very complex system. Despite many studies about its functioning, many secrets are remained. Development of the imaging methods provides a strong tool to study the brain functionality. Imaging methods are based on different measures of the brain activity and each of them has particular advantages and disadvantages. Integration of different imaging methods is considered in recent years to come to an optimum method with the advantages of the both. In this study, we use EEG and fMRI data to estimate the sources with high spatiotemporal resolution.

Simultaneous acquisition of the EEG and fMRI and their complementary characteristic are encouraging the integration of these two imaging methods. EEG has a high temporal resolution around ms while, fMRI has high spatial resolution around mm.

In this study, we propose a semi-blind source separation method which is called referencebased source separation (R-SS). It estimates uncorrelated sources which discriminate between two states of the observations. It uses the information about timing of two different states of the observations. The method can be used in any application which has a reference signal. Here, we evaluate this method on the epileptic intracranial EEG signal. We use a realistic simulation signal and also the clinical data of 5 epileptic patients who were seizure free after a resective surgery. R-SS method is based on general eigenvalue decomposition (GEVD) so it has an analytic solution and it is fast. The evaluations show that it is stable against the additive noise and it works well on the shorter time periods in comparison with its competitor proposed in [START_REF] Amini | Directed Differential Connectivity Graph of Interictal Epileptiform Discharges[END_REF] Consequently, R-SS does better that source separation, it does source extraction, i.e. it selects the relevant source to extract, based on a very simple reference. It means that it can classify the whole set of the sources to two subsets: reference and non-reference.

Regularization of a linear equation is used when we have underdetermined linear equation.

In these equations many solutions exist and the regularization term would constrain the solution space to be able to find a unique solution which is congruent with the reality. So the regularization term should be chosen according to the nature of the problem. General constraints which are used in the regularization are the � 2 and the � 1 norms. The � 2 norm tends to find correlated sources, while the � 1 norm tends to find sparse sources. The combination of these two norms also has been considered to estimate sparse correlated sources.

The � 0 norm is also considered instead to � 1 norm, and a fast and stable algorithm has been proposed to solve it. Here, we use the combination of the � 0 and the � 2 norms which estimates sparse correlated sources which is fast and stable. We use lexicography optimization method to minimize both � 0 and � 2 norms. The method is called � 2,0 regularization. EEG-fMRI integration method proposed in this study belongs to the fMRI-informed EEG group. The result of fMRI analysis is used as a guide to analyze the EEG data. The method uses a simple assumption about the EEG and fMRI relation which is the spatial correlation between their sources. It also assumes that the spatial maps of the sources are sparse, which is based on the fact that the brain activity is localized. These assumptions are the trivial physiological facts. However invisible sources in each modality may be missed in the integration method. In fact, fMRI-infromed EEG method emphasizes on signals which are located at places proposed by fMRI and for which the EEG signal has specific properties. So new sources either in fMRI or in EEG (but not common) cannot change the results as much as common sources.

In the proposed, the R-SS method is used to estimate some uncorrelated temporal sources which can be considered as the basis of the activity sources. This basis is used to define the source space. The EEG inverse problem is projected onto the source space to have a sparse linear equation in which the projection of the scalp EEG signal on the source space is a linear combination of the spatial map of the basis of the source space. In the source space the sparse decomposition method based on � 2,0 regularization. The block diagram of the proposed is shown in Figure 1.

The proposed integration method is evaluated on simulated and actual data sets. Experimental results prove that the proposed method has superiority to a similar method, multiple sparse priors (MSP), proposed in [START_REF] Friston | Multiple sparse priors for the m/eeg inverse problem[END_REF]. First, it estimates the spatial sources closer to the real activation sites, Second, it is more stable against noise and subject variations. Third, it is more stable against invisible sources.

Another advantage of the proposed method is that it uses the Pareto optimization method, and no threshold level is defined in the method.

Perspectives

This work has been partly developed in the framework of the project Challenges for Extraction and Separation of Sources (CHESS, 2012-ERC-AdG-320684), but it only proposed partial answers to the challenges identified in this project. In fact, although EEG-fMRI integration is studied much in recent years, some challenges still exist as open problems. Here, we want to look at the problem from a more general point of view. We will consider EEG-fMRI integration as a source separation/extraction of multimodal observations. Separation/extraction of sources are wide concepts in information sciences, since sensors provide signal mixing and an essential step consists in separating/extracting useful information from useless one, the noise. Here, we want to explain three challenges in this field with the application of EEG-fMRI integration.

Multimodal Source Separation

Due to technological progresses in microelectronics and computer science, the actual world is observed through many sensor arrays, leading to multimodal and multidimensional observations. Multidimensionality is due to using a set of sensors, e.g. microphone arrays in audio processing, antenna arrays in communication, hundreds of electrodes in electro-(EEG) or magnetoencephalography (MEG). Multimodality, which implies heterogeneity, is due to the fact that the same phenomenon is measured using different physical principles and different types of sensors, e.g., brain activity can be measured with EEG and fMRI, based on electric and BOLD effects, respectively.

However, designing efficient signal processing methods for multimodal data (at the signal level and not at the decision or higher levels) still remains an open problem of very high interest as pointed out in [200,[START_REF] Khaleghi | Multisensor data fusion: A review of the state-of-the-art[END_REF]. In fact, multimodal source separation has been addressed according to ad hoc methods; especially in biomedical signal processing, A few authors considered the problem of jointly analyzing EEG/fMRI data [START_REF] Calhoun | A Feature-Based Approach To Combine Multimodal Brain Imaging Data[END_REF][START_REF] Daunizeau | Symmetrical eventrelated EEG-fMRI information fusion in a variational Bayesian framework[END_REF]. Others considered joint analysis of EEG or fMRI of a group of patients [START_REF] Calhoun | Fusion of Multisubject Hemodynamic and Event-Related Potential Data Using Independent Component Analysis[END_REF] (multimodality is then related to multiple patients). In general, the common signals correspond to the redundancy property, while specific signals are related to complementarity. Sometimes, the two modalities have the same sampling frequency, and sometimes they have not (e.g. in EEG-fMRI integration), it means that samples have to be registered, or completed if one modality undersamples the signals. One can also wonder how a conflict, e. g. related to a wrong sensor output [START_REF] Moussaoui | On the decomposition of Mars hyperspectral data by ICA and Bayesian positive source separation[END_REF], can be modeled for being detected and cancelled. But, there is no general theoretical framework for modeling the typical properties and characteristics of multimodalities.

The necessity of investigating various strategies for analyzing multimodal data is felt strongly: one can jointly analyze the whole data, or analyze separately each data set and then study similarity between estimated sources of the various modalities. Do both strategies achieve the same performance? Are there conditions for which performance of one method outperforms performance of the other? In this multimodal framework, are there new results concerning identifiability1 and separability2 ? Especially in EEG-fMRI integration, the goal of multimodality analysis is not specified uniquely, because there exist mismatches between two modalities. Figure 8.1 shows the different spaces that we are interfaced in multimodal analysis. Are The common sources the most important part? Should the common sources be estimated with better resolution, while the uncommon sources also should be estimated? Or we should define another space called useful sources. Is estimation of the useful sources the goal of multimodal analysis?

Extraction of Sources in High-or Low-dimension Data

Let us consider brain activity like a function f (x, y, z, t), where x, y, z are spatial coordinates and t is the time. Then an ideal resolution is, spatially a few 10,000 of voxels (which is achieved with MRI), and temporally a bandwidth of a few 100Hz (which is achieved with EEG). 1) By merging EEG and MRI, is it possible to achieve high resolution both spatially and temporally ? 2) Practically, EEG is observed on a few hundreds of electrodes, so that estimation of brain activity (i.e. f (x, y, z, t) with high spatial resolution is a strongly underdetermined inverse problem. How to do for making the problem to have a solution ? In this PhD, we propose a specific approach based on the spatial sparsity of fMRI sources for solving the problem. But more generally, there are many theoretical issues. Here are a few of them: Are we sure that we can improve the solution? What kind of problems can we solve? Can we predict the performance gain? What are the properties of each modality which insure performance gain? EEG has a large number of temporal samples and fMRI has a large number of spatial samples.

In the case of high-dimension mixtures, in addition to the computational load, a first problem occurs -even most of the time ignored -when processing the whole (very large) data by averaging. Some sources are relevant locally, but not for the whole data. Especially, rare events will not be correctly modeled. For a more accurate approximation, even for rare events, one has to consider local models. However, the main issue is to define the segmentation in local segments or regions? A first idea is to use the approximation of whole data as a preliminary step for determining segments or regions where the approximation is poor and must be improved. A second idea is to use priors on the source of interest for segmenting the data. Priors could be related to either particular property of the sources of interest, e.g. non-stationarity, temporal coloration, positivity, sparsity, or to similarity with some reference signals. For a high-dimension multimodal data, these new ideas can be investigated for avoiding to process the whole data, and only selecting sources of interest, common to the various modalities and/or specific to each one.

For low dimension data, one is faced to underdetermined mixtures, i.e. mixtures with more sources than observations. In such cases, source separation problem is much more complex. In fact, even if the mixing mapping is known, without extra priors, the source restoration is an illposed problem with an infinite number of solutions [START_REF] Taleb | On underdetermined source separation[END_REF]. Source sparsity is one prior which leads to a general very powerful framework [START_REF] Gribonval | A survey of Sparse Component Analysis for Blind Source Separation: principles, perspectives, and new challenges[END_REF] for solving underdetermined mixtures, but assume all the sources are sparse. One possible basic idea is to use a deflation procedure. The first source is extracted based on a spatial filter, driven by one prior. Then, after removing its contribution in the whole data, we extract the second source, based on another prior, etc. We also have to define what priors can be used and to design related criteria. However, in addition to well-known error accumulation due to the deflation process, at a given step, the spatial filter defines a particular signal subspace which can contain useful and useless sources, or two different useful sources. Moreover, at each iteration of the deflation procedure, after removing the contribution of the new extracted source, the data rank decreases of 1, which is very annoying with a small number of sensors. For preserving the data rank, one possible direction of research is to associate to source separation methods, models of the sources (useful or not) so that each signal selected by the spatial filter can be post-processed [START_REF] Sameni | A Deflation Procedure for Subspace Decomposition[END_REF] using the model for only removing one source.

In this challenge, the main theoretical issues are the following: how to select the sources among a huge number when the data is very high-dimension, or low-dimension? What are the performance bounds with respect to analysis of the whole data? How can we preserve the data rank when the number of sensors is small? What priors and criteria can be used? Can we jointly extract sources of interest in multimodal data?

EEG source localization and sparse decomposition

The efficacy of modern Bayesian techniques for quantifying uncertainty and explicitly accounting for prior assumptions make them attractive candidates for source localization. However, it is not always transparent how these methods relate, nor how they can be extended to handle more challenging problems, nor which ones should be expected to perform better in various situations relevant to MEG/EEG source imaging. Moreover, many other source localization algorithms are developed in parallel, which can be similar to one of the previous methods or it can be considered as an especial case of a more general method. [START_REF] Wipf | A unified Bayesian framework for MEG/EEG source imaging[END_REF] proposed a Bayesian framework for EEG source localization and he shows that many other methods can also be compatible with the framework. For example optimization with constraints relates to Bayesian framework and especially to MAP solutions. In chapter 3 some examples (MNE, MCE, MSP, SBL) are included. To avoid many parallel works, any new source localization methods may be compared with Bayesian framework to see if it belongs to these techniques.

In this study, we proposed a new regularization method � 2 and � 0 norm. Investigating the relation between the new method and the Bayesian framework may help us to use the efficient algorithms in this context to get more accurate estimations with higher speed. Consider two random variables x and y with zero mean and stemming from a multi-normal distribution with the total covariance matrix:

C = � C xx C xy C yx C yy � = E � � x y � � x y � T � (B.1)
Now, suppose that the goal is to find the directions of maximum data correlation. Consider the linear combinations x = x T ŵx and y = y T ŵy of the two variables, respectively. This means that the function to be maximized is: An important property of canonical correlations is that they are invariant with respect to affine transformations of x and y. An affine transformation is given by a translation of the origin followed by a linear transformation. The translation of the origin of x ory has no effect on ρ since it leaves the covariance matrix C unaffected. Invariance with respect to scalings of x and y follows directly from equation B.2.

ρ = E [xy] � E [x 2 ] E [y 2 ] = E � ŵT x xy T ŵy � � E � ŵT x xx T ŵx � E � ŵT
Another important property is that the canonical correlations are uncorrelated for different solutions, i.e.

       E[xx] = E[w T x i xx T w x j ] = w T x i C xx w x j = 0 E[yy] = E[w T y i yy T w y j ] = w T y i C yy w y j = 0 f ori � = j, E[xy] = E[w T
x i xy T w y j ] = w T x i C xy w y j = 0 (B.9) according to equation 4.4 and 4.5.

B.2 Relation to SNR

The correlation is strongly related to signal to noise ratio (SNR), which is a more commonly used measure in signal processing. This relation is used later in this thesis. Consider a signal x and two noise signals η 1 and η 2 all having zero mean and all being uncorrelated with each other.

Let S = E[x 2 ] and N i = E[η 2 i ] be the energy of the signal and the noise signals respectively. Then the correlation between a (x + η 1 ) and b(x

+ η 2 ) is ρ = E[a(x + η 1 )b(x + η 2 )] � E[a 2 (x + η 1 ) 2 ]E[b 2 (x + η 2 ) 2 ] = x 2 � (E[x 2 ] + E[η 2 1 ])(E[x 2 ] + E[η 2 2 ]) = S � (S + N 1 )(S + N 2 ) (B.10)
Note that the amplfication factors a and b do not affect the correlation or the SNR.

B.2.1 Equal noise energies

In the special case where the noise energies are equal, i.e. N 1 = N 2 = N , Eq. B.10 can be written as

ρ = S S + N
This means that the SNR can be written as

S N = ρ 1 -ρ
Here, it should be noted that the noise affects the signal twice, so this relation between SNR and correlation is perhaps not so intuitive.

B.2.2 Correlation between a signal and the corrupted signal

Another special case is when N 1 = 0 and N 2 = N . Then, the correlation between a signal and a noise-corrupted version of that signal is

ρ = S � S(S + N )
In this case, the relation between SNR and correlation is

S N = ρ 2 1 -ρ 2

Méthode

En imagerie céreébrale, on suppose que l'activité est localisée sur une grille de plusieurs milliers de dipôles, de sorte que le modèle direct conduit à l'équation :

X = GJ + n (1)
Linear Model R-SS × Pareto � 2,0 IRMf Data Spatial Map SM 1 h(SM ) X = GJ X s = A s S J s = B s S S =    s � 1 . . . s � i *    A s = G s B s

Estimation de l'espace des sources

Nous utilisons la méthode R-SS introduite dans [1] pour extraire les sources temporelles (signaux) non corrélées qui ont la plus grande similitude avec un signal de référence. La méthode R-SS est utilisée lorsque l'on peut segmenter les signaux en intervalles de temps qui ne se chevauchent pas, chacun exclusivement associés à un état parmi deux états différents. Ici, nous considérons deux états, noté C 1 et C 2 , qui correspondent à un état de référence et un état de non-référence, respectivement. L'intervalle de référence est lié au signal d'intérêt tandis que l'intervalle de nonréférence est associé au fond, au contrôle ou à l'état de repos. Cette segmentation dépend de l'expérience ou du problème auquel nous sommes confrontés. R-SS utilise le concept de l'analyse en composantes indépendantes quasi-périodiques, voisin de la Common Spatial Pattern (CSP) qui peut être résolue par décomposition en valeurs propres généralisée (GEVD) [2,3]. En outre, cette méthode utilise les valeurs propres pour estimer le nombre i * de sources pertinentes, c'està-dire similaires à la source de référence. Les enregistrements EEG peuvent être modélisés par une combinaison linéaire de sources non corrélées :

X = AS (2)
où A est une matrice N × N de mélange, dont a j est la j-ème colonne, S est une matrice N × T contenant les sources non corrélées, dont s � i est la i-ème ligne. Sans perte de généralité, les sources temporelles sont normalisées ce qui signifie que SS � = I, où I est la matrice identité N × N , et S � désigne la transposée de S. Nous utilisons la méthode R-SS pour extraire les sources de i * liés à un état de référence. Par conséquent, l'équation (2) peut être réécrite :

X = [a 1 . . . a i * ] [s 1 . . . s i * ] � + [a i * +1 . . . a N ] [s i * +1 . . . s N ] � En notant A s = [a 1 . . . a i * ] la partie pertinente de A et a s = [a i * +1 . . . a N ] la partie restante, l'équation (2) s'écrit : X = X s + X s = A s [s 1 . . . s i * ] � + A s [s i * +1 . . . s N ] � (3)
où les dimensions de

A s et A s sont N × i * et N × (N i * ), respectivement.
Comme X et J sont liées par la relation linéaire (1), J peut aussi être estimée à partir de S selon :

J = BS (4)
où B est une matrice inconnue de mélange de taille M × N qui prédit J à partir de S.

Nous pouvons partitionner J en deux parties associées aux sources pertinentes et non pertinentes, comme nous l'avons fait pour X.

J = J s + J s = B s [s 1 . . . s i * ] � + B s [s i * +1 . . . s N ] � (5)
où les dimensions de

B s et B s sont M × i * et M × (N -i * ), respectivement.
Les sources temporelles non corrélées pertinentes estimées par la R-SS constituent les lignes de la matrice Ŝ = [ŝ 1 • • • ŝi * ] � et forment ce que nous appèllerons à partir de maintenant espace source.

Localisation spatiale

Projection dans l'espace des sources Dans cette étape, nous projetons le problème EEG direct dans l'espace source en multipliant le modèle direct par la matrice Ŝ. Le problème direct projetée dans l'espace source s'écrit donc :

X Ŝ� = GJ Ŝ� + n Ŝ� (6)
En utilisant (3) et (5) et le fait que les sources normalisées ne sont pas corrélées entre elles (S Ŝ ∼ I), ni corrélées avec le bruit (n Ŝ ∼ 0), nous pouvons écrire :

A s I + A ŝ0 = GB s I + GB ŝ0 + 0 (7)
où I est une matrice identité de taille i * × i * et 0 est une matrice nulle. Par conséquent, A s et B s sont liées par l'équation suivante :

A s = GB s (8)
k = arg min b s k �a s k -Gb s k � + λ 2 �W 2 b s K � + λ 1 �W 1 b s K � 0 où �.� et �.� 1 représentent les normes � 2 et � 1 , respectivement. Ce problème d'optimisation peut être écrit comme : minimize b s k �W 1 b s k � 0 subject to �a s k -Gb s k � + λ 2 �W 2 b s k � < δ.
L'optimisation � 2,0 peut être simplifiée en utilisant le changement de variables suivant [4]:

G * = (1 + λ 2 ) -1 2 � G √ λ 2 W 2 � a * k = � a s k 0 � b * k = � 1 + λ 2 b s k
Le problème d'optimisation peut alors être écrit comme :

minimize b * k �W 1 b * k � 0 subject to �a * k -G * b * k � < δ (9)
CONTENTS Un algorithme rapide a été proposé dans [5] pour résoudre (9). Nous avons utilisé cet algorithme pour estimer b s 1 , • • • , b s i * qui sont les colonnes de B s . Dans notre problème, a s k (vecteur de taille N × 1) est la k-ème colonne de A s , c'est-à-dire la projection sur la k-ème source, et les matrices de poids W 2 de et W 1 de taille N × N sont des matrices diagonales initialisées par l'inverse de la carte spatiale extraite des données IRMf. A cet effet, les cartes spatiales de l'IRMf sont normalisées par la valeur maximale. Il est noté SP M IRM f , qui est un vecteur de taille M × 1. Les éléments diagonaux de W 1 et W 2 sont initialisés selon la relation :

diag(W 1 ) i = diag(W 2 ) i = 1 (1 -α) + αSP M i IRM f (
s i ) � for i = 1, • • • , M subject to (b s i ) � ∈ P ⊂ � i * (11)
Il se compose alors de i * fonctions objectifs qui doivent être maximisées simultanément. (b s i ) � , la i -ième ligne de la matrice B s , contient les contributions des sources sur le sommet i. D'un point de vue géométrique, chaque (b s i ) � peut être considéré comme un point dans un espace de dimension i * . Nous classons l' espace de recherche P selon le concept de non-domination de Pareto [6] : un point est un membre de la couche non-dominée si ce point domine tous les autres, ou s'il n'y a pas d'autre point qui le domine. Un point de (b

s i ) � domine le point (b s k ) � , si ∀l, b s i (l) ≥ b s k (l) et ∃l * , b s i (l) > b s k (l * ) 1 .
L'ensemble de tous les points non dominés est appelée couche non-dominée, ce que l'on note D(P ) ⊂ P . Dans la suite, nous expliquons comment estimer D(P ) en utilisant l'algorithme d'optimisation de Pareto [6]. Prenons M vecteurs de décision de dimension i * , (b s i ) � dans l'espace de recherche P . La couche non-dominée est obtenue en utilisant l'algorithme d'optimisation de Pareto suivant [6]:

IRMf

Nous utilisons le logiciel neuRosim écrit en R [7] pour la simulation d'IRMf, notamment en raison de sa rapidité pour la génération de données [8]. Une fonction double-Gamma, qui est une combinaison linéaire de deux fonctions Gamma, est utilisée pour modéliser FRH [9,10]. ion.ucl.ac.uk/spm). Nous supposons que les dipôles sont sur la matière grise corticale, où chaque dipôle est perpendiculaire à la couche corticale. Nous avons mis les décours temporels de la tâche et du contrôle dans les lignes associées de la matrice J. Nous avons

Données réelles

Les données cliniques proviennent de ftp://ftp.mrc-cbu.cam.ac.uk/personal/rik.henson/ wakemandg_hensonrn/. Un bref résumé des caractéristiques de l'acquisition et des étapes de prétraitement est présenté dans ce section. Plus de détails peuvent être trouvés dans [12,13].

Participants

Les participants sont seize jeunes adultes en bonne santé (huit femmes et huit hommes). Le protocole de l'étude est approuvée par le conseil d'examen éthique local (référence CPREC 2005,08). Chaque sujet a effectué des expériences sur deux jours : une fois en MEG + EEG et une fois en l'IRMf + IRM.

Stimuli et tâches

Le paradigme est similaire à celui utilisé précédemment sous EEG, MEG et IRMf [12,13]. Une croix de fixation centrale (présentée avec une durée aléatoire de 400 -600 ms) est suivie d'un visage normal ou visage brouillé (présenté pour une durée aléatoire de 800 à 1000 ms), suivi par un cercle central de 1700 ms.

Prétraitement des EEG

Les étapes de pré-traitement des EEG sont les suivantes:

� Sélection d'une fenêtre de -500 à +1000 ms par rapport au début de la stimulation.

� Sous-échantillonnage à 250 Hz (en utilisant un filtre passe-bas anti-repliement avec une fréquence de coupure d'environ 100 Hz).

� Rejet des portions où l'EOG a dépassé 150µV .

Traitement IRMf

Les étapes de traitement sont effectuées avec SPM8 (http://www.fil.ion.ucl.ac.uk/spm).

Le modèle linéaire général (GLM) est utilisé pour modéliser la réponse BOLD et l'estimation du maximum de vraisemblance est utilisée pour estimer les paramètres. Une carte paramétrique statistique (SPM) de la F-statistique a été seuillée pour p < 0, 05 et les régions avec au moins 10 voxels contigus sont conservées comme régions actives.

Cette analyse produit des régions actives dans les régions fusiforme et frontale occipitale (FFA et OFA, respectivement).

Résultats

La méthode S4 est évaluée en utilisant des données simulées et cliniques. Le paramètre α (de l'équation (10.)) est égal à 0,9, sauf spécification contraire. Cette valeur est proposée dans [14] et elle est utilisé par d'autres auteurs dans les méthodes de groupe en EEG informée par IRMf.

Données simulées

L'efficacité de la S4 est évaluée à l'aide de simulations informatiques, et comparée à deux autres méthodes, mises en aeuvre dans SPM8 (http://www.fil.ion.ucl.ac.uk/spm). La première méthode est MSP (Multiple Sparse Prior) (MSP) [15]. Dans la seconde méthode, nous utilisons la carte spatiale IRMf comme information préalable dans MSP pour voir comment l'information IRMf aide à résoudre le problème inverse EEG.

Précision spatiale

Pour comparer les performances des méthodes, deux mesures de performance sont définies.

L'écart de localisation (LB) est défini comme la distance euclidienne entre le site moyen de

Données cliniques

Sources spatiales

Les sources spatiales pour l'analyse de groupe de 16 sujets sont montrées à la figure 9. Pour une analyse de groupe, les données de l'EEG ont été analysées pour chaque sujet et le logiciel SPM a été utilisé pour faire l'analyse de groupe. On utilise un t-test sur un échantillon et on seuille la t-statistiqueà p < 0, 05 et on conserve les régions pour lesquelles au moins il y a au moins 10 voxels contigës.

Stabilité

Pour évaluer la stabilité de la méthode S4, nous calculons les taux de faux positifs (FP) et vrais positifs (TP) en utilisant l'idée de [18]. A cet effet, on partitionne les données de l'EEG en deux groupes sans intersection, pour chacun de 8 sujets. Pour chaque groupe, la carte spatiale est extraite à l'aide de S4. Un groupe est supposé être la vérité de terrain, l'autre est le groupe de test. Pour évaluer les vrais positifs et les faux positifs, les écarts de localisation des régions dans l'image de test sont comparés à la l'écart de la localisation des régions des images du groupe vérité de terrain. Si une région dans les données de test est au moins 20 mm plus loin que de toutes les régions de la vérité de terrain, il est considéré comme un faux positif. Sinon, c'est un vrai positif. La distance de 20 mm est choisie car la meilleure résolution MSP dans la simulation

Données réelles

Les résultats de sources spatiales sont présentés à la figure 9. Les régions d'activation sont un sous-ensemble de celles obtenues par l'analyse EEG seule (sans IRMf). Ces sources sont situées dans FFA et OFA, ce qui est en accord avec les études précédentes concernant visages normaux contre visages brouillés [12,13]. La caractéristique la plus importante ici est que ces régions d'activation devraient présenter peu de variation spatiale entre les différents participants. Pour tester cette propriété, nous avons estimé les taux de faux positifs et de vrais positifs en utilisant l'algorithme proposé dans [18]. Les résultats montrent une grande stabilité ou un faible variation spatiale de la méthode sur les différents groupes de participants.

Conclusion

Nous Les avantages de la méthode S4 par rapport aux autres sont les suivantes.

� Utilisation de la norme � 0 conduit à une solution plus stable.

� La base de projection utilisée dans S4 conduit directement à la carte spatiale d'intérêt.

Figure 1 . 1 :

 11 Figure 1.1: Samples of EEG scalp cap.

Figure 1 . 2 :

 12 Figure 1.2: Comparison of spatial and temporal resolution of some imaging methods.

Figure 1 . 3 :

 13 Figure 1.3: Comparison of spatial and temporal resolution of some imaging methods.

Figure 1 . 4 :

 14 Figure 1.4: A schematic of the main EEG-fMRI integration methods. Taken from [50].

Figure 1 . 5 :

 15 Figure 1.5: Using the EEG features in the fMRI analysis.

Figure 1 . 6 :

 16 Figure 1.6: Using the fMRI spatial map in the ECD model.

Figure 1 . 7 :

 17 Figure 1.7: Data model for mCCA and jICA [89].

Figure 2 .

 2 Figure 2.1 shows different groups in a schematic. A review of these different methods has been presented in Chapter 1. Here, we just mentioned some advantages and disadvantages of these methods. Methods in EEG-informed fMRI analysis and other model-based methods are very efficient when a realistic model of EEG and fMRI relation is available. However, by now none of the proposed methods are considered as the best one and none of them can interpret the EEG and fMRI data relations completely. In EEG-informed fMRI analysis and symmetric feature-based methods, the feature which should be chosen for the method are not fixed and the best choice is still an open question. Symmetric data driven methods are not based on the real physiological relations between the EEG and fMRI data. The assumed EEG-fMRI model is very simple.Therefore, there is no guarantee for them to work well in all conditions. In fMRI-informed EEG analysis method, the least assumptions are taken for the relation between the EEG and fMRI data. Therefore, They are more general and they are used in different applications more frequently. However, the fMRI invisible sources and EEG invisible sources are the challenges for

Figure 2 . 1 :

 21 Figure 2.1: EEG-fMRI integration categorization.

3 .

 3 Part IV: the integrating method is explained in this part including two chapters: (a) Chapter 6: The theoretical basis of the proposed method is explained in this chapter. (b) Chapter 7: Here, the proposed method is evaluated by the simulated and actual EEG-fMRI data set. 4. Chapter 8: the concluding remarks, perspectives, and the list of related publications are brought in this chapter.

E

  n represents the noise, and M is the number of mesh vertices and, practically, N � M . Source-localization is the extraction of the spatial information of the brain activation from the EEG data. It would be the result of the inverse solution of the linear forward equation (3.1).

a. 9 ) 3 . 4 �W 1

 9341 linear model a = Gb, where a and b are N × 1 and M × 1 vectors, is as follows: b = arg min b �a -Gb� + λ 2 �W 2 b� + λ 1 |W 1 b| (3Weighted Smooth � 0 (WSL0) Weighted � 0 regularization problem for a linear model is formulated as follows: minimize b b� 0 subject to a = Gb (3.10)

Figure 3 . 1 :

 31 Figure 3.1: WSL0 algorithm.

Figure 3 . 2 :

 32 Figure 3.2: � 2,0 regularization algorithm.

First, G 1

 1 created a system Gb = a of n = 40 equations in m = 100 unknowns having two different sparse solutions b 1 and b 2 with �b 1 � 0 = 13 and �b 2 � 0 = 13 but the non-zero points are in two regions. For b 1 , non-zero points are chosen from a normal curve which has its maximum value in the middle non-zero point, and it is randomly placed in the first half of the time points. For b 2 , non-zero points are chosen randomly in the second half of the time points. This was done as follows: of size n = 40 by m/2 = 50 was randomly created with entries drawn independently from a standardized Gaussian distribution. Then b � 1 of length m/2 = 50 was created to have 13 connected non-zero entries (whose center locations were chosen randomly), and x = G 1 b � 1 was calculated. Similarly, G 2 of size n = 40 by m/2 -1 = 49, and b � 2 of length m/2 -1 = 49 with

Figure 3 . 1 :

 31 Figure 3.1: Source signal sample.

Figure 3 .

 3 Figure 3.1 shows a signal sample for b 1 and b 2 . Then, we chose the weighting vector w such that it is equal to α for non-zero elements of b 1 and it is equal to 1α for non-zero elements of b 2 , where 0 < α < 1. So, for α < 0.5, b 1 the first source should be estimated and for α > 0.5 the second source should be estimated.

Figure 3 .

 3 2 depicts the percentages of estimating b 1 , b 2 and failure rates versus α.

Figure 3 .

 3 Figure 3.2 shows that WSL0 estimates b 1 for α < 0.4 and estimates b 2 for α > 0.4, while � 2,0 estimates b 1 for for α < 0.5 and estimates b 2 for α > 0.8. It shows that � 2,0 has the tendency to estimate b 1 which is a connected source, while WSL0 has the tendency to estimate b 2 which is a spread source. The results of weighted elastic net (WEN) are also shown for λ 2 = 0.025 and λ 1 = 0.05. The results show that WEN tends to estimate the connected source, but its accuracy is less than the � 2,0 regularization method.

Figure 3 . 2 :

 32 Figure 3.2: The percentages of b 1 and b 2 estimates, and the percentages of failure versus α.

Figure 4

 4

  .1, three multi-lead depth electrodes (A, B, and C) are considered, and positioned parallel to each other. We consider eight brain sources: two epileptic sources are placed (i) between electrodes A and B (e 1 ), and (ii) close to electrode C (e 2 ), and six non-epileptic sources are randomly placed in the volume around the electrodes. The locations of the electrodes and sources are kept constant in simulations, and for simplicity are restricted to belong to a 2-D plane. The electrical activity of each brain source is represented by a current dipole [162]. The orientations of the six non-epileptic source dipoles are randomly chosen in the 2-D plane and kept constant in simulations, while the orientation of two epileptic source dipoles are assumed as one of the three possible orientations: a tangential orientation (i.e. along the electrode axis), a radial orientation (i.e. orthogonal to the electrode axis), and a 'mixed' orientation (i.e. with a 45 • -angle with the electrode axis). Three simulations (see Fig. 5.1) are performed where the orientations of two epileptic dipoles are assumed as (1) both orthogonal (D 0 ), (2) both tangential (D 1 ), and (3) both mixed (D 2 ). The time-varying magnitudes of the source dipole moments are assumed to represent the global neuronal activity in the source regions. They are obtained from a neural mass model (Modified version of the model proposed in [167]

Figure 5 . 1 :

 51 Figure 5.1: Schematic diagram of the simulated data. Five electrodes (A, B, C, D and E), six non-epileptic and two epileptic (e 1 and e 2 ) sources in three different orientations (D 0 , D 1 , D 2 ) are demonstrated. The squares indicate the electrode leads. Each electrode consists of 10 equally spaced recording contacts (3.5-mm inter-contact spacing). The thick and thin frames demonstrate the first and second layer nodes of Pareto, respectively.

Figure 5 . 2 :

 52 Figure 5.2: Schematic diagram of the simulated data. Three electrodes (A, B, and C), six nonepileptic and two epileptic (e 1 and e 2 ) sources in three different orientations (D 0 , D 1 , D 2 ) are demonstrated. The squares indicate the electrode leads. Each electrode consists of 10 equally spaced recording contacts (3.5-mm inter-contact spacing).

Figure 5 . 3 :

 53 Figure 5.3: Implantation scheme of the iEEG electrodes (a) Sagittal view; (b) Coronal view.Table 5.3: Type of epilepsy of the patients. L/R:left/right; T:temporal; midInsG: middle short gyrus of the insula. patients focal epilepsy P1 LT P2 LT P3 LT P4 RT P5 RmidInsG

3 .

 3 dis: the average of minimum distances (mm) between IED and vSOZ electrode leads. The smaller this value, on average the closer the set of IED regions to vSOZ. 4. ovp: the average percentage of the number of IED electrode leads which are in the neighborhood (< 1.5cm) of at least one of the vSOZ electrode leads. ovp = 100% shows that we obtained at least one IED electrode lead in proximity of each of vSOZ electrode leads.

. 6 :

 6 Quantitative comparison between �IED regions based on dDCG and IED regions based on R-SS. The optimum value in each row for each measure is in bold. Precision Sensitivity dis (mm) % ovp % ovp2 R-SS �IED R-SS �IED R-SS �IED R-SS �IED R-SS �IED P 1 100 100 100 80 1.6 8 100 67 100 70 P 2 100 100 25 25 0 0 100 100 85 85 P 3 100 100 100 67 3.2 4 100 100 100 100 P 4 100 100 80 80 2 11.5 100 67 79 79 P 5 100 100 100 100 5.3 5.3 100 100 100 100 mean 100 100 81 70 2.4 5.8 100 87 93 87

Figure 5 . 4 :

 54 Figure 5.4: Mean and standard deviation of F P e (solid) and F N e (dashed) for different percent of IED time intervals (left), and different error percentages of IED labeling (right) for five patients. From top to bottom, P 1 to P 5 .

  Electroencephalography (EEG) is a non-invasive technique that records scalp electrical activity generated by brain structures. The scalp electric potentials are directly proportional to neural activity currents[15]. Moreover, EEG reflects neural activity with high temporal resolution around millisecond. The EEG inverse problem looks for the localization and dynamics of cerebral activity sources which greatly helps to figure out the aspects of cognitive processes.

2 .

 2 Choose a new node (i = i + 1): (a) If any node in D(P ) dominates node i go to step 3. (b) Else add node i to D(P ) and remove any nodes of D(P ) that node i dominates.

  the electrodes in clinical data acquisition. The time courses of the sources are generated as an additive source problem in which two main sources are distinguished, namely (1) the activation caused by an experimental design, and (2) the background. The experiment is defined by two states, task and control trials. The details are as follows.

Figure 7 . 1 :

 71 Figure 7.1: Activation time courses.

Figure 7 . 2 :

 72 Figure 7.2: Stimuli and Tasks.

Figure 7 . 3 :

 73 Figure 7.3: fMRI processing steps.

Figure 7 . 4 :

 74 Figure 7.4: fMRI spatial map.

Figure 7 . 5 :

 75 Figure 7.5: Sample activation sites in simulated data (the figure is generated with xjView toolbox (http://www.alivelearn.net/xjview)).

Figure 7 . 6 :

 76 Figure 7.6: Mean of localization bias over 50 random activation sites for five methods (MSP, MSP+fMRI, SBL, SBL+fMRI, WEN0) is compared.

Figure 7 . 7 :

 77 Figure 7.7: Mean of source distribution index over 50 random activation sites for five methods (MSP, MSP+fMRI SBL, SBL+fMRI, WEN0) is compared.

( a )

 a Localization bias for three α values. (b) Source distribution index for three α values.

Figure 7 . 8 :

 78 Figure 7.8: Effect of α value on the proposed method in existence of an EEG extra source.

Figure 7 . 9 :

 79 Figure 7.9: Activation sites in the clinical data of 16 subjects by MSP (the figure is generated with xjView toolbox, http://www.alivelearn.net/xjview).

Figure 7 . 10 :

 710 Figure 7.10: Activation sites in the clinical data of 16 subjects by the proposed method (the figure is generated with xjView toolbox, http://www.alivelearn.net/xjview).

Figure 7 . 11 :

 711 Figure 7.11: Zoom of the activation sites in the clinical data of 16 subjects by the proposed method (the figure is generated with xjView toolbox, http://www.alivelearn.net/xjview).

Figure 7 . 12 :

 712 Figure 7.12: Histogram of TPR and FPR in 20 partiotioning test.

Figure 7 .

 7 Figure 7.8 and Table7.1 show that by decreasing the α value, the probability of the extra

Figure 8 . 1 :

 81 Figure 8.1: The concept of multimodality.

8. 3

 3 List of publications � Samareh Samadi, Hamid Soltanian-Zadeh, Symmetric analysis of EEG and fMRI using nonlinear models: a Simulation Study, SPARS'09, Saint-Malo, France. � Samareh Samadi, Ladan Amini, Hamid Soltanian-Zadeh, Christian Jutten, Identification False positive and false negative in terms of electrode leads Appendix B Canonical Correlation Analysis B.1 Definitions and Properties

(B. 4 )λ

 4 Setting the derivatives to zero gives the equation system� C xy ŵy = ρλ x C xx ŵxC yx ŵx = ρλ x C yy ŵy (x is the ratio between the standard deviation of y and the standard deviation of x and vice versa. The λs can be interpreted as scaling factors between the linear combinations. Rewritingequation system B.1 gives � C -1 xx C xy C -1 yy C yx ŵx = ρ 2 ŵx C -1 yy C yx C -1 xx C xy ŵy = ρ 2 ŵy (B.7)Hence, ŵx and ŵy are found as the eigenvectors of the matricesC -1 xx C xy C -1 yy C yx and C -1 yy C yx C -1xx C xy respectively. The corresponding eigenvalues ρ 2 are the squared canonical correlations. The eigenvectors corresponding to the largest eigenvalue ρ 2 1 are the vectors ŵx1 and ŵy1 that maximize the correlation between the canonical variates x 1 = x T ŵx 1 and y 1 = y T ŵy1 . written as� C xy ŵy = r µx µy C xx ŵx C yx ŵx = rµy µx C yy ŵy (B.8) which is recognized as equation for ρλ x = r µx µy and ρλ y = r µy µx . Solving for w x and w y in equation B.8, gives equation B.7 with r 2 = ρ 2 . This shows that the equations for the canonical correlations are obtained as the result of maximizing the Rayleigh quotient.

  (EEG) est une technique d'imagerie cérébrale non invasive importante, capable d'enregistrer l'activité neuronale avec une grande résolution temporelle (ms), mais avec une résolution spatiale faible. Le problème inverse en EEG est un problème difficile, fortement sous-déterminé : des contraintes ou des a priori sont nécessaires pour aboutir à une solution unique. Récemment, l'intégration de signaux EEG et d'imagerie par résonance magnétique fonctionnelle (IRMf) a été largement considérée. Les données EEG et IRMf relatives à une tâche donnée, reflètent les activités neuronales des mêmes régions. Nous pouvons donc supposer qu'il existe des cartes spatiales communes entre données EEG et IRMf. En conséquence, résoudre le problème inverse en EEG afin de trouver les cartes spatiales des sources EEG congruentes avec celles obtenues par l'analyse de signaux IRMf semble être une démarche réaliste. Le grand défi reste de déterminer la relation entre l'activité neuronale électrique (EEG) et l'activité hémodynamique (IRMf), qui n'est pas parfaitement connue à ce jour. La plupart des études actuelles reposent sur un modèle neurovasculaire simpliste par rapport à la réalité. Dans ce travail, nous utilisons des a priori et des faits simples et généraux, qui ne dépendent pas des données ou de l'expérience et sont toujours valides, comme contraintes pour résoudre le problème inverse en EEG. Ainsi, nous résolvons le problème inverse en EEG en estimant les sources spatiales parcimonieuses, qui présentent la plus forte corrélation avec les cartes spatiales obtenues par IRMf sur la même tâche. Pour trouver la représentation parcimonieuse du signal EEG, relative à une tâche donnée, on utilise une méthode (semi-aveugle) de séparation de sources avec référence (R-SS), qui extrait les sources dont la puissance est la plus corrélée à la tâche (référence). Cette méthode a été validée sur des simulations réalistes et sur des données réelles d'EEG intracrânienne (iEEG) de patients épileptiques. Cette représentation du signal EEG dans l'espace des sources liées à la tâche, qui sont en petit nombre, est ainsi parcimonieuse. En recherchant les fonctions d'activation d'IRMf similaires à ces sources, on déduit les cartes spatiales d'IRMf très précises de la tâche. Ces cartes fournissent une matrice de poids, qui impose que les voxels activés en IRMf doivent être plus importants que les autres voxels dans la résolution du problème inverse en EEG. Nous avons d'abord validé cette méthode sur des données simulées, puis sur des données réelles relatives à une expérience de reconnaissance de visages. Les résultats montrent en particulier que cette méthode est très robuste par rapport au bruit et à la variabilité inter-sujets. CONTENTS La méthode proposée est expliquée dans la section 2. Les données de validation sont présentées à la section 3. Les résultats sont portés à la section 4. Remarques, Discussion et conclusion occupent les sections 5 et 6, respectivement.

Figure 1 :

 1 Figure 1: Schéma de la méthode S4 proposée. � 2,0 est la méthode de régularisation, et h(SM ) = (1α) + αSM

  10) où α peut être estimée expérimentalement. La valeur α = 0 supprime l'effet de pondération de la carte spatiale des données IRMf et α = 1 contraint les résultats strictement à la carte spatiale fournie par l' IRMf. Puisque SP M IRM f fournit les valeurs maximales tandis que la localisation EEG est un problème de minimisation (9), l'inverse de la carte spatiale IRMf est utilisé comme matrice de pondération. En d'autres termes, une valeur faible pour un voxel donné implique qu'il n'est pas essentiel pour la parcimonie de la variable et l'algorithme mettra l'accent sur la diminution des valeurs des autres voxels. Nous résolvons (9) i * fois (c'est-à-dire pour chaque source pertinente extraite par R-SS), pourb = b s k et a = a s k , k = 1 • • • i * , quisont les colonnes de B s et a s , respectivement. Pour trouver la région d'activation ou le position des grands coefficients dans B s , nous utilisons la méthode d'optimisation de Pareto. Optimisation de Pareto B s est une matrice de taille M × i * dont les colonnes présentent les contributions des sources liées à chaque sommet de maillage. Toutes les i * sources sont impliquées dans l'activation neuronale en relation avec la référence. Par conséquent, elles sont toutes importantes, mais avec des niveaux différents : la contribution de la première source est supérieure à la seconde, et ainsi de suite. Par conséquent, si nous utilisons la somme ou la somme des carrés, les contributions des sources faibles seraient cachées par la première source. Pour surmonter ce problème et être capable de prendre en compte la contribution de chaque source, nous utilisons alors la méthode d'optimisation multi-objectif, appelée méthode de Pareto [6]. L'optimisation de Pareto trouve une solution optimale sans seuil, ce qui évite le choix de ce paramètre. Un problème d'optimisation multi-objectif, au sens de Pareto, a la forme suivante : maximize (b

1 .

 1 Initialiser D(P ) avec le premier point (i = 1) avec une valeur de (b 1 ) � . Cela peut être n'importe quel point. 2. Choisissez un nouveau point (i = i + 1): (a) Si un point dans D(P ) domine le point i passez à l'étape 3. (b) Sinon ajouter le point i à D(P ) et supprimer tous les points de D(P ) que le point i domine.3. Si i n'est pas égal à M aller à l'étape 2.

  Une modélisation de Monte Carlo (échantillonnage sur 50 emplacements aléatoires de sources) est utilisée pour comparer les différentes méthodes avec la méthode proposée. Pour chaque site d'activation deux ensembles de données EEG et IRMf sont créées. 50 sommets aléatoires sont sélectionnés en tant que centres de leurs régions actives. Les sources sont supposées assez diffuses, ce qui est modélisé par un filtre gaussien de lissage. Ici, nous simulons les données EEG et IRMf après l'étape de prétraitement. Cela signifie que, par exemple en IRMf, le décalage temporel inter-coupe n'est pas simulé et en EEG, la moyenne des essais après prétraitement est simulée. L'expérience a deux états : la tâche et le contrôle. Le début des états de tâches et de contrôle est défini comme le début des stimulations (données cliniques) avec le visage et les états brouillés, respectivement. Tous les paramètres des méthodes d'imagerie sont choisis comme dans le cas de données cliniques : positions des électrodes, fréquence d'échantillonnage en EEG et TR, dimensions de l'image dans l'IRMf, etc.

CONTENTS��

  figures 3 et 4. � Décours temporels sur les électrodes. Pour calculer les décours temporels sur les électrodes, nous utilisons le modèle linéaire direct (1). Pour estimer la matrice de gain, on utilise l' IRM structurelle d'un objet réel avec ses repères et des positions des électrodes. Pour ce faire, nous utilisons BEM avec 8196 naeuds développé dans SPM8 (http://www.fil.

Figure 2 :Figure 3 :Figure 4 :

 234 Figure 2: Sites d'activation pour les données simulées (la figure est générée avec la boîte à outils xjView (http://www.alivelearn.net/xjview)).

  est utilisée pour modéliser le problème direct. Les images IRM de chaque participant sont segmentées et spatialement normalisées selon un modèle du cerveau INM T1 dans l'espace Talairach. L'inverse de la transformation de normalisation est ensuite utilisée pour déformer un maillage cortical d'un cerveau de modèle de l'espace de l'INM vers l'espace IRM de chaque participant. Le maillage possède 8196 sommets (4098 par hémisphère) avec un espacement moyen d'environ 5 mm entre sommets. La matrice en champ proche est ensuite calculée pour les dipôles (orientés perpendiculairement à la maille) en tous les points du maillage cortical. Tout est fait avec le logiciel SPM.

Figure 5 :p �b p �r p � p �b p � ( 14 )

 514 Figure 5: Ecart de localisation pour les trois méthodes (MSP, MSP + IRMf, S4)

Figure 6 :Figure 7 :

 67 Figure 6: Indice de diffusion de sources pour les trois méthodes (MSP, MSP + IRMf, S4)

Figure 8 :

 8 Figure 8: Effet de la valeur de α sur SDI dans le cas d'une source supplémentaire EEG.

Figure 9 :

 9 Figure 9: Sites d'activation pour les données cliniques de 16 sujets (la figure est généré avec la boîte à outils xjView, http://www.alivelearn.net/xjview). Les voxels sombres indiquent un niveau d'activation plus élevé.

  avons proposé une nouvelle méthode d'intégration de données EEG et IRMf pour l'identification des régions actives du cerveau, avec haute résolution spatiale, à partir d'enregistrements EEG de scalp. En utilisant la GEVD, nous estimons le filtre spatial qui fournit des sources temporelles classées en fonction de leur similitude avec un signal de référence. En utilisant les sources temporelles ainsi extraites, le signal EEG a été projeté dans l'espace des sources pertinentes. La résolution du problème inverse de EEG dans cet espace des sources fournit des solutions stables et uniques sous la contrainte de parcimonie. Pour la décomposition parcimonieuse, nous avons utilisées � 0 en raison de ses avantages par rapport à � 1 . En outre, afin de maintenir les effets de groupe, nous avons utilisé l'idée du filet élastique (Elastic net) en combinant les normes � 2 et � 0 . Les données IRMf sont utilisées pour guider le problème inverse en EEG : l'inverse de la carte spatiale IRMf est utilisé comme matrice de pondération dans l'optimisation � 2,0 . La méthode a été appliquée sur des données simulées et sur des données réelles EEG et IRMf d'une expérience de perception de visages. Les résultats sur données simulées sont comparés avec deux autres méthodes : 1) MSP sans IRMf, et 2) MSP avec l'IRMf. Ils montrent que la carte spatiale estimée par S4 est beaucoup plus proche des sites d'activation réels et les résultats sont plus stables par rapport à autres. Les sites d'activation estimés sur les données cliniques sont situées dans les régions de la FFA et OFA qui sont en accord avec des études précédentes. En outre, les résultats sont stables pour les différents participants, comme le montre l'étude des FP et TP réalisée utilisant différents partitionnements des données.

Table 5 . 1 :

 51 IED electrodes for different SIR values in 2D electrodes.

Table 5 . 2 :

 52 IED electrodes 3D implantation for SIR=5dB.

	Electrode lead number
	D0	1st Layer 2nd Layer E6 D0 D1 D2 D3 D4 E5
	D1	1st Layer 2nd Layer D3 E3 D1 D2 E4 E5 E6 E7
	D2	1st Layer 2nd Layer A1 D2 E5 A0 C9

Table 5

 5

.

6

. In this Table, the IED regions estimated by

Table 5 . 4 :

 54 Parameters of the patients' iEEG. N : number of the bipolar channels; T : length of the original (non segmented) iEEG signal (minute); L l : number of IED or non-IED time intervals.

		P1 P2 P3 P4 P5 mean
	N	104 105 111 109 100 106
	T (minutes) 61	56	42	90	66	55.44
	L 1	298 614 223 160 223 304
	L 2	143 200 195 183 148 174

Table 5 . 5 :

 55 Comparison between the visually inspected SOZ (vSOZ) by the epileptologist, the leading IED regions (�IED regions) estimated based on dDCG and the IED regions estimated by R-SS method.

	P 1 antHC	postHC pHcG amyg mTP
	vSOZ ×

Table 7 . 1 :

 71 Detection percentage of the extra EEG source.

	α	0.7 0.8 0.9
	Detection (%) 54 20 10

  Also in this case, since ρ changes sign if w x or w y is rotated 180 , it is sufficient to find the positive values. This function cannot be written as a Rayleigh quotient. But also in this case, it can be shown that the critical points of this function coincide with the critical points of a Rayleigh quotient with proper choices of R 1 and R 2 . The partial derivatives of ρ with respect to w x and w y are:Proof. The partial derivative ρ with respect to w x is C xy w y (w T x C xx w x w T y C y yw y ) -1 2 C xx w x w T y C y yw y w T x C xx w x w T y C y yw y

	∂ρ ∂w x	=	(w T x C xx w x w T y C y yw y ) w T x C xx w x w T y C y yw y 1 2 C xy w y x = (w T -w T x C xx w x w T y C y yw y ) -1 2	� �	C xy w y -	w T x C xy w y w T x C xx w x	� � C xx w x
				= �w x � -1 ( ŵT x C xx ŵx	ŵT y C y y ŵy ) -1 2 �	C xy ŵy -	ŵT x C xy ŵy ŵT x C xx ŵx � C xx ŵx
							=	a �w x �	C xy ŵy -	ŵT x C xy ŵy ŵT x C xx ŵx	C xx ŵx	, a ≥ 0.
							y yy T ŵy	� =	�	ŵT x C xx ŵx ŵT x C xy ŵy ŵT y C yy ŵy	(B.2)
			          	∂ρ ∂w x ∂ρ ∂w y	= =	a �w x � a �w y �	� � C xy ŵy -C yx ŵx -	ŵT x C xy ŵy ŵT x C xx ŵx ŵT ŵT y C yy ŵy y C yx ŵx	C xx ŵx C yy ŵy	� �	(B.3)
	where a is a positive scalar.

  où X est une matrice N × T qui contient les mesures sur le scalp, N est le nombre de canaux (capteurs) et T est le nombre de points dans le temps, G N ×M est la matrice de gain, J M ×T

	pertinentes (en petit nombre) qui sont obtenues par application de la méthode de séparation de
	sources avec référence (R-SS) proposée dans [1]. R-SS est une méthode de séparation de sources
	semi-aveugle qui extrait les sources caractéristiques d'un état par rapport à un autre. R-SS est
	utilisée pour extraire les sources temporelles correspondant à la tâche d'intérêt. Les sources non
	corrélées séparées par R-SS peuvent être considérées comme constituant une base de l'activité
	neuronale associée à la tâche d'intérêt. Par projection dans l'espace source, nous reformulons le
	problème inverse en EEG pour obtenir une représentation parcimonieuse qui est la carte spatiale
	des sources.
	La carte spatiale extraite de l'analyse IRMf de groupe est utilisée comme matrice de pondéra-

contient les grandeurs dipolaires à différents instants, n représente le bruit, et M est le nombre de voxels avec N � M . Le problème inverse en EEG est fortement sous-déterminé (car N � M ) et ne conduit pas à une solution unique, car il y a moins d' observations que les variables inconnues (sources). Par conséquent, des a priori supplémentaires (contraintes) sont nécessaires pour régulariser le problème d'optimisation associé et trouver une solution unique. Ici, nous utilisons la parcimonie spatiale de l'activité cérébrale comme contrainte de régularisation. On suppose qu'à chaque instant, seules quelques régions du cerveau sont actives et suffisantes pour expliquer les mesures. Cette hypothèse est en accord avec la réalité. L'architecture de la méthode proposée est illustrée à la figure 1. Dans la suite, cette méthode fondée sur la parcimonie spatiale dans l'espace des sources, est appelée S4 pour Spatial Sparsity in the Source Space.

Les signaux EEG ne sont pas essentiellement parcimonieux, de sorte que pour exploiter la parcimonie, les données doit être projetées dans un autre espace ou bien une fonction parcimonieuse des données doit être utilisée dans le terme de pénalité. Pour obtenir plus de parcimonie dans le problème inverse EEG, nous projetons les données EEG dans l'espace des sources tion de l'algorithme de dé composition parcimonieuse. La méthode de Pareto est une méthode d'optimisation multi-objectif utilisée pour trouver la carte spatiale optimale, i.e. correspondant aux différentes sources pertinentes, à partir des différentes cartes spatiales fournies par l'IRMf. SPM 1 h()

  où la i-ème colonne de B s contient la projection du signal de dipôle de courant de tous les voxels sur la i-ème source. Par conséquent, les voxels actifs liés à i-ème source ont un grand coefficient dans la i-ème colonne de B s . En d'autres termes, chaque colonne de B s représente la carte spatiale correspondante des sources. L'équation (8) montre que A s , la projection des observations sur l'espace source, est une combinaison linéaire des cartes spatiales. Comme l'activité cérébrale est localisée, la plupart des voxels ne sont pas actifs et ont un coefficient faible (proche de zéro) dans B s , donc B s est une matrice creuse.Pour produire un modèle parcimonieux avec une bonne précision, tout en favorisant un regroupement spatial des voxels actifs, en combinant la norme � 2 avec la norme � 0 , nous avons proposé une nouvelle méthode, appelée optimisation � 2,0 .

	Optimisation � 2,0

L'optimisation � 2,0 pour le modèle linéaire a s k = Gb s k (b s k est un vecteur colonne de B s k de taille M × 1 et a s k est un vecteur colonne de A s k de taille N × 1) revient à calculer : bs

Table 1 :

 1 Pourcentage de détection de la source EEG supplémentaire.

	α	0.7 0.8 0.9
	Détection (%) 54 20 10

Please note that in the two inequalities, one of them is a strict inequality.

Here, signal relates to epileptic sources, and interference to non-epileptic (background) sources. The SIR is then calculated for the significant electrodes, A0 and C9, by computing the power ratio (in dB) of the signals related to the epileptic sources (e1 and e2) over all the signals related to the non-epileptic sources.

Please note that in the two inequalities, one of them is a strict inequality.

Identifiability means the ability to estimate the unknown mixing mapping or its parameters.

Separability means the ability to estimated -from the observed mixtures -signals all depending on a unique source

Noter que dans les deux inégalités, la seconde est une inégalité stricte.

Acknowledgments

I praise God for all His blessings. Throughout my PhD, I have received tremendous support from Professor Hamid SOLTANIAN-ZADEH, my supervisor. I have been very privileged to have him as an adviser who is not only enthusiastic and professional, but also friendly.

There is no doubt that I was lucky to have Professor Christian JUTTEN as my supervisor. He always listen to me attentively and provide invaluable comments. His distinctive and timely input often enabled me to gain a new perspective on the various problems I encountered. He made me full of hope for doing such a great work with all my constraints. I think of him as Robin Hood of my life. He was not only a great scientist with brilliant ideas but also he was so kind and supportive friend for me and my family.

I sincerely thank my committee members: the president of the jury Prof. Habib BENALI for evaluating our work and accepting the problem of video conference, the reviewers Prof. Fabrice WENDLING, who also supports our study with simulations and Mohammad Bagher SHAM-SOLLAHI for reading and evaluation of my thesis very carefully and their invaluable comments, and the examiners Dr. Gholam-Ali HOSSEIN-ZADEH for teaching me the alphabets of fMRI and also for his insightful feedbacks. I would like to thank Dr. Massoud BABAIE-ZADEH for teaching me the alphabets of research and his research results on sparse decomposition helped me much. I would like to thank Prof. Olivier DAVID and Prof. Monica BACIU for supporting us with simulations and actual data. I thank my husband Dr. Farhad MOOSAKHANI who supports me from deep in his heart. I love him forever and I would never forget his helps and patience in my PhD duration. Without my family and friends, I would have never completed this work. I owe them all a debt of gratitude for their never ending support. Mostly, I want to thank my father and my brothers, Behzad and Amir and their wives for their encouragement and humor. v vi Particular thanks go to Dr. Ladan AMINI for her great support. I also want to thank her family for the warm welcomes I received whenever I was invited to her home. I thank Meri and his husband Reza for all their helps. They were so friendly that I think of them as a part of my family. I also gratefully acknowledge Delphine COSANDIER-RIM ÉL É, Marcela PERRONE, Gaetan YVERT for providing the data, interpretations, and their helpful assistance in this study. I would like to thank all my friends who helped me more than I expect them. They are Mohammad NIKNAZAR, Hamid KHOSRAVI, Pouya ZAMANI, Jonas CHATEL-GOLDMAN, Ronald PHLYPO, Mr. YAZDANPOUR, Mr. HOSSEINI and his wife, Mitra.

I thank all of my friends in GIPSA Lab for the enjoyable memories and also unpleasant memories that taught me to be kind with others. Finally, I thank Hermin and Helia, my mother and my daughter, they were with me in all my difficult moments and encouraged me to do my best. I hope I can create the moments full of joy for them in future. CONTENTS � L'optimisation de Pareto nous aide à éviter les méthodes statistiques supplémentaires et le choix d'un seuil. L'inconvénient de S4 est que le résultat final dépend de la première étape de séparation de sources. Cependant, nous n'avons pas eu de fausse découverte dans les expériences de simulation. Dans les travaux futurs, S4 peut être appliquée aux données simultanées EEG et IRMf. Définir différentes matrices de pondération pour les normes � 2 et � 0 peut être étudiée. La valeur optimisée de α, qui contrôle la contribution de la carte spatiale IRMf dans la matrice de pondération, pourrait être définie en utilisant les données. Dans les résultats de la simulation, il a été montré que les régions manquées dans IRMf sont également omises dans S4. Pour surmonter ce problème, nous pourrions étudier l'idée suivante : après application de S4 pour localiser les sites d'activation communes de l'EEG et IRMf, nous pourrions utiliser le modèle EEG direct pour calculer le signal EEG de scalp des régions communes EEG et IRMf et le supprimer des observations. Nous pourrions également supprimer les sommets du maillage liées aux régions communes EEG et IRMf du modèle direct et utiliser des méthodes de localisation EEG pour trouver des sites d'activation supplémentaires dans EEG dans les données résiduelles après suppression des régions communes et de leurs signaux. Les sites d'activation découverts par cette méthode conduiraient à des régions qui pourraient être liées à la tâche ou au fond, ou un faux positif. L'idée doit être étudiée plus en détails afin de rechercher des sites d'activation fiables. Remerciements Ce travail a été partiellement financé par le projet CHESS, 2012-ERC-AdG-320684. Les auteurs tiennent à remercier Daniel Wakeman et Richard Henson pour fournir les ensembles de données réelles.

Bibliography

The disadvantage of the proposed method is that the final result is dependent on the first source separation step and any mistake in the first step will propagate to the rest. However, we did not have any false discovery in the simulation experiments.

In the future work, the proposed method may be applied to simultaneous EEG and fMRI data. Defining different weighting matrices for � 2 and � 0 terms may be investigated and also, the optimized value of α, which controls the contribution of the fMRI spatial map in the weighting matrix, can be defined using the data.

Chapter 7

Evaluation of the EEG-fMRI integration method

Introduction

In this chapter the integration method is evaluated with simulated data (Section 7.2.1) and the actual data (Section 7.2.2). Results are presented in Section 7.3. Discussion and concluding remarks are held in Sections 7.4 and 7.5.

Data

Simulated Data

Monte Carlo modeling (sampling over randomized 50 source locations) is used to compare different methods with the proposed methods. For each activation site a related set of EEG and fMRI data is created. 50 random vertices (voxels) are selected as the center of their active regions. The sources are assumed to be spread, so a smoothing Gaussian filter is used to model the spread sources.

Here, we simulate the EEG and fMRI data after the preprocessing step. It means that for example in fMRI, slice timing is not simulated and in EEG, the mean of the trials after preprocessing is simulated.

The experiment has two states: task and control. Onsets of task and control states are set as the onsets of the clinical data for face and scrambled states, respectively. All the parameters of imaging methods are chosen as the clinical data, like the electrode positions and sampling rate in EEG and TR and image dimension in fMRI.

fMRI

We use the package neuRosim in R [START_REF] Core | R: A Language and Environment for Statistical Computing[END_REF] for fMRI simulation, because the data generation in neuRosim is fairly fast [START_REF] Welvaert | {neuRosim}: An {R} Package for Generating {fMRI} Data[END_REF]. A double-gamma function, which is a linear combination of two gamma functions, is used to model HRF [START_REF] Friston | Eventrelated fMRI: characterizing differential responses[END_REF][START_REF] Glover | Deconvolution of impulse response in event-related BOLD fMRI[END_REF]. Fifty fMRI simulations are created with � Time courses of the electrodes. To compute the electrodes time courses, we use the forward linear model (6.1). To estimate the gain matrix, we use the structural MRI of a real subject with its fiducials and electrode positions. For this purpose, we use BEM with 8196 nodes developed in SPM8 (http://www.fil.ion.ucl.ac.uk/spm). We assume that the dipoles are on the cortical gray matter, where each dipole is perpendicular to the cortical sheet.

We put the task and control time courses in the related rows of the matrix J. We put the estimated activation sources in the predefined sites and use a Gaussian filter with full width at half maximum (FWHM) of 8 mm to estimate the activation of the neighboring mesh vertices. After computing the effects of activations on the electrode space, we add the background signal to the electrodes with different SNR values in the range of [-20, 5] dB and the step size of 5 dB.

Clinical Data

Clinical data is taken from ftp://ftp.mrc-cbu.cam.ac.uk/personal/rik.henson/wakemandg_ hensonrn/. A brief summary of the acquisition characteristics, and preprocessing steps are brought in this section. More details can be found in [START_REF] Henson | MEG and EEG data fusion: Simultaneous localisation of face-evoked responses[END_REF][START_REF] Henson | A Parametric Empirical Bayesian framework for the EEG/MEG inverse problem: generative models for multisubject and multimodal integration[END_REF][START_REF] Wakeman | Functional and structural connectivity in faceprocessing: MEG, EEG, fMRI, MRI and DWI data[END_REF].

Participants

The participants were sixteen healthy young adults (eight female). The study protocol was approved by the local ethics review board (CPREC reference 2005.08). Each subject performed two days experiments; once for concurrent MEG + EEG and once for fMRI + MRI.

Stimuli and Tasks

The paradigm was similar to that used previously under EEG, MEG, and fMRI [START_REF] Henson | MEG and EEG data fusion: Simultaneous localisation of face-evoked responses[END_REF][START_REF] Henson | A Parametric Empirical Bayesian framework for the EEG/MEG inverse problem: generative models for multisubject and multimodal integration[END_REF]. A central fixation cross (presented for a random duration of 400-600 ms) was followed by a face or scrambled face (presented for a random duration of 800-1000 ms), followed by a central circle for 1700 ms (See Figure 7.2).

EEG Acquisition

The EEG data was recorded simultaneously with MEG with a VectorView system (Elekta Neuromag, Helsinki, Finland). EEG data was recorded with 70 Ag-AgCl electrodes with sampling rate of 1.1 kHz. The low-pass filter with cut off frequency of 350 Hz was applied to the EEG data.

An elastic cap (EASYCAP GmbH, Herrsching-Breitbrunn, Germany) was used according to the extended 10-10 and a nose electrode was used as the recording reference.

Appendix A

False positive and false negative in terms of electrode leads F P e and F N e are calculated between members of IED # and GT . IED # is the set of estimated IED electrode leads of each recalculated trial and members of GT is assumed as the ground truth.

We calculate the Euclidean distance, d ij , between the ith member of IED # and jth member of GT , where i = 1, . . . , N � , j = 1, . . . , N . We threshold d ij giving b ij = 1 if d ij < th LF P , otherwise zero. LFP stands for local field potential. According to [208], the LFP recorded from each electrode lead can be related to a neural population within 0.5 -3 mm of the electrode tip. In [START_REF] Ulmer | FMRI: Basics and Clinical Applications[END_REF], an LFP coherence about 0.15 -0.35 (0 <coherence< 1) for approximately 3 -4 mm primary visual cortical distance was reported in the frequency range of 2 -60 Hz. In our recordings, the inter-distance between two adjacent electrode leads is 3.5 mm, so for the thresholds less than this distance, there is no neighbor electrode lead at least on a single electrode.

As such th LF P equal to 4 mm is chosen. F P e and F N e are calculated as:

UNIVERSITE DE GRENOBLE INSTITUT POLYTECHNIQUE DE GRENOBLE