Travelling waves in models of neural tissue: from localised structures to periodic waves

We consider travelling waves (fronts, pulses and periodics) in spatially extended one dimensional neural field models. We demonstrate for an excitatory field with linear adaptation that, in addition to an expected stable pulse solution, a stable anti-pulse can exist. Varying the adaptation strength we unravel the organizing centers of the bifurcation diagram for fronts and pulses, with a mixture of exact analysis for a Heaviside firing rate function and novel numerical schemes otherwise. These schemes, for non-local models with space-dependent delays, further allow for the construction and continuation of periodic waves. We use them to construct the dispersion curve – wave speed as a function of period – and find that they can be oscillatory and multi-valued, suggesting bistability of periodic waves. A kinematic theory predicts the onset of wave instabilities at stationary points in the dispersion curve, leading to period doubling behaviour, and is confirmed with direct numerical simulations. We end with a discussion of how the construction of dispersion curves may allow a useful classification scheme of neural field models for epileptic waves.


Background
The analysis of waves in models arising from the study of the nervous system has a long tradition.The seminal example is the development and numerical analysis of the model for action potential propagation in an excitable axonal fibre by Hodgkin and Huxley [1], and see [2] for an excellent review.This has been followed by more rigorous mathematical analysis using tools from geometric singular perturbation theory for the existence of pulses (homoclinics) and wave trains (periodics) [3], as well as the development of a stability theorem [4].However, the detailed properties of waves in detailed biophysical models is often best pursued with a mixture of both mathematical and numerical analysis.This is nicely exemplified by the work of Miller and Rinzel [5], who numerically construct the dispersion relation (between speed and period) for steadily propagating periodic wave trains, and then use a mathematical kinematic theory, which applies the dispersion relation, instantaneously, to individual pulses to predict how interspike time intervals can be distorted during propagation.Going beyond the single neuron scenario traveling waves in neurobiology are often studied at the tissue level, using neuroimaging methods such as electroencephalography (EEG).For example these waves can take the form of spindle waves seen at the onset of sleep [6], the propagation of synchronous discharge during an epileptic seizure [7] and waves of excitation associated with sensory processing [8].Such waves are a consequence of synaptic interactions and the intrinsic behaviour of local neuronal circuitry.The propagation speed of these waves is of the order cm s −1 , an order of magnitude slower than that of action potential propagation along axons.The class of computational models that are believed to support synaptic waves differ radically from classic models of waves in excitable systems.Most importantly, synaptic interactions are non-local (in space), involve communication (space-dependent) delays (arising from the finite propagation velocity of an action potential) and distributed delays (arising from neurotransmitter release and dendritic processing).In contrast models of axonal fibres are typically based on local partial differential equation models (that track the opening and closing of ion channels and model voltage transport along a fibre with a diffusive process).A common way to model the spatiotemporal evolution of coarse grained variables such as synaptic or firing rate activity in populations of neurons is through the use of a neural field model (in which space is continuous).Despite a growing body of analysis on synaptic waves in neural field models, focusing mainly on fronts and pulses (and see [9,10] for further discussion) surprisingly relatively little analysis has been performed on periodic waves.Rather in most neural field models of EEG only spatially homogeneous solutions and their instabilities to pattern states have been analysed e.g.[11][12][13].Interestingly in some of these models such instabilities have been shown to give rise to periodic travelling waves, e.g. in the work of Curtu and Ermentrout on neural field models with adaptation [14] and that of Venkov et al. [15] in models with short-range inhibition, long range excitation and axonal delays, though not subsequently analysed in the full nonlinear regime (far from bifurcation).This begs the question as to whether, like many models for action potential wave trains in axons, it is useful to analyse and classify neural field models in terms of their dispersion curves and the large period limit where one recovers the description of a travelling pulse.This is precisely the question we address in this paper, by considering a minimal one dimensional single population neural field model.Importantly the techniques we develop are readily applicable to more general neural field models, including those with multiple populations in higher dimensions and with various forms of feedback, such as spike frequency adaptation.
In Section "Neural tissue models" we discuss some simple neural field models with a focus on those with purely excitatory connectivity and linear adaptation as these are perhaps the simplest ones that support travelling waves.Direct simulations are used to show that not only do they manifest travelling pulses and periodic wave trains, but also a localised wave of decreased activity that we term an anti-pulse.The existence and stability of both pulses and anti-pulses is determined in Section "Travelling fronts and pulses" for the case of an idealised Heaviside firing rate function.Novel numerical techniques for handling more general firing rate functions are developed in Section "Numerical techniques for analysing neural fields with sigmoidal firing rates".Not only do these allow us to show that the analytical results obtained for Heaviside firing rates are similar to those for steep sigmoidal functions, they also open the door for the numerical study of periodic waves and their dispersion curves.These are computed in Section "Dispersion curves of wavetrains" and discussed in the context of a kinematic theory that has previously been used so effectively for excitable models of axons to predict and organise irregular

Neural tissue models
Neural field models describe the coarse-grained activity of neuronal populations, and make no attempt to model the detailed behaviour of single neurons.Rather they model the statistical properties of tissue activity that is more directly relevant to EEG.Due to the long-range axonal connectivities in cortex neural field models are typically formulated as non-local integral or integro-differential equations.For a recent review see [16].Compared to the analysis of local models (of ordinary or partial differential equation type) their analysis (mathematical and numerical) is not as thoroughly developed.Fortunately, in certain physiologically realistic regimes equivalent partial differential equation (PDE) models can be formulated.In particular this then allows the use of powerful techniques from nonlinear PDE theory to be brought to bear, especially the numerical analysis of travelling waves.
To illustrate this possibility we now formulate a simple continuum neural field model in one spatial dimension and present some simulations of travelling wave behaviour.
Neural field models of cortex typically assume a density of neurons at a point with inputs that arise from the delayed and weighted contribution of activity at other points in the tissue.Let us introduce the synaptic activity at a point x in the tissue at time t and call this u(x, t).The output from this activity will be taken to be a population firing rate of the form f (u), for some sigmoidal shape f with steepness parameter β > 0 and threshold θ.A simple neural field model then takes the symbolic form where ψ is the drive to the synapse and captures information about both anatomical connectivity patterns and the distribution of axonal delays, while the temporal differential operator Q describes synaptic filtering.For example to capture a post synaptic potential with a bi-exponential response of the form e −α 1 t − e −α 2 t for t > 0 (up to a normalisation factor α 1 α 2 /(α 2 − α 1 )) we would use: To generate a normalised α-function response, α 2 te −αt , we would simply set α 1 = α = α 2 , and to model a synapse with a simple exponential response αe −αt we would let α 2 → ∞ and set α 1 = α.The drive to the synapse takes the integral form where sets the domain of the tissue, w describes the anatomical connectivity and v represents the velocity of an action potential.Models of the type (2) have a long history in the field of neural tissue modelling and have been particularly useful in understanding the generation of EEG rhythms, as in the work of Nunez [17], Jirsa and Haken [18], and Liley et al. [11].Moreover, they are naturally extended to include adaptive processes, such as spike frequency adaptation [19] (i.e. the addition of a current that activates in the presence of high activity).These are known to favour the generation of travelling pulses over travelling fronts.In this case we would modify the drive according to ψ → ψ − a, and choose a simple dynamic model for a such as Here τ and κ are constants describing the timescale separation between activity and adaptation and how strong the (linear) adaptation becomes, respectively.More sophisticated models of metabolic processes whose combined effect is to modulate neuronal response could also be considered, though for the purposes of this paper we shall work with (5).
Figure 1 shows simulations of model (2) for two different values of the adaptation strength κ.We choose 1/v = 0, θ = 1/4, τ = 7 and a steep sigmoidal activation function (1) with β = 42.We use an equidistant spatial discretisation with periodic boundary conditions and compute ψ using fast Fourier transforms (FFT) (as for 1/v = 0 it has a spatial convolution structure).With the same initial conditions quite different behaviour emerges.For κ = 1.1 a travelling pulse appears, while for κ = 0.9, we observe a moving localised area with lower activity and high background activity.We will refer to this solution as a travelling anti-pulse.A solution of this type has previously been observed in [20], though not exhaustively analysed.In the next section we will do so, for a Heaviside firing rate function (β → ∞), and show that 2θ(κ + 1) = 1 determines a critical value where the dynamics change from pulse to anti-pulse.

Travelling fronts and pulses
The simulation in Figure 1 indicates that a stable travelling anti-pulse solution exists.We will consider the existence and stability for a Heaviside firing rate function, making use of an Evans function approach.This is a powerful tool for the stability analysis of nonlinear waves on unbounded domains.It was originally formulated by Evans [4] in the context of a stability theorem about excitable nerve axon equations of Hodgkin-Huxley type.The extension to integral models is far more recent, see [21] for a discussion.For neural field models with axonal delays it has previously been noted these do not typically induce any change of wave stability [22,23] (though it will affect the shape and speed of a wave).Since our simulations indicate similar effects for the anti-pulse, we will initially consider the case 1/v = 0 for simplicity.The discussion below is adapted from that presented in [22,23], to which we refer the reader for further details.

Existence of the anti-pulse
To show existence and stability it is convenient to first rewrite the integro-differential equation ( 2) in integral form.We define η(t) = e −t and η a = e −t/τ /τ with η(t) = η a (t) = 0 if t < 0. The model with adaptation can then be written in the integral form where * denotes temporal convolution (e.g. The construction of travelling waves, e.g.fronts and pulses, proceeds by introducing the comoving frame ξ = x − ct with c > 0. Travelling waves are then stationary solutions u(x, t) = u(ξ ) where with Introducing the Fourier transform, equation ( 6) can be arranged to the following form Exploiting the pole structure of ( 11), an inverse Fourier transform yields A travelling wave solution can then be written succinctly as For f (u) = H(u − θ), a Heaviside step function, the pulse solution can be found by specifying where the solution is above the threshold θ.An anti-pulse has a profile u(ξ ) above threshold everywhere except for ξ = x − ct ∈ [− , 0].So the profile is formally given by where we note that the width and speed c are yet unknown.These can be found by imposing the crossing conditions u(− ) = u(0) = θ.This yields the following two equations When we take the limit → ∞, we see that ( 15) and ( 16) which recovers a result in [24] for right-moving inactivating fronts (left-moving waves in their case).For completeness we note that the travelling pulse satisfies These are determined analogously as for the anti-pulse but with the requirement that the profile is above threshold only for ξ ∈ (− , 0) [24].The speed of activating fronts (moving to the right) can be obtained from ( 18) in the limit → ∞.The above equations apply to right-moving waves, i.e. c > 0.

Stability of the anti-pulse
We will now construct the Evans function.First, we linearise (6) around the travelling solution ū(ξ ) by writing u(ξ , t) = ū(ξ ) + z(ξ , t).In particular, we will look for exponentially decaying/growing separable solutions of the form z(ξ , t) = z(ξ )e λt .Collecting O(z) terms yields the eigenvalue equation For the Heaviside firing rate function we have f (u(r)) = δ(u(r) − θ)/|u (r)|, which appears inside the integral only and hence poses no difficulties.Since we have crossing points at ξ = 0 and ξ = − , we get z(ξ where The eigenvalue problem should be self-consistent at ξ = 0 and ξ = − giving the system of equations A nontrivial solution for z exists if E(λ) := det(A(λ) − I) = 0.The function E(λ) is the Evans function for the travelling anti-pulse.Note that the same construction applies to the travelling pulse.Solving ( 15) and ( 16) we find that c(κ) has a turning point, so that we have a fast and a slow branch, denoted by c + and c − respectively.Figure 2 shows the roots of the Evans function along the pulse and anti-pulse branches.There is a trivial eigenvalue λ = 0 due to translation invariance.The nontrivial eigenvalue is positive for the slow (anti-) pulse c − , and negative for the fast branch c + .Hence, the latter is stable.

A combined view on fronts and pulses
From equations ( 15) and ( 16), it is straightforward to determine the speed as a function of the adaptation strength κ, see Figure 3.For completeness, we also plot the speeds of the (right-moving) travelling anti-front, front and pulse.What is apparent is that all curves meet at the point C given by Here we have a codimension 2 heteroclinic cycle bifurcation.We have checked using simulations that varying κ for Gaussian connectivity induces a similar transition between pulse and anti-pulse behaviour.The unfolding of this bifurcation involves two heteroclinic and two homoclinic bifurcation curves.It implies that the anti-pulse exists generically.We have also plotted the profiles for κ = 0.65 and κ = 0.75 with τ = 7 and θ = 0.3, see Figure 4.This shows the difference between the slower and faster solutions.The afterovershoot in activity is due to diminished adaptation during the anti-pulse.When τ is decreased, two things happen.First, the profile develops a slow and strongly damped oscillation as λ ± become complex.Second, the fast branch of stable pulse solutions disappears, see Figure 3 (right) as the adaptation is fast and pulls down all activity quickly.This happens when the pulse and anti-pulse branch have become tangent to the inactivating front and front branches for → ∞, respectively, i.e. they have interchanged their position in the (c, κ) plane.

Numerical techniques for analysing neural fields with sigmoidal firing rates
Despite the elegant analysis that can be done for a Heaviside firing rate function this is a special case and the techniques developed above for existence and stability do not generalise.A natural alternative approach is to develop numerical schemes for investigating the persistence of fronts and pulses for more general sigmoidal firing rate functions, such as given by equation ( 1), as well as the inclusion of axonal delays where 1/v is finite.In this section we develop two complementary numerical schemes, both of which are easily embedded within a continuation framework to track travelling wave solution branches in parameter space.The first exploits shapes of synaptic connectivity for which an equivalent PDE model can be obtained, as in the Jirsa-Haken-Nunez brain wave equation [18], and the second is a more general approach that has no restriction on the choice of kernel shape.In both cases the challenge is to compute the synaptic drive (for a travelling wave) given by: in a computationally efficient and accurate fashion.

Numerical continuation I: Equivalent PDE
Although neural field models are inherently nonlocal, for special choices of connectivity an equivalent PDE model can be formulated, as described in [9].This has also been used for simulation and continuation in, for instance, [11,[25][26][27][28].In this case travelling wave and pulse solutions are readily found with standard continuation techniques.In illustration of this point consider an exponentially decaying connectivity w(x) = exp(−|x|)/2.
Then on an infinite domain, it can be shown [29] that ψ(x, t) given by (4) satisfies the PDE where Equation ( 25) is often referred to as the Jirsa-Haken-Nunez brain wave equation.We then look for waves u(x, t) = u(ξ ) with ξ = x − ct of the PDE system given by ( 2), ( 5) and (25).Under the replacement ∂ x → d ξ and ∂ t → −cd ξ , we obtain, for our explicit choices, the 4th order ordinary differential equation system: It is easily checked that the model has either one or three equilibria with (u, ψ, φ, a) = (u, f (u), 0, κu), where u is a solution of u = f (u)/(1 + κ).Orbits connecting one or two equilibria are homoclinic and heteroclinic solutions corresponding to travelling pulses and fronts, respectively.These are readily analysed using a numerical toolbox such as MATCONT [30].However, the equivalent PDE approach is limited to certain special choices of synaptic connectivity possessing a rational Fourier transform structure [25].

Numerical continuation II: Integral model
Many modelling studies employ Gaussian connectivity rather than a decaying exponential for which an equivalent (finite order) PDE cannot be derived.In such cases, we spatially discretise the neural field and do continuation of waves in the moving frame ξ = x − ct directly as in [31] (for a neural field model with refractoriness) and in [32] for stationary patterns (in a model without axonal delays).Here we will work with the choice w(x) = σ exp(−(σ x/2) 2 )/2 √ π, with σ > 0. We discretize u and a on an equidistant mesh and compute derivatives u ξ , a ξ with central finite differences.Typically, we use 2 N , N = 11 mesh points which gave stable numerical results.Next we use periodic boundary conditions and an integral phase condition, i.e. (u − u 0 ) udξ = 0 where u 0 is some reference solution, e.g. the previously computed orbit.Next we compute the spatial convolution using FFT.When we consider axonal delays, i.e. 0 < v < ∞, the convolution structure is not obvious.The crucial point is that the drive (for stationary solutions in the traveling wave frame) is a sum of convolutions that can each be computed efficiently using FFTs.To see this we break the integration of ψ(ξ ) into two parts, one over the left-half of the real line and the other over the right.Introducing a ± = c/v ± 1 allows us to write (24) as and in the last line we see that each term separately has a convolution structure.We conclude that ψ(ξ ) (with finite v) can be cast into a suitable form by an asymmetric scaling of the connectivity function.In the limit when the wavespeed equals the transmission speed, the scaled connectivity function w(z/a − ) approaches a delta distribution.Initial data is taken from a stable wave found with simulations.We then compute periodic solutions of the integral model with standard pseudo-arclength continuation with free parameters the speed c and the size of the spatial domain T.

Numerical results
The results from the analysis with a Heaviside activation function (for fronts and pulses) are recovered for smooth, but steep sigmoidal activation functions.For illustration we choose β = 42 and keep τ = 7 and θ = 0.30, see Figure 5.The diagram with the heteroclinic cycle persists including the values of κ where the fronts and pulse emanate from c = 0.The curve for the speed of the anti-pulse also exhibits the same limiting behaviour for parameters near the heteroclinic cycle bifurcation.However, the homoclinic bifurcation curve corresponding to the anti-pulse now emanates from a Zero-Hopf (ZH) bifurcation point.Its unfolding involves homoclinic and heteroclinic bifurcation curves that oscillate towards the ZH-point generically.The corresponding waves have small amplitude and are unstable.Note also that the critical value of κ corresponding to a saddle-node bifurcation in which the higher steady state disappears, is lower compared to the case of Heaviside activation function.Additionally, there is a Hopf bifurcation curve http://www.epjnonlinearbiomedphys.com/content/2/1/3 0.4 0.8 1. starting from (κ, c) ≈ (1/τ , 0) where the homoclinic curve corresponding to pulse solutions emanates as well.Fixing system parameters, we can follow a periodic orbit from the Hopf curve by continuation in c and the period T. Continuing further, the periodic orbits accumulate on a homoclinic curve corresponding to the anti-pulse or pulse, which ever occurs first.There is also a homoclinic orbit to the middle steady state (left purple curve in Figure 5), but this steady state is unstable itself and consequently the travelling wave as well.We do not consider this branch further.
Decreasing the slope parameter to β = 9.5 we obtain the bifurcation diagram shown in Figure 6 that is much more complex than for higher values of β.In essence, we find that the range of κ where three steady states coexist, shrinks.Consequently, also the range over which fronts exist shrinks, but slightly more.The travelling pulse, on the other hand, exists over a much wider range of adaptation strengths, even where only the lower steady state exists.Next travelling waves have become complex as the appearance of the LPC curve indicates and as discussed in the Section "Dispersion curves of wavetrains".The diagram contains an additional saddle-node curve and the Hopf curve no longer starts from c = 0, but turns.First we note that the heteroclinic cycle still persists.Also another Zero-Hopf point at the left saddle-node bifurcation curve is present together with a similar structure of homoclinic curves.One has a high speed and is associated to an unstable wave as explained above.The other homoclinic curve has several turns until at κ ≈ 0.5425 it approaches the heteroclinic cycle tangentially to a heteroclinic curve.These are the anti-pulse and inactivating front, respectively.The heteroclinic curve is very similar to the one computed for β = 42, but the branch terminates when it meets the Hopf curve.The other heteroclinic curve representing the activating front is even split into two branches.Emanating from c = 0 it meets the right saddle-node curve tangentially and then exists for 0.2848 < κ < 0.7647 between two Hopf curves for higher values of c.The upper branch is tangent to the homoclinic curve representing the travelling pulse.This homoclinic curve exists to the right of the saddle-node bifurcation where it turns and then undulates towards the left Zero-Hopf point.There is another homoclinic curve emanating from the left Zero-Hopf point with high speeds but similar to the other purple branches, it corresponds to unstable solutions.The stability of the front and pulse solutions is similar as for higher values of β.It is remarkable that many of the features obtained from the Heaviside analysis persist for not so steep sigmoidal firing rate functions (though the overall bifurcation structure can be more complex).
Travelling waves can be found by starting from a Hopf bifurcation and then to continue periodic orbits in two parameters, e.g. the period T plus the speed c or the strength κ.When the period goes to infinity, the wave approaches one of the homoclinic curves discussed above.In between, in the (c, T)-diagram the curve of periodic orbits makes several turns.These limit point of cycles (LPC) mark the increasing complexity of travelling wave solutions, as we will discuss below.Here the primary LPC's are indicated by a thick blue solid line.

Dispersion curves of wavetrains
As briefly discussed in the Background, dispersion curves for periodic wavetrains in axonal models have proved very useful for understanding the behaviour of more irregular wavetrains using a kinematic theory.Thus it is first useful to set the scene for dispersion curves in neural field models by reviewing this approach for a simple excitable fibre model with FitzHugh-Nagumo (FHN) dynamics.

Dispersion curves for a FHN model
We write the FHN model as where Despite its simplicity this model has a rich dynamical structure with a variety of travelling waves and pulses, as discussed and analysed in [33][34][35][36].The dispersion curve for periodic orbits (stationary in a co-moving frame) relates the speed of a wave to its period, giving c = c(T).In the limit of a large period one recovers the homoclinic orbit describing a solitary travelling pulse.http://www.epjnonlinearbiomedphys.com/content/2/1/3 For the FHN model the dispersion curve has a slow and fast branch, and a minimum period largely determined by the refractory time-scale.If the shape of the periodic wave (over one period) is pulse-like then one may invoke a kinematic theory to determine the stability of solution branches, as well as describe the evolution of inter-pulse intervals in more exotic wavetrain patterns (that can bifurcate from periodic orbits).For a review of this approach we refer the reader to [37].Importantly periodic waves are predicted to be stable if c (T) > 0, with wave bifurcations (found to be period-doubling bifurcations for the FHN model) occurring at stationary points in the dispersion curve where c (T) = 0.It is well known in the context of travelling waves in PDEs that the eigenvalues of the steady state corresponding to the homoclinic solution can strongly affect the shape of the dispersion curve [33].The stable and unstable eigenvalues closest to the imaginary axis are called leading.When the leading eigenvalues are real, the dispersion curve approaches the asymptotic wavespeed (namely that of the solitary pulse) monotonically.When some are complex, the dispersion curve displays oscillations if the complex eigenvalues are the closest to the imaginary axis.Moreover, in the FHN model the saddle-focus boundary can be determined analytically.One can determine the dispersion curve for parameter values on both sides of the saddle-focus boundary, as in [35,38], see Figure 7. Fixing b, the origin has three real eigenvalues for larger values of a, while for lower values of a the origin is a saddle-focus.This is reflected in the dispersion curve, where for large values of a the dispersion curve is monotone and the origin is a "simple" saddle as the eigenvalues show.For smaller values of a the origin is a saddle-focus and the dispersion curve is oscillatory.Note that the monotonicity and oscillations are also visible in the activity profile.

Dispersion curves for neural field models
We turn now to the construction of dispersion curves in neural field models for pulse-like periodic waves.Although we may do this explicitly for the case of a Heaviside firing rate function (extending the technique for fronts and pulses used in Section "Travelling fronts and pulses") we prefer instead to focus on the numerical construction of dispersion curves using the approaches developed in Section "Numerical techniques for analysing neural fields with sigmoidal firing rates" .Moreover, as for the FHN model, we pinpoint the saddle-focus boundary as it predicts where one finds complex dispersion curves.Based on this we show several bifurcation diagrams where we collect the primary backbone of waves and pulses.From these diagrams it is straightforward to understand the shape of dispersion curves for various values of the system parameters.
We start with high values of β keeping 1/v = 0 for now.We choose the adaptation strength κ such that there is both a fast and a slow pulse and compute the dispersion curve by solving the travelling wave equation for u(ξ ) for fixed speed c with the periodicity constraint that u(ξ ) = u(ξ + T), to determine the dispersion curve c = c(T), see Figure 8.In this case the dispersion curve is monotone, which reflects that the leading eigenvalues of the steady state are real.Also the profile approaches the steady state without oscillating.The faster branch is stable according to the kinematic theory, which agrees with direct simulations.
For lower values of β, the dispersion curve is qualitatively different.As the wavenumber grows, both the fast and slow branch oscillate toward the asymptotic value, see Figure 9.The oscillations are due to complex eigenvalues of the steady state.Also the profile spirals towards the steady state.The stability of the wavetrains is much more delicate as with every other turn in the dispersion curve, there is a change of stability.We have found near every extremum a period-doubled branch c 2 (T) which also oscillates as T → ∞ leading to period-doubled branches c 4 (T).This structure repeats itself, but we show only the main additional branches.Hence there exist complex stable wavetrains where one or more pulses are grouped together, see Figure 9.
The transition from saddle to (wild) saddle-focus occurs for decreasing values of β or τ or increasing κ and involves several codimension 2 homoclinic bifurcations.At the transition two scenarios are possible.First, two real leading eigenvalues can collide and become complex as shown in Figure 10.For τ = 7 we found that the complex eigenvalues were already closer to the imaginary axis for β < 15.78.This is one of the two Belyakov bifurcations [39] also underlying the transition in the FHN model [35,38].For most parameter choices, especially high τ this situation is typical.In the second scenario, for lower values of τ the transition from monotone to oscillatory dispersion curves occurs as the complex leading unstable eigenvalues are closer to the imaginary axis than the real stable one as illustrated in Figure 10 at β = 22.73 for fixed τ = 4.4.For values of β below this critical value we then find oscillatory dispersion curves.We remark that at β = 26.35 the unstable eigenspace of the equilibrium of the homoclinic orbit is three-dimensional, which, to the best of our knowledge, is the first observation of this codimension 2 homoclinic The spatial period has been divided by the wavenumber to illustrate that the additional branches start near the extrema.Inset shows that the green branch also emerges from an extremum.Right: The spatial arrangement of travelling pulses can be much more complex with one or more pulses grouped together.Colors correspond to the branches of the dispersion curve.
bifurcation.A detailed investigation is outside the scope of this paper, but we found no bifurcations emerging from this point.Note that when β is decreased further the complex pair of eigenvalues approaches the imaginary axis suggesting a Shilnikov-Hopf scenario, similarly as [35].
Moving parameters further, increasing κ and/or decreasing β, the dispersion curve oscillates more and more.At a certain moment, the fast branch is multi-valued.Here for given T the dispersion curve has two points with c (T) > 0 on the faster branch suggesting bistability of waves.Indeed, simulations show that there are two stable waves with a different profile and velocity for the same wavelength, see Figure 11 (left).Changing parameters even more, the dispersion curve breaks up into two disconnected sets creating gaps, see Figure 11 (right).Here the primary T-periodic wave does not exist for certain values of T marking the birth of ever more complex wavetrains.Eventually, for lower values of β the part with the homoclinics disappears in a saddle-node bifurcation.Here the system is oscillatory and the dispersion relation behaves as c = γ T for some constant γ .
Let us finally also mention the effects of model choices on dispersion curves for neural field models.Using the new numerical continuation scheme we have developed, we have computed the dispersion curves for a few model variants.We do this for parameters such that waves in all variants exist.First, and most importantly, we have not observed any qualitative changes when we change to Gaussian connectivity or alpha synapses or include delays.The quantitative aspects are still considerable though.The use of Gaussian connectivity can increase the speed, while transmission delays slow them down, see Figure 12 (left).Their combined effect shows that Gaussian connectivity still increases the speed.In addition, we have observed that as κ is increased, the waves first disappear for exponential connectivity and then still exist for Gaussian connectivity.For an alpha synapse the wavespeed decreases roughly twofold.Also the wavelength at which the first extremum in the dispersion curve occurs, is halved, see Figure 12 (right).

Discussion and conclusions
In this paper we have developed the numerical tools to explore and continue travelling wave solutions for non-local neural field models with space-dependent axonal delays.Moreover, we have validated our approach against the analytically tractable case of a Heaviside firing rate and shown how bifurcation diagrams of this special case are modified as one moves toward more physiologically realistic shallower sigmoidal firing rate shapes.Interestingly we have shown that as well as pulses and fronts expected for excitatory networks with inhibitory feedback, also anti-pulses are a robust travelling wave solution.Moreover, the bifurcation diagram for travelling localised states, i.e. fronts, pulses, and anti-pulses, is organised around a co-dimension 2 heteroclinic cycle bifurcation.Our main result, however, has been the numerical construction of dispersion curves.We have shown that they offer similar insight into the behaviour and instability of periodic travelling waves as originally found in their application to excitable reaction diffusion systems.Namely, that the eigen-spectrum of fixed points for a homoclinic orbit corresponding to a pulse can strongly affect the shape of the dispersion curve by complex leading eigenvalues arising through a saddle to saddle-focus transition leading to oscillations in the dispersion curve), as the LPC curve in Figure 6 indicates.Not only can this lead to wave instabilities at stationary points which predicted by a kinematic theory, are period doubling bifurcations, it also leads to multi-stability and waves of super-normal speed, as seen in Figures 5,6 and 11.Such super-normal speeds are greater than that of the isolated pulse and can be read off from the dispersion curve for large values of the period.Furthermore it is possible to see gaps in dispersion curve for the case when β is small.
Importantly our work opens up a novel way to classify the behaviours (and the consequences for spatio-temporal wave propagation) for a broad spectrum of neural field models that have been used in the modelling of EEG, such as the Liley model [11], and in particular those for epilepsy, such as in the work of Marten et al. [40] and Goodfellow et al. [13].Namely, we expect the main similarities, or differences, between these models to be captured through a comparison of their dispersion curves.We conjecture that this may clarify why coherent oscillations, as opposed to travelling waves, are found to be favoured in the model presented in [40].This, and related work on how to sculpt the shape of dispersion curves through the detailed form of the neural field model using multiple populations, higher-order synapse models, distributions of axonal speeds and so on, will be reported upon elsewhere.

Figure 4
Figure 4 The profile of the front and pulse solutions all travelling to the right with θ = 0.3, τ = 7 and Heaviside activation function.The dashed solutions are slower and unstable.(Top) Travelling activating front (Left) and pulse (Right) for κ = 0.75.(Bottom) Travelling inactivating front (Left) and anti-pulse (Right) for κ = 0.65.The solid curve ( = 9.346, c = 0.4858) corresponds to a stable anti-pulse and the dashed ( = 2.394, c = 0.243) to an unstable solution.

Figure 7
Figure 7 Dispersion curves for (left) a = 0.15 and (right) a = 0.01 and b = 0.0025.Stars/dots indicate eigenvalues of the origin for higher/lower asymptotic wavespeeds, respectively.The dispersion curve and the profile of u have an oscillatory tail for a = .01.