Claudia Fonte

Valentin Schmutz

Long time behavior of an age and leaky memory-structured neuronal population equation

published or not. The documents may come

Long time behavior of an age and leaky memory-structured neuronal population equation

Introduction

Multidimensional mean-field models in theoretical neuroscience are notoriously challenging to analyse [START_REF] Richardson | Dynamics of populations and networks of neurons with voltageactivated and calcium-activated currents[END_REF][START_REF] Vellmer | Theory of spike-train power spectra for multidimensional integrate-and-fire neurons[END_REF][START_REF] Beiran | Contrasting the effects of adaptation and synaptic filtering on the timescales of dynamics in recurrent networks[END_REF][START_REF] Samuel P Muscinelli | How single neuron properties shape chaotic dynamics and signal transmission in random neural networks[END_REF] but their study is a necessary step towards understanding how the multiple timescales present at the single neuron level [START_REF] Pozzorini | Temporal whitening by power-law adaptation in neocortical neurons[END_REF][START_REF] Teeter | Generalized leaky integrate-and-fire models classify multiple neuron types[END_REF] affect the dynamics of large networks of neurons.

One-dimensional mean-field population equations for spiking neurons with 'escape noise' [START_REF] Gerstner | Neuronal dynamics: From single neurons to networks and models of cognition[END_REF]Ch. 9] have been a subject of mathematical study since the seminal papers of Pakdaman, Perthame and Salort [START_REF] Khashayar Pakdaman | Dynamics of a structured neuron population[END_REF][START_REF] Khashayar Pakdaman | Relaxation and self-sustained oscillations in the time elapsed neuron network model[END_REF][START_REF] Khashayar Pakdaman | Adaptation and fatigue model for neuron networks and large time asymptotics in a nonlinear fragmentation equation[END_REF], providing rigorous foundations to earlier works in theoretical neuroscience [START_REF] Hugh | Excitatory and inhibitory interactions in localized populations of model neurons[END_REF][START_REF] Gerstner | Associative memory in a network of 'spiking' neurons[END_REF][START_REF] Gerstner | Time structure of the activity in neural network models[END_REF][START_REF] Gerstner | Population dynamics of spiking neurons: fast transients, asynchronous states, and locking[END_REF]. The fact that these population equations correspond to the mean-field limit of large networks of interacting neurons has been proved in [START_REF] De Masi | Hydrodynamic limit for interacting neurons[END_REF][START_REF] Fournier | On a toy model of interacting neurons[END_REF][START_REF] Chevallier | Mean-field limit of generalized Hawkes processes[END_REF]. All these population equations are derived from spiking neuron models that are of the 'renewal' type (with the exception of the model with fatigue of [START_REF] Khashayar Pakdaman | Adaptation and fatigue model for neuron networks and large time asymptotics in a nonlinear fragmentation equation[END_REF]), which means that they can capture the effect of neuronal refractoriness but neglect slower neuronal timescales, like that of spike frequency adaptation and short-term synaptic plasticity.

To take into account slow neuronal timescales, state-of-the-art phenomenological spiking neuron models are multidimensional [START_REF] Kobayashi | Made-to-order spiking neuron model equipped with a multi-timescale adaptive threshold[END_REF][START_REF] Teeter | Generalized leaky integrate-and-fire models classify multiple neuron types[END_REF] or kernel-based [START_REF] Truccolo | A point process framework for relating neural spiking activity to spiking history, neural ensemble, and extrinsic covariate effects[END_REF][START_REF] Pillow | Spatio-temporal correlations and visual signalling in a complete neuronal population[END_REF][START_REF] Pozzorini | Automated high-throughput characterization of single neurons by means of simplified spiking models[END_REF] (and see [START_REF] Gerstner | Neuronal dynamics: From single neurons to networks and models of cognition[END_REF]Ch. 6.4]). Recently, one of the author showed that, for a large class of multidimensional spiking neuron models, the mean-field limit can be characterized by an 'age and leaky memory'-structured evolution PDE in 1 + d spatial dimensions, d being the number of 'leaky memory' variables [START_REF] Schmutz | Mean-field limit of Age and Leaky memory dependent Hawkes processes[END_REF].

The aim of this work is to study the long time behavior of the 1 + 1 dimensional case (i.e. one age and one leaky memory variable) of the equation proposed in [START_REF] Schmutz | Mean-field limit of Age and Leaky memory dependent Hawkes processes[END_REF].

The age-and leaky memory-structured model

The population model we consider describes the evolution of a density ρ t over the state-space (a, m) ∈ R + × R * + , where a and m are the 'age' and 'leaky memory' variables of the neuron, and ρ t (a, m) represents the density of neurons in state (a, m) at time t.

The nonlinear evolution problem for the density ρ t , for the initial datum u 0 , writes

∂ t ρ t + ∇ • (bρ t ) = -f (a, m, εx t )ρ t , (1a)
ρ t (0, m) = p(t, γ -1 (m)), (1b)

x t = t 0 ∞ 0 ∞ 0 h(t -s, a, m)f (a, m, εx s )ρ s (a, m)dadmds, (1c)
ρ 0 = u 0 . (1d)
The dynamics of the model can be decomposed in three elements: the behavior of neurons between spikes, the spike-triggered jumps and the interaction between neurons.

Between spikes, neurons are transported along the vector field b(a, m) = (1, -λm), with λ > 0 (∇• denotes the divergence operator over the state-space).

Neurons spike at a rate f (a, m, εx t ), where f : R + × R * + × R → R + is the 'firing rate function'. When a neuron spikes, its age a is reset to 0 and its leaky memory variable m jumps to γ(m), where γ : R + → R * + is the 'jump mapping' and is assumed to be a strictly increasing C 1 -diffeomorphism. As a consequence, the border condition (1b) has a simple interpretation: the density of neurons in state (0, m) at time t is equal to the the marginal density of neurons which have their leaky memory variable in state γ -1 (m) and spike at time t. This quantity is denoted p(t, γ -1 (m)) and has the explicit expression p(t, γ -1 (m)) := ½ m>γ(0) (γ -1 ) ′ (m) ∞ 0 f (a, γ -1 (m), εx t )ρ t (a, γ -1 (m))da.

(1e)

The indicator function ½ m>γ(0) reflects the fact that m is always strictly positive and the term (γ -1 ) ′ (m) is necessary to guarantee the conservation of the total mass of neurons. Indeed, formally,

∂ t ρ t = ½ m>γ(0) (γ -1 ) ′ (m) ∞ 0
f (a, γ -1 (m), εx t )ρ t (a, γ -1 (m))dadm -f (a, m, εx t )ρ t = 0, by a change of variable. Lastly, neurons interact through the 'total postsynaptic potential' x t , which integrates the past spiking activity of the population, filtered by the 'interaction function' h : R + × R + × R * + → R. x t , weighted by the connection strength ε ∈ R, influences the firing rate f . If we write N (t) the mean firing rate

N (t) := ∞ 0 ∞ 0 f (a, m, εx t )ρ t (a, m)dadm,
and if we take h independent of a and m, then x t takes the form

x t = t 0 h(t -s)N (s)ds,
and h is simply a delay kernel, as in the age-structured model of [START_REF] Khashayar Pakdaman | Dynamics of a structured neuron population[END_REF]. In our formulation, h in Eq. (1c) allows to model more general interactions. For example, in Sec. 1.2.2, we show that by choosing h(t, a, m) = ĥ(t)(1 -m), we can include the effects of a classical short-term synaptic plasticity model [START_REF] Tsodyks | Neural networks with dynamic synapses[END_REF].

Motivation

The model [START_REF] Beiran | Contrasting the effects of adaptation and synaptic filtering on the timescales of dynamics in recurrent networks[END_REF] extends the age-structured model of Pakdaman, Perthame and Salort [START_REF] Khashayar Pakdaman | Dynamics of a structured neuron population[END_REF] by the addition of a leaky memory variable which can accumulate over spikes (as opposed to the age variable which is reset to 0 at each spike) and hence introduces a slow timescale in the population dynamics. Such a slow timescale is typically used to account for some form of fatigue mechanism, which can act on the spiking activity (spike frequency adaptation) or on synaptic transmission (short-term synaptic depression). Slow fatigue at the single neuron level can lead to nontrivial emergent behaviors at the population level, like population bursts [START_REF] Van Vreeswijk | Patterns of synchrony in neural networks with spike adaptation[END_REF][START_REF] Gigante | Diverse population-bursting modes of adapting spiking neurons[END_REF][START_REF] Gast | A mean-field description of bursting dynamics in spiking neural networks with short-term adaptation[END_REF] (see Fig. 1), which are not observed in the age-or voltage-structured models of [START_REF] Khashayar Pakdaman | Dynamics of a structured neuron population[END_REF] and [START_REF] De Masi | Hydrodynamic limit for interacting neurons[END_REF]. Even though some population equations have been successfully used in the computational neuroscience literature to study emergent behaviors in networks of neurons with fatigue [START_REF] Gigante | Diverse population-bursting modes of adapting spiking neurons[END_REF][START_REF] Gast | A mean-field description of bursting dynamics in spiking neural networks with short-term adaptation[END_REF], these population equations were obtained at the cost of a timescale separation approximation, making them inexact. In contrast, the model ( 1) is an exact mean-field limit [START_REF] Schmutz | Mean-field limit of Age and Leaky memory dependent Hawkes processes[END_REF].

This model is also directly related to standard models of spike frequency adaptation and shortterm synaptic depression.

Spike frequency adaptation

The recent spike history of a neuron can modulate its firing rate f , leading to spike frequency adaptation [START_REF] Benda | A universal model for spike-frequency adaptation[END_REF]. If h is independent of a and m and if γ(m) = m + Γ, for a fixed Γ>0, (1) becomes

∂ t ρ t + ∇ • (bρ t ) = -f (a, m, εx t )ρ t , (2a)
ρ t (0, m) = ½ m> Γ ∞ 0 f (a, m -Γ, εx t )ρ t (a, m -Γ)da, (2b)
x t = t 0 h(t -s) ∞ 0 ∞ 0 f (a, m, εx s )ρ s (a, m)dadmds, (2c)
ρ 0 = u 0 . (2d)
If η : R + → R is a bounded function such that lim a→+∞ η(a) = 0 (η is the 'refractory kernel' [16, Sec. 9.3]), we can define f more explicitly:

f (a, m, x) := f (η(a) -m + x), (2e)
where f : R → R + is typically a non-decreasing function. Since m makes jumps of size Γ > 0 at each spike and decays at exponential speed -λ between spikes, m accumulates over spikes, which decreases the firing rate f (Eq. (2e)), leading to spike frequency adaptation [START_REF] Benda | A universal model for spike-frequency adaptation[END_REF]. More specifically, Eq. ( 2) is a population equation for adaptive SRM 0 (Spike Response Model) neurons [START_REF] Jolivet | Predicting spike timing of neocortical pyramidal neurons by simple threshold models[END_REF][START_REF] Gerstner | Neuronal dynamics: From single neurons to networks and models of cognition[END_REF].

Populations of spiking neurons with spike frequency adaptation exhibit self-sustained population bursts when the connectivity strength is sufficiently strong [START_REF] Van Vreeswijk | Patterns of synchrony in neural networks with spike adaptation[END_REF][START_REF] Gigante | Diverse population-bursting modes of adapting spiking neurons[END_REF][START_REF] Gast | A mean-field description of bursting dynamics in spiking neural networks with short-term adaptation[END_REF]. In Fig. 1, we show simulations of (2) for two different connectivity strengths ε. For large ε, we observe self-sustained bursts, whereas for small ε, we observe relaxation to a stationary state.

Short-term synaptic depression

The recent spike history of a presynaptic neuron can modulate the synaptic transmission, leading to short-term synaptic plasticity [START_REF] Robert | Short-term synaptic plasticity[END_REF]. We will consider here the case of depressive synapses and use the model of [START_REF] Tsodyks | Neural networks with dynamic synapses[END_REF] (with a change of variable for convenience). In this case, the state-space is (a, m) ∈ R + ×]0, 1[. Taking f independent of m, h of the the form h(t, a, m) := ĥ(t)(1 -m) and γ(m) := 1 -υ + υm for a fixed υ ∈]0, 1[, (1) becomes

∂ t ρ t + ∇ • (bρ t ) = -f (a, m, εx t )ρ t , (3a)
ρ t (0, m) = ½ m>γ(0) 1 υ ∞ 0 f (a, εx t )ρ t (a, γ -1 (m))da, (3b)
x t = t 0 ĥ(t -s) 1 0 ∞ 0 (1 -m)f (a, εx s )ρ s (a, m)dadmds, (3c)
ρ 0 = u 0 . ( 3d
)
Note that the term 1 υ on the RHS of Eq. (3b) simply comes from the fact that |(γ -1 ) ′ (m)|= 1 υ , for all m ∈]0, 1[. Here, at each spike, m makes strictly positive jumps which size tends to 0 as m tends to 1 (γ(1) = 1) and decays at exponential speed -λ between spikes. If m is close to 1, synaptic transmission is weak because of the factor (1 -m) in Eq. (3c). This is an implementation of the purely depressive Tsodyks-Markram model [START_REF] Tsodyks | Neural networks with dynamic synapses[END_REF], a standard model of short-term synaptic depression.

As observed in [START_REF] Romani | Mean-field analysis of selective persistent activity in presence of short-term synaptic depression[END_REF], the stationary state of populations of neurons with short-term synaptic plasticity can be described by a simple formula, which we prove in Sec. 4.3. (2) showing self-sustained bursts for strong connectivity strength (ε ≫ 0) and relaxation to a stationary state for weak connectivity strength (small ε). Simulations of a network of 5 • 10 5 adaptive SRM 0 neurons approximating Eq. ( 2) with identical parameters (except for ε) and identical initial conditions. The raster plots below the plots for the time-evolution of the total postsynaptic potential x t represent the spikes of 100 randomly selected neurons.

Assumptions and main results

The main result of this work is the exponential stability of the age and leaky memory-structured equation (1) in the weak connectivity regime (Theorem 3). Before proving the exponential stability, we first establish the well-posedness of (1) in the appropriate function space (Theorem 1) and show that stationary solutions exist and are unique for sufficiently weak connectivity (Theorem 2).

Here, we study the weak solutions to (1) for an initial datum in

L 1 + := L 1 (R + × R * + , R + ) and write L 1 + (R * + ) := L 1 (R * + , R + ). Definition (Solutions). (ρ, x) ∈ C(R + , L 1 + ) × C(R +
) is a solution to (1), for the initial datum

u 0 ∈ L 1 + , if x t = t 0 ∞ 0 ∞ 0 h(t -s, a, m)f (a, m, εx s )ρ s (a, m)dadmds, ∀t ≥ 0, (4a)
and if for all

ϕ ∈ C ∞ c (R + × R + × R * + ), 0 = ∞ 0 ∞ 0 u 0 (a, m)ϕ(0, a, m)dadm + ∞ 0 ∞ 0 ∞ 0 ρ t (a, m) [∂ t + ∂ a -λm∂ m ]ϕ + (ϕ(t, 0, γ(m)) -ϕ(t, a, m))f (a, m, εx t ) dadmdt. (4b)
To prove the well-posedness of (1), we need some simple assumptions of the firing rate function f and the interaction function h:

Assumption 1. f is bounded and L f -Lipschitz, i.e. |f (a, m, x) -f (a * , m * , x * )|≤ L f (|a -a * |+|m -m * |+|x -x * |),
and h is bounded and continuous.

Since our goal is the use Harris' theorem to study the long time behavior of (1), the wellposedness in L 1 (which is treated in [START_REF] Schmutz | Mean-field limit of Age and Leaky memory dependent Hawkes processes[END_REF]) is not enough and we need the well-posedness in a weighted L 1 space (where the weight satisfies a Lyapunov condition [START_REF] Sean | Stability of Markovian processes III: Foster-Lyapunov criteria for continuous-time processes[END_REF]) with a global-in-time estimate in the weighted L 1 norm.

Using the weight function

w : R + × R + → [1, ∞), (a, m) → 1 + m,
we define the function space

L 1 + (w) := g ∈ L 1 (R + × R * + , R + ) g L 1 (w) := ∞ 0 ∞ 0 g(a, m)w(a, m)dadm < ∞ .
To obtain a global-in-time estimate in the L 1 + (w) norm, we further need that the jump mapping γ grows at most linearly. For convenience, we formulate this assumption in terms of the jump size function Γ:

Assumption 2. There exists Γ ∈ C 1 b (R + + , R * + ) such that γ(m) = m + Γ(m), ∀m ∈ R * + .
Theorem 1 (Well-posedness). Grant Assumption 1. For any initial datum u 0 ∈ L 1 + , there exists a unique weak solution (ρ, x) to (1). This solution satisfies (i) (L 1 -stability)

ρ t L 1 = u 0 L 1 , ∀t > 0, (ii) (Global bound in L 1 + (w)) if, in addition, Assumption 2 holds and u 0 ∈ L 1 + (w), then ∀t > 0, ρ t L 1 (w) ≤ u 0 L 1 (w) e -αt + b α (1 -e -αt ), (5)
for some constants α > 0 and b ∈ R.

In contrast to [START_REF] Schmutz | Mean-field limit of Age and Leaky memory dependent Hawkes processes[END_REF], the well-posedness proof presented here does not involve any probabilistic argument. The proof consists of two consecutive applications of Banach's fixed-point theorem, where a first fixed-point gives the unique solution to a linearized version of (1) which is then used in a second fixed-point treating the nonlinearity of [START_REF] Beiran | Contrasting the effects of adaptation and synaptic filtering on the timescales of dynamics in recurrent networks[END_REF].

The second step towards the exponential stability proof is the study of the existence and uniqueness of the stationary solutions to [START_REF] Beiran | Contrasting the effects of adaptation and synaptic filtering on the timescales of dynamics in recurrent networks[END_REF]. For this step, we require: Assumption 3. nothing (i) There exists ∆ abs > 0 and σ > 0 such that

f (a, m, x) ≥ σ, ∀(a, m, x) ∈ [∆ abs , +∞[×R * + × R.
(ii) There exists

C γ ∈]0, 1] such that C γ ≤ γ ′ ≤ 1. (iii) h(a, m) = ∞ 0 h(t, a, m)dt is bounded.
The first point of Assumption 3 sets a lower bound on the firing rate function f for any a ≥ ∆ abs and hence allows for an absolute refractory period ∆ abs > 0, i.e. a period of time following a spike during which f = 0 (which is an important neurodynamical feature [START_REF] Gerstner | Neuronal dynamics: From single neurons to networks and models of cognition[END_REF]Sec. 1.1]). This assumption is also used in [START_REF] José | Asymptotic behaviour of neuron population models structured by elapsed-time[END_REF].

In the second point of Assumption 3, the lower bound 0 < C γ ≤ γ ′ guarantees that γ is strictly increasing, which reflects the idea that m is a 'leaky memory' variable of the past neuronal activity. On the other hand, the upper bound γ ′ ≤ 1, which can be rewritten in terms of the jump size function Γ as Γ ′ ≤ 0, prevents the variable m from growing too fast and allows for a potential saturation of the memory, as in the example with short-term synaptic plasticity [START_REF] José | Asymptotic behaviour of neuron population models structured by elapsed-time[END_REF]. The third point of Assumption 3 is purely technical.

We emphasize that the two examples shown above, spike frequency adaptation (2) and shortterm synaptic depression (3), satisfy Assumption 3.

Theorem 2 (Stationary solutions). Grant Assumptions 1 and 3.

(i) There exists a stationary solution to (1).

(ii) There exists ε * > 0 such that for all ε ∈] -ε * , +ε * [, the stationary solution to (1) is unique.

Over the course of this work, we obtained the existence of the stationary solution by two different approaches. The first approach is based on the Doeblin-Harris technique [START_REF] Hairer | Yet another look at Harris' ergodic theorem for Markov chains[END_REF] and is similar to that of [START_REF] José | Asymptotic behaviour of neuron population models structured by elapsed-time[END_REF]. First, we show that when x t is fixed and time-invariant in (1) (neurons are non-interacting), the system satisfies a Harris condition -this constitutes a key result of this work -, and we can use Harris' theorem to get the stationary solution. Then, we use the Lipschitz continuity the stationary solutions with respect to the fixed x to prove the existence of a stationary solution for arbitrary connectivity strengths ε. Finally, for ε small enough, we also get the uniqueness of the stationary solution, by Banach's fixed-point theorem.

The second approach relies on the fact that the stationary solutions solve an integral equation, for which we can show that a solution exists by the Schauder's fixed-point theorem. In the process, we get several estimates on the stationary solutions, namely that they are continuous, bounded, and exponentially decaying in m. However, this approach does give uniqueness.

As mentioned above, the application of Harris' Theorem requires us to consider solutions in the weighted space L 1 (w). However, in the case where the state-space of the leaky memory variable m is bounded, the situation is simpler: we can use Doeblin's theorem in L 1 . The following assumption guarantees that m stays in a bounded state-space: Assumption 4. There exists G > 0 such that for all m ∈ R + , γ(m) < G.

Note that this assumption is satisfied in the example with short-term synaptic plasticity (3), with G = 1.

Finally, to study the exponential stability of (1), we need on exponential decay on h:

Assumption 5. There exists h, C h > 0 such that h(t, a, m) ≤ C h e -ht , ∀(t, a, m).
By a perturbative argument similarly as in [START_REF] Mischler | Weak and strong connectivity regimes for a general time elapsed neuron network model[END_REF], we obtain our main result:

Theorem 3 (Exponential stability in the weak connectivity regime). Grant Assumptions 1, 3 and 5. For any W > 0, there exists

ε * * W > 0 such that for ε ∈] -ε * * W , +ε * * W [, there exists C ≥ 1 and c W > 0 such that for all initial data u 0 ∈ L 1 + (w) with u 0 L 1 = 1 and u 0 L 1 (w) ≤ W , the solution (ρ, x) to (1) satisfies ρ t -ρ ∞ L 1 (w) + |x t -x ∞ |≤ Ce -c W t u 0 -ρ ∞ L 1 (w) + 1 , ∀t ≥ 0, (6)
where (ρ ∞ , x ∞ ) is the unique stationary solution given by Theorem 2 (ii).

If, in addition, we grant Assumption 4, then there exists ε * * > 0 such that for all ε ∈]-ε * * , +ε * * [, there exists C ′ ≥ 1 and c > 0 such that for all initial data

u 0 ∈ L 1 + with u 0 L 1 = 1, ρ t -ρ ∞ L 1 + |x t -x ∞ |≤ C ′ e -ct u 0 -ρ ∞ L 1 + 1 , ∀t ≥ 0. (7)
From the neuronal modeling point of view, this result is not surprising: when the connection strength is weak enough, neurons do not synchronize and the population activity converges to a stationary state. This was already proved for simpler one-dimensional models (see below) and the addition of a leaky memory variable carrying the effect of spike frequency adaptation or short-term synaptic plasticity does not change this behavior.

Methodology discussion

The asymptotic stability of the age-structured model of [START_REF] Khashayar Pakdaman | Dynamics of a structured neuron population[END_REF] in the weak connectivity regime has been studied by the same authors in [START_REF] Khashayar Pakdaman | Dynamics of a structured neuron population[END_REF][START_REF] Khashayar Pakdaman | Relaxation and self-sustained oscillations in the time elapsed neuron network model[END_REF] using entropy methods and assuming that f is a step-function, in [START_REF] Mischler | Relaxation in time elapsed neuron network models in the weak connectivity regime[END_REF][START_REF] Mischler | Weak and strong connectivity regimes for a general time elapsed neuron network model[END_REF] using spectral analysis of semigroups in Banach spaces and in [3] using Doeblin's theorem. For the treatment of the strong connectivity regime, we refer to [START_REF] Khashayar Pakdaman | Dynamics of a structured neuron population[END_REF][START_REF] Mischler | Weak and strong connectivity regimes for a general time elapsed neuron network model[END_REF].

On the closely related voltage-structured model [START_REF] De Masi | Hydrodynamic limit for interacting neurons[END_REF], [START_REF] Cormier | Long time behavior of a mean-field model of interacting neurons[END_REF] also proved asymptotic stability in the weak connectivity regime using Laplace transform techniques. For this model, the nonlinear stability of the stationary solutions has been recently studied in [START_REF] Cormier | A mean-field model of integrate-and-fire neurons: non-linear stability of the stationary solutions[END_REF] (see also [START_REF] Drogoul | Exponential stability of the stationary distribution of a mean field of spiking neural network[END_REF]) and can identify Hopf bifurcations [START_REF] Cormier | Hopf bifurcation in a mean-field model of spiking neurons[END_REF]. This type of analysis would be challenging in the case of (1) in general, as there is no representation formula for the stationary solutions when ε = 0.

Finally, Doeblin's theorem has also been used in [START_REF] Dumont | The mean-field equation of a leaky integrate-and-fire neural network: measure solutions and steady states[END_REF] in the case of the 'threshold crossing' neuronal population equation of [START_REF] Omurtag | On the simulation of large populations of neurons[END_REF]. Note also that closely related methods have been used by probabilists to study the ergodicity of single neuron models [START_REF] Höpfner | Ergodicity and limit theorems for degenerate diffusions with time periodic drift. application to a stochastic Hodgkin-Huxley model[END_REF][START_REF] Duarte | Stability, convergence to equilibrium and simulation of non-linear Hawkes processes with memory kernels given by the sum of Erlang kernels[END_REF].

Our approach combines strategies from [START_REF] Mischler | Relaxation in time elapsed neuron network models in the weak connectivity regime[END_REF] and [START_REF] José | Asymptotic behaviour of neuron population models structured by elapsed-time[END_REF]. On the one hand, our proof is based on the application of Harris' theorem for the linear problem, like in [START_REF] José | Asymptotic behaviour of neuron population models structured by elapsed-time[END_REF] 1 , which simplifies the proof of [START_REF] Mischler | Relaxation in time elapsed neuron network models in the weak connectivity regime[END_REF]. On the other hand, we use an argument from in [START_REF] Mischler | Relaxation in time elapsed neuron network models in the weak connectivity regime[END_REF] to deal with delay effects, which are not considered in [START_REF] José | Asymptotic behaviour of neuron population models structured by elapsed-time[END_REF]. Note that our model is two-dimensional (by the addition of the leaky memory variable), whereas the aforementioned works only considered one-dimensional models.

Plan of the paper

The proof of Theorem 1 (Well-posedness) is presented in Section 2. In Section 3, we prove the exponential stability of (1) in the non-interacting case ε = 0 using Harris' and Doeblin's theorems. The proof of Theorem 2 (Stationary solutions) is presented in Section 4 which is divided in three parts: in the first part, we present a proof which uses the exponential stability of the non-interacting case; in the second part, we present an alternative proof for the existence of stationary solutions which does not involve the Doeblin-Harris method; and in the last part, we present proof for the formula of [START_REF] Romani | Mean-field analysis of selective persistent activity in presence of short-term synaptic depression[END_REF] in the case of short-term synaptic plasticity (3). Finally, Section 5 is dedicated to the proof of Theorem 3 (Exponential stability in the weak connectivity regime).

Well-posedness

This section is dedicated to the proof of Theorem 1, which we decompose is several Lemmas. First, we verify the a priori L 1 -stability of the solutions to [START_REF] Beiran | Contrasting the effects of adaptation and synaptic filtering on the timescales of dynamics in recurrent networks[END_REF], an important technical result we use later in the proof. Then, we introduce a linearized version of ( 1) and show that it is well-posed by an application of Banach's fixed-point theorem. Another Banach's fixed-point argument is used to treat the nonlinearity of (1) and concludes the proof of the well-posedness of (1) in L 1 . Finally, we prove the global bound in L 1 + (w) of the solutions (point (ii) of Theorem 1), which we will need in order to apply Harris' theorem in the next Sections.

Lemma 1 (A priori L 1 -stability). Grant Assumption 1. If (ρ, x) is a weak solution to (1) for the initial datum u 0 ∈ L 1 + , then

ρ t L 1 = u 0 L 1 , ∀t > 0.
Proof. By a standard cut-off in time argument, we have that for all T > 0 and for all ϕ

∈ C ∞ c (R + × R + × R * + ), ∞ 0 ∞ 0 ρ T (a, m)ϕ(T, a, m)dadm - ∞ 0 ∞ 0 u 0 (a, m)ϕ(0, a, m)dadm = T 0 ∞ 0 ∞ 0 ρ t (a, m) [∂ t + ∂ a -λm∂ m ]ϕ + (ϕ(t, 0, γ(m)) -ϕ(t, a, m))f (a, m, εx t ) dadmdt.
Let χ be a function in

C ∞ c (R + × R * + , R + ) such that χ(a, m) = 1, for all a 2 + m 2 ≤ 1.
For all n ∈ N * , we write

ϕ n ∈ C ∞ (R + × R + × R * + )
the classical solution to the transport equation

∂ t ϕ n (t, a, m) + ∂ a ϕ n (t, a, m) -λm∂ m ϕ n (t, a, m) = 0, (8a)
ϕ n (0, a, m) = χ(a/n, m/n). (8b)
Because of the finite speed of propagation of the transport equation, for all n, there exists a function

ϕ n ∈ C ∞ c (R + × R + × R * + ) such that ϕ n (t, a, m) = ϕ n (t, a, m), for all (t, a, m) ∈ [0, T ] × R + × R * . Hence, for all n ∈ N * , ∞ 0 ∞ 0 ϕ n (T, a, m)ρ T (a, m)dadm - ∞ 0 ∞ 0 ϕ n (0, a, m)u 0 (a, m)dadm = T 0 ∞ 0 ∞ 0 ∂ t ϕ n +∂ a ϕ n -λm∂ m ϕ n + ϕ n (t, 0, γ(m)) -ϕ n (t, a, m) f (a, m, εx t ) ρ t (a, m)dadmdt.
As ϕ n is a solution to Eq. (8a) on time [0, T ], we get

∞ 0 ∞ 0 ϕ n (T, a, m)ρ T (a, m)dadm - ∞ 0 ∞ 0 ϕ n (0, a, m)u 0 (a, m)dadm = T 0 ∞ 0 ∞ 0 ϕ n (t, 0, γ(m)) -ϕ n (t, a, m) f (a, m, εx t ) ρ t (a, m)dadmdt. For all (t, a, m) ∈ [0, T ] × R + × R * + , ϕ n (t, a, m) ---→ n→∞ 1
, since the initial datum tends to 1 as n → ∞ (Eq. ( 8b)) and by finite speed of propagation. Thus, by dominated convergence, we get

∞ 0 ∞ 0 ρ T (a, m)dadm - ∞ 0 ∞ 0 u 0 (a, m)dadm = 0. (9)
Since ρ is nonnegative, this concludes the proof.

Lemma 1 will allow us to prove the well-posedness of (1) by the means of fixed-point arguments. Let us first introduce a linearized version of Eq. ( 1): for all x ∈ C(R + ), we consider the linear evolution problem

∂ t ρ t + ∇ • (bρ t ) = -f (a, m, εx t )ρ t , (10a)
ρ t (0, m) = ½ m>γ(0) (γ -1 ) ′ (m) ∞ 0 f (a, γ -1 (m), εx t )ρ t (a, γ -1 (m))da, (10b)
ρ 0 = u 0 . (10c)
We can see Eq. ( 10) as the Kolmogorov forward equation of a time-dependent Makrov process. Indeed, we can rewrite Eqs. (10a) and (10b) as

∂ρ t = L t ρ t (11)
where, for all suitable test function φ : R

+ × R * + → R, L * t φ(a, m) = b(a, m)∇ • φ(a, m) + [φ(0, γ(m)) -φ(a, m)]f (a, m, εx t ). (12)
L * t is the time-dependent generator of a piecewise deterministic Markov process with degenerate jumps.

The linearized equation [START_REF] Duarte | Stability, convergence to equilibrium and simulation of non-linear Hawkes processes with memory kernels given by the sum of Erlang kernels[END_REF] will play a special role in the following Sections and it therefore deserves its own Proposition: Proposition 1 (Well-posedness of the linearized equation ( 10)). Grant Assumption 1. For any initial datum u 0 ∈ L 1 + and any x ∈ C(R + ), there exists a unique weak solution ρ x ∈ C(R + , L 1 + ) to Eq. [START_REF] Duarte | Stability, convergence to equilibrium and simulation of non-linear Hawkes processes with memory kernels given by the sum of Erlang kernels[END_REF]. Furthermore, ρ x satisfies (i) For all t > 0 and for all m ∈ R * + ,

ρ x t (0, m) = ½ m>γ(0) (γ -1 ) ′ (m) ∞ 0 f (a, γ -1 (m), εx)ρ x t (a, γ -1 (m))da, ρ x t (a, m) =      u 0 (a -t, e λt m) exp λt - t 0 f (a -t + s, e λ(t-s) m, εx)ds if a ≥ t, ρ x t-a (0, e λa m) exp λa - t t-a f (a -t + s, e λ(t-s) m, εx)ds if 0 < a < t.
(ii) For all t > 0 and for all φ

∈ C ∞ c (R + × R * + ), ρ x t , φ = u 0 , φ + t 0 ρ x t , L * x φ ds. ( 14
)
Proof. Fix x ∈ C(R + ). For all p ∈ C(R + , L 1 + (R * + )
) and u 0 ∈ L 1 + , we know, from the standard theory of transport equations, that there is a unique weak solution to

∂ t ρ t + ∇ • (bρ t ) = -f (a, m, εx t )ρ t , ρ t (0, m) = p t (m), ρ 0 = u 0 ,
which we denote ρ x,p and is given by the representation formula,

ρ x,p t (a, m) :=      u 0 (a -t, e λt m) exp λt - t 0 f (a -t + s, e λ(t-s) m, εx s )ds if a ≥ t, p t-a (e λa m) exp λa - t t-a f (a -t + s, e λ(t-s) m, εx s )ds if 0 < a < t. ρ x,p is in C(R + , L 1 ), since ∀t ∈ R + , ρ x,p t L 1 ≤ u 0 L 1 + t 0 p s L 1 ds.
We have

½ m>γ(0) (γ -1 ) ′ (m) ∞ 0 f (a, γ -1 (m), x)ρ x,p t (a, γ -1 (m))da (t,m)∈R + ×R * + ∈ C(R + , L 1 + (R * + )) since ∀t ∈ R + , ∞ γ(0) (γ -1 ) ′ (m) ∞ 0 f (a, γ -1 (m), x)ρ x,p t (a, γ -1 (m))dadm ≤ f ∞ ρ x,p t L 1 ≤ f ∞ u 0 L 1 + t 0 p s L 1 ds .
Hence, we can define, for any T > 0, the operator Φ x T :

C([0, T ], L 1 + (R * + )) → C([0, T ], L 1 + (R * + )) p → ½ m>γ(0) (γ -1 ) ′ (m) ∞ 0 f (a, γ -1 (m), x)ρ x,p t (a, γ -1 (m))da (t,m)∈[0,T ]×R * + . For any p, q ∈ C([0, T ], L 1 + (R * + )), Φ x T (p) -Φ x T (q) C([0,T ],L 1 ) ≤ f ∞ sup t∈[0,T ] ρ x,p t -ρ x,q t L 1 ≤ f ∞ T 0 p s -q s L 1 ds ≤ T f ∞ p -q C([0,T ],L 1 ) . Therefore, if 0 < T < f -1 ∞ , Φ x
T is a contraction. By Banach's fixed-point theorem, there exists a unique ρ x ∈ C([0, T ], L 1 + ) solving Eq. [START_REF] Duarte | Stability, convergence to equilibrium and simulation of non-linear Hawkes processes with memory kernels given by the sum of Erlang kernels[END_REF]. Since the choice of the contracting T does not depend on the initial datum, we can iterate the above argument on successive time intervals of length T and conclude that there exists a unique ρ x ∈ C(R + , L 1 + ) solving Eq. ( 10) for which the formula (i) is satisfied. Then, (ii) follows from a standard cut-off in time argument. Now, we can prove the existence and uniqueness of a solution to the nonlinear problem (1) by the means of a second application of Banach's fixed-point theorem.

Proof of the well-posedness of (1) in L 1 . For any x ∈ C(R + ), we take the ρ x given by Proposition 1. We have

t 0 R + ×R * + h(t -s)f (εx s )ρ x s dadmds t∈R + ∈ C(R + ) since ∀t ∈ R + , t 0 R + ×R * + h(t -s)f (εx s )ρ x s dadmds ≤ h ∞ f ∞ t 0 ρ x s L 1 ds.
Hence, for any T > 0, we can define the operator

Ψ T : C([0, T ]) → C([0, T ]) x → t 0 R + ×R * + h(t -s)f (εx s )ρ x s dadmds t∈[0,T ]
.

For any x, y ∈ C([0, T ]), we have

Ψ T (x) -Ψ T (y) C([0,T ]) ≤ T h ∞ sup t∈[0,T ] R + ×R * + |f (εx t )ρ x t -f (εy t )ρ y t | dadm ≤ T h ∞ sup t∈[0,T ] εL f |x t -y t | ρ x t L 1 + f ∞ ρ x t -ρ y t L 1 . By Grönwall's lemma, ρ x t L 1 ≤ u 0 L 1 exp( f ∞ t), since ∀t ∈ [0, T ], ρ x t L 1 ≤ u 0 L 1 + f ∞ t 0 ρ x s L 1 ds.
On the other hand, we have, for all t ∈ [0, T ],

ρ x t -ρ y t L 1 ≤ f ∞ t 0 ∞ 0 ∞ 0 ρ x s (a, m) exp - t s f (u -s, e -λ(u-s) γ(m), εx u )du ρ y s (a, m) exp - t s f (u -s, e -λ(u-s) γ(m), εy u )du dadmds ≤ f ∞ t 0 ρ x s -ρ y s L 1 ds + tε f ∞ L f x -y C([0,T ]) t 0 ρ x s L 1 ds.
Hence, by Grönwall's lemma, for all t ∈ [0, T ],

ρ x t -ρ y t L 1 ≤ εL f u 0 L 1 exp( f ∞ t) -1 2 f ∞ x -y C([0,T ]) .
Gathering the bounds, we get

Ψ T (x) -Ψ T (y) C([0,T ]) ≤ T ε h ∞ L f u 0 L 1 exp( f ∞ T ) 1 + exp( f ∞ T ) x -y C([0,T ]) .
For T small enough, Ψ T is a contraction and, by Banach fixed-point theorem, has a unique fixedpoint. Thus, there exists a unique solution (ρ, x) ∈ C([0, T ], L 1 + ). Since by Lemma 1, ρ T L 1 = u 0 L 1 , we can iterate this argument on successive time intervals of length T and conclude that there exists a unique solution in C(R + , L 1 + ).

To conclude the proof of Theorem 1, it remains to show the estimate Eq. ( 5). Under Assumption 2, the weight function

w : R + × R + → [1, ∞),
(a, m) → 1 + m defined above, satisfies w(a, m) → ∞ when m → ∞ and the Lyapunov condition on m:

∃α > 0, b ≥ 0 such that L * t w ≤ -αw + b. (16)
Indeed, for all (t, a, m

) ∈ R + × R + × R * + , L * t w(a, m) = -λm + Γ(m)f (a, m, εx t ) ≤ -λw(a, m) + λ + Γ ∞ f ∞ .
Importantly, the constants α and b do not depend on x.

Lemma 2 (Global bound in L 1 + (w)). Grant Assumptions 1 and 2. If the initial datum u 0 is in

L 1 + (w), then ρ t ∈ L 1 + (w) for all t ≥ 0. Moreover, ∀t > 0, ρ t L 1 (w) ≤ u 0 L 1 (w) e -αt + b α (1 -e -αt ), (17)
where the constants α and b are taken from the Lyapunov condition [START_REF] Gerstner | Neuronal dynamics: From single neurons to networks and models of cognition[END_REF].

Proof. We divide the proof in two steps: first, we prove that the solution is stable in L 1 + (w) with a weaker and time dependent bound; then, we use this first bound to apply the dominated convergence theorem and obtain Eq. ( 17) by Grönwall's lemma.

Step 1. Fix any

T > 0. Let χ ∈ C ∞ c (R + , R + ) be a non-increasing function such that χ(x) = 1 if 0 ≤ x ≤ 1 and χ(x) = 0 if x > 2. For all n ∈ N * , let us write ϕ k (a)χ n (m) := χ(a/k)χ(m/n). For all n, k, wχ n ϕ k ∈ C ∞ c (R + × R + , R + ). Hence, ∀n ∈ N * , ρ T , wχ n ϕ k = u 0 , wχ n ϕ k + T 0 ρ t , L * x (wχ n ϕ k ) dt,
where

L * x (wχ n ϕ k ) = ∂ a (wχ n ϕ k ) -λm∂ m (wχ n ϕ k ) + w(γ(m))χ n (γ(m))ϕ k (0) -wχ n ϕ k f = wχ n 1 k χ ′ (a/k) + wϕ k 1 n χ ′ (m/n) -λmwχ n + (w(γ(m))χ n (γ(m))ϕ k -wχ n ϕ k )f.
From the L 1 -stability and the fact that both w∂ m χ n and wχ n are bounded, we take the limit k → ∞ with the dominated convergence theorem to obtain

ρ T , wχ n = u 0 , wχ n + T 0 ρ t , w 1 n χ ′ (m/n) -λmwχ n + (w(γ(m))χ n (γ(m)) -wχ n )f dt. ( 18
)
From the properties of χ and γ, we get

w(0, γ(m))χ n (γ(m)) ≤ w(0, m + Γ ∞ )χ n (γ(m)) ≤ (1 + Γ ∞ )w(a, m)χ n (m), 1 n wχ ′ (m/n) ≤ 1 + 2n n χ ′ ∞ .
Then, since f is bounded, there exists a constant C, which does not depend on n, such that ρ T , wχ n ≤ u 0 , wχ n + T 0 ρ t , C(wχ n (m) + 1) dt.

We can now apply Grönwall's lemma to obtain ρ T , wχ n ≤ max( u 0 , wχ n Ce Ct , C).

It follows from Fatou's lemma that ρ t = S x T u 0 ∈ L 1 + (w).

Step 2. To improve the previous estimate, we come back to [START_REF] Gigante | Diverse population-bursting modes of adapting spiking neurons[END_REF] and use dominated convergence in n (domination being guaranteed by the bound above) to get

ρ T , w = u 0 , w + T 0 ρ t , L * x w dt.
By the Lyapunov condition ( 16),

ρ T L 1 (w) ≤ u 0 L 1 (w) -α T 0 ρ t L 1 (w) + T b.
Finally, Eq. ( 17) is obtained by Grönwall's lemma.

Remark . Following the same steps as in the proof above, we can show that the bound Eq. ( 17) also holds for the linearized equation [START_REF] Duarte | Stability, convergence to equilibrium and simulation of non-linear Hawkes processes with memory kernels given by the sum of Erlang kernels[END_REF] and does not depend on x nor the constants α and b.

Exponential stability in the non-interacting case

If x ∈ C(R + ) in the linearized equation ( 10) is time-invariant, i.e. x ≡ x for some x ∈ R, then Eq. ( 10) can be seen as the dynamics of a non-interacting population of neurons. In this section, we prove the exponential stability in the non-interacting case using Harris' and Doeblin's theorems. This is the key result of this work and will allow us to prove the existence and uniqueness of the stationary solution to (1) (Section 4) and the exponential convergence to it (Section 5). For x ∈ R, u 0 ∈ L 1 , we denote ρ x t the unique solution to Eq. ( 10) for the initial datum u 0 and x ≡ x, given by Proposition 1. We write, using the semigroup notation,

S x t u 0 := ρ x t , ∀t ≥ 0. ( 19
)
To show that the Eq. ( 19) is exponentially stable we will use Harris' theorem in the general case or Doeblin's theorem if Assumption 4 is granted. In both cases, the main technical difficulty is to verify the Doeblin minoration condition (Lemma 3) as the jumps of the process described by Eq. ( 12) are degenerate and the model is two-dimensional.

Lemma 3 (Doeblin minoration condition). Grant Assumptions 1 and 3. Fix any x ∈ R. For all R > 0, there exists T > 0 and a positive non-zero measure ν such that

∀u 0 ∈ L 1 + , S x T u 0 ≥ ν R + ×]0,R] u 0 dadm. (20)
Proof. We proceed in two steps. First (Step 1), we choose a time T > 0 and a rectangle [0, ā] × [m, m] ⊂ R + ×R * + (with nonzero Lebesgue measure) and show that the density S x T u 0 ∈ L 1 has a lower bound on [0, ā] × [m, m] which depends on a Lebesgue integral in R 2 + involving u 0 . Then (Step 2), we perform a change of variable to express this lower bound in terms of R + ×]0,R] u 0 dadm. The proof only relies on the expression of S x t u 0 given by the method of characteristics (see Proposition 1) and this allows treating a typically probabilistic question -the Doeblin minoration condition -from a transport point of view. This is possible because S x t is the stochastic (mass-conservative) semigroup of a piecewise deterministic Markov process.

The constants ∆ abs , σ and C γ are taken from Assumption 3.

Step 1: Fix R > 0. Since γ(e -λ∆ abs γ(0)) > γ(0) and γ(e -λt γ(e -λ∆ abs R)) → γ(0) as t → ∞, there exists ā > 0 and T > ā + ∆ abs such that m =: γ(e -λ(T -ā-∆ abs ) γ(e -λ∆ abs R)) < e -λā γ(e -λ∆ abs γ(0)) =: m.

Eq. ( 21) has the following heuristic interpretation: if we see S x t as the stochastic semigroup of the piecewise deterministic Markov process defined by the generator Eq. ( 12), for any initial point (a 0 , m 0 ) ∈ R + ×]0, R] and any landing point (a, m) ∈ [0, ā] × [m, m] at time T , there is a 'possible' trajectory going from (a 0 , m 0 ) to (a, m), with exactly two jumps (spikes). Since the trajectories of the process are determined by the jump times, we will exploit the fact that these 'possible' trajectories correspond to jump times with strictly positive probability density. Below, we take a transport point of view on this probabilistic argument.

For all (a, m)

∈ [0, ā] × [m, m], (S x T u 0 )(a, m) ≥ ½ {a<T } (S x T -a u 0 )(0, e λa m) exp λa - T T -a f (a -T + s, e λ(T -s) m, x)ds ≥ ½ {a<T } e -f ∞ T e λa (S x T -a u 0 )(0, e λa m) ≥ ½ {a<T } e -f ∞ T σe λa (γ -1 ) ′ (e λa m) ∞ ∆ abs (S x T -a u 0 )(a ′ , γ -1 (e λa m))da ′ = ½ {a<T } e -f ∞ T σ d dm γ -1 (e λa m) ∞ ∆ abs (S x T -a u 0 )(a ′ , γ -1 (e λa m)) (⋆) da ′ .
Above, we went back in time to the last jump time T -a. Let us notice that γ -1 (e λa m) ≥ γ -1 (e λa m) > 0. We can therefore define

a * a,m := 1 λ log γ(0) -log γ -1 (e λa m) .
Note that a * a,m satisfies γ -1 (e λa * a,m γ -1 (e λa m)) = 0. In other words, a * a,m is the minimal time between the last and second last jumps for a trajectory landing at (a, m) at time T . We can easily verify that, by our choice of {T, ā, m, m}, ∆ abs ≤ a * a,m < T -a -∆ abs . This guarantees that it is possible to make two jumps in [0, T ] and land at (a, m) at time T while respecting the absolute refractoriness of the neuron (i.e. there needs to be a time interval ≥ ∆ abs between jumps). This allows us to go further back in time to the second last jump.

For all

a ′ ∈ [a * a,m , T -a -∆ abs ], (⋆) ≥ ½ {a ′ <T -a} e -f ∞ T σ (γ -1 ) ′ (e λa ′ γ -1 (e λa m)) e λa ′ ∞ ∆ abs (S x T -a-a ′ u 0 )(a ′′ , γ -1 (e λa ′ γ -1 (e λa m))) (⋆⋆) da ′′ .
Then, we can go further back to time 0 to get u 0 :

(⋆⋆) ≥ ½ {a ′′ ≥T -a-a ′ } e -f ∞ T e λ(T -a-a ′ ) u 0 (a ′′ -(T -a -a ′ ), e λ(T -a-a ′ ) γ -1 (e λa ′ γ -1 (e λa m))).
Putting all the lower bounds together, we get

(S x T u 0 )(a, m) ≥ ½ {a<T } e -3 f ∞ T σ 2 T -a-∆ abs a * a,m ∞ T -a-a ′ d dm e λ(T -a-a ′ ) γ -1 (e λa ′ γ -1 (e λa m)) u 0 (a ′′ -(T -a -a ′ ), e λ(T -a-a ′ ) γ -1 (e λa ′ γ -1 (e λa m)))da ′′ da ′ . Since γ ′ ≤ 1 (Assumption 3), d dm e λ(T -a-a ′ ) γ -1 (e λa ′ γ -1 (e λa m)) ≥ e λT .
Thus,

(S x T u 0 )(a, m) ≥ ½ {a<T } e (λ-3 f ∞ )T σ 2 T -a-∆ abs a * a,m ∞ 0 u 0 (a 0 , e λ(T -a-a ′ ) γ -1 (e λa ′ γ -1 (e λa m)))da 0 da ′ . ( 22
)
We have obtained that on [0, ā] × [m, m], the density (S x T u 0 ) is lower bounded by a constant depending on a Lebesgue integral on R 2 + involving u 0 .

Step 2: Now, we want express the lower bound Eq. ( 22) in terms of R + ×]0,R] u 0 dadm by a change of variable. Let us define the function ψ T a,m :

ψ T a,m : [a * a,m , T -a -∆ abs ] → R + , a ′ → e λ(T -a-a ′ ) γ -1 (e λa ′ γ -1 (e λa m)).
We verify that (ψ T a,m ) ′ > 0:

For all a ′ ∈ [a * a,m , T -a], (ψ T a,m ) ′ (a ′ ) = λe λ(T -a-a ′ ) (γ -1 ) ′ (e λa ′ γ -1 (e λa m))e λa ′ γ -1 (e λa m) -γ -1 (e λa ′ γ -1 (e λa m)) . ( 23
)
As Γ > 0 and γ ′ ≤ 1 (Assumption 3), we have

(ψ T a,m ) ′ (a ′ ) > λe λ(T -a-a ′ ) (γ -1 ) ′ (e λa ′ γ -1 (e λa m))e λa ′ γ -1 (e λa m) -e λa ′ γ -1 (e λa m) = λe λ(T -a) γ -1 (e λa m) (γ -1 ) ′ (e λa ′ γ -1 (e λa m)) ≥1 -1 ≥ 0.
Therefore, ψ T a,m is a strictly increasing

C 1 -diffeomorphism from [a * a,m , T -a-∆ abs ] to [ψ T a,m (a * a,m ), ψ T a,m (T - a -∆ abs )].
We can now rewrite Eq. ( 22):

(S x T u 0 )(a, m) ≥ e (λ-3 f ∞ )T σ 2 T -a-∆ abs a * a,m ∞ 0 u 0 (a 0 , ψ T a,m (a ′ ))da 0 da ′ = e (λ-3 f ∞ )T σ 2 ψ T a,m (T -a-∆ abs ) ψ T a,m (a * a,m ) ∞ 0 u 0 (a 0 , m 0 ) ((ψ T a,m ) -1 ) ′ (m 0 ) da 0 dm 0 .
Going back to Eq. ( 23), and using the fact that there exists

C γ such that C γ ≤ γ ′ ≤ 1 (Assumption 3), we have, for all a ′ ∈ [a * a,m , T -a -∆ abs ], (ψ T a,m ) ′ (a ′ ) ≤ λe λ(T -a-a ′ ) C -1 γ e λa ′ γ -1 (e λa m) ≤ λe λT C -1 γ m. Hence, (S x T u 0 )(a, m) ≥ e -3 f ∞ T σ 2 C γ λm ψ T a,m (T -a-∆ abs ) ψ T a,m (a * a,m ) ∞ 0 u 0 (a 0 , m 0 )da 0 dm 0 .
In addition, by our choice of {T, ā, m, m}, we have

ψ T a,m (a * a,m ) = 0, ψ T a,m (T -a -∆ abs ) = e λ∆ abs γ -1 (e λ(T -a-∆ abs ) γ -1 (e λa m)) > R. Therefore, (S x T u 0 )(a, m) ≥ e -3 f ∞ T σ 2 C γ λm R 0 ∞ 0 u 0 (a 0 , m 0 )da 0 dm 0 .
Since we have supposed that (a, m) ∈ [0, ā] × [m, m], this concludes the proof.

With the Lyapunov condition ( 16) and the Doeblin minoration condition [START_REF] Höpfner | Ergodicity and limit theorems for degenerate diffusions with time periodic drift. application to a stochastic Hodgkin-Huxley model[END_REF], we can apply a version of Harris' theorem: Theorem 4 (Harris). Grant Assumptions 1 and 3. For all x ∈ R, there exists a unique ρ

x ∞ ∈ L 1 + (w) with ρ x ∞ L 1 = 1 such that S x t ρ x ∞ = ρ x ∞
, for all t ≥ 0, and there exists K ≥ 1 and a > 0 such that for all initial data u 0 ∈ L 1 + (w) with u 0 L 1 = 1,

S x t u 0 -ρ x ∞ L 1 (w) ≤ Ke -at u 0 -ρ x ∞ L 1 (w)
, ∀t ≥ 0.

Furthermore, by Lemma 2, we have that

ρ x ∞ L 1 (w)
≤ b α , where the constants α and b are taken from the Lyapunov condition [START_REF] Gerstner | Neuronal dynamics: From single neurons to networks and models of cognition[END_REF].

Proof. This is a classic result which proof can be found in the work of Hairer and Mattingly [START_REF] Hairer | Yet another look at Harris' ergodic theorem for Markov chains[END_REF].

If, in addition, Assumption 4 holds, we can simply apply Doeblin's theorem: Theorem 5 (Doeblin). Grant Assumptions 1, 3 and 4. For all x ∈ R, there exists a unique

ρ x ∞ ∈ L 1 + with ρ x ∞ L 1 = 1 such that S x t ρ x ∞ = ρ x ∞
, for all t ≥ 0, and there exists K ≥ 1 and a > 0 such that for all initial data

u 0 ∈ L 1 + with u 0 L 1 = 1, S x t u 0 -ρ x ∞ L 1 ≤ Ke -at u 0 -ρ x ∞ L 1 , ∀t ≥ 0. (25)
Proof. See, for example, Theorem 2.3 in [START_REF] José | Asymptotic behaviour of neuron population models structured by elapsed-time[END_REF].

We say that ρ x ∞ is the invariant probability measure of the semigroup (S x t ) t∈R + . Note that both theorems imply the next corollary.

Corollary 1. Grant Assumptions 1 and 3. For all x ∈ R, there exists a unique ρ

x ∞ ∈ L 1 + (w) with ρ x ∞ L 1 = 1 solving ∂ a ρ x ∞ (a, m) -λ∂ m (mρ x ∞ (a, m)) = -f (a, m, x)ρ x ∞ (a, m), (26a)
ρ x ∞ (0, m) = ½ m>γ(0) (γ -1 ) ′ (m) ∞ 0 f (a, γ -1 (m), x)ρ x ∞ (a, γ -1 (m))da, (26b)
in the weak sense. Furthermore, we have that

ρ x ∞ ∈ C(R + , L 1 + (R * + )) ∩ L ∞ (R + , L 1 + (R * + )).

Stationary solutions for arbitrary connectivity strength

In this section, we study the stationary solutions to (1), namely the solution to 27) in the weak sense.

∂ a ρ ∞ (a, m) -λ∂ m (mρ ∞ (a, m)) = -f (a, m, εx ∞ )ρ ∞ (a, m), (27a)
ρ ∞ (0, m) = ½ m>γ(0) (γ -1 ) ′ (m) ∞ 0 f (a, γ -1 (m), εx ∞ )ρ ∞ (a, γ -1 (m))da, (27b)
x ∞ = ∞ 0 ∞ 0 h(a, m)f (a, m, εx ∞ )ρ ∞ (a, m)dadm. (27c) Definition. (ρ ∞ , x ∞ ) ∈ L 1 + (w) ∩ C(R + , L 1 + (R * + )) ∩ L ∞ (R + , L 1 + (R * + )) × R + is a stationary solution to (1) if ρ ∞ L 1 = 1 and if it solves Eq. (

Existence and uniqueness using the Doeblin-Harris method.

We present two important Lipschitz continuity results, which will allow us to prove the existence (and the uniqueness when ε is small) of stationary solutions. The following lemma plays the same role as Theorem 4.5 in [START_REF] José | Asymptotic behaviour of neuron population models structured by elapsed-time[END_REF]: Lemma 4 (Lipschitz continuity at finite T ). Grant Assumptions 1 and 3. For all initial data u 0 ∈ L 1 + (w) and for all T > 0, there exists a constant C T, u 0 L 1 (w) > 0 such that

∀ x 1 , x 2 ∈ R, S x 1 T u 0 -S x 2 T u 0 L 1 (w) ≤ C T, u 0 L 1 (w) | x 1 -x 2 |. (28)
Proof. For all t > 0,

S x 1 t u 0 -S x 2 t u 0 L 1 (w) = ∞ 0 ∞ t u 0 (a -t, e λt m) exp λt - t 0 f (a -t + s, e λ(t-s) m, x 1 )ds -u 0 (a -t, e λt m) exp λt - t 0 f (a -t + s, e λ(t-s) m, x 2 )ds w(a, m)dadm + ∞ 0 t 0 ρ x 1 t-a (0, e λa m) exp λa - t t-a f (a -t + s, e λ(t-s) m, x 1 )ds -ρ x 2 t-a (0, e λa m) exp λa - t t-a f (a -t + s, e λ(t-s) m, x 2 )ds w(a, m)dadm =: Q 1 + Q 2 . Q 1 = ∞ 0 ∞ 0 u 0 (a, m) exp - t 0 f (a + s, e -λs m, x 1 )ds -exp - t 0 f (a + s, e -λs m, x 2 )ds w(a + t, e -λt m)dadm ≤ ∞ 0 ∞ 0 u 0 (a, m) t 0
f (a + s, e -λs m, x 1 ) -f (a + s, e -λs m, x 2 ) ds w(a + t, e -λt m)dadm

≤ tL f | x 1 -x 2 | ∞ 0 ∞ 0 u 0 (a, m)w(a + t, e -λt m)dadm ≤ tL f u 0 L 1 (w) | x 1 -x 2 |,
where in the last inequality we used

w(a + t, e -λt m) ≤ w(a, m), ∀a ≥ 0, m ≥ 0. (29)
Q 2 = ∞ 0 t 0 ρ x 1 t-a (0, m) exp - t t-a f (a -t + s, e λ(t-s-a) m, x 1 )ds -ρ x 2 t-a (0, m) exp - t t-a
f (a -t + s, e λ(t-s-a) m, x 2 )ds w(a, e -λa m)dadm.

By changes of variables,

Q 2 = ∞ 0 t 0 ρ x 1 s (0, m) exp - t-s 0 f (u, e -λu m, x 1 )du -ρ x 2 s (0, m) exp - t-s 0 f (u, e -λu m, x 2 )du w(t -s, e -λ(t-s) m)dsdm ≤ ∞ 0 t 0 ρ x 1 s (0, m) exp - t-s 0 f (u, e -λu m, x 1 )du -exp - t-s 0 f (u, e -λu m, x 2 )du w(t -s, e -λ(t-s) m)dsdm + ∞ 0 t 0 ρ x 1 s (0, m) -ρ x 2 s (0, m) w(t -s, e -λ(t-s) m)dsdm =: Q 2,1 + Q 2,2 Q 2,1 ≤ t f ∞ L f | x 1 -x 2 | t 0 ∞ 0 ∞ 0 (γ -1 ) ′ (m) ρ x 1 s (a, γ -1 (m))w(t, m)dadmds ≤ t f ∞ L f | x 1 -x 2 | t 0 ∞ 0 ∞ 0 ρ x 1 s (a, m)w(t, m + Γ ∞ )dadmds ≤ t(1 + Γ ∞ ) f ∞ L f | x 1 -x 2 | t 0 ρ x 1 s L 1 (w) ds,
where in the last inequality we used

w(t, m + Γ ∞ ) = 1 + m + Γ ∞ ≤ (1 + Γ ∞ )w(a, m), ∀a ≥ 0, m ≥ 0. (30)
By Lemma 2,

Q 2,1 ≤ t 2 (1 + Γ ∞ ) f ∞ L f u 0 L 1 (w) + b α | x 1 -x 2 |. Q 2,2 ≤ f ∞ t 0 ∞ 0 ∞ 0 (γ -1 ) ′ (m) ρ x 1 s (a, γ -1 (m)) -ρ x 2 s (a, γ -1 (m)) w(t, m)dadmds ≤ f ∞ t 0 ∞ 0 ∞ 0 ρ x 1 s (a, m) -ρ x 2 s (a, m) w(t, m + Γ ∞ )dadmds ≤ (1 + Γ ∞ ) f ∞ t 0 S x 1 s u 0 -S x 2 s u 0 L 1 (w) ds,
where again, in the last inequality, we used Eq. ( 30). Fix T > 0. Gathering the bounds for

Q 1 , Q 2,1
and Q 2,2 we see that there exists constants C > 0 and

C ′ T, u 0 L 1 (w)
> 0 such that, for all t ∈ [0, T ],

S x 1 t u 0 -S x 2 t u 0 L 1 (w) ≤ C t 0 S x 1 s u 0 -S x 2 s u 0 L 1 (w) ds + tC ′ T, u 0 L 1 (w) | x 1 -x 2 |.
By Grönwall's lemma, for all t ∈ [0, T ],

S x 1 t u 0 -S x 2 t u 0 L 1 (w) ≤ C ′ T, u 0 L 1 (w) | x 1 -x 2 | C exp(Ct) -1 . (31)
Since Eq. ( 31) holds for all t ∈ [0, T ], this achieves the proof.

For any (ũ, x) ∈ L 1 + (]γ(0), +∞[) × R, consider the transport equation

∂ a ̺(a, m) -λ∂ m (m̺(a, m)) = -f (a, m, x)̺(a, m), ̺(0, m) = ũ(m).
It has a unique weak solution

ρ ũ,x ∞ ∈ C(R + , L 1 + (R * + )) ∩ L ∞ (R + , L 1 + (R * + ))
given by the method of characteristics: for all (a, m

) ∈ R + × R * + , ρ ũ,x ∞ (a, m) = ũ(e λa m) exp λa - a 0 f (s, e λ(a-s) m, x)ds . (32)
We can now define the operator

Φ := (Φ 1 , Φ 2 ) on L 1 + (]γ(0), +∞[) × R where, for all (ũ, x) ∈ L 1 + (]γ(0), +∞[) × R, Φ 1 (ũ, x)(m) := ½ m>γ(0) (γ -1 ) ′ (m) ∞ 0 f (a, γ -1 (m), x)ρ (ũ,x) ∞ (a, γ -1 (m))da, (33a)
Φ 2 (ũ, x) := ∞ 0 ∞ 0 h(a, m)f (a, m, x)ρ ũ,x ∞ (a, m)dadm. ( 33b
) (ρ ∞ , x ∞
) is a stationary solution if and only if it is a fixed-point of Φ. Whence, we get the a priori estimates:

Lemma 6. Grant Assumptions 1 and 3. There exists θ ∈]0, 1[ such that for all (ũ, x)

∈ L 1 + (]γ(0), +∞[)× R, (i) Φ 1 (ũ, x) L 1 = ũ L 1 . (ii) For all m ∈ R * + , |Φ 1 (ũ, x)(m)|≤ ½ m>γ(0) f ∞ λγ -1 (m) ũ L 1 . (iii) ∞ 0 Φ 1 (ũ, x)(m)mdm ≤ max ∞ 0 ũ(m)mdm, γ(0) 1 -θ ũ L 1 . (iv) For all β ∈]0, min(f ) λ [, ∞ γ(0) Φ 1 (ũ, x)(m) γ -1 (m) β dm ≤ f ∞ λγ(0) β min(f ) λ -β ũ L 1 . (v) Φ 2 (ũ, x) ≤ h ∞ ũ L 1 . Proof. (i) By changes of variables on m, Φ 1 (ũ, x) L 1 = ∞ 0 ∞ 0 f (a, m, x)ũ(e λa m) exp λa - a 0 f (s, e λ(a-s) m, x)ds dadm = ∞ 0 ũ(m) ∞ 0 f (a, e -λa m, x) exp - a 0 f (s, e -λs m, x)ds da =1 (by Assumption 3 (i)) dm. (ii) |Φ 1 (ũ, x)(m)| ≤ ½ m>γ(0) f ∞ ∞ 0 ũ(e λa γ -1 (m)) exp(λa)da = ½ m>γ(0) f ∞ λγ -1 (m) ∞ 0 ũ(e λa γ -1 (m))γ -1 (m)λ exp(λa)da = ½ m>γ(0) f ∞ λγ -1 (m) ∞ γ -1 (m) ũ(y)dy ≤ ũ L 1
, where for the last equality we used the change of variable y = e λa γ -1 (m).

(iii) Performing the same change of variable as for (i) and using the fact that γ(m) ≤ γ(0) + m, ∀m ∈ R + (since γ ′ ≤ 1), we have

≤ 1 -(1 -e -λǫ ) exp(-f ∞ ǫ) =: θ < 1. Whence, ∞ 0 Φ 1 (ũ, x)(m)mdm ≤ θ ∞ 0 ũ(m)mdm + γ(0) ũ L 1 .
To see that

∞ 0 Φ 1 (ũ, x)(m)mdm ≤ max ∞ 0 ũ(m)mdm, γ(0) 1 -θ ũ L 1 ,
we can distinguish three cases: if

∞ 0 ũ(m)mdm = ∞, the inequality is trivial; if γ(0) 1-θ ũ L 1 ≤ ∞ 0 ũ(m)mdm < +∞, then ∞ 0 Φ 1 (ũ, x)(m)mdm ≤ ∞ 0 ũ(m)mdm -(1 -θ) ∞ 0 ũ(m)mdm + γ(0) ũ L 1 ≤ ∞ 0 ũ(m)mdm;
and finally if

∞ 0 ũ(m)mdm < γ(0) 1-θ ũ L 1 , then ∞ 0 Φ 1 (ũ, x)(m)mdm ≤ θ γ(0) 1 -θ ũ L 1 + γ(0) ũ L 1 = γ(0) 1 -θ ũ L 1 . (iv) ∞ γ(0) Φ 1 (ũ, x)(m) γ -1 (m) β dm = ∞ 0 ∞ 0 1 m f (a, m, x)ũ(e λa m) exp λa - a 0 f (s, e λ(a-s) m, x) dadm ≤ f ∞ ∞ 0 ∞ 0 1 m β ũ(e λa m) exp λa -min(f )a dadm,
and the change of variable y = e λa m,

= f ∞ ∞ 0 ∞ m 1 λm 1+β ũ(y) exp -min(f ) 1 λ ln y m dydm = f ∞ λ ∞ 0 ∞ m
m min(f )/λ-1-β ũ(y)y -min(f )/λ dydm, using Fubini's theorem and since min(f

)/λ -β > 0, = f ∞ λ ∞ 0 ũ(y)y -min(f )/λ y 0 m min(f )/λ-1-β dm = y min(f )/λ-β min(f )/λ-β dy = f ∞ λ min(f ) λ -β ∞ 0 ũ(y)y -β dy.
Finally, it is easy to check that ∞ 0 ũ(y)y -β dy ≤ γ(0) -β ũ L 1 . (v) Use Eq. (33b) and see the proof of (i).

By these estimates, we see that there exists β, C 1 , C 2 , C 3 , C 4 > 0 such that the set C × B ∈ L 1 (]γ(0), +∞[) × R, where

C := u ∈ L 1 + (]γ(0), +∞[) u L 1 ≤ 1; u ≤ C 1 γ -1 (•) a.e.; ∞ 0 u(m)mdm ≤ C 2 ; ∞ γ(0) u(m) γ -1 (m) β dm ≤ C 3 and B := [-C 4 , +C 4 ]
, is stable by the operator Φ.

In order to apply Schauder's fixed-point theorem, we will need Lemma 7. Grant Assumptions 1 and 3. C is convex, closed and compact for the weak topology σ(L 1 , L ∞ ).

Proof. It is easy to verify that C is convex. Since C is convex, if suffices to show that it is strongly closed to show that it is weakly closed. Let u n be a sequence of elements of C which converge strongly to u ∈ L 1 (]γ(0), +∞[). By the strong convergence, u L 1 ≤ 1. We can extract a subsequence u n k such that u n k converges to u a.e. Taking the pointwise limit, we have that u ≤ C 1 γ -1 (•) a.e. Furthermore, by Fatou's lemma,

∞ γ(0) u(m)mdm ≤ lim inf k→+∞ ∞ γ(0) u n k (m)mdm ≤ C 2 and ∞ γ(0) u(m) γ -1 (m) β dm ≤ lim inf k→+∞ ∞ γ(0) u n k (m) γ -1 (m) β dm ≤ C 3 .
Hence, C is strongly closed.

To show that C is weakly compact, we will show that a.

sup u∈C u L 1 < ∞, b. ∀ǫ > 0, ∃R > 0 such that ∞ R u(m)dm < ǫ for all u ∈ C , c. C is equi-integrale, i.
R > 0, ∞ R u(m)dm ≤ 1 R ∞ 0 u(m)mdm ≤ C 2 R .
To show (c.), let us first observe that for all δ 1 > 0,

γ(0)+δ 1 γ(0) u(m)dm ≤ γ -1 (γ(0) + δ 1 ) β ∞ γ(0) u(m) γ -1 (m) β dm ≤ γ -1 (γ(0) + δ 1 ) β C 3 . For any ǫ > 0, let us choose δ 1 > 0 such that γ -1 (γ(0) + δ 1 ) β C 3 ≤ ǫ 2 . Then, for all Borel set A ⊂ R + with |A|≤ δ, A u(m)dm ≤ γ(0)+δ 1 γ(0) u(m)dm + A\[0,γ(0)+δ 1 ] u(m)dm ≤ ǫ 2 + δ C 1 γ -1 (γ(0) + δ 1 )
.

Hence, we can choose δ = min δ 1 , ǫγ -1 (γ(0)+δ 1 ) 2C 1 and (c.) is verified. By Dunford-Pettis theorem, C is weakly relatively compact. Finally, since C is weakly closed, C is weakly compact.

We can now give an alternative proof of the existence of stationary solutions to (1) for arbitrary connectivity strength ε: Proof of Theorem 2 (i). We verify that the operator Φ is weakly continuous: For any sequence

(u n , x n ) → (u, x) in C × R and for any ϕ ∈ L ∞ (R + ), (Φ 1 (u n , x n ) -Φ 1 (u, x))ϕ(m)dm ≤ Q n 1 + Q n 2 + Q n 3 ,
where

Q n 1 := ∞ 0 ∞ 0 (u n (e λa m) -u(e λa m))ϕ(γ(m))e λa f (a, m, x)e -a 0 f (τ,e λ(a-τ ) m,x)dτ dadm , Q n 2 := ϕ ∞ ∞ 0 ∞ 0 u n (e λa m)e λa |f (a, m, x) -f (a, m, x n )|e -a 0 f (τ,e λ(a-τ ) m,x)dτ dadm, Q n 3 := ϕ ∞ ∞ 0 ∞ 0 u n (e λa m)e λa f (a, m, x n ) e -a 0 f (τ,e λ(a-τ ) m,x)dτ -e -a 0 f (τ,e λ(a-τ ) m,xn)dτ dadm.
Making the change of variable, ydy = e λa mdm in Q 1 we get

Q n 1 = ∞ 0 (u n (y) -u(y)) ∞ 0 ϕ(γ(ye -λa ))f (a, ye -λa , x)e -a 0 f (τ,e -λτ y,x) dadm . Since u n converges to u in σ(L 1 , L ∞ ) and ∞ 0 ϕ(γ(ye -λa ))f (a, ye -λa , x)e -a 0 f (τ,e -λτ y,x)dτ da ≤ ϕ ∞ ∞ 0 f (a, ye -λa , x)e -a 0 f (τ,e -λτ y,x)dτ da = ϕ ∞ , Q n 1 converges to 0. On the other hand, since f is bounded and Lipschitz, Q n 2 , Q n 3 ≤ u n L 1 C|x n -x|≤ C|x n -x|. Whence, Φ 1 is a continuous operator with to respect the weak topology σ(L 1 , L ∞ ).
The continuity of Φ 2 is shown analogously, taking ϕ = h (h is a bounded). Since C is stable by Φ, convex and weakly compact (Lemma 7), we can apply Schauder's fixedpoint theorem to obtain the existence of a fixed-point, which gives the existence of a stationary solution.

Corollary 2. Grant Assumptions 1 and 3. If f is of class C k , then u(m) is a function of class C k for all m > γ(0). Consequently, the stationary solutions of (1)

are of class C k . Proof. If (u, x) is a fixed-point of Φ, then u(m) = ½ m>γ(0) (γ -1 ) ′ (m) ∞ 0 f (a, γ -1 (m), x)u(e λa γ -1 (m))
exp λa -a 0 f (s, e λ(a-s) γ -1 (m), x)ds da. [START_REF] Richardson | Dynamics of populations and networks of neurons with voltageactivated and calcium-activated currents[END_REF] Making the change of variable y = e λa γ -1 (m) in a, as in the estimate (ii) of Lemma 6, we obtain

u(m) = ½ m>γ(0) (γ -1 ) ′ (m) λγ -1 (m) ∞ γ -1 (m) f (g(y, m), y, x)u(y) exp - g(y,m) 0 f (s, e s y, x)ds dy, (35)
where g(y, m) = ln y λ(γ -1 (m)) . We conclude with a bootstrap argument: if u is L 1 , then the right hand side of Eq. ( 35) is a continuous function of m, meaning that u is continuous. But if u is continuous, then the right hand side is of class C 1 , etc. Corollary 3. Grant Assumptions 1 and 3. There exist a constant C > 0, such that the stationary solution ρ ∞ satisfies,

ρ ∞ (a, m) ≤ Ce -σ(a-∆ abs ) m . (36)
Proof. From the previous theorem it follows that there is C such that u(m) ≤ C/m, which together with [START_REF] Pozzorini | Automated high-throughput characterization of single neurons by means of simplified spiking models[END_REF] implies

ρ ∞ (a, m) ≤ C e -a 0 f (s,e λ(a-s) γ -1 (m),x)ds m . (37)
The estimate follows from Assumption 3 (i).

Formula in the case of short-term synaptic depression

In general, there is no explicit formula for the invariant probability measure solving Eq. ( 26). However, in the case of short-term synaptic depression Eq. ( 3), we can derive an explicit expression for the total postsynaptic potential

X(x) := ∞ 0 ĥ(t) 1 0 ∞ 0 (1 -m)f (a, x)ρ x ∞ (a, m)dadmdt, (38)
for any x ∈ R. This fact has been reported in the theoretical neuroscience literature [START_REF] Romani | Mean-field analysis of selective persistent activity in presence of short-term synaptic depression[END_REF]; we provide here a rigorous and analytic justification for it. For all x ∈ R, let us introduce the quantities

I x := ∞ 0 af (a, x) exp - a 0 f (s, x)ds da = ∞ 0 exp - a 0 f (s, x)ds da, P x(λ) := ∞ 0 e -λa f (a, x) exp - a 0 f (s, x)ds da.
I x can be interpreted as the mean inter-spike interval of a neuron receiving a constant input x. P x(λ) can be seen as the Laplace transform of the inter-spike interval distribution of that neuron, evaluated in λ.

Proposition 2. Grant Assumptions 1 and 3. For all x ∈ R,

X(x) = ∞ 0 ĥ(t)dt 1 I x 1 -P x(λ)
1 -υP x(λ) .

Proof. Using the method of characteristics (i.e. combining Eqs. (26b) and ( 32)), we have On the other hand, .

Finally, we have

X(x) = ∞ 0 ĥ(t)dt

Exponential stability in the weak connectivity regime

To study the long time behavior of Eq. (1) in the weak connectivity regime, we proceed by a perturbation of the non-interacting case [START_REF] Hairer | Yet another look at Harris' ergodic theorem for Markov chains[END_REF], taking x = εx ∞ , where x ∞ is given by the unique stationary solution to Eq. (1) when ε ∈] -ε * , +ε * [ (ε * is taken from Theorem 2 (ii)). In this section, we keep the small ε fixed and we work under Assumptions 1, 3 and 5. We roughly follow the same line of argument as [START_REF] Mischler | Weak and strong connectivity regimes for a general time elapsed neuron network model[END_REF]Sec. 5].

For convenience, we first rewrite (1) in a more formal and compact form:

∂ t ρ t = -∂ a ρ t + λ∂ m (mρ t ) -f (εx t )ρ t + δ a 0 (γ • Π) * f (εx t )ρ t , (39a)
x t = t 0 h(t -s)f (εx s )ρ s dadm ds, (39b)

ρ 0 = u 0 , (39c)
where δ a 0 indicates that (singular) mass enters in a = 02 , Π : (a, m) → m is the projection on m and * denotes the pushforward measure. To write Eq. ( 39) as an evolution equation, we introduce an auxiliary transport equation on R + × R + × R * + ∂ t ζ t = -∂ s ζ t + δ s 0 f (εx t )ρ t , ζ 0 = 0, Also, we have, for all t ≥ 0,

Z (1) t L 1 ≤ |ε|2L f ρ t L 1 |x t -x ∞ |≤ |ε|2L f ζt L 1 (µ) , (45a)
Z (1) t L 1 (w) ≤ |ε|2L f ρ t L 1 (w) |x t -x ∞ |≤ |ε|2L f u 0 L 1 (w)+ b α ζt L 1 (µ) , (45b)
Z (2) t L 1 (µ) ≤ |ε| h ∞ L f ρ t L 1 |x t -x ∞ |≤ |ε| h ∞ L f ζt L 1 (µ) , (45c)
where we have used Theorem 1 (ii) in the first line and Lemma 1 in the second.

Lemma 8. Grant Assumptions 1, 3 and 5 and take (ρ 0 , ζ0 ) as in Eq. [START_REF] Robert | Short-term synaptic plasticity[END_REF]. There exists K 1 ≥ 1 and a 1 > 0 such that, for all initial data u 0 ∈ L 1 + (w) with u 0 L 1 = 1, S Λ t (ρ 0 , ζ0 )

L 1 (w)×L 1 (µ)
≤ K 1 e -a 1 t (ρ 0 , ζ0 ) L 1 (w)×L 1 (µ) , ∀t ≥ 0.

If in addition, we grant Assumption 4, then there exists K 2 ≥ 1 and a 2 > 0 such that, for all initial data u 0 ∈ L 1 + with u 0 L 1 = 1, S Λ t (ρ 0 , ζ0 )

L 1 ×L 1 (µ)
≤ K 2 e -a 2 t (ρ 0 , ζ0 ) L 1 ×L 1 (µ) , ∀t ≥ 0.

(47)

Proof. We write (S Λ t (ρ 0 , ζ0 ) (1) , S Λ t (ρ 0 , ζ0 ) (2) ) := S Λ t (ρ 0 , ζ0 ) the first and second component of S Λ t (ρ 0 , ζ0 ). By Theorem 4, there exists K ≥ 0 and a > 0 such that,

S Λ t (ρ 0 , ζ0 ) (1)
L 1 (w)

≤ Ke -at ρ0 L 1 (w) , ∀t ≥ 0.

Then, S Λ t (ρ 0 , ζ0 ) (2) L 1 (µ) = t 0 f (εx ∞ )S Λ t-s (ρ 0 , ζ0 ) (1) L 1 C h e -hs ds

+ ∞ t ζ0 (s) L 1 C h e -hs ds ≤ C h f ∞ K t 0
e -a(t-s) e -hs ds ρ0 L 1 (w) + e -ht ζ0 L 1 (µ) .

Gathering the bounds on the two components and observing that the function t → t 0 e -a(t-s) e -hs ds decays exponentially, we can conclude that there exists K 1 ≥ 1 and a 1 > 0 such that Eq. (46) holds.

For Eq. (47) we use Theorem 5 and follow the same argument.

We can now prove our main result:

Proof of Theorem 3. By Duhamel's formula (44), Eq. ( 46) in Lemma 8 and the bounds Eqs. (45), for all t ≥ 0, (ρ t , ζt ) L 1 (w)×L 1 (µ) ≤ S Λ t (ρ 0 , ζ0 )

L 1 (w)×L 1 (µ) + t 0 S Λ t-s (Z (1) s , Z (2) s )

L 1 (w)×L 1 (µ) ds ≤ K 1 e -a 1 t (ρ 0 , ζ0 ) L 1 (w)×L 1 (µ) + K 1 t 0 e -a 1 (t-s) (Z (1) s , Z (2) s )

L 1 (w)×L 1 (µ) ds ≤ K 1 e -a 1 t (ρ 0 , ζ0 ) L 1 (w)×L 1 (µ) + |ε| CW t 0 e -a 1 (t-s) (ρ s , ζs ) L 1 (w)×L 1 (µ) ds =: Q(t),

Figure 1 :

 1 Figure 1: Simulations of Eq.(2) showing self-sustained bursts for strong connectivity strength (ε ≫ 0) and relaxation to a stationary state for weak connectivity strength (small ε). Simulations of a network of 5 • 10 5 adaptive SRM 0 neurons approximating Eq. (2) with identical parameters (except for ε) and identical initial conditions. The raster plots below the plots for the time-evolution of the total postsynaptic potential x t represent the spikes of 100 randomly selected neurons.

e

  -λa m)f (a, e -λa m, x) exp -a 0 f (s, e -λs m, x)ds dadm≤ -λa f (a, e -λa m, x) exp -a 0 f (s, e -λs m, x)ds da =:ϑ(m) dm + γ(0) ũ L 1 . There exists θ ∈]0, 1[ such that for all m ∈ R * + , ϑ(m) < 1: Fix ǫ > 0.ϑ(m) ≤ ǫ 0 f (a, e -λa m, x) exp -a 0 f (s, e -λs m, x)ds da + ∞ ǫ e -λǫ f (a, e -λa m, x) exp -a 0 f (s, e -λs m, x)ds da = 1 -(1 -e -λǫ ) ∞ ǫ f (a, e -λa m, x) exp -a 0 f (s, e -λs m, x)ds da = 1 -(1 -e -λǫ ) exp -ǫ 0 f (s, e -λs m, x)ds

  e. ∀ǫ > 0, ∃δ > 0 such that for all Borel set A ⊂ R + with |A|≤ δ and for all u ∈ C , A u(m)dm ≤ ǫ, and use Dunford-Pettis theorem. (a.) is clearly verified. (b.) is also verified since for all

0 ½

 0 e λa m<1 ρ x ∞ (0, e λa m) exp λa -

0 ½e

 0 e λa m<1 mf (a, x)ρ x ∞ (0, e λa m) exp λa --λa mf (a, x)ρ x ∞ (0, m) exp -, x)ρ x ∞ (a, m)dadm = P x(λ)(1 -υ) I x(1 -υP x(λ))

f 1 I x 1 - 1 -

 111 (a, x)ρ x ∞ (a, m)dadm -P x(λ)(1 -υ) (1 -υP x(λ)) P x(λ)1 -υP x(λ) .

To be precise, they use Doeblin's theorem.

δ a 0 should not be confused with the Dirac distribution δ0=a. Using δ0=a, by integration by parts of weak solutions, Eq. (39a) should write ∂tρt = -∂aρt + λ∂m(mρt) -f (εxt)ρt + δ0=a (γ • Π) * f (εxt)ρtρt(0, •) .

Acknowledgements

We thank Stéphane Mischler for supervising this work and Wulfram Gerstner for many useful discussions. This research has been funded by the European Union's Horizon 2020 research and innovation programme under the Marie Skłodowska-Curie grant agreement No 754362 and the Swiss National Science Foundation (grant no. 200020_184615).

Lemma 5 (Lipschitz continuity at T = ∞). Grant Assumptions 1 and 3. Writing ρ x ∞ ∈ L 1 + (w) the invariant probability measure given by Theorem 4 for any x ∈ R, the function

is Lipschitz and there exists C > 0 such that

Proof. Since f is lipschitz in x, we have, for any

from where we only need to bound the first term on the right hand side. We can use Theorem 4 and Lemma 4: for any T ∈ R + ,

where K and a are the exponential stability constants of Theorem 4. Choosing T such that Ke -aT = 1/2, we get

Gathering the bounds concludes the proof.

Theorem 6 (Stationary solutions). Grant Assumptions 1 and 3. We have (i) There exists a stationary solution to (1).

(ii) There exists ε * > 0 such that for all ε ∈] -ε * , +ε * [, the stationary solution to (1) is unique.

Proof. For all x ∈ R, let us write ρ x ∞ ∈ L 1 + (w) the unique invariant measure given by Theorem 4 and let us also take the function Υ from Lemma 5. By Corollary 1,

Hence, the study of the existence and the uniqueness of stationary solutions is reduced to the study of the existence and the uniqueness of the fixed-point of Υ.

Since for all

Then, the continuity of Υ guarantees the existence of a fixed-point, which proves (i).

To obtain (ii), we observe that the Lipschitz constant of Υ is |ε|C: if we take |ε|< ε * := C -1 , Υ is a contraction and we can apply Banach's fixed-point theorem to conclude.

Alternative proof for the existence using Schauder's fixed-point theorem

We include here an alternative proof for the existence of a stationary solution, which is interesting for two reasons: on the one hand, it does not rely on the Harris-Doeblin method, and on the other hand, it provides some estimates on the stationary solutions. which solution is given by the method of characteristics:

Using the auxiliary equation, Eq. ( 39) is equivalent to

where

Now, we write Eq. ( 40) as the sum of a linear equation and a perturbation:

t ), (42a)

where

Then, using Eq. ( 41), by the linearity of the operator Λ and writing ρt := ρ t -ρ ∞ and ζt := ζ t -ζ ∞ , we get

t , Z

t ), (43a)

Writing (S Λ t ) t∈R + the semigroup associated with the operator Λ, we have, by Duhamel's formula,

Let us define the weighted space

Note that, for all t ≥ 0,

where CK is a constant depending on W . We have, for all t ≥ 0,

Whence, by Grönwall's lemma,

For all t ≥ 0, we have

Whence, choosing ε * * W := a 1 CW ∧ ε * , we easily see that there exists C ≥ 1 and c W > 0 such that Eq. ( 6) holds.

For Eq. ( 7) we use Eq. (47) instead of Eq. ( 46) and follow the same argument.


V Bonnet

J Lagarde

P Fraisse

N Ramdani

S Ramdani

P Poignet

B Bardy

Modeling human postural coordination to improve the control of balance in humanoids

This paper models recent data in the field of postural coordination showing the existence of self-organized postural states, and transition between them, underlying supra-postural tracking movements. The proposed biomechanical model, capitalizing on stability and optimization criteria, captures the complex postural behaviors observed in humans and can be used to implement efficient balance control principles in humanoids. Experimental results on humanoids show the relevance of this work.

I. INTRODUCTION

The field of movement science, perception-action and motor control in humans provides a natural source of inspiration for researchers in humanoid robotics. Due to their anthropomorphic structure, humanoid robots often present dynamic similarities with humans that are interesting to study in both fields of research. The goal of the present paper is to model basic phenomena characterizing human postural coordination, which have motivated in depth analyses during the last decade in the field of human movement science. In 1985, Nasher and Mc Collum [START_REF] Nashner | The organization of postural movements: a formal basis and experimental synthesis[END_REF] observed two postural strategies when the whole body oscillates in upright position. These postural strategies are characterized by a large activity of the ankles during the the ankle strategy, and by the coordinative movements of the hips and the ankles in the hip strategy. Within the framework of coordination dynamics [START_REF] Kelso | Dynamics Patterns, The Self-Organization of brain and behavior[END_REF], Bardy et al. [START_REF] Bardy | Postural coordination modes considered as emergent phenomena[END_REF] proposed the use of a collective variable to describe in a simple way the complex biological couplings underlying the ordering in space and time of postural coordination. To identify this collective variables, Bardy et al. [START_REF] Bardy | Dynamics of Human Postural Transitions[END_REF] used a parametric visual tracking task (Fig. 1). They had standing participants moving back and forth in the sagittal plane in order to track the displacement of a the target, with the knees locked and the soles in permanent contact with the ground. This simple task allowed the observation of several self-organized properties of the postural system, such as phase transition, multistability, self-organization, critical fluctuations, hysteresis, and critical slowing down. The collective variable able to capture both fully and in a very compact way these properties was the phase difference between the ankle and the hip.

These experimental observations suggest that posture may be organized in terms of dynamical principles [START_REF] Kelso | Dynamics Patterns, The Self-Organization of brain and behavior[END_REF], and call for the existence of common principles governing pattern formation in complex biological systems. An attempt to reproduce these complex couplings should ameliorate the control of balance in humanoids. One first approach reaching for this global objective is to ensure the validity of the observed principles in humanoid robots.

A biomechanical model using a double inverted pendulum associated with an under-constraints optimization method was proposed by [START_REF] Martin | Optimization model predictions for postural coordination modes[END_REF] to reproduce the dynamical and mechanical aforementioned observations. Furthermore the objective function used in [START_REF] Martin | Optimization model predictions for postural coordination modes[END_REF] was also bio-inspired. [START_REF] Engelbrecht1 | Invariant characteristics of horizontalplane minimum-torque-change movements with one mechanical degree of freedom[END_REF] has shown that during stance the human body aims at minimizing the variations in energy during tracking tasks. However, this model cannot describe critical slowing down and hysteresis, two important phenomena observed in human behaviors and characterizing non linear systems. The present work has two distinct objectives:

1) develop a biomechanical model reproducing the human observations, 2) implement the model in humanoid robots. The implementation of the model permits the execution of pre-programmed motion even if balance cannot be maintained, contrary to humans. Humanoid robot sensors provide information to decompose and understand a given task perfectly.

For this purpose, the target tracking task used by [START_REF] Bardy | Postural coordination modes considered as emergent phenomena[END_REF] and [START_REF] Bardy | Dynamics of Human Postural Transitions[END_REF] was chosen, and the Martin et al.s optimization model [START_REF] Martin | Optimization model predictions for postural coordination modes[END_REF] for this task was examined in the context of humanoid robotics.

An analysis of the natural frequencies was performed and compared with the simulations resulting from the combination of an optimization and a biomechanical model. Finally, in order to validate the resulting model, the optimal joint amplitudes were executed in open loop on two different humanoid robots.

II. BIOMECHANICAL MODEL FOR THE TRACKING TASK

A. Biomechanical modeling

Here we replicate the simple biomechanical model proposed by Martin et al. [START_REF] Martin | Optimization model predictions for postural coordination modes[END_REF] for postural sway in the sagittal plane. The model (Fig. 2) is formed by three rigid segments including the feet, the lower limbs and the upper part of the body which are linked by the ankle and the hip joints (such joints are modeled as frictionless hinges). The entire sole of the feet were assumed to stay in contact with the floor. For the inverse dynamic model calculation, the Lagrange movement equations were used. We changed only the joint markers, with respect to the original Martin's model.

Γ 1 = θ1 (m 1 k 2 1 l 2 1 + l 2 1 m 2 + I 1 + I 2 ) (1)
+ θ2 (m 2 l 1 l 2 k 2 cos(θ 1 -θ 2 ) + I 2 ) + θ2 1 (m 2 l 1 l 2 k 2 sin(θ 1 -θ 2 )) + gl 1 sin(θ 1 )(m 1 k 1 + m 2 ) Γ 2 = θ2 (m 2 k 2 2 l 2 2 + I 2 ) (2)
+ θ1 (m 2 l 1 l 2 k 2 cos(θ 1 -θ 2 ) + I 2 ) + θ2 1 (m 2 l 1 l 2 k 2 sin(θ 1 -θ 2 )) + gm 2 k 2 l 2 sin(θ 2 )
where, Γ 1 , Γ 2 represent the torque at the ankle and hip respectively. θ 1 , θ1 , θ1 are position, velocity and acceleration of joint 1. Finally, the dynamic parameters of each segment are mass m, lenght l, center of mass position k, inertia I, and the gravity constant g .

Balance is described by the position of center of pressure (CoP) within the BoS the sagittal plane. It can be expressed as a function of the dynamic parameters (eq.3).

CoP

= (Γ 1 -F gx d + m 0 k 0 g) F gy (3)
where F gx is the horizontal ground reaction force and F gy the vertical ground reaction force. For the calculation of the ground reaction forces, Cahouet et al. [START_REF] Cahouet | Static optimal estimation of joint accelerations for inverse dynamics problem solution[END_REF] proposed to use Euler's equations.

B. Natural frequencies of the double inverted pendulum

Intrinsic properties of the dynamic model can explain some phenomena observed in experiments, around specific frequencies. Indeed, each mechanical system has natural frequencies at which minimal inputs can provide maximal reactions. This interesting fact has been studied in other research areas such as in humanoid structure construction [START_REF] Williamson | Robot Arm Control Exploiting Natural Dynamics[END_REF]. Natural frequencies are calculated via the linearization of the inverse dynamic model, taking only the homogeneous equations of the system (4).

θ1 α 1 + θ2 β 1 + θ 1 γ 1 = 0 θ1 α 2 + θ2 β 2 + θ 1 γ 2 = 0 (4)
and α, β, γ are coefficients that represent the dynamic parameters; they are defined by the following equations:

α 1 = (m 1 k 2 1 l 2 1 + l 2 1 m 2 + I 1 + I 2 ) β 1 = (m 2 l 1 l 2 k 2 + I 2 ) γ 1 = gl 1 (m 1 k 1 + m 2 ) α 2 = (m 2 l 1 l 2 k 2 + I 2 ) β 2 = (m 2 k 2 2 l 2 2 + I 2 ) γ 2 = gm 2 k 2 l 2
By considering the periodic solution of (eq. 4), one can solve the system for amplitudes. The associated characteristic equation has two real positive solutions which are the natural frequencies of the system.

ω 01 = -1 2(α 1 β 2 -β 1 α 2 ) -2(α 1 β 2 -β 1 α 2 )(-γ 1 β 2 -α 1 γ 2 + (γ 1 β 2 -α 1 γ 2 ) 2 + 4 γ 1 γ 2 β 1 α 2 ) 1/2 (5)
ω 02 = -1 2(α 1 β 2 -β 1 α 2 ) √ 2(α 1 β 2 -β 1 α 2 )(γ 1 β 2 +α 1 γ 2 + (γ 1 β 2 -α 1 γ 2 ) 2 + 4 γ 1 γ 2 β 1 α 2 ) 1/2 (6)
The numerical values of the natural frequencies for the two tested robots are:

• HRP2: ω 01 = 0.48 Hz and ω 02 = 1.05 Hz • HOAP3: ω 01 = 0.74 Hz and ω 02 = 1.65 Hz The system response shows two resonance peaks at the natural frequencies. This observation tends to confirm the privilege for certain frequencies to naturally decrease the joint power necessary to accomplish the movement, and to increase the stability of the system since CoP position is linked to the ankle torque (eq. 3).

C. Optimal amplitude calculation for the tracking task

In Bardy et al.'s experiments [START_REF] Bardy | Postural coordination modes considered as emergent phenomena[END_REF] [4], the head movement was assumed to be periodic with a constant distance between the head and the target. It thus makes sense to use a periodic solution for the system. Martin et al. proposed to decompose the joint positions in Fourier series.

θ 1 (t) = a 01 2 + N k=1 (a k1 cos(kωt) + b k1 sin(kωt)) (7)
where a k1 and b k1 are the ankle amplitude coefficients for the k th harmonic, ω the pulsation of the target motion and N the number of harmonics in Fourier series.

The system constraints are thus to keep a constant distance between the head and the sinusoidal target (eq.8) and to hold the orthostatic balance (eq.9), i.e, to maintain the CoP inside the base of support (BoS).

A t = l 1 (cos(θ 1 ( T t 2 )) -cos(θ 1 (0)))+ l 2 (cos(θ 1 ( T t 2 ) + θ 2 ( T t 2 )) -cos(θ(0) + θ 2 (0))) (8)
CoP min ≤ CoP ≤ CoP max (9)
where Tt 2 is the half period of the target motion, A t its amplitude, CoP max and CoP min the upper and the lower BoS limits, respectively. In the litterature, it is often assumed [START_REF] Engelbrecht1 | Invariant characteristics of horizontalplane minimum-torque-change movements with one mechanical degree of freedom[END_REF] that the human body tries to minimize the global torque variations. For this reason, (eq.10) is named the minimum torque change criterion.

J n = n j=1 d Γ 2 1 dt (t j ) + d Γ 2 2 dt (t j ) Δt ( 10
)
where n is the number of discrete time t j and Δt the time discretization step, and 10 is the criterion in discrete time.

The optimization problem boils down to find the optimal Fourier coefficients that respect these constraints.

III. SIMULATION OF THE TRACKING TASK

A. Simulation parameters

For the tracking task simulation, we used two different models of humanoid robot, HRP2 and HOAP3. HRP2 (1.54m, 58kg) has the mass and the length repartition similar to the Winter's anthropometric model [START_REF] Winter | Biomechanics and motor control of human movement third edition[END_REF]. HOAP3 however is not close to the anthropometric tables (70cm, 8.8kg). It is interesting to see how the optimization evolves with a nonstandard anthropometric model. The initial posture parameters all positioned the double inverted pendulum in an upright position. We took in (eq. 7) only the first Fourier harmonic, to keep a periodic solution. Simulations were performed with different input parameters.

The amplitude of the target varied between 20 and 100 mm, and the oscillation frequencies varied from 0.2Hz to 1.2Hz. The length of the BoS must also be a variable parameter, because in the model the detachment of the sole is not considered. The humanoid feet are rigid, wich can introduce a bouncing phenomenon. It is necessary to adapt the length of the BoS. The difference in the phase between the two joints was measured by using the Lissajous curves. For the model corresponding to HRP2 parameters, two distinctive stable modes appeared (Fig. 3a), similar to the outcome of the anthropometric model used by Martin: the inphase mode between the ankle and the hip at low frequencies and the anti-phase mode at high frequencies. Around the natural frequency ω 01 , and independently of the length of the BoS and A t , a minimum peak-to-peak variation of the CoP was visible (Fig. 3b) with a reduction of the total joint power used (Fig. 3c).

B. Simulation results

For the model corresponding to the HOAP3 parameters, two distinctive stable modes again appeared (Fig. 4a), but an antiphase mode was also the preferred mode at very low frequencies, probably due to the non-antropometric parameters. With this exception, all previous observations were confirmed. It is interesting to note that a modification of the inertia and CoM position of each segment, to approach the Winter's parameters [START_REF] Winter | Biomechanics and motor control of human movement third edition[END_REF], produced results closer to the HRP2 simulations.

(Fig. 5), (Fig. 6) illustrate Fourier's first optimal coefficients (i.e the amplitude) for fixed length of BoS and A t . These curves, constant before the transition, can be approximated to provide a direct relation between the frequency and the optimal amplitude.

Moreover the size of the constraints (BoS and A T ) influenced the frequency of phase transition. There is converging evidence from these results that the CoP governs the emergence of the coordination mode; when CoP approaches the BoS limit, the hip-ankle relative phase switches to a new value. The anti-phase mode thus appears to be the most stable mode in terms of balance, whereas the in phase mode appears to be the mode minimizing energy. These results also tend to show that the natural frequency is a privileged frequency at which the system needs minimum joint power and has a better equilibrium(eq. 3). The second natural frequency has no physical sense since the constraints are activated and change the model used for calculation of the natural frequencies. The second minimum on the joint power curves, located just after the transition, is due to the ankle torque. In fact, the torque produced by the in phase mode near the transition is important (in order to maintain balance), which is not the case after the transition. In conclusion, we observed a decrease in the peak-to-peak CoP movement during the anti-phase mode, a reduction in the total joint power for the in phase mode, and an increase in balance with a minimization of energy around the natural frequencies.

IV. EXPERIMENTAL RESULTS AT THE TRACKING TASK

To show the effectiveness of the proposed optimized biomechanical model, we implemented in HOAP3 and HRP2 the offline optimal joint amplitudes calculated for each frequency of the target motion. The robots were used without stabilizers and positioned initially in an upright position. The main features we expected to observe were the minimization of the joint energy during the in phase mode, the minimization of the CoP peak-to-peak displacement during the anti-phase mode, the shaping role of the CoP location in the phase transition, and the CoP and torque values around the natural frequencies.

A. Effects of the coordination mode on the CoP

To analyse the influence of the two coordination modes on the CoP position, we executed the same head motion at various frequencies. The sign of the hip coefficient was changed and consequently the coordination bifurcated from in phase to antiphase. On (Fig. 7) we can see clearly that the anti-phase mode reduced the CoP peak-to-peak displacement of the HOAP3 robot.

CoP evolution of HRP2 (Fig. 8) described the same phenomenon.

The anti-phase mode thus appears to improved the general equilibrium in both robots.

B. Effects of the coordination mode on the energy

In this part, only the HRP2 results are presented since the robot HOAP3 has no joint torque sensors. The methodology was identical to the one used for the CoP analysis.

The in-phase mode (Fig. 9a,b) dissipated less joint power than the anti-phase mode (Fig. 9c,d). This was due to the hip joint which increased after the phase transition. This shows clairly that the in-phase mode minimizes joint power.

C. The role of CoP position in phase transition

In section III, we evidenced that the CoP dynamics governs both the general equilibrium and the coordination mode adopted. We tried to evidence this effect with humanoid robots by observing when the CoP location crossed the BoS limits. Because however we did not take into account a possible takeoff of the sole from the ground, rebound behavior sometime appeared due to the rigidity of the humanoid feet. Thus, robots could fall because of these uncontrolled rebounds. We thus searched for a compromise between the length of the BoS and the amplitude of target motion. The experiment was performed with short lengths of BoS. For example, with HOAP3, a BoS length of 1cm together with a head amplitude of 5cm, produced a stable value of the in phase mode at 0.85Hz while the predicted value was 0.92Hz. As an other example using HRP2 (BoS = 2cm and A T = 4cm ), the experimental in phase stable value was around 0.55Hz while the predicted value was 0.69Hz. The difference between experimental and predicted results was more important for HRP2 than for HOAP3, probably because HRP2 has a compliance element located at is ankle level. This compliant element can add oscillations which are not considered in our model. The differences can also be due to an incomplete model: the missing friction parameter, among others, may reduce these differences. Above these values the only stable mode was the anti-phase mode. This confirms that the CoP is an important variable constraining the adopted coordination pattern.

D. Robot postural states around natural frequencies

Section III evidenced privileged frequencies, with minimal variations in peak-to-peak CoP and energy.

HRP2 sensors provide the joint torques as function of frequency. Our experiments did not reveal a direct effect of natural frequencies, but a slow down in the increasing of torques with frequency near the theoretical natural frequency. Further tests are necessary in order to fully understand the effects of natural frequencies.

For the CoP, the test is performed with HOAP3. A diminution of the peak-to-peak CoP displacements is visible (Fig. 10), but less important and not at the theoretical frequency. That is certainly due to the imprecisions of the model (section .IV-C). Nevertheless, the important fact is the presence of a privileged frequency.

V. DISCUSSION AND PERSPECTIVES By using humanoid robots we validated the results obtained by Bardy et al. [START_REF] Bardy | Dynamics of Human Postural Transitions[END_REF]. We showed experimentally, that the antiphase postural mode is more efficient for the global stability of the postural system at high frequencies, and that the in phase mode minimizes joint power at low frequencies. Phase transition occurred when CoP reached the limits of the BoS. This proves that the distance between the CoP and the limits of the BoS is a key variable which has to be taken into account by the nervous system. The natural frequency for in-phase mode appears like a privileged mode in terms of CoP displacements and in terms of joint power. Based on these first results, we aim at developing a ankle-hip controller which takes into account the natural frequencies. For example, we would like to generate adequate trajectories in order to maintain the frequency spectrum close to the natural frequency. However, the model analysed in this paper does not take into account all the dynamical phenomena observed in humans. For this reason, and also in order to improve the humanoid robot control, we plan to develop a more complete model of the human postural system, capitalizing on nonlinear coupled oscillators (ankles and hips). Haken Kelso and Bunz (HKB) have proposed a nonlinear coupled-oscillators system for bimanual coordinations [START_REF] Haken | A theorical model of phase transitions in human hand movements[END_REF]. The HKB model is able to produce bi-stability, phase transition, amplitude evolution with frequency increase. With few changes, such as those proposed by [START_REF] Fuchs | The HKB model revisited: How varying the degree of symmetry controls dynamics[END_REF], we plan to have a more complete picture of human postural coordination, in order to create a new controller for the ankles and the hips.

The use of nonlinear coupled oscillators associated with the results obtained by the biomechanical model will be proposed to create a new controller between the ankle and the hip. This controller will be able to automatically change the phase difference, the joint amplitude necessary to maintain balance, capitalizing on the natural frequency properties.

Fig. 1 .

 1 Fig. 1. Task and set-up used in Bardy et al. experiments. Participants faced a white target on a large screen and tracked its motion in the sagittal plane. The frequency of the target increased or decreased continuously.

Fig. 2 .

 2 Fig. 2. Biomechanical model in the sagittal plane used to describe human postural coordination

Fig. 3 .

 3 Fig.3. Typical simulations (a T = 4cm , BoS = 2cm) for the HRP2 model. In-phase (f < 0.67) and anti-phase (f > 0.67) coordination modes (a), CoP under the BoS limits for the in-phase mode with a minimum variation at ω 0 and on the BoS limits after the phase transition (b), the total net joint power (b) with a inflection close to ω 0 (c).

Fig. 4 .Fig. 5 .

 45 Fig.[START_REF] Bardy | Dynamics of Human Postural Transitions[END_REF]. Typical simulations (a T = 4cm, BoS = 2cm) for the HOAP3 model. Two coordination modes anti-phase (f < 0.47,f > 0.92 ) and inphase (0.47 < f < 0.92) (a), CoP with a minimum variation close to ω 0 (b), the total net joint power (b) with a inflection close to ω 0 (c), the two last curves (d) and (e) represent the optimal amplitude (first Fourier coefficients.) at each frequency

Fig. 6 .Fig. 7 .

 67 Fig.6. Optimized amplitudes for the HOAP3 model. The anti-phase mode was observed at very low frequencies (f < 0.47Hz) due to the specific HOAP3 parameters; then amplitudes followed the HRP2 amplitude pattern.

Fig. 8 .

 8 Fig. 8. CoP evolution during two periods of head oscillation for HRP2 (A T = 2cm, BoS = 2cm); the anti-phase mode reduces the CoP peak-topeak distance.

Fig. 9 .

 9 Fig.9. The joint power (W) is inferior for the in-phase mode (c,d) compared to the anti-phase mode (a,b). Ankle power is slightly the same for in phase mode (a) and anti-phase mode (d), hip power during anti-phase becomes larger in the movement.

Fig. 10 .

 10 Fig.10. Evolution of the experimental CoP peak-to-peak displacement as a function of frequency. A minimum variation observable at 0.62Hz (not depending of the BoS and the A T ) but less important than the theoretical frequency.

Fig. 11 .

 11 Fig. 11.HOAP3 experiments with in-phase coefficients at high frequency(f = 1Hz, BoS = 1cm, At = 5cm). The robot cannot maintain its balance and falls backward.

Fig. 12 .

 12 Fig. 12. HRP2 experiments with in-phase coefficients at high frequency (f = 0.6Hz, BoS = 2cm, At = 4cm). The robot cannot maintain its balance and falls backward.

ACKNOWLEDGMENT

Author Nacim Ramdani is on leave from CERTES EA 3481 Universite Paris 12 Val de Marne.