Distributed Nonlocal Feedback Delays May Destabilize Fronts in Neural Fields, Distributed Transmission Delays Do Not

The spread of activity in neural populations is a well-known phenomenon. To understand the propagation speed and the stability of stationary fronts in neural populations, the present work considers a neural field model that involves intracortical and cortico-cortical synaptic interactions. This includes distributions of axonal transmission speeds and nonlocal feedback delays as well as general classes of synaptic interactions. The work proves the spectral stability of standing and traveling fronts subject to general transmission speeds for large classes of spatial interactions and derives conditions for the front instabilities subjected to nonlocal feedback delays. Moreover, it turns out that the uniqueness of the stationary traveling fronts guarantees its exponential stability for vanishing feedback delay. Numerical simulations complement the analytical findings.

The spatially-extended neural network under study implies axonal connections with finite transmission speeds, which essentially leads to transmission delays. This delay depends strongly on the axonal branching architecture [2] and the degree of myelination of axonal branches [42]. For instance, unmyelinated axons exhibit a small transmission speed in the range of 0.1-1.0 m/s and occur mainly in short-range intracortical connections. In long-range axonal fibers such as cortico-cortical connections, the axons are myelinated leading to a faster transmission speed in the range of 1 m/s-100 m/s. Consequently, the resulting transmission delay between two spatial locations depend on the axonal paths and varies between 0.5 ms and 100 ms. Since these delay times are in the same range as time constants of synaptic responses of tens of milliseconds, effects of finite transmission delays on the spatio-temporal evolution of activity occur [6,30,46]. Although one may estimate such mean transmission delay times along axonal fibers, physiological studies point out that the transmission speed depends on the specific path the action potential takes, and hence varies in a single axonal branching structure from one neuron to another [42]. In addition, the detailed branching structure of neural tissue changes on a time scale of few months [45] or even few days [21], and hence changes the transmission speed. Consequently, it is not reasonable to study the effect of a single transmission speed, but rather a distribution of speeds. The present work considers such a distribution of transmission speeds and extends previous studies assuming a single speed only [13,30,37,38,46].
The model under study considers two types of axonal pathways. In intra-cortical connections, the lengths of axonal paths may vary due to the absence of fiber bundles of fixed length yielding a transmission delay proportional to the distance between two spatial locations. In contrast, cortico-cortical feedback connections may exhibit fiber bundles with fixed length yielding a constant feedback delay. By virtue of the distribution of transmission speeds, these two pathways are modeled by a distribution of transmission speeds and distributed feedback delays.
Most recent studies of extended neural networks considered either transmission speeds in intraarea connections or delay in feedback connections, although experimental findings indicate the presence of both connections [8], a distribution of transmission speeds [20] and feedback delays [8]. Only few previous studies consider both a single transmission delay in intraarea connections and the feedback delay [24,28,36,51]. The present work extends these studies by a detailed spectral stability analysis.
It is well known that transmission delays and feedback delays may destabilize stationary activity yielding oscillatory phenomena [4,7,30]. To better understand the spatio-temporal dynamics of neural populations, it is essential to comprehend the role of delays and their effects on the activity propagation. The present work undertakes the analysis of stationary front for the general case of large classes of delay distributions and spatial interactions and, therefore, aims to reveal answers for this problem.

Methods
This section introduces the model equation and gives the major previous results, which represent the basis of the novel results given in Sect. 3. To not expand the section too much, the discussion of the already published findings is kept short and concise.

The Model Equation
The model under study describes the activity in neural populations on a mesoscopic spatial and temporal scale with typical spatial range of 500 µm and temporal time constants of 5-10 ms. The corresponding spatial domain is coarse-grained and exhibits spatial patches. This structure reflects the macrocolumnar structure observed in primary sensory areas [19]. Such a neural field model allows the successful reproduction of electroencephalographic activity on the head [35] and the successful description of spiral waves in neural tissue [23]. More precisely, the neural field model describes a rate-coding neural population involving synaptic interactions whose spatio-temporal evolution obeys the following nonlinear scalar integro-differential equation: where u = u(x, t) represents the mean membrane potential of a spatial patch at position x and time t [4] and u o is the initial activity. This model neglects units since they can be considered by appropriate scaling of time and space [27]. The prefactors α ≥ 0 and β ≥ 0 are nonnegative constants and reflect the synaptic weights of intraarea connections and feedback connections, respectively, and α + β > 0. The function ξ ≥ 0 represents the probability density distribution of axonal transmission speeds. Similarly, the function η ≥ 0 represents the probability density function for feedback delays. We would like to point out that Eq. (1) describes the evolution of a minimal scalar model which, however, considers a large number of features. In spite of its reduced form, i.e., scalar and activity evolves in a one-dimensional space, and it promises to give insights into the effect of various delay types.
To learn more about the wave speed, at some point in the work, we will investigate the dependence of the wave speed of the traveling wave front on the transmission speeds and delays. To this end, the authors assume the axonal transmission speed distribution and the feedback delay distribution to a sum of two terms i.e., the presence of two axonal transmission speeds and two feedback delays. This choice reflects the presence of short intracortical connections showing a small transmission speed [31] and small constant feedback delay and long-range cortico-cortical connections, which exhibit large axonal transmission speeds and large feedback delays.
In addition, S = H (u − θ) denotes the transfer function of the model and is chosen to the Heaviside step function: H (u − θ) = 0 for all u < θ, H (0) = 1 2 , and H (u − θ) = 1 for all u > θ. This assumption is valid for identical neurons in the population [27]. Although this is a strong approximation, it gives first insights into the possible dynamics of homogeneous neural populations. The parameter θ is constant and represents the mean firing threshold of the neurons.
The spatial kernel functions K and W are real-valued and reflect two different nonlocal axonal connectivities between neurons and synapses: K represents intracortical axonal interaction in the neural population, which exhibits finite transmission speeds, whereas W denotes the axonal connectivity along fibers that leave the neural population and reenter it with a constant delay. These functions may be seen as weighted sums of probability density functions of axonal connections of subnetworks where the weights represent the synaptic strengths in each subnetwork [27]. We require that ξ = 0 at least in a small open interval (0, c 0 ), where c 0 > 0 is a positive constant. This assumption is reasonable since the resulting traveling wave front propagates with a wave speed μ 0 , 0 < μ 0 < c 0 that is equal to or smaller than the transmission speed c 0 due to physical reasoning. It is assumed that for two positive constants C > 0 and ρ > 0.
To investigate various superpositions of subnetworks, the paper considers the following three general classes of synaptic interactions: The synaptic feedback coupling W satisfies the same assumptions as K and it belongs to one of the three classes. Certainly, K and W are not necessarily in the same class. Moreover, K and W are not necessarily symmetric functions.

The Standing Wave Front
In the presence of a single firing threshold θ as in Eq. (1), a traveling wave front connects two constant states, i.e., one above and one below the threshold θ . Let β = 0 and choose u(x, t) = u(t) in Eq. (1), then with two exponentially stable stationary states U stationary−0 = 0 and U stationary−1 = α.
They are stable irrespective of the spatial interactions and the delay distributions. Suppose that α + β = 2θ and αK(0) + βW (0) > 0. Now let us consider a stationary standing wave front with profile U = U(x). Without loss of generality, suppose that the front crosses the threshold θ at the point Plugging such a solution back into the equation, noting that the solution is independent of time, we get Amari had studied solutions similar to (5) for general stationary states in his celebrated work [3]. Section 3.1 will elaborate in some detail the spectral stability of the stationary standing wave front given by (5). It has been shown in many previous studies [4,27] that certain spatial interactions may destabilize spatially constant stationary states subject to the nonlinear gain func-tion S = H (u − θ). In the present model, the functional derivative of the nonlinear gain function is given by δS[u(x, t)]/δu(x, t) = δ(u(x, t) − θ). In other words, spatial interactions at a certain spatial location x with activity u(x, t) contribute to the stability of a stationary state only if the location is close to the threshold θ .

The Traveling Wave Front and Its Stability
After the study of the existence of a standing wave front, this subsection focuses on the existence of a travelling wave front. Let α ≥ 0, β ≥ 0, θ > 0 be constants such that 0 < 2θ < α + β. As shown in the previous section, there are two constant solutions U 1 = 0 and U 2 = α + β to Eq. (1).

The Front Shape
Suppose that u(x, t) = U(x + μt) is a traveling wave front of Eq. (1), where μ > 0 represents the wave speed and z = x + μt represents a moving coordinate. Due to translation invariance, suppose that the traveling wave front satisfies the condi- Suppose that the front satisfies the boundary conditions lim z→−∞ U(z) = 0, lim z→∞ U(z) = α + β and lim z→±∞ U (z) = 0. For physical reasons, suppose that the wave speed satisfies the conditions 0 < μ < c 0 , where c 0 = sup{c > 0 : ξ = 0 on (0, c) and ξ ≥ 0 on (c, ∞)} is the smallest occurring transmission speed. Then the traveling wave front U = U(z) and the wave speed μ satisfy the equation After a series of change of variables (such as ω = y − μ c |z − y| and x = c c+s(z−ω)μ × (z − ω), etc.), this above equation becomes is the sign function, which is defined by s(x) = −1 for all x < 0, s(0) = 0 and s(x) = 1 for all x > 0. By using the integrating factor idea and integration by parts, we find the representation of the front The derivative of U = U(z) is given by These expressions are useful in a later part of the work.

The Wave Speed of the Traveling Wave Front
To compute the wave speed, setting z = 0 and U(0) = θ in (8) and making some simple change of variables yields This equation may be rewritten as where Define the sub speed index functions φ 21 = φ 21 (μ) and φ 22 There holds φ 2 (μ) = φ 21 (μ) + φ 22 (μ). We will see that there exists a unique positive solution (that is, the wave speed μ 0 ) to the equation φ(μ) = α+β 2 − θ . Interestingly, the speed index functions φ 1 and φ 2 allow us to express the slope of the front at the threshold θ by utilizing (9) and (10) These expressions will be useful in later discussions of the spectral stability of the traveling wave front. Moreover, for single transmission speed and single feedback delay (i.e., Thus, there is no distinction between intracortical and feedback interactions. This is obvious from Eq. (1) where the two integrals may be written as a single one.

The Spectral Stability of the Traveling Wave Front
Real neural structures exhibit a certain level of background noise, which may disturb the propagation of activity. Hence, it is important to study the stability of the stationary traveling wave front with respect to small perturbations. This kind of analysis has already been performed before for similar equations in [4,14,16,43]. Previous studies on integral and/or partial differential equations involving infinite delays have found criteria for the existence, uniqueness, and stability of traveling wave solutions [1,18]. Following the successful extrapolation approach for infinite delays in integral and/or partial differential equations based on finite delays [12,18] and recalling the close relationship of partial differential equations and integrodifferential equations of the type discussed here [25], we assume in the following that mathematical analysis based on finite delays are applicable. An additional discussion of this approximation may appear necessary, but is neglected since it would exceed the major aim of the present work. Motivated by a previous study of Coombes and Owen [10], we study the spectral stability of the traveling wave front of (1) by constructing Evans functions which are well known from the literature of integral and/or partial differential equations [43].
To study the spectral stability of the front, we must rewrite the equation under consideration in moving coordinate and linearize the new equation with respect to the traveling wave front to obtain a linear equation. Then we seek solutions of the form exp(λt)ψ(z) to separate time and the moving coordinate z and to obtain an eigenvalue problem with the complex eigenvalue λ and the eigenfunction ψ; see also the Appendix. Solving the eigenvalue problem yields the definition of the Evans function. It turns out that a complex number λ 0 is an eigenvalue of the eigenvalue problem if and only if λ 0 is a zero of the Evans function. Note that λ 0 = 0 is always an eigenvalue, reflecting the translation invariance of the travelling wave front. Previous studies have shown that the equation E(λ) = 0 determines the isolated spectrum {λ n : n = 1, 2, 3, . . .}. As before, the stationary front is unstable if there exists some eigenvalue λ 0 with positive real part Re λ 0 > 0 or if the neutral eigenvalue λ 0 = 0 is not simple.
The spectrum of the linear differential operator L(λ) consists of two parts: the essential spectrum and the isolated spectrum (point spectrum comprising the eigenvalues). Please see the Appendix for the definitions of the linear differential operators L(λ) and L 0 . Are there any other spectrum in Ω = {λ ∈ C : Re λ > −1} other than the essential spectrum and the isolated spectrum to the linear differential operator L(λ)? By assumption, the kernel functions K and W converge to zero exponentially fast as x → ±∞. Therefore, the operator [L(λ) − L 0 ](L 0 ) −1 is a compact operator in C 0 (R). The residual spectrum does not exist in our model.
The essential spectrum is easy to calculate by following the original ideas of John Evans [14]. The complex number λ ∈ C belongs to the essential spectrum if and only if λ is a complex number such that the solutions of the differential equation is bounded on R. The solution of this equation is given by ψ(z) = C exp(− λ+1 μ z). This solution is bounded on R if and only if λ = −1 + ir, for some real number r. Since the subsequent sections show that linear deviations from stationary traveling wave front obey differential equations of the type (17), we find that the essential spectrum contains those values of λ with Re λ = −1 < 0 only, i.e., the essential spectrum does not threaten the stability of the traveling wave front.
It remains to investigate whether the isolated spectrum (the eigenvalues) threatens the stability, i.e., whether there are eigenvalues λ 0 with positive real part Re λ 0 > 0, threatening the stability of the stationary front.
Define the Evans function E = E(λ) for the traveling wave front of Eq. (1) by where E 1 = E 1 (λ) and E 2 = E 2 (λ) are also called Evans functions, defined by in the open domain Ω = {λ ∈ C : Re λ > −1}. Here, the Evans function E = E(λ) is also called the stability index function and represents a sum of two single stability index functions E 1 = E 1 (λ) and E 2 = E 2 (λ), reflecting the contribution of transmission delay and constant delay, respectively.

Properties of the Evans functions:
A complex number λ 0 is an eigenvalue of the eigenvalue problem L(λ)ψ = λψ (details to be given in the Appendix) if and only if λ 0 is a zero of the Evans function. Moreover,

The Numerical Simulations
To complement the analytical study, subsequent sections show numerical integrations of the model equation (1). To this end, if not stated otherwise, the spatial kernel functions are chosen to with the real parameters s > r, ρ > 1, σ > 0. Equation (1) is integrated numerically in time by an Euler-forward method with time step t = 0.02, the spatial integral has been computed on a circular spatial grid of length L = 60 and 600 intervals by the Monte Carlo-type VEGAS-algorithm [22] with 5000 random draws. If not stated differently, the initial values of the neural activity have been chosen to the analytical stationary front U(x − cT ) perturbed by random values γ (x, t) taken from a uniform distribution γ (x, T ) ∈ [−0.025, 0.025] in the initial interval −12 ≤ T ≤ 0.

Results
This section presents the new findings on the effect of distributed delays in both standing wave front and traveling wave front. They extend previous results mentioned in the previous section.

The Standing Wave Front
Derivation of an eigenvalue problem. Recall that the standing wave front U = U(x) is a time independent solution of the nonlinear scalar integro-differential equation To study the spectral stability of the stationary standing wave front, we linearize this equation about the front U = U(x) and we find that v(x, t) = u(x, t) − U(x) obeys the differential equation (neglecting higher order terms) . Making a change of variable and using (7), we find the equation where αK(0) + βW (0) > 0. We understand that the contribution of the spatial interactions to the activity is maximum if (0)] are large. For kernel functions with a maximum value at the origin, the largest contribution is expected at x ≈ x 0 , i.e., close to the threshold θ . Hence, one expects a large change of the activity close to the threshold θ . This is different from other spatial interactions, such as gamma-distributed kernel functions [27], where the strongest contribution occurs away from the threshold. A standing wave front is translation invariant. We may let x 0 = 0. Suppose that v(x, t) = ψ(x) exp(λt) is a solution of the above equation, where λ is a complex constant and ψ = ψ(z) is a bounded continuous function defined on R. After canceling out exp(λt), we obtain the eigenvalue problem The spectral stability of the stationary front (5) is determined by the eigenvalues of this eigenvalue problem. Letting x = 0 in this equation and canceling out ψ(0) because it is not equal to zero, we get At first, we observe that the spectral stability of the standing wave front just depends on the spatial self- For all kernel functions K and W , for all probability density functions ξ and η, f is an increasing of λ. There exists a unique solution λ 0 = 0 to the equation f (λ) = 0. Moreover, λ 0 = 0 is a simple solution (that is, λ 0 = 0 is a simple eigenvalue). Therefore, the standing wave front is spectrally stable. Let us consider a very special case: absent intraarea connections (i.e., α = 0) and focus on nonlocal feedback connections with a single feedback delay τ 0 (i.e. η(τ ) = δ(τ − τ 0 ), where τ 0 > 0 is a positive constant). Then λ + 1 = exp(−λτ 0 ). That is (λ + 1) exp(λτ 0 ) = 1. Given the positive constant τ 0 > 0, it is very easy to see that there exists a unique solution λ 0 = 0 to (λ + 1) exp(λτ 0 ) = 1.
Moreover, Fig. 1 shows the space-time activity of oscillatory unstable standing front gained by a numerical simulation of Eq. (1) for a single delay (a) and a distribution of two delays (b). We observe oscillatory activity close to the threshold value θ consistent with the reasoning in Sect. 2.2.

Traveling Wave Front
This subsection investigates analytically the uniqueness of the stationary travelling wave front and the spectral stability of the front. Numerical simulations validate the analytical findings.

The Uniqueness
Taking a close look at the implicit equation (10), the question arises whether its solution, i.e., the wave speed μ 0 , is unique, or whether there are several possible traveling wave fronts with different wave speeds. Previous studies on propagating front in neural fields involving a single axonal finite transmission speed [11] have established the uniqueness of the traveling wave front. The present work extends these studies by considering distributed transmission speeds and distributed feedback delays.
If the speed index functions φ 1 (μ), φ 2 (μ), and φ(μ) defined in (11), (12), and (10) are monotonic in μ, then the wave speed of the traveling wave front is unique. By using rigorous mathematical analysis (to keep the paper from too long, details are not to be given here), we find that • For all kernel functions (K, W ) in classes (A) and (B), for all speed and delay distributions (ξ, η), for all μ ∈ (0, c 0 ), there hold the following estimates: • For all kernel functions (K, W ) in class (C), for all speed and delay distributions (ξ, η), for all μ ∈ (0, c 0 ), there hold the following estimates: for two positive constants μ * and μ * * , where 0 < μ * < c 0 and 0 < μ * * < c 0 .
Overall, there exists a unique wave speed and there is a unique traveling wave front to Eq. (1).

Dependence on Transmission Speed and Feedback Delays
Now we investigate the change of the wave speed while changing the speed and delay distributions assuming the distributions (2). For simplicity, we choose c 2 /c 1 =

Stability Analysis
Analysis of the eigenvalues of the eigenvalue problem. The introduction of the speed index functions and the stability index functions given in Sect. 2.3.3 reveals novel relationship between the uniqueness of the wave speed and the stability of the traveling wave front. Comparing the stability index functions given in Eqs. (19), (20) and the speed index functions in Eqs. (11), (12), we can obtain the following results.
(III) For both intracortical and feedback connections (i.e., α > 0, β > 0, and 0 < 2θ < α + β), the stability index function reads Similar to cases (I) and (II), we obtain the relationship By virtue of the monotonicity of the speed index function φ 1 (μ) + φ 22 (μ), γ = 1 is the only solution of (29), i.e., λ 0 = 0 is the only eigenvalue and the front is spectrally stable. Now let us consider distributed feedback delays, the stability index function (18) and Eqs. (10), (15) yield Since | exp(λx/μ 0 )| < 1 for all λ with Re λ > 0, the same reasoning as in (I) and (II) applies and the right side of Eq. We write λ = x + iy and E(λ) = E real (λ) + iE imag (λ). Note that both the real part E real (λ) and the imaginary part E imag (λ) of the Evans function E(λ) are real harmonic functions of x and y. Hence, they satisfy the mean value formula; see Evans [14][15][16][17]. As a result, E(λ) also satisfies the mean value formula. Thus, |E(λ)| cannot attain a local maximum inside any open domain of C. By using a strong maximum principle of |E(λ)| in Ω, we find that 0 < |E(λ)| < 1, for all λ ∈ C with Re λ > 0. The spectral stability follows immediately. To illustrate this finding, Fig. 3 shows numerical simulations of Eq. (1) and we observe stationary propagating front for kernel functions in the three classes (A), (B), and (C). The simulations consider both distributed axonal transmission speeds and distributed feedback delays. We observe that starting from a noise-perturbed stationary front, the activity approaches the smooth stationary propagating front after some time, i.e., the front is also exponentially stable in the presence of feedback delay. This numerical result for nonvanishing feedback connections complements the analytical result for vanishing feedback delay above.

Stability of the Standing Wave Front
The results in Sect. 3.1 show that transmission delays do not destabilize standing wave front. In contrast, nonlocal feedback delays may destabilize standing wave front for certain delays resulting to oscillatory activity. The corresponding analytical conditions for a single delay hold for feedback connections in class (C), i.e., for local inhibition-lateral excitation interactions. Since such synaptic couplings as well as feedback delays are omnipresent in the brain, oscillatory standing wave front may occur frequently in real neural structures. For distributed delays, such instability may be possible, but no analytical method is known up to today proving oscillatory instability. Nevertheless, numerical simulations of standing wave front subject to distributed feedback delays reveal oscillatory instability as well.

Wave Speed of the Traveling Wave Front
Section 3.2 shows that the speed index function is monotonically increasing, i.e., ∂φ ∂μ > 0 on the interval where the wave speed exist. Consequently, the wave speed is unique. Moreover, the same subsection shows analytically that the increase of the transmission speed, i.e., the decrease of the transmission delay, increases the wave speed. Similarly, the increase of the feedback delays decreases the wave speed. Summarizing these results, increasing delays slows down the wave speed of the traveling wave front.

Stability of the Travelling Wave Front
The results in Sect. 3.2 indicate that traveling wave front is spectrally stable in the absence of delayed feedback due to the uniqueness of the wave speed. To our best knowledge, this relationship between the uniqueness and the stability has not been found before, although claimed in previous studies [50,51].
Focusing on the effect of delayed feedback, the stability analysis in Sect. 3.2.3 reveals the stability of traveling wave front for distributed delays. To support the analytical results, numerical simulations were performed showing stable traveling wave front, cf. Fig. 3. Additional extensive numerical studies (not shown) on the effect of distributed feedback delays have found stable traveling wave front only.
Summarizing the latter results, distributed transmission delays do not destabilize traveling wave front, but feedback delays may induce oscillatory instability.

Summary and Outlook
We find that the standing wave front and the stationary traveling wave front involving distributed transmission speeds exhibit a unique stable traveling wave front u(x, t) = U(x + μ 0 t) to Eq. (1) given any pair of synaptic couplings (K, W ), probability density functions (ξ , η), synaptic rate constants (α, β), any threshold θ and if assumptions (3) hold. In addition, the standing wave front and the stationary traveling wave front are spectrally stable in the presence of small external perturbations due to the found relationship between uniqueness and stability. In contrast, the additional presence of nonlocal feedback delays may render the stationary front instability. We find that the delay-induced loss of stability is oscillatory. This is valid for the present model involving a single synaptic time scale. It is worth mentioning, however, that neural fields involving multiple synaptic time scales are sensitive to transmission delays [33].
Future work may further investigate the relation of uniqueness and stability in other wave phenomena, such as traveling pulse solutions or global waves. Moreover, it is interesting to further study the impact of nonlocal delayed feedback on these wave phenomena. These studies may permit deeper insight into the role of delayed feedback loops which are omnipresent in neural systems.
(2) If there exists an eigenvalue with positive real part to the eigenvalue problem L(λ)ψ = λψ, then we say that the traveling wave front is exponentially unstable. (3) If max{Re λ : λ = 0, λ ∈ σ (L(λ))} < 0, that is, there exists no nonzero eigenvalue to the eigenvalue problem L(λ)ψ = λψ in the right half plane {λ ∈ C : Re λ > 0}, but λ 0 = 0 is not a simple eigenvalue of the eigenvalue problem L(λ)ψ = λψ, then we say that the traveling wave front is algebraically unstable.
Therefore, we have finished the proof of the spectral stability.