Numerical resolution of the inverse source problem for EEG using the quasi-reversibility method

The paper concerns the numerical resolution of the inverse source problem for electroencephalography. We propose an approach which is able to take into account some heterogeneity properties (namely a varying electrical conductivity) of the head tissues, in particular of the skull layer. It combines two consecutive steps: (i) a data completion procedure from the scalp to the cortex using the quasi-reversibility method, (ii) a source estimation method from these cortical transmitted data within the brain (modeled as sphere), developed in the software tool FindSources3D. Numerical simulations in the case of the multi-layered spherical head model illustrate both the promising and limiting features of the approach.


Introduction
Electroencephalography (EEG) is a non-invasive cerebral exploration method and one of the most widespread clinical and functional brain imaging techniques.Measurements of the electric potential generated by normal or pathological brain activity are taken at electrodes placed at the surface of the scalp.They record in a passive way the voltage potential fluctuations between different cortical regions.EEG-monitoring can provide periodic and repeated evaluations at the bedside of the patient.The important goal of brain imaging using EEG is to localize cerebral sources within the brain that generate the measured EEG signals on the scalp (together with an appropriate current flux).EEG is known to have an excellent time resolution and is able to record neural events of one millisecond time order.It can be used in presurgical evaluation for refractory epilepsy [7,51].Further, it provides a tool for studying functional brain connectivity, as for instance the development of language networks in children [24].Note that EEG is the single available clinical process able to follow the premature neonate functional brain development [62].
From a mathematical point of view, source localization from incomplete overdetermined Cauchy-type data is an (ill-posed) inverse problem that aims to determine the location of electric brain sources and their characteristics (moments if the sources are dipolar).A large number of analytical and numerical methods for the source reconstruction exist in the literature (see e.g.[17,28], references therein and in the sequel of this section).They are intimately linked with the head model (geometry and electric conductivity distribution), and with the modeling of electric brain sources.The forward problem for EEG consists in finding the potential on the scalp for a given electric source located in the brain.Generally, multilayered head models are considered and the conductivity of each layer (brain, cerebrospinal fluid -CSF, skull, scalp) is assumed to be constant (whence piecewise constant in the overall head, (see [33] and references therein).On the one hand, the spherical head model has gathered much interest from the beginning of EEG source analysis since an explicit expression for the potential is available [49].On the other hand, realistic head models obtained from the segmentation of anatomical images using magnetic resonance imaging (MRI) are able to take into account the geometrical complexity of the different tissues [53].Concerning the source modeling, dipolar or distributed sources can be considered.In equivalent current dipole methods, sources are modeled by a few pointwise dipoles within the brain.They are characterized by their position(s) and moment(s).In distributed source models, the sources are assumed to be located at voxels of the grey/white matter interface and perpendicular to this surface, whence their moments are characterized by their amplitude only.At the frequencies of brain electromagnetic phenomena, the quasi-static approximation of Maxwell equations is valid [31,47] and the (scalar) electric potential solves the conductivity Poisson partial differential equation (PDE) in the head domain.For the numerical resolution of the forward problem, both boundary elements (see e.g.[27,38,50] and references therein) and 3D finite elements (see e.g.[60,61] and references therein) are commonly used.Spherical head models are still frequently considered to solve the forward problem.The accuracy of the EEG source reconstruction strongly depends on the accuracy of the associated forward model.
Source estimation in a domain from data measured at some distance from the source support is a widely addressed inverse problem, in particular for brain imaging purposes.In addition to EEG data on the scalp, available data may consist in measures of the magnetic field outside the head (magnetoencephalography, MEG, see [31]), or measures of the electric potential inside the head (sEEG).These devices record time varying signals, from which static values at a fixed time instant can be obtained, either by direct sampling or performing time-space signal separation through a singular value decomposition (SVD) and component analyses, see e.g.[21].
Depending on the source model, two big classes of methods implemented in software libraries devoted to the resolution of inverse EEG/MEG problems can be distinguished: -Minimum Norm Estimation (MNE) methods, where the source term is assumed to be distributed on the cortex surface (within the brain, which should be known from MRI), its amplitude (moments) being unknown [14,25,29,63].This leads to a number of degrees of freedom much bigger than the number of data.Solutions to this underdetermined linear inverse problem are obtained by norm minimization of the output error combined with a statistical modelization of the source distibution, see ‡ [26,42].
-Dipolar source estimation methods (of dipole fitting family) use a source model reduced to a few pointwise dipoles, make no assumption on their positions, leading to an overdetermined non-linear problem.A large variety of methods exists.Without being exhausive, we cite MUSIC-Multiple Signal Classification [21,44,46,48], beamformers, see [57] and dedicated softwares §, Bayesian methods, see [6,58,59] and references therein.
These different approaches assume a modeling of the head with a piecewise constant conductivity.In the present paper, we propose a numerical method for the EEG inverse source problem which is able to take into account the heterogeneity of the head tissues, in particular of the skull layer, yet under isotropy assumptions (i.e.scalar conductivity).The electric conductivity of the skull can hardly be assumed to be constant because of both hard and spongy bones, especially the presence of fontanels for neonates.Fontanels are tissues that are in the process of ossification.They therefore possess different (inhomogeneous) electric properties in comparison to an homogeneous skull.Furthermore, skull conductivity models vary between individuals [1,3,52] and through time for a single person.The inhomogeneity of the skull conductivity prevents the application of boundary element methods which are currently used for EEG source localization.Finite element method is required to model the skull volume itself.Few reconstruction methods for a varying conductivity are proposed.In [16], a minimization approach using the Broyden-Fletcher-Goldberg-Shanno (BFGS) iterative algorithm is performed with good reconstructions of dipolar sources (position and moment) for realistic geometries.But the sensitivity of the convergence with respect to the initial guess constitutes a disadvantage.Furthermore, the number of sources is assumed to be known a priori.To overcome these drawbacks, we develop a new method which consists of two steps.The approach firstly transmits the recorded data from the scalp to the cortex, ‡ And the dedicated software MNE https://mne.tools/stable.
and then applies a source reconstruction method in the brain.The data transmission issue from the scalp to the cortex amounts to Cauchy problems for Laplace-Poisson PDE whose solutions carry harmonicity properties.They remain ill-posed (additional assumptions are required to ensure existence, uniqueness, and stability of solutions) and have been studied in [5,12,23,31].The resolution of their regularized versions, of Tykhonov type, can be carried out using boundary element methods on the layers interfaces, with meshes as described in [27,47] or appropriate families of functions in simple geometrical settings [2,12,40].In non homogeneous situations that we want to consider, Cauchy transmission issues for the related conductivity PDE have to be solved through the whole corresponding volume, using FEM for realistic geometries.We choose to regularize the Cauchy problems by applying the quasi-reversibility (QR) method.This is a non-iterative approach which was introduced by Lattès and Lions in [39].It has been successfully adopted and validated for Laplace's equation (see e.g.[8,10,11,19,35]).Its variational setting is numerically interesting since finite element methods can be used.We propose to generalize it to a varying conductivity and to study its performance as a data transmission procedure from scalp to brain surface (cortical mapping).This step is central to take into account the heterogeneity of the overall head tissue.We provide theoretical results as well as a numerical study for spherical domains.Different configurations are tested: one single QR method for transmission through several -3 or 4 -homogeneous layers, with or without modeling of fontanels in the skull.We then combine it with a source reconstruction method within the brain.Among the different approaches, we focus on the one developed in the software tool FindSources3D (FS3D) which belongs to the (static) dipolefit class in a spherical representation of the head with piecewise homogeneous conductivity.In particular, the brain is assumed to be a sphere and the conductivity in the brain to be constant.FS3D makes use of an efficient analytic method that allows to retrieve both the number and the characteristics of dipolar sources (locations, moments, even in time correlated situations), see [12].It involves the computation of the spherical harmonic coefficients from the approximated Cauchy data on the cortex.Then, FS3D considers best rational approximation on families of 2-D planar cross-sections in order to locate singularities and to determine the expected quantity of brain electric sources.From those planar singularities, the sources are finally estimated, together with their moment, in a last clustering step.Through this process, FS3D is also able to recover time correlated sources, which is an important advantage with respect to other software tools addressing the same problem.The paper is organized as follows.In Section 2, we describe the forward and inverse problems for EEG.In Section 3, we present the quasi-reversibility method that transmits the data from the electrodes to the brain surface.In Section 4, the algorithm implemented within the software tool FindSources3D (FS3D) is precised.Section 5 is devoted to the numerical validation of both the quasi-reversibility method itself as a cortical mapping procedure and https://www-sop.inria.fr/apics/FindSources3D/en/.
of its coupling with FS3D source localisation method for solving the inverse EEG source method.Various simulations from synthetic data are presented.A conclusion is provided in Section 6.

The forward model
In the low frequency range under consideration in EEG measurements, the electromagnetic field satisfies the quasi-static Maxwell equations where the time derivatives are neglected [31,34].The head tissues are characterized by the electric permittivity ε > 0, magnetic permeability µ > 0 supposed to be constant and a conductivity σ with positive values.Under the assumption that the constitution laws of the medium are linear, the quasi-static Maxwell equations are given by and model the interaction between the electric field E, the magnetic induction B, and the current density J, that are R 3 -valued quantities, while the charge density ρ ∈ R .In the brain and in particular in the cortex, the synchronized effect among a multitude of neurons creates an intracellular current denoted by J s .The total current density J produced by cerebral activity splits into two terms (Ohm's law): the current term in the extracellular region of the brain being described by σE.It follows from ( 1)-(ii), that the electric field derives from a scalar potential u, i.e.E = ∇u .Now, consider a simply connected bounded domain Ω ⊂ R 3 with C 2 smooth boundary ∂Ω.
Taking the divergence of (1)-(iii) together with (2) yields the following transmission equation for the electric potential u in Ω (in the distributional sense, for distributions J s ): In the typical multi-layer head model, we distinguish a number of nested layers corresponding to the brain (containing gray and white matters), maybe to the cerebrospinal fluid (CSF), the skull, and scalp.Therefore, consider a partition of Ω into L open nested subdomains We also assume that Ω 0 coincides with the brain and is simply connected.We denote by Ω L = R 3 \ Ω the outside of Ω and by S i , i = 0, . . ., L − 1, the interface between Ω i and Ω i+1 .We assume that (S i ) i are closed C 2 surfaces.Let n i be the unit normal vector on S i oriented towards the exterior of Ω i .We denote by Γ = ∂Ω the exterior boundary of the whole domain Ω.This configuration includes the classical spherical model of three concentric spheres representing brain, skull and scalp (see Figure 1 left).Left: Spherical head model with L = 3 and one source (c q , p q ).Right: Fontanels and skull of a neonate.
The space W 1,∞ (Ω i ) coincides with the space of Lipschitz continuous functions on Ωi .This assumption on the conductivity is one of the main features of this work.It allows in particular to take into account the specific heterogeneity of the skull layer.Indeed, if the piecewise constant conductivity assumption is relevant for tissues of the brain and scalp (and CSF), it is not the case of the skull, whose conductivity actually varies.For instance, the skull's structure could be modeled by a spongiosa part (soft bone) located between inner and outer compacta parts (hard bone) [52].For the neonates, the presence of fontanels in the skull (see Figure 1 right) could be modeled also through a varying conductivity [3].
The source model of neural activity can be described by a sum of Q electric dipoles located in the brain (e.g.[31,56]).Each dipole is characterized by its position c q ∈ Ω 0 and its moment p q which is a vector of R 3 .The points (c q ) q are assumed to be mutually distinct and the vectors p q are non vanishing, i.e. c q = c s ∀q = s and p q = 0 ∀q, s ∈ {1, . . ., Q}.
The current density J s thus reads: where δ cq denotes the Dirac distribution supported at c q .The r.h.s. of (3) is then given by Since σ L = 0, no electric current can flow out of the head in Ω L and the electric potential u is then solution of the following boundary problem with an homogeneous Neumann condition For given sources (c q , p q ) ∈ Ω 0 × R 3 , q = 1, . . ., Q, and a known distribution of the conductivity σ in Ω, computing a solution u to problem ( 5) is the forward problem for EEG.Notice that (5) includes the following transmission conditions (jumps) on any interface A direct variational formulation of ( 5) in H 1 (Ω) is not possible since the source term F belongs to the Sobolev space H s (R 3 ) for s < −5/2.In order to overcome this difficulty, we apply the subtraction approach [3,23,64].It consists of decomposing the potential u, solution to ( 5), into a potential ũ which contains the singularities and a regular lifting w : With a view to define the singular potentials ũq , we assume that the brain conductivity σ 0 is constant in a neighborhood of each source (see [3] for details).More precisely we fix a family of disjoint open balls (V q ) q=1,...,Q such that V q Ω 0 , c q ∈ V q and σ 0|Vq ≡ α q = σ 0 (c q ) ∈ R, for any q ∈ {1, . . ., Q}. ( The singular potential ũq is solution of the following Poisson equation It can be obtained by convolution of the fundamental solution of the Laplace equation with 1 α q p q • ∇δ cq as follows The potential ũq has a singularity at the source point c q , but is smooth everywhere else.The auxiliary function w satisfies the following boundary value problem Since problem (11) involves a Neumann boundary condition, its solution is determined up to an additive constant only.We therefore introduce the following subspace of H 1 (Ω) The variational formulation for w then reads as follows: and H-elliptic.Furthermore, the linear form l(•) is continuous on H.The compatibility condition l(1) = 0 can be proven with the help of the solid angle formula (see e.g.[37]).Then, Lax-Milgram's theorem guarantees the uniqueness of the solution w ∈ H.We refer to [3,22] for any details.In [4, Prop. 1, Cor. 1], uniqueness and regularity properties of solutions to the forward problem ( 5) are established for smooth enough geometries and piecewise constant conductivities, that may extend to varying conductivities.They also use the above decomposition (7) and the strict ellipticity of the operator, together with the De Giorgi-Nash theorem in order to show that w ∈ H 1 (Ω) is in fact Hölder smooth in Ω.

The inverse source problem
We denote by Γ * ⊂ Γ an open subset of Γ with meas(Γ * ) = 0.The inverse problem for EEG consists in reconstructing the source F of the form (4) from a single Dirichlet measurement f = u | Γ * of a solution u to (5) on Γ * .More precisely, the aim is to determine the quantity Q, the locations (c q ) and moments (p q ), q = 1, . . ., Q, in the source F generating the data f .In experiments, a finite number of surface electrodes, on which voltages are recorded, are attached to the boundary.Consider L well-separated pointwise electrodes {e } L =1 on Γ * .Measurements are given by with The function f has to belong to (the trace space) H 1/2 (Γ * ) and is chosen from the measurements U as an appropriate interpolant of u at points e .Observe that EEG data U are subject to noise and errors, whence the choice of such an f that robustly approximates them is judicious: Concerning the overall inverse problem associated to (5), while the source term appears in divergence form in the right-hand side, its properties very much depend on assumption on the source distribution.In the present situation with pointwise dipoles, the uniqueness of the source term F from available data is established in [23] for a constant conductivity σ.The uniqueness result allows to know if the inverse problem is well posed in the following sense: if two measured potentials coincide on a nonempty subset Γ * , then, they are generated by the same source of the form (4). Stability considerations are also essential properties in view of a numerical resolution.In [16], uniqueness and stability results have been proven with a varying conductivity.The proof of the uniqueness property is based on the subtraction method (following (7)), and the next regularity result for the potential w, Lemma 1 (see [3] for details).
where a ij ∈ W 1,∞ (D),and (a ij ) is a symmetric matrix satisfying the uniform ellipticity condition for some constant λ > 0 i,j Remark 1.In the remarks following [55, Thm 1.4], it is stated that Lemma 2 holds true for u ∈ H 1 (D).Further, if the conductivity is constant in the layer Ω i , then a version of this lemma is available from [41] whenever the coefficients a ij in P do not depend on x.It also asserts that the unique continuation result holds true for u ∈ H 1 (D).
We now precise our procedure to numerically solve the inverse source problem for EEG through two steps (Sections 3, 4).

Data transmission using the quasi-reversibility method
The first step of the reconstruction method is called the cortical mapping procedure.It deals with the numerical resolution of a data transmission problem for elliptic equations.To this end, we propose to apply the quasi-reversibility method (called hereafter QR method) which provides a regularized solution of Cauchy problems in a bounded domain.

Consider the domain
), the Cauchy problem to be solved is the following: find u ∈ H 1 (ω) such that with a varying conductivity σ.The part Γ * is called the accessible part of the boundary Γ.The inaccessible part is given by Γ 1 = (Γ \ Γ * ) ∪ S 0 .The aim is to compute the Cauchy data (u + , ∂ n u + ) on the interior inaccessible part S 0 which is the cortex surface.Then, using the transmission conditions (6), we deduce u The Cauchy problem ( 14) is known to be ill-posed in the sense that it may not possess a solution for arbitrary data (f, g) in H 1/2 (Γ * ) × H −1/2 (Γ * ), and that the solution is not stable with respect to the data [30].On the other hand, if problem (14) admits a solution, it is unique.This is an immediate consequence of unique continuation results for elliptic equations [55] (see Lemma 2 above).

The quasi-reversibility method
The quasi-reversibility method replaces the ill-posed Cauchy problem with a well-posed variational problem including a small regularization parameter.To solve the Cauchy problem (14), we consider the mixed formulation of the quasi-reversibility method which is constructed and studied for Laplace's equation in [8,11] and we adapt it to the presence of a varying conductivity.To this end, let us introduce the spaces We introduce the linear form (•) defined on W by First, let us give the variational formulation of the problem (14).
Proof.Let u ∈ V .Using a test-function v ∈ V , and integrating on ω, it follows from the classical Green's theorem that Hence, if u ∈ V is a solution to (14), it actually belongs to u ∈ V f and satisfies the weak formulation (15 For regularization purposes, the idea is to solve the following weak mixed formulation: where δ is a regularization parameter.The regularized solution (u δ , λ δ ) tends to (u, 0) when δ tends to 0, where u is solution to (14), as stated in the Theorem 2 below.
Proof.We can apply [11,Thm 2.4].Indeed, the trace operator on Γ * , that is Problem ( 17) can be written as a classical variational formulation: find involving the continuous coercive bilinear form A δ defined on (V 0 × W ) 2 by and the linear form δ ((v, µ)) = −δ r, v V − (w) + a(r, w) given on V 0 × W . Existence and uniqueness follow from Lax-Milgram's Lemma for any boundary measurement f ∈ H 1/2 (Γ * ) and for all δ > 0. In addition, if there exists a solution u ∈ V f to the Cauchy problem ( 14), Lemma 1 and 2 ensure that u is unique.The convergence of the sequence (u δ , λ δ ) to (u, 0) follows as in [11] since u is the only solution to the variational formulation (15).
In practice, the data f are noisy and will, in general, not belong to the trace space H 1/2 (Γ * ).The boundary condition u = f on Γ * has thus to be enforced weakly.Assume that f ∈ L 2 (Γ * ).
In this configuration, a relaxed mixed formulation is prescribed in [11].For a varying conductivity, it reads: find where the subscript α = (δ, η) indicates that the solution of ( 18) depends on the regularization parameter δ > 0 and the relaxation parameter η > 0. The following theorem holds, with For any α = (δ, η) such that δ > 0 and η > 0, problem (18) admits a unique solution (u α , λ α ) ∈ V × W .In addition, we have in V × W for any fixed η > 0. Here, u is the unique solution of the Cauchy problem (14).Futhermore, we have the following estimates Proof.Problem ( 18) can be written in the following closed form where the bilinear form A α (•, •) and the linear form L α (•) are defined on V × W by The continuity of A α (•, •) and L α (•) is clearly obtained, and coercivity comes from the relation since the scalar product in W is the one of V .We thus apply Lax-Milgram's Theorem and prove existence and uniqueness of a solution of problem (18).The solution u ∈ V f to the Cauchy problem ( 14) satisfies (see Lemma 3): By subtracting (22) from the second equation of problem (18), we get Now, take v = u α − u in problem (18) and w = λ α in (23) and subtract the latter from the first one.Since we have u = f on Γ * , this yields the following relation which gives From (25), we deduce that u α V ≤ u V , i.e. the sequence (u α ) α is bounded with respect to the two parameters δ and η.We thus get the estimation (20), and the convergence of (λ α ) α to 0 whenever the regularization parameter δ tends to 0. In a same way, we get the estimation (21) on the boundary term.By consequence, the sequence ( = 0, hence for any fixed η > 0. Finally, let us prove the convergence of the sequence (u α ) α .Since this sequence is bounded in V , it admits a subsequence that weakly converges in V to a limit u ∈ V .Passing to the limit in the second equation in (18) yields a(u, w) = (w) , ∀w ∈ W , if δ → 0. On the one hand, we have from the weak convergence of (u α ) α in V .On the other hand, according to (21), Then, from (24) , we get Because (u α ) α weakly converges to u, the above inequality implies the strong convergence, at least for a subsequence.It follows from a standard argument that the whole sequence converges strongly to u.
Remark 4. In Theorem 3, the estimate (20) indicates a convergence order for the convergence of λ α to 0 independently from the choice of η.The estimate (21) suggests that large values of η should accelerate the convergence of the trace of u α on Γ * .
Numerical validation of the data transmission using the QR method will be addressed in Section 5.

Source reconstruction
We have computed approximated Cauchy data on the cortex S 0 , transmitted from the scalp through the skull using the QR method from Section 3. The next stage of the approach consists in solving the inverse source problem within the brain Ω 0 .For simplicity, we stick to the homogeneous brain situation where the conductivity σ 0 is constant within Ω 0 .Following the model described in Section 2.1, we get from ( 7) and (10) (with α q = σ 0 ) that the singular part ũ of the potential u, solution to (4), (5) in Ω 0 , is given by: ũ The inverse source problem in Ω 0 reads: from u and ∂ n u (or σ 0 ∂ n u) on S 0 , estimate the quantity Q ≥ 1, the source locations c q ∈ Ω 0 and their moments p q ∈ R 3 for q = 1, • • • , Q.
Assume that Ω 0 is a ball, hence S 0 = S 0 , the sphere centered at 0, of radius r 0 > 0.
In this spherical setting, we solve the inverse source problem using the following four-steps process, a method which is implemented in the FS3D software ¶ [12].These main bricks are summarized below.
Step 1: Anti-harmonic projection.From u and ∂ n u on S 0 , we get their expansions on the spherical harmonics and that of ũ, see [ , useful information about the locations of sources is obtained from the study of the singularities s d,p,q of f d,p in D d,p , illustrated by Figure 2 (we typically take D = 13, P = 20).
• For fixed p, d, each function f d,p has Q poles s d,p,q in D d,p of multiplicity 3, that are also branched singularities (of order 3/2) if Q ≥ 2; if Q = 1, f d,p is rational with a single triple pole in D d,p (recall that Q is the quantity of dipolar sources c q ).
• For a fixed direction d, the P × Q singularities s d,p,q for p = 1, • • • , P are all contained in Q planes, each of these planes containing the axis ∆ d and the source location c q (whence perpendicular to Π d,p ), in each of these Q planes (for fixed q), the P singularities s d,p,q for p = 1, • • • , P are linked by a curve (line) that passes through the actual source location c q , at maximum distance from the axis ∆ d .
In the following steps, we use these properties in order to estimate the source locations c q and their moments p q .
Step 3: singularities estimation, rational approximation on circles.For fixed d, p, the estimation of the singularities s d,p,q is performed through best quadratic rational approximation of f d,p on T d,p , see [12,Sec. 3.3,3.4].The approximation function takes the rational form π(z)/(z − α) 3 , parameterized by its triple pole α ∈ D d,p and its numerator π, a polynomial of degree less than 2 with complex-valued coefficients.¶ https://www-sop.inria.fr/apics/FindSources3D/en/.The underlying optimization problem concerns 4 complex-valued (or 8 real-valued) parameters and is solved by a standard non-linear minimization algorithm + .The initial guess is estimated either by brute grid evaluation, or by more sophisticated methods relying on the Hankel matrix of the Fourier coefficients of f d,p .Whenever Q = 1, in the single source situation, each function f d,p is rational with one triple pole of degree 3 and this approach allows to recover the triple pole s d,p,1 , at least from "exact" data ũ2 .A nice feature of this process is that it also works out in general for the multiple sources case where Q > 1, which it allows to estimate together with the singularities s d,p,q , at least when the sources are separated enough (at sufficiently different heights) in the slicing direction d.This is illustrated in Figure 2 (a), (b).On the contrary, the rational approximation error furnishes a criterion allowing to reject the badly identified poles.This algorithm is further quite stable and robust to noise on the data.
For a fixed direction d, we then get a collection of triple poles α d,p , p = 1, . . ., P and their arguments arg(α d,p ) ∈ (−π, π].This set is clustered to produce Q estimated angles θ d,q , q = 1, . . ., Q. Together with the direction ∆ d , each angle θ d,q defines a plane that approximately contains the source location c q .In each of these Q clusters, we then look for the maximum modulus ρ d,q of its elements and for the value p which achieves this maximum.Thus ρ d,q and θ d,q define a point ĉd,q which is an estimation of c q (related to the direction d).
Step 4: 3D clustering, final tuning of positions and moments.
Gathering all the D × Q estimations ĉd,q for a family of directions d = 1, . . ., D , form in Ω 0 a collection of points clouds around the Q sources locations c q .Whenever the sources are far enough one from the others in the chosen slicing directions, and the data not too noisy, we indeed obtained D × Q estimations ĉd,q of the source locations, given by the Q line extrema in D directions (Figure 2, (c)).These clouds are separated by a clustering algorithm which returns Q clusters.The respective mean values (barycenter) of these clusters are furnished as initial values for the estimated positions.Once the source locations are fixed, the data depend linearly on the moment values.A first estimation of these moments is obtained by minimizing the L 2norm of the residual error at the electrodes.A final non linear least-squares optimization algorithm + returns both moments and tuned positions (Figure 2, (d)).When true sources are known and their number correctly recovered * , relative position errors are computed.

Numerical results
In this section, we study the performance of the quasi-reversibility (QR) method for data transmission from scalp to cortex together with the source localisation in the cortex using FS3D.We apply a single QR step, with or without accounting for the presence of fontanel in the skull.We consider the three-layer spherical head model (i.e.L = 3, see Figure 1, left).The radii of the brain (Ω 0 ), skull (Ω 1 ) and scalp (Ω 2 ) are respectively 5cm, 5.4cm and 6cm.These sizes correspond to a cranial perimeter of 37.7cm which is approximatively the one of a newborn child.We add a 2mm cerebrospinal fluid (CSF) layer around the cortex in some cases (i.e.L = 4, see Sections 5.6 and 5.7).The new approach is tested for several choices of conductivity distribution σ.The conductivities are assumed to be constant in brain and scalp where σ 0 = σ 2 = 0.33S.m −1 and in CSF where σ CSF = 1.8S.m −1 .In the skull layer Ω 1 , two configurations are modeled and compared: (a) a constant conductivity σ 1 = 0.04S.m −1 , (b) the presence of the anterior fontanel (see Figure 1 right) and its ossification process by defining the following conductivity function We fix σ skull = 0.04S.m −1 and σ f = 0.3S.m −1 for both the neonatal skull conductivity and the fontanel conductivity [54].The function g, which is defined by g(x) = e −s(x 2 1 +x 2 2 ) , models the process of fontanel ossification.In the numerical experiments hereafter, s is set to 10 4 and thus the fontanel is limited to the region Ω f = {x ∈ Ω 1 : x 2 1 + x 2 2 ≤ r 2 , x 3 > 0} with r = 10mm (see Figure 14).

Generation of synthetic data
Synthetic data are generated by solving the forward problem (5) (see Section 2.1) for the neonatal head models (a) and (b).The subtraction method is implemented with FreeFem++ (see [32]).Discretization is performed with 3D Finite Elements of type P2.Error estimates are proved in [3].We use a tetrahedral mesh of the three-layer spherical domain Ω, composed of 617 167 elements.
The measurement region Γ * on scalp corresponds to the union of small patches around electrode positions.Three sets of electrode positions with 33, 61 and 257 electrodes are examined (see Figure 3).In practice, the mesh is refined around the electrode positions using the Mmg remeshing software [15].The surface mesh of the scalp also explicitly includes the boundary of small circular patches representing the electrodes around each electrode position, thanks to the level-set meshing capabilities of Mmg on triangular surface meshes (see Figure 4, left).Then, the synthetic EEG measurement on each electrode patch on the scalp is taken as the average value of the potential on the electrode patch (see Figure 4).The goal is to emulate the features of a real electrode measurement (single value measurement), while avoiding the inaccuracy of interpolating a pointwise measurement on a regular nonadapted finite element mesh.
In the following we will consider two sets of symmetric sources: • two distant sources c q d = (2.4,−0.8, 2.9) and c q −d = (−2.4,−0.8, 2.9).The respective moments of the sources are p q d = (0.9, 0.2, 0.3)J and p q −d = (−0.9,0.2, 0.3)J, with intensity J = 4.1e −5 A.m −2 .These two sources are far away from each other and from the fontanel, and lie at 1.15cm below the brain surface S 0 .
Numerical experiments in Sections 5.2 and 5.3 will be performed with one distant source (c q d , p q d ).We will use the two distant sources (c q d , p q d ) and (c q −d , p q −d ) in Section 5.4.Section 5.6 will present results with one close source (c qc , p qc ), while Sections 5.5 and 5.7 cover the two close sources case (c qc , p qc ) and (c q −c , p q −c ).Such configurations are relevant for EEG, since the electric activity for a specific task is generally activated on one or two regions of the brain, and here modeled by pointwise dipolar sources.

Numerical validation of the quasi-reversibility method
This section is dedicated to the numerical study of the data transmission step using the QR method (18).Convergence of the discretized solution to the exact solution is proved in [8,9] for Laplace's equation and can be generalized to a varying conductivity.One goal is to assess numerically the effect of the regularization and the relaxation parameters δ and η on the quality of the reconstruction of the potential for different noise levels.We propose a calibration of the parameters with respect to a piecewise constant conductivity.It consists of a QR step for L = 3 embedded homogeneous layers.In the noiseless case, we take δ = 10 −5 and we set η = 10 15 , which corresponds to an approach with no relaxation (equivalent to (17)).Here we consider one distant source (c q d , p q d ).The source term of the EEG forward problem ( 5) is thus given by F = p q d • ∇δ cq d and the data are generated in the way explained in Section 5.1.Numerical solutions of the mixed formulation (18) are obtained using P2 Lagrange finite elements.We implement the quasi-reversibility formulations with FreeFem++ (see [32]).We compare the reconstructed potential, denoted u + h,α , using the quasi-reversibility method with the one u + h computed with the forward EEG solver.We measure the error in the reconstructed Dirichlet and Neumann data at the cortex S 0 as .
First, the accessible part Γ * is taken as the entire "upper" half of the scalp S 2 where z > 0 instead of the electrode patches.On Figures 5 and 6 we report the relative L 2 errors e D,S 0 and e N,S 0 for δ varying between 10 −6 and 10 −3 and η varying between 10 −2 and 10 3 , for 5% and 20% noise respectively.For 5% noise, the smallest reconstruction error is achieved for δ = 10 −5 and η = 1, with e D,S 0 = 0.11 and e N,S 0 = 0.221.For 20% noise, we get the smallest reconstruction errors, e D,S 0 = 0.142 and e N,S 0 = 0.292, again for δ = 10 −5 and η = 1.Comparatively, in the noiseless case, the reconstruction error is below 10%.We get e D,S 0 = 0.041 and e N,S 0 = 0.084.Now, we consider the configuration with 257 electrodes (see Figure 3).Again, we use synthetic data (namely f = u h| Γ * and g = 0) and add 5% and 20% white Gaussian noise to the value of the potential at each electrode.We aim at finding the best (δ, η) combination for each of the two noise levels.Figures 7 and 8 show the relative L 2 errors e D,S 0 and e N,S 0 for 5% and 20% noise respectively.For 5% noise, the smallest reconstruction error is achieved for δ = 10 −4 and η = 100, with e D,S 0 = 0.217 and e N,S 0 = 0.418.For 20% noise, we get the smallest reconstruction errors, e D,S 0 = 0.298 and e N,S 0 = 0.537, for δ = 10 −4 and η = 10.Comparatively, in the noiseless case, the reconstruction error is e D,S 0 = 0.125 and N,S 0 = 0.265.These values of the regularization and relaxation parameters will be kept throughout the rest of the numerical experiments for noisy synthetic electrode data.As expected, the data transmission is less efficient when the noise level increases.Nevertheless, we will see that the quasi-reversibility method allows to keep the valuable information on the data at the cortex which are needed for the source reconstruction.

Two sources, no fontanel
Here we consider the symmetric distant sources (c d , p q d ) and (c −d , p q −d ).QR reconstructions on the cortex and source reconstruction with FS3D for noiseless data and noisy data with 5% and 20% white Gaussian noise on the value of the potential at each electrode are reported in Figure 12 for 257 electrodes.Here and throughout the rest of the paper, an asterisk ( * ) on the position error means that the position error corresponds to the best one-to-one pairings of true and reconstructed sources, even though the number of reconstructed sources is larger than the number of true sources.The QR method well captures qualitatively the two localized perturbations of the electric potential at the cortex corresponding to the two sources.Results are better for the Dirichlet trace of the reconstructed electric potential than for its Neumann trace.The coupling QR-FS3D provides very precise source localizations, even with noisy data.

Two close sources, no fontanel
We test now the challenging case of two symmetric close sources (c qc , p qc ) and (c q −c , p q −c ).
The data are computed for the case of 257 electrodes.QR reconstructions on the cortex and source reconstruction with FS3D for noiseless data and noisy data are shown in Figure 13.
Results are again concluding and the two sources are efficiently localized.Their orientations are also well retrieved.Even if the reconstruction errors after the data transmission are higher than for two distant sources, we can still distinguish the presence of two distinct electric activities at the cortex, which are associated with the two sources.The distinction is less easy for 20% noise.

One source, CSF fontanel
The above results attest efficiency of the method for a piecewise constant conductivity.
The is that the transmission step is realized in one step using a QR formulation, and not layer by layer as other approaches.Let us study now the effect of an heteregeneous skull layer.For these simulations, we consider the head model (b) including the main fontanel in the skull layer and a varying conductivity (27).We add also a 2mm cerebrospinal fluid (CSF) layer around the cortex.The CSF layer has high conductivity σ CSF = 1.8S.m −1 .Figure 14 shows the mesh of such a domain.We perform the same previous tests with one close source (c qc , p qc ) located close to the brain-CSF interface.The results are reported in Figures 15, 16 and validate the approach in the case of a more realistic conductivity.We compare what happens when not taking the fontanel into account in the QR reconstruction (columns 1-3), and when taking the fontanel into account (columns 4-6).In Figure 15, we consider measures furnished by 61 electrodes on the scalp, a more practical situation [45].There, taking the fontanel into account gives slightly better results.In Figure 16, with 257 electrodes, the accuracy difference between the two columns is less obvious in presence of noise, though the results remain pretty good.The data transmission method is able to treat complex head models and is complementary with the reconstruction approach of FS3D.Here again, at the cortex surface, the QR method keeps the location of the support of the Dirichlet and Neumann traces which is useful for solving the inverse source problem.

Two close sources, CSF and fontanel
Finally, we test the reconstruction of the two symmetric close sources (c qc , p qc ) and (c q −c , p q −c ) using synthetic data generated with the head model (b) and CSF. Figure 17 (61 electrodes) and Figure 18 (257 electrodes) compare the case without (columns 1-3) and with (columns 4-6) the fontanel in the QR reconstruction.The configuration with 257 electrodes provides the best results.For both 61 and 257 electrodes cases, the sources are efficiently reconstructed until 5% noise level.In the noiseless case, the presence of the fontanel leads to slightly more accurate reconstructions, since no parasite source is detected as in the head model without fontanel.For noisy cases, results are quite similar with or without the anterior fontanel.Nevertheless, visually, it seems that the information associated with the two close sources is better reconstructed by the QR process with the presence of the fontanel, until 5% noise level.Fontanel would play a role in the localization of close sources for neonates from EEG measurements.We summarize all the results of the section in Table 1.The method allows to treat a varying skull conductivity and is validated for such an heterogeneous medium.

Conclusion
This work presents results concerning EEG inverse source problems for spherical inhomogeneous head models, obtained from the joint use of two numerical methods, QR and FS3D.The QR method propagates external measures (taken at electrodes on the scalp) through inhomogeneous layers to an internal spherical surface, the boundary of a domain which contains unknown sources.The Cauchy data obtained on that surface are processed by FS3D to estimate dipolar sources.FS3D takes full advantage of the source locations being contained in a spherical domain.The QR transmission step could be coupled with other source reconstruction approaches.
A numerical study was performed for preliminary configurations (spherical head models, noisy synthetic data, a particular heterogeneity in the skull) which shows both the advantages of the method and the needed improvements.In particular, Neumann data on the cortex are not as accurately reconstructed as the Dirichlet trace in Section 5, which slightly affects the source localization in presence of noise and for several sources.Nevertheless, QR satisfactorily performs the transmission process, in one single step through the various tissues (scalp, skull, CSF).provides an interesting homogeneous or heteregeneous multilayer transmission algorithm which is easy to implement.Such an additional transmission step will also be analysed in order to handle more realistic head geometries.Indeed, the QR method itself can provide data transmission in meshed complex geometries.Mixed geometries the head, with general layers and an additional spherical domain included in the brain and enclosing unknown sources could be considered.
Another development will be processing of time varying records of EEG/MEG signal.FS3D is yet equiped for time/space separation of asynchronous sources.
A further interesting perspective would be to propose numerical (iterative) methods for a simultaneous reconstruction of tissue conductivities and electric brain sources.It would allow for patient-specific EEG [52].

Figure 1 .
Figure 1.Left: Spherical head model with L = 3 and one source (c q , p q ).Right: Fontanels and skull of a neonate.

Lemma 2 .
Let D ⊂ R 3 be a smooth bounded connected open domain.Let γ be a nonempty open set of the boundary ∂D of D. Consider the linear operator Consequently, we get u = f on Γ * and u is a solution of the weak Cauchy problem.The uniqueness of the solution yields u = u.

Figure 5 .
Figure 5. Relative L 2 error the QR reconstruction in the Dirichlet (left) and Neumann (right) data at the cortex in the half sphere case with 5% noise, for varying values of δ and η.

Figure 6 .
Figure 6.Relative L 2 error of the QR reconstruction in the Dirichlet (left) and Neumann (right) data at the cortex in the half sphere case with 20% noise, for varying values of δ and η.

Figure 7 .
Figure 7. Relative L 2 error of the QR reconstruction in the Dirichlet (left) and Neumann (right) data at the cortex in the 257 electrodes case with 5% noise, for varying values of δ and η.

Figure 8 .Figure 9 .Figure 10 .Figure 11 .
Figure 8. Relative L 2 error of the QR reconstruction in the Dirichlet (left) and Neumann (right) data at the cortex in the 257 electrodes case with 20% noise, for varying values of δ and η.

eFigure 12 .Figure 13 .
Figure 12.Two distant sources (c q d , p q d ) and (c q −d , p q −d ), head model (a): QR reconstruction on the cortex with 257 electrodes and source reconstruction with FS3D.From top to bottom: exact solution, then reconstruction with 0%, 5% and 20% noise.From left to right: Dirichlet trace, Neumann trace, source reconstruction with FS3D.The asterisk ( * ) means too many reconstructed sources, and position error with the best true/reconstructed source pairing is shown.

Figure 14 .Figure 15 .
Figure 14.Mesh of the three conductivity regions around the cortex corresponding to CSF, skull and scalp, including the main fontanel region with varying conductivity.

eFigure 16 .
Figure 16.One close source (c qc , p qc ), model (b) with CSF and fontanel: QR reconstruction on the cortex with 257 electrodes and source reconstruction with FS3D.Columns 1-3: without fontanel in the QR reconstruction; Columns 4-6: taking the fontanel into account.From top to bottom: exact solution, then reconstruction with 0%, 5% and 20% noise.From left to right: Dirichlet trace, Neumann trace, source reconstruction with FS3D.

eFigure 17 .Figure 18 .
Figure 17.Two close sources (c qc , p qc ) and (c q−c , p q−c ), model (b) with CSF and fontanel: QR reconstruction on the cortex with 61 electrodes and source reconstruction with FS3D.Columns 1-3: without fontanel in the QR reconstruction; Columns 4-6: taking the fontanel into account.From top to bottom: exact solution, then reconstruction with 0%, 5% and 20% noise.From left to right: Dirichlet trace, Neumann trace, source reconstruction with FS3D.The asterisk ( * ) means too many reconstructed sources, and position errors with the best true/reconstructed source pairings are shown.
Sources: configuration of the sources (see Section 5.1).Head model : constant conductivity model (a) or presence of the anterior fontanel (b) and CSF.For the rows "same w/o font.",data is generated with model (b) and CSF, but the anterior fontanel is not taken into account in the QR reconstruction step.#E : number of electrodes.e D,S0 , e N,S0 : error in the reconstructed Dirichlet and Neumann data at the cortex S 0 .e pos : relative error on the source position.The asterisk ( * ) means too many reconstructed sources, and position errors with the best true/reconstructed source pairings are shown.
12, Sec.3.1.2].Choose two positive integers D and P ; choose D directions (axes) ∆ d for d = 1, • • • , D; for each direction ∆ d , choose a family of P parallel slicing planes Π d,p for [12,1, • • • , P , perpendicular to ∆ d .This defines a family of circles T d,p contained in S 0 ∩Π d,p , boundaries of disk D d,p contained in Ω 0 ∩ Π d,p .On the circle T d,p , get the trace (restriction) of ũ2 , hence its complex-valued extension f d,p to the complex plane C.From(26), and[12, Sec.3.2]