Mean-field limit of interacting 2D nonlinear stochastic spiking neurons

TL;DR

This paper develops a mean-field model for a network of stochastic spiking neurons, derives its limit, and validates it through numerical simulations showing propagation of chaos and bifurcation phenomena.

Contribution

It introduces a new nonlinear stochastic model of neuron networks, derives its mean-field limit, and provides numerical methods to analyze its properties.

Findings

Propagation of chaos observed numerically

Convergence from microscopic to mean-field model

Evidence of Hopf bifurcation at high connectivity

Abstract

In this work, we propose a nonlinear stochastic model of a network of stochastic spiking neurons. We heuristically derive the mean-field limit of this system. We then design a Monte Carlo method for the simulation of the microscopic system, and a finite volume method (based on an upwind implicit scheme) for the mean-field model. The finite volume method respects numerical versions of the two main properties of the mean-field model, conservation and positivity, leading to existence and uniqueness of a numerical solution. As the size of the network tends to infinity, we numerically observe propagation of chaos and convergence from an individual description to a mean-field description. Numerical evidences for the existence of a Hopf bifurcation (synonym of synchronised activity) for a sufficiently high value of connectivity, are provided.

Tables

Table 1. Table 1: Model parameters

Parameter/function	CV test	Invariant distributions	Hopf test
$I\left(t\right)$	$2$	$-2,1$	$0$
${\tau }_w$	$1$	$2$	$13$
$b$	$1$	$1,0.05$	$0.011$
$\overline{v}$	$1.0$	$1.8$	$-1.5$
$\overline{w}$	$1.5$	$5.5,1.5$	$1.5$
$J$	$3.1$	$0$	$5$
$\lambda \left(v\right)$	$0.1+e^{\left(v-1\right)}$	$e^v$	$e^v$
$F\left(v\right)$	$e^v-5v$	$e^v-v$	$e^v-v$

Equations

$${\dot{v}}_i(t)$$

$${\dot{w}}_i(t)$$

$$F_1(v)$$

$$F_2(v)$$

$$F_3(v)$$

$$\displaystyle\bigl{(}v_{i}(t),w_{i}(t)\bigr{)}$$

$$\displaystyle\bigl{(}v_{j}(t),w_{j}(t)\bigr{)}$$

$$\{\begin{array}{rl}{v}_i^N(t)& =v_{i,0}^N+{\int }_0^t\tilde{\mathcal{V}}(v_i^N(s),w_i^N(s))\text{d}s\\ & +{\int }_0^t{\int }_0^{\infty }(\bar{v}-V_i^N(s^{-})){\mathbf{missing}1}_{\{z\le \lambda (v_i^N(s^{-}))\}}{\mathcal{N}}^i(\text{d}z,\text{d}s)\\ & +\frac{J}{N}\sum _{j\ne i}{\int }_0^t{\int }_0^{\infty }{\mathbf{missing}1}_{\{z\le \lambda (v_j^N(s^{-}))\}}{\mathcal{N}}^j(\text{d}z,\text{d}s),\\ w_i^N(t)& =w_{i,0}^N+{\int }_0^t\mathcal{W}(v_i^N(s),w_i^N(s))\text{d}s\\ & +{\int }_0^t{\int }_0^{\infty }\bar{w}{\mathbf{missing}1}_{\{z\le \lambda (v_i^N(s^{-}))\}}{\mathcal{N}}^i(\text{d}z,\text{d}s),\end{array}$$

$${\Psi }^N(t)\stackrel{def.}{=}\frac{1}{N}\sum _{j\ne i}{\int }_0^t{\int }_0^{\infty }{\mathbf{missing}1}_{\{z\le \lambda (v_j^N(s^{-}))\}}{\mathcal{N}}^j(\text{d}z,\text{d}s).$$

$${\Psi }^N(t)\underset{N\mathrm{large}}{\simeq }\Psi (t)\stackrel{def.}{=}\mathbb{E}{\int }_0^t{\int }_0^{\infty }{\mathbf{missing}1}_{\{z\le \lambda (v_1^N(s^{-}))\}}{\mathcal{N}}^1(\text{d}z,\text{d}s)$$

$$\Psi (t)=\mathbb{E}{\int }_0^t{\int }_0^{\infty }{\mathbf{missing}1}_{\{z\le \lambda (v_1^N(s^{-}))\}}\text{d}z\text{d}s={\int }_0^t\psi (s)\text{d}s\text{ with }\psi (t)\stackrel{def.}{=}\mathbb{E}(\lambda (v(t))).$$

$${\mu }^N(t,\text{d}v,\text{d}w)$$

$$\{\begin{array}{rl}v(t)& =v_0+{\int }_0^t\tilde{\mathcal{V}}(v(s),w(s))\text{d}s\\ & +{\int }_0^t{\int }_0^{\infty }(\bar{v}-v(s^{-})){\mathbf{missing}1}_{\{z\le \lambda (v(s^{-}))\}}\mathcal{N}(\text{d}z,\text{d}s)\\ & +J{\int }_0^t\mathbb{E}(\lambda (v(s)))\text{d}s,\\ w(t)& =w_0+{\int }_0^t\mathcal{W}(v(s),w(s))\text{d}s\\ & +{\int }_0^t{\int }_0^{\infty }\bar{w}{\mathbf{missing}1}_{\{z\le \lambda (v(s^{-}))\}}\mathcal{N}(\text{d}z,\text{d}s)\end{array}$$

$${\mathcal{L}}_{\mu (t)}\phi (v,w)={\mathcal{L}}_{\mu (t)}^{\text{t}}\phi (v,w)+{\mathcal{L}}^{\text{j}}\phi (v,w),$$

$${\mathcal{L}}_{\mu }^{\text{t}}\phi (v,w)$$

$${\mathcal{L}}^{\text{j}}\phi (v,w)$$

$$\displaystyle\qquad\qquad=\lambda(v)\bigl{(}\varphi(\bar{v},\bar{w}+w)-\varphi(v,w)\bigr{)}\,,$$

$${\mathcal{V}}_{\mu }(v,w)$$

$$\pi (v,w,\text{d}{v}^{\prime },\text{d}{w}^{\prime })$$

$$\frac{\text{d}}{\text{d}t}⟨\mu (t),\phi ⟩$$

$$⟨{\mathcal{L}}_{\mu }^{\ast }\mu ,\phi ⟩$$

$$=⟨\mu ,{\mathcal{V}}_{\mu }{\partial }_v\phi +\mathcal{W}{\partial }_w\phi ⟩+⟨{\delta }_{\bar{v}}\otimes {\int }_{\mathbb{R}}\lambda (v^{\prime }){\sigma }_{\bar{w}}\mu (t,\text{d}{v}^{\prime },\cdot ),\phi ⟩-⟨\mu ,\lambda \phi ⟩$$

$$=⟨{\mathcal{L}}_{\mu }^{\text{t}\ast }\mu +{\mathcal{L}}^{\text{j}\ast }\mu ,\phi ⟩,$$

$$\displaystyle\frac{\partial}{\partial t}\mu(t)=\mathcal{L}^{*}_{\mu(t)}\,\mu(t)=\bigl{(}\mathcal{L}^{\textrm{\tiny t}*}_{\mu(t)}+\mathcal{L}^{\textrm{\tiny j}*}\bigr{)}\,\mu(t)$$

$$\displaystyle\frac{\partial}{\partial t}\mu(t,v,w)+\frac{\partial}{\partial v}\bigl{(}\mu(t,v,w)\,\mathcal{V}_{\mu(t)}(v,w)\bigr{)}+\frac{\partial}{\partial w}\bigl{(}\mu(t,v,w)\,\mathcal{W}(v,w)\bigr{)}$$

$$=-\lambda (v)\mu (t,v,w)+{\delta }_{\bar{v}}(v){\int }_{\mathbb{R}}\lambda (v^{\prime }){\sigma }_{\bar{w}}\mu (t,v^{\prime },w)\text{d}{v}^{\prime },\mu (0)={\mu }_0$$

$$\Gamma$$

$${\Omega }_1$$

$$\displaystyle\bigl{\langle}\mu(t),\nabla\varphi\cdot F_{\mu(t)}\bigr{\rangle}=\textstyle\iint_{\mathbb{R}^{2}}\nabla\varphi\cdot\bigl{(}\mu(t)\,F_{\mu(t)}\bigr{)}$$

$$\displaystyle\quad=\textstyle\iint_{\Omega_{1}}\nabla\varphi\cdot\bigl{(}\mu(t)\,F_{\mu(t)}\bigr{)}+\iint_{\Omega_{2}}\nabla\varphi\cdot\bigl{(}\mu(t)\,F_{\mu(t)}\bigr{)}$$

$$\displaystyle\quad=\textstyle-\iint_{\Omega_{1}}\varphi\,\textrm{div}\bigl{(}\mu(t)\,F_{\mu(t)}\bigr{)}+\int_{\mathbb{R}}\varphi({\bar{v}},w)\,\mu(t,{\bar{v}}^{-},w)\,F_{\mu(t)}({\bar{v}},w)\cdot\bigl{(}\begin{smallmatrix}1\\ 0\end{smallmatrix}\bigr{)}\,{\textrm{d}}w$$

$$\displaystyle\textstyle\qquad-\iint_{\Omega_{2}}\varphi\,\textrm{div}\bigl{(}\mu(t)\,F_{\mu(t)}\bigr{)}+\int_{\mathbb{R}}\varphi({\bar{v}},w)\,\mu(t,{\bar{v}}^{+},w)\,F_{\mu(t)}({\bar{v}},w)\cdot\bigl{(}\begin{smallmatrix}-1\\ 0\end{smallmatrix}\bigr{)}\,{\textrm{d}}w$$

$$\displaystyle\quad=\textstyle-\iint_{\mathbb{R}^{2}}\varphi\,\textrm{div}\bigl{(}\mu(t)\,F_{\mu(t)}\bigr{)}-\int_{\mathbb{R}}\varphi({\bar{v}},w)\,\bigl{[}\!\!\bigl{[}\mu(t,\cdot,w)\bigr{]}\!\!\bigr{]}\,\mathcal{V}_{\mu(t)}({\bar{v}},w)\,{\textrm{d}}w$$

$$\bigl{[}\!\!\bigl{[}\mu(t,\cdot,w)\bigr{]}\!\!\bigr{]}\stackrel{{\scriptstyle\smash{\scriptscriptstyle\mathrm{def.}}}}{{=}}\mu(t,{\bar{v}}^{+},w)-\mu(t,{\bar{v}}^{-},w)$$

$$\displaystyle\partial_{t}\mu(t)+\partial_{v}\bigl{(}\mu(t)\,\mathcal{V}_{\mu(t)}\bigr{)}+\partial_{w}\bigl{(}\mu(t)\,\mathcal{W}\bigr{)}=-\lambda\,\mu(t)\,,$$

$$\text{with }\mu (t=0)={\mu }_0,on{\Omega }_i,i=1,2,$$

$$\displaystyle\mathcal{V}_{\mu(t)}({\bar{v}},w)\,\bigl{[}\!\!\bigl{[}\mu(t,\cdot,w)\bigr{]}\!\!\bigr{]}=\int_{\mathbb{R}}\lambda(v^{\prime})\,\sigma_{\bar{w}}\mu(t,v^{\prime},w)\,{\textrm{d}}v^{\prime}\,,\ \forall w\in\mathbb{R}\,.$$

$$(v_i^0,w_i^0)\sim {\mu }_0,$$

Full text

{AMS}

68Q25, 28Dxx, 65Cxx, 65Mxx, 65C05, 65C20

1 Introduction

Bridging the gap between microscopic and macroscopic descriptions of biological neural networks is one of the current challenges in neuroscience. Microscopic description of individual neurons has been widely studied, from the historical model of Louis Lapicque [28] to the Nobel winning model of Hodgkin and Huxley [22], reproducing reality in a very satisfying way, both qualitatively and quantitatively. On the other end of the spectrum, electroencephalography [2] and functional magnetic resonance imaging [34] provide insight in brain activity and connectivity at a macroscopic level. But how can we explain the macroscopic measurements from the microscopic description? A direct approach consists of considering a large model coupling single neuron models in a finite size network description. The complexity of this approach increases with the size of the system, leading to a high computational cost for simulations, and a difficult theoretical analysis. An alternative approach is to use the mean-field approximation, where the individual description is homogeneous stemming from exchangeability of the particles. An advantage of the mean-field description is the possibility to use the tools of PDEs theory for theoretical analysis. Also, in this approach, the complexity remains constant, independent of the number of neurons in the system.

The family of two-dimensional (2D) nonlinear spiking neuron model [20, 43, 44] is efficient at reproducing the majority of observed membrane potential behaviours. It is a good compromise between the simple but limited integrate-and-fire model, and the rich but highly complex Hodgkin-Huxley model. However, it has been shown that networks composed of such neurons, where the spiking process is deterministic, might lead to a blow up of the solution in finite time [10]. The model introduced recently in [31, 17], where the spiking event is stochastic, elegantly circumvents this problem. Additionally, stochastic spiking with a probability increasing with the membrane potential value is biologically relevant [12] and has been used [40] to model cortical columns. Yet, the neuron model they introduced in [31] is only 1D, limiting the dynamics to integrator neurons without bursting nor resonance for instance.

The use of 2D nonlinear spiking neurons in a network and the study of the mean-field approximation brings many difficulties which have been tackled only recently [11]. The first difficulty is that spikes are modelled as finite time explosions rather than threshold crossing. Note that although the singular behaviour of [10] might disappear for this neuron model, we focus on the framework of [31]. The second difficulty is that the rate function which governs the spiking mechanism is unbounded which renders the use of classical theory for piecewise deterministic Markov processes [9] more difficult. Note that a related network with scalar neurons, e.g. the quadratic IF neuron, is widely studied in the literature under the name of the Ott-Antonsen ansatz [36, 37, 38].

The study of similar noiseless case networks have been studied in the past [33, 32] under the first order moment closure approximation which effectively makes the system one-dimensional. However, this approximation implies that one loses information regarding the adaptation variable. Nevertheless, the authors were able to compute bifurcation diagrams thanks to this approximation and to predict bursting activity at the network level. Recently, several approximations have been re-visited [1] with similar aims. Note that network bursts have also been predicted recently based on separation of timescales approximation [16].

Please note that this mean-field approach for spiking neural networks has been extensively used in the past [5, 39, 35] and put on rigorous grounds in [7, 10]. The use of generalised integrate and fire neurons in mean-field models has been studied in [19, 31, 17]. The mean-field limits of networks of spiking neurons modelled by Hawkes processes has been studied by [8]. Several extensions have been provided namely by introducing space dependency [14] or dendritic compartments [23, 18].

From the simulation point of view, finite size network dynamics might be efficiently computed using a Monte Carlo approach and taking advantage of graphics processing unit (GPU), whereas conservative mean-field may be simulated using finite volume methods (FVM).

The PDEs obtained as mean-field limits of networks of 2D nonlinear spiking neuron are challenging to simulate (see [26]) and to study theoretically: the velocity field can be explosive (one may think about adaptive exponential integrate and fire models [4]), the reset introduces a spatial singularity in the form of a Dirac measure, and the network activity includes a nonlocal term. Qualitative properties of FVM are crucial in order to reliably approximate such equations (see [30] for a diffusive case), especially when the aim is to capture invariant measures.

In this article, we introduce a 2D nonlinear stochastic spiking neural network combining [43] and [31]. In Section 2, we present the model and heuristically derive its mean-field limit. Section 3 is devoted to numerical methods and their properties. Finally, in Section 4 we present numerical simulations of the model.

2 A detailed stochastic model and its mean-field approximation

2.1 Detailed stochastic model

We consider a population of $N$ identical neurons. Let $v_{i}(t)\in\mathbb{R}$ (resp. $w_{i}(t)\in\mathbb{R}$) denote the membrane potential (resp. the adaptation current) of neuron $i$ at time $t$. The orders of magnitude are typically millisecond for the time, millivolt for membrane potential and pico Ampère for the adaptation current. The adaptation current is not directly linked to a biological quantity, but represents internal processes such as ion channel dynamics or propagation of depolarisation.

The model we consider is a piecewise deterministic Markov process (PDMP). The flow is given by

[TABLE]

where $\tau_{w}$ and $b$ are real parameters, and $I$ corresponds to the external stimuli (input current). Several choices are possible for the nonlinearity $F$, corresponding to different classical models, such as

[TABLE]

Coupled with Eq. 1 and independent of each other, each neuron $i$ spikes at a given rate $\lambda(v^{i})>0$ depending on the membrane potential, and imposes a jump transition given by

[TABLE]

where $\bar{w}>0$. Note that is is possible to adapt, at minor cost, the algorithms that follows to the case $\bar{w}<0$.

The first part of the transition describes the reset following the emission of a spike by an isolated neuron while the second part models the excitatory interaction between neurons. In short, when a neuron spikes, it increases the membrane potentials of the post-synaptic neurons by the amount $\frac{J}{N}$. We focus on the case $J>0$ although $J<0$ can be treated analogously.

We assume that the mechanisms acting on each neuron are independent and identically distributed at initial time according to $\mu_{0}({\textrm{d}}v,{\textrm{d}}w)$; and that the spike events for every neuron $i$ and the initial condition of the system are mutually independent.

This stochastic approach was originally introduced in [20, 31] in the framework of 1D linear spiking neurons. Stochastic firing is biologically justified [12], and avoids a mathematical problem appearing within networks using deterministic firing [10]. In this work, we generalise this approach to the case of 2D nonlinear spiking neurons.

2.2 mean-field approximation

Let $(v^{N}_{i}(t),w^{N}_{i}(t))_{1\leq i\leq N}$ denote the dynamic Eq. 1–Eq. 2 described in the previous section. It can be rewritten as the following interacting particle system:

[TABLE]

where $\mathcal{N}^{i}$ are $N$ independent Poisson random measures with intensity measure ${\textrm{d}}z\times{\textrm{d}}s$ (Lebesgue measure on $\mathbb{R}^{2}_{+}$). We suppose that the initial conditions $(v^{N}_{i,0},w^{N}_{i,0})$ are independent with the same distribution $\mu_{0}$, for all $i$, and that they are also independent from the $\mathcal{N}^{i}$’s. Particles in (3) interact only through the term:

[TABLE]

Suppose that the interacting particle system (3) features a propagation of chaos property [42]; this will have to be demonstrated in a future study. This property, that we do not prove in the present work, roughly means that when $N$ is large, the particles tend to behave like independent particles with the same limit distribution. Hence, according to the law of large numbers, the (stochastic) interaction term $\Psi^{N}{t}$ converges to the following (deterministic) expression:

[TABLE]

as $\mathcal{N}^{1}$ is a Poisson random measure with intensity measure ${\textrm{d}}z\times{\textrm{d}}s$, the limit interaction term is:

[TABLE]

In other words, we have a mean-field limit, indeed the empirical distribution of the interacting particles:

[TABLE]

converges to a (deterministic) distribution $\mu(t,{\textrm{d}}v,{\textrm{d}}w)$ which represents the distribution of a “limit” particle $(v(t),w(t))$ described as:

[TABLE]

where $(v_{0},w_{0})\sim\mu_{0}$ and $\mathcal{N}({\textrm{d}}z,{\textrm{d}}s)$ is a random Poisson measure with intensity measure ${\textrm{d}}z\times{\textrm{d}}s$, $(v_{0},w_{0})$ and $\mathcal{N}({\textrm{d}}z,{\textrm{d}}s)$ are independent.

2.3 Derivation of the nonlinear PDE

The PDMP Eq. 5 is of McKean-Vlasov type as it is not simply an equation for $(v(t),w(t))$ per se as the right hand side of Eq. 5 depends on $(v(t),w(t))$ and on its law through the deterministic current $\psi(t)=\mathbb{E}(\lambda(v(t)))$. The infinitesimal generator of this process is:

[TABLE]

with

[TABLE]

and

[TABLE]

The generator is defined for all test functions $\varphi:\mathbb{R}^{2}\mapsto\mathbb{R}$ in $\mathcal{C}^{1}_{b}(\mathbb{R}^{2})$ (continuously differentiable in $(v,w)$ and bounded). The operator $\mathcal{L}^{\textrm{\tiny t}}_{\mu}$ corresponds to the “ODE” part of the McKean-Vlasov PDMP Eq. 5 and $\mathcal{L}^{\textrm{\tiny j}}$ to the “pure jump” part.

The evolution of the distribution $\mu(t)$ is given by the Kolmogorov forward equation as a weak PDE:

[TABLE]

for all test functions, where $\left\langle\mu,\varphi\right\rangle\stackrel{{\scriptstyle\smash{\scriptscriptstyle\mathrm{def.}}}}{{=}}\iint_{\mathbb{R}^{2}}\varphi(v,w)\,\mu({\textrm{d}}v,{\textrm{d}}w)$ is the usual duality bracket. Note that the adjoint $\mathcal{L}_{\mu}^{*}$ of $\mathcal{L}_{\mu}$ is given by:

[TABLE]

with $\sigma_{\bar{w}}\mu(t,{\textrm{d}}v,{\textrm{d}}w)\stackrel{{\scriptstyle\smash{\scriptscriptstyle\mathrm{def.}}}}{{=}}\mu(t,{\textrm{d}}v,{\textrm{d}}w-{\bar{w}})$, so that Eq. 8 is a weak form of the following strong form of the forward Kolmogorov equation:

[TABLE]

which reads:

[TABLE]

where the shift operator $\sigma_{\bar{w}}$ is now $\sigma_{\bar{w}}\mu(t,v,w)\stackrel{{\scriptstyle\smash{\scriptscriptstyle\mathrm{def.}}}}{{=}}\mu(t,v,w-{\bar{w}})$.

2.4 The mean-field as coupled transport equations

The PDE Eq. 10 can be represented as a transport PDE with a jump condition at the interface:

[TABLE]

Indeed, let:

[TABLE]

In Eq. 8 consider the term:

[TABLE]

where

[TABLE]

denotes the jump of $\mu(t,v,w)$ through the interface $\Gamma$. Using test functions $\varphi$ with support included in $\Omega_{1}\cup\Omega_{2}$ we get:

[TABLE]

where F_{\mu(t)}=\bigl{(}\mathcal{V}_{\mu(t)}\,,\,\mathcal{W}\bigr{)}; and with test functions $\varphi$ with support included in $\Gamma$ we get the interface condition:

[TABLE]

Hence the PDE Eq. 11 with the interface condition Eq. 12 is equivalent to Eq. 10. Note that the PDE Eq. 11 can be solved using the method of characteristics and that the boundary conditions are specified if and only if the characteristics enter the domain.

3 Numerical methods

3.1 Monte Carlo simulation of the particle system

We now present a Monte Carlo simulation procedure of the interacting particles system introduced in Section 2.1 and rewritten as Eq. 3. The Monte Carlo procedure is exact up to the time discretisation of the trajectories.

For the finite size network description, we take advantage of the all-to-all connectivity of the network to design a simple Monte Carlo algorithm. Let $v_{i}^{n}$ (resp. $w_{i}^{n}$) denote the membrane potential (resp. the adaptation current) of neuron $i$ at time $t^{n}$. After initialisation on each neuron $i$:

[TABLE]

the method consists of three steps. First, we use an Euler step for the ODE integration to predict the deterministic evolution on each neuron $i$:

[TABLE]

Second, we assume that the jump rate $\lambda$ does not evolve much in between the jump times. Hence, neuron $i$ fires with probability $\lambda(v_{i}^{*})\,\Delta t$:

[TABLE]

independently on each neuron $i$. Third, the sum of all the spikes:

[TABLE]

is distributed to each neuron separately:

[TABLE]

This method is stable under certain conditions. For the deterministic part of the PDMP, the classical Courant–Friedrichs–Lewy condition for forward Euler must hold. For the probabilistic part, the time step has to be small enough so that, for each neuron, at most one jump occurs during the time step. This depends on the parameters of the model. In practice, we use the time step for the flow part and plot at the same time as the jump frequency, and then reduce the time step if necessary.

The exact Monte Carlo methods that are developing more and more in the field of neuroscience [45, 29]. Here, in the case of a simple Euler scheme for the flow approximation, the proposed “almost” exact Monte Carlo method is very effectively simulated on graphics processing units (GPU).

3.2 Simulation of the mean-field PDE Eq. 10

As this PDE is defined on $\mathbb{R}^{2}$, a preliminary step is to consider a compact subset:

[TABLE]

large enough to contain the initial condition and the solution during the evolution.

The numerical method for the simulation of the macroscopic description Eq. 10 is based on a splitting strategy [41]: we solve alternatively the transport part, using a conservative finite volume method (FVM, see [47]), and the source part, using the approximate solution of an ODE system. We define a regular $N_{v}\times N_{w}$ grid on $\Omega$:

[TABLE]

with

[TABLE]

with $\Delta v>0$ and $\Delta w>0$ being the mesh size in the $v$ and $w$ direction respectively, and suppose $v_{\max}=v_{N_{v}}$, $w_{\max}=w_{N_{w}}$. Define also $v_{i\pm 1/2}=v_{i}\pm\Delta v/2$, $w_{j\pm 1/2}=v_{j}\pm\Delta w/2$, and:

[TABLE]

Cells are defined by

[TABLE]

for $1\leq i\leq N_{v}$, $1\leq j\leq N_{w}$.

Let $\Delta t>0$ denote the time step, and $t^{n}=n\,\Delta t$. We now propose a finite volume approximation $\mu^{n}=(\mu^{n}_{i,j})_{1\leq i\leq N_{v},1\leq j\leq N_{w}}$ of $\mu(t)$ at time $t^{n}$:

[TABLE]

with initial condition $\mu^{0}_{i,j}\simeq\left\langle\mu_{0},{\mathrm{\mathbf{missing}}1}_{\Omega_{i,j}}\right\rangle/(\Delta v\,\Delta w)$.

For the update $\mu^{n}$ to $\mu^{n+1}$, we adopt a splitting technique for Eq. 10 that will be detailed later. This technique alternates numerical approximations of the “transport” part:

[TABLE]

and of the “jump” part:

[TABLE]

We first detail the approximation of these two PDEs and then the splitting technique.

Discretisation of transport part

We describe the update $\mu\to\mu^{\prime}$ which corresponds to the numerical approximation of Eq. 13 on the interval $[0,\Delta t]$ with initial condition $\mu$ and final value at $\Delta t$ assigned to $\mu^{\prime}$. We adopt an upwind finite volume scheme on the structured mesh $\Omega^{\textrm{\tiny d}}$ for the space and a semi-implicit Euler scheme for the time:

[TABLE]

$\mathbf{D}(\mu)\,\mu^{\prime}$ is the approximation of divergence operator $\partial_{v}(\mu^{\prime}\,\mathcal{V}_{\mu})+\partial_{w}(\mu^{\prime}\,\mathcal{W})$ given by

[TABLE]

based on the upwind numerical fluxes defined by

[TABLE]

(see [47]) and the approximated vector field, using numerical integration, by:

[TABLE]

Transition (15) reads:

[TABLE]

with

[TABLE]

Note that $\mathbf{A}_{\Delta t}(\mu)$ and $\mathbf{D}(\mu)$ can be considered as $(N_{v}\times N_{w})^{2}$ square (sparse) matrices, but to avoid notational complexities, we will consider them as linear operators on the set of real functions defined on $\Omega^{\textrm{\tiny d}}$.

Concerning the boundary conditions, we impose null fluxes on the boundaries, leading to

[TABLE]

This semi-implicit approach allows to avoid the restrictive Courant–Friedrichs–Lewy condition that would appear due to the strong nonlinearity $F$ (possibly exponential). Also, the matrix of the linear system to solve at each iteration is sparse, because the nonlocal term is treated explicitly.

Remark 1**.**

*It is easy to see that $\mathbf{D}$ is composed of a matrix diagonal by bands with 5 bands at $-N_{v},-1,0,1,N_{v}$. This is a direct consequence of the formulas Eq. 16 based on a classical stencil scheme: each cell is (possibly) connected to its 4 direct neighbours plus itself. The same holds for $\mathbf{D}^{*}$. *

Discretisation of the jump part

We describe the update $\mu\to\mu^{\prime}$ that corresponds to the numerical approximation of Eq. 14 on the interval $[0,\Delta t]$ with initial condition $\mu$ and final value at $\Delta t$ assigned to $\mu^{\prime}$. Equation Eq. 14 is the forward Kolmogorov equation of a pure jump Markov process with infinitesimal generator $\mathcal{L}^{\textrm{\tiny j}}$. First, we restrict the latter operator to $\Omega$. The most evident way to achieve this is to consider: \lambda(v)\,\bigl{(}\varphi(\bar{v},(\bar{w}+w)\wedge w_{\max})-\varphi(v,w)\bigr{)}, which amounts to an accumulation point at $(\bar{v},w_{\max})$. This operator is then approximated by a jump process on the grid $\Omega^{\textrm{\tiny d}}$ by considering the generator:

[TABLE]

The associated forward Kolmogorov equation is:

[TABLE]

where $\tilde{\mathcal{L}}^{\textrm{\tiny j}}_{(i,j),(i^{\prime},j^{\prime})}=\tilde{\mathcal{L}}^{\textrm{\tiny j}}\phi(i,j)$ with $\phi(i,j)=\delta_{ii^{\prime}}\,\delta_{jj^{\prime}}$. Straightforward calculations lead to the following system of ODEs:

[TABLE]

Note that this last equation is:

[TABLE]

which indeed corresponds to the fact that $(\bar{i},N_{w})$ is an accumulation point.

Hence, the update $\mu\to\mu^{\prime}$ on interval $[0,\Delta t]$ consists of solving (20) with initial condition $\mu$ and to set $\mu^{\prime}=\mu(\Delta t)$. System (20) can be explicitly solved but not in a convenient way. As the previous step of the splitting method is of the first order, it is consistent to propose a simple first order approximation of (20).

Note that the solutions of the first two sets of equations in (20) are $\mu_{i^{\prime},j^{\prime}}(t)=e^{-\lambda(v_{i^{\prime}})\,t}\,\mu_{i^{\prime},j^{\prime}}(0)$ for all $(i^{\prime},j^{\prime})$ not in $\{\bar{i}\}\times\{\bar{j}+1,\dots,N_{w}\}$. For the other components of (20), we choose to make an approximation which features the same time Euler scheme as (20) and which respects both positivity and mass conservation properties, namely:

[TABLE]

defined by:

[TABLE]

Operator $\mathbf{B}_{\Delta t}$ is first order accurate in time. Indeed, the first order of the Taylor expansion of $\mathbf{B}_{\Delta t}\mu$ in $\Delta t$ corresponds to the Euler scheme for Eq. 20.

Remark 2**.**

*The case $\bar{w}<0$ can easily be adapted. *

Time step strategy

For the iteration $\mu^{n}\rightarrow\mu^{n+1}$ we use a symmetric Strang splitting method of the second order [21, p. 82]:

[TABLE]

coupled with an adaptation of the classical Euler-Richardson extrapolation strategy for the control the time step size [15, 48] (see Algorithm 1).

The choice of the indicator for the evaluation of the evolution of the solution has been made in order to accurately approximate the mean membrane potential $\left\langle\mu(t),\phi_{1}\right\rangle$ with $\phi_{1}(v,w)=v$. It is indeed a good measure of error as the membrane potential of each individual neuron explodes in finite time. Another reasonable choice would be to control the mean firing rate $\left\langle\mu(t),\lambda\right\rangle$. However, for problems with localised (in time) activity such as in Figure 3 Bottom, the indicator would be small during small network spiking activity independently of the underlying dynamics, that is why we do not use this.

Note that in contrast to the classic Euler-Richardson algorithm, we do not write:

[TABLE]

This solution gives a second order accuracy for a fixed time step, but the explicit formulation would break the positivity property of our scheme.

Remark 3**.**

*We did not prove the convergence of our approximate solution to the solution of the PDE Eq. 10 as $N_{v},N_{w}\to\infty$. This requires at least to prove that the PDE Eq. 10 is well posed and to provide some properties regarding its dynamics (e.g. a priori bounds…). Nevertheless, we can readily see that our scheme above is incomplete. We need to set $v_{\max},w_{\max}$ and $v_{\min},w_{\min}$ such that our approximation of the original jump process by the absorbing one does not affect the dynamics too much. The accumulating point $(\bar{i},w_{\max})$ (resp. $(\bar{i},w_{\min})$) in the case $\bar{w}>0$ (resp. $\bar{w}<0$) should be compensated by a redistribution of the mass by the drift if one wants to avoid concentration of mass that is not a feature of PDE Eq. 10. A simple way to achieve this is, for example, to chose $w_{\max}$ high enough so that $(\bar{i},w_{\max})$ is above the w-nullcline and the v-nullcline. This way, the flow is downward to the left at $(\bar{i},w_{\max})$ and any mass at this point will be re-injected into the dynamics. Then, we chose $v_{min}$ small enough, on the left of the $v$-nullcline, to ensure that the vector field is entrant in the domain. Finally, we chose $v_{max}$ large enough, on the right of the $v$-nullcline to capture the explosive behaviour. *

3.3 Properties

Define the mass of a discrete solution $\mu$ as

[TABLE]

We have designed the numerical schemes so that $\mathbf{A}_{\Delta t}(\mu^{n})$ and $\mathbf{B}_{\Delta t}$ are mass conservative. It follows that the general Algorithm 1 is mass conservative as well. Let us now focus on proving the positivity of the algorithm. Note that the same property was proved for a general mesh in [3] using a different method. It should be straightforward to adapt our proof to a general (non-regular) mesh. The main idea of our proof is to note that the adjoint of the discrete divergence $\mathbf{D}^{*}$ is diagonally dominant.

Theorem 3.1**.**

For $\Delta t>0$ and $\mu^{0}\geq 0$, we have the following properties:

there exists a unique solution $\nu$ to the equation:

[TABLE]

with null fluxes on the boundaries Eq. 18, 2. 2.

if $\mu$ is non-negative, then so is $\nu$.

Proof 3.2**.**

For simplicity, as we work at time fixed in this proof, we drop the dependency of $\mathbf{A}$ and $\mathbf{D}$ on $\mu^{0}$. The idea of the proof is to show that $I+\Delta t\,\mathbf{D}^{*}$ is diagonally strictly dominant hence invertible with positive diagonal and non-negative off-diagonal elements. To this end, we first have to identify the adjoint $\mathbf{D}^{*}$ of $\mathbf{D}$. For this, we consider the following quantity

[TABLE]

Using the boundary conditions Eq. 18, we get:

[TABLE]

Note that from Eq. 16, for any $\psi\in\mathbb{R}^{N_{v}-1}$:

[TABLE]

hence:

[TABLE]

Similarly:

[TABLE]

From Eq. 24 and Eq. 25, the diagonal element of $\mathbf{D}^{*}$ are non-negative, the off-diagonal elements are non-positive. Finally, the matrix $\mathbf{D}^{*}$ is diagonally dominant (not strictly). Indeed, each term in Eq. 24–Eq. 25 satisfies this property and the set of matrices satisfying this property is obviously convex.

*Let us now conclude. First $M=I+\Delta t\,\mathbf{D}^{*}$ is diagonally strictly dominant hence invertible. Second $M$ has a non-negative (resp. non-positive) diagonal (resp. off-diagonal), so up to a scaling, we can assume that $\max_{i}M_{ii}<1$, so that $P\stackrel{{\scriptstyle\smash{\scriptscriptstyle\mathrm{def.}}}}{{=}}I-M$ is non-negative and $\left\|P\right\|_{\infty}<1$. We can thus expand $M^{-1}=(I-P)^{-1}$ as an infinite sum $I+P+P^{2}+\cdots$ which is non-negative too. The proof is complete. *

Lemma 3.3**.**

$\mathbf{B}_{\Delta t}$* is non-negative in the sense that for all non-negative $\mu$, $\mathbf{B}_{\Delta t}\,\mu$ is also non-negative. *

Proof 3.4**.**

*From the expression Eq. 22 of $\mathbf{B}_{\Delta t}$, this is a consequence of $\lambda\,\Delta t$ being non-negative. *

Proposition 3.5**.**

*The full numerical scheme based on the splitting strategy consisting of formulas Eq. 17 and Eq. 22 is non-negative. This is the discrete version of the fact that the solution of Eq. 8 is a positive measure. *

Proof 3.6**.**

*This result is a consequence of the previous results on non-negativity of $\mathbf{A}_{\Delta t}(\mu)$ and $\mathbf{B}_{\Delta t}$. *

4 Numerical simulations

Let us define useful statistical quantifiers in order to analyse the results, mean membrane potential for both descriptions:

[TABLE]

The implementation is done in Julia language. For the Algorithm 1, the linear system associated to $\mathbf{A}_{\Delta t}(\mu^{n})$ is encoded in a sparse matrix at each iteration and a general linear solver is then called. GPU simulations of the finite size network were carried on a Nvidia Tesla K80 card. In all simulations, we used a Mersenne twister for the generation of random numbers.

4.1 Convergence and propagation of chaos

In this paragraph, we study the behaviour of the finite size network Eq. 1 as the number $N$ of neurons is increased. We use the default parameters (see Table 1) and we start from a Gaussian distribution:

[TABLE]

with $\sigma_{1}=\sigma_{2}=1$, $\mu_{1}=-1.3$, $\mu_{2}=2.28$.

On Fig. 1 (top) we present the results of a convergence test. We plot the means membrane potential $V_{N},V_{\infty}$ computed using two distinct approaches, microscopic and macroscopic. For the microscopic approach, we have used two different methods. The label PDMP1 refers to the method developed in this article whereas the label PDMP2 refers to the Julia library PiecewiseDeterministicMarkovProcesses.jl, based on a stiff solver for the flow and a specific algorithm. This second method, based on [46], is almost exact in that the errors come from the numerical flow in between jumps and the random generator. However, PDMP2 cannot be used for a network size larger than $\sim 10^{4}$ and that is why we use it as a reference for PDMP1. The three curves are in very good agreement. We note that the convergence is linear in $N$ towards the mean-field approximation, and that a plateau is reached for $N=10^{6}$ (machine zero is reached).

We say that a $N$-neurons system, such as Eq. 1, propagates chaos if the particles that compose the system become independent as the total number tends to infinity. This concept was originally introduced in [25] in the framework of kinetic theory. In order to numerically support this property, we adapt the algorithm developed in [6]. For each size of network $N$, we consider $M$ realisations. Each realisation consists of simulating the system Eq. 1 until a final time $T$, long enough such that the system reaches a steady state. Then, we randomly pick (uniform distribution) two neurons, let’s say neuron $i$ and $j$, and we keep them. At the end, we plot the 2D histogram in order to analyse the correlations. On Fig. 1 (bottom) we present the results of this test, with $N=10^{4}:10^{6}$ and $M=400$. A correlation is clearly visible for the case $N=10^{4}$, and seems to vanish as $N$ grows, as the cases $N=10^{5}$ seem to indicate. This would be a numerical argument in favour of propagation of chaos in the case of this model.

4.2 Equilibrium and bifurcation

The invariant distributions $\mu^{\textrm{\tiny inv}}$ of (10) are invariant distribution(s) of an isolated neuron, e.g. for $J=0$, but for the current $I+J\left\langle\mu^{\textrm{\tiny inv}},\lambda\right\rangle$. We thus have to look for the invariant distributions of an isolated neuron in order to study their existence in the nonlinear case. Two examples are shown in Figure 2 Left, in the case of two and no equilibria for the underlying vector field. These results are obtained by using the Algorithm 1 for very long times. It appears that the implementation respects the properties of positivity and conservation of mass. In the first case, one equilibrium is an attracting focus to which the dynamics are attracted. The isolated neuron is a PDMP with embedded jump chain $(\bar{v},w_{n})_{n}$ being a Markov chain. In the Right column, we compare the invariant distribution of the embedded chain with the (renormalised) quantity $\int_{\mathbb{R}}\lambda(v)\,\sigma_{\bar{w}}\,\mu^{\textrm{\tiny inv}}({\textrm{d}}v,\cdot)$. In the second case, the density presents several peaks on the reset line below the $v$ nullcline. This represents several consecutive spikes before a quiescent period during which the dynamics re-accumulate under the $v$-nullcline: this is typically a bursting behaviour.

Based on extensive simulations, we make the following conjecture in the case $J=0$.

*The isolated neuron with rate function $\lambda=\exp$ and $F(v)=e^{v}-v$ has a unique invariant distribution to which the (linear) dynamics converge exponentially fast. Furthermore, this invariant distribution has a density with respect to the Lebesgue measure. *

It seems that the conjecture should hold for $\lambda$ positive increasing and $F$ such that $\lambda/F$ is not integrable111to ensure that the network spikes at least once at $v=+\infty$, but we don’t have enough data to support this.

The existence of invariant distributions for the nonlinear equation (10) is much less trivial to study numerically. Based on the above conjecture, it should hold that there is a unique attracting invariant distribution when the connectivity $J$ is small.

We therefore seek to study the effect of connectivity strength by varying the coupling parameter $J$. Following the approach done in [13], we look for synchronised activity within the network. We use same initial condition and parameters as in the previous paragraph, except for the parameter $J$.

Results are displayed on Fig. 3. For values of $J$ smaller than $6.15$ (inset of Fig. 3 Top Left), the network is not synchronised, and tends to relax to an invariant distribution. Above a threshold value $J^{*}\approx 6.15$, we observe synchronised activity within the network characterised by a periodic solution to (10); see Fig. 4 for an example of a periodic orbit far from the “Hopf bifurcation point”. Figure 3 is a numerical evidence of the possible existence of a Hopf bifurcation for the network, based on the parameter $J$. Please note that we could not get the scaling behaviour [27] in $\propto\sqrt{J-J^{*}}$ for the amplitude of the periodic orbit (inset of Fig. 3 Top Left) possibly due to a subcritical bifurcation.

5 Conclusion

In this work, we have presented a new nonlinear stochastic model of a network of stochastic spiking neurons. This model naturally avoids blow-up solutions that may appear using “threshold crossing” based spiking like for the Integrate and Fire neuron model. The system admits a mean-field limit: we have derived the PDE heuristically, and numerical simulations have confirmed the heuristic. We have shown that this mean-field may be seen as a coupled transport equation, and could be entirely defined by deduction from the solution on the reset interface.

We have designed a Monte Carlo method to simulate the system of neurons. On the other hand, we have designed a reliable finite volume method in order to simulate the PDE. This numerical method is overall conservative, theoretically second order accurate in time and preserves positivity. As further work, it would be interesting to use mesh adaptation.

Taking advantage of those numerical tools, we have studied numerically the network. As the network size tends to infinity, numerical simulations tend to show propagation of chaos, and convergent behaviour from microscopic description Eq. 1 to macroscopic one Eq. 10. Two open theoretical questions remain on this point: does the particle system propagate chaos, and is it possible to rigorously prove the convergence from the particle system to the mean-field model?

By varying the strength of the connectivity, we have observed a Hopf bifurcation, signature of a synchronisation of the activity within the system. This work may be pursued in several ways: prove the existence and uniqueness of invariant distributions, and the theoretical and numerical study of bifurcations.

Acknowledgements

This research has received funding from the European Union’s Horizon 2020 Framework Programme for Research and Innovation under the Specific Grant Agreement No. 785907 (Human Brain Project SGA2).

We would like to thank the Julia community for their help online. We would like to thank Laurent Monasse for fruitful discussions.