where and are given constants, with a fixed parameter whereas is a scaling small parameter describing the intensity of local interactions between neurons. For all and for all , the connectivity kernel only depends on the relative distance between neurons and is given by
where . This scaling with respect to means that when goes to zero, space interactions are highly dominated by local ones compared to long range correlations. In the rest of this article, we assume that the connectivity kernel is nonnegative and rapidly vanishing at infinity, hence we introduce the following quantities,
which will play an important role later. A typical example for is a Gaussian function, or the indicator function in a compact set.
In , we proved that as the number of neurons goes to infinity and for , the set of neurons at time and position can be described by a distribution function solution to a mean-field equation,
with and given by
Here, we want to construct numerical solutions to (1.3)–(1.4) using particle methods, which consist in approximating the distribution function by a finite number of macro-particles. The trajectories of these particles are determined from the characteristic curves corresponding to the (1.3). Indeed, for any initial data
with finite second moments inand , the solution to (1.3)–(1.4) is uniquely defined as the push-forward of by the flow of the characteristic system of equations associated to (1.3)–(1.4), which can be written for and as
that is, for any test-function and ,
We also define for all and the following macroscopic quantities,
so that is the average neuron density in the network at time and location , and is the average pair membrane potential - adaptation variable. Therefore, we observe that may be written with respect to the macroscopic quantities and as
where denotes the standard convolution product in .
Before describing and analyzing a class of numerical methods for (1.3)–(1.4) in the presence of strong local space interactions (), we first briefly expound what may be expected from the continuous model in the limit .
On the one hand, by integrating (1.3) with respect to , we observe that for all
and moreover we suppose that it does not depend neither on , so that with
On the other hand, using (1.8), we observe that
Of course, this system is not closed since the right hand side of the equation on again depends on the distribution function . However, in the regime of strong local interactions , that is, in the limit , the singular term in indicates that the distribution function converges towards a Dirac distribution in centered in . Then applying a Taylor expansion of the solution , the right hand side of (1.11) gives rise to a diffusive operator for the spatial interactions at zeroth order with respect to . It yields that converges towards a limit pair satisfying the FHN reaction-diffusion system,
for more details on this asymptotic analysis.
We now come to our main concern in the present article and seek after a numerical method that is able to capture these expected asymptotic properties, even when numerical discretization parameters are kept independent of hence are not adapted to the stiffness degree of the space interactions. Our objective enters in the general framework of so-called Asymptotic Preserving (AP) schemes, first introduced and widely studied for dissipative systems as in [21, 23]. Yet, in opposition with collisional kinetic equations in hydrodynamic or diffusion limits, transport equations like (1.3) involve of course some stiffness in time but it is also crucial to take care of the space discretization in order to capture the correction terms of the non-local operator . By many respects this makes the identification of suitable schemes much more challenging.
In , the author proposed a numerical approximation to (1.3)–(1.4) using a standard particle method. However, as the parameter goes to , that is when the range of interactions between neurons shrinks and their amplitude grows, the time step and spatial grid size have to tend to zero too, hence the scheme cannot be consistent with the limit system (1.12) in the limit . In a different context [15, 16], F. Filbet & L. M. Rodrigues developed a particle method for the Vlasov-Poisson system with a strong external magnetic field, which is able to capture accurately the non stiff part of the evolution while allowing for coarse discretization parameters.
Here, we show how this approach may be extended to transport equations like (1.3) to deal with the time discretization. However, it is not sufficient since an appropriate space discretization technique is mandatory to capture the diffusive operator in (1.12) in the limit . In , the authors apply a spectral collocation method to provide numerical approximations of reaction-diffusion equations, with fractional spatial diffusion. Their method obviously can also be applied for local diffusions as in the FitzHugh-Nagumo reaction-diffusion system (1.12). On the other hand, the spectral collocation method also provides numerical approximations of differential equations with integral terms. For example, in [28, 11, 12, 13, 14] the authors use fast spectral methods for the non-local Boltzmann operator, which lead to compute the time evolution of Fourier coefficients of the solution instead of the solution itself. Therefore, this approach considerably simplifies the computation of the integral collision term and may be applied in our context. Moreover, we will show that a suitable formulation allows to perform a Taylor expansion of the solution in the Fourier space and to recover a consistent discretization of the macroscopic system (1.12) in the limit , which guarantee the asymptotic preserving property. Finally, another difficulty in our framework is to prove the convergence when vanishes of the nonlinear term in (1.5) involving the cubic function . The idea to circumvent this issue is to use, as in the continuous framework , the stiff term in (1.5), which stands for the interactions between neurons throughout the network to prove that the solution converges towards a Dirac mass in , that is all the membrane potential of the neurons at position are synchronized. Thus, it is possible to identify the asymptotic of the nonlinear term in (1.5). We will show that the particle approximation of the distribution in is particularly well suited to achieve this.
The rest of the paper is organized as follows. In Section 2, we present the particle method for the transport equation (1.3)–(1.4) and propose an apropriate time discretization technique in order to preserve the correct asymptotic when . Then, we provide first and second order schemes and verify the consistency when tends to zero. Finally, in Section 3, we present some numerical simulations to illustrate our results, and to study the dynamics of (1.3)–(1.4) with different different sets of parameters and different heterogeneous neuron densities.
2 A numerical scheme for the FitzHugh-Nagumo transport equation
This section is devoted to the construction of the numerical schemes for (1.3)- (1.4). We first focus on the discretization of the nonlocal operator in (1.4), for which we propose a spectral collocation method based on the discrete fast Fourier method. Then, we treat the transport equation (1.3) using a particle method for the microscopic variable and provide first and second order semi-implicit schemes for the time discretization. This algorithm is constructed in order to get a consistent approximation in the limit .
For sake of clarity, we drop the dependence with respect to on the distribution function and on the non-local operator .
2.1 Computation of the Non-local operator
We first look for an approximation of the operator given in (1.4). In view of applying a Fourier spectral method in space, we write as
Then we define a truncated operator in the following way.
Suppose that , where is the ball of radius centered at the origin and choose . Then, for any , is solution to
where for any ,
where denotes the characteristic function in the ball
denotes the characteristic function in the ball.
Actually the operator can be seen as convolution products between and the connectivity kernel , that is,
where is given by
In the sequel, we choose for simplicity such that , and consider a set of equidistant points with where is an even integer. An efficient strategy to approximate this nonlocal term is the spectral or spectral collocation methods [19, 28]. We suppose that the density and the macroscopic membrane potential are both known at the mesh points , then we compute an approximation of the Fourier coefficients for as,
and get a trigonometric polynomial
Therefore, we substitute this polynomials in (2.2), which yields a discrete operator given by
where is given by
and is the expansion coefficient depending on the connectivity kernel
Finally the approximation of the operator is provided by
Let us focus on the computation of the kernel modes for any fixed parameter . In the spirit of  for the Boltzmann equation, our purpose is to prove that these coefficients can be computed as one-dimensional integrals, so that we can store them in an array, but also to compute the asymptotic limit when in order to ensure that the scheme is consistent and stable when .
Using the change of variable , for and , we get:
Then, changing the variable into , we get
To complete the computation of the function , we have to study separately each possible value for the spatial dimension .
One-dimensional case: .
Since , it is straightforward to check that for any ,
hence we get:
Two-dimensional case: .
Let and . In this case, if we note , then using spherical coordinates, we have
where is the Bessel function of order , defined with
Consequently, we get
Three-dimensional case: .
Let and . Hence, if we note , and then using spherical coordinates, we get
where . Thus, the kernel mode is given by
Now let us investigate the asymptotic behavior of the discrete operator when . To this aim we set the space of trigonometric polynomial of degree in each direction, defined as 
equipped with the classical norm , which satisfies  for any
and for any and , we also have
Finally we define by the projection operator from to such that , for all .
Let and consider a connectivity kernel satisfying (1.2) with
Then, for all , there exists a positive constant , depending on , such that for all ,
Moreover for any trigonometric polynomial , we have
On the one hand, for any , we perform a Taylor expansion of at and using the assumptions (2.7) on , it yields
On the other hand, we have
Gathering these results and using (1.2), there exists a constant , depending on , such that
Then, we consider and for , we substitute the latter result in the expression (2.4) of , it yields for each
Thus, from the definition of (2.3), we know that and get
2.2 Particle/Spectral methods for (1.3)
We now consider the transport equation (1.3) and apply a standard particle method. This kind of numerical scheme was first introduced by Harlow  for the numerical computation of specific problems in fluid dynamics, and precisely mathematically studied later . Thus a large diversity of particle methods were developed for the simulation in fluid mechanics and plasma physics (see for instance [15, 16] and references therein). The method consists in approximating the solution to (1.3) with a sum of Dirac masses centered in a finite number of solutions of the characteristic system (1.5). These solutions stand for some particles characterized by a pair membrane potential-adaptation variable .
We approximate the solution to the transport equation (1.3) at each point ,
where , stands for the Dirac measure, and for any , is the solution of the spatially discretized characteristic system which can be written as follows, and
with a given initial data for and is the projection operator on . Moreover, we define the macroscopic potential at each point , as
2.3 Time discretization
The time discretization of the system (2.10) is the key point to get an Asymptotic-Preserving of the transport equation (1.3). We have to be especially careful about the stiff nonlocal terms in (1.5): on the one hand, we cannot use a fully explicit scheme, which does not provide an AP-scheme, and on the other hand, a fully implicit time discretization would be too costly because of the spectral collocation method for the nonlocal terms.
Therefore, our strategy consists in applying implicit-explicit numerical scheme, and to treat as an additional unknown of the system. In the following, we consider and for all , we set .
A first order semi-implicit scheme
We propose a first order semi-implicit scheme, that is for any time step and any particle index , we approximate solution to (2.10) by given by the following system
where denotes an approximation of the macroscopic membrane potential. Using the linearity of and the fact that , the system (2.12) yields that . Moreover, since the projection is linear, and according to its definition (2.3), we get that the right term in the first equation in (2.12) reads
On the one hand, let us emphasize that the stiff term, for , is treated implicitly but can be solved exactly whereas other terms, nonlinear with respect to , are considered explicitly. Formally speaking, when tends to zero, at each point , , the microscopic potential converges to .
On the other hand, the macroscopic membrane potential might be given by (2.11) from the values . Unfortunately, this approach would not give the correct asymptotic behavior of the macroscopic membrane potential when . Therefore, we consider as an additional variable solution to the following scheme
Observe here that the nonlinear term is computed implicitly from whereas the stiff term is now explicit.