# Dynamics of Periodic Waves in a Neural Field Model

Nikolai Bessonov, Anne Beuter, Sergei Trofimchuk, Vitaly Volpert

► To cite this version:

Nikolai Bessonov, Anne Beuter, Sergei Trofimchuk, Vitaly Volpert. Dynamics of Periodic Waves in a Neural Field Model. Mathematics, 2020, 8 (7), pp.1-17. 10.3390/math8071076. hal-03027337

HAL Id: hal-03027337

https://hal.science/hal-03027337v1

Submitted on 11 Mar 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
Article

# Dynamics of Periodic Waves in a Neural Field Model

Nikolai Bessonov<fnref n="1" />, Anne Beuter<fnref n="2" />,<fnref n="3" />, Sergei Trofimchuk<fnref n="4" /> and Vitaly Volpert<fnref n="5" />,<fnref n="6" />,<fnref n="7" />,*

Received: 12 March 2020; Accepted: 29 June 2020; Published: 2 July 2020

<figure><img src="image_2.png" /></figure>

check for
updates

**Abstract:** Periodic traveling waves are observed in various brain activities, including visual, motor, language, sleep, and so on. There are several neural field models describing periodic waves assuming nonlocal interaction, and possibly, inhibition, time delay or some other properties. In this work we study the influences of asymmetric connectivity functions and of time delay for symmetric connectivity functions on the emergence of periodic waves and their properties. Nonlinear wave dynamics are studied, including modulated and aperiodic waves. Multiplicity of waves for the same values of parameters is observed. External stimulation in order to restore wave propagation in a damaged tissue is discussed.

**Keywords:** neural field model; integro-differential equation; waves; brain stimulation

# 1. Introduction

## 1.1. Brain Activity and Periodic Travelling Waves

The brain displays a variety of highly nonlinear, complex dynamics across multiple spatial and temporal scales [1]. About 86 × 10⁹ neurons of the human brain entertain complex and fluctuating interactions. Understanding the dynamics of these interactions and the mechanisms underlying their control remains a technical and theoretical challenge. In other words, when healthy brain processes evolve toward abnormal and pathological states as a result of disease, degeneration or traumatic injury, how can a therapeutic intervention be used to reposition the control parameters and guide the dynamics back toward a healthy state? These brain processes are described today as interacting networks of nodes/hubs and edges which for the whole brain constitute the human connectome [2]. While the connectome is focused on anatomical connections, the dynamics of the networks are represented by functional connections (also called the dyname [1]). For example, brain functional connections described as a graph explore how signals are transmitted along neuroanatomical pathways and interact with local dynamics. Functional connections are often investigated via modeling [1]. One possibility is to use mathematical models to identify how an outside intervention such as neural electrical stimulation can modify local dynamics and how local dynamics will in turn affect other brain regions.

Cortical brain dynamics are investigated by means of periodic traveling waves (TW) characterized by their speed and frequency [3]. They describe the distribution of electric potential in the brain cortex. They are measured as a mean field potential (averaged macroscopic level). Propagating waves

<footnote n="1">¹ Institute of Problems of Mechanical Engineering, Russian Academy of Sciences, 199178 Saint Petersburg, Russia; nickbessonov1@gmail.com</footnote>

<footnote n="2">² Bordeaux INP, Avenue des Facultes, 33400 Talence, France; anne.beuter@wanadoo.fr</footnote>

<footnote n="3">³ CorStim SAS, 700 Avenue du Pic Saint Loup, 34090 Montpellier, France</footnote>

<footnote n="4">⁴ Instituto de Matematica y Fisica, Universidad de Talca, Casilla 747, Talca, Chile; trofimch@inst-mat.utalca.cl</footnote>

<footnote n="5">⁵ Institut Camille Jordan, UMR 5208 CNRS, University Lyon 1, 69622 Villeurbanne, France</footnote>

<footnote n="6">⁶ INRIA Team Dracula, INRIA Lyon La Doua, 69603 Villeurbanne, France</footnote>

<footnote n="7">⁷ Peoples' Friendship University of Russia (RUDN University), 6 Miklukho-Maklaya St, 117198 Moscow, Russia</footnote>

<footnote>* Correspondence: volpert@math.univ-lyon1.fr</footnote>
are observed during various types of brain activity. They provide subthreshold depolarization to individual neurons and increase their spiking probability. According to Muller et al. [3], TW “travel over spatial scales that range from the mesoscopic (single cortical areas and millimetres of cortex) to the macroscopic (global patterns of activity over several centimeters) and extend over temporal scales from tens to hundreds of milliseconds.” It has been proposed that TW mediate information transfer in the cortex.

Propagating waves increase the probability of neuron firing due to depolarization of neuronal membrane [4]. In [3] it was suggested that TW can be “spontaneously generated by recurrent circuits or evoked by external stimuli and travel along brain networks at multiple scales, transiently modulating spiking and excitability as they pass.” The phase relations between oscillations in different cortical regions produce the TW, and depending on the distance, axonal conduction delays can reach up to tens of milliseconds. These TW correlate with the subject’s performance, propagate in specific directions and synchronize distributed cortical networks that are communicating [5].

Botella-Soler et al. [6] identified for each subject the set of intracranial contacts that showed a larger percentage of detected events during slow wave cortical activity. They called these contacts “hubs” because the slow wave events in their travel through the cortical networks seemed to have a great probability of passing through the region close to the contact. Using probabilities, they were able to reconstruct a preferential propagation network for each subject. Slow waves have been reported to propagate across cortical areas at about 1m/s with multiple propagation paths and several points of origin. It seems that the slow waves have a preference to start in the prefrontal cortex and to end in posterior and temporal regions of the cortex. These waves appear to shape and strengthen neuronal networks.

Let us also note that time delays are intrinsic to the dynamics of brain networks, nodes and edges. Propagation speed along axons depends on axonal length and diameter. Conduction times along neural circuits also depend on degree of fiber myelination. The introduction of time delays in models can lead to significant changes in brain dynamics. They can be averaged or treated as distributed delays. In this paper we address the question of the effects of these delays in the dynamics of periodic cortical waves.

Thus, according to the biological observations, TW propagate in the cortex, activating and coordinating different parts of the brain. In this work, we study some of their properties. In the next section, we introduce the model. Then we present stability analysis which determines the conditions of wave appearance. Their nonlinear dynamics will be discussed in Section 3.

## 1.2. Neural Field Model

Neural field models were first introduced in [7]. Periodic traveling waves are described by several models (see [8–12] and Appendix A). In this work we consider one equation model with delay, and we discuss two mechanisms of the emergence of such waves, which were not sufficiently investigated previously. The first mechanism is related to the asymmetric connectivity functions [13], and the second one is determined by the delay in the response function. It is known that the loss of stability of the homogeneous in the space solution in this model does not lead to the bifurcation of periodic waves [8]. We show that they still appear for some larger values of time delay. We consider the one-dimensional neural field equation for the electric potential in the brain cortex written in the form

$$ \frac{\partial u}{\partial t} = D \frac{\partial^2 u}{\partial x^2} + W_a - W_i - \sigma u, \quad (1) $$

where D is the diffusion coefficient; $W_a$ and $W_i$ are given by the expressions

$$ W_a(x, t) = \int_{-\infty}^{\infty} \phi_a(x - y) S_a \left( u \left( y, t - \frac{|x - y|}{q_a} - \tau_a \right) \right) dy, \quad (2) $$
$$W_i(x, t) = \int_{-\infty}^{\infty} \phi_i(x - y) S_i \left( u \left( y, t - \frac{|x - y|}{q_i} - \tau_i \right) \right) dy \quad (3)$$

and characterize neuron activation and inhibition. These expressions describe the intensity of signal coming from all points y to the point x; $S_a(u)$ and $S_i(u)$ are smooth functions; $q_a$ and $q_i$ are the excitation speeds; $|x - y|/q_{a,i}$ is the time delay due to the excitation propagation from the point y to the point x; $\phi_a$ and $\phi_i$ are the connectivity functions,

$$\phi_a(r) = \begin{cases} a_1 e^{-b_1 r}, & r > 0 \\ a_3 e^{b_3 r}, & r < 0 \end{cases}, \quad \phi_i(r) = \begin{cases} a_2 e^{-b_2 r}, & r > 0 \\ a_4 e^{b_4 r}, & r < 0 \end{cases}, \quad (4)$$

where $a_i, b_i$ are some positive constants. Response functions $S_a(u)$ and $S_i(u)$ are non-negative, non-decreasing functions usually considered as sigmoid-type functions; $\tau_a$ and $\tau_i$ are time delays in neuron response to the activating and inhibitory signals. The last term on the right-hand side of Equation (1) describes signal decay with a decay rate $\sigma > 0$.

Neural field models are often considered without the diffusion term in the studies of both brain oscillations and traveling waves (see, e.g., [14–17]). In the case of small diffusion coefficients, its influence is not essential [18]. It can describe the ephaptic effect [19], ion diffusion and gap junction. The integral terms in Equation (1) characterize nonlocal neuron communication due to axons and dendrites. The speeds $q_a$ and $q_i$ of electric impulse propagation along axons are of the order 2–4 m/s, while the speed of wave propagation is one-two orders of magnitude less [20] (p. 213). Thus, we will consider a large speed limit $q_a = q_i \to \infty$:

$$\frac{\partial u}{\partial t} = D \frac{\partial^2 u}{\partial x^2} + \int_{-\infty}^{\infty} (\phi_a(x-y)S_a(u(y,t-\tau_a)) - \phi_i(x-y)S_i(u(y,t-\tau_i))) dy - \sigma u. \quad (5)$$

Up to the diffusion term, which is not very essential for small diffusion coefficients, this model is a particular case of the two equation model considered in [8]. If the kernels in the two equations are the same, then the system can be reduced to the single equation.

Different neural field models describe the propagation of periodic traveling waves (see Appendix A). In this work we will consider two other mechanisms of their emergence. One of them is related to asymmetric kernels $\phi_a$ and $\phi_i$, whose existence is confirmed by the experimental observations and used in theoretical considerations [13,21]. Another one is determined by the secondary bifurcation due to time delay for symmetric connectivity functions. Similarly to [8], we observe that the loss of stability of the homogeneous in space stationary solution leads to the appearance of periodic time oscillations independent of the space variable or of stationary periodic in space solutions. Periodic traveling waves bifurcate in the instability region, and they are unstable close to the bifurcation point. We will see that they can become stable under further change of parameters.

We studied the dynamics of the periodic waves, including their non-uniqueness for the same values of parameters. This property seems to us important, since the co-existence of waves with different frequencies is experimentally observed.

## 2. Stability

### 2.1. Linearization and Eigenvalues

In this section we will consider the Equation (5) on the interval $0 < x < L$ with periodic boundary conditions. We extend the function $u(x, t)$ by periodicity on the whole axis, $-\infty < x < \infty$, so that the integrals in Equation (5) are well defined. Let $u_0$ be a solution of the equation

$$\phi_a^* S_a(u) + \phi_i^* S_i(u) - \sigma u = 0,$$
where $\phi_a^* = \int_{-\infty}^{\infty} \phi_a(x)dx$, $\phi_i^* = \int_{-\infty}^{\infty} \phi_i(x)dx$. Then $u_0$ is a stationary solution of Equation (5). Linearizing this equation about $u_0$, we obtain the eigenvalue problem:

$$Dv'' + S'_{a}(u_0)e^{-\lambda\tau_a} \int_{-\infty}^{\infty} \phi_a(x-y)v(y)dy - S'_{i}(u_0)e^{-\lambda\tau_i} \int_{-\infty}^{\infty} \phi_i(x-y)v(y)dy - \sigma v = \lambda v. \quad (6)$$

Here, $v$ is a small perturbation of the stationary solution $u_0$. Applying the Fourier transform, we get

$$S'_{a}(u_0)e^{-\lambda\tau_a}\tilde{\phi}_a(\xi) - S'_{i}(u_0)e^{-\lambda\tau_i}\tilde{\phi}_i(\xi) - D\xi^2 - \sigma = \lambda, \quad (7)$$

where

$$\tilde{\phi}_a(\xi) = \frac{a_1 b_1}{b_1^2 + \xi^2} + \frac{a_3 b_3}{b_3^2 + \xi^2} + i\xi \left( \frac{a_1}{b_1^2 + \xi^2} - \frac{a_3}{b_3^2 + \xi^2} \right), \\ \tilde{\phi}_i(\xi) = \frac{a_2 b_2}{b_2^2 + \xi^2} + \frac{a_4 b_4}{b_4^2 + \xi^2} + i\xi \left( \frac{a_2}{b_2^2 + \xi^2} - \frac{a_4}{b_4^2 + \xi^2} \right)$$

are Fourier transforms of the functions $\phi_a$ and $\phi_i$, respectively. Let us note that Fourier transform of non-integrable functions is considered in the sense of generalized functions. Using Fourier series instead of Fourier transform allows one to use a classical function instead of generalized functions. However, the advantage of the Fourier transform is that the periodicity of solutions is not imposed.

Set $\lambda = iv$. Separating the real and imaginary parts in Equation (7), we obtain:

$$S'_{a}(u_0) \cos(\nu\tau_a) \left( \frac{a_1 b_1}{b_1^2 + \xi^2} + \frac{a_3 b_3}{b_3^2 + \xi^2} \right) + S'_{a}(u_0) \xi \sin(\nu\tau_a) \left( \frac{a_1}{b_1^2 + \xi^2} - \frac{a_3}{b_3^2 + \xi^2} \right) - S'_{i}(u_0) \cos(\nu\tau_i) \left( \frac{a_2 b_2}{b_2^2 + \xi^2} + \frac{a_4 b_4}{b_4^2 + \xi^2} \right) - S'_{i}(u_0) \xi \sin(\nu\tau_i) \left( \frac{a_2}{b_2^2 + \xi^2} - \frac{a_4}{b_4^2 + \xi^2} \right) - D\xi^2 = \sigma, \quad (8)$$

$$-S'_{a}(u_0) \sin(\nu\tau_a) \left( \frac{a_1 b_1}{b_1^2 + \xi^2} + \frac{a_3 b_3}{b_3^2 + \xi^2} \right) + S'_{a}(u_0) \xi \cos(\nu\tau_a) \left( \frac{a_1}{b_1^2 + \xi^2} - \frac{a_3}{b_3^2 + \xi^2} \right) + S'_{i}(u_0) \sin(\nu\tau_i) \left( \frac{a_2 b_2}{b_2^2 + \xi^2} + \frac{a_4 b_4}{b_4^2 + \xi^2} \right) - S'_{i}(u_0) \xi \cos(\nu\tau_i) \left( \frac{a_2}{b_2^2 + \xi^2} - \frac{a_4}{b_4^2 + \xi^2} \right) = v. \quad (9)$$

## 2.2. Symmetric Connectivity Functions with Time Delay

If the symmetry condition

$$a_1 = a_3, a_2 = a_4, b_1 = b_3, b_2 = b_4 \quad (10)$$

is satisfied, then Equations (8) and (9) are as follows:

$$S'_{a}(u_0) \cos(\nu\tau_a) \frac{a_1 b_1}{b_1^2 + \xi^2} - S'_{i}(u_0) \cos(\nu\tau_i) \frac{a_2 b_2}{b_2^2 + \xi^2} - D\xi^2/2 = \sigma/2, \quad (11)$$

$$-S'_{a}(u_0) \sin(\nu\tau_a) \frac{a_1 b_1}{b_1^2 + \xi^2} + S'_{i}(u_0) \sin(\nu\tau_i) \frac{a_2 b_2}{b_2^2 + \xi^2} = v/2. \quad (12)$$

We can express $\xi^2$ through $v$ from the last equation and substitute them into Equation (11). The resulting equation with respect to $v$ should be solved numerically or asymptotically. Since the calculations are sufficiently complex, we will consider here a simplified case where $\tau_a = 0$. Numerical simulations show that periodic waves can exist in this case (Section 3). Our aim here is to analyze their bifurcations from the homogeneous in space solution. Assuming that $v \ge 0$, we also get the conjugate solution $-v$.
If $\tau_a = 0$, then Equations (11) and (12) are as follows:

$$S'_{a}(u_{0}) \frac{a_{1}b_{1}}{b_{1}^{2} + \xi^{2}} - S'_{i}(u_{0}) \cos(v\tau_{i}) \frac{a_{2}b_{2}}{b_{2}^{2} + \xi^{2}} - D\xi^{2}/2 = \sigma/2, \quad (13)$$

$$S'_{i}(u_{0}) \sin(v\tau_{i}) \frac{a_{2}b_{2}}{b_{2}^{2} + \xi^{2}} = v/2. \quad (14)$$

By the stability boundary we will understand here the value $\sigma = \sigma_0$ such that system (13), (14) does not have solutions for $\sigma > \sigma_0$, and it has at least one solution for $\sigma < \sigma_0$. This stability boundary $\sigma_0$ is uniquely defined and it can be negative.

**Proposition 1.** *At the stability boundary, $\xi = 0$ or $v = 0$.*

**Proof.** Consider Equation (14) as an equation with respect to $v$ for a fixed $\xi$. It has a solution $v \neq 0$ if and only if

$$2\tau_i S'_i(u_0) \frac{a_2 b_2}{b_2^2 + \xi^2} > 1. \quad (15)$$

Suppose that this condition is satisfied for some $\xi_0 > 0$, and denote by $v(\xi)$ solution of Equation (14) in the vicinity of $\xi = \xi_0 > 0$, and $v_0 = v(\xi_0) > 0$. Assume that the corresponding value $\sigma_0$ belongs to the stability boundary.

Taking into account Equation (14), we can write Equation (13) as follows:

$$2S'_{a}(u_{0}) \frac{a_{1}b_{1}}{b_{1}^{2} + \xi^{2}} - v \cot(v\tau_{i}) - D\xi^{2} = \sigma. \quad (16)$$

Let $\xi_1$ be sufficiently close to $\xi_0$ such that condition (15) remains satisfied, and $\xi_1 < \xi_0$. Considering $v(\xi_0) < \pi/\tau_i$, we get $v(\xi_1) > v(\xi_0)$. Therefore,

$$v(\xi_1) \cot(v(\xi_1)\tau_i) < v(\xi_0) \cot(v(\xi_0)\tau_i).$$

Set

$$F(\xi) = 2S'_{a}(u_{0}) \frac{a_{1}b_{1}}{b_{1}^{2} + \xi^{2}} - v(\xi) \cot(v(\xi)\tau_{i}) - D\xi^{2}.$$

Thus, $F(\xi_1) > F(\xi_0)$. Hence, $\sigma_0 = F(\xi_0)$ cannot belong to the stability boundary since $\sigma_1 = F(\xi_1) > \sigma_0$, and system (13), (14) has a solution for $\sigma = \sigma_1$. This contradiction shows that $\xi_0 = 0$ when $v_0 \neq 0$.

Let $2\pi/\tau_i < v(\xi_0) < 5\pi/(2\tau_i)$. Then there is another solution $v_1(\xi_0)$ of Equation (14) such that $\pi/(2\tau_i) < v_1(\xi_0) < \pi/\tau_i$. Therefore, $\cos(v_1(\xi_0)\tau_i) < 0$, $\cos(v(\xi_0)\tau_i) > 0$, and

$$\sigma_0 = F(\xi_0) < 2S'_{a}(u_0) \frac{a_1 b_1}{b_1^2 + \xi_0^2} - v_1(\xi_0) \cot(v_1(\xi_0)\tau_i) - D\xi_0^2.$$

Hence, $\sigma_0$ cannot belong to the stability boundary.

Finally, in all other cases with $v(\xi_0) > 5\pi/(2\tau_i)$ we obtain a contradiction similarly to the two cases considered before. This contradiction shows that at the stability boundary, either $\xi = 0$ or $v = 0$. $\square$

Let us now determine the stability boundary with respect to time delay for a fixed $\sigma$. If $v = 0$, then from Equation (13) we get

$$\sigma = \frac{a_1}{b_1^2 + \xi^2} - \frac{a_2}{b_2^2 + \xi^2} - D\xi^2, \quad (17)$$
where $\alpha_1 = 2S'_a(u_0)a_1b_1$, $\alpha_2 = 2S'_i(u_0)a_2b_2$. Assuming that $v \neq 0$, we set $\beta_i = \alpha_i/(b_i^2 + \xi^2)$, $i = 1,2$. Then $\sin(v\tau_i) = v/\beta_2$,

$$ \beta_1 - \sqrt{\beta_2^2 - v^2} = \sigma, \quad v^2 = \beta_2^2 - (\beta_1 - \sigma)^2. $$

**Example 1.** Consider the values of parameters: $S'_a(u_0) = S'_i(u_0) = 20, a_1 = a_2 = 4, b_1 = 40, b_2 = 20, \alpha_1 = 6400, \alpha_2 = 3200, \sigma = 0.01$, Then we find

$$ \xi = 0, \beta_1 = 4, \beta_2 = 8, v^2 \approx 48, v \approx 6.93, \quad \tau_i = \frac{1}{v} \arcsin \frac{v}{\beta_2} = 0.151, $$

and

$$ \xi = \pi, \beta_1 = 3.975, \beta_2 = 7.805, v^2 \approx 45.123, v \approx 6.717, \quad \tau_i = \frac{1}{v} \arcsin \frac{v}{\beta_2} = 0.154 $$

Hence, the homogeneous in space oscillations appear for a smaller value of time delay than periodic traveling waves. Therefore, these waves are unstable in the vicinity of the bifurcation point. We will see in the next section that they become stable for larger values of time delay.

## 2.3. Asymmetric Connectivity Function without Time Delay

If conditions (10) do not hold, then eigenvalue (7) has a nonzero imaginary part, $\lambda(\xi) = \alpha(\xi) + iv(\xi)$,

$$ \alpha(\xi) = \frac{s_a a_1 b_1}{b_1^2 + \xi^2} + \frac{s_a a_3 b_3}{b_3^2 + \xi^2} - \frac{s_i a_2 b_2}{b_2^2 + \xi^2} - \frac{s_i a_4 b_4}{b_4^2 + \xi^2} - D\xi^2 - \sigma, $$

$$ v(\xi) = s_a \xi \left( \frac{a_1}{b_1^2 + \xi^2} - \frac{a_3}{b_3^2 + \xi^2} \right) - s_i \xi \left( \frac{a_2}{b_2^2 + \xi^2} - \frac{a_4}{b_4^2 + \xi^2} \right), \quad (18) $$

$s_a = S'_a(u_0), s_i = S'_i(u_0)$. Condition $\alpha(\xi) = 0$ determines the stability boundary,

$$ \sigma = \frac{s_a a_1 b_1}{b_1^2 + \xi^2} + \frac{s_a a_3 b_3}{b_3^2 + \xi^2} - \frac{s_i a_2 b_2}{b_2^2 + \xi^2} - \frac{s_i a_4 b_4}{b_4^2 + \xi^2} - D\xi^2 \equiv \Phi_2(\xi). \quad (19) $$

If $\sigma < \Phi_2(\xi)$ for some values of $\xi$, then the solution loses its stability due to a pair of complex conjugate eigenvalues $\lambda(\xi) = \alpha(\xi) \pm iv(\xi)$. For $\alpha = 0$ corresponding to the stability boundary, the bounded solution of the linearized equation can be written as follows:

$$ u(x,t) = e^{ivt} e^{i\xi x} + e^{-ivt} e^{-i\xi x} = \cos(vt + \xi x). $$

We find $\xi$ from equality $\alpha(\xi) = 0$ and $v$ from (18). The frequency of the periodic wave equals $\xi$ and its speed $c = -v/\xi$ can be determined from (18):

$$ c = -s_a \left( \frac{a_1}{b_1^2 + \xi^2} - \frac{a_3}{b_3^2 + \xi^2} \right) + s_i \left( \frac{a_2}{b_2^2 + \xi^2} - \frac{a_4}{b_4^2 + \xi^2} \right). \quad (20) $$

Different waves with the frequencies satisfying condition $\sigma < \Phi(\xi)$ can exist for the same values of parameters. Their speed can be increasing or decreasing functions of frequency according to (20).

## 3. Numerical Results

Numerical simulations of Equation (5) in a bounded interval $0 \le x \le L$ with periodic boundary conditions will be started in the case without time delay and continued with the case of time delay in the response functions $S_a(u) = S_b(u) = \arctan(hu)$, $h > 0$. We will finish this section with modeling of the stimulation of the damaged tissue in order to restore wave propagation.
The numerical method uses an implicit finite difference scheme with Thomas's algorithm for the inversion of the tridiagonal matrix. Initial conditions are considered either in the form of a piece-wise constant function or a sinus-like function with a given periodicity. The continuation method is used to follow the branches of solutions. In this case, the value of some parameter gradually changes. As an initial condition for the new value of parameter, we take the result of the simulation obtained for the previous value of parameter. The time step was usually taken as 0.05 and the space step 0.005. The accuracy of numerical simulations was controlled by decreasing the time and space steps. The typical values of parameters used in numerical simulations are presented in Table A1. Their specific values are given below.

### 3.1. Wave Propagation without Time Delay

In the case without time delay ($\tau_a = \tau_i = 0$) and with symmetric connectivity functions $\phi_a$ and $\phi_i$, where $a_1 = a_3, a_2 = a_4, b_1 = b_3, b_2 = b_4$, periodic in space stationary solutions bifurcate from the constant solution. Linear stability analysis shows that traveling waves with nonzero speed do not bifurcate in this case.

If the connectivity functions are not symmetric, then traveling waves with nonzero speed are observed (see Figures 1 and 2). Linear stability analysis allows the determination of the wave speed and frequency near the bifurcation point. Let us consider an example with the following values of parameters: $a_1 = 3, a_2 = 4, a_3 = 3, a_4 = 1, b_1 = 40, b_2 = 20, b_3 = 40, b_4 = 20, D = 10^{-3}, h = 20, L = 2$. Analysis of the function $\alpha(\xi)$ (Section 2.3) shows that it has a maximum at $\xi = 20.53$. It is positive for $\sigma < \sigma_c = 0.4$ resulting in the bifurcation of the periodic wave with the speed given by Equation (20). For the chosen values of parameters, we get $c = -0.54$. Numerical simulations show the emergence of a periodic wave with 7 periods in the interval $[0, L]$ and with the corresponding frequency $\xi = 21.98$. Since the number of periods is an integer, the frequency of the wave observed in numerical simulations is not precisely equal to the analytical value. The wave with the closest frequency emerges since the corresponding eigenvalue has the maximal real part in comparison with the waves with other periods. The speed of this wave $c = -0.52$ is close to the analytical value.

For the values of parameters far from the stability boundary, two types of solutions are observed in numerical simulations: periodic waves with a constant speed (Figure 1) and aperiodic waves with oscillating speed (Figure 2). In the first case, the wave has conventional form $w(x - ct)$, where $w(x)$ is a periodic in space function and $c$ is a constant. In the second case, the amplitude of spatial peaks and the wave speed oscillate.

<figure><img src="image_3.png" /><figcaption>**Figure 1.** A snapshot of periodic wave described by Equation (5) for $\tau_a = \tau_i = 0$ (left) and a zoom-in on the lower part of the graph (right). The values of parameters are as follows: $D = 10^{-4}$, $\sigma = 0.01$, $L = 2$, $a_1 = a_2 = 0.6$, $a_3 = a_4 = 4$, $b_1 = b_3 = 40$, $b_2 = b_4 = 20$. Here $S_a(u) = S_i(u) = \arctan(hu)$, and $h = 20$. Here and in all examples below, the stationary solution $u_0 = 0$.</figcaption></figure>
<figure><img src="image_4.png" /><figcaption>**Figure 2.** A snapshot of aperiodic wave described by Equation (5) for $\tau_a = \tau_i = 0$ (left) and a zoom-in on the lower part of the graph (right). Periodic (Figure 1) and aperiodic waves can exist for the same values of parameters. Convergence of solutions to one of them is determined by the initial conditions. The values of parameters are as follows: $D = 10^{-4}$, $\sigma = 0.01$, $L = 2$, $a_1 = a_2 = 0.6$, $a_3 = a_4 = 4$, $b_1 = b_3 = 40$, $b_2 = b_4 = 20$. Here $S_a(u) = S_i(u) = \arctan(hu)$, and $h = 20$.</figcaption></figure>

Different periodic regimes can co-exist for the same values of parameters. Figure 3 shows the dependence of the wave speed on the values $a_1(= a_2)$ for all other parameters fixed. There are three branches of solutions corresponding to different spatial frequencies. Thus, there are three different periodic waves for the same values of parameters with the speeds depending on their frequency.

<figure><img src="image_5.png" /><figcaption>**Figure 3.** The speed of periodic waves for different values of the parameters $a_1 = a_2$. Three curves correspond to different values of the spatial frequency: (1) solution with 12 periods in the interval $[0, L]$, (2) 13 periods, (3) 14 periods. Branching in curve 1 shows the transition to modulated oscillations with the maximal and minimal values of the oscillating speed. The values of parameters are as follows: $D = 10^{-4}$, $\sigma = 0.01$, $L = 2$, $a_3 = a_4 = 4$, $b_1 = b_3 = 40$, $b_2 = b_4 = 20$. Here $S_a(u) = S_i(u) = \arctan(hu)$, and $h = 20$.</figcaption></figure>

For a sufficiently small $a_1$, transition to modulated waves can occur (curve 1) with the amplitude and speed depending on time. Such solutions can be qualitatively approximated by the function $u(x, t) = (k_1 + \epsilon \sin(k_2x + k_3t)) \sin(k_4x + k_5t)$, where $k_i$ are some constants. For even smaller values of
$a_1$, transition to aperiodic oscillations is observed for all three branches of solutions. These oscillations can coexist with periodic or modulated periodic waves for the same values of parameters.

## 3.2. Time Delay and Symmetric Connectivity Functions

### 3.2.1. Initial Conditions

As it was discussed in the previous section, the loss of stability of the homogeneous in space stationary solution leads either to the homogeneous in space time oscillations or to the stationary periodic in space solutions. Consider time delay as a bifurcation parameter. If it exceeds a critical value, then time oscillations emerge. Periodic traveling waves bifurcate for a larger value of time delay, and they are unstable in the vicinity of the bifurcation point. Numerical simulations show that they can become stable under a further increase of $\tau$. Nevertheless, periodic time oscillations homogeneous in space remain stable. Therefore, we need to choose some particular initial conditions in order to get periodic traveling waves.

The simulations presented in this section were carried out with two types of initial conditions. In the first case, we consider the equation

$$ \frac{\partial u}{\partial t} = D \frac{\partial^2 u}{\partial x^2} + I(x, t) \qquad (21) $$

on the time interval $0 \le t \le T_0$. Here $I(x,t) = I_0 \cos(px + qt)$, and the initial condition $u(x,0) = 0$. The result of this simulation is considered as the initial condition for Equation (5). We set $I_0 = 0.5, q = 0.015, T_0 = 20$; the value of $p$ is taken 3, 6, 9 depending on the required space periodicity.

The second type of initial conditions is used in the continuation method. The results of the simulations of Equation (5) for some values of parameters are used as initial conditions for some other values of parameters.

### 3.2.2. Multiplicity of Waves and Parameter Dependence

An example of periodic traveling waves with different initial conditions and the same values of parameters are shown in Figure 4. We observe one, two and three-period waves generated by function $I(x,t)$ for $p = 3, 6$ and $9$. These waves have different speeds and amplitudes. The waves with larger wavelengths have higher speeds and amplitudes. If time delay is sufficiently small, then the waves become unstable, and the transition to periodic time oscillations independent of the space variable is observed. If the value $I_0$ is small enough, then the solution converges to the homogeneous time oscillations.

<figure><img src="image_6.png" /><figcaption>**Figure 4.** Three types of periodic waves observed for the same values of parameters and having different spatial frequency and speed: one-period wave with speed $-0.027$ (left); two-period wave with speed $-0.012$ (middle); three-period wave with speed $-0.0094$. The values of parameters are as follows: $D = 10^{-4}$, $\sigma = 0.01$, $L = 2$, $\tau_a = 0$, $\tau_i = 12$, $a_1 = a_2 = a_3 = a_4 = 4$, $b_1 = b_3 = 40$, $b_2 = b_4 = 20$. Here $S_a(u) = S_i(u) = \arctan(hu)$, and $h = 20$.</figcaption></figure>
Increase of time delay leads to the increase of the wave amplitude and to the decrease of its speed (Figure 5). Wave propagation is determined by the transmission of activating signal between the neurons. Since time delay retards this transmission, the speed decreases. This result was previously obtained analytically for monotone waves [22]. Since the activating signal propagates slower with time delay, the local response and the wave amplitude increase.

<figure><img src="image_7.png" /><figcaption>**Figure 5.** The absolute value of the speed (1a, 2a) and the amplitude (1b, 2b) for the one-period wave (1a, 1b) and two-period wave (2a, 2b) depending on $\tau_i$. Connected points correspond to stable solutions; separate points to unstable solutions. The latter lead to the appearance of stationary solutions periodic in space. The values of parameters are as follows: $D = 10^{-4}$, $\sigma = 0.01$, $\tau_a = 0$, $a_1 = a_2 = a_3 = a_4 = 4$, $b_1 = b_3 = 40$, $b_2 = b_4 = 20$. Here $S_a(u) = S_i(u) = \arctan(hu)$, and $h = 20$.</figcaption></figure>

Let us note that the critical value of time delay $\tau_i \approx 1.5$ found in numerical simulations and the value of the wave speed near the bifurcation point $c \approx -2$ correspond to the analytical values determined in the example in Section 2.2 (Figure 5, curves 1a, 1b). The analytical value of the speed $c = -\nu/\zeta = -2.13$ at the bifurcation point is slightly different from the numerical value. Since the periodic wave is unstable near the bifurcation point, we can determine its speed in numerical simulations only at some distance from the critical value $\tau_i$.

Figure 6 shows the dependence of periodic waves on the value $b_2$ (equal to $b_4$). The wave amplitude and speed decrease with the increase of $b_2$. There is a critical value $b_2 \approx 40$ for which the speed becomes zero, and a transition to another branch of solutions is observed. These are stationary, periodic in space solutions with growing amplitude as $b_2$ increases. It is interesting to note the existence of weakly oscillating time periodic solutions in a narrow interval between traveling waves and stationary solutions.

The results presented above were obtained for a single delay $\tau_i$ in the inhibition term while $\tau_a = 0$. If we fix $\tau_i = 1$ and increase $\tau_a$ beginning from $\tau_a = 0$, then the amplitude and the speed of traveling waves are not monotonous. The former first decreases, passes through the minimum and then increases, while the latter increases in the beginning and decreases for larger values of the delay (not shown).
<figure><img src="image_8.png" /><figcaption>**Figure 6.** Dependence of the amplitude of the observed regimes on $b_2(= b_4)$. Periodic traveling waves are located on branch 1; stationary solutions on branch 3. A short branch 2 contains weakly oscillating solutions. The values of parameters are as follows: $D = 10^{-4}$, $\sigma = 0.01$, $L = 2$, $\tau_a = 0$, $\tau_i = 1$, $a_1 = a_2 = a_3 = a_4 = 4$, $b_1 = b_3 = 40$. Here $S_a(u) = S_i(u) = \arctan(hu)$, and $h = 20$.</figcaption></figure>

### 3.3. Stimulation

Wave propagation can be different in a damaged tissue [23]. In this section we discuss its possible restoration by external stimulation. This question was considered in [24] for monotonous wave fronts.

#### 3.3.1. Exact Solution of the Stimulation Problem

Let us write Equation (5) for the normal tissue

$$ \frac{\partial u}{\partial t} = D \frac{\partial^2 u}{\partial x^2} + J(u) - \sigma u, \qquad (22) $$

where $J(u) = \int_{-\infty}^{\infty} \phi_a(x-y)S_a(u(y,t-\tau_a)) dy - \int_{-\infty}^{\infty} \phi_i(x-y)S_i(u(y,t-\tau_i)) dy$,

and a similar equation for the damaged tissue

$$ \frac{\partial v}{\partial t} = D \frac{\partial^2 v}{\partial x^2} + J^*(v) - \sigma v, \qquad (23) $$

where $J^*(v) = \int_{-\infty}^{\infty} \phi_a^*(x-y)S_a^*(v(y,t-\tau_a)) dy - \int_{-\infty}^{\infty} \phi_i^*(x-y)S_i^*(v(y,t-\tau_i)) dy$,

and * denotes the functions for the damaged tissue. In the presence of stimulation $I(x,t)$ this equation is as follows:

$$ \frac{\partial z}{\partial t} = D \frac{\partial^2 z}{\partial x^2} + J^*(z) - \sigma z + I(x,t). \qquad (24) $$

We will choose this function such that solution $z(x,t)$ of Equation (24) becomes equal solution $u(x,t)$ of Equation (22). We set $z(x,t) = u(x,t)$, and substitute $u(x,t)$ in Equation (24). Then

$$ I(x,t) = \frac{\partial u}{\partial t} - D \frac{\partial^2 u}{\partial x^2} - J^*(u) + \sigma u = J(u) - J^*(u). $$

Thus, the stimulation function

$$ I(x,t) = J(u) - J^*(u) $$
gives solution of the complete reconstruction problem.

### 3.3.2. Approximate Solution of the Stimulation Problem

The stimulation function suggested above uses the solution $u(x, t)$ which may not be known for the patient with the damaged brain tissue. In this case, some approximate solutions of the stimulation problem can be used. Consider the integral $J^*(v)$ in the following form:

$$J^*(v) = \int_{-\infty}^{\infty} W(x)W(y) (\phi_a(x-y)S_a(v(y,t-\tau_a)) - \phi_i(x-y)S_i(v(y,t-\tau_i))) dy,$$

where $W(x) = w_0 < 1$ for $x_1 \le x \le x_2$, $x_1, x_2 \in (0, L)$, and $W(x) = 1$ outside the interval $[x_1, x_2]$. Hence, the connectivity function decreases if $x$ or $y$ belongs to the damaged area $[x_1, x_2]$. If we set $w_0 = 0$, then the connectivity function vanishes if one of the two points $x$ or $y$ (neurons) belongs to the damaged area.

Without damage, the periodic traveling wave solution of Equation (22) and the integral $J(u)$ can be approximated by a cosine function (Figure 7). Hence, we will look for the stimulation function in the form approximating the periodic wave: $I(x, t) = I_0(x) \cos(px + qt)$. We set $p = 6, q = 1$ to approximate the frequency and the speed of the wave in the normal tissue. We set $I_0(x) = i_0$ for $x_1 \le x \le x_2$, and $I_0(x) = i_1$ outside the interval $[x_1, x_2]$. The results of numerical simulations are presented in Figure 7 with the comparison of the normal tissue, damaged tissue without stimulation and damaged tissue with stimulation. We observe a periodic wave in the normal tissue. The behavior of the solution is completely different in the damaged tissue. The solution is close to 0 at the damaged interval (green in the middle figure), and it oscillates periodically in time from both sides of this interval. Stimulation restores the wave propagation with the same frequency and speed. The choice of the stimulation amplitude $i_0 = 0.6$ and $i_1 = 0.1$ (for the example in the figure) is important. If we set $i_0 = 0.6, i_1 = 0$, the stimulation is not successful; there is no wave propagation. Thus, stimulation should be done not only inside the damaged interval but also around it.

<figure><img src="image_9.png" /><figcaption>**Figure 7.** Periodic traveling wave in the normal tissue and the integral $J(u)$ (left). Snapshots of solutions in the damaged tissue for two different moments of time (a) and (b) (middle). Solution with stimulation becomes close to the periodic wave (right). The values of parameters are as follows: $D = 10^{-4}$, $\sigma = 0.01$, $L = 2$, $\tau_a = 0$, $\tau_i = 1$, $a_1 = a_2 = a_3 = a_4 = 4$, $b_1 = b_3 = 40$, $b_2 = b_4 = 20$, $S_a(u) = S_i(u) = \arctan(hu)$, $h = 20$, $p = 6$, $q = 1$, $i_0 = 0.6$, $i_1 = 0.1$. Green interval $[0.5, 1.07]$ shows the damaged tissue, $w_0 = 0$.</figcaption></figure>
## 4. Discussion

Brain functioning is determined by large scale networks of epicenters (hubs) located in the cortex and connected by white matter fiber tracks. The structure of these networks depends on the particular type of the brain activity (motor, language and so on) and on the inter-individual variation. The epicenters exchange information by means of signaling along the cortex or along the white matter fibers. Apparently, this signaling occurs in the form of traveling waves. Periodic traveling waves are observed in thalamus, visual cortex, hippocampus and other parts of the brain. There is more and more evidence that they play the key roles in brain functioning. However, their exact roles and organization are not known. Their characteristics, such as speed, frequency and amplitude, can vary in wide limits (Appendix B), and it is not clear how they are initiated and stop, and how they are related to some particular types of activities. For example, periodic waves are observed in the beginning of produced speech, and they inverse their direction at the end [25]. The corresponding mechanisms governing these activities are not yet understood. These investigations are at the stage of accumulation of biological data and of elaboration of different models whose role and utility will become clearer with time.

Neural field models are widely used to study brain patterns, including stationary structures, pulses, wave fronts, and periodic waves. There are a number of models which describe periodic waves (Appendix A) since it is relatively easy to find these waves on the basis of linear stability analysis of the homogeneous in space stationary solution. The mechanisms of this instability can be related to a combination of inhibition, time delay, refractoriness (also a variant of time delay) and nonlocal interaction. In this work we study the influence of asymmetry of connectivity functions and of time delay in neuron response.

Neuron connectivity is provided by numerous axons whose density decrease approximately exponentially as a function of distance [26,27]. In modeling, the connectivity functions are usually considered to be symmetric; that is, axon connections between points x and y have the same density as between y and x. There is some evidence that connectivity can be asymmetric [13].

In this work we study periodic waves emerging due to the asymmetry of the connectivity function and due to time delay in the response function. Both of them correspond to the existing biological mechanisms. They allow us to study periodic waves in the minimal model which consists of a single integro-differential equation. Other models contain two equations (Appendix A). In the case of asymmetric connectivity function, periodic waves bifurcate due to the loss of stability of the homogeneous in space stationary solution. They become unstable, while bifurcating periodic waves are stable. This is different in the case of time delay. Periodic waves bifurcate from the unstable stationary solution because of space-independent oscillations. Increasing time delay leads to their stabilization.

Though the model contains quite many parameters, using some combination of parameters reduces their number. Next, there are some additional relations between the parameters which determine the stability boundary. Furthermore, there are some physical estimates, such as the decay rate of neuron connectivity function or the value of time delay. Finally, the wave speed and frequency, which depend on parameters, should be in some experimentally observed range. Altogether, these constraints determine some limited intervals of parameter variation (Appendix B).

### Cortex Damage and Stimulation

The network of hubs related to some brain function (connectome) can be damaged because of stroke or other factors leading to partial or complete loss of the corresponding function. Post-stroke patient recovery is often incomplete, and usually limited to six months after the accident. Various stimulation techniques are discussed in the literature, but their results are controversial [28].

Post-stroke brain damage can influence propagation of brain waves between the hubs of the connectome. Traveling waves in the cortex can have several functions, including activation of some of its parts. This activation facilitates firing of individual neurons [4]. TW reflect information originating
from right and left hemispheres, traveling short and long distances, and containing various time delays. In other words, TW coordinate multiple faster and slower speech events by preparing the arrival of signals traveling along white matter tracts to specific hubs. In particular, direct brain recording (ECOG) showed the presence of TW during consonant-vowel syllables' pronunciation [25]. Summarizing the results by Gross et al. [29] examining how brain waves help us make sense of speech in healthy subjects, Weaver [30] indicated the presence of different time scales in speech production from tens of milliseconds (phoneme) to hundreds of milliseconds (intonation). Today, the language networks have been identified with precision, including phonological processing, speech planning, language semantics, spatial cognition and other functions [31].

We address in this work the question about possible restoration of cortex waves by external stimulation. We show that appropriate choice of injected current allows the recovery of wave speed and frequency, possibly leading to a better communication between the hubs of the connectome. This proof of concept is the first step on the long road to a possible application of this approach to patient rehabilitation.

**Author Contributions:** Conceptualization, V.V.; Data curation, A.B.; Formal analysis, S.T.; Methodology, A.B. and V.V.; Software, N.B.; Writing original draft, V.V. All authors have read and agreed to the published version of the manuscript.

**Funding:** A. Beuter was funded by CorStim SAS. S. Trofimchuk was partially supported by FONDECYT (Chile), project 1190712. This research was also funded by "RUDN University Program 5-100" and the French-Russian project PRC2307.

**Conflicts of Interest:** The authors declare no conflict of interest.

## Appendix A. Periodic Waves in Different Models

### Two equations with time delay.

In the work [8] a single neuron population model (excitation)

$$ \tau \frac{\partial u}{\partial t} = \int_{-\infty}^{\infty} P(x-y) \psi(u(y, t-d)) du - u \quad (A1) $$

and a two population model (excitation and inhibition)

$$ \tau \frac{\partial u}{\partial t} = \int_{-\infty}^{\infty} (P_{11}(x-y)\psi_1(u(y, t-d)) - P_{12}(x-y)\psi_2(v(y, t-d))) dy - u \quad (A2) $$

$$ \tau \frac{\partial v}{\partial t} = \int_{-\infty}^{\infty} (P_{21}(x-y)\psi_1(u(y, t-d)) - P_{22}(x-y)\psi_2(v(y, t-d))) dy - v \quad (A3) $$

are considered. Here $P(x)$ and $P_{ij}(x)$ are symmetric positive functions. A particular example of step-wise constant functions is studied. Periodic traveling waves cannot exist for the first model. They are observed for the second model. The existence of pulse solutions in a similar model without delay was studied in [32].

### One equation with distributed speed and delay.

Equation with a distributed propagation speed and time delay is considered in [9]:

$$ \begin{aligned} L \left( \frac{\partial u}{\partial t} \right) = & \alpha \int_{0}^{\infty} g(v) \int_{-\infty}^{\infty} K(z) S(u(x+z, t - |z|/v)) dz dv + \\ & \beta \int_{0}^{\infty} f(\tau) \int_{-\infty}^{\infty} F(z) S(u(x+z, t - \tau)) dz d\tau, \end{aligned} \quad (A4) $$
where L is a second-order differential operator; the functions K(z) and F(z) include both activatory and inhibitory kernels. Different regimes are observed: periodic in time and independent of space, stationary periodic in space, periodic traveling waves.

The model with distributed time delay is considered in [10]

$$u(x,t) = \int_{-\infty}^{\infty} w(x-y)dy \int_{-\infty}^{t} \eta(t-s) f(u(y,s - |x-y|/v))ds.$$

The same regimes as above and oscillating Turing structures are observed.

## Neural field model with linear adaptation.

It is a two equation model

$$\tau \frac{\partial u}{\partial t} = -u - \beta v + \int_{D} w(x-y) F(u(y,t)) dy, \quad (A5)$$

$$\frac{1}{\alpha} \frac{\partial v}{\partial t} = u - v, \quad (A6)$$

where the second variable represents a linear adaptation (see [11] and the references therein). A large variety of waves and patterns are observed, including stationary periodic in space solutions, traveling waves, modulated traveling waves and stationary and oscillating bumps.

## Neural field model with refractoriness.

Periodic traveling waves are also found in one-equation model without inhibition term but with neuron refractoriness (time delay after firing) [12]:

$$\frac{1}{r} \frac{\partial u}{\partial t} = -u + \left(1 - \int_{t-1}^{t} u(x,s) ds\right) f(w \otimes u).$$

Here ⊗ denotes spatial convolution.

## Appendix B. The Values of Parameters

Periodic traveling waves behave as cos(βt + ζx) with the time frequency β, space frequency ζ and speed c = -β/ζ. With the interval length L = 2 cm, we get ζ = 2π/L ≈ 3 cm⁻¹.

Let us consider the example of the simulations in Figure 5 with τ = 1 and wave speed 0.3. If this value of τ corresponds to the characteristic time delay 10 ms ([9,33] Chapter 2.3.1), that is, the time unit in the simulation corresponds to 0.01 s, then c = 0.3 × 100 = 30 cm/s.

The time frequency β = cζ = 90 s⁻¹ belongs to the upper limit of the observed range. The value of the wave speed, and respectively, the time frequency linearly dependent on it can be decreased by the variation of parameters aᵢ and bᵢ.

Connectivity functions can be estimated from the data in [26,27]. It exponentially decreases with the rate of decrease in the interval 3–10 times at the at the distance 0.03 cm. This corresponds to the exponential exp(-μx) with μ in the range 30 ÷ 40 cm⁻¹.

Propagation speed measures vary depending of the methodology used. When macroscopic waves are recorded from EEG or from ECoG which have low spatial and high temporal resolutions, the propagation speeds varies between 1 and 10 m/s. As indicated by Muller et al. [3], these results are compatible with the range of axonal conduction speeds of myelinated white matter fibers in the cortex. However, when measuring mesoscopic waves’ propagation speed using local field potential (LFP) from multielectrode arrays (MEAs) or from optical imaging signals recorded with voltage-sensitive dyes (VSDs) having high spatial and temporal resolution, the propagation speeds varies from 0.1 to 0.8 m per second, consistent with the axonal conduction speed of the unmyelinated long-range horizontal fibers within the superficial layers of the cortex, as indicated by Muller et al. [3].
**Summary of parameters.**

Parameters of the model are summarized in the following table. Let us note that all parameters are dimensionless. They can be recalculated to dimensional values as in the example above.

<table><caption>**Table A1.** Typical values of parameters of the model.</caption><thead><tr><th>Parameter</th><th>Name</th><th>Unit</th><th>Typical Value</th></tr></thead><tbody><tr><td><i>D</i></td><td>diffusion coefficient</td><td>length<sup>2</sup>/time</td><td>10<sup>-4</sup></td></tr><tr><td><i>L</i></td><td>length of the interval</td><td>length</td><td>2</td></tr><tr><td><i>a</i><sub>1</sub>, <i>a</i><sub>2</sub>, <i>a</i><sub>3</sub>, <i>a</i><sub>4</sub></td><td>factors in connectivity functions</td><td>1/length</td><td>1 &divide; 5</td></tr><tr><td><i>b</i><sub>1</sub>, <i>b</i><sub>2</sub>, <i>b</i><sub>3</sub>, <i>b</i><sub>4</sub></td><td>exponents in connectivity functions</td><td>1/length</td><td>20 &divide; 40</td></tr><tr><td><i>S</i>'<sub>a</sub>(0), <i>S</i>'<sub>i</sub>(0)</td><td>growth rate of response functions</td><td>1/time</td><td>20</td></tr><tr><td><i>&tau;</i><sub>a</sub>, <i>&tau;</i><sub>i</sub></td><td>time delay in response functions</td><td>time</td><td>1 &divide; 10</td></tr><tr><td><i>&sigma;</i></td><td>potential decay rate</td><td>1/time</td><td>0.01</td></tr></tbody></table>

**References**

1. Kopell, N.; Gritton, H.J.; Whittington, M.A.; Kramer, M.A. Beyond The Connectome: The Dynome. *Neuron* **2014**, *83*, 1319–1328. [CrossRef] [PubMed]

2. Sporns, O. The human connectome: Origins and challenges. *NeuroImage* **2013**, *80*, 53–61. [CrossRef] [PubMed]

3. Muller, L.; Chavane, F.; Reynolds, J.; Sejnowski, T.J. Cortical travelling waves: Mechanisms and computational principles. *Nat. Rev. Neurosci.* **2018**, *19*, 255–268. [CrossRef] [PubMed]

4. Wu, J.Y.; Huang, X.; Zhang, C. Propagating Waves of Activity in the Neocortex: What They Are, What They Do. *Neuroscientist* **2008**, *14*, 487–502. [CrossRef]

5. Zhang, H.; Watrous, A.J.; Patel, A.; Jacobs, J. Theta and Alpha Oscillations Are Travelling Waves in the Human Neocortex. *Neuron* **2018**, *98*, 1269–1281. [CrossRef] [PubMed]

6. Botella-Soler, V.; Valderrama, M.; Crépon, B.; Navarro, V.; Quyen, M.L.V. Large-Scale Cortical Dynamics of Sleep Slow Waves. *PLoS ONE* **2012**, *7*, e30757. [CrossRef]

7. Wilson, H.R.; Cowan, J.D. A Mathematical Theory of the Functional Dynamics of Cortical and Thalamic Nervous Tissue. *Kybernetik* **1973**, *13*, 55–80. [CrossRef]

8. Senk, J.; Korvasova, K.; Schuecker, J.; Hagen, E.; Tetzlaff, T.; Diesmann, M.; Helias, M. Conditions for travelling waves in spiking neural networks. *arXiv* **2018**, arXiv:1801.06046v1.

9. Atay, F.M.; Hutt, A. Neural Fields with Distributed Transmission Speeds and Long-Range Feedback Delays. *SIAM J. Appl. Dyn. Syst.* **2006**, *5*, 670–698. [CrossRef]

10. Venkov, N.A.; Coombes, S.; Matthews, P.C. Dynamic instabilities in scalar neural field equations with space-dependent delays. *Physica D* **2007**, *232*, 1–15. [CrossRef]

11. Ermentrout, G.B.; Folias, S.E.; Kilpatrick, Z.P. Spatiotemporal Pattern Formation in Neural Fields with Linear Adaptation. In *Neural Fields*; Coombes, S., beim Graben, P., Potthast, R., Wright, J., Eds.; Springer: Berlin/Heidelberg, Germany, 2014; pp. 119–151.

12. Meijer, H.G.E.; Coombes, S. Travelling waves in a neural field model with refractoriness. *J. Math. Biol.* **2014**, *68*, 1249–1268. [CrossRef] [PubMed]

13. Pinotsis, D.A.; Hansen, E.; Friston, K.J.; Jirsa, V.K. Anatomical connectivity and the resting state activity of large cortical networks. *Neuroimage* **2013**, *65*, 127–138. [CrossRef] [PubMed]

14. Modolo, J.; Bhattacharya, B.; Edwards, R.; Campagnaud, J.; Legros, A.; Beuter, A. Using a virtual cortical module implementing a neural field model to modulate brain rhythms in Parkinson’s disease. *Front. Neurosci.* **2010**, *4*, 45. [CrossRef] [PubMed]

15. Ermentrout, B.; McLeod, J.B. Existence and uniqueness of travelling waves for a neural network. *Proc. R. Soc. Edinburgh* **1994**, *134A*, 1013–1022. [CrossRef]

16. Chen, Z.; Ermentrout, B.; Wang, X.J. Wave Propagation Mediated by GABAB Synapse and Rebound Excitation in an Inhibitory Network: A Reduced Model Approach. *J. Comput. Neurosci.* **1998**, *5*, 53–69. [CrossRef]
17. Amari, S. Dynamics of pattern formation in lateral-inhibition type neural fields. *Biol. Cybern.* **1977**, *27*, 77–87.
[CrossRef]

18. Bessonov, N.; Beuter, A.; Trofimchuk, S.; Volpert, V. Estimate of the travelling wave speed for an integro-differential equation. *Appl. Math. Lett.* **2019**, *88*, 103–110. [CrossRef]

19. Buzsaki, G.; Anastassiou, C.A.; Koch, C. The origin of extracellular fields and currents—EEG, ECoG, LFP and spikes. *Nat. Rev. Neurosci.* **2016**, *13*, 407–420. [CrossRef]

20. Pinto, D.J.; Ermentrout, G.B. Spatially structured activity in synapticalaty coupled neuronal networks: I. Travelling fronts and pulses. *SIAM J. Appl. Math.* **2001**, *62*, 206–225. [CrossRef]

21. Kilpatrick, Z.P. Coupling layers regularizes wave propagation in laminar stochastic neural fields. *Phys. Rev. E* **2014**, *89*, 022706. [CrossRef]

22. Moussaoui, A.; Volpert, V. Speed of wave propagation for a nonlocal reaction-diffusion equation. *Appl. Anal.* **2018**, *1*, 1–15. [CrossRef]

23. Rabiller, G.; He, J.-W.; Nishijima, Y.; Wong, A.; Liu, J. Perturbation of brain oscillations after ischemic stroke: a potential biomarker for post-stroke function and therapy. *Int. J. Mol. Sci.* **2015**, *16*, 25605–25640. [CrossRef]

24. Beuter, A.; Balossier, A.; Trofimchuk, S.; Volpert, V. Modelling of post-stroke stimulation of cortical tissue. *Math. Biosci.* **2018**, *305*, 146–159. [CrossRef] [PubMed]

25. Rapela, J. Travelling waves appear and disappear in unison with produced speech. *arXiv* **2018**, arXiv:1806.09559v1.

26. Pelt, J.; van Ooyen, A. Estimating neuronal connectivity from axonal and dendritic density fields. *Front. Comput. Neurosci.* **2013**, *7*, 160.

27. van Ooyen, A.; Carnell, A.; de Ridder, S.; Tarigan, B.; Mansvelder, H.D.; Bijma, F.; Gunst, M.D.; Pelt, J.V. Independently Outgrowing Neurons and Geometry- Based Synapse Formation Produce Networks with Realistic Synaptic Connectivity. *PLoS ONE* **2014**, *9*, e85858. [CrossRef]

28. Carter, A.R.; Connor, L.T.; Dromerick, A.W. Rehabilitation after stroke: Current state of the science. *Curr. Neurol. Neurosci. Rep.* **2010**, *10*, 158–166. [CrossRef]

29. Gross, J.; Hoogenboom, N.; Thut, G.; Schyns, P.; Panzeri, S.; Belin, P.; Garrod, S. Speech Rhythms and Multiplexed Oscillatory Sensory Coding in the Human Brain. *PLoS. Biol.* **2013**, *11*, e1001752. [CrossRef]

30. Weaver, J. How Brain Waves Help Us Make Sense of Speech. *PLoS Biol.* **2013**, *11*, e1001753. doi:10.1371/journal.pbio.1001753. [CrossRef]

31. Sarubbo, S.; Benedictis, A.D.; Merler, S.; Mandonnet, E.; Balbi, S.; Granieri, E.; Duffau, H. Towards a Functional Atlas of Human White Matter. *Hum. Brain Mapp.* **2015**. [CrossRef]

32. Pinto, D.J.; Ermentrout, G.B. Spatially structured activity in synapticalaty coupled neuronal networks: II. Lateral inhibition and standing pulses. *SIAM J. Appl. Math.* **2001**, *62*, 226–243. [CrossRef]

33. Gerstner, W.; Kistler, W.M.; Naud, R.; Paninski, L. *Neuronal Dynamics: From Single Neurons to Networks and Models of Cognition*; Cambridge University Press: Cambridge, UK, 2014.

<figure><img src="image_10.png" /></figure>

© 2020 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).