Neural field models with transmission delays and diffusion

A neural field models the large scale behaviour of large groups of neurons. We extend previous results for these models by including a diffusion term into the neural field, which models direct, electrical connections. We extend known and prove new sun-star calculus results for delay equations to be able to include diffusion and explicitly characterise the essential spectrum. For a certain class of connectivity functions in the neural field model, we are able to compute its spectral properties and the first Lyapunov coefficient of a Hopf bifurcation. By examining a numerical example, we find that the addition of diffusion suppresses non-synchronised steady-states while favouring synchronised oscillatory modes.


Introduction
In the study of neurological disease, non-invasive imaging techniques are often used to get an understanding of the structure and functioning of the brain on intermediate scales.As they give a course-grained view of the neuronal activity, mean-field models are a natural fit to describe the observed dynamics.In this paper we will use a neural field model with gap-junctions, electrical connections between neurons, which are thought to be related to observed synchronisation of neural tissue in Parkinson's disease [Schwab et al., 2014a,b].We will study the effect of gap junctions on the dynamics of the model.We will mainly focus on the stability of steady-states, periodic oscillations and the bifurcations which lead to a qualitative change in behaviour.
To properly address the difference in time-scales between gap-junctions and synaptic connections, we use a neural field with transmission delays for the synaptic connections.This leads to a complicated model which is infinite-dimensional and has spatially-distributed delays.The dynamical theory for such models is not readily available.In this paper, we address the analytic problems which arise from these abstract delay differential equations.
We use the sun-star calculus as the basic functional analytic tool to cast the equation in the variation-of-constants form.We exploit the results in Janssens [2019] that allow the linear part of the equation, without the delays, to be unbounded, as is the case for the diffusion operator.

Background
Neural field models try to bridge the gap between single neurons models, e.g.[Hodgkin and Huxley, 1952], and whole brain models, e.g.[Sanz Leon et al., 2013], by modeling the qualitative behaviour of large groups of neurons.In the seminal work of Wilson andCowan [1972, 1973], they modelled two populations of excitatory and inhibitory neurons and analysed the dynamical properties of the resulting model.A neural field uses spatial and temporal averaging of the membrane voltage of a population of neurons.The synaptic connections are modelled by a convolution of a connectivity kernel and a nonlinear activation function.
These models have been simplified by Amari [1977] by combining the excitatory and inhibitory populations into a single population and made more realistic by Nunez [1974] by including transmission delays.These delays arise from the finite propagation speed of action potentials across an axon and the delay due to dendritic integration.Faye and Faugeras [2010] investigated stability properties of stationary solutions using methods from functional analysis.Later, van Gils et al. [2013] put the neural fields into the functional analytic sun-star framework to compute normal form coefficients for Hopf and double Hopf bifurcations.Dijkstra et al. [2015] expanded their analysis to the Pitchfork-Hopf bifurcation.We will build on these papers by introducing gap-junctions into the neural field model and studying the resulting dynamics.
Gap-junctions are electrical connections between neurons, which directly exchange ions through a connexin-protein.This is in contrast to synaptic connections, where a potential is induced across the synapse by neurotransmitters.These gap-junctions are thought to be related to Parkinson's disease by synchronising neurons in the globus pallidus [Schwab et al., 2014a,b].Gap-junctions can be modelled as a simple diffusion process [Coombes et al., 2014].There have been some attempts to incorporate gap-junctions into networks of coupled neurons [Amitai et al., 2002;Laing, 2015], but to our knowledge not yet within a proper neural field model.
The first approach to develop a geometric theory for delay equations along the lines of ODEs was given by Hale [1971] who used formal adjoint operators.To avoid the formal adjoint, Diekmann et al. [1995] used the sun-star calculus for delay differential equations.This theory uses the duality pairing between a subspace of the dual space X and its dual X * , pronounced as "X sun star".This allows for a formulation of the variation-of-constants formula on the space X * from which the center manifold and bifurcation theorems can be derived.Recently, Janssens [2019] has begun expanding the sun-star calculus to the cases where the linear part, which contains no delays, is an unbounded operator.This case, with the additional assumption that this operator defines a compact semigroup, was examined earlier by Wu [2012] using the formal adjoint method to study reaction-diffusion systems with delays.One of the main benefits of the sun-star framework, instead of the formal adjoint method, is that it uses only standard methods from functional analysis and that it is possible to establish a perturbation theory for both compact and non-compact semigroups.For our neural field this means that we can use the same framework to study the model with and without diffusion.
In this paper we will build on the work of Janssens [2019] and prove the necessary theorems to use the sun-star calculus to study our neural field model with diffusion.We will then derive the spectrum and resolvent of a neural field with delays, diffusion and a connectivity kernel of a sum of exponentials.Finally we will compute the first Lyapunov coefficient of a Hopf-bifurcation and verify our results by simulating the full neural field numerically.The numerical analysis of such models was first attempted by R. Bellingacci in his Master Thesis at Utrecht University.

Modelling
In this section we will derive the neural field model with transmission delays and gap junctions.This is largely based on a derivation by Ermentrout and Cowan [1980].
We start with a collection of neurons i = 1, 2, 3, • • • and denote the (somatic) potential of neuron i at time t by u i (t) and its firing rate by f i (t).We assume there is a nonlinear dependence of f i on u i given by f i (t) = S i (u i (t)) We define Φ i,j (t) to be the postsynaptic potential appearing on postsynaptic cell i due to a single spike from presynaptic cell j.We assume a linear summation of the postsynaptic potentials, so the total potential received at the soma due to the synaptic connection between cell i and j can be modelled as where τ i,j is the delay due to the finite propagation speed of action potentials along an axon and other factors such as dendritic integration.We define Ψ i (t) to be the potential appearing in neuron i due to a gap-junction current I i,gap (t).The resulting model for u i becomes We can reduce this integral equation if we have a model for Φ and Ψ.For cell i, let us consider a passive membrane with a time constant 1/α i , a resistance R i and an injected postsynaptic current I i,j,syn (t) and similarly when a gap-junction current is injected If we now apply the Laplace transform L to equation (1), we get We assume that the synaptic dynamics are dominated by the time-scale of the membrane.This means we can reduce I i,j,syn (t) to w i,j δ(t), where δ is the Dirac-delta distribution and w i,j represents the strength of the synaptic connection, where a negative value corresponds to inhibition.Taking the inverse Laplace transform results in a system of differential equations 1 We want to model this network of cells by a neural field.Suppose we have a sequence of similar neurons i = 1, 2, • • • , M on the interval Ω = [−1, 1] and we model the gap-junctions as a simple resistor between adjacent neurons we arrive at the formula 1 α We will now take the limit as M → ∞ while scaling g by M 2 and w i,j by 1/M , to find our neural field model ∂u ∂t We haven't specified yet what happens with the gap-junctions at the boundary of our domain.It is natural to assume that no current leaks away at the boundaries, which corresponds to Neumann boundary conditions in the neural field: ∂u ∂x (t, ±1) = 0

Overview
The paper is divided into three parts, each of which can mostly be read independently.
In section 2, we construct the sun-star calculus for these equations and derive the variationof-constants formula.In particular we prove a novel characterisation for sun-reflexivity.Furthermore we consider linearisation, the corresponding spectrum and a normal form derivation for Hopf bifurcation of the nonlinear equations.In appendix A we elaborate on the case when the unbounded linear operator is the diffusion operator.We expect the reader to be familiar with the basics of the sun-star framework in Diekmann et al. [1995].
In section 3 we derive formulas for the eigenvalues and eigenvectors for a neural field with a connectivity defined by a sum of exponentials.We also explicitly construct the solution to the resolvent problem for this class of neural field models.
The section 4 we do a numerical study for a neural field model with specific parameter values.We compute the first Lyapunov coefficient for the Hopf bifurcation and investigate how it is influenced by the diffusion term.We will also investigate the emergence of periodic behaviour using numerical simulations of the neural field.
(ADDE).Next we study the linearisation and obtain results on the spectrum.Finally we construct a method for computing the first Lyapunov coefficient for a Hopf bifurcation of nonlinear equations.We build on the theory developed in Janssens [2019], who considers a class of abstract delay differential equations with a possibly unbounded linear part.
Consider two Banach spaces Y and X = C([−h, 0]; Y ) over R or C. Let S be a strongly continuous semigroup on Y with its generator B and let G : X → Y be a (nonlinear) globally Lipschitz-continuous operator.Note that the assumption that B is compact is not necessary, in contrast to what is assumed by Wu [2012].
Here u t ∈ X, where u t (θ) = u(t + θ) for t ≥ 0 and θ ∈ [−h, 0].On X we consider the strongly continuous semigroup T 0 defined by Here ϕ ∈ X, t ≥ 0 and θ ∈ [−h, 0].This semigroup is also called the shift-semigroup and is related to the problem for The solution of problem ( 6) is then given by v t := T 0 (t)ϕ.
Lemma 1. [Engel and Nagel, 1999, Theorem VI.6.1]The generator A 0 of the shiftsemigroup T 0 is given by We will interpret the (ADDE) as problem (6) with some nonlinear perturbation G : X → Y and use a variation-of-constants formula in X to obtain results about the perturbed problem, such as normal form coefficients for local bifurcations.As G maps X into Y , we would like to embed Y in a natural way into X.A naive approach would be to use a delta-function as an embedding.However, this embedding is not bounded, so the domain of A 0 would not be preserved under perturbation.This is indeed the case, as the rule for extending a function beyond its original domain, i.e. φ(0) = Bϕ(0), is incorporated in D(A 0 ).Hence adding a perturbation to the rule for extension changes the domain of the generator.A way out is to embed this problem into a larger space.A natural choice would be Y × X, where we have a continuous embedding : Y → Y × {0} and we can separate the extension and translation part of A 0 into Y × {0} and {0} × X respectively.
More formally we use the sun-star calculus as developed in Diekmann et al. [1995] to construct the space X * , which contains the space Y × X.We will first restrict the dual space X * to the sun space X , on which T * 0 is strongly continuous.Then taking the dual we obtain the dual space X * .It is convenient to present the relationship of the various spaces schematically in the following 'duality' diagram, see figure 1.

Characterisation of the sun-dual
We can represent X * , the dual space of X, as N BV ([0, h]; Y * ), the space of functions f : [0, h] → Y * of bounded variation on [0, h], normalised such that f (0) = 0 and f is right continuous on (0, h).This follows from a generalisation of the Riesz Representation Theorem for Y -valued functions proven by Gowurin [1936].The (complex valued) duality pairing between X and X * is given by the Riemann-Stieltjes integral, for ϕ ∈ X and f ∈ X * f, ϕ := and the action of A * 0 is given by A * 0 f = B * y * χ 0 +g, where χ 0 = 1 (0,h] , i.e the characteristic function of (0, h]. Proof.Let f ∈ D(A * 0 ) and ϕ ∈ D(A 0 ).Without loss of generality we can write A * 0 f = cχ 0 + g, where c ∈ Y * and g ∈ N BV ([0, h]; Y * ) and g(h) = 0. Using the integration by parts formulas for Riemann-Stieltjes integrals [Janssens, 2019, Proposition A.15, A.18, A.19] we obtain We will now want to use some limiting argument.However, the Riemann-Stieltjes integral lacks good convergence properties.In the scalar case, we could interpret this integral as a Lebesque-Stieltjes integral, which has better convergence properties.For a general Banach space Y and continuous integrands, the equivalent would be the Bartle integral [Singer, 1957].The Bartle integral has an equivalent theorem to the Lebesque Dominated Convergence Theorem.For uniformly bounded, pointwise converging sequences we can interchange the limit and the integral.[Bartle, 1956, Theorem 6].
For some 0 < s < t ≤ h and y ∈ Y , we may choose ( φn ) n∈N as a uniformly bounded sequence in X such that φn (0) = ϕ n (0) = 0 and it converges pointwise to y1 [−t,−s] , i.e. the characteristic function of [−t, −s].We then substitute ϕ for ϕ n in (9) Taking the limit as n → ∞, we get using the dominated convergence of the Bartle integral that Since y was arbitrary, we infer that Letting s ↓ 0, we obtain for t ∈ [0, h] where y * = lim s↓0 f (s).Now we substitute this formula for f into f, A 0 ϕ and use integration by parts and the fact that φ(0) = Bϕ(0) to find that We compare this to equation ( 9) Since ϕ(0) can be chosen arbitrary, y * , Bϕ(0) = c, ϕ(0) implies that c ∈ D(B * ) and c = B * y * .
Conversely let f be of the form in (8), then by the above computations we find that We can characterise the sun-dual X as the subspace of X * where T * 0 is strongly continuous or equivalently X = D(A * 0 ), where the closure is with respect to the norm on X * .The proof of the following theorem includes showing that T * 0 is strongly continuous on the set E (10), that D(A * 0 ) ⊆ E, and that E is closed.Theorem 3. [Janssens, 2019, Theorem 1 and Remark 4] The space X , the sun-dual of X with respect to T 0 , is given by the set is an isometric isomorphism.
From now on we will identify X with Y × L 1 ([0, h]; Y * ).The corresponding duality pairing between X and X is then given by ϕ , ϕ := y , ϕ(0) Now we can describe the action of T 0 and A 0 .
Definition 4. The strongly continuous semigroup T 1 on L 1 ([0, h]; Y * ) is defined as Theorem 5. [Janssens, 2019, Theorem 1] For the action of T 0 on X we have where the integral is the weak * Lebesque integral with values in Y .
Theorem 6.For the sun-dual of A 0 on X we have that Proof.By definition we have that and ιA 0 ϕ = A * 0 ιϕ .Let ϕ = (y , g) ∈ X such that ιϕ ∈ D(A * 0 ).Recall that the embedding ι is given by ( 11) . This is the case when B * y + g(0+) ∈ Y and we can write g as We conclude that A 0 ϕ = (B * y + g(0), ġ).

Characterisation of the sun-star space
We can represent X * , the dual of X , as Y * × (L 1 ([0, h]; Y * ) * , where Y * is the dual of Y .In general (L 1 ([0, h]; Y * ) * cannot be identified with L ∞ ([−h, 0]; Y * * ).However, the latter space can be embedded into the former.
Theorem 7. [Kreuter, 2015, Proposition 2.18 and 2.22] There exists an isometric embedding of The canonical embedding j : X → X * is defined as jϕ, ϕ = ϕ , ϕ .The continuous embedding : Y → X * is defined as = (j Y y, 0), where j Y is the canonical embedding of Y into Y * .It is possible to find an explicit representation of j.
Lemma 9.For ϕ ∈ X, jϕ = (j Y ϕ(0), ϕ).Moreover, j is a continuous embedding and Hence jϕ = (j Y ϕ(0), φ).The other statements are generally known to hold for the canonical embedding of X into X * , cf.Diekmann et al. [1995, Appendix II, Cor. 3.16, Prop. 3.17] As we don't have an explicit norm or measure on (L 1 ([0, h]; Y * ) * we cannot say anything in general about A * 0 .However, it is possible to find a representation of

In this case the action of
Then by Lemma 45 and Theorem 6, i.e. g(h) = 0, we can rewrite (17) as Taking g ≡ 0 we get that y * , B * y = γ, y for all y ∈ Y such that B * y ∈ Y by Theorem 6. Hence y ∈ D(B ), which implies that et al., 1986, Corollary 4.2], we find that y * ∈ Y * * .Furthermore by [Clément et al., 1986, Theorem 4.3] we have for all y * ∈ D(B * ) and y ∈ D(B * ).
Then again using Lemma 45 we get that for any (y , g) In this case, the action of A * 0 is given by A * 0 jϕ = (B * j Y ϕ(0), φ).Proof.This follows immediately from Theorem 10 and Lemma 9 Note that for A * 0 the rule for extension, φ(0) = Bϕ(0), is no longer included in the domain of A * 0 , but is represented in the action of A * 0 , which resolves the problem with A 0 stated at the beginning of this section.
The previous theorem allows us to formulate an equivalence between the sun-reflexivity of X, i.e.X = j(X) and the ordinary reflexivity of Theorem 12. X is sun-reflexive with respect to T 0 if and only if Y is reflexive.
Proof.Suppose that Y is reflexive.Then by Theorem 7 and Lemma 8, X * can be represented as Y * ×L ∞ ([−h, 0]; Y ) and hence the full domain of A * 0 is given by Theorem 10 = j Y (Y ).Next we note that C 1 functions are dense in the continuous functions and C 0 is closed with respect to the L ∞ -norm.Hence we conclude that Taking the norm closure of both sides, we conclude that proper subset of X , so X is not sun-reflexive.

Variation-of-constants formulation
As the space X * solves the problems mentioned in the beginning of this section.We can formulate a variation-of-constants formula for the (ADDE) as an abstract integral equation (AIE).
As the integrand of (AIE) takes values in X * , the integral is taken to be a weak * integral.It is possible to show that the integral maps to the range of j(X) and hence the (AIE) is well-defined.
The Banach Fixed Point Theorem in combination with the bound in ( 21) gives the existence of a unique global solution of (AIE).
Corollary 14. [Janssens, 2019, Corollary 9] Let G : X → Y be globally Lipschitz continuous.For every initial condition ϕ ∈ X there exists a unique solution v ∈ C(R + , X) such that u t = v(t) satisfies (AIE) for all t ≥ 0.
We would like to show that this unique solution of the (AIE) can be translated over to a (classical) solution of the (ADDE).However, this is in general not the case when B is unbounded.Therefore we recall a weaker solution concept from Wu [2012].
for all t ≥ 0 and u satisfies the (ADDE) Note that Definition 15 is quite restrictive as only specific initial conditions ϕ ∈ X are admissible.There is the following correspondence between classical and mild solutions of (ADDE) Lemma 17. [Wu, 2012, Theorem 2.1.4]A classical solution of (ADDE) is also a mild solution of (ADDE) Conversely when G has a globally Lipschitz continuous Fréchet derivative and is also a classical solution of (ADDE).
Note that Theorem 25 below implies that the conditions on G are equivalent to ϕ ∈ D(A).
It is possible to construct a one-to-one correspondence between solutions of (AIE) and mild solutions of (ADDE).
1. Suppose that u is a mild solution of (ADDE).Define v : R + → X by Then v is a solution of (AIE).
2. Suppose that v is a solution of (AIE).Define u Then u is a mild solution of (ADDE).
Corollary 19.Suppose G is a globally Lipschitz operator and it has a globally Lipschitz Fréchet derivative then for all , there exists a unique classical solution of the (ADDE).

Linearisation
We want to investigate the behaviour near a fixed point.We will show that for the linearised problem, we can perturb the semigroup T 0 with generator A 0 to a semigroup T with generator A. In the next section we will investigate the spectral properties of A.
Linearising equation (ADDE) near a fixed point u, which we take without loss of generality to be u ≡ 0, results in the linear problem (LINP).
As with the general nonlinear problem we can define an abstract integral equation.
where L := DG(0).Then due to Lemma 13 and Corollary 14 we can define the strongly continuous semigroup T (t)ϕ := u t when DG(0) is globally Lipschitz.
Lemma 20. [Janssens, 2019, Theorem 19] Let DG(0) be globally Lipschitz continuous, then there exists a unique strongly continuous semigroup T on X such that for all ϕ ∈ X and for all t ≥ 0.
The strongly continuous semigroup T has a generator A. We want to establish how the perturbed generator A relates to the original generator A 0 , which can be done using the sun-star framework.A technical detail which we need to check is that the sun dual space X is the same with respect to T and T 0 .
Lemma 21. [Janssens, 2019, Proposition 20] X is also the maximal subspace of strong continuity of the adjoint semigroup T * on X * .The adjoint generator A * is given by and the generator A of the T is given by Finally X is also the maximal subspace of strong continuity of the sun-star semigroup T .
One could think that we could extend this argument and show that D(A * ) = D(A * 0 ) and A * = A * 0 + Lj −1 .However, this is not the case when we lack sun-reflexivity, i.e.X = j(X).We can circumvent these problems by restricting the domain to j(X) Lemma 22. [Janssens, 2019, Proposition 22] It holds that and A * = A * 0 + Lj −1 on this subspace.We can extend the Corollary 11 for A * 0 to A * , which will be needed for the computation of normal form coefficients.
Corollary 23.For ϕ ∈ X the following statements are equivalent In this case, the action of A * is given by A * jϕ = (B * j Y ϕ(0) + j Y DG(0)ϕ, φ).
Proof.The statement on the domain follows immediately from Lemma 22 and Corollary 11.Furthermore we have that We are now able to state the result which relates A to A 0 .
Theorem 24.[Janssens, 2019, Corollary 23] For the generator A of the semigroup T we have that We can cast ( 27) in a form which can also be found in Engel and Nagel [1999, Theorem VI.6.1] by using Corollary 11.
Theorem 25.For the generator A of the semigroup T we have that ) and ϕ has an a.e.derivative φ ∈ L ∞ ([−h, 0]; Y ).Furthermore we have that be the strongly continuous semigroup generated by B .By Diekmann et al. [1995, Appendix II Proposition 3.17] we have that Finally, for the action of A we derive

Spectral properties
In this section we state some results on the spectrum of the operator A, notably its essential spectrum and a method for computing its eigenvalues.
For an operator A on X the resolvent set ρ(A) is the set of all z ∈ C such that the operator z − A has a bounded inverse.The resolvent operator R(z, A) : can be decomposed into the point spectrum σ p (A) and the essential spectrum σ ess (A).
Conversely, suppose that z ∈ ρ(A 0 ) en let y ∈ Y .Then the constant function ψ(θ As z is not an eigenvalue of A 0 , by the above reasoning it is not an eigenvalue of B and hence z − B is injective. So we conclude that σ(A) = σ(B) and σ ess (A 0 ) = σ ess (B).
If DG(0) is compact then we can make inferences on the essential spectrum of A from the spectrum of A 0 .
Proof.We will proof this by working in the dual space.This is possible as σ ess (A) = σ ess (A * ), which is a consequence of the properties of Fredholm operators.[Katō, 1995, Theorem IV.5.14]On X * , A * = A * 0 + L * due to Lemma 21.As is bounded, L = DG(0) is compact and so is its adjoint L * due to Schauder's theorem, [Katō, 1995, Theorem III.4.10].Hence A * is a compact perturbation of A * 0 .One of the defining properties of Weyl's essential spectrum is that it is invariant under compact perturbations, cf.[Katō, 1995, Theorem IV.5.35].

So we conclude that
For computation the eigenvalues we follow Engel and Nagel [1999].We introduce the family of operators K z : Y → Y , H z : X → X and W z : X → Y parametrized by z ∈ C, defined as . Using these we can define the characteristic operator Now we formulate the main theorem of this section which allows us to reduce the computation of the eigenvalues and eigenvectors in X to a computation on Y Theorem 28.[Engel and Nagel, 1999, Proposition VI.6.7]For every z ∈ C, ϕ ∈ R(z−A) if and only if ∆(z)q = W z ϕ has a solution q ∈ D(B).Moreover z ∈ ρ(A) if and only if this q is unique.In that case the resolvent is given by where θ ∈ [−h, 0] and ψ ∈ X.Finally, ψ ∈ D(A) is an eigenvector corresponding to λ ∈ σ p (A) if and only if ψ(θ) = e λθ q, where q ∈ D(B) is non-trivial and satisfies ∆(λ)q = 0

Hopf-bifurcation
We are interested in the nonlinear behaviour of (ADDE).In this section we develop techniques to compute the first Lyapunov coefficient of an Hopf bifurcation assuming that a sufficiently smooth center manifold exists.These techniques can be extended to other local bifurcations, but we will not address those here.In this section, we will follow the methods from van Gils et al. [2013].
Suppose that σ(A) contains a pair of simple purely imaginary eigenvalues λ = ±iω with ω > 0 and no other eigenvalues on the imaginary axis.Let ψ ∈ X be the corresponding eigenvector of A and ψ ∈ X be the corresponding eigenvector of A respectively, We normalise these vectors such that ψ , ψ = 1 (33) The center subspace X 0 is spanned by the basis Ψ = {ψ, ψ} of eigenvectors corresponding to the critical eigenvalues of A.Here ψ denotes the complex conjugate of ψ.If ζ ∈ X 0 then ζ = zψ + z ψ for some z ∈ C.
We assume that there exists a locally invariant smooth critical center manifold W c loc ⊂ X, which is tangent to X 0 .The critical center manifold has the formal expansion Due to Theorem 18 the (ADDE) and (AIE) formulations are equivalent.By weak * differentiation of (AIE) and exploiting the finite dimensionality of W c loc , one can show that a solution v ∈ C(R + ; X), v(t) = u t , of (AIE) satisfies the abstract ODE Where the nonlinearity R : X → Y is given by Let ζ(t) = z(t)ψ + z(t) ψ be the projection of v(t) onto the center subspace X 0 .The function z(t) satisfies a complex ODE which is smoothly equivalent to the Poincaré normal form ż = iωz + c 1 z|z| 2 + O(|z| 4 ) (37) where z, c 1 ∈ C. In polar coordinates, z = re iθ , this is orbitally equivalent to where l 1 is the first Lyapunov coefficient determined by the formula It is well known, e.g.Kuznetsov [2013], that in generic unfoldings of (38), l 1 < 0 implies that the bifurcation is supercritical and that a stable limit cycle exists near one of the branches.On the other hand, l 1 > 0 implies that the bifurcation is subcritical and that an unstable limit cycle exists near one of the branches.
The critical center manifold W c loc has the expansion (34) and due to the time-invariance of W c loc we have If we differentiate both sides with respect to time and use the abstract ODE (35) for the left-hand side, we obtain the homological equation We can substitute the expansion of the nonlinearity (36), the normal form (37), and the expansion of the critical center manifold (34) into the homological equation ( 41) to derive the normal form coefficients.If we equate coefficients of the corresponding powers of z and z we obtain the following equations Lemma 29 (Fredholm solvability).Let z / ∈ σ ess (A).Then z − A : D(A ) → X has closed range.In particular Proof.From the definition of the essential spectrum, R(z − A) is closed [Katō, 1995, Section IV.5.1], and R(z − A * ) is also closed by Banach's Closed Range Theorem [Katō, 1995, Theorem IV.5.13] Due to Banach's Closed Range Theorem, ϕ * is a solution of We now return to equations (42).As {0, 2iω} ⊂ ρ(A) = ρ(A ) we can use the resolvent of A * to solve the first two equations.However, iω ∈ σ(A) so for the last equation of ( 42) we need to use the theorem above.The corresponding eigenspace N (A * − λ) is spanned by ψ , so we can compute for the normal form coefficient by We are not yet able to compute the normal form coefficient explicitly as we don't have an explicit representation of ψ or a representation of the resolvent of A * .However, we resolve this by using spectral projections.
Let P and P * be the spectral projections on X and X * corresponding to some eigenvalue λ, respectively.Then P * ϕ * = νjψ for some ν ∈ C and ϕ * , ψ = ϕ * , P ψ = P * ϕ * , ψ = ν jψ, ψ = ν Hence we seek to determine ν.From the Dunford integral representation it follows that where C λ is a sufficiently small open disk centered at λ and ∂C λ its boundary.The element on the left in the pairing ( 44) is of the form ϕ * = y, y ∈ Y .In this case we can reduce R(z, A * )ϕ * to ∆ −1 (z)y due to the following theorem.
Now given that we can compute the resolvent ∆ −1 (z) and the Fréchet derivatives of G, we have a method to compute the first Lyapunov coefficient l 1 = 1 ω Re c 1 and the normal form coefficients h 20 and h 11 .
Due to Theorem 28 we have that λ is an eigenvalue and ψ an eigenvector if and only if ψ(θ) = qe λθ and q ∈ D(B) satisfies characteristic equation (CE).
Where in this case K z : Y → Y is a parametrized family of operators for z ∈ C defined as where c j (z) := S (0)αη j e −τ 0 z = 0 and k j (z) := µ j + z The case without diffusion, i.e. d = 0, has already been extensively studied in van Gils et al.
[2013] and Dijkstra et al. [2015].So in this section we will develop formula's for the eigenvalues, eigenvectors and resolvent with nontrivial diffusion, i.e. d > 0.
For the following section we adopt the notational convention that bold-faced variables correspond to vectors a = (a 1 • • • a n ) T where its length is clear from the context

Eigenvalues
So we are looking for non-trivial solutions q ∈ D(B) of As this is a mixed differential-integral equation, it is in general hard to solve.We will use the method of Dijkstra et al. [2015] to convert (CE) into a differential equation (ODE), which we can solve.Then substituting the general solution of (ODE) back into (CE) yields appropriate conditions on q.This is possible due to the following observations.
Proof.As q ∈ C 2 (Ω) and the range of K z is contained in C 3 (Ω) we have that Bq ∈ C 2 (Ω), which means that q ∈ C 4 (Ω).By induction, we conclude that q ∈ C ∞ (Ω).
Differentiating the kernel functions in the (CE) in the distributional sense yields for j ∈ So we define the differential operator L z j for j ∈ {1, • • • , N }.
For this operator L j we have that for Proof.We have that z ∈ L if and only if P z (k j (z)) = 0 for some j ∈ {1, • • • , N }.
∈ L we can rewrite P z (ρ m ) as We can divide out the product to conclude that for m ∈ {1, Next we find formula's for K z j cosh(ρ m (z)x) and K z j sinh(ρ m (z)x).To compute these integrals we split the interval [−1, 1] into the intervals [−1, x] and [x, 1].On these intervals e −k|x−x | is an C 1 function in x so we can compute the following anti-derivatives for these smooth branches.
Using these anti-derivatives, we can evaluate the integrals K z j cosh(ρ m (z)x) and K z j sinh(ρ m (z)x).For clarity we omit the dependence on z in the remainder of this section.
Now we are ready to substitute the general solution q of (ODE), ( 51), back into (CE).
Due to the characteristic equation ( 53) the first line in equation ( 55) vanishes.When z / ∈ L, cosh(k j x) and sinh(k j x) for j ∈ {1, • • • , N } are linearly independent.Hence the second line vanishes if and only if S z,even a = S z,odd b = 0, where matrices S z,even and S z,odd are defined as As q ∈ D(B), we also need to take the boundary conditions into account as To satisfy the boundary conditions, we augment the matrices S z,even and S z,odd as follows: Now we have square matrices S z,even , S z,odd ∈ C (N +1)×(N +1) .There exists a non-trivial solution q ∈ D(B) of (CE) if and only if det(S z,even ) = 0 or det(S z,odd ) = 0.
When det(S λ,even ) = 0, the corresponding eigenvector ψ ∈ X is given by Where a is a vector in the nullspace of S λ,even .
When det(S λ,odd ) = 0, the corresponding eigenvector ψ ∈ X is given by Where b is a vector in the nullspace of S λ,odd .
Proof.Let q ∈ D(B) be a solution of (CE) for some λ ∈ C. Then by Theorem 32 q ∈ C ∞ so it is also a solution of (ODE).
tions and the fact that (z − B)R(z, B)y = y, we have that We have that the above equation vanishes when all the terms within square brackets vanish.The term (64a) vanishes naturally due to characteristic equation in (53) as z / ∈ L.
As R(z, B) maps into D(B), the boundary conditions q (±1) = 0 reduces to We can split equation ( 65) into 3 sufficient equations Note that the equations (66b) and (66c) are equivalent to If we combine the equations (67) with terms in square brackets in (64c) and (64d) we get the matrix equations.

∂ ∂x
We see that in equation ( 69a) the sum should be constant.Using equation (66a) we see that this constant is zero.
We can rewrite these equations by introducing some matrices.We define the diagonal matrices Ĉ, Ŝ ∈ C(Ω, Proof. As σ p (A) and σ p (B) contain only isolated eigenvalues and ρ m (z) and det(P z (k i,j (z))) are analytic in z, the set S contains only isolated values.Hence such a C λ exists.Suppose λ is an even eigenvalue.As S ∩ C λ = ∅ and σ(A) ∩ C λ = {λ}, we have that the ∆ −1 (z)y is given by Theorem 36 for z ∈ C λ .We observe that all components of the resolvent are analytic for all z ∈ C λ expect for the constants of integration a c (z).This analyticity simplifies ( 81) to The reasoning for odd eigenvalues is similar.

Numerical Results
In this section we will examine a specific numerical example.We will compute eigenvalues and the first Lyapunov coefficient for a Hopf bifurcation and investigate the effect of varying the diffusion parameter d.
For J we choose the following substraction of two exponentials, as in Dijkstra et al. [2015] J This connectivity is a model of a population of excitatory neurons acting on a short distance combined with a population of inhibitory neurons acting on a longer distance, see figure 4 For the activation function S we choose the sigmoidal function As S is an odd function, S (0) = 0 and hence that D 2 G(0) ≡ 0. This simplifies the computation of first Lyapunov coefficient l 1 of (47) to We can compute this integral using Theorem 37 with y =1 2 D 3 G(0)(ψ, ψ, ψ).
We fix the following values for parameters α = 1 and τ 0 = 3 4 and use γ as the bifurcation parameter.We want to compare two cases: without diffusion, i.e. d = 0, and with diffusion, i.e. d > 0.
As one might already have observed, the diffusion has little effect on the Hopf bifurcation.
We observe more generally that the eigenvalues which are off the real axis are barely effected by the introduction of diffusion, while the eigenvalues on the real axis become more negative, see figure 3. 1 A possible explanation is that the eigenvector corresponding to the eigenvalue on the imaginary axis has very little spatial curvature, see figure 4. As diffusion penalizes curvature, its effect on this eigenvector would be small.1 without and with diffusion respectively.Note that with diffusion the eigenvector satisfies the boundary conditions at x = 1 and x = −1, while this is not the case without diffusion.

Discretisation
To obtain an approximate solution of (ADDE) we discretise the spatial domain Ω into an equidistant grid of n x points, x 1 , . . ., x n x , with a width of δ = 2 n x −1 .We discretise the integral operator G using the trapezoidal rule and the diffusion operator B using a central difference method and a reflection across the boundary for the boundary conditions.This results in a second order spatial discretisation.The discretisation of the (ADDE) for n ∈ {1, • • • , n x } and t ∈ R + becomes a set of delay equations (DDE), cf.[Faye and Faugeras, 2010] Here ξ m is defined as Now we are left with a set of n x ordinary delay differential equations which we solve with a standard DDE-solver.Note that the (DDE) is very similar to the discrete model (3) from which the (ADDE) is derived.Only the terms at the boundary are different due to the second order discretisation.

Simulations
We will now perform some simulations around the Hopf-bifurcation with diffusion.We set n x = 50 and take as initial conditions an odd function and an even function, For figure 5 we took γ = 3 and for figure 6 γ = 4.
For γ = 3, both initial conditions (92) converge to the trivial equilibrium.The odd initial condition converges monotonously to the trivial equilibrium, while the even initial condition converges to the trivial equilibrium in an oscillatory manner.For γ = 4, there are (at least) two non-trivial stable states.The odd initial condition converges to some non-trivial equilibrium and the even initial condition converges to some limit cycle, which is due to the Hopf bifurcation.This is similar to the results of [Dijkstra et al., 2015], where the non-trivial equilibrium arises from a Pitchfork bifurcation.The bi-stability is also exemplified in the eigenvalues, see figure 7, as we have a positive real eigenvalue and a pair of complex eigenvalues with a positive real component.
We have seen that increasing the value of d, decreases the eigenvalues on the real axis.This would imply that the non-trivial equilibrium becomes unstable or disappears, probably through a Pitchfork bifurcation.Indeed when we use the initial condition and compare the dynamics for d = 0.2 and d = 0.5 in figure 8.The initial condition converges to a non-trivial equilibrium when d = 0.2, but it converges to a limit cycle when d = 0.5.

Discussion
We have demonstrated that the sun-star calculus for delay equations is a natural setting to study the neural field models both with and without diffusion.We have proved the the necessary theorems to construct the sun-star calculus for abstract delay differential equations.In particular we proved a novel characterisation for sun-reflexivity in Theorem 12.The sun-star calculus provides a variation-of-constants formulation for the nonlinear problem and produces results on the spectral properties of the system, notably the essential spectrum.Following Diekmann et al. [1995], the variation-of-constants formulation can be used to make a center manifold reduction, however some technical details still need to be worked out.Given a center manifold reduction, we are able to compute the first Lyapunov coefficient for the Hopf bifurcation.This procedure can quite easily be extended to normal coefficients of other local bifurcations.
For certain specific connectivity functions we have derived analytical conditions when λ is an eigenvalue for a neural field with a connectivity function which is a sum of exponentials.
We have also constructed the corresponding eigenvectors and the resolvent.Numerical results show that the diffusion term does not cause oscillations to arise due to a Hopfbifurcation.However, in a bi-stable system, consisting of a steady-state and an oscillatory mode, increasing the diffusion leads to a system where only the oscillations are stable.
Gap junctions, modelled by the diffusion term in our neural field, are thought to be linked to synchronisation in Parkinson's disease, cf.[Schwab et al., 2014a].Further research could be undertaken to see whether the effects can be observed in a Neural Field Model with physiological values for the parameters.
We used a neural field model with a connectivity function which is a sum of exponentials.This connectivity function is commonly used to aggregate the effect of multiple different types of cells, e.g.excitatory and inhibitory neurons.However, introducing a diffusion term into this model leads to gap junctions between similar and different populations of neurons of the same strength.This may not be physiologically feasible.A way to circumvent this is to use a neural field model with multiple populations.In such a model, it is possible to introduce only gap junctions between neurons of the same population.
We have studied a neural field on a 1-dimensional closed domain.However, when modelling the neuronal activity in the cortex, it is common to use 2-dimensional domains, cf.[Coombes et al., 2014].For a rectangular domain, characterising the spectrum is still an open problem.On a spherical domain, Visser et al. [2017] have characterised the spectrum for a neural field with transmission delays and have computed normal form coefficients of Hopf and double Hopf-bifurcation.It seems possible to extend the analysis of that paper to include a diffusion term into that neural field model.Due to the general nature of the theoretical results of section 2, these results, including the sun-star framework, the variation of constants formulation and the essential spectrum, also hold for neural field models on arbitrary domains.

A.2 Sun-star calculus
We will now develop the sun-star calculus for the diffusion operator B. We can take d = 1 and α = 0 for this section, without loss of generality, as the sun-star calculus is invariant with respect to bounded perturbations of the generator of the semi-group.
As a consequence of the Riesz representation theorem, Y * can be represented as N BV (Ω), the functions of normalised bounded variation.The corresponding norm on N BV (Ω) is the total variation norm and the duality pairing is given by the Riemann-Stieltjes integral: We will now try to find a representation for B * .
Theorem 41.The dual space Y * can be represented as N BV (Ω).Furthermore, y * ∈ D(B * ) if and only if for x ∈ (−1, 1] (w ) (s) ds Note that the sun-dual Y is almost the same as in Diekmann et al. [1995, Theorem II.5.2],where it is taken with respect to the first derivative with the condition ẏ(0) = 0.However, in that case there was an extra condition in Y that functions g ∈ L 1 could be extended be zero for θ ≥ h.In our case with diffusion we have a fixed domain on which the diffusion takes place, so this condition is not present.

Figure 1 :
Figure 1: A schematic representation of the various spaces in sun-star calculus following Janssens [2019] to find an explicit representation of the dual operator A * 0 and its corresponding domain D(A * 0 ).Theorem 2. The domain of A * 0 is given by D(A * 0 ) := {f ∈ N BV ([0, h]; Y * )| there exists y * ∈ D(B * ) and g ∈ N BV ([0, h]; Y * ) ϕ has an a.e.derivative} We use that X is the closure of D(A * 0 ) with respect to the norm on X * .First the closure of D(B * ) with respect to the Y * -norm results in the space Y .By van Neerven [1990, Corollary 2.5] reflexivity implies sun-reflexivity, Y 42)They all have the form (z− A * )ϕ * = ψ * (43)Here z ∈ C and ψ * ∈ X * are given.When z ∈ ρ(A) then (43) has a unique solution.However, if z ∈ σ(A), then a solution ϕ * doesn't necessarily exist for all ψ * .The following lemma, which is equivalent to vanGils et al. [2013, Lemma 33], provides a condition for solvability.

Figure 3 :Figure 4 :
Figure 3: The eigenvalues of A at parameter values in table 1 of the Hopf bifurcation without and with diffusion respectively.

Figure 7 :
Figure 7: The eigenvalues of A for γ = 4 and d = 0.2