Lou Zonca

David Holcman

Analysis of bursting dynamics in a modified facilitation-depression model accounting for afterhyperpolarization

Analysis of bursting dynamics in a modified facilitation-depression model accounting for afterhyperpolarization

Lou Zonca, David Holcman

Introduction

Electrophysiological recordings of neuronal networks reveal periods of synchronous high-frequency activity called bursts separated by interbursts (quiet time periods). Bursting can either be due to intrinsic channel properties driven by Ca 2+ and/or voltage-gated channels, or by collective properties of the neuronal network [1]. Several models have been proposed to generate bursting, starting with the classical Wilson-Cowan oscillator, where two reciprocally coupled populations of excitatory and inhibitory neurons exhibit bursting [START_REF] Wilson | Excitatory and inhibitory interactions in localized populations of model neurons[END_REF][START_REF] Wilson | A mathematical theory of the functional dynamics of cortical and thalamic nervous tissue[END_REF]. Bursters are modeled as slow-fast dynamical systems using the Hodgkin-Huxley formalism, where the fast dynamics are responsible for the fast spiking and are modulated by the slow variables representing the mean voltage dynamics [START_REF] Izhikevich | Dynamical Systems in Neuroscience: the Geometry of Excitability and Bursting[END_REF]. The classical Hindmarsh-Rose model [START_REF] Hindmarsh | A model of neuronal bursting using three coupled first order differential equations[END_REF] implements such strategy with three variables: one for the membrane potential, one for the fast ion channels (fast subsystem) and one for the slow ion channels (slow subsystem). Following this model, different type of bursters have been developed, such as ones with low spike frequency at the beginning and the end of a burst [START_REF] Rinzel | Bursting oscillations in an excitable membrane model[END_REF]. Similarly, parabolic bursters exhibit fast-oscillation frequencies that vary along time with the burst [START_REF] Ermentrout | Parabolic bursting in an excitable system coupled with a slow oscillation[END_REF]. In periodic bursters, the slow variable, is oscillating periodically, between a set of stable attractors, during these transitions the fast variable exhibits spiking [START_REF] Golubitsky | Bursting in coupled cell systems[END_REF]. Bursters can be classified according to their topological bifurcation diagram, where the fast subsystem can lead to two main changes of the state space: 1) resting to spiking, when a stable equilibrium transitions to an attractive limit cycle, 2) spiking to resting for the opposite transition [START_REF] Coombes | Bursting: The Genesis Of Rhythm In The Nervous System[END_REF][START_REF] Zeldenrust | Neural coding with burstscurrent state and future perspectives[END_REF]. Bursts that emerge as a network property have been studied using different modeling approaches such as coupled integrate and fire neurons [START_REF] Brunel | Fast global oscillations in networks of integrate-and-fire neurons with low firing rates[END_REF][START_REF] Neltner | On synchrony of weakly coupled neurons at low firing rate[END_REF], improved recently by adding noise to connected Hodgkin-Huxley type neurons, to allow desynchronisation [START_REF] Chizhov | Conductance-based refractory density approach for a population of bursting neurons[END_REF]. Bursting can also depend on the balance between excitatory and inhibitory neurons: coupling excitatory neurons results in in-phase bursting within the network, whereas inhibitory coupling leads to anti-phase dynamics [START_REF] Shi | Burst synchronization of electrically and chemically coupled map-based neurons[END_REF]. Furthermore, time-delays [START_REF] Roxin | Role of delays in shaping spatiotemporal dynamics of neuronal activity in large networks[END_REF] play a crucial role in synchronisation, by generating coherent bursting in the Hindmarsh-Rose model, specifically when the time-delays are inversely proportional to the coupling strength [START_REF] Liang | Phase synchronization of inhibitory bursting neurons induced by distributed time delays in chemical coupling[END_REF]. Central Pattern generators such as the respiratory rhythm in the pre-botzinger complex [START_REF] Smith | Pre-botzinger complex: a brainstem region that may generate respiratory rhythm in mammals[END_REF][START_REF] Cui | Defining prebtzinger complex rhythm-and pattern-generating neural microcircuits in vivo[END_REF], mastication or oscillatory motor neurons [START_REF] Marder | Invertebrate central pattern generation moves along[END_REF] are involved in the genesis and maintenance of rhythmic patterns. Interestingly, several coupled pacemaker neurons receiving an excitatory input from tonic firing neurons can either lead to bursting, tonic spiking or resting depending on the values of the channel conductances and the neuronal coupling level [START_REF] Butera | Models of respiratory rhythm generation in the pre-btzinger complex. i. bursting pacemaker neurons[END_REF][START_REF] Butera | Models of respiratory rhythm generation in the pre-btzinger complex. ii. populations of coupled pacemaker neurons[END_REF][START_REF] Del Negro | Models of respiratory rhythm generation in the pre-btzinger complex. iii. experimental tests of model predictions[END_REF]. Rhythm generation based on network bursting also depends on the bursting frequency and the interburst intervals. Synaptic properties shape the genesis and maintenance of bursts [START_REF] Staley | Presynaptic modulation of ca3 network activity[END_REF][START_REF] Verderio | Astrocytes are required for the oscillatory activity in cultured hippocampal neurons[END_REF][START_REF] Cohen | Homeostatic presynaptic suppression of neuronal network bursts[END_REF]. Synaptic short-term plasticity modeled in the mean-field approximation, are based on facilitation, depression and network firing rate [START_REF] Tsodyks | The neural code between neocortical pyramidal neurons depends on neurotransmitter release probability[END_REF], leading to a three-dimensional dynamical system. Long interburst intervals have been generated by introducing a double depression model [START_REF] Guerrier | Robust network oscillations during mammalian respiratory rhythm generation driven by synaptic dynamics[END_REF]. Interestingly, different levels of facilitation and depression lead to various network dynamics [START_REF] Barak | Persistent activity in neural networks with dynamic synapses[END_REF] such as resting, bursting or spiking and, when noise is added, up and down state transitions [START_REF] Holcman | The emergence of up and down states in cortical networks[END_REF]. Such models were used to interpret bursting in small hippocampal neuronal islands [START_REF] Dao Duc | Bursting reverberation as a multiscale neuronal network process driven by synaptic depression-facilitation[END_REF] to show that the correlation between successive bursts and interbursts could result from synchronous depressing-facilitating synapses. However in all these models, the interburst phase has not attracted much attention. The interburst in hippocampal pyramidal neurons is shaped by various type of potassium and calcium ionic channels [START_REF] Mckiernan | Ca1 pyramidal cells have diverse biophysical properties, affected by development, experience, and aging[END_REF][START_REF] David | Calcium-activated afterhyperpolarizations regulate synchronization and timing of epileptiform bursts in hippocampal ca3 pyramidal neurons[END_REF][START_REF] Tzingounis | Hippocalcin gates the calcium activation of the slow afterhyperpolarization in hippocampal pyramidal cells[END_REF][START_REF] Tzingounis | Contribution of kcnq2 and kcnq3 to the medium and slow afterhyperpolarization currents[END_REF], leading to medium and slow hyperpolarizing currents in the cells, a phenomenon known as afterhyperpolarization (AHP). We develop here a facilitation-depression model that accounts for AHP in order to better describe these interburst intervals. The manuscript is organized as follows: in the first part, we introduce a new dynamical system where we have added the AHP to the facilitation-depression model. We show that this model perturbed by a stochastic noise on the voltage variable can produce bursting periods followed by interburst intervals. We then study the phase-space that reveals three critical points (one attractor and two saddles). To further characterize the distribution of burst durations, we study the distribution of exit points on the stable manifold delimiting the basin of attraction of the stable equilibrium. Finally, we find an analytical formula for the burst and AHP durations. We also study the influence of network connectivity as well as facilitation and depression parameters on burst and interburst using numerical simulations.

A generalized facilitation-depression model accounting for AHP 1.Model description

We recall here the depression-facilitation short-term synaptic plasticity, a mean-field model for a sufficiently connected ensemble of neurons which consists of a stochastic dynamical system made of three equations [START_REF] Tsodyks | The neural code between neocortical pyramidal neurons depends on neurotransmitter release probability[END_REF][START_REF] Dao Duc | Bursting reverberation as a multiscale neuronal network process driven by synaptic depression-facilitation[END_REF] for the mean voltage h, the depression y, and the synaptic facilitation x:

τ ḣ = -h + Jxyh + + √ τ σ ω ẋ = X -x t f + K(1 -x)h + (1) ẏ = 1 -y t r -Lxyh + ,
The population average firing rate is given by h + = max(h, 0), which is a linear threshold function of the synaptic current [START_REF] Holcman | The emergence of up and down states in cortical networks[END_REF]. The term Jxy reflects the combined effect of synaptic short-term dynamics on the network activity. The second equation describes facilitation, while the third one describes depression.

The mean number of connections (synapses) per neurons is accounted for by the parameter J [START_REF] Bart | Modeling the spontaneous activity of the auditory cortex[END_REF]. We previously distinguished [START_REF] Dao Duc | Bursting reverberation as a multiscale neuronal network process driven by synaptic depression-facilitation[END_REF] the parameters K and L which describe how the firing rate is transformed into molecular events that are changing the duration and the probability of vesicular release respectively. The time scales t f and t r define the recovery of a synapse from the network activity. Finally, ω is an additive Gaussian white noise and σ its amplitude, that represent fluctuations in the firing rate.

The model (1) does not account for long AHP periods, that could be due to potassium channels [START_REF] Mckiernan | Ca1 pyramidal cells have diverse biophysical properties, affected by development, experience, and aging[END_REF], leading to a refractory period. To account for AHP, we modified the classical depression-facilitation model by introducing two features: 1) A new equilibrium state representing hyperpolarization 2) two timescales, for the medium and slow recovery to steady state.

To implement these novel properties, we decomposed the burst in four steps: step 1 starts with burst initiation and ends when the depression y starts increasing again, where we consider that hyperpolarization is initiated.

Step 2 lasts until y grows above the threshold Y h again. During this phase, we change the time constant τ of h to τ mAHP and the resting value of h from T to T AHP < T so that system (1) becomes:

τ 0 ḣ = -(h -T 0 ) + Jxy(h -T 0 ) + + √ τ 0 σ ω ẋ = X -x τ f + K(1 -x)(h -T 0 ) + ẏ = 1 -y τ r -Lxy(h -T 0 ) + , (2)
where τ 0 = τ mAHP and T 0 = T AHP for y < Y h and ẏ > 0. These changes forces the voltage to hyperpolarize. In step 3 the depression y is still increasing, that is: ẏ > 0, Y AHP < y and h > H AHP .

During this phase, we change the time constant τ mAHP to τ 0 = τ sAHP and the resting value of h is set to its initial value T 0 = T . These modifications accounts for the slow recovery from hyperpolarization to the resting state, this phase end when y reaches a second threshold Y AHP and h reaches another threshold H AHP .

Step 4 models the resting state, where y > Y AHP and h ≥ H AHP , τ 0 = τ and T 0 = T . All parameters are defined in Table 1.

We use [START_REF] Tsodyks | The neural code between neocortical pyramidal neurons depends on neurotransmitter release probability[END_REF][START_REF] Barak | Persistent activity in neural networks with dynamic synapses[END_REF][START_REF] Holcman | The emergence of up and down states in cortical networks[END_REF][START_REF] Dao Duc | Bursting reverberation as a multiscale neuronal network process driven by synaptic depression-facilitation[END_REF] to determine the values of the parameters for the classical facilitation-depression part and the order of magnitudes reviewed in [START_REF] Mckiernan | Ca1 pyramidal cells have diverse biophysical properties, affected by development, experience, and aging[END_REF] for the new AHP parameters (T AHP , τ mAHP and τ sAHP , Table 1). We will determine in sections 2.5.1 and 2.6 the effect on bursting dynamics of varying these parameters. Numerical simulations of equations ( 2) with a sufficient level of noise exhibit spontaneous bursts in the voltage variable followed by AHP periods (fig. 1A-B, upper). We segmented the simulated time series into two phases: burst (fig. 1C, blue) and interburst, which is further segmented into a AHP (pink) and quiescent phase (QP, green). We recall that we define the quiescent phase as the period where the voltage fluctuates around its equilibrium h = 0. This segmentation allows us to obtain the distributions of burst, AHP and QP durations (fig. 1D).

Phase-space analysis

We describe here the phase-space of the deterministic part of system (2).

Equilibrium points

Attractor. The first equilibrium point A is given by h = 0, x = X, y = 1 and the Jacobian at this point is

J A =         -1 + JX τ 0 0 K(1 -X) - 1 τ f 0 LX 0 - 1 τ r         . (3)
The eigenvalues

λ 1 = -1 + JX τ , λ 2 = - 1 τ f and λ 3 = - 1 τ r
are real strictly negative, confirming the nature of the attractor (fig. 1B and 2A, yellow star). Note that with the parameters of Table 1, the system shows three orders of magnitude as λ 1 = -12.6, λ 2 = -1.1, λ 3 = -0.34. This shows that the dynamics near the attractor is very anisotrope, very restricted to the plan perpendicular to the eigenvector associated to the highest eigenvalue |λ 1 |.

Saddle-points. Another solution of ḣ = 0 is given by Jxy = 1. Then ẋ = 0 yields showing a burst trajectory. The trajectory is decomposed into a QP (green), a burst (blue) and an AHP (pink) phase.

X -x τ f + K(1 -x)(h -T -T 0 ) = 0 ⇔ h = T + T 0 + x -X τ f K(1 -x) , re-injecting h in ẏ = 0 gives 1 - 1 Jx τ r - L J X -x τ f K(1 -x) = 0 ⇔ (Jτ f K + Lτ r )x 2 -(τ f K(J + 1) + LXτ r )x + τ f K = 0, the discriminant of this equation is ∆ = (τ f K(J + 1) + LXτ r ) 2 -4(Jτ f K + Lτ r )τ f K > 0. ( 4
The phase-space is divided into 3 regions according to the AHP dynamics: 1) the medium dynamics of hyperpolarization τ 0 = τ mAHP & T 0 = T AHP under and right of the orange surface where the trajectory is highlighted (orange circles).

2)

The slow recovery dynamics (τ 0 = τ sAHP & T 0 = 0, region under the purple plan) where the trajectory is highlighted (purple triangles).

3) The fast dynamics (τ 0 = τ & T 0 = 0). C. Segmentation of the voltage time series in burst (blue) and interburst (AHP (pink) and QP (green)). D. Distribution of burst (left, blue), AHP (center, pink) and QP (right, green) durations, extracted from numerical simulations lasting 10 4 s.

Thus

x 1,2 = τ f K(J + 1) + LXτ r ± √ ∆ 2(Jτ f K + Lτ r ) y 1,2 = 1 Jx 1,2 h 1,2 = T + T 0 + x 1,2 -X τ f K(1 -x 1,2 )
.

(5)

The Jacobians of the system at these points are

J S 1,2 =         0 Jy 1,2 (h 1,2 -T -T 0 ) + τ 0 Jx 1,2 (h 1,2 -T -T 0 ) + τ 0 K(1 -x 1,2 ) - 1 τ f -K(h 1,2 -T -T 0 ) + 0 - L J -Ly 1,2 (h 1,2 -T -T 0 ) + - 1 τ r -Lx 1,2 (h 1,2 -T -T 0 ) +         (6)
With the parameter values of Table 1, y 1,2 > Y h and thus T 0 = 0. Moreover, ẏ|y 1,2 < 0 so τ 0 = τ . We compute numerically the eigenvalues of the matrices J S 1,2 . The first saddle point S 1 has one real strictly negative eigenvalue and two complex-conjugate eigenvalues with positive real-parts, S 1 is a saddle-focus (with a repulsive focus and a stable manifold of dimension 1, fig. 2B). The second saddle point S 2 has two real negative eigenvalues and one positive one, it is a saddle-point with a stable manifold of dimension two and unstable of dimension one (fig. 2C).

Boundary of the basin of attraction associated to the stable equilibrium A

To determine the boundary of the basin of attraction for the point A, we ran numerical simulations of the deterministic system 2 (no noise), where we sample the entire (h,x,y)-space for the initial points and monitored where the trajectories escaped from the basin of attraction, characterized by a long trajectory, which describes the bursting phase. The limit values for the initial points, where trajectories escape define the separatrix surface Σ (fig. 2, cyan surface). Note that this separatrix Σ is constructed with a precision ∆h = 0.01 for a normalized amplitude of h to 1, which is smaller than the spatial scale of the stochastic component of the simulation σ √ τ ∆t ≈ 0.07.

To characterize the range of bursting durations, we determined numerically the durations of the shortest (red) and longest (purple) trajectories starting in the upper neighborhood of the separatrix Σ and ending below h = 0 (fig. 2D). The extreme trajectories are determined when we sampled the initial condition in the discretized approximation of Σ by a grid (x k , y q ) = (k∆x, q∆y) ∈ [0, 1] 2 , where we use for numerical computation ∆x = ∆y = 0.025.

To better understand how the stochastic system (2) bursts, we studied the distribution of exit points around the basin of attraction of A. We ran simulations with initial point A and a fixed level of noise, for each burst, and we recorded the intersection point of the trajectory and the separatrix (exit point). We show this distribution of points on the separatrix in fig. 3. In section 1.3 we shall explicit it analytically.

Distribution of exit points

We determine now the distribution of exit points located on the separatrix Σ when the initial point is at the attractor A = (0, X, 1). In this region of the phase-space, the dynamics simplifies to the system [START_REF] Wilson | Excitatory and inhibitory interactions in localized populations of model neurons[END_REF]. A. Repulsive trajectories from the saddle point S 1 (pink) and S 2 (blue) with corresponding eigenvectors (dashed arrows). B. Inset around S 1 . Real (dashed red arrow) and imaginary (dashed green arrow, see inset) parts of the eigenvectors associated to the complex conjugate repulsive eigenvalues and attractive eigenvector (dashed blue arrow). C. Inset around S 2 . Attractive eigenvectors (dashed blue and black arrows), and repulsive one (dashed red arrow). D. Longest (purple) and shortest (red) bursting trajectories starting outside the basin of attraction. without AHP:

τ ḣ = -h + Jxyh + + √ τ σ ω ẋ = X -x τ f + K(1 -x)h + ẏ = 1 -y τ r -Lxyh + , (7)
which can be written in the matrix form

ṡ = B(s) + √ σ Ẇ (8)
where s = (h, x, y) T and

B(s) =         b 1 (s) = - h τ + Jxyh + τ b 2 (s) = X -x τ f + K(1 -x)h + b 3 (s) = 1 -y τ r -Lxyh +         (9)
and √ σ = diag σ τ , 0, 0 . The probability density function q(s) of exit points is obtained by conditioning that the trajectories of the process ( 8) are absorbed on Σ. It is solution of the Fokker-Planck renewal equation (FPE) [START_REF] Schuss | Theory and Applications of Stochastic Processes: An Analytical Approach[END_REF][START_REF] Schuss | Nonlinear filtering and optimal phase tracking[END_REF] -∂ ∂h

(Jxy -1)h τ q(s) - ∂ ∂x X -x τ f + K(1 -x)h q(s) - ∂ ∂y 1 -y τ r -Lxyh q(s) + σ 2τ ∂ 2 ∂h 2 q(s) = δ(s -s 0 ) (10) q(s ∈ Σ|s 0 ) = 0. ( 11
)
We use WKB approximation to search for a solution of equation [START_REF] Zeldenrust | Neural coding with burstscurrent state and future perspectives[END_REF] in the form

q(s|s 0 ) = Q σ (s)e - ψ(s) σ , (12)
where Q σ is a regular function with the formal expansion

Q σ (s) = ∞ i=0 Q i (s)σ i . ( 13
)
The function ψ satisfies the eikonal equation

(Jxy -1)h τ ∂ψ ∂h + X -x τ f + K(1 -x)h ∂ψ ∂x + 1 -y τ r -Lxyh ∂ψ ∂y + 1 2τ ∂ψ ∂h 2 = 0 (14)
We use the method of characteristics to solve the eikonal equation. Setting

p = ∇ψ =   p 1 p 2 p 3   , (15)
and using the classical notation

F (s, ψ, p) = b 1 (s)p 1 + b 2 (s)p 2 + b 3 (s)p 3 + 1 2τ p 2 1 , (16)
the characteristics are given by dh dt

= F p 1 = b 1 + 1 τ p 1 dx dt = F p 2 = b 2 dy dt = F p 3 = b 3 , (17)
dp 1 dt = -F h = - Jxy -1 τ p 1 -K(1 -x)p 2 + Lxyp 3 dp 2 dt = -F x = - Jyh τ p 1 + 1 τ f + Kh p 2 + Lyhp 3 dp 3 dt = -F y = - Jxh τ p 1 + 1 τ r + Lxh p 3 (18)
and

dψ dt = 1 2τ p 2 1 . (19)
We should solve ( 17)-( 19) starting at the attractor s 0 = A, however, this characteristic will be trapped at A. To avoid this difficulty, we follow the method proposed in [START_REF] Schuss | Nonlinear filtering and optimal phase tracking[END_REF] p.165-170, and we start from points located in a neighborhood V A of A. In V A , the solution of the eikonal equation has a quadratic approximation

ψ(s) = 1 2 s T Rs + o(|s| 2 ). ( 20
)
To find the matrix R, we use the linearized eikonal equation around the attractor A

(J A s) T • ∇ψ + 1 2τ p 2 1 = 0, (21)
where J A is the Jacobian defined in (3). This matrix equation does not have a unique solution, but we shall use the one given by ψ(s)

≈ (1 -JX)h 2 . ( 22
)
Choosing initial conditions on the contours ψ(s) = δ = 0.05, that is

h = ± δ 1 -JX ≈ 0.28, (23)
we computed the characteristics numerically (fig. 4A-B). To determine the exit points distribution, we now solve the transport equation ( 24)

1 -Jxy τ h ∂Q 0 ∂h + Q 0 + 1 τ f + Kh Q 0 - X -x τ f + K(1 -x)h ∂Q 0 ∂x + 1 τ r + Lxh Q 0 - 1 -y τ r -Lxyh ∂Q 0 ∂y - 1 τ ∂Q 0 ∂h ∂ψ ∂h - Q 0 2τ ∂ 2 ψ ∂h 2 = 0. ( 24
)
To find Q 0 , we follow the method from [START_REF] Schuss | Nonlinear filtering and optimal phase tracking[END_REF] p.172-175. We rewrite equation ( 24)

B • ∇Q 0 + 1 τ ∂Q 0 ∂h ∂ψ ∂h = -∇ • B + 1 2τ ∂ 2 ψ ∂h 2 Q 0 ( 25
)
where B is defined in [START_REF] Coombes | Bursting: The Genesis Of Rhythm In The Nervous System[END_REF]. Along the characteristics, ( 25) is

dQ 0 (s(t)) dt = ∇Q 0 (s(t)) • ds(t) dt = -∇ • B(s(t)) + 1 2τ ∂ 2 ψ(s(t)) ∂h 2 Q 0 (s(t)). ( 26
)
Our goal is to compute Q 0 on the separatrix and for that purpose, we need to evaluate ∂ 2 ψ(s(t)) ∂h 2 by differentiating the characteristics equations ( 17)-( 19) with respect to the initial point s 0 = s(0). Setting

s j (t) = ∂s(t) ∂s j 0 , p j (t) = ∂p(t) ∂s j 0 , ∂ 2 ψ(s(t)) ∂s i ∂s j = R i,j (t), (27)
we have R(t) = P (t)S(t) -1 , where P (t) (resp. S(t)) is the matrix with columns p j (t) (resp. s j (t)). The initial conditions are

s i j (0) = δ i,j , p i j (0) = ∂ 2 ψ(0) ∂s i ∂s j = R i,j . (28)
The dynamic has the form

ds 1 1 (t) dt = dh 1 dt = ∂b 1 ∂h + 1 τ ∂p 1 ∂h h 1 ds 1 2 (t) dt = dh 2 dt = ∂b 1 ∂h + 1 τ ∂p 1 ∂h h 2 ds 1 3 (t) dt = dh 3 dt = ∂b 1 ∂h + 1 τ ∂p 1 ∂h h 3 ds 2 1 (t) dt = dx 1 dt = ∂b 2 ∂x x 1 ds 2 2 (t) dt = dx 2 dt = ∂b 2 ∂x x 2 ds 2 3 (t) dt = dx 3 dt = ∂b 2 ∂x x 3 ds 3 1 (t) dt = dy 1 dt = ∂b 3 ∂y y 1 ds 3 2 (t) dt = dy 2 dt = ∂b 3 ∂y y 2 ds 3 3 (t) dt = dy 3 dt = ∂b 3 ∂y y 3 (29)
and because we are only interested in R 1,1 we only need to compute the first row of P (t), thus

dp 1 1 (t) dt = - Jxy -1 τ ∂p 1 ∂h -K(1 -x) ∂p 2 ∂h + Lxy ∂p 3 ∂h h 1 dp 1 2 (t) dt = - Jxy -1 τ ∂p 1 ∂h -K(1 -x) ∂p 2 ∂h + Lxy ∂p 3 ∂h h 2 dp 1 3 (t) dt = - Jxy -1 τ ∂p 1 ∂h -K(1 -x) ∂p 2 ∂h + Lxy ∂p 3 ∂h h 3 . (30)
In the limit t → ∞ the characteristic that hits the saddle point S 2 is tangent to the separatrix and

-∇ • B + 1 2τ ∂ 2 ψ ∂h 2 Q 0 → -∇ • B |S 2 ≈ 1.82. Indeed,
∂ 2 ψ ∂h 2 tends to 0 near the saddle point S 2 as shown in fig. 4C. Thus, near the saddle point, we have

dQ 0 (s(t)) dt = -(∇ • B |S 2 + o(1))Q 0 (s(t)). ( 31
)
The solution is approximated by

Q 0 (s(t)) = Q 0 (s(0))e -∇ • B |S 2 t(1 + o(1)) . (32)
Finally, the characteristic s(t) near the saddle point S 2 can be expressed with respect to the arc length s:

s(t) ≈ t 0 ṡ2 (u) 2 du, (33)
where s 2 is the dominant coordinate of s ∈ Σ in the eigenvectors basis of the jacobian J S 2 of system (2) at S 2 , (λ 1 ≈ -4.58 and λ 2 ≈ -0.25), thus locally

s(t) ≈ t 0 s 2 (0)e 2λ 2 u du, (34)
and

s(t) ≈ t 0 s 2 (0)e 2λ 2 u du = s 2 (0) e λ 2 t -1 λ 2 . ( 35
)
Finally, using [START_REF] David | Calcium-activated afterhyperpolarizations regulate synchronization and timing of epileptiform bursts in hippocampal ca3 pyramidal neurons[END_REF] and [START_REF] Bart | Modeling the spontaneous activity of the auditory cortex[END_REF], we obtain locally

Q 0 (s) = Q 0 (0)s - ∇ • B |S 2 λ 2 , (36)
where

- ∇ • B |S 2 λ 2 ≈ -7.23.
To connect the solution q of the FPE [START_REF] Zeldenrust | Neural coding with burstscurrent state and future perspectives[END_REF] to the distribution of exit points, we have to account for the boundary layer function q σ that has to be added to the transport solution in the form Q 0 q σ . This product satisfies the boundary condition [START_REF] Brunel | Fast global oscillations in networks of integrate-and-fire neurons with low firing rates[END_REF]. We do not compute here q σ as the computation follows the one of [START_REF] Schuss | Nonlinear filtering and optimal phase tracking[END_REF] p. 182-183 near the separatrix. It is a regular function of the form -

2 π ργ(s 1 ,s 2 )/ √ σ 0 e -η 2 /2 dη,
where ρ is the distance to the separatrix Σ in a neighborhood of S 2 and γ(s 1 , s 2 ) a regular function.

Finally, we recall that the exit point distribution per unit surface ds is given by

p Σ (s|s 0 ) = 1 N J(s|s 0 ) • ν(s)ds for s ∈ Σ (37)
where the probability flux is

J(s|s 0 ) =           Jxy -1 τ hq(s) - σ 2τ ∂q(s) ∂h X -x τ f + K(1 -x)h q(s) 1 -y τ r -Lxyh q(s)           , (38)
and the normalization constant is

N = Σ J(s|s 0 ) • ν(s)ds, (39)
where ν(s) is the unit normal vector at the point s. The flux is computed by differentiating expression [START_REF] Neltner | On synchrony of weakly coupled neurons at low firing rate[END_REF],

q(s|s 0 ) = q σ (s)Q 0 (s)e - ψ(s) σ . (40)
We obtain

J(s|s 0 ) • ν(s)ds = -N 2σ π q(s|s 0 )γ(s 1 , s 2 )ds = K 0 s- ∇ • B |S 2 λ 2 e - ψ(s) σ ds, (41)
where γ(s 1 , s 2 ) has been approximated by its value at s = 0. Furthermore, in the limit s → 0, s-

∇ • B |S 2 λ 2
tends to infinity, however it is compensated by e -ψ(s) σ which is small enough, as we observe numerically. We plotted the distribution of exit points in fig. 4D-E for K 0 = 1. Finally, we compare the distribution p Σ with the one obtained from the stochastic simulations of system (2) with the same level of noise (σ = 3). Both distributions are peaked, showing that the exit points are constrained in a small area of the separatrix. To conclude this part, our two different numerical methods confirm that the exit point distribution is peaked, thus the trajectories associated to the bursting periods are confined in a tubular neighborhood of a generic trajectory and thus the distribution of the bursting times is peaked, as observed in fig. 1D. 2 Computing analytically burst and AHP durations

Approximated equations

In this section we develop an approximation procedure to compute the mean bursting and hyperpolarization durations from the AHP facilitation-depression model [START_REF] Wilson | Excitatory and inhibitory interactions in localized populations of model neurons[END_REF]. The approximation procedure is based on the following considerations: because in the first phases of burst and AHP, the voltage h evolves much faster than the facilitation x and depression y, to compute the duration of the bursting phase, we will replace the dynamics of h in the depression and facilitation equations by a piecewise constant function H(t) (fig. 5). This approximation decouples the system (2), thus x and y can be computed.

We shall now specify the function H(t). In the bursting phase, it is constant equal to

H 1 for t ∈ [0; t 1 ]
where t 1 will be specified in section 2.4. In the hyperpolarization phase, H(t) = H 2 for t ∈ [0; t 2 ]. For t > t 2 (that will also be specified in section 2.4), we choose H(t) = 0 to account for the recovery phase.

H(t) =    H 1 , for t ∈ [0, t 1 ] H 2 , for t ∈]t 1 , t 2 ] 0 for t > t 2 . (42)
The approximated system of equations becomes:

τ 0 ḣ = -(h -T 0 (t)) + Jxy(h -T 0 (t)) + ẋ = X -x τ f + K(1 -x)H(t) ẏ = 1 -y τ r -LxyH(t) (43)
where the AHP is accounted for by changing the threshold and timescales as follows

T 0 (t) =    0 for t ∈ [0, t 1 ] T AHP for t ∈]t 1 , t 2 ] 0 for t > t 2 and τ 0 (t) =    τ for t ∈ [0, t 1 ] τ mAHP for t ∈]t 1 , t 2 ] τ sAHP for t > t 2 .
(44)

Computing the facilitation and depression dynamics in three phases

Phase 1

We integrate the facilitation and depression equations in (43). During the bursting phase (fig. 5, phase 1, blue) H(t) = H 1 . We use the following initial conditions: x(0) = X and y(0) = 1 (resting values). We obtain

x(t) = A 1 e -α 1 t + B 1 , (45)
where

α 1 = 1 τ f + KH 1 , A 1 = KH 1 (X -1) α 1 , B 1 = X τ f + KH 1 α 1 . (46)
Injecting expression (45) in the third equation of system (43), we obtain

ẏ = 1 -y τ r -L(A 1 e -α 1 t + B 1 )H 1 y t 1 t 2 H 2 H 1
1. Burst: fast timescale τ 0 =τ, no hyperpolarization: T 0 =0

2. AHP initiation: medium timescale τ 0 =τ mAHP , hyperpolarization: T 0 = T AHP <0

3. Recovery: slow timescale τ 0 =τ sAHP , no hyperpolarization: T 0 =0 The solution is

y(t) = C 1 + 1 τ r t 0 exp(f 1 (s))ds exp(-f 1 (t)), (47)
where the function

f 1 (t) = β 1 t - LA 1 H 1 α 1 e -α 1 t .
To approximate the integral t 0 exp(f 1 (s))ds, we use that f 1 is monotonic on the interval [0; t 1 ], thus using a Taylor expansion at order 1, we get

t 0 exp(f 1 (s))ds ≈ exp(f 1 (t)) t 0 exp(f 1 (t)(s -t))ds = exp(f 1 (t)) f 1 (t) (1 -exp(-tf 1 (t))). ( 48
)
Using expression (47), we obtain the approximation for t ∈ [0, t 1 ]

y(t) ≈ 1 τ r (1 -exp(-tf 1 (t))) f 1 (t) + C 1 exp(-f 1 (t)), (49)
where

β 1 = 1 τ r + B 1 LH 1 and C 1 = exp - LH 1 A 1 α 1 . ( 50
)
Phase 2

The second phase starts at t 1 where H(t) = H 2 , where the equations and the approximation are similar to the paragraph above. However we use the following initial conditions: x(t - 1 ) = x(t + 1 ) and y(t - 1 ) = y(t + 1 ). This yields for t ∈ [t 1 ; t 2 ],

x(t

) = A 2 e -α 2 t + B 2 , (51)
where

α 2 = 1 τ f + KH 2 , A 2 = (x(t - 1 ) -B 2 )e α 2 t 1 , B 2 = X τ f + KH 2 α 2 . ( 52
)
y(t) ≈ 1 τ r (1 -exp(-(t -t 1 )f 2 (t))) f 2 (t) + C 2 exp(-f 2 (t)),
where

f 2 (t) = β 2 t - LA 2 H 2 α 2 e -α 2 t , (53)
β 2 = 1 τ r + B 2 LH 2 and C 2 = y(t - 1 ) exp(f 2 (t 1 )). ( 54
)
Phase 3

The recovery phase starts at t 2 where H(t) = 0. We use the following initial conditions: x(t - 2 ) = x(t + 2 ) and y(t - 2 ) = y(t + 2 ). This yields for t ≥ t 2 to the solution

x(t) = X + (x(t - 2 ) -X) exp - t -t 2 τ f (55)
y(t) = 1 + (y(t - 2 ) -1) exp - t -t 2 τ r .

Computing the approximated voltage in the three phases

To compute the voltage, we will use the approximations for x and y described for the three phases in paragraph 2.2. Although the previous approximations might be drastic for x and y, we shall see that they provide a very good approximation for h. In addition, they allowed to decouple the system of equations and thus h can now be computed explicitly.

Phase 1

The first equation in system (43) is

τ ḣ = -h + Jxyh + (56)
where the initial condition is h(0) = H1 . A direct integration leads to

h(t) = H1 exp - t τ + J τ t 0 x(s)y(s)ds . (57)
To obtain an explicit dependency of the solution h with respect to the parameters, we will use expressions (45) and (49) for x and y respectively to compute the integral in expression (57). This calculation is detailed in appendix A. We note that H1 could be different from H 1 , indeed to guarantee that the facilitation and depression, that have slower dynamics compared to the voltage, are immediately in the bursting state we choose H 1 H1 (see Table 3).

Phase 2

In phase 2, we use equation (44) for T 0 = T AHP and τ 0 = τ mAHP so that

τ mAHP ḣ = -(h -T AHP ) + Jxy(h -T AHP ) + . ( 58
)
We use the initial condition h(t + 1 ) = h(t - 1 ), and obtain by a direct integration

h(t) = (h(t - 1 ) -T AHP ) exp - t -t 1 τ mAHP + J τ mAHP t t 1
x(s)y(s)ds + T AHP .

(59)

Similar to phase 1, we detail this calculation in appendix A.

Phase 3

Finally, when t > t 2 , h enters into its slow relaxation phase leading in equation ( 44) to T 0 = 0 and τ 0 = τ sAHP , and the initial condition h(t - 2 ) = h(t + 2 ). A direct integration of equation ( 43) leads to

h(t) = h(t - 2 ) exp - t -t 2 τ sAHP + J τ sAHP t t 2
x(s)y(s)ds .

Depression y

Figure 6: Analytical approximation (green) vs exact solution (dashed magenta) for h, x and y.

Identification of the termination times t 1 and t 2

End of phase 1

Following burst activation, medium and slow K + channels start to be activated forcing the voltage to hyperpolarize. To account for the overall changes in the voltage dynamics due to this K + channels activation, we change the recovery timescale τ 0 to τ mAHP (equation ( 44)) and H(t) to H 2 in (42) at time t 1 . In practice the hyperpolarization initiation is defined in the region where h is decreasing after reaching its maximum, as the first time t 1 where h(t 1 ) = h 0 (expression (57)), leading to equation

t 1 B 1 J -τ r β 1 τ r β 1 J - 1 τ r LH 1 α 1 ln      1 + LA 1 H 1 β 1 e -α 1 t 1 1 + LA 1 H 1 β 1      + e -α 1 t 1 -1 α 1 LA 1 H 1 B 1 τ r β 2 1 - e -2α 1 t 1 -1 2α 1 (LA 1 H 1 ) 2 B 1 τ r β 3 1 + e -(β 1 +LA 1 H 1 )t 1 -1 β 1 + LA 1 H 1 B 1 (1 -τ r β 1 ) τ r β 1 + e -(α 1 +β 1 +LA 1 H 1 )t 1 -1 α 1 + β 1 + LA 1 H 1 A 1 (-τ r β 2 1 + β 1 -LH 1 B 1 ) τ r β 2 1 + e -(2α 1 +β 1 +LA 1 H 1 )t 1 -1 2α 1 + β 1 + LA 1 H 1 - LA 2 1 H 1 τ r β 2 1 + (LA 1 H 1 ) 2 B 1 τ r β 3 1 + e -(3α 1 +β 1 +LA 1 H 1 )t 1 -1 3α 1 + β 1 + LA 1 H 1 L 2 A 3 1 H 2 1 τ 2 r β 3 1 = τ J ln h 0 H 1 . ( 61
)
This transcendental equation cannot be solved explicitly however, with the parameters from Table 1 and Table 3, and the order of t 1 , we can neglect the exponential terms in (61), leading to

t 1 = 1 Γ 1 -Γ 2 ln 1 1 + Γ 3 + Γ 4 α 1 + Γ 5 β 1 + LA 1 H 1 + Γ 6 α 1 + β 1 + LA 1 H 1 + Γ 7 2α 1 + β 1 + LA 1 H 1 + Γ 8 3α 1 + β 1 + LA 1 H 1 + Γ 9 2α 1 + τ J ln h 0 H1 , (62)
where

Γ 1 = B 1 J -τ r β 1 τ r β 1 J , Γ 2 = - 1 τ r LH 1 α 1 , Γ 3 = LA 1 H 1 β 1 , Γ 4 = LA 1 H 1 B 1 τ r β 2 1 , Γ 5 = B 1 (1 -τ r β 1 e LA 1 h 1 α 1 ) τ r β 1 , Γ 6 = A 1 (-τ r β 2 1 e LA 1 h 1 α 1 + β 1 -LH 1 B 1 ) τ r β 2 1 , Γ 7 = LA 2 1 H 1 τ r β 2 1 LH 1 B 1 β 1 -1 , Γ 8 = L 2 A 3 1 H 2 1 τ 2 r β 3 1 and Γ 9 = - (LA 1 H 1 ) 2 B 1 τ r β 3 1 .
(63) Using the parameter values from Table 1 and Table 3, we obtain t 1 ≈ 200 ms. This time suggests that the medium and slow K + channels start to be activated quite early following burst initiation.

End of phase 2

The second phase, dominated by the hyperpolarization, ends when the voltage reaches asymptotically its minimum. In practice we introduce a threshold h AHP so that when h(t 2 ) = h AHP (expression (59)), we switch into the third phase (see [START_REF] Schuss | Where is the exit point?[END_REF] and ( 44)). This leads to equation

(t 2 -t 1 ) B 2 J -τ r β 2 τ r β 2 J - 1 τ r LH 2 α 2 ln      1 + LA 2 H 2 β 2 e -α 2 t 1 e -α 2 (t 2 -t 1 ) 1 + LA 2 H 2 β 2 e -α 2 t 1      + e -α 2 (t 2 -t 1 ) -1 α 2 e -α 2 t 1 LA 2 H 2 B 2 τ r β 2 2 - e -2α 2 (t 2 -t 1 ) -1 2α 2 e -2α 2 t 1 (LA 2 H 2 ) 2 B 2 τ r β 3 2 + e -(β 2 +LA 2 H 2 )(t 2 -t 1 ) -1 β 2 + LA 2 H 2 e -(β 2 +LA 2 H 2 )t 1 B 2 (1 -C 2 e LA 2 H 2 α 2 τ r β 2 ) τ r β 2 e -(α 2 +β 2 +LA 2 H 2 )(t 2 -t 1 ) -1 α 2 + β 2 + LA 2 H 2 e -(α 2 +β 2 +LA 2 H 2 )t 1 A 2 (-C 2 e LA 2 H 2 α 2 τ r β 2 2 + β 2 -LH 2 B 2 ) τ r β 2 2 + e -(2α 2 +β 2 +LA 2 H 2 )(t 2 -t 1 ) -1 2α 2 + β 2 + LA 2 H 2 e -(2α 2 +β 2 +LA 2 H 2 )t 1 - LA 2 2 H 2 τ r β 2 2 + (LA 2 H 2 ) 2 B 2 τ r β 3 2 + e -(3α 2 +β 2 +LA 2 H 2 )(t 2 -t 1 ) -1 3α 2 + β 2 + LA 2 H 2 e -(3α 2 +β 2 +LA 2 H 2 )t 1 L 2 A 3 2 H 2 2 τ 2 r β 3 2 = τ mAHP J ln h AHP -T AHP h(t - 1 ) -T AHP .
Here all terms are of the same order thus we cannot neglect any of them. Since we just need to estimate the value of t 2 to calibrate our approximated model we solve numerically the following transcendental equation

Λ 1 (t 2 -t 1 ) + Λ 2 ln 1 + Λ 3 e -α 2 (t 2 -t 1 ) 1 + Λ 3 + Λ 4 e -α 2 (t 2 -t 1 ) -1 α 2 + Λ 5 e -(β 2 + LA 2 H 2 )(t 2 -t 1 ) -1 β 2 + LA 2 H 2 +Λ 6 e -(α 2 + β 2 + LA 2 H 2 )(t 2 -t 1 ) -1 α 2 + β 2 + LA 2 H 2 + Λ 7 e -(2α 2 + β 2 + LA 2 H 2 )(t 2 -t 1 ) -1 2α 2 + β 2 + LA 2 H 2 +Λ 8 e -(3α 2 + β 2 + LA 2 H 2 )(t 2 -t 1 ) -1 3α 2 + β 2 + LA 2 H 2 + Λ 9 e -2α 2 (t 2 -t 1 ) -1 2α 2 - τ mAHP J ln h AHP -T AHP h(t - 1 ) -T AHP = 0, (64)
where

Λ 1 = B 2 J -τ r β 2 τ r β 2 J , Λ 2 = - 1 τ r LH 2 α 2 , Λ 3 = LA 2 H2 β 2 e -α 2 t 1 , Λ 4 = LA 2 H 2 B 2 τ r β 2 2 e -α 2 t 1 , Λ 5 = B 2 (1 -τ r β 2 C 2 e LA 2 H 2 α 2 ) τ r β 2 e -(β 2 + LA 2 H2)t 1 , Λ 6 = A 2 (-C 2 e LA 2 H 2 α 2 τ r β 2 2 + β 2 -LH 2 B 2 ) τ r β 2 2 e -(α 2 + β 2 + LA 2 H 2 )t 1 , Λ 7 = LA 2 2 H 2 τ r β 2 2 LH 2 B 2 β 2 -1 e -(2α 2 + β 2 + LA 2 H 2 )t 1 , Λ 8 = L 2 A 3 2 H 2 2 τ 2 r β 3 2 e -(3α 2 + β 2 + LA 2 H 2 )t 1 and Λ 9 = - (LA 2 H 2 ) 2 B 2 τ r β 3 2 e -2α 2 t 1 . (65)
Using parameter values defined in Table 1 and Table 3, and the value of t 1 computed in the previous section we obtain t 2 ≈ 1.4 s. The obtained analytical approximation is plotted in fig. 6 (green) in comparison to the exact solution obtained using numerical simulations (dashed magenta).

Bursting and AHP durations Bursting duration

The burst duration is defined from the voltage jump at time t = 0 to h(t) = H 1 and ends when h(t i ) = 0 for the first time. In practice, we use expression (59) as in section 2.4 for the end of phase 2 however, here t i -t 1 is small enough to allow us to use Taylor expansions to second order leading to the quadratic equation

Λ(t i -t 1 ) 2 + Λ(t i -t 1 ) - τ mAHP J ln -T AHP h(t - 1 ) -T AHP = 0, (66)
where

Λ = Λ 2 Λ 3 α 2 2 2(1 + Λ 3 ) 2 + 1 2 α 2 Λ 4 + (β 2 + LA 2 H 2 )Λ 5 + (α 2 + β 2 + LA 2 H 2 )Λ 6 + (2α 2 + β 2 + LA 2 H 2 )Λ 7 +(3α 2 + β 2 + LA 2 H 2 )Λ 8 + 2α 2 Λ 9 ,
and

Λ = - Λ 2 Λ 3 α 2 1 + Λ 3 + Λ 1 -Λ 4 -Λ 5 -Λ 6 -Λ 7 -Λ 8 -Λ 9
We keep the positive root

t i = t 1 + -Λ -Λ 2 + 4 Λ τ mAHP J ln -T AHP h(t - 1 ) -T AHP 2 Λ . (67)
Using parameters from Table 1 and Table 3, we obtain t i ≈ 940 ms, which is comparable to the bursting times observed in experimental data [START_REF] Chever | Astroglial networks promote neuronal coordination[END_REF], and from our numerical simulations (fig. 1D).

AHP duration

The AHP starts at time t i computed above, however using expression (60) the termination time would be infinite. Thus, we introduce a threshold and define the end of AHP t e such as h(t e ) = . In practice, the value can be estimated from the amplitude of the voltage fluctuations at equilibrium. We obtain from expression (60)

- 1 J + X (t e -t 2 ) -τ r X(y(t - 2 ) -1)   e - t e -t 2 τ r -1   -τ f (x(t - 2 ) -X)    e - t e -t 2 τ f -1    - (y(t - 2 ) -1)(x(t - 2 ) -X)τ f τ r τ f + τ r    e -(t -t 2 ) τ f + τ r τ f τ r -1    = τ sAHP J ln h(t - 2 )
because t e -t 2 is large enough, we neglect the exponential terms so that

(t e -t 2 ) X - 1 J + τ r X(y(t - 2 ) -1) + τ f (x(t - 2 ) -X) + (x(t - 2 ) -X)(y(t - 2 ) -1)τ f τ r τ f + τ r = τ sAHP J ln h(t - 2 )
, leading to

t e = t 2 + τ sAHP J ln h 2 -τ r X(y(t - 2 ) -1) -τ f (x(t - 2 ) -X) - (x(t - 2 ) -X)(y(t - 2 ) -1)τ f τ r τ f + τ r J JX -1 ,
Using the parameter values from Table 1 and Table 3 we obtain t e ≈ 15.4 s and ∆ AHP ≈ 14.4 s, which is coherent with the durations obtained from the numerical simulations (fig. 1D), as well as classical AHP durations found in the literature [START_REF] Mckiernan | Ca1 pyramidal cells have diverse biophysical properties, affected by development, experience, and aging[END_REF].

Study of parameter influence on burst and AHP durations

To evaluate the influence of the main parameters on the bursting and AHP durations we plotted these times vs the recovery timescales τ mAHP and τ sAHP , the hyperpolarization level T AHP and the arbitrary thresholds h 0 , H1 , h AHP and . First, the burst duration that varies between 0.5 and 3s, is an increasing function of τ mAHP and does not depend much on T AHP in the range [-15; -40] (fig. 7A). In addition, the AHP duration increases with τ sAHP , but in a larger range from 9 to 35s. However, the hyperpolarization level T AHP has a larger influence on this duration (fig. 7B). To verify that the arbitrary thresholds that we use do not influence much the burst and AHP durations, we plotted them in fig. 7C-F with respect to the phase 1 termination threshold h 0 , the phase 2 termination threshold h AHP , the duration of phase 1 t 1 and the AHP termination threshold respectively. These figures show that there is almost no dependency with respect to H1 and T AHP , as well as h 0 and h AHP due to the effect of the logarithmic term.

Numerical analysis of burst and interburst durations: effect of J, K, L parameters

To study the influence of the network connectivity J on burst, AHP and QP durations, we ran numerical simulations of the stochastic system (2), where we varied J, as well as the facilitation and depression parameters K and L. To determine the time distributions of burst and interburst, we segmented the traces obtained for 5000 seconds simulations with a noise amplitude σ = 6 and computed the mean value of the bursts (fig. 8A), AHP (fig. 8B) and QP durations (fig. 8C). Interestingly, we observe two different regimes depending on the values of the parameters: no bursts (J < 3.05 for K = 0.047, L = 0.028; J < 3.2 for K = 0.037, L = 0.028; J < 3.5 for K = 0.027, L = 0.028; fig. 8 left column, or J < 3.7 for K = 0.037, L = 0.038 and J < 4.1 for K = 0.037, L = 0.048, right column) and bursts followed by AHP (for higher values of J). Surprisingly, in the bursting regime changing J does not influence the mean burst duration. However, AHP durations decreases as J increases. Finally, QP durations reach a peak at the transition value of J between the two regimes and then quickly decrease around QP ≈ 25s. We note that the mean burst durations obtained here are longer than the ones observed in fig. 1D, this is due to the fact that in these simulations, we used σ = 6 (vs σ = 3 for fig. 1D). Indeed, increasing the noise increases the mean burst duration because, at the beginning of the burst, the deterministic part of the trajectory is still perturbed by the noise component, leading to a longer trajectory when the noise level is higher.

To conclude, a sufficient connectivity level is necessary to generate bursting, however once the dynamics enter into this regime, increasing the level of neuronal connectivity does not change much the bursting times.

Phase 1 termination threshold h 0 Duration of phase 1 t

H 1 =350 H 1 =360 H 1 =370 H 1 =380 H 1 =390 H 1 =400

Conclusion and discussion

We present here a novel mean-field model of synaptic short-term plasticity for the voltage, depression and facilitation variables that now accounts for long AHP periods. This model generalizes the depressionfacilitation model introduced in [START_REF] Tsodyks | The neural code between neocortical pyramidal neurons depends on neurotransmitter release probability[END_REF] and developed in [START_REF] Holcman | The emergence of up and down states in cortical networks[END_REF][START_REF] Dao Duc | Bursting reverberation as a multiscale neuronal network process driven by synaptic depression-facilitation[END_REF][START_REF] Mongillo | Synaptic theory of working memory[END_REF][START_REF] Dao Duc | Synaptic dynamics and neuronal network connectivity are reflected in the distribution of times in up states[END_REF]. The AHP significantly increases the interburst duration by introducing a recovery phase after network bursting. When a Gaussian noise of small amplitude is added to the dynamics, it exhibits spontaneous bursts followed by AHP periods. We have studied here the distribution of bursts and of interbursts, decomposed in AHP and QP durations. Interestingly, we found that the distribution of bursts durations is quite concentrated (subsection 1.1). To explain this property, we studied the three-dimensional phase-space of the dynamical system (2), that contains one attractor and two saddle points. By computing numerically the two-dimensional stable manifold at one of the saddles, we found the distribution of exit points (on this manifold) when the initial point of the stochastic dynamics is located at the attractor. To compute this distribution we used two methods: 1) stochastic simulations, and 2) the method of characteristics to solve the FPE (10) in the limit of small noise. In both cases, we found a peaked distribution of exit points close to the saddle point, as predicted for two-dimensional stochastic systems [START_REF] Schuss | Nonlinear filtering and optimal phase tracking[END_REF][START_REF] Schuss | The exit distribution on the stochastic separatrix in kramers' exit problem[END_REF][START_REF] Schuss | Where is the exit point?[END_REF], summarized by expression [START_REF] Schuss | The exit distribution on the stochastic separatrix in kramers' exit problem[END_REF]. After the stochastic trajectories have crossed the separatrix, they follow an almost deterministic behavior, confirming that the distribution of exit points on the separatrix defines the spread of the distribution of burst durations. We also derived here analytical formulas (subsection 2.5) that reveal the influence of the parameters on burst and AHP durations. These computations can be used to calibrate the AHP parameters with respect to the expected values of burst and AHP durations, that could be measured experimentally. This model could thus be used to decipher the main mechanisms leading to changes in bursting and interburst dynamics, for example when the neuronal network is disrupted, during epilepsy or in the case of a glial network alteration [START_REF] Chever | Astroglial networks promote neuronal coordination[END_REF]. Classical bursting models describe accurately the burst phase [START_REF] Izhikevich | Dynamical Systems in Neuroscience: the Geometry of Excitability and Bursting[END_REF][START_REF] Ermentrout | Parabolic bursting in an excitable system coupled with a slow oscillation[END_REF][START_REF] Coombes | Bursting: The Genesis Of Rhythm In The Nervous System[END_REF][START_REF] Chizhov | Conductance-based refractory density approach for a population of bursting neurons[END_REF], but interburst is often considered as the continuation in the phase-space of the deterministic trajectories. Here the interburst phase is composed of a deterministic refractory period, the AHP, followed by the escape from an attractor due to noise (subsection 1.1). During successive bursts, trajectories are not reset at the attractor, but explore the basin of attraction. This exploration depends on the previous bursting trajectory. Thus, we expect a correlation between successive burst and interburst durations. This correlation may also depend on the amplitude of the voltage fluctuations. Finally, we predict that modifying the AHP duration could affect bursting, because it corresponds to a change in the attractor's position and dominates the effect of synaptic depression.

A Calculation details of burst and AHP durations A.1 Integral term of h in phase 1

To compute the integral in expression (57), we split it into two parts:

t 0 x(s)y(s)ds = C 1 t 0 (A 1 e -α 1 s + B 1 )e -f 1 (s) I ds + t 0 (A 1 e -α 1 s + B 1 )
1 -e -sf 1 (s) τ r f 1 (s) ds II .

We start by I:

I = C 1 A 1 t 0 e -(α 1 + β 1 )s + LA 1 H 1 α 1 e -α 1 s ds I A + C 1 B 1 t 0 e -β 1 s + LA 1 H 1 α 1 e -α 1 s ds I B .
Using a Taylor expansion at first order, e -α 1 s ≈ 1 -α 1 s, we obtain

I A (t) ≈ A 1 C 1 t 0 e LA 1 H 1 α 1 -α 1 + β 1 + LA 1 H1 s ds ≈ - A 1 e -(α 1 + β 1 + LA 1 H1)t -1 α 1 + β 1 + LA 1 H1
and

I B (t) ≈ - B 1 e -(β 1 + LA 1 H1)t -1 β 1 + LA 1 H1 .
Similarly, we write II = II A + II B , where

II A (t) = A 1 τ r t 0 e -α 1 s 1 -e -β 1 s -LA 1 H 1 se -α 1 s β 1 + LA 1 H 1 e -α 1 s ds ≈ A 1 τ r β 1           t 0 e -α 1 s 1 + LA 1 H 1 β 1 e -α 1 s ds (i) - t 0 e -(α 1 + β 1 + LA 1 H 1 )s 1 + LA 1 H 1 β 1 e -α 1 s ds (ii)           .
For (i), using the change of variable u = e -α 1 s , we obtain

(i) = - 1 α 1 e -α 1 t 1 du 1 + LA 1 H 1 β 1 1 u = - β 1 α 1 LA 1 H 1 ln      1 + LA 1 H 1 β 1 e -α 1 t 1 + LA 1 H 1 β 1     
For small s, se α 1 s ≈ s and using the condition Finally,

LA 1 H 1 β 1 < 1,
II A (t) ≈ - 1 τ r LH 1 α 1 ln      1 + LA 1 H 1 β 1 e -α 1 t 1 + LA 1 H 1 β 1      + A 1 τ r β 1 e -(α 1 + β 1 + LA 1 H 1 )t -1 α 1 + β 1 + LA 1 H 1 - LA 2 1 H 1 τ r β 2 1 e -(2α 1 + β 1 + LA 1 H 1 )t -1 2α 1 + β 1 + LA 1 H 1 + L 2 A 3 1 H 2 1 τ 2 r β 3 1 e -(3α 1 + β 1 + LA 1 H 1 )t -1 3α 1 + β 1 + LA 1 H 1 .
Similarly, we obtain the following expression for

A.2 Integral term of h in phase 2

Our goal is now to compute expression (59). We decompose it into four parts: All computations and approximations are similar except that we integrate between t 1 and t. We obtain

I A (t) ≈ - A 2 C

A.3 Integral term of h in phase 3

Similarly as in phases 1 and 2 we compute the integral in expression (60) and obtain

Figure 1 :

 1 Figure1: Exploration of the AHP-model. A. Time series for the mean voltage h (upper, normalized), the facilitation x (center) and the depression y (lower) simulated from eq. (2). B. Three dimensional phase-space of the AHP-model showing a burst trajectory. The trajectory is decomposed into a QP (green), a burst (blue) and an AHP (pink) phase. The phase-space is divided into 3 regions according to the AHP dynamics: 1) the medium dynamics of hyperpolarization τ 0 = τ mAHP & T 0 = T AHP under and right of the orange surface where the trajectory is highlighted (orange circles). 2) The slow recovery dynamics (τ 0 = τ sAHP & T 0 = 0, region under the purple plan) where the trajectory is highlighted (purple triangles). 3) The fast dynamics (τ 0 = τ & T 0 = 0). C. Segmentation of the voltage time series in burst (blue) and interburst (AHP (pink) and QP (green)). D. Distribution of burst (left, blue), AHP (center, pink) and QP (right, green) durations, extracted from numerical simulations lasting 10 4 s.

Figure 2 :

 2 Figure2: Phase-space of the dynamical system[START_REF] Wilson | Excitatory and inhibitory interactions in localized populations of model neurons[END_REF]. A. Repulsive trajectories from the saddle point S 1 (pink) and S 2 (blue) with corresponding eigenvectors (dashed arrows). B. Inset around S 1 . Real (dashed red arrow) and imaginary (dashed green arrow, see inset) parts of the eigenvectors associated to the complex conjugate repulsive eigenvalues and attractive eigenvector (dashed blue arrow). C. Inset around S 2 . Attractive eigenvectors (dashed blue and black arrows), and repulsive one (dashed red arrow). D. Longest (purple) and shortest (red) bursting trajectories starting outside the basin of attraction.

∑Figure 3 :

 3 Figure 3: Stochastic dynamics in the phase-space of fig.2. A. Exit points (black dots on Σ) of the stochastic system (2) from the basin of attraction (5000s numerical simulations, σ = 3) with longest (purple) and shortest (red) trajectories. B. Escaping stochastic trajectory (black), Inset around the attractor A (yellow star). C. Top view of A and center of mass C M of the exit points (red cross). D. Inset of C.

F 6 FFigure 4 :

 64 Figure 4: Exit probability. A-B. Characteristics crossing the separatrix Σ (the darker the line color is, the lower the value of ψ on Σ is) and distribution exit points obtained from numerical simulations (yellow); visualized with two different angles. C. R 1,1 (t) = ∂ 2 ψ ∂h 2 vs time along the characteristic, obtained numerically. D-E. PDF of the exit points p Σ = Q 0 e -ψ σ on the separatrix Σ compared to the distribution obtained from the stochastic simulations (white histogram); visualized with two different angles.

Figure 5 :

 5 Figure 5: Approximated voltage step function H(t).

Figure 7 :

 7 Figure

Figure 8 :

 8 Figure 8: Influence of the network connectivity J on bursting dynamics. A. (resp.B, C) Mean burst (resp. AHP, QP) duration in seconds from 5000s simulations for J varying from 2.95 to 5.25 and three values of K (left) and L (right) with a fixed noise level (σ = 6).

t t 1 x

 1 (s)y(s)ds = I A + I B + II A + II B .

Table 1 :

 1 Model parameters

		Parameters	Values
	τ	Fast time constant for h	0.05s
	τ mAHP Medium time constants for h 0.15s
	τ sAHP Slow time constants for h	5s
	J	Synaptic connectivity	4.21
	K	Facilitation rate	0.037Hz
	X	Facilitation resting value	0.08825
	L	Depression rate	0.028Hz
	τ r	Facilitation time rate	2.9s
	τ f	Depression time rate	0.9s
	T	Depolarization parameter	0
	σ	Noise amplitude	3
	T AHP	Undershoot threshold	-30

  Medium voltage recovery timescale τ mAHP (s) Slow voltage recovery timescale τ sAHP (s)

	A					B
	Burst duration t i (s)					AHP duration t e -t i (s)
	1 (s) 0.26 C					2 (s)	D
	0.14 0.18 0.22	200	400	600	800	End time of phase 2 t
	0.42 0.46 0.5 0.54 i (s) Burst duration t E					Phase 2 termination threshold h AHP -T AHP 11 12 13 14 e (s) End time of AHP t F
		0.16 0.18 0.2 0.22 0.24 Duration of phase 1 t 1 (s)		-9 -8 -7 -6 -5 -4 AHP termination threshold ε

  we expand the denominator to second -(α 1 + β 1 + LA 1 H 1 )s + β 1 + LA 1 H 1 )t -1 α 1 + β 1 + LA 1 H 1 + LA 1 H 1 β 1 e -(2α 1 + β 1 + LA 1 H 1 )t -1 2α 1 + β 1 + LA 1 H 1 + β 1 + LA 1 H 1 )t -1 3α 1 + β 1 + LA 1 H 1 .

	order to obtain
	(ii) ≈	0	t	e   1 -	LA 1 H 1 β 1	e -α 1 s +	LA 1 H 1 β 1	2	e -2α 1 s	  ds
				≈ -	e -(α 1 -	β 1 LA 1 H 1	2	e -(3α 1

  + LA 1 H 1 )t -1 β 1 + LA 1 H 1 -LA 1 H 1 β 1 e -(α 1 + β 1 + LA 1 H 1 )t -1 α 1 + β 1 + LA 1 H 1 + LA 1 H 1 β 1 2 e -(2α 1 + β 1 + LA 1 H 1 )t -1 2α 1 + β 1 + LA 1 H 1

	II B (t) ≈	B 1 τ r β 1	0	t	  1 -	LA 1 H 1 β 1	e -α 1 s +	β 1 LA 1 H 1	2
												
												 ds
	II B ≈	B 1 τ r β 1	  t +	LA 1 H 1 β 1	e -α 1 t -1 α 1	-	LA 1 H 1 β 1	2	e -2α 1 t -1 2α 1	+	e -(β 1

e -2α 1 s -e -(β 1 + LA 1 H 1 )s

+ LA 1 H 1 β 1 e -(α 1 + β 1 + LA 1 H 1 )s -LA 1 H 1 β 1 2 e -(2α 1 + β 1 + LA 1 H 1 )s   .

  I B (t) ≈ -B 2 C 2 e LA 2 H 2 α 2 (e -(β 2 + LA 2 H2)t -e -(β 2 + LA 2 H 2 )t 1 )β 2 + LA 2 H 2 II A (t) ≈ -1 τ r LH 2 α 2 ln -(α 2 + β 2 + LA 2 H 2 )t -e -(α 2 + β 2 + LA 2 H 2 )t 1 α 2 + β 2 + LA 2 H 2 + β 2 + LA 2 H 2 )t -e -(2α 2 + β 2 + LA 2 H 2 )t 1 2α 2 + β 2 + LA 1 H 2 + β 2 + LA 2 H 2 )t -e -(3α 2 + β 2 + LA 2 H 2 )t 1 3α 2 + β 2 + LA 1 H 2 -2α 2 t -e -2α 2 t 1 2α 2 + e -(β 2 + LA 2 H 2 )t -e -(β 2 + LA 2 H 2 )t 1 β 2 + LA 2 H 2 -LA 2 H 2 β 2 e -(α 2 + β 2 + LA 2 H 2 )t -e -(α 2 + β 2 + LA 2 H 2 )t 1 α 2 + β 2 + LA 2 H 2 + β 2 + LA 2 H 2 )t -e -(2α 2 + β 2 + LA 2 H 2 )t 1 2α 2 + β 2 + LA 2 H 2

		    	1 + 1 +	LA 2 H 2 β 2 LA 2 H 2 β 2 e -α 2 t 1 e -α 2 t	    
	+	A 2 τ r β 2
	-e -(2α 2 + LA 2 2 H 2 τ r β 2 2 L 2 A 3 2 H 2 2 τ 2 r β 3 2 e -(3α 2 II B (t) ≈ B 2 τ r β 2   t -t 1 + β 2 e + LA 2 H 2 β 2 e -α 2 t -e -α 2 t 1 α 2 -LA 2 H 2 β 2 2 LA 2 H 2 2 e -(2α 2

2 e LA 2 H 2 α 2 (e -(α 2 + β 2 + LA 2 H 2 )t -e -(α 2 + β 2 + LA 2 H 2 )t 1 ) α 2 + β 2 + LA 2 H 2 e   .

  -t 2 ) -τ r X(y(t2 -) -1)(e -t -t 2 τ r -1) -τ f (x(t2 -) -X)(e -t -t 2 τ f -1) -(y(t2 -) -1)(x(t2 -) -X)τ f τ r τ f + τ r (e -(t -t 2 ) τ f + τ r τ f τ f -1).A.4 Numerical values of the intermediate and approximation parameters

	t t 2 = X(t Parameters Values x(s)y(s)ds = t t 2 (X + X(y(t2 -) -1)e t 2 -s τ r + (x(t2 -) -X)e +(y(t2 -) -1)(x(t2 -) -X)e (t 2 -s)( 1 τ f + Γ 1 -0.24 Λ 1 -0.18 Γ 2 5.2.10 -6 Λ 2 11.47 Γ 3 -0.91 Λ 3 -0.077 Γ 4 1.4.10 -3 Λ 4 4.4.10 -3 Γ 5 -0.99 Λ 5 0.054 Γ 6 0.91 Λ 6 0.89 Γ 7 1.9.10 -6 Λ 7 0.07 Γ 8 4.0.10 -4 Λ 8 1.8.10 -3 Γ 9 1.3.10 -3 Λ 9 2.8.10 -3	t 2 -s τ f 1 τ r ) ds

Table 2 :

 2 Intermediate parameters

		Parameters	Values
	H 1	Approximation of h for x and y during phase 1 8000
	H1	Initial value of h	400
	H 2	Approximation of h for x and y during phase 2 -1
	h 0	End of phase 1 threshold	400
	h AHP End of phase 2 threshold	-29
		End of AHP threshold	-5
	t 1	End of phase 1 time	200ms
	t 2	End of phase 2 time	1.37s
	A 1	Approximation of x on phase 1 parameter	-0.91
	B 1	Approximation of x on phase 1 parameter	0.99
	C 1	Approximation of y on phase 1 parameter	1.98
	α 1	Approximation of x on phase 1 parameter	297Hz
	β 1	Approximation of y on phase 1 parameter	224Hz
	A 2	Approximation of x on phase 2 parameter	1.16
	B 2	Approximation of x on phase 2 parameter	0.06
	C 2	Approximation of y on phase 2 parameter	0.0017
	α 2	Approximation of x on phase 2 parameter	1.07Hz
	β 2	Approximation of y on phase 2 parameter	0.34Hz

Table 3 :

 3 Approximation parameters


Philippe Berthet

Abel Gomez

Laurence Solovei

Matthieu Gilbert

Sebastien Bommart

Philippe Gaudard

Ludovic Canaud

Pierre Alric

Charles-Henri Marty-Ané

Titanium Implant Failure After Chest Wall Osteosynthesis

niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.