Travelling waves in a neural field model with refractoriness

At one level of abstraction neural tissue can be regarded as a medium for turning local synaptic activity into output signals that propagate over large distances via axons to generate further synaptic activity that can cause reverberant activity in networks that possess a mixture of excitatory and inhibitory connections. This output is often taken to be a firing rate, and the mathematical form for the evolution equation of activity depends upon a spatial convolution of this rate with a fixed anatomical connectivity pattern. Such formulations often neglect the metabolic processes that would ultimately limit synaptic activity. Here we reinstate such a process, in the spirit of an original prescription by Wilson and Cowan (Biophys J 12:1–24, 1972), using a term that multiplies the usual spatial convolution with a moving time average of local activity over some refractory time-scale. This modulation can substantially affect network behaviour, and in particular give rise to periodic travelling waves in a purely excitatory network (with exponentially decaying anatomical connectivity), which in the absence of refractoriness would only support travelling fronts. We construct these solutions numerically as stationary periodic solutions in a co-moving frame (of both an equivalent delay differential model as well as the original delay integro-differential model). Continuation methods are used to obtain the dispersion curve for periodic travelling waves (speed as a function of period), and found to be reminiscent of those for spatially extended models of excitable tissue. A kinematic analysis (based on the dispersion curve) predicts the onset of wave instabilities, which are confirmed numerically.


Introduction
The continuum approximation of neural activity can be traced back to work of Beurle (1956), who built a model describing the proportion of active neurons per unit time in a given volume of randomly connected nervous tissue. A major limitation of this very early neural field model is its neglect of refractoriness or any process to mimic the metabolic restrictions placed on maintaining repetitive activity. It was Cowan (1972), (1973) who first developed neural field models with some notion of refractoriness. At the same time they also emphasised the importance of modelling neural population in terms of an excitatory subpopulation and an inhibitory subpopulation. Indeed over the years many studies of the Wilson-Cowan excitatory-inhibitory model have been made, with applications to problems in neuroscience ranging from the generation of electroencephalogram rhythms through to visual hallucinations, and see (Coombes et al. 2003) for a review. However, many of these subsequent studies drop the refractory term and focus more on the role of excitatory-inhibitory interactions in generating neural dynamics. Perhaps one exception to this is the work of Curtu and Ermentrout (2001), who have shown that the original Wilson-Cowan model with refractoriness can drive oscillations even in the absence of inhibition. Their work was done for a point model which begs the question as to whether refractoriness alone can allow for periodic waves to be generated in a spatially extended excitatory network. This is an especially intriguing issue given that neural field models with some form of inhibition or negative feedback, such as spike frequency adaptation, have traditionally been invoked to explain wave behaviour in cortex, including fronts, pulses, target waves and spirals (Ermentrout and McLeod 1993;Pinto and Ermentrout 2001;Huang et al. 2004).
In this paper we reinstate the original refractory term of Wilson and Cowan in a minimal neural field model describing a single population in one spatial dimension. This model is briefly reviewed in Sect. 2. In Sect. 3 we present a linear stability analysis, as well as a weakly nonlinear analysis, of the homogeneous steady state that predicts the onset of periodic travelling wave patterns in a purely excitatory network. This is confirmed by direct numerical simulations that show periodic travelling waves with profiles that appear as either single or multiple spikes of activity. A novel numerical continuation scheme is developed to track solution properties in a co-moving frame (speed, period, and profile shape) as a function of physiologically important system parameters (such as refractory time-scale, strength of anatomical connectivity, and firing threshold). These are obtained after recognising that the original model can be reformulated as a delay-differential equation for an exponentially decaying choice of anatomical weight distribution. The delay is set by the time-scale of the refractory process. In Sect. 4 we numerically construct the wave speed as a function of the wave period, to obtain the so-called dispersion curve. Here we avoid special case choices of the weight distribution and develop a numerical scheme that can handle the original delayed integro-differential model. The dispersion curve for an excitatory network with an exponentially decaying weight distribution is shown to have a shape reminiscent of that seen in the study of nonlinear reaction-diffusion systems, and in particular those arising in the study of an axon or active dendrite (Miller and Rinzel 1981). Using the dispersion curve we further develop a kinematic model that allows predictions about non-regular spike trains to be made, including period-doubling scenarios subsequently confirmed by direct numerical simulations. In addition, we establish an example of a homoclinic orbit of chaotic saddle-focus type in an infinite-dimensional system. Finally in Sect. 5 we present a brief discussion of the work in this paper.

The Wilson-Cowan model with refractoriness
Wilson and Cowan considered the spatio-temporal evolution of the activity of synaptically interacting excitatory and inhibitory neural sub-populations (Wilson and Cowan 1972). A recent review of their model can be found in (Coombes et al. 2013;Bressloff 2012). A common reduction of their original model, and one often employed as a minimal model of cortex, takes the form of a scalar integro-differential equation: Here u = u(x, t) ∈ R is a temporal coarse-grained variable describing the proportion of neural cells firing per unit time at position x ∈ R at the instant t ∈ R + . The symbol ⊗ represents spatial convolution, the function w describes an effective anatomical connectivity or weight distribution and is a function of the distance between two points, and τ is the relaxation time-scale. The nonlinear function f describes the expected proportion of neurons receiving at least threshold excitation per unit time, and is often taken to have a sigmoidal form. In major contrast to the original Wilson-Cowan equations refractory terms are not included in this model. To reinstate such terms in (1) we follow (Wilson and Cowan 1972), and more recently (Curtu and Ermentrout 2001), and model the fraction of cells in their absolute refractory period R by Since only a fraction (1 − z) of cells can be activated and actually contribute to any firing activity the model (1) is modified to For convenience we rescale time t → Rt and define r = R/τ to obtain the model that we shall work with for the remainder of this paper: As a choice of firing rate we shall take the sigmoid with threshold θ and steepness parameter β. For the choice of weight distribution we shall consider symmetric normalised kernels such that w(x) = w(|x|) and ∞ −∞ w(y)dy = 1.

Analysis of waves
Direct numerical simulations of a purely excitatory network, see below, show the possibility of periodic travelling waves. This is particularly interesting because these are not typically found in neural field models with pure excitation, though they are often encountered in the presence of some form of negative feedback, such as may arise with the inclusion of an inhibitory sub-population or a form of spike frequency adaptation, as reviewed in (Ermentrout 1998). If these patterns arise via the instability of the homogeneous steady state, then they can be predicted using a classic Turing instability analysis. Their analysis beyond the point of instability can be pursued with a weakly nonlinear analysis, to develop a set of amplitude equations (typically in the form of coupled complex Ginzburg-Landau equations), as in (Curtu and Ermentrout 2004;Venkov et al. 2007). However, this is only relevant close to the bifurcation point, and it is much more informative to gain an insight into the fully nonlinear properties of waves using numerical analysis. This has been pursued at length for many excitable systems, and especially for single neuron models of the axon or active dendrite with single (Miller and Rinzel 1981;Röder et al. 2007) or multi-pulse (Evans et al. 1982;Feroe 1982;Hastings 1982;Kuznetsov 1994;Lord and Coombes 2002) periodic waves. However, the study of periodic travelling waves has largely been ignored in the neural field community, which is surprising since this can inform a kinematic analysis [elegantly reviewed in (Keener and Sneyd 1998)] to predict instabilities to more exotic classes of travelling wave solution. We build on a Turing analysis and develop precisely this approach below.

Linear stability analysis of homogeneous solutions
A homogeneous fixed point with u(x, t) = u is given by the solution of the nonlinear algebraic equation The model displays up to three different fixed points depending on β and θ , see Fig. 1, which we denote by u i with u 1 < u 2 < u 3 when three fixed points exist.
Linearising around u and considering perturbations of the form e ikx e λt gives a dispersion relation for the pair (λ, k) The Turing bifurcation point is defined by the smallest non-zero wave number k c that satisfies Re (λ(k c )) = 0. It is said to be static if Im (λ(k c )) = 0 and dynamic if Im (λ(k c )) ≡ ω c = 0. A static bifurcation may then be identified with the tangential intersection of ω = ω(ν) and ν = 0 at ω = 0. Similarly a dynamic bifurcation is identified with a tangential intersection with ω = 0. Beyond a dynamic instability one would expect the emergence of a periodic travelling wave of the form e i(k c x+ω c t) (with speed ω c /k c ). For the sigmoidal function (5) we have that f = β f (1 − f ). For an exponential kernel w(x) = Se −S|x| /2, with S > 0, then W (k) = (1 + (k/S) 2 ) −1 , which has a maximum of one at the origin. Hence for λ ∈ R we have that For large r we see that E(λ, 0) is a decreasing function of λ so that a fixed point will be stable (to a static instability with k c = 0) if lim λ→0 E(λ, 0) < 0, or equivalently, β > β s where For λ ∈ C it is natural to write λ = ν + iω and equate real and imaginary parts of (7) to obtain two equations for ν and ω, which we write in the form The simultaneous solution of these two equations gives the pair (ν(k), ω(k)). For λ = iω the above pair of equations reduces to Since there are singularities at ω = 2nπ , these equations define a series of parametric curves T i = T (u i ) defined in the regions ω ∈ (2nπ, 2(n + 1)π ) for n ∈ Z. We see that there will be a solution if and only if It is clear that (13) defines a quadratic function in u and it turns out that T 2 does not exist, only T 1 and T 3 . In Fig. 2 we show only the branches with ω < 2π as we did not find any with larger ω. Using the observation that sin(ω)/ω < 1 and (1−cos ω)/ω 2 < 1/2 we see that a solution for Hence a dynamic instability will occur before a static instability when β < β s and r is sufficiently large. To determine whether the static instability gives rise to a travelling or a standing wave it is useful to perform a weakly nonlinear analysis.

Weakly nonlinear analysis: amplitude equations
A characteristic feature of the dynamics of systems beyond an instability is the slow growth of the dominant eigenmode, giving rise to the notion of a separation of scales. This observation is key in deriving the so-called amplitude equations.
In this approach information about the short-term behaviour of the system is discarded in favour of a description on some appropriately identified slow time-scale. By Taylor-expansion of the dispersion curve near its maximum one expects the scalings Re λ ∼ r − r c , k − k c ∼ √ r − r c , close to bifurcation, where r is the bifurcation parameter. Since the eigenvectors at the point of instability are of the type Z L e i(ω c t+k c x) + Z R e i(ω c t−k c x) + cc, for r > r c emergent patterns are described by an infinite sum of unstable modes (in a continuous band) of the form e ν 0 (r −r c )t e i(ω c t+k c x) e ik 0 √ r −r c x . Let us denote r = r c + 2 δ where is arbitrary and δ is a measure of the distance from the bifurcation point. Then, for small we can separate the dynamics into fast eigen-oscillations e i(ω c t+k c x) , and slow modulations of the form e ν 0 2 t e ik 0 x . If we set as further independent variables T = 2 t for the modulation time-scale and X = x for the long-wavelength spatial scale (at which the interactions between excited nearby modes become important) we may write the weakly nonlinear solution as It is known from the standard theory (Hoyle 2006) that weakly nonlinear solutions will exist in the form of either travelling waves (TWs), Z L = 0 or Z R = 0, or standing waves (SWs), Z L = Z R . In the Appendix we show that (ignoring spatial variation) the amplitude equations take the form where Λ, Φ and Ψ are known functions of system parameters. A linear stability analysis of the above amplitude equations generates the conditions for selection between TWs or SWs. If Re Λ and Re Φ have opposite sign, then a TW exists and for Re Φ < 0 and Re (Φ − Ψ ) > 0 it is stable. If Re Λ and Re (Φ + Ψ ) have opposite sign, then a SW exists and for Re (Φ + Ψ ) < 0 and Re (Φ − Ψ ) < 0 it is stable. Evaluation of the coefficients Φ and Ψ (see Appendix), yields Re Φ > 0 and Re (Φ + Ψ ) > 0. Therefore, for typical parameter values the Turing instability is subcritical and travelling and standing waves are unstable. Also, along T 3 , waves exist for δ > 0 if r > r s ≈ 20.3 and for δ < 0 if r < r s . Still, these can be used to start to track waves numerically. Note that a similar weakly nonlinear analysis for waves in a neural field model without refractoriness though with adaptation has been performed in (Curtu and Ermentrout 2004), and for axonal delays in (Venkov et al. 2007).
In the next part we will consider waves and how these can grow beyond a dynamic Turing bifurcation.

Excitability and waves
From the Turing analysis above we expect to see travelling waves for sufficiently large r and suitable θ . A similar observation, based on the numerical simulation of a lattice model with nearest-neighbour coupling has previously been made by Curtu and Ermentrout (2001). These authors further point out that for some range of β values This can be re-written in delay-differential form aṡ We graph the nullcline of the u-equation z = 1 − u/ f (u) together with the steady states constraint u = z and two trajectories in Fig. 3. These show that if an initial displacement from the steady state is sufficiently large, then the trajectory makes a large excursion. In other words, the system is excitable.
Next we want to show travelling waves for the full spatially extended system defined by (4) using direct numerical simulations. We evolve the state as follows. We use an equidistant spatial discretisation with 2 11 mesh points with periodic boundary conditions. We compute the spatial convolution using Fourier transforms. The history integral z(x, t) = t t−1 u(x, s)ds is calculated using a trapezoidal rule with 100 points. This gives a system that can be simulated with matlab's dde23-solver. Supplying the history as u(x, τ ) = .05 + 0.7 exp(−80(x − 1 − 0.63τ ) 2 ) we obtain a travelling wave, see Fig. 4 (left). Other initial history can also lead to travelling waves as long as the amplitude at one spot is sufficient to cause excitation of neighbouring tissue and the initial spot decays due to refractoriness. The figure suggests that there is enough space to fit in a second moving pulse and this is indeed possible, see Fig. 4 (right). Note that the time from the one pulse to the next is different than from the previous pulse. The travelling waves shown in Fig. 4 have nearly the same velocities namely c ≈ 0.6302 for the travelling one pulse, and c ≈ 0.6309 for the two pulses. The trajectories near the steady state in Fig. 3 also go some way to explaining why the two pulses can move slightly faster. For two pulses in one periodic domain the next pulse arrives when the system is less refractory, i.e. z is lower, than with only one pulse. We study the dependence of the speed on the wavenumber below.
If the domain for the wave becomes infinite, the travelling wave approaches a pulse which is a homoclinic orbit. The linear stability of the steady state in the moving frame ξ = x + ct classifies the homoclinic orbit. All travelling waves profiles can be constructed as stationary profiles of (4) in the travelling frame ξ = x + ct, namely as solutions of the dynamical system: Focusing now on the homoclinic orbit we consider small perturbations u(ξ ) =ū+v(ξ ) to obtain This has solutions of the form v = exp(λξ ) where λ is a solution of the transcendental equation Here we impose the condition |Re (λ)| < S to ensure convergence of the integral over R in (18). Fortunately, this includes the imaginary axis allowing stability analysis. If S → ∞ we recover the formula derived in (Curtu and Ermentrout 2001) for the Hopf bifurcation of the point model. Solving for the (single) steady state we findū ≈ 0.0554. We insert the wavespeed c = 0.6303 as observed in the simulations and numerically solve the eigenvalue equation (19) for many random starting values near the origin in the complex plane, see Fig. 5. We find a single positive real eigenvalue λ 1 ≈ 8.1902 and the other eigenvalues are complex pairs with negative real part. The leading stable eigenvalues are λ 2,3 ≈ −5.8021 ± 3.8026i (and satisfy the constraint |Re (λ)| < S).
We conclude thatū is a saddle-focus with saddle-quantity λ 1 /Re (λ 2,3 ) > 1. This implies the existence of N-homoclinic loops for all N = 1, 2, 3, . . .. In particular, we may expect that travelling waves with 3 pulses form an isolated branch and have a saddle-node bifurcation (Gonchenko et al. 1997), even though the system is infinitedimensional.
It is an open challenge to rigorously prove this (and extend the result from the finite dimensional setting). However, it is likely that Lin's method can be applied in this case, along the lines considered in (Lin 1990) for analysing the bifurcation of a unique periodic orbit from a homoclinic orbit to an equilibrium in a DDE.

Numerical continuation of waves
Here we start with the derivation of an equivalent DDE in a co-moving frame for the travelling waves. This is a standard approach and normally would allow us to study periodic orbits that relate to waves in the original system. However, we found that for the available numerical tools to work we needed to modify the equations artificially. Therefore we computed the waves in an alternative and novel way. Rather than using a PDE approach we inserted the co-moving frame directly and used the discretisation of that system as described below.

A delay differential equation for waves
We write z(x, t) = t t−1 u(x, s)ds and transform (4) to a delayed partial differential equation (DPDE) using Fourier techniques (see (Coombes et al. 2003) for more details) to obtain Next we insert the wave ansatz ξ = t + x/c and obtain the following delay differential equation (DDE) One can make the following observations about this equivalent DDE. First, suppose we had inserted the more common ansatz ξ = x − ct for right-going waves. We would then have obtained an advanced instead of a delayed term. Also, this time rescaling gives a constant delay which is numerically more stable. Second, for simulations one can also compute the history integral from the DDE for z in (21). However, the history of z, must satisfy the constraint So, if we know the history of z we actually know u(τ ) for τ ∈ (−2, 0]. We use system (21) to determine periodic orbits with numerical continuation. The continuation problem automatically specifies enough history so that this does not pose a problem. Solving (21) using knut (Roose and Szalai 2007) we did not achieve convergence. This can be understood from the steady state problem of system (16). This is illdefined as any constant may be added to z giving a continuum of solutions as we lose the constant of integration. The integral constraint yields z =ū and hence maximally three steady states exist. Adding a small additional term ε(u − z) to z ξ in (16) and (21) with ε = 10 −6 incorporates the integral constraint. The continuation results for both system (16) and (21) were then numerically stable and in agreement with final profiles obtained from simulations.

Direct continuation of the integral equation
As the DDE-approach only works for the modified system, we develop a novel numerical scheme to track periodic solutions of the integral equation in a co-moving frame. The novelty is that we do not introduce auxilary variables as for the DDE, but compute the (convolution) integrals directly using fast Fourier transforms (FFT). Working with the non-local model (17) directly also allows us to treat a more general class of weight distributions and not only those of exponential form (though we do not pursue this here).
Similarly as for the simulations of (4) we use an equidistant spatial grid ξ with N mesh points for the interval [0, Δ]. We employ central finite differences for the temporal derivative. We have one convolution with the connectivity w. For our periodic solutions u, this convolution can be computed by taking the FFT of w and u, multiplying element-wise and then applying the inverse FFT. It is sufficient to take the FFT of w and not its periodic summation as the connectivity w decays sufficiently fast so that w(−Δ/2) ≈ 0. Hence the circular convolution theorem can be employed. We observe that the integral over [ξ − c, ξ] can be seen as another convolution of u with g = H (ξ − c)H (−ξ) with H the Heaviside step function. It is not strictly necessary, but we have the Fourier transforms of g and w analytically and can evaluate them immediately without transforming these with a FFT. For simplicity and convergence, we use N = 2 13 . Next we use the periodic boundary condition u 0 = u N +1 to eliminate one equation. Then to break the translational invariance we add the integral phase condition (u − u 0 )udξ = 0 where u 0 is some reference solution; here we take the previously computed point. This results in N + 1 equations for N + 1 unknowns. Then we use the pseudo-arclength condition v, (x − x 0 ) = h, where v is the tangent vector, h is the step-size and the brackets indicate the standard inner-product (Meijer et al. 2009). We add the spatial period Δ as an additional parameter.

Initial data for the continuation
The numerical continuation of travelling waves needs an initial point sufficiently close to an actual solution branch. There are two ways to obtain such data. The first is to take parameters corresponding to a Turing instability. The initial profile is then u(ξ ) =ū + ε cos(ξ ). The second way is to use a simulation where u approaches a travelling wave. Then the final profile u p of u can be used as initial point for the continuation. This is sufficient for our novel method. For the DDE, we need the auxilary variables z, ψ, φ along the periodic orbit. For z we integrated u p with the trapezoidal rule over [ξ −c, ξ]. Convoluting u p with the connectivity w yields ψ and φ is obtained from numerical differentation of ψ. The initial parameters are the same as in the simulation.

Parameter dependence of the travelling waves
From Fig. 2 we find Turing instabilities for r = 13 at θ 1 = 0.3038 and θ 3 = 0.3018 with ω = 0.6229 and ω = 4.088, respectively. Starting the numerical continuation from T 3 and varying θ we find travelling waves, see Fig. 6. First they are unstable, but the branch turns at θ = 0.2747 where the travelling waves are stable until θ = 0.3458. In between, near θ = 0.305, there is bistability of two different waves where the branch turns twice. The difference in the profile is an additional local minimum, present for lower values of θ . Finally, the branch ends at T 1 . For r = 10 there is only one Turing instability for θ = 0.3046 with ω = 3.7941. Following this branch we encounter similar scenarios, but the branch ends by approaching a non-uniform stationary profile near θ = 0.3060 when the speed c vanishes. Also the amplitude shows that the wave emanates from these Turing points (at T 1,3 ). The amplitude of travelling waves near Turing points is also illustrated in Fig. 7. This shows for fixed θ = 0.3008 that the prediction of the amplitude equations and results of the numerical continuation match well. We were unable to observe these small amplitude waves in simulations close to the Turing points corroborating that these bifurcations are subcritical. We have also investigated the effect of other system parameters, see Fig. 8. For decreasing S the wave speed increases as waves more easily excite neighbouring tissue. Increasing r leads to faster stable waves. Figure 9 shows the dependence of the wave speed on the spatial period Δ. The direct continuation of the integral equation and those obtained using (21) agreed very well, i.e. to the first four digits. The symbols P(n) indicate the number of pulses within a period Δ and P(2 + 1) differs from P(3) in that the interpulse distance is different, see Fig. 11. The P(1) solution is expected to be stable (using a kinematic analysis) when c (Δ) > 0 A kinematic theory of wave propagation is one attempt to follow the progress of localised pulse shapes, within a periodic wave, at the expense of a detailed description of their shape (Rinzel and Maginu 1984). Suppose that a pulse has a well defined arrival time at some position x then we denote the arrival of the nth pulse at position x by T n (x). A periodic wave, of period Δ, is then completely specified by the set of ordinary differential equations

Dispersion curves and kinematic theory
with solution T n (x) = nΔ + x/c, where c = c(Δ) is the dispersion curve, such as that obtained numerically in Fig. 9. The kinematic formalism asserts that there is a description of irregular spike trains in the above form such that where T n (x) − T n−1 (x) is recognised as the instantaneous period of the wave train at position x. A steadily propagating wave train is stable if under the perturbation T n (x) → T n (x) + u n (x) the system converges to the unperturbed solution during propagation, or u n (x) → 0 as x → ∞. For the case of uniformly propagating periodic traveling waves of period Δ we insert the perturbed solution in (23), so that to first order in the u n Thus, a uniformly spaced, infinite wave train with period Δ is stable (within the kinematic approximation) if and only if c (Δ) > 0. Hence, for the dispersion curves of shown in Fig. 9 it would seem to a first approximation that it is always the faster of the two periodic branches that is stable. Note that where there are bumps in the dispersion curve defining so-called supernormal wave speeds (wave speeds are faster than the corresponding speed of the large period wave) then it is only the supernormal wave of smaller period that is stable. Corresponding conclusions can also be made about subnormal waves (waves of slower speed compared to the wave of large period) on the slower branch. Also shown in Fig. 9 are waves with multiple pulses per period Δ and we indicate the number n of pulses on a branch by P(n).
According to the kinematic prediction there is a change of stability at the stationary points of the dispersion curve, i.e. the extrema in Fig. 9. As the branch persists another solution branch must bifurcate from a stationary point. Therefore we expect these points to act as organising centers of the waves. Indeed with (21) we have verified that the 2-pulse solution starts from a period-doubling bifurcation very close to the highest stationary point. In addition we plotted the dispersion curve with (Δ/N ), i.e. period per pulse, see Fig. 10. This may be demonstrated on an infinite domain, but here we show the following simulations. Consider a periodic domain Δ = 4.4, where the P(2) dispersion curve has positive slope. If we consider two pulses at equal distance then this is equivalent to the P(1) branch at Δ = 2.2. Indeed, this is the result found in Fig. 4 (right), where we determined the wavespeed c ≈ 0.6310. Now we take a slightly displaced double one-pulse solution, i.e. one pulse starts at x = 1 and the other at x = 3.19. Initially these pulses travel as two separate pulses but then adjust their speeds, and inter-pulse time, to travel together at a slightly lower speed c = 0.6302, see Fig. 10.
The travelling wave solution with three pulses corroborated the chaotic saddlefocus scenario even further as it formed an isolated branch in a two parameter diagram as expected, Fig. 11 (left). Interestingly we found that the P(3) and P(2 + 1) branch were different in the inter-peak distance, see Fig. 11 (right). On the P(3) branch the three pulses travel together where the distance between the first and second and that of the second and third is nearly equal around 1.639, while on the P(2 + 1) branch the third pulse follows at a distance of 2.45 from the second. Importantly we have included an absolute refractory process as in the original work of Wilson and Cowan (1972) and shown how to analyse this using a mixture of linear (Turing) analysis and novel numerical techniques. Despite the long history and extensive study of this type of model, to the best of our knowledge this is the first analysis of moving N -pulses in a neural field model with refractoriness. Moreover, we have shown that the types of travelling pulse patterns in this class of neural field model can be captured with a reduced kinematic description. This highlights the importance of the shape of the dispersion curve and its usefulness in predicting the behaviour of more exotic travelling wave packets. Given that other variants of neural field models, such as those that include axonal delays (Venkov et al. 2007), synaptic depression (Kilpatrick and Bressloff 2010), and slow inhibitory feedback (Taylor and Baier 2011) are also known to support periodic travelling waves it is of interest to construct dispersion curves for these models and contrast their shapes (and in effect the types of wave that they would be able to support). This is a topic of ongoing research and will be reported upon elsewhere.
where * denotes complex conjugation. The adjoint of L with respect to this inner product can be found as where We observe that the kernel of L † is spanned by the same set of functions as the kernel of L. Hence the solvability condition takes the form e ±i(ω c t∓k c x) , Lu n = L † e ±i(ω c t∓k c x) , v n = 0, n ≥ 2.