M Ali Awan

Damien Masson
email: damien.masson@univ-paris-est.fr

Eduardo Tovar

Energy-aware Task Allocation onto Unrelated Heterogeneous Multicore Platform for Mixed Criticality Systems

Heterogeneous multicore platforms have become an attractive choice to deploy mixed criticality systems demanding diverse computational requirements. One of the major challenges is to efficiently harness the computational power of these multicore platforms while deploying mixed criticality applications. The problem is acerbated with an additional demand of energy efficiency. It is particularly relevant for the battery powered embedded systems. We propose a partitioning algorithm for unrelated heterogeneous multicore platforms to map mixed criticality applications that ensures the timeliness property and reduces the energy consumption.

I. INTRODUCTION

Modern real-time (RT) applications are becoming increasingly complex everyday. These applications may have different criticality levels. A criticality level corresponds to a level of assurance against failure. The highest-criticality functions are vital for the operations of the system and any violation of temporal constraint (fault) may lead to a disastrous consequences. However, fault of lower criticality functions may cause a tolerable decrease in the quality of service.

In the past, applications of different criticality were hosted on separate components. A recent trend in RT and embedded system domain to integrate the functionalities belonging to different criticality levels on a single hardware platform has paved a way towards mixed criticality systems. Mixed criticality systems are common in many domains such as avionics and automotive industry. For example, an unmanned aerial vehicle may have functionalities corresponding to flight critical (concerning the safety of flight) and mission critical (reconnaissance and surveillance) operations. It is vital for the integrity of the aerial vehicle to ensure the safe operation of a flight in all conditions while reconnaissance and surveillance are secondary functionalities.

In order to increase the level of assurance against failure for the high criticality functions, pessimistic estimates of the worst-case-execution time (WCET) are determined using static analysis. However, a mixed criticality system designed with the pessimistic estimates of WCET for the high criticality functions (tasks) will lead to inefficient utilisation of the resources and unnecessary over-provisioning of the system. This is driven by the fact that high criticality WCET estimates are rarely achieved in real execution. To efficiently exploit the available resources Steve Vestal [START_REF] Vestal | Preemptive scheduling of multi-criticality systems with varying degrees of execution time assurance[END_REF] proposed a mixed criticality system model. In this model, tasks have different estimates of WCET corresponding to different criticality levels. A system starts in a low criticality mode and the schedulability of the system is ensured assuming a WCET estimate corresponding to current low criticality mode. If any task exceeds its execution beyond the WCET estimate defined for the given criticality mode of operation, a system transitions to high criticality mode.

Modern multicore platforms provide sufficient computing capabilities to deploy mixed criticality systems on a single chip. Among different kinds of multicore architectures, heterogeneous multicore platforms composed of more than one heterogeneous processing units (cores) has gained popularity to perform specific tasks well and cheap. Various components (applications) in a mixed criticality system have diverse computing requirements that makes heterogeneous multicore platforms an attractive choice to opt. Despite all benefits, efficient mapping of mixed criticality applications on heterogeneous multicore platform is one of the major challenges faced by system designers.

The problem is acerbated with an additional limitation of energy supply. This is particularly relevant in battery powered and nomadic embedded systems. Even when the application is technically feasible upon the targeted platform in the sense that the platform can provide a sufficient computing capacity for the execution of the application, it has become unreasonable to expect to implement such a system without addressing the issue of reducing its energy consumption. The energy consumption of a heterogeneous platform can be reduced in two different ways. Firstly, a task should be mapped to a processor where it dissipates minimum dynamic power. Secondly, the energy-aware scheduling mechanism can be employed to further reduce the energy consumption on each core. This work explores the former technique.

We use partitioned scheduling to map the given set of tasks belonging to different criticality levels on a heterogeneous multicore platform such that available computational capacity of the hardware platform is efficiently exploited, and energy consumption is minimised while ensuring the timeliness properties of the system. This is an NP-hard problem in a strong sense as it is a special case of bin packing and hence, heuristics are the way forward.

There has been some efforts in the state-of-the-art [START_REF] Lakshmanan | Resource allocation in distributed mixed-criticality cyber-physical systems[END_REF]- [START_REF] Rodriguez | Multicriteria evaluation of partitioned edf-vd for mixed-criticality systems upon identical processors[END_REF] to propose new techniques and evaluate the potential of existing bin packing heuristic to perform the allocation with an objective to optimise the utilisation of the system. However, the state-of-the-art on the energy minimisation of mixed criticality systems is very limited. Broekaet et al. [START_REF] Florian Broekaert | Towards power-efficient mixed critical systems[END_REF] proposed system level power management approach that allocates energy budgets to virtual machines. In this techniques an overrunning virtual machine deducts its additional budget from budget of low priority virtual machine to be scheduled in future. Legout et al. [START_REF] Legout | Mixed-Criticality Multiprocessor Real-Time Systems: Energy Consumption vs Deadline Misses[END_REF] proposed a technique that provides trade-off between energy reduction and the number of deadline misses of low criticality tasks. Finally, Zhang et al. [START_REF] Zhang | Design optimization of security-sensitive mixed-criticality real-time embedded systems[END_REF] allocates energy budgets to different cores and perform allocation through genetic algorithm to extend the battery life time of the system. They assume that the energy consumption of a task is only a function of its execution time. We consider a realistic power model in which the energy consumption of a task depends on the characteristics of the core type.

II. SYSTEM MODEL

This work considers a partitioned unrelated heterogeneous multicore platform π def = [π 1 , π 2 , . . . , π M ] composed of M distant cores. We adopt a common mixed criticality task model, initially proposed by Steve Vestal [START_REF] Vestal | Preemptive scheduling of multi-criticality systems with varying degrees of execution time assurance[END_REF] and later extended by other researchers in the RT research community [START_REF] Burns | Mixed criticality systems-a review[END_REF]. We consider a dual criticality system in which a system executes in either low criticality mode (L-mode) or high criticality mode (H-mode) of operation. We assume a set of n independent sporadic tasks

τ def = {τ 1 , τ 2 , . . . , τ n }. Each task τ i def = T i , D i , κ i , → C i (H), → C i (L), → E i is characterised by its minimum inter-arrival time T i , deadline D i , criticality level κ i , a vector of WCET profile → C i (H) in H-mode of operation, a vector of WCET profile → C i (L) in L-mode of
operation and a vector of energy consumption profile

→ E i .
Each task τ i belongs to either high criticality level (H-task) or low criticality level (L-task), i.e., κ i ∈ {L, H}. The WCET profiles

→ C i (H) def = (C 1 i (H), C 2 i (H), . . . , C M i (H)) and → C i (L) def = (C 1 i (L), C 2 m (L), . . . , C M i (L)
) are the WCET estimates in the H-mode and L-mode of operation, respectively, on different cores indexed from 1 to M . Similarly, energy consumption profile

→ E i def = (E 1 i , E 2 m , . . . , E M i ) represents the energy consumption of task τ i on different cores in L-mode.
It is very common to assume in state-of-the-art that energy consumption of a task is a function of its execution time.

However, in reality, energy consumption on a certain processor depends also on the set of instructions it has to execute to perform the desired functionality. Various instructions use different parts of CPU, and hence may result in a different estimate of energy consumption. Therefore, two applications with identical execution time may consume different energy depending on the characteristics of the instructions used, and the number of cache misses involved. In this work, we consider a realistic power model in which the energy consumption of a task depends on the characteristics of the core type and hence, computed on each core using existing energy measurement technique [START_REF] Snowdon | Koala: A platform for os-level power management[END_REF]. This approach has the flexibility to incorporate the effect of other system resource usage on energy consumption such as memory and caches etc.

We assume that τ (L)

def = {τ i ∈ τ |κ i ∈ L} and τ (H) def = {τ i ∈ τ |κ i ∈
H} represent the subsets of lowand high-criticality tasks in τ , respectively. It is assumed that

∀τ i ∈ τ (H) ∧ ∀m ∈ [1, 2, . . . , M ], C m i (L) ≤ C m i (H) and ∀τ i ∈ τ (L) ∧ ∀m ∈ [1, 2, . . . , M ], C m i (H) = 0.
Tasks are not allowed to migrate after assignment as we assume partitioned scheduling. The schedulability analysis proposed by Ekberg and Yi [START_REF] Ekberg | Bounding and shaping the demand of generalized mixed-criticality sporadic task systems[END_REF] is used to test the feasibility of the system in both low and high criticality mode. This analysis shortens the deadlines of τ (H) in the L-mode to keep the schedule ahead of time which allows a safe transition from L to H-mode. The utilisation of task τ i in low and high criticality mode on processor π m are represented as

U m i (L) def = C m i (L) T i and U m i (H) def = C m i (H) T i , respectively.
A subset of task assigned to a core π m is represented as τ (π m ). The overall utilisation of tasks assigned to core π m in L and H-mode is defined as U m (L)

def = ∀τi∈τ (π m ) C m i (L) and U m (H) def = ∀τi∈τ (π m )
C m i (H), respectively. We assume an arbitrary deadline model in our system in which task deadline has not relation with its minimum inter-arrival time.

III. PRELIMINARIES

Initially, we define the following concepts needed to explain the proposed allocation heuristics.

Definition 1 (Energy density ED m i ): The energy density of a task τ i on a core π m is defined as

ED m i def = E m i T i
, where E m i corresponds to the energy consumption of task in the low criticality mode. This is a measure of the average power dissipation in the low criticality mode.

Definition 2 (Energy density difference EDD m i ): The energy density difference of a task τ i on a core π m is defined as

EDD m i def = min{ED h i : h = m ∧ ED h i ≥ ED m i } -ED m i .
This value corresponds to the difference between the next higher energy density of τ i on any other core and the energy density on the current core.

Definition 3 (Low criticality utilisation difference LU D m i ): The low criticality utilisation difference of a task τ i on a core π m is defined as

LU D m i def = min{U h i (L) : h = m ∧ U h i (L) ≥ U m i (L)} -U m i (L).
This value is the difference between the next higher low criticality utilisation of τ i on any other core and the low criticality utilisation on the current core.

Definition 4 (High criticality utilisation difference HU D m i ): The high criticality utilisation difference of a task τ i ∈ τ (H) on a core π m is defined as

HU D m i def = min{U h i (H) : h = m ∧ U h i (H) ≥ U m i (H)} -U m i (H).
This measure represents the difference of next higher high criticality utilisation of τ i ∈ τ (H) on any other core and the high criticality utilisation on the current core.

Initially, we compute ED m i and LU D m i for all tasks, and HU D m i for high criticality tasks (i.e., ∀τ i ∈ τ (H)) on each core. Afterwards, we generate the following three lists out of these values. i }. S EDD is a list that contains all the elements of EDD sorted in non-increasing order.

IV. PROPOSED APPROACH

In our proposed approach, we compute a set of feasible allocations and select the one that minimises the overall energy consumption of the system. Each feasible allocation ensures the schedulability of assigned tasks on each core both in LOW and HIGH criticality modes. The following method is used to collect the feasible allocations. A suffrage based task-to-core allocation scheme ILLED (proposed by Awan et al. [?] for single criticality systems) ranks the tasks based on a metric called density difference1 while performing allocation to their preferred core. The density difference of a task shows how much a system will lose in term of given criteria (utilisation or energy) if a task is not allocated to its preferred core. The ILLED algorithm ranks tasks based on density difference metric and performs the allocation to optimise the given criteria (utilisation or energy). This ranking plays an important role in the allocation process. The tasks ranked higher have the higher probability of getting mapped to their preferred core. The proposed allocation heuristics manipulates S LU D , S HU D and S EDD density difference lists to obtain a set of feasible task-to-core mappings. The ILLED algorithm signals a failure if the mapping is not schedulable2 in either high criticality or low criticality mode.

The pseudocode of the proposed approach is given in Algorithm 1. First of all, we perform the task-to-core mapping based on a list S EDD through ILLED algorithm. If the taskto-core mapping is successful, meaning it ensures the schedulability in both high criticality and low criticality modes, we are done. As the ordering of the tasks in S EDD is based on the energy density difference, therefore, ILLED generates a mapping to minimise the energy consumption. However, if a system is not schedulable in either of high or low criticality mode, we need to determine an order of a list that leads to a feasible task-to-core mapping ensuring schedulability in both modes.

In mixed criticality systems, it is very important to ensure the schedulability of the high criticality tasks in both low and high criticality modes. A conservative allocation approach gives preference to the high criticality tasks in the assignment process to ensures their schedulability in both modes. Such an allocation can be achieved through S HU D list. Assume, we divide the list S HU D into two parts S H and S L corresponding to high and low criticality tasks respectively, and apply the ILLED algorithm on S H followed by S L .

It is also desirable to ensure the schedulability of the low criticality tasks in a low criticality mode. In practice, mixed criticality system mostly stays in low criticality mode and occasionally transitions into a high criticality mode triggered due to an overrun or other external/internal events. Therefore, it is more beneficial to optimise the energy efficiency in the low criticality mode. An allocation based on a list S LU D ordering maximises the schedulability in the low criticality mode as it is sorted with respect to low criticality utilisation. However, there may be a possibility that system may not be feasible in high criticality mode when an allocation is performed with such an ordering (S LU D ). Therefore, we start with a list S LU D in our allocation process and gradually promote the high criticality tasks in our S LU D list corresponding to the ordering given in S HU D to generate all feasible allocations.

The reordering mechanism is explained as follows. We perform the allocation with respect to the given list and compute its energy consumption if the allocation is successful (schedulable both in low and high mode). Afterwards, we modify S LU D list, by promoting the entry of high criticality task by one position to get a new list. As modification in S LU D list is performed with respect to the list S HU D , therefore, we select the high criticality task corresponding to the entry on

(E[i]) 23: return A[index]
top of the S HU D list and promote its corresponding entry in S LU D by one position. If the task corresponding to the top most element in S HU D is also at the same level as in S LU D , we need to select the task corresponding to the second element of S HU D in S LU D and promote it by one position in S LU D unless it reaches the same position as in S HU D . As mentioned previously, after every modification we have to compute the energy consumption, if the allocation is successful. After all these promotions of high criticality tasks, S LU D will have the same ordering as given in S HU D . However, S LU D does not represent the same list as S HU D as high criticality tasks are ordered with respect to HU D m i in S HU D . Therefore, as a final step, we perform the allocation based on a list S HU D and compute the energy consumption on successful allocation. The allocation with a list S HU D is performed in a slightly different manner. We first divide the list S HU D into two lists S HT and S LT representing high and low criticality tasks. All elements in S HT and S LT are sorted in non-increasing order with respect to HU D m i and LU D m i respectively. First, we perform the allocation of S HT and then on top of this we perform the allocation of low criticality task from a list S LT . The energy consumption is computed on successful allocation and stored in a variable E HU D . Let E LU D be the minimum of the energy consumptions computed for all the iteration performed on a list S LU D . A minimum of E HU D and E LU D will give us our desired allocation that ensures the schedulability and reduces the total energy consumption.

Finally, we explain an adapted ILLED algorithm [?]. This algorithm assigns the tasks to their favourite core in the specified order of the provided list. If a task cannot be assigned to its preferred core, the current entry of the task in the provided list is updated with the value on next preferred core in the order. The entry corresponding to this task is added again in the provided list on its appropriate location based on its new updated value. On a successful assignment the task from the provided list is removed. This allocation scheme should continue unless the provided list is empty or a task cannot be assigned to any core. We used Ekberg and Yi's [START_REF] Ekberg | Bounding and shaping the demand of generalized mixed-criticality sporadic task systems[END_REF] test to check the schedulability of a core after mapping any task.

V. CONCLUSION AND FUTURE DIRECTIONS

This work presents the energy-aware task-to-core allocation for mixed criticality systems deployed on heterogeneous multicore platform considering a realistic power model. Initial assessment of this algorithm predicted upto 8.9% gains over the traditional first-fit bin packing heuristics in terms of energy consumption and scheduled task-set. In the future, it is intended to evaluate the performance of the proposed approach for different types of applications and hardware platforms. Another possible extension is to integrate the effect of leakage power dissipation in the proposed approach.

1 )= LU D min 1 , LU D min 2 , LU D min 3 ,

 1123 Assume a set LU D such that LU D def i }. S LU D is a list that contains a sequence of numbers of LU D sorted in a non-increasing order. 2) Assume a set HT def = HU D min j1 , HU D min j2 , . . . , HU D min j |τ (H)| , where HU D min jx def = min ∀m {HU D m jx } and τ jx ∈ τ (H). Similarly, assume another set LT def = LU D min k1 , LU D min k2 , . . . , LU D min k |τ (H)|, where τ kx ∈ τ (L). Let S HT and S LT represent the sequence of elements of set HT and LT , respectively, sorted in non-increasing order, then, S HU D is a list that concatenates S HT and S LT , i.e., S HU D def = S HT ||S LT , where || represents the concatenation sign.3) Assume that EDD is a set defined as EDD def

  Algorithm 1 Pseudocode of the proposed algorithm Input: τ, π Output: Assignment 1: Compute S EDD , S LU D and S HU D 2: Assume A[i] represent the i th assignment 3: A[0] := Perform the ILLED allocation with S EDD 4: if (A[0] is feasible) then return A[0] Perform the ILLED allocation with S HT followed by S LT , 19: if (A[i] is feasible) then

	5: end if
	6: i := 0
	7: while (true) do
	8:	A[i] := Perform the ILLED allocation with S LU D
	9:	if (A[i] is feasible) then
	10:	E[i] := Compute energy consumption of A[i]
	11:	i := i + 1
	12:	end if
	13:	Promote high criticality task's entry in S LU D
	14:	if (Promotion fails) then
	15:	break
	16:	end if
	17: end while
	18: A[i] := 20: E[i] := Compute energy consumption of A[i]
	21: end if
	22: index := Find index that gives min ∀i

The density difference is a generic term used to indicate any of the following quantities, i) low criticality utilisation difference, ii) high criticality utilsation difference or iii) energy density difference.

The ILLED algorithm is initially designed for single criticality system. To make it compatible with mixed criticality systems, we replaced the feasibility test of the ILLED algorithm with Ekberg and Yi's schedulability analysis[START_REF] Ekberg | Bounding and shaping the demand of generalized mixed-criticality sporadic task systems[END_REF].

This work was partially supported by National Funds through FCT/MEC (Portuguese Foundation for Science and Technology) and co-financed by ERDF (European Regional Development Fund) under the PT2020 Partnership, within project UID/CEC/04234/2013 (CISTER Research Centre); also by FCT/MEC and the EU ARTEMIS JU within project(s) ARTEMIS/0001/2013 -JU grant nr.621429 (EMC2) and ARTEMIS/0003/2012 -JU grant nr.333053 (CONCERTO).


Malal M Darbas

Mohamedou M Diallo

Abdellatif El Badia

S Lohrengel

SENSITIVITY ANALYSIS FOR THE EEG MODEL IN NEONATES WITH RESPECT TO VARIATIONS OF THE CONDUCTIVITY

whether they are published or not. The documents may come

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 regions on the scalp. The electrical activity recorded is the synchronous activity of a large number of neighboring neurons in the cerebral cortex beneath the skull. The measurements provide valuable information about the sources that are at the origin of pathological activities of the brain. In particular, EEG is one of the main diagnostic tests in presurgical evaluation for refractory epilepsy.

The accuracy of the EEG source reconstruction relies heavily on the accuracy of the associated forward model. EEG source reconstruction is an inverse problem that aims to identify the sources responsible of electrical brain activity from the knowledge of the measured potentials at the electrodes on the scalp. The EEG forward problem 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 [START_REF] De Munck | A fast method to compute the potential in the multisphere model[END_REF][START_REF] Zhang | A fast method to compute surface potentials generated by dipoles within multilayer anisotropic spheres[END_REF]. On the other hand, realistic head models obtained from segmentation of magnetic resonance imaging (MRI) are able to take into account the precise geometry of the different tissues. Several source models have been developed as e.g. partial integration, the St. Venant model or the subtraction approach [START_REF] Wolters | Numerical mathematics of the subtraction approach for the modeling of a current dipole in EEG source reconstruction using finite element head models[END_REF][START_REF] Bauer | Comparison Study for Whitney (Raviart-Thomas) -Type Source Models in Finite-Element-Method-Based EEG Forward Modeling[END_REF]. For the numerical resolution of the forward problem, both boundary elements and 3D finite elements are commonly used [START_REF] Hallez | Review on solving the forward problem in EEG source analysis[END_REF][START_REF] Grech | Review on solving the inverse problem in EEG source analysis[END_REF].

Inaccuracies in the EEG forward problem impact the precision of source localization in the inverse problem. In the head model for adults, the effect of skull anisotropy and inhomogeneity as well as the uncertainty in the tissue conductivities have been investigated [START_REF] Vorwerk | A guideline for head volume conductor modeling in EEG and MEG[END_REF]. From a mathematical point of view, different cost functions have been analyzed [START_REF] Badia | A stable recovering of dipole sources from partial boundary measurements[END_REF][START_REF] Farah | Problèmes inverses de sources et lien avec l'Electro-encéphalo-graphie[END_REF].

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. Two characteristics are inherent to the neonate. The first one is that neonates show higher skull conductivities than children or adults [START_REF] Gargiulo | The effect of fontanel on scalp EEG potentials in the neonate[END_REF]. 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 bone. The model for adults considering an homogeneous skull conductivity is a priori not appropriate in the context of neonates. It may be noticed also 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.

The present work aims at proposing an EEG model for neonates and its numerical validation for both the multi-layer spherical head model and a realistic neonate head model. Moreover, with regard to EEG source reconstruction, it is important to understand the influence of fontanel and skull conductivities on the forward problem.

Commonly, this is done via the computation of two error functionals, respectively the RDM (Relative Difference Measure) and MAG (MAGnification factor), that allow to compare the difference between two models. In this paper, we introduce an additional analysis tool to investigate the sensitivity of the potential solution of the forward 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. From the best of our knowledge, this is the first time that a mathematical sensitivity analysis has been performed in the context of EEG models.

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. The 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 conductivity.

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 [START_REF] Faugeras | The inverse EEG and MEG problems : The adjoint state approach I : The continuous case[END_REF]. In terms of the electric field E and the magnetic field H, this yields

∇ • (εE) = ρ, (2.1a) curl E = 0, (2.1b) curl H = J, (2.1c) ∇ • (µH) = 0. (2.1d)
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, (2.2)

J = σE + J s ,
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.

(2.

3)

E = ∇u.
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 Ω

(2.4) -∇ • (σ∇u) = ∇ • J s .
In the typical multi-layer head model, we distinguish three to five layers for the brain (containing gray and white matters, CSF), skull, and scalp. Therefore, consider a partition of Ω into L open nested subdomains (Ω i ) i=1,...,L (see Figure 2.1) such that

Ω = L i=1
Ω i and

Ω i ∩ Ω j = ∅ ∀i = j.
For i = 1, . . . , L -1, we denote by Γ i the interface between Ω i and Ω i+1 . We further denote by Γ ∞ = ∂Ω the exterior boundary of the whole domain Ω. Let n i be the unit normal vector on Γ i oriented towards the exterior of Ω i . We assume that (Γ i ) i are closed regular surfaces such that

Γ i ∩ Γ j = ∅ ∀i, j ∈ {1, . . . , L -1} ∪ {∞}.
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. In order to take into account the presence of the fontanels in neonates, we will consider in the sequel conductivities σ i that are functions of the position, i.e. σ i = σ i (x) in Ω i . We further assume that σ i has H 1 -regularity, (2.5) The source model of neural activity can be described by a sum of M electric dipoles located in the brain (e.g. [START_REF] Schneider | A multistage process for computing virtual dipole sources of EEG discharges from surface information[END_REF]). Each dipole is characterized by its position S m ⊂ Ω 1 and its moment q m which is a vector of R 3 . The current density J s thus reads

σ i ∈ H 1 (Ω i ) ∀i = 1, . . . , L.
J s = M m=1 q m δ Sm
where δ Sm denotes the delta distribution at S m . The right hand side of (2.4) is then given by (2.6)

F def = ∇ • J s = M m=1 q m • ∇δ Sm .
In the sequel, we assume that the conductivity is constant in a given neighborhood of each source.

More precisely we fix a family of open balls

(V m ) m=1,...,M such that V m ⊂⊂ Ω 1 , S m ∈ V m and
(2.7)

σ 1|Vm ≡ c m ∈ R
for any m ∈ {1, . . . , M }. Without loss of generality, we assume that the balls are non intersecting,

V m ∩ V p = ∅ if m = p.
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

(2.8) -∇ • (σ∇u) = F, in Ω, σ∂ n u = 0, on Γ ∞ ,
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

[u] |Γi = 0 on Γ i (i = 1, . . . , L -1), (2.9a) [σ∂ n u] |Γi = 0 on Γ i (i = 1, . . . , L -1). (2.9b)
Here and below, [f ] |Γi = f |Ωi -f |Ωi+1 denotes the jump of the quantity f across the interface Γ i . For given sources (S m , q m ) m and a known distribution of the conductivity σ, problem (2.8) is the forward EEG problem.

Existence and uniqueness result

The source term F given by (2.6) belongs to H s (R 3 ) for any s < -5/2 which prohibits a variational formulation of (2.8) in H 1 (Ω). To overcome this problem, we adapt an idea from [START_REF] Farah | Problèmes inverses de sources et lien avec l'Electro-encéphalo-graphie[END_REF]. This method was also introduced in [START_REF] Wolters | Numerical mathematics of the subtraction approach for the modeling of a current dipole in EEG source reconstruction using finite element head models[END_REF] as the subtraction approach for a single source.

3.1.

Lifting of the singularity. The idea of the subtraction method is to decompose the potential u, solution to (2.8), into a potential ũ which contains the singularity and a regular lifting w

(3.1) u = ũ + w, with ũ = M m=1 ũm .
The singular potential ũm is the solution of the Poisson equation in an unbounded homogeneous conductor of conductivity c m = σ 1 (S m ),

(3.2) -c m ∆ũ m = q m • ∇δ Sm in R 3 .
It is actually obtained by convolution of the fundamental solution of the Laplace equation with the right hand side 1 c m q m • ∇δ Sm and reads

(3.3) ũm (x) = 1 4πc m q m • (x -S m ) |x -S m | 3 , ∀x ∈ R 3 \ {S m }.
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 ∇

• (σ 1 ∇ũ m ) is well defined on Ω 1 by ∇ • (σ 1 ∇ũ m ) = ∇ • ((σ 1 -c m )∇ũ m ) + c m ∆ũ m .
Indeed, both terms on the right hand side of the above identity are well defined as distributions on Ω 1 , even if σ 1 has only H 1 -regularity, since σ 1 -c m vanishes identically on V m and ũm is regular outside ∪ m V m . Therefore, we get

-∇ • (σ 1 ∇w) = -∇ • (σ 1 ∇(u -ũ)) = F + M m=1 ∇ • ((σ 1 -c m )∇ũ m ) + c m ∆ũ m on Ω 1 .

It follows from the definition of ũm

that F = - M m=1 c m ∆ũ m . Next, taking into account the identity (3.4) ∇ • ((σ 1 -c m )∇ũ m ) = (σ 1 -c m )∆ũ m + ∇(σ 1 -c m ) • ∇ũ m
and noticing that the first term on the right hand side vanishes on Ω 1 , we get

-∇ • (σ 1 ∇w) = M m=1 ∇(σ 1 -c m ) • ∇ũ m = ∇σ 1 • ∇ũ.
On Ω i , i = 1, the potentials ũm are regular and we see from a direct computation that w is solution of the following problem

(3.5) -∇ • (σ i ∇w) = ∇σ i • ∇ũ, in Ω i (i = 1, . . . , L), σ∂ n w = -σ∂ n ũ, on Γ ∞ ,
with transmission conditions

[w] |Γi = 0 on Γ i , (3.6a) [σ∂ n w] |Γi = (σ i+1 -σ i )∂ n ũ on Γ i , (3.6b)
for i = 1, . . . , L -1. Recall that the right hand side of (3.5) vanishes in any neighborhood V m of the sources, since σ 1 is constant on V m , and that it is smooth in

Ω \ ∪ M m=1 V m . 3
σ i ∇w •∇v dx = L-1 i=1 ˆΓi (σ i+1 -σ i )∂ n ũv ds- ˆΓ∞ σ L ∂ n ũv ds+ L i=1 ˆΩi (∇σ i •∇ũ)v dx.
Since ∇(σ i -c m ) = ∇σ i , we can show with the help of (3.4) that (3.7) is equivalent to

(3.8) ˆΩ σ∇w • ∇v dx = M m=1 ˆΩ(c m -σ)∇ũ m • ∇v dx - ˆΓ∞ c m ∂ n ũm v ds
which is the variational formulation of the boundary value problem (3.9)

     -∇ • (σ∇w) = M m=1 ∇ • ((σ -c m )∇ũ m ) , in Ω, σ∂ n w = -σ∂ n ũ, on Γ ∞ .
In the following, we focus on formulation (3.8). We introduce the bilinear form a(•, •) defined on H 1 (Ω) × H 1 (Ω) by Proof. Let ṽ be the solution of -∆ṽ = q m • ∇δ Sm obtained by convolution of the source term with the fundamental solution G of the Laplacian in R 3 . We then get from differentiation rules for the convolution product

ṽ(x) = (G * (q m • ∇δ Sm )) (x) = q m • ∇(G * δ Sm )(x) = q m • ∇G(x -S m )
where ∇ denotes the gradient with respect to the variable x. The symmetry of the Green's function G yields further

q m • ∇G(x -S m ) = -q m • ∇ Sm G(x -S m ).
Therefore, we have

∂ n ṽ(x) = -q m • ∇ Sm (∂ n G(x -S m )) ,
for any x = S m . This implies

ˆΓ∞ ∂ n ṽ ds = -q m • ∇ Sm ˆΓ∞ ∂ n G(• -S m ) ds,
since S m ∈ Ω 1 and thus S m ∩ Γ ∞ = ∅. Now, notice that the formula of the solid angle (see e.g. [START_REF] Kress | Linear Integral Equations[END_REF]) states that

(3.14) ˆΓ∞ ∂ n G(x -y) ds(x) =    -1 if y ∈ Ω, -1 2 if y ∈ Γ ∞ , 0 otherwise.
Hence, the quantity ´Γ∞ ∂ n G(• -S m ) ds is constant equal to -1 and its gradient with respect to S m vanishes which implies ˆΓ∞ ∂ n ṽ ds = 0. Now, (3.13) follows since ũm = 1 c m ṽ.

A solution to (3.8) 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,

(3.15) V = v ∈ H 1 (Ω) ˆΩ v dx = 0 .
On V , the Poincaré-Wirtinger inequality holds true,

(3.16) v 0,Ω ≤ C P ∇u L 2 (Ω) ∀v ∈ V.
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.10) and (3.11), respectively. Then the variational problem

(3.17) Find w ∈ V such that a(w, v) = l(v), ∀v ∈ H 1 (Ω)
has exactly one solution w ∈ V . Moreover, the following estimate holds true,

(3.18) w H 1 (Ω) M m=1 ∇ũ m L 2 (Ω\Vm) + ∂ n ũm L 2 (Γ∞) .
Proof. It follows from standard arguments in variational theory that the bilinear form a(•, •) is continous on H 1 (Ω) × H 1 (Ω) and V -elliptic. Since c m -σ vanishes in V m , we have

(3.19) |l(v)| ≤ M m=1 ∇ũ m L 2 (Ω\Vm) + C T ∂ n ũm L 2 (Γ∞) v H 1 (Ω) ,
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

a(w, v) = l(v) ∀v ∈ V.
Next, let v belong to H 1 (Ω). We have v -v Ω ∈ V where v Ω = 1 |Ω| ´Ω v dx is the mean value of v. Since l(1) = 0 according to the compatibility condition (3.12) and ∇v = ∇(v -v Ω ), we get

l(v) = l(v -v Ω ) = a(w, v -v Ω ) = a(u, 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.

The following theorem states the global H 2 -regularity of the variational solution of (3.17) in the subdomains Ω i .

Theorem 2. In addition to the assumptions of Theorem 1, assume that σ i ∈ W 1,∞ (Ω i ) for any i = {1, . . . , L} and that Γ i is of class C 2 for i ∈ {1, . . . , L -1} ∪ {∞}. Let w ∈ H 1 (Ω) be the solution of the variational problem (3.17). Then we have w |Ωi ∈ H 2 (Ω i ) for any i = {1, . . . , L} and

(3.20) w H 2 (Ωi) M m=1 ∇ũ m H 1 (Ωi\Vm) + ∂ n ũm H 1 (Γ∞) .
The proof of Theorem 2 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

-∇ • (σ i ∇w) = f i with f i def = M m=1 ∇ • ((σ i -c m )∇ũ m = M m=1 ∇(σ i -c m ) • ∇ũ m .
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 . 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. Nevertheless, we give the full proof of global H 2 -regularity on the subdomains in Appendix B since most textbooks deal in general only with Dirichlet data on a single subdomain.

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. Mathematically, a rigorous way to describe sensitivity is given by Gâteaux differentiability which expresses a weak concept of derivative. Definition 3. Let F : X → Y be an application between two Banach spaces X and Y . Let U ⊂ X be an open set. The directional derivative D µ F (σ) of F at σ ∈ U in the direction µ ∈ X is defined as

D µ F (σ) = lim h→0 F (σ + hµ) -F (σ) h if the limit exists. If D µ F (σ) exists for any direction µ ∈ X and if the application µ → D µ F (σ) is linear continuous from X to Y , F is called Gâteaux differentiable at σ. Now, let (V m ) m be a fixed family of neighborhoods of the sources such that for any m, S m ∈ V m ⊂ Ω 1 and V m ∩ V m = ∅ if m = m .
We introduce the parameter space

P = σ ∈ L ∞ (Ω) σ |Vm ≡ const. ∀m = 1, . . . , M
as well as the (open) subset

P adm = {σ ∈ P | σ min < σ < σ max }
of admissible conductivities. Here, σ min and σ max are two fixed positive constants. According to Theorem 1, problem (3.8) 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

(4.1) ũm (x, σ) = 1 4πc m q m • (x -S m ) |x -S m | 3
is the solution of the Poisson equation

-c m ∆ũ m (•, σ) = q m • ∇δ Sm in R 3
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 have

(4.2) ũm (x, σ + hµ) = 1 4π(c m + hp m ) q m • (x -S m ) |x -S m | 3 ,
where p m = µ(S m ), and ũm (•, σ + hµ) is solution of the following perturbed Poisson equation,

-(c m + hp m )∆ũ m (•, σ + hµ) = q m • ∇δ Sm in R 3 .
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 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

(4.3) D µ ũm (•, σ) = - p m c m ũm (•, σ) with p m = µ(S m ) and c m = σ(S m ).
Proof. A straightforward computation of the differential quotient yields

(4.4) ũm (x, σ + hµ) -ũm (x, σ) h = - p m c m + hp m ũm (x, σ), ∀x = S m
and (4.3) follows. The right-hand side of (4.3) is obviously linear and continuous in µ since p m = µ(S m ).

Theorem 4. Let σ ∈ P adm such that σ + hµ ∈ P adm for any h ∈ [-h 0 , h 0 ] and any µ ∈ P with µ ∞ = 1. Then the solution w(•, σ) of (3.8) 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

w 1 ∈ V such that ˆΩ σ∇w 1 • ∇v dx = -ˆΩ µ∇w • ∇v dx (4.5) + M m=1 ˆΩ(c m -σ)∇ũ 1 m • ∇v dx - ˆΓ∞ c m ∂ n ũ1 m v dx + M m=1 ˆΩ(p m -µ)∇ũ m • ∇v dx - ˆΓ∞ p m ∂ n ũm v dx ,
for all v ∈ V . Here, we note ũ1 m = D µ ũm (•, σ) and ũm = ũm (•, σ). Proof. 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 (Ω) and satisfies the compatibility condition l(1) = 0. Indeed, both boundary integrals in (4.5) vanish for v ≡ 1 due to the assertion of Lemma 1, taking into account that ũ1 m = -pm cm ũm . In order to investigate the Gâteaux derivative of w with respect to σ, let µ ∈ P with µ ∞ = 1. From the assumptions, we have σ + hµ ∈ P adm for any h ∈ [-h 0 , h 0 ].

In the sequel we shall omit the dependence of w and ũm on the parameters σ and µ for better reading and set

ˆΩ(σ + hµ)∇w h • ∇v dx = M m=1 ( ˆΩ ((c m + hp m ) -(σ + hµ)) ∇ũ m,h • ∇v dx - ˆΓ∞ (c m + hp m )∂ n ũm,h v ds), and
(4.7) ˆΩ σ∇w • ∇v dx = M m=1 ˆΩ(c m -σ)∇ũ m • ∇v dx - ˆΓ∞ c m ∂ n ũm v ds , for all v ∈ V .
The assumption on σ and µ implies that σ + hµ > σ min for any h ∈ [-h 0 , h 0 ] and therefore the bilinear form in (4.6) is V -elliptic with a constant independent from h. The righthand side of (4.6) defines a continuous linear form on H 1 (Ω) and we get the following estimate from similar arguments as for (3.18) in Theorem 1:

(4.8) w h H 1 (Ω) M m=1 ∇ũ m,h L 2 (Ω\Vm) + ∂ n ũm,h L 2 (Γ∞) .
In order to identify the Gâteaux derivative of w, we introduce the differential quotients

w 1 h = w h -w h and ũ1 m,h = ũm,h -ũm h .
Subtracting (4.7) from (4.6) and dividing by h leads to

ˆΩ σ∇w 1 h • ∇v dx = -ˆΩ µ∇w h • ∇v dx (4.9) + M m=1 ˆΩ(c m -σ)∇ũ 1 m,h • ∇v dx - ˆΓ∞ c m ∂ n ũ1 m,h v ds + M m=1 ˆΩ(p m -µ)∇ũ m,h • ∇v dx - ˆΓ∞ p m ∂ n ũm,h v ds .
We compare the above formulation for the differential quotient w 1 h with the variational formulation (4.5) for w 1 :

ˆΩ σ∇ w 1 h -w 1 • ∇v dx = -ˆΩ µ∇ (w h -w) • ∇v dx (4.10) + M m=1 ˆΩ(c m -σ)∇ ũ1 m,h -ũ1 m • ∇v dx - ˆΓ∞ c m ∂ n ũ1 m,h -ũ1 m v ds + M m=1 ˆΩ(p m -µ)∇ (ũ m,h -ũm ) • ∇v dx - ˆΓ∞ p m ∂ n (ũ m,h -ũm ) v ds .
Notice that the integrals over Ω in the second and third term vanish on V m since σ |Vm ≡ c m and µ |Vm ≡ p m . Taking v = w 1 h -w 1 in (4.10), we get the following estimate from classical inequalities in variational theory,

w 1 h -w 1 H 1 (Ω) ∇(w h -w) L 2 (Ω) (4.11) + M m=1 ∇ ũ1 m,h -ũ1 m L 2 (Ω\Vm) + ∂ n ũ1 m,h -ũ1 m L 2 (Γ∞) + M m=1 ∇ (ũ m,h -ũm ) L 2 (Ω\Vm) + ∂ n (ũ m,h -ũm ) L 2 (Γ∞)
The second and third term in the right-hand side of (4.11) 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 (4.9) and we get

h -1 w h -w H 1 (Ω) ∇w h L 2 (Ω) + M m=1 ∇ũ 1 m,h L 2 (Ω\Vm) + ∂ n ũ1 m,h L 2 (Γ∞) (4.12) + M m=1 ∇ũ m,h L 2 (Ω\Vm) + ∂ n ũm,h L 2 (Γ∞) .
The right hand side in (4.12) is obviously bounded when h tends to zero. Indeed, we have lim h→0 ũ1 m,h = ũ1 m and lim h→0 ũm,h = ũm , and ∇w h is bounded in terms of ũm,h according to (4.8). Multiplying (4.12) by h shows that w h -w H 1 (Ω) is of order h. Consequently, (4.11) becomes

w 1 h -w 1 H 1 (Ω) h M m=1 ∇ũ m L 2 (Ω\Vm) + ∂ n ũm L 2 (Γ∞) .
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

w 1 H 1 (Ω) µ ∞ ∇w L 2 (Ω) + M m=1 |µ(S m )| ∇ũ m L 2 (Ω\Vm) + ∂ n ũm L 2 (Γ∞)
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 σ.

Lemma 2. Let ũm = ũm (•, σ) and ũm,h = ũm (•, σ + hµ) be given by (4.1) and (4.2), respectively. Under the assumptions of Theorem 4, the following estimates hold true:

∇ (ũ m,h -ũm ) L 2 (Ω\Vm) h||∇ũ m || L 2 (Ω\Vm) , (4.13a) ∂ n (ũ m,h -ũm ) L 2 (Γ∞) h||∂ n ũm || L 2 (Γ∞) . (4.13b) Further, let ũ1 m,h = ũm,h -ũm h
and denote by ũ1 m = D µ ũm (•, σ) the Gâteaux-derivative of ũm at σ in the direction µ. Under the assumptions of Theorem 4, the following estimates hold true:

∇ ũ1 m,h -ũ1 m L 2 (Ω\Vm) h||∇ũ m || L 2 (Ω\Vm) , (4.14a) ∂ n ũ1 m,h -ũ1 m L 2 (Γ∞) h||∇ũ m || L 2 (Γ∞) . (4.14b)
Proof. From the definition of ũm and ũm,h we get

(4.15) ũm,h (x) -ũm (x) = 1 4π 1 c m + hp m - 1 c m q m • (x -S m ) |x -S m | 3 = - hp m c m + hp m ũm (x).
Integration of the gradient of the above expression over Ω\V m yields (

ũ1 m,h -ũ1 m = - p m c m + hp m + p m c m ũm = h p 2 m c m (c m + hp m ) ũm .
The boundedness of |p m | and c m (resp. c m + hp m ) yields (4.14a) by integration of the gradient of (4.16) over Ω \ V m . (4.14b) follows in the same way.

Remark 5. Assume that the variation of the conductivity in a given direction µ occurs only in the skull Ω 2 . In this particular case, the singular potential ũ defined by (3.3) is independent from µ and we have

w(•, σ + hµ) -w(•, σ) = u(•, σ + hµ) -ũ(•, σ 1 ) -(u(•, σ) + ũ(•, σ 1 )) = u(•, σ + hµ) -u(•, σ).
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).

Finite Element formulation of the EEG problem

In this section, we state the discretization of problem (3.17) by standard finite elements of type P1.

5.1. Discretization and convergence analysis. Throughout this subsection, we assume that the regularity assumptions of Theorem 2 are fulfilled. Consider a family {T h } h of tetrahedral meshes satisfying the usual regularity assumptions (see e.g. [START_REF] Ciarlet | The Finite Element Method for Elliptic Problems[END_REF])). 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

T h,i = T ∈ T h N h ∩ T ⊂ Ω i ∀i = 1, . . . , L. Then, let Ω i,h = T ∈T h,i T and Γ i,h = N h ∩ Γ i .
The following conditions state that T h fits approximately Ω and its subdomains

L i=1 Ω i,h = Ω h (5.1a) Γ i,h = Ω i,h ∩ Ω i+1,h ∀i = 1, . . . , L -1. (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 σ i on Ω i,h . To this end, assume that for any i = 1, . . . , L there is a domain Ω i such that

Ω i ⊂ Ω i and Ω i,h ⊂ Ω i .
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, (5.2)

X h = v h ∈ C 0 (Ω h ) v h|T ∈ P 1 (T ) ∀T ∈ T h
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

(5.3) a h (w h , v h ) = l h (v h ) ∀v h ∈ V h with (5.4) a h (w h , v h ) = L i=1 ˆΩi,h σ i ∇w h • ∇v h dx and l h (v h ) = L i=1 ˆΩi,h F i • ∇v h + ˆΓ∞,h gv h ds
where

F i = M m=1 (c m -σ i )∇ũ m on Ω i and g = - M m=1 c m ∂ n ũm on Γ ∞,h .
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 function g 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) follows 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 6. For a family of discrete problems (5.3), the bilinear forms a h (•, •) are uniformly V h -coercive if there is a constant α independent of h such that (5.5)

a h (v h , v h ) ≥ α v h 2 H 1 (Ω h ) ∀v h ∈ V h .
Notice that the family (a h (•, •)) h of bilinear forms defined by (5.4) is uniformly V h -coercive 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.

The discrete bilinear form a h (•, •) is naturally defined for elements in X h . In order to get error estimates of the discretization error, we need to extend this definition to elements in

L i=1 H 1 ( Ω i ), ∀ w = ( wi ) i=1,...,L ∈ L i=1 H 1 ( Ω i ), a h ( w, v h ) def = L i=1 ˆΩi,h σ i ∇ wi • ∇v h dx.
Similarly, we aim to define the continuous bilinear form a(•, •) as well as the linear form l(•) for elements v h ∈ X h which are only defined on Ω h . Without loss of generality, we may assume that any tetrahedron T ∈ T h has at most one face on the external boundary Γ ∞,h . Now, let T be a "boundary" tetrahedron and denote by f 1 h its face situated on Γ ∞,h . Denote by f 2 h , . . . , f 4 h the faces of T that are not contained in Γ ∞,h . The face f 1 h is the approximation of a curved surface f 1 ⊂ Γ ∞ . Then, let T be the curved tetrahedron with faces f 1 , f 2 h , . . . , f 4 h . Now, consider v h ∈ X h . On T , v h coincides with a polynom p of degree less or equal than 1. Then let us define v h on T in a unique way by

v h| T = p.
This definition allows to define a(w h , v h ) and l(v h ) for elements in X h .

Theorem 7. 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,

(5.6) ∀K ∈ h T h , h h K ≤ C inv . Let w = ( wi ) i=1,...,L ∈ L i=1 H 2 ( Ω i
) be an extension of w, the solution of problem (3.17), such that wi|Ω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,

L i=1 wi -w h H 1 (Ω i,h ) h M m=1 ∇ũ m H 1 (Ω\Vm) + ∂ n ũm H 1 (Γ∞) (5.7) +h 3/2 M m=1 ∂ n ũm W 1,∞ (Γ ∞,h ) + ∇ũ m H 1 (Γ ∞,h ) . (5.8)
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 7 relies on the following 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. It follows in the same way as in [START_REF] Ciarlet | The Finite Element Method for Elliptic Problems[END_REF] taking into account the extensions of the conductivity σ on the different subdomains Ω i .

Proposition 2. Consider a family of discrete problems of the form (5.3) for which 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 (5.9)

L i=1 wi -w h H 1 (Ω i,h ) inf v h ∈V h L i=1 wi -v h H 1 (Ω i,h ) + sup ξ h ∈V h |a h ( w, ξ h ) -l h (ξ h )| ξ h H 1 (Ω h )
where w = ( wi ) i=1,...,L is an arbitrary element of

L i=1 H 1 ( Ω i ).
According to the regularity of the solution w of (3.17) on the subdomains Ω i , we can state the following estimate of the interpolation error. Proposition 3. Consider a family (X h ) h of finite element spaces of type P1 on regular triangulations (T h ) h . Then, the following estimate holds true for any w = ( wi ) i=1,...,L ∈

L i=1 H 2 ( Ω i ). (5.10) wi -Π i,h wi H 1 (Ω i,h ) h wi H 2 ( Ωi) .
The next proposition is concerned with the estimate of the consistency error

|a h ( w,ξ h )-l h (ξ h )| ξ h H 1 (Ω h ) : Proposition 4. Let w = ( wi ) i=1,...,L ∈ L i=1 H 2 ( Ω i
) be an extension of w, the solution of problem (3.17), such that wi|Ω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.

Assume that the mesh family (T h ) h satisfies the inverse assumption (5.6).

Let w h ∈ X h be the discrete solution of problem (5.3) where X h is the space of linear finite elements on the mesh T h . Then,

|a h ( w, v h ) -l h (v h )| v h H 1 (Ω h ) h M m=1 ∇ũ m H 1 (Ω\Vm) + ∂ n ũm H 1 (Γ∞) (5.11) +h 3/2 M m=1 ∂ n ũm W 1,∞ (Γ ∞,h ) + ∇ũ m H 1 (Γ ∞,h ) Proof. We have a h ( w, v h ) -l h (v h ) = L i=1 ˆΩi,h ∩Ωi (σ i ∇w i -F i ) • ∇v h dx + ˆΩi,h \Ωi σ i ∇ wi -F i • ∇v h dx - ˆΓ∞,h g v h ds = ˆΩ (σ∇w -F ) • ∇v h dx - L i=1 ˆΩi\Ωi,h (σ∇w -F ) • ∇v h dx + L i=1 ˆΩi,h \Ωi σ i ∇ wi -F i • ∇v h dx - ˆΓ∞,h g v h ds = L i=1 ˆΩi,h \Ωi σ i ∇ wi -F i • ∇v h dx - L i=1 ˆΩi\Ωi,h (σ∇w -F ) • ∇v h dx + ˆΓ∞ g v h ds - ˆΓ∞,h g v h ds,
taking into account the variational equality (3.17) for w. The first two terms in the above expression are of order h according to the following estimates (see e.g. [START_REF] Raviart | Introduction à l'Analyse Numérique des Equations aux Dérivées Partielles[END_REF]):

σ i ∇ wi -F i L 2 (Ω i,h \Ωi) h σ i ∇ wi -F i H 1 ( Ωi) , (5.12a) σ i ∇w i -F i L 2 (Ωi\Ω i,h ) h σ i ∇w i -F i H 1 (Ωi) .
(5.12b) Notice further that the right hand sides in (5.12) depend continuously on the data F and g if we assume that the extension operator from Ω i to Ω i is continuous. Then,

L i=1 ˆΩi,h \Ωi σ i ∇ wi -F i • ∇v h dx - ˆΩi\Ωi,h (σ∇w -F ) • ∇v h dx h M m=1 ( ∇ũ m H 1 (Ω\Vm) + ∂ n ũm H 1 (Γ∞) ) v h H 1 (Ω h ) .
Here, the H 1 -norm of v h on Ω has been majored by v h H 1 (Ω h ) up to a multiplicative constant which is possible whenever the mesh parameter h is small enough (see e.g. [START_REF] Hernández | Finite element Approximation of Spectral Problems with Neumann Boundary Conditions on curved domains[END_REF]).

In order to estimate the boundary integrals, let f h ⊂ Γ ∞,h be one of the boundary faces. Then f h is the first order approximation of a curved surface f ⊂ Γ ∞ defined by the nodes of f h . Let ϕ be the parametrization of class C 2 that maps f h onto f . ϕ(s) -s represents the interpolation error with polynomials of degree 1.

Then,

ˆΓ∞ g v h ds - ˆΓ∞,h g v h ds = f h ⊂Γ ∞,h ˆfh ((g v h )(ϕ(s)) -(g v h )(s)) ds + ˆfh (g -g)(s) v h (s) ds = f h ⊂Γ ∞,h ˆfh (∇ s (g v h ))(s) • (ϕ(s) -s) + o( ϕ(s) -s ) ds + ˆfh (g -g)(s) v h (s) ds f h ⊂Γ ∞,h v h ∇ s g L 2 (f h ) + g ∇ s v h L 2 (f h ) ϕ -Π f h ϕ L 2 (f h ) + g -g L 2 (f h ) v h L 2 (f h ) .
It follows from the regularity of g that g W 1,∞ (f h ) is bounded independently from h. We then get from standard estimates in polynomial vector spaces that

v h ∇ s g L 2 (f h ) + g ∇ s v h L 2 (f h ) g W 1,∞ (f h ) v h H 1 (f h ) h -1/2 g W 1,∞ (f h ) v h H 1 (T )
provided that the meshes satisfy the inverse assumption (5.6). We further deduce from the classical interpolation estimates that

ϕ -Π f h ϕ L 2 (f h ) h 2
with a constant depending on |ϕ| H 2 (f h ) that can be chosen independent from h. It remains to estimate the term g -

g L 2 (f h ) v h L 2 (f h ) .
Recall that the exterior normal on f is given by n = -∇ s ϕ since ϕ is the parametrization of f . Thus,

g -g = M m=1 c m ∇ũ m • ∇ s ϕ -n |f h . Now, ∇ s ϕ -n |f h
is the gradient of the interpolation error between ϕ and its P1-interpolate. This yields

g -g L 2 (f h ) h M m=1 ∇ũ m L 2 (f h )
with a constant depending on |ϕ| H 2 (f h ) but independent from h. Finally, it follows from the continuity of the trace operator on the reference element that

v h 2 L 2 (f h ) meas(f h ) meas(T ) h 2 v 2 H 1 (T ) h v 2 H 1 (T )
provided the inverse assumption (5.6) holds true. This yields the following estimate

|a h ( w, v h ) -l h (v h )| v h H 1 (Ω h ) h M m=1 ∇ũ m H 1 (Ω\Vm) + ∂ n ũm H 1 (Γ∞) +h 3/2 M m=1 ∂ n ũm W 1,∞ (Γ ∞,h ) + h 3/2 M m=1 ∇ũ m H 1 (Γ ∞,h )
which proves (5.11).

5.2.

Numerical validation in the multi-layer spherical model. In this section, we consider a three-layer spherical head model (see Fig. 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 [START_REF] Lew | ainen, Effects of sutures and fontanels on MEG and EEG source analysis in a realistic infant head model[END_REF] (unless indicated otherwise).

We consider a family of three tetrahedral meshes with decreasing mesh size h (cf. Table 5.2). 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.8), is deduced from the discrete solution w h and decomposition (3.1). All simulations are executed with the software FreeFem++ [START_REF] Hecht | FreeFem++ Manual[END_REF]. The different linear systems are solved with the iterative solver GMRES with a tolerance equal to 10 -6 .

Two criteria are commonly used in numerical validation of EEG models [START_REF] Wolters | Numerical mathematics of the subtraction approach for the modeling of a current dipole in EEG source reconstruction using finite element head models[END_REF]. The first, called the Relative Difference Measure (RDM), is computed as follows (5.13) RDM

:= u h u h L 2 (Γ ∞,h ) - u ref u ref L 2 (Γ ∞,h ) L 2 (Γ ∞,h )
.

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

(5.14) MAG := 1 - u h L 2 (Γ ∞,h ) u ref L 2 (Γ ∞,h ) .
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 [START_REF] De Munck | A fast method to compute the potential in the multisphere model[END_REF][START_REF] Zhang | A fast method to compute surface potentials generated by dipoles within multilayer anisotropic spheres[END_REF]. 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 values 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 Fig. 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. The inclusion of the main fontanel in the three-layer spherical model is performed by the definition of a region Ω f ⊂ Ω 2 . In Ω 2 , the conductivity is defined by

Mesh

σ 2 (x) = σ skull +(σ f -σ skull )g(x) with σ skull = 0.04S.m -1 , σ f = 0.3S.m -1 and a gaussian function g(x) = e -α(x 2 1 +x 2 2 )
. The parameter in the definition of g is set to α = 10 4 which amounts to saying that the fontanel is limited (up to 1.5%) to the region 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 ref computed on the very fine reference mesh M ref (cf. Table 5.2). We define global errors on the whole domain Ω by

Ω f := {x = (x 1 , x 2 , x 3 ) ∈ Ω 2 : x 2 1 +x 2 2 ≤ L 2 } with L = 20mm (see Fig.
u h -u ref 2 H 1 (Ω) def = w h -w ref 2 H 1 (Ω) .
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. Left: one single source S = (0, 0, 40mm), q = (0, 0, J). Right: two sources S 1 = (0, 0, 10mm), S 2 = (0, 10mm, 0) with moments q 1 = (0, 0, J), q 2 = (0, J, 0).

Intensity J = 10 -6 A.m -2 .
We report in Fig. 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. [START_REF] Wolters | Numerical mathematics of the subtraction approach for the modeling of a current dipole in EEG source reconstruction using finite element head models[END_REF]). Here, δ denotes the distance of the sources from the interface brain/skull,

δ = max m=1:M dist(S m , Γ 1 )
which is related to the eccentricity e by e = 1 -δ r 1 .

Numerical discussion of the EEG models in neonates

This section is concerned with the comparison of different EEG models in neonates. We investigate the impact of considering the fontanels in neonate brain modeling through their conductivity σ f , as well as the effect of the skull conductivity σ skull . We use a panel of conductivities which are found in literature (e.g. [START_REF] Roche-Labarbe | High-resolution electroencephalography and source localization in neonates[END_REF][START_REF] Lew | ainen, Effects of sutures and fontanels on MEG and EEG source analysis in a realistic infant head model[END_REF][START_REF] Gargiulo | The effect of fontanel on scalp EEG potentials in the neonate[END_REF]). 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 values (σ f , σ skull ) on the measured electric potential at the scalp, we compute the factors RDM and MAG (cf. formula (5.13) and (5.14)) with respect to u h and u ref which are, respectively, the solutions to the fontanels and no-fontanel models.

We consider a realistic head model of a healthy newborn obtained from MR and CT images of the Amiens' hospital database (see Fig. 6.1). The segmented tissues correspond to brain (Ω 1 ), cerebrospinal fluid (Ω 2 ), skull (Ω 3 ), fontanels (Ω 3,f ⊂ Ω 3 ) and scalp (Ω 4 ). The diameter of the computational domain is about 12cm. The mesh properties of the model are reported in Table 6. 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 . Fig. 6.2 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 impact the EEG measurements in neonates through the ratio σ f /σ skull of tissue conductivities. Table 2. Four-layer realistic head model

Mesh Nodes Tetrahedra

Numerical sensitivity analysis

This section is devoted to study numerically how a slight variation of the conductivity affects EEG measurements. More precisely, we are interested in analyzing conductivity perturbations in neonatal skull including fontanels. We consider the realistic head model presented in Section 6 and fix the respective fontanel and skull conductivities to σ f = 0.3S.m -1 and σ skull = 0.04S.m -1 .

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 7.1 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 farther source position. 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 7.2 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 Fig. 7.3). 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 of the fontanels 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 Left: moment q = (0, J, J). Right: moment q = (J, J, 0). 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 impact seems to be localized near the fontanels and depends on the eccentricity and orientation of the dipolar sources.

Concluding remarks

In this paper, we propose a rigorous mathematical framework for an EEG model for neonates that is able to take into account the presence and ossification of fontanels and their specific electrical properties. To this end, the different tissues are assumed to be inhomogeneous and tissue conductivities are represented by functions of the position. Sources of electrical activity in the brain are represented as current dipoles which leads to an elliptic problem with a Dirac source term on the right-hand side. 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 tissue conductivities. Convergence estimates are obtained for standard finite elements of type P1.

Firstly, the numerical validation has been performed 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, we have studied the influence of the fontanel and skull conductivities on the electrical potential values at the scalp for a realistic neonatal head model. 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.

The analysis of the EEG forward model is an essential preliminary to the resolution of the corresponding inverse source problem. In the head model of adults, theoretical and numerical results exist [START_REF] Badia | A stable recovering of dipole sources from partial boundary measurements[END_REF][START_REF] Farah | Problèmes inverses de sources et lien avec l'Electro-encéphalo-graphie[END_REF]. compléter la biblio. Abdellatif. For the neonatal head model, identifiability and stability results for the inverse EEG source problem have been obtained parallel to the present work [START_REF] Diallo | An inverse dipole EEG source problem in neonates[END_REF]. 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.

(

2) Let v ∈ L p (Ω) with 1 < p < ∞. Assume that there is a constant C > 0 such that ∆ h k v L p (Ω ) ≤ C for any Ω ⊂⊂ Ω satisfying |h| < dist(Ω , ∂Ω). Then, the weak derivative ∂ k v exists and ∂ k v L p (Ω) ≤ C.

Proof. (of Theorem 2)

Recall that w ∈ H 1 (Ω) is the variational solution of the following problem,

(B.2) ˆΩ σ∇w • ∇v dx = ˆΩ F • ∇v dx + ˆΓ∞ gv ds ∀v ∈ H 1 (Ω)
where

F def = M m=1 (c m -σ)∇ũ m and g def = - M m=1 c m ∂ n ũm .
Interior regularity of w in the subdomains Ω i follows from classical results for elliptic equations with variable coefficients (see e.g. [START_REF] Gilbarg | Elliptic Partial Differential Equations of Second Order[END_REF]). Indeed, the regularity of the conductivity, σ i ∈ W 1,∞ (Ω i ), fits with the assumptions in [START_REF] Gilbarg | Elliptic Partial Differential Equations of Second Order[END_REF] since W 1,∞ (Ω i ) can be identified with the space of uniformly Lipschitz continuous functions on Ω i . Therefore, on any subdomain Ω i , w is the weak solution of the partial differential equation -∇ • (σ i ∇w) = f where the right hand side

f i = M m=1 ∇ • ((σ i -c m )∇ũ m ) = M m=1 (∇(σ i -c m ) • ∇ũ m + (σ -c m )∆ũ m ) = M m=1 ∇(σ i -c m ) • ∇ũ m belongs to L 2 (Ω i ) for any i = 1, . . . , L since ∇(σ i -c m ) ∈ L ∞ (Ω i )
and ũm is regular on the support of σ i -c m . Hence, for any Ω i ⊂⊂ Ω i , we have w ∈ H 2 (Ω i ) and the estimate

(B.3) w H 2 (Ω i ) w H 1 (Ωi) + ∇ũ L 2 (Ωi)
holds true with a constant depending on σ and on the distance d = dist(Ω i , ∂Ω i ).

In order to prove H 2 -regularity up to the boundary of the subdomains, we adapt the technique from [START_REF] Gilbarg | Elliptic Partial Differential Equations of Second Order[END_REF] to Neumann boundary conditions. Let x 0 ∈ Γ ∞ be a point situated on the exterior boundary of Ω. Let Ψ be a C 2 -diffeomorphism satisfying Condition 1 and let

B + = B R (x 0 ) ∩ Ω, D = Ψ(B R (x 0 )), D + = Ψ(B + ), and T = Ψ(B R (x 0 ) ∩ Γ ∞ ). Now, consider the vector space H 1 0,∂B + \Γ∞ (B + ) of functions in H 1 0 (B R (x 0 )) restricted to B + . Any function v ∈ H 1 0,∂B + \Γ∞ (B + ) satisfies γ 0 v = 0 on ∂B + \ Γ ∞ ,
ˆD+ 3 i,j=1 âij (ξ)∂ j ŵ(ξ)∂ i v dξ = ˆD+ F (ξ) • ∇v(ξ) dξ + ˆT∞ ĝ( ξ)v( ξ) d ξ, ∀v ∈ H 1 0,∂D + \T (D + ).
Here, the coefficients âij (ξ) ∈ W According to (B.5), we get

ˆD+ 3 i,j=1 âij (ξ + he k )∆ h k ∂ j ŵ(ξ)∂ i v(ξ) dξ = ˆD+ ∆ h k F • ∇v dξ + ˆT ∆ h k ĝv d ξ - ˆD+ 3 i,j=1 ∆ h k âij (ξ)∂ j ŵj (ξ)∂ i v(ξ) dξ ≤ ∆ h k F L 2 (D + ) ∇v L 2 (D + ) + ∆ h k ĝ L 2 (T ) v L 2 (T ) + 3 i,j=1 ∆ h k âij L ∞ (D + ) ∂ j ŵ L 2 (D + ) ∂ i v L 2 (D + )   ∂ k F L 2 (D + ) + ∂ k ĝ L 2 (T ) + i,j=1 ∂ k âij L ∞ (D + ) ŵ H 1 (D + )   ∇v L 2 (D + )
where the last estimate follows from (B.1) in Proposition 5. Notice also that v From the ellipticity property (B.10) and the boundedness of the coefficients âij in W

∂ k F L 2 (D + ) + ∂ k ĝ L 2 (T ) + ŵ H 1 (D + ) η∇(∆ h k ŵ) L 2 (D + ) + ∂ k F L 2 (D + ) + ∂ k ĝ L 2 (T ) + ŵ H 1 (D + ) ŵ H 1 (D + ) F H 1 (D + ) + ĝ H 1 (T ) η∇(∆ h k ŵ) L 2 (D + ) + F H 1 (D + ) + ĝ H 1 (T ) 2 (B.7)
according to the standard a-priori estimate for the solution of (B.4), ŵ F L 2 (D + ) + ĝ L 2 (T ) .

Notice that the constant in (B.7) now depends on âij and η and their derivatives.

Young's inequality 2ab ≤ εa 2 + 1 ε b 2 applied to the first term of (B.7) yields the following estimate for η∇(∆ h k ŵ) L 2 (D + ) , η∇(∆ h k ŵ) L 2 (D + ) ≤ C(â ij , ∇η) F H 1 (D + ) + ĝ H 1 (T ) .

According to Proposition 5.2, this proves that ∂ k ŵ belongs to H 1 (D + ) for k = 1 and k = 2. The remaining second derivative ∂ 2 33 ŵ can be estimated directly from the partial differential equation

- 3 i,j=1 ∂ i (â ij ∂ j ŵ) = -∇ • F ∈ L 2 (D + ).
Hence, ŵ belongs to the space H 2 (D + ) and satifies the estimate

∇ ŵ H 1 (D + ) ≤ C(â ij , η) F H 1 (D + ) + ĝ H 1 (T ) .
Going back to the original domain B + ⊂ Ω L , we get w = ŵ • Ψ ∈ H 2 (B ρ ∩ Ω L ) for any ρ < R.

Finally, the estimate

(B.8) w H 2 (Bρ∩Ω L ) ≤ C(σ, Ψ) ∇ũ H 1 (Bρ∩Ω L ) + ∂ n ũ H 1 (Bρ∩Γ∞)
follows from (B.13) and (B.14). For a point x 0 situated on an interface Γ i , the same technique applies with ĝ = 0. This yields w ∈ H 2 (B ρ ∩ Ω i ) for any ρ < R. Since w ∈ H 2 loc (Ω i ) by interior regularity, we infer that w ∈ H 2 (Ω i ). Choosing a finite number of points x

Figure 2 . 1 .

 21 Figure 2.1. Three-layer head model.

( 3 .Lemma 1 .

 31 10) a(w, v) = ˆΩ σ∇w • ∇v dx, as well as the linear form (3.11) l(v) = M m=1 ˆΩ(c m -σ)∇ũ m • ∇v dx -ˆΓ∞ c m ∂ n ũm v ds defined for any v ∈ H 1 (Ω). Since (3.9) is a problem with Neumann boundary condition, the variational formulation (3.8) allows a solution only under the compatibility condition (3.12) l(1) = M m=1 ˆΓ∞ c m ∂ n ũm (x) ds = 0. Condition (3.12) follows from the following lemma by summing over m. Let ũm be the solution of equation (3.2) given by (3.3). Then (3.13) ˆΓ∞ c m ∂ n ũm ds = 0 ∀m = 1, . . . , M.

  w def = w(•, σ) and w h def = w(•, σ + hµ) as well as ũm def = ũm (•, σ) and ũm,h def = ũm (•, σ + hµ) Now, recall that w h and w are the respective solutions in V of (4.6)

Figure 5 . 1 .

 51 Figure 5.1. Behavior of factors RDM and MAG with respect to the eccentricity of the dipole. Different mesh sizes (finest mesh M 3 ). Neonatal three-layer spherical head model without fontanels. Exact reference solution.

  5.2). Hence, the Gaussian g allows to model the process of fontanel ossification and leads to a continuous function σ 2 . Quite different values may be found in literature for both the fontanel conductivity and the neonatal skull conductivity. A possible parameter set of tissue conductivities and the resulting EEG models in neonates are discussed in Section 6.

Figure 5 . 2 .

 52 Figure 5.2. A spherical head model with the main fontanel.

Figure 5 . 3 .

 53 Figure 5.3. Errors in H 1 -norm with respect to the mesh size h in logarithm scale. Three-layer spherical head model with the anterior fontanel (Gaussian behavior for the fontanel conductivity). Numerical reference solution computed on M ref .Left: one single source S = (0, 0, 40mm), q = (0, 0, J). Right: two sources S 1 = (0, 0, 10mm), S 2 = (0, 10mm, 0) with moments q 1 = (0, 0, J), q 2 = (0, J, 0). Intensity J = 10 -6 A.m -2 .

Figure 5 . 4 .

 54 Figure 5.4. Behavior of factors RDM and MAG with respect to the eccentricity dipole position for different meshes. Three-layer spherical model with the anterior fontanel (Gaussian behavior for the fontanel conductivity). Numerical reference solution computed on M ref .

Figure 6 . 1 .

 61 Figure 6.1. Realistic head model of a neonate. Left: skull and fontanels. Right: mesh of the fontanels.

  Boundary faces h min [m] h max [m]

Figure 6 . 2 .

 62 Figure 6.2. Variations of factors RDM and MAG with respect to different conductivities (σ f , σ skull ). Four-layer realistic head model. Reference solution computed with the the model without fontanels.

Figure 7 . 1 .

 71 Figure 7.1. Sensitivity of the electric potential on the scalp with respect to eccentricity. Distance source-interface brain/CSF ≈ 5mm (left) and ≈ 15mm (right).

Figure 7 . 2 .

 72 Figure 7.2. Sensitivity of the electric potential on the scalp with respect to orientation. Distance source-interface brain/CSF ≈ 15mm.Left: moment q = (0, J, J). Right: moment q = (J, J, 0).

Figure 7 . 3 .

 73 Figure 7.3. Sensitivity of the electric potential on the scalp for a deep source.

  and may be extended by zero on Ω to an admissible test function in (B.2). For v ∈ H 1 0,∂B + \Γ∞ (B + ), the variational equation (B.2) may be transformed into a similar equation on D + in the following way. Let ŵ = w•Ψ -1 ∈ H 1 (D + ) and v = v • Ψ -1 ∈ H 1 0,∂D + \T (D + ) with v ∈ H 1 0,∂B + \Γ∞ (B + ). According to Lemma 3, ŵ satisfies the equation (B.4)

  1,∞ (D + ) satisfy the coercivity condition (B.10). The functionsF : D + → R 3 and ĝ : T → R are respectively of class W 1,∞ (D + ) and C 1 . Now, let v ∈ H 1 0 ( D) with supp(v) ⊂⊂ D. Then v|D + belongs to the space H 1 0,∂D + \T (D + ). Let |h| < dist(supp(v), ∂ D). For k = 1 or k = 2, we consider the differential quotient of v in the direction e k , ∆ h k v(ξ) = v(ξ + he k ) -v(ξ) h . Notice that ∆ h k v belongs to H 1 0,∂D + \T (D + ) since k = 3 and thus ξ ± he k ∈ T for any ξ = ( ξ, 0) ∈ T with dist(ξ, ∂T ) > |h|. Then, âij ∂ j ŵ ∂ i ∆ -h k v dξ = -ˆD+ F • ∇ ∆ -h k v dξ -ˆT ĝ∆ -h k v d ξ = ˆD+ ∆ h k F • ∇v dξ + ˆT ∆ h k ĝv d ξ, (B.5)noticing that ∆ -h k v defines an admissible test function in (B.4). On the other hand, we have âij (ξ + he k )∆ h k ∂ j ŵ(ξ) + ∆ h k âij (ξ)∂ j ŵ(ξ) .

  L 2 (T ) ∇v L 2 (D + ) due to the trace theorem and Poincaré's inequality since v = 0 on ∂D + \ T . All terms of the right hand side of (B.6) are bounded according to the regularity of the data F , ĝ and âij . Now, let η ∈ C 1 0 ( D) be a cut-off function with 0 ≤ η ≤ 1 and η ≡ 1 on D ⊂⊂ D. Let |h| < dist(supp(η), ∂ D) and take k = 1 or k = 2. Then, the function v = η 2 ∆ h k ŵ ∈ H 1 0,∂D + \T (D + ) is an admissible test function in (B.4).

0Lemma 3 .

 3 and appropriate balls B ρ (x (j) 0 ) that recover L-1 i=1 ∪Γ ∞ , we obtain (3.20) from (B.3) and (B.8). The following Lemma describes the effect of the transformation Ψ from the physical domain B + to a straightened domain D + on the (localized) variational formulation (B.4). Let x 0 ∈ ∂Ω i , i ∈ {1, . . . , L}. Let Ψ : B → D be a C 2 -diffeomorphism defined on a ball B = B(x 0 ) and satisfying Condition 1. For B R (x 0 ) ⊂⊂ B, set B + = B R (x 0 ) ∩ Ω i , D = Ψ(B R (x 0 )), D + = Ψ(B + ), and T = Ψ(B R (x 0 ) ∩ ∂Ω i ). For any function v ∈ H 1 (B + ), we denote by v the function in H 1 (D + ) defined by v = v • Ψ -1 . Under the assumptions of Theorem 2, we have (B.9)ˆB+ σ(x)∇w(x) • ∇v(x) dx = ˆD+ 3 i,j=1 âij (ξ)∂ j ŵ(ξ)∂ i v(ξ) dξ where âij ∈ W 1,∞ (D + ) and âij (ξ)d j d i ∀d ∈ R 3 , ∀ξ ∈ D +with a constant λ min > 0 independent from ξ. Furthermore, there are functionsF ∈ W 1,∞ (D + ) and ĝ ∈ C 1 (T ) such that (B.11) M m=1 ˆB+ (c m -σ)∇ũ m • ∇v = ˆD+ F (ξ) • ∇v dξ and (B.12) -M m=1 ˆBR (x0)∩Γ∞ c m ∂ n ũm v ds = ˆT ĝ( ξ)v( ξ) d ξ where ξ = ( ξ, ξ 3 ) with ξ ∈ T . F and ĝ satisfy the following estimates in terms of the data ũ = M m=1 ũm , (B.13) F H 1 (D + ) ≤ C(Ψ, σ) ∇ũ H 1 (B R (x0)∩Ωi)

  4.13a) since |p m | ≤ µ ∞ = 1 and |c m + hp m | ≥ σ min . (4.13b) follows by integration of the normal derivative.

	Next, recall that ũ1 m = -pm cm ũm according to Proposition 1. Together with identity (4.15), we
	thus get
	(4.16)

Table 1 .

 1 Definition of meshes (neonatal three-layer spherical head model)

		Nodes	Tetrahedra Boundary nodes h min [m] h max [m]
	M 1	102 540	594 907	16 936	8.16 10 -4 4.81 10 -3
	M 2	302 140	1 855 005	23 339	6.35 10 -4 3.07 10 -3
	M 3	596 197	3 632 996	54 290	4.1 10 -4 2.46 10 -3
	M ref 2 754 393 17 263 316	124 847	2.5 10 -4 1.51 10 -3

  1,∞ (D + ), F L 2 (D + ) + ∂ k ĝ L 2 (T ) + ŵ H 1 (D + ) ∇v L 2 (D + ) +2 sup D |∇η| η∇(∆ h k ŵ) L 2 (D + ) ∆ h k ŵ L 2 (D + )with a constant depending on âij . But∇v L 2 (D + ) η∇(∆ h k ŵ) L 2 (D + ) + 2 sup D |∇η| ŵ H 1 (D + ) and ∆ h k ŵ L 2 (D + ) is bounded by ŵ H 1 (D + ) . Therefore, η∇(∆ h k ŵ) 2 L 2 (D + )

	we get
	ˆD+	3
	η∇(∆ h k ŵ) 2 L 2 (D + )	i,j=1	âij (ξ + he k )∂ j ∆ h k ŵ	η 2 ∂ i ∆ h k ŵ ∂i v-2η(∂iη)∆ h k	ŵ dξ
	∂ k

Acknowledgements. We are very grateful to H. Azizollahi of the group GRAMFC INSERM U1105 at Amiens' hospital who provided us with the mesh for the realistic head model.

Appendix A

For completeness, we give the details for the proof of equivalence between the variational formulations (3.7) and (3.8):

In this section we prove global H 2 -regularity for the variational solution of problem (3.17) on each subdomain Ω i . According to the regularity assumptions on Ω i , any appropriate neighborhood N of a point x 0 ∈ ∂Ω i can be straightened by a regular mapping satisfying the following conditions:

In the special case where i = 1, we further assume that B R (x 0 ) ∩ {S m } = ∅ for any m = 1, . . . , M .

The proof relies on estimates of the differential quotient of a function v in the direction of the canonical basis vector e k , k = 1, . . . , 3, i.e.

For the sake of completeness, we recall the following properties of ∆ h k : Proposition 5.

(1) Let v ∈ W 1,p (Ω). Then ∆ h k v ∈ W 1,p (Ω ) for any Ω ⊂⊂ Ω satisfying |h| < dist(Ω , ∂Ω), and the following estimate holds true:

and

Proof. The substitution rule for integrals yields (B.9) with

The regularity of âij is now given by the C 0,1 -regularity of σ on B + since DΨ and DΨ -1 are of class C 1 . For any ξ ∈ D + , the matrix A(ξ (B.11) may be obtained in a similar way. First, let F (x) = M m=1 (c m -σ)∇ũ m on B + . The regularity of F = (F 1 , F 2 , F 3 ) is given by the regularity of σ on B + since ũm is C ∞ on supp(c m -σ). Therefore, F k ∈ W 1,∞ (B + ) (or equivalently C 0,1 (B + )). Applying the substitution rule yields (B.11) with

Finally, define the function g = -

|T is the parametrization of the boundary B R (x 0 ) ∩ Γ ∞ and n = ∇ϕ. Since, ∇ũ m is C ∞ on Γ ∞ , the regularity of the function g is given by the regularity of ∇ϕ, i.e. C 1 . Once again, we obtain (B.12) from the substitution rule with