Numerical Solution of Stochastic Neural Fields with Delays

TL;DR

This paper presents a novel numerical algorithm for solving stochastic neural field equations with delays, demonstrating how noise influences stationary solutions under spatially heterogeneous inputs.

Contribution

The paper introduces a new numerical method for stochastic neural fields with delays, enabling analysis of noise effects on stationary solutions.

Findings

Noise affects the stability of stationary solutions.

Spatial heterogeneity influences neural field dynamics.

The algorithm successfully simulates stochastic neural fields.

Abstract

We introduce a new numerical algorithm for solving the stochastic neural field equation (NFE) with delays. Using this algorithm we have obtained some numerical results which illustrate the effect of noise in the dynamical behaviour of stationary solutions of the NFE, in the presence of spatially heterogeneous external inputs.

Equations

$$\frac{\partial }{\partial t}V(x,t)=I(x,t)-\alpha V(x,t)+{\int }_{\Omega }K(|x-y|)S(V(y,t-\tau (x,y)))dy,$$

$$t\in [0,T],x\in \Omega \subset \mathbb{R},$$

$$d{U}_t(x)=\left(I(x,t)-\alpha U_t(x)+{\int }_{\Omega }K(|x-y|)S(U_t(y))dy\right)dt+ϵd{W}_t(x),$$

$$d{U}_t(x)=\left(I(x,t)-\alpha U_t(x)+{\int }_{\Omega }K(|x-y|)S(U_{t-\tau }(y))dy\right)dt+ϵd{W}_t(x),$$

$$U_t(x)=U_0(x,t),t\in [-{\tau }_{max},0],x\in \Omega ,$$

$$U_t(x)=\sum _{k=0}^{\infty }{u}_t^k{v}_k(x),$$

$$(d{U}_t,v_i)=\left[(I(x,t),v_i)-\alpha (U_t,v_i)+\left({\int }_{\Omega }K(|x-y|)S(U_{t-\tau }(y))dy,v_i\right)\right]dt+ϵ(d{W}_t,v_i).$$

$$d{W}_t(x)=\sum _{k=0}^{\infty }{v}_k(x){\lambda }_{k}d{\beta }_t^k,$$

$$E{W}_t(x)W_s(y)=\min (t,s)\frac{1}{2\xi }\exp \left(\frac{-\pi }{4}\frac{|x-y{|}^2}{{\xi }^2}\right),$$

$${\lambda }_k^2=\exp \left(-\frac{{\xi }^2{k}^2}{4\pi }\right).$$

$$d{u}_t^i=\left[(I(x,t),v_i)-\alpha u_t^i+(KS{)}^i({\bar{u}}_{t-\tau })\right]dt+ϵ{\lambda }_{i}d{\beta }_t^i,$$

$$(KS{)}^i({\bar{u}}_t)={\int }_{\Omega }{v}_i(x)\left({\int }_{\Omega }K(|x-y|)S\left(\sum _{k=1}^{\infty }{u}_{t-\tau }^k{v}_k(y)\right)dy\right)dx.$$

$$(KS{)}^i({\bar{u}}_{t-\tau })={\int }_{\Omega }{v}_i(x)K(|x-y|)\left({\int }_{\Omega }S\left(\sum _{k=0}^{\infty }{u}_{t-\tau }^k{v}_k(y)\right)dy\right)dx.$$

$$U_t^N(x)=\sum _{k=0}^{N-1}{u}_t^{k,N}{v}_k(x).$$

$$d{u}_t^{i,N}=\left[(I(x,t),v_i)-\alpha u_t^{i,N}+(KS{)}^{i,N}({\bar{u}}_{t-\tau })\right]dt+ϵ{\lambda }_{i}d{\beta }_t^i,$$

$$(KS{)}^{i,N}({\bar{u}}_{t-\tau })=h^2\sum _{j=0}^N{v}_i(x)\left(\sum _{l=1}^{N}K(|x_l-x_j|)S\left(\sum _{k=1}^N{u}_{t-\tau }^k{v}_k(x_j)\right)\right)$$

$$u_j^{k,N}\approx u_{t_j}^{k,N}.$$

$$u_{j+1}^{i,N}=u_j^{i,N}+h_t\left[(I(x_i,t_j),v_i)-\alpha u_{j+1}^{i,N}+(KS{)}^{i,N}({\bar{u}}_{t_j-\tau })\right]+\sqrt{h_t}ϵ{\lambda }_i{w}_i,$$

$$u_{j+1}^{i,N}=\frac{u_j^{i,N}+h_t\left[(I(x_i,t_j),v_i)+(KS{)}^{i,N}({\bar{u}}_{t_j})\right]+\sqrt{h_t}ϵ{\lambda }_i{w}_i}{1+\alpha h_t}.$$

$$t_{j-d-1}\le t_j-\tau \le t_{j-d}$$

$$h_t{\delta }_t=\tau -d{h}_t.$$

$$u_{t_{j}-\tau}\approx\left\{\begin{array}[]{c}u(t_{j-d}),\quad\makebox{if}\,\delta_{t}<0.5,\\ u(t_{j-d-1}),\quad\makebox{if}\,\delta_{t}\geq 0.5.\end{array}\right.$$

$$v_k(x)=\exp (ikx),k=0,1,..,N.$$

$$K(x)=2\exp (-0.08x)\left(0.08\sin (\pi x/10)+\cos (\pi x/10)\right).$$

$$I(x)=-3.39967+8\exp \left(-\frac{x^2}{18}\right).$$

$$U_{max,max}(t)=\underset{s\in \{1,\dots ,100\}}{\max }\underset{i\in \{1,\dots ,100\}}{\max }u(s,x_i,t);$$

$$U_{min,max}(t)=\underset{s\in \{1,\dots ,100\}}{\min }\underset{i\in \{1,\dots ,100\}}{\max }u(s,x_i,t);$$

$$U_{max,min}(t)=\underset{s\in \{1,100\}}{\max }\underset{i\in \{1,\dots ,100\}}{\min }u(s,x_i,t);$$

$$U_{min,min}(t)=\underset{s\in \{1,\dots ,100\}}{\min }\underset{i\in \{1,\dots ,100\}}{\min }u(s,x_i,t);$$

$$E(u(x,t))\approx \frac{1}{100}\sum _{s=1}^{100}u(s,x,t);$$

Full text

1 Introduction

Neural Field Equations (NFE) are a powerful tool for analysing the dynamical behaviour of populations of neurons. The analysis of such dynamical mechanisms is crucially important for understanding a wide range of neurobiological phenomena [2]. In this work we will be concerned with the NFE in the form

[TABLE]

[TABLE]

where $V(x,t)$ (the unknown function) denotes the membrane potential in point $x$ at time $t$; $I(x,t)$ represents the external sources of excitation; $S$ is the dependence between the firing rate of the neurons and their membrane potentials (sigmoidal or Heaviside function); and $K(|x-y|)$ gives the connectivity between neurons at $x$ and $y$, $\alpha$ is a constant (related to the decay rate of the potential); $\tau(x,y)>0$ is a delay, depending on the spatial variables (it results from the finite propagation speed of nervous stimulus in the brain).

Equation (1) (without delay) was introduced first by Wilson and Cowan [15], and then by Amari [1], to describe excitatory and inhibitory interactions in populations of neurons.

Intensive studies of Hopf bifurcations occuring in neural fields have been carried out in the last decade. In [2] the authors investigate the occurence of spatially localized oscillations (or breathers) in two-dimensional neural fields with excitatory and inhibitory interactions; in [13] the authors obtain sufficient conditions for the stability of stationary solutions of neural field equations; the dependency of the stationary solutions of NFE with respect to the stiffness of the nonlinearity and the contrast of external inputs is studied in [14].

Though the above mentioned results were obtained analytically, numerical simulations play a fundamental role in studying brain dynamics. Thus, the availability of efficient numerical methods is an important ingredient for improving the understanding of neural processes. Concerning equation (1), numerical approximations were obtained in [3]. The computational method applies quadrature rule in space to reduce the problem to a system of delay differential equations, which is then solved by a standard algorithm for this kind of equations. A more efficient approach was recently proposed in [6] [7], where the authors introduce a new approach to deal with the convolution kernel of the equation and use Fast Fourier Transforms to reduce significantly the computational effort required by numerical integration. Recently, a new numerical method for the approximation of two-dimensional neural fields has been introduced, based on an implicit second order scheme for the integration in time and using Chebyshev interpolation to reduce the dimensions of the matrices [11]. Some applications of this algorithm to Neuroscience problems have been discussed in [12].

As in other sciences, in Neurobiology it is well-known that better consistency with some phenomena can be provided if the effects of random processes in the system are taken into account. In the recent work of Kühn and Riedler [9], the authors study the effect of additive noise in Neural Field Equations. With this purpose they introduce the stochastic integro-differential equation

[TABLE]

where $t\in[0,T]$, $x\in\Omega\subset\mathbb{R}^{n}$, $W_{t}$ is a Q-Wiener process.

The main goal of the present work is to analyse the effect of noise in certain neural fields, which allow different types of stationary solutions. In this case we consider the following modification of equation (2)

[TABLE]

where, as in the deterministic case, $\tau$ is a delay, depending on the distance $|x-y|$. Equation (3) is completed with an initial condition of the form

[TABLE]

where $U_{0}(x,t)$ is some given stochastic process, $\tau_{max}$ is the maximum value of the delay ( $\tau_{max}=|\Omega|/v$), where $v$ is the propagation speed of the signals. We assume that $U_{t}(x)$ satisfies periodic boundary conditions in space. We will consider domains of the form $\Omega=[-l,l]$, including the limit case when $l\rightarrow\infty$.

2 Numerical Approximation

To construct a numerical approximation of the solution of (2) in the one-dimensional case, we begin by expanding the solution $U_{t}(x)$ using the Karhunen-Loeve formula:

[TABLE]

where $v_{k}$ are the eigenfunctions of the covariance operator of the noise in (2), which form an orthogonal system (their explicit form is indicated below). To derive a formula for the coefficients $u_{t}^{k}$ we take the inner product of equation (2) with the basis functions $v_{i}$:

[TABLE]

We expand $dW_{t}$ as

[TABLE]

where the functions $\beta_{t}^{k}$ form a system of independent white noises in time and $\lambda_{k}$ are the eigenvalues of the covariance operator of the noise. As an important particular case, we consider the one described in [9], p.7. In this case the correlation function satisfies

[TABLE]

where $\xi$ is a parameter modeling the spatial correlation length. In this case, if $\xi<<2\pi$, the eigenvalues of the covariance operator satisfy

[TABLE]

By substituting (5) into (6), taking into account (7) and the ortogonality of the system $v_{k}$, we obtain

[TABLE]

where $(KS)^{i}(\bar{u}_{t})$ denotes the nonlinear term of the system:

[TABLE]

Due to the symmetry of the kernel, the last integral may be rewritten as

[TABLE]

When using the Galerkin method, we define an approximate solution by truncating the series expansion (5)

[TABLE]

Then the coefficients $u_{t}^{k,N}$ satisfy the following nonlinear system of stochastic delay differential equations:

[TABLE]

where $(KS)^{i,N}(\bar{u}_{t})$ is given by

[TABLE]

$i=0,...,N-1$. In this case we are introducing in $[-L,L]$ a set of $N+1$ equidistant gridpoints $x_{j}=-l+j*h$, $j=0,...,N$, where $h=2l/N$, and using the rectangular rule to evaluate the integral in (10).

Since the problem has been reduced to the system (12) (a system of nonlinear stochastic delay differential equations), we can apply the Euler-Maruyama method to the solution of this system. For the discretisation in time, we introduce on the interval $[0,T]$ a uniform mesh with step size $h_{t}$, such that $t_{j}=jh_{t}$, $j=0,1,...,n$. Then the solution $u_{t}^{k,N}$ of (12) will be approximated by a vector $(u_{1}^{k,N},u_{2}^{k,N},...,u_{n}^{k,N})$, where

[TABLE]

In these notations, the Euler-Maruyama method may be written as

[TABLE]

where $w_{i}$ is a random variable with normal distribution ($w_{i}=N(0,1)$), $j=0,....,n$, $i=0,...,N-1$. In the right-hand side of (14) we have written $\alpha u_{j+1}^{i,N}$, meaning that we are using a semi-implicit version of the Euler-Maruyama method.

Further we can rewrite equations (14) in the form

[TABLE]

The meaning of the inner product in (15) will be explained below. In order to compute $u_{t_{j}-\tau}^{k}$ we should take into account that $\tau=\frac{|x_{k_{1}}-x_{k_{2}}|}{v}$ (the time spent by the signal to travel between $x_{k_{1}}$ and $x_{k_{2}}$). In general $\tau$ may not be a multiple of $h_{t}$. Let $d$ and $\delta_{t}$ be the integer and the fractional part of $\frac{\tau}{h_{t}}$. In this case, we have

[TABLE]

and

[TABLE]

The needed value of the solution $u_{t_{j}-\tau}$ in (14) is then approximated by

[TABLE]

Concerning the choice of the basis functions, we consider a set of orthogonal functions, similar to the onde described in [9], page 7. More precisely, we define

[TABLE]

Note that with this choice of the basis functions the inner product in (15) and the sums in (13) can be interpreted as the Discrete Fourier Transform (DFT). In particular, the set of inner products $(I(x,t_{j}),v_{i})$ , $i=1,...,N$ may be seen as the DFT of the vector $I_{N}$, which contains the values of the function $I(x,t)$ at the grid points $x_{k}=-l+kh$, $k=1,..,N$. In this case $t$ is fixed ($t=t_{j}$). Therefore, these inner products can be evaluated efficiently by the Fast Fourier Transform (FFT).

3 Numerical Examples

3.1 Noise effects on stationary multi-bump solutions

We have applied our algorithm to analyse the effect of noise on the formation of multi-bump solutions in dynamic neural fields, in the presence of space dependent external stimulli. In [10], the authors have investigated the formation of regions of high activity (bumps) in neural fields, which can be switched ’on’ and ’off’ by transient stimulli. The stability analysis of such patterns in the deterministic case was carried out in [4] and in [5]. The effects of noise on such stationary pulse solutions was studied in [8]. In this case, the firing rate function $S(x)$ is the Heaviside function; the connectivity kernel is given by

[TABLE]

In [5] it was proved that such oscillatory connectivity kernels support the formation of stationary stable multi-bump solutions, induced by external inputs. In this case the external input has the form

[TABLE]

This specific example was considered in [4], p. 37. The parameters of the numerical approximation are $N=100,h=1,l=50;n=200,h_{t}=0.02.$ We start by considering the deterministic case. Using our code with $\epsilon=0$, we are able to reproduce three kinds of stationary solutions, which are displayed in Fig.1 and Fig.2.

As the first experiment with the stochastic case we have performed a simulation with 100 paths, with noise level $\epsilon=0.01$ starting with the initial condition $U_{0}(x,t)\equiv 0$, over the time interval $t\in[0,4]$. Our aim was to investigate the evolution of the paths of the stochastic equation and their relation with the stationary multi-bump solutions observed in the deterministic case. We remark that each solution $V(x,t)$ of the deterministic equation (1), as $t$ tends to infinity, tends to a certain stationary solution $u(x)$, which is characterized by a certain maximal value $u_{max}=\max_{x\in[-l,l]}u(x)$ (which is usually located at the origin) and a minimal value $u_{min}=\min_{x\in[-l,l]}u(x)$ (which is attained at two symmetric points). Usually a solution may have other local minima and maxima, whose absolute values do not exceed those of $u_{min}$ and $u_{max}$. The greater is the number of bumps, the higher is the value of $u_{max}$ and the lower is the value of $u_{min}$. Note that $u_{max}$ and $u_{min}$ are respectively the limits, as $t\rightarrow\infty$, of $u_{max}(t)$ and $u_{min}(t)$, that is, $u_{max}(t)=\max_{x\in[-l,l]}V(x,t)$ and $u_{min}(t)=\min_{x\in[-l,l]}V(x,t)$. We have used these properties to analyse the paths of the stochastic equation. Let $u(s,x,t)$ denote the approximate value of $U_{t}(x)$, given by the $s$-th path. We use the following notations:

[TABLE]

[TABLE]

[TABLE]

[TABLE]

We consider the following approximations of mathematical expectations:

[TABLE]

[TABLE]

[TABLE]

In Fig. 3 we can see the time evolution of the solution maximum (left) and the time evolution of the solution minimum (right). In the first case the graphs of $U_{max,max}(t)$, $U_{min,max}(t)$ and $E_{max}(t)$ are displayed; the second graphic shows values of $U_{max,min}(t)$, $U_{min,min}(t)$ and $E_{min}(t)$ . These figures show that the average value of the maximum increases with time, while the average value of the minimum decreases; as a result the average amplitude of the oscillations (height of the bumps) increases with time. Moreover, as it could be expected the dispersion (maximal difference between values of different paths) also increases with time. As it happens with $u_{min}(t)$ and $u_{max}(t)$ in the deterministic case, the values $E_{min}(t)$ and $E_{max}(t)$ stabilize as $t$ increases. For $j=200(t_{j}=4)$, their values are close to the ones of $u_{min}$ and $u_{max}$, in the case of a deterministic one-bump solution.

Fig.4 shows the distribution (at $t=4$) of $u_{max}$ (left) and $u_{min}$ (right), at $t=4$, for the diferent paths. We see that the maxima are concentrated on the range $[15.8,16.6]$, which corresponds to the stationary one-bump solution (see Fig. 1), and $[20,21.2]$, which corresponds to three-bump and five-bump solutions (see Fig. 2). The minima are concentrated on the interval $[-8.2,-7.4]$, which corresponds to the stationary one-bump solution, and $[-14.,-12.5]$, which corresponds to three-bump and five-bump solutions of the deterministic equation. This suggests that under the effect of not very strong noise, the most probable values of the solution of the stochastic equation are close to the ones of the stationary solutions of the deterministic equation.

As a second numerical experiment, we have performed a simulation with 100 paths, with noise level $\epsilon=0.01$, over the time interval $t\in[0,4]$, taking as initial condition $U_{0}(x,t)$ a stationary one-bump solution (of the type represented in Fig.1).

In Fig. 5 we again observe the evolution of $U_{max,max}(t)$, $U_{min,max}(t)$ and $E_{max}(t)$ (left-hand side); the values of values of $U_{max,min}(t)$, $U_{min,min}(t)$ and $E_{min}(t)$ (right-hand side). As in the previous case (Fig. 3) we observe that the average amplitude of the solutions oscillations increases with time; however for this level of noise ($\epsilon=0.01$), the amplitude of oscillations tends to stabilize and the mean value of the stochastic solution remains close to the one-bump deterministic solution.

In Fig. 6, the graph of the average solution $E(u(x,t))$ at $t=4$ is displayed, between the graphs containing the minimal and maximal values (of all the paths); these graphs correspond to the simulation starting with the stationary one-bump solution and with noise level $\epsilon=0.01$.

In Fig. 7 we can see the distribution (at $t=4$) of the values of $u_{max}$ (left) and $u_{min}$ (right), for the same initial values and noise level. We see that the maxima are concentrated on the range $[15.8,16.6]$ and the minima are concentrated on the range $[-9.4,-8.3]$. This means that with level noise $\epsilon=0.01$ when the initial value of the simulation is a deterministic one-bump solution, the most probable values of the stochastic solution are concentrated near the ones of this stationary solution.

As a third numerical experiment, we have performed a simulation with 100 paths, with level noise $\epsilon=0.05$, over the time interval $t\in[0,4]$, taking again as initial condition $U_{0}(x,t)$ a stationary one-bump solution.

The evolution of the maximum and minimum of the solutions is displayed in Fig. 8; in each case we can see the graphs of the average value of the 100 paths ($E_{max}(t)$ and $E_{min}(t)$, respectively). In Fig. 9, the graph of the average solution at $t=4$ is plotted, between the graphs containing the minimal and maximal values (of all the paths). As in the previous cases we observe that the average amplitude of the solutions oscillations increases with time; moreover, for this level of noise ($\epsilon=0.05$), the amplitude of oscillations increases so much that some paths behave like three-bump solutions or five-bump solutions.

In Fig. 10 the distribution (at $t=4$) of the maxima and minima (respectively) of the diferent paths are displayed. We see that most of the maxima are on the range $[16,25]$ and most of the minima are on the range $[-19,-9]$. This reflects the fact that with noise level $\epsilon=0.05$ many paths starting from a stationary one-bump solution can after some time take values that are characteristic of three-bump or five-bump deterministic solutions.

The algorithm was implemented in Mathematica [16] and the computations were performed in a PC with a 1.7Ghz processor and 8 Gb of installed memory (RAM). The computation of the numerical examples presented in this section (with 100 paths) takes about 1 hour.

4 Conclusions

In this paper we have introduced a new numerical algorithm for the approximation of the stochastic neural field equation with delay. This numerical algorithm uses the Galerkin method and is inspired in the Kühn and Riedler’s approach [9]. The choice of the basis functions and grid points allows the use of the Fast Fourier Transform to perform summations, which improves significantly the efficiency of the algorithm.

To test the algorithm we have applied it to the numerical solution of a neural field, in the presence of external stimulli, where stationary one-bump and multi-bump solutions are known to exist in the deterministic case. The numerical results suggest that for a low level of noise the trajectories of the stochastic equation are concentrated near the stationary solutions of the deterministic one. In particular, if the initial condition is the trivial solution, in the deterministic case the solution tends with time to a one-bump deterministic solution. But in the presence of not very strong noise the trajectories of the stochastic equation split into several classes, each of them close to a different stationary solution of the deterministic equation.

In conclusion, we can say that our results are consistent with the study of Kilpatrick and Ermentrout [8], where the authors conclude that upon breaking the translation symmetry of a neural field, by introducing spatially heterogeneous input or synapses, bumps in the stochastic neural field can become temporarily fixed to a finite number of locations.