Hanna Becker

Laurent Albera
email: laurent.albera@univ-rennes1.fr

Pierre Comon

Martin Haardt

Gwénaël Birot

Fabrice Wendling

Martine Gavaret

Christian-George G Bénar

Isabelle Merlet

Ltsi

EEG extended source localization: Tensor-based vs. conventional methods

Keywords: EEG, Tensor decomposition, Space-Time-Frequency analysis, Space-Time-Wave-Vector analysis, Distributed source localization

come   L'archive ouverte pluridisciplinaire

Introduction

In the context of drug-resistant partial epilepsy, scalp Electroencephalographic (EEG) recordings constitute the first step in the pre-surgical evaluation of patients, and can be used to delineate the regions from where epileptic paroxysms arise. To this end, source localization algorithms have been applied to interictal spikes that can be frequently observed on the scalp with a high Signal-to-Noise-Ratio (SNR) compared to ictal discharges [START_REF] Ebersole | Noninvasive localization of epileptogenic foci by EEG source modeling[END_REF][START_REF] Merlet | Dipole modeling of interictal and ictal EEG and MEG[END_REF][START_REF] Michel | EEG source imaging[END_REF][START_REF] Gavaret | Electric source imaging in frontal lobe epilepsy[END_REF][START_REF] Plummer | EEG source localization in focal epilepsy: Where are we now?[END_REF]. A number of techniques for the localization of several equivalent dipoles each of which models an entire active brain region have been proposed (see [START_REF] Grech | Review on solving the inverse problem in EEG source analysis[END_REF] for a review). To accurately localize extended sources, for which the equivalent dipole model is not appropriate, one can further resort to a multipole source model as suggested in [START_REF] Jerbi | Localization of realistic cortical activity in MEG using current multipoles[END_REF]. Now, if both the localization of the distributed sources and their spatial extent need to be identified such as in the context of drug-resistant partial epilepsy, a distributed source model can be employed where each source is characterized by a large number of dipoles and a complex geometry. The determination of these distributed sources constitutes an ill-posed problem for which various types of methods based on different assumptions on the source activity have been developed. One class of methods which is referred to as minimum norm solutions comprising techniques like the original minimum norm method introduced in [START_REF] Hämäläinen | Interpreting magnetic fields of the brain: Estimates of current distributions[END_REF], LORETA (low resolution brain electromagnetic tomography) [START_REF] Pascual-Marqui | Low resolution electromagnetic tomograph: A new method for localizing elecrical activity in the brain[END_REF], LAURA (local autoregressive average) [START_REF] Peralta Menendez | Electrical neuroimaging based on biophysical constraints[END_REF], sLORETA (standardized low resolution brain electromagnetic tomography) [START_REF] Pascual-Marqui | Standardized low resolution brain electromagnetic tomography (sLORETA): technical details[END_REF], and dSPM (dynamic statistical parametric mapping) [START_REF] Dale | Dynamic statistical parametric mapping: combining fMRI and MEG for high-resolution imaging of cortical activity[END_REF] reconstructs the source activity at every point of a given grid in a tomographic way based on constraints such as, for example, maximal smoothness in the case of LORETA. However, these methods are reported to yield blurred results [START_REF] Xu | Lp norm iterative sparse solution for EEG source localization[END_REF]. Several techniques have been developed to improve the spatial resolution of the minimum norm solutions, but they tend to result in too sparse source estimates [START_REF] Xu | Lp norm iterative sparse solution for EEG source localization[END_REF]. Moreover, except for sLORETA, the minimum norm solutions are biased (for a single dipole) even in the absence of noise [START_REF] Sekihara | Localization bias and spatial resolution of adaptive and non-adaptive spatial filters for MEG source reconstruction[END_REF]. A comparison of several distributed source localization methods including minimum norm and maximum entropy of the mean (MEM) approaches on realistic simulated data of extended epileptic sources can be found in [START_REF] Grova | Evaluation of EEG localization methods using realistic simulations of interictal spikes[END_REF].

Another source localization approach includes statistical techniques such as MUSIClike methods [START_REF] Mosher | Multiple dipole modeling and localization from spatio-temporal MEG data[END_REF][START_REF] Mosher | Source localization using recursively applied and projected (RAP) MUSIC[END_REF][START_REF] Albera | Brain source localization using a fourth-order deflation scheme[END_REF]. The proposed algorithms in this field are generally conceived for equivalent dipole localization. However, recently, the 2q-ExSo-MUSIC (2q-th order extended source multiple signal classification) algorithm [START_REF] Birot | Localisation of extended brain sources from EEG/MEG: the ExSo-MUSIC approach[END_REF] has been developed for the localization of distributed sources that summarize a large number of dipoles with highly synchronized activities. This model describes well the epileptic paroxysms that are observable by EEG measurements and that often involve large cortical regions as has been shown by several studies [START_REF] Cosandier-Rimélé | The neuronal sources of EEG: Modeling of simultaneous scalp and intracerebral recordings in epilepsy[END_REF][START_REF] Gavaret | Electric source imaging in frontal lobe epilepsy[END_REF][START_REF] Tao | Intracranial EEG substrates of scalp EEG interictal spikes[END_REF][START_REF] Merlet | Reliability of dipole models of epileptic spikes[END_REF][START_REF] Ebersole | Magnetoencephalography/magnetic source imaging in the assessment of patients with epilepsy[END_REF]. Similarly, studies with simultaneous electrocorticography (ECoG) and magnetoencephalography (MEG) recordings [START_REF] Mikuni | Simultaneous recording of epileptiform discharges by MEG and subdural electrodes in temporal lobe epilepsy[END_REF][START_REF] Oishi | Epileptic spikes: magnetoencephalography versus simultaneous electrocorticography[END_REF][START_REF] Shigeto | Feasibility and limitations of magnetoencephalographic detection of epileptic discharges: simultaneous recording of magnetic fields and electrocorticography[END_REF] have demonstrated that a certain cortical area needs to be active in order to observe epileptic spikes at the surface MEG sensors. From an ensemble of circular-shaped source regions (or "disks"), that comprises the entire cortex, 2q-ExSo-MUSIC selects the disks which optimize a metric computed from the 2q-th order signal subspace of the data. Yet this method makes it difficult and computationally very demanding to localize several correlated, simultaneously active source regions, which we subsequently refer to as patches.

To simplify the problem of localizing several potentially correlated patches in the presence of background activity, it is desirable to apply a preprocessing technique to the measurement data that separates the simultaneously active source regions into different components and reduces the noise. In the past, several authors have studied the application of the Canonical Polyadic (CP) decomposition to Space-Time-Frequency (STF)transformed EEG data, which is obtained by computing a wavelet transform [START_REF] Miwakeichi | Decomposing eeg data into space-time-frequency components using parallel factor analysis[END_REF][START_REF] Morup | Parallel factor analysis as an exploratory tool for wavelet transformed event-related EEG[END_REF]De Vos et al., 2007a,b;[START_REF] Deburchgraeve | Neonatal seizure localization using parafac decomposition[END_REF] or the Wigner-Ville distribution [START_REF] Weis | Mutli-dimensional space-time-frequency component analysis of event related EEG data using closed-form PARAFAC[END_REF] over the time dimension of the measurements. Under certain conditions on the signals, this method provides separate space, time, and frequency characteristics for each source region and therefore allows us to treat each patch individually in the the next steps.

In [START_REF] Becker | Multiway space-timewave-vector analysis for source localization and extraction[END_REF][START_REF] Becker | Multi-way space-timewave-vector analysis for EEG source separation[END_REF], we have introduced an alternative method, which is based on a local spatial Fourier transform of the EEG measurements. This leads to a Space-Time-Wave-Vector (STWV) tensor, that can also be decomposed using the CP model. An advantage of this approach compared to the STF analysis consists in its robustness to correlated source activities. This is of particular interest when patches with identical, but shortly delayed, source activities have to be identified, i.e., in the context of interictal spikes when spreading of epileptic spikes is suspected between two regions. Furthermore, we have addressed the topic of extended source localization after the tensor-based preprocessing by proposing the Binary Coefficient Matching Pursuit (BCMP) algorithm [START_REF] Becker | Multi-way space-timewave-vector analysis for EEG source separation[END_REF]. This technique iteratively selects dipoles from a predefined grid to model a larger source region. It was shown to accurately recover sources of some extent in the context of a spherical head model and of simplistic simulated EEG signals.

The next step in the validation process of these tensor methods and in particular of the STWV analysis consists in evaluating the algorithms on more realistic simulated data and in testing them on real EEG measurements, which is one of the objectives of the present paper. Additionally, we analyze the theoretical basis of the tensor-based methods to understand the underlying mechanisms and conditions that are necessary for the STF and STWV techniques to work. To our knowledge, this has not been studied yet.

In order to demonstrate the interest of tensor-based preprocessing for extended source localization applications, we conduct an extensive analysis of the performance of the tensor-based methods in comparison to conventional methods, both in terms of computational complexity and source localization accuracy. The latter is examined using realistic simulated data that is generated for a realistic head model and extended sources whose signals are modeled with a physiologically relevant model [START_REF] Cosandier-Rimélé | A physiologically plausible spatio-temporal model for EEG signals recorded with intracerebral electrodes in human partial epilepsy[END_REF].

Since the BCMP algorithm fails to accurately recover the patches of the realistic model, we propose a new algorithm for extended source localization which we refer to as the disk algorithm (DA). This method uses an optimization strategy that is inspired by the ExSo-MUSIC source localization procedure [START_REF] Birot | Localisation of extended brain sources from EEG/MEG: the ExSo-MUSIC approach[END_REF], but is based on a different metric for the disk selection since it identifies the disks that best describe the spatial characteristics of only one source or one part of a source. Furthermore, we compare the performance of this algorithm to that of 4-ExSo-MUSIC as a representative of subspacebased approaches for extended source localization and sLORETA and cortical LORETA (cLORETA) as representatives of the family of minimum norm solutions. We study the influence of the patch size and of the SNR related to the patch anatomical location and determine the limits of the examined methods regarding the localization of deep and closely-spaced patches.

Finally, we apply the different source localization methods to actual EEG interictal spikes, obtained from a patient with frontal lobe epilepsy, and compare the results to the epileptogenic brain regions involved during interictal spikes as identified from intracerebral stereotactic electroencephalography (SEEG) recordings of the same patient. This paper is organized as follows: in Section 2, we present the data model on which the examined methods are based and review the concept of tensor-based preprocessing as well as of the STF and STWV analyses. Furthermore, we describe the newly proposed disk algorithm as well as the other examined extended source localization algorithms and analyze their computational complexity. Finally, we provide information about the simulation setup. Section 3 then presents the results of the conducted simulation study and the analysis of actual EEG data. A discussion of the results and their interpretation is given in Section 4, which concludes the paper. A theoretical analysis of the STF and STWV methods and sufficient conditions under which these methods yield exact results are made available in the online supplementary material.

The following notation is used throughout the paper: bold italic uppercase letters denote tensors, e.g., T , bold uppercase letters denote matrices, e.g., A, bold lowercase letters denote column vectors, e.g., a, and plain font denotes scalars, e.g., X i,j,k , T i,j or a i . Moreover, (•) T denotes a transposition and (•) + stands for the Moore-Penrose pseudoinverse.

Methods

Data model

It is generally admitted that the electrical activity that can be observed at the surface of the scalp originates mostly from pyramidal cells that are located in the gray matter.

These cells are arranged in parallel and exhibit an orientation that is perpendicular to the cortex. In order to obtain a signal of sufficient amplitude to be measurable at the surface of the scalp, a certain number of simultaneously active neuronal populations is required.

These populations can be modeled by a grid of dipole sources located at the surface of the cortex [START_REF] Dale | Improved localization of cortical activity by combining EEG and MEG with MRI cortical surface reconstruction: a linear approach[END_REF]. The electric potential data that is recorded at the N r electrodes of an EEG sensor array for N t time samples then constitutes the superposition of all dipole signals contained in the signal matrix S ∈ R N d ×Nt that are transmitted to the surface of the scalp and some instrumentation noise X i :

X = GS + X i (1)
where G ∈ R Nr×N d is the lead field matrix that describes the propagation in the head volume conductor and depends on spatial parameters of the head, such as the geometry of the brain, skull, and scalp as well as their conductivities, and N d denotes the number of grid dipoles. For a given head model and source space, the lead field matrix can be computed numerically (see also Section 2.5.1).

In the context of epilepsy, the regions of interest can be modeled by extended sources that can be described as the union of (one or) several non-necessarily contiguous areas of the cortex (so-called patches) with highly correlated source activities [START_REF] Cosandier-Rimélé | A physiologically plausible spatio-temporal model for EEG signals recorded with intracerebral electrodes in human partial epilepsy[END_REF][START_REF] Birot | Localisation of extended brain sources from EEG/MEG: the ExSo-MUSIC approach[END_REF]. Assuming that artifacts have already been removed from the data as can be achieved by methods as described, for instance, in (Comon and Jutten, 2010;[START_REF] Albera | ICA-based EEG demonising: a comparative analysis of fifteen methods[END_REF], all dipoles that do not belong to an extended source can be considered to emit background activity. Consequently, in order to distinguish between the extended sources, that we are looking for, and the noisy background activity, we can rewrite the data model (1) in the following way:

X = R r=1 kr∈Ωr g kr s T kr + l / ∈∪ R r=1 Ωr g l s T l + X i (2) X = X + X b + X i = X + N. (3)
Here Ω r is the set of indices of the dipoles that belong to the r-th extended source, g k is the lead field vector of the k-th dipole, and s k is the associated signal vector that corresponds to the k-th row vector of S. The matrix X comprises the data generated by the signals of interest whereas the noise matrix N summarizes background activity X b and instrumentation noise X i .

As the activity is highly synchronized within one extended source, we can assume that all signals emitted by dipoles of the same extended source are quasi-equal, i.e., s kr ≈ s r .

Consequently, we can approximate the data matrix by

X ≈ R r=1 h r s T r + N = HS + N (4)
where the spatial mixing vector h r of the patch r corresponds to the linear combination of the lead field vectors of all active grid dipoles:

h r = kr∈Ωr g kr = Gc r (5)
which can also be written as the product of the lead field matrix G and a coefficient vector c r whose k-th element, k = 1, . . . , N d , is given by

c k,r =    1 if k ∈ Ω r 0 otherwise.

Tensor-based preprocessing

To separate several simultaneously active patches and to reduce the background activity, tensor-based preprocessing can be applied. The idea of this method consists in exploiting the structure of three-way data which is obtained by applying a transform to the two-dimensional measurements. Under the hypothesis that the resulting data, which depends on three variables, is trilinear, the tensor can be decomposed in a unique way (under mild conditions) up to scale and permutation ambiguities into separate characteristics for each variable with the help of the CP decomposition (also sometimes referred to as Parallel Factor Analysis (PARAFAC)). It is thus possible to get an accurate estimate of the spatial mixing matrix or the signal matrix without imposing statistical constraints on sources (unlike Independent Component Analysis (ICA) (Comon and Jutten, 2010)),

which may be physiologically difficult to justify.

CP decomposition

Each element of a third order tensor X of size I × J × K can be written in the form:

X i,j,k = P p=1 a i,p b j,p c k,p (6)
which is generally called a polyadic decomposition of the tensor X. If P is the smallest integer for which equality (6) holds, P corresponds to the rank of the tensor. In this case, the trilinear representation ( 6) is referred to as the canonical polyadic (CP) decomposition of the tensor X [START_REF] Comon | Tensor decompositions, alternating least squares and other tales[END_REF]. The variables a i,p , b j,p , and c k,p are elements of three matrices A ∈ C I×P , B ∈ C J×P , and C ∈ C K×P , respectively, called loading matrices. There is almost surely a finite number of decompositions of the tensor X into the three loading matrices A, B and C, up to scale and permutation indeterminacies, if P < IJK I+J+K-2 [START_REF] Comon | Tensor decompositions, alternating least squares and other tales[END_REF]; tighter bounds also give sufficient conditions for uniqueness such as the condition based on Kruskal rank that has been derived in [START_REF] Kruskal | Three-way arrays: Rank and uniqueness of trilinear decompositions[END_REF][START_REF] Sidiropoulos | On the uniqueness of multilinear decompositions of n-way arrays[END_REF].

The CP decomposition can be computed using Alternating Least Squares (ALS) [START_REF] Bro | Multi-way analysis in the food industry: Models, algorithms and applications[END_REF], efficient gradient-based algorithms [START_REF] Phan | Low complexity damped Gauss-Newton algorithms for CANDECOMP/PARAFAC[END_REF][START_REF] Sorber | Optimization-based algorithms for tensor decompositions: canonical polyadic decomposition, decomposition in rank-(l r , l r , 1) terms and a new generalization[END_REF], or semialgebraic methods as presented in [START_REF] Luciani | Semi-algebraic canonical decomposition of multi-way arrays and joint eigenvalue decomposition[END_REF]Römer andHaardt, 2008, 2013). Because of its good performance and robustness to collinear factors, overestimation of the number of CP components, and initialization, we employ the DIAG (Direct Algorithm for canonical polyadic decomposition) algorithm (Luciani andAlbera, 2011, 2014).

To get a 3-dimensional data tensor which can be treated by the CP decomposition, one can either compute a transform over time of the electric potential measurements, which leads to the STF analysis, or a transform over space, yielding STWV data. These methods will be described in the subsequent sections.

STF analysis

An often used technique for the time-frequency analysis of EEG data consists in applying a wavelet transform to the time signals {x(r, t)} of the different channels [START_REF] Miwakeichi | Decomposing eeg data into space-time-frequency components using parallel factor analysis[END_REF][START_REF] Morup | Parallel factor analysis as an exploratory tool for wavelet transformed event-related EEG[END_REF]De Vos et al., 2007a,b;[START_REF] Deburchgraeve | Neonatal seizure localization using parafac decomposition[END_REF]. The resulting three-way data can then be stored into the data tensor

W (r, t, f ) = ∞ -∞ x(r, τ )ψ(a, τ, t)dτ . (7)
The frequency f can be estimated from the scale a of the wavelet ψ(a, τ, t) by f = f c /(aT

)
where f c is the center frequency of the wavelet and T is the interval between time samples.

In order to decompose the tensor W using the CP decomposition, we assume that for each extended source, the time and frequency variables separate, leading to a trilinear tensor. This is approximately the case under the hypothesis of oscillatory signals. The tensor can then be decomposed as:

W (r i , t j , f k ) ≈ R p=1 a(r i ; p)b(t j ; p)c(f k ; p) (8)
where r i , t j , and f k represent the sampled space, time, and frequency variables and a(r i ; p), b(t j ; p), and c(f k ; p) denote elements of the loading matrices A, B, and C indicating the space, time, and frequency characteristics, respectively. The number of components R corresponds to the number of extended sources.

The loading matrix A containing the spatial characteristics generally constitutes a good estimate for the spatial mixing matrix H.

STWV analysis

If a local spatial Fourier transform is calculated within a certain region on the scalp, selected by the spherical window function w(r ′r) centered at sensor position r (see Section 2.5.2 for more details), the STWV tensor

F (r, t, k) = ∞ -∞ w(r ′ -r)x(r ′ , t)e jk T r ′ dr ′ (9)
is obtained [START_REF] Becker | Multi-way space-timewave-vector analysis for EEG source separation[END_REF]. Here, the third variable k is the wave vector.

Under the assumption that the space and wave vector variables separate for each extended source, which is approximately the case for superficial sources, the tensor F can be approximated by the CP model and be decomposed into space, time, and wave vector characteristics a(r i ; p), b(t j ; p), and c(k l ; p):

F (r i , t j , k l ) ≈ R p=1 a(r i ; p)b(t j ; p)c(k l ; p). (10)
In the case of the STWV analysis, the temporal characteristics Ŝ = B constitute a good approximation of the signal matrix S. An estimate Ĥ for the lead field matrix H can thus be obtained from the pseudo-inverse Ŝ+ of the estimated signal matrix Ŝ and the data matrix X:

Ĥ = X Ŝ+ . (11)

Trilinear approximation

Even though the STF analysis has been widely used, up to now only intuitive conditions such as oscillatory signals that presumably lead to trilinear data have been provided.

But no theoretical validation that justifies the application of the CP decomposition to the STF data tensor has been performed and the mechanisms underlying the STF method are still insufficiently explored. The same is true for the STWV technique. Therefore, we

have analyzed what happens when applying the DIAG algorithm to STF or STWV data, that are not exactly trilinear, and have clarified under which conditions this procedure yields exact results. For the interested reader, our findings are made available in the online supplementary material and discussed with respect to the results obtained for two scenarios that are considered in Section 3.2.

Extended source localization

Based on Equation ( 5), the objective of extended source localization consists in recovering the coefficient vectors c r , r = 1, . . . , R, for all patches. This can be achieved based on an estimate Ĥ = ĥ1 , . . . , ĥR of the spatial mixing matrix H, which permits to identify the grid dipoles for each patch separately. In principle, any localization algorithm that acts on a vector of spatial measurements can be employed to this end. Here, we introduce the disk algorithm (DA) that uses an optimization strategy inspired by the ExSo-MUSIC approach [START_REF] Birot | Localisation of extended brain sources from EEG/MEG: the ExSo-MUSIC approach[END_REF] but with a different metric built from the spatial mixing matrix estimated by the STF or STWV analysis. In addition, we consider surface grid based versions of the sLORETA [START_REF] Pascual-Marqui | Standardized low resolution brain electromagnetic tomography (sLORETA): technical details[END_REF] and LORETA [START_REF] Pascual-Marqui | Low resolution electromagnetic tomograph: A new method for localizing elecrical activity in the brain[END_REF] algorithms, which belong to the class of minimum norm estimates, also called generalized inverse solutions, and which are applied to a data vector x. Furthermore, we employ a statistical subspace-based approach, 4-ExSo-MUSIC [START_REF] Birot | Localisation of extended brain sources from EEG/MEG: the ExSo-MUSIC approach[END_REF], that is applied to the data matrix X. The sLORETA, cLORETA, and 4-ExSo-MUSIC methods are used as references to evaluate the results of the tensor-based DA. Subsequently, we give a short description of the proposed source localization algorithm, that is based on the tensor decomposition, and the three reference methods.

Disk algorithm (DA)

To localize extended sources based on the estimated spatial mixing matrix obtained by the STWV and STF analyses, we present in this section the newly developed disk algorithm. In addition to the hypotheses that are made on the structure of the STF or STWV tensor in order to separate the sources, this approach assumes that the measurements can be described by extended sources with a piece-wise constant spatial distribution according to data model ( 4) and ( 5). The concept underlying the disk algorithm consists in recovering the extended source from a number of small circular-shaped patches of grid dipoles, the disks. Each extended source is then obtained by the union of several disks.

For each grid dipole, several disks composed of the 0 to D max -1 nearest dipoles and the current grid dipole as central point are determined. Assuming equal coefficients for all grid dipoles belonging to the extended source, the spatial mixing vector of each disk is computed by adding the lead field vectors of all grid dipoles belonging to the disk. The reconstructed spatial mixing vector d k of the k-th disk, k = 1, . . . , D max N d , can then be obtained from

d k = Gĉ k (12)
where the elements of the coefficient vector ĉk that are associated with dipoles belonging to the disk are equal to 1 whereas all other elements are 0. To determine which disks of the parameter space best describe the measurements, the spatial mixing vectors d k of all disks are then compared to the estimated spatial mixing vector ĥr of the r-th source with the help of the following metric, which is based on the normalized inner product:

M iprod (ĉ r ) = - ĥT r Gĉ r 2 ĉT r G T Gĉ r . ( 13
)
This metric is minimized with respect to the coefficient vector ĉr . The norm of the spatial mixing vector ĥr is fixed. The extended source localization procedure then selects a certain number of disks which are associated with the smallest values of this metric.

sLORETA

Minimum norm solutions are based on the assumption that among all possible source configurations that explain the measurements, the source distribution with minimal energy is the most likely. In this context, the idea of the conventional sLORETA algorithm [START_REF] Pascual-Marqui | Standardized low resolution brain electromagnetic tomography (sLORETA): technical details[END_REF] consists in standardizing the minimum norm solution for the current density distribution, computed on a voxel grid inside the brain, by weighting the current dipoles by their estimated variances. This method can be adapted to our data model by replacing the voxel grid by a surface grid that considers only dipoles on the cortical surface with an orientation perpendicular to this surface. Due to the assumption of radial orientation, we solve the inverse problem for one coefficient per dipole which represents its amplitude instead of recovering three components for each dipole in order to determine both dipole strength and orientation as the conventional sLORETA algorithm does. Thanks to an inversion lemma, the solution of the optimization problem arg min

c ||x -Gc|| 2 + λ||c|| 2 ,
where λ is a regularization parameter, can be written as

ĉ = Qx. (14)
Here, more details. This leads to unbiased estimates for single dipole sources [START_REF] Sekihara | Localization bias and spatial resolution of adaptive and non-adaptive spatial filters for MEG source reconstruction[END_REF]. The i-th dipole coefficient, i = 1 . . . N d , can then be computed as

Q = G T GG T + λI -1 ∈ R N d
1 ĉi,r = (K i,i ) -1 2 Q i,: x (15)
The resulting coefficient vector ĉ contains continuous values for the dipole strengths. For the localization of extended sources, all the grid dipoles whose coefficients exceed a certain threshold are associated to the extended source.

cLORETA

The objective of the LORETA algorithm [START_REF] Pascual-Marqui | Low resolution electromagnetic tomograph: A new method for localizing elecrical activity in the brain[END_REF] consists in identifying a maximally smooth source distribution. Similarly to sLORETA, it has originally been conceived for volumetric source localization, but can be modified to work on a surface grid [START_REF] Wagner | Smooth reconstruction of cortical sources from EEG and MEG recordings[END_REF]. In this case, it is also referred to as cLORETA. To obtain a smooth source distribution, cLORETA solves the following optimization problem:

min c ||x -Gc|| 2 + λ||LWc|| 2 (16)
where L corresponds to the surface Laplacian matrix whose elements are given by [START_REF] Wagner | Smooth reconstruction of cortical sources from EEG and MEG recordings[END_REF])

L i,j =          -1 j d i,j j 1 d i,j for i = j 1 d i,j j d i,j
if i and j are indices of adjacent dipoles 0 otherwise where d i,j denotes the distance between dipoles i and j, i, j = 1, . . . , N d , and W is a diagonal weight matrix with W i,i = ||g i || -1 . The solution to ( 16) can be computed as

ĉ = WL T LW -1 G T G WL T LW -1 G T + λI -1
x.

(17)

In order to obtain an estimate of the extended source, the resulting coefficient vector ĉ is thresholded and only the grid dipoles with coefficients exceeding the threshold are retained to characterize the extended source.

4-ExSo-MUSIC

Contrary to the other source localization algorithms used in this paper, which are deterministic, the 4-ExSo-MUSIC algorithm [START_REF] Birot | Localisation of extended brain sources from EEG/MEG: the ExSo-MUSIC approach[END_REF] interprets the N t measurement vectors x contained in the data matrix X as N t realizations of a random vector.

The algorithm exploits the estimated Fourth Order (FO) cumulants of the latter random vector, which are contained in the quadricovariance matrix C 4 . In order to resort to the FO statistics of the data, which offer an asymptotic robustness with respect to Gaussian noise and a better resolution than second order (SO) statistics [START_REF] Chevalier | High-resolution direction finding from higher order statistics: the 2q-MUSIC algorithm[END_REF][START_REF] Birot | Localisation of extended brain sources from EEG/MEG: the ExSo-MUSIC approach[END_REF], the FO marginal cumulants of the epileptic source signals are assumed to be non-zero. The 4-ExSo-MUSIC algorithm is then based on the FO signal subspace, which is computed by an eigenvalue decomposition (EVD) of the quadricovariance matrix:

C 4 = E s E n   Λ s 0 0 Λ n     E T s E T n   . (18)
Here, the matrices E s ∈ R N 2 ×P and E n ∈ R N 2 ×(N 2 -P ) span the signal and noise subspaces where P corresponds to the dimension of the signal subspace and Λ s ∈ R P ×P and P ) are diagonal matrices containing the signal and noise eigenvalues.

Λ n ∈ R (N 2 -P )×(N 2 -
The source localization is then accomplished based on data model ( 4) and the hypothesis of extended sources with a piece-wise constant spatial distribution by performing an exhaustive search over a parameter space that is composed of disks which are constructed as described in Section 2.3.1. Nevertheless, the metric which is optimized in order to find the appropriate union of disks is different of that used in the DA algorithm. More particularly, the 4-ExSo-MUSIC algorithm identifies the disks whose spatial mixing vectors d k = Gĉ k are closest to the signal subspace, leading to the MUSIC metric

M MUSIC (ĉ) = (d ⊗ d) T E s E T s (d ⊗ d) (d ⊗ d) T (d ⊗ d) = (ĉ ⊗ ĉ) T (G ⊗ G) T E s E T s (G ⊗ G)(ĉ ⊗ ĉ) (ĉ ⊗ ĉ) T (G ⊗ G) T (G ⊗ G)(ĉ ⊗ ĉ) , (19)
which is maximized with respect to the coefficient vector ĉ. Here, d ⊗ d denotes the Kronecker product of the spatial mixing vector d with itself, yielding the FO spatial mixing vector. Since 4-ExSo-MUSIC does not separate the sources, all source regions are recovered simultaneously and are characterized by the union of all disks for which the metric ( 19) is above a certain threshold value.

Study of the computational complexity

In order to compare the different extended source localization methods relative to their computational cost, we have analyzed the number of real-valued multiplications that are required for each step of the different source localization processes, comprising tensor construction, tensor decomposition, cumulant and spectrum estimation, and source localization algorithms. The results are summarized in a table that is available in the online supplementary material.

Simulations

Data generation

In order to generate physiological data, we employ a realistic head model with three compartments representing the brain, the skull, and the scalp, whose surfaces are obtained from a normal MRI. The data is generated for N r = 91 EEG sensors. The source space is defined by the triangularized inner cortical surface (grey matter / white matter interface),

where a grid dipole is placed at the centroid of each of the triangles. The grid consists of 19626 triangles (9698 for the left hemisphere and 9928 for the right hemisphere) and on average, each triangle describes 5 mm 2 of the cortical surface. The lead field vectors contained in the matrix G ∈ R 91×19626 are then computed numerically for all grid dipoles using a Boundary Element Method (ASA, ANT, Enschede, Netherlands).

For the generation of extended sources, we consider a number of patches each of which consists of 100 adjacent grid dipoles corresponding to a cortical area of approximately 5 cm 2 . Using a model of coupled neuronal populations as previously described [START_REF] Wendling | Relevance of nonlinear lumped-parameter models in the analysis of depth-EEG epileptic signals[END_REF][START_REF] Cosandier-Rimélé | A physiologically plausible spatio-temporal model for EEG signals recorded with intracerebral electrodes in human partial epilepsy[END_REF][START_REF] Cosandier-Rimélé | Computational modeling of epileptic acitivity: from cortical sources to EEG signals[END_REF],

Tensor-based preprocessing

To separate several simultaneously active patches and to reduce the noise on the estimated lead field matrix, we employ tensor-based preprocessing. Here, we consider the application of both the STF and the STWV analysis. To this end, we construct the STF data tensor by computing a wavelet transform over time of the electric potential data.

Due to the application of a real-valued Morlet-wavelet, the resulting tensor W is also real-valued and of size N r × N t × N f where N r , N t , and N f denote the number of sensors, time samples, and frequency samples, respectively. Here, we consider N f = 50 frequency samples and frequencies ranging from 0.5 to 20 Hz.

To obtain STWV transformed data, we select sensor data with the help of a spherical window function that is centered at different electrode positions on the scalp and attribute a weight 1 to sensors located within the sphere and 0 to all other sensors. Then we compute a discrete non-uniform spatial short term Fourier transform of the selected data. This leads to the complex-valued tensor F ∈ C N ′ r ×Nt×N k , where N ′ r is the effective number of sensors. To ensure meaningful results, the local Fourier transform is only computed for N ′ r manually selected sensors which are deemed to be surrounded by a sufficient number of electrodes, leading to an effective number of sensors which is strictly less than N r .

Here, we choose a radius of 3.3 cm for the spherical window function, and select 71 electrodes, surrounded by 7 to 12 sensors within the spherical window, for the local spatial Fourier transform. The variable N k = 63 denotes the number of wave vector

samples k = [k 1 , k 2 , k 3 ] T ,
which are chosen arbitrarily to contain all possible combinations of k 1 , k 2 , k 3 ∈ {0, ±1, ±2} such that there are no two wave vectors k 1 and k 2 for which

k 1 = -k 2 .
The tensors are then decomposed using the DIAG algorithm (Luciani andAlbera, 2011, 2014). In case of the STWV tensor, it is followed by 10 ALS iterations to ensure that the loading matrix B that contains the temporal characteristics of the signal is real-valued. The number of CP components is chosen such that it equals the number of patches. Finally, the spatial mixing matrix is estimated according to the STF and STWV methods.

Extended source localization

The patches are localized based on the spatial mixing vectors obtained by the STF and STWV preprocessing using DA. Depending on the tensor method used for the preprocessing, we subsequently refer to these techniques as STF-DA and STWV-DA. If several patches are simultaneously active and if the corresponding STF and STWV tensors have an approximately trilinear structure, the tensor-based preprocessing separates the patches. The localization is thus performed for each patch separately. Furthermore, we employ sLORETA, cLORETA, and 4-ExSo-MUSIC for extended source localization to compare the results of the tensor-based methods to other approaches. sLORETA and cLORETA are applied to spatially prewhitened data2 as we found that this generally improves the results. Since these two methods do not take into account the temporal information, they are applied to the time sample that exhibits the highest variance over all EEG channels, corresponding to the maximum of the epileptiform spike. The FO statistics employed in the 4-ExSo-MUSIC algorithm are estimated from the same data that is used for the construction of the STF and STWV tensors, comprising N t = 200 time samples. For both DA and 4-ExSo-MUSIC, we consider disks that are composed of up to D max = 100 grid dipoles.

Evaluation

For the comparison of the different extended source localization methods, we use the receiver operating curves (ROC), which represent a plot of the true positive fraction (TPF) over the false positive fraction (FPF) of the estimated grid dipoles. If I denotes the set of extended source dipoles which belong to the original source, Î is the set of correctly identified extended source dipoles, and J represents the set of all grid dipoles, then the two measures are defined as follows:

TPF = #{I ∩ Î} #{I} (20) FPF = #{ Î} -#{I ∩ Î} #{J } -#{I} . ( 21
)
Here, #{I} denotes the number of elements in the set I. Different TPF and FPF values are achieved by varying the threshold values for the extended source localization algorithms. The ROC curves are generally plotted for a FPF ranging from 0 % (no dipoles that are falsely associated to the patch) to 6 %, which corresponds to approximately 60 cm 2 of cortex that is erroneously associated to the patch.

In the case of several patches, the STF and STWV preprocessing separate the patches, yielding a separate ROC curve for each patch. In order to compare these results with the solutions of the other methods, we combine the individual ROC curves by multiplying the FPF by the number of separate patches and computing the arithmetic mean of the TPF values for each FPF. This procedure assumes that the grid dipoles erroneously identified for each source are distinct, such that the FPFs add up. Thus, the resulting ROC curve represent the worst case overall performance for these methods.

Real data

Real EEG data were acquired with a 62-channel measurement system using the common average reference with a sampling rate of 1000 Hz. Our analysis is based on 9 interictal spikes that were selected from the recordings. We have considered data segments comprising the ascending and descending parts of the spikes as well as parts of the following wave. The selected time intervals are marked in Figure 8. A realistic head model was built by segmenting the patient's MRI using the BrainVisa software [START_REF] Rivière | A freely available anatomis/brainVISA package for structural morphometry of the cortical sulci[END_REF]. The lead field matrix was then computed for a cortical mesh with 20003 vertices using Brainstorm [START_REF] Tadel | Brainstorm: A userfriendly application for MEG/EEG analysis[END_REF] and OpenMEEG [START_REF] Gramfort | OpenMEEG: opensource software for quasi static bioelectromagnetics[END_REF][START_REF] Kybic | A common formalism for the integral formulations of the forward EEG problem[END_REF]. In this case, each vertex of the mesh corresponded to one grid dipole.

The sources were localized using STWV-DA, STF-DA, and 4-ExSo-MUSIC as well as using sLORETA and cLORETA. Contrary to STWV-DA, STF-DA, and 4-ExSo-MUSIC, which exploit the data of the whole time interval, sLORETA and cLORETA were applied to three time points corresponding to the first (negative) peak, the second (positive) peak, and the wave. The STWV and STF tensors were constructed and decomposed as described in Section 2.5.2. We analyzed the results obtained for R = 1, R = 2, and R = 3 CP components in accordance with the number of sources that could be expected according to the SEEG recordings that were available for the same patient. For the source localization using STWV-DA, STF-DA, and 4-ExSo-MUSIC, we employed a maximal disk size of 200 dipoles. The number of disks or grid dipoles to consider, which determines the size of the identified patch, was chosen such that the goodness-of-fit (GOF) value

GOF = || X -X rec || F || X|| F (22)
was minimal. Here, X rec corresponds to the data matrix that is reconstructed from the estimated source configuration. In case of STF-DA, STWV-DA, and 4-ExSo-MUSIC, X rec = R r=1 dr sT r where dr denotes the reconstructed spatial mixing vector for the combination of a certain number of disks for the r-th component and sr denotes the corresponding patch signal that can be computed as S = D+ X with D = d1 , . . . , dR . For sLORETA and cLORETA, X rec = ds T where d = d∈ Î g d corresponds to the sum of the lead field vectors of the considered grid dipoles, which are characterized by the set Î of dipole indices and which are identified by thresholding the coefficient vector ĉ. The corresponding patch signal is computed as s = d+ X.

The source localization results can be evaluated based on the findings of the SEEG, which give a strong hypothesis on the actual source regions of the epileptic activity. To this end, in Figure 8, we marked by small spheres the positions of the three SEEG electrodes for which the highest amount of epileptic spikes were automatically detected [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] during SEEG recordings. Note that the automatic detection was based on an independent evaluation of the recordings of each SEEG electrode which means that the epileptic activity at the three identified sites could be independent or concomitant.

A more detailed analysis of the SEEG recordings showed that in some cases, an epileptic spike at the anterior SEEG electrode may be associated with an epileptic spike at the central SEEG electrode, delayed by about 20 ms, and a spike at the posterior SEEG site delayed by 70 ms with respect to the spike at the anterior site. This suggests that epileptic activity is propagated from the anterior SEEG electrode to the posterior SEEG electrode.

Results

Single patch scenarios

Analysis of the patch position

Since superficial sources exhibit more focal distributions of the electric potential than deep sources, this may favor the source localization procedure. Furthermore, the signals emanating from deep sources lead to smaller amplitudes at the sensor level than those of superficial sources and therefore correspond to smaller SNRs for the same background activity. It is thus significant to determine the influence of the patch position on the localization accuracy of the different source localization methods. To this end, in the first simulation, we consider 9 different patches with various locations (see Figure 2) on the left hemisphere, leading to SNR ranging from 3.9 dB to 11 dB where the SNR was computed as SNR = ||X|| 2 ||X b || 2 . We evaluate the performance of the source localization algorithms based on the ROC curves, which are plotted in Figure 2 along with images of the recovered patches for a fixed FPF of approximately 0.2 % and for the source localization method that yields the best results in each case.

As can be seen in Figure 2, for all simulation scenarios other than patch BasTe, STF-DA and STWV-DA outperform sLORETA and cLORETA. For patch BasTe, sLORETA leads to slightly better results than STF-DA and STWV-DA, whereas the latter two clearly yield more accurate results than cLORETA. While the performance of sLORETA and, in particular, cLORETA approaches that achieved by STF-DA and STWV-DA for some scenarios, e.g., patches InfPa and SupOcc, the results of sLORETA and cLORETA are rather poor for other scenarios, e.g., patches MidTe and ParaHipp. Finally, the 4-ExSo-MUSIC algorithm generally yields comparable results to STF-DA and STWV-DA. For the examined single patch scenarios, there is hardly any difference between the results of the STF and the STWV analysis. Concerning the patch location, the simulation results show that the examined source localization methods generally yield good results for superficial patches (except for patches PreC (precentral) and MidTe (mid temporal)). However, they show some difficulties for accurately recovering deep patches such as patches BasTe (basal temporal) and ParaHipp (para hippocampal).

Influence of the patch size

Besides the patch location, another factor that influences the result of the source localization is the patch size. This is why, in the second simulation, we consider a scenario where this parameter is varied over 10, 25, 50, 100, 200, 300, and 400 triangles that composed the patch corresponding to an area of 0.5, 1.25, 2.5, 5, 10, 15, and 20 cm 2 of cortex. Figure 3 shows the resulting patches of different sizes and the TPF score of the examined source localization methods for a fixed FPF of 0.2 %. Except for sLORETA, all source localization methods display the best performance for intermediate patch sizes (around 100 or 200 patch dipoles corresponding to an area of 5 to 10 cm 2 and the green zone in Figure 3 (left)) with decreasing accuracy for both small and very large patches.

Two patch scenarios

One of the main interests of the tensor-based methods consists in their ability to separate several simultaneously active patches. This point is therefore studied in the following simulations. Please note that a discussion of the results obtained for two of the considered scenarios with respect to the theoretical analysis of the STF and STWV analyses can be found in the supplementary material "Trilinear approximation".

Influence of the patch distance

An important factor for the distinction of two patches is their distance, especially for the STWV analysis, which exploits the difference between spatial distributions of the electric potential for each patch. To determine the influence of the patch distance on the source localization results, we consider in the following three configurations of two superficial patches with large, medium, and small distances amounting to approximately 13.5 cm, 9 cm, and 5 cm, respectively. The corresponding ROC curves as well as the original and recovered patches are plotted in Figure 4. For all three patch configurations, STWV-DA clearly outperforms all other approaches. The localization accuracy (in terms of the TPF) diminishes with decreasing source distance. For two close sources, at an FPF of about 0.2 %, only one patch has been partly recovered by STWV-DA while at larger distances, for the same FPF, both patches have been almost completely identified (cf.

Figure 4).

Influence of the depth of the patch

In the previous simulations for two patches, we considered only superficial patches.

However, as already stated in Section 3.1.1, the depth of a patch plays an important role in the outcome of the source localization process. To determine its impact on the source separation and localization of two patches, we conduct a simulation study with three patch The recovered patches and resulting ROC curves are displayed in Figure 5 and show that all tested source localization algorithms have great difficulties in identifying both patches.

Both the STF and the STWV analysis fail to accurately separate the sources. The 4-ExSo-MUSIC algorithm thus features the best performance for the first two scenarios.

However, in the first scenario, it only permits to recover the superficial patch and in the second scenario it first identifies a false patch. In the third scenario involving two deep patches, sLORETA yields the best source localization result. It recovers parts of both patches, but does not permit to identify the true patch forms and extents.

Influence of the number of CP components

The number of CP components identified in the decomposition of the STF and STWV tensors should be chosen according to the number of extended sources. However, in practice, the number of sources is unknown and has to be estimated from the measurements.

While the estimation of the number of sources is out of the scope of the present paper, we analyze in this section the sensitivity of the STF and STWV based source localization methods to the number of CP components used in the tensor decomposition. To this end, we consider a scenario with a single patch, InfPa, and a scenario consisting of the two patches InfFr and InfPa. Then, we decompose the STF and STWV tensors using one CP component and using two CP components. For both cases, we perform source localization using STF-DA and STWV-DA. The resulting ROC curves are shown in Figure 6. For STF-DA, the results that are achieved with one or two CP components are the same, indicating that the spatial mixing vectors obtained for both components must be almost identical. With a TPF close to 100% for an FPF of about 1%, the results obtained by STF-DA are good for the single patch scenario, but poor for two patches where the TPF does not exceed 50% for an FPF smaller than or equal to 6%, which suggests that only one patch is localized. For STWV-DA, with a 1-component CP decomposition, one obtains the same results as with STF-DA for both scenarios. With a 2-component CP decomposition, on the other hand, the results of STWV-DA are worse than those obtained for one component in the single patch case, but considerably better than those obtained with one component in the two patch case. This shows that the correct choice of the number of CP components is important to achieve accurate results with STWV-DA.

Computational Complexity

In Section 2.4, we provided a number of formulas for the computational complexity of the tensor-based preprocessing methods and the source localization algorithms. However, these formulas are difficult to interpret due to their dependence on a large number of parameters. To give an idea of the computational complexity that can be encountered in practice for the compared methods, we compute the number of real-valued multiplications for fixed values of parameters that were used for the computer simulations for two patches presented above and vary only the parameters associated with the data, namely the number of sensors and the number of time samples. As Figure 7 shows, the 4-ExSo-MUSIC algorithm clearly exhibits the highest computational cost, which augments rapidly with increasing number of sensors, but hardly changes for the examined numbers of time samples, since the cost of the EVD of the cumulant matrix depends only on the number of sensors and dominates the cost of all other operations. The cLORETA algorithm has the lowest computational complexity, followed by sLORETA. For both methods, the computational complexity does not depend on the number of time samples since the algorithms are applied to a single data vector. The number of real-valued multiplications required for STF and STWV analyses are comparable for about 200 time samples. For identical tensor sizes, the STWV analysis is computationally more expensive than the STF technique since the STWV tensor is complex-valued. However, for an increasing number of time samples, the calculation of the wavelet transform becomes computationally more expensive, which explains the increased computational burden of STF-DA compared to STWV-DA for large time samples in Figure 7 (right).

Real data results

To validate the tensor-based methods and in particular the STWV analysis, we applied these techniques to real EEG measurements that were recorded for a patient suffering from epilepsy. For comparison, we also analyzed the data using 4-ExSo-MUSIC, sLORETA, and cLORETA. In Figure 8, the patch dipoles that were identified with the examined algorithms for the 9 different spikes are marked in color. The colorscale varies from 0 to 9 depending on the number of spikes for which a dipole was determined to be active.

STWV-DA, STF-DA, 4-ExSo-MUSIC, and sLORETA identified patches mostly on the left hemisphere, at locations that were close to the two or three SEEG contacts that recorded the highest amount of spikes. However, all tested source localization methods also identified patches on the right hemisphere. This was particularly the case for cLORETA. In case of STWV-DA, the patches identified for most of the spikes were located in between the two posterior or the two frontal SEEG contacts, in the superior frontal gyrus or in the mesial areas of the frontal lobe. Some isolated patches were located in the posterior of the right hemisphere. Comparing the results of STWV-DA obtained for R = 1, R = 2, and R = 3 CP components, we observe that for R = 1, only one small patch, located in between the two frontal SEEG contacts, was identified. For R = 2, the recovered patches were larger and located close to the three SEEG contacts, whereas for R = 3, the majority of patches was localized between the two posterior SEEG contacts. The patches localized by STF-DA were mostly located between the two frontal SEEG contacts and were comparable for the three tested tensor ranks. However, for R = 3, in some cases, patches were also identified at the equivalent position on the right hemisphere. Moreover, a small number of patches were localized in the vicinity of the third SEEG contact in the pre-central gyrus. The 4-ExSo-MUSIC algorithm identified patches in between the frontal and the posterior SEEG contact, with most patches located close to the central SEEG contact. The patch size slightly increased for larger ranks of the signal subspace.

Otherwise, the different ranks of the signal subspace lead to similar results. As for STF-DA, patches were also localized at the equivalent positions on the right hemisphere. The patches localized by sLORETA are globally more anterior, while cLORETA identified patches all over the frontal parts of the left and right hemispheres.

Discussion and conclusion

We have compared the results of tensor-based methods to conventional techniques for extended source localization. We reviewed two types of tensor methods, namely STF and STWV analyses, and proposed an algorithm, DA, for source localization based on the results of the tensor decomposition. These methods have been evaluated with respect to three different aspects: (i) their theoretical foundation, (ii) their computational complexity, and (iii) their source localization performance on simulated and real data.

Analysis of the tensor-based methods

Whereas previous references on STF and STWV analyses have provided only intuitive conditions concerning their applicability in practical situations, in the supplementary material of this paper, we have derived sufficient conditions under which the application of the STF and STWV techniques is justified. Although these mathematical conditions are very restrictive and it is difficult to translate them into physiological conditions, which can be verified in practice, we will subsequently discuss several points which can be deduced from the identified conditions and which influence the functioning of the STF and STWV analyses:

Source strengths For a correct separation, the sources should lead to surface measurements of comparable strengths. Sources with significant differences in amplitude or Time-frequency or space-wave-vector characteristics The STF analysis assumes the time-frequency content of each source to be of rank 1 and the STWV analysis is based on the hypothesis of a rank-1 space-wave vector content of each source.

In practice, it is generally sufficient if the singular values of the time-frequency or space-wave-vector matrix of each source decrease quickly.

In fact, as the simulations conducted in Section 3 showed, despite their approximate nature, the tensor-based methods, STF-DA and STWV-DA, lead to good results for a number of scenarios. The usefulness of the STWV method became particularly apparent for several simultaneously active superficial patches, in which case it accurately separated the patches and exhibited a considerable performance gain compared to all other methods. Nevertheless, it is important to retain that even though tensor-based methods can provide a powerful tool for EEG preprocessing and source separation as we demonstrated in Sections 3.1 and 3.2.1 of this paper, they may fail even in ideal situations where there is neither noise nor artifacts and all patch signals are identical. This is particularly the case for the STF analysis when patches emit highly correlated activities (spreading activity from one region to another one) or for the STWV analysis when deep sources are expected (see results of simulations in Section 3.2.2).

Another source of errors for the tensor-based techniques stems from the imperfect synchronization of the signals that are emitted by the dipoles of one patch. Both the STF and the STWV analysis are based on the model ( 4), which approximates model (1) by assuming the same signal for all dipoles within a patch. If the activity of the patch dipoles is not sufficiently synchronous, this model is incorrect and leads to perturbations of the estimated spatial mixing vectors and thereby of the source localization results. The same problem also applies to the 4-ExSo-MUSIC algorithm.

A remaining problem of the tensor methods consists in the estimation of the number of active patches, which we assumed to be known in this paper. As the analysis of STF-DA and STWV-DA with respect to the number of identified CP components (see Section 3.3) showed, this parameter has a high impact on the results of the STWV analysis. An inappropriate number of components may cause the STWV analysis to fail by separating patches into several components or by mixing different patches in one component. For the STF analysis, the number of CP components seemingly did not have an impact on the results. But this insensitivity may be explained by the inability of the STF analysis to identify components that can be associated to different patches because the timefrequency content of the simulated patch activities is nearly identical. In practice, the number of patches has to be determined from the measurements, which is a difficult task, especially in the context of delayed signals for simultaneously active patches. Similarly, for 4-ExSo-MUSIC, one has to estimate the dimension of the FO signal subspace, which raises the same difficulties.

Computational complexity analysis

Most evaluations of source localization algorithms are based exclusively on performance measures that compare the estimated source distributions to the ground truth whereas considerations concerning the computational complexity of different approaches are generally neglected. Nevertheless, the computational cost of a source localization method should also be taken into account as it influences the usefulness of the method in practice. In this paper, we have therefore included a comparison of the different techniques with respect to their computational complexity. Our analysis has revealed that for a few hundred time samples generally used in the STF or STWV analyses, the cost of 4-ExSo-MUSIC is considerably higher than that of STF-DA and STWV-DA. In particular, it increases strongly with the number of sensors, leading to an intractable number of operations for large sensor numbers. On the other hand, the computational cost of 4-ExSo-MUSIC augments only linearly with the number of time samples, which makes it an efficient algorithm for a small number of sensors and a high number (several thousands) of time samples.

The cLORETA algorithm has the smallest computational complexity among the examined source localization methods, followed by sLORETA. Their cost can be mostly associated with the treatment of the lead field matrix (generalized inverse computation) that has to be accomplished prior to the actual source localization. However, this step needs to be performed only once per patient and the actual source localization can then be performed many times at a much smaller cost. Yet as the simulations showed, the cLORETA and sLORETA principles for source localization do not yield as accurate results obtained with simulated data as 4-ExSo-MUSIC, STF-DA, and STWV-DA.

Simulation results

Contrary to previous studies of tensor-based approaches [START_REF] Miwakeichi | Decomposing eeg data into space-time-frequency components using parallel factor analysis[END_REF][START_REF] Morup | Parallel factor analysis as an exploratory tool for wavelet transformed event-related EEG[END_REF]De Vos et al., 2007a,b;[START_REF] Deburchgraeve | Neonatal seizure localization using parafac decomposition[END_REF][START_REF] Becker | Multi-way space-timewave-vector analysis for EEG source separation[END_REF], which have mainly focused on source separation and equivalent dipole localization, we have conducted our analysis using a realistic head model and have aimed at identifying the spatial extent of the sources. The localization of extended sources has also been addressed in several other publications [START_REF] Limpiti | Cortical patch basis model for spatially extended neural activity[END_REF][START_REF] Friston | Multiple sparse priors for the M/EEG inverse problem[END_REF][START_REF] Ding | Reconstructing cortical current density by exploring sparseness in the transform domain[END_REF][START_REF] Bolstad | Space-time event sparse penalization for magneto-/electroencephalography[END_REF][START_REF] Birot | Localisation of extended brain sources from EEG/MEG: the ExSo-MUSIC approach[END_REF][START_REF] Chowdhury | MEG source localization of spatially extended generators for epileptic activity: comparing entropic and hierarchical Bayesian approaches[END_REF], based on different algorithms including scanning methods, MEM, Bayesian and sparse approaches. While these algorithms are generally shown to yield satisfactory results for scenarios with only one extended source, difficulties in localizing several patches, in particular in the case of correlated activities, have been reported in several studies. In this paper, we have analyzed the use of tensor-based methods to overcome this problem by separating the sources prior to the localization. In this context, an important contribution of the paper consists in the proposition of the disk algorithm, which permits us to accurately localize extended sources based on the results of the tensor-based preprocessing step. This method is based on a parameterization of the extended source, similar to [START_REF] Limpiti | Cortical patch basis model for spatially extended neural activity[END_REF][START_REF] Birot | Localisation of extended brain sources from EEG/MEG: the ExSo-MUSIC approach[END_REF], but utilizes a different metric to identify the elements of the source space that best describe the measurements. To evaluate the performance of the proposed STF-DA and STWV-DA techniques in comparison to sLORETA, cLORETA, and 4-ExSo-MUSIC, a simulation setup inspired by [START_REF] Grova | Evaluation of EEG localization methods using realistic simulations of interictal spikes[END_REF] is used.

Our simulation study has shown that for single patch scenarios, 4-ExSo-MUSIC yields comparable results to those of STF-DA and STWV-DA and that these methods outperform sLORETA and cLORETA. As no source separation is performed, the spatial mixing vector is almost perfectly recovered by both the STF and the STWV analysis and these methods yield approximately the same results. In this case, the errors of the patch estimates are mostly due to DA. Similar to previous studies conducted in the context of MEG source localization [START_REF] Chowdhury | MEG source localization of spatially extended generators for epileptic activity: comparing entropic and hierarchical Bayesian approaches[END_REF], we have observed that the extended source localization algorithms do not work as well for deep sources as for superficial sources. This shows that the improvement of these methods for the localization of deep sources is an important research track for future work.

For small patches, the SNR is low, which may explain the poor performance of the source localization algorithms. However, among all considered methods, STF-DA and STWV-DA lead to good results except for very small patches (< 2.5 cm 2 ). On the other hand, for large patches, the limited size of 100 grid dipoles for the disks that we employed in the 4-ExSo-MUSIC technique and the disk algorithm is not adequate anymore and impedes the source localization. This may explain the better performance of cLORETA compared to STF-DA. STWV-DA, and 4-ExSo-MUSIC for large patches and shows that the maximal disk size has to be chosen carefully according to the expected size of the source region to obtain accurate source localization results. To avoid this problem, one could systematically use a high maximal disk size, but it has to be kept in mind that this may considerably increase the computational complexity.

For scenarios with two superficial patches, STWV-DA exhibited the best performance. This can be explained by the fact that the STWV analysis correctly separates the spatial mixing vectors of the two patches and therefore permits to localize each patch individually. Due to the highly correlated signals of the two patches, which differ only by a small time delay, the STF analysis fails to separate the patches, therefore impeding source localization. 4-ExSo-MUSIC needs to localize both patches simultaneously, which does not work as well as the localization of a single patch and thus does not yield as accurate results as STWV-DA. As for the single patch scenarios, cLORETA and sLORETA generally do not permit to achieve as accurate results as STWV-DA. However, if we consider scenarios with one deep patch and one superficial patch or two deep patches, the tensorbased techniques feature a poor performance because they do not accurately separate the two patches. In the presence of one superficial and one deep patch, this is mainly due to the different strengths of the patch signals recorded at the surface. Furthermore, the STF technique fails because of the highly correlated signals of the patches and the STWV method struggles with the wide-spread distribution of the electric potential of deep sources.

Real data

We have applied the STF-DA and STWV-DA algorithms to actual EEG measurements of an epileptic patient and have localized patches that show a good correspondence to the positions of the SEEG electrodes detecting frequent interictal epileptic activity. We employed tensor decompositions with R = 1, R = 2, and R = 3 components because we expected up to three active patches. In all cases, we obtained results that are concordant with the sites identified from the SEEG recordings, yielding patches that are close to two or three of the marked SEEG contacts for most spikes. In some cases, the determined patches can be farther away from the marked SEEG electrodes and include regions on the right hemisphere. This could be due to lower SNRs for the single spikes or to propagation phenomena, which occur during the spike and wave complex of the analyzed epileptic spikes. Nevertheless, it is difficult to consider these results as a "false" localization since, in the absence of simultaneous SEEG/EEG recordings, the involvement of these remote regions cannot be ruled out. For STF-DA, we did not observe significant discrepancies between the results obtained for different tensor ranks. Due to propagation effects, the source signals can be expected to be highly correlated and the STF analysis is therefore unlikely to separate the sources, distinguishing rather different components of one source. This would explain the insensibility of the results to the employed number of CP components, which has also been observed in the simulation study. For STWV-DA, the best results are achieved for R = 2, in which case patches are localized in proximity to the three marked SEEG contacts. For R = 3, the source localization results are slightly less concordant with the sites identified from the SEEG recordings and for R = 1, patches are only identified close to the two frontal SEEG contacts. This suggests that the STWV-DA method is able to separate only two sources. If there is a higher number of active patches, their signals are probably too correlated or their amplitudes too different to enable their separation using the STWV analysis, leading to worse results for higher tensor ranks.

Altogether, we deduce that the tensor-based methods are also well suited for the analysis of real data, as is the 4-ExSo-MUSIC algorithm, for which we obtained similar results.

The results of these three methods are more concordant with the identified SEEG contacts than those of sLORETA and cLORETA, which frequently identified patches on the contralateral hemisphere.

Conclusion

On the whole, we conclude that the disk algorithm applied to the results of the tensorbased preprocessing and in particular of the STWV analysis yields the best results in terms of computational complexity and localization accuracy in many cases. However, the tensor-based methods should be used carefully as they do not provide accurate results under certain conditions such as correlated sources in case of the STF analysis and deep sources for the STWV analysis. Hence their suitability for the analysis of the data at hand should be examined, which can, for example, be accomplished based on prior knowledge about the underlying sources or using other source localization methods for comparison.

  ×Nr is the Tikhonov-regularized inverse matrix of G and I denotes the identity matrix of size N r × N r . Next, the covariance matrix of ĉ takes the form: K = G T (GG T + λI) -1 G. The idea of sLORETA consists in weighting solution ĉ by the inverse of the diagonal matrix (Diag{K}) 1 2 ; see[START_REF] Palmero-Soler | swLORETA: a novel approach to robust source localization and synchronization tomography[END_REF] for

Figure 1 :

 1 Figure 1: Example of highly correlated signals for the dipoles of one patch.

Figure 2 :

 2 Figure 2: ROC curves and recovered patches for the 9 examined single patch scenarios and an FPF of approximately 0.2 % (0.13 % for patch InfPa). Triangles belonging to the original patch are marked in red, correctly identified triangles are dark red and erroneously identified triangles are yellow. L (left), R (right), F (frontal), and O (occipital) indicate the orientations of the figures. Patches BasTe and ParaHipp are viewed from the bottom. The SNR for each scenario is indicated beneath the estimated patch.

Figure 3 :

 3 Figure 3: Different examined patch sizes (left; frontal view) and TPF as a function of the number of patch dipoles for a fixed FPF of 0.2 % (right). The colorbar on the x-axis of the Figure on the right associates the patch dipoles with the triangles of the same color in the Figure on the left. The results of the different methods are coded as follows: Plain blue line -sLORETA, plain beige line with crosses (x) -cLORETA, plain green line with circles -STF-DA, dashed red line with stars -STWV-DA, dashed dotted black line with crosses (+) -4-ExSo-MUSIC.

Figure 4 :

 4 Figure 4: ROC curves and recovered patches for scenarios consisting of patches SupFr and SupOcc (large distance), patches InfFr and InfPa (medium distance), and patches InfFr and SupFr (small distance). Triangles belonging to the original patches are marked in red, correctly identified triangles are dark red and erroneously identified triangles are yellow. L (left), R (right), F (frontal), and O (occipital) indicate the orientations of the figures. Patches SupFr and SupOcc are viewed from the top. The SNR for each scenario is indicated beneath the estimated patch.

Figure 5 :

 5 Figure 5: ROC curves and original and recovered patches for two patch scenarios with varying patch depths: patches InfFr (superficial) and Cing (deep), patches MidTe (superficial) and ParaHipp (deep), and patches BasTe (deep) and ParaHipp (deep). Triangles belonging to the original patches are marked in red, correctly identified triangles are dark red and erroneously identified triangles are yellow. L (left) and R (right) indicate the orientations of the figures. Patches BasTe and ParaHipp are viewed from the bottom. The SNR for each scenario is indicated beneath the estimated patch.

Figure 6 :

 6 Figure 6: ROC curves obtained for STF-DA and STWV-DA based on a CP decomposition with one or two components for a single patch and a two-patch scenario.

Figure 7 :

 7 Figure 7: Computational complexity of the different source localization methods depending on the number of sensors (left) for N t = 200 time samples and as function of the number of time samples (right) for 91 electrodes.

Figure 8 :

 8 Figure 8: Results of STWV-DA and STF-DA for R = 1, R = 2, and R = 3 and of 4-ExSo-MUSIC for P = 1, P = 3, and P = 6 (corresponding to 1, 2, and 3 correlated sources) as well as of sLORETA and cLORETA for the different spikes. The patch dipoles are colored according to the number of spikes (from 0 to 9) for which they were identified. Small blue spheres indicate the positions of the three SEEG electrodes for which the highest amount of epileptic spikes was automatically detected during SEEG recordings. We also illustrated the spike intervals considered for STWV-DA, STF-DA, and 4-ExSo-MUSIC and the time points considered for cLORETA and sLORETA (displayed on the time signal recorded by electrode AF7).

  combinations of deep and superficial sources often lead to the idenfication of the source with the highest surface amplitude only. Correlation of the source time signals The time signals of the different sources should not be too correlated. Low correlations facilitate the source separation. Correlation of the spatial mixing vectors For close sources, the spatial mixing vectors are highly correlated, making the source separation difficult. Distant sources, on the other hand, lead to a limited spatial correlation and favor the source separation.

In equation (15), the Matlab notation is used to denote the i-th row of a matrix.

The spatial prewhitening is accomplished by multiplying the data and the lead field matrix by the prewhitening matrix P = K + , where K is the square root of the covariance matrix of the background activity.

Acknowledgement H. Becker was supported by Conseil Régional PACA and by CNRS France. The work of P. Comon was partly funded by the European Research Council under the EC 7th framework programme FP7/2007-2013 Grant Agreement no. 320594. Furthermore, we acknowledge the support of Programme ANR 2010 BLAN 0309 01 (project MULTI-MODEL).


est destinée au dépôt

Influence of orography and canopy conditions on friction

velocity observed during frontal events using doppler sodar observations

V. Kotroni, Christine Amory-Mazaudier

To cite this version:

V. Kotroni, Christine Amory-Mazaudier. Influence of orography and canopy conditions on friction velocity observed during frontal events using doppler sodar observations. Journal of Applied Meteorology, 1993, 32 (3), pp.1-16. hal-01011946


Jean-Claude Andre

Jean-Paul Goutorbe

Alain Perrier

François Becker

Pierre Bessemoulin

Philippe Bougeault

Yves Brunet

Wilfried Brutsaert

Tobby Carlson

Richard Cuenca

Evaporation over land-surfaces: First results from HAPEX-MOBILHY special observing period

come


Amal Tahri

Laurence Duchien

Jacques Pulou

Jacques Pulou Software Evolution

Software Evolution Multi-View : From the Smart Home to the Cloud

Challenges Motivations

Future Work Ongoing Work



  Choosing FM [1] formalism to manage variability and heterogeneity.  Extending FM by addressing Localization Constraints e.g., Colocation, Dislocation of software components.  Proposing a DSL to express the deployment constraints. Software Evolution Multi-View : From the Smart Home to the Cloud Smart Home characteristics  Heterogeneity: different elements e.g., devices, services and technologies.  Variability: different characteristics e.g., the amount of offered resources for a given element.  Volatility: changing set of applications deployed on a changing set of elements Where to deploy new applications? In the Smart Home? On the Cloud? Or distributed on both?  How to deploy them? Only in one place? By distributing the applications between the Cloud and the Smart Home? 1. Balance the computational load between the Smart Home and the Cloud. 2. Reduce the acquisition cost for users by limiting the embedded systems resources. 3. Reduce the scalability cost when application set grows. Context 1. Optimization of the deployment of a new componentbased application onto a distributed environment (Smart Home + Cloud) despite the heterogeneity and the variability. 2. Run Time self-adaptation of the application deployment to volatility. 1. React to the volatility by run time self-adaption of the application deployment using code offloading and migration between the Smart Home and the Cloud 2. Build a self-* architecture supporting the selfconfiguration, self-optimization and self-adaptation of the application deployment[2]. nom@orange.com; * prenom.nom@inria.fr} [1]Benavides, David, Sergio Segura, and Antonio Ruiz-Cortés, « Automated Analysis of feature models 20 years later : A Literature view » Information System 35.6 (2010) : 615-636. [2]Movahedi, Zeinab et al. « A Survey of Autonomic Network Architectures and Evaluation Criteria » Communications Surveys & Tutorials, IEEE 14-2(2012):464-490.