Mireille Bossy
email: mireille.bossy@inria.fr

Joaquín Fontbona
email: fontbona@dim.uchile.cl.

Héctor Olivero
email: hector.olivero@uv.cl.

Synchronization of stochastic mean field networks of Hodgkin-Huxley neurons with noisy channels

Keywords: Hodgkin-Huxley neurons, synchronization of neuron networks, mean-field limits, propagation of chaos, stochastic differential equations AMS Subject classification: 60H99, 60K35, 82C22, 82C32, 92B20, 92B25

In this work we are interested in a mathematical model of the collective behavior of a fully connected network of finitely many neurons, when their number and when time go to infinity. We assume that every neuron follows a stochastic version of the Hodgkin-Huxley model, and that pairs of neurons interact through both electrical and chemical synapses, the global connectivity being of mean field type. When the leak conductance is strictly positive, we prove that if the initial voltages are uniformly bounded and the electrical interaction between neurons is strong enough, then, uniformly in the number of neurons, the whole system synchronizes exponentially fast as time goes to infinity, up to some error controlled by (and vanishing with) the channels noise level. Moreover, we prove that if the random initial condition is exchangeable, on every bounded time interval the propagation of chaos property for this system holds (regardless of the interaction intensities). Combining these results, we deduce that the nonlinear McKean-Vlasov equation describing an infinite network of such neurons concentrates, as time goes to infinity, around the dynamics of a single Hodgkin-Huxley neuron with chemical neurotransmitter channels. Our results are illustrated and complemented with numerical simulations.

Introduction

The dynamics of a neuron's voltage is the result of the passage of ions through its membrane. This ion flux takes place through specific proteins which act as gated channels. According to the Hodgkin-Huxley model of a nerve neuron [START_REF] Hodgkin | A quantitative description of membrane current and its application to conduction and excitation in nerve[END_REF], the coupled behavior of the voltage of the neuron V t , with the proportions m t , h t and n t of open channels of the different ions involved in this process (respectively activation Sodium channels, deactivation Sodium channels and activation Potassium channels), can be described by the following system of ordinary differential equations:

V t = V 0 + t 0 F (V s , m s , n s , h s )ds x t = x 0 + t 0 ρ x (V s )(1 -x s ) -ζ x (V s )x s ds (1.1)
where, here and in the sequel, x generically represents the m, n, h components and F : R×[0, 1] 4 → R, defined by

F (V, m, n, h) = I -g K n 4 (V -V K ) -g Na m 3 h(V -V Na ) -g L (V -V L ), (1.2)
represents the effect on the voltage of the ionic channels and of an external current I (assumed constant for simplicity). The rate functions ρ x and ζ x , originally considered in [START_REF] Hodgkin | A quantitative description of membrane current and its application to conduction and excitation in nerve[END_REF], have some generic form given in (2.3) and (2.4) below; see also Figure 2 and Tables 1 and2 for their shape and for biologically meaningful values of the parameters. We refer the reader to Ermentrout and Terman [17] for a concise discussion on the Hodgkin-Huxley (HH in the sequel) model and its deduction, as well as for general background on mathematical models. From a mathematical point of view, system (1.1) defines a rich dynamical system, the properties of which has been extensively studied. As an example of its various possible behaviors, Figure 3 below illustrates different possible responses of system (1.1) to the value of the input current I (no oscillations, oscillations of various types, damping) all other parameters of the model being fixed. See e.g. [17] and Izhikevich [START_REF] Izhikevich | Dynamical Systems in Neuroscience[END_REF], and references therein for detailed accounts on dynamical properties of (1.1) and related neuron models. Lower dimensional dynamics have also been proposed as simpler alternatives to (1.1), the most important ones being the FitzHugh-Nagumo model (FitzHugh [19], Nagumo et al. [42]) and the Morris-Lecar model (Morris and Lecar [41]). These are able to reproduce some of the dynamical features of the HH system (1.1) and are easier to study from the mathematical point of view, but they are less realistic regarding some of its relevant features.

A different approach to model the electric activity of neurons are integrate-and-fire models, introduced in Lapicque [36]. In these models the electric potential evolves according to some ordinary differential equation until it reaches a certain fixed threshold; the neuron then emits a potential spike and the voltage is reset to some reference value, from which its evolution restarts following the same dynamics. We refer the reader to Burkitt [12], [13] for a review of this class of models.

In the last decade, there has been an increasing interest of the mathematical and computational neuroscience communities in understanding the role of stochasticity in neurons' dynamics, as well as in mathematical models for it. We refer the reader to Goldwyn et al. [START_REF] Goldwyn | Stochastic differential equation models for ion channel noise in hodgkin-huxley neurons[END_REF] and to for a discussion on different ways in which randomness might be introduced in the HH model, their biological interpretation and their pertinence. See also [17, Chapter 10] for general background on this issue. Classically, random models arise in the form of finite Markov chains describing a discrete number of open gates which approximate the ion channel dynamics, or by directly introducing Gaussian additive or multiplicative white noise (that is, a Brownian motion or a stochastic integral with respect to it) in the voltage or ion channels dynamics in (1.1). More recently, hybrid (also called piecewise deterministic) Markov processes have also been proposed as microscopic counterparts of the HH or other deterministic models. In this setting, the channel variables are replaced by discrete continuous time processes whose jump rates depend on the voltage, while keeping a continuous description for the latter, see Austin [START_REF] Austin | The emergence of the deterministic hodgkin-huxley equations as a limit from the underlying stochastic ion-channel mechanism[END_REF], Pakdaman et al. [START_REF] Pakdaman | Fluid limit theorems for stochastic hybrid systems with aplication to neuron models[END_REF] and references therein. We also refer the reader to Bossy et al. [START_REF] Bossy | Clarification and complement to "mean-field description and propagation of chaos in networks of hodgkin-huxley and fitzhugh-nagumo neurons[END_REF], Dangerfield et al. [START_REF] Dangerfield | Modeling ion channel dynamics through reflected stochastic differential equations[END_REF], Wainrib [51] and Sacerdote and Giraudo [START_REF] Sacerdote | Stochastic Integrate and Fire Models: A Review on Mathematical Methods and Their Applications[END_REF] for further discussion on stochastic models in this context, the latter one in the case of integrate-and-fire models.

In the present work we consider stochastic versions of the HH model (1.1) which arise as diffusive scaling limits of hybrid models of the type studied in [START_REF] Austin | The emergence of the deterministic hodgkin-huxley equations as a limit from the underlying stochastic ion-channel mechanism[END_REF] and [START_REF] Pakdaman | Fluid limit theorems for stochastic hybrid systems with aplication to neuron models[END_REF]. More precisely, we are interested in networks of N such neurons in mean field interaction, which can be described by a system of stochastic differential equations of the form:

V (i) t = V (i) 0 + t 0 F (V (i) s , m (i) s , n (i) s , h (i) s )ds - t 0 1 N N j=1 J E (V (i) s -V (j) s ) - 1 N N j=1 J Ch y (j) s (V (i) s -V rev )ds, x (i) t = x (i) 0 + t 0 b x (V (i) s , x (i) s )ds + t 0 σ x (V (i) s , x (i) s )dW x,i s , x = m, n, h, y.
Here, (W x,i : i = 1, . . . , N, x = m, n, h, y) are one dimensional Brownian motions and the interaction between neurons account for the effect of electrical and chemical synapses (the biological interpretation of the interactions terms and in particular of the variables y (i) is given in next section). We refer to equation (2.1) below for the explicit form of the system we will consider, and for Hypothesis 2.1 for our assumptions on its coefficients. The collective behavior of neurons, and the way it emerges from their individual features and synaptic activity, is indeed a central question in neuroscience. In particular, considerable efforts have been devoted to understanding synchronization of neurons, an ubiquitous phenomenon seemingly related to the generation of rhythms (such as the respiratory one or the heartbeat) but also to more complex neurologic functionalities. For instance, at the brain level, synchronization has been connected to memory formation, see Axmacher et al. [3], but also to disorders such as epileptic seizures, see Jiruska et al. [START_REF] Jiruska | Synchronization and desynchronization in epilepsy: controversies and hypotheses[END_REF]. Since the neuroscience literature on this topic is huge, it is not our intention to thoroughly comment on it here, and we refer the reader to [17,Chapters 8,[START_REF] Bossy | Modeling the wind circulation around mills with a lagrangian stochastic approach[END_REF], [START_REF] Izhikevich | Dynamical Systems in Neuroscience[END_REF]Chapter 10] for further discussion on biological roles of neuron synchronization, and mathematical approaches to it. For a broad perspective on synchronization, we also refer the reader to Pikovsky et al. [START_REF] Pikovsky | Synchronization: a universal concept in nonlinear sciences[END_REF].

In this paper we are interested in synchronization due to a strong enough coupling between neurons. This phenomenon differs from synchronization owed to common noise addressed e.g. in Marella and Ermentrout [38] (where uncoupled oscillators subject to a common noise are observed to get synchronized) or in Pikovsky [46] (where synchronization results from the action of a random forcing term). In our case, noise is unshared by the interacting neurons and does therefore not contribute to their synchronization; indeed, it actually prevents the perfect asymptotic synchronization of the network. (This phenomenon might be compared to noise-induced deviations from stable cycles in noisy oscillators, see [5,6] for a large deviations approach to that problem.) Moreover, we will understand and quantify synchronization in terms of the empirical variance of the system of neurons and of its dissipation, an approach which does not rely on the stability properties of individual neuron dynamics nor, in particular, on the existence of some oscillatory limiting behavior.

In the case of interacting oscillators, a central mathematical tool is phase reduction. Introduced by Kuramoto [35], it is based on the idea that stable periodic solutions of a nonlinear oscillator can be parametrized by its phase in the limit cycle. Kuramoto's model has proved useful to understand synchronization mechanisms of simple coupled oscillators (see [17] and [START_REF] Izhikevich | Dynamical Systems in Neuroscience[END_REF] ), or for ensembles of population of neurons with intrinsic and extrinsic noise (see Bressloff and Ming Lai [START_REF] Bressloff | Stochastic synchronization of neuronal populations with intrinsic and extrinsic noise[END_REF] and the references therein), or even in the limiting case of infinite oscillators with noise and mean field interaction (see Bertini et al. [7] and the discussion in [51, Chapter 4]). However, to our knowledge, applications of these ideas to networks of HH-type neurons have so far been restricted to small deterministic networks and "weak coupling" regimes (see e.g. Hansel and Mato [START_REF] Hansel | Patterns of synchrony in a heterogeneous hodgkin-huxley neural network with weak coupling[END_REF] and Hansel et al. [START_REF] Hansel | Phase dynamics for weakly coupled hodgkin-huxley neurons[END_REF] ). We refer the reader to Ostojic et al. [START_REF] Ostojic | Synchronization properties of networks of electrically coupled neurons in the presence of noise and heterogeneities[END_REF] for synchronization results in the case of integrate-and-fire networks.

A related question is the asymptotic behavior of networks when the number of neurons tends to infinity. In that sense, networks of N neurons in mean field interaction, in which every neuron experiences a pairwise interaction of strength-order 1/N with each other, provide a mathematically tractable (though not completely realistic) framework to address this question. Indeed, in a mean field network, the evolutions of finitely many neurons are expected to become independent as N goes to infinity, a property known as propagation of chaos. In the case of exchangeable particles, this is equivalent to the convergence of the dynamics of the empirical law of the system to some deterministic flow of probability laws, typically described by a nonlinear McKean-Vlasov partial differential equation (also termed "mean field" equation in this context); see Méléard [39], Sznitman [49] for background on propagation of chaos. We refer the reader to e.g. Faugeras et al. [START_REF] Faugeras | A constructive mean-field analysis of multi population neural networks with random synaptic weights and stochastic inputs[END_REF] for formal derivations of mean field equations for multi type population networks of integrate-and-fire neurons, and respectively to Delarue et al. [START_REF] Delarue | Global solvability of a networked integrateand-fire model of mckean-vlasov type[END_REF] and Perthame and Salort [START_REF] Perthame | On a voltage-conductance kinetic system for integrate and fire neural networks[END_REF] for probabilistic and PDE approaches to the global solvability of that equation (which can in principle have explosive solutions) when the interaction is small. See also Fournier and Löcherbach [START_REF] Fournier | On a toy model of interacting neurons[END_REF] for further recent results on propagation of chaos for integrate and fire models. The propagation of chaos for mean field networks of neurons described by stochastic differential equations, including stochastic, multi type HH and FitzHugh-Nagumo networks, has been addressed in Baladron et al. [4], and then rigorously established in Bossy et al. [START_REF] Bossy | Clarification and complement to "mean-field description and propagation of chaos in networks of hodgkin-huxley and fitzhugh-nagumo neurons[END_REF]. We also refer the reader to Mischler et al. [START_REF] Mischler | On a kinetic fitzhugh-nagumo model of neuronal network[END_REF] for the mean field description of a network of FitzHugh-Nagumo neurons.

In this present work, we establish that, under strong enough electrical connectivity of the network (i.e. large enough J E ), the N neurons get synchronized, up to an error proportional to the channels' noise level σ 2 , at an exponential rate which is independent of N . Moreover, we exhibit a deterministic single-neuron dynamics which is "mimicked", as time goes to infinity, by every neuron of the system (2.1), over short enough moving time-windows, up to an error that vanishes with σ 2 and N -1 . As far as we know, this is a first mathematical result which establishes the synchronization of large networks of neurons. We also establish the propagation of chaos for system (2.1), or its convergence to solutions to a McKean-Valsov equation, for arbitrary parameters of the model (and for slightly more general coefficients than in [START_REF] Bossy | Clarification and complement to "mean-field description and propagation of chaos in networks of hodgkin-huxley and fitzhugh-nagumo neurons[END_REF] in the single population case). This allows us to transfer our synchronization results to the limiting PDE, which can be understood as the description of an infinite network of neurons. Our theoretical results will be complemented with simulations, which in particular point out that synchronization phenomena might also hold for small electrical connectivity and even for pure chemical connectivity (J E = 0).

The remainder of the paper is organized as follows. In Section 2, we detail the model we consider and state precisely our main results. Section 3 is devoted to numerical experiments, both to illustrate our theoretical statements and to explore related phenomena in mean field settings not covered by them (like multi-type neuron populations or chemical-only synapses). In Section 4 we a discuss possible improvements and extensions of our results, and some open questions. The mathematical proofs of our results are given in the Appendix sections.

Model and main results

We start by briefly recalling how chemical and electrical synapses in networks of neurons are modeled (we follow [17, Chapter 7] which we also refer to for further background on synaptic channels).

In chemical synapses, a neurotransmitter is released to the intercellular media (technically the synaptic cleft), from a pre-synaptic neuron to the post-synaptic one through synaptic channels, which are voltage-gated just as ion channels are. With each pre-synaptic neuron we can thus associate a new variable y in [0, 1] which represents its proportion of open synaptic channels at each time. The dynamics of this variable can be modeled in a similar way as those of ion channels, that is, in terms of certain rate functions ρ y and ζ y depending on the membrane potential V of that same neuron, and on some parameters (see (2.5)). The choice of these parameters determines the characteristic (inhibitory or excitatory) of the chemical synapse. Hence, in a fully connected network of N similar neurons, chemical synapses coming from a pre-synaptic neuron j should induce on the voltage V (i) of the post-synaptic neuron i an instantaneous variation at time t of

- J Ch N y (j) t (V (i) t -V rev ),
where y

(j)
t is the proportion of open synaptic channels of neuron j, J Ch ≥ 0 is a constant representing the chemical conductance of the network and V rev is a reference potential. The factor 1 N is introduced in order that the contribution of each incoming synapse to the neuron i has similar weight, which corresponds to a global interaction of mean field type.

On the other hand, the interior of one neuron can be directly connected with another neuron's one through an intercellular channel called gap junction, which allows the constant flow of ions between them, as a result of their possibly different potentials. We thus may assume that pre-synaptic neuron j contributes to the variation of the voltage of post-synaptic neuron i by the amount

- J E N (V (i) s -V (j) s ),
where J E ≥ 0 is the electrical conductance (that can be thought of as a measure of the connectivity of the network) and the factor 1 N appears by similar reasons as before. Connections of this type are termed electrical synapses and are less frequent than chemical ones; on the other hand, they transmit information faster. (See also Hormuzdi et al. [START_REF] Hormuzdi | Electrical synapses: a dynamic signaling system that shapes the activity of neuronal networks[END_REF] for a deeper discussion on electrical synapses.) In all the sequel, for each fixed N we consider a stochastic process X = (X (1) , . . . , X (N ) ) valued in (R 5 ) N , with coordinates

X (i) t = (V (i) t , m (i) t , n (i) t , h (i) t , y (i) t )
given for i = 1, . . . , N and t ≥ 0 by the solution of the system of stochastic differential equations:

V (i) t = V (i) 0 + t 0 F (V (i) s , m (i) s , n (i) s , h (i) s )ds - t 0 1 N N j=1 J E (V (i) s -V (j) s ) - 1 N N j=1 J Ch y (j) s (V (i) s -V rev )ds, x (i) t = x (i) 0 + t 0 ρ x (V (i) s )(1 -x (i) s ) -ζ x (V (i) s )x (i) s ds + t 0 σ x (V (i) s , x (i) s )dW x,i s , (2.1)
where (W x,i : i ∈ N, x = m, n, h, y) are independent one dimensional Brownian motions independent of X 0 and F is defined in (1.2). Notice that, for notational simplicity, the dependence of system (2.1) on N is omitted. Throughout this work, we will additionally make the following assumptions on system (2.1):

Hypothesis 2.1.

1. For x = m, n, h and y, ρ x and ζ x are strictly positive, locally Lipschitz continuous functions defined on R.

2. For x = m, n, h and y, functions σ x : R 2 → R are given by

σ x (v, z) = σ |ρ x (v)(1 -z) + ζ x (v)z|χ(z), (2.2)
with χ : R → [0, 1] a Lipschitz continuous function with support contained in [0, 1] and σ ≥ 0.

3. One has (m

(i) 0 , n (i) 0 , h (i) 0 , y (i)
0 ) ∈ [0, 1] 4 a.s. These assumptions cover, for parameters a x r , a x d > 0, functions of the form

ρ x (V ) = a x r (V -V x r ) 1 -exp (-λ x r (V -V x r )) , ζ x (V ) = a x d exp (-λ x d (V -V x d )) , (2.3)
for x = m, n, and

ρ h (V ) = a h r exp -λ h r (V -V h r ) , ζ h (V ) = a h d 1 + exp -λ h d (V -V h d ) , (2.4)
considered in the original HH model [START_REF] Hodgkin | A quantitative description of membrane current and its application to conduction and excitation in nerve[END_REF], as well as functions

ρ y (V ) = a y r T max 1 + exp (-λ(V -V T )) , ζ y (V ) = a y d (2.5)
associated with synaptic channels in [17, Chapter 7]. Diffusion coefficients σ x defined in terms of the functions ρ x and ζ x as in (2.2) have been considered in [START_REF] Bossy | Clarification and complement to "mean-field description and propagation of chaos in networks of hodgkin-huxley and fitzhugh-nagumo neurons[END_REF], and arise naturally in diffusive scaling limits of the hybrid models studied in [START_REF] Pakdaman | Fluid limit theorems for stochastic hybrid systems with aplication to neuron models[END_REF].

Observe that for functions ρ h , ζ m and ζ n as above, the coefficients of system (2.1) do not satisfy classic conditions for wellposedness. This well-posedness will be proved in Lemma A.1 below, relying on results in [START_REF] Bossy | Clarification and complement to "mean-field description and propagation of chaos in networks of hodgkin-huxley and fitzhugh-nagumo neurons[END_REF] that ensure that under Hypothesis 2.1 the processes (m

(i) t , n (i) t , h (i) t , y (i) t )
remain in [0, 1] 4 . Notice that the absolute value in (2.2) can then be removed.

Synchronization

By synchronization we will understand the dissipation of the empirical variance of the network (2.1) as time goes by. In Figure 1 we show two extreme situations in this regard. On the left, for a small interaction parameter J E and noise σ = 0 we observe a chaotic behavior resulting in an empirical variance of constant order in time. On the right, for large J E and σ = 0 we observe the fast emergence of a coherent evolution implying the dissipation of the empirical variance. Our main result, Theorem 2.3 will provide a quantitative picture of this behavior with respect to noise level σ and the size of the network N , for large enough connectivity J E . Our results require the following additional assumption: (a) Small interaction parameter and noise. Hypothesis 2.2. Hypothesis 2.1 holds and moreover: 4. The parameter g L in (1.2) (termed "leak conductance") is strictly positive.

5. There exists a constant V max 0 > 0 not depending on N such that sup i=1,...,N

|V (i) 0 | ≤ V max 0 a.s.
We also need to introduce notation for some empirical means, namely

V N t = 1 N N i=1 V (i) t , XN t = 1 N N i=1 X (i) t and xN t = 1 N N i=1 x (i) t for x = m, n, h, y.
Moreover, for each N ≥ 1 and t 1 ≥ 0 we let ( X N,t1 t

: t ≥ t 1 ) = ( V N t , m N t , n N t , h N t , y N t ) : t ≥ t 1 denote
the solution of the ordinary differential equation

V N t = V N t1 + t t1 F ( V N s , m N s , n N s , h N s ) -J Ch y N s ( V N s -V rev )ds, x N t = x N t1 + t t1 ρ x ( V N s )(1 -x N s ) -ζ x ( V N s ) x N s ds, x = m, n, h, y.
(2.6) with random initial condition X N,t1 t1 = XN t1 . We are now is position to state our main result about synchronization of the system (2.1): Theorem 2.3. Assume Hypothesis 2.2 holds and that (X

(1) 0 , . . . , X (N ) 0
) is an exchangeable random vector. a) Synchronization. There exist constants J 0 E > 0, C 0 ζ,ρ > 0 and λ 0 > 0 not depending on N ≥ 1, σ ≥ 0 or X 0 , and there exists a time t 0 ≥ 0 not depending on N ≥ 1 or σ ≥ 0, such that for each J E > J 0 E the solution X of (2.1) satisfies, for every t ≥ t 0 and each i ∈ {1, . . . , N }:

E |X (i) t -XN t | 2 ≤ E |X (i) t0 -Xt0 | 2 e -λ 0 (t-t0) + σ 2 C 0 ζ,ρ λ 0 . (2.7)
In particular, lim sup t→∞ E |X

(i) t -XN t | 2 ≤ σ 2 C 0 ζ,ρ λ 0 . b) Synchronized dynamics. Assume J E > J 0 E .
Then, there are constant K 0 , K 0 > 0 depending only on the parameters of the voltage dynamics in (1.2) and, for each δ ≥ 0, constants K δ , K δ > 0 depending on the coefficients in (2.1) and on δ (increasingly) but not on N , such that for every t 1 ≥ t 0 and each i ∈ {1, . . . , N }:

sup t1≤t≤t1+δ E |X (i) t -X t1,N t | 2 ≤ K 0 ∧ 2 K 0 e -λ0(t1-t0) + σ 2 C 0 ζ,ρ λ 0 (1 + δK δ ) + δK δ σ 2 N C 0 ζ,ρ .
(2.8) Some remarks on this result are in order:

Remark 2.4.

(i) The constants C 0 ζ,ρ and t 0 depend explicitly on the coefficients of the system, with the latter possibly depending also on V max 0 . On the other hand, bounds for λ 0 and J 0 E which do not depend on the initial data can also be obtained. The remaining constants are explicit and do not depend on the initial condition. No constant is claimed to be optimal.

(ii) The bound K 0 in Theorem 2.3 b) is deduced from global bounds (which we establish) on the voltage processes and their average, and its role is only to prevent the r.h.s. from growing arbitrarily with δ. The estimate becomes informative as t 1 → ∞ for small enough σ 2 > 0, δ > 0 and N -1 , and for any δ > 0 and N if σ 2 = 0.

(iii) Aside from the assumption that g L > 0, Theorem 2.3 holds regardless of the values of the parameters of voltage dynamics F in (1.2), and in particular if the input current I is replaced by a uniformly bounded function or (suitably measurable) process. (iv) The exchangeability assumption can be removed at the price of adding inside the expectations in (2.7) and (2.8) the averages over i.

Mean field limit

We next address the question of the behavior of system (2.1) as N → ∞. We need to introduce additional notation:

• We denote by P(R × [0, 1] 4 ) the space of Borel probability measures on R × [0, 1] 4 endowed with the weak topology, and by P 2 (R × [0, 1] 4 ) its subspace of probability measures with finite second moment, endowed with the Wasserstein distance W 2 . That is, for all

µ 1 , µ 2 ∈ P 2 (R × [0, 1] 4 ), W 2 (µ 1 , µ 2 ) = inf µ∈Π(µ1,µ2) R 5 |r 1 -r 2 | 2 µ(dr 1 , dr 2 ),
with Π(µ 1 , µ 2 ) the set of probability measures on (R×[0, 1] 4 ) 2 with first and second marginals equal to µ 1 and µ 2 respectively. It is well known that the infimum is attained and that W 2 defines a complete metric on P 2 (R × [0, 1] 4 ) inducing the weak topology, strengthened with the convergence of second moments (see Villani [START_REF] Villani | Optimal transport, old and new, volume 338 of Grundlehren der Mathematischen Wissenschaften[END_REF] for the relevant properties of Wasserstein metrics). • Elements of R × [0, 1] 4 describing the state space of a single neuron's dynamics will be written in the form

(v, u) = (v, (u m , u n , u h , u y )), with u = (u m , u n , u h , u y ) ∈ [0, 1] 4 . • We introduce the function Φ : (R × [0, 1]) × (R × [0, 1] 4 ) → R given by Φ(w, z, v, u) = F (v, u m , u n , u h ) -J E (v -w) -J Ch z(v -V rev )
and, for each channel type x = m, n, h, y, we define functions b x , a x : R ×

[0, 1] 4 → R by b x (v, u) :=ρ x (v)(1 -u x ) -ζ x (v)u x and a x (v, u) :=(ρ x (v)(1 -u x ) + ζ x (v)u x )χ(u x ) (that is, σ 2 x (v, u x ) = σ 2 a x (v, u)). • Given µ ∈ P 2 (R × [0, 1] 4 ) we write µ V = R 5 vµ(dv, du m , du n , du h , du y ) ∈ R, µ x = R 5 u x µ(dv, , du m , du n , du h , du y ) ∈ [0, 1] for x = m, n, h, y and µ = ( µ V , ( µ x ) x=m,n,h,y ) ∈ R × [0, 1] 4 .
• Finally, with δ x denoting the Dirac mass at x ∈ R × [0, 1] 4 , we write

µ N t := 1 N N i=1 δ X (i) t ∈ P 2 (R × [0, 1] 4 ) (2.9)
for the empirical measure of system (2.1) at time t ≥ 0.

Theorem 2.5. Assume Hypothesis 2.2 and moreover that for all N ≥ 1, (X

0 , . . . , X (N ) 0

) are i.i.d. random vectors with (compactly supported) common law µ 0 ∈ P(R × [0, 1] 4 ) not depending on N . a) For each T > 0, the process (µ N t : t ∈ [0, T ]) converges in law on C([0, T ]; P 2 (R × [0, 1] 4 )), when N tends to ∞, to a deterministic and uniquely determined flow of probability measures (µ t : t ∈ [0, T ]) having uniformly bounded compact support. Moreover (µ t : t ≥ 0) in C(R + ; P 2 (R × [0, 1] 4 )) is a global solution (in the sense of distribution) of the non linear McKean-Vlasov Fokker Planck equation

∂ t µ t = ∂ v Φ( µ V t , µ y t , •, •)µ t + x=m,n,h,y 1 2 σ 2 ∂ 2 uxux (a x µ t ) -∂ ux (b x µ t ) (2.10)
with initial condition µ 0 .

b) There is a constant C(T ) > 0 depending on V max 0 , T > 0 and on the coefficients of system (2.1), but not on N , such that

sup t∈[0,T ] E W 2 2 (µ N t , µ t ) ≤ C(T )N -2/5 . (2.11) c) If additionally functions ρ x and ζ x are of class C 2 (R), (or of class C 1 (R) when σ = 0), (µ t : t ≥ 0)
given in part a) is the unique weak solution of (2.10) with initial condition µ 0 which has supports bounded uniformly in time.

Remark 2.6. i) We have not been able to prove uniqueness of weak (measure) solutions to (2.10) in full generality. However, the global weak solution of (2.10) given by Theorem 2.5 a) is uniquely determined. ii) Classically (see [39], [49]), convergence in law of µ N t to µ t for fixed t ≥ 0 implies the asymptotic independence as N → ∞ of any subfamily (X

(1) t , . . . , X (k) t ) of fixed size k ≤ N of system (2.
1) (the propagation of chaos property). Somewhat counterintuitively, this is not incompatible with part a) of Theorem 2.3, even when σ 2 = 0. iii) Parts a) and b) of Theorem 2.5 also hold for general exchangeable µ 0 -chaotic initial conditions (X

(1) 0 , . . . , X (N ) 0 ) (that is, such that µ N 0 converges in law to µ 0 on P 2 (R × [0, 1] 4 ))
, in which case one must add a term of the form CE W 2 2 (µ N 0 , µ 0 ) on the right hand side of (2.11). iv) The first assertion in Theorem 2.5 would be standard if the coefficients in (2.1) were globally Lipschitz. Under the key Hypothesis 2.2 we will be able to reduce the proof to the Lipschitz case. Moreover, this assumption will allow us to take full advantage of the estimates for empirical measures of i.i.d. samples proved in Fournier and Guillin [20], from where the convergence rate (2.11) will stem.

Equation (2.10) can be interpreted as the dynamical description of a system of infinitely many HH neurons in mean field interaction. Thanks to Theorem 2.5 and to the uniformity in N of the results in Theorem 2.3, we can now finally transfer our synchronization results to this infinite dimensional setting. For each

t 1 ≥ 0, define ( X t1,∞ t : t ≥ t 1 ) = ( V ∞ t , m ∞ t , n ∞ t , h ∞ t , y ∞ t ) : t ≥ t 1
as the solution of the ordinary differential equation

V ∞ t = V ∞ t1 + t t1 F ( V ∞ s , m ∞ s , n ∞ s , h ∞ s ) -J Ch y ∞ s ( V ∞ s -V rev )ds, x ∞ t = x ∞ t1 + t t1 ρ x ( V ∞ s )(1 -x ∞ s ) -ζ x ( V ∞ s ) x ∞ s ds, x = m, n, h, y , (2.12
)

with deterministic initial condition X t1,∞ t1 = µ t1 ,
where

(µ t : t ≥ 0) ∈ C([0, ∞), P 2 (R × [0, 1] 4 ))
is the global weak solution of (2.10) with initial condition µ 0 given by Theorem 2.5 a). We have:

Corollary 2.7. Under the assumptions of Theorem 2.5 and for the same constants as in Theorem 2.3, whenever

J E > J 0 E we have: a) For every t ≥ t 0 , W 2 2 (µ t , δ µt ) ≤ W 2 2 (µ t0 , δ µt 0 )e -λ 0 (t-t0) + σ 2 C 0 ζ,ρ λ 0 . (2.13) In particular, lim sup t→∞ W 2 2 (µ t , δ µt ) ≤ σ 2 C 0 ζ,ρ
λ 0 . b) For every t 1 ≥ t 0 and δ ≥ 0 we have:

sup t1≤t≤t1+δ W 2 2 (µ t , δ X t 1 ,∞ t ) ≤ K 0 ∧ 2 K 0 e -λ0(t1-t0) + σ 2 C 0 ζ,ρ λ 0 (1 + δK δ ) .
We next present some numerical simulations which illustrate the validity of our theoretical results (at least from a qualitative point of view) and moreover we explore the behavior of system (2.1) when several neurons subpopulations are considered and when only chemical interaction is present. Furthermore, in view of the numerical experiments, we discuss some of the limitations and possible extensions of our theoretical results.

Numerical Experiments

Inspired in Bossy et al. [START_REF] Bossy | Modeling the wind circulation around mills with a lagrangian stochastic approach[END_REF], we have implemented numerical simulations of system (2.1) by means of an Exponential Projective Euler Scheme (EPES) which we next describe.

For a given time horizon T > 0 and a natural number M , we consider the time grid {t 0 = 0, t 1 = T /M, t 2 = 2T /M, . . . , t k = kT /M, . . . , t M = T }. As initial condition for each neuron in the system we consider independent random variables, uniformly distributed on [-100, 100]×[0, 1] 4 . Given the value of the system at t k , the value for V (i) t k+1 is computed as the exact solution to the ODE

V (i) t = V (i) t k + t t k F ( V (i) s , m (i) t k , n (i) t k , h (i) t k )ds - t 0 1 N N j=1 J E ( V (i) s -V (j) t k ) - 1 N N j=1 J Ch y (j) t k ( V (i) s -V rev )ds.
which is indeed a linear ODE since F is linear in V . To compute x (i)

t k+1 we first solve the SDE x(i) t = x (i) t k + t t k ρ x ( V (i) t k )(1 -x(i) s ) -ζ x ( V (i) t k )x (i) s ds + t t k σ x ( V (i) t k , x (i) t k )dW x,i s , x = m, n, h, y,
which corresponds to an Ornstein-Uhlenbeck process, so that x(i) t k+1 can be exactly simulated. However, since conditionally on x

(i) tj , j ≤ k the law of x(i) t k+1 is Gaussian, {x (i) t k+1 / ∈ [0, 1]} happens with
positive probability, so we are led to define x (i)

t k+1 by projecting x(i) t k+1 onto [0, 1]
, that is:

x (i) t k+1 =      0, x(i) t k+1 ∈ (-∞, 0) x(i) t k+1 , x(i) t k+1 ∈ [0, 1] 1, x(i) t k+1 ∈ (1, +∞).
In Appendix E we prove the convergence in L 2 -norm of the EPES applied to (2.1). We also provide the rate of convergence which is 1/2 as for the classical Euler scheme.

In our simulations we have used as cut-off function (see Hypothesis 2.1-2)

χ(u) = 0.1 exp -0.5 1-(2u-1) 2 u ∈ (0, 1) 0 u / ∈ (0, 1),
whereas, the specific values of the constants we have used are given in Table 1, and the rate functions ρ x and ζ x are given in Table 2 and shown in Figure 2. Although our results hold irrespectively of the value of input current I, we have taken in all simulations I = 25, in which case a noiseless single neuron with the chosen parameters has a limiting regime of sustained oscillations, see Figure 3.

gNa 120 [mS/cm 3 ] gK 36 [mS/cm 3 ] gL 0.3 [mS/cm 3 ] VNa 50 [mV] VK -77 [mV] VL -54.4 [mV]
Table 1: Values for the constants in the function for F . Taken from [17, p.23]

Numerical experiments illustrating our theoretical results

Our first numerical experiments illustrate the results of part a) in Theorem 2.3. In Figure 4 we show one trajectory of the system (2.1) under purely electrical interaction, for different sizes of the network and levels of noise. The first row shows the trajectories of a network of 10 neurons for σ = 0, σ = 0.5 and σ = 1. From the second to the fourth row, the trajectories of networks of 100, 1000 and 10000 neurons, respectively, are shown. The scale of all plots is the same. We observe that the qualitative behavior in terms of σ is the same for all rows: as expected, the noiseless network ultimately reaches perfect synchronization, whereas for σ > 0 the trajectories of the neurons lie in a band whose width Channel type Rate Functions for deactivation Sodium channels increases with σ, as predicted by Theorem 2.3, a). Moreover, the speed at which synchronizations takes place does not depend on the size of the network nor on the level of noise.

ρx(V ) ζx(V ) Sodium (Na) Activation Channels m 0.1(V + 40) 1 -exp - V + 40 10 4 exp - V + 65 18
Sodium (Na) Deactivation Channels h 0.07 exp - V + 65 20 1 1 + exp - V + 35 10
Potassium (K) Activation Channels n 0.01(V + 55) 1 -exp - V + 55 10 0.125 exp - V + 65 80 Neurotransmitter Channels y 5 1 + exp (-0.2(V -2.0)) 0.18
h(V) h(V) (b) Deactivation Sodium Channels.
In our second experiment, we estimate the expected value of the empirical variance of a network of various sizes and for different levels of noise. More precisely we estimate the mean of

SV t = 1 N N i=1 V (i) t -V N t 2 , Sx t = 1 N N i=1 x (i) t -xN t 2 , x = m, n, h.
over 50000 Monte Carlo replica for each value of σ ∈ {0.1, 0.5, 1} and N ∈ {10, 100, 1000, 10000} (we now use σ = 0.1 instead of σ = 0 since in the latter case the obtained plot quickly becomes flat). The computation of the empirical variance for each time step and replica was done using the corrected two-phase algorithm to avoid catastrophic cancellations (see [14]). The results of this experiment appear in Figure 5, where a different variable is presented in each row, from top to bottom: voltage (V ), Sodium activation channels (m), Potassium channels (n) and Sodium deactivation channels (h). Each column corresponds to a different level of noise, increasing from σ = 0.1 on the left, to σ = 0.5 in the middle and to σ = 1 on the right. In each subfigure we show the dissipation of the expected value of the empirical variance for networks of 10, 100, 1000 and 10000 neurons. Just as for one trajectory of the system, we observe again a quick synchronization, now measured in terms of the average dispersion over many trajectories, at speed which does not depend on the noise or the size of the network, with the heights of the peaks increasing with σ. Notice that double peaks are expected from Figure 5 already: even a small dispersion of the phase among different neurons can induce a high dispersion of their voltages and channels right before and after a potential spike is emitted. This dispersion increases with N , but tends to stabilize as N becomes large (notice that the red and green lines in Figure 5 are indistinguishable), consistently with Theorem 2.5.

From this last observation, it is also interesting to point out that the maximum variance over timewindows of fixed length δ > 0 which drift to infinity cannot decrease for every possible value of δ, unless σ 2 = 0. Indeed, in the noisy case the voltage of significantly many neurons can in principle differ from the voltage of the underlying one-neuron dynamics, over time-windows larger than its period, by as much as the whole asymptotic range of the voltages dynamics. Thus, albeit not sharp, the estimates in part b) of Theorem 2.3 and Corollary 2.7 are qualitatively correct.

Beyond Theorem 2.3

We next carry out two type of experiments in situations not covered by our theoretical results on synchronization. In the first one we study the behavior of a more realistic network with several subpopulations of neurons. The dynamics of the i-th neuron in the subpopulation of type α ∈ P , with P denoting the set of subpopulations, is given by

V (i) t = V (i) 0 + t 0 F α (V (i) s , m (i) s , n (i) s , h (i) s )ds - t 0 γ∈P 1 N γ Nγ j=1 J αγ E (V (i) s -V (j) s ) - γ∈P 1 N γ Nγ j=1 J αγ Ch y (j) s (V (i) s -V αγ rev )ds, x (i) t = x (i) 0 + t 0 ρ α x (V (i) s )(1 -x (i) s ) -ζ α x (V (i) s )x (i) s ds + t 0 σ α x (V (i) s , x (i) s )dW x,i s , x = m, n, h, y , t ≥ 0 ,
where N γ is the number of neurons in subpopulation γ. We notice that in this case the electric and chemical conductivity parameters are |P | × |P | matrices. Propagation of chaos for such systems as N → ∞ was proved in [START_REF] Bossy | Clarification and complement to "mean-field description and propagation of chaos in networks of hodgkin-huxley and fitzhugh-nagumo neurons[END_REF] (though under slightly more stringent assumptions on the coefficients).

In Figure 6 we show one trajectory of a network of 100 neurons with two subpopulations, each of them with 50 neurons. On the left (plot (a)), we consider the two subpopulations with different levels of noise and different input current for each of them (different F ), meanwhile the electrical conductance matrix J αγ E is homogeneous, with all the components equal to 1. We observe that the whole network gets synchronized. We believe that Theorem 2.3 can be easily extended to this case (or, more generally, when inf α,γ∈P J αγ E is big enough). In the middle (plot (b)), we observe that, if in addition to considering different subpopulations, the matrix J E is not homogeneous (taking in some entries strictly smaller values than the largest value 1), then synchronization can be observed in each subpopulation but not globally. More precisely, in these examples the two populations synchronize out-of-phase. Finally, on the right (plot (c)), we observe no evidence of synchronization when neurons in one population with same F are electrically connected with two small but different values for J E . This is in line with our theoretical result that thresholds the synchronization of the dynamics for a big enough inf α,γ∈P J αγ E , even if Theorem 2.3 gives only a sufficient condition on J E . According to Kopell and Ermentrout [34], "a small amount of electrical conductance can increase the degree of synchronization far more that a much larger increase in inhibitory conductance". This is consistent with what we have observed in our numerical experiments. The effect of the electrical interaction is stronger and faster than the effect of the chemical interaction. Therefore to appreciate the effect of the chemical synapses, the second type of experiment we have performed concerns only chemical synapses, that is, J E = 0. In Figure 7 we show one trajectory of a network of 100 neurons interacting through inhibitory chemical synapses (in this case, we choose V rev = -75 according to [17, p.163]). We observe an anti-phase synchronization that persists in time (see in plot (a) the transition to the stationary regime in plot (b)), in which two clusters of simultaneously firing neurons emerge. We note that the relative sizes of the two clusters is random and might change from one simulation to another one.

On the other hand, Figure 8 shows one trajectory of a network of 100 neurons interacting through excitatory chemical synapses (with V rev = 0, see [17, p.161]). Some kind of synchronization, similar to the case of purely electrical synapses (see e.g. Figure 6(a)), emerges also here, although the shape of the oscillations is different and the frequency of the spikes is smaller. The numerical experiments presented in Section 3.1 show that, from a quantitative point of view, our theoretical results should still allow for considerable improvements. Indeed, our simulations indicate that the actual global bounds on the voltage, the critical interaction strength above which synchronization happens and the asymptotic discrepancy from synchronization are much smaller than suggested from rough estimates of the bounds in our theoretical results and their proofs. In turn, the actual exponential synchronization rate seems to be much higher. Also, our theoretical results treat the values of the voltage and channels variables jointly, although they are of considerably different orders of magnitude (i.e. channel variables and variances are negligible compared to the voltages ones). The numerical experiments show in turn that the theoretically described behavior happens at each variable's scale.

We must also emphasize that the anti-or out-of-phase-synchronization responses put in evidence in the numerical experiments presented in Section 3.2 are not well captured by the empirical variance criteria proposed in our Theorem 2.3. In those cases, a natural, though challenging strategy would be to extend the phase reduction approach and results developed e.g. in [START_REF] Hansel | Phase dynamics for weakly coupled hodgkin-huxley neurons[END_REF] and [START_REF] Hansel | Patterns of synchrony in a heterogeneous hodgkin-huxley neural network with weak coupling[END_REF] for finite deterministic networks of HH neurons in order to obtain rigorous synchronization results, regardless of the network size, and then in the mean field limit. Another interesting but also challenging question is studying the existence of stationary measures for the McKean-Vlasov limit equation in relation with some characteristic behavior of the system and its possible synchronization regimes, in the vein of recent works of Bertini et al. [8], Giacomin et al. [START_REF] Giacomin | Coherence stability and effect of random natural frequencies in populations of coupled oscillators[END_REF] and Luçon and Poquet [START_REF] Bertini | Dynamical aspects of mean field plane rotators and the kuramoto model[END_REF]. These questions are left for future works.

A Basic properties of the model (2.1)

We start establishing three basic facts about the system of stochastic differential equations (2.1): its (strong) global well-posedness, the fact that the open channels proportion processes stay (as required) in [0, 1] and, finally, and explicit global bound for the voltage processes in terms of a bound for the initial values.

Lemma A.1. Assume Hypothesis 2.1. Then, strong existence and pathwise uniqueness holds for system (2.1). Moreover, a.s. for all t ≥ 0 and every i = 1, . . . , N we have (m

(i) t , n (i) t , h (i) t , y (i) t ) ∈ [0, 1] 4 .
In particular, the absolute value in (2.2) can be removed.

Proof. It is enough to prove the result for deterministic initial data so we assume this is the case. Take M > 0 fixed, and for j = 1, 3, 4 define truncation functions p j M on R by

p j M (x) =    x j x ∈ [-M, M ] M j x ∈ (M, ∞) (-M ) j x ∈ (-∞, -M ). Let X (M ) := (X (1,M ) , . . . , X (N,M ) ) with X (i,M ) t = (V (i,M ) t , m (i,M ) t , n (i,M ) t , h (i,M ) t , y (i,M ) t
), i = 1, . . . , N be defined by

V (i,M ) t = V (i,M ) 0 + t 0 F M (V (i,M ) s , m (i,M ) s , n (i,M ) s , h (i,M ) s ) - 1 N N j=1 J E (V (i,M ) s -V (j,M ) s ) - 1 N N i=1 J Ch p 1 M (y (j,M ) s )(p 1 M (V (i,M ) s ) -V rev )ds, x (i,M ) t = x (i,M ) 0 + t 0 ρ x (p 1 M (V (i,M ) s ))(1 -p 1 M (x (i,M ) s )) -ζ x (p 1 M (V (i,M ) s ))p 1 M (x (i,M ) s )ds + t 0 σ x (p 1 M (V (i,M ) s ), x (i,M ) s
)dW x,i s , x = m, n, h, y, (A.1) where

F M (v, m, n, h) = I-g K p 4 M (n)(p 1 M (v)-V K )-g Na p 3 M (m)p 1 M (h)(p 1 M (v)-V Na )-g L (v-V L ). (A.2)
Is is immediate that the drift coefficients in system (A.1) are Lipschitz continuous. This is less clear in the case of the diffusion coefficients, so we check this point next. Notice that

max (v,u)∈R×[0,1] ρ x (p 1 M (v))(1 -u) + ζ x (p 1 M (v))u ≤ S M := max v∈[-M,M ] ρ x (v) + ζ x (v) < ∞
whereas, thanks to point 2) in Hypothesis 2.1,

min (v,u)∈R×[0,1] ρ x (p 1 M (v))(1 -u) + ζ x (p 1 M (v))u ≥ δ M := min v∈[-M,M ] {ρ x (v), ζ x (v)} > 0.
Therefore, one can find a bounded Lipschitz continuous function g x : R → R + such that g x (s) = √ s on (δ M /2, 2S M ) and rewrite the diffusion coefficients in (A.1) as

σ x (p 1 M (v), u) = σg x (|ρ x (p 1 M (v))(1 -u) + ζ x (p 1 M (v))u|)χ(u).
It is then easily seen that

|σ x (p 1 M (v), u) -σ x (p 1 M (v ), u )| ≤ C M (|u -u | + |v -v |) for some C M > 0 in each of the three cases (u, u ) ∈ [0, 1] 2 , (u, u ) ∈ ([0, 1] 2 ) c and (u, u ) ∈ [0, 1] × [0, 1] c
for any v, v ∈ R. Thus, global pathwise well-posedness for system (A.1) holds.

Thanks to the second assumption in point 4) of Hypothesis 2.1 and the fact that

σ x (v, u) = 0 for (v, u) ∈ R × (0, 1) c and [ρ x (v)(1 -u) -ζ x (v)u][1 (-∞,0] (u) -1 [1,-∞) (u)] ≥ 0 for (v, u) ∈ R 2 ,
we can more apply Proposition 3.3 in [START_REF] Bossy | Clarification and complement to "mean-field description and propagation of chaos in networks of hodgkin-huxley and fitzhugh-nagumo neurons[END_REF] to get that x (1,M ) , . . . , x (N,M ) are confined in [0, 1] for all time (notice that the proof of that result still works if Hypothesis 2.1 i) therein that χ be compactly supported in (0, 1) is replaced by χ being supported in [0, 1]).

We can now use standard arguments to deduce global existence and pathwise uniqueness of a solution to system (2.1). Indeed, setting θ M = inf{t ≥ 0 : |X (M ) t | ≥ M }, using the global Lipschitz character of its coefficients together with Itô calculus and Gronwall's lemma we get for every M > M that a.s. for all t ≥ 0, X (M )

t∧θ M = X (M )
t∧θ M . This implies that θ M > θ M a.s. and allows us to unambiguously define a process X solving (2.1) on the random interval [0, θ), with [START_REF] Ambrosio | Gradient flows: in metric spaces and in the space of probability measures[END_REF] 3 . Using this control on the right hand side of the equations for V (1,M ) , . . . , V (N,M ) in (A.1) and Gronwall's lemma we get

θ := sup M >0 θ M , by X t = X (M ) t for all t ∈ [0, θ M ]. On the other hand, since |p 1 M (z)| ≤ |z| for all z ∈ R, for two constants C 1 , C 2 > 0 not depending on M > 0 we have |F M (v, m, n, h)| ≤ C 1 + C 2 |v| for every (v, m, n, h) ∈ R × [0,
E [|X t∧θ M |] ≤ C(t),
for some constant C(t) > 0 not depending on M . This yields M P [θ M < t] ≤ C(t), whence P [θ < ∞] = 0 letting M and then t ∞ . The statement follows.

Remark A.2.

i) The arguments given in the previous proof also show that each of the functions

σ x is locally Lipschitz on R × [0, 1].
ii) The same proof also works for some extensions of our model. For instance, if independent Brownian motions are added to each of the voltage processes.

We next show that under the additional Hypothesis 2.2, each of the voltage processes is bounded uniformly in time and in N . Below and in all the sequel we denote

V max t,∞ := max i=1,...,N sup s∈[t,∞) |V (i) s |.
We also set

R max := max r,s,u∈[0,1] |I + g Na V Na r + g K V K s + g L V L + J Ch V rev u|.
Proposition A.3. Under Hypothesis 2.2, for every N ≥ 1 and t ≥ 0 we have a.s.

V N t ≤ V max 0 e -gLt + 2R max g L (1 -e -gLt ) and V max t,∞ ≤ V * t := 4R max g L + 2V max 0 e -gLt . (A.3)
As a consequence for every N ≥ 1, there exists at least one invariant law µ N ∞ for the solution to (2.1), namely there exists a solution (X t , t ≥ 0) to (2.1) such that X t has law µ N ∞ for all t ≥ 0 as soon as X 0 has law µ N ∞ . Moreover, this invariant measure is exchangeable.

Remark A.4. i) The bound V * t on V max t,∞ is in general not optimal. For instance, if V max 0 < 2Rmax g L , one can choose V max 0 > 2Rmax
g L +JE and get from the last identity in (A.5) that V max 0,∞ ≤ 4Rmax g L < V * 0 . However, in order to state a synchronization result that holds for a general class of initial conditions V 0 , the fact that the bound V * t does not depend on the electrical connectivity J E and that V * ∞ := lim t→∞ V * t = 4Rmax g L does not depend on the initial condition will be crucial. See point i) in Remark B.4 for a related discussion.

ii) If point 2) of Hypothesis 2.2 does not hold, by slightly modifying the arguments of Lemma A.3 we still can get the a.s. bound

V (i) t ≤ 4R max g L + 2 |V 0 | √ N e -gLt
implying a uniform in N bound for E(V max t,∞ ) if for instance all the random variables V i,N 0 , i = 1, . . . , N , N ≥ 1 are equal in law and have finite second moment. However, we have not been able to fully extend our results to such a framework.

iii) The same arguments also show that a bound like (A.3) holds with V max t,∞ replaced by

V max t,∞ := max i=1,...,N sup s∈[t,∞) | V (i) s |.
That is, the voltages obtained with the EPE scheme are also uniformly bounded.

In the proof of Proposition A.3 and later, we will make use of the the following version of Gronwall's lemma (see e.g. Ambrosio et al. [1, p. 88]).

Lemma A.5. Let θ : [0, +∞) → R be a locally absolutely continuous function and a, b ∈ L 1 loc ([0, +∞)) be given functions satisfying, for λ ∈ R,

d dt θ 2 (t) + 2λθ 2 (t) ≤ a(t) + 2b(t)θ(t) for L ∞ -a.e. t > 0.
Then for every T > 0 we have Proof of Proposition A.3. Setting

e
R (i) s := I + g Na V Na m (i) s 3 h (i) s + g K V K n (i) s 4 + g L V L + J Ch V rev ȳN s , and
A (i) s := g Na m (i) s 3 h (i) s + g K n (i) s 4 + g L + J Ch ȳN s ,
the dynamics of the potential can be written as

V (i) t = V (i) 0 + t 0 R (i) s -A (i) s V (i) s -J E V (i) s + J E V N s ds.
Therefore, we get

V (i) t 2 = V (i) 0 2 + 2 t 0 R (i) s V (i) s -A (i) s V (i) s 2 -J E V (i) s 2 + J E V N s V (i) s ds. (A.4)
and

|V t | 2 = |V 0 | 2 + 2 t 0 N i=1 R (i) s V (i) s -A (i) s V (i) s 2 -J E |V s | 2 + N J E ( V N s ) 2 ds.
Notice that

( V N s ) 2 = 1 N 2 N i,j=1 V (i) s V (j) s ≤ 1 2N 2 N i,j=1 (V (i) s ) 2 + (V (j) s ) 2 = 1 N |V s | 2 ,
which yields

d dt |V t | 2 + 2g L |V t | 2 ≤ 2|R t ||V t |.
By Lemma A.5 we deduce that

|V t | ≤ |V 0 |e -gLt + 2e -gLt t 0 e gLs |R s |ds. Since |R (i) s | ≤ R max , we then get V N t ≤ |V t | √ N ≤ |V 0 | √ N e -gLt + 2R max g L (1 -e -gLt ) ≤ V max 0 e -gLt + 2R max g L (1 -e -gLt ).
which is the first desired inequality. Using this in (A.4) yields

d dt (V (i) t ) 2 + 2(g L + J E )(V (i) t ) 2 ≤ 2|R (i) t + J E V N t ||V (i) t | ≤ 2 R max + J E V max 0 - 2R max g L e -gLt + 2J E R max g L |V (i) t |.
Applying once again Lemma A.5, we obtain

V (i) t ≤ V max 0 e -(gL+JE)t + 2e -(gL+JE)t t 0 e (gL+JE)s R max g L + 2J E g L + J E V max 0 - 2R max g L e -gLs ds = V max 0 e -(gL+JE)t + 2R max g L + 2J E g L (g L + J E ) (1 -e -(gL+JE)t ) + 2 V max 0 - 2R max g L e -gLt -e -(gL+JE)t = 2R max g L g L + 2J E g L + J E + 2 V max 0 - 2R max g L e -gLt + 2R max g L + J E -V max 0 e -(gL+JE)t ≤ 4R max g L + 2V max 0 e -gLt = V * t (A.5)
which implies the asserted bounds on V max t,∞ . Let us now deduce the existence of an invariant distribution which is exchangeable. Let P N t denote the semigroup associated to the solution of (2.1), that is for each X ∈ (R × [0, 1] 4 ) N and B Borel set of (R × [0, 1] 4 ) N , P N t (X , B) = P (X t ∈ B|X 0 = X ) . Consider also the probability measure R N T (λ) on (R × [0, 1] 4 ) N , defined for any law λ as

R N T (λ)(B) = (R×[0,1] 4 ) N 1 T T 0 P N t (X , B)dt λ(dX ).
Since the voltage component is uniformly bounded in time, by (A.5), the solution to (2.1) lies in the compact set

([-4 Rmax gL -2V max 0 , 4 Rmax gL + 2V max 0 ] × [0, 1] 4
) N , and then for any (T M ) ∞, and any λ with compact support, the sequence (R N T M (λ), M ≥ 0) is tight and has a subsequence weakly converging to some probability measure µ N ∞ . According to Krylov-Bogoliubov Theorem, µ N ∞ is invariant for P N t . Let us now choose and exchangeable initial law λ. For any measurable and bounded function ψ, the identity

(R×[0,1] 4 ) N P N t (X , dy)ψ(y)λ(dX ) = (R×[0,1] 4 ) N P N t (X , dy)(ψ • π)(y)λ(dX ).
for any N -permutation π of the coordinates follows directly from the exchangeable structure of the system of equations (2.1). Therefore, R N T M (λ) is exchangeable for any T M , and the corresponding µ N ∞ is exchangeable as the weak limit of exchangeable measures.

B Synchronization: proof of Theorem 2.3 a)

In the sequel, for any locally bounded real function f on R and each R > 0 we will write

f ∞,R := sup v∈[-R,R] |f (v)|.
We will repeatedly use a simple control of the increments of the function F , stated in next lemma for convenience: Lemma B.1. We have

(F (V 1 , m 1 , n 1 , h 1 ) -F (V 2 , m 2 , n 2 , h 2 )) (V 1 -V 2 ) ≤ -g L (V 1 -V 2 ) 2 + 4g K |V 2 -V K ||n 1 -n 2 ||V 1 -V 2 | + 3g Na |V 2 -V Na ||m 1 -m 2 ||V 1 -V 2 | + g Na |V 2 -V Na ||h 1 -h 2 ||V 1 -V 2 |. (B.1) for every m i , n i , h i ∈ [0, 1] and V i ∈ R, i = 1, 2.

Proof. Since

x 4 -y 4 = (x 2 + y 2 )(x + y)(x -y), and

x 3 u -y 3 v = u(x 2 + xy + y 2 )(x -y) + y 3 (u -v),
we get

F (V 1 , m 1 , n 1 , h 1 ) -F (V 2 , m 2 , n 2 , h 2 ) = -(g K n 4 1 + g Na m 3 1 h 1 + g L )(V 1 -V 2 ) -g K (V 2 -V K )(n 2 1 + n 2 2 )(n 1 + n 2 )(n 1 -n 2 ) -g Na (V 2 -V Na )h 1 (m 2 1 + m 1 m 2 + m 2 2 )(m 1 -m 2 ) -g Na (V 2 -V Na )m 3 2 (h 1 -h 2 )
and the asserted bound follows.

The following result is the core of the proof of Theorem 2.3:

Proposition B.2.
For each V * > 0, there are constants J * E > 0 and λ * > 0 not depending on N nor on σ such that for each J E > J * E and any solution X of (2.1) satisfying V max 0,∞ ≤ V * , one has

E |X (i) t -X (j) t | 2 ≤ E |X (i) 0 -X (j) 0 | 2 e -λ * t + σ 2 2C * ζ,ρ λ * ∀ t ≥ 0,
for all i, j ∈ {1, . . . , N }, where

C * ζ,ρ = x=m,n,h,y ρ x ∨ ζ x ∞,V * < ∞. Proof. Let us write ∆V t = V (i) t -V (j) t and ∆x t = x (i) t -x (j)
t . Thanks to the bound (B.1), we have

(∆V t ) 2 = (∆V 0 ) 2 + 2 t 0 [F (V (i) s , m (i) s , n (i) s , h (i) s ) -F (V (j) s , m (j) s , n (j) s , h (j) s )]∆V s ds - t 0 (2J E + 2J Ch ȳN s )(∆V s ) 2 ds ≤ (∆V 0 ) 2 + t 0 8g K |V (j) s -V K ||∆n s ||∆V s | + 6g Na |V (j) s -V Na ||∆m s ||∆V s |ds + t 0 2g Na |V (j) s -V Na ||∆h s ||∆V s |ds - t 0 (2g L + 2J E + 2J Ch ȳN s )(∆V s ) 2 ds ≤ (∆V 0 ) 2 + t 0 ε m (∆m s ) 2 + ε n (∆n s ) 2 + ε h (∆h s ) 2 ds - t 0 2g L + 2J E + 2J Ch ȳN s - 9M 2 Na ε m - 16M 2 K ε n - M 2 Na ε h (∆V s ) 2 ds,
where we have used Young's inequality: ab ≤ ε x a 2 + b 2 4εx for x = m, n, h, y, with ε x > 0 to be chosen later, and where we have set

M Na = g Na max v∈[-V * ,V * ] |v -V Na |, M K = g K max v∈[-V * ,V * ] |v -V K |.
On the other hand, for the channel types x = m, n, h, y, we have

E (∆x t ) 2 = E (∆x 0 ) 2 + 2 t 0 E (1 -x (i) t )(ρ x (V (i) t ) -ρ x (V (j) t ))∆x s ds -2 t 0 E x (i) t (ζ x (V (i) t ) -ζ x (V (j) t ))∆x s ds -2 t 0 E ρ x (V (j) s ) + ζ x (V (j) s ) (∆x s ) 2 ds + t 0 E σ 2 x (V (j) s , x (j) s ) + σ 2 x (V (i) s , x (i) s ) ds.
By our assumptions, for all t ≥ 0 we have fror k = i, j,

σ 2 x (V (k) t , x t ) ≤ σ 2 (1 -x (k) t )ρ x (V (k) t ) + x (k) t ζ x (V (k) t ) ≤ σ 2 ρ x ∨ ζ x ∞,V * .
Using Young's inequality in the same way as before yields

E (∆x t ) 2 ≤ E (∆x 0 ) 2 + t 0 E (L * ρx + L * ζx ) 2 ε x (∆V s ) 2 ds -(2η x -ε x ) t 0 E (∆x s ) 2 ds + 2t σ 2 ρ x ∨ ζ x ∞,V * ,
where L * f denotes the Lipschitz constant on [-V * , V * ] of a locally Lipschitz function f , and where

η * x := inf v∈[-V * ,V * ] {ρ x (v) + ζ x (v)} > 0.
Adding up, we get

E |X (i) t -X (j) t | 2 ≤ E |X (i) 0 -X (j) 0 | 2 - t 0 E 2g L + 2J E - 9M 2 Na + (L * ρm + L * ζm ) 2 ε m - 16M 2 K + (L * ρn + L * ζn ) 2 ε n - M 2 Na + (L * ρ h + L * ζ h ) 2 ε h - (L * ρy + L * ζy ) 2 ε y (∆V s ) 2 ds -(2η m -2ε m ) t 0 E (∆m s ) 2 ds -(2η n -2ε n ) t 0 E (∆n s ) 2 ds -(2η h -2ε h ) t 0 E (∆h s ) 2 ds -(2η y -ε y ) t 0 E (∆y s ) 2 ds + 2t σ 2 C * ζ,ρ .
Define now λ * as the optimal value of the problem max J,εm,εn,ε h ,εy>0

Ψ(J, ε m , ε n , ε h , ε y ) ,
where

Ψ(J, ε m , ε n , ε h , ε y ) := min 2g L +2J- 12M 2 Na + (L * ρ h + L * ζ h ) 2 ε m - 16M 2 K + (L * ρn + L * ζn ) 2 ε n - M 2 Na + (L * ρ h + L * ζ h ) 2 ε h - (L * ρy + L * ζy ) 2 ε y , 2η * m -2ε m , 2η * n -2ε n , 2η * h -2ε h , 2η * y -ε y .
Notice that λ * is strictly positive since Ψ(J, ε m , ε n , ε h , ε y ) can be made so by taking small enough ε x > 0 for x = m, n, h, y and then large enough J > 0. Calling J * E the smallest J > 0 such that (J, ε m , ε n , ε h , ε y ) ∈ arg max Ψ, it follows that for every

J E > J * E , E |X (i) t -X (j) t | 2 ≤ E |X (i) 0 -X (j) 0 | 2 -λ * t 0 E |X (i) s -X (j) s | 2 + 2t σ 2 C ζ,ρ .
Applying Lemma A.5, we obtain

E |X (i) t -X (j) t | 2 ≤ e -λ * t 2 E |X (i) 0 -X (j) 0 | 2 + (e λ * t -1) λ * 2σ 2 C * ζ,ρ 1/2
, and the desired result.

The next result removes the dependance of the previous one on the bound V * , at the price of ensuring exponentially fast synchronization only from some time instant t 0 ≥ 0 on. It will then be easy to deduce part a) of Theorem 2.3.

Theorem B.3.

There are constants J 0 E > 0 and λ 0 > 0 not depending on N ≥ 1, on σ ≥ 0 nor on the initial data, and t 0 ≥ 0 not depending on N ≥ 1 nor on σ ≥ 0, such that for each J E > J 0 E the solution X of (2.1) satisfies, for every t ≥ t 0 ,

E |X (i) t -X (j) t | 2 ≤ E |X (i) t0 -X (j) t0 | 2 e -λ 0 (t-t0) + σ 2 2C 0 ζ,ρ λ 0 , ∀ i, j ∈ {1, . . . , N }, where C 0 ζ,ρ := x=m,n,h,y ρ x ∨ ζ x ∞, 5Rmax g L < ∞ .
Proof. Fix 0 ∈ (0, 1), take t 0 ≥ 0 such that 2V max 0 e -gLt0 ≤ 0 Rmax g L and, conditionally on the sigmafield generated by (X s : s ≤ t 0 ), apply Proposition B.2 to the shifted process X := (X t+t0 :

t ≥ 0) with V * = V * t0 ≤ (4 + 0 ) Rmax g L ≤ 5 Rmax g L .
The proof is then achieved taking expectation in the obtained inequality.

We can now finish the proof of Theorem 2.3. a). Here and in the sequel we denote by SV t and Sx t the empirical variance of voltages and x type channels at time t, respectively:

SV t = 1 N N i=1 V (i) t -V N t 2 and Sx t = 1 N N i=1 x (i) t -xN t 2 .
Proof of Theorem 2.3. a). Applying in the conclusion of Theorem B.3 the elementary identity

1 N 2 N i,j=1 (α i -α j ) 2 = 2 N N k=1 (α k -ᾱN ) 2 for every α 1 , . . . , α N ∈ R, with ᾱN = 1 N N i=1 α i we get E   SV t + x=mn,n,h,y Sx t   ≤ E   SV t0 + x=mn,n,h,y Sx t0   e -λ 0 (t-t0) + σ 2 C 0 ζ,ρ λ 0
in the general case. If, additionally, exchangeability of the initial condition is assumed, the path law of system (2.1) is exchangeable, by pathwise uniqueness. The asserted inequality follows.

Remark B.4.

i) Theorems 2.3 and B.3 show that, for large enough J E , synchronization of the network (2.1) always occurs, as long as the initial voltage V 0 is bounded, but regardless of its actual values. More precisely, the time t 0 > 0 which depends on V max 0 , on Rmax g L and on some arbitrary choice of the parameter 0 > 0, but not on J E , is one possible time after which we can grant that the voltage trajectories stay in some fixed interval not depending on V 0 . Then, after t 0 and if J E was chosen large enough, synchronization occurs at least at the exponential rate λ 0 which depends on coefficients of the system (2.1) but no longer on the initial data. In turn, for large enough J E , Proposition B.2 ensures synchronization from t 0 = 0 on but only if V max 0 is small enough.

ii) Notice that the function Ψ in the proof of Proposition B.2 (and hence the constant λ * therein) increases when its parameter V * decreases, whereas C * ζ,ρ decreases when V * does. Therefore, letting 0 → 0 (or t 0 → ∞) yields the best (by this approach) bounds for the lim sup in Theorem 2.3. Moreover, the largest possible exponential rate λ 0 > 0 and the smallest possible interaction strength J 0 E ≥ 0 that can be obtained (but not necessarily attained) in Theorems 2.3 and B.3 by our approach are λ * and J * E corresponding to V * = 4Rmax g L . These choices are certainly not optimal in general.

C Synchronized dynamics: proof of Theorem 2.3 b)

Our next goal is to prove part b) of Theorem 2.3.

Remark C.1. Proceeding in a similar way as in the proof of Proposition A.3 one checks that the process (2.6)

satisfies d dt | V t | 2 2 + 2g L | V t | 2 2 ≤ 2R max | V t |, which now yields, for any t ≥ t 1 , | V t | ≤ | V t1 |e -gL(t-t1) + 2R max g L (1 -e -gL(t-t1) ).
Applying on

V N t1 = V t1 the first bound in Lemma A.3 we get that | V t | ≤ V max 0 e -gLt + 2Rmax gL for every t ≥ t 1 . Thus, if t 0 ≥ 0 is chosen as in Theorem B.3, we deduce that max sup s∈[t1,∞) | Vs |, sup s∈[t1,∞) | V s | ≤ V * t0 2 ≤ 5R max 2g L . (C.1)
We first prove Proposition C.2. Let t 0 be as in Theorem B.3 and δ > 0. There are constants K 1,δ , K 2,δ > 0 increasingly depending on δ > 0, but not depending on N nor on the initial condition, such that for each t 1 ≥ t 0 ,

E sup t1≤t≤t1+δ | XN,t1 t -X N t | 2 ≤ (V * t0 ) 2 + 4 e -λ 0 (t1-t0) + σ 2 C 0 ζ,ρ λ 0 δK 1,δ + δK 2,δ σ 2 N C 0 ζ,ρ . (C.2)
Proof. For notational simplicity we write in the proof X N t := X N,t1 t . Notice that the average process satisfies the dynamics

V N t = V N t1 + t t1 1 N N i=1 F (V (i) s , m (i) s , n (i) s , h (i) s ) -J Ch ȳN s ( V N s -V rev )ds xN t = xN t1 + 1 N N j=1 t t1 ρ x (V (j) s )(1 -x (j) s ) -ζ x (V (j) s )x (j) s + 1 N N i=1 t t1 σ x (V (j) s , x (j) s )dW x,j s .
Therefore, after some manipulations, we get that

V N t -V t 2 = t t1 1 N N i=1 F (V (i) s , m (i) s , n (i) s , h (i) s ) -F ( V N s , mN s , nN s , hN s ) 2 ds + 2 t t1 F ( V N s , mN s , nN s , hN s ) -F ( V N s , m N s , n N s , h N s ) ( V N s -V N s )ds + t t1 (1 -2J Ch ȳN s )( V N s -V N s ) 2 + 2J Ch ( V N s -V rev )(ȳ N s -y N s )( V N s -V N s )ds = I 1 + I 2 + I 3 .
By Jensen's inequality and the bound (C.1) we have

I 1 ≤ t t1 1 N N i=1 F (V (i) s , m (i) s , n (i) s , h (i) s ) -F ( V N s , mN s , nN s , hN s ) 2 ds ≤ t t1 1 N N i=1 4 g K (n (i) s ) 4 + g Na (m (i) s ) 3 h (i) s + g L (V (i) s -V N s ) 2 ds + t t1 1 N N i=1 4 g Na ( V N s -V Na )h (i) s (m (i) s ) 2 + m (i) s mN s + ( mN s ) 2 (m (i) s -mN s ) 2 ds + t t1 1 N N i=1 4 g Na ( V N s -V K ) (n (i) s ) 2 + (n N s ) 2 n (i) s + nN s (n (i) s -nN s ) 2 ds + t t1 1 N N i=1 4 g Na ( V N s -V Na )( mN s ) 3 (h (i) s -hN s ) 2 ds ≤ K 1 V t t1 SV s + Sm s + Sn s + Sh s ds, with K 1 V explicitly depending on sup v∈[-5Rmax 2g L , 5Rmax 2g L
] max{|v -V Na |, |v -V K |}, g K and g Na . Meanwhile, using (B.1) we get

I 2 ≤ t t1 -2g L ( V N s -V N s ) 2 + 4g K | V N s -V K | ( V N s -V N s ) 2 + (n N s -n N s ) 2 + 3g Na | V N s -V Na | ( V N s -V N s ) 2 + ( mN s -m N s ) 2 + g Na | V N s -V Na | ( V N s -V N s ) 2 + ( hN s -h N s ) 2 ds ≤ K 2 V t t1 ( V N s -V N s ) 2 + (n N s -n N s ) 2 + ( mN s -m N s ) 2 + ( hN s -h N s ) 2 ds,
with K 2 V also depending on those quantities and on g L . By similar arguments, we get

I 3 ≤ K 3 V t t1 ( V N s -V N s ) 2 + (ȳ N s -y N s ) 2 ds for some K 3 V depending on J Ch and on sup v∈[-5Rmax 2g L , 5Rmax 2g L
] |v -V rev |. We thus get:

V N t -V t 2 ≤ K 1 V t t1 SV s + Sm s + Sn s + Sh s ds + KV t t1 | XN s -X N s | 2 ds
for some explicit KV a.s., from where

E sup t1≤s≤t V N s -V s 2 ≤ K 1 V t t1 E SV s + Sm s + Sn s + Sh s ds+ KV t t1 E sup t1≤u≤s | XN u -X N u | 2 ds. (C.3)
On the other hand, for x type channels we get

xN t -x N t = 1 N N j=1 t t1 ρ x (V (j) s )(1 -x (j) s ) -ζ x (V (j) s )x (j) s -ρ x ( V N s )(1 -xN s ) + ζ x ( V N s )x N s ds + t t1 ρ x ( V N s )(1 -xN s ) -ζ x ( V N s )x N s -ρ x ( V N s )(1 -x N s ) + ζ x ( V N s ) x N s ds + 1 N N j=1 t t1 σ x (V (j) s , x (j) s )dW x,j s .
For t ∈ (t 1 , t 1 + δ) we deduce:

(x N t -x N t ) 2 ≤ 3δ t t1 1 N N j=1 ρ x (V (j) s )(1 -x (j) s ) -ζ x (V (j) s )x (j) s -ρ x ( V N s )(1 -xN s ) + ζ x ( V N s )x N s 2 ds + 3δ t 0 ρ x ( V N s )(1 -xN s ) -ζ x ( V N s )x N s -ρ x ( V N s )(1 -x N s ) + ζ x ( V N s ) x N s 2 ds + 3 1 N N i=1 t t1 σ x (V (j) s , x (j) s )dW x,j s 2 .
The previous yields,

E sup t1≤s≤t (x N s -x N s ) 2 ≤ 3δ t t1 E   1 N N j=1 ρ x (V (j) s )(1 -x (j) s ) -ζ x (V (j) s )x (j) s -ρ x ( V N s )(1 -xN s ) + ζ x ( V N s )x N s 2   ds + 3δ t 0 E ρ x ( V N s )(1 -xN s ) -ζ x ( V N s )x N s -ρ x ( V N s )(1 -x N s ) + ζ x ( V N s ) x N s 2 ds + 3E   sup t1≤s≤t 1 N N i=1 s t1 σ x (V (j) u , x (j) u )dW x,j u 2   = I 1 + I 2 + I 3 .
Denoting by L f,R a Lipschitz constant of a function f on [-R, R] and using standard arguments, we get that

I 1 ≤ K x δ(L 2 ρx, 5Rmax 2g L + L 2 ρx+ζx, 5Rmax 2g L ) t t1
E SV s + Sx s ds and that

I 2 ≤ K x δ(L 2 ρx, 5Rmax 2g L + L 2 ρx+ζx, 5Rmax 2g
L ) t t1 E ( V N s -V N s ) 2 + (x N s -x N s ) 2 ds
for all t ∈ (t 1 , t 1 + δ). By Doob's inequality, we moreover obtain

I 3 ≤ 3 • 4E 1 N 2 N i=1 t t1 σ 2 x (V (j) s , x (j) s )ds ≤ 12σ 2 δ N ρ x ∨ ζ x ∞, 5Rmax 2g L .
Summarizing, for the x-type channel we have shown that for all t ∈ (t 1 , t 1 + δ),

E sup t1≤s≤t (x N s -x N s ) 2 ≤ δK x t t1 E SV s + Sx s + sup t1≤u≤s | XN u -X N u | 2 ds+ 12σ 2 δ N ρ x ∨ζ x ∞, 5Rmax 2g L
(C.4) for some constants K x > 0. Putting together (C.3) and (C.4) we get for all t ∈ (t 1 , t 1 + δ) and some constants

K 1 , K 2 > 0, E sup t1≤s≤t | XN s -X N s | 2 ≤ (1 + δ)K 1 t1+δ t1 E   SV s + x=m,n,h,y Sx s   ds + 12σ 2 δ N C 0 ζ,ρ + (1 + δ)K 2 t t1 E sup t1≤u≤s | XN u -X N u | 2 ds,
from where, using Gronwall's inequality, we deduce:

E sup t1≤t≤t1+δ | XN t -X N t | 2 ≤ e K2(1+δ)   K 1 (1 + δ) t1+δ t1 E   SV s + x=m,n,h,y Sx s   ds + 12σ 2 δ N C 0 ζ,ρ   .
We can now use Theorem 2.3 to bound the integral on the r.h.s. With K 1,δ = e K2(1+δ) K 1 (1 + δ) and K 1,δ = 12e K2(1+δ) we get, for all t 1 ≥ t 0 , that

E sup t1≤t≤t1+δ | XN t -X N t | 2 ≤ E   SV t0 + x=m,n,h,y Sx t0   1 λ 0 (1 -e -λ 0 δ )e -λ 0 (t1-t0) K 1,δ + σ 2 C 0 ζ,ρ λ 0 δK 1,δ + δK 2,δ σ 2 N C 0 ζ,ρ ≤ (V * t0 ) 2 + 4 e -λ 0 (t1-t0) + σ 2 C 0 ζ,ρ λ 0 δK 1,δ + δK 2,δ σ 2 N C 0 ζ,ρ since SV t0 ≤ (V * t0 ) 2 .
Proof of Theorem 2.3. b). Notice on hand that, for each t ≥ t 1 , we always have the bounds

|X (i) t -X t1,N t | 2 ≤ 2 SV t + 2| V t -V N t | 2 + 4 ≤ 4(V * t0 ) 2 + 4 ≤ K 0 := 4 5R max g L 2 + 4,
thanks to (C.1) and that V * t0 ≤ 5Rmax g L . On the other hand, combining Proposition C.2 with Theorem 2.3. a) we get for every t

∈ [t 1 , t 1 + δ] that E |X (i) t -X t1,N t | 2 ≤ 2 K 0 e -λ0(t-t0) + σ 2 C 0 ζ,ρ λ 0 (1 + δK 1,δ ) + δK 2,δ σ 2 N C 0 ζ,ρ with K 0 = 5Rmax g L 2
+ 4. The statement follows.

D Propagation of Chaos and synchronization for the McKean-Vlasov limit: proofs of Theorem 2.5 and Corollary 2.7

We first address the asymptotic behavior of the flow of empirical measures (2.9) when N → ∞ and the proof of Theorem 2.5. In particular, we will prove the propagation of chaos property for system (2.1). Following the classic pathwise approach developed in [49] and [39], we first establish:

Theorem D.1. Under the assumptions of Theorem 2.5, we have: a) Let W x , x = m, n, h, y be independent standard Brownian motions and (V 0 , m 0 , n 0 , h 0 , y 0 ) an independent random vector with law µ 0 . There is existence and uniqueness, pathwise and in law, of a solution X = ( V t , m t , n t , h t , y t , t ≥ 0) to the nonlinear stochastic differential equation (in the sense of McKean) with values in R × [0, 1] 4 : c) For each T > 0, let

V t = V 0 + t 0 F ( V s , m s , n s , h s )ds - t 0 J E ( V s -E[ V s ])ds - t 0 J Ch E[ y s ]( V s -V rev )ds, x t = x 0 + t 0 ρ x ( V s )(1 -x s ) -ζ x ( V s ) x s ds + t 0 σ x ( V s , x s )dW x s , x = m, n, h, y (D.1) such that for all t ≥ 0, | V t | ≤ 4R max /g L + 2V max 0 e -
X (i) = ( V (i) t , m (i) t , n (i) t , h (i) t , y (i)
t ) : t ∈ [0, T ] , i = 1, . . . , N be independent copies of the nonlinear process (D.1) each of them driven by the same Brownian motions (W x,i , x = m, n, h, y) and with same initial conditions X (i) 0 = X (i) 0 as the Nparticle system (2.1). Then, there is a constant C(T ) > 0 such that for every N ≥ 1 and i ∈ {1, . . . , N },

E sup 0≤t≤T |X (i) t -X (i) t | 2 ≤ C(T ) N .
Proof. The statements a), b) and c) would be standard if the coefficients in each of the N components of (2.1) were replaced by globally Lipschitz functions of X (i) s and X (j) s , see Theorems 2.2 and 2.3 in [39]. In particular, with functions p j M and F M defined for fixed M > 0 as in Lemma A.1, for any T > 0 there is existence and uniqueness, pathwise and in law, of a solution to the nonlinear stochastic differential equation on [0, T ]:

V M t = V 0 + t 0 F M ( V M s , m M s , n M s , h M s )ds - t 0 J E ( V M s -E[ V M s ])ds - t 0 J Ch E[p 1 M ( y s )](p 1 M ( V M s ) -V rev )ds, x M t = x 0 + t 0 ρ x (p 1 M ( V M s )(1 -p 1 M ( x M s )) -ζ x (p 2 M ( V M s ))p 1 M ( x M s )ds + t 0 σ x (p 1 M ( V M s ), x M s )dW x s , x = m, n, h, y. (D.2) Moreover, letting X (i,M ) = ( V (i,M ) t , m (i,M ) t , n (i,M ) t , h (i,M ) t , y (i,M ) t ) : t ∈ [0, T ] , i = 1, . . . , N
be independent copies of the nonlinear process (D.2) driven by the same Brownian motions (W x,i , x = m, n, h, y) and with same initial conditions X (i) 0 = X (i) 0 as the system (X (1,M ) , . . . , X (N,M ) ) defined in (A.1), we obtain that

E sup 0≤t≤T |X (i,M ) t -X (i,M ) t | 2 ≤ C M (T ) N
for every N ≥ 1 and i ∈ {1, . . . , N }, and some constant C M (T ) > 0.

We notice now that, by Proposition A.3, for M > 0 large enough the system (X (1) , . . . , X (N ) ) is also a solution to the system of equations (A.1). Pathwise uniqueness of the latter yields for all such M > 0 that (X (1) , . . . , X (N ) ) = (X (1,M ) , . . . , X (N,M ) ) on [0, T ], from where

E sup 0≤t≤T |X (i) t -X (i,M ) t | 2 ≤ C M (T ) N (D.3)
for every N ≥ 1 and i ∈ {1, . . . , N }. Furthermore, for any M > 0

P sup 0≤t≤T x (i,M ) t ≥ M + ε ≤ P sup 0≤t≤T x (i,M ) t ≥ M + 2C M (T ) N ε 2
Taking M = 1, letting N → ∞ and then ε → 0 we deduce that x (i,M ) t ≤ 1 a.s. for every t ∈ [0.T ] and i ∈ N. In a similar way, x

(i,M ) t ≥ 0 and | Ṽ (i,M ) t | ≤ V max t,∞ hold
a.s. for every t ∈ [0, T ] and i ∈ N. This implies that for M > 0 large enough but fixed, a solution to (D.2) also solves (D.1), and proves the existence part in a).

We show now that any solution have uniform in time bounded compact support, from which uniqueness in part a) will immediately follow. We shall first consider a solution (U t , q m t , q n t , q h t , q y t ) of (D.1) with explosion time ξ, and we will show that it coincides with

( V M t , m M t , n M t , h M t , y M t , t ≥ 0) for a M big enough. For M > 1, we define τ M = inf{t ≥ 0 : max{|U t |, |q m t |, |q n t |, |q h t |, |q y t |} ≥ M }.
Then we observe that the coefficients of (D.1) applied to (U t , q m t , q n t , q h t , q y t , 0 ≤ t ≤ τ M ) coincide with the truncated coefficients of (D.2) and thanks to the uniqueness property for (D.2) we conclude that almost surely (U, q m , q n , q h , q y ) t∧τ

M = ( V M , m M , n M , h M , y M ) t∧τ M .
In particular, we observe that q x t∧τ M ∈ [0, 1] for x = m, n, h, y, and that τ M = inf{t ≥ 0 :

|U t | ≥ M } for M > 1. Moreover the second order moment E[U 2
t∧τ M ] is uniformly bounded in M , since

U 2 t∧τ M = V 2 0 + 2 t∧τ M 0 U s F (U s , q m s , q n s , q h s )ds -2 t∧τ M 0 J E U s (U s -E[U s ])ds -2 t∧τ M 0 J Ch E[q y s ]U s (U s -V rev )ds,
from where, it is easy to show that

E(U 2 t∧τ M ) ≤ C 1 + C 2 t 0 E(U 2 s∧τ M ),
and therefore, thanks to Gronwall's inequality

E(U 2 t∧τ M ) ≤ C 1 e C2t .
On the other hand E(U 2 t∧τ M ) = E(U 2 t 1 τ M >t ) + M 2 P(τ M ≤ t) and then we can conclude for all t ≥ 0 and all M ≥ 1

P(τ M ≤ t) ≤ C 1 e C2t M 2 . Since τ M
ξ, we conclude that for all t P(ξ ≤ t) = 0, from where it follows that ξ is almost surely infinite. Now, since (U t , q m t , q n t , q h t , q y t ) has no explosion, we apply Proposition 3.3 in [START_REF] Bossy | Clarification and complement to "mean-field description and propagation of chaos in networks of hodgkin-huxley and fitzhugh-nagumo neurons[END_REF] to get that almost surely q x t ∈ [0, 1] for any t > 0. Using this, we derive a more precise bound for the second order moment:

E(U 2 t ) ≤ E(V 2 0 ) + 2 t 0 E(R 2 s ) E(U 2 s ) -g L E(U 2 s )ds,
where as in the proof of Proposition A.3,

R s ≤ R max := max a,b,c∈[0,1] |I + g L V L + g K V K a + g Na V Na b + J Ch V rev c|.
Applying one more time Lemma A.5 we conclude

E(U 2 t ) ≤ E(V 2 0 )e -gLt + 2R max g L (1 -e -gLt ).
Thus, the second moment of any solution of (D.1) is uniformly bounded in time. Moreover, since the initial condition V 0 is bounded, proceeding exactly as in the proof of Proposition A.3 we obtain that

|U t | ≤ 4R max g L + 2V max 0 e -gLt ,
with the same bound V max 0 for V 0 . In conclusion, solutions of (D.1) are non explosive, even more they are uniformly bounded in time. Choosing M > 4R max /g L + 2V max 0 , we get τ M = ∞ almost surely, and for any t ≥ 0, (U, q m , q n , q h , q y ) t = ( V M , m M , n M , h M , y M ) t .

Hence equation (D.1) has a unique solution.

Part b) derives from a direct and easy application of the Ito's formula to compute

E[ψ( X t )] = R×[0,1] 4 ψ(x)µ t (dx)
for a C ∞ c test function ψ, thanks to the fact that the Lebesgue integrals on the right hand side of the Itô formula will be all bounded, since the supports of the laws (µ t : t ≥ 0) are contained in some compact set, and by continuity of coefficients.

Part c) is immediate taking large enough M in (D.3).

We are now in position to prove Proof of of Theorem 2.5. a) We write C T := C([0, T ], R × [0, 1] 4 ). Part c) of Theorem D.1 implies that for each T > 0 and k ≥ 1 the convergence Law(X (1) , . . . X (k) ) → µ ⊗k with µ = Law( X(1) ) holds on the space C k T as N → ∞. By Proposition 2.2. in [49] or Proposition 4.2. in [39], this implies that the empirical measure

µ N := 1 N N i=1 δ X (i) • ∈ P(C T ),
with P(C T ) denoting the space of probability measures on C T endowed with the weak topology, converges in law to the (deterministic) probability measure µ. The first assertion of the theorem follows then from the fact that the mapping associating with ν ∈ P(C T ) its flow (ν t : t ∈ [0, T ]) ∈ C([0, T ]; P(R×[0, 1] 4 )) of one-dimensional time-marginals laws is continuous, together with part b) of Theorem D.1 (notice that C([0, T ]; P(R × [0, 1] 4 )) can be replaced by C([0, T ]; P 2 (R × [0, 1] 4 )) since all the random measures involved have a common compact support).

b) We observe first that for each t ≥ 0 one has

E W 2 2 (µ N t , µ t ) ≤ 2E W 2 2 (µ N t , μN t ) + 2E W 2 2 (μ N t , µ t ) ,
where μN t is the empirical measure of any random i.i.d. sample of the law µ t constructed in the same probability as µ N t . Taking μN

t := 1 N N i=1 δ X(i) t with X (i)
t , i = 1, . . . , N the processes defined in part c) of Theorem D.1 we get for every t ∈ [0, T ] that

E W 2 2 (µ N t , µ t ) ≤ 2 C(T ) N + 2E W 2 2 (μ N t , µ t ) .
On the other hand, we have sup t∈[0,T ] ( |z| q µ t (dz)) 1/q < ∞ for each q ≥ 1, using for instance the bound obtained at the end of the proof of Theorem D.1. We can therefore apply Theorem 1 in [20] with p = 2, d = 5 and a sufficiently large q > 2, to get that E W 2 2 (μ N t , µ t ) ≤ CN -2/5 . The second assertion thus follows. c) In order to prove uniqueness for the McKean-Vlasov equation (2.10), we adapt to our setting a generic argument going back at least to Gärtner [23]. Assume for a while that for each compactly supported ν 0

∈ P(R × [0, 1] 4 )) and (ν * t : t ∈ [0, T ]) ∈ C([0, T ], P 2 (R × [0, 1] 4 )) the linear Fokker- Planck equation ∂ t ν t = ∂ v Φ( (ν * t ) V , (ν * t ) y , •, •)ν t + x=m,n,h,y 1 2 σ 2 ∂ 2 uxux (a x ν t ) -∂ ux (b x ν t ) (D.4)
has at most one weak solution with supports bounded uniformly in t ∈ [0, T ]. By similar arguments as in Lemma A.1, strong well-posedness holds for the stochastic differential equation:

V * t = V * 0 + t 0 F (V * s , m * s , n * s , h * s )ds - t 0 J E (V * s -(ν * s ) V )ds - t 0 J Ch (ν * s ) y (V * s -V rev )ds, x * t = x * 0 + t 0 ρ x (V * s )(1 -x * s ) -ζ x (V * s )x * s ds + t 0 σ x (V * s , x * s )dW x s , x = m,
: t ∈ [0, T ]) = (µ t : t ∈ [0, T ]). This yields, for all t ∈ [0, T ], that µ t = law(V * t , m * t , n * t , h * t , y * t )
, for the process defined as in (D.5), with ν * s = µ s for all s ∈ [0, T ]. In other words, this process solves the nonlinear stochastic differential equation (D.1). From Theorem D.1 we conclude that (µ t : t ∈ [0, T ]) = (law( Xt ) : t ∈ [0, T ]), that is, there is uniqueness of solutions in C([0, T ], P(R × [0, 1] 4 )) of (2.10) having uniformly bounded support.

Hence, in order to conclude the proof of Theorem 2.5 it is enough to show that, given functions α, β ∈ C([0, T ], R) and ν 0 ∈ P 2 (R × [0, 1] 4 ) there is at most one solution (ν t : t ∈ [0, T ]) ∈ C([0, T ], P(R × [0, 1] 4 )) with support bounded uniformly in [0, T ], to the distribution formulation of equation (D.4)

ψ(t, v, u)ν t (dv, du) = ψ(0, v, u)ν 0 (dv, du) - t 0 Φ(α s , β s , v, u)∂ v ψ(s, v, u) + ∂ s + x=m,n,h,y 1 2 σ 2 a x ∂ 2 uxux + b x ∂ ux ψ(s, v, u) ν s (dv, du) ds (D.6)
for all t ∈ [0, T ] and for an extended class of test function

ψ ∈ C 1,1,2 b ([0, T ] × R × [0, 1] 4
). Let ρ x and ζ x denote compactly supported functions coinciding with ρ x and ζ x on some compact set K ⊂ R containing the supports of the measures ν V t for t ∈ [0, T ], and define σ x , a x and b x in terms of them in a similar way as σ x , a x and b x were defined in terms of ρ x and ζ x . For a given t > 0, consider the following Cauchy problem in R 5 : for all (s, v, u)

∈ [0, t) × R × R 4 , ∂ s -Φ(α s , β s , v, u)∂ v + x=m,n,h,y 1 2 σ 2 a x ∂ 2 uxux + b x ∂ ux f t (s, v, u) =0, f t (t, v, u) = ψ(v, u). (D.7)
By the Feynman-Kac formula (see e.g. Karatzas and Shreve [START_REF] Karatzas | Brownian Motion and Stochastic Calculus[END_REF]), if a solution

f t ∈ C b ([0, t]×R 5 )∩ C 1,1,2 b ([0, t) × R × R 4 ) exists, then it is given by f t (s, v, u) := E(ψ(X s,v,u t )) (D.8) where (X s,v,u r := (V r , m r , n r , h r , y r ) : r ∈ [s, t]
) is the unique (pathwise and in law) solution in [s, t] of the stochastic differential equation:

V r = v + r s F (V θ , m θ , n θ , h θ )dθ - r s J E (V θ -α θ )ds - t s J Ch β θ (V θ -V rev )dθ, x r = u x + r s ρ x (V θ )(1 -x θ ) -ζ x (V θ )x θ dθ + r s σ x (V θ , x θ )dW x θ , x = m,
= f t in (D.6) yields ψ(v, u)ν t (dv, du) = E(ψ(X s,v,u t ))ν 0 (dv, du)
for all ψ ∈ C 2 0 (R 5 ), which uniquely determines ν t . Notice that when σ = 0, the required regularity for φ and for f turns from C 1,1,2 to C 1,1,1 and the Feymann Kac formula in the argument can be replaced by the characteristics formula. The proof of part c) is complete. d) This is immediate from parts b) and d) of Theorem D.1 .

Proof of Corollary 2.7. Recall first that, for any ν ∈ P 2 (R × [0, 1] 4 ) and w ∈ R × [0, 1] 4 ), one has

W 2 2 (ν, δ w ) = |z -w| 2 ν(dz).
Moreover, for every t ≥ t 1 and N ≥ 1 it holds by exchangeability that:

E |X (i) t -X t1,N t | 2 = E W 2 2 (µ N t , δ X t 1 ,N t ) .
Therefore, it is enough to prove that, for any

t 1 ≥ 0, sup t1≤t≤t1+δ E W 2 2 (µ t , δ X t 1 ,∞ t ) -W 2 2 (µ N t , δ X t 1 ,N t ) → 0
as N → ∞. Given t ≥ t 1 and N ≥ 1, let π N t (dz, dz ) be a coupling between µ t and µ N t . Then, for some constant C > 0 not depending on t ≥ t 1 nor on N ≥ 1, we have

W 2 2 (µ t , δ X t 1 ,∞ t ) -W 2 2 (µ N t , δ X t 1 ,N t ) = π N t (dz, dz ) |z -X t1,∞ t | 2 -|z -X t1,N t | 2 ≤ C |z -z |π N t (dz, dz ) + | X t1,∞ t -X t1,N t |
since the supports of µ t and µ N t and the processes X t1,∞ t and X t1,N t are uniformly bounded in t ≥ t 1 and N . The latter property also allows us to write the dynamics in (2.6) and (2.12) using globally Lipschitz coefficients. Thanks to Gronwall's lemma this yields the estimates

sup t1≤t≤t1+δ | X t1,∞ t -X t1,N t | ≤ C δ | X t1,∞ t1 -X t1,N t1 | = C δ | µ t1 -µ N t1 | ≤ C δ |z -z |π N t1 (dz, dz ) for some constant C δ > 0 not depending on N . Since |z-z |π N t (dz, dz ) ≤ |z -z | 2 π N t (dz, dz ) 1/2 ,
by taking the above couplings to be optimal for W 2 , we get the estimate

sup t1≤t≤t1+δ E W 2 2 (µ t , δ X t 1 ,∞ t ) -W 2 2 (µ N t , δ X t 1 ,N t ) ≤ C sup t1≤t≤t1+δ E 1/2 W 2 2 (µ t , µ N t )
for some C > 0. We conclude thanks to Theorem 2.5.

E Strong Convergence Rate Result for the Exponential Projective Euler Scheme (EPES)

The main object of this section is to prove the convergence of the numerical scheme presented in Section 3 to the model (2.1) and establish the following rate of convergence Proposition E.1. Asumme Hypothesis 2.2, if χ(x) = O(x(1 -x)), then there exists a constant C depending on the parameters of the system, but independent of ∆t, such that for any i = 1, . . . , N :

E V (i) t -V (i) t 2 + x=m,n,h,y E |x (i) t -x (i) t | 2 ≤ C∆t.
We decompose the proof of this proposition in several preliminary results. The next result follows from the uniform bound for V (i) t

(see iii) in Remark A.4) and some standard arguments on local approximation of SDEs, so we omit the proof.

Lemma E.2. Under Hypothesis 2.2, there exists a constant C depending on the parameters of the system, but independent of ∆t such that sup i=1,...,N E V

(i) t -V (i) η(t) 2 ≤ C∆t 2 , sup i=1,...,N E x(i) t -x (i) η(t) 2 ≤ C∆t.
Next we establish a the key step in the convergence of the scheme, namely that, with extremely high probability, the processes x (i) and x(i) coincide.

Lemma E.3. Asumme Hypothesis 2.2, if χ(x) = O(x(1 -x)), then there exists a constant C depending on the parameters of the system, but independent of ∆t, such that sup i=1,...,N x=m,n,h,y

P

x(i) t / ∈ [0, 1] ≤ exp -C ∆t .

Remark E.4. It is not difficult to see that

E x(i) t -x (i) t 2 = E x(i) t -x (i) t 2 1 {x (i) t / ∈[0,1]} ≤ sup j=1,...,N E 2(x j t ) 2 + 1 P x(i) t / ∈ [0, 1] ≤ K exp - C 2∆t .
Notice that the RHS above tends to zero faster than any power of ∆t when ∆t → 0.

Proof of Lemma E.3. We first notice that conditional to F η(t) , x(i) corresponds to an Ornstein-Uhlenbeck process, therefore its law is Gaussian with known conditional mean and conditional variance given by

E η(t) x(i) t = x (i) η(t) exp -(ρ x + ζ x ) ( V (i) η(t) )(t -η(t)) + ρ x ( V (i) η(t) ) (ρ x + ζ x ) ( V (i) η(t) ) 1 -exp -(ρ x + ζ x ) ( V (i) η(t) )(t -η(t)) , Var η(t)
x(i)

t = σ 2 x ( V (i) η(t) , x (i)
η(t) ) 2 (ρ x + ζ x ) ( V (i) η(t) ) 1 -exp -2 (ρ x + ζ x ) ( V (i) η(t) )(t -η(t)) .
Observe that the conditional variance is strictly positive if t > η(t), x η(t) = 1 the diffusions coefficient vanish, and in that case the solution to the ODE for x(i) remains in [0, 1] almost surely, we can restrict ourselves to the case x (i) η(t) ∈ (0, 1). Using the inequality for Gaussian concentration, conditional to F η(t) , we have

P η(t)
x(i) On the other hand, E η(t)

x(i) t is a weighted mean between to quantities in [0, 1], therefore .

E η(t) x(i) t ≥ x (i) η(t) ∧ ρ x ( V (i) η(t) ) (ρ x + ζ x ) ( V (i) η(t) ) , hence P η(t) x(i) t ≤ 0 ≤ 1 2 exp   -ρ 2 x ( V (i) η(t) ) (ρ x + ζ x ) 2 ( V (i) η(t) )σ 2 x ( V (i) η(t) , x (i)
(E.1) To bound the first exponential in the right-hand side of the last inequality, we notice that since the process V (i) is uniformly bounded and σ is bounded, we can easily exhibit a constant C 1 > 0 independent of i such that ρ 2

x ( V (i)

η(t) ) (ρ x + ζ x ) 2 ( V (i) η(t) )σ 2 x ( V (i) η(t) , x (i)
η(t) ) ≥ C 1 .
For the second term in the right-hand side of (E.1), since x 2 /χ(x) 2 is bounded from below in (0, 1), there exists C 2 > 0, such that

x (i) η(t) 2 σ 2 x ( V (i) η(t) , x (i)
η(t) ) ≥ C 2 ,
from which we conclude

P x(i) t ≤ 0 ≤ exp - C 1 ∧ C 2 ∆t .
An analogous computation shows

P x(i) t ≥ 1 ≤ exp - C 1 ∧ C 2 ∆t .
The last preliminary step in the proof of Proposition E. Then there exists a constant C depending on the parameters of the system, but independent of ∆t, such that u(t) ≤ u(η(t)) + C∆t 2 e C∆t . (E.2)

Proof. Thanks to the boundedness of the processes, drift and diffusion coefficients, that is b x and σ x , behave like Lipschitz functions, just as in the proof of Lemma A.1. Then, thanks to Itô formula and pivoting with in drift and diffusion with the point ( V

(i) s , x(i) s ) E (x (i) t - x(i) t ) 2 ≤ E (x (i) η(t) -x (i) η(t) ) 2 + 2 t η(t) E (x (i) s -x(i) s ) b x (V (i) s , x (i) s ) -b x ( V (i) s , x(i) s ) ds + 2 t η(t) E (x (i) s -x(i) s ) b x ( V (i) s , x(i) s ) -b x ( V (i) η(t) , x(i) s ) ds + t η(t) 2E σ x (V (i) s , x (i) s ) -σ x ( V (i) s , x(i) s ) 2 + 2E σ x ( V (i) s , x(i) s ) -σ x ( V (i) η(t) , x (i)
η(t) ) 2 ds.
from where the Lipchitz property of the coefficients, Lemma E.2 to bound the terms involving the local error and some classical arguments lead to

E (x (i) t - x(i) t ) 2 ≤ E (x (i) η(t) - x(i) η(t) ) 2 + C t η(t) E (x (i) s -x(i) s ) 2 + E (V (i) s -V (i) s ) 2 ds + C∆t 2 .
On the other hand, for the voltage error we obtain first the a.s. bound

V (i) t -V (i) t 2 ≤ V (i) η(t) -V (
|V (i) s -V (i) s | 2 + (V (i) s + V rev ) 2 1 N N j=1
|y (j) s -y Thanks to the exchangeability of the particles, it follows that

E   1 N N j=1 |y (j) s -y (j) η(t) | 2   = E y (i) s -y (i) η(t) 2 , E   1 N N j=1 |V (j) s -V (j) η(t) | 2   = E V (i) s -V (i) η(t) 2 ,
and then, since the processes are uniformly bounded, we get that

E V (i) t -V (i) t 2 ≤ E V (i) η(t) -V (i) η(t) 2 + C t η(t) E V (i) s -V (i) s 2 ds+ + C t η(t) E V (i) s -V (i) η(t) 2 + x=m,n,h,y E |x (i) s -x (i) η(t) | 2 ds ≤ E V (i) η(t) -V (i) η(t) 2 + C t η(t) E |V (i) s -V (i) s | 2 ds+ + C t η(t) E V (i) s -V (i) s 2 + C∆t 2 + x=m,n,h,y E |x (i) s -x(i) s | 2 + C∆tds ≤ E V (i) η(t) -V (i) η(t) 2 + C t η(t) E V (i) s -V (i) s 2 + x=m,n,h,y E |x (i) s -x(i) s | 2 ds + C∆t 2 .
We can summarize the previous computations as

u(t) ≤ u(η(t)) + C t η(t)
u(s)ds + C∆t 2 , from where we conclude thanks to Gronwall's inequality.

Proof of Proposition E.1 . From the previous Lemma, denoting u k = u(t k ) we obtain the following recurrence relationship:

u k+1 ≤ Au k + B, u 0 = 0, A = e C∆t , B = C∆t 2 e C∆t .
Iterating this inequality, it is easy to conclude that

u k+1 ≤ B A k+1 -1 A -1 ≤ C∆t 2 e C∆t k+1
e C∆t -1 = C∆t 2 e Ct k+1 e C∆t -1 .

But when ∆t → 0, we have e C∆t -1 ∼ C∆t, and therefore u k+1 ≤ C∆t. Inserting this in (E.2), we conclude

E V (i) t -V (i) t 2 + x=m,n,h,y E |x (i) t -x (i) t | 2 ≤ E V (i) t -V (i) t 2 + x E |x (i) t - x(i) t | 2 + E |x (i) t -x (i) t | 2 ≤ C∆t + x P x(i) t / ∈ [0, 1] ,
from where the statement follows, applying Lemma E.3.

  .(b) Large interaction parameter and no noise.

Figure 1 :

 1 Figure 1: Two typical situations of the evolution of the network (2.1).

  for Chemical Neurotransmitter channelsy(V) y(V) (d) Synaptic Channels.

Figure 2 :

 2 Figure 2: Charateristic plot of rate functions ρx and ζx.

  .(c) I = 100.

  . (d) I = 200.

Figure 3 :

 3 Figure 3: Reponses of the model (1.1) depending on the input current I: no oscillations if I = 0 (a); large amplitude and low frequency oscillations if I = 10 (b); small amplitude and high frequency oscillations if I = 100 (c); damped oscillations if I = 200 (d).

  .(a) N = 10, σ = 0.

  .(b) N = 10, σ = 0.5.

  .(c) N = 10, σ = 1.

  (d) N = 100, σ = 0.

  (e) N = 100, σ = 0.5.

  (f) N = 100, σ = 1. (g) N = 1000, σ = 0. (h) N = 1000, σ = 0.5. (i) N = 1000, σ = 1. (j) N = 10000, σ = 0. (k) N = 10000, σ = 0.5. (l) N = 10000, σ = 1.

Figure 4 :

 4 Figure 4: Synchronization of a network under pure electrical interaction for different network sizes N and noise levels σ.

  of the Empirical Variance of the Voltage for =0.1 N=10 N=100 N=1000 N=10000 (a) V , σ = 0.1.

  of the Empirical Variance of the Voltage for =0.5 N=10 N=100 N=1000 N=10000 (b) V , σ = 0.5.

  of the Empirical Variance of the Voltage for =1 N=10 N=100 N=1000 N=10000 (c) V , σ = 1.

  Value of the Empirical Variance for m for =0.1 N=10 N=100 N=1000 N=10000 (d) m, σ = 0.1.

  Value of the Empirical Variance for m for =0.5 N=10 N=100 N=1000 N=10000(e) m, σ = 0.5.

  Value of the Empirical Variance for m for =1 N=10 N=100 N=1000 N=10000 (f) m, σ = 1.

  Value of the Empirical Variance for n for =0.1 N=10 N=100 N=1000 N=10000 (g) n, σ = 0.1. Value of the Empirical Variance for n for =0.5 N=10 N=100 N=1000 N=10000 (h) n, σ = 0.5. Value of the Empirical Variance for n for =1 N=10 N=100 N=1000 N=10000 (i) n, σ = 1.

  Value of the Empirical Variance for h for =0.1 N=10 N=100 N=1000 N=10000 (j) h, σ = 0.1.

  Value of the Empirical Variance for h for =0.5 N=10 N=100 N=1000 N=10000(k) h, σ = 0.5.

  Value of the Empirical Variance for h for =1 N=10 N=100 N=1000 N=10000 (l) h, σ = 1.

Figure 5 :

 5 Figure 5: Dissipation of the empirical variance for different level of noise. First row: expected empirical variance of the voltage; second to fourth rows: expected empirical variance of the Sodium activation channels (m), the Potassium channels (n) and the Sodium deactivation channels (h) respectively. JCh = 0 and JE = 1 in all simulations.

  Populations are different, the interaction is homogeneous.

  Populations are different, the interaction is heterogeneous.

  Populations are equal, the interaction is heterogeneous.

Figure 6 :

 6 Figure 6: Trajectories for network of 100 neuron with two subpopulations.

  Voltage of 100 neurons after 900 [ms].

Figure 7 :

 7 Figure 7: Trajectories for network of 100 neuron with inhibitory chemical synapses

  (a) Voltage of 100 neurons at the beginning of the simulation.

  Voltage of 100 neurons after 900 [ms].

Figure 8 :

 8 Figure 8: Trajectory of a network of 100 neurons with excitatory chemical synapses

e

  λT |θ(T )| ≤ θ 2 (0) + sup t∈[0,T ] λt |b(t)|dt.

  gLt almost surely. b) (µ t := law( Xt ) : t ≥ 0) is a weak solution globally defined in C((0, +∞]; P 2 (R × [0, 1] 4 )) of the McKean-Vlasov equation (2.10).

  ) | 2 ds.

Table 2 :

 2 Rate functions for the dynamics of the channels. Taken from [17, p.23] for the Sodium and Potassium channels, and from page [17, pp.160,163] for the neurotransmitter channel.

  It follows using Itô's formula that a unique weak solution to equation (D.4) with uniformly bounded supports does exist, and is given by ν t = law(V * Now, any solution (µ t : t ∈ [0, T ]) in C([0, T ], P(R × [0, 1] 4 )) of (2.10) with uniformly bounded supports also solves the linear equation (D.4) with (ν * t

	n, h, y,
	(D.5)
	with (V * 0 , m * 0 , n * 0 , h * 0 , y t , m * t , n * t , h * t , y *

* 0 ) independent of the Brownian motions W x and with law ν 0 . Moreover, one can check that x * t ∈ [0, 1] a.s. for all t ∈ [0, T ] and that the process (V * t : t ∈ [0, T ]) is bounded. t ) for all t ∈ [0, T ].

  n, h, y.Moreover, for v chosen in some fixed compact set, this solution is bounded independently of s ∈ [0, t], and one has x r ∈ [0, 1] for all r ∈ [s, t]. Hence, under the assumption that σ > 0, ρ x and ζ x are of class C 2 (R), one can moreover prove, following the lines ofFriedman [22, p.124], that the function f t defined by (D.8) actually is of class C 1,1,2

b ([0, t) × R × R 4

) and solves the Cauchy problem (D.7). Putting ψ

  Since for t small enough1 -exp -2 (ρ x + ζ x ) ( V )(t -η(t)) ≤ 2 (ρ x + ζ x ) ( Vand t -η(t) ≤ ∆t, we can bound the conditional variance, and then it follows

	t ≤ 0 = P η(t)	   	x(i) t -E η(t) Var η(t)	x(i) t x(i) t	≤	-E η(t) Var η(t)	x(i) t x(i) t	    ≤	1 2	exp	  -	E η(t) Var η(t)	2 x(i) x(i) t t	   .
			(i) η(t) (i) η(t) )(t -η(t)),
	P η(t)	x(i) t ≤ 0 ≤	1 2	exp	  -	E η(t) x ( V σ 2 (i) η(t) , x x(i) t	2

  1 is the following Lemma E.5. Under hypotheses of Proposition E.1, consider u(t) := E V

	(i) t	-V	(i) t	2	+	E |x	(i) t -	x(i)
					x=m,n,h,y

t | 2 .

Partially supported by Nucleo Milenio NC 130062, Beca Chile Postdoctorado and Fondecyt Postdoctoral Grant 3180777.

This research was partially supported by the supercomputing infrastructure of NLHPC (ECM-02), Conicyt.