Abstract
We derive the meanfield 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 HodgkinHuxley model or by one of its simplified version, the FitzHughNagumo 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 meanfield 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 McKeanVlasov equations or nonlocal partial differential equations resembling the McKeanVlasovFokkerPlanck equations. We prove the wellposedness of the McKeanVlasov equations, i.e. the existence and uniqueness of a solution. We also show the results of some numerical experiments that indicate that the meanfield equations are a good representation of the mean activity of a finite size network, even for modest sizes. These experiments also indicate that the McKeanVlasovFokkerPlanck equations may be a good way to understand the meanfield dynamics through, e.g. a bifurcation analysis.
Mathematics Subject Classification (2000): 60F99, 60B10, 92B20, 82C32, 82C80, 35Q80.
Keywords:
meanfield limits; propagation of chaos; stochastic differential equations; McKeanVlasov equations; FokkerPlanck equations; neural networks; neural assemblies; HodgkinHuxley neurons; FitzHughNagumo neurons1 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 largescale neural assemblies has been a great endeavor in the past decades. One of the main interests of largescale 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 barshaped 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 longstanding 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 [69] and of course Sznitman [1012]. 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. Finitesize 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 Lyapunovlike 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 meanfield 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 meanfield 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 meanfield 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 meanfield limit when varying parameters can be understood by studying the meanfield FokkerPlanck equation.
We start by reviewing such models in the ‘Spiking conductancebased models’ section to motivate the present study. It is in the ‘Meanfield equations for conductancebased 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 conductancebased models. The detailed proof of the second main theorem, that of the convergence of the network equations to the meanfield limit, is given in the Appendix. In the ‘Numerical simulations’ section, we begin to address the difficult problem of the numerical simulation of the meanfield equations and show some results indicating that they may be an efficient way of representing the mean activity of a finitesize 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 conductancebased models
This section sets the stage for our results. We review in the ‘HodgkinHuxley model’ section the HodgkinHuxley model equations in the case where both the membrane potential and the ion channel equations include noise. We then proceed in the ‘The FitzHughNagumo model’ section with the FitzHughNagumo 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 meanfield limits in the ‘Meanfield equations for conductancebased models’ section. Before we jump into this, we conclude in the ‘Meanfield methods in computational neuroscience: a quick overview’ section with a brief overview of the meanfield 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 [1517]. 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 integrateandfire 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 HodgkinHuxley model
One of the most important models in computational neuroscience is the HodgkinHuxley 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: voltagegated persistent current with four activation gates, voltagegated transient current with three activation gates and one inactivation gate, and Ohmic leak current, , which is carried mostly by chloride ions (). In this paper, we only use the spaceclamped HodgkinHuxley model which we slightly generalize to a stochastic setting in order to better take into account the variability of the parameters. The advantages of this model are numerous, and one of the most prominent aspects in its favor is its correspondence with the most widely accepted formalism to describe the dynamics of the nerve cell membrane. A very extensive literature can also be found about the mathematical properties of this system, and it is now quite well understood.
The basic electrical relation between the membrane potential and the currents is simply:
where is an external current. The detailed expressions for , and can be found in several textbooks, e.g. [17,20]:
where (respectively, ) is the maximum conductance of the potassium (respectively, the sodium) channel; is the conductance of the Ohmic channel; and n (respectively, m) is the activation variable for (respectively, for Na). There are four (respectively, three) activation gates for the (respectively, the Na) current which accounts for the power 4 (respectively, 3) in the expression of (respectively ). h is the inactivation variable for Na. These activation/deactivation variables, denoted by in what follows, represent a proportion (they vary between 0 and 1) of open gates. The proportions of open channels are given by the functions and . The proportions of open gates can be computed through a Markov chain modeling assuming the gates to open with rate (the dependence in V accounts for the voltagegating of the gate) and to close with rate . These processes can be shown to converge, under standard assumptions, towards the following ordinary differential equations:
The functions and are smooth functions whose exact values can be found in several textbooks such as the ones cited above. Note that half of these six functions are unbounded when the voltage goes to −∞, being of the form , with and as two positive constants. Since these functions have been fitted to experimental data corresponding to values of the membrane potential between roughly −100 and 100 mVs, it is clear that extremely large in magnitude and negative values of this variable do not have any physiological meaning. We can therefore safely, smoothly perturb these functions so that they are upperbounded by some large (but finite) positive number for these values of the membrane potential. Hence, the functions and are bounded and Lipschitz continuous for . A more precise model taking into account the finite number of channels through the Langevin approximation results in the stochastic differential equation^{a}
where and are independent standard Brownian motions. is a function that vanishes outside . This guarantees that the solution remains a proportion, i.e. lies between 0 and 1 for all times. We define
In order to complete our stochastic HodgkinHuxley neuron model, we assume that the external current is the sum of a deterministic part, noted as , and a stochastic part, a white noise with variance built from a standard Brownian motion independent of and . Considering the current produced by the income of ion through these channels, we end up with the following system of stochastic differential equations:
This is a stochastic version of the HodgkinHuxley model. The functions and are bounded and Lipschitz continuous (see discussion above). The functions n, m and h are bounded between 0 and 1; hence, the functions and are Lipschitz continuous.
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 . The system of ordinary differential equations (ODEs) has been solved using a RungeKutta scheme of order 4 with an integration time step . In Figure 2, we show the same time evolution when noise is added to the channel variables and the membrane potential.
Fig. 1. Solution of the noiseless HodgkinHuxley 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 HodgkinHuxley 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 (see Equation 2), while for the noise in the ion channels, we have used the following χ function (see Equation 1):
with and for all the ion channels. The system of SDEs has been integrated using the EulerMaruyama scheme with .
Because the HodgkinHuxley model is rather complicated and highdimensional, many reductions have been proposed, in particular to two dimensions instead of four. These reduced models include the famous FitzHughNagumo and MorrisLecar models. These two models are twodimensional approximations of the original HodgkinHuxley 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 FitzHughNagumo model which we now introduce.
2.2 The FitzHughNagumo model
In order to reduce the dimension of the HodgkinHuxley model, FitzHugh [15,16,21] introduced a simplified twodimensional 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 voltagelike variable V having a cubic nonlinearity and a slower recovery variable w. It can be written as:
where is a cubic polynomial in V which we choose, without loss of generality, to be . The parameter models the input current the neuron receives; the parameters a, and describe the kinetics of the recovery variable w. As in the case of the HodgkinHuxley model, the current is assumed to be the sum of a deterministic part, noted I, and a stochastic white noise accounting for the randomness of the environment. The stochastic FitzHughNagumo equation is deduced from Equation 4 and reads:
Note that because the function is not globally Lipschitz continuous (only locally), the wellposedness 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 lefthand side of the figure shows the case with no noise while the righthand side shows the case where noise of intensity (see Equation 5) has been added.
Fig. 3. Time evolution of the membrane potential and the adaptation variable in the FitzHughNagumo model. Left: without noise. Right: with noise. See text.
The deterministic model has been solved with a RungeKutta method of order 4, while the stochastic model, with the EulerMaruyama scheme. In both cases, we have used an integration time step .
2.3 Partial conclusion
We have reviewed two main models of spaceclamped single neurons: the HodgkinHuxley and FitzHughNagumo 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 wellposed. As pointed out above, this analysis extends to a large number of reduced versions of the HodgkinHuxley 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, , we note , as the population it belongs to. We note as the number of neurons in population . Since we want to be as close to biology as possible while keeping the possibility of a mathematical analysis of the resulting model, we consider two types of simplified, but realistic, synapses: chemical and electrical or gap junctions. The following material concerning synapses is standard and can be found in textbooks [20]. The new, and we think important, twist is to add noise to our models. To unify notations, in what follows, i is the index of a postsynaptic neuron belonging to population , and 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 with a voltage difference:
The synaptic reversal potentials are approximately constant within each population: . The conductance is the product of the maximum conductance with a function that denotes the fraction of open channels and depends only upon the presynaptic neuron j:
The function is often modelled [20] as satisfying the following ordinary differential equation:
The positive constants and characterize the rise and decay rates, respectively, of the synaptic conductance. Their values depend only on the population of the presynaptic neuron j, i.e. and , but may vary significantly from one population to the next. For example, gammaaminobutyric acid synapses are slow to activate and slow to turn off while the reverse is true for and AMPA synapses [20]. denotes the concentration of the transmitter released into the synaptic cleft by a presynaptic spike. We assume that the function is sigmoidal and that its exact form depends only upon the population of the neuron j. Its expression is given by (see, e.g. [20]):
Destexhe et al. [23] give some typical values of the parameters , and .
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 HodgkinHuxley model (Equation 2), we assume that the proportion of active channels is actually governed by a stochastic differential equation with diffusion coefficient depending only on the population γ 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 are assumed to be independent from one neuron to the next.
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 is the maximum conductance.
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 . Because the synaptic transmission through a synapse is affected by the nature of the environment, the maximum conductances are affected by dynamical random variations (we do not take into account such phenomena as plasticity). What kind of models can we consider for these random variations?
The simplest idea is to assume that the maximum conductances are independent diffusion processes with mean and standard deviation , i.e. that depend only on the populations. The quantities , being conductances, are positive. We write the following equation:
where the , , , are NPindependent zero mean unit variance white noise processes derived from NPindependent standard Brownian motions , i.e. , which we also assume to be independent of all the previously defined Brownian motions. The main advantage of this dynamics is its simplicity. Its main disadvantage is that if we increase the noise level , the probability that becomes negative increases also: this would result in a negative conductance!
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 CoxIngersollRoss model [25] given by:
Note that the righthand side only depends upon the population . Let be the initial condition, it is known [25] that
This shows that if the initial condition is equal to the mean , the mean of the process is constant over time and equal to . Otherwise, if the initial condition is of the same sign as , i.e. positive, then the long term mean is and the process is guaranteed not to touch 0 if the condition holds [25]. Note that the long term variance is .
2.5 Putting everything together
We are ready to write the equations of a network of HodgkinHuxley or FitzHughNagumo 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, , and that the same is true for the variance, . We further assume that the N Brownian motions are Nindependent 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 FitzHughNagumo neurons
We assume that the parameters , and in Equation 5 of the adaptation variable of neuron 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 ith neuron in a fully connected network of FitzHughNagumo neurons with chemical synapses is determined by the variables that satisfy the following set of 3N stochastic differential equations:
is given by Equation 8; , by Equation 9; and , , are Nindependent Brownian processes that model noise in the process of transmitter release into the synaptic clefts.
Signpreserving 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 , . We obtain
which is a set of stochastic differential equations.
2.5.2 Network of HodgkinHuxley neurons
We provide a similar description in the case of the HodgkinHuxley neurons. We assume that the functions and , , that appear in Equation 2 only depend upon .
Simple maximum conductance variation. If we assume that the maximum conductances fluctuate according to Equation 11, the state of the ith neuron in a fully connected network of HodgkinHuxley neurons with chemical synapses is therefore determined by the variables that satisfy the following set of 5N stochastic differential equations:
Signpreserving maximum conductance variation. If we assume that the maximum conductances fluctuate according to Equation 12, we use the same idea as in the FitzHughNagumo case of augmenting the state space of each individual neuron and obtain the following set of stochastic differential equations:
2.5.3 Partial conclusion
Equations 14 to 17 have a quite similar structure. They are wellposed, i.e. given any initial condition, and any time , they have a unique solution on which is squareintegrable. A little bit of care has to be taken when choosing these initial conditions for some of the parameters, i.e. 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 ‘Meanfield equations for conductancebased 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 FitzHughNagumo network. Looking at Equation 14, we define the threedimensional state vector of neuron i to be . Let us now define , , by
It appears that the intrinsic dynamics of the neuron i is conveniently described by the equation
We next define the functions , for , by
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 signpreserving variation of the maximum conductances, still for the FitzHughNagumo neurons. The state of each neuron is now P+3dimensional: we define . We then define the functions , , by
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 , for , by
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 HodgkinHuxley neurons.
Let us note d as the positive integer equal to the dimension of the state space in Equation 18 () or 19 () or in the corresponding cases for the HodgkinHuxley model ( and ). The reader will easily check that the following four assumptions hold for both models:
(H1) Locally Lipschitz dynamics: For all , the functions and are uniformly locally Lipschitz continuous with respect to the second variable. In detail, for all , there exists independent of such that for all , the ball of of radius U:
(H2) Locally Lipschitz interactions: For all , the functions and are locally Lipschitz continuous. In detail, for all , there exists such that for all , we have:
(H3) Linear growth of the interactions: There exists a such that
(H4) Monotone growth of the dynamics: We assume that and satisfy the following monotonous condition for all :
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 wellposed:
Proposition 1Letbe a fixed time. Ifon, for, Equations 18 and 19 together with an initial condition, of squareintegrable 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 implies that Equations 2 and 5, describing the stochastic FitzHughNagumo and HodgkinHuxley neurons, are wellposed.
We are interested in the behavior of the solutions of these equations as the number of neurons tends to infinity. This problem has been longstanding 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 Meanfield methods in computational neuroscience: a quick overview
Obtaining the equations of evolution of the effective meanfield 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 largescale networks have relied on continuum limits ever since the seminal work of Amari, and Wilson and Cowan [2730]. Such models represent the activity of the network by macroscopic variables, e.g. the populationaveraged 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 spatiotemporal pattern formation in spatially extended models (see, e.g. [3336]). 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 integrateandfire neuron showed that under certain conditions, neurons in large assemblies end up firing asynchronously, producing null correlations [3941]. 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 meanfield 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 finitesize effects for asynchronous states is generally not reduced to the study of mean firing rates and can include higher order moments of firing activity [4244]. 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 singlecell firings [4951]. How the variability at the singleneuron 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 avalanchesize distributions has been studied both via models and experiments [5255]. 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,5659]. Most of these approaches are proved correct in some parameter regions using statistical physics tools such as path integrals and VanKampen expansions, and their analysis often involve a moment expansion and truncation. Using a different approach, a static meanfield study of multipopulation 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 FokkerPlanck formalism, the stability of an asynchronous state in the network. Later on, Gerstner in [61] built a new approach to characterize the meanfield 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 spinglass 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 meanfield models. The meanfield 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, nonMarkov solution, i.e. with infinite memory, for each initial condition.
Brunel and Hakim considered a network of integrateandfire 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 FokkerPlanck equation. A similar formalism was used in [44] which, when complemented with selfconsistency equations, resulted in the dynamical description of the meanfield 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 conductancebased 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 nonzero 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 HodgkinHuxley 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 meanfield 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 meanfield equations. The specific study of their solutions is therefore a crucial step, which will be developed in forthcoming papers.
3 Meanfield equations for conductancebased 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 Ninteracting neurons, each of which satisfying a nonlinear stochastic differential equation of the type described in the ‘Spiking conductancebased models’ section. In the remainder of this section, we work in a complete probability space satisfying the usual conditions and endow with a filtration .
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 , denoted by , increases as the network size increases and moreover that the asymptotic proportion of neurons in population α is nontrivial, i.e. as 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 in and has an intrinsic dynamics governed by a drift function and a diffusion matrix assumed uniformly locally Lipschitz continuous with respect to the second variable. For a neuron i in population α, the dynamics of the ddimensional process governing the evolution of the membrane potential and additional variables (adaptation, ionic concentrations), when there is no interaction, is governed by the equation:
Moreover, we assume, as it is the case for all the models described in the ‘Spiking conductancebased 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 , the neuron i belongs to and the populations γ of the presynaptic neurons. These functions, , are scaled by the coefficients , so the maximal interaction is independent of the size of the network (in particular, neither diverging nor vanishing as N goes to infinity).
As discussed in the ‘Spiking conductancebased 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 independent Brownian motions of dimension δ 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 and which are adapted to the filtration .
The resulting equation for the ith neuron, including the noisy interactions, reads:
Note that this implies that and have the same law whenever , given identically distributed initial conditions.
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 McKeanVlasov 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,911]. 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 meanfield 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 meanfield equation:
where is a process independent of that has the same law, and denotes the expectation under the law of . In other words, the meanfield equation can be written, denoting by the law of (hence, also of ):
In these equations, , for , are independent, standard, mdimensional Brownian motions. Let us point out the fact that the righthand side of Equations 22 and 23 depends on the law of the solution; this fact is sometimes referred to as ‘the process is attracted by its own law’. This equation is also classically written as the McKeanVlasovFokkerPlanck equation on the probability distribution p of the solution. This equation which we use in the ‘Numerical simulations’ section can be easily derived from Equation 22. Let , , be the probability density at time t of the solution to Equation 22 (this is equivalent to ), then we have:
with
The P equations (Equation 24) yield the probability densities of the solutions of the meanfield equations (Equation 22). Because of the propagation of chaos result, the are statistically independent, but their probability functions are clearly functionally dependent.
We now spend some time on notations in order to obtain a somewhat more compact form of Equation 22. We define to be the dPdimensional process . We similarly define f, g, b and β as the concatenations of functions , , and , respectively. In details, , and . The term of noisy synaptic interactions requires a more careful treatment. We define and , and the product ⊙ of an element and an element as the element of :
We obtain the equivalent compact meanfield 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 wellposedness of the meanfield equation (Equation 22). The second is about the convergence of the solutions of the network equations to those of the meanfield 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 2Under assumptions (H1) to (H4), there exists a unique solution to the meanfield equation (Equation 22) onfor any.
Let us denote by the set of probability distributions on the set continuous functions , and the space of squareintegrable processes. Let (respectively, ) also be a family of P (respectively, )independent, m (respectively δ)dimensional, adapted standard Brownian motions on . Let us also note as the (random) initial condition of the meanfield equation. We introduce the map Φ acting on stochastic processes and defined by:
We have introduced in the previous formula the process with the same law as and independent of . There is a trivial identification between the solutions of the meanfield equation (Equation 22) and the fixed points of the map Φ: any fixed point of Φ provides a solution for Equation 22, and conversely, any solution of Equation 22 is a fixed point of Φ.
The following lemma is useful to prove the theorem:
Lemma 3Letbe a squareintegrable random variable. LetXbe a solution of the meanfield equation (Equation 22) with initial condition. Under assumptions (H3) and (H4), there exists a constantdepending on the parameters of the system and on the horizonT, such that:
Proof Using the Itô formula for , we have:
where is a stochastic integral, hence with a null expectation, .
This expression involves the term . Because of assumption (H3), we clearly have:
It also involves the term which, because of assumption (H4), is upperbounded by . Finally, assumption (H3) again allows us to upperbound the term by .
Finally, we obtain
Using Gronwall’s inequality, we deduce the boundedness of the solutions of the meanfield equations. □
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 secondorder moment.
Existence. Let be a given stochastic process, and define the sequence of probability distributions on defined by induction by . Define also a sequence of processes , , independent of the sequence of processes and having the same law. We note this as ‘X and Z i.i.d.’ below. We stop the processes at the time the first hitting time of the norm of to the constant value U. For convenience, we will make an abuse of notation in the proof and denote . This implies that belongs to , the ball of radius U centered at the origin in , for all times .
Using the notations introduced for Equation 25, we decompose the difference as follows:
and find an upperbound for by finding upperbounds for the corresponding norms of the four terms , , and . Applying the discrete CauchySchwartz inequality, we have:
and treat each term separately. The upperbounds for the first two terms are obtained using the CauchySchwartz inequality, those of the last two terms using the BurkholderDavisGundy martingale moment inequality.
The term is easily controlled using the CauchySchwarz inequality and the use of assumption (H1):
Taking the sup of both sides of the last inequality, we obtain
from which follows the fact that
The term is controlled using the CauchySchwartz inequality, assumption (H2), and the fact that the processes 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 is controlled using the fact that it is a martingale and applying the BurkholderDavisGundy martingale moment inequality and assumption (H1):
The term is also controlled using the fact that it is a martingale and applying the BurkholderDavisGundy martingale moment inequality and assumption (H2):
Putting all of these together, we get:
From the relation with , we get by an immediate recursion:
and is finite because the processes are bounded. The BienayméTchebychev inequality and Equation 27 now give
and this upper bound is the term of a convergent series. The BorelCantelli lemma stems that for almost any , there exists a positive integer (ω denotes an element of the probability space Ω) such that
and hence
It follows that with probability 1, the partial sums:
are uniformly (in ) convergent. Denote the thus defined limit by . It is clearly continuous and adapted. On the other hand, the inequality (Equation 27) shows that for every fixed t, the sequence is a Cauchy sequence in . Lemma 3 shows that .
It is easy to show using routine methods that indeed satisfies Equation 22.
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 and are globally Lipchitz continuous; hence, the previous proof shows that there exists a unique solution to equations (Equation 22) associated with the truncated functions. This solution satisfies the equation
Let us now define the stopping time as
It is easy to show that
implying that the sequence of stopping times is increasing. Using Lemma 3 which implies that the solution to Equation 22 is almost surely bounded, for almost all , there exists such that for all . Now, define , . Because of Equation 29, we have , and it follows from Equation 28 that
and letting , we have shown the existence of solution to Equation 22 which, by Lemma 3, is squareintegrable.
Uniqueness. Assume that X and Y are two solutions of the meanfield equations (Equation 22). From Lemma 3, we know that both solutions are in . Moreover, using the bound Equation 26, we directly obtain the inequality:
which, by Gronwall’s theorem, directly implies that
which ends the proof. □
We have proved the wellposedness of the meanfield equations. It remains to show that the solutions to the network equations converge to the solutions of the meanfield equations. This is what is achieved in the next theorem.
Theorem 4Under assumptions (H1) to (H4), the following holds true:
• Convergence^{c}: For each neuroniof populationα, the law of the multidimensional processconverges towards the law of the solution of the meanfield equation related to populationα, namely.
• Propagation of chaos: For any, and anyktuple, the law of the processconverges towards^{d}, i.e. the asymptotic processes have the law of the solution of the meanfield 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 meanfield equations (Equation 22) and of their variants (such as the McKeanVlasovFokkerPlanck 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 FitzHughNagumo 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 nonlocal. By this, we mean that in the case of the McKeanVlasov 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 FokkerPlanck 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 McKeanVlasov equations
The fact that the McKeanVlasov 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 EulerMaruyama method a number of times with different initial conditions, and averages the trajectories of the solutions over the number of simulations.
In detail, let and . The discretetime dynamics implemented in the stochastic numerical simulations consists of simulating M times a Ppopulation discretetime process , solution of the recursion, for i in population α:
where and are independent d and δdimensional standard normal random variables. The initial conditions , , are drawn independently from the same law within each population for each Monte Carlo simulation . One then chooses one neuron in each population . If the size N of the population is large enough, Theorem 4 states that the law, noted as , of should be close to that of the solution of the meanfield equations for . Hence, in effect, simulating the network is a good approximation (see below) of the simulation of the meanfield or McKeanVlasov equations [68,69]. An approximation of can be obtained from the Monte Carlo simulations by quantizing the phase space and incrementing the count of each bin whenever the trajectory of the neuron at time t falls into that particular bin. The resulting histogram can then be compared to the solution of the McKeanVlasovFokkerPlanck 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) and the solution of the meanfield equations (Equation 22) is of order , the first term being related to the error made by approximating the solution of the network of size N, by an EulerMaruyama method, and the second term, to the convergence of towards the meanfield equation when considering globally Lipschitz continuous dynamics (see proof of Theorem 4 in the Appendix). In our case, as shown before, the dynamics is only locally Lipschitz continuous. Finding efficient and provably convergent numerical schemes to approximate the solutions of such stochastic differential equations is an area of active research. There exist proofs that some schemes are divergent [70] or convergent [71] for some types of drift and diffusion coefficients. Since our equations are not included in either case, we conjecture convergence since we did not observe any divergence and leave the proof for future work.
4.2 Numerical simulations of the McKeanVlasovFokkerPlanck equation
For solving the McKeanVlasovFokkerPlanck 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 , of the probability density function of population α 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 solution of (Equation 24) at the sample points. The computation of the integral terms that appear in the McKeanVlasovFokkerPlanck equation is achieved through a recursive scheme, the NewtonCotes method of order 6 [74]. The dimensionality of the space being large and numerical errors increasing with the dimensionality of the integrand, such precise integration schemes are necessary. For an arbitrary real function f to be integrated between the values and , this numerical scheme reads:
where Δx is the integration step, and is chosen to be an integer multiple of 5.
The discretization of the derivatives with respect to the phase space parameters is done through the following fourthorder central difference scheme:
for the firstorder derivatives, and
for the secondorder derivatives (see [75]).
Finally, we have used a RungeKutta 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 meanfield equations
We illustrate these ideas with the example of a network of 100 FitzHughNagumo 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 points defined by
where , , , , and define the volume in which we solve the network equations and estimate the histogram defined in the ‘Numerical simulations of the McKeanVlasov equations’ section, while ΔV, Δw and Δy are the quantization steps in each dimension of the phase space. For the simulation of the McKeanVlasovFokkerPlanck 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 McKeanVlasovFokkerPlanck equation with the method of lines is the product (in our case, we chose ). This can become fairly large if we increase the precision of the phase space discretization. Moreover, increasing the precision of the simulation in the phase space, in order to ensure the numerical stability of the method of lines, requires to decrease the time step Δ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 lefthand parts of Figures 4 and 5, we have used one population of 100 excitatory FitzHughNagumo neurons connected with chemical synapses. We performed Monte Carlo simulations of the network equations (Equation 14) with the EulerMaruyama method in order to approximate the probability density. The model for the time variation of the synaptic weights is the simple model. The p.d.f. of the initial condition is Gaussian and reads
Fig. 4. Joint probability distribution. computed with the Monte Carlo algorithm for the network equations (Equation 14) (left) compared with the solution of the McKeanVlasovFokkerPlanck equation (Equation 24) (right), sampled at four times . Parameters are given in Table 1, with a current corresponding to a stable limit cycle. Initial conditions (first column of Table 1) are concentrated inside this limit cycle. The two distributions are similar and centered around the limit cycle with two peaks (see text).
Fig. 5. Joint probability distribution. computed with the Monte Carlo algorithm for the network equations (Equation 14) (left) compared with the solution of the McKeanVlasovFokkerPlanck equation (Equation 24) (right), sampled at four times . Parameters are given in Table 1, with a current corresponding to a stable limit cycle. Initial conditions (first column of Table 1) are concentrated inside this limit cycle. The two distributions are similar and centered around the limit cycle with two peaks (see text).
The parameters are given in the first column of Table 1. In this table, the parameter is the time at which we stop the computation of the trajectories in the case of the network equations and the computation of the solution of the McKeanVlasovFokkerPlanck equation in the case of the meanfield equations. The sequence indicates that we compute the solutions at those four time instants corresponding to the four rows of Figures 4 and 5. The phase space has been quantized with the parameters shown in the second column of the same table to solve the McKeanVlasovFokkerPlanck equation. This quantization has also been used to build the histograms that represent the marginal probability densities with respect to the pairs and of coordinates of the state vector of a particular neuron. These histograms have then been interpolated to build the surfaces shown in the lefthand side of Figures 4 and 5. The parameters of the FitzHughNagumo model are the same for each neuron of the population: they are shown in the third column of Table 1.
Table 1. Parameters used in the simulations of the neural network and for solving the McKeanVlasovFokkerPlanck 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 and , the probability that the maximum conductances change sign is very small. Finally, the parameters of the chemical synapses are shown in the sixth column. The parameters Γ and Λ are those of the χ function (Equation 3). The solutions are computed over an interval of time units with a time sampling of for the network and for the McKeanVlasovFokkerPlanck equation. The rest of the parameters are the typical values for the FitzHughNagumo equations.
The marginals estimated from the trajectories of the network solutions are then compared to those obtained from the numerical solution of the McKeanVlasovFokkerPlanck 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 for the external current (this value corresponds to the existence of a stable limit cycle for the isolated FitzHughNagumo neuron), and the initial conditions have the values , and ; therefore, the initial points of the trajectories in the phase space are concentrated inside the limit cycle. We therefore expect that the solutions of the neural network and the McKeanVlasovFokkerPlanck equation will concentrate their mass around the limit cycle. This is what is observed in Figures 4 and 5, where the simulation of the neural network (lefthand side) is in very good agreement with the results of the simulation of the McKeanVlasovFokkerPlanck equation (righthand side). Note that the densities display two peaks. These two peaks correspond to the fact that depending upon the position of the initial condition with respect to the nullclines of the FitzHughNagumo equations, the points in the phase space follow two different classes of trajectories, as shown in Figure 6. The two peaks then rotate along the limit cycle in the space (see also the ‘Numerical simulations with GPUs’ section).
Fig. 6. Projection of 100 trajectories in the (top left), (top right) and (bottom) planes. The limit cycle is especially visible in the projection (red curves). The initial conditions split the trajectories into two classes corresponding to the two peaks shown in Figures 4 and 5. The parameters are the same as those used to generate these two pictures.
Figures 4 and 5 show a qualitative similarity between the marginal probability density functions obtained by simulating the network and those obtained by solving the FokkerPlanck equation corresponding to the meanfield equations. To make this more quantitative, we computed the KullbackLeibler divergence between the two distributions.
We performed Monte Carlo simulations of the network equations up to for increasing values of the network size N. As shown in Figure 7, the KullbackLeibler divergence does decrease with increasing values of N, thereby confirming the fact that even for relatively small values of N, the average behavior of the network is well represented by the meanfield system described by the McKeanVlasovFokkerPlanck equation.
Fig. 7. Variation of the KullbackLeibler divergence. Variation of the KullbackLeibler divergence between the marginal probability density function estimated from the network equations and computed from the McKeanVlasovFokkerPlanck equation as a function of the network size. We have performed Monte Carlo simulations of the network equations up to time .
4.4 Numerical simulations with GPUs
Unfortunately, the algorithm for solving the McKeanVlasovFokkerPlanck equation described in the previous section is computationally very expensive. In fact, when the number of points in the discretized grid of the phase space is big, i.e. when the discretization steps ΔV, Δw and Δy are small, we also need to keep Δt small enough in order to guarantee the stability of the algorithm. This implies that the number of equations that must be solved has to be large and moreover that they must be solved with a small time step if we want to keep the numerical errors small. This will inevitably slow down the simulations. We have dealt with this problem by using a more powerful hardware, the graphical processing units (GPUs).
We have changed the RungeKutta scheme of order 2 used for the simulations shown in the ‘Numerical simulations of the McKeanVlasovFokkerPlanck equation’ section and adopted a more accurate RungeKutta scheme of order 4. This was done because with the more powerful machine, each computation of the righthand side of the equation is faster, making it possible to use four calls per time step instead of two in the previous method. Hence, the parallel hardware allowed us to use a more accurate method.
One of the purposes of the numerical study is to get a feeling for how the different parameters, in particular those related to the sources of noise, influence the solutions of the McKeanVlasovFokkerPlanck equation. This is meant to prepare the ground for the study of the bifurcation of these solutions with respect to these parameters, as was done in [76] in a different context. For this preliminary study, we varied the input current I and the parameter controlling the intensity of the noise on the membrane potential in Equations 14. The McKeanVlasovFokkerPlanck equation writes in this case:^{e}
The simulations were run with the χ function (Equation 3); the initial condition described by Equation 31 and the parameters are shown in Table 2. These parameters are similar to those used in the previous numerical simulations, but they differ in the size of the grid which is larger in this case.
Table 2. Parameters used in the simulations of the McKeanVlasovFokkerPlanck equation on GPUs
Four snapshots of the solution are shown in Figure 8 (corresponding to the values and of the external input current and of the standard deviation of the noise on the membrane potential), and three are shown in Figure 9 (corresponding to the values and ). In the figures, the left column corresponds to the values of the marginal , and the right column corresponds to the values of the marginal . Both are necessary to get an idea of the shape of the full distribution . The first row of Figure 8 shows the initial conditions. They are the same for the results shown in Figure 9. The second, third and fourth rows of Figure 8 show the time instants , and at convergence (the time units differ from those of the previous section, but it is irrelevant to this discussion). The three rows of Figure 9 show the time instants , and at convergence. In both cases, the solution appears to converge to a stationary distribution whose mass is distributed over a ‘blurred’ version of the limit cycle of the isolated neuron. The ‘blurriness’ increases with the variance of the noise. The four movies for these two cases are available as Additional files 1, 2, 3 and 4.
Fig. 8. Marginals of the solutions to the McKeanVlasovFokkerPlanck equation. Marginals with respect to the V and w variables (left) and to the V and y variables (right) of the solution of the McKeanVlasovFokkerPlanck equation. The first row shows the initial condition; the second, the marginals at time 30.0; the third, the marginals at time 50.0; and the fourth, the stationary (large time) solutions. The input current I is equal to 0.4 and . These are screenshots at different times of movies available as Additional files 1 and 2.
Fig. 9. Marginals of the solutions to the McKeanVlasovFokkerPlanck equation. Marginals with respect to the V and w variables (left) and to the V and y variables (right) of the solution of the McKeanVlasovFokkerPlanck equation. The first row shows the marginals at time 30.0, the second the marginals at time 50.0 and the third the stationary (large time) solutions. The input current I is equal to 0.7 and . These are screenshots at different times of movies available as Additional files 3 and 4.
The results shown in Figures 8 and 9 and in Additional files 1, 2, 3 and 4 were obtained using two machines, each with seven nVidia Tesla C2050 cards, six 2.66 GHz dualXeon X5650 processors and 72G of ram. The communication inside each machine was done using the lpthreads library and between machines using MPI calls. The mean execution time per time step using the parameters already described is 0.05 s.
The reader interested in more details in the numerical implementations and in the gains that can be achieved by the use of GPUs can consult [77].
In Figure 10, we show a solution to the McKeanVlasovFokkerPlanck equation which is qualitatively quite different from the solutions shown in Figures 8 and 9: The stationary solution is concentrated at a point in space. This is an indication that perhaps, between the values −0.8 and 0.4 of the input current, the solutions to the McKeanVlasovFokkerPlanck equation have bifurcated. The numerical tools we have developed may be a way to build an intuition to guide a rigorous analysis of these phenomena.
Fig. 10. Marginals of the solutions to the McKeanVlasovFokkerPlanck equation at convergence. Marginals with respect to the V and w variables (left) and to the V and y variables (right) of the solution of the McKeanVlasovFokkerPlanck equation at convergence. The parameters are those in Table 1 except for the input current I which is equal to −0.8, and . Compare with the last row of Figure 9 (see text).
5 Discussion and conclusion
In this article, we addressed the problem of the limit in law of networks of biologically inspired neurons as the number of neurons tends to infinity. We emphasized the necessity of dealing with biologically inspired models and discussed at length the type of models relevant to this study. We chose to address the case conductancebased network models that are a relevant description of the neuronal activity. Mathematical results on the analysis of these diffusion processes in interaction resulted to the replacement of a set of NPddimensional coupled equations (the network equations) in the limit of large Ns by Pddimensional meanfield equations describing the global behavior of the network. However, the price to pay for this reduction was the fact that the resulting meanfield equations are nonstandard stochastic differential equations, similar to the McKeanVlasov equations. These can be expressed either as implicit equations on the law of the solution or, in terms of probability density function through the McKeanVlasovFokkerPlanck equations, as a nonlinear, nonlocal partial differential equation. These equations are, in general, hard to study theoretically.
Besides the fact that we explicitly model real spiking neurons, the mathematical part of our work differs from that of previous authors such as McKean, Tanaka and Sznitman (see the ‘Introduction’ section) because we are considering several populations with the effect that the analysis is significantly more complicated. Our hypotheses are also more general, e.g. the drift and diffusion functions are nontrivial and satisfy the general condition (H4) which is more general than the usual linear growth condition. Also, they are only assumed locally (and not globally) Lipschitz continuous to be able to deal, for example, with the FitzHughNagumo model. A locally Lipschitz continuous case was recently addressed in a different context for a model of swarming in [67].
Proofs of our results, for somewhat stronger hypotheses than ours and in special cases, are scattered in the literature, as briefly reviewed in the ‘Introduction’ and ‘Setting of the problem’ sections. Our main contribution is that we provide a complete, selfsufficient proof in a fairly general case by gathering all the ingredients that are required for our neuroscience applications. In particular, the case of the FitzHughNagumo model where the drift function does not satisfy the linear growth condition involves a generalization of previous works using the more general growth condition (H4).
The simulation of these equations can itself be very costly. We, hence, addressed in the ‘Numerical simulations’ section numerical methods to compute the solutions of these equations, in the probabilistic framework, using the convergence result of the network equations to the meanfield limit and standard integration methods of differential equations or in the FokkerPlanck framework. The simulations performed for different values of the external input current parameter and one of the parameters controlling the noise allowed us to show that the spatiotemporal shape of the probability density function describing the solution of the McKeanVlasovFokkerPlanck equation was sensitive to the variations of these parameters, as shown in Figures 8 and 9. However, we did not address the full characterization of the dynamics of the solutions in the present article. This appears to be a complex question that will be the subject of future work. It is known that for different McKeanVlasov equations, stationary solutions of these equations do not necessarily exist and, when they do, are not necessarily unique (see [78]). A very particular case of these equations was treated in [76] where the authors consider that the function is linear, is constant and . This model, known as the firingrate model, is shown in that paper to have the Gaussian solutions when the initial data is Gaussian, and the dynamics of the solutions can be exactly reduced to a set of 2Pcoupled ordinary differential equations governing the mean and the standard deviation of the solution. Under these assumptions, a complete study of the solutions is possible, and the dependence upon the parameters can be understood through bifurcation analysis. The authors show that intrinsic noise levels govern the dynamics, creating or destroying fixed points and periodic orbits.
The meanfield description has also deep theoretical implications in neuroscience. Indeed, it points towards the fact that neurons encode their responses to stimuli through probability distributions. This type of coding was evoked by several authors [47], and the meanfield approach shows that under some mild conditions, this phenomenon arises: all neurons belonging to a particular population can be seen as independent realizations of the same process, governed by the meanfield equation. The relevance of this phenomenon is reinforced by the fact that it has recently been observed experimentally that neurons had correlation levels significantly below what had been previously reported [13]. This independence has deep implications on the efficiency of neural coding which the propagation of chaos theory accounts for. To illustrate this phenomenon, we have performed the following simulations. Considering a network of 2, 10 and 100 FitzHughNagumo neurons, we have simulated 2,000 times the network equations over some time interval . We have picked at random a pair of neurons and computed the time variation of the crosscorrelation of the values of their state variables. The results are shown in Figure 11. It appears that the propagation of chaos is observable for relatively small values of the number of neurons in the network, thus indicating once more that the theory developed in this paper in the limit case of an infinite number of neurons is quite robust to finitesize effects.^{f}
Fig. 11. Variations over time of the crosscorrelation of variables of several FitzHughNagumo neurons in a network. Top left: 2 neurons. Top right: 10 neurons. Bottom: 100 neurons. The crosscorrelation decreases steadily with the number of neurons in the network.
The present study develops theoretical arguments to derive the meanfield equations resulting from the activity of large neuron ensembles. However, the rigorous and formal approach developed here does not allow direct characterization of brain states. The paper, however, opens the way to rigorous analysis of the dynamics of large neuron ensembles through derivations of different quantities that may be relevant. A first approach could be to derive the equations of the successive moments of the solutions. Truncating this expansion would yield systems of ordinary differential equations that can give approximate information on the solution. However, the choice of the number of moments taken into account is still an open question that can raise several deep questions [46].
Appendix 1: Proof of Theorem 4
In this appendix, we prove the convergence of the network equations towards the meanfield equations (Equation 22) and of the propagation of chaos property. The proof follows standard proofs in the domain as generally done, in particular by Tanaka or Sznitman [6,10], adapted to our particular case where we consider a nonzero drift function and a time and spacedependent diffusion function. It is based on the very powerful coupling argument, which identifies the almost sure limit of the process as the number of neurons tends to infinity, as popularized by Sznitman in [12], but whose idea dates back from the 1970s (for instance, Dobrushin uses it in [5]). This process is exactly the solution of the meanfield equation driven by the same Brownian motion as and with the same initial condition random variable. In our case, this leads us to introduce the sequence of independent stochastic processes having the same law as , , solution of the meanfield equation:
with initial condition , the initial condition of the neuron i in the network, which was assumed to be independent and identically distributed. and are the Brownian motions involved in the network equation (Equation 21). As described previously, is a process independent of that has the same law. Denoting, as described previously, the probability distribution of solution of the meanfield equation (Equation 22) by , the law of the collection of processes for some fixed , namely , is shown to be the limit of the process solution of the network equations (Equation 21) as N goes to infinity.
We recall, for completeness, Theorem 4:
Theorem 4Under assumptions (H1) to (H4), the following holds true:
• Convergence: For each neuroniof populationα, the law of the multidimensional processconverges towards the law of the solution of the meanfield equation related to populationα, namely.
• Propagation of chaos: For any, and anykuplet, the law of the processconverges towards, i.e. the asymptotic processes have the law of the solution of the meanfield equations and are all independent.
Proof
On our way, we also prove that
which implies, in particular, convergence in law of the process towards solution of the meanfield equations (Equation 22).
The proof basically consists of thoroughly analyzing the difference between the two processes as N tends to infinity. The difference is the sum of eight terms (we dropped the index N for the sake of simplicity of notations) denoted by through :
It is important to note that the probability distribution of these terms does not depend on the neuron i. We are interested in the limit, as N goes to infinity, of the quantity . We decompose this expression into the sum of the eight terms involved in Equation 35 using Hölder’s inequality and upperbound each term separately. The terms and are treated exactly as in the proof of Theorem 2. We start by assuming that f and g are uniformly globally K Lipschitz continuous with respect to the second variable. The locally Lipschitz case is treated in the same manner as done in the proof of Theorem 2 (1) by stopping the process at time , (2) by using the Lipschitz continuity of f and g in the ball of radius U and (3) by a truncation argument and using the almost sure boundedness of the solutions extending the convergence to the locally Lipschitz case.
As seen previously, we have:
Therefore,
Hence, we have:
Therefore,
The terms and are treated in the same fashion, but instead of using the CauchySchwartz inequality, the BurkholderDavisGundy martingale moment inequality are used. For , in detail,
We are left with the problem of controlling the terms and that involve sums of processes with bounded second moment, thanks to Proposition 3 and assumption (H3). We have:
and using the BurkholderDavisGundy martingale moment inequality,
Each of these two expressions involves an expectation which we write:
All the terms of the sum corresponding to indexes j and k such that the three conditions , and are satisfied are null since in that case, , , and are independent and have the same law for .^{g} In effect, denoting the measure of their common law by , we have:
expanding further and renaming the second z variable to y in the last term, we obtain:
which is indeed equal to 0 by the Fubini theorem.
Therefore, there are no more than nonnull terms in the sum, and all the terms have the same value (that depends on Θ), which is bounded by Lemma 3 and assumption (H3). We denote the supremum of these values for across all possible pairs of populations by , and the smallest value of the by . We have shown that
Finally, we have:
for some positive constants and . Using Gronwall’s inequality, we obtain:
for some positive constant . The righthand side of this inequality tends to zero as N goes to infinity proving the propagation of chaos property. In order to show a convergence with speed as stated in the theorem, we use the fact:
and the righthand side of the inequality is bounded for all Ns because of the hypothesis for . This ends the proof. □
Electronic Supplementary Material
Additional file 1. Time evolution of the marginal of the solution to the McKeanVlasovFokkerPlanck equation. The four images in the left part of Figure 8 are four snapshots of this movie taken at time 0 (initial condition), time 30, time 50 and at a large enough time for the solution to be stationary. The input current is equal to 0.4, and the standard deviation of the membrane potential noise, to 0.27. (AVI 2.0 MB)
Format: AVI Size: 5.3MB Download file  Watch movie
Additional file 2. Time evolution of the marginal of the solution to the McKeanVlasovFokkerPlanck equation. The four images in the right part of Figure 8 are four snapshots of this movie taken at time 0 (initial condition), time 30, time 50 and at a large enough time for the solution to be stationary. The input current is equal to 0.4, and the standard deviation of the membrane potential noise, to 0.27. (AVI 1.5 MB)
Format: AVI Size: 4.2MB Download file  Watch movie
Additional file 3. Time evolution of the marginal of the solution to the McKeanVlasovFokkerPlanck equation. The three images in the left part of Figure 9 are three snapshots of this movie taken at time 30, time 50 and at a large enough time for the solution to be stationary. The input current is equal to 0.7, and the standard deviation of the membrane potential noise, to 0.45. (AVI 3.0 MB)
Format: AVI Size: 3.5MB Download file  Watch movie
Additional file 4. Time evolution of the marginal of the solution to the McKeanVlasovFokkerPlanck equation. The three images in the right part of Figure 9 are three snapshots of this movie taken at time 30, time 50, and at a large enough time for the solution to be stationary. The input current is equal to 0.7, and the standard deviation of the membrane potential noise, to 0.45. (AVI 2.2 MB)
Format: AVI Size: 3.4MB Download file  Watch movie
Competing interests
The authors declare that they have no competing interests.
Authors’ contributions
JB and DF developed the code for solving the stochastic differential equations, the McKeanVlasov equations and the McKeanVlasovFokkerPlanck equations. They ran the numerical experiments and generated all the figures. DF derived some of the McKeanVlasov equations in a heuristic fashion. OF and JT developed the models, proved the theorems and wrote the paper. All authors read and approved the final manuscript.
Acknowledgements
This work was partially supported by the ERC grant #227747 NerVi, the FACETSITN MarieCurie Initial Training Network #237955 and the IP project BrainScaleS #269921.
End notes

More precisely, as shown in [79,80], the convergence is to a larger  13dimensional  system with an invariant fourdimensional manifold on which the solution lives given appropriate initial conditions. See also [81].

As we will see in the proof, most properties are valid as soon as tends to infinity as N goes to infinity for all , the previous assumption will allow quantifying the speed of convergence towards the asymptotic regime.

The type of convergence is specified in the proof given in the Appendix.

We have included a small noise (controlled by the parameter ) on the adaptation variable w. This does not change the previous analysis, in particular proposition 1, but makes the McKeanVlasovFokkerPlanck equation wellposed in a cube of the state space with 0 boundary value, see e.g. [82].

Note that we did not estimate the correlation within larger networks since, as predicted by Theorem 4, it will be smaller and smaller, requiring an increasingly large number of Monte Carlo simulations.

Note that and as soon as . In the case where , it is easy to check that when j (respectively, k) is equal to i, all terms such that (respectively, ) are equal to 0.
References

McKean H: A class of Markov processes associated with nonlinear parabolic equations.
Proc Natl Acad Sci USA 1966, 56(6):19071911. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

McKean H: Propagation of chaos for a class of nonlinear parabolic equations. In Stochastic Differential Equations. Air Force Office Sci. Res., Arlington; 1967:4157. [Lecture Series in Differential Equations 7]

Braun W, Hepp K: The Vlasov dynamics and its fluctuations in the 1/n limit of interacting classical particles.
Commun Math Phys 1977, 56(2):101113. Publisher Full Text

Dawson D: Critical dynamics and fluctuations for a meanfield model of cooperative behavior.
J Stat Phys 1983, 31:2985. Publisher Full Text

Dobrushin RL: Prescribing a system of random variables by conditional distributions.
Theory Probab Appl 1970, 15:458486. Publisher Full Text

Tanaka H: Probabilistic treatment of the Boltzmann equation of Maxwellian molecules.

Tanaka H, Hitsuda M: Central limit theorem for a simple diffusion model of interacting particles.

Tanaka H: Some probabilistic problems in the spatially homogeneous Boltzmann equation. In Theory and Application of Random Fields Lecture Notes in Control and Information Sciences. Edited by Kallianpur G. Springer, Berlin; 1983:258267.

Tanaka H: Limit theorems for certain diffusion processes with interaction. In Stochastic Analysis. NorthHolland, Amsterdam; 1984:469488. [NorthHolland Mathematical Library 32]

Sznitman A: Nonlinear reflecting diffusion process, and the propagation of chaos and fluctuations associated.
J Funct Anal 1984, 56(3):311336. Publisher Full Text

Sznitman A: A propagation of chaos result for Burgers’ equation.
Probab Theory Relat Fields 1986, 71(4):581613. Publisher Full Text

Sznitman AS: Topics in propagation of chaos. In Ecole d’Eté de Probabilités de SaintFlour XIX 1989. Edited by Burkholder D, Pardoux E, Sznitman AS. Springer, Berlin; 1991:165251. [Lecture Notes in Math. 1464] PubMed Abstract  Publisher Full Text

Ecker A, Berens P, Keliris G, Bethge M, Logothetis N, Tolias A: Decorrelated neuronal firing in cortical microcircuits.
Science 2010, 327(5965):584. PubMed Abstract  Publisher Full Text

Hodgkin A, Huxley A: A quantitative description of membrane current and its application to conduction and excitation in nerve.
J Physiol 1952, 117:500544. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Fitzhugh R: Theoretical effect of temperature on threshold in the HodgkinHuxley nerve model.
J Gen Physiol 1966, 49(5):9891005. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

FitzHugh R: Mathematical models of excitation and propagation in nerve. In Biological Engineering. Edited by Schwan HP. McGrawHill Book Co., New York; 1969:185.

Izhikevich EM: Dynamical Systems in Neuroscience: The Geometry of Excitability and Bursting. MIT Press, Cambridge; 2007.

Lapicque L: Recherches quantitatifs sur l’excitation des nerfs traitee comme une polarisation.
J. Physiol. Paris 1907, 9:620635. PubMed Abstract

Tuckwell HC: Introduction to Theoretical Neurobiology. Cambridge University Press, Cambridge; 1988.

Ermentrout GB, Terman D: Foundations of Mathematical Neuroscience.. Springer, Berlin; 2010. [Interdisciplinary Applied Mathematics]

FitzHugh R: Mathematical models of threshold phenomena in the nerve membrane.

Nagumo J, Arimoto S, Yoshizawa S: An active pulse transmission line simulating nerve axon.

Destexhe A, Mainen Z, Sejnowski T: Synthesis of models for excitable membranes, synaptic transmission and neuromodulation using a common kinetic formalism.
J Comput Neurosci 1994, 1(3):195230. PubMed Abstract  Publisher Full Text

Kandel ER, Schwartz JH, Jessel TM: Principles of Neural Science. 4th edition. McGrawHill, New York; 2000.

Cox JC, Ingersoll JC Jr, Ross SA: A theory of the term structure of interest rates.
Econometrica 1985, 53(2):385407. Publisher Full Text

Mao X: Stochastic Differential Equations and Applications. 2nd edition. Horwood, Chichester; 2008.

Amari S: Characteristics of random nets of analog neuronlike elements.

Amari S: Dynamics of pattern formation in lateralinhibition type neural fields.
Biol Cybern 1977, 27(2):7787. PubMed Abstract  Publisher Full Text

Wilson H, Cowan J: Excitatory and inhibitory interactions in localized populations of model neurons.
Biophys J 1972, 12:124. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Wilson H, Cowan J: A mathematical theory of the functional dynamics of cortical and thalamic nervous tissue.

Hammerstein A: Nichtlineare Integralgleichungen nebst Anwendungen.
Acta Math 1930, 54:117176. Publisher Full Text

Faugeras O, Grimbert F, Slotine JJ: Absolute stability and complete synchronization in a class of neural fields models.
SIAM J Appl Math 2008, 61:205250. Publisher Full Text

Coombes S, Owen MR: Bumps, breathers, and waves in a neural network with spike frequency adaptation.
Phys Rev Lett 2005., 94(14)
Article ID 148102

Ermentrout B: Neural networks as spatiotemporal patternforming systems.
Rep Prog Phys 1998, 61:353430. Publisher Full Text

Ermentrout G, Cowan J: Temporal oscillations in neuronal nets.
J Math Biol 1979, 7(3):265280. PubMed Abstract  Publisher Full Text

Laing C, Troy W, Gutkin B, Ermentrout G: Multiple bumps in a neuronal model of working memory.
SIAM J Appl Math 2002, 63:6297. Publisher Full Text

Chossat P, Faugeras O: Hyperbolic planforms in relation to visual edges and textures perception.
PLoS Comput Biol 2009., 5(12)
Article ID e1000625

Veltz R, Faugeras O: Local/global analysis of the stationary solutions of some neural field equations.
SIAM J Appl Dyn Syst 2010, 9(3):954998. Publisher Full Text

Abbott L, Van Vreeswijk C: Asynchronous states in networks of pulsecoupled neuron.
Phys Rev 1993, 48:14831490. Publisher Full Text

Amit D, Brunel N: Model of global spontaneous activity and local structured delay activity during delay periods in the cerebral cortex.
Cereb Cortex 1997, 7:237252. PubMed Abstract  Publisher Full Text

Brunel N, Hakim V: Fast global oscillations in networks of integrateandfire neurons with low firing rates.
Neural Comput 1999, 11:16211671. PubMed Abstract  Publisher Full Text

Brunel N: Dynamics of sparsely connected networks of excitatory and inhibitory spiking neurons.
J Comput Neurosci 2000, 8:183208. PubMed Abstract  Publisher Full Text

El Boustani S, Destexhe A: A master equation formalism for macroscopic modeling of asynchronous irregular activity states.
Neural Comput 2009, 21:46100. PubMed Abstract  Publisher Full Text

Mattia M, Del Giudice P: Population dynamics of interacting spiking neurons.
Phys Rev E, Stat Nonlinear Soft Matter Phys 2002., 66(5)
Article ID 51917

Cai D, Tao L, Shelley M, McLaughlin DW: An effective kinetic representation of fluctuationdriven neuronal networks with application to simple and complex cells in visual cortex.
Proc Natl Acad Sci USA 2004, 101(20):77577762. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Ly C, Tranchina D: Critical analysis of dimension reduction by a moment closure method in a population density approach to neural network modeling.
Neural Comput 2007, 19(8):20322092. PubMed Abstract  Publisher Full Text

Rolls ET, Deco G: The Noisy Brain: Stochastic Dynamics as a Principle of Brain Function. Oxford University Press, Oxford; 2010.

Softky WR, Koch C: The highly irregular firing of cortical cells is inconsistent with temporal integration of random EPSPs.
J Neurosci 1993, 13:334350. PubMed Abstract  Publisher Full Text

Brunel N, Latham P: Firing rate of noisy quadratic integrateandfire neurons.
Neural Comput 2003, 15:22812306. PubMed Abstract  Publisher Full Text

Plesser HE: Aspects of signal processing in noisy neurons. PhD thesis. GeorgAugustUniversität; 1999.

Touboul J, Faugeras O: First hitting time of double integral processes to curved boundaries.
Adv Appl Probab 2008, 40(2):501528. Publisher Full Text

Beggs JM, Plenz D: Neuronal avalanches are diverse and precise activity patterns that are stable for many hours in cortical slice cultures.
J Neurosci 2004, 24(22):52165229. PubMed Abstract  Publisher Full Text

Benayoun M, Cowan JD, van Drongelen W, Wallace E: Avalanches in a stochastic model of spiking neurons.
PLoS Comput Biol 2010., 6(7)
Article ID e1000846

Levina A, Herrmann JM, Geisel T: Phase transitions towards criticality in a neural system with adaptive interactions.
Phys Rev Lett 2009., 102(11)
Article ID 118110

Touboul J, Destexhe A: Can powerlaw scaling and neuronal avalanches arise from stochastic dynamics?
PLoS ONE 2010., 5(2)
Article ID e8982

Bressloff P: Stochastic neural field theory and the systemsize expansion.

Buice MA, Cowan JD: Fieldtheoretic approach to fluctuation effects in neural networks.
Phys Rev E, Stat Nonlinear Soft Matter Phys 2007., 75(5)
Article ID 051919

Buice M, Cowan J, Chow C: Systematic fluctuation expansion for neural network activity equations.
Neural Comput 2010, 22(2):377426. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Ohira T, Cowan J: Masterequation approach to stochastic neurodynamics.
Phys Rev E, Stat Nonlinear Soft Matter Phys 1993, 48(3):22592266. Publisher Full Text

Treves A: Meanfield analysis of neuronal spike dynamics.
Network 1993, 4(3):259284. Publisher Full Text

Gerstner W: Time structure of the activity in neural network models.
Phys Rev E, Stat Nonlinear Soft Matter Phys 1995, 51:738758. Publisher Full Text

Faugeras O, Touboul J, Cessac B: A constructive meanfield analysis of multipopulation neural networks with random synaptic weights and stochastic inputs.
Front Comput Neurosci 2009.
doi:10.3389/neuro.10.001.2009

Guionnet A: Averaged and quenched propagation of chaos for spin glass dynamics.
Probab Theory Relat Fields 1997, 109(2):183215. Publisher Full Text

Chizhov AV, Graham LJ: Population model of hippocampal pyramidal neurons, linking to refractory density approach to conductancebased neurons.
Phys Rev E, Stat Nonlinear Soft Matter Phys 2007., 75
Article ID 011924

Sompolinsky H, Crisanti A, Sommers H: Chaos in random neural networks.
Phys Rev Lett 1988, 61(3):259262. PubMed Abstract  Publisher Full Text

Sompolinsky H, Zippelius A: Relaxational dynamics of the EdwardsAnderson model and the meanfield theory of spinglasses.
Phys Rev B, Condens Matter Mater Phys 1982, 25(11):68606875. Publisher Full Text

Bolley F, Cañizo JA, Carrillo JA: Stochastic meanfield limit: nonLipschitz forces and swarming.
Math Models Methods Appl Sci 2011, 21(11):21792210. Publisher Full Text

Talay D, Vaillant O: A stochastic particle method with random weights for the computation of statistical solutions of McKeanVlasov equations.

Bossy M, Talay D: A stochastic particle method for the McKeanVlasov and the Burgers equation.
Math Comput 1997, 66(217):157192. Publisher Full Text

Hutzenthaler M, Jentzen A, Kloeden P: Strong and weak divergence in finite time of Euler’s method for stochastic differential equations with nonglobally Lipschitz continuous coefficients.
Proc R Soc, Math Phys Eng Sci 2011, 467(2130):15631576. Publisher Full Text

Hutzenthaler M, Jentzen A: Convergence of the stochastic Euler scheme for locally Lipschitz coefficients.
Found Comput Math 2011, 11(6):657706. Publisher Full Text

Schiesser W: The Numerical Method of Lines: Integration of Partial Differential Equations. Academic Press, San Diego; 1991.

Schiesser WE, Griffiths GW: A Compendium of Partial Differential Equation Models: Method of Lines Analysis with Matlab. 1st edition. Cambridge University Press, New York; 2009.

Ueberhuber CW: Numerical Computation 2: Methods, Software, and Analysis. Springer, Berlin; 1997.

Morton KW, Mayers DF: Numerical Solution of Partial Differential Equations: An Introduction. Cambridge University Press, Cambridge; 2005.

Touboul J, Hermann G, Faugeras O: Noiseinduced behaviors in neural mean field dynamics.
SIAM J Appl Dyn Syst 2012, 11(1):4981. Publisher Full Text

Baladron J, Fasoli D, Faugeras O: Three applications of GPU computing in neuroscience. Comput Sci Eng 2012, 14:4047.

Herrmann S, Tugaut J: Nonuniqueness of stationary measures for selfstabilizing processes.
Stoch Process Appl 2010, 120(7):12151246. Publisher Full Text

Pakdaman K, Thieullen M, Wainrib G: Fluid limit theorems for stochastic hybrid systems with application to neuron models.
Adv Appl Probab 2010, 42(3):761794. Publisher Full Text

Wainrib G: Randomness in neurons: a multiscale probabilistic analysis. PhD thesis. Ecole Polytechnique; 2010.

Goldwyn JH, Imennov NS, Famulare M, SheaBrown E: Stochastic differential equation models for ion channel noise in HodgkinHuxley neurons.
Phys Rev E, Stat Nonlinear Soft Matter Phys 2011., 83(4)
Article ID 041908

Evans LC: Partial Differential Equations. American Mathematical Society, Providence; 1998. [Graduate Studies in Mathematics 19]