### Abstract

We derive the mean-field equations arising as the limit of a network of interacting spiking neurons, as the number of neurons goes to infinity. The neurons belong to a fixed number of populations and are represented either by the Hodgkin-Huxley model or by one of its simplified version, the FitzHugh-Nagumo model. The synapses between neurons are either electrical or chemical. The network is assumed to be fully connected. The maximum conductances vary randomly. Under the condition that all neurons’ initial conditions are drawn independently from the same law that depends only on the population they belong to, we prove that a propagation of chaos phenomenon takes place, namely that in the mean-field limit, any finite number of neurons become independent and, within each population, have the same probability distribution. This probability distribution is a solution of a set of implicit equations, either nonlinear stochastic differential equations resembling the McKean-Vlasov equations or non-local partial differential equations resembling the McKean-Vlasov-Fokker-Planck equations. We prove the well-posedness of the McKean-Vlasov equations, i.e. the existence and uniqueness of a solution. We also show the results of some numerical experiments that indicate that the mean-field equations are a good representation of the mean activity of a finite size network, even for modest sizes. These experiments also indicate that the McKean-Vlasov-Fokker-Planck equations may be a good way to understand the mean-field dynamics through, e.g. a bifurcation analysis.

**Mathematics Subject Classification (2000): **
60F99, 60B10, 92B20, 82C32, 82C80, 35Q80.

##### Keywords:

mean-field limits; propagation of chaos; stochastic differential equations; McKean-Vlasov equations; Fokker-Planck equations; neural networks; neural assemblies; Hodgkin-Huxley neurons; FitzHugh-Nagumo neurons### 1 Introduction

Cortical activity displays highly complex behaviors which are often characterized by the presence of noise. Reliable responses to specific stimuli often arise at the level of population assemblies (cortical areas or cortical columns) featuring a very large number of neuronal cells, each of these presenting a highly nonlinear behavior, that are interconnected in a very intricate fashion. Understanding the global behavior of large-scale neural assemblies has been a great endeavor in the past decades. One of the main interests of large-scale modeling is characterizing brain functions, which most imaging techniques are recording. Moreover, anatomical data recorded in the cortex reveal the existence of structures, such as the cortical columns, with a diameter of about 50 μm to 1 mm, containing the order of 100 to 100,000 neurons belonging to a few different types. These columns have specific functions; for example, in the human visual area V1, they respond to preferential orientations of bar-shaped visual stimuli. In this case, information processing does not occur at the scale of individual neurons but rather corresponds to an activity integrating the individual dynamics of many interacting neurons and resulting in a mesoscopic signal arising through averaging effects, and this effectively depends on a few effective control parameters. This vision, inherited from statistical physics, requires that the space scale be large enough to include sufficiently many neurons and small enough so that the region considered is homogeneous. This is, in effect, the case of the cortical columns.

In the field of mathematics, studying the limits of systems of particle systems in interaction has been a long-standing problem and presents many technical difficulties. One of the questions addressed in mathematics was to characterize the limit of the probability distribution of an infinite set of interacting diffusion processes, and the fluctuations around the limit for a finite number of processes. The first breakthroughs to find answers to this question are due to Henry McKean (see, e.g. [1,2]). It was then investigated in various contexts by a large number of authors such as Braun and Hepp [3], Dawson [4] and Dobrushin [5], and most of the theory was achieved by Tanaka and collaborators [6-9] and of course Sznitman [10-12]. When considering that all particles (in our case, neurons) have the same, independent initial condition, they are mathematically proved using stochastic theory (the Wasserstein distance, large deviation techniques) that in the limit where the number of particles tends to infinity, any finite number of particles behaves independently of the other ones, and they all present the same probability distribution, which satisfies a nonlinear Markov equation. Finite-size fluctuations around the limit are derived in a general case in [10]. Most of these models use a standard hypothesis of global Lipschitz continuity and linear growth condition of the drift and diffusion coefficients of the diffusions, as well as the Lipschitz continuity of the interaction function. Extensions to discontinuous càdlàg processes including singular interactions (through a local time process) were developed in [11]. Problems involving singular interaction variables (e.g. nonsmooth functions) are also widely studied in the field, but are not relevant in our case.

In the present article, we apply this mathematical approach to the problem of interacting neurons arising in neuroscience. To this end, we extend the theory to encompass a wider class of models. This implies the use of locally (instead of globally) Lipschitz coefficients and of a Lyapunov-like growth condition replacing the customary linear growth assumption for some of the functions appearing in the equations. The contributions of this article are fourfold:

1. We derive, in a rigorous manner, the mean-field equations resulting from the interaction of infinitely many neurons in the case of widely accepted models of spiking neurons and synapses.

2. We prove a propagation of chaos property which shows that in the mean-field limit, the neurons become independent, in agreement with some recent experimental work [13] and with the idea that the brain processes information in a somewhat optimal way.

3. We show, numerically, that the mean-field limit is a good approximation of the mean activity of the network even for fairly small sizes of neuronal populations.

4. We suggest, numerically, that the changes in the dynamics of the mean-field limit when varying parameters can be understood by studying the mean-field Fokker-Planck equation.

We start by reviewing such models in the ‘Spiking conductance-based models’ section to motivate the present study. It is in the ‘Mean-field equations for conductance-based models’ section that we provide the limit equations describing the behaviors of an infinite number of interacting neurons and state and prove the existence and uniqueness of solutions in the case of conductance-based models. The detailed proof of the second main theorem, that of the convergence of the network equations to the mean-field limit, is given in the Appendix. In the ‘Numerical simulations’ section, we begin to address the difficult problem of the numerical simulation of the mean-field equations and show some results indicating that they may be an efficient way of representing the mean activity of a finite-size network as well as to study the changes in the dynamics when varying biological parameters. The final ‘Discussion and conclusion’ section focuses on the conclusions of our mathematical and numerical results and raises some important questions for future work.

### 2 Spiking conductance-based models

This section sets the stage for our results. We review in the ‘Hodgkin-Huxley model’ section the Hodgkin-Huxley model equations in the case where both the membrane potential and the ion channel equations include noise. We then proceed in the ‘The FitzHugh-Nagumo model’ section with the FitzHugh-Nagumo equations in the case where the membrane potential equation includes noise. We next discuss in the ‘Models of synapses and maximum conductances’ section the connectivity models of networks of such neurons, starting with the synapses, electrical and chemical, and finishing with several stochastic models of the synaptic weights. In the ‘Putting everything together’ section, we write the network equations in the various cases considered in the previous section and express them in a general abstract mathematical form that is the one used for stating and proving the results about the mean-field limits in the ‘Mean-field equations for conductance-based models’ section. Before we jump into this, we conclude in the ‘Mean-field methods in computational neuroscience: a quick overview’ section with a brief overview of the mean-field methods popular in computational neuroscience.

From the mathematical point of view, each neuron is a complex system, whose dynamics is often described by a set of stochastic nonlinear differential equations. Such models aim at reproducing the biophysics of ion channels governing the membrane potential and therefore the spike emission. This is the case of the classical model of Hodgkin and Huxley [14] and of its reductions [15-17]. Simpler models use discontinuous processes mimicking the spike emission by modeling the membrane voltage and considering that spikes are emitted when it reaches a given threshold. These are called integrate-and-fire models [18,19] and will not be addressed here. The models of large networks we deal with here therefore consist of systems of coupled nonlinear diffusion processes.

#### 2.1 Hodgkin-Huxley model

One of the most important models in computational neuroscience is the Hodgkin-Huxley
model. Using pioneering experimental techniques of that time, Hodgkin and Huxley [14] determined that the activity of the giant squid axon is controlled by three major
currents: voltage-gated persistent

The basic electrical relation between the membrane potential and the currents is simply:

where

where
*n* (respectively, *m*) is the activation variable for
*h* is the inactivation variable for Na. These activation/deactivation variables, denoted
by
*V* accounts for the voltage-gating of the gate) and to close with rate

The functions
^{a}

where

In order to complete our stochastic Hodgkin-Huxley neuron model, we assume that the
external current

This is a stochastic version of the Hodgkin-Huxley model. The functions
*n*, *m* and *h* are bounded between 0 and 1; hence, the functions

To illustrate the model, we show in Figure 1 the time evolution of the three ion channel variables *n*, *m* and *h* as well as that of the membrane potential *V* for a constant input

**Fig. 1.** Solution of the noiseless Hodgkin-Huxley model. *Left*: time evolution of the three ion channel variables *n*, *m* and *h*. *Right*: corresponding time evolution of the membrane potential. Parameters are given in
the text.

**Fig. 2.** Noisy Hodgkin-Huxley model. *Left*: time evolution of the three ion channel variables *n*, *m* and *h*. *Right*: corresponding time evolution of the membrane potential. Parameters are given in
the text.

For the membrane potential, we have used
*χ* function (see Equation 1):

with

Because the Hodgkin-Huxley model is rather complicated and high-dimensional, many reductions have been proposed, in particular to two dimensions instead of four. These reduced models include the famous FitzHugh-Nagumo and Morris-Lecar models. These two models are two-dimensional approximations of the original Hodgkin-Huxley model based on quantitative observations of the time scale of the dynamics of each variable and identification of variables. Most reduced models still comply with the Lipschitz and linear growth conditions ensuring the existence and uniqueness of a solution, except for the FitzHugh-Nagumo model which we now introduce.

#### 2.2 The FitzHugh-Nagumo model

In order to reduce the dimension of the Hodgkin-Huxley model, FitzHugh [15,16,21] introduced a simplified two-dimensional model. The motivation was to isolate conceptually
essential mathematical features yielding excitation and transmission properties from
the analysis of the biophysics of sodium and potassium flows. Nagumo and collaborators [22] followed up with an electrical system reproducing the dynamics of this model and
studied its properties. The model consists of two equations, one governing a voltage-like
variable *V* having a cubic nonlinearity and a slower recovery variable *w*. It can be written as:

where
*V* which we choose, without loss of generality, to be
*a*,
*w*. As in the case of the Hodgkin-Huxley model, the current
*I*, and a stochastic white noise accounting for the randomness of the environment. The
stochastic FitzHugh-Nagumo equation is deduced from Equation 4 and reads:

Note that because the function
*g*lobally Lipschitz continuous (only locally), the well-posedness of the stochastic
differential equation (Equation 5) does not follow immediately from the standard theorem
which assumes the global Lipschitz continuity of the drift and diffusion coefficients.
This question is settled below by Proposition 1.

We show in Figure 3 the time evolution of the adaptation variable and the membrane potential in the case
where the input *I* is constant and equal to 0.7. The left-hand side of the figure shows the case with
no noise while the right-hand side shows the case where noise of intensity

**Fig. 3.** Time evolution of the membrane potential and the adaptation variable in the FitzHugh-Nagumo
model. *Left*: without noise. *Right*: with noise. See text.

The deterministic model has been solved with a Runge-Kutta method of order 4, while
the stochastic model, with the Euler-Maruyama scheme. In both cases, we have used
an integration time step

#### 2.3 Partial conclusion

We have reviewed two main models of space-clamped single neurons: the Hodgkin-Huxley and FitzHugh-Nagumo models. These models are stochastic, including various sources of noise: external and internal. The noise sources are supposed to be independent Brownian processes. We have shown that the resulting stochastic differential Equations 2 and 5 were well-posed. As pointed out above, this analysis extends to a large number of reduced versions of the Hodgkin-Huxley such as those that can be found in the book [17].

#### 2.4 Models of synapses and maximum conductances

We now study the situation in which several of these neurons are connected to one
another forming a network, which we will assume to be fully connected. Let *N* be the total number of neurons. These neurons belong to *P* populations, e.g. pyramidal cells or interneurons. If the index of a neuron is *i*,
*i* is the index of a postsynaptic neuron belonging to population
*j* is the index of a presynaptic neuron to neuron *i* belonging to population

#### 2.4.1 Chemical synapses

The principle of functioning of chemical synapses is based on the release of a neurotransmitter
in the presynaptic neuron synaptic button, which binds to specific receptors on the
postsynaptic cell. This process, similar to the currents described in the Hodgkin
and Huxley model, is governed by the value of the cell membrane potential. We use
the model described in [20,23], which features a quite realistic biophysical representation of the processes at
work in the spike transmission and is consistent with the previous formalism used
to describe the conductances of other ion channels. The model emulates the fact that
following the arrival of an action potential at the presynaptic terminal, a neurotransmitter
is released in the synaptic cleft and binds to the postsynaptic receptor with a first
order kinetic scheme. Let *j* be a presynaptic neuron to the postynaptic neuron *i*. The synaptic current induced by the synapse from *j* to *i* can be modelled as the product of a conductance

The synaptic reversal potentials
*j*:

The function

The positive constants
*j*, i.e.
*j*. Its expression is given by (see, e.g. [20]):

Destexhe et al. [23] give some typical values of the parameters

Because of the dynamics of ion channels and of their finite number, similar to the
channel noise models derived through the Langevin approximation in the Hodgkin-Huxley
model (Equation 2), we assume that the proportion of active channels is actually governed
by a stochastic differential equation with diffusion coefficient
*γ* of *j* of the form (Equation 1):

In detail, we have

Remember that the form of the diffusion term guarantees that the solutions to this
equation with appropriate initial conditions stay between 0 and 1. The Brownian motions

#### 2.4.2 Electrical synapses

The electrical synapse transmission is rapid and stereotyped and is mainly used to
send simple depolarizing signals for systems requiring the fastest possible response.
At the location of an electrical synapse, the separation between two neurons is very
small (≈3.5 nm). This narrow gap is bridged by the *gap junction channels*, specialized protein structures that conduct the flow of ionic current from the presynaptic
to the postsynaptic cell (see, e.g. [24]).

Electrical synapses thus work by allowing ionic current to flow passively through the gap junction pores from one neuron to another. The usual source of this current is the potential difference generated locally by the action potential. Without the need for receptors to recognize chemical messengers, signaling at electrical synapses is more rapid than that which occurs across chemical synapses, the predominant kind of junctions between neurons. The relative speed of electrical synapses also allows for many neurons to fire synchronously.

We model the current for this type of synapse as

where

#### 2.4.3 The maximum conductances

As shown in Equations 6, 7 and 10, we model the current going through the synapse
connecting neuron *j* to neuron *i* as being proportional to the maximum conductance

The simplest idea is to assume that the maximum conductances are independent diffusion
processes with mean

where the
*NP*-independent zero mean unit variance white noise processes derived from *NP*-independent standard Brownian motions

One way to alleviate this problem is to modify the dynamics (Equation 11) to a slightly more complicated one whose solutions do not change sign, such as for instance, the Cox-Ingersoll-Ross model [25] given by:

Note that the right-hand side only depends upon the population

This shows that if the initial condition

#### 2.5 Putting everything together

We are ready to write the equations of a network of Hodgkin-Huxley or FitzHugh-Nagumo
neurons and study their properties and their limit, if any, when the number of neurons
becomes large. The external current for neuron *i* has been modelled as the sum of a deterministic part and a stochastic part:

We will assume that the deterministic part is the same for all neurons in the same
population,
*N* Brownian motions
*N*-independent Brownian motions and independent of all the other Brownian motions defined
in the model. In other words,

We only cover the case of chemical synapses and leave it to the reader to derive the equations in the simpler case of gap junctions.

#### 2.5.1 Network of FitzHugh-Nagumo neurons

We assume that the parameters
*i* are only functions of the population

*Simple maximum conductance variation.* If we assume that the maximum conductances fluctuate according to Equation 11, the
state of the *i*th neuron in a fully connected network of FitzHugh-Nagumo neurons with chemical synapses
is determined by the variables
*N* stochastic differential equations:

*N*-independent Brownian processes that model noise in the process of transmitter release
into the synaptic clefts.

*Sign-preserving maximum conductance variation.* If we assume that the maximum conductances fluctuate according to Equation 12, the
situation is slightly more complicated. In effect, the state space of the neuron *i* has to be augmented by the *P* maximum conductances

which is a set of

#### 2.5.2 Network of Hodgkin-Huxley neurons

We provide a similar description in the case of the Hodgkin-Huxley neurons. We assume
that the functions

*Simple maximum conductance variation.* If we assume that the maximum conductances fluctuate according to Equation 11, the
state of the *i*th neuron in a fully connected network of Hodgkin-Huxley neurons with chemical synapses
is therefore determined by the variables
*N* stochastic differential equations:

*Sign-preserving maximum conductance variation.* If we assume that the maximum conductances fluctuate according to Equation 12, we
use the same idea as in the FitzHugh-Nagumo case of augmenting the state space of
each individual neuron and obtain the following set of

#### 2.5.3 Partial conclusion

Equations 14 to 17 have a quite similar structure. They are well-posed, i.e. given
any initial condition, and any time
*n*, *m* and *h*, which take values between 0 and 1, and the maximum conductances when one wants to
preserve their signs.

In order to prepare the grounds for the ‘Mean-field equations for conductance-based
models’ section, we explore a bit more the aforementioned common structure. Let us
first consider the case of the simple maximum conductance variations for the FitzHugh-Nagumo
network. Looking at Equation 14, we define the three-dimensional state vector of neuron
*i* to be

Let us next define

It appears that the intrinsic dynamics of the neuron *i* is conveniently described by the equation

We next define the functions

and the function

It appears that the full dynamics of the neuron *i*, corresponding to Equation 14, can be described compactly by

Let us now move to the case of the sign-preserving variation of the maximum conductances,
still for the FitzHugh-Nagumo neurons. The state of each neuron is now *P*+3-dimensional: we define

and the functions

It appears that the intrinsic dynamics of the neuron *i* isolated from the other neurons is conveniently described by the equation

Let us finally define the functions

It appears that the full dynamics of the neuron *i*, corresponding to Equation 15 can be described compactly by

We let the reader apply the same machinery to the network of Hodgkin-Huxley neurons.

Let us note *d* as the positive integer equal to the dimension of the state space in Equation 18
(

(H1) *Locally Lipschitz dynamics*: For all
*U*:

(H2) *Locally Lipschitz interactions*: For all

(H3) *Linear growth of the interactions*: There exists a

(H4) *Monotone growth of the dynamics*: We assume that

These assumptions are central to the proofs of Theorems 2 and 4.

They imply the following proposition stating that the system of stochastic differential equations (Equation 19) is well-posed:

**Proposition 1***Let*
*be a fixed time*. *If*
*on*
*for*
*Equations *18 *and* 19 *together with an initial condition*
*of square*-*integrable random variables*, *have a unique strong solution which belongs to*

*Proof* The proof uses Theorem 3.5 in chapter 2 in [26] whose conditions are easily shown to follow from hypotheses 2.5.3 to (H2). □

The case

We are interested in the behavior of the solutions of these equations as the number of neurons tends to infinity. This problem has been long-standing in neuroscience, arousing the interest of many researchers in different domains. We discuss the different approaches developed in the field in the next subsection.

#### 2.6 Mean-field methods in computational neuroscience: a quick overview

Obtaining the equations of evolution of the effective mean-field from microscopic dynamics is a very complex problem. Many approximate solutions have been provided, mostly based on the statistical physics literature.

Many models describing the emergent behavior arising from the interaction of neurons in large-scale networks have relied on continuum limits ever since the seminal work of Amari, and Wilson and Cowan [27-30]. Such models represent the activity of the network by macroscopic variables, e.g. the population-averaged firing rate, which are generally assumed to be deterministic. When the spatial dimension is not taken into account in the equations, they are referred to as neural masses, otherwise as neural fields. The equations that relate these variables are ordinary differential equations for neural masses and integrodifferential equations for neural fields. In the second case, they fall in a category studied in [31] or can be seen as ordinary differential equations defined on specific functional spaces [32]. Many analytical and numerical results have been derived from these equations and related to cortical phenomena, for instance, for the problem of spatio-temporal pattern formation in spatially extended models (see, e.g. [33-36]). The use of bifurcation theory has also proven to be quite powerful [37,38]. Despite all its qualities, this approach implicitly makes the assumption that the effect of noise vanishes at the mesoscopic and macroscopic scales and hence that the behavior of such populations of neurons is deterministic.

A different approach has been to study regimes where the activity is uncorrelated. A number of computational studies on the integrate-and-fire neuron showed that under certain conditions, neurons in large assemblies end up firing asynchronously, producing null correlations [39-41]. In these regimes, the correlations in the firing activity decrease towards zero in the limit where the number of neurons tends to infinity. The emergent global activity of the population in this limit is deterministic and evolves according to a mean-field firing rate equation. However, according to the theory, these states only exist in the limit where the number of neurons is infinite, thereby raising the question of how the finiteness of the number of neurons impacts the existence and behavior of asynchronous states. The study of finite-size effects for asynchronous states is generally not reduced to the study of mean firing rates and can include higher order moments of firing activity [42-44]. In order to go beyond asynchronous states and take into account the stochastic nature of the firing and understand how this activity scales as the network size increases, different approaches have been developed, such as the population density method and related approaches [45]. Most of these approaches involve expansions in terms of the moments of the corresponding random variables, and the moment hierarchy needs to be truncated which is not a simple task that can raise a number of difficult technical issues (see, e.g. [46]).

However, increasingly many researchers now believe that the different intrinsic or extrinsic noise sources are part of the neuronal signal, and rather than being a pure disturbing effect related to the intrinsically noisy biological substrate of the neural system, they suggest that noise conveys information that can be an important principle of brain function [47]. At the level of a single cell, various studies have shown that the firing statistics are highly stochastic with probability distributions close to the Poisson distributions [48], and several computational studies confirmed the stochastic nature of single-cell firings [49-51]. How the variability at the single-neuron level affects the dynamics of cortical networks is less well established. Theoretically, the interaction of a large number of neurons that fire stochastic spike trains can naturally produce correlations in the firing activity of the population. For instance, power laws in the scaling of avalanche-size distributions has been studied both via models and experiments [52-55]. In these regimes, the randomness plays a central role.

In order to study the effect of the stochastic nature of the firing in large networks,
many authors strived to introduce randomness in a tractable form. Some of the models
proposed in the area are based on the definition of a Markov chain governing the firing
dynamics of the neurons in the network, where the transition probability satisfies
a differential equation, the *master equation*. Seminal works of the application of such modeling for neuroscience date back to
the early 1990s and have been recently developed by several authors [43,56-59]. Most of these approaches are proved correct in some parameter regions using statistical
physics tools such as path integrals and Van-Kampen expansions, and their analysis
often involve a moment expansion and truncation. Using a different approach, a static
mean-field study of multi-population network activity was developed by Treves in [60]. This author did not consider external inputs but incorporated dynamical synaptic
currents and adaptation effects. His analysis was completed in [39], where the authors proved, using a Fokker-Planck formalism, the stability of an asynchronous
state in the network. Later on, Gerstner in [61] built a new approach to characterize the mean-field dynamics for the spike response
model, via the introduction of suitable kernels propagating the collective activity
of a neural population in time. Another approach is based on the use of large deviation
techniques to study large networks of neurons [62]. This approach is inspired by the work on spin-glass dynamics, e.g. [63]. It takes into account the randomness of the maximum conductances and the noise at
various levels. The individual neuron models are rate models, hence already mean-field
models. The mean-field equations are not rigorously derived from the network equations
in the limit of an infinite number of neurons, but they are shown to have a unique,
non-Markov solution, i.e. with infinite memory, for each initial condition.

Brunel and Hakim considered a network of integrate-and-fire neurons connected with constant maximum conductances [41]. In the case of sparse connectivity, stationarity, and in a regime where individual neurons emit spikes at a low rate, they were able to analytically study the dynamics of the network and to show that it exhibits a sharp transition between a stationary regime and a regime of fast collective oscillations weakly synchronized. Their approach was based on a perturbative analysis of the Fokker-Planck equation. A similar formalism was used in [44] which, when complemented with self-consistency equations, resulted in the dynamical description of the mean-field equations of the network and was extended to a multi population network. Finally, Chizhov and Graham [64] have recently proposed a new method based on a population density approach allowing to characterize the mesoscopic behavior of neuron populations in conductance-based models.

Let us finish this very short and incomplete survey by mentioning the work of Sompolinsky and colleagues. Assuming a linear intrinsic dynamics for the individual neurons described by a rate model and random centered maximum conductances for the connections, they showed [65,66] that the system undergoes a phase transition between two different stationary regimes: a ‘trivial’ regime where the system has a unique null and uncorrelated solution, and a ‘chaotic’ regime in which the firing rate converges towards a non-zero value and correlations stabilize on a specific curve which they were able to characterize.

All these approaches have in common that they are not based on the most widely accepted microscopic dynamics (such as the ones represented by the Hodgkin-Huxley equations or some of their simplifications) and/or involve approximations or moment closures. Our approach is distinct in that it aims at deriving rigorously and without approximations the mean-field equations of populations of neurons whose individual neurons are described by biological, if not correct at least plausible, representations. The price to pay is the complexity of the resulting mean-field equations. The specific study of their solutions is therefore a crucial step, which will be developed in forthcoming papers.

### 3 Mean-field equations for conductance-based models

In this section, we give a general formulation of the neural network models introduced
in the previous section and use it in a probabilistic framework to address the problem
of the asymptotic behavior of the networks, as the number of neurons *N* goes to infinity. In other words, we derive the limit in law of *N*-interacting neurons, each of which satisfying a nonlinear stochastic differential
equation of the type described in the ‘Spiking conductance-based models’ section.
In the remainder of this section, we work in a complete probability space

#### 3.1 Setting of the problem

We recall that the neurons in the network fall into different populations *P*. The populations differ through the intrinsic properties of their neurons and the
input they receive. We assume that the number of neurons in each population
*α* is nontrivial, i.e.
*N* goes to infinity^{b}.

We use the notations introduced in the ‘Partial conclusion’ section, and the reader should refer to this section to give a concrete meaning to the rather abstract (but required by the mathematics) setting that we now establish.

Each neuron *i* in population *α* is described by a state vector noted as
*i* in population *α*, the dynamics of the *d*-dimensional process

Moreover, we assume, as it is the case for all the models described in the ‘Spiking conductance-based models’ section, that the solutions of this stochastic differential equation exist for all time.

When included in the network, these processes interact with those of all the other
neurons through a set of continuous functions that only depend on the population
*i* belongs to and the populations *γ* of the presynaptic neurons. These functions,
*N* goes to infinity).

As discussed in the ‘Spiking conductance-based models’ section, due to the stochastic
nature of ionic currents and the noise effects linked with the discrete nature of
charge carriers, the maximum conductances are perturbed dynamically through the
*δ* that were previously introduced. The interaction between the neurons and the noise
term is represented by the function

In order to introduce the stochastic current and stochastic maximum conductances,
we define two independent sequences of independent *m*- and *δ*-dimensional Brownian motions noted as

The resulting equation for the *i*th neuron, including the noisy interactions, reads:

Note that this implies that

These equations are similar to the equations studied in another context by a number of mathematicians, among which are McKean, Tanaka and Sznitman (see the ‘Introduction’ section), in that they involve a very large number of particles (here, particles are neurons) in interaction. Motivated by the study of the McKean-Vlasov equations, these authors studied special cases of equations (Equation 21). This theory, referred to as the kinetic theory, is chiefly interested in the study of the thermodynamics questions. They show the property that in the limit where the number of particles tends to infinity, provided that the initial state of each particle is drawn independently from the same law, each particle behaves independently and has the same law, which is given by an implicit stochastic equation. They also evaluate the fluctuations around this limit under diverse conditions [1,2,6,7,9-11]. Some extensions to biological problems where the drift term is not globally Lipschitz but satisfies the monotone growth condition (Equation 20) were studied in [67]. This is the approach we undertake here.

#### 3.2 Convergence of the network equations to the mean-field equations and properties of those equations

We now show that the same type of phenomena that were predicted for systems of interacting
particles happen in networks of neurons. In detail, we prove that in the limit of
large populations, the network displays the property of propagation of chaos. This
means that any finite number of diffusion processes become independent, and all neurons
belonging to a given population *α* have asymptotically the same probability distribution, which is the solution of the
following mean-field equation:

where

In these equations,
*m*-dimensional Brownian motions. Let us point out the fact that the right-hand side
of Equations 22 and 23 depends on the law of the solution; this fact is sometimes
referred to as ‘the process
*p* of the solution. This equation which we use in the ‘Numerical simulations’ section
can be easily derived from Equation 22. Let
*t* of the solution

where the

with

The *P* equations (Equation 24) yield the probability densities of the solutions

We now spend some time on notations in order to obtain a somewhat more compact form
of Equation 22. We define
*dP*-dimensional process
*f*, *g*, *b* and *β* as the concatenations of functions

We obtain the equivalent compact mean-field equation:

Equations 22 and 24 are implicit equations on the law of

We now state the main theoretical results of the paper as two theorems. The first theorem is about the well-posedness of the mean-field equation (Equation 22). The second is about the convergence of the solutions of the network equations to those of the mean-field equations. Since the proof of the second theorem involves similar ideas to those used in the proof of the first, it is given in the Appendix.

**Theorem 2***Under assumptions* (H1) *to* (H4), *there exists a unique solution to the mean*-*field equation* (*Equation *22) *on*
*for any*

Let us denote by
*P* (respectively,
*m* (respectively *δ*)-dimensional, adapted standard Brownian motions on

We have introduced in the previous formula the process

The following lemma is useful to prove the theorem:

**Lemma 3***Let*
*be a square*-*integrable random variable*. *Let**X**be a solution of the mean*-*field equation* (*Equation *22) *with initial condition*
*Under assumptions *(H3) *and* (H4), *there exists a constant*
*depending on the parameters of the system and on the horizon**T*, *such that*:

*Proof* Using the Itô formula for

where

This expression involves the term

It also involves the term

Finally, we obtain

Using Gronwall’s inequality, we deduce the

This lemma puts us in a position to prove the existence and uniqueness theorem:

*Proof* We start by showing the existence of solutions and then prove the uniqueness property.
We recall that by the application of Lemma 3, the solutions will all have bounded
second-order moment.

*Existence*. Let
*X* and *Z* i.i.d.’ below. We stop the processes at the time
*U*. For convenience, we will make an abuse of notation in the proof and denote
*U* centered at the origin in

Using the notations introduced for Equation 25, we decompose the difference

and find an upperbound for

and treat each term separately. The upperbounds for the first two terms are obtained using the Cauchy-Schwartz inequality, those of the last two terms using the Burkholder-Davis-Gundy martingale moment inequality.

The term

Taking the sup of both sides of the last inequality, we obtain

from which follows the fact that

The term
*X* and *Z* are independent with the same law:

Taking the sup of both sides of the last inequality, we obtain

from which follows the fact that

The term

The term

Putting all of these together, we get:

From the relation

and

and this upper bound is the term of a convergent series. The Borel-Cantelli lemma
stems that for almost any
*ω* denotes an element of the probability space Ω) such that

and hence

It follows that with probability 1, the partial sums:

are uniformly (in
*t*, the sequence

It is easy to show using routine methods that

To complete the proof, we use a standard truncation property. This method replaces
the function *f* by the truncated function:

and similarly for *g*. The functions

Let us now define the stopping time as

It is easy to show that

implying that the sequence of stopping times

and letting

*Uniqueness*. Assume that *X* and *Y* are two solutions of the mean-field equations (Equation 22). From Lemma 3, we know
that both solutions are in

which, by Gronwall’s theorem, directly implies that

which ends the proof. □

We have proved the well-posedness of the mean-field equations. It remains to show that the solutions to the network equations converge to the solutions of the mean-field equations. This is what is achieved in the next theorem.

**Theorem 4***Under assumptions* (H1) *to* (H4), *the following holds true*:

• Convergence^{c}: *For each neuron**i**of population**α*, *the law of the multidimensional process*
*converges towards the law of the solution of the mean*-*field equation related to population**α*, *namely*

• Propagation of chaos: *For any*
*and any**k*-*tuple*
*the law of the process*
*converges towards*^{d}
*i*.*e*. *the asymptotic processes have the law of the solution of the mean*-*field equations and are all independent*.

This theorem has important implications in neuroscience that we discuss in the ‘Discussion and conclusion’ section. Its proof is given in the Appendix.

### 4 Numerical simulations

At this point, we have provided a compact description of the activity of the network when the number of neurons tends to infinity. However, the structure of the solutions of these equations is complicated to understand from the implicit mean-field equations (Equation 22) and of their variants (such as the McKean-Vlasov-Fokker-Planck equations (Equation 24)). In this section, we present some classical ways to numerically approximate the solutions to these equations and give some indications about the rate of convergence and the accuracy of the simulation. These numerical schemes allow us to compute and visualize the solutions. We then compare the results of the two schemes for a network of FitzHugh-Nagumo neurons belonging to a single population and show their good agreement.

The main difficulty one faces when developing numerical schemes for Equations 22 and 24 is that they are non-local. By this, we mean that in the case of the McKean-Vlasov equations, they contain the expectation of a certain function under the law of the solution to the equations (see Equation 22). In the case of the corresponding Fokker-Planck equation, it contains integrals of the probability density functions which is a solution to the equation (see Equation 24).

#### 4.1 Numerical simulations of the McKean-Vlasov equations

The fact that the McKean-Vlasov equations involve an expectation of a certain function
under the law of the solution of the equation makes them particularly hard to simulate
directly. One is often reduced to use Monte Carlo simulations to compute this expectation,
which amounts to simulating the solution of the network equations themselves (see [68]). This is the method we used. In its simplest fashion, it consists of a Monte Carlo
simulation where one numerically solves the *N* network equations (Equation 21) with the classical Euler-Maruyama method a number
of times with different initial conditions, and averages the trajectories of the solutions
over the number of simulations.

In detail, let
*M* times a *P*-population discrete-time process
*i* in population *α*:

where
*d*- and *δ*-dimensional standard normal random variables. The initial conditions
*N* of the population is large enough, Theorem 4 states that the law, noted as
*t* falls into that particular bin. The resulting histogram can then be compared to the
solution of the McKean-Vlasov-Fokker-Planck equation (Equation 24) corresponding to
population *α* whose numerical solution is described next.

The mean square error between the solution of the numerical recursion (Equation 30)
*N*,

#### 4.2 Numerical simulations of the McKean-Vlasov-Fokker-Planck equation

For solving the McKean-Vlasov-Fokker-Planck equation (Equation 24), we have used the
*method of lines* [72,73]. Its basic idea is to discretize the phase space and to keep the time continuous.
In this way, the values
*α* at each sample point *X* of the phase space are the solutions of *P* ODEs where the independent variable is the time. Each sample point in the phase space
generates *P* ODEs, resulting in a system of coupled ODEs. The solutions to this system yield the
values of the probability density functions
*f* to be integrated between the values

where Δ*x* is the integration step, and

The discretization of the derivatives with respect to the phase space parameters is done through the following fourth-order central difference scheme:

for the first-order derivatives, and

for the second-order derivatives (see [75]).

Finally, we have used a Runge-Kutta method of order 2 (RK2) for the numerical integration
of the resulting system of ODEs. This method is of the explicit kind for ordinary
differential equations, and it is described by the following *Butcher tableau*:

#### 4.3 Comparison between the solutions to the network and the mean-field equations

We illustrate these ideas with the example of a network of 100 FitzHugh-Nagumo neurons
belonging to one, excitatory, population. We also use chemical synapses with the variation
of the weights described by (Equation 11). We choose a finite volume, outside of which
we assume that the probability density function (p.d.f.) is zero. We then discretize
this volume with

where
*V*, Δ*w* and Δ*y* are the quantization steps in each dimension of the phase space. For the simulation
of the McKean-Vlasov-Fokker-Planck equation, instead, we use Dirichlet boundary conditions
and assume the probability and its partial derivatives to be 0 on the boundary and
outside the volume.

In general, the total number of coupled ODEs that we have to solve for the McKean-Vlasov-Fokker-Planck
equation with the method of lines is the product
*t* used in the RK2 scheme. This can strongly impact the efficiency of the numerical
method (see the ‘Numerical simulations with GPUs’ section).

In the simulations shown in the left-hand parts of Figures 4 and 5, we have used one population of 100 excitatory FitzHugh-Nagumo neurons connected
with chemical synapses. We performed

**Fig. 4.** Joint probability distribution.
*left*) compared with the solution of the McKean-Vlasov-Fokker-Planck equation (Equation 24)
(*right*), sampled at four times

**Fig. 5.** Joint probability distribution.
*left*) compared with the solution of the McKean-Vlasov-Fokker-Planck equation (Equation 24)
(*right*), sampled at four times

The parameters are given in the first column of Table 1. In this table, the parameter

**Table 1.** Parameters used in the simulations of the neural network and for solving the McKean-Vlasov-Fokker-Planck
equation

The parameters for the noisy model of maximum conductances of Equation 11 are shown
in the fourth column of the table. For these values of
*χ* function (Equation 3). The solutions are computed over an interval of

The marginals estimated from the trajectories of the network solutions are then compared to those obtained from the numerical solution of the McKean-Vlasov-Fokker-Planck equation (see Figures 4 and 5 right), using the method of lines explained above and starting from the same initial conditions (Equation 31) as the neural network.

We have used the value

**Fig. 6.** Projection of 100 trajectories in the
*top left*),
*top right*) and