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

Marion Darbas, Juliette Leblond, Jean-Paul Marmorat, Pierre-Henri Tournier

► To cite this version:

Marion Darbas, Juliette Leblond, Jean-Paul Marmorat, Pierre-Henri Tournier. Numerical resolution of the inverse source problem for EEG using the quasi-reversibility method. Inverse Problems, 2023, 39 (11), pp.115003. 10.1088/1361-6420/acf9c6. hal-03880526v2

HAL Id: hal-03880526

https://inria.hal.science/hal-03880526v2

Submitted on 4 Sep 2023

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Numerical resolution of the inverse source problem for EEG using the quasi-reversibility method.

M. Darbas<fnref n="1" />, J. Leblond<fnref n="2" />, J.-P. Marmorat<fnref n="3" />, and P.-H. Tournier<fnref n="4" />

E-mail: darbas@math.univ-paris13.fr

**Abstract.** The paper concerns the numerical resolution of the inverse source problem for electroencephalography (EEG). 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 (QR) method, (ii) a source estimation method from these cortical transmitted data within the brain (modeled as sphere), developed in the software tool FindSources3D (FS3D). Numerical simulations in the case of the multi-layered spherical head model illustrate both the promising and limiting features of the approach.

*Keywords:* Inverse problems, elliptic PDE, quasi-reversibility method, data transmission, source estimation, brain imaging (EEG).

Submitted to: *Inverse Problems*

## 1. 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

<footnote n="1">¹ LAGA CNRS UMR 7539, Institut Galilée, 99 avenue Jean-Baptiste Clément, 93430 Villetaneuse, France</footnote>

<footnote n="2">² Team FACTAS, Centre Inria d'Université Côte d'Azur, BP 93, 06902 Sophia Antipolis Cedex, France</footnote>

<footnote n="3">³ CMA MinesParis, 1 rue Claude Daunesse, CS 10 207, 06904 Sophia Antipolis Cedex, France</footnote>

<footnote n="4">⁴ Sorbonne Université, CNRS, Université Paris Cité, Inria, Laboratoire Jacques-Louis Lions (LJLL), F-75005 Paris, France</footnote>
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, multi-
layered 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 distribution, 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 exhaustive, 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,

<footnote>‡ And the dedicated software MNE https://mne.tools/stable.</footnote>

<footnote>§ The software EEGLAB (signal processing features), https://sccn.ucsd.edu/eeglab; see also Brainstorm (source estimation using weighted minimum-norm estimates, beamforming), https://neuroimage.usc.edu/brainstorm, FieldTrip, https://www.fieldtriptoolbox.org.</footnote>
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.

## 2. The forward and inverse problems for EEG

### 2.1. 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 $\epsilon > 0$, magnetic permeability $\mu > 0$ supposed to be constant and a conductivity $\sigma$ with positive values. Under the assumption that the constitution laws of the medium are linear, the quasi-static Maxwell equations are given by

$$ (i) \nabla \cdot (\varepsilon \mathbf{E}) = \rho, \quad (ii) \operatorname{curl} \mathbf{E} = 0, \quad (iii) \operatorname{curl} \mathbf{H} = \mathbf{J}, \quad (iv) \nabla \cdot (\mu \mathbf{H}) = 0, \quad (1) $$

and model the interaction between the electric field $\mathbf{E}$, the magnetic induction $\mathbf{B}$, and the current density $\mathbf{J}$, that are $\mathbb{R}^3$-valued quantities, while the charge density $\rho \in \mathbb{R}$. In the brain and in particular in the cortex, the synchronized effect among a multitude of neurons creates an intracellular current denoted by $\mathbf{J}^s$. The total current density $\mathbf{J}$ produced by cerebral activity splits into two terms (Ohm's law):

$$ \mathbf{J} = \sigma \mathbf{E} + \mathbf{J}^s, \quad (2) $$

the current term in the extracellular region of the brain being described by $\sigma \mathbf{E}$. It follows from (1)-(ii), that the electric field derives from a scalar potential $u$, i.e. $\mathbf{E} = \nabla u$.

Now, consider a simply connected bounded domain $\Omega \subset \mathbb{R}^3$ with $C^2$ smooth boundary $\partial \Omega$. Taking the divergence of (1)-(iii) together with (2) yields the following transmission equation for the electric potential $u$ in $\Omega$ (in the distributional sense, for distributions $\mathbf{J}^s$):

$$ -\nabla \cdot (\sigma \nabla u) = \nabla \cdot \mathbf{J}^s. \quad (3) $$

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 $\Omega$ into $L$ open nested subdomains $(\Omega_i)_{i=0,...,L-1}$ such that

$$ \bar{\Omega} = \bigcup_{i=0}^{L-1} \bar{\Omega}_i \text{ while } \Omega_i \cap \Omega_j = \emptyset, \forall i \neq j, \text{ and } \bar{\Omega}_i \cap \bar{\Omega}_j = \emptyset, \forall j \notin \{i-1, i, i+1\}. $$

We also assume that $\Omega_0$ coincides with the brain and is simply connected. We denote by $\Omega_L = \mathbb{R}^3 \setminus \bar{\Omega}$ the outside of $\Omega$ and by $S_i, i = 0, ..., L-1$, the interface between $\Omega_i$ and
$\Omega_{i+1}$. We assume that $(S_i)_i$ are closed $\mathcal{C}^2$ surfaces. Let $\mathbf{n}_i$ be the unit normal vector on $S_i$ oriented towards the exterior of $\Omega_i$. We denote by $\Gamma = \partial\Omega$ the exterior boundary of the whole domain $\Omega$. This configuration includes the classical spherical model of three concentric spheres representing brain, skull and scalp (see Figure 1 left).

<figure><img src="image_2.png" /><figcaption>**Figure 1.** Left: Spherical head model with $L = 3$ and one source $(\mathbf{c}_q, \mathbf{p}_q)$. Right: Fontanels and skull of a neonate.</figcaption></figure>

Let $\sigma_i = \sigma_{|\Omega_i}$ ($i = 0, \dots, L$) denote the conductivity of the subdomain $\Omega_i$. We set $\sigma_L = 0$ in $\Omega_L$, since the conductivity of the air is zero. We assume that $\sigma_i$ belongs to $W^{1,\infty}(\Omega_i)$ such that

$$0 < c \le \sigma_i(\boldsymbol{x}) \le C < +\infty, \quad \boldsymbol{x} \in \Omega_i, \quad i = 0, \dots, L-1.$$

The space $W^{1,\infty}(\Omega_i)$ coincides with the space of Lipschitz continuous functions on $\bar{\Omega}_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 $\mathbf{c}_q \in \Omega_0$ and its moment $\mathbf{p}_q$ which is a vector of $\mathbb{R}^3$. The points $(\mathbf{c}_q)_q$ are assumed to be mutually distinct and the vectors $\mathbf{p}_q$ are non vanishing, i.e.

$$\mathbf{c}_q \neq \mathbf{c}_s \quad \forall q \neq s \quad \text{and} \quad \mathbf{p}_q \neq 0 \quad \forall q, s \in \{1, \dots, Q\}.$$

The current density $\mathbf{J}^s$ thus reads:

$$\mathbf{J}^s = \sum_{q=1}^{Q} \mathbf{p}_q \delta_{c_q}$$
where $\delta_{c_q}$ denotes the Dirac distribution supported at $c_q$. The r.h.s. of (3) is then given by

$$F = \nabla \cdot \mathbf{J}^s = \sum_{q=1}^{Q} \mathbf{p}_q \cdot \nabla \delta_{c_q}. \qquad (4)$$

Since $\sigma_L = 0$, no electric current can flow out of the head in $\Omega_L$ and the electric potential $u$ is then solution of the following boundary problem with an homogeneous Neumann condition

$$-\nabla \cdot (\sigma \nabla u) = F \text{ in } \Omega, \sigma \partial_n u = 0 \text{ on } \Gamma. \qquad (5)$$

For given sources $(c_q, p_q) \in \Omega_0 \times \mathbb{R}^3$, $q = 1, \dots, Q$, and a known distribution of the conductivity $\sigma$ in $\Omega$, 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 $S_i$, $i = 0, \dots, L-1$,

$$0 = [u]|_{S_i} = u_{|S_i}^{-} - u_{|S_i}^{+}, \quad 0 = [\sigma \partial_n u]|_{S_i} = \sigma_i \partial_n u_{|S_i}^{-} - \sigma_{i+1} \partial_n u_{|S_i}^{+}. \qquad (6)$$

A direct variational formulation of (5) in $H^1(\Omega)$ is not possible since the source term $F$ belongs to the Sobolev space $H^s(\mathbb{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 $\tilde{u}$ which contains the singularities and a regular lifting $w$:

$$u = \tilde{u} + w, \text{ with } \tilde{u} = \sum_{q=1}^{Q} \tilde{u}_q. \qquad (7)$$

With a view to define the singular potentials $\tilde{u}_q$, we assume that the brain conductivity $\sigma_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,\dots,Q}$ such that $V_q \subsetneq \Omega_0$, $c_q \in V_q$ and

$$\sigma_{0|V_q} \equiv \alpha_q = \sigma_0(c_q) \in \mathbb{R}, \text{ for any } q \in \{1, \dots, Q\}. \qquad (8)$$

The singular potential $\tilde{u}_q$ is solution of the following Poisson equation

$$-\alpha_q \Delta \tilde{u}_q = \mathbf{p}_q \cdot \nabla \delta_{c_q} \text{ in } \mathbb{R}^3. \qquad (9)$$

It can be obtained by convolution of the fundamental solution of the Laplace equation with $\frac{1}{\alpha_q} \mathbf{p}_q \cdot \nabla \delta_{c_q}$ as follows

$$\tilde{u}_q(\mathbf{x}) = \frac{1}{4\pi\alpha_q} \mathbf{p}_q \cdot \frac{(\mathbf{x} - \mathbf{c}_q)}{|\mathbf{x} - \mathbf{c}_q|^3}, \quad \forall \mathbf{x} \in \mathbb{R}^3 \setminus \{\mathbf{c}_q\}. \qquad (10)$$
The potential $\tilde{u}_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

$$ \left\{ \begin{array}{l} -\nabla \cdot (\sigma \nabla w) = \sum_{q=1}^{Q} \nabla \cdot ((\sigma - \alpha_q) \nabla \tilde{u}_q) \quad \text{in } \Omega, \\ \\ \sigma \partial_n w = -\sigma \sum_{q=1}^{Q} \partial_n \tilde{u}_q \quad \text{on } \Gamma. \end{array} \right. \qquad (11) $$

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(\Omega)$

$$ H = \left\{ v \in H^1(\Omega) \mid \int_{\Omega} v \, dx = 0 \right\}. $$

The variational formulation for $w$ then reads as follows:

$$ \text{find } w \in H \text{ such that } b(w, v) = l(v) \quad \forall v \in H^1(\Omega), \qquad (12) $$

where

$$ b(w, v) = \int_{\Omega} \sigma \nabla w \cdot \nabla v \, dx $$

is a symmetric linear form on $H^1(\Omega) \times H^1(\Omega)$ and

$$ l(v) = \sum_{q=1}^{Q} \int_{\Omega} (\alpha_q - \sigma) \nabla \tilde{u}_q \cdot \nabla v \, dx - \sum_{q=1}^{Q} \int_{\Gamma} \alpha_q \partial_n \tilde{u}_q v \, ds $$

is a linear form on $H^1(\Omega)$. The bilinear form $b(\cdot, \cdot)$ is clearly continuous on $H^1(\Omega) \times H^1(\Omega)$ and $H$-elliptic. Furthermore, the linear form $l(\cdot)$ 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 \in 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 \in H^1(\Omega)$ is in fact Hölder smooth in $\Omega$.

## 2.2. The inverse source problem

We denote by $\Gamma^* \subset \Gamma$ an open subset of $\Gamma$ with meas($\Gamma^*) \neq 0$. The inverse problem for EEG consists in reconstructing the source $F$ of the form (4) from a single Dirichlet measurement $f = u|_{\Gamma^*}$ of a solution $u$ to (5) on $\Gamma^*$. More precisely, the aim is to determine the quantity
$Q$, the locations ($c_q$) and moments ($p_q$), $q = 1, \dots, 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 $\mathcal{L}$ well-separated pointwise electrodes $\{e_\ell\}_{\ell=1}^L$ on $\Gamma^*$. Measurements are given by

$$u(e_\ell) = f(e_\ell) = U_\ell, \quad \ell = 1, \dots, \mathcal{L}, \qquad (13)$$

with $U = (U_\ell)_{\ell=1}^L \in \mathbb{R}^L$. The function $f$ has to belong to (the trace space) $H^{1/2}(\Gamma^*)$ and is chosen from the measurements $U_\ell$ as an appropriate interpolant of $u$ at points $e_\ell$. Observe that EEG data $U_\ell$ are subject to noise and errors, whence the choice of such an $f$ that robustly approximates them is judicious:

$$f(e_\ell) \approx U_\ell, \quad |f(e_\ell) - U_\ell| \to 0 \text{ when the error level } \to 0, \quad \ell = 1, \dots, \mathcal{L}.$$

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 $\sigma$. 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 $\Gamma^*$, 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).

**Lemma 1.** Assume that $S_i$ is of class $C^2$ for any $i \in \{0, \dots, L-1\}$. Let $\sigma \in L^\infty(\Omega)$ be such that $0 < \sigma_{\min} \le \sigma(\mathbf{x}) \le \sigma_{\max}$ for almost any $\mathbf{x} \in \Omega$, where $\sigma_{\min}$ and $\sigma_{\max}$ are two given positive constants. In addition, assume $\sigma_i \in W^{1,\infty}(\Omega_i)$, $i = 0, \dots, L-1$, satisfying condition (8). Let $w \in H^1(\Omega)$ be the solution of the variational problem (12). Then we have $w|_{\Omega_i} \in H^2(\Omega_i)$ for any $i = \{0, \dots, L-1\}$.

The proof of the uniqueness also involves the unique continuation result for elliptic equations recalled in Lemma 2 (from [55, Thm 1.4]) with $a_{ii} = \sigma$ and $a_{ij} = 0$ for $i \ne j$.

**Lemma 2.** Let $D \subset \mathbb{R}^3$ be a smooth bounded connected open domain. Let $\gamma$ be a nonempty open set of the boundary $\partial D$ of $D$. Consider the linear operator

$$Pu = \sum_{i,j=1}^{3} \frac{\partial}{\partial x_i} a_{ij} \frac{\partial u}{\partial x_j},$$
where $a_{ij} \in W^{1,\infty}(D)$, and $(a_{ij})$ is a symmetric matrix satisfying the uniform ellipticity condition for some constant $\lambda > 0$

$$ \sum_{i,j} a_{ij}(x) \xi_i \xi_j \geq \lambda |\xi|^2 \quad \forall x \in D, \forall \xi \in \mathbb{R}^3. $$

If $u \in H^2(D)$ satisfies $Pu = 0$ in $D$ and $u|_{\gamma} = \partial_n u|_{\gamma} = 0$, then $u = 0$ in $D$.

**Remark 1.** In the remarks following [55, Thm 1.4], it is stated that Lemma 2 holds true for $u \in H^1(D)$. Further, if the conductivity is constant in the layer $\Omega_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 \in H^1(D)$.

We now precise our procedure to numerically solve the inverse source problem for EEG through two steps (Sections 3, 4).

### 3. 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.

#### 3.1. Presentation of the Cauchy problem

Consider the domain $\omega := \cup_{i=1}^{L-1} \Omega_i = \Omega \setminus \bar{\Omega}_0$. From boundary data $(f,g) \in H^{1/2}(\Gamma^*) \times H^{-1/2}(\Gamma^*)$, the Cauchy problem to be solved is the following: find $u \in H^1(\omega)$ such that

$$ \left\{ \begin{array}{ll} \nabla \cdot (\sigma \nabla u) = 0 & \text{in } \omega, \\ u = f & \text{on } \Gamma^*, \\ \sigma \partial_n u = g & \text{on } \Gamma^*, \end{array} \right. \qquad (14) $$

with a varying conductivity $\sigma$. The part $\Gamma^*$ is called the accessible part of the boundary $\Gamma$. The inaccessible part is given by $\Gamma_1 = (\Gamma \setminus \bar{\Gamma^*}) \cup S_0$. The aim is to compute the Cauchy data $(u^+, \partial_n u^+)$ on the interior inaccessible part $S_0$ which is the cortex surface. Then, using the transmission conditions (6), we deduce $u_{|S_0}^- = u_{|S_0}^+$ and $\sigma_0 \partial_n u_{|S_0}^- = \sigma_1 \partial_n u_{|S_0}^+$.

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}(\Gamma^*) \times H^{-1/2}(\Gamma^*)$, 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).
### 3.2. 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

$$V = H^1(\omega), V_f = \{u \in V; u = f \text{ on } \Gamma^*\}, \text{ and } W = \{\mu \in V; \mu = 0 \text{ on } \Gamma_1\}.$$

Let $a(\cdot, \cdot)$ be the bilinear form involved in (15) defined on $V \times V$ by

$$a(u, v) := \int_{\omega} \sigma \nabla u \cdot \nabla v d\boldsymbol{x}, \text{ and let } \langle u, v \rangle_V = \int_{\omega} \nabla u \cdot \nabla v d\boldsymbol{x}.$$

We introduce the linear form $\ell(\cdot)$ defined on $W$ by

$$\ell(v) = \int_{\Gamma^*} g v_{|_{\Gamma^*}} d\boldsymbol{x}.$$

First, let us give the variational formulation of the problem (14).

**Lemma 3.** Assume that $\sigma$ satisfies the hypothesis of Lemma 1. For $f \in H^{1/2}(\Gamma^*)$, if the function $u \in V$ is a solution to (14), then $u \in V_f$ and for all $v \in W$,

$$a(u, v) = \ell(v). \qquad (15)$$

*Proof.* Let $u \in V$. Using a test-function $v \in V$, and integrating on $\omega$, it follows from the classical Green's theorem that

$$\int_{\omega} \nabla \cdot (\sigma \nabla u) v d\boldsymbol{x} = -a(u, v) + \int_{\partial \omega} \sigma \partial_{n} u v ds. \qquad (16)$$

Hence, if $u \in V$ is a solution to (14), it actually belongs to $u \in V_f$ and satisfies the weak formulation (15) since $\partial\omega = \Gamma^* \cup \Gamma_1$, $\sigma\partial_n u|_{\Gamma^*} = g$ and $v \in W$ ($v|_{\Gamma_1} = 0$). $\square$

For regularization purposes, the idea is to solve the following weak mixed formulation: for small $\delta > 0$, find $(u_\delta, \lambda_\delta) \in V_f \times W$ such that

$$\left\{ \begin{array}{ll} \delta \langle u_\delta, v \rangle_V + a(v, \lambda_\delta) & = 0, \quad \forall v \in V_0, \\ a(u_\delta, w) - \langle \lambda_\delta, v \rangle_V & = \ell(w), \quad \forall w \in W, \end{array} \right. \qquad (17)$$

where $\delta$ is a regularization parameter. The regularized solution $(u_\delta, \lambda_\delta)$ tends to $(u, 0)$ when $\delta$ tends to 0, where $u$ is solution to (14), as stated in the Theorem 2 below.
**Theorem 2.** Let $f \in H^{1/2}(\Gamma^*)$. For any $\delta > 0$, the mixed formulation (17) admits a unique solution $(u_\delta, \lambda_\delta) \in V_f \times W$. In addition, we have

$$ \lim_{\delta \to 0} (u_\delta, \lambda_\delta) = (u, 0) $$

in $V_f \times W$, where $u \in V$ is the unique solution of the initial Cauchy problem (14).

*Proof.* We can apply [11, Thm 2.4]. Indeed, the trace operator on $\Gamma^*$, that is $u \in V \mapsto u|_{\Gamma^*} \in H^{1/2}(\Gamma^*)$, is onto. Let us introduce some $r \in V$ such that $r|_{\Gamma^*} = f$. Let us set $\tilde{u}_\delta = u_\delta - r$. Problem (17) can be written as a classical variational formulation: find $(\tilde{u}_\delta, \lambda_\delta) \in V_0 \times W$ such that for all $(v, w) \in V_0 \times W$

$$ A_\delta((u, \lambda), (v, w)) = \ell_\delta((v, w)) $$

involving the continuous coercive bilinear form $A_\delta$ defined on $(V_0 \times W)^2$ by

$$ A_\delta((u, \lambda), (v, w)) = \delta\langle u, v \rangle_V + a(v, \lambda) - a(u, w) + \langle \lambda, w \rangle_V, $$

and the linear form $\ell_\delta((v, \mu)) = -\delta\langle r, v \rangle_V - \ell(w) + a(r, w)$ given on $V_0 \times W$. Existence and uniqueness follow from Lax-Milgram's Lemma for any boundary measurement $f \in H^{1/2}(\Gamma^*)$ and for all $\delta > 0$. In addition, if there exists a solution $u \in V_f$ to the Cauchy problem (14), Lemma 1 and 2 ensure that $u$ is unique. The convergence of the sequence $(u_\delta, \lambda_\delta)$ to $(u, 0)$ follows as in [11] since $u$ is the only solution to the variational formulation (15). $\square$

In practice, the data $f$ are noisy and will, in general, not belong to the trace space $H^{1/2}(\Gamma^*)$. The boundary condition $u = f$ on $\Gamma^*$ has thus to be enforced weakly. Assume that $f \in L^2(\Gamma^*)$.

In this configuration, a relaxed mixed formulation is prescribed in [11]. For a varying conductivity, it reads: find $(u_\alpha, \lambda_\alpha) \in V \times W$ such that

$$ \left\{ \begin{array}{ll} \delta\langle u_\alpha, v \rangle_V + \eta^2 \langle u_{\alpha|_{\Gamma^*}}, v|_{\Gamma^*} \rangle_{L^2(\Gamma^*)} + a(v, \lambda_\alpha) & = \eta^2 \langle f, v|_{\Gamma^*} \rangle_{L^2(\Gamma^*)}, \quad \forall v \in V, \\ a(u_\alpha, w) - \langle \lambda_\alpha, w \rangle_V & = \ell(w), \quad \forall w \in W, \end{array} \right. \quad (18) $$

where the subscript $\alpha = (\delta, \eta)$ indicates that the solution of (18) depends on the regularization parameter $\delta > 0$ and the relaxation parameter $\eta > 0$. The following theorem holds, with $\| \cdot \|_V = \sqrt{\langle \cdot, \cdot \rangle_V}$.

**Theorem 3.** Let $f \in L^2(\Gamma^*)$. For any $\alpha = (\delta, \eta)$ such that $\delta > 0$ and $\eta > 0$, problem (18) admits a unique solution $(u_\alpha, \lambda_\alpha) \in V \times W$. In addition, we have

$$ \lim_{\delta \to 0} (u_\alpha, \lambda_\alpha) = (u, 0) \quad (19) $$
in $V \times W$ for any fixed $\eta > 0$. Here, $u$ is the unique solution of the Cauchy problem (14). Furthermore, we have the following estimates

$$ \|\lambda_{\alpha}\|_{V} \leq \sqrt{\delta}\|u\|_{V} \quad \forall \eta > 0, \qquad (20) $$

$$ \|u_{\alpha|_{\Gamma^*}} - f\|_{L^2(\Gamma^*)} \leq \frac{\sqrt{\delta}}{\eta} \|u\|_V. \qquad (21) $$

*Proof.* Problem (18) can be written in the following closed form

$$ \left\{ \begin{array}{l} \text{Find } (u_{\alpha}, \lambda_{\alpha}) \in V \times W \text{ such that} \\ A_{\alpha}((u_{\alpha}, \lambda_{\alpha}), (v, w)) = L_{\alpha}((v, w)) \quad \forall (v, w) \in V \times W, \end{array} \right. $$

where the bilinear form $A_{\alpha}(\cdot, \cdot)$ and the linear form $L_{\alpha}(\cdot)$ are defined on $V \times W$ by

$$ A_{\alpha}((u, \lambda), (v, w)) = \delta\langle u, v \rangle_V + \eta^2 \langle u_{|_{\Gamma^*}}, v_{|_{\Gamma^*}} \rangle_{L^2(\Gamma^*)} + a(v, \lambda) - a(u, w) + \langle \lambda, w \rangle_V, $$

and

$$ L_{\alpha}((v, w)) = \eta^2 \langle f, v_{|_{\Gamma^*}} \rangle_{L^2(\Gamma^*)} - \ell(w). $$

The continuity of $A_{\alpha}(\cdot, \cdot)$ and $L_{\alpha}(\cdot)$ is clearly obtained, and coercivity comes from the relation

$$ A_{\alpha}((u, \lambda), (u, \lambda)) = \delta\|u\|_{V}^{2} + \eta^{2}\|u|_{\Gamma^{*}}\|_{L^{2}(\Gamma^{*})} + \|\lambda\|_{V}^{2} \geq \min(\delta, 1)\|(u, \lambda)\|_{V \times W}^{2} $$

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 \in V_f$ to the Cauchy problem (14) satisfies (see Lemma 3):

$$ a(u, w) = \ell(w), \quad \forall w \in W. \qquad (22) $$

By subtracting (22) from the second equation of problem (18), we get

$$ a(u_{\alpha} - u, w) - \langle \lambda_{\alpha}, w \rangle_{V} = 0, \quad \forall w \in W. \qquad (23) $$

Now, take $v = u_{\alpha} - u$ in problem (18) and $w = \lambda_{\alpha}$ in (23) and subtract the latter from the first one. Since we have $u = f$ on $\Gamma^*$, this yields the following relation

$$ \delta(u_{\alpha}, u_{\alpha} - u)_{V} + \eta^{2} \|u_{\alpha|_{\Gamma^{*}}} - f\|_{L^{2}(\Gamma^{*})}^{2} + \|\lambda_{\alpha}\|_{V}^{2} = 0 \qquad (24) $$

which gives

$$ \delta \|u_{\alpha}\|_{V}^{2} + \eta^{2} \|u_{\alpha|_{\Gamma^{*}}} - f\|_{L^{2}(\Gamma^{*})}^{2} + \|\lambda_{\alpha}\|_{V}^{2} \leq \delta \|u_{\alpha}\|_{V} \|u\|_{V}. \qquad (25) $$

From (25), we deduce that $\|u_{\alpha}\|_{V} \leq \|u\|_{V}$, i.e. the sequence $(u_{\alpha})_{\alpha}$ is bounded with respect to the two parameters $\delta$ and $\eta$. We thus get the estimation (20), and the convergence of $(\lambda_{\alpha})_{\alpha}$ to 0 whenever the regularization parameter $\delta$ tends to 0. In a same way, we get the estimation (21) on the boundary term. By consequence, the sequence $(u_{\alpha|_{\Gamma^*}})_{\alpha}$ tends to $f$
if $\lim_{\delta \to 0} (\frac{\sqrt{\delta}}{\eta}) = 0$, hence for any fixed $\eta > 0$. Finally, let us prove the convergence of the sequence $(u_\alpha)_\alpha$. Since this sequence is bounded in $V$, it admits a subsequence that weakly converges in $V$ to a limit $\bar{u} \in V$. Passing to the limit in the second equation in (18) yields

$$a(\bar{u}, w) = \ell(w), \quad \forall w \in W,$$

if $\delta \to 0$. On the one hand, we have

$$\langle u_{\alpha|_{\Gamma^*}}, \varphi \rangle_{L^2(\Gamma^*)} \to \langle \bar{u}_{|_{\Gamma^*}}, \varphi \rangle_{L^2(\Gamma^*)}, \quad \forall \varphi \in L^2(\Gamma^*),$$

from the weak convergence of $(u_\alpha)_\alpha$ in $V$. On the other hand, according to (21), $u_{\alpha|_{\Gamma^*}} \to f$ strongly in $L^2(\Gamma^*)$ when $\lim_{\delta \to 0} (\frac{\sqrt{\delta}}{\eta}) = 0$. Consequently, we get $\bar{u} = f$ on $\Gamma^*$ and $\bar{u}$ is a solution of the weak Cauchy problem. The uniqueness of the solution yields $\bar{u} = u$.

Then, from (24), we get

$$\begin{aligned}
\|u_{\alpha} - u\|_{V}^{2} &= \langle u_{\alpha}, u_{\alpha} - u \rangle_{V} - \langle u, u_{\alpha} - u \rangle_{V} \\
&= -\frac{\eta^{2}}{\delta} \|u_{\alpha|_{\Gamma^{*}}} - f\|_{L^{2}(\Gamma^{*})} - \frac{1}{\delta} \|\lambda_{\alpha}\|_{V}^{2} - \langle u, u_{\alpha} - u \rangle_{V} \\
&\leq -\langle u, u_{\alpha} \rangle_{V} + \|u\|_{V}^{2}.
\end{aligned}$$

Because $(u_\alpha)_\alpha$ 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$. $\square$

*Remark 4.* In Theorem 3, the estimate (20) indicates a convergence order for the convergence of $\lambda_\alpha$ to 0 independently from the choice of $\eta$. The estimate (21) suggests that large values of $\eta$ should accelerate the convergence of the trace of $u_\alpha$ on $\Gamma^*$.

Numerical validation of the data transmission using the QR method will be addressed in Section 5.

## 4. 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 $\Omega_0$. For simplicity, we stick to the homogeneous brain situation where the conductivity $\sigma_0$ is constant within $\Omega_0$. Following the model described in Section 2.1, we get from (7) and (10) (with $\alpha_q = \sigma_0$) that the singular part $\tilde{u}$ of the potential $u$, solution to (4), (5) in $\Omega_0$, is given by:

$$\tilde{u}(\boldsymbol{x}) = \frac{1}{4\pi\sigma_0} \sum_{q=1}^{Q} \boldsymbol{p}_q \cdot \frac{(\boldsymbol{x} - \boldsymbol{c}_q)}{|\boldsymbol{x} - \boldsymbol{c}_q|^3}, \quad \forall \boldsymbol{x} \in \mathbb{R}^3 \setminus \{\boldsymbol{c}_q\}. \qquad (26)$$
The inverse source problem in $\Omega_0$ reads: from $u$ and $\partial_n u$ (or $\sigma_0 \partial_n u$) on $S_0$, estimate the quantity $Q \ge 1$, the source locations $c_q \in \Omega_0$ and their moments $p_q \in \mathbb{R}^3$ for $q = 1, \dots, Q$. Assume that $\Omega_0$ is a ball, hence $S_0 = \mathbb{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 $\partial_n u$ on $S_0$, we get their expansions on the spherical harmonics and that of $\tilde{u}$, see [12, Sec. 3.1.2].

**Step 2: 2D slicing.** Choose two positive integers $D$ and $P$; choose $D$ directions (axes) $\Delta_d$ for $d = 1, \dots, D$; for each direction $\Delta_d$, choose a family of $P$ parallel *slicing* planes $\Pi_{d,p}$ for $p = 1, \dots, P$, perpendicular to $\Delta_d$. This defines a family of circles $T_{d,p}$ contained in $S_0 \cap \Pi_{d,p}$, boundaries of disk $D_{d,p}$ contained in $\Omega_0 \cap \Pi_{d,p}$. On the circle $T_{d,p}$, get the trace (restriction) of $\tilde{u}^2$, hence its complex-valued extension $f_{d,p}$ to the complex plane $\mathbb{C}$.

From (26), and [12, Sec. 3.2], 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 \ge 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 \times Q$ singularities $s_{d,p,q}$ for $p = 1, \dots, P$ are all contained in $Q$ planes, each of these planes containing the axis $\Delta_d$ and the source location $c_q$ (whence perpendicular to $\Pi_{d,p}$),

- in each of these $Q$ planes (for fixed $q$), the $P$ singularities $s_{d,p,q}$ for $p = 1, \dots, P$ are linked by a curve (line) that passes through the actual source location $c_q$, at maximum distance from the axis $\Delta_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 $\pi(z)/(z-\alpha)^3$, parameterized by its triple pole $\alpha \in D_{d,p}$ and its numerator $\pi$, a polynomial of degree less than 2 with complex-valued coefficients.

¶ https://www-sop.inria.fr/apics/FindSources3D/en/.
<figure><img src="image_3.png" /><figcaption>**Figure 2.** From left to right. (a), (b), Steps 2, 3: one direction `((0, x₃))`, *Q* = 2 sources (blue), theoretical singularities (green), estimated poles (red); (c), Step 3: theoretical singularity lines cross at sources locations; (d), Step 4: clustering of estimated positions (red), last least-square locations and moments tuning (blue).</figcaption></figure>

The underlying optimization problem concerns 4 complex-valued (or 8 real-valued) parameters and is solved by a standard non-linear minimization algorithm<fnref n="1" />. 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 $\tilde{u}^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 $\alpha_{d,p}$, $p = 1, \dots, P$ and their arguments $\arg(\alpha_{d,p}) \in (-\pi, \pi]$. This set is clustered to produce *Q* estimated angles $\theta_{d,q}$, $q = 1, \dots, Q$. Together with the direction $\Delta_d$, each angle $\theta_{d,q}$ defines a plane that approximately contains the source location $\mathbf{c}_q$. In each of these *Q* clusters, we then look for the maximum modulus $\rho_{d,q}$ of its elements and for the value *p* which achieves this maximum. Thus $\rho_{d,q}$ and $\theta_{d,q}$ define a point $\hat{\mathbf{c}}_{d,q}$ which is an estimation of $\mathbf{c}_q$ (related to the direction *d*).

*Step 4: 3D clustering, final tuning of positions and moments.* Gathering all the $D \times Q$ estimations $\hat{\mathbf{c}}_{d,q}$ for a family of directions $d = 1, \dots, D$, form in $\Omega_0$ a collection of points clouds around the *Q* sources locations $\mathbf{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 \times Q$ estimations $\hat{\mathbf{c}}_{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

<footnote n="1">¹ Matlab optimization toolboxes: `fmincon`, `lsqlin`.</footnote>
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²-norm 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.

## 5. 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 ($\Omega_0$), skull ($\Omega_1$) and scalp ($\Omega_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 $\sigma$. The conductivities are assumed to be constant in brain and scalp where $\sigma_0 = \sigma_2 = 0.33S.m^{-1}$ and in CSF where $\sigma_{CSF} = 1.8S.m^{-1}$. In the skull layer $\Omega_1$, two configurations are modeled and compared: (a) a constant conductivity $\sigma_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

$$ \sigma_1(\mathbf{x}) = \sigma_{skull} + (\sigma_f - \sigma_{skull})g(\mathbf{x})\mathbf{1}_{\Omega_f}(\mathbf{x}), \quad \mathbf{x} = (x_1, x_2, x_3) \in \Omega_1. \quad (27) $$

We fix $\sigma_{skull} = 0.04S.m^{-1}$ and $\sigma_f = 0.3S.m^{-1}$ for both the neonatal skull conductivity and the fontanel conductivity [54]. The function $g$, which is defined by $g(\mathbf{x}) = e^{-s(x_1^2+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 $\Omega_f = \{\mathbf{x} \in \Omega_1 : x_1^2+x_2^2 \le r^2, x_3 > 0\}$ with $r = 10\text{mm}$ (see Figure 14).

### 5.1. 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 $\Omega$, composed of 617 167 elements.

<footnote>* This may not always be the case for several sources or noisy data; see Figures 12, 13, 17, 18 and Table 1.</footnote>
The measurement region Γ* on the 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 non-adapted finite element mesh.

In the following we will consider two sets of symmetric sources:

* two *distant* sources $\mathbf{c}_{qd} = (2.4, -0.8, 2.9)$ and $\mathbf{c}_{q-d} = (-2.4, -0.8, 2.9)$. The respective moments of the sources are $\mathbf{p}_{qd} = (0.9, 0.2, 0.3)J$ and $\mathbf{p}_{q-d} = (-0.9, 0.2, 0.3)J$, with intensity $J = 4.1e^{-5}\text{A}\cdot\text{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$.

* two *close* sources $\mathbf{c}_{qc} = (1, -1, 4)$ and $\mathbf{c}_{q-c} = (-1, -1, 4)$ of respective moments $\mathbf{p}_{qc} = (1.2, 1.5, 3.8)J$ and $\mathbf{p}_{q-c} = (-1.2, 1.5, 3.8)J$, with intensity $J = 1e^{-5}\text{A}\cdot\text{m}^{-2}$. These two sources are close to each other and lie below the fontanel, at 0.76cm below $S_0$.

Numerical experiments in Sections 5.2 and 5.3 will be performed with one distant source ($\mathbf{c}_{qd}, \mathbf{p}_{qd}$). We will use the two distant sources ($\mathbf{c}_{qd}, \mathbf{p}_{qd}$) and ($\mathbf{c}_{q-d}, \mathbf{p}_{q-d}$) in Section 5.4. Section 5.6 will present results with one close source ($\mathbf{c}_{qc}, \mathbf{p}_{qc}$), while Sections 5.5 and 5.7 cover the two close sources case ($\mathbf{c}_{qc}, \mathbf{p}_{qc}$) and ($\mathbf{c}_{q-c}, \mathbf{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.

<figure><img src="image_4.png" /><figcaption>**Figure 3.** Electrode positions (side view). From left to right: 33, 61 and 257 electrodes test cases.</figcaption></figure>
<figure><img src="image_5.png" /><figcaption>Figure 4. From left to right: mesh adapted to the electrode patches, exact synthetic data on scalp, average value of the potential on each electrode patch used as input data for QR.</figcaption></figure>

## 5.2. 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 $\delta$ and $\eta$ 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 $\delta = 10^{-5}$ and we set $\eta = 10^{15}$, which corresponds to an approach with no relaxation (equivalent to (17)). Here we consider one distant source ($c_{qd}, p_{qd}$). The source term of the EEG forward problem (5) is thus given by $F = p_{qd} \cdot \nabla \delta_{qd}$ 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,\alpha}^{+}$, 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

$$e_{D,S_0} = \frac{\|u_{h,\alpha|S_0}^{+} - u_{|S_0}^{+}\|_{L^2(S_0)}}{\|u_{|S_0}^{+}\|_{L^2(S_0)}}, \quad e_{N,S_0} = \frac{\|\partial_n u_{\alpha|S_0}^{+} - \partial_n u_{|S_0}^{+}\|_{L^2(S_0)}}{\|\partial_n u_{|S_0}^{+}\|_{L^2(S_0)}}.$$

First, the accessible part $\Gamma^{*}$ 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 $\delta$ varying between $10^{-6}$ and $10^{-3}$ and $\eta$ varying between $10^{-2}$ and $10^{3}$, for 5% and 20% noise respectively. For 5% noise, the smallest reconstruction error is achieved for $\delta = 10^{-5}$ and $\eta = 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 $\delta = 10^{-5}$ and $\eta = 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$.
<figure><img src="image_6.png" /><figcaption>**Figure 5.** 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 5% noise, for varying values of $\delta$ and $\eta$.</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>**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 $\delta$ and $\eta$.</figcaption></figure>

Now, we consider the configuration with 257 electrodes (see Figure 3). Again, we use synthetic data (namely $f = u_{h|_{\Gamma^*}}$ 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 $(\delta, \eta)$ 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 $\delta = 10^{-4}$ and $\eta = 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 $\delta = 10^{-4}$ and $\eta = 10$. Comparatively, in the noiseless case, the reconstruction error is $e_{D,S_0} = 0.125$ and
$e_{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.

<figure><img src="image_8.png" /><figcaption>**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 $\delta$ and $\eta$.</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>**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 $\delta$ and $\eta$.</figcaption></figure>
5.3. One source, no fontanel

We study the whole process in the case of the head model (a) and one distant dipolar source ($c_{qd}, p_{qd}$) as in the previous section. QR reconstructions on the cortex along with source reconstruction with FS3D for noiseless data and noisy data with 5% and 20% noise on the value of the potential at each electrode are shown in Figures 9, 10 and 11 for 33, 61 and 257 electrodes respectively. We report only the relative error on the source position $e_{pos}$. The first row of each figure (here and in the remaining of Section 5) corresponds to the reference data simulated by solving the forward problem and the true source localisation. The characteristics of the source (location and orientation) are very accurately retrieved in each configuration. Best results are obtained for 257 electrodes. The impact of the electric brain source is visible on a specific area of the cortex (first row). Such a behaviour of the Cauchy data is still observed after the transmission process. Thus, knowing this information, the source reconstruction is efficient.
<figure><img src="image_10.png" /><figcaption>**Figure 9.** One distant source ($c_{q_d}, p_{q_d}$) and head model (a): QR reconstruction on the cortex with 33 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.</figcaption></figure>
<figure><img src="image_11.png" /><figcaption>**Figure 10.** One distant source ($c_{q_d}, p_{q_d}$) and head model (a): QR reconstruction on the cortex with 61 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.</figcaption></figure>
<figure><img src="image_12.png" /><figcaption>**Figure 11.** One distant source ($c_{q_d}, p_{q_d}$) and 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.</figcaption></figure>
### 5.4. Two distant sources, no fontanel

Here we consider the two symmetric distant sources ($c_{q_d}, p_{q_d}$) and ($c_{q_{-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.

### 5.5. Two close sources, no fontanel

We test now the challenging case of two symmetric close sources ($c_{q_c}, p_{q_c}$) 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.
<figure><img src="image_13.png" /><figcaption>**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.</figcaption></figure>
<figure><img src="image_14.png" /><figcaption>**Figure 13.** Two close sources ($c_{q_c}, p_{q_c}$) and ($c_{q_{-c}}, p_{q_{-c}}$) and 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 errors with the best true/reconstructed source pairings are shown.</figcaption></figure>
5.6. One source, CSF and fontanel

The above results attest the efficiency of the method for a piecewise constant conductivity.
The advantage is that the transmission step is easily realized in one step using a QR
formulation, and not layer by layer as other approaches.

Let us study now the effect of an heterogeneous 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 σ<sub>CSF</sub> = 1.8S.m<sup>-1</sup>. Figure 14 shows the mesh of such a
domain. We perform the same previous tests with one close source (c<sub>qc</sub>, p<sub>qc</sub>) 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.

<figure><img src="image_15.png" /><figcaption>**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.</figcaption></figure>
<figure><img src="image_16.png" /><figcaption>**Figure 15.** One close source ($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.</figcaption></figure>

5.7. *Two close sources, CSF and fontanel*

Finally, we test the reconstruction of the two symmetric close sources ($c_{q_c}, p_{q_c}$) 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
<figure><img src="image_17.png" /><figcaption>**Figure 16.** One close source ($c_{q_c}, p_{q_c}$), 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.</figcaption></figure>

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.
<figure><img src="image_18.png" /><figcaption>**Figure 17.** Two close sources ($c_{q_c}, p_{q_c}$) 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.</figcaption></figure>
<figure><img src="image_19.png" /><figcaption>**Figure 18.** Two close sources ($c_{q_c}, p_{q_c}$) and ($c_{q_{-c}}, p_{q_{-c}}$), 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. The asterisk (*) means too many reconstructed sources, and position errors with the best true/reconstructed source pairings are shown.</figcaption></figure>
<table><caption>**Table 1.** Summary of reconstruction errors for each case. *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,S_0}$, $e_{N,S_0}$: 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.</caption><thead><tr><th rowspan="2">Sources</th><th rowspan="2">Head model</th><th rowspan="2">#E</th><th colspan="3">0% noise</th><th colspan="3">5% noise</th><th colspan="3">20% noise</th></tr><tr><th>e<sub>D,S<sub>0</sub></sub></th><th>e<sub>N,S<sub>0</sub></sub></th><th>e<sub>pos</sub></th><th>e<sub>D,S<sub>0</sub></sub></th><th>e<sub>N,S<sub>0</sub></sub></th><th>e<sub>pos</sub></th><th>e<sub>D,S<sub>0</sub></sub></th><th>e<sub>N,S<sub>0</sub></sub></th><th>e<sub>pos</sub></th></tr></thead><tbody><tr><td>1 distant</td><td>(a)</td><td>33</td><td>0.387</td><td>0.638</td><td>0.0064</td><td>0.428</td><td>0.68</td><td>0.0068</td><td>0.578</td><td>0.817</td><td>0.011</td></tr><tr><td>1 distant</td><td>(a)</td><td>61</td><td>0.245</td><td>0.487</td><td>0.0026</td><td>0.281</td><td>0.512</td><td>0.0039</td><td>0.356</td><td>0.59</td><td>0.0057</td></tr><tr><td>1 distant</td><td>(a)</td><td>257</td><td>0.125</td><td>0.265</td><td>0.0011</td><td>0.217</td><td>0.418</td><td>0.0023</td><td>0.298</td><td>0.537</td><td>0.0031</td></tr><tr><td>2 distant</td><td>(a)</td><td>257</td><td>0.128</td><td>0.246</td><td>0.0013</td><td>0.206</td><td>0.356</td><td>0.0026</td><td>0.263</td><td>0.439</td><td>0.0073*</td></tr><tr><td>2 close</td><td>(a)</td><td>257</td><td>0.213</td><td>0.382</td><td>0.0027</td><td>0.371</td><td>0.606</td><td>0.0066*</td><td>0.485</td><td>0.734</td><td>0.0101*</td></tr><tr><td>1 close</td><td>(b) + CSF</td><td>61</td><td>0.299</td><td>0.896</td><td>0.0038</td><td>0.446</td><td>0.971</td><td>0.0066</td><td>0.597</td><td>1.02</td><td>0.0103</td></tr><tr><td>1 close</td><td>same w/o font.</td><td>61</td><td>0.333</td><td>0.869</td><td>0.0039</td><td>0.435</td><td>0.933</td><td>0.0068</td><td>0.581</td><td>0.989</td><td>0.0106</td></tr><tr><td>1 close</td><td>(b) + CSF</td><td>257</td><td>0.264</td><td>0.878</td><td>0.0032</td><td>0.367</td><td>0.938</td><td>0.0049</td><td>0.448</td><td>0.965</td><td>0.0068</td></tr><tr><td>1 close</td><td>same w/o font.</td><td>257</td><td>0.279</td><td>0.762</td><td>0.0035</td><td>0.346</td><td>0.878</td><td>0.0048</td><td>0.437</td><td>0.934</td><td>0.0068</td></tr><tr><td>2 close</td><td>(b) + CSF</td><td>61</td><td>0.272</td><td>0.904</td><td>0.0053</td><td>0.393</td><td>0.969</td><td>0.0073*</td><td>0.522</td><td>1.02</td><td>0.0171*</td></tr><tr><td>2 close</td><td>same w/o font.</td><td>61</td><td>0.34</td><td>0.905</td><td>0.0063*</td><td>0.389</td><td>0.948</td><td>0.0103*</td><td>0.499</td><td>0.99</td><td>0.0182*</td></tr><tr><td>2 close</td><td>(b) + CSF</td><td>257</td><td>0.218</td><td>0.796</td><td>0.004</td><td>0.338</td><td>0.916</td><td>0.0071</td><td>0.431</td><td>0.968</td><td>0.0122*</td></tr><tr><td>2 close</td><td>same w/o font.</td><td>257</td><td>0.31</td><td>0.811</td><td>0.0035*</td><td>0.327</td><td>0.897</td><td>0.0076</td><td>0.403</td><td>0.942</td><td>0.0115*</td></tr></tbody></table>

## 6. 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). It provides an interesting homogeneous or heterogeneous multi-layer 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 of 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 the processing of time varying records of EEG/MEG signal. FS3D is yet equipped 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].

## References

[1] M. ANTONAKAKIS, S. SCHRADER, U. AYDIN, A. KHAN, J. GROSS, M. ZERVAKIS, S. RAMPP, C.H. WOLTERS, *Inter-Subject Variability of Skull Conductivity and Thickness in Calibrated Realistic Head Models*, Neuroimage, **223** (2020), 1173-1175.

[2] B. ATFEH, L. BARATCHART, J. LEBLOND, J.R. PARTINGTON *Bounded extremal and Cauchy-Laplace problems on the sphere and shell*, J. Fourier Analysis and Applications, **16**(2) (2010), 177-203.

[3] H. AZIZOLLAHI, M. DARBAS, M. DIALLO, A. EL BADIA, S. LOHRENGEL, *EEG in neonates: forward modeling and sensitivity analysis with respect to variations of the conductivity*, Mathematical Biosciences and Engineering, **15**(4) (2018), 905-932.

[4] L. BARATCHART, A. BEN ABDA, F. BEN HASSEN, J. LEBLOND, *Recovery of pointwise sources or small inclusions in 2D domains and rational approximation*, Inverse Problems, **21** (2005), 51-74.

[5] L. BARATCHART, L. BOURGEOIS, J. LEBLOND, *Uniqueness results for inverse Robin problems with bounded coefficients*, J. Functional Analysis, **270**(1) (2016), 2508-2542.

[6] C.G. BÉNAR, R.N. GUNN, C. GROVA, B. CHAMPAGNE, J. GOTMAN, *Statistical Maps for EEG Dipolar Source Localization*, IEEE Trans. on Biomed. Eng., **52**(3) (2005).

[7] S.R. BENBADIS, S. BENICZKY, E. BERTRAM, S. MACIVER, S. L. MOSHÉ, *The role of EEG in patients with suspected epilepsy*, Epileptic Disord **22**(2) (2020), 143-155.

[8] L. BOURGEOIS, *A mixed formulation of quasi-reversibility to solve the Cauchy problem for Laplace's equation*, Inverse Problems, **21** (2005), 1087-1104.

[9] L. BOURGEOIS, L. CHESNEL, *On quasi-reversibility solutions to the Cauchy problem for the Laplace equation: regularity and error estimates*, Math. Mod. Num. Anal., **54** (2020), 493-529.

[10] L. BOURGEOIS, J. DARDÉ, *A duality-based method of quasi-reversibility to solve the Cauchy problem in the presence of noisy data*, Inverse Problems, **26**(9) (2010), 095016.

[11] L. BOURGEOIS, A. RECOQUILLAY, *A mixed formulation of the Tikhonov regularization and its application to inverse PDE problems*, Mathematical Modelling and Numerical Analysis, **52**(1) (2018), 123-145.

[12] M. CLERC, J. LEBLOND, J.-P. MARMORAT, T. PAPADOPOULO, *Source localization in EEG using rational approximation on plane sections*, Inverse Problems **28** (2012), 055018.

[13] M. CLERC, J. LEBLOND, J.-P. MARMORAT, C. PAPAGEORGAKIS, *Uniqueness result for an inverse conductivity recovery problem with application to EEG*, Rendiconti dell'Istituto di Matematica dell'Università di Trieste. An International Journal of Mathematics (2016), 48.
[14] F. COSTA, H. BATATIA, T. OBERLIN, C. D'GIANO, J.-Y. TOURNERET, *Bayesian EEG source localization using a structured sparsity prior*, Neuroimage, **144** (2017), 142-152.

[15] C. DAPOGNY, C. DOBRZYNSKI, P. FREY, *Three-dimensional adaptive domain remeshing, implicit domain meshing, and applications to free and moving boundary problems*, Journal of Computational Physics **262** (2014), 358-378.

[16] M. DARBAS, M. DIALLO, A. EL BADIA, S. LOHRENGEL, *An inverse dipole source problem in inhomogeneous media: application to the EEG source localization in neonates*, Journal of Inverse and Ill-Posed Problems, **27**(2) (2019), 255-282.

[17] M. DARBAS, S. LOHRENGEL, *Review of mathematical modeling of electroencephalography (EEG)*, Jahresbericht der Deutschen Mathematiker-Vereinigung, **121**(1) (2019), 3-39.

[18] J. DARDÉ, *Méthodes de quasi-réversibilité et de lignes de niveau appliquées aux problèmes inverses elliptiques*, Thèse de doctorat, Université Paris-Diderot, 2010.

[19] J. DARDÉ, *Iterated quasi-reversibility method applied to elliptic and parabolic data completion problems*, Inverse Problems and Imaging, **10**(2), (2016), 379-407.

[20] R. DAUTRAY, J.-L. LIONS, *Mathematical Analysis and Numerical Methods for Science and Technology*, Vol. 1, Springer, 2000.

[21] A. DELORME, S. MAKEIG, *EEGLAB: an open source toolbox for analysis of single-trial EEG dynamics including independent component analysis*, J. Neuroscience Methods, **134**(1), (2004).

[22] M. M. DIALLO, *Problème inverse de sources en Electro-Encéphalo-Graphie chez le nouveau-né*, Thèse de doctorat, Université de Picardie Jules Verne, France, 2017.

[23] A. EL BADIA, T. HA DUONG, *An inverse source problem in potential analysis*, Inverse Problems, **16** (2000), 651-663.

[24] I. GAUDET, A. HÜSSER, P. VANNASING, A. GALLAGHER, *Functional Brain Connectivity of Language Functions in Children Revealed by EEG and MEG: A Systematic Review*, Front. Hum. Neurosci **14**:62 (2020).

[25] A. GRAMFORT, M. KOWALSKI, M. HÄMÄLÄINEN, *Mixed-norm estimates for the M/EEG inverse problem using accelerated gradient methods*, Phys. Med. Biol. **57**(7) (2012), 1937-1961.

[26] A. GRAMFORT, M. LUESSI, E. LARSON, D. ENGEMANN, D. STROHMEIER, C. BRODBECK, L. PARKKONEN, M. HÄMÄLÄINEN, *MNE software for processing MEG and EEG data*, Neuroimage, **86**(10) (2013).

[27] A. GRAMFORT, T. PAPADOPOULO, E. OLIVI, M. CLERC, *Forward field computation with OpenMEEG*, Computational Intelligence and Neuroscience, 2011(1687-5265):923703.

[28] R. GRECH, T. CASSAR, J. MUSCAT, K.P. CAMILLERI, S.G.FABRI, M. ZERVAKIS, P. XANTHOPOULOS, V. SAKKALIS, B. VANRUMSTE, *Review on solving the inverse problem in EEG source analysis*, J. Neuro. Rehabil. **5**(25) (2008), 1-33.

[29] C. GROVA, J. DAUNIZEAU, J.M. LINA, C.G. BÉNAR, H. BENALI, J. GOTMAN, *Evaluation of EEG localization methods using realistic simulations of interictal spikes Neuroimage*, **29**(3) (2006), 734-753.

[30] J. HADAMARD, *Sur les probèmes aux dérivées partielles et leur signification physique*, Princeton University Bulletin, **13**(4), (1902), 49-52.

[31] M. HÄMÄLÄINEN, R. HARI, J. ILMONIEMI, J. KNUUTILA, O.V. LOUNASMAA, *Magnetoencephalography-theory, instrumentation, and applications to noninvasive studies of the working human brain*, Rev. Mod. Phys., **65** (1993), 413-497.

[32] F. HECHT, *New Development in FreeFem++*, Journal of Numerical Mathematics **20**(3-4) (2012), 251-265.

[33] R.J. ILMONIEMI, J. SARVAS, *Brain Signals: Physics and Mathematics of MEG and EEG*, Ch. 3, MIT Press (2019).

[34] J. D. JACKSON, *Classical Electrodynamics*, Third edition, John Wiley & Sons, 1998.

[35] M. V. KLIBANOV, F. SANTOSA, *A computational quasi-reversibility method for cauchy problems for*
*Laplace's equation*, SIAM J. Appl. Math., **51** (1991), 1653-1675.

[36] V.A. KOZLOV, V.G. MAZYA, A.V. FOMIN, *An iterative method for solving the Cauchy problem for elliptic equation*, Comput. Math. Phys., **31** (1991), 45-52.

[37] R. KRESS, *Linear Integral Equations*, Second edition, Applied Mathematical Sciences 82, Spinger-Verlag, 1999.

[38] J. KYBIC, M. CLERC, T. ABBOUD, O. FAUGERAS, R. KERIVEN, T. PAPADOPOULO, *A common formalism for the integral formulations of the forward EEG problem*, IEEE Trans. on Med. Imaging, **24** (2005), 12-28.

[39] R. LATTÈS, J.-L. LIONS, *Méthode de Quasi-réversibilité et Applications*, Dunod, Paris, 1967.

[40] J. LEBLOND, C. PADURET, S. RIGAT, M. ZGHAL, *Source localization in ellipsoids by best meromorphic approximation in planar sections*, Inverse Problems, **24**(3) (2008), 035017.

[41] J. LE ROUSSEAU, G. LEBEAU, *On Carleman estimates for elliptic and parabolic operators. Applications to unique continuation and control of parabolic equations*, ESAIM: Control, Optimisation and Calculus of Variations, **18**(3) (2012), 712-747.

[42] F.-H. LIN, J.W. BELLIVEAU, A.M. DALE, M. HÄMÄLÄINEN, *Distributed current estimates using cortical orientation constraints*, Human Brain Mapp., **27**(1) (2006).

[43] B. LÜTKENHÖNER, *Dipole source localization by means of maximum likelihood estimation. I. Theory and simulations*, Electroencephalography and Clinical Neurophysiology, **106**(4) (1998), 314-321.

[44] N. MÄKELÄ, M. STENROOS, J. SARVAS, R.J. ILMONIEMI, *Truncated rap-music (trap-music) for meg and EEG source localization*, Neuroimage, **167** (2018), 73-83.

[45] G. MALFILÂTRE, L. MONY, D. HASAERTS, P. VIGNOLO-DIARD, M.-D. LAMBLIN, E. BOUREL-PONCHEL, *Technical recommendations and interpretation guidelines for electroencephalography for premature and full-term newborns*, Neurophysiologie Clinique, **51**(1) (2021), 35-60.

[46] J.C. MOSHER, R.M. LEAHY, *Recursive MUSIC: a framework for EEG and MEG source localization*, IEEE Trans. IEEE Trans. on Biomed. Eng., **45**(11) (1998).

[47] J.C. MOSHER, R.M. LEAHY, P.S. LEWIS, *EEG and MEG: forward solutions for inverse methods*, IEEE Trans. on Biomed. Eng., **46** (1999), 245-259

[48] J.C. MOSHER, R.M. LEAHY, *Source localization using recursively applied and projected (RAP) MUSIC*, IEEE Trans. on Signal Proc., **47**(2) (1999), 332-340.

[49] J. C. DE MUNCK, M. J. PETERS, *A fast method to compute the potential in the multisphere model*, IEEE Trans. Biomed. Eng. **40** (1993), 1166-1174.

[50] M. NEMAIRE, *Inverse potential problems, with applications to quasi-static electromagnetics*, PhD Thesis, Université Côte d'Azur, France, 2023.

[51] D. NIKODIJEVIC, N. BANEVA-DOLNENEC, D. PETROVSKA-CVETKOVSKA, D. CAPAROSKA, *Refractory Epilepsy-MRI, EEG and CT scan, a Correlative Clinical Study*, Open Access Maced J Med Sci. **4**(1) (2016), 98-101.

[52] C. PAPAGEORGAKIS, *Patient specific conductivity models: characterization of the skull bones*, PhD Thesis, Université Côte d'Azur, France, 2017.

[53] J. PESCATORE, I. BLOCH, S. BAILLET, L. GARNERO, *F.E.M. tetrahedral mesh of head tissues from M.R.I. under geometrical and topological constraints for applications in E.E.G. and M.E.G.*, Neuroimage **13**(6) (2001), 218.

[54] N. ROCHE-LABARBE, A. AARABI, G. KONGOLO, C. GONDRY-JOUET, M. DÜMPELMANN, R. GREBE, F. WALLOIS, *High-resolution electroencephalography and source localization in neonates*, Hum. Brain Mapp., **29** (2008), 167-176.

[55] M. SALO, *Unique continuation for elliptic equations*, Notes, Fall, University of Jyväskylä (2014).

[56] M. SCHNEIDER, *A multistage process for computing virtual dipole sources of EEG discharges from surface information*, IEEE Trans. on Biomed. Eng., **19** (1972), 1-12.

[57] K. SEKIHARA, S. NAGARAJAN, D. POEPPEL, A. MARANTZ, Y. MIYASHITA, *Application of an MEG*
*Eigenspace Beamformer to reconstructing spatio-temporal activities of neural sources*, Human Brain Map **15** (2002) 199–215.

[58] S. SOMMARIVA, A. SORRENTINO, *Sequential Monte Carlo samplers for semi-linear inverse problems and application to magnetoencephalography*, Inverse Problems, **30** (2014).

[59] A. SORRENTINO, L. PARKKONEN, A. PASCARELLA, C. CAMPI, M. PIANA, *Dynamical MEG source modeling with multi-target Bayesian filtering*, Hum. Brain Mapp., **30** (2009), 1911–1921.

[60] J. VORWERK, *New Finite Element Methods to Solve the EEG/MEG Forward Problem*, Ph.D. thesis, University of Münster, Münster, Germany (2016).

[61] J. VORWERK, M. CLERC, M. BURGER, C. WOLTERS, *Comparison of boundary element and finite element approaches to the EEG forward problem*, Biomedizinische Technik. Biomedical engineering, **57** (2012).

[62] F. WALLOIS, L. ROUTIER, C. HEBERLÉ, M. MAHMOUDZADEH, E. BOUREL-PONCHEL, S. MOGHIMI, *Back to basics: the neuronal substrates and mechanisms that underlie the electroencephalogram in premature neonates*, Neurophysiol. Clin. **51**(1)(2021), 5–33.

[63] D. WIPF, S. NAGARAJAN, *A unified Bayesian framework for MEG/EEG source imaging*, Neuroimage, **44**(3) (2009), 947–966.

[64] C.H. WOLTERS, H. KÖSTLER, C. MÖLLER, J. HÄRDTLEIN, L. GRASEDYCK, W. HACKBUSCH, *Numerical mathematics of the subtraction approach for the modeling of a current dipole in EEG source reconstruction using finite element head models*, SIAM J. Sci. Comput., **30**(1) (2007), 24–45.