Clarification and Complement to “Mean-Field Description and Propagation of Chaos in Networks of Hodgkin–Huxley and FitzHugh–Nagumo Neurons”

In this note, we clarify the well-posedness of the limit equations to the mean-field N-neuron models proposed in (Baladron et al. in J. Math. Neurosci. 2:10, 2012) and we prove the associated propagation of chaos property. We also complete the modeling issue in (Baladron et al. in J. Math. Neurosci. 2:10, 2012) by discussing the well-posedness of the stochastic differential equations which govern the behavior of the ion channels and the amount of available neurotransmitters.


Introduction
The paper of Baladron et al. [1] studies quite general networks of neurons and aims to prove that these networks propagate chaos in the sense originally developed by Sznitman [2] after the seminal work of Kac on mean-field limits and McKean's work [3] on diffusion processes propagating chaos. As observed by the authors, the membrane potentials of the neurons in the networks they consider are described by interacting stochastic particle dynamics. The coefficients of these McKean-Vlasov systems are not globally Lipschitz. Therefore the classical results of the propagation of chaos theory do not directly apply and a specific analysis needs to be performed. The main theorems (existence and uniqueness of the limit system when the number of particles B O. Faugeras olivier.faugeras@inria.fr 1 tends to infinity, propagation of chaos property) are stated under a fairly general hypothesis on the coefficients. Unfortunately the proof in [1], pp. 24-25, involves an erroneous management of hitting times in combination with a truncation technique, and the limit system may not be well defined under the too general hypothesis used by the authors. Indeed, the following equation, where φ is a bounded and locally Lipschitz function and Z 0 is a random variable, satisfies the hypothesis made in [1], pp. 15-16: However, Scheutzow exhibited examples of a function φ and initial condition Z 0 for which many solutions do exist: see Counterexample 2 in [4] and the remark which follows it. 1 This note restricts the neuron model to the much used variants of the FitzHugh-Nagumo and Hodgkin-Huxley models. Our objective is two-fold: first, we discuss a modeling issue on the diffusion coefficients of the equations describing the proportions of open and closed channels that guarantees that these variables do not escape from the interval [0, 1]. This was not completely achieved in [1] and can be seen as a complement to this paper.
Second, we give a rigorous proof of the propagation of the chaos property.

The Models
In this section we present and discuss the stochastic models considered by Baladron et al. [1] for the electrical activity of p populations of neurons. Each population has a label α and N α elements. We denote by P the set of the p population labels and by N := α∈P N α the total number of neurons. Given the neuron i in a population α, the stochastic time evolution of the membrane potential is denoted by (V i t ). In the case of the Hodgkin-Huxley model, the sodium and potassium activation variables, which represent proportions of open gates along the neuron i are, respectively, denoted by (n i t ), (m i t ). The sodium inactivation variable, which is also a proportion of open gates, is denoted by (h i t ). In the case of the FitzHugh-Nagumo model, the recovery variable is denoted by (w i t ). Both models feature synaptic variables (y i t ) which represent the proportion of available neurotransmitters at the synapses of neuron i.
The synaptic connections between neurons are assumed to be chemical in [1]. We make the same assumption here. This implies that the synaptic current I syn ij from the presynaptic neuron j in population γ to the postsynaptic neuron i in population α can be written where V αγ is the synaptic reversal potential of the j → i synapse, assumed to be approximately constant across populations, and g ij (t) the electric conductance of that synapse. Hence the postsynaptic neuron i belongs to population α and the presynaptic neuron j to population γ . This conductance is the product of the maximal conductance, noted J ij t , of the synapse by the proportion y j t of neurotransmitters available at neuron j . Conductances are positive quantities.
The processes (V i t , n i t , m i t , h i t , w i t , y i t ) are defined by means of the stochastic differential systems (1) or (4) in the N -neuron model section below. The mean-field limit processes are defined in (9) and (5). Well-posedness of those systems is postponed to Sect. 4.

The N -Neuron Model
The variants of the FitzHugh-Nagumo and Hodgkin-Huxley dynamics proposed in Baladron et al. [1] to model neuron networks are all of the two types below; the only differences concern the algebraic expressions of the function F α and the fact that the FitzHugh-Nagumo model does not depend on the variables (n i t , m i t , h i t ) but on the recovery variable (w i t ) only. Conversely, the Hodgkin-Huxley model does not depend on w i t . In what follows we denote by q the vector (n, m, h) of R 3 in the case of the Hodgkin-Huxley model, and the real (w) in the case of the FitzHugh-Nagumo model. We also note ( Given a neuron i, the number p(i) = α denotes the type of the population it belongs to.
Equations (1) and (4) below are those studied in [1]. They correspond to two different models for the maximum conductances. The first one does not respect the positivity constraint while the second one guarantees that these quantities stay positive. In these equations, all quantities which are not time indexed are constant parameters.

Simple Maximum Conductance Variation
For i and j such that p(i) = α and p(j ) = γ , the model assumes that J ij t fluctuates around a mean J αγ according to a white noise with standard deviation σ J αγ : For (B iγ t ) a family of independent Brownian motions, independent of the Brownian family (W i t ), the equations describing the dynamics of the state vector of neuron i in population α are coupled with The reader may wonder about the reason for the square root term and the function χ in the diffusion coefficient of the SDE for the processes x i and y i . The square root arises from the fact that this SDE is a Langevin approximation to a stochastic hybrid, or piecewise deterministic, model of the ion channels. There is a finite (albeit large) number of such ion channels and each of them can be modeled as jump Markov processes coupled to the membrane potential. Altogether ion channels and membrane potentials are described by a piecewise deterministic Markov process which, as shown for example in [6], can be approximated by the solution to the SDE shown in (1) and (4). Hypothesis 2.1(i) below on the function χ implies that the processes x i and y i are valued in [0, 1]: see Sect. 4.

The Mean-Field Limit Models
When making the N α s tend to infinity, the linear structure of the above N -neuron models w.r.t. the (y i

Simple Maximum Conductance Variation
For all α in P, coupled with or where again q α t stands for (w α t ) in the FitzHugh-Nagumo model and for (n α t , m α t , h α t ) in the Hodgkin-Huxley model.
Notice that the diffusion coefficients of the (y α t ) play no role in the above meanfield dynamics since one readily sees that

Sign-Preserving Maximum Conductance Variation
With the same notation, for all α in P, for all γ ∈ P, coupled with (6) or (7). As in the simple maximum conductance variation case, the diffusion coefficients of the (y α t ) play no role in the above mean-field dynamics and E[y α t ] is given by (8).

Hypotheses
Our hypotheses on the coefficients of the neuron models are the following.
We also assume that V α 0 and J αγ 0 admit moments of any order. Moreover, the support of the law of y α 0 belongs to [0, 1], as well as the support of the laws of n α 0 , m α 0 , h α 0 , for all α in P. The support of the law of J αγ 0 belongs to (0, +∞).

Remark 2.2
For each neuron population α, the function S α represents the concentrations of the transmitter released into the synaptic cleft by a presynaptic spike.
Our hypothesis on the support of the function χ is essential to force the proportion processes (y i t ), (n i t ), (m i t ), (h i t ) to live in the interval (0, 1). In the case of the FitzHugh-Nagumo model, for all α we have Finally, the i.i.d. hypothesis in part (iv) is only used in Sect. 4 where it allows simplifications, but it can be relaxed to initial chaos by classical arguments in the propagation of chaos literature.

Remark 2.3
We notice that a one-sided Lipschitz condition also appears in the work by Luçon and Stannat [7] for stochastic particle systems in random environments in which they model one-population networks of FitzHugh-Nagumo neurons. However, their model does not include synaptic interactions as ours does. This has led us in particular to consider square root diffusion coefficients in the dynamics of the synaptic variables which, as shown below, requires specific arguments to prove that the particle systems are well posed and propagate chaos.

Remark 2.4
The boundedness of the functions ρ x and ζ x is a technical hypothesis which simplifies the analysis but can be relaxed, provided that the above models have solutions which do not explode in finite time. However, this comfortable hypothesis is quite reasonable for neuron models since the membrane potentials take values between −100 mV and 100 mV. It is therefore implicitly understood that Lipschitz functions which reasonably fit data within the interval [−100, 100] are extended to bounded Lipschitz functions on the entire real line.

SDEs in Rectangular Cylinders
In the above N -neuron and limit models one requires that, for all i, α and x = n, m, h, the concentration processes (x i t ), (x α t ), (y i t ) and (y α t ) are well defined and take values in the interval [0, 1]. Each one of these processes is one-dimensional but not Markov since the coefficients of their dynamics depend on (V t ) and thus on all the components of the systems (4) and (1). In this context, classical arguments for one-dimensional Markov diffusion processes such as Feller's tests or comparison theorems cannot be used to show that the processes do not exit from [0, 1]. We thus need to develop an ad hoc technique. Instead of focusing on the above neuron models, we rather introduce a more general setting.
Consider the stochastic differential equation where We aim to exhibit conditions on the coefficients of (11) which imply that the process (X t , U t ) takes values in the infinite rectangular cylinder Remark 3.1 Many stochastic models of the type (11) which arise in physics need to satisfy the constraint that the process (X t ) is valued in the hypercube, say, [0, 1] k . The algebraic expressions of the coefficients derived from physical laws may be 'naturally' defined only for x in [0, 1] k . However, one typically can construct continuous extensions of these coefficients on the whole Euclidean space. These extensions may be arbitrarily chosen, provided that they satisfy Hypothesis 3.2 and that Eq. Assume the following.

Hypothesis 3.2
The locally Lipschitz continuous functions b i , A i , β and Γ are such that, on some filtered probability space (Ω, F, P, (F t , t ≥ 0)), there exists a weak solution to (11) which does not explode in finite time. In addition: The following argument implies that we may limit ourselves to the case k = 1. Set U := (X (2) , . . . , X (k) , U). Then, for obviously defined new coefficients β 1 , A , etc., and a R r+k−1 -valued Brownian motion W * one has dX (1) If we can prove that X (1) takes values in [0, 1], then the same arguments would show that all the other components enjoy the same property. We therefore consider the system (11) with k = 1.
Proof We limit ourselves to the proof that 0 ≤ X (1) t for all t ≥ 0 a.s. We can use very similar arguments to get the other inequality.

The Models Are Well-Posed Diffusions in Rectangular Cylinders and Propagate Chaos
In this section, we check that the particle systems and the mean-field limit systems are well posed, and that the components of the processes (y i t ), (x i t ), (y α t ), (x α t ) take values in [0, 1]. Then we prove that the particle systems propagate chaos toward the law of the limit systems (5) and (9).
Our situation differs from the above mentioned Scheutzow's counterexamples [4] in the fact that the interaction kernel is globally Lipschitz and the functions F α are one-sided Lipschitz (they are not only locally Lipschitz). These features of the neuronal models under consideration protect one from non-uniqueness of the solutions. Proof Observe that the coefficients of (1) and (4) are locally Lipschitz. This is obvious for the drift coefficients. In view of the assumption on the function χ , the diffusion coefficients obviously are locally Lipschitz at all point (v, x) (respectively, (v, y)) with x or y in R \ [0, 1]; this also holds true at all the other points since, for all λ 1 > 0 and λ 2 > 0, the value of z such that λ 1 (1 − z) + λ 2 z = 0 never belongs to [0, 1], so that the arguments of the square roots in the diffusion coefficients are strictly positive when x (respectively, y) belongs to [0, 1]. It results from the preceding observation that solutions to (1) and (4) exist and are pathwise unique up to their respective explosion times: see, e.g., Protter [9], Chap. V, Sect. 7, Theorem 38. Set ξ i Using the one-sided Lipschitz condition (10) and Itô's formula, it is an easy exercise to see that E|ξ i T | 2 is finite for all T > 0, from which it readily follows that E sup 0≤t≤T |ξ i t | 2 is finite for all T > 0. Therefore the explosion times of (1) and (4)

Well-Posedness of the N -Neuron Models
for x = n, m, h, and where V t is some real-valued continuous process. Hypothesis 3.2(ii) is satisfied by the drift coefficients of (12) and (13): The desired result follows, by applying Proposition 3.3.

Remark 4.2
The preceding discussion shows that one can get rid of the absolute values in the diffusion coefficients of all the models in Sect. 2.

Well-Posedness of the Mean-Field Limit Models
For the next proposition we slightly reinforce the hypotheses on the functions ρ x and ζ x . The boundedness from below by a strictly positive constant is justified from a biological point of view (see the discussion in [1], Sect. 2.1, p. 5).

Proposition 4.3 Under Hypothesis 2.1 and the coercivity condition
the systems (5) and (9), complemented with (6) or (7) Proof Existence and pathwise uniqueness are obtained by slightly extending the fixed point method developed by Sznitman [2] for particle systems with bounded Lipschitz coefficients. We essentially combine arguments already available in the literature (e.g. see [7] and references therein) and therefore only emphasize the additional minor arguments required by the above neuron models. As exactly the same arguments can be used to treat Eqs. (5) and (9), we limit ourselves to consider the second one.
We start with the following observation. Given the Brownian motions (B αγ , α, γ ∈ P) and the constant J αγ , the processes (J αγ t , α, γ ∈ P) are unique pathwise solutions according to the Yamada and Watanabe theorem (see, e.g., Karatzas and Shreve [8], Chap. 5, Theorem 2.13). Let ϕ(t) be an arbitrary continuous function. Consider the system obtained by substituting ϕ(t) for E[y γ t ] into (9). Similar arguments as in the proof of Proposition 4.1 show that this new system has a unique pathwise solution. Now, we denote by = p 2 + 6p the dimension of the state space of the process Let (r t ) := v α t , j αγ t , ψ α t , w α t , z α t ; α, γ ∈ P be the canonical process of C(0, T ; R ). Let C T be the subspace of C(0, T ; R ) consisting of the paths of the canonical process such that (ψ α t , z α t ) takes values in [0, 1] 4 for all t, α ∈ P.
Equip the space M T of probability measures on C T with the standard Wasserstein (2) metric: where Λ(π 1 , π 2 ) denotes the set of all coupling measures of π 1 and π 2 . Given π in M T , set y (π),γ t dθ ds for all γ ∈ P.
Let us construct a contraction map Φ on M T as follows. For all π in M T , Φ(π) is the probability law of the process , and apply Itô's formula to ( R t ) 2 . In order to deduce that there exists a positive constant K T uniform w.r.t. π 1 and π 2 such that it suffices to use classical arguments, plus the following ingredients: • The one-sided Lipschitz condition (10); • the fact that y (π),α t is uniformly bounded w.r.t. π in M T and t in [0, T ]; • the additional coercivity condition (14) implies (the same remarks apply to the diffusion coefficients of (y (π 1 ),α t ) and (y (π 2 ),α t )); • the existence of a positive K T uniform w.r.t. π 1 and π 2 such that, for all α ∈ P, Classical arguments allow one to deduce from (15) that, for some possibly new positive constant K T , One then finds that Φ is a contraction map (see Sznitman [2]), from which the desired existence and pathwise uniqueness result follows for (9). It remains to again use Proposition 3.3 to get the last part in the statement.

Convergence
In this part, we analyze the convergence of the N -neurons system given in (4) to the mean-field limit described in (9). The convergence of the model (1) to the solution of (5) results from a straightforward adaptation of the following proposition and of its proof. Moreover, as in the proof of Proposition 4.3 we use again Remark 4.4 to shorten the presentation.
The above L 2 (Ω)-estimate obviously implies that the law of any subsystem of size k with p(i α ) = α, converges to the law P ⊗k when the N α tend to ∞. In other words, the reordered system R i α ,N t , α ∈ P , p(i α ) = α, i α ∈ {1, . . . , N} is P-chaotic.
Proof of Proposition 4.5 A short discussion of our methodology. We only present the main ideas of the proof which follows and adapts Sznitman [2] or Méléard [10]. To help the reader follow the lengthy calculations, we start with an explanation of the main differences between our problem where some of the coefficients of our stochastic differential equations are not globally Lipschitz continuous and the classical Lipschitz coefficients framework. In a nutshell, we are dealing with a particle system of the generic form where the Brownian motions W i are independent, and the functions φ, f , and σ are such that one gets existence and strong uniqueness of a solution with finite moments, as well as the existence and strong uniqueness of the mean-field limit Now, let ( X i t ) be independent copies of (X t ) driven by the Brownian motions W i . Under strong enough Lipschitz hypotheses on φ, f and σ , one typically obtains, for some C > 0 and all 0 ≤ t ≤ T , Using independence arguments one readily gets Using Gronwall's lemma, one deduces that This method fails when one of f , σ or φ is not globally Lipschitz. In our case the drift f is not globally Lipschitz, but of the form (see Eqs. (20), The fact that the first part F α of the drift is only one-sided Lipschitz is easy to overcome. To handle the second part we make use of the following three properties: Notice that because in our case the function f is not globally Lipschitz, the convergence rate for E[sup 0≤t≤T α∈P |R i α ,N t − R i α t | 2 ] is of the order of 1/ √ N , as indicated by the inequality (18), instead of 1/N in the Lipschitz case.
Details of our proof. We now turn to the proof of inequality (18).
By the symmetry of the coupled systems, we can fix the index set (1 α , α ∈ P) and rewrite the SDEs (4) and (9) in the following condensed form: for all α ∈ P, where we have introduced the empirical measure μ N , and the time-marginal law P s = P • ( R 1 α s , α ∈ P) −1 .
We denote by r the canonical variable on R , that we decompose in the following set of p coordinates on R p+6 : r := r α , α ∈ P = v α , j αγ ; γ ∈ P , y α , w α , x α ; α ∈ P .
According to Remark 4.4, the diffusion coefficient is defined as and is Lipschitz on the state subspace of the process (V α t , y α t , w α t , x α t ). The drift coefficient is defined as is one-sided Lipschitz in the sense of (10) in Hypothesis 2.1(iii), and k is defined as follows. For any probability measure μ on R l , Remark 4.6 Notice that the characteristic function 1 {η=γ } is unnecessary in the above definition but, combined with the hypothesis (17), which fixes the constants {c γ ; γ ∈ P}, it has the advantage of matching the notations in Eqs. (1) and (4), which helps identifying the interaction kernel in the limit equations. The measure μ(d(r η ; η ∈ P)) is on R whose state coordinates are here labeled in P.
In all the sequel C is a positive constant which may change from line to line and is independent of N and 0 ≤ t ≤ T , but it may depend on T .
Step 1. We prove that the processes V i t have bounded moments of any positive order. We start with reminding the reader that the CIR processes (J iγ ) have bounded moments of any positive order when their initial conditions enjoy the same property (see e.g. Lemma 2.1 in Alfonsi [11]). In view of the Hypotheses 2.1(i) and (iv), one can show that the same is true for the processes (V i ) and ( V i ) by using the following argument. Apply the Itô formula to (V i ) 2q , q ≥ 1, till time τ n = inf{t ≥ 0; |V i t | ≥ n}, and take expectations to get We then observe that ds is negative and can be ignored. It then remains to use Hypothesis 2.1 and classical arguments to deduce that E[(V i t ) n ] is finite for all n ≥ 1.
Step 2. A first bound for the random variables Because of the polynomial form of the non-Lipschitz part of the drift, it is not a good idea to introduce the expectation too early in the calculation of the bound for where H is an unbounded random variable correlated with R 1 α ,N t . We therefore postpone taking expectations to Step 3.
Apply Itô's formula to |R 1 α ,N t − R 1 α t | 2 . We obtain where s is a martingale. By Itô isometry and the result in Step 1 above, sup 0≤t≤T E|M 1 α ,N t | 2 ≤ C. Applying the Lipschitz and one-sided-Lipschitz properties for σ and b, we obtain Now, we are interested in (k[s, We introduce the empirical measure of the coupling system μ N We consider in turn the two terms in the right-hand side of (24). First, from the definition of k in (22) we get Since the (J 1 α γ t ) and the (y i t , i = 1, . . . , N) are positive, the first term in the righthand side is negative. We bound the second term by using Young's inequality: Next we consider the second contribution coming from the right-hand side of (24). By Young's inequality where for any j = 1 α , where we have set k(r α , r η ) = γ ∈P (v α − V αγ )j αγ 1 c γ 1 {η=γ } y η , and the symbol σ stands for sigma algebra (which must not be confused with the above diffusion coefficient). Thus Combining (23) with the last inequalities (25), (24), and (26), By Gronwall's lemma and integration by parts, we have where for all t ∈ [0, T ], since (M 1 α ,N t ) is a martingale, We now set Defining the drift and diffusion for processes y i by Notice that the processes (Z t ), defined by Ce C(t−s) Z s ds.
Step 3. The bound for E[sup 0≤t≤T α∈P |R i α ,N t − R i α t | 2 ]. Combining the last inequality (28) with (27), we have is such that sup 0≤t≤T E(γ t ) 2 ≤ C N , since (δy s ) 2 ≤ 1 a.s. Taking the expectation of the last inequality, we get by applying again Gronwall's lemma in the case of a non-decreasing remainder. Coming back to (27), we get where again is such that sup 0≤t≤T E(β t ) 2 ≤ C N . Using (17), this ends the proof of the proposition.

Conclusion
In this note we have set the work published in [1] on a totally rigorous footing. In doing so we also have shed some new light on the way to incorporate noise in the ion channels equations for the Hodgkin-Huxley model and in the amount of neurotransmitters at the synapses in both the Hodgkin-Huxley and the FitzHugh-Nagumo models.
The techniques in this paper could be extended to a more generic form of interaction kernel k[r; μ] in (22). Notice also that the Hypothesis 2.1(iii) should allow one to prove the convergence in time to equilibrium of the mean-field limits.