EEG in neonates: Forward modeling and sensitivity analysis with respect to variations of the conductivity.

The paper is devoted to the analysis of electroencephalography (EEG) in neonates. The goal is to investigate the impact of fontanels on EEG measurements, i.e. on the values of the electric potential on the scalp. In order to answer this clinical issue, a complete mathematical study (modeling, existence and uniqueness result, realistic simulations) is carried out. A model for the forward problem in EEG source localization is proposed. The model is able to take into account the presence and ossification process of fontanels which are characterized by a variable conductivity. From a mathematical point of view, the model consists in solving an elliptic problem with a singular source term in an inhomogeneous medium. A subtraction approach is used to deal with the singularity in the source term, and existence and uniqueness results are proved for the continuous problem. Discretization is performed with 3D Finite Elements of type P1 and error estimates are proved in the energy norm (H¹-norm). Numerical simulations for a three-layer spherical model as well as for a realistic neonatal head model including or not the fontanels have been obtained and corroborate the theoretical results. A mathematical tool related to the concept of Gâteau derivatives is introduced which is able to measure the sensitivity of the electric potential with respect to small variations in the fontanel conductivity. This study attests that the presence of fontanels in neonates does have an impact on EEG measurements.


(Communicated by Eugene Kashdan)
Abstract. The paper is devoted to the analysis of electroencephalography (EEG) in neonates. The goal is to investigate the impact of fontanels on EEG measurements, i.e. on the values of the electric potential on the scalp. In order to answer this clinical issue, a complete mathematical study (modeling, existence and uniqueness result, realistic simulations) is carried out. A model for the forward problem in EEG source localization is proposed. The model is able to take into account the presence and ossification process of fontanels which are characterized by a variable conductivity. From a mathematical point of view, the model consists in solving an elliptic problem with a singular source term in an inhomogeneous medium. A subtraction approach is used to deal with the singularity in the source term, and existence and uniqueness results are proved for the continuous problem. Discretization is performed with 3D Finite Elements of type P1 and error estimates are proved in the energy norm (H 1 -norm). Numerical simulations for a three-layer spherical model as well as for a realistic neonatal head model including or not the fontanels have been 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 quasistatic Maxwell equations where the time derivatives are neglected [22,17]. In terms of the electric field E and the magnetic field H, this yields Here, ρ is the charge density, ε and µ are, respectively, the electric permittivity and magnetic permeability, and J is the electric current density. Following Ohm's law, the current density splits into two terms, where J s is the density of the impressed neural currents and σ 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.
Now, consider a bounded regular domain Ω ⊂ R 3 with boundary ∂Ω. Taking the divergence of (2.1c) together with Ohm's law (2.2) yields the following transmission equation for the electric potential u in Ω 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 ∈ Ω 1 and its moment q m which is a vector of R 3 . The points (S m ) m are assumed to be mutually distinct and the vectors q m are non vanishing, i.e. S m = S p ∀m = p and q m = 0 ∀m, p ∈ {1, ..., M }. (2.5) The current density J s thus reads where δ Sm denotes the delta distribution at S m . The right hand side of (2.4) is then given by 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 Ω into L open subdomains (Ω i ) i=1,...,L (see Figure 2.1) such that We assume the subdomains to be nested which means that Then, for i = 1, . . . , L − 1, we denote by Γ i the interface between Ω i and Ω i+1 , and assume that (Γ i ) i are closed regular surfaces. Let n i be the unit normal vector on Γ i oriented towards the exterior of Ω i . Notice that in a more general setting, we could allow interfaces Γ ij between two arbitrary subdomains Ω i and Ω j provided the interfaces are not intersecting. We finally denote by Γ ∞ = ∂Ω the exterior boundary of the whole domain Ω. This configuration includes the classical spherical model of three concentric spheres representing brain, skull and scalp (see Figure 2.1). Now, let σ i def = σ |Ωi (i = 1, . . . , L) denote the conductivity of the subdomain Ω i . In the head model of an adult, the different layers are assumed to be homogeneous and isotropic, and therefore each σ 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 σ i that are functions of the position, i.e. σ i = σ i (x) in Ω i . where the source term F is given by (2.6). Since F vanishes identically in a neighborhood of the interfaces Γ i , u satisfies the transmission conditions 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 (Ω) is not possible since the source term F given by (2.6) belongs to H s (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ũ = −G F where G is the fundamental solution of the Laplacian. Then, the function w = u −ũ satisfies −σ i ∆w = 0 in the subdomains Ω i and non-homogeneous transmission conditions with regular right hand sides on the interfaces Γ 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 σ 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ũ which contains the singularity and a regular lifting w u =ũ + w, withũ = where c m = σ 1 (S m ). It can be obtained by convolution of the fundamental solution of the Laplace equation with the right hand side 1 c m q m · ∇δ Sm as follows We see that the potentialũ 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 ∈ {1, . . . , M }, the quantity ∇ · (σ∇ũ m ) is well defined on Ω by Indeed, both terms on the right hand side of the above identity are well defined as distributions on Ω since σ − c m vanishes identically on V m andũ m is regular outside ∪ m V m . Therefore, we get −∇ · (σ∇w) = −∇ · (σ∇(u −ũ)) = F + Problem (3.5) can be reformulated on each subdomain Ω i provided the conductivity σ is piecewise regular. To this end, assume that σ i ∈ W 1,∞ (Ω i ). We may notice that the first term on the right hand side of the following identity vanishes, since ∆ũ m is a distribution with a single point support according to (3.3) and vanishes outside the neighborhood V m , We thus get On Ω i , i = 1, the potentialsũ m are regular and we see from a direct computation that w is solution of the following problem with transmission conditions Recall that the right hand side of (3.7) is well defined and vanishes in any neighborhood V m of the sources, since σ 1 is constant on V m .

Variational formulation.
In this section, we give a variational formulation of problem (3.7)-(3.8) for the auxiliary function w. Let v ∈ H 1 (Ω). Multiplying (3.7) by v |Ωi , integrating over Ω i and summing up, we obtain the following formulation from Green's formula and the boundary and transmission conditions, Since ∇(σ i − c m ) = ∇σ i , we can show with the help of (3.6) that (3.9) is equivalent toˆΩ 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(·, ·) Notice that a(·, ·) and l(·) are well defined on H 1 (Ω) whenever the conductivity σ belongs to L ∞ (Ω) 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 Condition (3.13) follows from the following classical lemma given in [16] that we recall for the convenience of the reader. 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 (Ω) which does not contain any constant other than zero, on which the Poincaré-Wirtinger inequality holds true, In the sequel, we write a b if there is a constant C > 0 independent from the quantities a and b such that a ≤ Cb. Theorem 1. Let σ ∈ L ∞ (Ω) be such that 0 < σ min ≤ σ(x) ≤ σ max for almost any x ∈ Ω, where σ min and σ max are two given positive constants. Assume further that σ |Vm is constant for any m = 1, . . . , M and denote by c m the value of σ on V m ⊂ Ω 1 . Let the bilinear form a(·, ·) and the linear form l(·) be given by (3.11) and (3.12), respectively. Then the variational problem has exactly one solution w ∈ V . Moreover, the following estimate holds true, Proof. It follows from standard arguments in variational theory that the bilinear form a(·, ·) is continous on where C T is the continuity constant of the trace operator from H 1 (Ω) to H 1/2 (Γ ∞ ). This proves that the linear form is continuous on H 1 (Ω) and thus on V . The Lax-Milgram theorem then yields existence and uniqueness of a function w ∈ V such that |Ω|´Ω v dx is the mean value of v. Since l(1) = 0 according to the compatibility condition (3.13) and 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(·) in the variational formulation (3.17) is well defined if and only if σ 1 is constant in a neighborhood of the source positions S m . Indeed, if this would not be the case, the term (c m − σ 1 )∇ũ m does no longer belong to (L 2 (Ω 1 )) 3 due to the singularity of ∇ũ m that behaves as 1 |x − S m | 3 , even if σ 1 is regular enough such that the value c m = σ 1 (S m ) makes sense.
The following theorem states the global H 2 -regularity of the variational solution of (3.17) in the subdomains Ω i . (3.20) The proof of Theorem 3 relies on standard techniques for elliptic partial differential equations. Indeed, we may notice that on each Ω i , the variational solution w satisfies the partial differential equation According to the assumptions on σ i , the function f i belongs to L 2 (Ω i ). Hence, classical arguments for partial differential equations with variable coefficients apply and yield interior regularity in each Ω i , se for example [7,19]. H 2 -regularity up to the boundary of the subdomains Ω i follows since the boundary Γ ∞ and the interfaces Γ i as well as the Neumann data σ∂ nũ 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 Definition 4. Let F : X → Y be an application between two Banach spaces X and Y . Let if the limit exists. If D µ F (σ) exists for any direction µ ∈ X and if the application Now, let (V m ) m be a fixed family of neighborhoods of the sources such that for as well as the (open) subset of admissible conductivities. Here, σ min and σ max are two fixed positive constants. According to Theorem 1, problem (3.10) with conductivity σ ∈ P adm admits a unique solution w(·, σ) ∈ V . The aim of this section is to prove differentiability of w with respect to σ and to identify its (Gâteaux) derivative in a given direction µ.
Since w depends on the singular potentialũ, we first analyze the derivative ofũ with respect to σ. To this end, recall that is the solution of the Poisson equation where c m = σ(S m ). Now, consider an arbitrary direction µ ∈ P and assume that σ + hµ belongs to P adm for small values of h. By definition, we havẽ The following proposition states thatũ m is Gâteaux differentiable at an interior point σ ∈ P and identifies its Gâteaux derivative: Proposition 1. Let σ ∈ P adm and h 0 > 0 such that σ + hµ ∈ P adm for any h ∈ [−h 0 , h 0 ] and any µ ∈ P with µ ∞ = 1. Then,ũ m (·, σ) is Gâteaux differentiable at σ and the Gâteaux derivative D µũm (·, σ) ofũ m at σ in the direction µ reads with p m = µ(S m ) and c m = σ(S m ).

Proof.
A straightforward computation of the differential quotient yields and (4.3) follows. The right-hand side of (4.3) is obviously linear and continuous in µ since p m = µ(S m ).
Theorem 5. Let σ ∈ P adm and h 0 > 0 such that σ + hµ ∈ P adm for any h ∈ [−h 0 , h 0 ] and any µ ∈ P with µ ∞ = 1. Then the solution w(·, σ) of (3.10) is Gâteaux differentiable with respect to σ and the Gâteaux derivative of w at σ in the direction µ ∈ P is the unique solution of the following variational problem : find . 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.
Thus, the Gâteaux derivative u 1 def = u (σ; µ) of the electric potential u at σ in the direction µ 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 Ω i are assumed to be regular. The discrete problem will thus be set on discretized domains Ω 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 {T h } h of tetrahedral meshes satisfying the usual regularity assumptions (see e.g. [8])). For any T ∈ T h , let h T be its diameter. Then h = max T ∈T h h T is the mesh parameter of T h . For any h, we denote by N h the set of the nodes of T h . We further introduce the discrete domain Ω h = T ∈T h T and its boundary Γ ∞,h = ∂Ω h . Notice that Ω h does not fit exactly the domain Ω and its subdomains Ω i since the latter are assumed to be at least of class C 2 . To this end, define the subtriangulation of T h related to Ω i by Then, let The following conditions state that T h fits approximately Ω and its subdomains 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 σ i on Ω i,h . To this end, assume that for any i = 1, . . . , L there is a domain Ω i such that Without loss of generality, we may assume that Ω i ∩ Ω i+2 = ∅ as well as Ω 2 ∩V m = ∅ for a given family of neighborhoods (V m ) m of the sources S m . Then, denote by σ i an extension of σ i on Ω i such that σ i ∈ W 1,∞ ( Ω i ) and σ i (x) ≥ σ min for almost every x ∈ Ω i . On T h , we introduce the standard vector space of Lagrange finite elements of type P1, where 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 ∩ L 2 0 (Ω h ) of P1 finite elements with zero mean value on Ω h .
Then the discrete problem reads: find w h ∈ V h such that Notice thatũ m is defined in a unique way on any extension Ω i outside the neighborhood V m of the source S m . Nevertheless, the functiong differs from the original data g = − M m=1 c m ∂ nũm on Γ ∞ since the normal vectors on Γ ∞ and Γ ∞,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 (·, ·) is clearly coercive on V h due to Poincaré-Wirtinger's inequality on Ω 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
Notice that the family (a h (·, ·)) h of bilinear forms defined by (5.4) is uniformly V hcoercive since the extensionsσ i are uniformly bounded from below by σ min and Ω h is included in the domain Ω def = L i=1 Ω i for any h. Therefore, Poincaré-Wirtinger's inequality holds true on Ω h with a constant independent from h. Theorem 8. Consider a regular family of meshes (T h ) h fitting the geometry of Ω and its subdomains in the sense of (5.1). Assume further that (T h ) h satisfies the following inverse assumption with a constant C inv > 0 independent from h, be an extension of w, the solution of problem (3.17), such thatw i|Ωi = w on Ω i for i = 1, . . . , L. Further, letσ = ( σ i ) i=1,...,L ∈ L i=1 W 1,∞ ( Ω i ) be an extension of σ such that for i = 1, . . . , L σ i|Ωi = σ i on Ω i and σ i (x) ≥ σ min almost everywhere.
Consider the family of discrete problems (5.3) and assume that the associated bilinear forms a h (·, ·) are uniformly V h -coercive. For any h, let w h ∈ V h be the solution of the discrete problem (5.3).
Then, the following error estimate holds true, The error estimate for the approximate solution u h follows immediately from (5.7). Indeed, u h is defined by u h = w h +ũ 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 σ on the different subdomains Ω i . Then, the interpolation error behaves as 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 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 σ 1 = σ 3 = 0.33S.m −1 for the brain and the scalp and σ 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 .
The second is the magnification factor (MAG) which is defined by The RDM and MAG are error functionals with respect to a reference solution u ref .
Obviously, an ideal model leads to RDM = 0 and 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 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 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 − dist(S, Γ 1,h )/r 1 for a dipole located at the point S. The results are reported in  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 Ω 2 , (5.10) The parameter σ skull > 0 represents the (constant) conductivity of the skull outside the fontanel, wheras σ f > 0 denotes the maximal conductivity inside the fontanel. Hereafter, we take σ skull = 0.04S.m −1 and σ 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 depending on the parameter α > 0. In the numerical experiments hereafter, α is set to 10 4 which amounts to saying that the fontanel is limited (up to 1.5%) to the region Ω f : Figure 5.2). Formula (5.10) leads to a regular function σ 2 and fits with the assumptions in Section 3.  Table 1). We define global errors on the whole domain Ω by 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 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 q 1 = (0, 0, J) and q 2 = (0, J, 0). The numerical results corroborate the convergence estimates of Section 5.1 that predict a theoretical convergence rate of τ = 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 σ 1 is constant on Ω 1 shows that the error in the H 1 -norm can be majored by δ −5/2 up to a constant independent from h and δ (cf. [40]). Here, δ denotes the distance of the sources from the interface brain/skull, δ = max 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 fullterm 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 (Ω 1 ), cerebrospinal fluid (Ω 2 ), skull (Ω 3 ), fontanels (Ω 3,f ⊂ Ω 3 ) and scalp (Ω 4 ). The T2weighted 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 × 0.47 × 0.47 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 × 0.35 × 0.63 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 (σ f , σ skull ) in the parameter set corresponds an EEG model, whereas the choice σ f = σ skull yields the model without fontanels. In order to quantify the influence of the values (σ f , σ 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 ref which are, respectively, the solutions to the fontanels and no-fontanel models.
The different conductivity values are fixed to σ 1 = σ 4 = 0.33S.m −1 for brain and scalp and σ 2 = 1.8S.m −1 for CSF. We compute the RDM and MAG factors for different couples (σ f , σ skull ) which are varying around a reference conductivity of σ f = 0.3S.m −1 and σ skull = 0.04S.m −1 . The uncertainty in head tissue conductivities is estimated at ±25% which leads to σ f ∈ [0.225, 0.375] and σ skull ∈ [0.03, 0.05]. A single dipole source is placed in Ω 1 at a distance of 5mm of the interface brain/CSF. Its moment is (0, 0, J) with intensity J = 2 √ 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 σ f /σ skull . The difference of electrical potential between models with or without fontanels is maximal (RDM=10% and MAG=15%) when the ratio σ f /σ 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 σ f /σ 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.  Table 2. Four-layer realistic head model 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.

Mesh Nodes Tetrahedra
The sensitivity u 1 of the potential in a given direction µ 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 Ω f , µ = 1 Ω f , as direction. Figure 6.4 compares the sensitivity of two sources of same moment q = (0, J, J) with J = 2 √ 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 15mm from the interface brain/CSF. We compare the sensitivity for the moments q 1 = (J, J, 0) and q 2 = (J, 0, J) with J = 2 √ 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.   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. estimate from classical inequalities in variational theory, 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 1 h by w h −w = hw 1 h . 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 This proves the strong convergence of the sequence (w 1 h ) h to w 1 in H 1 (Ω). It remains to show that D µ w belongs to L(L ∞ , V ). For fixed µ, D µ w(·, σ) is defined by the solution of (4.5). But the right-hand side of (4.5) is linear in µ. This is obvious for the first and the third term since w andũ m are independent from µ and p m = µ(S m ) is linear in µ. The second term depends on µ only via the Gâteaux derivativeũ 1 m ofũ m . According to Proposition 1, we haveũ 1 m = − µ(Sm) σ(Sm)ũ m which is a linear expression in µ.
In order to prove that the linear application µ → D µ w is continuous from P to V , we estimate w 1 with the help of formulation (4.5). Taking v = w 1 , we get using again thatũ 1 m = − µ(Sm) σ(Sm)ũ m . This yields the continuity of D µ w with respect to µ and proves that w(·, σ) is Gâteaux differentiable with respect to the conductivity σ.