1.1. Finite particle system

Let us consider a network consisting of K large populations of neurons, where the number of neurons in the kth population is denoted by $N_k$ and the total number of neurons in the network is $N = N_1 + \dots + N_K$ . Let $Z^{k, n}_t$ represent the number of spikes of the nth neuron belonging to the kth population during the time interval [0, t]. The family of (simple) counting processes $\{(Z^{k, n}_t)_{t\geq 0} ,\, 1\leq k \leq K, \, 1\leq n \leq N_k \}$ is characterized by the intensity processes $(\lambda^{k, n}(t))_{t\geq 0}$ , which are defined through the relation

\begin{equation*}\mathbb{P}\big(Z^{k, n}_t \text{ has a jump in } (t, t+ dt]|\mathcal{F}_t\big) = \lambda^{k, n}(t) dt,\end{equation*}

where $\mathcal{F}_t$ contains the information about the processes $(Z^{k, n}_t)_{t\geq 0}$ up to time t. The interested reader is referred to [Reference Delattre, Fournier and Hoffmann14] for the existence and formal definition of these processes. The mean-field framework considered here corresponds to intensities $\lambda^{k, n}(t)$ given by

(1) \begin{equation}\lambda^{k, n}(t) = f_k\left(\sum_{l=1}^K \frac{1}{N_l} \sum_{1\leq m \leq N_l} \int_{(0, t)} h_{kl}(t-s) dZ^{l, m}_s \right),\end{equation}

where $\{h_{kl}\;:\;\mathbb{R}_+ \to \mathbb{R} \}$ is a family of synaptic weight functions, which model the influence of population l on population k. When $h_{kl}$ is positive, the influence is excitatory. Conversely, when $h_{kl}$ is negative, the influence is inhibitory. The function $f_k\;:\;\mathbb{R} \to \mathbb{R}_+ $ is the spiking rate function of population k. The expression ‘mean-field framework’ refers to the fact that the intensity $\lambda^{k, n}(t)$ depends on the whole system only through the ‘mean-field’ behaviour of each population, namely $\frac{1}{N_l} \sum_{1\leq m \leq N_l} dZ^{l, m}_s$ . Furthermore, letting $N\to \infty$ (while keeping K fixed), we assume that $N_k/N\to p_k>0$ for all k.

Throughout the paper we assume that the functions $f_k$ satisfy the following conditions:

(A)The spiking rate functions $f_k$ are strictly positive, Lipschitz-continuous, non-decreasing, and such that $0< f_k \leq f^{\max}_k \text{ for } k = 1,\dots,K$ , where $f^{\max}_k$ is a finite constant.

In this paper, Erlang-type memory kernels and a cyclic feedback system of interactions are considered. This means that for each k, population k is influenced only by population $k+1$ , where we identify $K+1$ with 1. In this case, all the memory kernels are null except the ones given by

(2) \begin{equation}h_{kk+1}(t) = c_k e^{-\nu_k t}\frac{t^{\eta_k}}{\eta_k!},\end{equation}

where $c_k = \pm 1$ . This constant determines whether the population has an inhibitory ( $c_k=-1$ ) or excitatory ( $c_k=+1$ ) effect. The parameter $\eta_k\geq 1$ is an integer determining the memory order for the interaction function from population $k+1$ to population k.

The parameters $\eta_k$ and $\nu_k$ determine, intuitively, the typical delay of interaction and its time width. The delay of the influence of population $k+1$ on population k attains its maximum $\eta_{k+1}/\nu_{k+1}$ units back in time, and its mean is $(\eta_{k+1}+1)/\nu_{k+1}$ . The larger this ratio is, the more important ‘old’ events are. When the ratio is fixed (equal to $\tau$ ), but both $\eta_k$ and $\nu_k$ tend to infinity, $h_{kk+1}$ tends to a Dirac mass in $\tau$ . This means that only one specific moment in time is important. The interested reader is referred to [Reference Ditlevsen and Löcherbach15] and [Reference Löcherbach27] for more details.

In this paper we are interested in the processes $\{(\bar{X}^{k,1}_t)_{t\geq 0},\, 1\leq k\leq K\}$ , which are the arguments of the function $f_k$ in Equation (1) and are defined by

(3) \begin{equation}\bar{X}^{k,1}_t = \frac{1}{N_{k+1}} \sum_{1\leq m \leq N_{k+1}} \int_{(0, t)} h_{kk+1}(t-s) dZ^{k+1, m}_s.\end{equation}

When the memory kernels are given in the form (2), the processes defined in (3) can be obtained as marginals of the process $(\bar{X}_t)_{t\geq 0} = \{ (\bar{X}^{k,j}_t)_{t\geq 0},\, 1\leq k\leq K, 1\leq j\leq \eta_k+1\}$ which solves the following system of dimension $\kappa = \sum_{k=1}^K (\eta_k+1)$ :

(4) \begin{equation}\begin{cases}d\bar X^{k,j}_t = \left[-\nu_k \bar X^{k,j}_t+\bar X^{k,j+1}_t\right]dt \quad \text{for } j= 1,\dots,\eta_k, \\[5pt] d\bar X^{k,\eta_k+1}_t = -\nu_k \bar X^{k,\eta_k+1}_tdt+c_k d\bar{Z}^{k+1}_t,\\[5pt] \bar{X}_0 = x_0 \in \mathbb{R}^\kappa,\end{cases}\end{equation}

where

\begin{equation*}\bar{Z}^{k+1}_t = \frac{1}{N_{k+1}}\sum_{n=1}^{N_{k+1}} Z^{k+1,n}_t,\end{equation*}

with each $Z^{k+1,n}_t $ jumping at rate $f(\bar{X}^{k+1,1}_{t-})$ ; see [Reference Ditlevsen and Löcherbach15] for more insight. This type of equation is sometimes called a Markovian cascade in the literature [Reference Ditlevsen and Löcherbach15]. Throughout the paper we refer to the solution of (4) as the PDMP.

The process $(\bar{X}_t)_{t\geq 0}$ summarizes and averages the influence of the past events. This process, along with the firing rate functions $f_k$ , determines the dynamics of $(Z^{k,n}_t)_{t\geq 0}$ , described by its intensity (1).

From a modelling point of view, the process $(\bar{X}^{k,1}_t)_{t\geq 0}$ can be roughly regarded as the voltage membrane potential of any neuron in population k. Then, the probability that a neuron will emit a spike is given as a function of its membrane potential. To summarize, the processes with coordinates (k, 1) defined by (3) describe the membrane potential in each population, whereas the other coordinates represent higher levels of memory for the process.

Note that the model presented so far starts with empty memory. The right-hand side of (1) is equal to $f_k(0)$ at time $t=0$ or equivalently $x_0=0_\kappa$ . However, one could easily generalize this to any initial condition $x_0$ in $\mathbb{R}^\kappa$ , as is done in the rest of the paper. Moreover, the interested reader is referred to [Reference Duarte, Löcherbach and Ost16], where a more general model is studied numerically and theoretically for $K=1$ population.

1.2. Notation

Now we focus on the case of two interacting populations of neurons ( $K = 2$ ), consisting of $N_1$ and $N_2$ neurons, respectively. Taking $K=2$ allows for an investigation of the interactions between populations of different sizes while avoiding heavy notation. Throughout the paper the following notation is used: $\mathbb{O}_{n \times m}$ denotes an $(n\times m)$ -dimensional zero matrix, and $0_n$ denotes an n-dimensional zero vector. It is convenient to rewrite the system (4) in the matrix–vector form

(5) \begin{equation}d\bar{X}_t = A \bar{X}_t dt + \Gamma \:d\bar{Z}_t,\qquad \bar{X}_0 = x_0\in \mathbb{R}^\kappa,\end{equation}

with

• • $A \in \mathbb{R}^{\kappa\times\kappa}$ defined as

  (6) \begin{equation}A = \left(\begin{matrix}A_{\nu_1} & \mathbb{O}_{(\eta_1+1) \times (\eta_2+1)} \\[5pt]\mathbb{O}_{(\eta_2+1) \times (\eta_1+1)} & A_{\nu_2}\end{matrix}\right),\end{equation}
  where $A_{\nu_k}$ is an $(\eta_k+1)\times (\eta_k+1)$ tri-diagonal matrix with lower diagonal equal to $0_{\eta_k}$ , diagonal equal to $(\!-\!\nu_k,\dots, -\nu_k)$ and upper diagonal equal to $(1,\dots, 1)$ ,

• • $\Gamma \in \mathbb{R}^{\kappa\times 2}$ having zero coefficients everywhere, except for $\Gamma_{\eta_1+1,2} = c_1$ and $\Gamma_{\kappa,1} = c_2$ ,

• • and $ \bar{Z}_t = \left( \bar{Z}^1_t, \bar{Z}^2_t \right)^T$ .

Throughout the paper the following convention is made. The coordinates of a generic vector x in $\mathbb{R}^\kappa$ are denoted by either $(x_i)_{i=1,\dots,\kappa}$ or $(x^{k,j})_{k=1,2;\, j=1,\dots,\eta_k+1}$ , with the relations $i = j$ if $k=1$ and $i = \eta_1+1+j$ if $k=2$ . The second notation is usually preferred since each population is easily identified by the index k. For some generic function $g\;:\;\mathbb{R}^\kappa \to \mathbb{R}^\kappa$ , the upper indices are used as follows: $\left( g(x)\right)^{k,j}$ . Moreover, it is sometimes more natural to consider some generic $\mathbb{R}^\kappa$ -valued process $x_t$ population-wise. Thus, it is split into two components $x^1_t = (x^{1,1}_t,\dots,x^{1,\eta_1+1}_t)\in \mathbb{R}^{\eta_1+1} $ and $x^2_t = (x^{2,1}_t,\dots,x^{2,\eta_2+1}_t)\in \mathbb{R}^{\eta_2+1}$ , such that $x_t = (x^1_t, x^2_t)^T \in \mathbb{R}^{\kappa}$ .

2.1. Strong error bound between the limiting stochastic diffusion and the PDMP

Any error bound of the diffusion approximation is determined by two facts, namely the approximation of a compensated Poisson process by a Brownian motion and the approximation of $N_k$ by $p_kN$ . We get rid of the second approximation by considering the SDE (7) with $p_k = N_k/N$ and denote the solution of this equation by Y. By choosing different notation we stress the fact that, in contrast to X, the solution Y depends on the exact number of neurons $N_k$ and not on its proportion, obtained in the mean-field limit. The same convention is used in [Reference Ditlevsen and Löcherbach15], where the following weak error bound is proved.

Theorem 1. ([Reference Ditlevsen and Löcherbach15], Theorem 4.) Grant Assumption (A), and suppose that all spiking functions $f_k$ belong to the space $C^5_b$ of bounded functions having bounded derivatives up to order 5. Then there exists a constant C, depending only on $f_1$ , $f_2$ and the bounds on their derivatives, such that for all $ \varphi \in C^4_b(\mathbb{R}^\kappa, \mathbb{R})$ and $t\geq 0$ ,

(11) \begin{equation}\sup_{x\in \mathbb{R}^k} \left| \mathbb{E}_x\varphi(\bar{X}_t) - \mathbb{E}_x\varphi(Y_t) \right| \leq Ct \frac{\|\varphi \|_{4,\infty}}{N^2},\end{equation}

where $\mathbb{E}_x$ denotes the conditional expectation given that $\bar{X}_0=Y_0=x$ .

In the following, we strengthen the above result, allowing for a comparison of trajectories of the PDMP and the stochastic diffusion.

Theorem 2. (Strong error bound.) Grant Assumption (A), and let $||\cdot||_{\infty}$ denote the sup norm on $\mathbb{R}^\kappa$ . For all $N>0$ , a solution $\bar{X}$ of (5) and a solution Y of (7) (with $p_k=N_k/N$ ) can be constructed on the same probability space, such that there exists a constant $C>0$ such that, for all $T>0$ ,

(12) \begin{equation}\sup_{t\leq T} \| \bar{X}_t - Y_t \|_{\infty} \leq \Theta_N e^{CT} \frac{\log(N)}{N}\end{equation}

almost surely, where $\Theta_N$ is a random variable with exponential moments whose distribution does not depend on N. In particular,

(13) \begin{equation}\mathbb{E}\left[ \sup_{t\leq T} \| \bar{X}_t - Y_t \|_{\infty} \right] \leq C e^{CT} \frac{\log(N)}{N}.\end{equation}

The proof of Theorem 2 is mainly inspired by [Reference Kurtz23] and relies on two main ingredients, a strong coupling between the standard Poisson process and the Brownian motion and a sharp result on the modulus of continuity for the Brownian motion. All the material is postponed to Appendix A.1.

When comparing (11) and (13), one notices that there is an exchange between the expectation sign and the absolute value. There are two prices to pay for such an exchange: first, a slower convergence rate with respect to N; second, a faster divergence rate with respect to t (the exponential term comes from a Grönwall-type argument). In the following remark we make precise the bound on the error which is caused by directly using the parameter $p_k$ instead of $N_k/N$ .

Remark 1. Let Y denote a solution of (7) (with parameter $p_k$ equal to $N_k/N$ ), and let X denote a solution of (7) with fixed values $p_k$ . Following the proof of Theorem 2, one can show that

\begin{equation*}\sup_{t\leq T} \| X_t - Y_t \| \leq \Theta_N e^{CT}\left( \frac{\log(N)}{N} + \max_k\left\{ \frac{1}{\sqrt{p_k N}} \left( 1 - \sqrt{p_k N/N_k}\right) \right\}\right),\end{equation*}

so that the strong error bound stated in the theorem also holds for the non-modified SDE if $\sqrt{p_k N/N_k} - 1$ is of order $N^{-1/2}$ or of faster order.

Fortunately, for any fixed N, setting $N_1 = \lfloor p_1N \rfloor$ and $N_2 = \lceil p_2 N \rceil$ ensures that $\sqrt{p_k N/N_k} - 1$ is of order $N^{-1} < N^{-1/2}$ , which grants that Theorem 2 holds for the SDE (7).

Since the SDE (7) transforms into the ODE (10) as N goes to infinity, the strong error bound can be used to prove the convergence of the PDMP to the solution of the ODE. However, this is beyond the scope of this paper.

2.2. Properties of the stochastic diffusion

The solution process $(X_t)_{t\geq 0}$ of SDE (7) is positive Harris recurrent with invariant measure $\pi$ which is of full support (see [Reference Löcherbach27]). This means that the trajectories visit all sets in the support of the invariant measure infinitely often almost surely. More precisely, for any initial condition $x_0$ and measurable set A such that $\pi(A)>0$ , $\limsup_{t\to +\infty} \mathbf{1}_{A}(X_t) = 1$ almost surely. Furthermore, by following the arguments in [Reference Mattingly, Stuart and Higham30], the technical results proven in [Reference Löcherbach27] can be used to prove geometric ergodicity of $(X_t)_{t\geq 0}$ as stated in Proposition 1 below.

In order to state the geometric ergodicity of $(X_t)_{t\geq 0}$ , let us first specify the Lyapunov function $G\;:\;\mathbb{R}^\kappa\to \mathbb{R}$ introduced in [Reference Ditlevsen and Löcherbach15]:

(14) \begin{equation}G(x) = \sum_{k=1}^K \sum_{j=1}^{\eta_k+1}\frac{j}{\nu^{j-1}_k}J(x^{k,j}),\end{equation}

where J is some smooth approximation of the absolute value. In particular, $J(x) = |x|$ for all $|x|\geq 1$ , and $\max\{ |J^\prime(x)|,|J^{\prime\prime}(x)|\} \leq J_c$ for all x, for some finite constant $J_c$ .

Proposition 1. (Geometric ergodicity.) Grant Assumption (A). Then the solution of the SDE (7) has a unique invariant measure $\pi$ on $\mathbb{R}^\kappa$ . For all initial conditions $x_0$ and all $m\geq 1$ , there exist $C = C(m) >0$ and $\lambda = \lambda(m)>0$ such that, for all measurable functions $g\;:\;\mathbb{R}^\kappa\to \mathbb{R}$ such that $|g|\leq G^m$ ,

\begin{equation*}\quad \left| \mathbb{E} g(X_t) - \pi(g)\right|\leq C G(x_0)^m e^{-\lambda t} \qquad \forall t \geq 0 .\end{equation*}

Proof. The proof closely follows that of Theorem 3.2 in [Reference Mattingly, Stuart and Higham30] and is based on Lyapunov and minorization conditions (the latter are implied by the existence of a smooth transition density and the irreducibility of the space).

1. (i) First we use the fact that G is a Lyapunov function for X [Reference Ditlevsen and Löcherbach15, Proposition 5], i.e., that there exist $\alpha, \beta >0$ such that

   \begin{equation*}\mathcal{A}^XG(x) \leq -\alpha G(x) + \beta,\end{equation*}
   where $ \mathcal{A}^XG(x)$ is the infinitesimal generator of (7).

2. (ii) Then we note that, from any initial condition $x_0$ , for any time $T>0$ and any open set O, the probability that $X_T$ belongs to O is positive. This is ensured by the controllability of the system (7) (see Theorem 4 in [Reference Löcherbach27]).

3. (iii) Finally, we note that the process $(X_t)_{t\geq 0}$ possesses a smooth transition density. Its existence is ensured by verifying the Hörmander condition, which is done in Proposition 7 of [Reference Ditlevsen and Löcherbach15].

The rest of the proof follows that of [Reference Mattingly, Stuart and Higham30, Theorem 3.2]: apply [Reference Mattingly, Stuart and Higham30, Theorem 2.5] to some discrete-time sampling of the process and conclude by interpolation.

Also note that the rank of the diffusion matrix $\sigma\sigma^T$ is smaller than the dimension of the system (7). This means that the system is not elliptic. However, the specific cascade structure of the drift ensures that the noise is propagated through the whole system via the drift term, so that the diffusion is hypoelliptic in the sense of stochastic calculus of variations [Reference Delarue and Menozzi13, Reference Malliavin and Thalmaier28]. We also note that the SDE (7) is semi-linear, with a linear term given by the matrix (6). Thus, its solution can be written in the form of a convolution equation (see, among others, [Reference Mao29, Section 3]).

Proposition 2. The unique solution of (7) satisfies

(15) \begin{equation}X_t = e^{At}x_0 + \int_0^t e^{A(t-s)} B(X_s) ds + \frac{1}{\sqrt{N}} \int_0^t e^{A(t-s)} \sigma(X_s) dW_s.\end{equation}

Proof. Consider the process $Y_t = e^{-At}X_t$ . By Itô’s formula one obtains

\begin{eqnarray*}d\left(e^{-At}X_t\right) &\,=\,& \left( -Ae^{-At}X_t + e^{-At} \left( AX_t + B(X_t)\right)\right)dt + \frac{e^{-At} }{\sqrt{N}} \sigma(X_t) dW_t \\[5pt] &\,=\,&e^{-At} B(X_t)dt + \frac{e^{-At}}{\sqrt{N}} \sigma(X_t) dW_t.\end{eqnarray*}

Integrating both parts yields

\begin{equation*}e^{-At}X_t = x_0 + \int_0^t e^{-As} B(X_s) ds + \frac{1}{\sqrt{N}} \int_0^t e^{-As} \sigma(X_s) dW_s.\end{equation*}

Multiplying the expression by $e^{At}$ gives the result.

Note that from this form, it is straightforward to see that the diffusion term is of full rank. Intuitively, this ensures the hypoellipticity. Systems of type (15) are called stochastic Volterra equations [Reference Abi Jaber, Larsson and Pulido1].

Now we focus on first and second moment bounds. The following results are needed, in particular, to ensure the accuracy of the approximation scheme in Section 3. In the following remark we provide some purely computational results in order to ease the further analysis.

Remark 2. Because of the block structure of the matrix A introduced in (6), its matrix exponential $e^{At}$ can be computed as

\begin{equation*}e^{At} = \left(\begin{matrix}e^{A_{\nu_1}t} & \quad\mathbb{O}_{(\eta_1+1) \times (\eta_2+1)} \\[8pt] \mathbb{O}_{(\eta_2+1) \times (\eta_1+1)} & \quad e^{A_{\nu_2}t}\end{matrix}\right),\end{equation*}

where $e^{A_{\nu_k}t}$ , $k=1,2$ , is an $(\eta_k+1)\times (\eta_k+1)$ upper-triangular matrix given by

(16) $${e^{{A_{{\nu _k}}}t}} = {e^{ - {\nu _k}t}}\left( {\matrix{ 1 & {\quad t} & {\quad {{{t^2}} \over 2}} & {\quad \ldots } & {\quad {{{t^{{\eta _k}}}} \over {{\eta _k}!}}} \cr 0 & {\quad 1} & {\quad t} & {\quad \ldots } & {\quad {{{t^{{\eta _k} - 1}}} \over {({\eta _k} - 1)!}}} \cr \vdots & {\quad \ddots } & {\quad \ddots } & {\quad \ddots } & {\quad \vdots } \cr \vdots & {\quad \vdots } & {\quad \ddots } & {\quad \ddots } & {\quad \vdots } \cr 0 & {\quad 0} & {\quad 0} & {\quad {\rm{ }} \ldots } & {\quad 1} \cr } } \right).$$

In further computations we will often use the vectors $e^{At} X_{s}$ . The elements of $e^{At} X_{s}$ are given by the formula

(17) \begin{equation}\left(e^{At} X_{s}\right)^{k,j} = e^{-\nu_kt}\sum_{m=j }^{\eta_k+1}\frac{t^{m-j}}{(m-j)!}X^{k,m}_s.\end{equation}

Theorem 3. (First moment bounds of the diffusion process.) Grant Assumption (A). The following bounds hold for the components of $\mathbb{E}[X_t]$ :

\begin{equation*}\mathcal{I}^{k,j}_{\min} \leq \mathbb{E}[X^{k,j}_t] \leq \mathcal{I}^{k,j}_{\max},\end{equation*}

where

\begin{align*}\mathcal{I}^{k,j}_{\min} &= \left(e^{At} x_{0}\right)^{k,j} + \left[1-e^{-t\nu_k}\sum_{l=0}^{\eta_k+1-j}\frac{(t\nu_k)^l}{l!} \right]\min\left\{0, \frac{c_k f_{k+1}^{\max}}{\nu_k^{(\eta_k+2-j)} } \right\}, \\[5pt]\mathcal{I}^{k,j}_{\max} &= \left(e^{At} x_{0}\right)^{k,j} + \left[1-e^{-t\nu_k}\sum_{l=0}^{\eta_k+1-j}\frac{(t\nu_k)^l}{l!} \right]\max\left\{0, \frac{c_k f_{k+1}^{\max}}{\nu_k^{(\eta_k+2-j)} } \right\}.\end{align*}

Proof of Theorem 3. From Proposition 2 and Remark 2, it follows that the convolution-based representation of the kth population is given by

\begin{align*}X^{k}_t = (e^{At}x_0)^k+\int_0^t e^{A_{\nu_k}(t-s)} B^k(X^{k+1}_s) ds + \frac{1}{\sqrt{N}} \int_0^t e^{A_{\nu_k}(t-s)} \sigma^k(X^{k+1}_s) dW_s. \end{align*}

Consequently, the jth components are given by

\begin{eqnarray*}X^{k,j}_t &\,=\,& \underbrace{\left(e^{At} x_{0}\right)^{k,j}}_{\,{:}\,{\raise-1.5pt{=}}\,T_1(t)}+\underbrace{\int_0^t c_k f_{k+1}(X^{{k+1},1}_s) \frac{e^{-{\nu_k}(t-s)}}{(\eta_k+1-j)!} (t-s)^{\eta_k+1-j} ds}_{\,{:}\,{\raise-1.5pt{=}}\,T_2(t)} \\[5pt] &+& \underbrace{\frac{1}{\sqrt{N}} \int_0^t \frac{c_k}{\sqrt{p_{k+1}}}\sqrt{f_{k+1}(X^{{k+1},1}_s)} \frac{e^{-{\nu_{k}}(t-s)}}{(\eta_{k}+1-j)!} (t-s)^{\eta_{k}+1-j} dW^{k+1}_s}_{\,{:}\,{\raise-1.5pt{=}}\,T_3(t)}.\end{eqnarray*}

Note that $\mathbb{E}[T_1(t)]=T_1(t)$ and $\mathbb{E}[T_3(t)]=0$ . It remains to consider $T_2(t)$ . The fact that the intensity function is bounded by $0<f_{k+1}\leq f_{k+1}^{\max}$ implies that

\begin{align*}\min\{0,c_k\} \frac{f_{k+1}^{\max}}{(\eta_k+1-j)!} I^{k,j} \leq \mathbb{E}[T^2(t)] \leq \max\{0,c_k\} \frac{f_{k+1}^{\max}}{(\eta_k+1-j)!} I^{k,j} ,\end{align*}

where

\begin{equation*}I^{k,j} = \int_0^t e^{-{\nu_k}(t-s)} (t-s)^{\eta_k+1-j} ds.\end{equation*}

Now, let us consider the integral $I^{k,j}$ :

\begin{equation*}\int_0^t e^{-{\nu_k}(t-s)} (t-s)^{\eta_k+1-j} ds = t^{\eta_k+1 - j } \int_0^t e^{-{\nu_kt}\frac{t-s}{t}} \left(\frac{t-s}{t}\right)^{\eta_k+1-j} ds.\end{equation*}

Setting $z = \frac{t-s}{t}$ yields

\begin{equation*}t^{\eta_k+2 - j } \int_0^1e^{-{\nu_kt}z} z^{\eta_k+1-j} dz = \frac{(\eta_k+1-j)!}{\nu_k^{(\eta_k+2-j)}}\left[1-e^{-\nu_kt}\sum_{l=0}^{\eta_k+1-j}\frac{(\nu_kt)^l}{l!} \right].\end{equation*}

This gives the result.

Remark 3. Recalling (17) and using that

\begin{equation*}\lim\limits_{t\to\infty} e^{-\nu_kt}\sum_{l=0}^{\eta_k+1-j}\frac{(t\nu_k)^l}{l!}=0,\end{equation*}

it follows from Theorem 3 that

\begin{equation*}\min\left\{0, \frac{c_k f_{k+1}^{\max}}{\nu_k^{(\eta_k+2-j)} } \right\} \leq \lim_{t\to\infty} \mathbb{E}[X^{k,j}_t] \leq \max\left\{0, \frac{c_k f_{k+1}^{\max}}{\nu_k^{(\eta_k+2-j)} } \right\}.\end{equation*}

The derived moment bounds give some intuition on how the system behaves in the long run. Remarkably, depending on whether $c_k$ is positive or negative, the trajectories of $(X_t^k)_{t \geq 0}$ are on average bounded by 0 from below or above, respectively. This is in agreement with the fact that the sign of $c_k$ defines whether the corresponding neural population is excitatory ( $c_k = +1$ ) or inhibitory ( $c_k = -1$ ). Moreover, we may immediately see the effect of increasing the memory order $\eta_k$ , depending on the constant $\nu_k$ . When $\nu_k = 1$ , the bounds for all j components are determined entirely by $c_k$ and the bounds of the intensity functions. When $\nu_k < 1$ and $\eta_k \to \infty$ , the first components, presenting the current state of the process, tend to infinity. Similarly, for $\nu_k > 1$ , the trajectories are attracted to 0. Finally, note that the first moment bounds do not depend on the number of neurons in the system.

Theorem 4. (Second moment bounds of the diffusion process.) Grant Assumption (A). The following bounds hold for $\mathbb{E}[(X^{k,j}_t)^2]$ :

\begin{eqnarray*}\mathbb{E}[(X^{k,j}_t)^2] &\leq& \left( \left(e^{At} x_{0}\right)^{k,j} \right)^2 + 2 \left(e^{At} x_{0}\right)^{k,j}\max\left\{0,\frac{c_k f^{\max}_{k+1}}{(\eta_k+1-j)!} I^{k,j}_1(t)\right\} \\[5pt] &+&\;f^{\max}_{k+1}\left(\frac{c_k }{(\eta_k+1-j)!}\right)^2\left(\sqrt{f^{\max}_{k+1}} \: I^{k,j}_1(t) + \sqrt{\frac{I^{k,j}_2(t)}{N\cdot p_{k+1}}} \right)^2,\end{eqnarray*}

where $I^{k,j}_u(t), \: u = 1,2$ , are defined as

\begin{eqnarray*}I^{k,j}_u(t) &\,{:}\,{\raise-1.5pt{=}}\,& \int_0^t e^{-u\nu_k(t-s)}(t-s)^{u(\eta_k+1-j)}ds \\[5pt] &\,=\,& \frac{(u(\eta_k + 1 - j))!}{(u\nu_k)^{u(\eta_k+1-j)+1}} \left[ 1-e^{-ut\nu_k} \sum_{l=0}^{u(\eta_k+1-j)}\frac{(ut\nu_k)^l}{l!} \right].\end{eqnarray*}

The proof of Theorem 4 is similar to that of Theorem 3 and is postponed to Appendix A.2.

Remark 4. Theorem 4 gives the following asymptotic bounds:

\begin{align*}\lim_{t\to\infty} \mathbb{E}[(X^{k,j}_t)^2] \leq \;f^{\max}_{k+1}\left(\frac{c_k }{(\eta_k+1-j)!}\right)^2 \left(\sqrt{f^{\max}_{k+1}} C^{k,j}_1 + \sqrt{\frac{C^{k,j}_2}{N\cdot p_{k+1}}}\right)^2,\end{align*}

where

\begin{equation*}C^{k,j}_u\,{:}\,{\raise-1.5pt{=}}\,\lim_{t\to\infty}I^{k,j}_u(t)=\frac{(u(\eta_k + 1 - j))!}{(u\nu_k)^{u(\eta_k+1-j)+1}}.\end{equation*}

Note that for $N \to \infty$ , the bound obtained in Theorem 4 equals the square of the bound for the first moment, derived in Theorem 3. This is in agreement with the fact that the stochastic system (7) transforms into an ODE as N increases [Reference Ditlevsen and Löcherbach15]. In other words, its diffusion coefficient tends to $\mathbb{O}_{\kappa\times 2}$ as N tends to infinity.

In Figure 1, the first and second moment bounds, derived in Theorem 3 and Theorem 4, respectively, are illustrated. In the left panel, we plot four sample trajectories (solid lines) of an inhibitory population $X^1$ and their lower first moment bounds (dashed lines). The main variable $X^{1,1}$ and its lower moment bound are depicted in black. The remaining three trajectories $X^{1,j}$ , $j=2,3,4$ , are auxiliary variables. They (and their corresponding bounds) are depicted in different shades of grey. We see that the trajectories can exceed the theoretical bounds, especially when the effect of noise is large. On average, the trajectories stay within the bounds. In the right panel, we plot the square of the first three components of $X^1$ (and their second moment bounds), omitting the fourth in order to stay within an easily interpretable scale. We conclude that the bounds are fairly precise for the parameter setting under consideration.

Figure 1. First (left panel) and second (right panel) moment bounds with respective trajectories of the inhibitory population $k=1$ . The rate function $f_2$ is given in Section 5. The parameters are $\eta_1=3$ , $\nu_1=2$ , $N=20$ , and $p_2=1/2$ .

3.1. Strong convergence in the mean-square sense

We now focus on the convergence in the mean-square sense and show that the numerical solutions obtained via the splitting approach converge to the process as the time step $\Delta\to 0$ with order 1. The frequently applied Euler–Maruyama scheme usually converges with mean-square order $1/{2}$ if the noise is multiplicative [Reference Kloeden, Platen and Schurz22, Reference Milstein and Tretyakov33], as is the case for the system (7). In the following result, thanks to the specific structure of the noise component, we show that the Euler–Maruyama scheme coincides with the Milstein scheme, which is known to converge with mean-square order 1. This result is then used to establish the convergence order of the splitting scheme.

Theorem 5. (Mean-square convergence of the splitting scheme.) Grant Assumption (A). Let $\tilde{X}_{t_i}$ denote the numerical method defined by (20) at time point $t_i$ and starting from $x_0$ . Then $\tilde{X}_{t_i}$ is mean-square convergent with order 1; i.e., there exists a constant $C>0$ such that

\begin{equation*}\left(\mathbb{E}\left[ \left\|X_{t_{i}} - \tilde X_{t_{i}}\right\|^2 \right]\right)^{\frac{1}{2}} \leq C\Delta\end{equation*}

for all time points $t_i$ , $i=1,\ldots,i_{\text{max}}$ , where $\|\cdot \|$ denotes the Euclidean norm.

Proof of Theorem 5. Let us denote by $\tilde{X}^{EM}$ a numerical solution of the SDE (7) obtained via the Euler–Maruyama method; that is,

(22) \begin{equation}\tilde X^{EM}_{t_{i+1}} = \tilde X^{EM}_{t_i} + \Delta\left( A\tilde X^{EM}_{t_i}+B(\tilde X^{EM}_{t_i}) \right) + \frac{\sqrt{\Delta}}{\sqrt N} \sigma(\tilde X^{EM}_{t_i})\xi_i.\end{equation}

First, we show that the Euler–Maruyama method, when applied to the system (7), coincides with the Milstein scheme, which is known to converge with mean-square order 1. To do so, we define the vector x by $x=(x^1,\ldots,x^{\kappa})$ , where $\kappa=\eta_1+\eta_2+2$ . Further, we recall that the jth component, $j=1,\ldots,\kappa$ , of the Milstein scheme differs from the jth component of the Euler–Maruyama scheme (22) only by the following additional term:

\begin{equation*}\sum_{m_1, m_2 = 1}^{2}\sum_{l=1}^{\kappa} \sigma^{l,m_1} \frac{\partial \sigma^{j,m_2}}{\partial x^l}I(m_1,m_2),\end{equation*}

where $\sigma^{j,m}$ denotes the value of the element at the jth row and the mth column of the diffusion matrix $\sigma$ at time $t_i$ , and

\begin{equation*}I(m_1,m_2)\,{:}\,{\raise-1.5pt{=}}\,\int_{t_i}^{t_{i+1}} \int_{t_i}^{{s_2}} dW_{s_1}^{m_1} dW_{s_2}^{m_2}.\end{equation*}

Now note that the term ${\partial \sigma^{j,m_2}}/{\partial x^l}$ is different from 0 only for $j=\eta_1+1$ , $m_2=1$ , $l=\eta_1+2$ and for $j=\eta_1+\eta_2+2$ , $m_1=2$ , $l=1$ . However, $\sigma^{l,m_1}$ equals 0 for those values of l. Thus, the above double sum equals 0 and the Euler–Maruyama method coincides with the Milstein scheme. This implies that

(23) \begin{equation}\| X_{t_i}-\tilde{X}^{EM}_{t_i} \|_{L^2} \leq C \Delta,\end{equation}

where $\|\cdot \|_{L^2}\,{:}\,{\raise-1.5pt{=}}\,\left( \mathbb{E}[ \|\cdot \|^2 ] \right)^{1/2}$ denotes the $L^2$ norm and C is some generic constant. For the second part, we provide a proof similar to the one presented in [Reference Milstein and Tretyakov32]. Applying the triangle inequality yields that

\begin{equation*}\| X_{t_i}-\tilde{X}_{t_i} \|_{L^2} \leq \| X_{t_i}-\tilde{X}^{EM}_{t_i} \|_{L^2} + \| \tilde{X}^{EM}_{t_i}-\tilde{X}_{t_i} \|_{L^2}.\end{equation*}

Given $X_{t_i}\,{:}\,{\raise-1.5pt{=}}\,x$ , let us denote by $\tilde{X}^{EM}_{t_{i+1}}(x,t_i)$ and $\tilde{X}_{t_{i+1}}(x,t_i)$ the one-step approximation of the Euler–Maruyama and splitting scheme, respectively. For instance, $\tilde{X}^{EM}_{t_{i+1}}(x,t_i)$ is given by Equation (22) with $\tilde X^{EM}_{t_i}$ replaced by x. By the definition of the matrix exponent, i.e. $e^{A\Delta} \,{:}\,{\raise-1.5pt{=}}\, I + \Delta A + \frac{\Delta^2}{2} A^2 + O(\Delta^3)$ , and by recalling (20), we obtain that

\begin{align*}\tilde X^{EM}_{t_{i+1}}(x,t_i)- \tilde X_{t_{i+1}}(x,t_i) &= x + \Delta A x + \Delta B( x) + \frac{\sqrt{\Delta}}{\sqrt N} \sigma(x)\xi_i \\[1pt] &- e^{A\Delta}\left( x + \Delta B( x) + \frac{\sqrt{\Delta}}{\sqrt N} \sigma(x)\xi_i \right) \\[1pt] &= x + \Delta Ax + \Delta B( x) + \frac{\sqrt{\Delta}}{\sqrt N} \sigma(x)\xi_i \\[1pt] &- x - \Delta B( x) - \frac{\sqrt{\Delta}}{\sqrt N} \sigma(x)\xi_i \\[1pt] &-\Delta Ax - \Delta^2 A B( x) - \frac{\Delta^{\frac{3}{2}}}{\sqrt N} A\;\sigma(x)\xi_i + R \\[1pt] &= - \Delta^2 A B( x) - \frac{\Delta^{\frac{3}{2}}}{\sqrt N} A\;\sigma(x)\xi_i + R,\end{align*}

where R is a corresponding stochastic remainder. Consequently, we get that

\begin{align*}\left\| \mathbb{E}\left[\tilde X^{EM}_{t_{i+1}}(x,t_i)- \tilde X_{t_{i+1}}(x,t_i) \right] \right\|=O(\Delta^2) , \\[5pt] \quad \left\|\tilde X^{EM}_{t_{i+1}}(x,t_i)- \tilde X_{t_{i+1}}(x,t_i) \right\|_{L^2} = O( \Delta^{\frac{3}{2}}).\end{align*}

Recalling (23), the result follows from the fundamental theorem on the mean-square order of convergence; see Theorem 1.1. in [Reference Milstein and Tretyakov33].

Theorem 5 states that as $\Delta\to 0$ , the approximated solution $(\tilde X_{t_i})_{i=0,\ldots,i_{\text{max}}}$ converges to the true process $(X_t)_{t\geq 0}$ in the mean-square sense. In practice, however, fixed time steps $\Delta>0$ are required. Thus, there is not yet any guarantee that the constructed numerical solutions share the same properties as the true solution of (7). For these reasons, we additionally study the ability of $(\tilde X_{t_i})_{i=0,\ldots,i_{\text{max}}}$ to preserve the derived first and second moment bounds and the geometric ergodicity of the SDE (7).

Note also that, unlike the case of ODE systems [Reference Hairer, Lubich and Wanner19], for stochastic equations the theoretical order of convergence usually cannot be increased by using the Strang composition instead of the Lie–Trotter approach. In practice, however, the Strang splitting often performs better than the Lie–Trotter method; see e.g. [Reference Ableidinger, Buckwar and Hinterleitner3, Reference Buckwar, Tamborrino and Tubikanec7]. This is also confirmed by our numerical experiments in Section 5.

3.2. Moment bounds of the approximated process

We are now interested in studying the qualitative properties of the splitting schemes for fixed time steps $\Delta>0$ . We start by illustrating that the constructed splitting schemes preserve the convolution-based structure of the model derived in Proposition 2. Using the one-step approximation (20) and performing back iteration yields

(24) \begin{equation}\tilde X_{t_i} = e^{At_i} x_0 + \Delta \sum_{l=1}^{i} e^{At_l}B(\tilde X_{t_{i-l}}) + \frac{\sqrt{\Delta}}{\sqrt N}\sum_{l=1}^{i} e^{At_l} \sigma(\tilde X_{t_{i-l}})\xi_{i-l}.\end{equation}

Note that the first term on the right side of (15) is preserved exactly. Moreover, the sums in (24) correspond to approximations of the integrals in (15) using the left-point rectangle rule. The expression (24) allows one to derive moment bounds for the numerical process in a similar fashion to that presented for the continuous process in the previous section.

Theorem 6. (First moment bounds of the approximated process.) Grant Assumption (A). The following bounds hold for the components of $\mathbb{E}[\tilde X_{t_i}]$ :

\begin{equation*}\tilde{\mathcal{I}}^{k,j}_{\min} \leq \mathbb{E}[\tilde X^{k,j}_{t_i}] \leq \tilde{\mathcal{I}}^{k,j}_{\max},\end{equation*}

where

\begin{align*}\tilde{\mathcal{I}}^{k,j}_{\min} &= \left(e^{A{t_i}} x_0 \right)^{k,j} + \Delta \sum_{l=0}^{i}e^{-{\nu_k}t_l} t_l^{\eta_k+1-j}\min\left\{0, \frac{c_k f_{k+1}^{\max}}{(\eta_k+1-j)! } \right\}, \\[5pt]\tilde{\mathcal{I}}^{k,j}_{\max} &= \left(e^{At_i} x_0\right)^{k,j} + \Delta \sum_{l=0}^{i}e^{-{\nu_k}t_l} t_l^{\eta_k+1-j}\max\left\{0, \frac{c_k f_{k+1}^{\max}}{(\eta_k+1-j)! } \right\}.\end{align*}

Proof of Theorem 6. From Remark 2 and (24), it follows that

\begin{align*}\tilde{X}^{k}_{t_i} = (e^{At_i}x_0)^k+\Delta \sum_{l=1}^{i} e^{A_{\nu_k}t_l} B^k(\tilde{X}^{k+1}_{t_{i-l}}) + \frac{\sqrt{\Delta}}{\sqrt{N}} \sum_{l=1}^{i} e^{A_{\nu_k}t_l} \sigma^k(\tilde{X}^{k+1}_{t_{i-l}}) \xi_{i-l}.\end{align*}

Consequently, the jth components are given by

\begin{align*}\tilde{X}^{k,j}_{t_i} &= \underbrace{(e^{At_i}x_0)^{k,j}}_{=T_1(t_i)}+\underbrace{ \frac{1}{(\eta_k+1-j)!} \Delta \sum_{l=1}^{i} c_k f_{k+1}(\tilde{X}^{{k+1},1}_{t_{i-l}}) e^{-{\nu_k}t_l} t_l^{\eta_k+1-j}}_{\,{:}\,{\raise-1.5pt{=}}\,\tilde{T}_2(t)} \\[5pt] &+ \underbrace{\frac{1}{(\eta_k+1-j)!} \frac{\sqrt{\Delta}}{\sqrt{N}} \sum_{l=1}^{i} \frac{c_k}{\sqrt{p_{k+1}}}\sqrt{f_{k+1}(\tilde{X}^{{k+1},1}_{t_{i-l}})} e^{-{\nu_k}t_l} t_l^{\eta_k+1-j} \xi^{k+1}_{i-l}}_{\,{:}\,{\raise-1.5pt{=}}\,\tilde{T}_3(t_i)}.\end{align*}

Note that $\mathbb{E}[T_1(t_i)]=T_1(t_i)$ and $\mathbb{E}[\tilde{T}_3(t)]=0$ . The fact that the intensity function is bounded by $0<f_{k+1}\leq f_{k+1}^{\max}$ implies the result.

Note that the bounds obtained in Theorem 6 equal those derived in Theorem 3, up to replacing the integrals (calculated in the proof of Theorem 3) by left Riemann sums. The accuracy of this approximation depends on the step size $\Delta$ . Under reasonably small choices of $\Delta$ , the bounds are preserved accurately for all $t_i$ . This is illustrated in the left panel of Figure 2, where we plot the first moment bound of the process (main variable of an excitatory population) and that of the approximated process, derived in Theorem 3 and Theorem 6, respectively. Different choices of $\nu_k$ are compared, and $\Delta=0.1$ is used for the bound of the approximated process.

Figure 2. First (left panel) and second (right panel) moment bounds of the excitatory population $k=2$ for different values of $\nu_2$ . The moment bounds for the diffusion are in solid lines and the moment bounds for the splitting scheme are in dashed lines. The bound of the rate function is fixed to $f_1^{\max} = 1$ . The parameters are $\eta_2=3$ , $N=100$ , $p_{1}=1/2$ , and the time step $\Delta=0.1$ is used.

The following corollary gives intuition for the long-time behaviour of the bounds.

Proof. (i) The zero bound is trivial. Considering

\begin{align*}\lim_{i\to\infty}\sum_{l=0}^{i}e^{-{\nu_k}t_l} t_l^{\kappa^{k,j}} &= \lim_{i\to\infty}\sum_{l=0}^{i}e^{-{\nu_k}l\Delta } (l\Delta)^{\kappa^{k,j}} \\[5pt] &= \lim_{i\to\infty}\Delta^{\kappa^{k,j}}\sum_{l=0}^{i}e^{-{\nu_k}l\Delta} l^{\kappa^{k,j}} = \Delta^{\kappa^{k,j}} Li_{-\kappa^{k,j}}\left(e^{-\nu_k\Delta} \right)\end{align*}

gives the result. The explicit form of the function is given in [Reference Wood43].

(ii) Let us rewrite once again the expression included in the limit:

\begin{align*}\left({-}1 \right)^{\kappa^{k,j}+1}& \lim_{\Delta\to 0} \frac{(\nu_k\Delta)^{\kappa^{k,j}+1}}{\nu^{\kappa^{k,j}+1}_k\kappa^{k,j}!} \sum_{m=0}^{\kappa^{k,j}} m!\: S(\kappa^{k,j}+1 ,m+1)\left(\frac{-1}{1-e^{-\nu_k\Delta}} \right)^{m+1} \\[5pt] &= \lim_{\Delta\to 0} \left[\frac{S(\kappa^{k,j}+1 ,\kappa^{k,j}+1)\kappa^{k,j}!}{\nu^{\kappa^{k,j}+1}_k\kappa^{k,j}!} \left(\frac{\nu_k\Delta}{1-e^{-\nu_k\Delta}} \right)^{\kappa^{k,j}+1} \right. \\[5pt] &- \left. \Delta \frac{S(\kappa^{k,j}+1 ,\kappa^{k,j})(\kappa^{k,j}-1)!}{\nu^{\kappa^{k,j}}_k\kappa^{k,j}!} \left(\frac{\nu_k\Delta}{1-e^{-\nu_k\Delta}} \right)^{\kappa^{k,j}} + O(\Delta^2) \right].\end{align*}

Note that

\begin{equation*}\lim_{\Delta\to 0} \left(\frac{\nu_k\Delta}{1-e^{-\nu_k\Delta}} \right) = 1.\end{equation*}

This implies that in the limit ${1}/{\nu^{\kappa^{k,j}+1}_k}$ is the only remaining term, since the rest converges to 0 as $\Delta\to 0$ . This gives the result.

In the first part of Corollary 1, the sums in Theorem 6 are calculated explicitly as $i \to \infty$ . This limit is described by polylogarithm functions. Note that the zero bounds derived in Remark 3, i.e. the upper bounds for the inhibitory population ( $c_k = -1$ ) and the lower bounds for the excitatory population ( $c_k = +1$ ), are preserved exactly by the splitting scheme for all times $t_i$ and for all choices of $\Delta>0$ . Moreover, the lower bounds for the inhibitory population and the upper bounds for the excitatory population are preserved accurately as $i \to \infty$ , provided that $\Delta$ is reasonably small. Indeed, as $i\to\infty$ and $\Delta\to 0$ (second part of Corollary 1), the bounds coincide with the ones obtained in Remark 3.

Theorem 7. (Second moment bounds of the approximated process.) Grant Assumption (A). Each component of $\mathbb{E}[(\tilde X^k_t)^2]$ is bounded by

\begin{align*}\mathbb{E}[(\tilde X^{k,j}_t)^2] &\leq \left( \left(e^{At} \ x_0\right)^{k,j} \right)^2 + 2 \left(e^{At} x_0\right)^{k,j}\max\left\{0,\frac{c_k f^{\max}_{k+1}}{(\eta_k+1-j)!} \tilde I^{k,j}_1(t)\right\} \\[5pt] &+f^{\max}_{k+1}\left(\frac{c_k }{(\eta_k+1-j)!}\right)^2\left(\sqrt{f^{\max}_{k+1}} \: \tilde I^{k,j}_1(t) + \sqrt{\frac{\tilde I^{k,j}_2(t)}{N\cdot p_{k+1}}}\right)^2,\end{align*}

where $\tilde I^{k,j}_u(t), \: u = 1,2$ , are defined as

\begin{equation*}\tilde I^{k,j}_u(t) \,{:}\,{\raise-1.5pt{=}}\, \Delta\sum_{l=0}^i e^{-u\nu_kt_l}t_l^{u(\eta_k+1-j)}.\end{equation*}

Proof of Theorem 7. The proof repeats the proof of Theorem 4, up to replacing integrals $I^{k,j}_u(t)$ by sums $\tilde I^{k,j}_u(t)$ .

As before, the second moment bounds obtained for the splitting scheme equal those derived for the true process in Theorem 4, except that the integrals are replaced by the corresponding Riemann sums. Using the same arguments as in the proof of Corollary 1, we conclude that the second moment bounds are also preserved accurately by the splitting scheme for reasonable choices of the time step $\Delta$ . A comparison of the theoretical and discrete second moment bounds is provided in the right panel of Figure 2.

3.3. Geometric ergodicity of the approximated process

Finally, our aim is to prove that the splitting scheme preserves the ergodic property of the underlying process, in the spirit of [Reference Ableidinger, Buckwar and Hinterleitner3, Reference Mattingly, Stuart and Higham30], providing a discrete analogue of Proposition 1. The main step is to establish a discrete Lyapunov condition for the approximated solution $(\tilde X_{t_i})_{i=0,\ldots,i_{\text{max}}}$ . This is granted by the following lemma.

Lemma 1. (Lyapunov condition for the approximated process.) Grant Assumption (A). The functional $\tilde G$ given by

\begin{equation*}\tilde G(x) = \sum_{k=1}^2\sum_{j=1}^{\eta_k+1}\frac{j}{\nu_k^{j-1}}\left|x^{k,j}\right|\end{equation*}

is a Lyapunov function for $\tilde X$ ; i.e., there exist constants $\alpha \in [0,1)$ and $\beta \geq 0$ such that

\begin{equation*}\mathbb{E}\left[\tilde G(\tilde X_{t_{i+1}})\vert \tilde X_{t_i} \right] \leq \alpha \tilde G(\tilde X_{t_i})+\beta.\end{equation*}

Proof. We bound the approximated solution obtained via (20) from above by a sum of three terms, thanks to the triangle inequality:

\begin{align*}\tilde G( \tilde X_{t_{i+1}}) &= \tilde G\left(e^{A\Delta} \tilde X_{t_i} + \Delta e^{A\Delta} B( \tilde X_{t_i}) + \frac{\sqrt{\Delta}}{\sqrt N} e^{A\Delta} \sigma( \tilde X_{t_i})\xi_i\right) \\[5pt] &\leq\underbrace{\tilde G\left(e^{A\Delta} \tilde X_{t_i}\right)}_{T_1} + \underbrace{\Delta \tilde G\left(e^{A\Delta} B( \tilde X_{t_i})\right)}_{T_2} + \underbrace{\frac{\sqrt{\Delta}}{\sqrt N} \tilde G\left(e^{A\Delta} \sigma( \tilde X_{t_i})\xi_i\right)}_{T_3}.\end{align*}

Note that the term $T_2$ , as well as the expectation of $T_3$ , is bounded by a constant depending on $f^{\max}_k$ , so that $\mathbb{E}[T_2\vert \tilde X_{t_i}]+\mathbb{E}[T_3\vert \tilde X_{t_i}]\leq \beta$ , and $\beta > 0$ since we consider the absolute value. Further, using the formulas (16)–(17), we can expand $T_1$ as follows:

\begin{eqnarray*}&&\tilde G\left(e^{A\Delta} \tilde X_{t_i}\right) = \sum_{k=1}^2e^{-\nu_k\Delta}\sum_{j=1}^{\eta_k+1}\frac{j}{\nu_k^{j-1}}\left|\sum_{m=j}^{\eta_k+1} \frac{\Delta^{m-j}}{(m-j)!} \tilde X^{k,m}_{t_i}\right|\\[5pt] &&\leq \sum_{k=1}^2e^{-\nu_k\Delta}\sum_{j=1}^{\eta_k+1}\frac{j}{\nu_k^{j-1}}\sum_{m=j}^{\eta_k+1} \frac{\Delta^{m-j}}{(m-j)!}\left| \tilde X^{k,m}_{t_i}\right| \end{eqnarray*}

\begin{eqnarray*} &\,=\,& \sum_{k=1}^2e^{-\nu_k\Delta}\left[\sum_{j=1}^{\eta_k+1}\frac{j}{\nu_k^{j-1}} \left| \tilde X^{k,j}_{t_i}\right|+\Delta \sum_{j=2}^{\eta_k+1}\frac{j-1}{\nu_k^{j-2}} \left| \tilde X^{k,j}_{t_i}\right|+\dots + \frac{\Delta^{\eta_k}}{\eta_k!} \left| \tilde X^{k,\eta_k+1}_{t_i}\right|\right] \\[5pt] &\,=\,& \sum_{k=1}^2e^{-\nu_k\Delta}\left[\sum_{j=1}^{\eta_k+1}\frac{j}{\nu_k^{j-1}} \left| \tilde X^{k,j}_{t_i}\right|+\Delta\nu_k \sum_{j=2}^{\eta_k+1}\frac{j-1}{\nu_k^{j-1}} \left| \tilde X^{k,j}_{t_i}\right|+\dots + \frac{(\nu_k\Delta)^{\eta_k}}{\eta_k!} \frac{1}{\nu_k^{\eta_k}}\left| \tilde X^{k,\eta_k+1}_{t_i}\right|\right] .\end{eqnarray*}

Note that, since $\nu_k>0$ , for all $m\geq 1$ it holds that

\begin{equation*}\sum_{j=m}^{\eta_k+1}\frac{(j-m+1)}{\nu_k^j}\left| \tilde X^{k,j}_{t_i}\right| \leq \sum_{j=1}^{\eta_k+1}\frac{j}{\nu_k^j}\left| \tilde X^{k,j}_{t_i}\right| = \tilde G\left( \tilde X^k_{t_i}\right).\end{equation*}

Thus, we have

\begin{equation*}\tilde G\left(e^{A\Delta} \tilde X_{t_i}\right) \leq \sum_{k=1}^2\left(e^{-\nu_k\Delta}\sum_{j=0}^{\eta_k}\frac{(\nu_k\Delta)^{j}}{j!} \right) \tilde G\left( \tilde X^k_{t_i}\right) .\end{equation*}

Let

\begin{equation*}\alpha = \max_k \left(e^{-\nu_k\Delta}\sum_{j=0}^{\eta_k}\frac{(\nu_k\Delta)^{j}}{j!} \right).\end{equation*}

Since $\eta_k$ is finite, we get $\alpha < 1$ , which implies the result.

Note that the statement of Lemma 1 holds without any assumption on the time step $\Delta$ . Also, the Lyapunov function is the same as for the continuous process up to smoothing the absolute value (see (14)). With a discrete Lyapunov condition established, the ergodicity is conditioned on two further technical steps. First, the transition probability of two (or more) consecutive steps, given by the recursive relation (20), must have a smooth transition density. This fact is granted by the hypoellipticity of the scheme (see Remark 5).

Second, the irreducibility condition must hold. This means that any point $y\in \mathbb{R}^\kappa$ can be reached from any starting point $x\in \mathbb{R}^\kappa$ in a fixed number of steps. In other words, we need a discrete-time analogue of Theorem 4 in [Reference Löcherbach27], granting the controllability of the SDE (7). This is provided by the following lemma, which is proved in Appendix A.3.

Lemma 2. (Irreducibility condition.) Grant Assumption (A). Let $\eta^* = \max_k\{\eta_k\}$ . Then, for all $ x, y \in \mathbb{R}^\kappa$ , there exists some sequence of two-dimensional vectors $(\xi_i)_{i=1,\dots,\eta^*+1}$ such that

\begin{equation*}y = \underbrace{\left(\psi_\Delta[\xi_{\eta^*+1}] \circ \dots \circ \psi_\Delta[\xi_{1}] \right)}_{\eta^*+1} (x),\end{equation*}

where $\psi_\Delta$ denotes one step of the scheme defined by (20), and the notation $[\cdot]$ is introduced to stress the dependency on the vectors $(\xi_i)_{i=1,\dots,\eta^*+1}$ .

Lemmas 1 and 2, combined with the hypoellipticity of the scheme, gives the following result, which is analogous to Theorem 7.3 in [Reference Mattingly, Stuart and Higham30].

Theorem 8. (Geometric ergodicity.) Grant Assumption (A). Then the process $(\tilde X_{t_i})_{i=0,\dots,i_{\max}}$ has a unique invariant measure $\pi^\Delta$ on $\mathbb{R}^\kappa$ . For all initial conditions $x_0$ and all $m\geq 1$ , there exist $\tilde C = C(m, \Delta) >0$ and $\tilde\lambda = \tilde\lambda(m,\Delta)>0$ such that, for all measurable functions $g\;:\;\mathbb{R}^\kappa\to \mathbb{R}$ such that $|g|\leq \tilde G^m$ ,

\begin{equation*}\forall i=0,\dots,i_{\max}, \qquad \left| \mathbb{E} g(\tilde X_{t_i}) - \pi^\Delta(g)\right|\leq \tilde C \tilde G(x_0)^m e^{-\tilde \lambda t_i}.\end{equation*}

4.1. Choice of an upper bound for the intensity

If $\bar{Z}_t=0$ , i.e. in the absence of any spike, it follows from (5) that $\bar{X}$ evolves as a linear ODE with matrix A, so that $\bar{X}_t = e^{At}x_0$ . In particular, for all neurons $n=1,\dots, N_k$ , it follows that

(25) \begin{equation}\lambda^{k,n}_t = f_k((e^{At}x_0)^{k,1}).\end{equation}

One possible choice for the dominating intensity $\tilde{\lambda}$ in the thinning procedure is to provide an upper bound of (25) which holds for all $t\geq 0$ . A straightforward candidate for such a bound is provided in the following lemma.

Lemma 3. For any $x\in \mathbb{R}^\kappa$ , let $\Phi_k(x) = \sup_{t\geq 0} (e^{At}x)^{k,1}$ . Then

\begin{equation*}\Phi_k(x) \leq \tilde{\Phi}_k(x) = \max_{j=1,\dots, \eta_k +1} \left\{ 0, \frac{x^{k,j}}{\nu_k^{j-1}} \right\}.\end{equation*}

Proof. The explicit expression for $(e^{At}x)^{k,1}$ is given in (17); that is,

\begin{equation*}(e^{At}x)^{k,1} = e^{-\nu_k t} \left( x^{k,1} + tx^{k,2} +\dots + \frac{t^{\eta_k}}{\eta_k !} x^{k,\eta_k+1} \right).\end{equation*}

Setting $y_j = x^{k,j}/(\nu_k)^{j-1}$ , one gets

\begin{equation*}(e^{At}x)^{k,1} = e^{-\nu_k t} \left( y_1 + t\nu_k y_2 +\dots + \frac{(t\nu_k)^{\eta_k}}{\eta_k !} y_{\eta_k+1} \right) \leq \max_k\{0,y_k\} e^{-\nu_k t} g(t).\end{equation*}

The result follows from the fact that $g(t) = 1+t\nu_k + \dots + (t\nu_k)^{\eta_k}/\eta_k! \leq e^{\nu_k t}$ .

Remark 6. Another possible choice of a uniform bound, similar to the one given in Lemma 3, is provided in [Reference Duarte, Löcherbach and Ost16]. Their method, adapted to our case, gives

\begin{equation*}\Phi_k(x) \leq e \max\left\{ 1, \left(\frac{\eta_k}{e\nu_k}\right)^{\eta_k} \right\} \max_{j}\{x^{k,j}\},\end{equation*}

which is larger, and thus less efficient, than the bound proposed in Lemma 3.

Since the functions $f_k$ are non-decreasing, the upper bound of $(e^{At}x)^{k,1}$ given in Lemma 3 provides the bound

(26) \begin{equation}\tilde{f}_k(x) = f_k( \tilde{\Phi}_k(x))\end{equation}

on the intensity. However, there is no guarantee that this bound is sharp. In most practical cases (especially when the functions $f_k$ are increasing fast), the procedure rejects a vast majority of the simulated points. Hence, a more efficient approach, based on the computation of the critical points of the function $(e^{At}x)^{k,1}$ , is proposed. Furthermore, instead of considering a bound for any $t>0$ , we choose a fixed time step $\tilde\Delta>0$ (such that one spike is likely to occur in the interval $[0,\tilde\Delta]$ ) and compute $\Phi_k^{\tilde\Delta}(x) = \sup_{0\leq t\leq {\tilde\Delta}} (e^{At}x)^{k,1}$ instead of $\Phi_k(x)$ . This choice has no impact on the precision of the simulation. It only influences the sharpness of the bound used in the method and thus its computational efficiency.

Lemma 4. For any $x\in \mathbb{R}^\kappa$ , it holds that

\begin{equation*}\Phi_k^{\tilde\Delta}(x) = \max_{0<t_c <{\tilde\Delta}} \left\{x^{k,1}, \big(e^{At_c}x\big)^{k,1}, \big(e^{A{\tilde\Delta}}x\big)^{k,1}\right\},\end{equation*}

where the maximum is taken over the critical points $t_c$ of $t\mapsto (e^{At}x)^{k,1}$ , which are the solutions of the equation

\begin{equation*}\big(\!-\!\nu_k x^{k,1} + x^{k,2}\big) + \dots + \big(\!-\!\nu_k x^{k,\eta_k} + x^{k,\eta_k+1}\big) \frac{(t_c)^{\eta_k-1}}{(\eta_k-1)!} + \big(\!-\!\nu_k x^{k,\eta_k+1}\big) \frac{(t_c)^{\eta_k}}{(\eta_k)!} = 0.\end{equation*}

Proof. The result follows from the computation of the time derivative of $ (e^{At}x)^{k,1}$ .

The critical points in Lemma 4 are given by polynomial roots, which can be accurately computed numerically. In most practical cases, the computational cost of the polynomial roots is compensated by the efficiency gained in the rejection method. Finally, let us define the upper-bound intensity function by

(27) \begin{equation}\tilde{f}_k^{\tilde\Delta}(x) = f_k(\Phi_k^{\tilde\Delta}(x)).\end{equation}

Note that when the population is inhibitory ( $c_k = -1$ ), the naive upper bound $\tilde{f}_k$ is constant with respect to time, because all the coordinates of $\bar{X}^1$ are always negative and the bound given by Lemma 3 is 0. Thus, $\tilde{f}_k\equiv f_k(0)$ . Of course, such a bound is not sharp in general. However, it is interesting to see how the two upper bounds $\tilde{f}_k$ and $\tilde{f}^{\tilde\Delta}_k$ behave for a particular realization of the intensity process for excitatory populations. Figure 3 gives a comparison of the paths of $\tilde{f}_2$ and $\tilde{f}^{\tilde\Delta}_2$ for the excitatory population (with $\tilde\Delta\equiv 1$ ). We observe that both bounds are fairly precise when the potential $\bar{X}^2_t$ (or, respectively, the intensity process) is decreasing. On these intervals the differences between the three trajectories are negligible. However, the accuracy of $\tilde{f}_2$ drops drastically on the intervals where the intensity grows. In general, the higher the amplitude of the oscillations, the poorer the performance of the naive bound. This is particularly visible when illustrated on systems with high memory order ( $\eta_k=3$ or 6). For $\eta_k=1$ both bounds perform well; however, $\tilde{f}^{\tilde\Delta}_k$ is clearly closer to the true process. The influence of the bound ( $\tilde{f}_2$ or $\tilde{f}^{\tilde\Delta}_2$ ) on the execution time is discussed in Section 5.

Figure 3. Intensity and intensity bounds for the second population (excitatory), $t\in[20, 100]$ . Red solid line: true intensity $\lambda^2_t$ . Black dashed line: $\tilde{f}_2^{\tilde\Delta}(\bar X_t)$ . Grey dashed line: $\tilde{f}_2(\bar X_t)$ . Intensity functions $f_1$ and $f_2$ are given in (28); $\nu_1=\nu_2 = 0.9$ , $N_1 = N_2 = 50$ .

4.2. Simulation algorithm

Now let us detail the recursive procedure, which is summarized in Algorithm 1. We choose a discrete time step ${\tilde\Delta}$ and a stopping time $t_{\max}$ , and we fix the initial values $t_0 = 0$ and $\bar{X}_0 = x_0\in \mathbb{R}^\kappa$ . Let us assume that the procedure’s current step is i, with current time $t_i$ and potential value $\bar{X}_i$ . Let us explain how $t_{i+1}$ and $\bar{X}_{i+1}$ are obtained. One simulates two independent exponential variables $\tau_1$ and $\tau_2$ with respective parameters $N_k\tilde{f}_k^{\tilde\Delta}(x)$ (one for each population). They represent the waiting times to the next spikes of the dominant process $\tilde{Z}$ for each respective population. Then two cases may occur:

1. 1. If $\min\{\tau_1, \tau_2\}> {\tilde\Delta}$ , no spike occurs in the interval $[t_i, t_i+{\tilde\Delta}]$ . We update $t_{i+1} = t_i+{\tilde\Delta}$ and $\bar{X}_{i+1} = e^{A{\tilde\Delta}}\bar{X}_i$ .

2. 2. If $\tau = \min\{\tau_1, \tau_2\}\leq {\tilde\Delta}$ , then the dominating point process $\tilde{Z}$ emits a spike at time $t^* = t_i+\tau$ . Let us denote by $k^*$ the population with the smallest waiting time; that is, $\tau = \tau_{k^*}$ . It remains to decide whether $t^*$ is also a spiking time for the process Z. If not, this point is discarded. We draw a uniform variable U on [0,Reference Abi Jaber, Larsson and Pulido1] and define the threshold R:

   \begin{equation*}R \,{:}\,{\raise-1.5pt{=}}\, \frac{f_{k^*}\left(e^{A\tau}\bar{X}_{i} \right)}{\tilde{f}_{k^*}^{\tilde\Delta}(\bar{X}_{i})}; \qquad R \in [0,1] \text{ by the definition of }\tilde{f}_k^{\tilde\Delta}(x).\end{equation*}

   • • If $U\geq R$ , then $t^*$ is discarded, i.e., no spike occurs in the interval $[t_i, t^*]$ . We update $t_{i+1} = t^*$ and $\bar X_{i+1}=e^{A\tau}\bar X_{i}$ .

   • • If $U<R$ , then $t^*$ is kept, i.e., we add $t^*$ to the list of one neuron of population $k^*$ chosen uniformly at random. We update $t_{i+1} = t^*$ and $\bar X_{i+1} = e^{A\tau}\bar X_{i} + \Gamma \unicode{x1D7D9}(k^*)$ , where $\unicode{x1D7D9}(k^*) = (\unicode{x1D7D9}_{k^* = 1}, \unicode{x1D7D9}_{k^* = 2})^T$ .

Algorithm 1: Simulation of model (1) with K = 2 populations.

Finally, the execution is stopped once $t_{i}\geq t_{\max}$ , i.e. once the time horizon of interest is reached. As output the algorithm returns a list of the spiking times of each neuron and the values of the processes $\bar{X}$ and $\lambda$ at the spiking times. At this stage it is clear why it is important to have a sharp upper bound. The closer the threshold R is to 1, the fewer points are rejected.

Algorithm 1 is most efficient when every iteration of the while loop enters the condition (2). Of course, that ideal case does not occur in practice. When lowering the value of ${\tilde\Delta}$ , the number of loops satisfying the condition (3) decreases, because the dominating intensity $\tilde{\lambda}$ is getting smaller. On the other hand, the number of loops fulfilling the condition (1) increases, because the exponentially distributed times have greater chances of being larger than ${\tilde\Delta}$ . The calibration of $\tilde{\Delta}$ is a difficult problem which is not addressed here. In practice, it is observed that the execution time is not very sensitive to the value of $\tilde{\Delta}$ . The main bottleneck of the thinning algorithm is the sharpness of the intensity bound. When the intensity functions are exponential, the execution time is more than halved with the bound of Lemma 4 compared to the bound of Lemma 3. This is illustrated in the right panel of Figure 8.

5.1. Comparison of the Euler–Maruyama method and the splitting schemes

In this section we are interested in comparing the performance of the splitting schemes with that of the frequently applied Euler–Maruyama method, for varying time steps $\Delta$ . The parameter values $\nu_1 = \nu_2 = 1$ , $\eta_1 = 3$ , $\eta_2 = 2$ , $N_1 = N_2 = 50$ are used, and the dimension of the system is thus $\kappa = 7$ . Except for the density and mean-square convergence plots, we consider the time interval [0,100]. Unless stated otherwise, we plot the variables $X^{k,1}_t$ for $k=1,2$ in black and the remaining $\eta_1+\eta_2$ auxiliary memory variables in grey.

5.1.1 Illustration of the mean-square convergence order

We start our study by comparing the convergence rates of the Euler–Maruyama method and the Lie–Trotter (20) and Strang (21) splitting schemes. The root-mean-square error (RMSE), approximating the left side of the equation in Theorem 5, is defined as

\begin{equation*}\text{RMSE}(\Delta)\,{:}\,{\raise-1.5pt{=}}\, \left( \frac{1}{M} \sum_{l=1}^{M} \| X^{(l)}_{t^*} - \tilde{X}^{(l)}_{t^*} \|^2 \right)^{1/2},\end{equation*}

where $X^{(l)}_{t^*}$ and $\tilde{X}^{(l)}_{t^*}$ denote the values at a fixed time $t^*$ of the lth simulated trajectory of the true process and its numerical approximation, respectively. The integer M is the total number of simulated differences. The value of the true process $X^{(l)}_{t^*}$ is obtained from the Euler–Maruyama scheme, using the small time step $\Delta=10^{-5}$ . The number of simulations is fixed to $M=10^3$ and $t^*=1$ .

We report the RMSE in Figure 4, where the x-axis corresponds to the logarithm (base 10) of the time step $\Delta$ and the y-axis corresponds to the logarithm (base 10) of the RMSE. The theoretical rate of convergence obtained in Theorem 5 (all considered schemes converge with order 1) is confirmed empirically. While the Lie–Trotter splitting and the Euler–Maruyama scheme show a similar RMSE for varying $\Delta$ , the RMSE obtained for the Strang splitting method is significantly smaller for all $\Delta$ under consideration, implying higher efficiency for that scheme. We stress, however, that from the fact that the rate of convergence is the same, it does not follow that they share the same qualitative properties when the step size $\Delta$ is fixed.

Figure 4. Mean-square order of convergence. The reference solution is obtained with the Euler–Maruyama method and the small time step $\Delta=10^{-5}$ . The numerical solutions are calculated for $\Delta=10^{-4},10^{-3},10^{-2},10^{-1}$ . The time $t^*$ is 1, and $M=10^3$ .

5.1.2 Illustration of the qualitative properties of the splitting schemes

We now illustrate how the proposed splitting schemes preserve the structure (e.g. the moments and the underlying invariant distribution) of the process, even for large values of $\Delta$ , while the Euler–Maruyama method may fail to do so. We start by studying sample trajectories (see Figure 5). All three methods yield comparable performance when $\Delta=0.01$ . For $\Delta = 0.5$ , the Euler–Maruyama scheme preserves the oscillations, but does not preserve the amplitude. The behaviour of the inhibitory population is less accurately approximated than that of the excitatory one. This problem worsens as $\Delta$ increases further to $0.7$ . An interesting observation is that, for time steps $\Delta$ that are not ‘small enough’, the Euler–Maruyama scheme may not preserve the mean of the process (see also Figure 7). Indeed, it has been observed that the Euler–Maruyama method preserves the first but not the second moment (see e.g. [Reference Ableidinger, Buckwar and Hinterleitner3, Reference Higham and Strømmen Melbø20]). In other words, the amplitude of the oscillations grows, but the mean is unchanged. In our case, however, since the trajectories are bounded by 0 from below or above (depending on the sign of $c_k$ ), the increased amplitude also introduces a bias in the first moment. Thus, the Euler–Maruyama approximation of the system (7) preserves neither the first nor the second moment. In contrast, the Lie–Trotter and Strang splitting schemes show comparably good performance. However, the Lie–Trotter splitting is less accurate in reproducing the delay between the current state of the process (black line) and the memory variables (grey lines) in the beginning of the interval, where the amplitude of the oscillations is large (see also Figure 7).

Figure 5. Sample trajectories of the system, simulated with the Euler–Maruyama scheme (top), the Lie–Trotter splitting scheme (middle), and the Strang splitting scheme (bottom), for varying $\Delta$ .

The difference between the schemes becomes clearer as we look at the phase portrait of the system (Figure 6). We observe again that both splitting schemes yield satisfactory approximations (for all $\Delta$ under consideration), the Strang approach slightly outperforming the Lie–Trotter method. In contrast, the phase portrait obtained with the Euler–Maruyama approximation fails to reproduce the behaviour of the process for $\Delta = 0.5$ or $\Delta = 0.7$ .

Figure 6. Phase portrait of the main variables, simulated with the Euler–Maruyama scheme (top), the Lie–Trotter splitting scheme (middle), and the Strang splitting scheme (bottom), for varying $\Delta$ and $x_0=(0,0,-3.5,-4,0,1.3,1.1)$ .

Similar conclusions can be drawn from Figure 7, where we visualize the marginal densities of the process. Each visualized density is estimated with a standard kernel density estimator, based on a simulated long-time trajectory ( $T=10^5)$ of each variable of the process. We observe again that the Euler–Maruyama method may not preserve the mean of the process (red dashed vertical lines). Moreover, the Euler–Maruyama scheme may even suggest a transition from a unimodal to a bimodal density as $\Delta$ increases.

Figure 7. Empirical density of the system, simulated with the Euler–Maruyama scheme (top), the Lie–Trotter splitting scheme (middle), and the Strang splitting scheme (bottom), for varying $\Delta$ and $T = 10^5$ . The red dashed vertical lines denote the mean of the main variables.

5.2. Comparison of the PDMP and the diffusion

Now we are interested in comparing the PDMP process $\bar{X}$ , simulated with the thinning algorithm detailed in Section 4, and the stochastic diffusion X, simulated with the property-preserving Strang splitting scheme introduced in Section 3. We simulate the trajectories of the stochastic diffusion process with the Strang splitting scheme, since it showed the best performance in the previous section.

5.2.1 Execution time

As a first step we consider the execution time. We compare the numerical cost of the simulation of the process $\bar{X}$ with two different intensity bounds (based on Lemmas 3 and 4) to that of the simulation of the diffusion X with the Strang splitting scheme.

We set $t_{\max}=100$ and vary the total number of neurons, taking $N = 20, 50, 100, 150, 200$ and $N_1 = N_2$ . In the case of the diffusion simulation, the parameter N does not influence the computational cost. Thus, we report the execution time for the diffusion simulation only for $N=200$ , taking $\Delta = 0.1$ , and report it as a reference value. The time step $\tilde\Delta$ for the thinning procedure is defined in an adaptive way within the while loop of Algorithm 1. In each step we use the last computed value of the intensities $\lambda_k$ and set $\tilde\Delta$ equal to $(N_1\lambda_1+N_2\lambda_2)^{-1}$ . This choice takes into account the scaling with respect to the number of neurons and the dynamics of the intensities. For instance, $\bar{X}^{2,1}$ roughly belongs to [0, Reference Ableidinger and Buckwar2] (see Figure 10), so that the intensity of population roughly belongs to [Reference Abi Jaber, Larsson and Pulido1, Reference Buckwar, Tamborrino and Tubikanec7] (with the intensity functions defined in (28)).

In Figure 8, two different sets of intensity functions, ‘truncated linear’ ones and ‘exponential’ ones, are studied. The plots show, in seconds, the mean execution times (over 100 realizations) required to simulate the process on the interval $[0,t_{\max}]$ , using the bounds $\tilde{f}(x)$ and $\tilde{f}^\Delta(x)$ (defined in (26) and (27), respectively).

Figure 8. Mean execution time for the PDMP (solid line) and the stochastic diffusion (dashed line) simulation for $t_{\max}=100$ over 100 realizations. Left panel: $f_1(x) = f_2(x) = \min\{1+x\unicode{x1D7D9}_{[x>0]}, 10\}$ . Right panel: $f_1$ and $f_2$ are given by Equation (28). The rest of the parameters are given in the beginning of Section 5.

Note that there is almost no difference in the performance of the algorithm with different bounds in the linear case (left panel of Figure 8). That means that the bound obtained in Lemma 3 is sharp enough. Note also that since $f^{\max}_k = 10$ , there occur only a few spikes and the process is simulated very fast. However, in the case of an exponential intensity (right panel of Figure 8), the execution time drastically increases. The process is simulated at least twice as fast with the local bound. The main reason is that the local bound $\tilde{f}^\Delta_k(x)$ rejects around 2% of spikes, while the bound $\tilde{f}_k(x)$ rejects around 90%. In general, we can conclude that the execution time depends linearly on the number of neurons for both the local and the general bound. Disregarding the bound chosen, neither algorithm can compete with the time required for simulating the diffusion. For $\Delta = 0.1$ and $T=100$ , the average running time with the exponential firing rate function is equal to $0.598$ s (with standard deviation $0.12$ s). For the linear one it is $0.597$ s (with standard deviation $0.15$ s). Thus, the execution time for the diffusion approximation does not depend on the firing rates.

Figure 9. Mean execution time for the PDMP and the stochastic diffusion (dashed line) simulation for $t_{\max}=50$ over 100 realizations. Left panel: $f_1(x) = f_2(x) = \min\{a+x\unicode{x1D7D9}_{[x>0]}, 40\}$ . Right panel: $f_1(x) = f_2(x) = \min\{a+e^x, 40\}$ , where the parameter a varies between $0.02$ and 1 with a step of $0.02$ . The rest of the parameters are given in the beginning of Section 5; $N_1 = N_2 = 25$ .

Obviously, the execution time of the thinning algorithm depends on the number of spikes of the system per time unit, which is well approximated by the temporal mean intensity (formally defined by $T^{-1}\int_0^T \lambda^{k,n}_s ds$ ). Hence, for a fixed number of neurons N, a modification of the intensity functions $f_k$ changes the execution time. In Figure 9, we report the mean execution time over 100 repetitions of the experiment, plotted against the temporal mean value of the inhibitory population for a fixed number of neurons $N_1=N_2=25$ . For simplicity, we take the same intensity function for both populations, but we add a baseline intensity parameter a in order to control the temporal mean intensity: more precisely, $f_1(x) = f_2(x) = \min\{a+x\unicode{x1D7D9}_{[x>0]}, 40\}$ for the truncated linear case, and $f_1(x) = f_2(x) = \min\{a+e^x, 40\}$ for the truncated exponential case. As a increases, the temporal mean intensity increases and so does the execution time: the latter even increases linearly with respect to the mean intensity, whatever bound method is used for the algorithm. Comparing the two plots in Figure 9, it appears that the truncated exponential case is computationally more expensive (the values and the slopes are larger for both bounds). Once again, this comes from the fact that the proportion of rejected spikes is larger in the truncated exponential case. In both cases, the dependency on the temporal mean value is linear (which further confirms the results presented in Figure 8). Finally, notice that the execution time for the stochastic diffusion approximation does not depend on the mean intensity.

Figure 10. Sample trajectories of the PDMP and the diffusion for varying N (excitatory population). Solid line: main variable of $\bar{X}$ . Dashed line: main variable of X (simulated with the Strang splitting scheme, using $\Delta = 0.01$ ).

A summary of the performances of both frameworks (stochastic diffusion and PDMP), with respect to the parameters of the model, is given below.

• • In both cases, the execution time increases as the dimension of the system grows, i.e. as $\eta_k$ increases.

• • For the stochastic diffusion, the execution time depends, in a linear manner, on the step size $\Delta$ .

• • For the PDMP, the execution time mainly depends, in a linear manner, on the temporal mean intensities of the two populations. They in turn depend linearly on the number N of neurons. They also depend, in a complex nonlinear manner, on the parameters $\nu_k$ , $\eta_k$ , and $f_k$ .

• • The computational cost of the simulation of the stochastic diffusion does not depend on the number of neurons N or the temporal mean intensity. Hence, it should be preferred to the simulation of the PDMP unless the number of neurons N or the temporal mean intensity is small.

5.2.2 Qualitative properties

It remains to determine whether the stochastic diffusion can really capture the behaviour of the underlying PDMP. To get an intuitive idea of how different processes behave when the number of neurons changes, we look at some sample trajectories. We take one realization of the PDMP and the stochastic diffusion process on a time interval of length $T = 300$ and plot them in Figure 10, cutting the short interval of time in the very beginning in order to observe the process in its oscillatory regime. For simplicity, we focus only on the second (excitatory) population. The trajectories in the top panel are simulated with $N_2=10$ , those in the middle panel with $N_2=50$ , and those in the bottom panel with $N_2=100$ .

Let us mention that Figure 10 is not an illustration of Theorem 2. Indeed, the trajectories are not coupled in such a way that (12) is satisfied. To our knowledge, there is no such result in the literature, and the coupling involved in the proof of Theorem 2 is not explicit. However, the figure illustrates the fact that the fluctuations of both trajectories vanish as N goes to infinity and that both converge to the solution of the ODE (10).

As a final step, we are interested in the long-time behaviour of the processes. We simulate both processes ( $\bar{X}$ and X) on a long-time interval, taking $T = 10^5$ , and report the respective marginal empirical densities in Figure 11. The densities of the PDMP are plotted with solid lines and those of the diffusion with dashed lines. Even for small N, the difference between the densities is negligible and their means are almost overlapping. As the number of neurons N increases, we observe that the empirical densities converge to some compactly supported distribution. Note that the mean-field limit is given by the ODE (10) as illustrated in Figure 10. Thus we expect that the support of the limit distribution is given by the amplitude of the solution of the ODE.

Figure 11. Empirical density of PDMP and the diffusion for $N=20$ (left) and $N=100$ (right). Solid line: $\bar{X}$ . Dashed line: X (simulated with the splitting scheme, $\Delta = 0.1$ , $T=10^5$ ). The red solid and dashed vertical lines denote the mean of the respective main variables.