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

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.

## Authors

• 1 publication
• 7 publications
03/18/2020

### Convergence Analysis of Asymptotic Preserving Schemes for Strongly Mangnetized Plasmas

The present paper is devoted to the convergence analysis of a class of a...
03/18/2020

### Convergence Analysis of Asymptotic Preserving Schemes for Strongly Magnetized Plasmas

The present paper is devoted to the convergence analysis of a class of a...
10/29/2021

### An Assessment of Solvers for Algebraically Stabilized Discretizations of Convection-Diffusion-Reaction Equations

We consider flux-corrected finite element discretizations of 3D convecti...
01/20/2022

### Numerical simulation of singularity propagation modeled by linear convection equations with spatially heterogeneous nonlocal interactions

We study the propagation of singularities in solutions of linear convect...
06/12/2020

### Asymptotic preserving IMEX-DG-S schemes for linear kinetic transport equations based on Schur complement

We consider a linear kinetic transport equation under a diffusive scalin...
03/29/2021

### Translating Numerical Concepts for PDEs into Neural Architectures

We investigate what can be learned from translating numerical algorithms...
01/02/2020

### Asymptotically compatible reproducing kernel collocation and meshfree integration for the peridynamic Navier equation

In this work, we study the reproducing kernel (RK) collocation method fo...
##### This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.

## 1 Introduction

The FitzHugh–Nagumo (FHN) system [17, 27]

, 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

neurons interacting with each other, where each one is labeled by , and endowed with a parameter for standing for the constant spatial position in the network. The FHN system accounts for the variations of the membrane potential of a neuron coupled to an auxiliary variable called the adaptation variable. It can be written as follows for all ,

 (1.1)

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

 Ψε(∥y∥):=1εdΨ(∥y∥ε),y∈Rd,

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,

 (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 [8], 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,

 (1.3)

with and given by

 ⎧⎪ ⎪ ⎪⎨⎪ ⎪ ⎪⎩Kε[fε](t,x,v)=1ε2∫Rd∫R2Ψε(∥x−x′∥)(v′−v)fε(t,x′,dv′,dw′)dx′,A(v,w)=τ(v−γw). (1.4)

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 in

and , 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

 ⎧⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪⎨⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪⎩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 the flow , hence the solution to (1.3)–(1.4) is given by

 fε(t,x,.)=Φt,x#fε0(x,.), (1.6)

that is, for any test-function and ,

 ∫Bφ(v,w)fε(t,x,dv,dw)=∫Φ−1t,x(B)φ∘Φt,xfε0(x,dv,dw).

We also define for all and the following macroscopic quantities,

 ρε⎛⎜⎝1VεWε⎞⎟⎠(t,x):=∫R2⎛⎜⎝1vw⎞⎟⎠fε(t,x,dv,dw), (1.7)

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

 Kε[fε](.,v)=1ε2[Ψε⋆(ρεVε)−Ψε⋆ρεv], (1.8)

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

 ρε(t,x)=ρε0(x),x∈Rd

and moreover we suppose that it does not depend neither on , so that with

 ρ0≥0,ρ0∈L∞(Rd). (1.9)

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

 ∫R2Kε[fε](t,x,v)fε(t,x,dvdw)=ρ0ε2[Ψε⋆(ρ0Vε)−(Ψε⋆ρ0)Vε]. (1.10)

Hence, multiplying (1.3) by (resp. ) and integrating with respect to and using (1.10), we get a time evolution equation for the macroscopic quantities as

 ⎧⎪ ⎪ ⎪⎨⎪ ⎪ ⎪⎩∂t(ρ0Vε)−ρ0ε2[Ψε⋆(ρ0Vε)−(Ψε⋆ρ0)Vε]=∫R2N(v)fε(.,dv,dw)−ρ0Wε,∂t(ρ0Wε)=ρ0A(Vε,Wε). (1.11)

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 [9], 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,

 ⎧⎪⎨⎪⎩ρ0(∂tV−¯¯¯σ[Δ(ρ0V)−VΔρ0]−N(V)+W)=0,ρ0(∂tW−A(V,W))=0, (1.12)

where is defined in (1.2). We refer to [9]

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 [8], 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 [5], 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 [9], 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

 K[f](t,x,v) = 1ε2∫RdΨε(∥y∥)ρ0(x−y)(V(t,x−y)−v)dy.

Then we define a truncated operator in the following way.

###### Lemma 2.1.

Suppose that , where is the ball of radius centered at the origin and choose . Then, for any , is solution to

 ∂tf+∂v[f(N(v)−w+KS[f])]+∂w[fA(v,w)]=0,

where for any ,

 KS[f](t,x,v)=χB(0,S)ε2(x)∫B(0,2S)Ψε(∥y∥)ρ0(x−y)(V(t,x−y)−v)dx′, (2.1)

where

denotes the characteristic function in the ball

.

###### Proof.

On the one hand, since and for all , the density , we get that for any , the transport equation (1.3) can be written as

 ∂tf+∂v[f(N(v)−w+χB(0,S)K[f])]+∂w[fA(v,w)]=0.

Then it is enough to consider only . On the other hand, the domain of integration of the operator is such that

 ∥y∥≤∥x∥+∥y−x∥≤2S,

hence for any ,

 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 . ∎

Actually the operator can be seen as convolution products between and the connectivity kernel , that is,

 KS[f](t,x,v)=1ε2(LS[ρ0V](t,x)−vLS[ρ0](x)),

where is given by

 LS[u]=Ψε⋆u,u∈{ρ0,ρ0V}. (2.2)

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,

 ˜u(t,k):=1ndx∑j∈Jnxu(t,xj)e−ik⋅xj,k∈Jnx.

and get a trigonometric polynomial

 unx(t,x):=∑k∈Jnx˜u(t,k)eik⋅x,u∈{ρ0,ρ0V}.

Therefore, we substitute this polynomials in (2.2), which yields a discrete operator given by

 LSnx[u]:=∑k∈Jnx˜LS[u](t,k)eik⋅x, (2.3)

where is given by

 ˜LS[u](t,k)=(2π)dˆΨε(k)˜u(t,k) (2.4)

and is the expansion coefficient depending on the connectivity kernel

 (2.5)

Finally the approximation of the operator is provided by

 KSnx[f](t,x,v)=1ε2(LSnx[ρ0V](t,x)−vLSnx[ρ0](x)). (2.6)

Let us focus on the computation of the kernel modes for any fixed parameter . In the spirit of [28] 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:

 ˆΨε(k)=1(2π)d∫π0Ψε(r)rd−1I(k,r)dr,

where

Then, changing the variable into , we get

 ˆΨε(k)=1(2π)d∫π/ε0Ψ(s)sd−1I(k,εs)ds.

To complete the computation of the function , we have to study separately each possible value for the spatial dimension .

##### One-dimensional case: d=1.

Since , it is straightforward to check that for any ,

 I(k,r)=2cos(r|k|),

hence we get:

 ˆΨε(k)=1π∫π/ε0Ψ(s)cos(εs|k|)ds.
##### Two-dimensional case: d=2.

Let and . In this case, if we note , then using spherical coordinates, we have

 =2∫π0cos(r∥k∥sinθ)dθ =2πJ0(r∥k∥),

where is the Bessel function of order , defined with

 J0:x∈R↦1π∫π0cos(xsinθ)dθ=∞∑l=0(−1)l(l!)2(x2)2l.

Consequently, we get

 ˆΨε(k)=12π∫π/ε0Ψ(s)sJ0(εs∥k∥)ds.
##### Three-dimensional case: d=3.

Let and . Hence, if we note , and then using spherical coordinates, we get

 =2π∫1−1exp(i∥q∥μ)dμ =4πSinc(r∥k∥),

where . Thus, the kernel mode is given by

 ˆΨε(k)=12π2∫π/ε0Ψ(s)|s|2Sinc% (εs∥k∥)ds.

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 [6]

 Snx=span{eik⋅x,−nx/2≤ks≤nx/2−1,s=1,…,d},

equipped with the classical norm , which satisfies [6] for any

 ∥u∥2L2=(2πnx)d∑j∈Jnx|u(xj)|2

and for any and , we also have

 ∫Tu(x)¯¯¯v(x)dx=(2πnx)d∑j∈Jnxu(xj)¯¯¯v(xj).

Finally we define by the projection operator from to such that , for all .

###### Proposition 2.2.

Let and consider a connectivity kernel satisfying (1.2) with

 ∫RdΨ(∥y∥)∥y∥4dy<∞. (2.7)

Then, for all , there exists a positive constant , depending on , such that for all ,

 ∣∣(2π)dˆΨε(k)−¯¯¯¯Ψ+¯¯¯σε2∥k∥2∣∣≤C(∥k∥4+1)ε4. (2.8)

Moreover for any trigonometric polynomial , we have

 ∥∥LSnx[u]−¯¯¯¯Ψu−¯¯¯σε2Δu∥∥L2≤Cε4(∥Δ2u∥L2+∥u∥L2). (2.9)
###### Proof.

On the one hand, for any , we perform a Taylor expansion of at and using the assumptions (2.7) on , it yields

 ∣∣∣(2π)dˆΨε(k)−∫π/ε0Ψ(s)sd−1ds−ε2∥k∥2∫π/ε0Ψ(s)sd+1ds∣∣∣≤∥k∥4ε4∫Rd∥y∥4Ψ(∥y∥)dy.

On the other hand, we have

 ∫∞π/εΨ(s)sd−1ds+ε2∫∞π/εΨ(s)sd+1ds≤ε4(1π4+1π2)∫Rd∥y∥4Ψ(∥y∥)dy.

Gathering these results and using (1.2), there exists a constant , depending on , such that

 ∣∣(2π)dˆΨε(k)−¯¯¯¯Ψ+¯¯¯σε2∥k∥2∣∣≤C(∥k∥4+1)ε4.

Then, we consider and for , we substitute the latter result in the expression (2.4) of , it yields for each

 ∣∣˜LS[u](k)−(¯¯¯¯Ψ+¯¯¯σε2∥k∥2)˜u(k)∣∣≤Cε2(∥k∥4+1)|˜u(k)|.

Thus, from the definition of (2.3), we know that and get

 ∥LSnx[u]−¯¯¯¯Ψu−¯¯¯σε2Δu∥L2 = ⎛⎝∑k∈Jnx∣∣˜LS[u](k)−(¯¯¯¯Ψ+¯¯¯σε2∥k∥2)˜u(k)∣∣2⎞⎠1/2, ≤ Cε4(∥Δ2u∥L2+∥u∥L2).

### 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 [18] for the numerical computation of specific problems in fluid dynamics, and precisely mathematically studied later [30]. 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 ,

 fM(t,x,dv,dw):=ρ0(x)MM∑p=1δVp(t,x)(dv)⊗δWp(t,x)(dw),

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

 ⎧⎪ ⎪ ⎪ ⎪ ⎪ ⎪⎨⎪ ⎪ ⎪ ⎪ ⎪ ⎪⎩dVpdt=Inx(N(Vp)+KSnx[fM](Vp))−Wp,dWpdt=A(Vp,Wp), (2.10)

with a given initial data for and is the projection operator on . Moreover, we define the macroscopic potential at each point , as

 ⎧⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪⎨⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪⎩ρ0=∫R2fM(t,x,dv,dw),ρ0VM(t,x):=∫R2vfM(t,x,dv,dw)=1MM∑p=1ρ0(x)Vp(t,x),ρ0WM(t,x):=∫R2wfM(t,x,dv,dw)=1MM∑p=1ρ0(x)Wp(t,x). (2.11)

From these macroscopic quantities, it is then possible to compute the discrete operator given in (2.6), where (2.10)–(2.11) are solved at each mesh point .

### 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

 ⎧⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪⎨⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪⎩Vn+1p−VnpΔt=Inx(N(Vnp)+1ε2[LSnx[ρ0VnM]−Vn+1pLSnx[ρ0]])−Wnp,Wn+1p−WnpΔt=A(Vn+1p,Wnp), (2.12)

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

 Inx(N(Vnp)+1ε2[LSnx[ρ0VnM]−Vn+1pLSnx[ρ0]])−Wnp=Inx(N(Vnp))+1ε2[LSnx[ρ0VnM]−Inx(Vn+1pLSnx[ρ0])]−Wnp.

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

 Vn+1M−VnMΔt=1Mp=1M∑Inx(N(Vn+1p))+1ε2[LSnx[ρ0VnM]−Inx(VnMLSnx[ρ0])]−WnM. (2.13)

Observe here that the nonlinear term is computed implicitly from whereas the stiff term is now explicit.

Now, we define a numerical parameter as , where and let us show the consistency of the numerical scheme (2.12)–(2.13) in the limit for a fixed numerical parameter .

###### Proposition 2.3 (Consistency in the limit ε→0 for fixed numerical parameters h).

Let be a fixed parameter and consider a connectivity kernel satisfying (1.2), (2.7) and a neuron density satisfying (1.9) at each grid point , . For all , and , let us assume that the triplet given by (2.12)–(2.13) is uniformly bounded with respect to . Then we define

 Wε,nM=1MM∑p=1Wε,np