1. Introduction

Electroencephalography (EEG) is a non-invasive functional brain imaging technique. EEG measures the electrical activity of the brain recorded by electrodes on the scalp and, more precisely, the voltage potential fluctuations between different cortical regions on the scalp. The recorded electrical activity is the synchronous activity of a large number of neighboring well-oriented neurons in the cerebral cortex beneath the skull.

The important goal of functional brain imaging using EEG is to localize cerebral sources generating measured EEG signals. The measurements provide valuable information about the sources that are at the origin of physological and pathological activities of the brain. In particular, EEG is one of the main diagnostic tests in presurgical evaluation for refractory epilepsy. Mathematical models are required to relate neural sources to EEG measurements. The accuracy of the EEG source reconstruction relies heavily on the accuracy of the associated forward model [1] which consists in computing the potential on the scalp for a given electrical source located in the brain.

On the one hand, the spherical multi-layer head model has gathered much interest from the beginning of EEG source analysis since an asymptotic formula for the potential is available [32,42]. On the other hand, realistic head models obtained from segmentation of magnetic resonance imaging (MRI) are able to take into account the geometrical complexity of the different tissues and their specific electrical characteristics. Several source models have been developed as e.g. partial integration, the St. Venant model or the subtraction approach [12,40,5,31,2]. For the numerical resolution of the forward problem, both boundary elements and 3D finite elements are commonly used [21,20,26,27]. In the head model for adults, the effect of anisotropy as well as the uncertainty in the tissue conductivities have been investigated [39]. All the aforementioned results are concerned with head models for adults or elderly children. In this paper, we are interested in an EEG model for neonates. The study is motivated by clinical questions. Few studies of neonates [35,18] exist due to the difficulties in creating realistic head models in neonates. Two characteristics are inherent to the neonate. The first one is that neonates show higher skull conductivities than children or adults [33]. The second one is the presence of fontanels in the skull. Fontanels are in the process of ossification and possess different electrical properties in comparison to the skull (cf. Figure 1.1). Clinicians wonder if it would be realistic to use the same EEG forward modeling for adults, early children and neonates. The underlying question is thus to quantify the impact of fontanels and the effects of uncertainty in the skull and fontanel conductivity on the values of the electric potential on the scalp. This investigation could be used directly for the development of realistic source localization methods from neonatal EEG.

From this motivation, the present work aims at proposing a mathematical EEG model for neonates and its numerical validation for both the multi-layer spherical head model and a realistic neonate head model. It may be noticed that the inhomogeneity of the skull conductivity prevents the application of boundary element methods which are currently used in commercial software for EEG source localization.

In order to answer clinical questions, we compare the proposed EEG model (with fontanels) and the one for adults (without fontanels). In neurophysiology studies, comparisons are commonly done via the computation of two error functionals, respectively the RDM (Relative Difference Measure) and the MAG (MAGnification factor). In this paper, we introduce an additional analysis tool to investigate the sensitivity of the potential solution of the problem with respect to the variations in the skull and fontanel conductivities. From a mathematical point of view, sensitivity is the directional derivative of the solution with respect to conductivity. It allows to analyze small variations in the tissue conductivities which currently occur from one patient to the other. Sensitivity analysis with respect to the conductivity can be employed to gather information about optimal design of the electrodes for EEG head caps [4,38]. In the present paper, sensitivity analysis is at the service of studying the effect of uncertainty in the knowledge of skull and fontanel conductivity on EEG measurements. Indeed, these parameters are not well known in the special case of neonates, and values of different orders of magnitude are proposed in the literature. We define in a rigorous way the mathematical setting of the sensitivity equation. Furthermore, the numerical discussion is performed on a realistic head model [3] obtained recently by the Inserm U1105 at Amiens' hospital. This model takes into account the fontanels and is more precise than previous ones [18].

The paper is organized as follows. In Section 2, we derive the EEG forward model. In Section 3, we present the subtraction approach to deal with the singularity of the source term and prove an existence and uniqueness result of a weak solution. In Section 4, we present a sensitivity analysis of the potential with respect to the conductivity. Section 5 is devoted to the numerical part: discretization, convergence analysis and various simulations. We discuss the validation of the EEG model, the impact of fontanels and the sensitivity of the potential with respect to a perturbation of the fontanel and/or skull conductivity.

---

2. The EEG forward problem in neonates

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 [22,17]. In terms of the electric field $\mathbf{E}$ and the magnetic field $\mathbf{H}$, this yields

$\nabla \cdot (\varepsilon\mathbf{E}) = \rho, \label{maxwell1}$

(2.1a)

$\text{curl}\mathbf{E} = 0,$

(2.1b)

$\text{curl}\mathbf{H} = \mathbf{J},$

(2.1c)

$\nabla \cdot (\mu\mathbf{H}) = 0.$

(2.1d)

Here, $\rho$ is the charge density, $\varepsilon$ and $\mu$ are, respectively, the electric permittivity and magnetic permeability, and $\mathbf{J}$ is the electric current density. Following Ohm's law, the current density splits into two terms,

$\mathbf{J} = \sigma\mathbf{E} + \mathbf{J}^s,$

(2.2)

where $\mathbf{J}^s$ is the density of the impressed neural currents and $\sigma$ denotes the conductivity distribution in the human head. It follows from (2.1b), that the electric field derives from a scalar potential $u$, i.e.

$\mathbf{E} = \nabla u.$

(2.3)

Now, consider a bounded regular domain $\Omega\subset\mathbb{R}^3$ with boundary $\partial\Omega$. Taking the divergence of (2.1c) together with Ohm's law (2.2) yields the following transmission equation for the electric potential $u$ in $\Omega$

$\label{e:elliptic} -\nabla \cdot (\sigma\nabla u) = \nabla \cdot {\boldsymbol{J}}^s.$

(2.4)

The source model of neural activity can be described by a sum of $M$ electric dipoles located in the brain (e.g. [22,37]). Each dipole is characterized by its position $S_m\in \Omega_1$ and its moment $\boldsymbol{q}_m$ which is a vector of $\mathbb{R}^3$. The points $(S_m)_m$ are assumed to be mutually distinct and the vectors $\bf{q_m}$ are non vanishing, i.e.

$S_m \neq S_p\, \, \, \, \forall m\neq p \, \, \, \, \mbox{and}\, \, \, \, {\bf{q_m} \neq 0} \, \, \, \, \forall m, p \in \{1, ..., M\}.$

(2.5)

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

${{\boldsymbol{J}}^s} = \sum\limits_{m = 1}^M {{{\boldsymbol{q}}_m}} {\delta _{{S_m}}}$

where $\delta_{S_m}$ denotes the delta distribution at $S_m$. The right hand side of (2.4) is then given by

$F \stackrel{\rm def}{ = } \nabla \cdot \boldsymbol{J}^s = \sum\limits_{m = 1}^M \boldsymbol{q}_m\cdot\nabla \delta_{S_m}.$

(2.6)

In the typical multi-layer head model, we distinguish three to five layers for the brain (containing gray and white matters, cerebrospinal fluid (CSF)), skull, and scalp. Therefore, consider a partition of $\Omega$ into $L$ open subdomains $(\Omega_i)_{i = 1, \ldots, L}$ (see Figure 2.1) such that

$\overline{\Omega} = \bigcup\limits_{i = 1}^L \overline{\Omega}_i\ \mbox{and} \\ \Omega_i\cap \Omega_j = \emptyset\ \forall i\neq j.$

We assume the subdomains to be nested which means that

$\overline{\Omega_i}\cap\overline{\Omega_j} = \emptyset\ \forall j\neq i-1, i, i+1.$

Then, for $i = 1, \ldots, L-1$, we denote by $\Gamma_i$ the interface between $\Omega_i$ and $\Omega_{i+1}$, and assume that $(\Gamma_i)_i$ are closed regular surfaces. Let $\boldsymbol{n}_i$ be the unit normal vector on $\Gamma_i$ oriented towards the exterior of $\Omega_i$. Notice that in a more general setting, we could allow interfaces $\Gamma_{ij}$ between two arbitrary subdomains $\Omega_i$ and $\Omega_j$ provided the interfaces are not intersecting. We finally denote by ${\Gamma _\infty } = \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 2.1). Now, let $\sigma_i \stackrel{\rm def}{ = } \sigma_{|\Omega_i}$ $(i = 1, \ldots, L)$ denote the conductivity of the subdomain $\Omega_i$. In the head model of an adult, the different layers are assumed to be homogeneous and isotropic, and therefore each $\sigma_i$ is a positive constant. For the EEG neonates' model, we propose to model the presence of fontanels in the skull via a variable conductivity. Their impact will thus be measured by the sensitivity of the model with respect to the conductivity. We will consider in the sequel conductivities $\sigma_i$ that are functions of the position, i.e. $\sigma_i = \sigma_i(\boldsymbol{x})$ in $\Omega_i$.

By considering that no electric current can flow out of the skull, the electric potential $u$ is then solution of the following boundary problem with homogeneous Neumann condition

$\label{e:bvp} \left\{ \begin{array}{rcll} -\nabla \cdot (\sigma\nabla u) = F,&\mbox{in}\ \Omega, \\ \sigma \partial_{\boldsymbol{n}} u = 0,&\mbox{on}\ {\Gamma _\infty }, \end{array} \right.$

(2.7)

where the source term $F$ is given by (2.6). Since $F$ vanishes identically in a neighborhood of the interfaces $\Gamma_i$, $u$ satisfies the transmission conditions

$[u]_{|\Gamma_i} = 0 \ \,\,\,\mbox{on}\,\,\, \Gamma_i\ (i = 1, \ldots, L-1), \label{e:transmission1}$

(2.8a)

$[\sigma\partial_{\boldsymbol{n}} u]_{|\Gamma_i} = 0\ \,\,\,\mbox{on}\,\,\, \Gamma_i\ (i = 1, \ldots, L-1).$

(2.8b)

Here and below, if $f_i = f_{|\Omega_{i}}, i = 1, \ldots, L$, $[f]_{|\Gamma_i} = f_{i|\Gamma_i} - f_{i+1|\Gamma_i}$ denotes the jump of the quantity $f$ across the interface $\Gamma_i$. For given sources $(S_m, \boldsymbol{q}_m)_m$ and a known distribution of the conductivity $\sigma$, problem (2.7) is the forward EEG problem.

---

3. Existence and uniqueness result

Problem (2.7) enters within the framework of standard elliptic problems, and its analysis is essentially covered by the classical results in variational theory. However, a careful formulation of these results (existence, uniqueness and regularity) in the context of a piecewise regular conductivity is of great importance. Indeed, the interest of studying the direct problem is twofold: on the one hand, these results are necessary to the sensitivity analysis with respect to the conductivity (Section 4), on the other, they allow the validation of the implementation of the method (Section 5.2).

Notice that a direct variational formulation of (2.7) in $H^1(\Omega)$ is not possible since the source term $F$ given by (2.6) belongs to $H^s(\mathbb{R}^3)$ for $s<-5/2$. In order to overcome this problem, we adapt an idea that has been presented in [12] in the case of piecewise constant conductivities. This method has also been introduced in [16,40,2] as the subtraction approach in the case of a single source and can be roughly described as follows: let $\widetilde u = - G\star F$ where $G$ is the fundamental solution of the Laplacian. Then, the function $w = u-\widetilde u$ satisfies $-\sigma_i\Delta w = 0$ in the sub-domains $\Omega_i$ and non-homogeneous transmission conditions with regular right hand sides on the interfaces $\Gamma_i$. Thus, $w$ is solution of a variational problem. Here, we aim to address the question of non-homogeneous conductivities. In addition to classical regularity assumptions on the functions $\sigma_i$, we will state below necessary and sufficient conditions on the behavior of these functions near the source points that allow to adapt the subtraction approach to this new configuration. We deal also with the general case of multiple dipolar sources.

---

3.1. Lifting of the singularity.

The principle of the subtraction method is to decompose the potential $u$, solution to (2.7), into a potential $\widetilde u$ which contains the singularity and a regular lifting $w$

$\label{changev} u = \widetilde u + w, \ \mbox{with}\ \widetilde u = \sum\limits_{m = 1}^M \widetilde u_m.$

(3.1)

In order to define the singular potentials $\widetilde u_m$, we assume in the sequel of the paper that the conductivity is constant in a neighborhood of each source. More precisely we fix a family of open balls $(\mathcal{V}_m)_{m = 1, \ldots, M}$ such that $\mathcal{V}_m\subset\subset \Omega_1$, $S_m\in\mathcal{V}_m$ and

$\sigma_{1|\mathcal{V}_m} \equiv c_m \in \mathbb{R}, \,\,\,\mbox{for any}\,\,\, m\in\{1, \ldots, M\}.$

(3.2)

Furthermore, since the locations $S_k$ are mutually distinct (2.5), one can assume that the balls $\mathcal{V}_m$ are non intersecting, $\mathcal{V}_m\cap \mathcal{V}_{p} = \emptyset$ if $m\neq p$.

The singular potential $\widetilde u_m$ is then solution of the following Poisson equation

$\label{poisson} -c_m\Delta\widetilde u_m = \boldsymbol{q}_m\cdot\nabla \delta_{S_m}\ {\rm in}\ \mathbb{R}^3,$

(3.3)

where $c_m = \sigma_1(S_m)$. It can be obtained by convolution of the fundamental solution of the Laplace equation with the right hand side $\dfrac{1}{c_m}\boldsymbol{q}_m\cdot\nabla\delta_{S_m}$ as follows

$\label{potsing} \widetilde u_m (\boldsymbol{x}) = \frac{1}{4\pi c_m}\boldsymbol{q}_m\cdot\frac{(\boldsymbol{x}-S_m)}{|\boldsymbol{x}-S_m|^3}, \ \forall\boldsymbol{x}\in\mathbb{R}^3\setminus \{S_m\}.$

(3.4)

We see that the potential $\widetilde u$ has a singularity at each source point $S_m$, but is smooth everywhere else.

In order to identify the problem satisfied by $w$, notice that for any $m\in\{1, \ldots, M\}$, the quantity $\nabla \cdot (\sigma\nabla \widetilde u_m)$ is well defined on $\Omega$ by

$\nabla \cdot (\sigma\nabla \widetilde u_m) = \nabla \cdot \left((\sigma-c_m)\nabla \widetilde u_m \right) + c_m\Delta \widetilde u_m.$

Indeed, both terms on the right hand side of the above identity are well defined as distributions on $\Omega$ since $\sigma - c_m$ vanishes identically on $\mathcal{V}_m$ and $\widetilde u_m$ is regular outside $\displaystyle \cup_m \mathcal{V}_m$. Therefore, we get

$-\nabla \cdot (\sigma\nabla w ) = -\nabla \cdot \left(\sigma\nabla(u-\widetilde u) \right) = F + \sum\limits_{m = 1}^M \nabla \cdot \left((\sigma-c_m)\nabla \widetilde u_m \right) + c_m\Delta\widetilde u_m\ \mbox{on}\ \Omega.$

It follows from the definition of $\widetilde u_m$ that $F = -\displaystyle\sum_{m = 1}^M c_m\Delta\widetilde u_m$, and $w$ is thus formally solution of the following boundary value problem

$\label{e:edp_Vsig} \left\{ \begin{array}{rcll} -\nabla \cdot (\sigma\nabla w) = \displaystyle\sum\limits_{m = 1}^M \nabla \cdot ((\sigma-c_m)\nabla \widetilde u_m)&\mbox{in}\ \Omega, \\ \sigma \partial_{\boldsymbol{n}} w = -\sigma\partial_{\boldsymbol{n}}\widetilde u,&\mbox{on}\ {\Gamma _\infty }, \end{array} \right.$

(3.5)

Problem (3.5) can be reformulated on each subdomain $\Omega_i$ provided the conductivity $\sigma$ is piecewise regular. To this end, assume that $\sigma_i\in W^{1, \infty}(\Omega_i)$. We may notice that the first term on the right hand side of the following identity vanishes, since $\Delta\widetilde u_m$ is a distribution with a single point support according to (3.3) and vanishes outside the neighborhood $\mathcal{V}_m$,

$\nabla \cdot \left( (\sigma_1-c_m)\nabla \widetilde u_m \right) = (\sigma_1-c_m)\Delta \widetilde u_m + \nabla(\sigma_1 -c_m)\cdot\nabla \widetilde u_m.$

(3.6)

We thus get

$-\nabla \cdot (\sigma_1\nabla w ) = \sum\limits_{m = 1}^M \nabla(\sigma_1 -c_m)\cdot\nabla \widetilde u_m = \nabla\sigma_1\cdot\nabla\widetilde u.$

On $\Omega_i$, $i \neq 1$, the potentials $\widetilde u_m$ are regular and we see from a direct computation that $w$ is solution of the following problem

$\label{e:bvpreform} \left\{ \begin{array}{rcll} -\nabla \cdot (\sigma_i\nabla w) = \nabla \sigma_i\cdot\nabla \widetilde u,&\mbox{in}\ \Omega_i (i = 1, \ldots, L), \\ \sigma \partial_{\boldsymbol{n}} w = -\sigma\partial_{\boldsymbol{n}}\widetilde u,&\mbox{on}\ {\Gamma _\infty }, \end{array} \right.$

(3.7)

with transmission conditions

$[w]_{|\Gamma_i} = 0 \ \mbox{on}\,\,\, \Gamma_i, \label{e:transmission1reform},$

(3.8a)

$[\sigma\partial_{\boldsymbol{n}} w]_{|\Gamma_i} = (\sigma_{i+1}-\sigma_i)\partial_{\boldsymbol{n}}\widetilde u\ \mbox{on}\,\,\, \Gamma_i,$

(3.8b)

for $i = 1, \ldots, L-1$. Recall that the right hand side of (3.7) is well defined and vanishes in any neighborhood $\mathcal{V}_m$ of the sources, since $\sigma_1$ is constant on $\mathcal{V}_m$.

---

3.2. Variational formulation

In this section, we give a variational formulation of problem (3.7)-(3.8) for the auxiliary function $w$. Let $v\in H^1(\Omega)$. Multiplying (3.7) by $v_{|\Omega_i}$, integrating over $\Omega_i$ and summing up, we obtain the following formulation from Green's formula and the boundary and transmission conditions,

$\begin{align} &\sum\limits_{i = 1}^L\int_{\Omega_i}\sigma_i\nabla w\cdot\nabla v\, d\boldsymbol{x}\nonumber\\ = &\sum\limits_{i = 1}^{L-1}\int_{\Gamma_i}(\sigma_{i+1}-\sigma_i)\partial_{\boldsymbol{n}}\widetilde u v\, \, \, ds-\int_{{\Gamma _\infty }}\sigma_L\partial_{\boldsymbol{n}}\widetilde u v \, ds + \sum\limits_{i = 1}^L\int_{\Omega_i}(\nabla \sigma_i\cdot\nabla \widetilde u)v\, d\boldsymbol{x}.\label{Psig} \end{align}$

(3.9)

Since $\nabla (\sigma_i-c_m) = \nabla \sigma_i$, we can show with the help of (3.6) that (3.9) is equivalent to

$\int_{\Omega}\sigma\nabla w\cdot\nabla v\, d\boldsymbol{x} = \sum\limits_{m = 1}^M\left(\int_{\Omega}(c_m-\sigma)\nabla \widetilde u_m\cdot\nabla v\, d\boldsymbol{x}-\int_{{\Gamma _\infty }}c_m\partial_{\boldsymbol{n}}\widetilde u_m v \, ds\right)$

(3.10)

which is the variational formulation of the boundary value problem (3.5). In the following, we focus on formulation (3.10). We introduce the bilinear form $a(\cdot, \cdot)$

$\label{bilinearform} a(w, v) = \int_{\Omega}\sigma\nabla w\cdot\nabla v\, d\boldsymbol{x},$

(3.11)

as well as the linear form

$\label{linearform} l(v) = \sum\limits_{m = 1}^M\left(\int_{\Omega}(c_m-\sigma)\nabla \widetilde u_m\cdot\nabla v\, d\boldsymbol{x}-\int_{{\Gamma _\infty }}c_m\partial_{\boldsymbol{n}}\widetilde u_m v \, ds\right).$

(3.12)

Notice that $a(\cdot, \cdot)$ and $l(\cdot)$ are well defined on $H^1(\Omega)$ whenever the conductivity $\sigma$ belongs to $L^{\infty}(\Omega)$ and satisfies the assumption (3.2). Since (3.5) is a problem with Neumann boundary condition, the variational formulation (3.10) allows a solution only under the compatibility condition

$\label{compatcondition} l(1) = \sum\limits_{m = 1}^M \int_{{\Gamma _\infty }}c_m\partial_{\boldsymbol{n}}\widetilde u_m (x) \, ds = 0.$

(3.13)

Condition (3.13) follows from the following classical lemma given in [16] that we recall for the convenience of the reader.

Lemma 1.   Let $\widetilde u_m$ be the solution of equation (3.3) given by (3.4). Then

$\int_{{\Gamma _\infty }}c_m\partial_{\boldsymbol{n}}\widetilde u_m \, ds = 0\ \forall m = 1, \ldots, M.$

(3.14)

Note that a solution to (3.10) is unique only up to an additive constant. To this end, we introduce the following subspace of $H^1(\Omega)$ which does not contain any constant other than zero,

$V = \left\{ {{ v\in H^1(\Omega)}\left| {{\int_{\Omega} v\, d\boldsymbol{x} = 0 }} \right., } \right\}$

(3.15)

on which the Poincaré-Wirtinger inequality holds true,

$\| v \|_{0, \Omega} \le C_P \|\nabla u \|_{L^2(\Omega)}\ \forall v\in V.$

(3.16)

In the sequel, we write $a\lesssim b$ if there is a constant $C>0$ independent from the quantities $a$ and $b$ such that $a\le Cb$.

Theorem 1. Let $\sigma\in L^{\infty}(\Omega)$ be such that $0<\sigma_{\min}\leq\sigma(\boldsymbol{x})\leq\sigma_{\max}$ for almost any $\boldsymbol{x}\in\Omega$, where $\sigma_{\min}$ and $\sigma_{\max}$ are two given positive constants. Assume further that $\sigma_{|\mathcal{V}_m}$ is constant for any $m = 1, \ldots, M$ and denote by $c_m$ the value of $\sigma$ on $\mathcal{V}_m\subset\Omega_1$. Let the bilinear form $a(\cdot, \cdot)$ and the linear form $l(\cdot)$ be given by (3.11) and (3.12), respectively. Then the variational problem

$Find \,\,\,\,\,\,w\in V \,\,\,\,\,\,such\,\,\,\,that\,\,\,\,\,\,a(w, v) = l(v), \,\,\,\forall v\in H^1(\Omega)$

(3.17)

has exactly one solution $w\in V$. Moreover, the following estimate holds true,

$\| w \|_{H^1(\Omega)} \lesssim \sum\limits_{m = 1}^M \left( \| \nabla \widetilde u_m \|_{L^2(\Omega\setminus\mathcal{V}_m)} + \| \partial_\boldsymbol{n}\widetilde u_m \|_{L^2({\Gamma _\infty })} \right).$

(3.18)

Proof. It follows from standard arguments in variational theory that the bilinear form $a(\cdot, \cdot)$ is continous on $H^1(\Omega)\times H^1(\Omega)$ and $V$-elliptic. Since $c_m-\sigma$ vanishes in $\mathcal{V}_m$, we have

$|l(v)|\lesssim \sum\limits_{m = 1}^M\left(\|\nabla \widetilde u_m\|_{L^2(\Omega\backslash\mathcal{V}_m)}+C_T\|\partial_\boldsymbol{n}\widetilde u_m\|_{L^2({\Gamma _\infty })}\right)\|v\|_{H^1(\Omega)},$

(3.19)

where $C_T$ is the continuity constant of the trace operator from $H^1(\Omega)$ to $H^{1/2}({\Gamma _\infty })$. This proves that the linear form is continuous on $H^1(\Omega)$ and thus on $V$. The Lax-Milgram theorem then yields existence and uniqueness of a function $w\in V$ such that

$a(w, v) = l(v)\ \forall v\in V.$

Next, let $v$ belong to $H^1(\Omega)$. We have $v-v_\Omega\in V$ where $v_\Omega = \frac{1}{|\Omega|}\int_\Omega v\, d\boldsymbol{x}$ is the mean value of $v$. Since $l(1) = 0$ according to the compatibility condition (3.13) and $\nabla v = \nabla (v-v_\Omega)$, we get

$l(v) = l(v-v_\Omega) = a(w, v-v_\Omega) = a(w, v)$

which proves that $w$ is the unique solution of problem (3.17). Finally, estimate (3.18) follows from the coercivity of the bilinear form together with estimate (3.19) for the linear form $l$.

Remark 2.  Notice that the linear form $l(\cdot)$ in the variational formulation (3.17) is well defined if and only if $\sigma_1$ is constant in a neighborhood of the source positions $S_m$. Indeed, if this would not be the case, the term $(c_m-\sigma_1)\nabla\widetilde u_m$ does no longer belong to $(L^2(\Omega_1))^3$ due to the singularity of $\nabla\widetilde u_m$ that behaves as $\displaystyle\frac{1}{|\boldsymbol{x}-S_m|^3}$, even if $\sigma_1$ is regular enough such that the value $c_m = \sigma_1(S_m)$ makes sense.

The following theorem states the global $H^2$-regularity of the variational solution of (3.17) in the subdomains $\Omega_i$.

Theorem 3.   In addition to the assumptions of Theorem 1, assume that $\sigma_i\in W^{1, \infty}(\Omega_i)$ for any $i = \{1, \ldots, L\}$ and that $\Gamma_i$ is of class $\mathcal{C}^2$ for $i\in\{ 1, \ldots, L-1\} \cup \{\infty \}$. Let $w\in H^1(\Omega)$ be the solution of the variational problem (3.17). Then we have $w_{|\Omega_i}\in H^2(\Omega_i)$ for any $i = \{1, \ldots, L\}$ and

$\| w \|_{H^2(\Omega_i)} \lesssim \sum\limits_{m = 1}^M \left( \| \nabla\widetilde u_m \|_{H^1(\Omega_i\setminus\mathcal{V}_m)} + \|\partial_\boldsymbol{n}\widetilde u_m\|_{H^1({\Gamma _\infty })}\right).$

(3.20)

The proof of Theorem 3 relies on standard techniques for elliptic partial differential equations. Indeed, we may notice that on each $\Omega_i$, the variational solution $w$ satisfies the partial differential equation

$-\nabla\cdot( \sigma_i\nabla w) = f_i$

with

$f_i \stackrel{\rm def}{ = } \sum\limits_{m = 1}^M \nabla\cdot( (\sigma_i - c_m)\nabla\widetilde u_m) = \sum\limits_{m = 1}^M \nabla(\sigma_i-c_m)\cdot\nabla\widetilde u_m.$

According to the assumptions on $\sigma_i$, the function $f_i$ belongs to $L^2(\Omega_i)$. Hence, classical arguments for partial differential equations with variable coefficients apply and yield interior regularity in each $\Omega_i$, se for example [7,19]. $H^2$-regularity up to the boundary of the subdomains $\Omega_i$ follows since the boundary ${\Gamma _\infty }$ and the interfaces $\Gamma_i$ as well as the Neumann data $\sigma\partial_n\widetilde u$ are regular.

---

4. Sensitivity analysis with respect to a perturbation of the conductivity

Sensitivity indicates the behavior of the potential when there is a slight variation of physical parameters. Here, we are interested in the sensitivity with respect to conductivity. This permits to measure the effects of uncertainty in the skull and fontanel conductivity on the model. Mathematically, a rigorous way to describe sensitivity is given by Gâteaux differentiability which expresses a weak concept of derivative.

Definition 4.  Let $F:X\to Y$ be an application between two Banach spaces $X$ and $Y$. Let $U\subset X$ be an open set. The directional derivative $D_\mu F(\sigma)$ of $F$ at $\sigma\in U$ in the direction $\mu\in X$ is defined as

$D_\mu F(\sigma) = \lim\limits_{h\to 0} \frac{F(\sigma+h\mu) -F(\sigma)}{h}$

if the limit exists. If $D_\mu F(\sigma)$ exists for any direction $\mu\in X$ and if the application $\mu \mapsto D_{\mu} F(\sigma)$ is linear continuous from $X$ to $Y$, $F$ is called Gâteaux differentiable at $\sigma$.

Now, let $(\mathcal{V}_m)_m$ be a fixed family of neighborhoods of the sources such that for any $m$, $S_m\in\mathcal{V}_m\subset\Omega_1$ and $\mathcal{V}_m\cap\mathcal{V}_{m^\prime} = \emptyset$ if $m\neq m^\prime$. We introduce the parameter space

$\mathcal{P} = \left\{ {{\sigma\in L^\infty(\Omega)}|{\sigma_{|\mathcal{V}_m}\equiv {\rm const.}\, \forall m = 1, \ldots, M}} \right\}$

as well as the (open) subset

$\mathcal{P}_\text{adm} = \left\{ {{\sigma\in \mathcal{P}}|{ \sigma_{\min} < \sigma < \sigma_{\max} }} \right\}$

of admissible conductivities. Here, $\sigma_{\min}$ and $\sigma_{\max}$ are two fixed positive constants. According to Theorem 1, problem (3.10) with conductivity $\sigma\in\mathcal{P}_\text{adm}$ admits a unique solution $w(\cdot, \sigma)\in V$. The aim of this section is to prove differentiability of $w$ with respect to $\sigma$ and to identify its (Gâteaux) derivative in a given direction $\mu$. Since $w$ depends on the singular potential $\widetilde u$, we first analyze the derivative of $\widetilde u$ with respect to $\sigma$. To this end, recall that

$\widetilde u_m(\boldsymbol{x}, \sigma) = \frac{1}{4\pi c_m}\boldsymbol{q}_m\cdot\frac{(\boldsymbol{x}-S_m)}{|\boldsymbol{x}-S_m|^3}$

(4.1)

is the solution of the Poisson equation

$-c_m\Delta\widetilde u_m(\cdot, \sigma) = \boldsymbol{q}_m\cdot\nabla \delta_{S_m} \ {\rm in}\ \mathbb{R}^3$

where $c_m = \sigma(S_m)$. Now, consider an arbitrary direction $\mu\in\mathcal{P}$ and assume that $\sigma + h\mu$ belongs to $\mathcal{P}_\text{adm}$ for small values of $h$. By definition, we have

$\widetilde u_m(\boldsymbol{x}, \sigma+ h\mu) = \frac{1}{4\pi(c_m+h p_m)} \boldsymbol{q}_m\cdot\frac{(\boldsymbol{x}-S_m)}{|\boldsymbol{x}-S_m|^3},$

(4.2)

where $p_m = \mu(S_m)$, and $\widetilde u_m(\cdot, \sigma+ h\mu)$ is solution of the following perturbed Poisson equation,

$-(c_m+h p_m)\Delta\widetilde u_m(\cdot, \sigma+h \mu) = \boldsymbol{q}_m\cdot\nabla \delta_{S_m}\ {\rm in}\ \mathbb{R}^3.$

The following proposition states that $\widetilde u_m$ is Gâteaux differentiable at an interior point $\sigma\in\mathcal{P}$ and identifies its Gâteaux derivative:

Proposition 1.   Let $\sigma\in\mathcal{P}_\text{adm}$ and $h_0>0$ such that $\sigma+h\mu\in\mathcal{P}_\text{adm}$ for any $h\in [-h_0, h_0]$ and any $\mu\in\mathcal{P}$ with $\|\mu\|_{\infty} = 1$. Then, $\widetilde u_m(\cdot, \sigma)$ is Gâteaux differentiable at $\sigma$ and the Gâteaux derivative $D_\mu \widetilde u_m(\cdot, \sigma)$ of $\widetilde u_m$ at $\sigma$ in the direction $\mu$ reads

$\label{Dumtilde} D_\mu \widetilde u_m(\cdot, \sigma) = -\frac{p_m}{c_m} \widetilde u_m(\cdot, \sigma)$

(4.3)

with $p_m = \mu(S_m)$ and $c_m = \sigma(S_m)$.

Proof. A straightforward computation of the differential quotient yields

$\label{um1tilde} \frac{\widetilde u_m(\boldsymbol{x}, \sigma+h\mu) - \widetilde u_m(\boldsymbol{x}, \sigma)}{h} = -\frac{p_m}{c_m+hp_m} \widetilde u_m(\boldsymbol{x}, \sigma), \ \forall \boldsymbol{x} \neq S_m$

(4.4)

and (4.3) follows. The right-hand side of (4.3) is obviously linear and continuous in $\mu$ since $p_m = \mu(S_m)$.

Theorem 5.   Let $\sigma\in\mathcal{P}_\text{adm}$ and $h_0>0$ such that $\sigma+h\mu\in\mathcal{P}_\text{adm}$ for any $h\in [-h_0, h_0]$ and any $\mu\in\mathcal{P}$ with $\|\mu\|_{\infty} = 1$. Then the solution $w(\cdot, \sigma)$ of (3.10) is Gâteaux differentiable with respect to $\sigma$ and the Gâteaux derivative of $w$ at $\sigma$ in the direction $\mu\in\mathcal{P}$ is the unique solution of the following variational problem : find $w^1\in V$ such that

$\int_{\Omega}\sigma\nabla w^1\cdot\nabla v\, d\boldsymbol{x} = -\int_{\Omega}\mu\nabla w\cdot\nabla v\, d\boldsymbol{x} \label{Vsens}\\ \,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,+\sum\limits_{m = 1}^M\left(\int_{\Omega}(c_m-\sigma)\nabla \widetilde u_m^1\cdot\nabla v\, d\boldsymbol{x} -\int_{{\Gamma _\infty }}c_m\partial_{\boldsymbol{n}}\widetilde u_m^1v\, d\boldsymbol{x}\right)\nonumber \\\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\, +\sum\limits_{m = 1}^M\left(\int_{\Omega}(p_m-\mu)\nabla \widetilde u_m \cdot\nabla v\, d\boldsymbol{x} -\int_{{\Gamma _\infty }}p_m\partial_{\boldsymbol{n}}\widetilde u_m v\, d\boldsymbol{x}\right), \nonumber$

(4.5)

for all $v\in V$. Here, we note $\widetilde u_m^1 = D_\mu\widetilde u_m(\cdot, \sigma)$ and $\widetilde u_m = \widetilde u_m(\cdot, \sigma)$.

The sensitivity equation (4.5) is obtained by varying the conductivity parameter $\sigma$. Formally, the approach can be explained as follows (assume for sake of simplicity that there is a single source point $S_m$): let $w$ be the solution of the unperturbed equation

$-\nabla \cdot(\sigma\nabla w) = \nabla \cdot((\sigma-c_m)\nabla \widetilde u_m)$

and consider a perturbation of $\sigma$ of the form $\sigma + h\mu$ for fixed $\mu$. Denote by $w_h$ the solution of the perturbed equation, i.e.

$-\nabla \cdot((\sigma+h\mu)\nabla w_h) = \nabla \cdot((\sigma+h\mu-c_m-hp_m)\nabla \widetilde u_{m, h})$

where $p_m = \mu(S_m)$ and $\widetilde u_{m, h}$ is the solution of the perturbed Poisson equation. Subtracting the above equations and dividing by the scalar $h$ yields

$-\nabla \cdot \left(\sigma\nabla \left( \frac{w_h-w}{h} \right)\right) = \nabla \cdot(\mu\nabla w_h) + \nabla \cdot\left((\sigma-c_m)\nabla \left( \frac{\widetilde u_{m, h}-\widetilde u_m}{h} \right)\right)$

$+ \nabla \cdot((\mu-p_m)\nabla \widetilde u_{m, h}).$

At the limit $h\to 0$, we see that $w^1 = \lim_{h\to 0}\frac{w_h-w}{h}$ satisfies the following equation

$-\nabla \cdot(\sigma\nabla w^1) = \nabla \cdot(\mu\nabla w) + \nabla \cdot((\sigma-c_m)\nabla \widetilde u_m^1) + \nabla \cdot((\mu-p_m)\widetilde u_m).$

The boundary equation can be obtained in a similar way which yields finally the variational formulation (4.5). The full proof of Theorem 5 with a mathematical justification of the sense in which the limit has to be taken is given in Appendix A.

Remark 6.   Assume that the variation of the conductivity in a given direction $\mu$ occurs only in the skull $\Omega_2$. In this particular case, the singular potential $\tilde{u}$ defined by (3.4) is independent from $\mu$ and we have

$w(\cdot, \sigma+h\mu) - w(\cdot, \sigma) = u(\cdot, \sigma+h\mu) - \tilde{u}(\cdot, \sigma_1) - (u(\cdot, \sigma) + \tilde{u}(\cdot, \sigma_1)) = u(\cdot, \sigma+h\mu) -u(\cdot, \sigma).$

Thus, the Gâteaux derivative $u^1\stackrel{\rm def}{ = }u^\prime(\sigma;\mu)$ of the electric potential $u$ at $\sigma$ in the direction $\mu$ coincides with $w^1$, the solution to (4.5).

---

5. Finite Element formulation of the EEG forward problem

In this section, we address the discretization of problem (3.17). Boundary Element Methods (BEMs) have been very popular for a long time since they allow to reduce the three-dimensional problem to only boundary integral equations on the interfaces between the different tissues, see, for example [26,27,28] and references therein. However, BEMs are restricted to a piecewise constant conductivity representing in EEG the standard model of homogeneous tissues in an adult's head. Here, because of fontanels, BEMs are thus prohibited, and we propose a discretization by means of three-dimensional Lagrange Finite Elements (FE) which are able to handle complex geometries.

Problem (3.5), even if it seems to be rather standard, does not fit rigorously the standard assumptions in finite element analysis. For the convenience of the reader who could be less familiar with advanced finite element techniques, we address in this section the following points.

The first one is concerned with the discretization of the computational domain. With regard to the application in EEG and conforming to the assumptions in Section 2, the different subdomains $\Omega_i$ are assumed to be regular. The discrete problem will thus be set on discretized domains $\Omega_{i, h}$, and this approximation of the geometry has to be taken into account in the error analysis. This item has been addressed frequently in the finite element community and isoparametric finite elements have been shown to keep the error estimates optimal in the case of a single curved domain [8].

Further, we have to deal with the singularity of the source term. In [40], error estimates have been stated for the subtraction approach in the case of a single subdomain and constant conductivity.

The aim of this section is to give a precise description of the discrete problem as well as a rigorous error analysis for the complex setting of the multi-layer model with variable and piecewise regular conductivity. From a practical point of view, this is an important issue for the validation of the numerical implementation.

---

5.1. Discretization and convergence analysis

Throughout this subsection, we assume that the regularity assumptions of Theorem 3 are fulfilled. Consider a family $\{\mathcal{T}_h\}_h$ of tetrahedral meshes satisfying the usual regularity assumptions (see e.g. [8])). For any $T\in\mathcal{T}_h$, let $h_T$ be its diameter. Then $h = \max_{T\in\mathcal{T}_h} h_T$ is the mesh parameter of $\mathcal{T}_h$. For any $h$, we denote by $\mathcal{N}_h$ the set of the nodes of $\mathcal{T}_h$. We further introduce the discrete domain $\Omega_h = \bigcup_{T\in\mathcal{T}_h} T$ and its boundary ${\Gamma _\infty }h = \partial\Omega_h$. Notice that $\Omega_h$ does not fit exactly the domain $\Omega$ and its subdomains $\Omega_i$ since the latter are assumed to be at least of class $\mathcal{C}^2$. To this end, define the subtriangulation of $\mathcal{T}_h$ related to $\Omega_i$ by

$\mathcal{T}_{h, i} = \left\{ {{T\in\mathcal{T}_h}{\mathcal{N}_h\cap T \subset\overline{\Omega}_i}\ \forall i = 1, \ldots, L. } \right\}$

Then, let

$\Omega_{i, h} = \bigcup\limits_{T\in\mathcal{T}_{h, i}} T\ \mbox{and}\ \Gamma_{i, h} = \mathcal{N}_h\cap\Gamma_i.$

The following conditions state that $\mathcal{T}_h$ fits approximately $\Omega$ and its subdomains

$\bigcup\limits_{i = 1}^L \Omega_{i, h} = \Omega_h$

(5.1a)

$\Gamma_{i, h} = \Omega_{i, h}\cap\Omega_{i+1, h}\ \forall i = 1, \ldots, L-1. \label{hyp:mesh2}$

(5.1b)

These assumptions guarantee that no element has nodes in the interior of two different subdomains. In order to formulate the discrete problem, we need to extend the functions $\sigma_i$ on $\Omega_{i, h}$. To this end, assume that for any $i = 1, \ldots, L$ there is a domain $\widetilde \Omega_{i}$ such that

$\overline{\Omega_i}\subset \widetilde \Omega_{i} \ \mbox{and}\ \Omega_{i, h}\subset \widetilde \Omega_{i} .$

Without loss of generality, we may assume that $\widetilde \Omega_{i} \cap\widetilde{\Omega}_{i+2} = \emptyset$ as well as $\widetilde{\Omega}_2\cap\mathcal{V}_m = \emptyset$ for a given family of neighborhoods $(\mathcal{V}_m)_m$ of the sources $S_m$. Then, denote by $\widetilde \sigma$ an extension of $\sigma_i$ on $\widetilde \Omega_{i}$ such that $\widetilde \sigma \in W^{1, \infty}(\widetilde \Omega_{i} )$ and $\widetilde \sigma (\boldsymbol{x})\ge \sigma_{\rm min}$ for almost every $\boldsymbol{x}\in\widetilde{\Omega}_i$.

On $\mathcal{T}_h$, we introduce the standard vector space of Lagrange finite elements of type P1,

$X_h = \left\{ {{v_h\in\mathcal{C}^0(\overline{\Omega_h})}{v_{h|T}\in\mathbb{P}_1(T)\ \forall T\in\mathcal{T}_h}} \right\}$

(5.2)

where $\mathbb{P}_1(T)$ denotes the space of polynomials of degree less or equal than 1 on $T$. We further define the discretization space $V_h = X_h\cap L^2_0(\Omega_h)$ of P1 finite elements with zero mean value on $\Omega_h$.

Then the discrete problem reads: find $w_h\in V_h$ such that

$a_h(w_h, v_h) = l_h(v_h)\ \forall v_h\in V_h$

(5.3)

with

$a_h(w_h, v_h) = \sum\limits_{i = 1}^L \int_{\Omega_{i, h}} \widetilde \sigma \nabla w_h\cdot\nabla v_h\, d\boldsymbol{x}\ \mbox{and}\ l_h(v_h) = \sum\limits_{i = 1}^L \int_{\Omega_{i, h}} \tilde{\boldsymbol{F}_i}\cdot\nabla v_h + \int_{{\Gamma _\infty }h} \tilde{g} v_h \, ds$

(5.4)

where

$\tilde{\boldsymbol{F}_i} = \sum\limits_{m = 1}^M (c_m-\widetilde \sigma )\nabla \widetilde u_m \ \mbox{on}\ \widetilde \Omega_{i} \ \mbox{and}\ \tilde{g} = - \sum\limits_{m = 1}^M c_m\partial_\boldsymbol{n}\widetilde u_m \ \mbox{on}\ {\Gamma _\infty }h.$

Notice that $\widetilde u_m$ is defined in a unique way on any extension $\widetilde{\Omega}_i$ outside the neighborhood $\mathcal{V}_m$ of the source $S_m$. Nevertheless, the function $\tilde{g}$ differs from the original data $g = - \sum_{m = 1}^M c_m\partial_\boldsymbol{n}\widetilde u_m$ on ${\Gamma _\infty }$ since the normal vectors on ${\Gamma _\infty }$ and ${\Gamma _\infty }h$ are not the same.

As for the continuous problem, existence and uniqueness of the solution of (5.3) follow from Lax-Milgram's theorem since $a_h(\cdot, \cdot)$ is clearly coercive on $V_h$ due to Poincaré-Wirtinger's inequality on $\Omega_h$. Notice that the compatibility condition $l_h(1) = 0$ can be proved as in Lemma 1.

Definition 7.   For a family of discrete problems (5.3), the bilinear forms $a_h(\cdot, \cdot)$ are uniformly $V_h$-coercive if there is a constant $\alpha$ independent of $h$ such that

$a_h(v_h, v_h) \ge \alpha \| v_h \|_{H^1(\Omega_h)}^2\ \forall v_h\in V_h.$

(5.5)

Notice that the family $(a_h(\cdot, \cdot))_h$ of bilinear forms defined by (5.4) is uniformly $V_h$-coercive since the extensions $\tilde{\sigma}_i$ are uniformly bounded from below by $\sigma_{\rm min}$ and $\Omega_h$ is included in the domain $\widetilde{\Omega}\stackrel{\rm def}{ = } \bigcup_{i = 1}^L \widetilde{\Omega}_i$ for any $h$. Therefore, Poincaré-Wirtinger's inequality holds true on $\Omega_h$ with a constant independent from $h$.

Theorem 8.   Consider a regular family of meshes $(\mathcal{T}_h)_h$ fitting the geometry of $\Omega$ and its subdomains in the sense of (5.1). Assume further that $(\mathcal{T}_h)_h$ satisfies the following inverse assumption with a constant $C_{\rm inv}>0$ independent from $h$,

$\forall K\in \bigcup\limits_h \mathcal{T}_h, \ \frac{h}{h_K} \le C_{\rm inv}.$

(5.6)

Let $\tilde{w} = (\tilde{w}_i)_{i = 1, \ldots, L}\in \prod_{i = 1}^L H^2(\widetilde \Omega_{i} )$ be an extension of $w$, the solution of problem (3.17), such that $\tilde{w}_{i|\Omega_i} = w$ on $\Omega_i$ for $i = 1, \ldots, L$. Further, let $\tilde{\sigma} = (\widetilde \sigma )_{i = 1, \ldots, L}\in\prod_{i = 1}^L W^{1, \infty}(\widetilde \Omega_{i} )$ be an extension of $\sigma$ such that for $i = 1, \ldots, L$ $\widetilde{\sigma}_{i|\Omega_i} = \sigma_i$ on $\Omega_i$ and $\widetilde \sigma (\boldsymbol{x})\ge \sigma_{\rm min}$ almost everywhere.

Consider the family of discrete problems (5.3) and assume that the associated bilinear forms $a_h(\cdot, \cdot)$ are uniformly $V_h$-coercive. For any $h$, let $w_h\in V_h$ be the solution of the discrete problem (5.3).

Then, the following error estimate holds true,

$\sum\limits_{i = 1}^L \| \tilde{w}_i-w_h \|_{H^1(\Omega_{i, h})} \lesssim h \sum\limits_{m = 1}^M \left( \| \nabla \widetilde u_m\|_{H^1(\Omega\setminus\mathcal{V}_m)} + \| \partial_\boldsymbol{n}\widetilde u_m\|_{H^1({\Gamma _\infty })} \right) \\ \, \, \, \, \, \, \, \, \, + h^{3/2} \sum\limits_{m = 1}^M \left( \| \partial_\boldsymbol{n}\widetilde u_m \|_{W^{1, \infty}({\Gamma _\infty }h)} + \| \nabla\widetilde u_m \|_{H^1({\Gamma _\infty }h)} \right). \nonumber$

(5.7)

The error estimate for the approximate solution $u_h$ follows immediately from (5.7). Indeed, $u_h$ is defined by $u_h = w_h + \widetilde u$ and therefore, the error $u_h - u$ coincides with $w_h-w$ which can be estimated by (5.7).

The result of Theorem 8 relies on an abstract error estimate which states that the discretization error may be estimated by the interpolation error with respect to $X_h$ and the consistency error due to the approximation of the domain of interest.

Such an abstract error estimate has been obtained in [8] in the case of one single domain and the result carries over to our configuration of multiple subdomains taking into account the extensions of the conductivity $\sigma$ on the different subdomains $\widetilde \Omega_{i}$. Then, the interpolation error behaves as $\mathcal{O}(h)$ according to the piecewise regularity of the solution $w$ on the subdomains. The consistency error depends on the order of the geometrical approximation of the subdomains. It should be at least of the same order as the interpolation error. Following our choice of P1 Finite Elements, a polygonal discretization of the computational domain is sufficient and leads to the second term in the right hand side of (5.7) which behaves as $\mathcal{O}(h^{3/2})$. Isoparametric finite elements will fit in the case of higher order approximation.

---

5.2. Numerical validation in the multi-layer spherical model

In this paragraph, we address the numerical validation of the subtraction approach. From our motivation of studying EEG in neonates, special attention is paid to the presence of fontanels in the skull. As mentioned above, fontanels are tissues that are in the process of ossification. They are taken into account in our model through the definition of an appropriate variable skull conductivity.

Let us consider a three-layer spherical head model (see Figure 2.1) representing the brain, skull and scalp with respective radii $r_1 = 50mm$, $r_2 = 54mm$ and $r_3 = 60mm$. These dimensions correspond to a cranial perimeter of 37.7cm which is approximatively the one of a newborn child. The adopted conductivity values are $\sigma_1 = \sigma_3 = 0.33S.m^{-1}$ for the brain and the scalp and $\sigma_2 = 0.04S.m^{-1}$ for the skull [30] (unless indicated otherwise).

We consider a family of three tetrahedral meshes with decreasing mesh size $h$ (cf. Table 1). The discretization of problem (3.17) is realized as explained in Section 5.1. The approximation $u_h$ of the electric potential $u$, solution to the EEG model (2.7), is deduced from the discrete solution $w_h$ and decomposition (3.1). All simulations are executed with the software FreeFem++ [23]. The different linear systems are solved with the iterative solver GMRES (with no restart) up to a tolerance of $10^{-6}$.

Two criteria are commonly used in the numerical validation of EEG models [40]. The first, called the Relative Difference Measure (RDM), is computed as follows

$\label{RDM} \mbox{RDM} : = \left\| \dfrac{u_h}{\| u_h\|_{L^2(\Gamma_{\!\infty, h})}} - \dfrac{u_{\tiny{\mbox{ref}}} }{\| u_{\tiny{\mbox{ref}}} \|_{L^2(\Gamma_{\!\infty, h})}} \right\|_{L^2(\Gamma_{\!\infty, h})}.$

(5.8)

The second is the magnification factor (MAG) which is defined by

$\label{MAG} \mbox{MAG}: = \left| 1 - \dfrac{\| u_h\|_{L^2(\Gamma_{\!\infty, h})}}{\|u_{\tiny{\mbox{ref}}}\|_{L^2(\Gamma_{\!\infty, h})}} \right|.$

(5.9)

The RDM and MAG are error functionals with respect to a reference solution $u_{\tiny{\mbox{ref}}}$. Obviously, an ideal model leads to $\mbox{RDM} = 0$ and $\mbox{MAG} = 0$. In the case where the different tissue conductivities of the multi-layer spherical model are homogeneous, the surface potential at the scalp can be expressed as an infinite series [32,42]. We thus take $u_{\tiny \mbox{ref}}$ to be a truncated series of the exact solution which can be easily computed. We calculate the RDM and MAG for different source positions and mesh sizes. The moment of the source is $\mathbf{q} = (0, 0, J)$ with intensity $J = 10^{-6}A.m^{-2}$. The dipolar source varies along the z-axis between $10mm$ and $49mm$ where the latter position corresponds to an eccentricity of 0.98. Recall that the source eccentricity measures the relative distance to the interface brain/skull and is defined by $1-{\rm dist}(S, \Gamma_{1, h})/r_1$ for a dipole located at the point $S$. The results are reported in Figure 5.1. The accuracy of the method is very satisfying over all dipole positions and all meshes. One sees that the RDM keeps below a value of 0.9% for all the tested meshes and even under 0.5% for the two finest meshes. The amplification factor MAG keeps under 0.7% for all meshes up to an eccentricity of 0.98. Globally, both the RDM and the MAG decrease as the mesh gets finer. This validates the subtraction method in the spherical head model without fontanels.

Next, we take into account the main fontanel, i.e. the anterior fontanel situated between the frontal and parietal bones (cf Figure 1.1). The inclusion of the main fontanel in the three-layer spherical model is performed by the following function defined in the subdomain $\Omega_2$,

$\sigma_2(\boldsymbol{x}) = \sigma_{skull}+(\sigma_{\!f}-\sigma_{skull})g(\boldsymbol{x}).$

(5.10)

The parameter $\sigma_{skull}>0$ represents the (constant) conductivity of the skull outside the fontanel, wheras $\sigma_{\!f}>0$ denotes the maximal conductivity inside the fontanel. Hereafter, we take $\sigma_{skull} = 0.04S.m^{-1}$ and $\sigma_{\!f} = 0.3S.m^{-1}$ as a typical parameter set for both the fontanel conductivity and the neonatal skull conductivity that can be found in the literature [35]. The function $g$ allows to model the process of fontanel ossification and is defined by

$g(\boldsymbol{x}) = e^{-\alpha (x_1^2 + x_2^2)}$

(5.11)

depending on the parameter $\alpha>0$. In the numerical experiments hereafter, $\alpha$ is set to $10^4$ which amounts to saying that the fontanel is limited (up to 1.5%) to the region $\Omega_f: = \{\boldsymbol{x} = (x_1, x_2, x_3)\in\Omega_2\ :\ x_1^2+x_2^2\le L^2 \}$ with $L = 20mm$ (see Figure 5.2). Formula (5.10) leads to a regular function $\sigma_2$ and fits with the assumptions in Section 3.

The numerical validation is performed for different configurations. Notice that no analytical solution is available for the spherical model with fontanels. Numerical solutions are therefore compared with a numerical reference one, $u_{\tiny{\mbox{ref}}}$, computed on the very fine reference mesh $M_{\rm ref}$ (cf. Table 1). We define global errors on the whole domain $\Omega$ by

$\| u_h - u_{\tiny{\mbox{ref}}} \|^2_{H^1(\Omega)} \stackrel{\tiny \mbox{def}}{ = } \| w_h - w_{\tiny{\mbox{ref}}} \|^2_{H^1(\Omega)}.$

Figure 5.3 shows two convergence curves in logarithmic scale of the relative error in the $H^1$-norm. The graph on the left corresponds to one dipolar source at position $S = (0, 0, 40mm)$ and moment $\mathbf{q} = (0, 0, J)$ with intensity $J = 10^{-6} A.m^{-2}$. The eccentricity of the source is thus equal to $0.8$. The graph on the right shows the convergence rate for two deep sources with an eccentricity $e = 0.2$ located respectively at $S_1 = (0, 0, 10mm)$ and $S_2 = (0, 10mm, 0)$ with moments $\mathbf{q}_1 = (0, 0, J)$ and $\mathbf{q}_2 = (0, J, 0)$. The numerical results corroborate the convergence estimates of Section 5.1 that predict a theoretical convergence rate of $\tau = 1$. One can also notice that the numerical convergence rate does depend neither on the eccentricity nor on the number of sources even if the errors are larger for sources located near to the brain/skull interface.

We report in Figure 5.4 the RDM and MAG coefficients for different source positions and mesh sizes. Conclusions are the same as those obtained for the spherical head model without fontanels. The factors RDM and MAG keep under 1.5% and 0.5% respectively for all meshes and eccentricities, and decrease with the mesh size. This validates the subtraction method in the case of the spherical head model with the anterior fontanel.

One notices in Figure 5.1 and Figure 5.4 that the error increases with the eccentricity. Indeed, a careful analysis of the right hand side in the error estimate (5.7) in the case where $\sigma_1$ is constant on $\Omega_1$ shows that the error in the $H^1$-norm can be majored by $\delta^{-5/2}$ up to a constant independent from $h$ and $\delta$ (cf. [40]). Here, $\delta$ denotes the distance of the sources from the interface brain/skull,

$\delta = \max\limits_{m = 1:M} {\rm dist}(S_m, \Gamma_1)$

which is related to the eccentricity $e$ by

$e = 1 - \frac{\delta}{r_1}.$

This deterioration in the precision of the simulated data is inherent to the numerical method and dependent on the mesh that is used in the simulations.

---

6. Numerical discussion of the EEG model in neonates

In this section, we investigate numerically the effects of fontanels and their conductivity on the electric potential measured at the scalp with following question in mind: is it important or not to consider fontanels in the EEG model in neonates?

---

6.1. A realistic head model

We consider a realistic head model of a healthy full-term newborn obtained from coregistration of MR and CT images of the Amiens' hospital database (see Figure 6.1). The diameter of the computational domain is about 12cm. Obtaining a realistic model was necessary for our numerical study. Nevertheless, neonatal cerebral image segmentation is a challenging task due to inherent difficulties such as low signal to noise ratio and partial volume effects. Segmentation of adult MR images has been extensively addressed by many researchers, and several methods for automatic segmentation have been developed in the case of healthy adults. However, due to insufficient anatomical similarity between neonatal and adult's images, these methods generally fail to segment accurately cerebral images of neonates. Extracting cranial bones and fontanels from neonatal MR images is also difficult and even unreliable in some cases. In this study, one has used the coregistration of MR and CT images of one healthy male neonate of 42 weeks gestational age at the time of the scan to construct a realistic volume conductor head model containing five different newborn's head compartments: brain ($\Omega_1$), cerebrospinal fluid ($\Omega_2$), skull ($\Omega_3$), fontanels ($\Omega_{3, f} \subset \Omega_{3}$) and scalp ($\Omega_4$). The T2-weighted MR image was acquired by a 1.5T GE MR scanner using pulse sequences of TR = 4500ms and TE = 149ms at original resolution of $1\times 0.47\times 0.47\ {\rm mm}^3$. The CT image was recorded in a separate session using a LightSpeed 16, GE medical system with a spatial resolution of $0.35\times 0.35\times 0.63\ {\rm mm}^3$. In a preprocessing step, the images were resampled to an isotropic 0.94 $mm^3$ resolution and interactively clipped to remove excess signal and extra parts of the neck. The CT and MR images were co-registered in SPM8 (Wellcome Trust Center for Neuroimaging, UCL, London) by means of a sequential approach [36]. The brain and CSF were segmented with help of the atlas based automatic segmentation in SPM8 from the MR image. The cranial bones were extracted from the CT image using a simple thresholding. Subsequently, the fontanel and sutures were delimited manually as gaps between the extracted bones. The scalp was delineated by removing the already segmented compartment from the MR image. A set of manual and automatic corrections were applied to the segmented tissues in order to eliminate unavoidable errors in neonatal cerebral image segmentation. The Iso2mesh toolbox [15] was used to construct a tetrahedral mesh of 1.1 mm resolution containing 590 878 elements and 108 669 nodes. The mesh properties of the model are reported in Table 2. Figure 6.2 shows the different compartments of the head model as well as their three-dimensional reconstruction. The resulting neonatal head model includes more precise information on neonatal skull geometry (especially concerning the size and the exact position of fontanels) than previous models [18]. It also considers the effects of the temporal fontanels.

---

6.2. Effects of fontanels

We compare different EEG forward models. We use a panel of conductivities which are found in literature (e.g. [35,30,18,3]). To each couple $(\sigma_{\!f}, \sigma_{skull})$ in the parameter set corresponds an EEG model, whereas the choice $\sigma_{\!f} = \sigma_{skull}$ yields the model without fontanels. In order to quantify the influence of the values $(\sigma_{\!f}, \sigma_{skull})$ on the measured electric potential at the scalp, we compute the factors RDM and MAG (cf. formula (5.8) and (5.9)) with respect to $u_h$ and $u_{\tiny{\mbox{ref}}}$ which are, respectively, the solutions to the fontanels and no-fontanel models.

The different conductivity values are fixed to $\sigma_1 = \sigma_4 = 0.33S.m^{-1}$ for brain and scalp and $\sigma_2 = 1.8S.m^{-1}$ for CSF. We compute the RDM and MAG factors for different couples $(\sigma_{\!f}, \sigma_{skull})$ which are varying around a reference conductivity of $\sigma_{\!f} = 0.3S.m^{-1}$ and $\sigma_{skull} = 0.04S.m^{-1}$. The uncertainty in head tissue conductivities is estimated at $\pm 25%$ which leads to $\sigma_{\!f} \in [0.225, 0.375]$ and $\sigma_{skull} \in [0.03, 0.05]$. A single dipole source is placed in $\Omega_1$ at a distance of $5$mm of the interface brain/CSF. Its moment is $(0, 0, J)$ with intensity $J = 2\sqrt{2} \ 10^{-6} A.m^{-2}$. Figure 6.3 represents the isolines of RDM and MAG for varying fontanel and skull conductivities. The factors RDM and MAG behave like functions of the ratio $\sigma_{\!f}/\sigma_{skull}$. The difference of electrical potential between models with or without fontanels is maximal (RDM = 10% and MAG = 15%) when the ratio $\sigma_{\!f}/\sigma_{skull}$ is the highest (equal to $0.375/0.03 = 12.5$). At a fixed conductivity of the fontanels, RDM and MAG decrease as the skull conductivity increases, but keep significant even at the lowest ratio (RDM = 6% and MAG = 8%). This numerical study shows that the presence of fontanels impacts the EEG measurements in neonates through the ratio $\sigma_{\!f}/\sigma_{skull}$ of tissue conductivities. The most pronounced effect occurs for higher ratio. These observations indicate that it would be judicious to include fontanels in the EEG forward model. It also says that uncertainty in conductivity plays an important role. The following numerical sensitivity analysis gives further insight into this issue.

---

6.3. Numerical sensitivity analysis

We study numerically how a slight variation of the conductivity affects EEG measurements. We are interested in analyzing conductivity perturbations in neonatal skull including fontanels.

The sensitivity $u^1$ of the potential in a given direction $\mu$ is computed as the solution of equation (4.5). Here, we consider a perturbation of the fontanel conductivity which amounts to take the characteristic function on $\Omega_{\!f}$, $\mu = 1_{\Omega_{\!f}}$, as direction.

Figure 6.4 compares the sensitivity of two sources of same moment $\mathbf{q} = (0, J, J)$ with $J = 2\sqrt{2} \ 10^{-6} A.m^{-2}$ that are localized at different positions. The moment is directed to the anterior fontanel and the sources are situated at a distance of, respectively, 5mm and 15mm from the interface brain/CSF. In both cases, the sensitivity is localized to an area above the anterior fontanel, but a diffusion effect appears in the case of the deeper source. Next, we consider a source that is located at a distance of $15$mm from the interface brain/CSF. We compare the sensitivity for the moments $\mathbf{q}_1 = (J, J, 0)$ and $\mathbf{q}_2 = (J, 0, J)$ with $J = 2\sqrt{2} \ 10^{-6} A.m^{-2}$. Figure 6.5 shows that the area of maximal sensitivity depends on the direction of the moment, but is still localized near the anterior fontanel. In order to illustrate the impact of fontanels on the measured potential independently from source location and orientation, we chose to compute the sensitivity for a deep source far away from the anterior fontanel (see Figure 6.6). Even if the absolute value of the quantity is about one order of magnitude smaller than in the near-interface configuration, one may notice that the impact is still significant near the fontanels.

The sensitivity analysis of this section provides information on those scalp areas on which the electrical potential is affected by small variations of the fontanel's conductivity. These variations appear naturally from one patient to the other and even for the same patient at different ages. The results corroborate the conclusions of the previous section that the fontanel and skull conductivity values impact the EEG forward solution. This effect seems to be localized near the fontanels and depends on the eccentricity and orientation of the dipolar sources.

---

7. Concluding remarks

In this paper, we have proposed an EEG forward problem for neonates, that amounts to solve an elliptic problem with a singular source term in a an inhomogeneous medium. The considered domain is made up of disjoint coated regions, the source term is a finite linear combination of dipoles and the conductivities are functions of the position in each region. This last assumption allows to take into account the presence of fontanels in the skull and their ossification process. In this context, a rigorous mathematical framework has been proposed. In order to deal with the singularity of the source term, we apply a subtraction method and give a proof for existence and uniqueness of the weak solution in the context of variable conductivities. Convergence estimates are obtained for standard finite elements of type P1.

Numerical validations were performed. Firstly, the numerical validation has been carried out for the three-layer spherical head model with or without fontanels. For a given source term, the error of the regular part of the solution has been computed in the $H^1$-norm on the whole computational domain for different meshes. The numerical convergence rate coincides with the theoretical results for any tested configuration. Classical error functionals as RDM and MAG factors have been computed for sources with different eccentricities. Conforming to theoretical results, the error increases if the source is placed near the interface between the brain and the neighboring layer.

Nevertheless, in all tested configurations the error keeps lower than 1.5% even for the coarsest mesh and the most eccentric source position. Higher quadrature rules could help to overcome the influence of the singular behavior of the source term.

Secondly, with the motivation of answering clinical questions, we have studied the influence of fontanels and uncertainty in fontanel and skull conductivities on the electrical potential values at the scalp using a realistic neonatal head model. This model is provided by Inserm U1105 (Amiens' hospital) database. The RDM and MAG factors between models with or without fontanels depend on the ratio of the fontanel conductivity over the skull conductivity. The difference is more pronounced for a higher ratio. With regard to EEG source reconstruction, these numerical observations attest the importance of considering the presence of the fontanels in the EEG forward model for neonates. It shows also that uncertain conductivity values impact the EEG forward solution. These conclusions are confirmed by a mathematical and numerical sensitivity analysis. Furthermore, this study provides useful information on areas where the electrical potential is the most sensitive to a variation of conductivity. The support of the sensitivity function depends on the source characteristics (position and moment) and is localized to an area above the fontanels.

Our findings are comparable with those of Azizollahi et al. [3] who studied the influence of different tissue conductivities including the fontanels. Their approach was based on a model of distributed sources. We are actually working together with the group of GRAMFC at Amiens' hospital for a detailed study of these new aspects.

The analysis of the EEG forward model is an essential preliminary step to the resolution of the corresponding inverse source problem. In the head model of adults, theoretical and numerical results exist [9,13,14,16,29,10,4]. For the neonatal head model, identifiability and stability results for the inverse EEG source problem have been obtained in parallel to the present work [11]. A numerical study that should help to understand the impact of the presence of fontanels on the source reconstruction in neonates, especially for epileptogenic sources, is underway.

---

Acknowledgments

The present work has been realized as part of the MIFAC project receiving financial support from the region Picardie (now Hauts-de-France).

---

Appendix A.

Proof of Theorem 5.  Let $w^1\in V$ be the solution of the following variational problem,

$\int_{\Omega}\sigma\nabla w^1\cdot\nabla v\, d\boldsymbol{x} = -\int_{\Omega}\mu\nabla w\cdot\nabla v\, d\boldsymbol{x} \label{Vsens2}\\ \, \, \, \, \, \, +\sum\limits_{m = 1}^M\left(\int_{\Omega}(c_m-\sigma)\nabla \widetilde u_m^1\cdot\nabla v\, d\boldsymbol{x} -\int_{{\Gamma _\infty }}c_m\partial_{\boldsymbol{n}}\widetilde u_m^1v\, d\boldsymbol{x}\right)\nonumber \\ \, \, \, \, \, \, +\sum\limits_{m = 1}^M\left(\int_{\Omega}(p_m-\mu)\nabla \widetilde u_m \cdot\nabla v\, d\boldsymbol{x} -\int_{{\Gamma _\infty }}p_m\partial_{\boldsymbol{n}}\widetilde u_m v\, d\boldsymbol{x}\right), \nonumber$

(A.1)

for all $v\in V$, where $\widetilde u_m^1 = D_\mu\widetilde u_m(\cdot, \sigma)$ and $\widetilde u_m = \widetilde u_m(\cdot, \sigma)$. We aim to prove that $w^1$ is the Gâteaux-derivative with respect to the conductivity $\sigma$ of the variational solution $w$ of problem (3.17).

First of all, notice that $w^1$ is well defined by (4.5) according to Lax-Milgram's theorem since the right hand side of (4.5) clearly defines a continuous linear form on $H^1(\Omega)$ and satisfies the compatibility condition $l(1) = 0$. Indeed, both boundary integrals in (4.5) vanish for $v\equiv 1$ due to the assertion of Lemma 1, taking into account that $\widetilde u_m^1 = -\frac{p_m}{c_m}\widetilde u_m$.

Then, in order to investigate the Gâteaux derivative of $w$ with respect to $\sigma$, we consider $\mu\in\mathcal{P}$ with $\|\mu\|_{\infty} = 1$ satisfying, from assumptions, $\sigma+h\mu\in\mathcal{P}_\text{adm}$ for any $h\in [-h_0, h_0]$.

In the sequel we shall omit the dependence of $w$ and $\widetilde u_m$ on the parameters $\sigma$ and $\mu$ for better reading and set

$w \stackrel{\rm def}{ = } w(\cdot, \sigma) \ {\rm and}\ w_h \stackrel{\rm def}{ = } w(\cdot, \sigma+h\mu)$

as well as

$\widetilde u_m \stackrel{\rm def}{ = } \widetilde u_m(\cdot, \sigma)\ {\rm and}\ \widetilde u_{m, h}\stackrel{\rm def}{ = } \widetilde u_m(\cdot, \sigma+h\mu)$

Now, recall that $w_h$ and $w$ are the respective solutions in $V$ of

$\label{Vsigh} \begin{array}{lll} \displaystyle \int_{\Omega}(\sigma+h\mu)\nabla w_h\cdot\nabla v\, d\boldsymbol{x}& = &\displaystyle \sum\limits_{m = 1}^M (\int_{\Omega}\left((c_m+h p_m)-(\sigma+h\mu)\right)\nabla \widetilde u_{m, h}\cdot\nabla v\, d\boldsymbol{x}\\ &-& \displaystyle \int_{{\Gamma _\infty }} (c_m+h p_m) \partial_{\boldsymbol{n}}\widetilde u_{m, h} v \, ds), \end{array}$

(A.2)

and

$\int_{\Omega} \sigma\nabla w\cdot\nabla v\, d\boldsymbol{x} = \sum\limits_{m = 1}^M\left(\int_{\Omega}(c_m- \sigma)\nabla \widetilde u_m \cdot\nabla v\, d\boldsymbol{x} -\int_{{\Gamma _\infty }}c_m\partial_{\boldsymbol{n}}\widetilde u_m v \, ds \right),$

(A.3)

for all $v \in V$. The assumption on $\sigma$ and $\mu$ implies that $\sigma+h\mu > \sigma_{\rm min}$ for any $h\in [-h_0, h_0]$ and therefore the bilinear form in (A.2) is $V$-elliptic with a constant independent from $h$. The right-hand side of (A.2) defines a continuous linear form on $H^1(\Omega)$ and we get the following estimate from similar arguments as for (3.18) in Theorem 1:

$\| w_h \|_{H^1(\Omega)} \lesssim \sum\limits_{m = 1}^M \left( \| \nabla \widetilde u_{m, h} \|_{L^2(\Omega\setminus\mathcal{V}_m)} + \| \partial_\boldsymbol{n}\widetilde u_{m, h}\|_{L^2({\Gamma _\infty })} \right).$

(A.4)

In order to identify the Gâteaux derivative of $w$, we introduce the differential quotients

$w_h^1 = \, \frac{w_h - w}{h} \mbox { and } \, \widetilde u^1_{m, h} = \frac{\widetilde u_{m, h} - \widetilde u_m}{h}.$

Subtracting (A.3) from (A.2) and dividing by $h$ leads to

$\label{limwh1} \int_{\Omega} \sigma\nabla w_h^1 \cdot\nabla v\, d\boldsymbol{x} = - \int_{\Omega} \mu \nabla w_h \cdot\nabla v\, d\boldsymbol{x} \\ \, \, \, \, \, \, + \sum\limits_{m = 1}^M \left(\int_{\Omega} (c_m- \sigma)\nabla \widetilde u^1_{m, h} \cdot\nabla v\, d\boldsymbol{x} - \int_{{\Gamma _\infty }} c_m\partial_{\boldsymbol{n}} \widetilde u^1_{m, h} \, v \, ds \right) \nonumber\\ \, \, \, \, \, \, + \sum\limits_{m = 1}^M \left(\int_{\Omega} (p_m- \mu) \nabla \widetilde u_{m, h}\cdot\nabla v\, d\boldsymbol{x} - \int_{{\Gamma _\infty }} p_m \partial_{\boldsymbol{n}} \widetilde u_{m, h} v \, ds\right). \nonumber$

(A.5)

We compare the above formulation for the differential quotient $w_h^1$ with the variational formulation (4.5) for $w^1$:

$\label{wh1_w1} \int_{\Omega} \sigma\nabla \left( w_h^1 -w^1 \right) \cdot\nabla v\, d\boldsymbol{x} = - \int_{\Omega} \mu \nabla \left( w_h-w \right) \cdot\nabla v\, d\boldsymbol{x}\\ \, \, \, \, \, \, \, \, \, \, \, + \sum\limits_{m = 1}^M \left(\int_{\Omega} (c_m- \sigma)\nabla \left( \widetilde u^1_{m, h} -\widetilde u^1_m \right) \cdot\nabla v\, d\boldsymbol{x} - \int_{{\Gamma _\infty }} c_m\partial_{\boldsymbol{n}} \left( \widetilde u^1_{m, h}-\widetilde u^1_m\right) \, v \, ds \right) \nonumber\\ \, \, \, \, \, \, \, \, \, \, \, + \sum\limits_{m = 1}^M \left(\int_{\Omega} (p_m- \mu) \nabla \left( \widetilde u_{m, h} -\widetilde u_m\right)\cdot\nabla v\, d\boldsymbol{x} - \int_{{\Gamma _\infty }} p_m \partial_{\boldsymbol{n}} \left( \widetilde u_{m, h} -\widetilde u_m\right)\, v \, ds\right). \nonumber$

(A.6)

Notice that the integrals over $\Omega$ in the second and third term vanish on $\mathcal{V}_m$ since $\sigma_{\vert \mathcal{V}_m}\equiv c_m$ and $\mu_{\vert \mathcal{V}_m}\equiv p_m$. Taking $v = w_h^1 -w^1$ in (A.6), we get the following estimate from classical inequalities in variational theory,

$\label{estim_wh1w1} \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \| w_h^1 -w^1\|_{H^1(\Omega)}\\ \, \, \, \, \lesssim \| \nabla ( w_h-w )\|_{L^2(\Omega)}\nonumber \\ \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, + \sum\limits_{m = 1}^M \left( \left\|\nabla \left( \widetilde u^1_{m, h} -\widetilde u^1_m \right) \right\|_{L^2(\Omega\setminus\mathcal{V}_m)} + \left\| \partial_\boldsymbol{n} \left( \widetilde u^1_{m, h} -\widetilde u^1_m \right) \right\|_{L^2({\Gamma _\infty })} \right) \nonumber\\ \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, + \sum\limits_{m = 1}^M \left( \left\|\nabla \left( \widetilde u_{m, h} -\widetilde u_m \right) \right\|_{L^2(\Omega\setminus\mathcal{V}_m)} + \left\| \partial_\boldsymbol{n} \left( \widetilde u_{m, h} -\widetilde u_m \right) \right\|_{L^2({\Gamma _\infty })} \right) \nonumber$

(A.7)

The second and third term in the right-hand side of (A.7) are of order $h$ and can be majored by the results of Lemma 2 hereafter. In order to estimate the first term, notice that $w_h-w$ is related to the differential quotient $w_h^1$ by $w_h-w = hw_h^1$. Hence, $h^{-1}(w_h-w)$ can be estimated by the norm of the linear form on the right-hand side of (A.5) and we get

$\label{estim_whw} \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, h^{-1} \| w_h -w\|_{H^1(\Omega)} \nonumber\\ \, \, \, \, \, \, \, \, \, \, \, \lesssim \|\nabla w_h\|_{L^2(\Omega)} + \sum\limits_{m = 1}^M \left( \left\|\nabla \widetilde u^1_{m, h} \right\|_{L^2(\Omega\setminus\mathcal{V}_m)} + \left\| \partial_\boldsymbol{n} \widetilde u^1_{m, h} \right\|_{L^2({\Gamma _\infty })} \right) \\ \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, + \sum\limits_{m = 1}^M \left( \left\|\nabla \widetilde u_{m, h} \right\|_{L^2(\Omega\setminus\mathcal{V}_m)} + \left\| \partial_\boldsymbol{n} \widetilde u_{m, h} \right\|_{L^2({\Gamma _\infty })} \right). \nonumber$

(A.8)

The right hand side in (A.8) is obviously bounded when $h$ tends to zero. Indeed, we have $\lim_{h\to 0} \widetilde u_{m, h}^1 = \widetilde u_m^1$ and $\lim_{h\to 0}\widetilde u_{m, h} = \widetilde u_m$, and $\nabla w_h$ is bounded in terms of $\widetilde u_{m, h}$ according to (A.4). Multiplying (A.8) by $h$ shows that $\| w_h -w\|_{H^1(\Omega)}$ is of order $h$. Consequently, (A.7) becomes

$\| w_h^1 -w^1\|_{H^1(\Omega)} \lesssim h \sum\limits_{m = 1}^M \left( \left\|\nabla \widetilde u_m \right\|_{L^2(\Omega\setminus\mathcal{V}_m)} + \left\| \partial_\boldsymbol{n} \widetilde u_m \right\|_{L^2({\Gamma _\infty })} \right) .$

This proves the strong convergence of the sequence $(w_h^1)_h$ to $w^1$ in $H^1(\Omega)$.

It remains to show that $D_\mu w$ belongs to $\mathcal{L}(L^{\infty}, V)$. For fixed $\mu$, $D_\mu w(\cdot, \sigma)$ is defined by the solution of (4.5). But the right-hand side of (4.5) is linear in $\mu$. This is obvious for the first and the third term since $w$ and $\widetilde u_m$ are independent from $\mu$ and $p_m = \mu(S_m)$ is linear in $\mu$. The second term depends on $\mu$ only via the Gâteaux derivative $\widetilde u_m^1$ of $\widetilde u_m$. According to Proposition 1, we have $\widetilde u_m^1 = -\frac{\mu(S_m)}{\sigma(S_m)}\widetilde u_m$ which is a linear expression in $\mu$.

In order to prove that the linear application $\mu\mapsto D_\mu w$ is continuous from $\mathcal{P}$ to $V$, we estimate $w^1$ with the help of formulation (4.5). Taking $v = w^1$, we get

$\|w^1\|_{H^1(\Omega)} \lesssim \|\mu\|_{\infty}\|\nabla w\|_{L^2(\Omega)} +\sum\limits_{m = 1}^M |\mu(S_m)| \left( \|\nabla \widetilde u_m\|_{L^2(\Omega\setminus\mathcal{V}_m)} + \|\partial_{\boldsymbol{n}}\widetilde u_m\|_{L^2({\Gamma _\infty })} \right)$

using again that $\widetilde u_m^1 = -\frac{\mu(S_m)}{\sigma(S_m)}\widetilde u_m$. This yields the continuity of $D_\mu w$ with respect to $\mu$ and proves that $w(\cdot, \sigma)$ is Gâteaux differentiable with respect to the conductivity $\sigma$.

Lemma 2. Let $\widetilde u_m = \widetilde u_m(\cdot, \sigma)$ and $\widetilde u_{m, h} = \widetilde u_m(\cdot, \sigma+h\mu)$ be given by (4.1) and (4.2), respectively. Under the assumptions of Theorem 5, the following estimates hold true:

$\left \|\nabla \left( \widetilde u_{m, h}-\widetilde u_m\right)\right \|_{L^2(\Omega\setminus\mathcal{V}_m)} \lesssim h||\nabla \widetilde u_m||_{L^2(\Omega\setminus\mathcal{V}_m)}, \label{estim1:a}\\ $

(A.9a)

$\left \|\partial_{\boldsymbol{n}}\left( \widetilde u_{m, h}-\widetilde u_m\right)\right \|_{L^2({\Gamma _\infty })} \lesssim h||\partial_{\boldsymbol{n}}\widetilde u_m||_{L^2({\Gamma _\infty })}. \label{estim1:b}$

(A.9b)

Further, let $\widetilde u_{m, h}^1 = \frac{\widetilde u_{m, h}-\widetilde u_m}{h}$ and denote by $\widetilde u_m^1 = D_\mu\widetilde u_m(\cdot, \sigma)$ the Gâteaux-derivative of $\widetilde u_m$ at $\sigma$ in the direction $\mu$. Under the assumptions of Theorem 5, the following estimates hold true:

$\left \|\nabla \left( \widetilde u_{m, h}^1-\widetilde u_m^1 \right)\right \|_{L^2(\Omega\setminus\mathcal{V}_m)} \lesssim h||\nabla \widetilde u_m||_{L^2(\Omega\setminus\mathcal{V}_m)}, \label{estim2:a}\\ $

(A.10a)

$\left \|\partial_{\boldsymbol{n}}\left(\widetilde u_{m, h}^1-\widetilde u_m^1\right)\right \|_{L^2({\Gamma _\infty })} \lesssim h||\nabla \widetilde u_m||_{L^2({\Gamma _\infty })}.$

(A.10b)

Proof. From the definition of $\widetilde u_m$ and $\widetilde u_{m, h}$ we get

$\widetilde u_{m, h}(\boldsymbol{x}) - \widetilde u_m(\boldsymbol{x}) = \frac{1}{4\pi} \left( \frac{1}{c_m+hp_m} - \frac{1}{c_m} \right) \frac{\boldsymbol{q}_m\cdot(\boldsymbol{x}-S_m)}{|\boldsymbol{x}-S_m|^3} = -\frac{hp_m}{c_m+hp_m} \widetilde u_m(\boldsymbol{x}).$

(A.11)

Integration of the gradient of the above expression over $\Omega\setminus\mathcal{V}_m$ yields (A.9a) since $|p_m| \le \| \mu\|_\infty = 1$ and $|c_m+ hp_m| \ge \sigma_{\rm min}$. (A.9b) follows by integration of the normal derivative.

Next, recall that $\widetilde u_m^1 = -\frac{p_m}{c_m}\widetilde u_m$ according to Proposition 1. Together with identity (A.11), we thus get

$\widetilde u_{m, h}^1 - \widetilde u_m^1 = \left( -\frac{p_m}{c_m+hp_m} + \frac{p_m}{c_m} \right) \widetilde u_m = h \frac{p_m^2}{c_m(c_m+hp_m)}\widetilde u_m.$

(A.12)

The boundedness of $|p_m|$ and $c_m$ (resp. $c_m+hp_m$) yields (A.10a) by integration of the gradient of (A.12) over $\Omega\setminus\mathcal{V}_m$. (A.10b) follows in the same way.

---