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 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 is a random variable, satisfies the hypothesis made in [1], pp. 15–16: However, Scheutzow exhibited examples of a function ϕ and initial condition 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 . 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 elements. We denote by the set of the p̅ population labels and by the total number of neurons.

Given the neuron i in a population α, the stochastic time evolution of the membrane potential is denoted by . 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 , . The sodium inactivation variable, which is also a proportion of open gates, is denoted by . In the case of the FitzHugh–Nagumo model, the recovery variable is denoted by . Both models feature synaptic variables 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 from the presynaptic neuron j in population γ to the postsynaptic neuron i in population α can be written where is the synaptic reversal potential of the synapse, assumed to be approximately constant across populations, and 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 , of the synapse by the proportion of neurotransmitters available at neuron j. Conductances are positive quantities.

The processes 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 and the fact that the FitzHugh–Nagumo model does not depend on the variables but on the recovery variable only. Conversely, the Hodgkin–Huxley model does not depend on . In what follows we denote by q the vector of in the case of the Hodgkin–Huxley model, and the real in the case of the FitzHugh–Nagumo model. We also note , independent Brownian motions ().

Given a neuron i, the number 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 and , the model assumes that fluctuates around a mean according to a white noise with standard deviation : For a family of independent Brownian motions, independent of the Brownian family , the equations describing the dynamics of the state vector of neuron i in population α are coupled with or 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 and . 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 and are valued in : see Sect. 4.

Sign-Preserving Maximum Conductance Variation

coupled with (2) or (3).

The Mean-Field Limit Models

When making the s tend to infinity, the linear structure of the above N-neuron models w.r.t. the , the linear structure of the dynamics of the , and the mutual independence of the Brownian motions , , lead one to mean-field dynamics. The limit processes are solutions to the SDEs (5) and (9) where denote independent Brownian motions.

For all α in , coupled with or where again stands for in the FitzHugh–Nagumo model and for in the Hodgkin–Huxley model.

Notice that the diffusion coefficients of the play no role in the above mean-field dynamics since one readily sees that

With the same notation, for all α in , coupled with (6) or (7).

As in the simple maximum conductance variation case, the diffusion coefficients of the play no role in the above mean-field dynamics and is given by (8).

Hypotheses

Our hypotheses on the coefficients of the neuron models are the following.

Hypothesis 2.1

On the ion channel models. The function χ is bounded Lipschitz continuous with compact support included in the interval .

The functions , are strictly positive, Lipschitz, and bounded.

On the chemical synapse model. The functions are of the sigmoid type, that is, with suitable positive parameters C, λ, δ depending on α. The parameters , are also positive.

On the membrane model. The drift terms are continuous, one-sided Lipschitz w.r.t. v and Lipschitz w.r.t. q. More precisely, there exist a positive real number L and a positive map such that for all , for all q, in or , v, in ,

The initial conditions. , , , , , , are i.i.d. random variables with the same law as , , , , , , , when . We also assume that and admit moments of any order.

Moreover, the support of the law of belongs to , as well as the support of the laws of , , , for all α in . The support of the law of belongs to .

Remark 2.2

For each neuron population α, the function 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 , , , to live in the interval .

In the case of the FitzHugh–Nagumo model, for all α we have , so that we may choose and .

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 and 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 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 , the concentration processes , , and are well defined and take values in the interval . Each one of these processes is one-dimensional but not Markov since the coefficients of their dynamics depend on 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 . 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 is -valued, is -valued, is a -valued -Brownian motion, and is a -valued Brownian motion.

We aim to exhibit conditions on the coefficients of (11) which imply that the process 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 is valued in the hypercube, say, . The algebraic expressions of the coefficients derived from physical laws may be ‘naturally’ defined only for x in . 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. (11) has a weak solution which does not explode in finite time. In Sect. 4 we develop this argument to show that all systems in Sect. 2 are well posed and that the concentration processes are all valued in .

Assume the following.

Hypothesis 3.2

The locally Lipschitz continuous functions , , β and Γ are such that, on some filtered probability space , there exists a weak solution to (11) which does not explode in finite time. In addition:

The functions , , satisfy, for all u in ,

The functions , , satisfy, for all u in and in ,

The following argument implies that we may limit ourselves to the case . Set . Then, for obviously defined new coefficients , , etc., and a -valued Brownian motion one has If we can prove that takes values in , then the same arguments would show that all the other components enjoy the same property. We therefore consider the system (11) with .

Proposition 3.3

Suppose thata.s. Under Hypothesis3.2

Proof

We limit ourselves to the proof that for all a.s. We can use very similar arguments to get the other inequality.

Let be a positive decreasing function of class with on and on . Let be the first time the process exits from the interval . As on , and Itô’s formula leads to As is decreasing, Hypothesis 3.2(ii) and the definition of imply that . As for all , one has . Finally, is a martingale. Therefore Fatou’s lemma implies Now consider a family of functions which pointwise converge to and satisfy . Lebesgue’s dominated convergence theorem implies In other words, the process is a modification of the null process. As they both are continuous processes, they are indistinguishable (see, e.g., Karatzas and Shreve [8], Chap. 1, Pb. 1.5). □

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 , , , take values in . 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 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.

Well-Posedness of the N-Neuron Models

Proposition 4.1

Under Hypothesis2.1the systems (1) and (4) have unique pathwise solutions on all time interval. In addition, the components of the processes, , , take values in.

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 (respectively, ) with x or y in ; this also holds true at all the other points since, for all and , the value of z such that never belongs to , so that the arguments of the square roots in the diffusion coefficients are strictly positive when x (respectively, y) belongs to .

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

Using the one-sided Lipschitz condition (10) and Itô’s formula, it is an easy exercise to see that is finite for all , from which it readily follows that is finite for all . Therefore the explosion times of (1) and (4) are infinite a.s.

Let us now check that the coordinates of take values in . Their dynamics is of the type for , and where 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 and . 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 Hypothesis2.1and the coercivity conditionthe systems (5) and (9), complemented with (6) or (7), have unique pathwise solutions on all time interval. In addition, all the components of the processtake values in.

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 and the constant , the processes are unique pathwise solutions according to the Yamada and Watanabe theorem (see, e.g., Karatzas and Shreve [8], Chap. 5, Theorem 2.13). Let be an arbitrary continuous function. Consider the system obtained by substituting for 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 the dimension of the state space of the process

Remark 4.4

Notice that we have lumped together the two models we are focusing on, i.e. FitzHugh–Nagumo and Hodgkin–Huxley, since the mathematical tools for handling them are similar: Hypothesis 2.1(iii) does not hold for the classic Hodgkin–Huxley model. However, well-posedness and propagation of chaos can be obtained in this case by adding standard localisation arguments to our proofs.

Let be the canonical process of . Let be the subspace of consisting of the paths of the canonical process such that takes values in for all t, .

Equip the space of probability measures on with the standard Wasserstein (2) metric: where denotes the set of all coupling measures of and .

Given π in , set

Let us construct a contraction map Φ on as follows. For all π in , is the probability law of the process solution to with , or . Notice that the probability law of the process is supported in .

Set , and apply Itô’s formula to . In order to deduce that there exists a positive constant uniform w.r.t. and such that it suffices to use classical arguments, plus the following ingredients: Classical arguments allow one to deduce from (15) that, for some possibly new positive constant , One then finds that Φ is a contraction map (see Sznitman [2]), from which the desired existence and pathwise uniqueness result follows for (9).

The one-sided Lipschitz condition (10);

the fact that is uniformly bounded w.r.t. π in and t in ;

the additional coercivity condition (14) implies As and take values in , all the terms of the type are bounded from above by (the same remarks apply to the diffusion coefficients of and );

the existence of a positive uniform w.r.t. and such that, for all ,

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.

Let denote the solution of (9), with law on . Let , the solution of the N-neuron system (4). Considering the family of Brownian motions in (4), and the set of i.i.d. initial random variables , we introduce a coupling between the and a set of N processes such that for all , all indices i such that are such that are independent copies of . More precisely, for each , such that , is the solution of coupled with or for , and starting at .

Under the hypotheses of Proposition 4.3 we have the following result as regards the propagation of chaos.

Proposition 4.5

Assume that for all, the proportionof neurons in populationγis a nonzero constant independent ofNand denoted: Then there exists a constantsuch that, for allsatisfying the assumption (17), for all set ofp̅indicesamongwith, the vector process, with one component in each population, satisfies

The above -estimate obviously implies that the law of any subsystem of size k with , converges to the law when the tend to ∞. In other words, the reordered system is -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 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 be independent copies of driven by the Brownian motions . Under strong enough Lipschitz hypotheses on ϕ, f and σ, one typically obtains, for some and all , 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), (21), (22)) The fact that the first part 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 is of the order of , as indicated by the inequality (18), instead of in the Lipschitz case.

the processes are positive,

the processes belong to ,

in the dynamics of , the term acts as a mean reverting term, stabilizing the moments of .

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 and rewrite the SDEs (4) and (9) in the following condensed form: for all , where we have introduced the empirical measure , the Brownian motion , and the time-marginal law .

We denote by r the canonical variable on , that we decompose in the following set of p̅ coordinates on : According to Remark 4.4, the diffusion coefficient is defined as and is Lipschitz on the state subspace of the process . The drift coefficient is defined as where 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 ,

Remark 4.6

Notice that the characteristic function is unnecessary in the above definition but, combined with the hypothesis (17), which fixes the constants , 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 is on whose state coordinates are here labeled in .

In all the sequel C is a positive constant which may change from line to line and is independent of N and , but it may depend on T.

Step 1. We prove that the processeshave bounded moments of any positive order. We start with reminding the reader that the CIR processes 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 and by using the following argument. Apply the Itô formula to , , till time , and take expectations to get We then observe that is negative and can be ignored. It then remains to use Hypothesis 2.1 and classical arguments to deduce that is finite for all .

Step 2. A first bound for the random variablesand. 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 or , since calculations with expectation lead to terms of the type , where H is an unbounded random variable correlated with . We therefore postpone taking expectations to Step 3.

Apply Itô’s formula to . We obtain where is a martingale. By Itô isometry and the result in Step 1 above, . Applying the Lipschitz and one-sided-Lipschitz properties for σ and b, we obtain Now, we are interested in . We introduce the empirical measure of the coupling system and write 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 and the are positive, the first term in the right-hand 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 is such that . Indeed, as the are i.i.d. with law , is the conditional expectation for any , where we have set , 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 , since is a martingale,

We now set Defining the drift and diffusion for processes by we have Notice that the processes , defined by is such that . Since and are Lipschitz on , we get Combining again Gronwall’s lemma and integration by parts we obtain Step 3. The bound for. Combining the last inequality (28) with (27), we have where is such that , since 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 . 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 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.