Mean field limits for nonlinear spatially extended hawkes processes with exponential memory kernels

TL;DR

This paper rigorously derives neural field equations as mean field limits of large spatially extended systems of nonlinear Hawkes processes, providing a mathematical foundation for neural modeling with exponential memory kernels.

Contribution

It introduces a rigorous derivation of neural field equations from large systems of nonlinear Hawkes processes with spatial extension and exponential memory kernels.

Findings

Proves convergence of neuron dynamics to a neural field equation as the number of neurons grows.

Establishes a mean field limit for spatially extended nonlinear Hawkes processes.

Provides a mathematical foundation for neural modeling using mean field theory.

Abstract

We consider spatially extended systems of interacting nonlinear Hawkes processes modeling large systems of neurons placed in Rd and study the associated mean field limits. As the total number of neurons tends to infinity, we prove that the evolution of a typical neuron, attached to a given spatial position, can be described by a nonlinear limit differential equation driven by a Poisson random measure. The limit process is described by a neural field equation. As a consequence, we provide a rigorous derivation of the neural field equation based on a thorough mean field analysis.

Equations

$$P\left(Z_i^{(N)}\text{mbox}hasajumpin]t,t+dt]|{\mathcal{F}}_t\right)={\lambda }_i^{(N)}(t)dt,$$

$$U_i^{(N)}(t):=e^{-\alpha t}{u}_0(x_i)+\frac{1}{N}\sum _{j=1}^{N}w\left(x_j,x_i\right){\int }_{]0,t]}{e}^{-\alpha (t-s)}d{Z}_j^{(N)}(s).$$

$$\{\begin{array}{l}\frac{\partial u(t,x)}{\partial t}=-\alpha u(t,x)+{\int }_{{\mathbb{R}}^d}w(y,x)f(u(t,y))\rho (dy),\\ u(0,x)=u_0(x).\end{array}$$

$$W_p({\nu }_1,{\nu }_2)=\underset{\pi \in \Pi ({\nu }_1,{\nu }_2)}{\inf }{\left({\int }_S{\int }_{S}d(x,y{)}^p\pi (dx,dy)\right)}^{1/p},$$

$$W_p({\nu }_1,{\nu }_2)=\inf \left\{E[d(X,Y{)}^p{]}^{1/p}:\mathcal{L}(X)={\nu }_1\text{mbox}and\mathcal{L}(Y)={\nu }_2\right\}.$$

$$\Vert h{\Vert }_{Lip}=\underset{x\ne y}{\sup }\frac{|h(x)-h(y)|}{d(x,y)}.$$

$$W_1({\nu }_1,{\nu }_2)=\underset{h\in Li{p}_1(S)\cap L^1\left(S,d|{\nu }_1-{\nu }_2|\right)}{\sup }\left\{{\int }_{S}h(x)({\nu }_1-{\nu }_2)(dx)\right\}.$$

$$U_i^{(N)}(t)=e^{-\alpha t}{u}_0(x_i)+\frac{1}{N}\sum _{j=1}^{N}w\left(x_j,x_i\right){\int }_{]0,t]}{e}^{-\alpha (t-s)}d{Z}_j^{(N)}(s).$$

$$d{U}_i^{(N)}(t)=-\alpha U_i^{(N)}(t)dt+\frac{1}{N}\sum _{j=1}^{N}w\left(x_j,x_i\right)d{Z}_j^{(N)}(t).$$

$$Z_i^{(N)}(t)={\int }_0^t{\int }_0^{\infty }1{}_{\left\{z\le f\left(U_i^{(N)}(s-)\right)\right\}}{\Pi }_i(dz,ds),$$

$$\frac{1}{N}\sum _{i=1}^{N}E\left[(Z_i^{(N)}(T))\right]\le T\left(f(0)+L_f\Vert u_0{\Vert }_{{\mathbb{R}}^d,\infty }\right)\exp \left\{T{L}_f\underset{j}{\sup }\Vert w_{x_j}{\Vert }_{L^1({\mu }^{(N)})}\right\},$$

$$\begin{gathered} \frac{1}{N}\sum _{i=1}^{N}E\left[(Z_i^{(N)}(T){)}^2\right]\le \exp \left\{T\left(1+4L_f^2{\left(\underset{j}{\sup }\Vert w_{x_j}{\Vert }_{L^2({\mu }^{(N)})}\right)}^2\right)\right\}\times \\ \times \left[T\left(f(0)+L_f\Vert u_0{\Vert }_{{\mathbb{R}}^d,\infty }\right)\exp \left\{T{L}_f\underset{j}{\sup }\Vert w_{x_j}{\Vert }_{L^1({\mu }^{(N)})}\right\}+2Tf(0{)}^2+4L_f^{2}T\Vert u_0{\Vert }_{{\mathbb{R}}^d,\infty }^2\right], \end{gathered}$$

$$|w(y,x)-w(y^{\prime },x^{\prime })|\le L_w(\Vert x-x^{\prime }\Vert +\Vert y-y^{\prime }\Vert ).$$

$$\underset{y\in {\mathbb{R}}^d}{\sup }\int |w(y,x){|}^2\rho (dx)<\infty ,$$

$$\underset{x\in {\mathbb{R}}^d}{\sup }\int |w(y,x){|}^2\rho (dy)<\infty .$$

$$(\max \{0,|w_{0^d}(x)|-L_w\Vert y\Vert \}{)}^2\le |w^x(y){|}^2\le 2\left(\Vert w_{0^d}{\Vert }_{{\mathbb{R}}^d,\infty }^2+L_w^2\Vert y{\Vert }^2\right).$$

$$\frac{|w_{0^{d}}(x)|^{2}}{4}\rho\left(B\Big{(}0^{d},\frac{|w_{0^{d}}(x)|}{2L_{w}}\Big{)}\right)\\ \leq\int|w^{x}(y)|^{2}\rho(dy)\leq 2\left(\|w_{0^{d}}\|^{2}_{{\mathbb{R}}^{d},\infty}+L^{2}_{w}\int\|y\|^{2}\rho(dy)\right),$$

$$\underset{x\in {\mathbb{R}}^d}{\sup }\int |w^x(y){|}^2\rho (dy)\le 2\left(\Vert w_{0^d}{\Vert }_{{\mathbb{R}}^d,\infty }^2+L_w^2\int \Vert y{\Vert }^2\rho (dy)\right)$$

$$\frac{\Vert w_{0^d}{\Vert }_{{\mathbb{R}}^d,\infty }^2}{4}\rho (B(0^d,\Vert w_{0^d}{\Vert }_{{\mathbb{R}}^d,\infty }/2L_w))\le \underset{x\in {\mathbb{R}}^d}{\sup }\int |w^x(y){|}^2\rho (dy).$$

$$P_{[0,T]}^{(N,N)}(d\eta ,dx)=\frac{1}{N}\sum _{i=1}^N{\delta }_{\left((Z_i^{(N)}(t){)}_{0\le t\le T},x_i\right)}(d\eta ,dx).$$

$$d_S(\eta ,\xi )=\underset{\varphi \in I}{\inf }\left\{\Vert \varphi {\Vert }_{[0,T],\ast }\vee \Vert \eta -\xi (\varphi ){\Vert }_{[0,T],\infty }\right\},$$

$$\Vert \varphi {\Vert }_{[0,T],\ast }=\underset{0\le s<t\le T}{\sup }\log \left(\frac{\varphi (t)-\varphi (s)}{t-s}\right).$$

$$d_{KR}(P,\tilde{P})=\underset{g\in Li{p}_1(D([0,T]\times \mathbb{N})\times {\mathbb{R}}^d)}{\sup }E\left[|⟨g,P-\tilde{P}⟩|\right],$$

$$d_{KR}(P_{[0,T]}^{(N,N)},P_{[0,T]})\to 0,$$

$${\int }_{{\mathbb{R}}^d}{\int }_{D([0,T],\mathbb{N})}[{\eta }_T^2+{\eta }_T]P_{[0,T]}(d\eta ,dx)<\infty .$$

$$P_{[0,T]}(d\eta ,dx)=P_{[0,T]}(d\eta |x)\rho (dx),$$

$$\lambda (t,x)=f\left(e^{-\alpha t}{u}_0(x)+{\int }_{{\mathbb{R}}^d}w(y,x){\int }_0^t{e}^{-\alpha (t-s)}\lambda (s,y)ds\rho (dy)\right).$$

$${\int }_{{\mathbb{R}}^d}\left[{\left({\int }_0^T\lambda (t,x)dt\right)}^2+{\int }_0^T\lambda (t,x)dt\right]\rho (dx)<\infty ,$$

$$|\lambda (t,x)-\lambda (t,y)|\le C\Vert x-y\Vert .$$

$$\lambda \mapsto F(\lambda )(t,x)=f\left(e^{-\alpha t}{u}_0(x)+{\int }_{{\mathbb{R}}^d}w(y,x){\int }_0^t{e}^{-\alpha (t-s)}\lambda (s,y)ds\rho (dy)\right),$$

Full text

60G55,

60G57,

60K35,

Hawkes processes,

Spatial mean field,

Propagation of Chaos,

Neural field equation,

Coupling,

, ,

and

t1

1 Introduction

The aim of this paper is to present a microscopic model describing a large network of spatially structured interacting neurons, and to study its large population limits. Each neuron is placed in ${\mathbb{R}}^{d}$. Its activity is represented by a point process accounting for the successive times at which the neuron emits an action potential, commonly referred to as a spike. The firing intensity of a neuron depends on the past history of the neuron. Moreover, this intensity is affected by the activity of other neurons in the network. Neurons interact mostly through chemical synapses. A spike of a pre-synaptic neuron leads to a change in the membrane potential of the post-synaptic neuron (namely an increase if the synapse is excitatory or a decrease if the synapse is inhibitory), possibly after some delay. In neurophysiological terms this is called synaptic integration. Thus, excitatory inputs from the neurons in the network increase the firing intensity, and inhibitory inputs decrease it. Hawkes processes provide good models for this synaptic integration mechanism by the structure of their intensity processes, see (1.1) below. We refer to [8], [9], [15] and to [19] for the use of Hawkes processes in neuronal modeling. For an overview of point processes used as stochastic models for interacting neurons both in discrete and in continuous time and related issues, see also [12].

In this paper, we study spatially structured systems of interacting Hawkes processes representing the time occurrences of action potentials of neurons. Each neuron is characterized by its spike train, and the whole system is described by the multivariate counting process $(Z^{(N)}_{1}(t),\ldots,Z^{(N)}_{N}(t))_{t\geq 0}$. Here, the integer $N\geq 1$ stands for the size of the neuronal network and $Z^{(N)}_{i}(t)$ represents the number of spikes of the $i$th neuron in the network during the time interval $[0,t].$ This neuron is placed in a position $x_{i}\in{\mathbb{R}}^{d},$ and we assume that the empirical distribution of the positions $\mu^{(N)}(dx):=N^{-1}\sum_{i=1}^{N}\delta_{x_{i}}(dx)$ converges111Convergence in the sense of the Wasserstein $W_{2}$-distance is considered in Scenario ($S_{2}$) below. to some probability measure $\rho(dx)$ on ${\mathbb{R}}^{d}$ as $N\to\infty$.

The multivariate counting process $(Z^{(N)}_{1}(t),\ldots,Z^{(N)}_{N}(t))_{t\geq 0}$ is characterized by its intensity process $(\lambda^{(N)}_{1}(t),\ldots,\lambda^{(N)}_{N}(t))_{t\geq 0}$ (informally) defined through the relation

[TABLE]

where ${\cal F}_{t}=\sigma(Z^{(N)}_{i}(s),\,s\leq t,{1\leq i\leq N}).$ We work with a spatially structured network of neurons in which $\lambda^{(N)}_{i}(t)$ is given by $\lambda^{(N)}_{i}(t)=f(U^{(N)}_{i}(t-))$ with

[TABLE]

Here, $f:{\mathbb{R}}\to{\mathbb{R}}_{+}$ is the firing rate function of each neuron and $w:{\mathbb{R}}^{d}\times{\mathbb{R}}^{d}\to{\mathbb{R}}$ is the matrix of *synaptic strengths *; for each $i,j\in\{1,\ldots,N\}$, the value $w(x_{j},x_{i})$ models the influence of neuron $j$ (located in position $x_{j}$) on neuron $i$ (in position $x_{i}$). The parameter $\alpha\geq 0$ is the leakage rate. Moreover, $u_{0}(x_{i})$ is the initial input to the membrane potential of neuron $i.$

Equation (1.1) has the typical form of the intensity of a multivariate nonlinear Hawkes process, going back to [17] and [16] who introduced Hawkes processes in a univariate and linear framework. We refer to [4] for the stability properties of multivariate nonlinear Hawkes processes, and to [10] and [6] for the study of Hawkes processes in high dimensions.

In this paper, we study the limit behavior of the system $(Z^{(N)}_{1}(t),\ldots,Z^{(N)}_{1}(t))_{t\geq 0}$ as $N\to\infty$. Our main result states that – under suitable regularity assumptions on the parameters $u_{0},w$ and $f$ – the system can be approximated by a system of inhomogeneous independent Poisson processes $(\bar{Z}_{x}(t))_{t\geq 0}$ associated with positions $x\in{\mathbb{R}}^{d}$ which can informally be described as follows. In the limit system, the spatial positions of the neurons are distributed according to the probability measure $\rho(dx)$. Given a position $x\in{\mathbb{R}}^{d}$, the law of the attached process $(\bar{Z}_{x}(t))_{t\geq 0}$ is the law of an inhomogeneous Poisson process having intensity given by $(\lambda(t,x))_{t\geq 0}$. Here $\lambda(t,x)=f(u(t,x))$ and $u(t,x)$ solves the scalar neural field equation

[TABLE]

Such scalar neural field equations (or neural field models) have been studied extensively in the literature, see e.g. [21], [22] and [1]. They constitute an important example of spatially structured neuronal networks with nonlocal interactions, see [5] for a recent and comprehensive review. Let us cite a remark made by Paul Bressloff, on page 15 of [5] :

There does not currently exist a multi-scale analysis of conductance-based neural networks that provides a rigorous derivation of neural field equations, although some progress has been made in this direction […]. One crucial step in the derivation of neural field equations presented here was the assumption of slowly varying synaptic currents, which is related to the assumption that there is not significant coherent activity at the level of individual spikes. This allowed us to treat the output of a neuron (or population of neurons) as an instantaneous firing rate. A more rigorous derivation would need to incorporate the mean field analysis of local populations of stochastic spiking neurons into a larger scale cortical model, and to carry out a systematic form of coarse graining or homogenization in order to generate a continuum neural field model.

Our model is not a conductance-based neural network. Nevertheless, with the present work, to the best of our knowledge, we present a first rigorous derivation of the well-known neural field equation as mean field limit of spatially structured Hawkes processes.

The paper is organized as follows. In Section 2 we introduce the model and provide an important a priori result, stated in Proposition 3, on the expected number of spikes of a typical neuron in the finite size system. In Section 3 we present our main results, Theorems 1 and 2, on the convergence of spatially extended nonlinear Hawkes processes towards the neural field equation (1.2). This convergence is expressed in terms of the empirical measure of the spike counting processes associated with the neurons as well as the empirical measure correspondent to their position. Therefore, we work with probability measures on the space $D([0,T],{\mathbb{N}})\times{\mathbb{R}}^{d}$ and a convenient distance defined on this space which is introduced in (3.8). The main ingredient of our approach is the unique existence of the limit process (or equivalently, of its intensity process), which is stated in Proposition 5. Once the unique existence of the limit process is granted, our proof makes use of a suitable coupling technique for jump processes which has already been applied in [10], [6] and [11], together with a control of the Wasserstein distance of order $2$ between the empirical distribution of the positions $\mu^{(N)}(dx)$ and the limit distribution $\rho(dx).$

As a by product of these convergence results, we obtain Corollaries 1 and 2, the former being closely connected to the classical propagation of chaos property and the latter stating the convergence of the process $(U^{(N)}_{1}(t),\ldots,U^{(N)}_{N}(t))_{t\geq 0}$ towards to the solution of the neural field equation.

In Section 4 are given the main technical estimates we use. Sections 5 and 6 are devoted to the proofs of our two main results, Theorem 1 and 2, together with Corollary 2. In Section 7 we discuss briefly the assumptions that are imposed on the parameters of the model by our approach. Finally, some auxiliary results are postponed to A.

2 General Notation, Model Definition and First Results

2.1 General notation

Let $(S,d)$ be a Polish metric space and $\mathcal{S}$ be the Borel $\sigma$-algebra of $S$. The supremum norm of any real-valued $\mathcal{S}$-measurable function $h$ defined on $S$ will be denoted by $\|h\|_{S,\infty}:=\sup_{x\in S}|h(x)|$. We will often write $\|h\|_{\infty}$ instead of $\|h\|_{S,\infty}$ when there is no risk of ambiguity. For any real-valued function $h(y,x)$ defined on $S\times S$ and $x\in S$, the $x$-section of $h$ is the function $h^{x}$ defined on $S$ by $h^{x}(y)=h(y,x)$ for all $y\in S$. Similarly, for any $y\in S$, the $y$-section of $h$ is the function $h_{y}$ defined on $S$ by $h_{y}(x)=h(y,x)$ for all $x\in S.$ The space of all continuous functions from $(S,d)$ to $({\mathbb{R}}_{+},|\cdot|)$ will be denoted by $C(S,{\mathbb{R}}_{+}).$ For any measure $\nu$ on $(S,\mathcal{S})$ and $\mathcal{S}$-measurable function $h:S\to{\mathbb{R}}$, we shall write $\left\langle h,\nu\right\rangle=\int_{S}h(x)\nu(dx)$ when convenient. For any $p\geq 1,$ we shall write $L^{p}(S,\nu)$ to denote the space of $\mathcal{S}$-measurable functions $h:S\to{\mathbb{R}}$ such that $\|h\|_{L^{p}(\nu)}:=(\int|h(x)|^{p}d\nu(x))^{1/p}<\infty$.

For two probability measures $\nu_{1}$ and $\nu_{2}$ on $(S,\mathcal{S})$, the Wasserstein distance of order $p$ between $\nu_{1}$ and $\nu_{2}$ associated with the metric $d$ is defined as

[TABLE]

where $\pi$ varies over the set $\Pi(\nu_{1},\nu_{2})$ of all probability measures on the product space $S\times S$ with marginals $\nu_{1}$ and $\nu_{2}$. Notice that the Wasserstein distance of order $p$ between $\nu_{1}$ and $\nu_{2}$ can be rewritten as the infimum of $E[d(X,Y)^{p}]^{1/p}$ over all possible couplings $(X,Y)$ of the random elements $X$ and $Y$ distributed according to $\nu_{1}$ and $\nu_{2}$ respectively, i.e.

[TABLE]

Let $Lip(S)$ denote the space of all real-valued Lipschitz functions on $S$ and

[TABLE]

We write $Lip_{1}(S)$ to denote the subset of $Lip(S)$ such that $\|h\|_{Lip}\leq 1.$ When $p=1,$ the Kantorovich-Rubinstein duality provides another useful representation for $W_{1}(\nu_{1},\nu_{2})$, namely

[TABLE]

Furthermore, the value of the supremum above is not changed if we impose the extra condition that $h$ is bounded, see e.g. Theorem 1.14 of [20].

2.2 The model and preliminary remarks

Throughout this paper we work on a filtered probability space $(\Omega,\mathcal{F},(\mathcal{F}_{t})_{t\geq 0},Q)$ which is rich enough such that all following processes may be defined on it. We consider a system of interacting nonlinear Hawkes processes which is spatially structured. In the sequel, the integer $N\geq 1$ will denote the number of processes in the system. Each process models the behavior of a specific neuron. All neurons are associated with a given spatial position belonging to ${\mathbb{R}}^{d}$. These spatial positions are denoted by $x_{1},\dots,x_{N}.$ In the following, ${\mathbb{R}}^{d}$ will be equipped with a fixed norm $\|\cdot\|.$ The positions $x_{1},\dots,x_{N}$ are assumed to be fixed in this section (unlike Section 5 where the positions are assumed to be random).

We now describe the dynamics of the $N$ Hawkes processes associated with these positions.

Definition 1**.**

Let $f:{\mathbb{R}}\to{\mathbb{R}}_{+}$, $w:{\mathbb{R}}^{d}\times{\mathbb{R}}^{d}\to{\mathbb{R}}$ and $u_{0}:{\mathbb{R}}^{d}\to{\mathbb{R}}$ be measurable functions and $\alpha\geq 0$ be a fixed parameter. A $(\mathcal{F}_{t})_{t\geq 0}$-adapted multivariate counting process $(Z_{1}^{(N)}(t),\ldots,Z_{1}^{(N)}(t))_{t\geq 0}$ defined on $(\Omega,\mathcal{F},(\mathcal{F}_{t})_{t\geq 0},Q)$ is said to be a multivariate Hawkes process with parameters $(N,f,w,u_{0},\alpha),$ associated with the positions $(x_{1},\ldots,x_{N}),$ if

$Q-$almost surely, for all pairs $i,j\in\{1,\dots,N\}$ with $i\neq j$, the counting processes $(Z^{(N)}_{i}(t))_{t\geq 0}$ and $(Z^{(N)}_{j}(t))_{t\geq 0}$ never jump simultaneously. 2. 2.

For each $i\in\{1,\dots,N\}$ and $t\geq 0$, the compensator of $Z^{(N)}_{i}(t)$ is given by $\int_{0}^{t}\lambda^{(N)}_{i}(s)ds$ where $(\lambda^{(N)}_{i}(t))_{t\geq 0}$ is the non-negative $(\mathcal{F}_{t})_{t\geq 0}-$progressively measurable process defined, for all $t\geq 0$, by $\lambda^{(N)}_{i}(t)=f(U^{(N)}_{i}(t-))$ with

[TABLE]

Remark 1**.**

Notice that for each $i\in\{1,\ldots,N\}$, the process $(U^{(N)}_{i}(t))_{t\geq 0}$ satisfies the following stochastic differential equation

[TABLE]

The functions $f:{\mathbb{R}}\to{\mathbb{R}}_{+}$ and $w:{\mathbb{R}}^{d}\times{\mathbb{R}}^{d}\to{\mathbb{R}}$ are called spike rate function and matrix of synaptic strengths respectively. The parameter $\alpha\geq 0$ is called the leakage rate. For each neuron $i\in\{1,\dots,N\}$, $u_{0}(x_{i})$ is interpreted as an initial input to the spike rate of neuron $i$ (its value depends on the position $x_{i}$ of the neuron). 222Without too much effort, the initial input $u_{0}(x_{i})$ could be replaced by a random input of the form $\frac{1}{N}\sum_{j=1}^{N}U_{i,j}$, where the random variables $U_{i,1},\ldots,U_{i,N}$ are i.i.d. distributed according to some probability measure $\nu(x_{i},du)$ defined on ${\mathbb{R}}^{d}$. In the limit $N\to+\infty$, we have the correspondence $\lim_{N\to+\infty}N^{-1}\sum_{j=1}^{N}U_{i,j}=\int u\nu(x_{i},du)=u_{0}(x_{i})$.

An alternative definition of multivariate Hawkes processes which will be used later on is the following.

Definition 2**.**

Let $(\Pi_{i}(dz,ds))_{1\leq i\leq N}$ be a sequence of i.i.d. Poisson random measures with intensity measure $dsdz$ on ${\mathbb{R}}_{+}\times{\mathbb{R}}_{+}$. A $(\mathcal{F}_{t})_{t\geq 0}$-adapted multivariate counting process $(Z_{1}^{(N)}(t),\ldots,Z_{N}^{(N)}(t))_{t\geq 0}$ defined on $(\Omega,\mathcal{F},(\mathcal{F}_{t})_{t\geq 0},Q)$ is said to be a multivariate Hawkes process with parameters $(N,f,w,u_{0},\alpha),$ associated with the positions $(x_{1},\ldots,x_{N}),$ if $Q-$almost surely, for all $t\geq 0$ and $i\in\{1,\ldots,N\}$,

[TABLE]

where $U^{N}_{i}(s)$ is given in (2.1).

Proposition 1**.**

Definitions 1 and 2 are equivalent.

We refer the reader to Proposition 3 of [10] for a proof of Proposition 1. In what follows we will work under

Assumption 1**.**

The function $f$ is Lipschitz continuous with Lipschitz norm $L_{f}>0.$

Assumption 2**.**

The initial condition $u_{0}$ is Lipschitz continuous with Lipschitz norm $L_{u_{0}}>0$ and bounded, i.e., $\|u_{0}\|_{{\mathbb{R}}^{d},\infty}<\infty.$

Proposition 2**.**

Under Assumption 1, there exists a path-wise unique multivariate Hawkes process with parameters $(N,f,w,u_{0},\alpha)$ such that $t\mapsto\sup_{1\leq i\leq N}E(Z_{i}^{(N)}(t))$ is locally bounded.

The proof of Proposition 2 relies on classical Picard iteration arguments and can be found (in a more general framework) in Theorem 6 of [10].

The next result provides an upper bound for the expected number of jumps in a finite time interval $[0,T]$, for each fixed $T>0$, which will be crucial later.

Proposition 3**.**

Under Assumptions 1 and 2, for each $N\geq 1$ and $T>0$, the following inequalities hold.

[TABLE]

and

[TABLE]

where $\mu^{(N)}(dx)=N^{-1}\sum_{i=1}^{N}\delta_{x_{i}}(dx)$ is the empirical distribution associated with the fixed spatial positions $x_{1},\ldots,x_{N}\in\mathbb{R}^{d}.$

The proof of Proposition 3 will be given in A.1.

3 Convergence of Spatially Extended Hawkes processes

In this section, we present two convergence results for the empirical process of nonlinear spatially extended Hawkes processes, our main results. We fix a probability measure $\rho(dx)$ on $({\mathbb{R}}^{d},{\cal B}({\mathbb{R}}^{d}))$; it is the expected limit of the empirical distribution $\mu^{(N)}(dx)=N^{-1}\sum_{i=1}^{N}\delta_{x_{i}}(dx).$ The following additional set of assumptions will be required as well.

Assumption 3**.**

The measure $\rho(dx)$ admits exponential moments, i.e., there exists a parameter $\beta>0$ such that $\mathcal{E}_{\beta}:=\int e^{\beta\|x\|}\rho(dx)<\infty$.

Assumption 4**.**

The matrix of synaptic strengths $w$ satisfies the following Lipschitz and boundedness conditions.

1. There exists a constant $L_{w}>0$ such that for all $x,x^{\prime},y,y^{\prime}\in{\mathbb{R}}^{d},$

[TABLE]

2. There exist $x_{0}$ and $y_{0}$ in $\mathbb{R}^{d}$ such that $||w(y_{0},\cdot)||_{{\mathbb{R}}^{d},\infty}<+\infty$ and $||w(\cdot,x_{0})||_{{\mathbb{R}}^{d},\infty}<+\infty$.

Remark 2**.**

Notice that the vectors $x_{0}$ and $y_{0}$ can be taken as the origin $0^{d}$ in $\mathbb{R}^{d}$. This follows directly from the Lipschitz continuity of $w$. Moreover, under Assumption 3, notice that Assumption 4-2. is equivalent to the fact that the $y$-sections $w_{y}$ and the $x$-sections $w^{x}$ of $w$ are uniformly square integrable with respect to $\rho(dx)$, i.e.,

[TABLE]

and

[TABLE]

In fact, (3.2) and (3.3) are the assumptions that naturally appear through the proofs (see e.g. (5.2) below where we need (3.2), and (4.10) where we need (3.3)).

Proof of the equivalence of Ass. 4-2. with (3.2) and (3.3), under Ass. 3 and (3.1).

We will only show that $\|w_{0^{d}}\|_{{\mathbb{R}}^{d},\infty}$ is finite if and only if (3.3) holds, the other case is treated similarly. Observe that the Lipschitz continuity of $w$ implies that for each $x,y\in{\mathbb{R}}^{d}$,

[TABLE]

As a consequence of the inequality above, it follows that

[TABLE]

where for each $z\in{\mathbb{R}}^{d}$ and $r>0$, $B(z,r)$ is the ball of radius $r$ centered at $z$. Now, taking the supremum with respect to $x,$ we deduce that

[TABLE]

and

[TABLE]

Under Assumption 3, we know that $\int\|y\|^{2}\rho(dy)$ is finite. Therefore, (3.4) shows that $\|w_{0^{d}}\|^{2}_{{\mathbb{R}}^{d},\infty}<\infty$ implies (3.3). On the other hand, (3.5) shows that (3.3) implies $\|w_{0^{d}}\|^{2}_{{\mathbb{R}}^{d},\infty}<\infty$ if $\rho(B(0^{d},\|w_{0^{d}}\|_{{\mathbb{R}}^{d},\infty}/2L_{w}))>0.$ Finally, if $\rho(B(0^{d},\|w_{0^{d}}\|_{{\mathbb{R}}^{d},\infty}/2L_{w}))=0,$ then it trivially holds that $\|w_{0^{d}}\|^{2}_{{\mathbb{R}}^{d},\infty}<\infty$ since $L_{w}>0$ and since $\rho(B(0^{d},\infty))=1.$ ∎

Let us now introduce, for each $T>0$ fixed, the empirical measure

[TABLE]

This measure $P^{(N,N)}_{[0,T]}$ is a random probability measure on the space $D([0,T],{\mathbb{N}})\times{\mathbb{R}}^{d},$ where $D([0,T],{\mathbb{N}})$ is the space of càdlàg functions defined on the interval $[0,T]$ taking values in ${\mathbb{N}}$. The measure is random since it depends on the realization of the $N$ counting processes $Z_{1}^{(N)},\ldots,Z_{N}^{(N)},$ defined on $(\Omega,\mathcal{F},(\mathcal{F}_{t})_{t\geq 0},Q)$.

Under two frameworks described below, our two main results state that $P^{(N,N)}_{[0,T]}$ converges as $N\to\infty$ to a deterministic probability measure $P_{[0,T]}$ defined on the space $D([0,T],{\mathbb{N}})\times{\mathbb{R}}^{d}$. This convergence is stated with respect to a convenient distance between random probability measures $P$ and $\tilde{P}$ on the space $D([0,T],{\mathbb{N}})\times{\mathbb{R}}^{d}$ which is introduced now.

For càdlàd functions $\eta,\xi\in D([0,T],{\mathbb{N}})$ we consider the distance $d_{S}(\eta,\xi)$ defined by

[TABLE]

where $I$ is the set of non-decreasing functions $\phi:[0,T]\to[0,T]$ satisfying $\phi(0)=0$ and $\phi(T)=T$ and where for any function $\phi\in I$ the norm $\|\phi\|_{[0,T],*}$ is defined as

[TABLE]

The metric $d_{S}(\cdot,\cdot)$ is equivalent to the classical Skorokhod distance. More importantly the metric space $(D([0,T],{\mathbb{N}}),d_{S})$ is Polish, see for instance [2].

Finally, for any random probability measures $P$ and $\tilde{P}$ on $D([0,T],{\mathbb{N}})\times{\mathbb{R}}^{d},$ we define the Kantorovich-Rubinstein like distance between $P$ and $\tilde{P}$ as

[TABLE]

where we recall that $\langle g,P-\tilde{P}\rangle=\int_{{\mathbb{R}}^{d}}\int_{D([0,T],{\mathbb{N}})}g(\eta,x)(P-\tilde{P})(d\eta,dx)$. Here the expectation is taken with respect to the probability measure $Q$ on $(\Omega,\mathcal{F},(\mathcal{F}_{t})_{t\geq 0})$, that is with respect to the randomness present in the jumps of the process.

Our convergence results are valid under two scenarios.

Definition 3**.**

Consider the two following assumptions:

[TABLE]

[TABLE]

Remark 3**.**

Under Scenario ($S_{1}$), the random positions $x_{1},\dots,x_{N}$ are interpreted as a random environment for our dynamics. This random environment is supposed to be independent of the Poisson random measures $\Pi_{i}(dz,ds),i\geq 1.$

In Scenario ($S_{2}$), the bound $N^{-1/d^{\prime}}$ is reasonable compared to the generic optimal quantization rate $N^{-1/d}$ [13, Theorem 6.2]. This bound is used afterwards to control the contribution of the spatial approximation in our spatial mean field approximation context. The construction of a sequence satisfying the claimed bound is described in Section 6.

In Theorems 1 and 2 below we will prove that as $N\to\infty$,

[TABLE]

almost surely with respect to the random environment under Scenario ($S_{1}$).

Let us now discuss briefly the properties that this limit $P_{[0,T]}$ should a priori satisfy. First, we obtain the following integrability property.

Proposition 4**.**

Under Assumptions 1, 2, 4 and either Scenario ($S_{1}$) or ($S_{2}$), the limit measure $P_{[0,T]}$ a priori satisfies that

[TABLE]

Proposition 4 is the limit version of Proposition 3, and its proof is postponed to A.2.

Let us now precisely define the limit measure $P_{[0,T]}.$ Firstly, consider any real-valued smooth test function $(\eta,x)\mapsto g(\eta,x)\equiv g(x)$ defined on $D([0,T],{\mathbb{N}})\times{\mathbb{R}}^{d}$ which does not depend on the variable $\eta$. Evaluating the integral of $g$ with respect to the probability measure $P^{(N,N)}_{[0,T]}$, and then letting $N\to\infty$, one deduces that the second marginal of $P_{[0,T]}$ on ${\mathbb{R}}^{d}$ must be equal to the probability measure $\rho(dx)$ (since $\mu^{(N)}(dx)$ converges to $\rho(dx)$).

Since $(D([0,T],{\mathbb{N}}),d_{S})$ is Polish, it follows from the Disintegration Theorem that $P_{[0,T]}$ can be rewritten as

[TABLE]

where $P_{[0,T]}(d\eta|x)$ denotes the conditional distribution of $P_{[0,T]}$ given the position $x\in{\mathbb{R}}^{d}$. Here, $x\mapsto P_{[0,T]}(d\eta|x)$ is Borel-measurable in the sense that $x\mapsto P_{[0,T]}(A|x)$ is Borel-measurable for any measurable subset $A\subset D([0,T],{\mathbb{N}}).$ From a heuristic which is explained below, the conditional distribution $P_{[0,T]}(d\eta|x)$ (for each $x\in\operatorname{Supp}(\rho)$, the support of $\rho$) turns out to be the law of a inhomogeneous Poisson point process with intensity process $(\lambda(t,x))_{0\leq t\leq T}$ where $\lambda=(\lambda(t,x),{0\leq t\leq T,x\in{\mathbb{R}}^{d}})$ is solution of the nonlinear equation

[TABLE]

The heuristic relies on the following argument: at the limit, we expect that the firing rate at time $t$ of the neurons near location $y$ should be approximately equal to $\lambda(t,y)$. Taking this into account, Equation (3.12) is the limit version of the interaction structure of our system described in Definition 1. In particular, the empirical mean with respect to the positions, i.e. the integral with respect to $\mu^{(N)}(dx)$, is replaced by an integral with respect to $\rho(dx)$.

Rewriting (3.10) in terms of the intensities, we obtain moreover that a priori,

[TABLE]

for each fixed $T>0.$ Existence and uniqueness of solutions for the nonlinear equation (3.12) is now ensured by

Proposition 5**.**

*Under Assumptions 1, 2 and 4, for any solution $\lambda$ of the equation (3.12) such that $t\mapsto\int_{{\mathbb{R}}^{d}}\int_{0}^{t}\lambda(s,y)ds\rho(dy)$ and $t\mapsto\int_{{\mathbb{R}}^{d}}\left(\int_{0}^{t}\lambda(s,y)ds\right)^{2}\rho(dy)$ are locally bounded, the following assertions hold.

1. For any $T>0,$ $\lambda\in C([0,T]\times{\mathbb{R}}^{d},{\mathbb{R}}_{+})$ and $\|\lambda\|_{{\mathbb{R}}^{d}\times[0,T],\infty}<\infty.$

$\lambda$ is Lipschitz-continuous in the space variable, that is, there exists a positive constant $C=C(f,u_{0},w,\alpha,\|\lambda\|_{[0,T]\times{\mathbb{R}}^{d},\infty},T)$ such that for all $x,y\in{\mathbb{R}}^{d}$ and for all $t\leq T,$*

[TABLE]

Furthermore, if we consider the map $F$ from $C([0,T]\times{\mathbb{R}}^{d},{\mathbb{R}}_{+})$ to itself defined by

[TABLE]

then for any $\lambda,\tilde{\lambda}\in C([0,T]\times{\mathbb{R}}^{d},{\mathbb{R}}_{+})$, the following inequality holds

[TABLE]

The proof of Proposition 5 is postponed to A.3. The inequality (3.14) together with a classical fixed point argument imply not only the existence but also the uniqueness of a solution of the equation (3.12).

We are now in position to state the two main results of the present paper.

Theorem 1**.**

Under Assumptions 1–4 and Scenario ($S_{1}$), there exists a positive constant $C=C(T,f,w,u_{0},\alpha,\beta)$ and a random variable $N_{0}$ depending only on the realization of the random positions $X_{1},\ldots,X_{N},\ldots,$ such that for all $N\geq N_{0}$,

[TABLE]

Moreover, if $\|w\|_{{\mathbb{R}}^{d}\times{\mathbb{R}}^{d},\infty}<\infty,$ then $N_{0}=1$ is valid. Furthermore, for any fixed $d^{\prime}>d,$ it holds that

[TABLE]

eventually almost surely as $N\to\infty.$

Theorem 2**.**

Under Assumptions 1–4 and Scenario ($S_{2}$), for each $T>0$, there exists a positive constant $C=C(T,f,w,u_{0},\alpha,\beta)$ such that for all $N\in{\mathbb{N}}$,

[TABLE]

Furthermore, for any fixed $d^{\prime}>2\vee d,$ it holds that

[TABLE]

eventually as $N\to\infty$.

The proofs of Theorems 1 and 2 are respectively postponed to Sections 5 and 6. Here are given two corollaries that are valid under either Scenario ($S_{1}$) or Scenario ($S_{2}$).

Corollary 1**.**

1. For any bounded functions $g,\tilde{g}\in Lip_{1}$,

[TABLE]

2. Suppose moreover that $\rho$ admits a $\mathcal{C}^{1}$ density, denoted by $f_{\rho}$, with respect to the Lebesgue measure. Then for any $1$-Lipschitz bounded functions $\varphi,\tilde{\varphi}:D([0,T],{\mathbb{N}})\to[-1,1]$ and any $x,\tilde{x}\in\mathbb{R}^{d}$ such that $f_{\rho}(x)\neq 0$ and $f_{\rho}(\tilde{x})\neq 0$, there exist kernel functions $\Phi_{N}(z),\tilde{\Phi}_{N}(z)$ with $\Phi_{N}(z)f_{\rho}(z)dz\stackrel{{\scriptstyle w}}{{\to}}\delta_{x}(dz)$ and $\tilde{\Phi}_{N}(z)f_{\rho}(z)dz\stackrel{{\scriptstyle w}}{{\to}}\delta_{\tilde{x}}(dz),$ as $N\to\infty,$ such that for

[TABLE]

we have

[TABLE]

Remark 4**.**

Equation (3.19) has to be compared to the property of propagation of chaos for standard or multi-class mean field approximations. Thanks to a suitable spatial scaling, which is contained in the test functions $g_{N}$ and $\tilde{g}_{N}$ and is explicit in the proof, we find that the activity near position $x$ is asymptotically independent of the activity near position $\tilde{x}$. Relating this result with multi-class propagation of chaos as defined in [14] for instance, let us mention that:

• •

if $x=\tilde{x}$, one recovers the chaoticity within a class,

• •

if $x\neq\tilde{x}$, one recovers the chaoticity between two different classes.

The functions $\Phi_{N}$ and $\tilde{\Phi}_{N}$ are approximations to the identity with spatial scaling $N^{p(d)}$. The optimal scaling obviously depends on the scenario under study. Under Scenario ($S_{1}$), we need $p(d)<[(4+d)(2d+1)]^{-1}$ whereas ($S_{2}$) the weaker condition $p(d)<[(2\vee d)(2d+1)]^{-1}$ is needed (the condition is weaker since the convergence is faster).

Proof of Corollary 1.

For any bounded functions $g,\tilde{g}\in Lip_{1}$, we apply the triangle inequality to obtain the following bound

[TABLE]

from which we deduce the first statement.

Next, to construct suitable sequences $(g_{N})_{N\geq 1}$ and $(\tilde{g}_{N})_{N\geq 1}$, let us first give some control for the density $f_{\rho}$. Since $f_{\rho}$ is $\mathcal{C}^{1}$, let $r>0$, $\varepsilon>0$ and $M>0$ be such that for all $y$ in $B(x,r)\cup B(\tilde{x},r)$, $f_{\rho}(y)\geq\varepsilon$ (recall that $f_{\rho}(x)\neq 0,$ $f_{\rho}(\tilde{x})\neq 0$) and $\|\nabla f_{\rho}(y)\|\leq M$.

Then, let $\Phi:\mathbb{R}^{d}\to\mathbb{R}$ be a mollifier, that is a compactly supported smooth function such that $\int\Phi(y)dy=1$. Let us define $\Phi_{N}$ and $\tilde{\Phi}_{N}$ by

[TABLE]

where $p(d)>0$ gives the spatial scaling of the approximations $\Phi_{N}$ and $\tilde{\Phi}_{N}$ (conditions on $p(d)$ raising to convergence are given in Remark 4 above and proved below). Since $\Phi$ is compactly supported, there exists $R>0$ such that $\operatorname{Supp}(\Phi)\subset B(0,R)$ which in particular implies that $\operatorname{Supp}(\Phi_{N})\subset B(x,RN^{-p(d)})$ and $\operatorname{Supp}(\tilde{\Phi}_{N})\subset B(\tilde{x},RN^{-p(d)})$. Until the end of the proof, we will assume that $N$ is large enough so that $RN^{-p(d)}\leq r$. Furthermore, we have

[TABLE]

and, by applying the quotient rule,

[TABLE]

We are now in position to define suitable sequences $(g_{N})_{N\geq 1}$ and $(\tilde{g}_{N})_{N\geq 1}$ by

[TABLE]

Obviously, the functions $\beta_{N}^{-1}g_{N}$ and $\beta_{N}^{-1}\tilde{g}_{N}$ belong to $Lip_{1}$ and $\max(\|g_{N}\|_{\infty},\|\tilde{g}_{N}\|_{\infty})\leq\alpha_{N}$.

On the one hand, applying the inequality obtained in the first step to $\beta_{N}^{-1}g_{N}$ and $\beta_{N}^{-1}\tilde{g}_{N}$, we deduce that

[TABLE]

is upperbounded by $2\beta_{N}\alpha_{N}d_{KR}\left(P^{(N,N)}_{[0,T]},P_{[0,T]}\right)$.

On the other hand, we have

[TABLE]

thanks to the continuity of $y\mapsto E[\varphi(\bar{Z}_{y})],$ which is a consequence of (4.13) proven below. Therefore,

[TABLE]

Gathering the steps above, we deduce (3.19) provided that $\beta_{N}\alpha_{N}d_{KR}\left(P^{(N,N)}_{[0,T]},P_{[0,T]}\right)$ goes to [math], which holds true if $p(d)<[(4+d)(2d+1)]^{-1}$ under Scenario ($S_{1}$) (apply Theorem 1) or if $p(d)<[(2\vee d)(2d+1)]^{-1}$ under Scenario ($S_{2}$) (apply Theorem 2). ∎

We close this section with the following observation. If for each $t\geq 0$ and $x\in{\mathbb{R}}^{d}$ we call

[TABLE]

then clearly $\lambda(t,x)=f(u(t,x))$ and $u(t,x)$ satisfies the scalar neural field equation

[TABLE]

Writing $U^{(N)}(t,x_{i}):=U^{(N)}_{i}(t),$ where $U^{(N)}_{i}(t)$ has been defined in (2.1) above, we obtain the convergence of $U^{(N)}(t,x_{i})$ to the solution $u(t,x)$ of the neural field equation in the following sense.

Corollary 2**.**

Under the conditions of either Theorem 1 or Theorem 2, we have that

[TABLE]

for any $T>0,$ where expectation is taken with respect to the randomness present in the jumps of the process.

The proof of this corollary goes along the lines of the proof of Theorem 1 and Theorem 2, in Sections 5 and 6 below.

4 Estimating $d_{KR}(P^{(N,N)}_{[0,T]},P_{[0,T]})$ for fixed positions $x_{1},\ldots,x_{N}$

Assume that the following quantities are given and fixed:

• •

the number of neurons $N$,

• •

the positions of the neurons $x_{1},\dots,x_{N}$.

Hereafter the following empirical measures will be used (the first and the last one are just reminders of (3.6) and (3.11) respectively),

[TABLE]

To estimate $d_{KR}(P^{(N,N)}_{[0,T]},P_{[0,T]})$ we shall proceed as follows. We will first show that $P^{(N,N)}_{[0,T]}$ and $P^{(\infty,N)}_{[0,T]}$ are close to each other by using a suitable coupling. The rate of convergence of such a coupling is a balance between the variance coming from the $N$ particles and the bias induced by the replacement of $\rho(dx)$ by $\mu^{(N)}(dx)$. Next, it will be shown that the $d_{KR}$-distance between $P^{(\infty,N)}_{[0,T]}$ and $P_{[0,T]}$ is controlled in terms of the Wasserstein distance between $\mu^{(N)}(dx)$ and $\rho(dx)$.

4.1 Estimating the $d_{KR}$-distance between $P^{(N,N)}_{[0,T]}$ and $P^{(\infty,N)}_{[0,T]}$

The aim of this subsection is to upper-bound $d_{KR}(P^{(N,N)}_{[0,T]},P^{(\infty,N)}_{[0,T]})$ when the positions $x_{1},\ldots,x_{N}\in{\mathbb{R}}^{d}$ are fixed.

Theorem 3**.**

Under Assumptions 1, 2 and 4, for each $N$ in ${\mathbb{N}}$ and $T>0$ there exists a constant $C=C(\alpha,f,w,T,u_{0})>0$ such that for a fixed choice $x_{1},\ldots,x_{N}\in{\mathbb{R}}^{d}$ of positions,

[TABLE]

Proof.

Fix a test function $g\in Lip_{1}$ and observe that by definition

[TABLE]

We will introduce in Equation (4.5) below a suitable coupling between the processes $(Z^{(N)}_{i}(t))_{0\leq t\leq T}$, $i=1,\dots,N$ and the processes $(\bar{Z}_{i}(t))_{0\leq t\leq T}$, $i\geq 1$, the latter being independent and distributed according to $P_{[0,T]}(d\eta|x_{i})$. In particular, we can decompose the equation above as

[TABLE]

with

[TABLE]

where $G_{i}:=g\left((\bar{Z}_{i}(t))_{0\leq t\leq T},x_{i}\right)$ for each $i\in\{1,\ldots,N\}$. To obtain the upper bound $A^{N}(T)$ we have used the $1$-Lipschitz continuity of $g$ and the inequality $d_{S}(\eta,\xi)\leq\sup_{t\leq T}|\eta(t)-\xi(t)|$ which is valid for all $\eta,\xi\in D([0,T],{\mathbb{N}}).$

Thus it suffices to obtain upper bounds for the expected values of $A^{N}(T)$ and $B^{N}(T)$. We start studying $B^{N}(T)$. By the Cauchy-Schwarz inequality and the independence of the $\bar{Z}_{i}$’s, it follows that

[TABLE]

so that we only need to control the variance of each $G_{i}$. Now, let $(\tilde{Z}_{i}(t))_{0\leq t\leq T}$ be an independent copy of $(\bar{Z}_{i}(t))_{0\leq t\leq T}$ and set $\tilde{G}_{i}=g\left((\tilde{Z}_{i}(t))_{0\leq t\leq T},x_{i}\right)$. In what follows, the expectation $\tilde{E}$ is taken with respect to $\tilde{G}_{i}$. Thus, by applying Jensen’s inequality we deduce that

[TABLE]

Then, since $d_{S}(\eta,\xi)\leq\sup_{0\leq t\leq T}|\eta(t)-\xi(t)|$ for all $\eta,\xi\in D([0,T],{\mathbb{N}})$ and both processes $(\bar{Z}_{i}(t))_{0\leq t\leq T}$ and $(\tilde{Z}_{i}(t))_{0\leq t\leq T}$ are increasing, the $1$-Lipschitz continuity of $g$ implies that

[TABLE]

Therefore, by applying once more Jensen’s inequality we obtain that

[TABLE]

Collecting all the estimates we then conclude that

[TABLE]

Now, noticing that $\bar{Z}_{i}(T)$ is a Poisson random variable with rate $\int_{0}^{T}\lambda(t,x_{i})dt,$

[TABLE]

Hence, we have just shown that

[TABLE]

Since clearly $(u+v)^{1/2}\leq u^{1/2}+v^{1/2}$ for all $u,v\geq 0$, it follows that

[TABLE]

We shall next deal with $A^{N}(T)$. Let us now introduce the coupling we consider here. Let $(\Pi_{i}(dz,ds))_{i\geq 1}$ be a sequence of i.i.d. Poisson random measures with intensity measure $dsdz$ on ${\mathbb{R}}_{+}\times{\mathbb{R}}_{+}$. By Proposition 1 the process $(Z^{(N)}_{i}(t))_{t\geq 0,i=1,\dots,N}$ defined for each $t\geq 0$ and $i\in\{1,\ldots,N\}$ by

[TABLE]

is also a multivariate nonlinear Hawkes process with parameters $(N,f,w,u_{0},\alpha),$ and the processes $(\bar{Z}_{i}(t))_{t\geq 0,i=1,\dots,N}$ defined for each $i\in\{1,\ldots,N\}$ and $t\geq 0$ as

[TABLE]

are independent and such that $(\bar{Z}_{i}(t))_{t\geq 0}$ is distributed according to $P(\cdot|x_{i})$ for each $i\in\{1,\ldots,N\}$. Now, for each $i$ and $t\geq 0$, let us define the following quantity

[TABLE]

Using successively that $f$ is Lipschitz and the triangle inequality we deduce that

[TABLE]

where $U^{(N)}(s,x_{i}):=U^{(N)}_{i}(s),$ with $U^{(N)}_{i}(s)$ as in (2.1), and

[TABLE]

Notice that since $e^{-\alpha(s-h)}\leq 1$ for $0\leq h\leq s$, we have

[TABLE]

which in turn implies that

[TABLE]

Now, write $W_{j}(t):=w\left(x_{j},x_{i}\right)\int_{[0,t[}e^{-\alpha(t-s)}d\bar{Z}_{j}(s)$. By the triangle inequality, we have

[TABLE]

Then the Cauchy-Schwarz inequality and the independence of the $W_{j}$’s implies that

[TABLE]

Now, $\textrm{Var}(W_{j}(t))=w(x_{j},x_{i})^{2}\int_{0}^{t}e^{-2\alpha(t-s)}\lambda(s,x_{j})ds,$ whence

[TABLE]

Applying Jensen’s inequality twice, the right-hand side of the inequality above can be upper-bounded by

[TABLE]

Since $\int_{0}^{T}\int_{0}^{t}e^{-2\alpha(t-s)}\lambda(s,x_{j})dsdt=\int_{0}^{T}\lambda(t,x_{j})(2\alpha)^{-1}(1-e^{-2\alpha(T-t)})dt\leq T\int_{0}^{T}\lambda(s,x_{j})dt$ (recall that $1-e^{-v}\leq v$ for all $v\geq 0$), we deduce that

[TABLE]

Finally we shall deal with $H^{(N)}_{i}(T)$. Proceeding similarly as above, we have

[TABLE]

and therefore it follows that

[TABLE]

Thus, defining $\delta(t)=N^{-1}\sum_{i=1}^{N}E[\Delta_{i}(t)]$ and then gathering the steps above gives

[TABLE]

Now, Gronwall’s lemma implies

[TABLE]

To deduce (4.2), it suffices to observe that $E[A^{N}(T)]\leq\delta(T)$ and the following control proven below: there exists a constant $C=C(\alpha,f,w,T,u_{0})>0$ such that

[TABLE]

Indeed, let $(\tilde{\Omega},\tilde{\cal A},\tilde{P})$ be a probability space on which are defined random variables $Y_{1}^{(N)}$ and $Y_{2}^{(N)}$ such that their joint law under $\tilde{P}$ is the optimal coupling achieving the Wasserstein distance $W_{2}(\mu^{(N)},\rho)$ of order $2$ between $\mu^{(N)}$ and $\rho.$ Then for all $t\leq T,$ using Hölder’s inequality,

[TABLE]

By Assumption 4,

[TABLE]

Moreover, using (3.13),

[TABLE]

implying (4.9). Finally, observe that together with (4.6), we obtain the same control for $E\left(\int_{\mathbb{R}}\int_{0}^{T}|U^{(N)}(t,x)-u(t,x)|dt\mu^{(N)}(dx)\right).$ ∎

4.2 Estimating the $d_{KR}$-distance between $P^{(\infty,N)}_{[0,T]}$ and $P_{[0,T]}$

In this subsection we give an upper bound for $d_{KR}(P^{(\infty,N)}_{[0,T]},P_{[0,T]})$ in terms of the Wasserstein distance between the empirical distribution $\mu^{(N)}(dx)$ and the limit measure $\rho(dx)$.

Proposition 6**.**

Grant Assumptions 1–4. For each $N\geq 1$ and $T>0$, there exists a positive constant $C=C(f,u_{0},w,\alpha,T)$ such that for any choice of $x_{1},\ldots,x_{N}\in{\mathbb{R}}^{d}$ the following inequality holds

[TABLE]

where $W_{1}(\mu^{(N)},\rho)$ is the Wasserstein distance between $\mu^{(N)}(dx)$ and $\rho(dx)$ associated with the metric $d(x,y)=\|x-y\|$ for $x,y\in{\mathbb{R}}^{d}.$

Proof.

Notice that for deterministic probability measures $P$ and $\tilde{P}$ which are defined on $D([0,T],{\mathbb{N}})\times{\mathbb{R}}^{d}$, the distance $d_{KR}\left(P,\tilde{P}\right)$ reduces to

[TABLE]

Now, fix a test function $g\in Lip_{1}.$ We take any coupling $W^{(N)}(dx,dy)$ of $\mu^{(N)}(dx)$ and $\rho(dy).$ Given $x$ and $y$ in $\mathbb{R}^{d}$, we use the canonical coupling of $\bar{Z}_{x}(t)$ and $\bar{Z}_{y}(t)$. That is, we pose

[TABLE]

and

[TABLE]

where we use the same Poisson random measure $\Pi(dz,ds)$ for the construction of $\bar{Z}_{x}(t)$ and of $\bar{Z}_{y}(t).$ Then

[TABLE]

since $\mu^{(N)}({\mathbb{R}}^{d})=\rho({\mathbb{R}}^{d})=1.$ By the Lipschitz continuity of $g,$

[TABLE]

Yet, as a consequence of the canonical coupling, it follows that

[TABLE]

where we used (3.13). Finally, the assertion follows from the definition of $W_{1}$. ∎

Before going to the proofs of Theorems 1 and 2, let us sum up the results obtained above and then present the scheme of the remaining proofs.

Remark 5**.**

Combining Theorem 3 and Proposition 6 together with the inequalities $W_{1}(\mu^{(N)},\rho)\leq W_{2}(\mu^{(N)},\rho)$ and $\|\cdot\|_{L^{1}(\mu^{(N)})}\leq\|\cdot\|_{L^{2}(\mu^{(N)})}$, it follows that to conclude the proofs of Theorems it suffices to control as $N\to+\infty$:

1- the Wasserstein distance $W_{2}(\mu^{(N)},\rho)$; 2- the supremum $\sup_{j}\|w_{x_{j}}\|_{L^{2}(\mu^{(N)})}$.

To treat the second point of the remark above, we use the following technical lemma in order to use (3.2).

Lemma 1**.**

Grant Assumption 4. Let $\mu$, $\rho$ be two probability measures. Assume that $\rho$ satisfies Assumption 3 and that $\mu$ is supported in $B(0^{d},r)$ for some $r>0$. Then, for any $\beta^{\prime}<\beta$, there exists a constant $C=C(w,\beta,\beta^{\prime})$ such that

[TABLE]

Proof.

The proof is based on a truncation argument. Let us define the auxiliary measure $\rho_{r}(dx)$ by

[TABLE]

which is the truncated version of $\rho(dx)$ to which we add a Dirac mass at $0^{d}$.

The Lipschitz continuity of $w$ implies that for any $y,x,z\in B(0^{d},r),$

[TABLE]

that is, $x\in B(0^{d},r)\mapsto(w_{y}(x))^{2}$ is Lipschitz with explicit constant, where the positive constant $C_{w}$ depends only on $w(0,0)$ and $L_{w}$. We hence deduce that for all $y\in B(0^{d},r),$

[TABLE]

thanks to the Kantorovich-Rubinstein duality. By the triangle inequality it follows then

[TABLE]

By the canonical coupling, we have

[TABLE]

Now, since for any $\beta^{\prime}<\beta$, there exists a constant $C>0$ depending only on $\beta$ and $\beta^{\prime}$ such that $\|x\|e^{\beta^{\prime}\|x\|}\leq Ce^{\beta\|x\|}$ we infer that $\int_{\|x\|>r}\|x\|\rho(dx)\leq C\mathcal{E}_{\beta}e^{-\beta^{\prime}r}$ which implies that $W_{1}(\rho,\rho_{r})\leq C\mathcal{E}_{\beta}e^{-\beta^{\prime}r}$.

Now, with $M_{r}=|w(0,0)|+L_{w}r$, we have for any $y\in B(0^{d},r)$,

[TABLE]

On the one hand, for any $\beta^{\prime}<\beta$, $\int_{||x||>r}||x||^{2}\rho(dx)\leq\int_{||x||>r}||x||^{2}e^{\beta^{\prime}||x||}\rho(dx)e^{-\beta^{\prime}r}$ and there exists a constant $C$ that only depends on $\beta$ and $\beta^{\prime}$ such that $||x||^{2}e^{\beta^{\prime}||x||}\leq Ce^{\beta||x||}$. Hence $\int_{||x||>r}||x||^{2}\rho(dx)\leq C\mathcal{E}_{\beta}e^{-\beta^{\prime}r}$. On the other hand, $1-\rho(B(0^{d},r))=\int_{||x||>r}\rho(dx)\leq\mathcal{E}_{\beta}e^{-\beta r}$. The same argument applies and gives the existence of a constant $C$ such that $\left(1-\rho(B(0^{d},r))\right)M_{r}^{2}\leq C\mathcal{E}_{\beta}e^{-\beta^{\prime}r}$. Finally, (4.14) follows from triangular inequality. ∎

5 Proof of Theorem 1

In this section we give the proof of Theorem 1 using the estimates obtained in the previous section. We work under the additional assumption that the positions $x_{1},\dots,x_{N}$ are realizations of i.i.d. random variables $X_{1},\ldots,X_{N},$ distributed according to $\rho.$

On the one hand, to control the Wasserstein distance, Theorem 1.6  Item (i) of [3] gives that for any fixed $d^{\prime}>d$ there exist constants $K$ and $N_{0}$ depending only on $d$, $\beta$ and $\mathcal{E}_{\beta}$ such that

[TABLE]

for any $0<\varepsilon\leq 1$ and $N\geq N_{0}\max\{1,\varepsilon^{-(4+d^{\prime})}\}$. As a consequence, for any fixed $d^{\prime}>d$ it follows that if $\varepsilon_{N}=O(N^{-\frac{1}{4+d^{\prime}}})$ then

[TABLE]

so that Borel-Cantelli’s lemma implies

[TABLE]

eventually almost surely.

On the other hand, define $R_{N}=N^{\gamma}$ where the constant $\gamma>0$ will be specified later and observe that by Markov’s inequality, $P\left(\cup_{i=1}^{N}\left\{\|X_{i}\|>R_{N}\right\}\right)\leq NP(\|X_{1}\|>R_{N})\leq\mathcal{E}_{\beta}e^{-\beta R_{N}(1-o(1))}$. As a consequence, $\sum_{N=1}^{\infty}P\left(\cup_{i=1}^{N}\left\{\|X_{i}\|>R_{N}\right\}\right)<\infty$ and by Borel-Cantelli’s lemma we deduce that for almost all realizations $X_{1},X_{2},\ldots$ there exists a $N_{0}$ depending on that realization such that for all $N\geq N_{0}$, $\mu^{(N)}$ is supported in $B(0^{d},R_{N})$.

Taking $d^{\prime}=d+1$, $0<\gamma<1/(10+2d)$ and applying Lemma 1 gives that almost surely,

[TABLE]

and in particular,

[TABLE]

Finally, the first two assertions of Theorem 1 follow immediately from (4.2) and (4.11) combined with (5.2) – or $\sup_{j}\|w_{x_{j}}\|_{L^{2}(\mu^{(N)})}\leq||w||_{\infty}$ in case $w$ is bounded; (3.16) is a consequence of (5.1).

6 Proof of Theorem 2

In this section we give the proof of Theorem 2 using the estimates obtained in Section 4. The main issue here is to construct a sequence of empirical distributions $\mu^{(N)}$ for which the Wasserstein distance $W_{2}(\mu^{(N)},\rho)$ is controlled.

We first use a truncation argument to reduce to the case where $\rho$ is compactly supported. Let $r>0$ be a truncation level to be chosen later. Let us define the measure $\rho_{r}$ by

[TABLE]

By the canonical coupling, we have

[TABLE]

Yet, if $r$ is large enough,

[TABLE]

hence, by Assumption 3, there exists $\beta>0$ such that

[TABLE]

Let us now describe how we construct an empirical distribution made of $N$ points adapted to any probability measure $\rho_{r}$ supported in $B(0^{d},r)$. For simplicity of our construction, we assume that $\mathbb{R}^{d}$ is endowed with the distance induced by the $\ell^{\infty}$ norm.

The construction is iterative, so we first explain how each step works. Let $\nu$ be a measure having support in the cube $B(0^{d},r)=[-r,r]^{d}$ and denote by $|\nu|\geq 1/N$ its mass. Let us prove that

[TABLE]

Assume that the $\nu$-mass of any cube of radius equal to $r\lfloor\left(N|\nu|\right)^{1/d}\rfloor^{-1}$ is less than $1/N$. There exists a covering of the cube $[-r,r]^{d}$ into $\lfloor\left(N|\nu|\right)^{1/d}\rfloor^{d}$ disjoint smaller cubes, each one of radius equal to $r\lfloor\left(N|\nu|\right)^{1/d}\rfloor^{-1}$. This implies

[TABLE]

yielding a contradiction. Then, treat separately the cases $\left(N|\nu|\right)^{1/d}\geq 2$ and $\left(N|\nu|\right)^{1/d}<2$ to prove $r\lfloor\left(N|\nu|\right)^{1/d}\rfloor^{-1}\leq 2r(N|\nu|)^{-1/d}$.

Applying the iterative step above to the probability measure $\rho_{r}$ gives the existence of a cube $C_{N}$ such that $\rho_{r}(C_{N})\geq 1/N$ and $\operatorname{Diam}(C_{N})\leq 4rN^{-1/d}$. Then, we define the measure

[TABLE]

Its mass is $1/N$. Applying the iterative step to $\tilde{\rho}_{r}=\rho_{r}-\rho_{r}^{N}$ (its mass is $(N-1)/N$) gives a cube $C_{N-1}$ such that $\tilde{\rho}_{r}(C_{N-1})\geq 1/N$ and $\operatorname{Diam}(C_{N-1})\leq 4r(N-1)^{-1/d}$. Similarly we define $\rho_{r}^{N-1}:=\frac{N^{-1}}{\tilde{\rho}_{r}(C_{N-1})}\tilde{\rho}_{r}\mathbf{1}_{C_{N-1}}$. In brief, applying $N$ times the iterative step gives a sequence of cubes $C_{1},\dots,C_{N}$ and associated measures $\rho_{r}^{1},\dots,\rho_{r}^{N}$ such that for all $k$,

[TABLE]

$\rho_{r}^{k}$ is a measure of mass $1/N$ supported in $C_{k}$, and $\rho_{r}=\sum_{k=1}^{N}\rho_{r}^{k}$.

For each $k$, let $x_{k}$ denote the center of $C_{k}$ and let $\mu^{(N)}=N^{-1}\sum_{k=1}^{N}\delta_{x_{k}}$ denote the associated empirical distribution. To control $W_{2}(\mu^{(N)},\rho_{r})$, we use the canonical coupling $\pi(dx,dy)=\sum_{k=1}^{N}\left(N^{-1}\delta_{x_{k}}\right)\otimes\rho_{r}^{k}$. Hence,

[TABLE]

If $d=1$, then $\sum_{k=1}^{+\infty}k^{-2/d}=\pi^{2}/6$ so $W_{2}(\mu^{(N)},\rho_{r})^{2}\leq g_{1}(r,N):=(4\pi^{2}r/6)N^{-1/2}$. If $d\geq 2$, by Hölder’s inequality333This is not optimal (see [7] for a refined version of this quantization argument).,

[TABLE]

so that

[TABLE]

We now chose a truncation level that depends on $N$, namely $r_{N}=N^{\varepsilon}$ for some $\varepsilon>0$. Combining (6.1) and the results above with the triangular inequality, we have, for $N$ large enough,

[TABLE]

Hence, for any $d^{\prime}>2\vee d$, there exist $K,N_{0}$ such that for all $N\geq N_{0}$,

[TABLE]

Taking $d^{\prime}=d+2$, $0<\varepsilon<1/(d+2)$ and applying Lemma 1 gives that

[TABLE]

and in particular,

[TABLE]

Finally, the first assertion of Theorem 2 follows immediately from (4.2) and (4.11) combined with (6.6); (3.18) is a consequence of (6.5).

Proof of Corollary 2.

Corollary 2 follows under both scenarios from (4.2) in Theorem 3, together with the arguments used to conclude the proofs of Theorem 1 and 2. ∎

7 Final discussion

In the previous sections, some technical assumptions have been imposed regarding the parameters of our model. These parameters are: the spike rate function $f$, the initial condition $u_{0}$ and the matrix of synaptic strengths $w.$ Here, we discuss these assumptions with respect to standard choices appearing e.g. in [5].

There are three main choices for the function $f$: a sigmoid-like function, a piecewise linear function or a Heaviside function (see page 6 of [5]). The first two choices obviously fit our Lipschitz condition (Assumption 1). A Heaviside function does of course not satisfy the Lipschitz condition. However, a Heaviside nonlinearity is less realistic and is mainly studied for purely mathematical reasons (to obtain explicit computations).

A typical choice for the function $u_{0}$ is a Gaussian kernel (see page 38 of [5]). It can describe an initial bump of the neural activity at some location of the cortex.

Usually, a (homogeneity) simplification is made concerning the function $w$: $w(y,x)$ is assumed to depend on $\|x-y\|$ only. Under this simplification, a common choice is the so-called Mexican hat function (see page 41 of [5]). Nevertheless, inhomogeneous neural fields where the previous simplification is dropped are also studied (see section 3.5 of [5] for instance). As a consequence of the modeling used in the present article, the interaction strength felt by neurons at position $x$ coming from neurons in the vicinity $dy$ of $y$ is given by $w(y,x)\rho(dy)$. Hence, inhomogeneity in neural networks can be considered in two ways:

• •

with an inhomogeneous matrix $w$ (whereas the limit spatial distribution $\rho$ is homogeneous),

• •

with an homogeneous $w$ but an inhomogeneous distribution $\rho$.

Let us discuss the spatial distribution $\rho$. Two standard choices are a uniform distribution over a bounded set (see [18]), or a finite sum of Dirac masses (in that case, the present paper is highly related to [11]). In these cases, $\rho$ is compactly supported and therefore satisfies Assumption 3. For unbounded domains, a typical choice is a Gaussian distribution satisfying Assumption 3 as well. Finally, let us notice that uniform distributions over unbounded domains (such as the Lebesgue measure on $\mathbb{R}^{d}$) are also considered for the study of the neural field equation. Such a model (which cannot correspond to a probability distribution on the positions) does not fit our assumptions and therefore cannot be obtained as the limit of a microscopic description following our approach.

Finally, the rate of convergence obtained in the present paper depends on the modeling scenario: the positions of the neurons are random or they are deterministic. Notice that our approach gives better rates of convergence in the deterministic framework (rate in $N^{-\frac{1}{2\vee d^{\prime}}}$ for any $d^{\prime}>d$) than in the random one (rate in $N^{-\frac{1}{4+d^{\prime}}}$ for any $d^{\prime}>d$).

Acknowledgements

This research has been conducted as part of the project Labex MME-DII (ANR11-LBX-0023-01); it is part of USP project Mathematics, computation, language and the brain, FAPESP project Research, Innovation and Dissemination Center for Neuromathematics (grant 2013/07699-0), CNPq projects Stochastic modeling of the brain activity (grant 480108/2012-9) and Plasticity in the brain after a brachial plexus lesion (grant 478537/2012-3). AD and GO are fully supported by a FAPESP fellowship (grants 2016/17791-9 and 2016/17789-4 respectively).

Appendix A Remaining Mathematical Proofs

A.1 Proof of Proposition 3

Using the inequality $f(u)\leq L_{f}|u|+f(0)$ valid for all $u\in{\mathbb{R}}$ (which follows from the Lipschitz continuity of $f$), we have that for each $1\leq i\leq N$,

[TABLE]

Thus, using that $e^{-\alpha(t-s)}\leq 1$ for all $0\leq s\leq t$, we obtain that

[TABLE]

Then, denoting $\beta(t)=N^{-1}\sum_{i=1}^{N}E\left[Z^{(N)}_{i}(t)\right]$ for each $t\geq 0$, it follows that

[TABLE]

Proposition 2 implies that $t\mapsto\beta(t)$ is locally bounded so that the first inequality stated in Proposition 3 follows from Gronwall’s inequality. We now turn to the control of the second moment of $Z^{(N)}_{i}(t).$ We first work with the stopped processes $Z_{i}^{(N)}(\cdot\wedge\tau_{K}),$ where $\tau_{K}=\inf\{t\geq 0:\sum_{i=1}^{N}Z_{i}^{(N)}(t)\geq K\},$ for some fixed truncation level $K>0.$ By Itô’s formula,

[TABLE]

(A.1) implies that

[TABLE]

Denoting $\gamma_{K}(t)=N^{-1}\sum_{i=1}^{N}E\left[(Z^{(N)}_{i}(t\wedge\tau_{K}))^{2}\right]$ for each $t\geq 0,$ it follows that

[TABLE]

This implies, applying once more Gronwall’s inequality, that

[TABLE]

We now obtain the result by letting $K\to\infty.$

A.2 Proof of Proposition 4

Taking the test function $g(\eta,x):=\eta_{T}$ which belongs to $Lip_{1},$ we obtain first for any $N,$

[TABLE]

and then, letting $N\to\infty$ and using (3.9),

[TABLE]

Yet

[TABLE]

by Proposition 3 together with (5.2) or (6.6) (depending on the chosen scenario).

To prove the second assertion, fix a truncation level $K>0$ and let $\Phi_{K}:{\mathbb{R}}_{+}\to{\mathbb{R}}_{+}$ be a smooth bounded function such that $\Phi_{K}(x)=x^{2}$ for all $x\leq K,$ $\Phi_{K}(x)\leq x^{2}$ for all $x\geq 0$ and such that $\|\Phi_{K}\|_{Lip}:=\sup_{x\neq y}\frac{|\Phi_{K}(x)-\Phi_{K}(y)|}{|x-y|}\leq 4K.$ Put then $g(\eta):=\frac{1}{4K}\Phi_{K}(\eta_{T}),$ by construction, this function belongs to $Lip_{1}.$ As before, we obtain that

[TABLE]

which implies, multiplying $g(\eta)$ by $4K,$ that

[TABLE]

where we have used that $\Phi_{K}(x)\leq x^{2}$ to obtain the rhs which does not depend on $K$ any more. As in the first step of the proof, the rhs of the above inequality is finite, thanks to Proposition 3 together with (5.2) or (6.6) (depending on the chosen scenario). Therefore, letting now $K\to\infty$ in the lhs of the above equation, the assertion follows.

A.3 Proof of Proposition 5

Using the inequality $f(u)\leq L_{f}|u|+f(0)$ valid for all $u\in{\mathbb{R}}$, one gets for each $0\leq t\leq T$,

[TABLE]

We use that $e^{-\alpha(t-s)}\leq 1$ for all $0\leq s\leq t,$ that $u_{0}$ is bounded and apply Hölder’s inequality to the last term (with respect to $\rho(dy)$) to obtain the upper bound

[TABLE]

Since $t\to\int_{{\mathbb{R}}^{d}}(\int_{0}^{t}\lambda(s,y)ds)^{2}\rho(dy)$ is locally bounded by assumption and since by Assumption 4 together with Remark 2, $x\mapsto\|w^{x}\|_{L^{2}(\rho)}$ is bounded, one obtains from the inequality above that $\lambda(t,x)$ is bounded on $[0,T]\times{\mathbb{R}}^{d},$ for any fixed $T>0.$

We now prove the continuity of the function $\lambda(t,x)$. To that end, take $t,t^{\prime}\in[0,T]$ with $t\leq t^{\prime}\leq T$ and $x,y$ in $\mathbb{R}^{d}$. On the one hand, by using successively the Lipschitz continuity of $f$, the triangle inequality and the boundedness of $u_{0},$ we deduce that

[TABLE]

Write $\|\lambda\|_{[0,T]\times{\mathbb{R}}^{d},\infty}:=\sup_{x\in{\mathbb{R}}^{d},t\leq T}\lambda(t,x)$ which is bounded thanks to the first step of the proof. It follows from the inequality above that

[TABLE]

On the other hand, the Lipschitz continuity of $f$ implies that

[TABLE]

Thus, using the Lipschitz-continuity of $u_{0},$ the Lipschitz-continuity of $w$ and the boundedness of $\lambda$, we deduce from the inequality above that

[TABLE]

Inequality (A.5) proves (3.13) and, together with (A.4), proves the continuity of $\lambda$.

Therefore it remains to establish (3.14). For that sake, observe that the Lipschitz continuity of $f$ implies that

[TABLE]

whence

[TABLE]

which implies the result.