Stability of the stationary solutions of neural field equations with propagation delays

In this paper, we consider neural field equations with space-dependent delays. Neural fields are continuous assemblies of mesoscopic models arising when modeling macroscopic parts of the brain. They are modeled by nonlinear integro-differential equations. We rigorously prove, for the first time to our knowledge, sufficient conditions for the stability of their stationary solutions. We use two methods 1) the computation of the eigenvalues of the linear operator defined by the linearized equations and 2) the formulation of the problem as a fixed point problem. The first method involves tools of functional analysis and yields a new estimate of the semigroup of the previous linear operator using the eigenvalues of its infinitesimal generator. It yields a sufficient condition for stability which is independent of the characteristics of the delays. The second method allows us to find new sufficient conditions for the stability of stationary solutions which depend upon the values of the delays. These conditions are very easy to evaluate numerically. We illustrate the conservativeness of the bounds with a comparison with numerical simulation.


Introduction
Neural fields equations first appeared as a spatial-continuous extension of Hopfield networks with the seminal works of Wilson and Cowan, Amari [1,2]. These networks describe the mean activity of neural populations by nonlinear integral equations and play an important role in the modeling of various cortical areas including the visual cortex. They have been modified to take into account several relevant biological mechanisms like spike-frequency adaptation [3,4], the tuning properties of some populations [5] or the spatial organization of the populations of neurons [6]. In this work we focus on the role of the delays coming from the finite-velocity of signals in axons, dendrites or the time of synaptic transmission [7,8]. It turns out that delayed neural fields equations feature some interesting mathematical difficulties. The main question we address in the sequel is that of determining, once the stationary states of a non-delayed neural field equation are well-understood, what changes, if any, are caused by the introduction of propagation delays? We think this question is important since non-delayed neural field equations are pretty well understood by now, at least in terms of their stationary solutions, but the same is not true for their delayed versions which in many cases are better models closer to experimental findings. A lot of work has been done concerning the role of delays in waves propagation or in the linear stability of stationary states but except in [9] the method used reduces to the computation of the eigenvalues (which we call characteristic values) of the linearized equation in some analytically convenient cases (see [10]). Some results are known in the case of a finite number of neurons [11,12] and in the case of a few number of distinct delays [13,14]: the dynamical portrait is highly intricated even in the case of two neurons with delayed connections.
The purpose of this article is to propose a solid mathematical framework to characterize the dynamical properties of neural field systems with propagation delays and to show that it allows us to find sufficient delay-dependent bounds for the linear stability of the stationary states. This is a step in the direction of answering the question of how much delays can be introduced in a neural field model without destabilization. As a consequence one can infer in some cases without much extra work, from the analysis of a neural field model without propagation delays, the changes caused by the finite propagation times of signals. This framework also allows us to prove a linear stability principle to study the bifurcations of the solutions when varying the nonlinear gain and the propagation times.
The paper is organized as follows: in Section 2 we describe our model of delayed neural field, state our assumptions and prove that the resulting equations are wellposed and enjoy a unique bounded solution for all times. In Section 3 we give two different methods for expressing the linear stability of stationary cortical states, that is, of the time independent solutions of these equations. The first one, Section 3.1, is computationally intensive but accurate. The second one, Section 3.2, is much lighter in terms of computation but unfortunately leads to somewhat coarse approximations. Readers not interested in the theoretical and analytical developments can go directly to the summary of this section. We illustrate these abstract results in Section 4 by applying them to a detailed study of a simple but illuminating example.

The model
We consider the following neural field equations defined over an open bounded piece of cortex and/or feature space ⊂ R d . They describe the dynamics of the mean membrane potential of each of p neural populations.
J ij (r,r)S σ j V j t − τ ij (r,r),r − h j dr We give an interpretation of the various parameters and functions that appear in (1). is a finite piece of cortex and/or feature space and is represented as an open bounded set of R d . The vectors r andr represent points in .
The function S : R → (0, 1) is the normalized sigmoid function: It describes the relation between the firing rate ν i of population i as a function of the membrane potential, for example, V i : We note V the pdimensional vector (V 1 , . . . , V p ).
The p functions I ext i , i = 1, . . . , p, represent external currents from other cortical areas. We note I ext the p-dimensional vector (I ext 1 , . . . , I ext p ). The p × p matrix of functions J = {J ij } i,j =1,...,p represents the connectivity between populations i and j , see below.
The p real values h i , i = 1, . . . , p, determine the threshold of activity for each population, that is, the value of the membrane potential corresponding to 50% of the maximal activity.
The p real positive values σ i , i = 1, . . . , p, determine the slopes of the sigmoids at the origin.
Finally the p real positive values l i , i = 1, . . . , p, determine the speed at which each membrane potential decreases exponentially toward its rest value.
A difference with other studies is the intrinsic dynamics of the population given by the linear response of chemical synapses. In [9,15], ( d dt + l i ) is replaced by ( d dt + l i ) 2 to use the alpha function synaptic response. We use ( d dt + l i ) for simplicity although our analysis applies to more general intrinsic dynamics, see Proposition 3.10 in Section 3.1.3.
For the sake of generality, the propagation delays are not assumed to be identical for all populations, hence they are described by a matrix τ (r,r) whose element τ ij (r,r) is the propagation delay between population j atr and population i at r. The reason for this assumption is that it is still unclear from physiology if propagation delays are independent of the populations. We assume for technical reasons that τ is continuous, that is, τ ∈ C 0 ( 2 , R p×p + ). Moreover biological data indicate that τ is not a symmetric function (that is, τ ij (r,r) = τ ji (r, r)), thus no assumption is made about this symmetry unless otherwise stated.
In order to compute the righthand side of (1), we need to know the voltage V on some interval [−T , 0]. The value of T is obtained by considering the maximal delay: Hence we choose T = τ m .

The propagation-delay function
What are the possible choices for the propagation-delay function τ (r,r)? There are few papers dealing with this subject. Our analysis is built upon [16]. The authors of this paper study, inter alia, the relationship between the path length along axons from soma to synaptic buttons versus the Euclidean distance to the soma. They observe a linear relationship with a slope close to one. If we neglect the dendritic arbor, this means that if a neuron located at r is connected to another neuron located atr, the path length of this connection is very close to r −r , in other words, axons are straight lines. According to this, we will choose in the following: where c is the inverse of the propagation speed.

Mathematical framework
A convenient functional setting for the non-delayed neural field equations (see [17][18][19]) is to use the space F = L 2 ( , R p ) which is a Hilbert space endowed with the usual inner product: To give a meaning to (1), we define the history space which is the Banach phase space associated with equation (3) below. Using the notation V t (θ ) = V(t + θ), θ ∈ [−τ m , 0], we write (1) as: where is the linear continuous operator satisfying (the notation | · | is defined in Definition A.2 of Appendix A) |L 1 | ≤ J L 2 ( 2 ,R p×p ) . Notice that most of the papers on this subject assume infinite, hence requiring τ m = ∞. This raises difficult mathematical questions which we do not have to worry about, unlike [9,15,[20][21][22][23][24]. We first recall the following proposition whose proof appears in [25].
Notice that this result gives existence on R + , finite-time explosion is impossible for this delayed differential equation. Nevertheless, a particular solution could grow indefinitely, we now prove that this cannot happen.

Boundedness of solutions
A valid model of neural networks should only feature bounded membrane potentials. We find a bounded attracting set in the spirit of our previous work with non-delayed neural mass equations. The proof is almost the same as in [19] but some care has to be taken because of the delays.

Theorem 2.2 All the trajectories of the equation (3) are ultimately bounded by the same constant
We note l = min i=1,...,p l i and from Lemma B.2 (see Appendix B.1): show that the open ball of F of center 0 and radius R, B R , is stable under the dynamics of equation (3). We know that V(t) is defined for all t ≥ 0s and that f < 0 on ∂B R , the boundary of B R . We consider three cases for the initial condition V 0 .
If V 0 C < R and set T = sup{t|∀s ∈ [0, t], V(s) ∈ B R }. Suppose that T ∈ R, then V(T ) is defined and belongs to B R , the closure of B R , because B R is closed, in effect to ∂B R . We also have d Thus we deduce that for ε > 0 and small enough, V(T + ε) ∈ B R which contradicts the definition of T . Thus T / ∈ R and B R is stable.
Finally we consider the case V 0 ∈ B R . Suppose that ∀t > 0, V(t) / ∈B R , then ∀t > 0, d dt V 2 F ≤ −2δ, thus V(t) F is monotonically decreasing and reaches the value of R in finite time when V(t) reaches ∂B R . This contradicts our assumption.

Stability results
When studying a dynamical system, a good starting point is to look for invariant sets. Theorem 2.2 provides such an invariant set but it is a very large one, not sufficient to convey a good understanding of the system. Other invariant sets (included in the previous one) are stationary points. Notice that delayed and non-delayed equations share exactly the same stationary solutions, also called persistent states. We can therefore make good use of the harvest of results that are available about these persistent states which we note V f . Note that in most papers dealing with persistent states, the authors compute one of them and are satisfied with the study of the local dynamics around this particular stationary solution. Very few authors (we are aware only of [19,26]) address the problem of the computation of the whole set of persistent states. Despite these efforts they have yet been unable to get a complete grasp of the global dynamics. To summarize, in order to understand the impact of the propagation delays on the solutions of the neural field equations, it is necessary to know all their stationary solutions and the dynamics in the region where these stationary solutions lie. Unfortunately such knowledge is currently not available. Hence we must be content with studying the local dynamics around each persistent state (computed, for example, with the tools of [19]) with and without propagation delays. This is already, we think, a significant step forward toward understanding delayed neural field equations.
From now on we note V f a persistent state of (3) and study its stability.
We can identify at least three ways to do this: 1. to derive a Lyapunov functional, 2. to use a fixed point approach, 3. to determine the spectrum of the infinitesimal generator associated to the linearized equation.
Previous results concerning stability bounds in delayed neural mass equations are 'absolute' results that do not involve the delays: they provide a sufficient condition, independent of the delays, for the stability of the fixed point (see [15,[20][21][22]). The bound they find is similar to our second bound in Proposition 3.13. They 'proved' it by showing that if the condition was satisfied, the eigenvalues of the infinitesimal generator of the semi-group of the linearized equation had negative real parts. This is not sufficient because a more complete analysis of the spectrum (for example, the essential part) is necessary as shown below in order to proof that the semi-group is exponentially bounded. In our case we prove this assertion in the case of a bounded cortex (see Section 3.1). To our knowledge it is still unknown whether this is true in the case of an infinite cortex.
These authors also provide a delay-dependent sufficient condition to guarantee that no oscillatory instabilities can appear, that is, they give a condition that forbids the existence of solutions of the form e i(k·r+ωt) . However, this result does not give any information regarding stability of the stationary solution.
We use the second method cited above, the fixed point method, to prove a more general result which takes into account the delay terms. We also use both the second and the third method above, the spectral method, to prove the delay-independent bound from [15,[20][21][22]. We then evaluate the conservativeness of these two sufficient conditions. Note that the delay-independent bound has been correctly derived in [25] using the first method, the Lyapunov method. It might be of interest to explore its potential to derive a delay-dependent bound.
We write the linearized version of (3) as follows. We choose a persistent state V f and perform the change of variable where the linear operatorL 1 is given by ⎧ ⎨ It is also convenient to define the following operator:

Principle of linear stability analysis via characteristic values
We derive the stability of the persistent state V f (see [19]) for the equation (1) or equivalently (3) using the spectral properties of the infinitesimal generator. We prove that if the eigenvalues of the infinitesimal generator of the righthand side of (4) are in the left part of the complex plane, the stationary state U = 0 is asymptotically stable for equation (4). This result is difficult to prove because the spectrum (the main definitions for the spectrum of a linear operator are recalled in Appendix A) of the infinitesimal generator neither reduces to the point spectrum (set of eigenvalues of finite multiplicity) nor is contained in a cone of the complex plane C (such an operator is said to be sectorial). The 'principle of linear stability' is the fact that the linear stability of U is inherited by the state V f for the nonlinear equations (1) or (3). This result is stated in the Corollaries 3.7 and 3.8. Following [27][28][29][30][31], we note (T(t)) t≥0 the strongly continuous semigroup of (4) on C (see Definition A.3 in Appendix A) and A its infinitesimal generator. By definition, if U is the solution of (4) we have U t = T(t)φ. In order to prove the linear stability, we need to find a condition on the spectrum (A) of A which ensures that T(t) → 0 as t → ∞.
Such a 'principle' of linear stability was derived in [29,30]. Their assumptions implied that (A) was a pure point spectrum (it contained only eigenvalues) with the effect of simplifying the study of the linear stability because, in this case, one can link estimates of the semigroup T to the spectrum of A. This is not the case here (see Proposition 3.4).
When the spectrum of the infinitesimal generator does not only contain eigenvalues, we can use the result in [27, Chapter 4, Theorem 3.10 and Corollary 3.12] for eventually norm continuous semigroups (see Definition A.4 in Appendix A) which links the growth bound of the semigroup to the spectrum of A: Thus, U is uniformly exponentially stable for (4) if and only if We prove in Lemma 3.6 (see below) that (T(t)) t≥0 is eventually norm continuous. Let us start by computing the spectrum of A.

Computation of the spectrum of A
In this section we use L 1 forL 1 for simplicity.
We now apply results from the theory of delay equations in Banach spaces (see [27,28,31]) which give the expression of the infinitesimal generator Aφ =φ as well as its domain of definition The spectrum (A) consists of those λ ∈ C such that the operator (λ) of L(F ) defined by (λ) = λId + L 0 − J(λ) is non-invertible. We use the following definition:

Definition 3.2 (Characteristic values (CV))
The characteristic values of A are the λs such that (λ) has a kernel which is not reduced to 0, that is, is not injective.
It is easy to see that the CV are the eigenvalues of A. There are various ways to compute the spectrum of an operator in infinite dimensions. They are related to how the spectrum is partitioned (for example, continuous spectrum, point spectrum. . .). In the case of operators which are compact perturbations of the identity such as Fredholm operators, which is the case here, there is no continuous spectrum. Hence the most convenient way for us is to compute the point spectrum and the essential spectrum (see Appendix A). This is what we achieve next.

Remark 1
In finite dimension (that is, dim F < ∞), the spectrum of A consists only of CV. We show that this is not the case here.
Notice that most papers dealing with delayed neural field equations only compute the CV and numerically assess the linear stability (see [9,24,33]).
We now show that we can link the spectral properties of A to the spectral properties of L λ . This is important since the latter operator is easier to handle because it acts on a Hilbert space. We start with the following lemma (see [34] for similar results in a different setting).
Proof Let us define the following operator.
Let us now prove the lemma. We already know that Lemma 3.3 is the key to obtain (A). Note that it is true regardless of the form of L and could be applied to other types of delays in neural field equations. We now prove the important following proposition.

Proposition 3.4 A satisfies the following properties:
. We apply [35, Theorem IV. 5.26]. It shows that the essential spectrum does not change under compact perturbation. As J(λ) ∈ L(F ) is compact, we find ess (−L 0 + J(λ)) = ess (−L 0 ). As an example, Figure 1 shows the first 200 eigenvalues computed for a very simple model one-dimensional model. We notice that they accumulate at λ = −1 which is the essential spectrum. These eigenvalues have been computed using TraceDDE, [36], a very efficient method for computing the CVs.

Let us show that ess
Last but not least, we can prove that the CVs are almost all, that is, except for possibly a finite number of them, located on the left part of the complex plane. This indicates that the unstable manifold is always finite dimensional for the models we are considering here.
But Hence, for λ large enough 1 / ∈ P ((λId + L 0 ) −1 J(λ)), which holds by the spectral radius inequality. This relationship states that the CVs λ satisfying λ > −l are located in a bounded set of the right part of C; given that the CV are isolated, there is a finite number of them.

Stability results from the characteristic values
We start with a lemma stating regularity for (T(t)) t≥0 : Lemma 3. 6 The semigroup (T(t)) t≥0 of (4) is norm continuous on C for t > τ m .
Proof We first notice that −L 0 generates a norm continuous semigroup (in fact a group) S(t) = e −tL 0 on F and thatL 1 is continuous from C to F . The lemma follows directly from [27, Theorem VI.6.6].
Using the spectrum computed in Proposition 3.4, the previous lemma and the formula (5), we can state the asymptotic stability of the linear equation (4). Notice that because of Corollary 3.5, the supremum in (5) is in fact a max.
It is however possible (note that a regularity condition has to be verified but this is done easily in our case) to extend (see [34]) the semigroup T(t) to the space F × L 2 ([−τ m , 0], F ). We noteT(t) this extension which has the same spectrum as T(t). Indeed, we can consider integral solutions of (4) with initial condition U 0 in L 2 ([−τ m , 0], F ). However, as L 0 U 0 (0) has no meaning because φ → φ(0) is not continuous in L 2 ([−τ m , 0], F ), the linear problem (4) is not well-posed in this space. This is why we have to extend the state space in order to make the linear operator in (4) continuous. Hence the correct state space is F × L 2 ([−τ m , 0], F ) and any function φ ∈ C is represented by the vector (φ(0), φ). The variation of constant formula becomes: where π 2 is the projector on the second component. Now we choose ω = − max p (A)/2 > 0 and the spectral mapping theorem implies that there exists M > 0 such that |T(t) | C ≤ Me −ωt and and concludes the proof.
Finally, we can use the CVs to derive a sufficient stability result.
Proof Suppose that a CV λ of positive real part exists, this gives a vector in the Kernel of (λ). Using straightforward estimates, it implies that min i l i ≤ J · DS(V f ) F , a contradiction.

Generalization of the model
In the description of our model, we have pointed out a possible generalization. It concerns the linear response of the chemical synapses, that is, the lefthand side ( d dt + l i ) of (1). It can be replaced by a polynomial in d dt , namely P i ( d dt ), where the zeros of the polynomials P i have negative real parts. Indeed, in this case, when J is small, the network is stable. We obtain a diagonal matrix P( d dt ) such that P(0) = L 0 and change the initial condition (as in the theory of ODEs) while the history space becomes  (1) writes J ij (r,r)S σ j V j t − τ ij (r,r),r − θ j dr Introducing the classical variable where −L 0 is the Vandermonde (we put a minus sign in order to have a formulation very close to 1) matrix associated to P and (L 1 ) k,l=1,...,d s = (δ k=d s ,l=1 L 1 ) k,l=1,...,d s , ]. It appears that equation (7) has the same structure as (1): L 0 , L 1 , are bounded linear operators; we can conclude that there is a unique solution to (6). The linearized equation around a persistent states yields a strongly continuous semigroup T (t) which is eventually continuous. Hence the stability is given by the sign of max (A) where A is the infinitesimal generator of T (t). It is then routine to show that λ ∈ (A) ⇔ (λ) ≡ P(λ) − J(λ) non-invertible. This indicates that the essential spectrum ess (A) of A is equal to i Root(P i ) which is located in the left side of the complex plane. Thus the point spectrum is enough to characterize the linear stability: Proposition 3.10 If max p (A) < 0 the persistent solution V f of (6) is asymptotically stable.
Using the same proof as in [20], one can prove that max

Principle of linear stability analysis via fixed point theory
The idea behind this method (see [37]) is to write (4) as an integral equation. This integral equation is then interpreted as a fixed point problem. We already know that this problem has a unique solution in C 0 . However, by looking at the definition of the (Lyapunov) stability, we can express the stability as the existence of a solution of the fixed point problem in a smaller space S ⊂ C 0 . The existence of a solution in S gives the unique solution in C 0 . Hence, the method is to provide conditions for the fixed point problem to have a solution in S; in the two cases presented below, we use the Picard fixed point theorem to obtain these conditions. Usually this method gives conditions on the averaged quantities arising in (4) whereas a Lyapunov method would give conditions on the sign of the same quantities. There is no method to be preferred, rather both of them should be applied to obtain the best bounds.
In order to be able to derive our bounds we make the further assumption that there exists a β > 0 such that: Note that the notation τ −β represents the matrix of elements 1/τ β ij .
We rewrite (4) in two different integral forms to which we apply the fixed point method. The first integral form is obtained by a straightforward use the variation-ofparameters formula. It reads The second integral form is less obvious. Let us define Z(r, t) = drJ(r,r) t t−τ (r,r) ds U(r, s).
Note the slight abuse of notation, namely (J(r,r) t t−τ (r,r) ds U(r, s)) i = jJ ij (r, r) t t−τ ij (r,r) ds U j (r, s).

τ m ,t] U(s) F . This shows that ∀t, Z(t) ∈ F .
Hence we propose the second integral form: We have the following lemma.

Proof
The idea is to write the linearized equation as: By the variation-of-parameters formula we have: We then use an integration by parts: Using the two integral formulations of (4) we obtain sufficient conditions of stability, as stated in the following proposition: Proposition 3.13 If one of the following two conditions is satisfied: Proof We start with the first condition. The problem (4) is equivalent to solving the fixed point equation U = P 2 U for an initial condition φ ∈ C 0 . Let us define B = C 0 ([−τ m , ∞), F ) with the supremum norm written · ∞,F , as well as We define P 2 on S φ . For all ψ ∈ S φ we have P 2 ψ ∈ B and (P 2 ψ)(0) = φ(0). We want to show that P 2 S φ ⊂ S φ . We prove two properties.

P 2 ψ tends to zero at infinity.
Choose ψ ∈ S φ . Using Corollary B.3, we have Z(t) → 0 as t → ∞. Let 0 < T < t, we also have For the first term we write: Similarly, for the second term we write

Now for a given ε > 0 we choose T large enough so that α sup s∈[T −τ m ,∞) ψ(s)
Putting all this together, for all t > t * : From (9), it follows that P 2 ψ → 0 when t → ∞.
Since P 2 ψ is continuous and has a limit when t → ∞ it is bounded and therefore Using (9) for all ψ 1 , ψ 2 ∈ S φ we have We conclude from Picard theorem that the operator P 2 has a unique fixed point in S φ . There remains to link this fixed point to the definition of stability and first show that where U(t, φ) is the solution of (4).
Let us choose ε > 0 and M ≥ 1 such that |e (J−L 0 )t | F ≤ M. M exists because, by hypothesis, max (J − L 0 ) < 0. We then choose δ satisfying and φ ∈ C such that φ C ≤ δ. Next define We already know that P 2 is a contraction on S φ,ε (which is a complete space). The last thing to check is P 2 S φ,ε ⊂ S φ,ε , that is ∀ψ ∈ S φ,ε , P 2 ψ ∞,F < ε. Using Lemma B.3 in Appendix B.2: Thus P 2 has a unique fixed point U φ,ε in S φ,ε ∀φ, ε which is the solution of the linear delayed differential equation, that is, → 0 in C, we have proved the asymptotic stability for the linearized equation.
The proof of the second property is straightforward. If 0 is asymptotically stable for (4) all the CV are negative and Corollary 3.8 indicates that V f is asymptotically stable for (3).
The second condition says that The asymptotic stability follows using the same arguments as in the case of P 2 .
We next simplify the first condition of the previous proposition to make it more amenable to numerics.
Proof This corollary follows immediately from the following upperbound of the integral Then if there exists α < 1, β > 0 such that τ 3 2 +β m J τ β L 2 ( 2 ,R p×p ) (1 + M ε ε J − L 0 L 2 ( 2 ,R p×p ) ) ≤ α, it implies that condi-tion 1 in Proposition 3.13 is satisfied, from which the asymptotic stability of V f follows.
Notice that ε > 0 is equivalent to max (−L 0 +J) < 0. The previous corollary is useful in at least the following cases: • IfJ − L 0 is diagonalizable, with associated eigenvalues/eigenvectors: λ n ∈ C, e n ∈ F , thenJ − L 0 = n e λ n t e n ⊗ e n and |e (J−L 0 )t | F ≤ e −t max n λ n = e t max (−L 0 +J) . • If L 0 = l 0 Id and the range ofJ is finite dimensional:J(r, r ) = N k,l=1 J kl e k (r) ⊗ e l (r ) where (e k ) k∈N is an orthonormal basis of F , then e (J−L 0 )t = e −l 0 ·Id·t eJ t and |e (J−L 0 )t | F ≤ e −l 0 t |eJ t | F . Let us write J = (J kl ) k,l=1,...,N the matrix associated toJ (see above). Then eJ t is also a compact operator with finite

Remark 3 If we suppose that we have higher order time derivatives as in Section 3.1.3, we can write the linearized equation aṡ
Suppose that L 0 is diagonalizable then . Also notice that J =L 1 | F , |L 1 | (C) ds ≤ |L 1 | C . Then using the same functionals as in the proof of Proposition 3.13, we can find two bounds for the stability of a stationary state V f : To conclude, we have found an easy-to-compute formula for the stability of the persistent state V f . It can indeed be cumbersome to compute the CVs of neural field equations for different parameters in order to find the region of stability whereas the evaluation of the conditions in Corollary 3.14 is very easy numerically.
The conditions in Proposition 3.13 and Corollary 3.14 define a set of parameters for which V f is stable. Notice that these conditions are only sufficient conditions: if they are violated, V f may still remain stable. In order to find out whether the persistent state is destabilized we have to look at the characteristic values. Condition 1 in Proposition 3.13 indicates that if V f is a stable point for the non-delayed equation (see [18]) it is also stable for the delayed-equation. Thus, according to this condition, it is not possible to destabilize a stable persistent state by the introduction of small delays, which is indeed meaningful from the biological viewpoint. Moreover this condition gives an indication of the amount of delay one can introduce without changing the stability.
Condition 2 is not very useful as it is independent of the delays: no matter what they are, the stable point V f will remain stable. Also, if this condition is satisfied there is a unique stationary solution (see [18]) and the dynamics is trivial, that is, converging to the unique stationary point.

Summary of the different bounds and conclusion
The next proposition summarizes the results we have obtained in Proposition 3.13 and Corollary 3.14 for the stability of a stationary solution.

Proposition 3.15
If one of the following conditions is satisfied: The only general results known so far for the stability of the stationary solutions are those of Atay and Hutt (see, for example, [20]): they found a bound similar to condition 2 in Proposition 3.15 by using the CVs, but no proof of stability was given. Their condition involves the L 1 -norm of the connectivity function J and it was derived using the CVs in the same way as we did in the previous section. Thus our contribution with respect to condition 2 is that, once it is satisfied, the stationary solution is asymptotically stable: up until now this was numerically inferred on the basis of the CVs. We have proved it in two ways, first by using the CVs, and second by using the fixed point method which has the advantage of making the proof essentially trivial.
Condition 1 is of interest, because it allows one to find the minimal propagation delay that does not destabilize. Notice that this bound, though very easy to compute, overestimates the minimal speed. As mentioned above, the bounds in condition 1 are sufficient conditions for the stability of the stationary state V f . In order to evaluate the conservativeness of these bounds, we need to compare them to the stability predicted by the CVs. This is done in the next section.

Numerical application: neural fields on a ring
In order to evaluate the conservativeness of the bounds derived above we compute the CVs in a numerical example. This can be done in two ways: • Solve numerically the nonlinear equation satisfied by the CVs. This is possible when one has an explicit expression for the eigenvectors and periodic boundary conditions. It is the method used in [9].
• Discretize the history space C in order to obtain a matrix version A N of the linear operator A: the CVs are approximated by the eigenvalues of A N . Following the scheme of [36], it can be shown that the convergence of the eigenvalues of A N to the CVs is in O( 1 N N ) for a suitable discretization of C. One drawback of this method is the size of A N which can be very large in the case of several neuron populations in a two-dimensional cortex. A recent improvement (see [38]), based on a clever factorization of A N , allows a much faster computation of the CVs: this is the scheme we have been using.
The Matlab program used to compute the righthand side of (1) uses a Cpp code that can be run on multiprocessors (with the OpenMP library) to speed up computations. It uses a trapezoidal rule to compute the integral. The time stepper dde23 of Matlab is also used.
In order to make the computation of the eigenvectors very straightforward, we study a network on a ring, but notice that all the tools (analytical/numerical) presented here also apply to a generic cortex. We reduce our study to scalar neural fields ⊂ R and one neuronal population, p = 1. With this in mind the connectivity is chosen to be homogeneous J (x, y) = J (x − y) with J even. To respect topology, we assume the same for the propagation delay function τ (x, y).
We therefore consider the scalar equation with axonal delays defined on = (− π 2 , π 2 ) with periodic boundary conditions. Hence F = L 2 π ( , R) and J is also π -periodic.
where the sigmoid S 0 satisfies S 0 (0) = 0. Remember that (13) has a Lyapunov functional when c = 0 and that all trajectories are bounded. The trajectories of the non-delayed form of (13) are heteroclinic orbits and no non-constant periodic orbit is possible.
We are looking at the local dynamics near the trivial solution V f = 0. Thus we study how the CVs vary as functions of the nonlinear gain σ and the 'maximum' delay c. From the periodicity assumption, the eigenvectors of (λ) are the functions cos(nx), sin(nx) which leads to the characteristic equation for the eigenvalues λ: where J is the Fourier Transform of J and s 1 ≡ S 0 (0). This nonlinear scalar equation is solved with the Matlab Toolbox TraceDDE (see [36]). Recall that the eigenvectors of A are given by the functions θ → e λθ cos(nx), θ → e λθ sin(nx) ∈ C where λ is a solution of (14). A bifurcation point is a pair (c, σ ) for which equations (14) have a solution with zero real part. Bifurcations are important because they signal a change in stability, a set of parameters ensuring stability is enclosed (if bounded) by bifurcation curves. Notice that if σ 0 is a bifurcation point in the case c = 0, it remains a bifurcation point for the delayed system ∀c, hence ∀c, σ = σ 0 , 0 ∈ (A). This is why there is a bifurcation line σ = σ 0 in the bifurcation diagrams that are shown later. The bifurcation diagram depends on the choice of the delay function τ . As explained in Section 2.1, we use τ (x, y) = |x − y| π , where the lower index π indicates that it is a π -periodic function. The bifurcation diagram with respect to the parameters (c, σ ) is shown in the righthand part of Figure 2 in the case when the connectivity J is equal to J (x) = (−1 + 1.5 cos(2x)) 2 π . The two bounds derived in Section 3.3 are also shown. Note that the delay-dependent bound is computed using the fact that J ≡ DS(0)J = s 1 J is self-adjoint. They are clearly very conservative. The lefthand part of the same figure shows the delay function τ .
The first bound gives the minimal velocity 1/c below which the stationary state might be unstable, in this case, even for smaller speed, the state is stable as one can see from the CV boundary. Notice that in the parameter domain defined by the 2 conditions bound.1. and bound.2., the dynamic is very simple: it is characterized by a unique and asymptotically stable stationary state, V f = 0.
In Figure 3 we show the dynamics for different parameters corresponding to the points labelled 1, 2 and 3 in the righthand part of Figure 2 for a random (in space) and constant (in time) initial condition φ (see (1)). When the parameter values are below the bound computed with the CV, the dynamics converge to the stable stationary state V f = 0. Along the Pitchfork line labelled (P) in the righthand part of Figure 2, there is a static bifurcation leading to the birth of new stable stationary states, this is shown in the middle part of Figure 3. The Hopf curve labelled (H) in the righthand part of Figure 2 indicates the transition to oscillatory behaviors as one can see in the righthand part of Figure 3. Note that we have not proved that the Hopf curve was indeed a Hopf bifurcation curve, we have just inferred numerically from the CVs that the eigenvalues satisfy the usual conditions for the Hopf bifurcation.
Notice that the graph of the CVs shown in the righthand part of Figure 2 features some interesting points, for example, the Fold-Hopf point at the intersection of the   . 3 Plot of the solution of (13) for different parameters corresponding to the points shown as 1, 2 and 3 in the righthand part of Figure 2 for a random (in space) and constant (in time) initial condition, see text. The horizontal axis corresponds to space, the range is (− π 2 , π 2 ). The vertical axis represents time.
Pitchfork line and the Hopf curve. It is also possible that the multiplicity of the 0 eigenvalue could change on the Pitchfork line (P) to yield a Bogdanov-Takens point. These numerical simulations reveal that the Lyapunov function derived in [39] is likely to be incorrect. Indeed, if such a function existed, as its value decreases along trajectories, it must be constant on any periodic orbit which is not possible. However the third plot in Figure 3 strongly suggests that we have found an oscillatory trajectory produced by a Hopf bifurcation (which we did not prove mathematically): this oscillatory trajectory converges to a periodic orbit which contradicts the existence of a Lyapunov functional such as the one proposed in [39].
Let us comment on the tightness of the delay-dependent bound: as shown in Proposition 3.13, this bound involves the maximum delay value τ m and the norm J τ β L 2 ( 2 ,R p×p ) , hence the specific shape of the delay function, that is, τ (r,r) = c r −r 2 , is not completely taken into account in the bound. We can imagine many different delay functions with the same values for τ m and J τ β L 2 ( 2 ,R p×p ) that will cause possibly large changes to the dynamical portrait. For example, in the previous numerical application the singularity σ = σ 0 , corresponding to the fact that 0 ∈ p (A), is independent of the details of the shape of the delay function: however for specific delay functions, the multiplicity of this 0-eigenvalue could change as in the Bogdanov-Takens bifurcation which involves changes in the dynamical portrait compared to the pitchfork bifurcation. Similarly, an additional purely imaginary eigenvalue could emerge (as for c ≈ 3.7 in the numerical example) leading to a Fold-Hopf bifurcation. These instabilities depend on the expression of the delay function (and the connectivity function as well). These reasons explain why the bound in Proposition 3.13 is not very tight.
This suggests another way to attack the problem of the stability of fixed points: one could look for connectivity functionsJ which have the following property: for all delay function τ , the linearized equation (4) does not possess 'unstable solutions', that is, for all delay function τ , p (A) < 0. In the literature (see [40,41]), this is termed as the all-delay stability or the delay-independent stability. These remain questions for future work.

Conclusion
We have developed a theoretical framework for the study of neural field equations with propagation delays. This has allowed us to prove the existence, uniqueness, and the boundedness of the solutions to these equations under fairly general hypotheses.
We have then studied the stability of the stationary solutions of these equations. We have proved that the CVs are sufficient to characterize the linear stability of the stationary states. This was done using the semigroups theory (see [27]).
By formulating the stability of the stationary solutions as a fixed point problem we have found delay-dependent sufficient conditions. These conditions involve all the parameters in the delayed neural field equations, the connectivity function, the nonlinear gain and the delay function. Albeit seemingly very conservative they are useful in order to avoid the numerically intensive computation of the CV.
From the numerical viewpoint we have used two algorithms [36,38] to compute the eigenvalues of the linearized problem in order to evaluate the conservativeness of our conditions. A potential application is the study of the bifurcations of the delayed neural field equations.
By providing easy-to-compute sufficient conditions to quantify the impact of the delays on neural field equations we hope that our work will improve the study of models of cortical areas in which the propagation delays have so far been somewhat been neglected due to a partial lack of theory.

Appendix A: Operators and their spectra
We recall and gather in this appendix a number of definitions, results and hypotheses that are used in the body of the article to make it more self-sufficient.
Definition A.1 An operator T ∈ L(E, F ), E, F being Banach spaces, is closed if its graph is closed in the direct sum E ⊕ F . Definition A. 2 We note |J | F the operator norm of a bounded operator J ∈ L(F , F ), that is, It is known, see, for example, [35], that |J | F ≤ J L 2 ( 2 ,R p×p ) .

Definition A.3 A semigroup (T(t)) t≥0 on a Banach space E is strongly continuous
if ∀x ∈ E, t → T (t)x is continuous from R + to E.

Definition A.4 A semigroup (T(t)) t≥0 on a Banach space E is norm continuous if t → T (t) is continuous from R + to L(E).
It is said eventually norm continuous if it is norm continuous for t > t 0 ≥ 0.

Definition A.6 A closed operator T ∈ L(E) of a Banach space E is semi-Fredholm if dim N (T ) or codim R(T ) is finite and R(T ) is closed in E.
Definition A.7 If T ∈ L(E) is a closed operator of a Banach space E the essential spectrum ess (T ) is the set of λs in C such that λId − T is not semi-Fredholm, that is, either R(λId − T ) is not closed or R(λId − T ) is closed but dim N (λId − T ) = codim R(λId − T ) = ∞.
Proof By the Cauchy-Schwarz inequality and Lemma B.1:

B.2 Stability
In this section we prove Lemma B.3 which is central in establishing the first sufficient condition in Proposition 3.13.

Lemma B.3
Let β > 0 be such that τ −β ∈ L 2 ( 2 , R p×p )). Then we have the following bound: and if we set y i (r) = j dr J ij (r,r)