Grégory Dumont
email: gregory.dumont@ens.fr

Jacques Henry
email: jacques.henry@inria.fr

Carmen Oana Tarniceriu
email: tarniceriuoana@yahoo.co.uk

Theoretical connections between mathematical neuronal models corresponding to different expressions of noise

Keywords: Neural noise, Noisy Leaky Integrate-and-Fire model, Escape rate, Fokker-Planck equation, Age structured model

Identifying the right tools to express the stochastic aspects of neural activity has proven to be one of the biggest challenges in computational neuroscience. Even if there is no definitive answer to this issue, the most common procedure to express this randomness is the use of stochastic models. In accordance with the origin of variability, the sources of randomness are classified as intrinsic or extrinsic and give rise to distinct mathematical frameworks to track down the dynamics of the cell. While the external variability is generally treated by the use of a Wiener process in models such as the Integrate-and-Fire model, the internal variability is mostly expressed via a random firing process. In this paper, we investigate how those distinct expressions of variability can be related. To do so, we examine the probability density functions to the corresponding stochastic models and investigate in what way they can be mapped one to another via integral transforms. Our theoretical findings offer a new insight view into the particular categories of variability and it confirms that, despite their contrasting nature, the mathematical formalization of internal and external variability are strikingly similar.

Introduction

The presence of variability in the neural activity is well documented nowadays [START_REF] Aa Faisal | Noise in the nervous system[END_REF]. In vivo as well as in vitro experiments stand for the evidence of irregular behavior of neuronal activity. For instance, spike trains of individual cortical neurons in vivo are highly irregular [START_REF] Shadlen | Noise, neural codes and cortical organization[END_REF][START_REF] Softky | The highly irregular firing of cortical cells is inconsistent with temporal integration of random epsps[END_REF], and, apart from randomness observed in spontaneous neuronal activity, recordings of in vitro experiments for input stimuli without any temporal structure exhibited irregular behavior of neural activity [START_REF] Gerstner | Spiking neuron models[END_REF]. It is commonly accepted now that the random influences over the neuronal firing activity to be designated as noise [START_REF] Longtin | Neuronal noise[END_REF]. Among the most notable sources of noise, few are usually reminded : thermal noise [START_REF] Manwany | Detecting and estimating signals in noisy cable sstructure. i: Neuronal noise sources[END_REF], the effect of signal transmission in a network [START_REF] Manwany | Detecting and estimating signals in noisy cable sstructure. i: Neuronal noise sources[END_REF], randomness of excitatory and inhibitory connections [START_REF] Brunel | Fast global oscillations in networks of integrateand-fire neurons with low firing rates[END_REF], global network effects or the finite number of ionic channels on a neuronal patch [START_REF] White | Channel noise in neurons[END_REF].

According to the location where the noise is generated, it has become a common procedure to classify these sources of noise as extrinsic or intrinsic [START_REF] Gerstner | Spiking neuron models[END_REF]. While intrinsic noise usually refers to random features underlying the firing process, therefore generated at cell level, the extrinsic noise is generally attributed to the network effects and signal transmission over the cell.

An important step toward the recent development in theoretical neuroscience consists in a deep understanding of the essence of this variability. However, its mathematical formalization is still an open problem; it has been a subject of intense research, (see for instance [START_REF] Gerstner | Spiking neuron models[END_REF] for a discussion on this issue) and many recent papers are trying to suitably mathematically model these effects. A typical way to mathematically model random processes is by the use of stochastic differential equations. Nevertheless, the particular form of the noise terms and their incorporation into stochastic neuron models are still subject of controversy. Among others, two approaches gain more visibility in the last decade; each of them corresponds to different treatment of noise.

The external influences over the transmembrane potential of a single neuron is usually modeled in the form of a Gaussian white noise and gives rise to an Ornstein Uhlenbeck process. In the probability theory, the representation of a random variable can be done by the use of the so called probability density function (pdf), which describes the likelihood for the stochastic process to take on a specific value. Since a stochastic equation of Langevin type can be translated into a Fokker-Planck (FP) equation, this allows the representation of the external noise category as a diffusion process [START_REF] Gardiner | Handbook of Stochastic Method for Physics[END_REF]. In particular, the FP equation [START_REF] Gardiner | Handbook of Stochastic Method for Physics[END_REF] for an OU process belongs to the most prominent models in the literature.

It should be stressed that the use of pdf concept in the field of mathematical neuroscience has already a long history, as it can be seen in [START_REF] Wilbur | A theoretical basis for large coefficient of variation and[END_REF][START_REF] Abbott | Asynchronous states in networks of pulse-coupled oscillators[END_REF], and it has lead to revealing new insights into phenomena related to neuronal behaviors. The resulted mathematical formalism is in particular pertinent for the simulation of large sparsely connected populations of neurons in terms of a population density function [START_REF] Omurtag | On the simulation of large population of neurons[END_REF][START_REF] Nykamp | A population density appraoch that facilitates large-scale modeling of neural networks : analysis and an application to orientation tuning[END_REF][START_REF] Dumont | A density model for a population of theta neurons[END_REF]. The formulation of the dynamics of neural networks in terms of population densities made possible mathematical descriptions of phenomena that emerge at population level, such as oscillations and synchrony caused by recurrent excitation [START_REF] Dumont | Population density models of integrate-and-fire neurons with jumps, well-posedness[END_REF][START_REF] Dumont | Synchronization of an excitatory integrate-and-fire neural network[END_REF][START_REF] María | Analysis of nonlinear noisy integrate & fire neuron models: blow-up and steady states[END_REF], by delayed inhibition feedback [START_REF] Brunel | Fast global oscillations in networks of integrateand-fire neurons with low firing rates[END_REF], by both recurrent excitation and inhibition [START_REF] Brunel | Dynamics of sparsely connected networks of excitatory and inhibitory spiking neurons[END_REF], and by gap junction [START_REF] Ostojic | Synchronization properties of networks of electrically coupled neurons in the presence of noise and heterogeneities[END_REF], the emergence of neural cascade [START_REF] Newhall | Cascadeinduced synchrony in stochastically-driven neuronal networks[END_REF][START_REF] Newhall | Dynamics of current-based, poisson driven, integrate-and-fire neuronal networks[END_REF], and of self criticality with synaptic adaptation [START_REF] Daniel Millman | Self-organized criticality occurs in non-conservative neuronal networks during 'up' states[END_REF]. For connections between models corresponding to probability respectively population densities, we refer to [START_REF] Nykamp | A population density appraoch that facilitates large-scale modeling of neural networks : analysis and an application to orientation tuning[END_REF][START_REF] Knight | Dynamical models of interacting neuron populations in visual cortex[END_REF][START_REF] Knight | Dynamics of encoding in neuron populations: Some general mathematical features[END_REF][START_REF] Brunel | Fast global oscillations in networks of integrateand-fire neurons with low firing rates[END_REF].

To account for the intrinsic variability, neuronal models incorporating a stochastic mechanism of firing have been considered. In this framework, the state of the neuron follows a deterministic trajectory and each spike occurs with a certain probability given in the form of an escape rate or stochastic intensity of firing. This assumption lead to the introduction of models where only the time passed since the last spike influences the dynamics. The associated pdf to such a model has the form of a so-called age-structured (AS) system [START_REF] Gerstner | Spiking neuron models[END_REF], [START_REF] Khashayar Pakdaman | Relaxation and selfsustained oscillations in the time elapsed neuron network model[END_REF]. Such a process is a special form of renewal processes, which are a particular category of point processes having the particularity that the times between two successive events are independent and identically distributed. The main assumption upon which the process is built is that the whole history that happened before the last event (firing time) is forgotten.

It is therefore suitable to consider such a process for the case where neurons do not show strong adaptation.

Although the two above reminded models are usually thought to express different treatments of noise, it has been shown though that they do generate similar statistics of spike trains [START_REF] Gerstner | Spiking neuron models[END_REF]. An approximation method has been presented in [START_REF] Plesser | Noise in integrate-and-fire neurons: from stochastic input to escape rates[END_REF] to explain this similarity.

However, no exact relation between the solutions of these models has been yet proven. We did make a first step in this direction in [START_REF] Dumont | Noisy threshold in neuronal models: connections with the noisy leaky integrate-and-fire model[END_REF] by mapping the solution to the AS system into the solution to the FP equation. The relation that we proved is only partial since the reverse mapping of the solutions still lacked. To completely tackle this problem, we have investigated therefore the possibility of giving an analytical transform of the AS system into the FP equation. Our theoretical findings that we present here highlight an unforeseen relationship between the FP equation and the AS (of von-Foester-McKendrick type) system, which allows, in particular, to transfer qualitative properties of the solutions one-to-another. This explains not only the similarities of the observed statistical activity of both models but also help us rise a unitary vision over the different formalisms used to describe different noise expressions acting over a neural cell.

In this paper we do not intend to discuss the issue of which sources of noise should be taken into account when modeling the neuronal dynamics [START_REF] Longtin | Neuronal noise[END_REF], nor to debate the appropriateness of models used to express different variability. Instead, we aim to highlight a case where two different implementation of noise can be directly linked by a mathematical relation.

The paper is structured as follows: We will start by reminding the two approaches of noise implementation and the main characteristics of the corresponding models in the first section. As the kernels of our integral transforms depend on the solutions to the forward and backward Chapman-Kolmogorov (CK) equation for a probability density function on an inter-spike interval, we will roughly remind in the second section the derivation and the meaning of these two systems. Our main results are given in the third section, and we end this paper by some conclusive remarks and discussing possible extensions of the present work. Finally, all along the paper, we present numerical simulations to illustrate the models presented in it.

Two standard approaches in expressing neuronal variability

Due to the influence of noise, the evolution in time of the state of a neuron is described by a suitably chosen stochastic process. In computational neuroscience, a first step toward the description of neural variability was made in [START_REF] Stein | Some models of neuronal variability[END_REF] and [START_REF] Gerstein | Random walk models for the spike activity of a single neuron[END_REF]. In this section, we remind the reader two standard mathematical formalizations of neural variability in use nowadays.

First we introduce the model that incorporates an extrinsic source of noise in the form of the noisy leaky integrate-and-fire model, and then present its associated FP equation. Next, we turn to a model that expresses the intrinsic noise via a stochastic firing intensity which will lead to characterizing the evolution of the state of a neuron in the form of an AS system.

Extrinsic noise. Fokker-Planck formalism

As a common procedure to handle the extrinsic variability, noise terms are explicitly added to the mathematical models that describe the evolution of the state of the nervous cell. In this way, the evolution of the state of a neuron is viewed as a random process which is described by a stochastic differential equation (SDE). Throughout this paper, we will illustrate all the theoretical considerations related to the extrinsic variability for the specific case of the noisy leaky integrate-and-fire (NLIF) model [START_REF] Izhikevich | Dynamical Systems in Neuroscience[END_REF]. The integrate-and-fire model describes a point process (see Fig. 1) and is largely used because it combines a relative mathematical simplicity with the capacity to realistically reproduce observed neuronal behaviors [START_REF] Gerstner | How good are neuron models[END_REF]. The model idealizes the neuron as a simple electrical circuit consisting in a capacitor in parallel with a resistor driven by a given current and was first introduced by Lapique in 1907, see [START_REF] Brunel | Lapicque's 1907 paper: from frogs to integrate-and-fire[END_REF][START_REF] Abbott | Lapique's introduction of the integrate-and-fire model neuron (1907)[END_REF] for historical considerations and [START_REF] Brunel | Quantitative investigations of electrical nerve excitation treated as polarization[END_REF] for a recent English translation of Lapique's original paper.

To be more specific, the model describes the subthreshold dynamics of a single neuron membrane's potential and a reset mechanism to account for the onset of an action potential:

A spike occurs whenever a given threshold V T is reached by the membrane potential variable V . Once the firing event occurs, the membrane potential is right away reset to a given value

Neuron model

Figure 1: Schematic representation of a point process model. The model is referred as to be a point process since it takes into account only the time events and the complex shape of the action potential is not explicitly modeled. Point process models mostly focus on the input/ouput relationship, i.e. the relationship between the input spike trains the cell receives via synaptic afferent and its response to it. The model is said stochastic when the input or the firing process is random.

V R . In the subthreshold regime, the membrane potential's dynamics is given by

τ d dt V (t) = -g(V (t) -V L ) + η(t),
where V (t) is the membrane potential at time t, τ is the membrane capacitance, g -the leak conductance, V L -the reversal potential and η(t) -a gaussian white noise, see [START_REF] Burkitt | A review of the integrate-and-fire neuron model: I. homogeneous synaptic input[END_REF] for a recent review and see [START_REF] Izhikevich | Dynamical Systems in Neuroscience[END_REF] for other spiking models. In what follows, we will use a normalized version of the above equation, i.e. we define µ as the bias current and v the membrane's potential which will be given by

110 µ = V L V T , v = V V T , v r = V R V T .
After re-scaling the time in units of the membrane constant g/τ , the normalized model reads

   d dt v(t) = µ -v(t) + ξ(t) If v > 1 then v = v r . (1)
ξ(t) is again a Gaussian white noise stochastic process with intensity σ: The NLIF model was introduced in [START_REF] Knight | Dynamics of encoding in a population of neurons[END_REF] and generalizations of it can be found in the more recent work [START_REF] Gerstner | Spiking neuron models[END_REF]. The first equation in ( 1) is a Langevin equation that contains a deterministic part expressed by the drift term µ -v(t) and a stochastic part in the form of the noise term ξ(t). The second line in (1) describes the onset of an action potential and the reset mechanism.

ξ(t) = 0, ξ(t)ξ(t ) = σδ(t -t ).
One popular way to deal with a SDE is to write down the associated Fokker-Planck equation for the associated probability density function (pdf). In the case of the NLIF model, the associated pdf p(t, v) express the likelihood of finding the membrane potential 120 at a given time t in a value v. Starting with the SDE (1) of Langevin type, the interested reader will found in [START_REF] Bressloff | Stochastic models of intra-cellular transport[END_REF] a rigorous derivation of the associated Fokker-Planck equation. For our specific case, the FP equation takes the following form:

∂ ∂t p(t, v) + Drift ∂ ∂v [(µ -v)p(t, v)] - Diffusion σ 2 2 ∂ 2 ∂v 2 p(t, v) = Reset δ(v -v r )r(t) . (2)
The equation above expresses three different processes: a drift process due to the determin-

p(t, 1) = 0, (3)
which expresses the firing process, and a reflecting boundary condition

lim v→-∞ (-µ + v)p(t, v) + σ 2 2 ∂ ∂v p(t, v) = 0, (4)
which states that there is no flux passing through this boundary.

The firing rate r(t) is defined as the flux at the threshold:

130 r(t) = - σ 2 2 ∂ ∂v p(t, 1). ( 5
)
To uniquely determine a solution, an initial condition is given:

p(0, v) = p 0 (v). (6)
Using the boundary conditions and the expression of r(t) given by ( 5), one can easily check the conservation property of the equation ( 2) by directly integrating it, so that, if the initial condition satisfies

1 -∞ p 0 (v) dv = 1, (7)
then the solution to ( 2)-( 6) necessarily satisfies the normalization condition

1 -∞ p(t, v) dv = 1. ( 8
)
We present in Fig 2 a simulation of the FP model ( 2)-( 6). The numerical results illustrate the stochastic process (1) and the time evolution of its associated probability density. In Fig 2, the time evolution of the density is depicted in the first panel while the black line only gives one particular realization of the Langevin type equation ( 1). Note the effect of noise on the dynamics of the membrane potential: For this particular neuron, the firing events occur whenever the membrane potential reaches the threshold, which is highlighted by the presence of a red dot in our simulation. The firing activity is also represented, the red curve corresponding to the FP equation ( 2)-( 6) and the blue curve to the stochastic process (1).

To get a better understanding of the time evolution of this probability density, we present in Fig. 3 different snapshots of the simulation. Under the drift and the diffusion effects, the density function gives a non zero flux at the threshold, and this flux is reset to v r according to the reset process. This effect can be seen clearly in the third panel of the simulation presented in Fig. 3. Asymptotically, the solution reaches a stationary density, which is shown in the last panel of Fig. 3.

Intrinsic noise. McKendrick-von Foerster formalism

To account for the intrinsic variability and the randomness in the firing process, a typical procedure is to consider deterministic trajectories and to assume that a firing event can occur at any time according to a certain probability rate [START_REF] Gerstner | Time structure of the activity in neural network models[END_REF]. This rate is often called escape rate or stochastic intensity of firing [START_REF] Gerstner | Spiking neuron models[END_REF] and will be denoted by S throughout this paper. Therefore the implementation of intrinsic noise lets the deterministic trajectories, expressed by one of the known single neuron models, unaffected, but instead influences the firing time, which is no longer deterministic. In this setting, a neuron may fire even though the formal threshold has not been reached yet or may stay quiescent after the formal threshold has been crossed.

Usually the expression of S depends on the momentary distance of the membrane potential to the formal threshold of firing, and therefore the process can be fully characterized by the amount of time passed by since the last action potential.

A random variable described in this way is known as a renewal process which is a class of stochastic processes characterized by the fact that the whole history of events that had occurred before the last one can be forgotten [START_REF] Cox | Renewal Theory[END_REF]. Note that the use of the renewal theory in neuroscience is justified by the fact that the neuron is not supposed to exhibit strong adaptation, see [START_REF] Deger | Fluctuations and information filtering in coupled populations of spiking neurons with adaptation[END_REF] for the case of adapting neurons. In [START_REF] Gerstner | Spiking neuron models[END_REF], [START_REF] Plesser | Noise in integrate-and-fire neurons: from stochastic input to escape rates[END_REF] there is treated the case of time dependent input which leads to the use of the more general case of non-stationary renewal systems. We will restrict in this paper to the stationary renewal theory since we only consider the case of time-independent stimulus, and therefore, the only variable on which the escape rate will depend is the age a, which stands for the time elapsed since the last spike. In such a framework, the neuron follows a stochastic process where the probability of surviving up to a certain age, P (a), is given by

P (a) = exp{- a 0 S(s) ds}, (9)
where the specific choice of S determines the properties of the renewal process. We will come back with details about the form of S(a) in the following sections.

To better understand the renewal process, one can focus on the probability density function that describes the relative likelihood for this random variable to take on a given value.

This description leads to a so-called age-structured system [START_REF] Gerstner | Spiking neuron models[END_REF] which consists in a partial differential equation with non-local boundary condition that is famous in the field of population dynamics, and which, in our specific case, is in the form of the McKendrick-von Foerster model [START_REF] Foerster | Some remarks on changing populations[END_REF]. Denoting by n(t, a) the probability density for a neuron to have at time t the

) 10
Because once a neuron triggers a spike its age is reset to zero, the natural boundary condition to be considered is

Reset n(t, 0) = r(t), (11)
where r(t) is the firing rate and is given by

r(t) = +∞ 0 S(a)n(t, a) da. ( 12
)
To completely describe the dynamics of n, an initial distribution is assumed known:

185 n(0, a) = n 0 (a). ( 13
)
From now on we will call solution of the age structured (AS) system, the well defined solution of [START_REF] Cox | Renewal Theory[END_REF] with initial condition [START_REF] Dumont | Synchronization of an excitatory integrate-and-fire neural network[END_REF] and boundary condition given by ( 11)- [START_REF] Dumont | Population density models of integrate-and-fire neurons with jumps, well-posedness[END_REF]. Using the boundary condition and the expression of r(t) given by ( 12), one can check easily the conservation property of the equation ( 2) by integrating it, so that if the initial condition satisfies

) 14
We present in Fig. 4 a simulation of the escape rate model and of its associated probability density function. Note that the main distinction on the stochastic aspect is that the noise does not act on the trajectory but only on the initiation of action potential. Again, we have made a comparison between the stochastic process (black curve) and the evolution of the density function (red curve). The simulation starts with a Gaussian as initial condition (the first panel of Fig. 5). Under the influence of the drift term, the density function advances in age, which is clearly seen in the upper plots of Fig. 5. After the spiking process, the age of the neuron is reset to zero. The effect is well perceived in the lower panels of Fig. 5. As expected from the model, the density function converges to an equilibrium state.

Chapman-Kolmogorov equation

The models introduced above, ( 2)-( 6) respectively ( 10)-( 13), have been shown to exhibit the same statistical activity [START_REF] Gerstner | Spiking neuron models[END_REF], even though the stochastic processes that they illustrate are conceptually different. To explain this behavior, we have been looking in which way the solutions of the two models can be analytically related.

The integral transforms that will be defined later on and which give this connection make use of the solutions to backward and forward Chapman-Kolmogorov (CK) equations for the stochastic process described by the NLIF model defined on an inter-spike interval.

Before tackling this issue, we will first remind few theoretical considerations about the CK equation, how the forward and backward systems are derived from it and the interpretations of the solutions to both systems. For a complete derivation and analysis of it, we refer to the classic work [START_REF] Gardiner | Handbook of Stochastic Method for Physics[END_REF], but we also refer to [START_REF] Bressloff | Stochastic models of intra-cellular transport[END_REF] for discussions about the forward and backward CK equations in biological processes.

Markov property

The Chapman-Kolmogorov equation is formulated in terms of conditional probabilities and it is built on the main assumption that the stochastic process in question satisfies the Markov property: The future states of the system depend only on the present state. Defining the conditional probability p(t, v|s, w) as the probability of a neuron that has not discharged up to time t to be in the state v at time t given that it started at time s from w, the CK equation reads:

p(t, v|s, w) = p(t, v|t , v )p(t , v |s, w) dv . (15)
Roughly speaking, the CK equation simply says that, for a Markov process, the transition from (s, w) to (t, v) is made in two steps: first the system moves from w to v at an intermediate time t and then from v to v at time t; the transition probability between the two states is calculated next by integrating over all possible intermediate states.

Starting from the integral CK equation ( 15), the forward CK equation, also known as the FP equation, as well as the backward CK equation are derived under the assumption of sample paths continuity of the stochastic process ( [START_REF] Gardiner | Handbook of Stochastic Method for Physics[END_REF]). While the FP equation is obtained by considering the evolution with respect to present time t, the backward CK equation is derived when considering the time development of the same conditional probability with respect to initial times.

To keep a certain level of simplicity, we shall not refer in this paper to the general form of both equations, and we refer to [START_REF] Gardiner | Handbook of Stochastic Method for Physics[END_REF] for the general case as well as for the derivation and interpretations in particular cases. We remind that, since we have considered in this paper the constant stimulus case, we deal with a stationary stochastic process. A stationary process has by definition the same statistics for any time translations, which implies that the associated joint probabilities densities are time invariant too; due to the relation between the conditional and joint probabilities and also to the Markov property, this allows expressing the conditional probabilities in the CK equation in terms of time differences. In particular, the FP equation as well as the backward CK equations can be written as in the following.

Forward and backward Chapman-Kolmogorov equations

We shall consider first the conditional probability for a neuron that started at the last firing time s from the corresponding state v r to be at time t in a state v; since in the following we shall keep the initial variables fix, and use the property of stationary processes, we can write down the associated FP equation in terms of a joint probability. We deal therefore with the following FP equation for the probability density function

ϕ(a, v) := p(t -s, v|0, v r )
for a neuron to have at age a := t -s the membrane's potential value v, where we have skipped from the notation the variables which are kept fix:

∂ ∂a ϕ(a, v) + ∂ ∂v [(µ -v)ϕ(a, v)] - σ 2 2 ∂ 2 ∂v 2 ϕ(a, v) = 0. ( 16
)
The initial condition is naturally given by

ϕ(0, v) = δ(v -v r ). (17)
Also, due to the firing process at threshold value, an absorbing boundary condition is imposed:

ϕ(a, 1) = 0. (18)
A reflecting boundary is also required:

v→-∞ (v -µ)ϕ(a, v) + σ 2 2 ∂ ∂v ϕ(a, v) = 0. ( 19
)
Note that, in contrast to the model ( 2)-( 6), here there is no conservation property required;

it is suitable therefore to think of ϕ as a probability density on an inter-spike interval, thus the reset term in (2) does not appear in this case. Nevertheless, the flux at the threshold is a quantity of big interest since it gives the inter-spike interval distribution (ISI) function,

250
that was introduced generically in the previous section, for a neuron that started at age zero from the reset value v r

ISI(a) = - σ 2 2 ∂ ∂v ϕ(a, 1). ( 20
)
The central notion of the renewal theory is the interval distribution of the events; it has the role of predicting the probability of the next event to occur in the immediate time interval.

In neuroscience context, the popularized notion is the inter-spike interval distribution. There 255 is a tight connection between the ISI distribution and the survivor function P due to their interpretations. Starting with [START_REF] Gerstner | Spiking neuron models[END_REF], one can define next the survivor function P as the probability to survive up to age a for a neuron that started at age zero from the position v r as

P(a) = 1 -∞ ϕ(a, v) dv. ( 21
)
Note that the relation reminded above that takes place between these two functions, i.e.

ISI(a) = - ∂ ∂a 1 -∞ ϕ(a, v) dv, (22)
is verified by integrating ( 16) with respect to v on the whole potential values interval and using the boundary conditions. Also, one may check directly that

1 -∞ ϕ(a, v) dv = 1 - a 0 ISI(s) ds.
Both these probabilities are defined here as they were introduced in [START_REF] Gerstner | Spiking neuron models[END_REF].

In the case of backward CK equation, the conditional probability for a neuron that started at age zero from a potential value v to survive up to the age a is considered,

ψ(a, v) := 1 -∞ p(t -s, w|0, v) dw.
We did use the same name for the time variable in the definition of ψ, i.e. age, for the sake of fluency, although inhere a represents simply the time passed since the initial time considered,

a = t -s,
and not the time passed since the last spike. Note that this probability is the same with the probability that at age a the present state w for a neuron that started at age zero from v did not reach yet the threshold. With this respect, the backward CK equation gives the evolution with respect to initial states:

∂ ∂a ψ(a, v) -(µ -v) ∂ ∂v ψ(a, v) - σ 2 2 ∂ 2 ∂v 2 ψ(a, v) = 0. ( 23
)
Obviously the suitable initial condition is

ψ(0, v) = 1, (24)
and, due to the same interpretation, a neuron that reaches the threshold value at age a has the probability of survival

ψ(a, 1) = 0. ( 25
)
The other boundary condition,

lim v→-∞ ∂ψ(a, v) ∂v = 0, ( 26
)
is obtained due to the corresponding boundary condition to the FP equation ( 19) and to the duality of the functions ϕ and ψ, see [START_REF] Gardiner | Handbook of Stochastic Method for Physics[END_REF] for details about the choice of boundary conditions 275 for forward as well as for backward CK equation.

Note that, in contrast with the survivor function introduced above as in [START_REF] Plesser | Noise in integrate-and-fire neurons: from stochastic input to escape rates[END_REF], ψ expresses the probability of survival for a neuron that started at age zero from a potential value v;

then, obviously, the following relation must hold:

1 -∞ ϕ(a, w) dw = ψ(a, v r ), (27)
equality that checks out immediately due to the definitions of both functions as solutions to equations ( 16) and ( 23) and using integration by parts:

0 = a 0 1 -∞ ∂ψ(t, v) ∂t ϕ(a -t, v)dvdt - a 0 1 -∞ ψ(t, v) ∂ϕ ∂a (a -t, v)dvdt = ψ(a, v r ) - 1 -∞ ϕ(a, v)dv.
Keeping in mind the definition of the survival function as

S(a) = - ∂ ∂a 1 -∞ ϕ(a, w) dw 1 -∞ ϕ(a, w) dw = ISI(a) 1 - a 0 ISI(s) ds (28)
which has the interpretation that, in order to emit a spike, the neuron has to stay quiescent in the interval (0, a) and then fire at age a, see [START_REF] Gerstner | Spiking neuron models[END_REF]. In Fig. 7, numerical simulations of the age-dependent death rate is presented. Let us notice that S clearly defines a positive function that converges towards a constant.

Analytical links between the FP and AS systems

Having set the theoretical framework, we are now ready to introduce our results regarding the connections between the two systems that correspond to different mathematical treatments of variability.We stress that both transformations given here have been obtained for an escape rate defined by the use of the ISI function defined by [START_REF] Izhikevich | Dynamical Systems in Neuroscience[END_REF].

More precisely, we will define integral transforms that map one-to-another the solutions to the two models considered, i.e. ( 2)-( 6) respectively ( 10)-( 13). One single important assumption over the initial probability densities of the systems is necessary for obtaining the results below, that is:

p 0 (v) = ∞ 0 ϕ(a, v) 1 -∞ ϕ(a, w) dw n 0 (a) da. ( 29
)
As we stressed before, the FP equation does not give any information about the firing times, merely it can indicate the condition for which a spike is triggered. Our relation simply assumes that the repartition of the potentials at a given initial time corresponds to an initial probability density of ages. It is assumed in this way that the neuron has fired at least once.

4.1.

From escape rate to diffusive noise.

In [START_REF] Dumont | Noisy threshold in neuronal models: connections with the noisy leaky integrate-and-fire model[END_REF], we have already studied the analytical transformation of the solution to the model expressing the internal variability into the solution corresponding to the external one.

Namely, we have shown that, once the relation ( 29) takes place, then for all times t, a similar relation between the density of potentials respectively ages takes place. The result states:

Theorem 1. Let p a solution to ( 2)-( 6) and n a solution to ( 10)-( 13), and p 0 (v) and n 0 (a)

two corresponding initial densities, respectively. Then if p 0 and n 0 satisfy

p 0 (v) = +∞ 0 ϕ(a, v) 1 -∞ ϕ(a, w) dw n 0 (a) da, (30)
the following relation holds true for any time t > 0:

p(t, v) = +∞ 0 ϕ(a, v) 1 -∞ ϕ(a, w) dw n(t, a) da. (31)
Here, ϕ(a, v) is the solution to ( 16)- [START_REF] Gerstner | Time structure of the activity in neural network models[END_REF].

In our above quoted work, the exact meaning of the solutions involved in [START_REF] Nykamp | A population density appraoch that facilitates large-scale modeling of neural networks : analysis and an application to orientation tuning[END_REF] and, consequently, of the integral (31) is given. A similar relation between the respective steady states of the models ( 2)-( 6), ( 10)-( 13) has been also proven.

Note that, the probabilistic meaning of the integral transform given in Theorem 1 can be interpreted using Bayes' rule. Since ϕ(a, v) is the probability density for a neuron to be at age a and at potential v, the kernel of the transform can be interpreted as the probability density for a neuron to be at potential value v given that it survived up to age a. Then, the product of this kernel with the solution n(t, a), which denotes the probability density at time t in state a, integrated over the all possible states a, gives indeed the probability density to be at time t in the state v.

The proof of Theorem 1. as well as additional results can be found in our previous work [START_REF] Dumont | Noisy threshold in neuronal models: connections with the noisy leaky integrate-and-fire model[END_REF].

From diffusive noise to escape rate

Transforming the solution to the FP model into the solution to the AS system is a little bit trickier. Before stating our result, let us make few comments about the features of this problem. Note first that the very nature of the age a contains all the information about time that is needed to properly define the integral transform [START_REF] Nykamp | A population density appraoch that facilitates large-scale modeling of neural networks : analysis and an application to orientation tuning[END_REF]. On the contrary, to define an inverse transform, one faces the problem of having a kernel that must depend on time. A second important aspect about the nature of the AS formalism is that the variable a also entails information about the last firing moment. Indeed, attributing an age to a neuron presupposes that the considered neuron has already initiated an action potential. From our perspective, the membrane potential variable v does not carry out such information.

The main result of the paper is the following:

Theorem 2.
Let p be a solution to ( 2)-( 6) and n be a solution to ( 10)-( 13), so that the compatibility condition (29) between the corresponding initial states p 0 and n 0 takes place.

Then, for any t such that 0 ≤ a < t, the following relation holds true:

n(t, a) = - ∂ ∂a 1 -∞ ψ(a, v)p(t -a, v) dv. (32)
Moreover, for any a ≥ t, a global relation between the corresponding probabilities at time t takes place:

+∞ t n(t, a) da = 1 -∞ ψ(t, v)p 0 (v)dv. (33)
Here, ψ(•, v) is the solution of the dual problem ( 23)-( 26).

Remark 1. All the solutions in Theorem 2 are taken in distributional sense. In this paper, we intend to maintain a certain level of simplicity by not going into details about the functional spaces in which the solutions are considered. We do refer for rigorous definitions of the solutions to the system (2)-( 6) to [START_REF] María | Analysis of nonlinear noisy integrate & fire neuron models: blow-up and steady states[END_REF]. Also, a more general case of the system (10)- [START_REF] Dumont | Synchronization of an excitatory integrate-and-fire neural network[END_REF] has been considered in [START_REF] Pakdaman | Dynamics of a structured neuron population[END_REF]. In what follows, therefore, we preferred to keep the computations at a formal level so that to preserve the intuitional meaning of the functions involved.

Remark 2. Note that relation [START_REF] Ostojic | Synchronization properties of networks of electrically coupled neurons in the presence of noise and heterogeneities[END_REF] does only provide a global relation between n(t, a) for a > t and p 0 (v), but not an exact expression of n(t, a) as a function of p 0 (v); the explicit 345 form of n(t, a) is known only if n 0 (a) is considered known.

Proof We start our proof by noting that the explicit solution of the AS system knowing r, can be calculated and is given by

n(t, a) =    P(a) P(a-t) n 0 (a -t), t ≤ a, P(a)r(t -a), t > a. (34)
During the first step of the proof we will denote by n(t, a) the quantity ∂ ∂a

We will prove that it satisfies [START_REF] Cox | Renewal Theory[END_REF] with the boundary condition [START_REF] Deger | Fluctuations and information filtering in coupled populations of spiking neurons with adaptation[END_REF]. Let us consider the

integral 350 1 -∞ ψ(a, v)p(t -a, v) dv.
A straightforward computation gives:

∂ ∂a 1 -∞ ψ(a, v)p(t -a, v) dv = 1 -∞ p(t -a, v) ∂ ∂a ψ(a, v) dv + 1 -∞ ψ(a, v) ∂ ∂a p(t -a, v) dv.
Due to the fact that ψ is solution to [START_REF] Knight | Dynamics of encoding in neuron populations: Some general mathematical features[END_REF], the first term in the right hand side can be written as

1 -∞ p(t -a, v) ∂ ∂a ψ(a, v) dv = 1 -∞ p(t -a, v) (µ -v) ∂ ∂v ψ(a, v) + σ 2 2 ∂ 2 ∂v 2 ψ(a, v) dv.
Integrating by parts both terms in the right hand side we get:

1 -∞ p(t -a, v) ∂ ∂a ψ(a, v) dv = σ 2 2 p(t -a, v) ∂ ∂v ψ(a, v) 1 -∞ - 1 -∞ σ 2 2 ∂ ∂v ψ(a, v) ∂ ∂v p(t -a, v) dv + [(µ -v)p(t -a, v)ψ(a, v)] 1 -∞ - 1 -∞ ψ(a, v) ∂ ∂v [(µ -v)p(t -a, v)] dv.
Applying again integration by parts in the second integral in the right hand side of the last expression and using the absorbing boundary condition for p and the boundary condition for ψ, it follows that:

1 -∞ p(t -a, v) ∂ ∂a ψ(a, v) dv = [(µ -v)p(t -a, v)ψ(a, v)]| -∞ - σ 2 2 ψ(a, v) ∂ ∂a p(t -a, v) 1 -∞ + 1 -∞ ψ(a, v) σ 2 2 ∂ 2 ∂v 2 p(t -a, v) - ∂ ∂v [(µ -v)p(t -a, v)] dv.
Using again the absorbing boundary condition and the reflecting boundary condition for the flux of the problem in p, we finally get that the first term is indeed

1 -∞ p(t -a, v) ∂ ∂a ψ(a, v) dv = 1 -∞ ψ(a, v) σ 2 2 ∂ 2 ∂v 2 p(t -a, v) - ∂ ∂v [(µ -v)p(t -a, v)] dv.
On the other hand, the second term in the equation can be written equivalently as

1 -∞ ψ(a, v) ∂ ∂a p(t -a, v) dv = 1 -∞ ψ(a, v) - ∂ ∂t p(t -a, v) dv, which implies that 1 -∞ ψ(a, v) ∂ ∂a p(t -a, v) dv = 1 -∞ ψ(a, v) ∂ ∂v [(µ -v)p(t -a, v)] - σ 2 2 ∂ 2 ∂v 2 p(t -a, v) -δ(v -v r )r(t -a) dv,
where r is the firing rate from the FP model. Thus, by adding the two terms, we finally get something expected

- ∂ ∂a 1 -∞ ψ(a, v)p(t -a, v) dv = ψ(a, v r )r(t -a),
expression that can also be written, due to [START_REF] Manwany | Detecting and estimating signals in noisy cable sstructure. i: Neuronal noise sources[END_REF], as

- ∂ ∂a 1 -∞ ψ(a, v)p(t -a, v) dv = r(t -a) 1 -∞
ϕ(a, v) dv = P(a)r(t -a).

Since, as reminded at the beginning of the proof, the right-hand side of the above relation expresses the solution of [START_REF] Cox | Renewal Theory[END_REF] with the boundary condition [START_REF] Deger | Fluctuations and information filtering in coupled populations of spiking neurons with adaptation[END_REF] for t > a, the first step of the proof is complete.

At this point, it remains to be proven that the relations ( 33) and ( 11)-( 12) take place.

We consider the region a ≥ t and we start by reminding that integration of (10) with initial condition [START_REF] Dumont | Synchronization of an excitatory integrate-and-fire neural network[END_REF] provides: P(a) P(a -t) n 0 (a -t), a ≥ t.

With this respect, given an initial density of ages, the solution for a ≥ t is completely determined by the initial density n 0 and the knowledge of ϕ, since, as reminded,

P(a) = 1 -∞ ϕ(a, v) da.
We wish though to relate the solution n for this case to the initial density of potentials, and to do so, we will see that, as anticipated, the sole relation ( 30) is sufficient.

355

Note for the beginning that the functions ψ and ϕ are adjoint; let us consider next the integral

1 -∞ ψ(t, v)ϕ(a -t, v) dv. (35)
One can show that the above integral does not depend explicitly of t by computing

∂ ∂t 1 -∞ ψ(t, v)ϕ(a -t, v) dv = 1 -∞ ∂ ∂t ψ(t, v)ϕ(a -t, v) + ψ(t, v) ∂ ∂t ϕ(a -t, v) dv. ( 36
)
Since, due to the fact that ψ is solution to backward CK equation, one can write down

1 -∞ ∂ ∂t ψ(t, v)ϕ(a -t, v) dv = 1 -∞ ϕ(a -t, v) (µ -v) ∂ ∂v ψ((t, v) + σ 2 2 ∂ 2 ∂v 2 ψ(t, v) dv,
integrating by parts in the right-hand side and using the boundary conditions, it follows that

360 1 -∞ ∂ ∂t ψ((t, v)ϕ(a -t, v) dv = 1 -∞ ψ(t, v) ∂ ∂v [-(µ -v)ϕ(a -t, v)] + σ 2 2 ∂ 2 ∂v 2 ϕ(a -t, v) dv = - 1 -∞ ψ(t, v) ∂ ∂t ϕ(a -t, v) dv.
Replacing therefore the last expression in [START_REF] Plesser | Noise in integrate-and-fire neurons: from stochastic input to escape rates[END_REF], we get indeed that

∂ ∂t 1 -∞ ψ(t, v)ϕ(a -t, v) dv = 0,
which implies that the above integral is only a function of a on the interval [0, a]. We obtain therefore, due to initial conditions

ϕ(0, v) = δ(v -v r ), ψ(0, v) = 1, that 1 -∞ ψ(t, v)ϕ(a -t, v) dv t=0 = 1 -∞ ϕ(a, v) dv = 1 -∞ ψ(t, v)ϕ(a -t, v) dv t=a .
Going back to the expression of the solution for this case, one can write down then equivalently

n(t, a) = n 0 (a -t) P(a -t) 1 -∞ ψ(t, v)ϕ(a -t, v) dv, a ≥ t,
and, integrating the last relation over [t, ∞):

365 ∞ t n(t, a) da = ∞ t n 0 (a -t) P(a -t) 1 -∞ ψ(t, v)ϕ(a -t, v) dv da.
Making the change of variable a -→ a = a -t and changing the order of integration, we finally get In order to prove [START_REF] Deger | Fluctuations and information filtering in coupled populations of spiking neurons with adaptation[END_REF] where r(t) is given by ( 12), we have first to prove the conservation property for n(t, a):

∞ 0 n(t, a) da = 1, ∀t > 0. ( 37
)
This is easily obtained by noting that

370 ∞ 0 n(t, a) da = t 0 n(t, a) da + ∞ t n(t, a) da = t 0 - ∂ ∂a 1 -∞ ψ(a, v)p(t -a, v) dv da + 1 -∞ ψ(t, v)p 0 (v) dv = 1 - 1 -∞ ψ(t, v)p 0 (v) dv + 1 -∞ ψ(t, v)p 0 (v) dv = 1.
Let us notice now that the integral in the right hand side of ( 12) should be split in two parts corresponding to the two branches of the solution (34):

(a -t) n 0 (a -t) ∞ t + ∞ t P (a) ∂ ∂a n 0 (a -t) P (a -t) da = P (t)n 0 (0) + ∞ t P (a) ∂ ∂t n 0 (a -t) P (a -t) dt = - ∂ ∂t ∞ t P (a) P (a -t) n 0 (a -t) da.
Adding up the two integrals, we get then that:

∞ 0 S(a)n(t, a) da = r(t) - ∂ ∂t ∞ 0 n(t, a) da,
where we used the expressions of the solution n(t, a) for a < t and a > t, respectively.

But since we do have, by the relation [START_REF] Shadlen | Noise, neural codes and cortical organization[END_REF], that the conservation property for n(t, a) takes place, it follows that Remark 3. In [START_REF] Dumont | Noisy threshold in neuronal models: connections with the noisy leaky integrate-and-fire model[END_REF] it was shown that, even if the relation between the initial states (30) is not satisfied, nevertheless, ( 31) is satisfied asymptotically as t goes to infinity. In a similar way, Theorem 2 shows that even if n(t, a) is well defined only for 0 < a < t, the relation [START_REF] Ostojic | Synchronization properties of networks of electrically coupled neurons in the presence of noise and heterogeneities[END_REF] assures that the influence of the undetermined part of the solution for a > t goes to zero as t goes to infinity.

Remark 4. Given the interpretations of all the functions involved in the relation [START_REF] Omurtag | On the simulation of large population of neurons[END_REF], one can interpret the result probabilistically in the following way: The probability density for a neuron to have the age a at time t is given by the flux generated by the probability density at potential v and time t -a that survived up to age a. The left-hand side in [START_REF] Ostojic | Synchronization properties of networks of electrically coupled neurons in the presence of noise and heterogeneities[END_REF] has the interpretation of the probability that the cell has at time t an age a belonging to [t, +∞).

With this respect, the relation simply says that, under assumption (29), the integral over v of the probability density at time 0 which survived up to t without firing gives indeed the probability that a > t.

Conclusions

How to give a good analytical treatment to neural variability? Theories abound and provide very different ways to deal with the stochastic aspect of nervous cells. Among them, probably the most popular approaches are the NLIF model with its associated FP equation and the escape rate model with its corresponding AS system. The FP equation is commonly used to represent the dynamics of the probability density function with respect to the membrane potential and has the advantage that the parameters of the system have bio-physical interpretations and could be measured/identified by experimentalists. Anyway, finding the membrane's potential density is rarely of use in practice. Models that give the evolution of a density with respect to the time elapsed since the last spike time (AS systems) have been considered in the literature: [START_REF] Gerstner | Spiking neuron models[END_REF][START_REF] Pakdaman | Dynamics of a structured neuron population[END_REF][START_REF] Khashayar Pakdaman | Relaxation and selfsustained oscillations in the time elapsed neuron network model[END_REF]. These models though, in reverse, are not linked to the bio-physical properties of the cell since they rely only on the assumption that a neuron fires with the probability given by the escape rate. The use of escape rate instead of consideration of diffusive noise has another obvious practical advantage: while FP equation ( 2) along with its boundary conditions and the source term given by ( 5) is difficult to be handled from a mathematical point of view, the AS model ( 10) -( 13) can be easily integrated and leads to an analytical solution. This model is well known in the age-structured systems theory and a large amount of papers have treated similar systems.

We mention that qualitative results for a similar model in neuroscience context have been obtained in [START_REF] Pakdaman | Dynamics of a structured neuron population[END_REF][START_REF] Khashayar Pakdaman | Relaxation and selfsustained oscillations in the time elapsed neuron network model[END_REF].

It has been shown in [START_REF] Plesser | Noise in integrate-and-fire neurons: from stochastic input to escape rates[END_REF] that the NLIF model can be mapped approximately onto an escape-rate model. In a recent study [START_REF] Dumont | Noisy threshold in neuronal models: connections with the noisy leaky integrate-and-fire model[END_REF], we have shown an analytical connection between the two models. We have proven there the existence of an exact analytical transform of the solution of the FP system into the solution to the AS system. To our knowledge, such a result has not been proven before. The present paper is intended as a sequel step toward the completion of the analytical link between these two known models.

The importance of our analytical connection consists in the fact that it gives a way to rely the density of ages to those of membrane potentials and underline in which way these densities depend on each-other; it is therefore, a first attempt to attribute age to potentials densities, information that is not carried out in the solution of the FP equation. Reversely, a more practical aspect of our transforms consists in the fact that it allows a simpler analysis of the properties of the solutions by noticing in which way their properties transfer one to another via the integral transforms proposed here. Finally, from a biological point of view, our theoretical result shows that the two mathematical formalization of variability are similar.

The usual way to treat the intrinsic noise is by considering a hazard function that depends on the momentary distance between the deterministic trajectory of the membrane potential and the formal spiking threshold [START_REF] Gerstner | Spiking neuron models[END_REF]. The link between an escape rate given in a general form as a function of this distance and a corresponding diffusion process associated to one of the single neuron stochastic models is not straightforward at all. In such a situation, there is no obvious way to expect any links between the hazard rate and the effect of a noisy stimulus. In this paper, we have shown that the particular choice of a hazard function built on the ISI function given by [START_REF] Izhikevich | Dynamical Systems in Neuroscience[END_REF] gives the desired link. Which means that to make the link possible, the escape rate has to be built via the properties of the FP equation, which characterizes the noisy stimulation.

We have to stress though that the results obtained here have been proven in the case of time independent stimulus. The case we considered is known in the framework of renewal systems as a stationary process. A possible extension of the present work for the case of time dependant stimulus remains thus for us an open issue to be investigated. Such consideration would allow us to discuss the case of interconnected neurons and thus get a better understanding of neural networks' dynamics. This is a current work in progress.

The case of coupled neurons puts additional problems since the corresponding ISI function will not be only age-dependent, but also time dependent. In particular, the corresponding AS system has a time dependent death rate, S(t, a), instead of S(a) as in our present manuscript.

Furthermore, in the case of a coupled population, our approach does not generalize so easily since the FP equation as well as the AS model will become nonlinear. However, our strategy still works and this case is in working progress. Finally, we would like to mention that our strategy is not adapted for the case of a stochastic process with different statistical properties.

Figure 2 :

 2 Figure 2: Simulation of the stochastic process (1) and of its associated FP equation (2). A) Evolution in time of the density probability; the brightness is proportional to the probability density and the black line illustrates one realization of the process. B) Raster plot depicting the precise spiking time of the cell over different trials. The red dots correspond to the particular realization of panel A. C) The firing rate given by the FP equation, (5), in red and by many realizations of the Langevin equation (1). A gaussian was taken as initial condition; the parameters of the simulation are: v r = 0.3, µ = 20, σ = 0.4.

Figure 3 :

 3 Figure 3: Simulation of the stochastic process (1) and of its associated FP system (2)-(5). A gaussian was taken as initial condition; the parameters of the simulation are: v r = 0.3, µ = 20, σ = 0.4. The plots show the evolution in time of the solution respectively t = 0, t = 0.1, t = 0.3, t = 0.5, t = 0.7, t = 7 for the respective panel A-B-C-D-E-F.

Figure 4 :Figure 5 :

 45 Figure4: Simulation of the stochastic process[START_REF] María | Analysis of nonlinear noisy integrate & fire neuron models: blow-up and steady states[END_REF] as well as of its associated AS equation (10)-[START_REF] Dumont | Synchronization of an excitatory integrate-and-fire neural network[END_REF]. A) Evolution in time of the density probability; the brightness is proportional to the probability density and the black line illustrates one realization of the process. B) Raster plot depicting the precise spiking time of the cell over different trials. The red dots correspond to the particular realization of panel A. C) The firing rate given by the AS equation[START_REF] Cox | Renewal Theory[END_REF] in red and by many realizations of[START_REF] María | Analysis of nonlinear noisy integrate & fire neuron models: blow-up and steady states[END_REF]. A gaussian was taken as initial condition. The function S used in the simulations is given by[START_REF] Daniel Millman | Self-organized criticality occurs in non-conservative neuronal networks during 'up' states[END_REF] with parameters v r = 0.3, µ = 20, σ = 0.4.

+∞ 0 n 0

 00 (a) da = 1, the solution satisfies the normalization condition +∞ 0 n(t, a) da = 1.(

Figure 6 :

 6 Figure 6: Simulation of the stochastic process and of its associated FP equation (16). A) Evolution in time of the density probability;the brightness is proportional to the probability density and the black line illustrates one realization of the process. B) Raster plot depicting the precise spiking time of the cell over different trials. The red dot corresponds to the particular realization of panel A. C) The ISI obtained from the FP equation (16) in red, and by many realizations of the Langevin equation (1). The parameters of the simulation are: v r = 0.3, µ = 20, σ = 0.4.

  lim

Figure 7 :

 7 Figure 7: Illustration of the ISI and the spiking rate function S. This illustration is obtained via a numerical simulation of the function S(a) given by (28) in the panel A and its corresponding ISI given by (22) in the B panel. The parameters of the simulation are v r = 0.7, µ = 5, σ = 0.1

  , v)p 0 (v) dv.

0 - 0 P 0 P 0 P

 0000 n(t, a) da, where in the first integral we have the solution given by the transform for t > a, i.e.n(t, a) = -v)p(t -a, v) dv = P (a)r(t -a).Taking the first integral, one can write down t P (a)P (a) n(t, a) da = -t (a)r(t -a) da = -P (t)r(0) + P (0)r(t) + t (a) ∂ ∂a r(t -a) da = r(t) -P (t)r(0) -(a)r(t -a) da.Exactly in the same way, by noting that P (a -t) n 0 (a -t) da, performing integrations by parts we arrive to: ∞ t -P (a) P (a) n(t, a) da = -P (a) P

380 ∞ 0 S

 0 (a)n(t, a) da = r(t) = n(t, 0), which ends the proof.

bimodality in interspike interval distributions. J. Theor. Biol., 105:345-368, 1983.


Assyr Abdulle
email: assyr.abdulle@epfl.ch

Martin E Huber
email: martin.huber@epfl.ch

Simon Lemaire
email: simon.lemaire@epfl.ch

An optimization-based numerical method for diffusion problems with sign-changing coefficients

A new optimization-based numerical method is proposed for the solution of diffusion problems with sign-changing conductivity coefficients. In contrast to existing approaches, our method does not rely on the discretization of a stabilized equation and the convergence of the scheme can be proved without any symmetry assumption on the mesh near the interface where the conductivity sign changes.

Résumé

Une méthode d'optimisation pour des problèmes de diffusion avec changement de signe. Nous proposons une nouvelle méthode, basée sur la résolution d'un problème de minimisation, pour l'approximation de problèmes de diffusion avec changement de signe. Cette approche, qui tire profit d'une reformulation du modèle initial sous la forme d'un problème de transmission, ne repose pas sur la discrétisation d'une équation stabilisée, et la convergence de la méthode est obtenue sans hypothèse de symétrie du maillage dans un voisinage de l'interface où la conductivité change de signe.

Version française abrégée

Dans cette note, nous introduisons une méthode d'optimisation pour l'approximation numérique de problèmes de diffusion dont la conductivité change de signe dans le domaine. La résolution numérique efficace de ce genre de problèmes est importante pour de nombreuses applications (e.g., super-lentilles, invisibilité), mais les méthodes existantes ne sont pour l'instant pas satisfaisantes. Dans [START_REF] Chesnel | T-coercivity and continuous Galerkin methods: application to transmission problems with sign-changing coefficients[END_REF], les deux approches envisagées reposent (i) sur la discrétisation d'une équation stabilisée, pour laquelle les taux de convergence obtenus sont sous-optimaux, ou (ii) sur des hypothèses de symétrie du maillage autour de l'interface où la conductivité change de signe, exigences pouvant s'avérer très contraignantes pour des interfaces générales (voir [START_REF] Bonnet-Ben Dhia | Mesh requirements for the finite element approximation of problems with sign-changing coefficients[END_REF]) ou en 3D. La méthode numérique introduite ici, qui utilise une reformulation du modèle initial en un problème de transmission, ne repose pas sur l'ajout de dissipation à l'équation, et nous montrons sa convergence pour des problèmes elliptiques présentant un changement de signe sans aucune hypothèse de symétrie sur le maillage. Nous notons que cette méthode numérique a pour la première fois été introduite par Gunzburger et al. [START_REF] Gunzburger | An optimization-based domain decomposition method for partial differential equations[END_REF] (voir aussi [START_REF] Gunzburger | Solution of elliptic partial differential equations by an optimization-based domain decomposition method[END_REF]), dans un contexte de décomposition de domaine pour des équations elliptiques classiques, sans preuve de convergence. L'application de cet algorithme à des problèmes elliptiques présentant un changement de signe est introduite dans cette note, et la convergence de la méthode est démontrée.

Introduction

Partial differential equations with sign-changing coefficients play an increasingly important role in the modeling of metamaterials, with applications ranging from superlensing to cloaking. In this paper, we consider in an open bounded polytopal domain Ω Ă R d , d P t2, 3u such that Ω " Ω e Y Ω i , with Ω e , Ω i disjoint open polytopal subsets of Ω with nonzero measure, the following sign-changing diffusion problem

´divps K∇ũq " f in Ω, ũ " 0 on BΩ, (1)
where K is a symmetric, uniformly elliptic and bounded matrix-valued field, and s : Ω Ñ R is a sign function such that s e :" s |Ωe " 1 and s i :" s |Ωi " ´1. The subscripts 'e' and 'i' respectively stand for the exterior (for which we assume that BΩ e X BΩ has nonzero measure) and the interior subdomains, and we denote their interface Γ " BΩ e X BΩ i . We assume that f P L 2 pΩq and study the following weak formulation of Problem (1): Find ũ P H 1 0 pΩq such that apũ, ϕq :" ps K∇ũ, ∇ϕq Ω " pf, ϕq Ω @ϕ P H 1 0 pΩq.

Several approaches have been developed to study the well-posedness of Problem [START_REF] Bernardi | Discrétisations variationnelles de problèmes aux limites elliptiques[END_REF]. We first mention the T-coercivity theory of Bonnet-Ben Dhia et al. [START_REF] Bonnet-Ben Dhia | T-coercivity for scalar interface problems between dielectrics and metamaterials[END_REF]. In this context, well-posedness (in the classical Hadamard sense or in the Fredholm sense) of Problem ( 2) is equivalent to finding a linear, bounded and bijective operator T : H 1 0 pΩq Ñ H 1 0 pΩq such that the bilinear form ap¨, T¨q is coercive (for classical well-posedness) or weakly coercive (for the Fredholm case). The operator T plays the role of an explicit inf ´sup operator. However, proving T-coercivity can be difficult for complex geometries, especially in 3D. More recently, another viewpoint has come out that consists in studying the "limit" as δ Ñ 0 `of the well-posed stabilized problem: Find ũδ P H 1 0 pΩq such that a δ pũ δ , ϕq :" apũ δ , ϕq ´iδpK∇ũ δ , ∇ϕq Ωi " pf, ϕq Ω @ϕ P H 1 0 pΩq,

where H 1 0 pΩq is complex-valued. This approach has been developed in the context of Helmholtz equations by Nguyen [START_REF] Nguyen | Asymptotic behavior of solutions to the Helmholtz equations with sign-changing coefficients[END_REF][START_REF] Nguyen | Limiting absorption principle and well-posedness for the Helmholtz equation with sign-changing coefficients[END_REF], and makes use of transport operators called reflections.

Different strategies have been studied by Chesnel and Ciarlet Jr. in [START_REF] Chesnel | T-coercivity and continuous Galerkin methods: application to transmission problems with sign-changing coefficients[END_REF] for the numerical approximation of Problem (2) (assuming that it is well-posed in the classical Hadamard sense). The first approach is based on a simplicial discretization T h of Ω, that respects the interface Γ and the construction of a conforming finite element space V 0 pT h q that is stable by the operator T. Well-posedness and optimal convergence rates can be shown for this approach. In practice, T-stability is achieved by means of symmetric meshes near the interface, whose construction is a nontrivial task for complicated interfaces (see [START_REF] Bonnet-Ben Dhia | Mesh requirements for the finite element approximation of problems with sign-changing coefficients[END_REF]) or 3D problems. On general meshes, two main approaches have been investigated by Chesnel and Ciarlet Jr. The first one consists in building a (mesh-dependent) operator T h , such that the bilinear form a is T h -coercive on V 0 pT h q (see, e.g., [START_REF] Nicaise | A posteriori error estimates for a finite element approximation of transmission problems with sign-changing coefficients[END_REF]). This kind of approach is limited by the fact that, in general, well-posedness cannot be proved for the whole range of admissible coefficients. The second idea consists in discretizing the stabilized equation [START_REF] Bonnet-Ben Dhia | Mesh requirements for the finite element approximation of problems with sign-changing coefficients[END_REF], and in scaling the dissipation δ as a function of h. However, and as studied in [START_REF] Chesnel | T-coercivity and continuous Galerkin methods: application to transmission problems with sign-changing coefficients[END_REF], this approach leads to suboptimal convergence rates.

In this note, we aim at proposing a new discretization approach, that is based on a reformulation of Problem (1) as a transmission problem, whose solution is obtained numerically from an optimization procedure. The approach we propose is proved to converge, and does not rely on any symmetry assumption on the mesh. As opposed to T-coercivity theory, we do not need to assume well-posedness of Problem (2) in the Hadamard sense (i.e. for any f P H ´1pΩq) to prove convergence of our numerical method. We only need to assume that, for the given f P L 2 pΩq we consider, the solution to Problem (2) exists and is unique. This allows us to treat cases that cannot be analyzed using the T-coercivity theory (cf. Remark 3). The numerical method we consider here has been introduced by Gunzburger et al. [START_REF] Gunzburger | An optimization-based domain decomposition method for partial differential equations[END_REF] (see also [START_REF] Gunzburger | Solution of elliptic partial differential equations by an optimization-based domain decomposition method[END_REF]), in the context of domain decomposition for classical elliptic equations. However, up to now, the convergence of this method has never been proved. We provide the first proof of convergence of such a numerical method, with application to sign-changing elliptic equations, for which this approach seems particularly promising.

Transmission problem and numerical method

We provide in this section, under suitable regularity assumptions, an alternative characterization of the solution to Problem (2), that will be our starting point for the design of the numerical method.

Transmission problem

Recall that f P L 2 pΩq. Hence, s K∇ũ P Hpdiv; Ωq. Denoting, for α P te, iu, by n α the unit normal vector to Γ pointing out of Ω α , and introducing g such that g :"

" pK∇ũq |Ωe ¨ne ı |Γ , (4)
we have that g P H ´1{2 pΓq. Here, H ´1{2 pΓq denotes the dual of H 1 {2 pBΩ i q when Γ " BΩ i or of H 1 {2 00 pΓq otherwise. For α P te, iu, and for any g P H ´1{2 pΓq, we consider in the subdomain Ω α the problem

s α pK∇u α pgq, ∇ϕ α q Ωα " pf, ϕ α q Ωα `sα xg, ϕ α y Γ @ϕ α P H 1 0zΓ pΩ α q, (5)
where H 1 0zΓ pΩ α q is the space of functions in H 1 pΩ α q that vanish on BΩ α zΓ.

In Ω e , Problem (5) always admits a unique solution u e pgq P H 1 0zΓ pΩ e q, as the measure of BΩ e X BΩ is nonzero. In Ω i , Problem (5) also admits a unique solution u i pgq P H 1 0zΓ pΩ i q if the measure of BΩ i X BΩ is nonzero. Otherwise, the problem in Ω i is purely Neumann and we assume that the flux g P H ´1{2 pΓq satisfies pf, 1q Ωi ´xg, 1y Γ " 0 to ensure that Problem (5) admits a solution, that is unique up to an additive constant. We fix the constant by imposing pu i pgq, 1q Γ " pu e pgq, 1q Γ . Finally, for g P H ´1{2 pΓq, we denote by upgq the function such that upgq |Ωα :" u α pgq, α P te, iu. Proposition 2.1 (Characterization of the solution to Problem (2)) We assume that (2) admits a unique solution ũ P H 1 0 pΩq, and that there exists s ą 1 2 such that ũ|Ωα P H 1`s pΩ α q for α P te, iu. Then, ' the flux g defined in (4) belongs to L 2 pΓq, and satisfies g " ´"p´K∇ũq |Ωi ¨ni ı |Γ ;

' upgq " ũ P H 1 0 pΩq; ' almost everywhere on the interface Γ, u e pgq " u i pgq;

' the problem inf gPL 2 pΓq }u e pgq ´ui pgq} 2 0,Γ admits g as its unique solution.

Remark 1 (Assumptions on ũ) The existence (and uniqueness) of a solution ũ P H 1 0 pΩq to Problem (2) is satisfied in practice in a large variety of situations (cf., e.g., [START_REF] Nguyen | Limiting absorption principle and well-posedness for the Helmholtz equation with sign-changing coefficients[END_REF]). The assumption ũ|Ωα P H 1`s pΩ α q for α P te, iu, s ą 1 2 , made in Proposition 2.1 is convenient for the analysis, as it enables to work in L 2 pΓq instead of working in H ´1{2 pΓq. This theoretical assumption is quite strict, but not mandatory in practice for the method to be applicable. In a forthcoming work [START_REF] Abdulle | An optimization-based method for sign-changing PDEs[END_REF], we will consider a test-case for which this regularity assumption is violated, and show that numerical convergence can still be observed.

Minimization problem and numerical method

We consider a family of simplicial conformal discretizations tT h u h of Ω, that respects the interface Γ and is shape-regular in the sense of Ciarlet [START_REF] Ciarlet | The finite element method for elliptic problems[END_REF]. The subscript h stands for the meshsize, i.e. the maximum diameter of all the simplices in T h . We denote by Γ h the set of faces of the mesh T h belonging to Γ. For k P N ‹ and α P te, iu, we introduce the space V k pT h,α q :" v P H 1 pΩ α q | v |T P P k d pT q, @T P T h,α ( , where P k d pT q is the space of d-variate polynomial functions of total degree less or equal to k in T P T h . We also introduce the subspace V k 0zΓ pT h,α q :" V k pT h,α q X H 1 0zΓ pΩ α q. For α P te, iu, and for any g P H ´1{2 pΓq, we consider, in the subdomain Ω α , the following conforming approximation of Problem (5):

s α pK∇u h,α pgq, ∇ϕ h,α q Ωα " pf, ϕ h,α q Ωα `sα xg, ϕ h,α y Γ @ϕ h,α P V k 0zΓ pT h,α q. ( 6
)
As in the continuous case, Problem (6) always admits a unique solution u h,e pgq P V k 0zΓ pT h,e q in Ω e , and u h,i pgq P V k 0zΓ pT h,i q in Ω i when the measure of BΩ i X BΩ is nonzero. For the case of a pure Neumann problem in Ω i , we assume that pf, 1q Ωi ´xg, 1y Γ " 0, ensuring existence of a solution and uniqueness up to an additive constant, that we fix imposing pu h,i pgq, 1q Γ " pu h,e pgq, 1q Γ . For g P H ´1{2 pΓq, we introduce u h pgq such that u h pgq |Ωα :" u h,α pgq, α P te, iu.

To define the minimization problem, we introduce F k pΓ h q :" q P L 2 pΓq | q |F P P k d´1 pF q, @F P Γ h ( when the measure of BΩ i X BΩ is nonzero, and we add the constraint pq, 1q Γ " pf, 1q Ωi otherwise. We then define the functional J h : F k pΓ h q Ñ R `such that J h pg h q :" }u h,e pg h q ´uh,i pg h q} 2 0,Γ `λphq}g h } 2 0,Γ , where λ : R ‹ `Ñ R ‹ `is a function such that lim hÑ0 λphq " 0, and we consider the minimization problem inf

g h PF k pΓ h q J h pg h q. (7)
The continuous function J h is coercive on F k pΓ h q, and hence admits at least one minimizer: there exists gh P F k pΓ h q such that J h pg h q ď J h pg h q for any g h P F k pΓ h q.

The approximation of ũ we then consider is u h pg h q. The following result is a direct consequence of (8) and of the properties of the L 2 -orthogonal projector from L 2 pΓq onto F k pΓ h q.

Lemma 2.1 (Estimate on J h pg h q) Assume that g defined in (4) is in L 2 pΓq. Then, the following holds:

J h pg h q ď }u h,e pgq ´uh,i pgq}

2 0,Γ `λphq}g} 2 0,Γ . (9)

Convergence

From now on, we write A À B when there exists a constant c ą 0, possibly depending on K, k, and on the geometry, but independent of h, and of ũ, g, such that A ď cB. Lemma 3.1 (Approximation properties) Suppose that the assumptions of Proposition 2.1 are fulfilled. Then, letting p " minpk, sq, there exists δ P `0,

for α P te, iu and with δ m α P p0, 1s only depending on the geometry of Ω α . By setting δ " min

´δm e 2 , δ m i 2 ¯,
we obtain the stated result. We next treat the case of a pure Neumann problem in Ω i . Due to the choice pu h,i pgq, 1q Γ " pu h,e pgq, 1q Γ to fix the constant, a straightforward Aubin-Nitsche duality argument fails. We therefore consider an auxiliary function u h,i pgq, defined as the unique discrete solution of the same Neumann problem in Ω i , but for which we fix the constant as in the continuous problem, namely pu h,i pgq, 1q Ωi " pũ, 1q Ωi . Then, there exists δ n i P p0, 1s, only depending on the geometry of Ω i , such that

}u h,i pgq ´ũ} 0,Ωi `h δ n i 2 }u h,i pgq ´ũ} 0,Γ `hδ n i }∇pu h,i pgq ´ũq} 0,Ωi À h p`δ n i |ũ| 1`p,Ωi . (11)
Note that, a priori, δ n i ‰ δ m i since regularity results for pure Neumann or mixed problems are different in general. Next, we observe that u h,i pgq ´uh,i pgq " c h , for a constant c h P R. We then write |Γ| c h " pu h,i pgq ´uh,i pgq, 1q Γ " pu h,e pgq ´uh,i pgq, 1q Γ " pu h,e pgq ´ũ, 1q Γ `pũ ´uh,i pgq, 1q Γ , which yields, using [START_REF] Nguyen | Negative index materials and their applications: recent mathematics progress[END_REF] and [START_REF] Nguyen | Asymptotic behavior of solutions to the Helmholtz equations with sign-changing coefficients[END_REF],

|c h | ď |Γ| ´1{2 ´}u h,e pgq ´ũ} 0,Γ `}u h,i pgq ´ũ} 0,Γ ¯À h p`δ m e 2 |ũ| 1`p,Ωe `hp`δ n i 2 |ũ| 1`p,Ωi . (12)
We finally obtain from [START_REF] Nguyen | Limiting absorption principle and well-posedness for the Helmholtz equation with sign-changing coefficients[END_REF] }u h,i pgq ´ũ} 0,Γ ď }u h,e pgq ´ũ} 0,Γ `2}u h,i pgq ´ũ} 0,Γ À h p`δ We note that, as a consequence of Lemma 3.1, it is always possible to choose λphq in [START_REF] Ciarlet | The finite element method for elliptic problems[END_REF] such that }u h,e pgq ´uh,i pgq} 2 0,Γ " O pλphqq or }u h,e pgq ´uh,i pgq} 2 0,Γ " o pλphqq. We then have the Theorem 3.2 (Convergence of the method) Suppose that the assumptions of Proposition 2.1 are fulfilled. Denote u h pg h q by ũh . Then, up to a choice of λphq in [START_REF] Ciarlet | The finite element method for elliptic problems[END_REF] such that }u h,e pgq ´uh,i pgq} 2 0,Γ " O pλphqq, there holds as h Ñ 0:

gh á g in L 2 pΓq, ∇ h ũh á ∇ũ in L 2 pΩq d , ũh Ñ ũ in L 2 pΩq. (13)
If we further choose λphq such that }u h,e pgq ´uh,i pgq} 2 0,Γ " o pλphqq, then we have

gh Ñ g in L 2 pΓq, ∇ h ũh Ñ ∇ũ in L 2 pΩq d , ũh Ñ ũ in L 2 pΩq. (14)
Proof We begin by remarking that, for α P te, iu,

u h,α pg h q " u h,α pgq `vh,α , (15)
where

v h,α P V k 0zΓ pT h,α q satisfies pK∇v h,α , ∇ϕ h,α q Ωα " ppg h ´gq , ϕ h,α q Γ @ϕ h,α P V k 0zΓ pT h,α q. (16)
From ( 16), we readily infer

}∇v h,α } 2 0,Ωα À }g h ´g} 0,Γ }v h,α } 0,Γ . (17)
Besides, the following estimate can be obtained from (9):

}g h } 2 0,Γ ď }u h,e pgq ´uh,i pgq} 2 0,Γ λphq `}g} 2 0,Γ . (18)
As soon as }u h,e pgq ´uh,i pgq} 2 0,Γ {λphq ď c, the combination of ( 18), (17), and of a continuous trace inequality (as well as of an appropriate Poincaré inequality [5, Equation (1.1)] when the measure of BΩ i X BΩ is zero) enables to infer that there exist g0 P L 2 pΓq, and v α P H 1 0zΓ pΩ α q, α P te, iu, so that the following convergences hold as h Ñ 0, up to a subsequence (retaining the same notation):

gh á g0 in L 2 pΓq, ∇v h,α á ∇v α in L 2 pΩ α q d , v h,α Ñ v α in L 2 pΩ α q, v h,α á v α in L 2 pΓq. (19)
Now, from (9) we infer }u h,e pg h q ´uh,i pg h q} 2 0,Γ ď }u h,e pgq ´uh,i pgq}

2 0,Γ `λphq}g} 2 0,Γ ď ´c `}g} 2 0,Γ ¯λphq. (20)
Owing to the fact that λphq tends to zero as h vanishes, we deduce from (20) that }u h,e pg h q ´uh,i pg h q} 0,Γ Ñ 0 as h Ñ 0. Combining this last result with (15) and Lemma 3.1, we get that }v h,e ´vh,i } 0,Γ Ñ 0 as h Ñ 0 and hence, owing to (19), that v e " v i almost everywhere on Γ. Passing to the limit in ( 16) using ( 19) and a strongly converging interpolant for test functions, it can be shown that v P H 1 0 pΩq such that v |Ωα :" v α for α P te, iu satisfies ps K∇v, ∇ϕq Ω " 0 for all ϕ P H 1 0 pΩq, which implies, from the well-posedness of Problem [START_REF] Bernardi | Discrétisations variationnelles de problèmes aux limites elliptiques[END_REF], that v " 0 and g0 " g (this last result is inferred considering the limit equation in a subdomain, and using a density argument along with the fact that g0 , g P L 2 pΓq). The uniqueness of the limits implies that the whole sequences converge in (19). Collecting these last results, (15), and Lemma 3.1, we prove [START_REF] Nicaise | A posteriori error estimates for a finite element approximation of transmission problems with sign-changing coefficients[END_REF]. If we further choose λphq such that }u h,e pgq ´uh,i pgq} 2 0,Γ {λphq Ñ 0 as h Ñ 0, we get from (18) and from the weak convergence of gh towards g that

}g} 0,Γ ď lim hÑ0 inf }g h } 0,Γ " lim hÑ0 sup }g h } 0,Γ ď }g} 0,Γ , (21)
which states the strong convergence of gh towards g in L 2 pΓq. Testing (16) with v h,α , and using (21) combined with the weak convergence result for v h,α in L 2 pΓq of (19), we finally get the strong convergence of ∇v h,α towards 0 in L 2 pΩ α q d for α P te, iu, proving (14) and concluding the proof. l

Numerical validation

We consider the 2D symmetric cavity test-case of [START_REF] Chesnel | T-coercivity and continuous Galerkin methods: application to transmission problems with sign-changing coefficients[END_REF], for which Ω :" p´1, 1q ˆp0, 1q and Ω e :" p´1, 0qˆp0, 1q (it follows that Ω i " p0, 1qˆp0, 1q and Γ " t0uˆr0, 1s). For this test-case, both the exterior and the interior problems feature mixed Dirichlet-Neumann boundary conditions. For the particular geometry considered here, elliptic regularity holds in both subdomains (see, e.g., [START_REF] Bernardi | Discrétisations variationnelles de problèmes aux limites elliptiques[END_REF]Remark I.3.6]), so that δ " 1 2 in Lemma 3.1. The tensor K is taken isotropic, with constant value in each subdomain (k |Ωe " k e ą 0 and k |Ωi " k i ą 0). For this particular setting, it is known [START_REF] Chesnel | T-coercivity and continuous Galerkin methods: application to transmission problems with sign-changing coefficients[END_REF] and is associated to the right-hand side f " ´divps k∇ũq P L 2 pΩq. It can be easily seen that, for α P te, iu, ũ|Ωα P H 1`l pΩ α q for any l ą 0, meaning that p " k in Lemma 3.1.

As in [START_REF] Chesnel | T-coercivity and continuous Galerkin methods: application to transmission problems with sign-changing coefficients[END_REF], we choose k e " 1 and k i " 1.001, so that ν " ´1.001. We run the computations on a sequence of non-symmetric (with respect to the interface Γ), unstructured meshes. Results are depicted on Figure 1 for linear pk " 1q and second-order pk " 2q FEM, for our approach (solid blue), and for classical FEM (dashed black) applied to (2). For our approach, the parameter λphq is chosen as Oph β q with βďp2k`1q, as required by Theorem 3.2 and Lemma 3.1. Here we choose λphq " 0.002 h 3 for k " 1 and λphq " 0.002 h 4.2 for k " 2. We first observe that all convergence plots for our approach are strictly monotone, as opposed to classical FEM. In the L 2 -norm, for h sufficiently small, we observe a slightly super-convergent behavior for our approach for both k " 1 and k " 2. In the H 1 -norm, for both orders, our method seems to reach the expected convergence rates for h sufficiently small, and clearly outperforms classical FEM. The choice of a small multiplicative coefficient 0.002 in the function λ is guided by the fact that the norm of g is big for such a contrast. Hence, for coarse meshes for which convergence is far from being reached, one has to give less weight to the second term of the functional J h in the minimization process, to have a chance to recover a correct approximation of g. This weight is useless when h is sufficiently small. Remark 2 (Stabilization parameter) For coarse meshes, one has to tune the parameter λphq in a nontrivial, contrast-depending way. This has to do with the regularization we use. If the L 2 -orthogonal projection Π k h g of g onto F k pΓ h q was known, then we would stabilize the functional J h with µphq Remark 3 It has to be noted that the symmetric cavity test-case can be analyzed by means of T-coercivity.

However, it cannot be on general meshes, which is the main outcome of our approach. Furthermore, there are also many cases (even in 2D) that cannot be analyzed, at the continuous level, using the T-coercivity approach. This is for example the case of the cloaking device studied in [10, Section 4]. In that case, the operator associated to the problem, viewed as an operator from H 1 0 pΩq to H ´1pΩq is not Fredholm, as its range is not closed. T-coercivity is hence inapplicable. However, it is known that, for compatible righthand sides f P L 2 pΩq, the solution exists and is unique in H 1 0 pΩq [START_REF] Nguyen | Negative index materials and their applications: recent mathematics progress[END_REF]. The convergence of our numerical method for compatible loadings as well as numerical experiments for such problems will be investigated in a forthcoming work [START_REF] Abdulle | An optimization-based method for sign-changing PDEs[END_REF].

2

  |ũ| 1`p,Ωi , which enables to conclude the proof, setting in that case δ " min ´δm

instead of λphq}g h } 2 0

 2 ,Γ , yielding, in view of (18), an optimal estimate on

  Proof We first assume that we have, both in Ω e and Ω i , mixed Dirichlet-Neumann boundary conditions.

	Then, classical arguments (particularly, Aubin-Nitsche duality argument, and a multiplicative continuous
	trace inequality) yield
	}u h,α pgq ´ũ} 0,Ωα	`h δ m α 2 }u h,α pgq ´ũ} 0,Γ	`hδ m α }∇pu h,α pgq ´ũq} 0,Ωα À h p`δ m α |ũ| 1`p,Ωα ,
					1 2	‰	such that
			}∇ ¯,
	}u h,e pgq ´uh,i pgq} 0,Γ ď	ÿ	}u h,α pgq ´ũ} 0,Γ À h p`δ ´|ũ| 1`p,Ωe `|ũ| 1`p,Ωi	¯.
			α Pte, iu

h pu h pgq ´ũq} 0,Ω À h p ´|ũ| 1`p,Ωe `|ũ| 1`p,Ωi

  that Problem (2) is wellposed in the classical Hadamard sense if and only if the contrast ν :" ´ki ke is different from ´1. Figure 1. Relative errors on Ω vs. meshsize for the symmetric cavity test-case with ν " ´1.001.

			10 0
			10 ´1				10 ´1
			10 ´2
							10 ´3
			10 ´3
			10 ´4				10 ´5
			10 ´5			k " 1			k " 1
			10 ´6			k " 2	10 ´7		k " 2
			10 ´3	10 ´2	10	´1	10 ´3	10 ´2	10	´1
	(a) Member of the mesh sequence		(b) Relative H 1 -norm error vs. h	(c) Relative L 2 -norm error vs. h
									The exact
	solution ũ P H 1 0 pΩq we consider is
	ũpx, yq :"	$ ' ' & ' ' %	ˆpx `1q k e k e ´ki px ´1q sinpπyq 2 ´2k e ´ki px `1q ˙sinpπyq k e ´ki	in Ω e , in Ω i ,