A space-time finite element method for neural field equations with transmission delays

TL;DR

This paper introduces a novel space-time finite element method using discontinuous Galerkin in time for neural field equations with transmission delays, enabling more flexible and accurate simulations of neural dynamics.

Contribution

It develops and analyzes a new dGcG-FEM algorithm for neural delay equations, addressing limitations of existing methods and including an a-priori error analysis.

Findings

Effective handling of space-dependent delays in neural models

Successful application to models with inhomogeneous kernels

Provides a detailed algorithm and error analysis

Abstract

We present and analyze a new space-time finite element method for the solution of neural field equations with transmission delays. The numerical treatment of these systems is rare in the literature and currently has several restrictions on the spatial domain and the functions involved, such as connectivity and delay functions. The use of a space-time discretization, with basis functions that are discontinuous in time and continuous in space (dGcG-FEM), is a natural way to deal with space-dependent delays, which is important for many neural field applications. In this article we provide a detailed description of a space-time dGcG-FEM algorithm for neural delay equations, including an a-priori error analysis. We demonstrate the application of the dGcG-FEM algorithm on several neural field models, including problems with an inhomogeneous kernel.

Equations

$$\frac{\partial V_i}{\partial t}(t,r)=-{\alpha }_i{V}_i(t,r)+\sum _{j=1}^p{\int }_{\Omega }{J}_{ij}(r,r^{\prime },t)S_j(V_j(t-{\tau }_{ij}(r,r^{\prime }),r^{\prime }))d{r}^{\prime },$$

$$\frac{\partial u}{\partial t}(t,x)=-\alpha u(t,x)+{\int }_{\Omega }J(x,r)S(u(t-\tau (x,r),r))dr,\alpha >0.$$

$$0<{\tau }_{max}=\underset{(x,r)\in \bar{\Omega }\times \bar{\Omega }}{\sup }\tau (x,r)<\infty .$$

$$\Vert \phi {\Vert }_X=\underset{s\in [-{\tau }_{max},0]}{\sup }\Vert \phi (s,\cdot ){\Vert }_Y,$$

$$\Vert J{\Vert }_C=\underset{(x,r)\in \bar{\Omega }\times \bar{\Omega }}{\sup }|J(x,r)|.$$

$$u_t(s)=u(t+s)\in C(\bar{\Omega }),s\in [-{\tau }_{max},0],t\ge 0.$$

$$G(\phi )(x)={\int }_{\Omega }J(x,r)S\left(\phi (-\tau (x,r),r)\right)dr.$$

$$\frac{\partial u}{\partial t}(t)=-\alpha u(t)+G(u_t),$$

$$\displaystyle\Bigl{(}\frac{\partial u}{\partial t}(t)+\alpha u(t),v\Bigr{)}-\left(G(u_{t}),v\right)=0,$$

$$u(s)=u_0(s),$$

$$\bigl{(}G(u_{t}),v\bigr{)}=\int_{\Omega}G(u_{t})(x)v(x)dx=\int_{\Omega}\int_{\Omega}J(x,r)S\left(u_{t}(-\tau(x,r),r)\right)dr\,v(x)dx.$$

$${\bar{\mathcal{T}}}_h=\left\{K_j:\bigcup _{j=1}^M{\bar{K}}_j={\bar{\Omega }}_h,K_j\cap K_i=\varnothing \text{if}i\ne j\right\}.$$

$${\mathcal{T}}_h^n=\left\{\mathcal{K}=G_{\mathcal{K}}^n(\hat{\mathcal{K}}):K\in {\bar{\mathcal{T}}}_h\right\},$$

$$V_h^n$$

$$V_h=\{u\in L^2(\mathcal{E}):u{\mid }_{{\mathcal{E}}_n}\in V_h^n,n=1,2,\dots ,N\}.$$

$$\displaystyle\sum_{n=1}^{N}\sum_{\mathcal{K}\in\mathcal{T}_{h}^{n}}\left[\Bigl{(}\frac{\partial u_{h}}{\partial t}+\alpha u_{h},v_{h}\Bigr{)}_{\mathcal{K}}-\int_{\mathcal{K}}\left[\int_{\Omega}J(x,r)S\left(u_{h}(t-\tau(x,r),r)\right)d\,r\right]v_{h}(t,x)dx\,dt\right]$$

$$+\sum _{n=2}^N\left([u_h{]}_{n-1},v_h^{n-1,+}\right)+\left(u_h^{0,+},v_h^{0,+}\right)=\left(u_0(0),v_h^{0,+}\right)$$

$$\hat{G}\psi ={\int }_{\Omega }J(\cdot ,r)S\left(\psi (-\tau (\cdot ,r),r)\right)dr,\psi \in \hat{X}.$$

$${\int }_{{\mathcal{K}}^n}\left(-u_h\frac{\partial v_h}{\partial t}+\alpha u_h{v}_h\right)dxdt+{\int }_{K(t_n)}{u}_h^{n,-}{v}_h^{n,-}dx$$

$$-{\int }_{{\mathcal{K}}^n}\hat{G}\left({u_h}_t\right)(x)v_{h}dxdt={\int }_{K(t_{n-1})}{u}_h^{n-1,-}{v}_h^{n-1,+}dx,$$

$$u_h(t,x){\mid }_{\mathcal{K}}=\sum _{m=1}^{N_p}{\hat{u}}_m^{\mathcal{K}}{\psi }_m^{\mathcal{K}}(t,x)$$

$${\int }_{\mathcal{K}}{\psi }_i^{\mathcal{K}}(t,x)\left({\int }_{\Omega }J(x,r)S\left(u_h(t-\tau (x,r),r)\right)dr\right)dxdt$$

$$\displaystyle={\int_{\mathcal{K}}}\psi_{i}^{\mathcal{K}}(t,x)\Bigl{(}\sum_{L\in\bar{\mathcal{T}}_{h}}\int_{L}J(x,r)S\Bigl{(}\sum_{m=1}^{N_{p}}\hat{u}_{m}^{L}\psi_{m}^{L}(t-\tau(x,r),r)\Bigr{)}dr\Bigr{)}dxdt.$$

$$\displaystyle\sum_{j=1}^{N_{p}}\left\{\hat{u}_{j}^{\mathcal{K}}\int_{\mathcal{K}}\Bigl{(}-\psi_{j}^{\mathcal{K}}(t,x)\frac{\partial}{\partial t}\psi_{i}^{\mathcal{K}}(t,x)+\alpha\psi_{j}^{\mathcal{K}}(t,x)\psi_{i}^{\mathcal{K}}(t,x)\Bigr{)}dxdt\right.$$

$$+{\hat{u}}_j^{\mathcal{K}}{\int }_{K(t_n^{-})}{\psi }_j^{\mathcal{K}}(t_n^{-},x){\psi }_i^{\mathcal{K}}(t_n^{-},x)dx\}$$

$$\displaystyle-\sum_{L\in\bar{\mathcal{T}}_{h}}{\int_{\mathcal{K}}}\psi_{i}^{\mathcal{K}}(t,x)\Bigl{[}\int_{L}J(x,r)S\Bigl{(}\sum_{m=1}^{N_{p}}\hat{u}_{m}^{L}\psi_{m}^{L}(t-\tau(x,r),r)\Bigr{)}dr\Bigr{]}dxdt$$

$$=\sum _{j=1}^{N_p}{\hat{u}}_j^{\mathcal{K},n-1}{\int }_{K(t_{n-1}^{+})}{\psi }_j^{\mathcal{K}}(t_{n-1}^{-},x){\psi }_i^{\mathcal{K}}(t_{n-1}^{+},x)dx,$$

$$S_k=\{w:[0,T]\to Y:w{\mid }_{I_n}=\sum _{j=0}^q{\phi }_j{t}^j,{\phi }_j\in Y,\forall n\ge 1\},$$

$$T_k:C\left([0,T],Y\right)\to S_k$$

$$\tilde{u}(t_{n-1})=u(t_{n-1}),$$

$${\int }_{I_t}$$

$$\Vert \tilde{u}(s)-u(s)\Vert \le C_I{k}_t^{q+1/2}{\left({\int }_{I_t}\Vert {\partial }_s^{q+1}u(s,\cdot ){\Vert }^{2}ds\right)}^{1/2},\text{ for }s\in I_t,$$

$$W_h=\{v\in C(\Omega ):v{\mid }_K\circ G_K^n\in {\hat{\mathcal{P}}}_r(\hat{K}),\forall K\in {\bar{\mathcal{T}}}_h\},$$

$$P_h:Y\to W_h$$

Full text

Key words. Neural fields, transmission delays, discontinuous Galerkin, finite element methods, space-time methods

AMS subject classifications. 65M60, 65M15, 65R20, 37M05, 92C20

1 Introduction

The motivation of this work is the need for numerical methods that can accurately and efficiently discretize delayed integro-differential equations originating from neural field models, in particular when the delay in the system is space dependent. Only a few studies considered so far the numerical treatment of neural field systems, see [8], [10], [11] and references therein. In [8], the authors used special types of delay and connectivity functions in order to reduce the spatial discretization to a large system of delay differential equations with constant time delays. This system was then solved with the Matlab solver dde23. In [10] a new numerical scheme was introduced that includes a convolution structure and hence allows the implementation of fast numerical algorithms. In both studies the connectivity kernel depends on the distance between two spatial locations. This choice has been shown to model successfully neural activity known from experiments, it introduces, however, also a limitation to the applicability of the presented techniques.

Here we propose the use of space-time finite element methods using discontinuous basis functions in time and continuous basis functions in space (dGcG-FEM), which are well established to solve ordinary and partial differential equations, e.g. [5], [6], [7], [9], [12], [13]. The novelty of this work is the successful application of the space-time dGcG-method to the neural field equations. The motivation of our choice is that the time-discontinuous Galerkin method has good long-time accuracy, [6], [12]. Moreover, the use of a space-time discretization is a natural way to deal also with the space-dependent delays. As it will be discussed later, there is no need in a space-time method to interpolate the solution from previous time levels. The space-time dGcG-method was successfully used for stiff systems and is well suited for mesh adaptation, which is of great importance when local changes in the solution are of interest. Further benefits are that we do not need to make restrictions, neither to the functions involved in the system, such as the connectivity kernel or the delay function, nor to the dimension or shape of the spatial domain.

In this article we present a novel space-time dGcG-method for delay differential equations. We provide a theoretical analysis of the stability and order of accuracy of the numerical discretization and demonstrate its application on a number of neural field problems. We focus on the design and an a-priori error analysis of the space-time dGcG-FEM for nonlinear neural field equations with space dependent delay.

The outline of this article is as follows. In the introductory Section 2 we recall a mathematical model for neural fields. In Section 3 we introduce the space-time dGcG-FEM method. The main difficulty is the treatment of the delay term in the neural field equations, which is investigated in detail in Section 3.2. An a-priori error analysis of the space-time discretization is given in Section 4. Next, we show in Section 5 some numerical simulations for the neural field equations in one spatial dimension with one population. These examples are taken from literature, [8], [14], where both analytical and numerical results are known for comparison. We demonstrate some further computational benefits of the space-time dGcG-FEM by introducing an inhomogeneous kernel in the delay term in Section 5.4. The numerical algorithms presented in [8] and [10], are not suitable for the treatment of local inhomogeneities.

In consecutive papers we will show computations on more complicated spatial domains and extend the model to more populations in the neural field system.

2 Neural fields with space dependent delays

The mathematical model for neural fields with space-dependent delays is as follows. Consider $p$ populations consisting of neurons distributed over a bounded, connected and open domain $\Omega\subset\mathbb{R}^{d},$ $d=1,2,3.$ For each $i,$ the variable $V_{i}(t,r)$ is the membrane potential at time $t,$ averaged over those neurons in the $i$th population positioned at $r\in\Omega.$ These potentials are assumed to evolve according to the following system of integro-differential equations

[TABLE]

for $i=1,\dots,p.$ The intrinsic dynamics exhibits exponential decay to the baseline level $0,$ as $\alpha_{i}>0.$ The propagation delays $\tau_{ij}(r,r^{\prime})$ measure the time it takes for a signal sent by a type-$j$ neuron located at position $r^{\prime}$ to reach a type-$i$ neuron located at position $r.$ The function $J_{ij}(r,r^{\prime},t)$ represents the connection strength between population $j$ at location $r^{\prime}$ and population $i$ at location $r$ at time $t.$ The firing rate functions are $S_{j}.$ For the definition and interpretation of these functions we refer to [15]. Some examples will be given in later sections.

Throughout this paper we consider a single population, $p=1,$ in a bounded domain $\Omega\subset\mathbb{R}^{d},$ on a time interval $[t_{0},T),$ with $T>t_{0}$ the final time,

[TABLE]

Note that we will only deal with autonomous systems. Therefore we assume from here on that the connectivity does not depend on time. We assume that the following hypotheses are satisfied for the functions involved in the system, (as in [14]): the connectivity kernel $J\in C(\bar{\Omega}\times\bar{\Omega}),$ the firing rate function $S\in C^{\infty}(\mathbb{R})$ and its $k$th derivative is bounded for every $k\in\mathbb{N}_{0},$ the delay function $\tau\in C(\bar{\Omega}\times\bar{\Omega})$ is non-negative.

Without loss of generality, we take $t_{0}=0.$ From the assumption on the delay function $\tau,$ we may set

[TABLE]

Note that when $\tau_{max}=0,$ the delay function $\tau(x,r)=0$ for all $(x,r)\in\bar{\Omega}\times\bar{\Omega},$ and in this case (2) reduces to an integro-differential equation without delay. As we will see later, our numerical method can handle this case as well.

Let $Y=C(\bar{\Omega})$ and set $X=C\left([-\tau_{max},0];Y\right).$ For $\varphi\in X,$ $s\in[-\tau_{max},0]$ and for $x\in\Omega$ we write $\varphi(s)(x)=\varphi(s,x),$ and its norm is given by

[TABLE]

where $\|\varphi(s,\cdot)\|_{Y}=\sup_{x\in\Omega}|\varphi(s,x)|.$ From the assumption on the connectivity kernel, it follows that it is bounded in the following norm

[TABLE]

We use the traditional notation for the state of the system at time $t$

[TABLE]

Define the nonlinear operator $G:X\to Y$ by

[TABLE]

Then the neural field equation (2) can be written as a delay differential equation (DDE) as

[TABLE]

where the solution is an element of $C([-\tau_{max},\infty);Y)\cap C^{1}([0,\infty);Y).$ Similarly, we have the state of the solution at time $t$ defined as $u_{t}(s)(x)=u(t+s,x),$ $s\in[-\tau_{max},0],$ $t\geq 0,$ $x\in\Omega.$ It was shown in [14] that under the above assumptions on the connectivity, the firing rate function and delay, the operator $G$ is well-defined and it satisfies a global Lipschitz condition.

Note that the assumptions on the firing rate function $S$ imposed in [14] were needed for further analysis of the neural field equations. For the numerical analysis presented in this paper it is sufficient to assume that $S$ is Lipschitz continuous.

3 The discontinuous Galerkin finite element

method

The starting point of our numerical discretization is the weak formulation. The numerical method is investigated for the nonlinear equation (4), which may be written in variational form as: Find $u\in C^{1}\left([0,T),Y\right)\cap C\left([-\tau_{max},T),Y\right)$ such that

[TABLE]

where $(\cdot,\cdot)$ is the usual $L^{2}(\Omega)$ inner product. Here the delay contribution is expressed as

[TABLE]

Note that for any $t>0,$ all functions in the inner product are elements of $Y=C(\bar{\Omega}),$ which is a dense subset of $L^{2}(\Omega),$ hence the inner product is well-defined.

3.1 The space-time dGcG-FEM discretization

Consider the neural field equations in the domain $\Omega.$ We will not distinguish between space and time variables and consider directly the space $\mathbb{R}^{d+1},$ where $d$ is the number of space dimensions.

Let $\mathcal{E}\subset\mathbb{R}^{d+1}$ be an open, bounded space-time domain in which a point has coordinates $(t,x)\in\mathbb{R}^{d+1},$ with ${x}\in\mathbb{R}^{d}$ the position vector and time $t.$ First, partition the time interval $\bar{I}=[0,T]$ using the time levels $0=t_{0}<t_{1}<\ldots<t_{N}=T$ and denote by $I_{n}=(t_{n-1},t_{n}]$ the $n$-th time interval of length $k_{n}=t_{n}-t_{n-1}.$ A space-time slab is defined as $\mathcal{E}^{n}=I_{n}\times\Omega.$ Second, we approximate the spatial domain $\Omega$ with $\Omega_{h}$ using a tessellation of non-overlapping hexahedral elements (line elements in 1D, quadrilaterals in 2D, etc.)

[TABLE]

The domain approximation is such that $\Omega_{h}\to\Omega$ as $h\to 0,$ where $h$ is the radius of the smallest sphere containing each element $K_{j}\in\bar{\mathcal{T}}_{h}.$ The space-time elements $\mathcal{K}_{j}^{n}$ are now obtained as $\mathcal{K}_{j}^{n}=(t_{n-1},t_{n})\times K_{j}$. The space-time tessellation is defined as

[TABLE]

where $G^{n}_{\mathcal{K}}$ denotes the mapping from the space-time reference element $\hat{\mathcal{K}}=(-1,1)^{d+1}$ to the space-time element in physical space $\mathcal{K},$ see Fig. 1. The tessellation $\mathcal{T}_{h}$ of the whole discrete space-time domain is $\mathcal{T}_{h}=\cup_{n=1}^{N}\mathcal{T}_{h}^{n}.$

The space-time FEM discretization is obtained by approximating the test and trial functions in each space-time element in the tessellation $\mathcal{K}^{n}\in\mathcal{T}_{h}^{n}$ with polynomial expansions that are assumed to be continuous within each space-time slab, but discontinuous across the interfaces of the space-time slabs, namely at times $t_{0},t_{1},\ldots,t_{N}.$

The finite element space associated with the tessellation $\mathcal{T}_{h}^{n}$ is defined as:

[TABLE]

where $\hat{\mathcal{P}}_{q}(-1,1)$ and $\hat{\mathcal{P}}_{r}(\hat{K}),$ respectively, represent $q$th-order polynomials on $(-1,1)$ and $r$th-order tensor product polynomials in the reference element $\hat{K}=(-1,1)^{d}$. Finally, define

[TABLE]

Note that the functions in $V_{h}$ are allowed to be discontinuous at the nodes of the partition of the time interval. We will use the notations $u^{n,\pm}=\lim_{s\to 0^{\pm}}u(t_{n}+s).$ Moreover, since $0\not\in I_{1},$ we specify $u^{0,-}=u_{0}(0).$

The space-time dGcG-FEM method applied to problem (5)-(6) can be formulated as: find $u_{h}\in V_{h}$ such that

[TABLE]

holds for all $v_{h}\in V_{h}$ and where $u_{h}^{0,-}=u_{0}(0).$ Here $[u_{h}]_{n}=u_{h}^{n,+}-u_{h}^{n,-}$ denotes the jump of $u_{h}$ at $t_{n}$ and $(\cdot,\cdot)_{\mathcal{K}}$ is the $L^{2}(\mathcal{K})$-inner product on a space-time element. The jumps were added to the weak formulation to ensure weak continuity between time slabs, since the basis functions in dGcG-FEM discretizations are discontinuous at the space-time slab boundary.

Note that throughout this paper the FEM solution will be denoted by $u_{h},$ which should not be confused with the state of the system notation introduced in Section 2. Moreover, it is important to remark that, for $u_{h}\in V_{h},$ the segments ${u_{h}}_{t},$ $t>0$ are not necessarily continuous, but piecewise continuous on $[-\tau_{max},0].$ Denoting the space of piecewise continuous functions on $[-\tau_{max},0]$ by $\hat{X}=PC\left([-\tau_{max},0];Y\right),$ we define the operator $\hat{G}:\hat{X}\to Y$ as

[TABLE]

Then the nonlinear integral operator in (3.1) is equal to $\hat{G}({u_{h}}_{t}).$

The weak formulation (3.1) can be transformed into an integrated-by-parts form, and since we added the jump term at each time level, it is possible to drop the summation over the space-time slabs. Moreover, after integration by parts, (3.1) can be decoupled into a sequence of local problems by choosing test functions that have support only in a single space-time slab $\mathcal{E}^{n}$. Hence we can solve the problem successively, i.e., using the known value $u_{h}(t_{n-1}^{-})$ from the previous space-time slab. The weak formulation for the dGcG-FEM discretization of the neural field equation is the following:

Find $u_{h}\in V_{h}^{n},$ such that for all $v_{h}\in V_{h}^{n}$ the variational equation is satisfied:

[TABLE]

with $\mathcal{K}^{n}\in\mathcal{T}_{h}^{n}$ for $n=1,\dots,N$.

Note here that the delay term may use values from space-time slabs where the solution was computed previously, but also from the current space-time slab, depending on the magnitude of the delay function compared to the time step. This problem will be discussed later in detail.

3.2 How to treat the delay term?

In this section we discuss the dGcG-FEM approximation of the delay term in the weak formulation (3.1). Introduce the approximation

[TABLE]

into (3.1) and set the test function $v_{h}(t,x)\mid_{\mathcal{K}}=\psi_{i}^{\mathcal{K}}(t,x),$ $i\in\{1,\dots,N_{p}\},$ with $N_{p}$ the number of degrees of freedom in element $\mathcal{K}$ and $\psi_{i}^{\mathcal{K}}$ standard Lagrange tensor product basis functions. The delay term becomes

[TABLE]

All integrals in the weak formulation are evaluated using Gaussian quadrature rules. Let us fix a quadrature point $(t_{q},x_{q})\in\mathcal{K}^{n}$ in a space-time element and let $\tau_{max}=\max_{(x,r)\in\bar{\Omega}\times\bar{\Omega}}\tau(x,r),$ as before. To compute the integral over a space element $L$ in (3.2), consider a space quadrature point $r_{qs}\in\Omega,$ and distinguish three cases for the time delay $t_{q}-\tau(x_{q},r_{qs})$, see Figure 2:

Case 1. If $-\tau_{max}\leq t_{q}-\tau(x_{q},r_{qs})\leq 0,$ then the solution at this time level is given by the initial solution, i.e., $u_{h}(t_{q}-\tau(x_{q},r_{qs}),r_{qs})=u_{0}(t_{q}-\tau(x_{q},r_{qs}),r_{qs}).$

Case 2. When $t_{q}-\tau(x_{q},r_{qs})\geq t_{n-1},$ then the delay term (3.2) is implicit since we remain in the same space-time slab $\mathcal{E}^{n},$ where the solution is unknown. Hence, when the delay time is small enough compared to the time step, an additional Newton method needs to be incorporated for the solution of the nonlinear system.

If we introduce the finite element approximations for $u_{h}$ and $v_{h}$ also into the other terms in the weak formulation (3.1), then we obtain for all $\mathcal{K}\in\mathcal{T}_{h}^{n}$

[TABLE]

where $\hat{u}_{j}^{\mathcal{K},n-1}$ are the coefficients of space-time element $\mathcal{K}$ in the space-time slab $\mathcal{E}^{n-1}.$

Case 3. When $0\leq t_{q}-\tau(x_{q},r_{qs})<t_{n-1},$ then the delay term is explicit since we go back to a previous space-time slab, where the FEM solution is already computed.

4 Error analysis

In this section we give an a-priori error analysis for the space-time dGcG-method (3.2). In the error analysis we will use a slightly modified version of the temporal interpolation functions defined in Proposition 4.1, [12]. First, define the space

[TABLE]

with $I_{n}=(t_{n-1},t_{n}]$ and $|I_{n}|=k_{n}.$ Note that these functions are allowed to be discontinuous at the nodes of the partition of the time interval, but continuous from the left in each subinterval $I_{n},$ i.e., $w(t_{n})=\lim_{t\to t_{n}^{-}}w(t).$ For the restriction of the functions in $S_{k}$ to $I_{n},$ we use the notation $S_{k}^{n}.$ Define the temporal polynomial interpolant

[TABLE]

as follows, see also [12].

Proposition 4.1**.**

Let $\tilde{u}=T_{k}u\in S_{k}$ be the time-interpolant of $u\in C\left([0,T],Y\right)\cap H^{q+1}\left([0,T],Y\right),$ $q\geq 0,$ with the following properties:

[TABLE]

The interpolation error then can be estimated as

[TABLE]

where $\partial_{s}^{q+1}$ denotes the $(q+1)$-th order derivative w.r.t. time, $k_{t}=|I_{t}|$ and the norm $\|\cdot\|=\|\cdot\|_{L^{2}(\Omega)}$ hereafter.

Observe that $\tilde{u}$ interpolates exactly at the nodes and the interpolation error is orthogonal to polynomials of degree at most $q-1.$ For constant polynomials ($q=0$) condition (17) is not used.

Next, define the spatial interpolant. Let $W_{h}$ be the space of tensor product polynomials of degree up to $r\geq 0$ on each space element $K_{j},$ i.e.,

[TABLE]

where $G^{n}_{K}$ denotes the mapping from the reference element $\hat{K}=(-1,1)^{d}$ to the element $K\in\bar{\mathcal{T}}_{h}$ in physical space. Let

[TABLE]

be the $L^{2}$-projection to the (spatial) finite element space, defined as $\left(P_{h}v,w_{h}\right)=\left(v,w_{h}\right)$ for all $w_{h}\in W_{h}.$ We use the standard interpolation estimate in space (see e.g. [2], [3])

[TABLE]

where $\|\cdot\|_{r+1}=\|\cdot\|_{H^{r+1}(\Omega)},$ $h$ denotes the maximal space element diameter as before, and the constant $C$ is independent of $h$ and $v.$

In the error analysis we also need the interpolation of the initial segment of the solution. Let the given initial function be $u_{0}\in X\cap H^{q+1}\left([-\tau_{max},0];H^{r+1}(\Omega)\right)$ for some $q,r\geq 0.$ Use a partition of the interval $[-\tau_{max},0]$ into $M$ subintervals $J_{i}$ of length $k_{i},$ respectively. On each $J_{i}$ we use the same temporal interpolation $\tilde{u}_{0}=T_{k}u_{0}$ of $u_{0},$ as introduced in Proposition 4.1. Then for all $s\in J_{i}$ we have

[TABLE]

where we use that the operator norm of the Lagrange interpolation $P_{h}$ is bounded, see [2], [4].

We will also need an estimate of the integral of the interpolation error on the partition of the initial segment. There exists $C>0$ generic constant (independent of the solution and mesh size), such that

[TABLE]

where we denoted the norm

[TABLE]

Next, we state the main result of the a-priori error analysis of the dGcG discretization (3.1) for the neural field equations.

Theorem 4.1**.**

Let $u\in C^{1}\left([0,T);Y\right)\cap H^{q+1}\left([0,T];H^{r+1}(\Omega)\right)$ be the solution of (4) for some $q,r\geq 0$ and with initial state $u_{0}\in X\cap H^{q+1}\left([-\tau_{max},0];H^{r+1}(\Omega)\right),$ and let $u_{h}\in V_{h}^{n}$ be the solution of (3.1). Then

[TABLE]

holds for $t_{N}\geq 0,$ $N$ the number of time slabs, where $C$ is a positive constant independent of the time step $k_{n}=t_{n}-t_{n-1}$ and the maximal space element diameter $h$. Here $m(n)\leq N-1,$ $1\leq n\leq N,$ is the multiplicity how many times we visited the interval $I_{n}$ due to the delay term.

Proof.

Let us decompose the error of the numerical discretization into the sum

[TABLE]

with $\theta$ the discretization error and $\rho$ the interpolation error. When $t\in[-\tau_{max},0],$ we only have the interpolation error of the given initial solution $u_{0},$ that is, $\theta(t,x)=0$ and $\rho(t,x)=P_{h}\tilde{u}_{0}(t)(x)-u_{0}(t,x).$ From here on, we suppress the spatial dependence where it is clear from the context. Since $\tilde{u}$ interpolates exactly at the nodes $t=t_{n-1},$ we have that

[TABLE]

holds for all $n\geq 1.$ Here the constant $C$ is independent of $h,$ see e.g. [2]. When we are in the interior of a time interval $I_{j}$, we decompose $\rho$ to be able to use the bound on the interpolation error in time and space, respectively, as in (4)

[TABLE]

for any $t\in I_{j}$ and $j=1,\dots,N.$ It is, therefore, sufficient to bound $\theta^{N}=\theta(t_{N}).$ Since both $u_{h}$ and $u$ satisfy the weak formulation (3.1) with $G$ and $\hat{G},$ respectively, we obtain that for all $v\in V_{h}^{n}$

[TABLE]

The variational equation (4) holds for any partition of the time interval $I_{n}$, hence the following equation is also valid for any $t\in(t_{n-1},t_{n}]$

[TABLE]

where $I_{t}=(t_{n-1},t].$ Using the assumptions on the interpolant, some terms in (4) will cancel, i.e., for all $t\in I_{n}$

[TABLE]

Let $v=2\theta\in S_{k}^{n}$ in (4). Then for each $I_{n}$ and $t\in I_{n}$ the following holds

[TABLE]

This may be further written as

[TABLE]

Using the Schwarz inequality and the inequality $2ab\leq\epsilon^{2}a^{2}+\frac{1}{\epsilon^{2}}b^{2}$ we obtain

[TABLE]

Since the nonlinearity $S$ is Lipschitz continuous with some Lipschitz constant $C_{S},$ we can estimate the nonlinear term as

[TABLE]

Let us estimate the first term on the right hand side of (4) as

[TABLE]

where we used the Schwarz inequality in each estimation step and $|\Omega|=\text{vol}(\Omega)$. Next, since $0<\tau(x,r)\leq\tau_{max},$ and $J(x,r)\leq\|J\|_{C},$ for all $(x,r)\in\bar{\Omega}\times\bar{\Omega},$ the following estimate is valid

[TABLE]

Hence we can further estimate (4) as

[TABLE]

Similarly as in (4) and (4), for the last term in (4) we obtain

[TABLE]

After introducing the above estimates into (4) we obtain that for all $t\in I_{n},$

[TABLE]

is valid for all $n\geq 1.$ Divide by $1-\epsilon^{2},$ where $0<\epsilon<1,$ and denote by

[TABLE]

Then, inequality (4) can be written as

[TABLE]

where $\eta(t)=\|\theta(t)\|^{2}.$ Apply Grönwall’s inequality to (39) to obtain

[TABLE]

When $t=t_{n},$

[TABLE]

where

[TABLE]

Note that the only time-dependent term in $\omega_{n}(t)$ is $\|\rho(t)\|^{2}.$ Hence, the integral term in (41) can be estimated as

[TABLE]

Therefore, we obtain for $n\geq 1$ that

[TABLE]

Let us recall that

[TABLE]

and observe that the right hand side of (44) can be estimated by the bound of the interpolation error and the bound of the integral of $\eta(t)$ over earlier time intervals, i.e., for $t\leq t_{n-1}.$ Hence we can write

[TABLE]

where $C_{i},$ $i=1,\dots,4$ depend on the parameters $\alpha,\beta,\epsilon,\|J\|_{C},$ $|\Omega|$ and $k_{n},$ such that $C_{i}=O(1)$ as $k_{n}\to 0.$

By integrating (40) we obtain the following general formula

[TABLE]

where we used that $\eta(s)=\|\theta(s)\|^{2}=0$ for $s\in[-\tau_{max},0]$ and (4) in the last inequality. Here $m=m(n)\leq n-1$ is the index of the interval $I_{m}$ for which $t_{n-1}-\tau_{max}\in I_{m}.$

As we can see, the integral of $\eta(s)$ can be bounded by the integral of $\|\rho(s)\|,$ hence in (4) we have

[TABLE]

We can use (4) to bound the integral of $\|\rho(s)\|$ as follows

[TABLE]

For $n=1,$ combining (4) with (4), (4), (4), (4) and (49) and using that $\eta(0)=0,$ we find that there exists a generic constant $C$, independent of the time step $k_{1}$ and the spatial mesh size $h,$ such that

[TABLE]

For $n=2,$ using again (4) and then (4), (4), (4) and (4), we find that there is a constant $C$, such that

[TABLE]

where $m(i)$ and $m(j)$ are the multiplicity how many times we visited the interval $J_{i}$ and $I_{j},$ respectively, in the integral of $\|\rho(s)\|$ over the delay interval. If $\tau_{max}$ is large compared to the time step, then $m(i)$ is consequently also larger.

We can repeat this procedure for the subsequent time intervals, which completes the proof of the theorem. ∎

5 Numerical simulations

In this section we present applications of the FEM discretization to the neural field equations, starting with delay differential equations with constant delay.

5.1 DDE with constant delay

Here we study the numerical solution of equations of the form

[TABLE]

with $\tau>0$ a constant delay and $f:\mathbb{R}^{2}\to\mathbb{R}$ linear, given by $f\left(u(t),u(t-\tau)\right)=-\alpha u(t)+u(t-\tau).$

To verify our results on the error analysis, we compare the time-discontinuous Galerkin FEM solution (dG(1)) using linear basis functions, with the exact solution computed for some delay intervals. Let the history function be $u_{0}(s)=-s,$ $s\in[-\tau,0]$ and $\tau=2.$ Figure 3 illustrates the solution when $\alpha=1,$ for which we know that it converges to a non-zero steady state. We set $k_{n}=k$ for all $n$ and distinguish two cases. First, when $\tau/k$ is not an integer, then the dG(1) method is second order accurate, which is consistent with our result on the error estimate. When $\tau/k$ is, however, integer then we observe a higher order accuracy of order three. Figure 3 shows both cases.

The numerical integration of delay differential equations is very sensitive to jump discontinuities in the solution or in its derivatives. Such discontinuity points are referred in the literature as breaking points, [1]. In case of constant delay, the breaking points are $\xi_{n}=n\tau$ for $n=1,2,\dots.$ The best procedure to guarantee the required accuracy is to include these breaking points in the set of mesh points. In our example, the derivative of the solution has discontinuity at $t=0.$ When the breaking points are also mesh points, i.e., when $\tau/k$ is integer then the error in the discontinuous Galerkin method is of order $O(k^{q+2}),$ which is of superconvergent order.

5.2 Integro-differential equations

One important result is the successful treatment of the fully implicit case, i.e., when the delay is zero. Hence, consider the integro-differential equation, obtained by removing the delay term in (2) and adding a given, sufficiently smooth, source term $g$

[TABLE]

with initial condition $u(0,x)=u_{0}(x).$ In our numerical simulation, we further simplify this equation by taking $J(x,r)=1,$ and $S\left(u(t,r)\right)=u(t,r)$ linear.

As a first example, we take $g=0,$ $\alpha=1$ and $u_{0}(x)=x.$ The exact solution of (53) is $u(t,x)=xe^{-t},$ which converges to zero as $t\to\infty$ for every $x\in\Omega.$ The time interval is divided equidistantly with time step $k_{n}=k.$ With this example we want to demonstrate that the time accuracy is not destroyed when we add a spatial integral term. The dGcG-FEM solution using linear basis functions, both in space and time, and the time accuracy for this example are plotted in Figure 4. We observe that the error in the dGcG-FEM method is of superconvergent order.

In the second example, we study the time accuracy when the solution of (53) is periodic in time. Take $g(t,x)=x\cos t+\alpha x\sin t,$ $\Omega=[-1,1].$ Then the exact solution $u(t,x)=x\sin t$ satisfies the initial condition $u(0,x)=0.$ We compute the error of the solution at several time levels in a period and observe the same temporal accuracy as in the first example.

5.3 The neural field equations

In this section we demonstrate the dGcG(1) method for an example analyzed in [14], both analytically and numerically. Consider the single population model (2), when the space is 1-dimensional. Space and time are rescaled such that $\bar{\Omega}=[-1,1]$ and the propagation speed is 1. This yields

[TABLE]

In this case, equation (2) becomes

[TABLE]

The connectivity and activation functions are, respectively,

[TABLE]

and

[TABLE]

Hopf bifurcations play an important role in the analysis of neural field equations. By choosing the steepness parameter $\sigma$ of the activation function as bifurcation parameter, we can simulate, using the dGcG(1) scheme, the space-time evolution of the solution beyond a Hopf bifurcation. As in [14], we choose the parameters $\alpha=1$ and $\sigma=6$ in the activation function (57) and the delay $\tau_{0}=1.$ In this simulation the connectivity function has a bi-exponential form

[TABLE]

with $\hat{c}_{1}=3.0,\ \hat{c}_{2}=-5.5,\ \mu_{1}=0.5,\ \mu_{2}=1.0.$ Figure 5 shows the time evolution of the system and Figure 6 is a surface plot of the numerical solution.

The initial function for this simulation is $u(t,x)=\epsilon=0.01,$ $t\in[-\tau_{max},0].$ Note that, because the size of the delay is relatively large compared to the time step, we do not need to linearize the system to solve the algebraic equations with a Newton method.

5.4 Neural fields with spatial inhomogeneity

Consider the neural field equation (55) with the locally changed connectivity

[TABLE]

where $J$ is given in (58) with the same parameters and $\tilde{\Omega}\subset\Omega.$ The activation function is given in (57) with the bifurcation parameter $\sigma=4,$ chosen below the threshold for Hopf bifurcation to occur in the homogeneous case, see [14]. In Figures 7, 8 and 9, we compare the solution of the system with homogeneous kernel, with the solution where we have locally changed the connectivity, specifically in one element, i.e., $\tilde{\Omega}=K\in\bar{\mathcal{T}}_{h}.$ Our simulations show that while the solution converges to a steady state in the homogeneous case, in the inhomogeneous case the solution becomes periodic ($\omega=15$). This is a new phenomenon observed in the one dimensional case. It requires, however, further bifurcation analysis in the two-parameter space $(\sigma,\omega).$

6 Concluding remarks

In this article we have presented a new space-time dGcG-FEM to solve delay integro-differential equations with space dependent delays. The main result is an a-priori error estimate of the space-time dGcG method, which also shows that the method is numerically stable. We demonstrated that by using a dGcG method we can handle general connectivity, synaptic activation and delay functions, and do not need to make any restriction on spatial dimension or shape of the domain. This makes it possible to extend our model to more general domains as well as more populations in the system, which are particularly interesting for our applications.

Acknowledgments

The first author was supported by the Hungarian Scientific Research Fund, Grant No. K109782. The ELI-ALPS project (GOP-1.1.1.-12/B-2012-0001, GINOP-2.3.6-15-2015-00001) is supported by the European Union and co-financed by the European Regional Development Fund.