Joachim Crevat
email: joachim.crevat@math.univ-toulouse.fr

Francis Filbet
email: francis.filbet@math.univ-toulouse.fr

Asymptotic preserving schemes for the FitzHugh-Nagumo transport equation with strong local interactions

Keywords:

This paper is devoted to the numerical approximation of the spatially extended FitzHugh-Nagumo transport equation with strong local interactions based on a particle method. In this regime, the time step can be subject to stability constraints related to the interaction kernel. To avoid this limitation, our approach is based on higher-order implicit-explicit numerical schemes. Thus, when the magnitude of the interactions becomes large, this method provides a consistent discretization of the macroscopic reaction-diffusion FitzHugh-Nagumo system. We carry out some theoretical proofs and perform several numerical experiments that establish a solid validation of the method and its underlying concepts.

Introduction

The FitzHugh-Nagumo (FHN) system [START_REF] Fitzhugh | Impulses and physiological sates in theoretical models of nerve membrane[END_REF], [START_REF] Nagumo | An active pulse transmission line simulating nerve axon[END_REF], models the pulse transmission in animal nerve axons and allows to describe complicated interactions of neurons in large neural networks. More precisely, we consider a network composed of n ∈ N neurons interacting with each other, where each one is labeled by i ∈ {1, ..., n}, and endowed with a parameter x i ∈ R d for d ∈ {1, 2, 3} standing for the constant spatial position in the network. The FHN system accounts for the variations of the membrane potential v i of a neuron coupled to an auxiliary variable w i called the adaptation variable. It can be written as follows for all i ∈ {1, ..., n},

             dv i dt = N (v i ) -w i + 1 n ε 2 n j=1 Ψ ε ( x i -x j ) (v j -v i ), dw i dt = τ (v i -γ w i ), (1.1)
where τ ≥ 0 and γ > 0 are given constants, N (v) = v (1 -v) (v -θ) with θ ∈ (0, 1) a fixed parameter whereas ε > 0 is a scaling small parameter describing the intensity of local interactions between neurons. For all ε > 0 and for all x ∈ R d , the connectivity kernel Ψ only depends on the relative distance between neurons and is given by

Ψ ε ( y ) := 1 ε d Ψ y ε , y ∈ R d ,
where Ψ : R + → R + . 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,

       Ψ := R d Ψ( y ) dy > 0, σ := 1 2 R d Ψ( y ) y 2 dy > 0, (1.2)
which will play an important role later. A typical example for Ψ is a Gaussian function, or the indicator function in a compact set.

In [START_REF] Crevat | Mean-field limit of a spatially-extended FitzHugh-Nagumo neural network[END_REF], we proved that as the number of neurons n goes to infinity and for Ψ ∈ Lip b (R + ), the set of neurons at time t > 0 and position x ∈ R d can be described by a distribution function f ε (t, x, .) solution of a mean-field equation,

     ∂ t f ε + ∂ v [f ε (N (v) -w + K ε [f ε ])] + ∂ w [f ε A(v, w)] = 0, f ε (t = 0, x, .) = f ε 0 (x, .), (1.3)
with K ε [f ε ] and A given by

       K ε [f ε ](t, x, v) = 1 ε 2 R d R 2 Ψ ε x -x (v -v) f ε (t, x , dv , dw )dx , A(v, w) = τ (v -γ w) .
(1.4)

Here, we want to construct numerical solutions of (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 f ε 0 with finite second moments in x ∈ R d and (v, w) ∈ R 2 , the solution of (1.3)-(1.4) is uniquely determined as the push-forward of f ε 0 by the flow of the characteristic system of equations associated to (1.3)-(1.4), which can be written for (t, x) ∈ R + × R d and (v, w) ∈ R 2 as

                 dV ε dt = N (V ε ) -W ε + K ε [f ε ] (t, x, V ε ) , dW ε dt = A (V ε , W ε ) , V ε (0) = v, W ε (0) = w.
(1.5)

Then we denote by Φ t,x the flow (v, w) ∈ R 2 → Φ t,x (v, w) = (V ε , W ε )(t, x, v, w), with Φ t,x (v, w) ∈ R 2 , hence the solution of (1.3)- (1.4) is given by f ε (t, x, .) = Φ t,x #f ε 0 (x, .), (1.6) where # is our notation for a push-forward, that is, for any test-function ϕ and B ⊂ R 2 , B ϕ(v, w) f ε (t, x, dv, dw) =

Φ -1 t,x (B)
ϕ • Φ t,x f ε 0 (x, dv, dw).

We also define for all (t, x) ∈ R + × R d and ε > 0 the following macroscopic quantities,

ρ ε    1 V ε W ε    (t, x) := R 2    1 v w    f ε (t, x, dv, dw), (1.7)
so that ρ ε (t, x) is the average neuron density in the network at time t and location x, and (V ε , W ε ) is the average pair membrane potential -adaptation variable. Therefore, we observe that K ε [f ε ] may be written with respect to the macroscopic quantities ρ ε and ρ ε V ε as

K ε [f ε ](., v) = 1 ε 2 [Ψ ε (ρ ε V ε ) -Ψ ε ρ ε v ] , (1.8)
where denotes the standard convolution product in x.

An important issue in the numerical simulation of (1.5) is that when the parameter ε is small, the numerical error of a classical time explicit scheme may become large. For instance with an explicit Euler scheme, the error may behave as O( ∆t ε 2 ), where ∆t is the time step, hence the scheme is not appropriate for ε 1. Here we want to design a numerical scheme which is less sensitive to this parameter ε > 0 in order to keep a control on the numerical error when ε 1.

Before describing and analyzing a class of numerical methods for (1.3)- (1.4) in the presence of strong local space interactions (ε 1), we first briefly expound what may be expected from the continuous model in the limit ε → 0.

On the one hand, by integrating (1.3) with respect to (v, w) ∈ R 2 , we observe that for all t ≥ 0

ρ ε (t, x) = ρ ε 0 (x), x ∈ R d
and moreover we suppose that it does not depend neither on ε, so that ρ ε (t, .) = ρ 0 with

ρ 0 ≥ 0, ρ 0 ∈ L ∞ (R d ).
(1.9)

On the other hand, using (1.8), we observe that

R 2 K ε [f ε ](t, x, v) f ε (t, x, dv dw) = ρ 0 ε 2 [ Ψ ε (ρ 0 V ε ) -(Ψ ε ρ 0 ) V ε ] . (1.10)
Hence, multiplying (1.3) by v (resp. w) and integrating with respect to (v, w) ∈ R 2 and using (1.10), we get a time evolution equation for the macroscopic quantities (ρ 0 V ε , ρ 0 W ε ) as

       ∂ t (ρ 0 V ε ) - ρ 0 ε 2 [Ψ ε (ρ 0 V ε ) -(Ψ ε ρ 0 )V ε ] = R 2 N (v)f ε (., dv, dw) -ρ 0 W ε , ∂ t (ρ 0 W ε ) = ρ 0 A (V ε , W ε ) . (1.11)
Of course, this system is not closed since the right hand side of the equation on ρ 0 V ε again depends on the distribution function f ε . However, in the regime of strong local interactions [START_REF] Crevat | Diffusive limit of a spatially-extended kinetic FitzHugh-Nagumo model[END_REF], that is, in the limit ε → 0, the singular term in ε -2 indicates that the distribution function f ε converges towards a Dirac distribution in v centered in V ε . Then applying a Taylor expansion of the solution V ε , 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 (ρ 0 V ε , ρ 0 W ε ) converges towards a limit pair (ρ 0 V, ρ 0 W ) satisfying the FHN reaction-diffusion system,

     ρ 0 ( ∂ t V -σ [∆ (ρ 0 V ) -V ∆ρ 0 ] -N (V ) + W ) = 0, ρ 0 ( ∂ t W -A(V, W ) ) = 0, (1.12)
where σ is defined in (1.2). We refer to [START_REF] Crevat | Diffusive limit of a spatially-extended kinetic FitzHugh-Nagumo model[END_REF] 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 [START_REF] Jin | Efficient asymptotic-preserving (AP) schemes for some multiscale kinetic equations[END_REF], [START_REF] Klar | An asymptotic-induced scheme for nonstationary transport equations in the diffusive limit[END_REF]. 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 K ε [f ε ]. By many respects this makes the identification of suitable schemes much more challenging.

One of the interest of the study of AP schemes is to numerically determine a rate of convergence of the transport equation (1.3) as the parameter ε goes to 0. Thus, we can compare this numerical rate of convergence with what we derived in the continuous framework in [START_REF] Crevat | Diffusive limit of a spatially-extended kinetic FitzHugh-Nagumo model[END_REF].

In [START_REF] Crevat | Mean-field limit of a spatially-extended FitzHugh-Nagumo neural network[END_REF], the author proposed a numerical approximation to (1.3)-(1.4) using a standard particle method. However, as the parameter ε goes to 0, 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 ε → 0. In a different context [START_REF] Filbet | Asymptotically stable particle-in-cell methods for the Vlasov-Poisson system with a strong external magnetic field[END_REF], [START_REF] Filbet | Asymptotically preserving particle-in-cell methods for inhomogeneous strongly magnetized plasmas[END_REF], 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 ε → 0. In [START_REF] Bueno-Orovio | Fourier spectral methods for fractional-in-space reactiondiffusion equations[END_REF], 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 [START_REF] Pareschi | Numerical Solution of the Boltzmann Equation I: Spectrally Accurate Approximation of the Collision Operator[END_REF], [START_REF] Filbet | High order numerical methods for the space non-homogeneous Boltzmann equation[END_REF], [START_REF] Filbet | Solving the Boltzmann equation in N log 2 N[END_REF], [START_REF] Filbet | A numerical scheme for the quantum Boltzmann equation with stiff collision terms[END_REF] and [START_REF] Filbet | On steady-state preserving spectral methods for homogeneous Boltzmann equations[END_REF], 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 ε → 0, 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 N . The idea to circumvent this issue is to use, as in the continuous framework [START_REF] Crevat | Diffusive limit of a spatially-extended kinetic FitzHugh-Nagumo model[END_REF], the stiff term in (1.5), which stands for the interactions between neurons throughout the network to prove that the solution f ε converges towards a Dirac mass in v, that is all the membrane potential of the neurons at position x 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 (v, w) ∈ R 2 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 appropriante time discretization technique in order to preserve the correct asymptotic when ε 1. 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.

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 K ε [f ε ] 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 (v, w) ∈ R 2 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 ε → 0.

For sake of clarity, we drop the dependence with respect to ε on the distribution function f ε and on the non-local operator K ε [f ε ].

Computation of the Non-local operator

We first look for an approximation of the operator K[f ] given in (1.4). In view of applying a Fourier spectral method in space, we write K[f ] as

K[f ](t, x, v) = 1 ε 2 R d Ψ ε ( y ) ρ 0 (x -y) (V (t, x -y) -v) dy .
Then we define a truncated operator K S [f ] in the following way.

Lemma 2.1. Suppose that Supp (ρ 0 ) ⊂ B(0, S), where B(0, S) is the ball of radius S > 0 centered at the origin and choose ε ∈ (0, 1). Then, for any

(t, x, v, w) ∈ R + × B(0, S) × R 2 , f is solution of ∂ t f + ∂ v f N (v) -w + K S [f ] + ∂ w [f A(v, w)] = 0,
where for any (x, v) ∈ B(0, S) × R,

K S [f ](t, x, v) = χ B(0,S) ε 2 (x) B(0,2S) Ψ ε ( y ) ρ 0 (x -y) (V (t, x -y) -v) dy, (2.1)
where χ B(0,S) denotes the characteristic function in the ball B(0, S).

Proof. On the one hand, since Supp (ρ 0 ) ⊂ B(0, S) and for all t ≥ 0, the density ρ(t) = ρ 0 , we get that for any x ∈ R d , the transport equation (1.3) can be written as

∂ t f + ∂ v f N (v) -w + χ B(0,S) K[f ] + ∂ w [f A(v, w)] = 0 .
Then it is enough to consider only x ∈ B(0, S). On the other hand, the domain of integration of the operator

K[f ] is such that y ≤ x + y -x ≤ 2 S, hence for any (t, x, v) ∈ R + × B(0, S) × R, K[f ](t, x, v) = 1 ε 2 B(0,2S) Ψ ε ( y ) ρ 0 (x -y) (V (t, x -y) -v) dy.
Thus, we define the truncated operator (2.1) as

ε 2 K S [f ] = χ B(0,S) K[f ].
Actually the operator K S [f ] can be seen as convolution products between (ρ 0 , ρ 0 V ) and the connectivity kernel Ψ ε , that is,

K S [f ](t, x, v) = 1 ε 2 L S [ρ 0 V ](t, x) -v L S [ρ 0 ](x) ,
where L S is given by

L S [u] = Ψ ε u, u ∈ {ρ 0 , ρ 0 V }. (2.2)
In the sequel, we choose for simplicity S = π/2 such that B(0, S) ⊂ T := [-π, π] d , and consider a set of equidistant points (x j ) j∈Jn x ⊂ T with J nx := -n x /2, n x /2 -1 d where n x is an even integer. An efficient strategy to approximate this nonlocal term is the spectral or spectral collocation methods [START_REF] Hesthaven | Spectral methods for time-dependent problems[END_REF][START_REF] Pareschi | Numerical Solution of the Boltzmann Equation I: Spectrally Accurate Approximation of the Collision Operator[END_REF]. We suppose that the density ρ 0 and the macroscopic membrane potential V are both known at the mesh points (x j ) j∈Jn x , then we compute an approximation of the Fourier coefficients for u ∈ {ρ 0 , ρ 0 V } as,

u(t, k) := 1 n d x j∈Jn x u(t, x j ) e -i k•x j , k ∈ J nx .
and get a trigonometric polynomial

u nx (t, x) := k∈Jn x u(t, k) e i k•x , u ∈ {ρ 0 , ρ 0 V } .
Therefore, we substitute this polynomials in (2.2), which yields a discrete operator L S nx given by

L S nx [u] := k∈Jn x L S [u](t, k) e i k•x , (2.3)
where L S [u] is given by

L S [u](t, k) = (2π) d Ψ ε (k) u(t, k) (2.4)
and Ψ ε is the expansion coefficient depending on the connectivity kernel

Ψ ε (k) = 1 (2π) d T Ψ ε ( x ) e -i k•x dx.
(2.5)

Finally the approximation K S nx [f ] of the operator K S [f ] is provided by K S nx [f ](t, x, v) = 1 ε 2 L S nx [ρ 0 V ](t, x) -v L S nx [ρ 0 ](x) . (2.6)
Let us focus on the computation of the kernel modes ( Ψ ε (k)) k∈Jn x for any fixed parameter ε > 0. In the spirit of [START_REF] Pareschi | Numerical Solution of the Boltzmann Equation I: Spectrally Accurate Approximation of the Collision Operator[END_REF] 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 ε → 0 in order to ensure that the scheme is consistent and stable when ε 1.

Using the change of variable x = r ω, for r ≥ 0 and ω ∈ S d-1 , we get:

Ψ ε (k) = 1 (2π) d π 0 Ψ ε (r) r d-1 I(k, r) dr, where I(k, r) := S d-1 exp(-i r k • ω) dω.
Then, changing the variable r into s = r/ε, we get

Ψ ε (k) = 1 (2π) d π/ε 0 Ψ(s) s d-1 I(k, ε s) ds.
To complete the computation of the function I, we have to study separately each possible value for the spatial dimension d ∈ {1, 2, 3}.

One-dimensional case: d = 1. Since S 0 = {-1, 1}, it is straightforward to check that for any k ∈ J nx , I(k, r) = 2 cos(r |k|), hence we get:

Ψ ε (k) = 1 π π/ε 0 Ψ(s) cos(ε s |k|) ds.
Two-dimensional case: d = 2. Let r ≥ 0 and k ∈ J nx . In this case, setting q = -r k, then using spherical coordinates, we have

I(k, r) = S 1 exp (i q • ω) dω = 2π 0 exp (i r k cos θ) dθ = 2 π 0 cos (r k sin θ) dθ = 2π J 0 (r k ),
where J 0 is the Bessel function of order 0, defined with

J 0 : x ∈ R → 1 π π 0 cos (x sin θ) dθ = ∞ l=0 (-1) l (l!) 2 x 2 2 l
.

Consequently, we get

Ψ ε (k) = 1 2π π/ε 0 Ψ (s) s J 0 (ε s k ) ds.
Three-dimensional case: d = 3. Let r ≥ 0 and k ∈ J nx . Hence, setting q = -r k, and then using spherical coordinates, we get

I(k, r) = S 2 exp (i q • ω) dω = 2π π 0 exp (i q cos(θ)) sin θ dθ = 2π 1 -1 exp (i q µ) dµ = 4π Sinc (r k ) ,
where Sinc(x) := sin(x)/x. Thus, the kernel mode Ψ ε (k) is given by

Ψ ε (k) = 1 2π 2 π/ε 0 Ψ (s) |s| 2 Sinc (ε s k ) ds.
Now let us investigate the asymptotic behavior of the discrete operator L S nx [u] when ε 1. To this approach we set S nx the space of trigonometric polynomial of degree n x /2 in each direction, defined as [START_REF] Canuto | Spectral Methods in Fluid Dynamics[END_REF]

S nx = span e ik•x , -n x /2 ≤ k j ≤ n x /2 -1, j = 1, . . . , d , equipped with the classical L 2 norm . L 2 , which satisfies for any u ∈ S nx u 2 L 2 = 2π n x d j∈Jn x |u(x j )| 2
and for any u and v ∈ S nx , we also have

T u(x) v(x) dx = 2π n x d j∈Jn x u(x j ) v(x j ).
Finally we define by I nx the projection operator from C(T) to S nx such that I nx u(x j ) = u(x j ), for all j ∈ J nx .

Proposition 2.2. Let d ∈ {1, 2, 3} and consider a connectivity kernel Ψ satisfying (1.2) with

R d Ψ( y ) y 4 dy < ∞. (2.7)
Then, for all k ∈ J nx , there exists a positive constant C > 0, depending on Ψ, such that for all ε > 0,

(2π) d Ψ ε (k) -Ψ + σ ε 2 k 2 ≤ C k 4 + 1 ε 4 . (2.8)
Moreover for any trigonometric polynomial u ∈ S nx , we have

L S nx [u] -Ψ u -σ ε 2 ∆u L 2 ≤ C ε 4 ∆ 2 u L 2 + u L 2 .
(2.9)

Proof. On the one hand, for any k ∈ J nx , we perform a Taylor expansion of I(k, .) at r = 0 and using the assumptions (2.7) on Ψ, it yields

(2π) d Ψ ε (k) - π/ε 0 Ψ(s)s d-1 ds -ε 2 k 2 π/ε 0 Ψ(s)s d+1 ds ≤ k 4 ε 4 R d y 4 Ψ( y )dy.
On the other hand, we have

∞ π/ε Ψ(s)s d-1 ds + ε 2 ∞ π/ε Ψ(s)s d+1 ds ≤ ε 4 1 π 4 + 1 π 2 R d y 4 Ψ( y )dy.
Gathering these results and using (1.2), there exists a constant C > 0, depending on Ψ, such that

(2π) d Ψ ε (k) -Ψ + σ ε 2 k 2 ≤ C ( k 4 + 1) ε 4 .
Then, we consider u ∈ S nx and for k ∈ J nx , we substitute the latter result in the expression (2.4) of L S [u](k), it yields for each

L S [u](k) -Ψ + σ ε 2 k 2 u(k) ≤ C ε 2 ( k 4 + 1) | u(k)|.
Thus, from the definition of (2.3), we know that L S nx [u] ∈ S nx and get

L S nx [u] -Ψ u -σ ε 2 ∆u L 2 =   k∈Jn x L S [u](k) -Ψ + σ ε 2 k 2 u(k) 2   1/2 , ≤ C ε 4 ( ∆ 2 u L 2 + u L 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 [START_REF] Harlow | The particle-in-cell computing method for fluid dynamics[END_REF] for the numerical computation of specific problems in fluid dynamics, and precisely mathematically studied later [START_REF] Raviart | An analysis of particle methods[END_REF]. Thus a large diversity of particle methods were developed for the simulation in fluid mechanics and plasma physics (see for instance [START_REF] Filbet | Asymptotically stable particle-in-cell methods for the Vlasov-Poisson system with a strong external magnetic field[END_REF], [START_REF] Filbet | Asymptotically preserving particle-in-cell methods for inhomogeneous strongly magnetized plasmas[END_REF] and references therein). The method consists in approximating the solution f 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 (v, w) ∈ R 2 .

In our case, since the transport equation (1.3) involves a high dimensional space (t, x, v, w) ∈ R + × R d × R 2 , a particle method seems to be the most natural approach. Moreover, although the cost of the simulation increases with the number of particles considered, this kind of method has already shown its efficiency to describe complex dynamics in plasma physics and fluid dynamics.

We approximate the solution f to the transport equation ( 1.3) at each point x ∈ T,

f M (t, x, dv, dw) := ρ 0 (x) M M p=1 δ Vp(t,x) (dv) ⊗ δ Wp(t,x) (dw),
where M ∈ N * , δ stands for the Dirac measure, and for any t ≥ 0, (V p , W p )(t) ∈ S nx is the solution of the spatially discretized characteristic system which can be written as follows, x ∈ T and 1

≤ p ≤ M          dV p dt = I nx N (V p ) + K S nx [f M ](V p ) -W p , dW p dt = A (V p , W p ) , (2.10)
with a given initial data (V 0 p , W 0 p ) ∈ S nx for 1 ≤ p ≤ M and I nx is the projection operator on S nx . Moreover, we define the macroscopic potential V M at each point (t, x) ∈ R + × T, as

                             ρ 0 = R 2 f M (t, x, dv, dw) , ρ 0 V M (t, x) := R 2 v f M (t, x, dv, dw) = 1 M M p=1 ρ 0 (x) V p (t, x) , ρ 0 W M (t, x) := R 2 w f M (t, x, dv, dw) = 1 M M p=1 ρ 0 (x) W p (t, x) .
(2.11)

From these macroscopic quantities, it is then possible to compute the discrete operator K S nx [f M ] given in (2.6), where (2.10)- (2.11) are solved at each mesh point (x j ) j∈Jn x .

Time discretization

The time discretization of (2.10) is the key point to get an asymptotic preserving scheme. The basic idea is do develop numerical methods that preserve the asymptotic limits (ε → 0) from the microscopic to the macroscopic models in the discrete setting. Contrary to multi-physics domain decomposition methods, the asymptotic preserving schemes only solve the microscopic equations avoiding the coupling of different models. This approach generates automatically macroscopic solvers when, in the asymptotic regime, the small time and space scales are not resolved numerically. This idea can be illustrated in Figure 2.1. Suppose, we start with a microscopic model P ε , which depends on a parameter ε, characterizing the small scale. As ε → 0, the model is approximated by a macroscopic model P 0 , which is independent of ε. We want to design a discretization P ε h of P ε , where h is the numerical parameter (mesh size and time step). If the asymptotic limit P ε h , as ε → 0 (with a fixed h), exists, then it is denoted by P 0 h . Furthermore when P 0 h is a stable and consistent approximation of P 0 , then the scheme P ε h is called asymptotic preserving. Error on an asymptotic preserving scheme is obtained from the following argument. Typically, we consider P ε and P 0 , corresponding for instance to (1.3)-(1.4) and its asymptotic model (1.12), we expect formally [START_REF] Crevat | Diffusive limit of a spatially-extended kinetic FitzHugh-Nagumo model[END_REF],

P ε P ε h P 0 P 0 h h → 0 ε → 0 h → 0 ε → 0
E 0 (ε) := P ε -P 0 = O(ε 2 ).
(2.12)

Then we assume that P ε h is an r-order approximation of P ε for a fixed ε > 0. Due to the presence of the small parameter ε > 0, a classical numerical analysis typically gives the following error estimates

E h (ε) := P ε h -P ε = O( h r /ε s ), s > 0, (2.13)
which blows-up when ε 1. Hence, the main issue of the asymptotic preserving analysis is to establish the discrete counterpart of the asymptotic analysis (2.12), that is, for a fixed h,

E h (ε) := P ε h -P 0 h = O(ε 2 ) (2.14)
and in the limit ε → 0,

E h (0) := P 0 h -P 0 = O( h r ). (2.15)
Clearly, if we add up the error estimates (2.12), (2.14) and (2.15), by the triangle inequality, we have

E h (ε) ≤ E 0 (ε) + E h (ε) + E h (0) = O(ε 2 + h r ). (2.16)
By comparing the two error estimates (2.13) and (2.16), it yields

E h (ε) ≤ C min h r ε s , ε 2 + h r ,
showing that when ε 1, the error does not blow-up. This formal argument applies to any asymptotic preserving schemes, although a rigorous proof will be problem dependent, based on the regularity of the solution to P ε and the specific scheme P ε h . Now, the aim is to apply this strategy to (2.10), which corresponds to the characteristic curves of (1.3)- (1.4). We have to be especially careful about the stiff nonlocal terms in (1.5), where the small parameter ε > 0 appears. We cannot use a fully explicit scheme, which does not provide an AP-scheme unless ∆t = O(ε 2 ), whereas 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 V M as an additional unknown of the system. In the following, we consider ∆t > 0 and for all n ∈ N, we set t n = n ∆t.

In this section, we propose a first and a second order time discretization scheme and prove some consistency properties when h is fixed and ε → 0 (see Proposition 2.3 and 2.5 ) corresponding to the error estimate E h (ε) in (2.14) and when ε is fixed and h → 0 (see Lemma 2.4 and 2.6) corresponding to the error estimates

E h (ε) in (2.

13).

A first order semi-implicit scheme

We propose a first order semi-implicit scheme, that is for any time step n ∈ N and any particle index 1 ≤ p ≤ M , we approximate (V p (t n ), W p (t n )) solution of (2.10) by (V n p , W n p ) ∈ S nx × S nx given by the following system

           V n+1 p -V n p ∆t -I nx N (V n p ) + 1 ε 2 L S nx [ρ 0 V n M ] -V n+1 p L S nx [ρ 0 ] + W n p = 0 , W n+1 p -W n p ∆t -A V n+1 p , W n p = 0 , (2.17)
where V n M denotes an approximation of the macroscopic membrane potential. Using the linearity of A and the fact that (V n p , W n p ) ∈ S nx × S nx , the system (2.17) yields that (V n+1 p , W n+1 p ) ∈ S nx × S nx . Moreover, since the projection I nx is linear, and L S nx [ρ 0 V n M ] ∈ S nx according to its definition (2.3), we get that the right term in the first equation in (2.17) reads

I nx N (V n p ) + 1 ε 2 L S nx [ρ 0 V n M ] -V n+1 p L S nx [ρ 0 ] -W n p = I nx N (V n p ) + 1 ε 2 L S nx [ρ 0 V n M ] -I nx V n+1 p L S nx [ρ 0 ] -W n p .
On the one hand, let us emphasize that the stiff term, for ε 1, is treated implicitly but can be solved exactly whereas other terms, nonlinear with respect to V p , are considered explicitly. Formally speaking, when ε tends to zero, at each point

x j ∈ T, j ∈ J nx , the microscopic potential V n+1 p converges to L S nx [ρ 0 V n M ]/L S nx [ρ 0 ]
. On the other hand, the macroscopic membrane potential V n M might be given by (2.11) from the values (V n p ) 1≤p≤M . Unfortunately, this approach would not give the correct asymptotic behavior of the macroscopic membrane potential when ε → 0. Indeed, as ε goes to 0, the first equation in (2.17) formally gives for all p and n:

V n+1 p ∼ V n M ,
which is the expected limit, but the non linear term N (V n p ) does not converge to N (V n M ). Therefore, we consider V n M ∈ S nx as an additional variable solution of the following scheme

V n+1 M -V n M ∆t - 1 M M p=1 I nx N (V n+1 p ) (2.18) - 1 ε 2 L S nx [ρ 0 V n M ] -I nx V n M L S nx [ρ 0 ] + W n M = 0 .
Observe here that the nonlinear term is computed implicitly from (V n+1 p ) 1≤p≤M whereas the stiff term is now explicit. Now, we define a numerical parameter h ∈ R 3 as h = (∆t, ∆x, 1/M ), where ∆x = 2π/n x and let us show the consistency of the numerical scheme (2.17)-(2.18) in the limit as ε → 0 for a fixed numerical parameter h. Proposition 2.3 (Consistency when ε → 0). Let h be a fixed parameter and consider a connectivity kernel Ψ : R + → R + satisfying (1.2), (2.7) and a neuron density ρ 0 ∈ S nx satisfying (1.9) at each grid point x j , j ∈ J nx . For all ε > 0, p ∈ {1, . . . , M } and n ∈ N, let us assume that the triplet (V ε,n p , W ε,n p , V ε,n M ) given by (2.17)-(2.18) is uniformly bounded with respect to ε > 0. Then we define

W ε,n M = 1 M M p=1 W ε,n p and for all j ∈ J nx , (V ε,n M , W ε,n M )(x j ) converges to (V n M , W n M )(x j ), as ε goes to 0, solution of                      V n+1 M -V n M ∆t = I nx N V n M -W n M + σ ∆I nx ρ 0 V n M -I nx V n M ∆ρ 0 , W n+1 M -W n M ∆t = A V n M , W n M .
(2.19)

Proof. For any p ∈ {1, . . . , M } and n ≥ 0, we denote by V ε,n p , W ε,n p , V ε,n M ε>0 the solution of (2.17)-(2.18) computed at the grid points (x j ) j∈Jn x . Since this sequence, abusively labeled by ε, is uniformly bounded, there exists a sub-sequence, still labeled in the same manner, which converges to

V n p , W n p , V n M when ε → 0.
On the one hand using the scheme (2.17) on V n+1 p , we may write

ε 2 V ε,n+1 p -V ε,n p ∆t = ε 2 I nx N (V ε,n p ) -ε 2 W ε,n p + L S ε,nx [ρ 0 V ε,n M ] -I nx V ε,n+1 p L S ε,nx [ρ 0 ] ,
and pass to the limit with respect to ε, it yields that for any j ∈ J nx ,

L S ε,nx [ρ 0 V ε,n M ](x j ) -I nx V ε,n+1 p L S ε,nx [ρ 0 ] (x j ) = L S ε,nx [ρ 0 V ε,n M ](x j ) -V ε,n+1 p (x i ) L S ε,nx [ρ 0 ](x j ) -→ ε→0 0.
Then, applying Proposition 2.2 to ρ 0 ∈ S nx , we have L S ε,nx [ρ 0 ] -Ψ ρ 0 L 2 → 0, when ε goes to 0, that is, for any j ∈ J nx L S ε,nx [ρ 0 ](x j ) -Ψ ρ 0 (x j ) -→ ε→0 0 .

Furthermore, applying again Proposition 2.2 to I nx (ρ 0 V ε,n M ) ∈ S nx , we also get

L S ε,nx [ρ 0 V ε,n M ](x j ) -Ψ ρ 0 V n M (x j ) -→ ε→0 0 ,
hence for any j ∈ J nx and p ∈ {1, . . . , M }, the limit V n+1 p (x j ) does no depend on p and is given by

V n+1 p (x j ) =    V n M (x j ), if ρ 0 (x j ) > 0, 0, else. Now we consider W ε,n M given by W ε,n M = 1 M M p=1 W ε,n p
and apply the second relation in (2.17), it gives by linearity of A,

W ε,n+1 M -W ε,n M ∆t = A   1 M M p=1 V ε,n+1 p , W ε,n M   ,
Passing to the limit ε → 0, we get an equation on the limit W n M given by

W n+1 M -W n M ∆t = A V n M , W n M .
On the other hand, we start from (2.18) and again apply Proposition 2.2, it yields that

V ε,n M L S ε,nx [ρ 0 ] -Ψρ 0 ε 2 -→ ε→0 σ V n M ∆ρ 0 , whereas L S ε,nx [ρ 0 V ε,n M ] -Ψ I nx ρ 0 V ε,n M ε 2 -→ ε→0 σ ∆I nx ρ 0 V n M .
Gathering these latter results, we get that when ε goes to zero,

L S ε,nx [ρ 0 V ε,n M ] -I nx V ε,n M L S ε,nx [ρ 0 ] ε 2 -→ ε→0 σ ∆I nx ρ 0 V n M -I nx V n M ∆ρ 0 . Therefore the limit V n+1 M is solution of V n+1 M -V n M ∆t = I nx N V n M -W n M + σ ∆I nx (ρ V n M ) -I nx ∆ρ 0 V n M . Finally, since the limit point (V n M , W n M ) is uniquely determined, actually all the sequence (V ε,n M , W ε,n M ) ε>0 converges.
Now, let us investigate the consistency error of the numerical scheme (2.17)-(2.18) as the parameter h goes to 0, for a fixed parameter ε like E h (ε) in (2.13). Let us note (V ε , W ε ) the flow of the characteristic system (1.5), and V ε the macroscopic potential as defined in (1.7). From the numerical scheme (2.17)-(2.18), we write the system in the form for all n ∈ N:

F 1 V ε,n p , V ε,n+1 p , W ε,n p , W ε,n+1 p , V ε,n M , V ε,n+1 M = 0,
where F 1 : R 6 → R 3 . We define the consistency error E n h (ε) as

E n h (ε) := F 1 V ε,n p , V ε,n+1 p , W ε,n p , W ε,n+1 p , V ε,n M , V ε,n+1 M L ∞ , (2.20)
where

. L ∞ is the classical L ∞ norm.
Lemma 2.4 (Consistency in the limit h → 0). Let 0 < ε < 1 be a fixed parameter and consider Ψ satisfying (1.2), and ρ 0 satisfying (1.9) and such that Supp(ρ 0 ) ⊂ B(0, S) where S > 0. We suppose that f ε the solution of (1.3)-(1.4) is differentiable twice with respect to time and there exists a constant C T > 0, independent of ε, such that in [0, T ] × B(0, S),

R 2 v 4 + w 4 f ε (., dv, dw) L ∞ + V ε L ∞ ≤ C T .
Consider the scheme (2.17)-(2.18) and the consistency error E n h (ε) in (2.20). Then, there exists another positive constant C > 0, independent of h and ε, such that for all n ∈ {0, . . . , [T /∆t]},

E n h (ε) ≤ C 1 + 1 ε 2 ∆x 3/2 ε 4 + ∆t + 1 M .
Proof. According to (1.3) and (1.5), we get that there exists a positive constant independent of ε such that for all t ∈ [0, T ] and z

= (x, v, w) ∈ R d+2 , ∂ 2 ∂t 2 V ε (t, z) + ∂ 2 ∂t 2 W ε (t, z) + ∂ 2 ∂t 2 (ρ 0 V ε )(t, x) ≤ C ε 4 .
Thus, from a Taylor expansion of the solution to (1.5) and the first equation of (1.11), the consistency error can be bounded as

E n h (ε) ≤ T 1 + T 2 + T 3 + T 4 + C ∆t ε 4 . with the non stiff terms                T 1 := A V ε (t n+1 ), W ε (t n ) -A (V ε (t n ), W ε (t n )) , T 2 := N (V ε (t n )) -I nx (N (V ε (t n ))) , T 3 := ρ 0 (x) M M p=1 I nx N (V ε (t n+1 )) - R 2 N (V ε (t n )) f ε 0 (dv , dw ) ,
and the stiff term

T 4 := 1 ε 2 I nx L S nx [ρ 0 ] V ε (t n+1 ) -Ψ ε [ρ 0 ] V ε (t n ) + 1 ε 2 I nx L S nx [ρ 0 V ε (t n )] -Ψ ε [ρ 0 V ε ] + 1 ε 2 I nx L S nx [ρ 0 ] V ε (t n ) -Ψ ε [ρ 0 ] V ε (t n ) .
First of all, for the first term, since second and fourth order moments of f ε in v and w are uniformly bounded with respect to ε, we directly get that

T 1 ≤ C ∆t,
for some positive constant C > 0. Then, in order to treat the second term, we use the estimate from Theorem 2.12 in [START_REF] Hesthaven | Spectral methods for time-dependent problems[END_REF] which yields that for all U ∈ H 2 (B(0, S)), there exists a constant C > 0 such that

U -I nx (U ) L ∞ (B(0,S)) ≤ C U H 2 (B(0,S)) ∆x 3/2 .
Since the initial data are regular enough, we get that there exists a positive constant independent of ε such that for all t ∈ [0, T ] and

z = (x, v, w) ∈ B(0, S) × R 2 , ∇ 2 x V ε (t, z) + ∇ 2 x W ε (t, z) + ∇ 2 x (ρ 0 V ε )(t, x) ≤ C ε 4 .
This leads to the estimate

T 2 ≤ C ∆x 3/2 ε 4 ,
for some positive constant C > 0. As for the third term, we decompose it as follows:

T 3 ≤ ρ 0 M M p=1 |I nx N (V ε (t n+1 ) -N (V ε (t n+1 ))| + ρ 0 M M p=1 |N (V ε (t n+1 ) -N (V ε (t n )| + ρ 0 M M p=1 N (V ε (t n ) - R 2 N (V ε (t n )) f ε 0 (dv dw ) .
The first term can be treated as previously with the estimate from [START_REF] Hesthaven | Spectral methods for time-dependent problems[END_REF]. Then, the second term is of order ∆t. Finally, the third term corresponds to an approximation of the integral term with the rectangle rule, hence it is of order 1/M . Consequently, we get that

T 3 ≤ C ∆x 3/2 ε 4 + ∆t + 1 M ,
for some positive constant C. As for the final term T 4 , we notice that for all U ∈ C(B(0, S)),

L S nx [U ] = I nx [Ψ ε U ] .
Therefore, we want to use the estimate form [START_REF] Hesthaven | Spectral methods for time-dependent problems[END_REF] here again. This gives us that there exists a positive constant C such that

T 4 ≤ C ε 2 ∆x 3/2 ε 4 + ∆t + 1 M .
Consequently, gathering the previous results, we get that there exists a positive constant C independent of h and ε such that if h is small enough,

E n h (ε) ≤ C 1 + 1 ε 2 ∆x 3/2 ε 4 + ∆t + 1 M .
The lack of uniform bounds, with respect to ε and h, on the numerical solution does not allow us to complete the rigorous analysis of the asymptotic preserving scheme (2.17)-(2.18). Furthermore, Proposition 2.3 only gives a convergence result without any error estimate as E h (ε) in (2.14). However, Proposition 2.3 indicates that in the limit ε → 0, the numerical scheme (2.17)-(2.18) becomes a first order explicit time approximation with respect to ∆t of the reaction-diffusion system (1.12), hence it constitutes together with Lemma 2.4 a first direction to provide a proof of the asymptotic preserving property (2.16).

A second order implicit-explicit Runge-Kutta scheme Now let us adapt the previous strategy to a second order implicit-explicit Runge-Kutta scheme for the system (2.10). We propose a combination of Heun's method for the explicit part, and an A-stable second

Numerical simulations

In this section, we provide examples of numerical computations to validate and compare the different time discretization schemes (2.17)-(2.18) and (2.21)-(2.25) introduced in the previous section.

First of all, we focus on the order of accuracy when ε is fixed and the numerical parameter h goes to zero. Then we study the behaviour of the numerical solutions for a fixed h and in the limit ε → 0, to show the convergence towards the solutions of the approximations (2.19) and (2.26)-(2.28) of the reaction-diffusion system (1.12).

Then, we display some simulations of the behaviour of a solution of (1.3)-(1.4) with a heterogeneous neuron density, and finally, we show some two-dimensional dynamics.

Throughout this section, except for first subsection, we fix the parameter of the nonlinearity N to θ = 0.1 and the other constants to τ = 0.005 and γ = 5, expect in the first subsection. This framework corresponds to the "excitable" regime of the well-known FHN reaction-diffusion system (3.5). Therefore, the system only admits one steady state which is the stable fixed point 0, and according to [START_REF] Carpenter | A geometric approach to singular perturbation problems with application to nerve impulse equations[END_REF], τ is small enough so that the solution of (3.5) exhibits slow/fast dynamics like traveling pulses.

Moreover, as for the connectivity kernel, we use the following truncated gaussian function

Ψ( z ) = 1 (2 π σ 0 ) d/2 exp - z 2 2 σ 0 , (3.1)
with σ 0 = 0.005 such that we have in (1.2), Ψ = 1 and σ = σ 0 2 .

Order of accuracy in the numerical parameters

In this subsection, we aim to verify the order of accuracy of our numerical methods proposed in Section 2 with respect to the numerical parameters h = (∆t, ∆x, 1/M ), when it goes to zero. We consider a simplified version of the nonlocal transport equation (1.3) with N (v) = -αv and τ = 0, that is, for t > 0 and

x ∈ R      ∂ t f ε + ∂ v (f ε [-α v -w + K ε [f ε ]]) = 0, f ε | t=0 (x, v, w) = δ V 0 (x) (v) ⊗ δ 0 (w), (3.2)
with V 0 given by

V 0 (x) = exp -100 |x| 2 , x ∈ R.
Consequently, in this configuration, we have ρ 0 ≡ 1, and the solution of (1.3)-(1.4) is given by f ε = δ V ε (v) ⊗ δ 0 (w) where V ε is the unique solution of the following reaction-diffusion equation for t > 0 and Thus, the parameter α > 0 determines the rate of convergence of V ε towards the stable state 0. Since the term N (V ε ) is now linear, the macroscopic equation on V ε is also linear (even if the equation on f ε is not) and we can exhibit an explicit solution using a Fourier transform in space. It yields that,

x ∈ R,        ∂ t V ε - 1 ε 2 Ψ ε V ε -Ψ V ε = -α V ε , V ε (0, x) = V 0 (x), x ∈ R. h L 2 error
V ε (t, ξ) = V 0 (ξ) exp -α + 1 ε 2 Ψ ε (ξ) -Ψ t .
where we choose the parameter α = 0.001, and the scaling parameter ε = 1. The domain in space is taken to be (-1, 1). We compute an approximation of the error on the macroscopic quantity V ε at each time step ). A linear regression yields that these numerical methods seems to be respectively first and second order in h. Therefore, with this parametrization, the order of accuracy corresponds to the one given by the time discretization, whereas the error due to the spectral discretization is negligible.

E n = V ε,n M -V ε (t n ) L 2 , n = 0, . . . N T , with N T = [T /∆t].

Order of accuracy in ε

We again consider the transport equation (1.3)-(1.4), with the initial data

f 0 (x, v, w) = δ V 0 (x) (v) ⊗ δ W 0 (x) (w), (3.3)
with V 0 = χ [-1,1] and W 0 ≡ 0.

In this configuration, we get ρ 0 ≡ 1 and the solution of the transport equation (1.3) is again a Dirac mass in (v, w) centered in (V ε , W ε ), solution of the nonlocal reaction-diffusion system for t > 0 and

x ∈ R d ,        ∂ t V ε - 1 ε 2 Ψ ε V ε -Ψ V ε = N (V ε ) -W ε , ∂ t W ε = A (V ε , W ε ) . (3.4)
The purpose is now to study the asymptotic when the scaling parameter ε goes to 0. It is expected that the macroscopic quantities (V ε , W ε ) converge towards the solution of the reaction-diffusion FHN system (1.12), which reads as follows when ρ 0 ≡ 1, for t > 0 and x ∈ R,

     ∂ t V -σ ∂ 2 x V = N (V ) -W, ∂ t W = τ (V -γ W ) . (3.5)
To investigate this asymptotic, we compute an approximation of the relative entropy given at any time t > 0 as

D ε (t) := R ρ 0 (x) |V ε (t, x) -V (t, x)| 2 + |W ε (t, x) -W (t, x)| 2 dx 1/2 , (3.6)
as ε goes to 0.

Here again, we approach D ε (t) with a rectangle rule. In [START_REF] Crevat | Diffusive limit of a spatially-extended kinetic FitzHugh-Nagumo model[END_REF], it is proven that for any t > 0, D ε (t) tends to 0 as ε goes to 0 with a rate of convergence larger than 2/7. However, when ρ 0 ≡ 1 and for compactly supported f ε , the rate of convergence is formally equal to 2.

Furthermore, since the solution of the transport equation is a Dirac mass in (v, w) ∈ R 2 , we take M = 1. Then, we choose ∆t = 0.01 and n x = 512 for the time and space discretization.

In Figure 3.1, we show the spatio-temporal profile of the mean membrane potential V ε computed from f ε the solution of the transport equation (1.3) for ε = 3.25 (panel (a)), ε = 3 (panel (b)) and ε = 1 (panel (c)). It shows that depending on the value of ε, the solution V ε presents dramatically different dynamics.

If ε is too large compared to the width of the considered interval, as in the case (a), two symmetric waves start to propagate, but quickly disappear, and then V ε converges to 0 everywhere as time goes on. On the contrary, for smaller values of ε as in the cases (b) and (c), that is ε ≤ 3, the function V ε has the shape of two symmetric counter-propagating traveling pulses. This is typically the kind of slow/fast dynamics expected for the solution of (3.5) according to [START_REF] Carpenter | A geometric approach to singular perturbation problems with application to nerve impulse equations[END_REF] with this set of parameters. Moreover, it seems that the speed of propagation of these waves decreases as ε grows, since in the case (b), the speed of propagation of these traveling pulses is slightly less than in the case (c). Then, we display in Table 3 and4, the numerical approximations of D ε (t) at fixed time t = 250 for several values of ε for the first order (left table) and the second order (right table) numerical schemes. Since the behavior of V ε is too different from its limit for smaller values of ε, we display linear regressions only from the line corresponding to ε = 1. These linear regressions yield that D ε (t) seems to be approximately of order two in ε for both numerical schemes, which corresponds to the one obtained by formal computations for the continuous problem [START_REF] Crevat | Diffusive limit of a spatially-extended kinetic FitzHugh-Nagumo model[END_REF].

Notice that the second order numerical scheme (2.21)-(2.25) represents a negligible improvement for the speed of convergence of D ε (t) as ε goes to 0. A key issue in numerical analysis is to perform a similar study on the discrete solution as the one we performed on the continuous problem [START_REF] Crevat | Diffusive limit of a spatially-extended kinetic FitzHugh-Nagumo model[END_REF] in order to establish the asymptotic preserving property of the scheme.

Heterogeneous neuron density

In the spirit of [START_REF] Bressloff | Traveling fronts and wave propagation failure in an inhomogeneous neural network[END_REF][START_REF] Bressloff | From invasion to extinction in heterogeneous neural fields[END_REF], the study of propagating waves in neural networks with spatial heterogeneities seems to be a fruitful topic. This subsection is therefore devoted to the illustration of the behaviour of the solution of the numerical scheme (2.17)-(2.18) with a non constant neuron density function ρ 0 . We choose the initial datum

f 0 (x, v, w) = ρ 0 (x) χ A v -V 0 (x) 10 χ A w -W 0 (x) 100 ,
with A = (-1/2, 1/2) where the density ρ 0 is a smooth approximation of 1 -χ B(0,6) and (V 0 , W 0 ) is chosen as

V 0 (x) = 1 if x 1 ∈ (-14, -13), 0 else, W 0 (x) = 0.1 if x 2 ≤ -14, 0 else. (3.7)
The domain in space is taken to be (-15, 15) 2 , discretized using n x = 512 points in each spatial coordinate and M = 50 particles per cell. It is expected that a wave will propagate initially from the left hand side in the homogeneous density of neurons. Then in the center of the domain, the density becomes inhomogeneous, which will perturb the wave propagation front. In Figure 3.2, we propose different scenario depending on the scaling parameter ε > 0. We display the profile of the solution V ε at time t = 300, 500 and 700 for ε = 5, 2 and 10 -2 . Clearly, the amplitude of the scaling parameter ε > 0 has an influence on the shape of the pulse but also on the speed of propagation.

First of all, the scrolling wave does not propagate through the ball B(0, 6), since the neuron density is too weak. Then, we can observe that as ε grows small, the speed of propagation and the width of the scroll wave increase. Thus, the heterogeneity does not have exactly the same effect. For ε = 5 and ε = 2 for example, the width of the gap in the neuron density is too large compared to the width of the traveling pulse. Therefore, the scroll wave breaks at its middle, and then recomposes once the heterogeneity is passed. Then, for smaller values of ε, as ε = 0.01, the traveling pulse starts to wrap the area where it cannot propagate before breaking and recomposing.

(a) t = 300 (b) t = 500 (c) t = 700

Rotating spiral waves

A spiral wave in the broadest sense is a rotating wave traveling outward from a center. Such spiral waves have been observed in many biological systems [START_REF] Winfree | The geometry of biological time[END_REF], [START_REF] Murray | Mathematical biology II: spatial models and biomedical applications[END_REF], such as mammalian cerebral cortex [START_REF] Huang | Spiral wave dynamics in neocortex[END_REF]. Although circular waves were predicted from early models of cortical activity [START_REF] Beurle | Properties of a mass of cells capable of regenerating pulses[END_REF], true spiral wave formation has been already obtained in numerical simulations of reaction-diffusion systems such as the Wilson-Cowan system [START_REF] Bueno-Orovio | Fourier spectral methods for fractional-in-space reactiondiffusion equations[END_REF][START_REF] Wilson | Excitatory and inhibitory interactions in localized populations of model neurons[END_REF].

In this section, we present numerical evidence for stable spiral waves considering the transport equation (1.3)-(1.4). We choose the initial datum [START_REF] Bueno-Orovio | Fourier spectral methods for fractional-in-space reactiondiffusion equations[END_REF] f

0 (x, v, w) = ρ 0 (x) χ A v -V 0 (x) 10 χ A w -W 0 (x) 100 ,
with A = (-1/2, 1/2) where the density ρ 0 is a smooth approximation of the characteristic function on the disk centered in 0 with radius 12, whereas (V 0 , W 0 ) is chosen as We first perform several computations changing the value of the scaling parameter and report in Figure 3.3, the profile of the numerical solution V ε obtained using the second order scheme (2.21)-(2.23) at the final time of the simulation t = 800. On the one hand, when ε ≥ 6, we observe that the initial wave first propagates into the domain, then it is damped and the solution converges to the stable steady state (V, W ) = (0, 0) times goes on (see Figure 3.3 (a) at time t = 800). On the other hand, when ε becomes smaller ε ∈ (4, 6), the solution evolves in a different manner. Indeed, the initial wave propagates into the physical domain where ρ 0 > 0, and a spiral wave appears at time t 20, where a traveling pulse emerges and propagates from the bottom left quarter of the domain, towards the bottom right quarter, which creates a rotating spiral wave at larger time. For these values of ε, the shape of the solution is very sensitive to ε (see for instance (b) and (c) in Figure 3.3). Finally, when ε ≤ 4, a spiral wave appears and it seems that the solution is not anymore sensitive to ε. In Figure 3.4, we report the numerical results for ε = 0.5 at different time t ∈ (0, 600). It illustrates how the spiral wave is generated from the initial data: a traveling pulse appears and begins to rotate clockwise, while the waves propagate up to the edge of the region where ρ 0 > 0. Moreover, it seems that once the spiral wave has appeared, its speed of rotation remains constant (see in (e) and (f) in Figure 3.4). Furthermore, in Figure 3.5, we report a zoom in the region where the traveling pulse appears. We observe that the center of the spiral moves and oscillates around a point. Finally in Figure 3.6, we propose the time evolution of the solution V ε at different points x = (-6, 3), x = (-8, 4) and x = (-8, 2). Close to the point x = (-6, 3), around which the spiral oscillates, time oscillations appear with an amplitude between -0.1 and 0.6 whereas in the neighboring points, different oscillations appear with a larger amplitude. Observe that at x = (-8, 4) and x = (-8, 2), the time oscillations look the same but are shifted.

V 0 (x) = 1 if x 1 ≤ -6 and x 2 ∈ (0, 3), 0 else, W 0 (x) = 0.1 if x 2 ≥ 3, 0 else. ( 3

Conclusion

In the present paper we have proposed a class of semi-implicit time discretization techniques for particle simulations to (1.3)-(1.4) coupled with a spectral collocation method for the space discretization. The main feature of our approach is to guarantee the accuracy and stability on slow scale variables even when the amplitude of local interactions becomes large, thus allowing a capture of the correct behavior with a large time step with respect to ε > 0. Even on large time simulations the obtained numerical schemes also provide an acceptable accuracy on the membrane potential when ε 1, whereas fast scales are automatically filtered when the time step is large compared to ε 2 .

As a theoretical validation we have proved that under some stability assumptions on numerical approximations, the slow part of the approximation converges when ε → 0 to the solution of a limiting scheme for the asymptotic evolution, that preserves the initial order of accuracy. Yet a full proof of uniform accuracy remains to be carried out in the frame of the continuous case [START_REF] Crevat | Diffusive limit of a spatially-extended kinetic FitzHugh-Nagumo model[END_REF]. The main challenge is to rigorously study the stability of the numerical solution for an appropriate norm and to find a bound uniformly with respect to ε and h.

Figure 2 . 1 :

 21 Figure 2.1: Asymptotic preserving diagram is performed to evaluate uniform error estimates P ε -P ε h with respect to ε.

= 1 Figure 3 . 1 :

 131 Figure 3.1: Order of accuracy in ε → 0: spatio-temporal profile of V ε (t, x) for (a) ε = 3.25, (b) ε = 3 and (c) ε = 1.

Figure 3 . 2 :

 32 Figure 3.2: Heterogeneous neuron density : plot of the solution V ε at different time t = 300, 500 and 700 for ε = 5 (top), ε = 2 (middle) and ε = 10 -2 (bottom).

. 8 ) 2 Figure 3 . 3 :

 8233 Figure 3.3: Rotating spiral waves : plot of the solution V ε at time t = 800 for different values of ε > 0.

Figure 3 . 4 :

 34 Figure 3.4: Rotating spiral waves : plot of the solution V ε for ε = 0.5 at different time t ∈ [0, 600].

Figure 3 . 5 :

 35 Figure 3.5: Rotating spiral waves : zoom on the solution V ε for ε = 0.5 at different time t ∈ [0, 400] around the point where the traveling pulse emerges.

Figure 3 . 6 :

 36 Figure 3.6: Rotating spiral waves : time evolution of the solution V ε for ε = 0.5 at different points around the where the traveling pulse emerges.

Table 1 :

 1 Order of accuracy in h → 0: evaluation of the numerical error at fixed time T = 10 of the numerical schemes (2.17)-(2.18).

		for (2.17)-(2.18)	Order
	1.e-01	5.48e-04	XXX
	5.e-02	2.73e-04	1.63
	2.e-02	1.09e-04	1.00
	1.e-02	5.47e-05	1.00
	5.e-03	2.73e-05	1.00
	2.e-03	1.09e-05	1.00
	1.e-03	5.47e-06	1.00
	5.e-04	2.73e-06	1.00
	h	L 2 error for (2.21)-(2.25)	Order
	1.e-01	1.23e-07	XXX
	5.e-02	3.56e-08	2.69
	2.e-02	8.35e-09	2.02
	1.e-02	2.07e-08	2.01
	5.e-03	5.01e-09	2.01
	2.e-03	1.23e-09	2.01
	1.e-03	2.95e-10	2.01
	5.e-04	2.95e-10	2.00

Table 2 :

 2 Order of accuracy in h → 0: evaluation of the numerical error at fixed time T = 10 of the numerical schemes (2.21)-(2.25).

  In Table 1 and 2, we report the numerical error for different values of h at fixed time T = 10 for the numerical schemes (2.17)-(2.18) (first table) and (2.21)-(2.25) (second table

Table 3 :

 3 Order of accuracy in ε → 0: approximation of D ε (t) at fixed time t = 250 with the first order scheme (2.17)-(2.18).

	ε	D ε (t) with (2.17)-(2.18)	Order
	5	1.21
	2	1.73	XX
	1	9.16e-01	XX
	5.e-01	2.60e-01	1.82
	2.e-01	4.17e-02	1.92
	1.e-01	1.04e-02	1.95
	5.e-02	2.60e-03	1.97
	2.e-02	4.17e-04	1.98
	1.e-02	1.04e-04	1.98
	5.e-03	2.62e-05	1.99
	2.e-03	4.24e-06	1.99
	1.e-03	8.65e-07	2.00
	ε	D ε (t) with (2.21)-(2.25)	Order
	5	1.21
	2	1.73	XX
	1	9.13e-01	XX
	5.e-01	2.59e-01	1.83
	2.e-02	4.15e-02	1.93
	1.e-01	1.04e-02	1.95
	5.e-02	2.59e-03	1.97
	2.e-02	4.15e-04	1.98
	1.e-02	1.03e-04	1.98
	5.e-03	2.59e-04	1.99
	2.e-03	6.94e-05	1.99
	1.e-03	1.74e-05	1.99

Table 4 :

 4 Order of accuracy in ε → 0: approximation of D ε (t) at fixed time t = 250 with the second order scheme (2.21)-(2.25).

Acknowledgements

The authors acknowledge support from ANITI (Artificial and Natural Intelligence Toulouse Institute) Research Chair and the project ChaMaNe (ANR-19-CE40-0024).

order singly diagonally implicit Runge-Kutta (SDIRK) method for the implicit part. According to the classification from [START_REF] Boscarino | High order semi-implicit schemes for time dependent partial differential equations[END_REF], we call it H-SDIRK2 (2,2,2). For all n ∈ N and p ∈ {1 . . . M }, we apply a first stage,

Hence we compute the additional variable V

M ∈ S nx solution of the following scheme

Then, we set

and compute the second stage with a semi-implicit step on (V

p , W

Moreover, V

M ∈ S nx is given by

M -W n M and W

(1)

W (1) p .

Finally, we get the numerical solution at time t n+1 through

(2.25)

Now, we prove an analogous result to Proposition 2.3 for the numerical scheme (2.21)-(2.25) as the parameter ε goes to 0 with a fixed numerical parameter h ∈ R 3 given by h = (∆t, ∆x, 1/M ), where ∆x = 2π/n x , ∆t > 0 and M ∈ N * .

Proposition 2.5 (Consistency when ε → 0). Let h to be fixed and consider a connectivity kernel Ψ :

2), (2.7) and a neuron density ρ 0 ∈ S nx satisfying (1.9) at each grid point x j , j ∈ J nx . For all ε > 0, p ∈ {1, . . . , M } and n ∈ N, let us assume that the triplet (V ε,n p , W ε,n p , V ε,n M ) given by (2.21)-(2.25) is uniformly bounded with respect to ε > 0. Then we define

where the second stage is given by

M ,

where V (1) M = 2V

(1)

M = 2W

(1)

The next time step is given by

M -W n M .

(2.28)

Proof. We fix a time step ∆t > 0, a set of equidistant points (x j ) j∈Jn x ⊂ T and p ∈ {1, . . . , M }. Then we denote by V ε,n p , W ε,n p , V ε,n M ε>0 the solution of (2.21)- (2.25). Up to a sub-sequence,

n M when ε → 0, hence we proceed exactly as in Proposition 2.3 and set

M .

Thus, we prove that V

M , W

(1) M corresponds to the solution of the first stage (2.26) and we have

Furthermore, we treat the second stage in the same manner for any j ∈ J nx and p ∈ {1, . . . , M }, the limit V

p (x j ) does not depend on p and is given by

Passing to the limit as ε → 0 in (2.24) and in the second equation in (2.23), it yields that (V

M , W

M ) satisfies (2.27) and finally (2.28).

Let us notice that the present strategy can be applied to a large class of second order schemes and can also be extended to a third order semi-implicit scheme. We refer to [START_REF] Boscarino | High order semi-implicit schemes for time dependent partial differential equations[END_REF] for the detailed description of the schemes. Now, let us investigate the consistency with respect to the numerical parameter h. Lemma 2.6 (Consistency h → 0). Let 0 < ε < 1 be a fixed parameter and consider Ψ satisfying (1.2), and ρ 0 satisfying (1.9) and such that Supp(ρ 0 ) ⊂ B(0, S) where S > 0. We suppose that f ε the solution of (1.3)-(1.4) is differentiable twice with respect to time and there exists a constant C T > 0, independent of ε, such that in [0, T ] × B(0, S),

Consider the scheme (2.21)-(2.25) and the consistency error E n h in (2.20). Then, there exists another positive constant C > 0, independent of h and ε, such that for all n ∈ {0, . . . , [T /∆t]},

Proof. The proof uses the same tools as in the proof of Lemma 2.4.

As in the previous section, Proposition 2.5 indicates that in the limit ε → 0, the numerical scheme (2.21)-(2.25) becomes a second order explicit time approximation with respect to ∆t of the reaction-diffusion system (1.12). Applying Lemma 2.6, we may conjecture that the asymptotic preserving property (2.16) is satisfied.

In the following section, we provide some numerical evidences on this issue.