M Darbas

S Lohrengel

B Sulis

Forward and inverse source problems for time-dependent electroencephalography

Keywords: time-dependent electroencephalography, dipolar sources, neuronal model, inverse source problem, moving source points. 2020 Mathematics Subject Classification: 92C50, 92-10, 35R30, 65N21

A new mathematical model for time-dependent electroencephalography (EEG) is developed and analysed. Evolution with time is introduced into the standard EEG model by considering dipolar sources with time-dependent moments and source positions. Dimensional analysis shows the validity of the quasi-stationary approximation for all tissues of the humain head. Non-linear systems of differential equations based on gating particles are used to model the postsynaptic current at the neuron level which, in turn, yields the dipolar source term of the boundary value problem. The wellposedness of the forward time-dependent EEG problem is proved by the subtraction approach for moments with L 2 -regularity in time and continuous source trajectories. Numerical results explain the pipeline from the simulation of the postsynaptic current up to the potential recorded at the electrodes in a 2D circular configuration and on the three-dimensional realistic head model of a neonate. The inverse source problem is formulated with the help of a time-dependent non-linear measurement operator and identifiability and stability results are proven. It is numerically solved by the Minimum Norm Estimate and the computation of the involved Lead Field Matrix is explained for the particular case of the subtraction approach. The reconstruction of the trajectory of a moving source point with time-dependent moment illustrates the approach for the inverse problem in the 2D configuration.

Introduction

Electroencephalography (EEG) is one of the most widespread functional brain imaging techniques. Measurements of the electric potential generated by normal or pathological brain activity are taken at electrodes placed at the surface of the scalp. They record in a passive and non-invasive way the voltage potential fluctuations between different cortical regions. EEG-monitoring can be done at the bedside of the patient which makes the technique particularly appropriate for neonates and premature babies. The important goal of brain imaging using EEG is to localize cerebral sources generating measured EEG signals. EEG is known to have an excellent resolution in time and is able to record neural events in order of one millisecond. It is one of the main tools in presurgical evaluation for refractive epilepsy since the measurements contain valuable information about the 1 INTRODUCTION 2 localization of active brain areas. Neonatal EEG is used to assess seizure recognition and classification and to make epilepsy syndrome diagnoses. It provides also prognostic information for other brain dysfunctions (e.g [START_REF] Koskela | Prognostic value of neonatal EEG following therapeutic hypothermia in survivors of hypoxic-ischemic encephalopathy[END_REF]).

The development of multimodal analysis, coupling EEG with other imaging modalities, is likely to improve electrical source analysis since it gives additional information and could allow to detect activity which is less visible in EEG. An example is the study of [START_REF] Roche-Labarbe | NIRS-measured oxy-and deoxyhemoglobin changes associated with EEG spike-andwave discharges in children[END_REF] where hemodynamic changes among 10s before an epileptic seizure have been observed. Different strategies can be considered, from coupling at the image level allowing to constrain further the solution of the EEG inverse problem, up to the physiological level where the evolution in time of the significative variables as cerebral blood volume, for instance, are described by differential equations. We are particularly interested in the coupling of EEG with optical imaging, taking advantage of the possible coregistration of EEG and NIRS (Near InfraRed Spectroscopy) which has already been operated by collaborators from GRAMFC INSERM UMR-S 1105 (Amiens' hospital) [START_REF] Nourhashemi | Neurovascular coupling in the developing neonatal brain at rest[END_REF] for neonates and premature babies. The present work is a first step in the mathematical modeling of such a coupling. Indeed, a time-dependent model for EEG is required in order to take into account different time scales and gaps between signals related to the same event. The time-dependent EEG model and its mathematical analysis is an essential brick towards a complete model for neurovascular coupling which describes the relationship between electrical brain activity and the hemodynamic response. To this end, we propose to get deeper insight in the link between sources of (normal or pathological) electrical brain activity and the resulting EEG pattern over a given time interval.

From a mathematical point of view, source localization is an inverse problem that relies heavily on the chosen forward model for the electromagnetic field [START_REF] Akalin Acar | Effects of Forward Model Errors on EEG Source Localization[END_REF]. Commonly, the following elliptic equation is considered for the forward problem

-∇ • (σ∇u) = ∇ • j p in Ω,
σ∂ n u = 0 on ∂Ω.

(

Here, the unknown u represents the electric potential in the head occupying the region Ω, σ is the electrical conductivity of the head tissues and j p denotes the primary current produced at the neural synapses.

As mentioned above, EEG measures phenomena that are time-dependent by nature [START_REF] Ponton | Spatio-temporal source modeling of evoked potentials to acoustic and cochlear implant stimulation[END_REF][START_REF] Makeig | Dynamic brain sources of visual evoked responses[END_REF]. This is true at the source level as well as for the electric potential: as an example, we cite the recent study of [START_REF] Rezaei | Reconstructing subcortical and cortical somatosensory activty via the RAMUS inverse source analysis technique using median nerve SEP data[END_REF] where somatosensory evoked potentials have been analyzed in a time window between 14 and 30 ms after stimulation. Whereas the earlier components of the signal correspond to deeper activity in the brainstem, later peaks reveal both cortical and subcortical activity.

But the forward model (1.1) is static. In this paper, we thus aim to answer the question: which is the right mathematical model for time-dependent EEG? Usually, (1.1) is derived from the full Maxwell equations by neglecting the time-derivative of the electric and magnetic field. Here, we conduct a dimensional analysis in the spirit of [START_REF] Rapetti | On quasi-static models hidden in Maxwell's equations[END_REF] in order to justify the absence of time-derivatives in an EEG model with time-dependent source terms. We then propose a modeling for the dynamics of the primary current j p at the neuron level. We have in mind certain situations (e.g. epileptic seizure or obstructive sleep apnea) happening during a short period of time. At these moments, the current j p can be net and strong. Furthermore, these episodes are generated and concentrated in a relatively small area of the brain.

We provide a rigorous mathematical formulation of the forward problem in the presence of dipolar sources with moving source points and dipolar moments with low regularity in time. To this end, we generalize the subtraction approach which is a well-known technique for the static problem (1.1) [START_REF] Badia | An inverse source problem in potential analysis[END_REF][START_REF] Wolters | Numerical mathematics of the subtraction method for the modelling of a current dipole in EEG source reconstruction using finite element head models[END_REF][START_REF] Azizollahi | EEG in neonates: forward modeling and sensitivity analysis with respect to variations of the conductivity[END_REF].

With regard to inverse source analysis, biomedical software generally brings already the possibility of solving the inverse problem step by step in time when the forward problem is formulated with the help of a lead-field matrix [START_REF] Grech | Review on solving the inverse problem in EEG source analysis[END_REF]. Spatio-temporal source analysis considers fixed source positions but time-varying amplitude of the moments. [START_REF] Scherg | Fundamentals of dipole source potential analysis[END_REF][START_REF] Ponton | Spatio-temporal source modeling of evoked potentials to acoustic and cochlear implant stimulation[END_REF], for example, use time-dependent scalp wave forms at the electrodes which are fitted to the recorded data . We can also cite probabilistic spatio-temporal models which embodies temporal and spatial priors in the source localization [START_REF] Trujillo-Barreto | Bayesian M/EEG source reconstruction with spatio-temporal priors[END_REF][START_REF] Hashemi | Efficient Hierarchical Bayesian Inference for Spatio-temporal Regression Models in Neuroimaging[END_REF]. In these approaches, the lead-field matrix is computed using the numerical resolution of the static forward problem (1.1). Here, we adopt the PDE's point of view based on a time-dependent forward problem to address the question of uniqueness for the inverse problem in the case of moving dipolar sources with moments of low regularity in time.

The paper is organized as follows. In Section 2, we derive a time-dependent EEG model. In Section 3, we study the forward problem. We address an existence and uniqueness result as well as numerical simulations of EEG measurements. Section 4 is devoted to the mathematical analysis of the time-dependent inverse source problem. We prove identifiability and stability results. Finally, numerical reconstruction of moving sources with time-dependent moments is performed in Section 5.

Modeling time-dependent EEG

Dimensional analysis of Maxwell's equations

We study an electromagnetic phenomenon in a space-time domain with characteristic length ℓ and time scale τ . The continuous medium is characterized by its electric permittivity ε, magnetic permeability µ and conductivity σ. Maxwell's equations describe the propagation of an electromagnetic field in this medium. Under the assumption that the constitution laws of the medium are linear, they are given by

∇ • B = 0, ∇ × B = µ(ε∂ t E + J), ∇ × E = -∂ t B, ∇ • E = ρ ε , ( 2.1)
and model the interaction between the electric field E, the magnetic induction B, the charge density ρ and the current density J. In a conducting medium with conductivity σ, Ohm's law states that the free current density is given by J = σE. In the presence of impressed current sources, a source term j p should be added to J. The static limit of Maxwell's equations is well understood in the case where all fields and sources are time-independent. As soon as there is some time-dependence, however, the full Maxwell system should apply unless mathematical analysis shows that some terms can be neglected. Electroencephalography measures the cerebral activity during the observation. This activity is clearly variable in time, meaning that time-dependent source terms should be considered. In this context, dimensional analysis is an interesting tool to reduce the set of Maxwell's equations and get simplified models under particular hypotheses. The idea is to observe the different physical quantities in terms of their units or dimensions. Smallness of some of them allows to neglect certain coupling between electric and magnetic fields [START_REF] Rapetti | On quasi-static models hidden in Maxwell's equations[END_REF]. To this end, we introduce the velocity of the system with modulus v = ℓ/τ . The light celerity in the medium with electromagnetic parameters ε and µ is given by c = 1/ √ εµ. Here, we are interested in limit configurations where v ≪ c.

For going further into details in the analysis, a first step is to express the quantities ℓ and c in terms of the electromagnetic parameters ε, µ and σ. Following [START_REF] Rapetti | On quasi-static models hidden in Maxwell's equations[END_REF], we introduce the quantity τ e = ε/σ which has dimension of time and can be interpreted as the electric charge diffusion time. Similarly, the quantity ℓcµσ can be shown to be dimensionless and can be written as the quotient of τ m = µσℓ 2 and τ em = ℓ/c which both have dimension of time. Whereas τ m can be associated with the time during which the electic field penetrates into the conducting medium, the quantity τ em corresponds to the time required for fields to propagate as an electromagnetic wave over a distance ℓ at the speed c. Comparing the order of magnitude of these characteristic times for the electromagnetic parameters of the different head tissues gives the range of validity of simplified models for Maxwell's equations in EEG.

We next proceed to the scaling of Maxwell's equations. To this end, set E(t, x) = eE ′ (t ′ , x ′ ) and B(t, x) = bB ′ (t ′ , x ′ ) where e, b are reference quantities and E ′ , B ′ are dimensionless quantities of order O(1) depending on t ′ = t/τ and x ′ = x/ℓ. In the sequel, prime notation corresponds to dimensionless operators or variables. Thus, Faraday's and Ampère's laws become respectively

∇ × E = -∂ t B ⇔ e ℓ ∇ ′ × E ′ = - b τ ∂ t ′ B ′ ⇔ ∇ ′ × E ′ = - ℓ τ b e ∂ t ′ B ′ , ( 2.2)
and

∇ × B = µ(ε∂ t E + σE) ⇔ b ℓ ∇ ′ × B ′ = µ(ε e τ ∂ t ′ E ′ + σeE ′ ) ⇔ ∇ ′ × B ′ = ℓ τ εµe b ∂ t ′ E ′ + µσℓe b E ′ . (2.3)
Using the definition of the characteristic times in (2.2) and (2.3), we get

∇ ′ × E ′ = - τ em τ cb e ∂ t ′ B ′ , ( 2.4)
and

∇ ′ × B ′ = τ em τ e cb ∂ t ′ E ′ + τ em τ e e cb E ′ . (2.5)
We are interested in the low frequency range where the characteristic dimensions of the system are such that τ em ≪ τ, τ e ≪ τ, and τ m ≪ τ.

(

Then, we cannot have simultaneously

τ em τ cb e = O(1) and τ em τ e cb = O(1)
and at least one of the time derivatives in (2.4) or (2.5) has to be neglected. We aim to determine the limit model from the only order between the characteristic time scales since these are available in a medium with given parameters ε, µ and σ for a configuration with known length ℓ and duration τ . One may notice that τ 2 em = τ e τ m which amounts to saying that any order between τ em and τ e induces an order of τ m with respect to τ em .

We distinguish the following three cases:

1. τ em ≪ τ e : According to assumption (2.6) on τ e and the relation between τ m , τ e and τ em , this implies the order

τ m ≪ τ em ≪ τ e ≪ τ.
From τ e ≪ τ , we deduce that τ em τ ≪ τ em τ e and the first term in the right hand side of (2.5) is small compared to the second one. The displacement current can thus be neglected and we get Ampere's law

∇ × B = µσE. (2.7)
Since the left hand side of (2.5) is of order 1, we further have

cb e ∼ τ em τ e ≪ 1.
Hence, the right hand side in (2.4) can be neglected, too, which yields

∇ × E = 0. ( 2
∇ • B = 0, ∇ × B = µJ, ∇ × E = 0, ∇ • E = ρ ε , ( 2.9)
called the quasi-stationary model.

Quasi-stationary approximation for EEG

The region of interest in EEG is the human head and the different tissues within it. In mathematical modeling of brain activity, one assumes that the permeability of head tissues is close to the one of free space. The parameters of the medium are given by ε = ε r ε 0 and µ = µ 0 where ε 0 = 8.85 × 10 -12 F.m -1 and µ 0 = 4π × 10 -7 H.m -1 are respective vacuum permittivity and permeability. The relative permittivity ε r and the conductivity σ depend on the tissue, but also on the frequency range. The typical frequency range of the electromagnetic fields in the head is inferior to 100 Hz [START_REF] Hämäläinen | Magnetoencephalography: theory, instrumentation, and applications to noninvasive studies of the working human brain[END_REF][START_REF] Faugeras | The inverse EEG and MEG problems: The adjoint state approach I: The continuous case Inria version 1[END_REF][START_REF] Gramfort | Mapping, timing and tracking cortical activations with MEG and EEG: Methods and application to human vision[END_REF] and we thus take τ ≈ 0.01s as the characteristic time scale. We consider a multilayer head model which distinguishes between white and grey matter and takes into account the cerebrospinal fluid (CSF) as well as skull and scalp. Uncertainty about the electric parameters should be taken into account since in vivo measurements in the human body are in general not available, especially for neonates and premature babies. In Table 1, we take the values for the electric permittivity from the IT'IS data base [23] which are given for a specific frequency from a parametrization fit based on the dispersion relation owing to [START_REF] Gabriel | Compilation of the Dielectric Properties of Body Tissues at RF and Microwave Frequencies[END_REF]. The conductivity values are taken from [START_REF] Lew | Effects of sutures and fontanels on MEG and EEG source analysis in a realistic infant head model[END_REF][START_REF] Azizollahi | Effects of Uncertainty in Head Tissue Conductivity and Complexity on EEG Forward Modeling in Neonates[END_REF] for neonates. We recall the formulae for the three characteristic times

τ e = ε/σ, τ m = µσℓ 2 , τ em = √ τ e τ m .
The characteristic length scale ℓ should be given by the thickness of the tissue in consideration. Notice, however, that τ e is independent from the characteristic length ℓ whereas τ m ∼ ℓ 2 and τ em ∼ ℓ. Hence, if the quasi-static approximation is valid for a given ℓ, it is also valid for any smaller length scale. We thus take for ℓ the dimensions of the neonatal head, i.e. ℓ ≈ 0.12 m. We then deduce the values of the three characteristic times. We have τ m ≪ τ em ≪ τ e ≪ τ for white matter, grey matter and skull. In CSF, the order is reversed and in the scalp, the three characteristic times are approximately of the same order with τ m ∼ τ em ∼ τ e ≪ τ . These situations correspond to the different cases studied in the previous dimensional analysis. Thus, for modeling EEG, the quasistationary approximation (2.9) of the full Maxwell's equations can be considered as a valid model. We focus on values for neonates and premature babies, but the conclusions of the dimensional analysis hold true for adults as well. In that case, the dimension of the head and the skull conductivity are replaced with ℓ ≈ 0.25 m and σ = 0.0043 S/m respectively. This yields an electric charge diffusion time τ e ≈ 1.3 × 10 -5 s which is still small compared to τ = 0.01s.

Equation ∇ × E = 0 shows that the electric field E derives from a scalar electric potential, denoted by u, and is written as

E(t, x) = -∇u(t, x), (t, x) ∈ (0, T ) × Ω, (2.10)
with T > 0 an observation time and Ω ⊂ R d (d = 2, 3) a head model. The currents in the extracellular region of the brain can be described as a homogeneous conductor, namely J = σE. In the brain and in particular in the cortex, the synchronized effect among a multitude of neurons creates an intracellular current denoted by j p . The current density J produced by cerebral activity thus splits into two terms

J = σE + j p . ( 2
= ∇ • (∇ × B) = µ∇ • (j p -σ∇u), (2.12)
which gives the following elliptic equation

∇ • (σ∇u) = ∇ • j p . (2.13)
In order to numerically reproduce spatially localized phenomena which are variable in time, we propose a time-dependent model for the source term j p in the form of a set of M electric dipoles located in the brain

j p (t, •) = M m=1 q m (t)δ(• -S m (t)). (2.14)
Here, S m (t) ∈ Ω and q m (t) ∈ R d are, respectively, the position and the moment of the m-th source at time t which is situated in the subdomain of Ω that models the brain or, more precisely, the grey matter. The right hand side of (2.13) then reads

F (t, •) := ∇ • j p (t, •) = M m=1 q m (t) • ∇δ(• -S m (t)). (2.15)
Assuming that the conductivity of air is zero and that no electric current can flow out of the scalp, the electric potential u is then solution of the following boundary problem with homogeneous Neumann condition

∇ • (σ∇u) = F in (0, T ) × Ω, σ∂ n u = 0 on (0, T ) × ∂Ω. (2.16)
Notice that (2.16) includes the following transmission conditions at any interface Γ p between two subdomains Ω p and Ω q

[u] |Γp = [σ∂ n u] |Γp = 0 on (0, T ) × Γ p .
(2.17)

Here, [f ] |Γp = f |Ωq -f |Ωp denotes the jump across the interface Γ p of the function f defined on Ω. Neurons are responsible for transmitting information between brain cells, through electrical impulses and chemical signals. The neuron is composed by three parts, the soma or cell body, the dendrites and the axon. When a stimulus occurs, the dendrites transmit the message to the soma which generates an action potential. The action potential is in turn sent through the axon and arrives at the presynaptic neuronal membrane which is separated from the postsynaptic neuron by the synaptic cleft. The action potential induces the release of neurotransmitters into the synaptic cleft which pass the information to the postsynaptic neuron and create the postsynaptic potential. The action potential can be described for example by the Hodgkin-Huxley model [START_REF] Hodgkin | A quantitative description of membrane current and its application to conduction and excitation in nerve[END_REF] which represents the nervous system as an electrical circuit. It was the first quantitative model that described the activation of membrane channels and was originally used to calculate action potentials in the squid giant axon. Since then, different models have been developed [START_REF] Fitzhugh | Impulses and Physiological States in Theoretical Models of Nerve Membrane[END_REF][START_REF] Hindmarsh | A model of neuronal bursting using three coupled first order differential equations[END_REF], but the Hodgkin-Huxley formalism remains the basic reference for mathematical modeling of the interdependence between the membrane potential and the different ionic currents. Basically, the neurons are traversed by sodium (Na+) and potassium (K+) currents which travel across voltage gated ionic channels. Each channel has a conductance that can be expressed as a maximum conductance, respectively ḡNa and ḡK , multiplied by a factor between 0 and 1. The Hodgkin-Huxley model considers that the potassium channel is composed by four independent components, closed or open. Thus, if n represents the proportion of open potassium channels, the probability to have these four components open is n 4 . This leads to a differential equation

The synaptic current

dn dt = α n (V )(1 -n) -β n (V )n,
where the opening and closing rates α n and β n depend on the membrane potential V. Then, the amount of ions that cross the membrane is given by

I K (t) = n 4 ḡK (V (t) -E K (t)),
with

E
dt = α m (V )(1 -m) -β m (V )m, dh dt = α h (V )(1 -h) -β h (V )h, (2.18)
with opening and closing rates α m (resp. α h ) and β m (resp. β h ) depending again on the membrane potential V . The sodium current is then given by

I N a (t) = m 3 hḡ N a (t)(V (t) -E N a (t)).
Kirchhoff's current law stipulates that the sum of the different ionic currents and the capacitive current C dV dt equals the applied exterior current (or stimulus) I app ,

I app (t) = C dV dt (t) + I N a (t) + I K (t) + I L (t),
where C is the specific membrane capacitance and I L denotes the leak current which depends linearly on the membrane potential V . The Hodgkin-Huxley model is the resulting set of non-linear differential equations

                             -C dV dt = m 3 hḡ N a (V -E N a ) + n 4 ḡK (V -E K ) + ḡL (V -E L ) -I app dn dt = α n (V )(1 -n) -β n (V )n dm dt = α m (V )(1 -m) -β m (V )m dh dt = α h (V )(1 -h) -β h (V )h (2.19)
where the expressions for the rates α(V ) and β(V ) of the different gating particles are obtained by experimental fitting and can be found in literature [START_REF] Hodgkin | A quantitative description of membrane current and its application to conduction and excitation in nerve[END_REF]. The resolution of problem (2.19) gives the dynamics of the membrane potential V and allows the numerical simulation of action potentials (see Section 3.2).

In response to an action potential (i.e. a presynaptic spike), neurotransmitters are released at the presynaptic terminal. These neurotransmitters then move into the synaptic cleft and bind with receptors in the postsynaptic neuron, opening ion channels. This movement of ions across the neuronal membrane generates a postsynaptic current I syn which is given by

I syn (t) = g(t)(V post (t) -V rev ). (2.20)
Here, the conductance g(t) depends on the presynaptic neuron, V post is the potential accross the postsynaptic membrane, and V rev is a (constant) reversal potential the value of which depends on the neurotransmitter.

There are different types of neurotransmitters and the model will depend on which neurotransmitter we consider. Here, we focus on a single channel model which is valid e.g. for AMPA or GABA A neurotransmitters which are, respectively, excitatory or inhibitory neurotransmitters. As before, we express the conductance as a maximum conductance ḡ and a factor between 0 and 1 modeling the proportion of open channels. Thus, we have g(t) = ḡs(t) in (2.20) where s is obtained by solving the differential equation [START_REF] Destexhe | Synthesis of models for excitable membranes, synaptic transmission and neuromodulation using a common kinetic formalism[END_REF]

ds dt = K 1 [T ](1 -s) -K 2 s, ( 2.21)
with K 1 , K 2 two constants that depend on the neurotransmitter and [T ] the concentration of neurotransmitters in the synaptic cleft. Here, we use a stationary relationship between [T ] and V pre from [START_REF] Destexhe | Synthesis of models for excitable membranes, synaptic transmission and neuromodulation using a common kinetic formalism[END_REF] according to which [T ] is given by

[T ](V pre ) = T max 1 + exp(-(V pre -V T )/K p ) , (2.22)
where T max is the maximum concentration of transmitters in the synaptic cleft, V T the value at which the concentration is halfed, K p models the steepness and V pre = V is the presynaptic action potential, solution to (2.19). Finally, we solve the differential equation

-C dV post dt = ḡs(t)(V post -V rev ) (2.23)
to get the postsynaptic potential V post . The corresponding postsynaptic current is then given by the right hand side of (2.23). Whereas action potentials only occur during a few milliseconds (1-2ms), postsynaptic potentials persist for a much longer period (20-40ms) and are dominating. It is commonly admitted that EEG signals result from the sum of postsynaptic currents generated by a large amount (≈ 10 6 ) of synchronized pyramidal neurons with similar orientiations. The above analysis describes potentials and currents at the neuron level.

If the active area of the brain at time t is localized at a position S m (t), the dipolar current source in the time-dependent forward EEG problem (2.16) can thus be defined with help of the postsynaptic current I syn by

q m (t) • δ Sm(t) = 10 6 I syn (t)u m (t) • δ Sm(t) (2.24)
where u m (t) ∈ R 3 is a unit vector describing the orientation along the dendrites of the synchronized neurons at time t. The solution u of (2.16) then models the spatial distribution in the head of the postsynaptic potential generated by the neuronal current I syn at S m over the time interval (0, T ).

The forward problem in time-dependent EEG

Existence and uniqueness result

In this section, we address the resolution of the forward problem (2.16). Let us give first some notations and hypotheses. Mathematically, a head model can be described as follows.

Let Ω ⊂ R d be a bounded simply connected domain with regular boundary Γ := ∂Ω and consider a partition of Ω into P open subdomains (Ω p ) p=1,...,P , such that

Ω = P p=1
Ωp and Ω p ∩ Ω q = ∅ ∀p ̸ = q.

Subdomains Ω p describe the different tissues of the head. In the case of concentric subdomains as in Figure 2, we denote by Γ p the interface between the subdomains Ω p and Ω p+1 and assume that (Γ p ) p are closed regular surfaces. Let n p be the unit normal vector to Γ p from Ω p to Ω p+1 .

We make the following assumptions on the moments, sources, and the conductivity:

S m (t) ̸ = S k (t), ∀m ̸ = k.
(H 3 ) All sources are located in the same subdomain Ω p 0 , p 0 ∈ {1, . . . , P }, and there is a convex domain

D ⊂⊂ Ω p 0 , such that M m=1 {S m (t), t ∈ [0, T ]} ⊂ D, ( 3.1)
and the conductivity is constant on

D: σ |D = σ D for a constant σ D > 0.
Due to the lack of regularity of the source term F , a direct variational formulation of (2.16) in H 1 (Ω) is not possible. A possibility is to apply the subtraction approach [START_REF] Badia | An inverse source problem in potential analysis[END_REF][START_REF] Azizollahi | EEG in neonates: forward modeling and sensitivity analysis with respect to variations of the conductivity[END_REF][START_REF] Wolters | Numerical mathematics of the subtraction method for the modelling of a current dipole in EEG source reconstruction using finite element head models[END_REF]. It consists of decomposing the potential u into a potential ũ which contains the singularity and a regular lifting w:

u = ũ + w on (0, T ) × Ω.
But in the time-dependent setting with moments q m belonging to L 2 (0, T ) d , the boundary value problem (2.16) with source term (2.15) can not be understood at fixed time t. The definition of the singular portential ũ thus needs to be done with caution. To this end, we introduce for m ∈ {1, . . . , M } and i ∈ {1, . . . , d}, the canonical source term

F (i) m : ∀t ∈ [0, T ], F (i) m (t, •) = e (i) • ∇δ(• -S m (t)) (3.2)
where e (i) denotes the i-th canonical basis vector of R d . Notice that

F (i) m is well defined at any time step t since S m is continuous on [0, T ]. The associated singular potential ũ(i) m (t, •), i = 1, . . . , d, is solution of the following Poisson equation ∀t ∈ [0, T ], σ D ∆ũ (i) m (t, •) = F (i) m (t, •) in R d .
Thus, ũ(i) m (t, •) is obtained by convolution in the space variable x of the fundamental solution of the Laplace equation with the right hand side

1 σ D e (i) • ∇δ(• -S m (t)) which leads to ũ(i) m (t, x) = 1 2 d-1 πσ D e (i) • x -S m (t) |x -S m (t)| d , ∀t ∈ [0, T ], ∀x ∈ R d \ {S m (t)}. (3.3)
At fixed time t ∈ [0, T ], consider the following boundary value problem with regular right hand side

   -∇ • (σ∇w (i) m (t, •)) = ∇ • ((σ -σ D )∇ũ (i) m (t, •))
in Ω,

σ∂ n w (i) m (t, •) = -σ∂ n ũ(i) m (t, •) on Γ. (3.4)
We prove here below that (3.4) admits a variational solution in an appropriated vector space. To this end, define the bilinear form a(•,

•) on H 1 (Ω) × H 1 (Ω) by a(w, v) = Ω σ∇w • ∇v dx.
For fixed t ∈ [0, T ], denote by l (i) m (t; •) the following linear form defined for v ∈ H 1 (Ω),

l (i) m (t; v) = Ω\ D(σ D -σ)∇ũ (i) m (t, •) • ∇v dx - Γ σ D ∂ n ũ(i) m (t, •)v ds.
Since problem (3.4) involves a Neumann boundary condition, its solution is determined up to an additive constant only. We therefore introduce the subspace of H 1 (Ω) of functions with vanishing mean value,

V = H 1 (Ω) ∩ L 2 0 (Ω), (3.5)
where

L 2 0 (Ω) = v ∈ L 2 (Ω) Ω v dx = 0 .
Then, the following proposition holds true:

Proposition 3.1. Let σ ∈ L ∞ (Ω) be such that 0 < σ min ≤ σ ≤ σ max a.e.
on Ω with constants σ min and σ max . Let S m ∈ C 0 (0, T ) d for any m ∈ {1, . . . , M } and assume that (H 2 ) and (H 3 ) are satisfied. Then, for m ∈ {1, . . . , M } and i ∈ {1, . . . , d}, the following problem admits a unique solution:

   Find w (i) m ∈ C 0 ([0, T ]; V ) such that a(w (i) m (t, •), v) = l (i) m (t; v) ∀v ∈ H 1 (Ω), ∀t ∈ [0, T ]. (3.6)
The solution of (3.6) satisfies (3.4) in a weak sense.

Proof. According to the assumptions on σ, the bilinear form a(•, •) is continuous and coercive on V × V . For the linear form, consider v ∈ H 1 (Ω) and t ∈ [0, T ]. We have

|l (i) m (t; v)| ≤ 2∥σ∥ L ∞ (Ω) ∥∇ũ (i) m (t, •)∥ L 2 (Ω\ D) ∥v∥ H 1 (Ω) + C Γ σ D ∥∂ n ũ(i) m (t, •)∥ L 2 (Γ) ∥v∥ H 1 (Ω)
where C Γ > 0 is the constant of the trace operator on H 1 (Ω). The potential ũ(i) m (t, •) is regular outside the domain D which allows to conclude that l

(i) m (t; •) is continuous on H 1 (Ω) and satisfies |l (i) m (t; v)| ≤ C∥v∥ H 1 (Ω) with a constant C = 2∥σ∥ L ∞ (Ω) max t∈[0,T ] ∥∇ũ (i) m (t, •)∥ L 2 (Ω\ D) + C Γ σ D max t∈[0,T ] ∥∂ n ũ(i) m (t, •)∥ L 2 (Γ)
which is independent from t. The compatibility condition l (i) m (t; 1) = 0 can be proved as in [START_REF] Azizollahi | EEG in neonates: forward modeling and sensitivity analysis with respect to variations of the conductivity[END_REF] with the help of the solid angle formula. Then, Lax-Milgram's theorem guarantees the uniqueness of the solution at fixed time t. We further get

∥w (i) m (t, •)∥ V ≤ C α ∀t ∈ [0, T ]
where α > 0 is the coercivity constant of a(•, •). This proves that w

(i) m ∈ C 0 ([0, T ]; V ). ■
Now, the time-dependent EEG problem (2.16) with source term (2.15) admits a solution in the following sense: Theorem 3.2. Let σ ∈ L ∞ (Ω) be such that 0 < σ min ≤ σ ≤ σ max a.e. on Ω with constants σ min and σ max . Assume (H 1 ) -(H 3 ) to be true. For i ∈ {1, . . . , d} and m ∈ {1, . . . , M }, let w (i) m be the solution of (3.6) where ũ(i) m is defined by (3.3). Finally, denote by q (i)

m the i-th component of the moment q m ∈ L 2 (0, T ) d . Then, u = M m=1 d i=1 q (i) m (ũ (i) m + w (i) m ) (3.7)
is solution of the time-dependent EEG problem

∇ • (σ∇u) = F in (0, T ) × Ω, σ∂ n u = 0 on (0, T ) × Γ,
with source term

F (t, •) = M m=1 q m (t) • ∇δ(• -S m (t)).
Proof. First notice that ∇ • (σ∇(ũ

(i) m + w (i) m )
) is well defined at any time t since

∇ • (σ∇(ũ (i) m + w (i) m )) = ∇ • ((σ -σ D )∇ũ (i) m ) + σ D ∆ũ (i) m + ∇ • (σ∇w (i) m ) = F (i) m .
We further have σ∂ n (ũ

(i) m + w (i) m ) = 0 on [0, T ] × Γ by construction. Now, let u be defined by (3.7). Then, ∇ • (σ∇u) = M m=1 d i=1 q (i) m ∇ • (σ∇(ũ (i) m + w (i) m )) = M m=1 d i=1 q (i) m F (i) m = F in (0, T ) × Ω
and

σ∂ n u = M m=1 d i=1 q (i) m ∂ n (ũ (i) m + w (i) m ) = 0 on (0, T ) × Γ. ■

Numerical simulation of the potential at the electrodes

In experiments, a finite number of surface electrodes are attached to the boundary Γ. Consider L well-separated pointwise electrodes {e ℓ } L ℓ=1 on Γ. Measurements are given by

u(e ℓ ) = U ℓ , ℓ = 1, . . . , L, (3.8)
with U = (U ℓ ) L ℓ=1 ∈ R L . As we have seen in the study of the forward problem, due to the Neumann boundary condition, the solution is determined up to an additive constant only. In theory, we can consider solutions in the subspace of H 1 (Ω) with vanishing mean value to have the uniqueness of the solution. In practice, the absolute potential cannot be measured by EEG and the measured voltage is the difference between the potential at electrode e ℓ and a reference potential set to zero at a reference electrode e 0 . Numerically, this is an alternative to fix the constant of the Neumann problem which guarantees the uniqueness of the solution. We present numerical simulations of the postsynaptic potential U recorded at the electrodes by solving the forward problem (2.16) with source term F using the subtraction approach. According to the previous section, we propose the computation of the source term F through the following steps: firstly, solve the Hodgkin-Huxley equations numerically to obtain the dynamics of the membrane potential V and simulate action potentials (spikes), then use the potential V to calculate the concentration of neurotransmitters according to (2.22) (GABA_A in the examples hereafter with constants taken from [START_REF] Destexhe | Synthesis of models for excitable membranes, synaptic transmission and neuromodulation using a common kinetic formalism[END_REF]) and finally deduce the postsynaptic current I syn from (2.20). Then, F is given by (2.24) for chosen orientation and position. In the simulations hereafter, we choose as stimulus a constant applied current equal to 9mA over a time domain of 100 ms. The Hodgkin-Huxley equations are solved by a variable-time-step, variable-order solver based on backward differentiation formulas of orders 1 to 5 and yield the time-course of the membrane potential shown on Figure 3. This choice produces regular spiking. Each spike has a characteristic time of about 1-2 ms. With this presynaptic potential, we solve differential equations (2.21) and (2.23) to get the postsynaptic current shown on Figure 4. According to values in literature [START_REF] Sterratt | Principles of Computational Modelling in Neuroscience[END_REF], the maximal conductance ḡ in (2.20) is at most 1 nS, voltage is expressed in mV, such that the postsynaptic current I syn is of the order of 10 -12 A. Now, if we consider that about a million neurons are simultaneously active and oriented in the same direction, the resulting current is of order 10 -6 A which is compatible with the values taken in [START_REF] Azizollahi | EEG in neonates: forward modeling and sensitivity analysis with respect to variations of the conductivity[END_REF]. Depending on the distance between the source point and the surface of the scalp, the recorded potential has a magnitude of several mV as observed in experimental EEG patterns. At any discrete time step t k , we solve problem (2.16) by the subtraction approach and use the finite element method with Lagrange-P 1 elements implemented in FreeFem++ [START_REF] Hecht | New development in FreeFem++[END_REF] to compute the regular potentials w (i) m (t k , •), i ∈ {1, . . . , d}, solution to the variational formulation (3.6). Then, the potential is computed using (3.7). Since the singularity of the potential u is taken account exactly in the subtraction approach, the discretization error is given by the error on the regular potential. We refer to [START_REF] Azizollahi | EEG in neonates: forward modeling and sensitivity analysis with respect to variations of the conductivity[END_REF] for a detailed error analysis. One may notice that the integrals of the linear forms l (i) m are well defined since the integration domain does not contain the singular point S m (t). However, it has been mentioned in [START_REF] Wolters | Numerical mathematics of the subtraction method for the modelling of a current dipole in EEG source reconstruction using finite element head models[END_REF] that the discretization error at a fixed mesh size is of order O(1/δ) where δ represents the distance of S m (t k ) to the nearest interface Γ p . Numerical verification of the code showed convergence rates with respect to the mesh size of 1.2 and 1.98 in the H 1 -and L 2 -norm, respectively. Visualization confirmed that the discretization error is maximal near the projection of the source position S m onto the nearest interface.

Hereafter, we provide numerical simulations of the postsynaptic potential in a circular multi-layer configuration in 2D and on a realistic head model of a healthy fullterm newborn obtained from coregistration of MR and CT images of the Amiens' hospital database (courtesy GRAMFC, INSERM U1105, Amiens, France (H. Azizollahi [START_REF] Azizollahi | Effects of Uncertainty in Head Tissue Conductivity and Complexity on EEG Forward Modeling in Neonates[END_REF])). The models distinguish between five tissues: white (WM) and grey matter (GM), CSF, skull and scalp for the circular 2D model and brain, CSF, skull, fontanels and scalp for the realistic head model.

Results for the 2D circular configuration.

Figure 6 represents the voltage measured at eight electrodes in the 2D circular configuration. One clearly distinguishes the three peaks which arise in the postsynaptic current (see Figure 4). The position of the electrodes and the source point as well as the width of the different layers of the model are given in Figure 5 and Table 3. In the present setting, the dipole is directed towards electrode e 5 , and e 1 is the reference electrode where the potential u is set to zero. Then, the amplitude of the measured voltage is maximal at e 5 and negative since the vectors q and (e 5 -S) are in opposite directions. Table 3: Width of the different tissues from interior to exterior layers. 2D circular configuration.

Figure 6: Postsynaptic potential at 8 electrodes between 0 and 100 ms. Circular 2D configuration.

Results for the realistic head model.

In this paragraph, we show numerical simulations of the postsynaptic potential for the realistic head model. Figure 7 (left) visualizes the distribution of four of eight electrodes where we record the brain activity. The source is located in the grey matter, close to the interface with the CSF and directed towards the surface (Figure 7, middle). Figure 7 (right) gives a snapshot of the simulated potential on the scalp at t = 10ms. In Figure 8, we represent the voltage at eight electrodes on the surface of the realistic head model over the time interval (0, T ). The reference electrode is situated outside the set of the represented electrodes. One clearly sees that electrodes localized in the scalp area towards which the source is oriented, record the most significant information (red line corresponding to e 2 in Figure 8).

Mathematical analysis of the inverse source problem

The inverse source problem we are interested in, consists in reconstructing the timedependent source term F of the form (2.15) from a measurement g := u |(0,T )×Γ * where Γ * ⊂ Γ with meas(Γ * ) ̸ = 0. More precisely, assuming the number M of sources to be known, the aim is to determine the locations (S m (t)) m=1,...,M and the moments (q m (t)) m=1,...,M of the source F as functions of time. Notice that numerical methods exist that are able to retrieve the number of sources independently from the source locations and its moments [START_REF] Badia | An inverse source problem in potential analysis[END_REF][START_REF] Nara | A projective method for an inverse source problem of the Poisson equation[END_REF]. These methods apply at fixed time steps t ∈ [0, T ]. We define the following non-linear mapping Λ, called measurement operator,

Λ : L 2 (0, T ) d × C 0 (0, T ) d M -→ L 2 (0, T ; H 1 2 (Γ * )) φ = (q m , S m ) m=1,..,M -→ u |(0,T )×Γ * , (4.1)
where u is the solution of the direct problem

∇ • (σ∇u) = F in (0, T ) × Ω, σ∂ n u = 0 on (0, T ) × Γ, (4.2)
with time-dependent source term

F (t, •) = M m=1 q m (t) • ∇δ(• -S m (t)). (4.
3)

The inverse problem then reads as follows: Given a measurement g ∈ L 2 (0, T ; H

1 2 (Γ * )), retrieve a source configuration φ ∈ L 2 (0, T ) d × C 0 (0, T ) d M such that Λ(φ) = g. (4.4)
In this section, we prove two essential properties for this inverse source problem: stability and identifiability. Notice that we do not address in this study the question of existence of a solution since it practically never arises when the inverse problem corresponds to a realistic situation as it is the case in EEG source reconstruction. Indeed, lack of existence indicates rather a wrong formulation of the inverse problem. In §4.2, we give an identifiability result if both the moments and the source points in configuration φ are continuous in time. In this particular case, identifiability follows directly from the corresponding result for static source terms [START_REF] Darbas | An inverse dipole source problem in inhomogeneous medi: application to the EEG source localization in neonates[END_REF]. Stability is proven in §4.3, again for source configurations that are continuous in time. The essential point is here to prove that the stability constant is time-independent. The main result of this section is given in §4.4 where identifiability is proven for moments that have only L 2 -regularity.

Differentiability of the measurement operator

The aim of this section is to prove regularity results for the measurement operator Λ defined in (4.1). To this end, let φ = (q m , S m ) 1≤m≤M ∈ (L 2 (0, T ) d × C 0 (0, T ) d ) M be a time-dependent source configuration, and consider

ψ = (p m , T m ) 1≤m≤M ∈ (L 2 (0, T ) d × C 0 (0, T ) d ) M such that S m (t) + T m (t) ⊂⊂ D ∀t ∈ [0, T ], m ∈ {1, . . . , M }. (4.5)
Let us define by u 1 ψ the solution of the following boundary value problem

∇ • (σ∇u 1 ψ ) = F 1 ψ in (0, T ) × Ω, σ∂ n u 1 ψ = 0 on (0, T ) × Γ, (4.6)
with the source term

F 1 ψ (t, •) = M m=1 p m (t) • ∇δ(• -S m (t)) -q m (t) • H(δ(• -S m (t))T m (t) (4.7)
where H(•) denotes the Hessian matrix with respect to the space variable.

The following proposition holds.

Proposition 4.1. The mapping Λ (4.1) is Fréchet differentiable and its derivative at φ, denoted by Λ ′ (φ), is given by

Λ ′ (φ) : L 2 (0, T ) d × C 0 (0, T ) d M -→ L 2 (0, T ; H 1 2 (Γ * )) ψ -→ (u 1 ψ ) |(0,T )×Γ * , (4.8)
where u 1 ψ is the solution of (4.6). The proof of Proposition 4.1 is an adaptation of Proposition 4.1 in [START_REF] Darbas | An inverse dipole source problem in inhomogeneous medi: application to the EEG source localization in neonates[END_REF] to a timedependent setting. It relies essentially on the subtraction approach for instationnary source terms developped in §3.1. We sketch the main steps of the proof and refer to [START_REF] Darbas | An inverse dipole source problem in inhomogeneous medi: application to the EEG source localization in neonates[END_REF] for any details.

Proof. Let φ ∈ (L 2 (0, T ) d × C 0 (0, T ) d ) M be a time-dependent source configuration and denote by ψ = (p m , T m ) 1≤m≤M any increment in (L 2 (0, T ) d × C 0 (0, T ) d ) M such that (4.5) holds true. Let ∥•∥ L 2 ×C 0 denote the canonical norm of (L 2 (0, T ) d × C 0 (0, T ) d ) M .
In order to prove that lim

∥ψ∥ L 2 ×C 0 →0 ∥Λ(φ + ψ) -Λ(φ) -Λ ′ (φ)ψ∥ L 2 (0,T ;H 1 2 (Γ * ))
∥ψ∥ L 2 ×C 0 = 0, (4.9)

denote by F ψ and F the source terms of type (4.3) corresponding, respectively, to source configurations φ + ψ and φ. Let u ψ and u be the associated solutions. Source term F 1 ψ and solution u 1 ψ corresponding to the term Λ ′ (φ)ψ have been defined in (4.7) and (4.6), respectively. Both u ψ and u can be rigorously defined by the subtraction approach in §3.1. The subtraction approach also applies to problem (4.6). Indeed, this has been proven in [START_REF] Darbas | An inverse dipole source problem in inhomogeneous medi: application to the EEG source localization in neonates[END_REF] for a static configuration and can be generalized as in §3.1 to the time-dependent case by decomposition of the moments p m (•) and q m (•) on the canonical basis of R d . Hence, (4.9) reads lim

∥ψ∥ L 2 ×C 0 →0 ∥u ψ -u -u 1 ψ ∥ L 2 (0,T ;H 1 2 (Γ * )) ∥ψ∥ L 2 ×C 0 = 0. (4.10)
Now, the function e ψ := u ψ -u -u 1 ψ can be split into a singular part ẽψ and a variationnal correction term belonging to H 1 (Ω). The singular part is obtained as before by convolution of F ψ -F -F 1 ψ with the fundamental solution G of the Laplacian operator. Taylor expansion of the first order derivatives of G(• -S m -T m ) allows to show that

∥ẽ ψ ∥ L 2 (0,T ;H 1 2 (Γ * ))
is of order O(∥ψ∥ 2 L 2 ×C 0 ). The regular part of e ψ is solution to a variational problem which involves ẽψ in a linear way and is thus also of order O(∥ψ∥ 2 L 2 ×C 0 ). This completes the proof. ■

Identifiability result in C 0

For source configurations φ that are continuous in time, i.e. q m ∈ C 0 (0, T ) d and S m ∈ C 0 (0, T ) d (m = 1, . . . , M ), we can generalize several results that have been obtained in [START_REF] Darbas | An inverse dipole source problem in inhomogeneous medi: application to the EEG source localization in neonates[END_REF] in the static configuration. To this end, consider the following (static) measurement operator Λ s

Λ s : R d × D M -→ H 1 2 (Γ * ) (q m , S m ) 1≤m≤M -→ u s|Γ * (4.11)
where (q m , S m ) 1≤m≤M is a static source configuration, i.e. q m ∈ R d and S m ∈ D ⊂ R d , and the measurement u s|Γ * is obtained by solving the following forward problem

       -∇ • (σ∇u s ) = M m=1 q m • ∇δ Sm in Ω, σ∂ n u s = 0 on Γ. (4.12)
Notice that the classical subtraction approach for the stationary EEG problem applies to (4.12). The link with the time-dependent measurement operator Λ is the following: Assume that φ = (q m , S m ) 1≤m≤M with q m ∈ C 0 (0, T ) d and S m ∈ C 0 (0, T ) d is a regular time-dependent source configuration. Then,

Λ s (φ(t)) = Λ(φ)(t, •) in H 1/2 (Γ * ) ∀t ∈ [0, T ]. (4.13)
The identifiability result for source configurations with regular moments and source points now reads as follows Theorem 4.2. Assume that the boundary Γ and the interfaces (Γ p ) p are of class C 2 . Let σ ∈ L ∞ (Ω) be such that 0 < σ min ≤ σ ≤ σ max a.e. on Ω with constants σ min and σ max .

In addition, assume that σ |Ωp ∈ W 1,∞ (Ω p ) for any p = 1, . . . , P .

Let g ℓ , ℓ = 1, 2 be two measurements in L 2 (0, T ; H 1 2 (Γ * )) originating respectively from source configurations φ ℓ = (q ℓ m , S ℓ m ) 1≤m≤M . Assume that q ℓ m , S ℓ m belong to C 0 (0, T ) d for any m = 1, . . . , M and ℓ = 1, 2, and let (H 2 ) and (H 3 ) hold true for φ ℓ , ℓ = 1, 2.

Then, if

g 1 = g 2 on [0, T ] × Γ * , φ 1 (t) ≡ φ 2 (t) up to a permutation π t of the integers (1, . . . , M ) for all t ∈ [0, T ].
Theorem 4.2 generalizes the static identifiability result of [START_REF] Darbas | An inverse dipole source problem in inhomogeneous medi: application to the EEG source localization in neonates[END_REF] to time-dependent source terms with regular moments and source points. The proof relies on the unique continuation principle for elliptic equations. Its application requires that the regular part w of the solution of the direct problem (4.2) belongs to H 2 on each subdomain Ω p which is true under the regularity assumptions on the conductivity σ and the boundary (resp. interfaces) of Ω. Notice that the C 0 -regularity assumption on the moments and source points guarantees that the arguments of the proof in [START_REF] Darbas | An inverse dipole source problem in inhomogeneous medi: application to the EEG source localization in neonates[END_REF] can be applied at any fixed time t ∈ [0, T ].

Stability result in C 0

In order to formulate a stability result, we need to fix a distance between two source configurations φ ℓ = (q ℓ m , S ℓ m ) 1≤m≤M , ℓ = 1, 2. Assume that φ 1 and φ 2 have the same number of sources M and belong to C 0 (0, T ) d × C 0 (0, T ) d M . Let t ∈ [0, T ] be fixed.

Setting S ℓ (t) = (S ℓ m (t)) 1≤m≤M , ℓ = 1, 2, the Hausdorff distance between S 1 (t) and S 2 (t) is defined by

d H (S 1 (t), S 2 (t)) = max max 1≤i≤M min 1≤j≤M ∥S 1 i (t) -S 2 j (t)∥, max 1≤i≤M min 1≤j≤M ∥S 1 j (t) -S 2 i (t)∥
where ∥ • ∥ designates the Euclidean norm in R d . It has been shown in [START_REF] Badia | Stability estimates for an inverse source problem of Helmholtz's equation from single Cauchy data at a fixed frequency[END_REF] that there exists a permutation π t of the integers {1, . . . , M } such that

d H (S 1 (t), S 2 (t)) = max 1≤m≤M ∥S 1 m (t) -S 2 πt(m) (t)∥. (4.14)
Notice that in the time-dependent setting, the permutation π t may depend on t. We then define the distance between φ 1 (t) and φ 2 (t) by

d π (φ 1 (t), φ 2 (t)) = max 1≤m≤M ∥q 1 m (t) -q 2 πt(m) (t)∥ + max 1≤m≤M ∥S 1 m (t) -S 2 πt(m) (t)∥ (4.15)
where π t is the permutation such that (4.14) holds true. Finally, we introduce the following notation for the distance between configurations φ 1 and φ 2 .

∥φ 1 -φ 2 ∥ π := T 0 |d π (φ 1 (t), φ 2 (t))| 2 dt 1 2 . ( 4.16)
Next, fix two positive constants 0 < a 1 < a 2 and denote by Q the set of all vectors q ∈ R d such that a 1 ≤ ∥q∥ ≤ a 2 . The following stability result holds true.

Theorem 4.3. Let φ ℓ (ℓ = 1, 2) be two configurations belonging to C 0 (0, T ) d × C 0 (0, T ) d M
. Assume that (H2) and (H3) are satisfied and that

a 1 ≤ ∥q ℓ m (t)∥ ≤ a 2 , (ℓ = 1, 2), ∀t ∈ [0, T ], ∀m = 1, . . . , M. (4.17)
Then, there exists a constant C > 0, independent from t and φ ℓ , such that

∥φ 1 -φ 2 ∥ π ≤ C T 0 ∥Λ(φ 1 )(t, •) -Λ(φ 2 )(t, •)∥ 2 H 1/2 (Γ * ) dt. (4.18)
Proof. We first notice that according to assumption (H 3 ), the trajectories of source points {S ℓ (t), t ∈ [0, T ]} belong to the convex and compact set D ⊂⊂ Ω ℓ 0 which is at a distance α > 0 from ∂Ω ℓ 0 . At a fixed time step t ∈ [0, T ], we thus can apply the stability result for the operator Λ s from [START_REF] Darbas | An inverse dipole source problem in inhomogeneous medi: application to the EEG source localization in neonates[END_REF] which reads in the present context

d π (φ 1 (t), φ 2 (t)) ≤ c∥Λ s (φ 1 (t)) -Λ s (φ 2 (t))∥ H 1/2 (Γ * ) , ∀φ 1 (t), φ 2 (t) ∈ Q M × D M .
In the above estimate, the constant c > 0 is independent from t since the estimate holds true uniformly for all source configurations which belong to the domain Q M × D M . Therefore, c only depends on Q and D, which are independent from t, as well as on the static operator Λ s . The stability estimate (4.18) then follows from integration over t and taking into account identity (4.13). ■

Identifiability result in L 2

We prove now a uniqueness theorem for source configurations with moments in L 2 (0, T ) d which ensures that the inverse source problem is well posed in the following sense: if two measured potentials coincide on a non empty set Γ * during a time period (0,T), they originate from the same source.

Theorem 4.4. Assume that (Γ p ) p and Γ are of class C 2 for any p ∈ {1, . . . , P -1}. Let σ ∈ L ∞ (Ω) be such that 0 < σ min ≤ σ(x) ≤ σ max for almost any x ∈ Ω, where σ min and σ max are two given positive constants. In addition, assume σ p ∈ W 1,∞ (Ω p ), p = 1, . . . , P . Let F be a time-dependent source of the form (2.15) whose moments and locations satisfy hypotheses (H1), (H2) and (H3). Let Γ * be a subset of Γ such that meas(Γ * ) ̸ = 0, and let u be the solution of problem (2.16). Then the locations (S m ) 1≤m≤M and the moments (q m ) 1≤m≤M are uniquely determined up to a permutation of the integers {1, . . . , M } from a single measurement g =: u | (0,T )×Γ * .

Proof. Let

φ ℓ = (q ℓ m , S ℓ m ) 1≤m≤M ∈ (L 2 (0, T ; R d ) × C 0 (0, T ; R d )) M , ℓ = 1, 2,
be two source configurations corresponding to respective sources F ℓ , ℓ = 1, 2, of the form (2.15) such that Λ(φ 1 ) = Λ(φ 2 ) in L 2 (0, T ; H We shall prove that

∥φ 1 -φ 2 ∥ π = 0, (4.20)
with ∥ • ∥ π defined by (4.16). The density of C 0 (0, T ) in L 2 (0, T ) implies the existence of sequences

φ ℓ n := (q ℓ m,n , S ℓ m,n ) ∈ (C 0 (0, T ) × C 0 (0, T )) M , ℓ = 1, 2, satisfying lim n→∞ ∥φ ℓ n -φ ℓ ∥ π = 0, ℓ = 1, 2. (4.21)
We have

∥φ 1 -φ 2 ∥ π ≤ ∥φ 1 -φ 1 n ∥ π + ∥φ 1 n -φ 2 n ∥ π + ∥φ 2 n -φ 2 ∥ π . (4.22)
Since the first and the third term tend to zero according to (4.21), it is sufficient to prove

lim n→∞ ∥φ 1 n -φ 2 n ∥ π = 0 (4.23)
for obtaining (4.20). The stability result for source configurations with regular moments (see Theorem 4.3) gives

∥φ 1 n -φ 2 n ∥ 2 π ≤ C T 0 ∥Λ(φ 1 n )(t, •) -Λ(φ 2 n )(t, •)∥ 2 H 1 2 (Γ * ) dt ≤ 2C T 0 ∥Λ(φ 1 n )(t, •) -Λ(φ 1 )(t, •)∥ 2 H 1 2 (Γ * ) dt +2C T 0 ∥Λ(φ 1 )(t, •) -Λ(φ 2 (t))∥ 2 H 1 2 (Γ * ) dt +2C T 0 ∥Λ(φ 2 )(t, •) -Λ(φ 2 n )(t, •)∥ 2 H 1 2 (Γ * )
dt.

(4.24)

From (4.19) we get

T 0 ∥Λ(φ 1 )(t, •) -Λ(φ 2 )(t, •)∥ 2 H 1 2 (Γ * ) dt = 0.
According to Proposition 4.1, the operator Λ is Fréchet differentiable and thus continuous.

Together with the continuity of the norm, we obtain

lim n→∞ ∥Λ(φ ℓ n ) -Λ(φ ℓ )∥ 2 H 1 2 (Γ * ) = 0, in L 2 (0, T ), ℓ = 1, 2, (4.25)
which implies lim n→∞ ∥φ 1 n -φ n ∥ π = 0 by (4.24). From (4.22), we then deduce

∥φ 1 -φ 2 ∥ π = 0. ( 4.26)
■

Numerical reconstruction

In this section, we illustrate source reconstruction for moving source points and timedependent moments. A broad number of performing inversion methods have been developed in electric source imaging at fixed time steps: among others, we cite LORETA, sLORETA [START_REF] Pascual-Marqui | Standardized low resolution brain electromagentic tomography (sLORETA): technical details[END_REF], MUSIC [START_REF] Mosher | Recursive MUSIC: a framework for EEG and MEG source localization[END_REF], RAMUS [START_REF] Rezaei | Reconstructing subcortical and cortical somatosensory activty via the RAMUS inverse source analysis technique using median nerve SEP data[END_REF] as well as [START_REF] Pascarella | An invivo validation of ESI methods with focal sources[END_REF] for a comparative study. Here, we use MNE at given time steps in order to illustrate the concept of reconstructing timedependent moments and moving source points in the 2D circular configuration, since it is one of the most simplest methods which yields satisfying results when minimal a priori information about the sources is available [START_REF] Hämäläinen | Interpreting magnetic fields of the brain: minimum norm estimates[END_REF].

Lead Field Matrix and Minimum Norm Estimate

We now make precise the numerical method that will be used to solve the time-dependent inverse source problem 4.4. Let (t k ) 1≤k≤K be a sequence of K time steps belonging to the interval [0, T ]. Let φ = (q m , S m ) 1≤m≤M be a time-dependent source configuration defined on [0, T ] such that (q m (t k ), S m (t k )) is well defined for any m, 1 ≤ m ≤ M and any k, 1 ≤ k ≤ K. Let e 1 , . . . , e L ∈ Γ denote the positions of L electrodes situated at the surface of the scalp. For the potential u, solution to the forward EEG problem (2.16) with source configuration φ, U k EEG = ( u(t k , e 1 ) . . . u(t k , e L ) ) t ∈ R L denotes the measurement vector at time step t k . Notice that we have L + 1 electrodes, taking into account the reference electrode situated at e 0 , but only L measurements at a fixed time step.

The resolution of the inverse problem is based on a Lead Field Matrix approach [17, 31] which we will briefly explain now. Consider a source space S = {S 1 , . . . , S R } ⊂ D where D is the domain of hypothesis (H3) in §3. In the present study, we will take for S the interior nodes of Ω p 0 . For any source point S m ∈ S and any coordinate 1 ≤ i ≤ d, we recall that the canonical source term F (i) m is defined by

F (i) m = e (i) • ∇δ(• -S m )
where e (i) denotes the i-th vector of the canonical basis of R d . Notice that in this approach, the source points S m do not depend on t. Let u (i) m be the associated potential, solution of the static EEG problem (4.12) with source term F (i) m . Then, the measurement vector U k EEG is given as the following matrix-vector product

U k EEG =          u 1 (e 1 ) u 2 (e 1 ) • • • u R (e 1 ) u 1 (e 2 ) u 2 (e 2 ) • • • u R (e 2 ) . . . . . . . . . . . . u 1 (e L ) u 2 (e L ) • • • u R (e L )                   q 1 (t k ) q 2 (t k ) . . . q R (t k )          , ( 5.1)
where u m (e s ) = u

(i) m (e s ) 1≤i≤d and q m (t k ) = q (i) m (t k ) 1≤i≤d
are vectors in R d . In this approach, a given source configuration is represented by the amplitudes and directions of the moments (q m ) 1≤m≤R at all source points of the source space. For example, a single source localized at a point S at time t k corresponds to a vector q k with non-zero entries only for nodes in S which are located near to S, the amplitudes of the source being distributed over all nodes in S. The matrix in (5.1) is called Lead Field Matrix and denoted by L. One may notice that L is of size L × (dR) where d = 2 or d = 3 is the dimension of the computational domain and requires the resolution of dR forward problems. Since the source space S contains in general a huge number of sources, the computational cost of the Lead Field Matrix can be prohibitive unless its implementation is optimized. We perform optimization by adapting the principle from [START_REF] Wolters | Efficient computation of lead field bases and influence matrix for the FEM-based EEG and MEG inverse problem[END_REF] to the subtraction approach. Recall that the subtraction method is based on the splitting of the discrete potential u h into a singular potential ũ which is known explictly and a regular (discrete) potential w h which is discretized by finite elements. Let N be the number of nodes of the whole finite element mesh and denote by W k ∈ R N the nodal vector associated with w h (t k , •). W k is the solution of a linear system of the form

KW k = J k ,
where the right hand side J k depends on the singular potential ũ. Since we only need the values of w h at the electrodes, we apply a restriction matrix R ∈ M L,N and define

W k EEG = RW k = RK -1 J k .
Here, R ℓj = 1 if and only if the j-th node of the mesh coincides with the ℓ-th electrode. Then, in order to compute the matrix B = RK -1 ∈ M L,N advantageously, notice that

B = RK -1 ⇔ B t = K -t R t ⇔ K t B t = R t ⇔ KB t = R t
owing to the symmetry of the matrix K. Consequently, the ℓ-th column of B t is the solution of a linear system with matrix K and right hand side R t •,ℓ where R t •,ℓ denotes the ℓ-th column of R t . This algorithm only needs the resolution of L linear systems K, L ≪ (dR). Moreover, all systems have the same matrix K.

In order to retrieve the lead field matrix from matrix B, it remains to take into account the linear dependence of ũ and w h on the moment vector q dR) . Indeed, according to (3.7) we have

k = q 1 (t k ) • • • q R (t k ) t ∈ R (
ũ(t k , •) = M m=1 d i=1 q (i) m (t k )ũ (i) m outside S and w(t k , •) = M m=1 d i=1 q (i) m (t k )w (i) m ,
where ũ(i) m and w (i) m are independent from the moments q m . At the finite element level, we thus have

U k = R( Ũ k + W k ) = B(K Ũ k + KW k ) = B(K Ũ k + J k ).
In order to be well defined, the coefficients of Ũ k ∈ R N satisfy Ũ k j = 0 if C j is a node belonging to the source space S. Notice that J k = KW k depends linearly on Ũ k since the right hand side of the variational formulation (3.6) involves integrals of ũ(i) m only on Ω \ D and on Γ. This yields the factorization of the lead field matrix into

L = B(K + K)A (5.2) where K ∈ M N (R) is a matrix such that J k = K Ũ k and A ∈ M N,(dR) satisfies Ũ k = Aq k .
The computational cost of the Lead Field Matrix from relation (5.2) is higher than for the source models presented in [START_REF] Wolters | Efficient computation of lead field bases and influence matrix for the FEM-based EEG and MEG inverse problem[END_REF]. It can however be optimized by taking into account the sparse structure of the matrices K and K and the explicit knowledge of the coefficients of A in terms of the mesh nodes. Notice also that the Lead Field Matrix will only be computed once for a given head model and source space. Now, for 1 ≤ k ≤ K, the measurement vector U k EEG is related to the moment vector

q k by U k EEG = Lq k + η k , ( 5.3)
where η k ∈ R (dR) represents additive white noise. The Minimum Norm Estimate then corresponds to the minimizer of the following least square problem,

             For 1 ≤ k ≤ K, find q k, * ∈ R (dR) such that q k, * = arg min q∈R (dR) ∥Lq -U k EEG ∥ 2 + α∥q∥ 2 .
(5.4)

Here, α > 0 is the regularization parameter of classical Tikhonov regularization.

Moving source points

In this section, we present the numerical reconstruction of the trajectory and amplitude of a moving point source in the 2D circular configuration. The synthetic measurements at the electrodes are generated by solving the forward problem for a given time-depending source configuration by the subtraction approach as explained in §3. We perturb these measurements by a Gaussian noise. When we solve the inverse problem with the MNE, we obtain an estimation of the moment q m at every node in the source space. In order to determine the position of the estimated source, we decide to consider the source with the highest amplitude. Thus, the position of a source at a time t k is given by

S k, * = arg max m=1,...,R |q k, * m |, (5.5)
where q k, * denotes the solution of (5.4). We then compute the error on [0, T ] in the discrete ℓ 2 -norm by

E S = 1 K K k=1 ∥S k, * -S k ∥ 2 .
(5.6) We consider a single source which moves on a trajectory of the form S(t) = (r cos(tt 0 ), r sin(t-t 0 )) in the grey matter layer, with a moment normal to the interface. Moreover, we fix 32 electrodes to reconstruct the moving source (see Figure 9) and we perturb these measures with 5% noise. The dimension of the source space S is R = 1028. It is important to note that the exact source trajectory is parametrized in the subdomain Ω p 0 , but the reconstruction method just allows to find sources in the source space, i.e. on a node. We thus cannot expect to find the exact position. As we can see on the Figure 9, the estimation of the position is close to the real trajectory of the moving source. In fact, the mean error between both lines is 0.24 cm, which is in agreement with the spatial resolution of the order of 1 cm of general distributed methods. An important point which is crucial in the resolution of inverse problems, is the choice of an adequate regularization parameter. There are various ways to find a suitable parameter. Here, the L-curve method (see e.g. [START_REF] Hansen | the use of the L-curve in the regularization of discrete ill-posed problems[END_REF]) yields α ≈ 10 5 . This value is in good agreement with the optimal regularization value which minimizes the error to the exact source position.

As we can see on Figure 10, the estimated moment has the same shape as the exact moment but there is a factor 100 between the amplitudes. The magnitude of the moments indeed depends on the regularization parameter α. Varying α allows to improve the approximation of the moments' magnitude at the price of increasing the error in position. Further study of parameter selection methods, based e.g. on Bayesian Hierarchical Modeling as in [START_REF] Rezaei | Reconstructing subcortical and cortical somatosensory activty via the RAMUS inverse source analysis technique using median nerve SEP data[END_REF], should be able to improve the reconstruction.

Conclusion

In this paper, we formulate the EEG problem regarding evolution in time and give a rigorous mathematical setting for both the direct and the inverse source problem. Sources of electrical brain activity are modeled by time-dependent moments and source positions. The subtraction approach which has been studied in detail for the static case in [START_REF] Azizollahi | EEG in neonates: forward modeling and sensitivity analysis with respect to variations of the conductivity[END_REF], is generalized to moments with low (L 2 ) regularity in time and moving source points, and allows discretization in space by Lagrange Finite Elements with optimal convergence order. The validity of the quasi-stationary approximation of Maxwell's equations which is obtained by neglecting all time-derivatives despite a time-dependent source term, is verified for all tissues of the brain with the help of dimensional analysis. A full numerical pipeline is presented that simulates the electric potential at the electrodes over a given time interval. The electric potential is generated by a time-dependent postsynaptic current at the neuron level which, in turn, results from the resolution of a non-linear system of differential equations. The latter models the electrical and chemical processes at the synapses. Numerical results in a circular 2D configuration as well as on the realistic head model of a full-term neonate illustrate the method for the direct problem. For a fixed and time-independent number of sources, the inverse problem is analyzed mathematically by studying the properties of a non-linear measurement operator. We prove identifiability for moments with L 2 -regularity and continuous moving sources. Numerically, the inverse problem is formulated as a linear least square problem that involves a Lead Field Matrix. The optimized computation of the latter matrix is explained in the case where the subtraction approach is used to solve the forward problem. As a proof of concept, we end with the numerical reconstruction of a moving source point with time-dependent moment in a 2D academic configuration. Extensions to three dimensional cases are currently underway.

In clinical examination, EEG monitoring is by nature time-dependent. Numerical source reconstruction is in general done with data extracted from the EEG pattern at a specific time that has been considered as significative with regard to the brain activity under study. Spatial-temporal source analysis takes advantage of EEG measurements over a given time interval, assuming the source point stationary. To the best of the authors' knowledge, no time-dependent mathematical model for EEG has been written or analyzed up to now. Another contribution of the present paper is to generate time-dependent dipolar source terms by solving classical neuronal models. Concerning the mathematical study of the forward EEG problem, we choose to adapt the subtraction approach since it yields a rigorous way to define its solution. It further allows to analyze precisely the discretization error in terms of the mesh size when Lagrange Finite Elements are used. For both the forward and the inverse problem, we pay special attention to the case of moments with L 2 -regularity in time which do not allow to write the problem at a fixed time step.

Our time-dependent EEG model is the first step towards a mathematical description of neurovascular coupling for neonates and premature babies. Such models based on a large set of ordinary differential equations exist at the cellular level [START_REF] Aubert | A Model of the Coupling between Brain Electrical Activity, Metabolism, and Hemodynamics: Application to the Interpretation of Functional Neuroimaging[END_REF]. We actually work on linking the time-dependent models of the two involved imaging modalities, EEG and optical imaging [START_REF] Arridge | Optical tomography: forward and inverse problems[END_REF][START_REF] Lohrengel | A homogenized cerebrospinal fluid model for diffuse optical tomography in the neonatal head[END_REF]. This requires a mathematical description of the relationship between neuronal activities and the corresponding hemodynamic responses, which are localized in space and time. The modeling of temporal and spatial correlations is a very challenging task. The resulting complete model will be useful for the numerical resolution of the inverse problem for coregistered EEG-NIRS, namely the reconstruction of both brain electrical sources and optical parameters.

Figure 1 :

 1 Figure 1: Schematic of a synapse

  K the resting potential of the potassium ion. The case of the sodium ion current is a little different. In fact, the channels can be active and open, inactive and open, or closed, and two gating particles, m and h, are needed, representing respectively the fraction of active open and inactive open channels. The three possible states lead to a couple of differential equations dm

Figure 2 :

 2 Figure 2: Three-layer head model

Figure 3 :

 3 Figure 3: Dynamics of the membrane potential V and simulation of action potentials (spikes).

Figure 4 :

 4 Figure 4: Postsynaptic current I syn .

Figure 5 :

 5 Figure 5: Electrode positions and source point.

Figure 7 :

 7 Figure 7: Electrode positions (left). Frontal cut and characteristics of the source (middle). Potential at the scalp at time t = 10ms (right). Realistic head model.

Figure 8 :

 8 Figure 8: Postsynaptic potential at 8 electrodes between 0 and 100 ms. Realistic head model.

Figure 9 :

 9 Figure 9: Distribution of the 32 electrodes (left), Estimated (blue line) vs. exact trajectory (red line) (right).

Figure 10 :

 10 Figure 10: Exact moment: x component (top left), y component (top right). Estimated moment: x component (bottom left), y component (bottom right). Ordinate's unit is [A].

  .8) 2. τ e ≪ τ em : We get the order τ e ≪ τ em ≪ τ m ≪ τ.

	As before, Ampere's law (2.7) follows from τ e ≪ τ and	τ em τ e	e cb	∼ 1 implies
	τ em τ	cb e	∼	τ em τ	τ em τ e	=	τ m τ	≪ 1.
	This yields again (2.8).
	3. τ e ∼ τ em : In this case, the three characteristic times are of the same order and we
	have
		τ m ∼ τ em ∼ τ e ≪ τ.
	Ampere's law (2.7) follows as above from τ e ≪ τ where the quantities e and cb are
	now of comparable order: e ∼ cb. Consequently,
		τ em τ	cb e	∼	τ em τ	≪ 1
	and (2.8) follows.
	Summing up, in a frequency range satisfying (2.6), we get the following approximation
	of Maxwell's equations

Table 1 :

 1 Electric parameter set at 100 Hz[23,[START_REF] Lew | Effects of sutures and fontanels on MEG and EEG source analysis in a realistic infant head model[END_REF] and characteristic times.

		White matter Grey matter	CSF	Skull	Scalp
	ε r	1.67 × 10 6	3.91 × 10 6	1.1 × 10 2	5.85 × 10 3	1.14 × 10 3
	σ [S/m]	0.14	0.33	1.8	0.04	0.33
	τ m [s]	2.53 × 10 -9	5.97 × 10 -9 3.25 × 10 -8 7.2 × 10 -10 6.0 × 10 -9
	τ em [s]	5 × 10 -7	8 × 10 -7	4.19 × 10 -9 3.1 × 10 -8 1.35 × 10 -8
	τ e [s]	1.05 × 10 -5	1.05 × 10 -4 5.4 × 10 -10 1.3 × 10 -6	3.1 × 10 -8

Table 2 :

 2 Table 2 summarizes the mesh parameters. The conductivity values of the different tissues are taken from Table 1. For the realistic head model, brain conductivity is set to 0.33 [S/m] and fontanels are considered with σ = 0.3 [S/m]. Model Nodes Tetrahedra h min [m] h max [m] diameter [m] Characteristics of the 2D circular and the realistic head model.

	M 2D	8232	16184	8.9 10 -4 2.2 10 -3	0.12
	M real 108 669	590 878	3.4 10 -4 14 10 -3	0.12

Acknowledgments

We are grateful to H. Azizollahi and GRAMFC INSERM UMR-S 1105 for providing us with the mesh of the realistic head model of a neonate.