The importance of different timings of excitatory and inhibitory pathways in neural field models.

In this paper, we consider a neural field model comprised of two distinct populations of neurons, excitatory and inhibitory, for which both the velocities of action potential propagation and the time courses of synaptic processing are different. Using recently-developed techniques, we construct the Evans function characterising the stability of both stationary and travelling wave solutions, under the assumption that the firing rate function is the Heaviside step. We find that these differences in timing for the two populations can cause instabilities of these solutions, leading to, for example, stationary breathers. We also analyse “anti-pulses”, a novel type of pattern for which all but a small interval of the domain (in moving coordinates) is active. These results extend previous work on neural fields with space-dependent delays, and demonstrate the importance of considering the effects of the different time-courses of excitatory and inhibitory neural activity.


Introduction
Models of neural fields have been studied extensively over the last few decades (Hansel & Sompolinsky 1998;Pinto & Ermentrout 2001a, 2001bCoombes et al. 2003;Laing & Troy 2003;Coombes & Owen 2004;Coombes 2005), in the hope that they will provide information about possible macroscopic patterns of activity in the cortex. These patterns are on a much larger spatial scale than that of individual neurons, so the models take a continuum limit in which space is continuous and the mean firing rate of neurons is the variable of interest. Among some of the patterns of interest are stationary, spatiallylocalised "bumps" of activity (Pinto & Ermentrout 2001b;Laing & Troy 2003;Blomquist et al. 2005). These are thought to be involved in working memory and feature selectivity in the visual system (Hansel & Sompolinsky 1998). Travelling waves of activity are also of interest. Experimentally, these can be induced by stimulating pharmacologically treated neural tissue slices (Pinto & Ermentrout 2001a;Huang et al. 2004); they have also been observed in several different areas of the cortex of awake animals, often when no stimulation 152 C. Laing & S. Coombes is present (Ermentrout & Kleinfeld 2001). An important question is: given a neural model, what sorts of spatiotemporal patterns can exist and which features of the model are important in determining the existence and stability of these patterns? In this paper, we answer these questions concentrating on delays and differences in timing relating to neural processing.
It is well-known that different types of synaptic events have quite different time-courses. For example, inhibitory GABA A synapses decay with a time constant of approximately 15 ms, while excitatory NMDA synapses have a time constant of approximately 80 ms (Koch 1999). Also, the speed of propagation of an action potential along an axon depends on the diameter of the axon, as well as whether it is myelinated or not (Koch 1999). Measured axonal conduction velocities in the cortex can differ by a factor of 10, depending on the type of neuron examined, and can be as low as 1 m/s (Swadlow 1990). These velocities are also use-dependent and can change over time (Swadlow 1985). Interestingly, anaesthetics have been shown to increase conduction velocity in myelinated fibres -the consequences of this are discussed in (Swindale 2003). In this paper, we study a neural field model that explicitly includes differences in timing (in both the speed of synaptic processing and the conduction speed) for two populations of neurons, excitatory and inhibitory.
Many authors modelling neural systems have assumed instantaneous communication between different parts of the domain (Hansel & Sompolinsky 1998;Pinto & Ermentrout 2001b;Laing & Troy 2003;), but recently the effects of delays due to finite conduction speeds and synaptic processing have been investigated. Some of these more recent studies have involved rate models, in which the firing rate of neurons is the variable of interest (Pinto & Ermentrout 2001a;Coombes et al. 2003;Hutt et al. 2003;Coombes & Owen 2004;Hutt 2004) while several others have investigated networks of spiking neurons (Bressloff 2000;Golomb & Ermentrout 2000, 2002. In all of these papers only one type of neuron is considered, and hence there is only one delay (Golomb & Ermentrout (2002) consider both excitatory and inhibitory neurons, but only one conduction velocity). In this paper, we study the more realistic case of a network of excitatory and inhibitory neurons, each of which have their own associated conduction velocities (and therefore delays), and their own time-constants for synaptic processing. We would like to know whether this extra level of complexity brings with it any new phenomena that cannot occur in simpler models.
The structure of the paper is as follows. In the next section, we discuss the model integral equations and show that for certain choices of connectivity functions, they can also be formulated as PDEs. Then, we discuss stationary (i.e., time-independent) localised solutions, followed by travelling solutions. In both cases we use an Evans function approach (Coombes & Owen 2004) to determine the stability of solutions and the conditions necessary for both drift and breathing bifurcations. Thereafter, we examine the collision of two fronts, which leads us to study anti-pulses. We also briefly consider networks in which excitatory connections have longer spatial range than inhibitory ones. We finish with a discussion in the final section.

The model
We analyze a neural field model with synaptic activity u = u e (x, t) − u i (x, t), x ∈ R, t ∈ R + , where u a , a ∈ {e, i }, is given by the integral equation Different timings in neural field models 153 where "e" and "i" label excitatory and inhibitory populations, respectively. Here, the symbol * represents a temporal convolution in the sense that The function η a (t) (with η a (t) = 0 for t < 0) represents a synaptic filter, whilst w a (x) = w a (|x|) is a synaptic footprint describing the anatomy of network connections and v a is the speed of axonal transmission for population a. The function f represents the firing rate of a single neuron. For the rest of this paper we shall take often choosing e = 1 and i = . Thus is a measure of the strength of the inhibitory population relative to that of the excitatory. It should be noted that Golomb & Ermentrout (2000) found some qualitatively different results in systems with exponential coupling (as above) as opposed to systems with other forms of coupling (Gaussian, square), and the same may happen for the system studied here.
The interpretation of this model is that it is the value of the synaptic activity, u(x, t), that determines whether neurons at position x are firing at time t. If they are, action potentials from population a then propagate at speed v a to other positions in the network, where they are synaptically filtered by other neurons in population a with function η a . Thus, we have assumed that interactions of this form within a population are much stronger than those between populations. The two populations only interact through the difference, u e − u i , driving both populations. This model is novel in its consideration of different conduction velocities for the two types of neuron. In the discussion we discuss a more general model, in which we have different firing rate functions for both populations and the possibility of synaptic interactions between the two populations of neurons, not just within them.
This non-local description of neural tissue neglects local delays (say within a hypercolumn of spatial scale 1 mm), though incorporates delays due to the finite velocity of action potential propagation between distinct cortical regions. For conduction velocities in the range 1.5-7 m/s (typical of white matter axons) these are expected to be significant over scales ranging from a single cortical area (of spatial scale 10 mm) up to the scale of inter-hemispherical collosal connections. Anatomical surveys show that 80% of the synapses of long-range lateral connections connect directly between pyramidal cells, which are thought to make excitatory synapses only (McGuire et al. 1991). The other 20% of the connections target inhibitory interneurons which in turn contact the pyramidal cells, and thus represent inhibitory connections. Even though the inhibitory connections are outnumbered, the net effect of having two distinct delayed pathways has often been ignored in modelling studies.

A PDE description
For the particular choice of synaptic footprint Equation 4 it is possible to obtain an equivalent PDE description of the integral Equation 2, using ideas developed by Jirsa & Haken (1996). To see this we write where and we use the notation ρ(x, t) = f • u(x, t). Introducing Fourier transforms of the following form allows us to write It is straightforward to show that the Fourier transform of Equation 6 is where We have using Equations 9 and 10 that where ω a = v a /σ a . We may now write Equation 8 as which upon inverse Fourier transforming gives the PDE: If we choose the synaptic response η a (t) to be the exponential: η a (t) = α a (t) exp (−α a t), where (t) is the Heaviside function, defined by (t) = 1 for t ≥ 0 and zero otherwise, we can also write the integral Equation 1 as the differential equation In numerical simulations of the model we integrate Equations 13 and 14 using finite difference approximations to the spatial derivatives and work with the choice where h can be thought of as a threshold and β as a gain parameter. Note that if we set v e = v i = v and η e = η i = η, i.e., we remove the differences in timings for the two neural populations, our system reduces to where w(y) ≡ w e (y) − w i (y), which is typically of "Mexican hat" shape. This equation has been studied elsewhere (Pinto & Ermentrout 2001a;Coombes et al. 2003;Coombes & Owen 2004), and provides a useful comparison. We now consider stationary spatially-localised patterns, or "bumps".

Stationary bump solutions
For general time-independent solutions, Equation 14 gives u a (x, t) = ψ a (x, t), and u(x, t) can be replaced by q (x), where q (x) satisfies where w(x) = w e (x)−w i (x). Throughout the rest of this paper, we shall focus on the particular choice of a Heaviside firing rate function, f (u) = (u − h). In this case, once the position of the threshold crossings are known, the explicit dependence on q of the right hand side of Equation 18 is removed. Apart from allowing an explicit construction of travelling bumps and waves this choice also allows for a direct calculation of wave stability via the construction of an Evans function (Coombes & Owen 2004). Although often chosen for mathematical reasons, the Heaviside function may be regarded as a natural restriction of sigmoidal functions (such as Equation 15) to the regime of high gain (β → ∞). Importantly, numerical simulations show that many of the qualitative properties of solutions in the high gain limit are retained with decreasing gain (Ermentrout & McLeod 1993;Pinto & Ermentrout 2001a;Coombes et al. 2003).

Existence
One-bump solutions with q (x) ≥ h for 0 < x < are given by Note that the system 18 is translationally invariant, so we can choose one end point of the bump to be at x = 0. The width of the bump is then determined by the self-consistent solution of q (0) = h = q ( ), which gives the equation: In Figure 1, we show a typical plot of bump width as a function of the threshold parameter h, for particular values of the other parameters. We note that the differences in timings for the two populations do not affect the existence of stationary patterns, only their stability, as seen below.

Stability via the Evans function
We can find the stability of stationary one-bump solutions by constructing the associated Evans function. Evans functions were first developed to study the stability of travelling waves in PDEs (Evans 1975). In essence the Evans function is an analytic tool whose zeros correspond to eigenvalues of the linearised problem obtained after considering perturbations around a travelling wave solution. Moreover, the order of the zero and the multiplicity of the eigenvalue match. A recent review of their use in determining the stability of travelling pulses in dissipative systems can be found in (Kapitula 2005). This theory has recently been extended to cover travelling waves in systems with nonlocal interactions, such as studied here (Coombes & Owen 2004;Kapitula et al. 2004;Pinto et al. 2005).
Generalising results in Coombes & Owen (2004) we may construct the Evans function of the stationary one-bump solution as where For our choice of an exponential synaptic response η a (t) = α a (t) exp (−α a t), we have simply thatη Note that we can obtain q (0) by using The discrete spectrum for the one-bump solution is then given by the zeros of the Evans function E(λ) = 0, so that a solution is stable if the spectrum only resides in the left hand complex plane (i.e. Re λ < 0). Note that a zero eigenvalue (satisfying E(0) = 0) is always expected due to translational invariance of the one-bump solution (with the corresponding eigenfunction q (x)). One natural way to find the zeros of E(λ) is to write λ = ν + i ω and plot the zero contours of Re E(λ) and Im E(λ) in the (ν, ω) plane. The Evans function is zero where the contours intersect. Note that it is sufficient for us to determine the location of the isolated spectrum for wave stability, since the systems under study in this paper are such that the real part of the continuous spectrum has a uniformly negative upper bound (Coombes & Owen 2004).

Examples
We now α e = α i = 1 and v i = 1 and decrease v e , we find that the wide stationary bump becomes unstable due to a single eigenvalue passing from the left half plane to the right half plane along the real axis. This is shown in Figure 2, where we plot contours of the real and imaginary parts of the Evans function over part of the complex (λ) plane for v e = 0.25 and v e = 0.15. The Evans function has zeros where the two contours cross. Note that we always have a zero at the origin, as expected. This instability manifests itself as a transition to a moving, or drifting, bump, as shown in Figure 3. Here we run the system with v e = 0.25 up until t = 100, when we switch to v e = 0.15, and add a small random perturbation to the solution. We see that the previously stable stationary bump is replaced by a stable moving bump. If the position of a bump of activity is used to encode some aspect of an action to be performed in the future (Goldman-Rakic 1995), this instability is undesirable. For simulations of the full system we use Equation 15 as the firing rate function with β = 150. This moving bump will be discussed later.
Conversely, by holding v e constant and decreasing v i we can make the wide stationary bump in Figure 1 go unstable through a Hopf bifurcation, in which a pair of complex eigenvalues pass through the imaginary axis. This is shown in Figure 4, where the Evans function is represented for v i = 0.4 (left) and v i = 0.2 (right). We see that between these values, a complex Figure 3. Instability of a stationary bump due to an eigenvalue passing through zero. At t = 100, v e was changed from 0.25 to 0.15 (and a small random perturbation was added). u e is plotted, red is high, blue is low. Other parameters are σ e = 1, The domain was discretised with 200 spatial points and periodic boundary conditions were used. pair of eigenvalues moves into the right half plane. The instability is shown in Figure 5, where we switch v i from 1 to 0.2 halfway through the simulation. We see an oscillatory instability develop, but it appears to be subcritical, and the system moves to the all-off state. Clearly this is also undesirable in the context of working memory. These instabilities are summarised in Figure 6, where we plot the curves of drift instabilities and Hopf bifurcations in the v e − v i plane. Note that we have followed only the curve of Hopf bifurcations corresponding to the imaginary roots with smallest imaginary part. Although we have not proven it, numerical computation of the Evans function at various points in the v e − v i plane suggests that the stationary bump is always stable in the wedge between the two curves.
We can also find supercritical Hopf bifurcations if we relax the restriction that α i = α e = 1, i.e., if we allow different time-courses for the synaptic processing in the two populations. In Figure 7 we show the effect of decreasing v e from 0.8 to 0.5, where α i = 1.8 and α e = 3. Oscillations appear, and their amplitude rapidly saturates. Presumably, by changing other Figure 5. Instability of a stationary bump due to a subcritical Hopf bifurcation. At t = 70 we switched v i from 1 to 0.2 (and a small random perturbation was added). Other parameters are v e = 1, α i = α e = 1, h = 0.1, σ e = 1, σ i = 2 and = 1. u e is plotted, with red high and blue low. 200 spatial points were used. Figure 6. The curve of drift instabilities, defined by E (0) = 0 (solid), and the curve of Hopf bifurcations, for which E(λ) has a conjugate pair of purely imaginary roots (dashed). Other parameters are α e = α i = 1, h = 0.1, σ e = 1, σ i = 2 and = 1. The stationary bump appears to be stable in the wedge between the two curves and unstable otherwise.

Discussion
It is possible to show that bifurcations of a stationary bump with an eigenvalue passing through zero are not possible for Equations 1-2 if v e = v i , α e = α i , and σ e < σ i (as here) i.e., this is a novel instability due to the differences in timings for the two populations. To show this, one calculates E(λ) and obtains where v e = v i = v, α e = α i = α, and we have used the fact that w(0) > w( ). Clearly E(0) = 0. If w( ) < 0, we also have that E (λ) > 0 for all λ > 0, and hence there can be no other positive real roots of E(λ). The condition w( ) < 0 means that we are on the upper branch in Figure 1. Similarly, it is possible to show that there are no purely imaginary roots of E(λ). We substitute λ = i ω into Equation 25, where ω ∈ R, ω = 0 and set the real and imaginary parts equal to zero, obtaining where θ = 2ω /v. Assuming that w( ) = 0, solving (27)  we see that it cannot simultaneously be satisfied, i.e., E(λ) has no purely imaginary roots. If w( ) = 0 the only root of E(λ) is λ = 0, which corresponds to the saddle-node bifurcation in Figure 1. Thus both types of instabilities are a result of the differences in timings for the two populations.
Bifurcations of the types discussed in this section have recently been observed in a system with one neural population but in which the threshold is a dynamic variable (Coombes & Owen 2005). There, the authors also saw a stationary bump start to move as an eigenvalue moved through zero, and stationary breathers caused by a supercritical Hopf bifurcation. Breathers have also been observed in neural field equations by Bressloff et al. (2003) and , but in those papers the authors made the domain inhomogeneous, inducing a bump to occur over a spatially-localised input. During normal awake operation the cortex continuously receives inhomogeneous inputs, so the response of a neural model with such inputs is of interest. However, under general anaesthesia the cortex will receive less input (and conduction velocities may be increased (Swindale 2003)) so it is also of interest to study the existence and stability of patterns in this situation. Blomquist et al. (2005) recently studied a two-layer neural field model that can be best thought of as a generalisation of that studied by Pinto & Ermentrout (2001b). In contrast with Pinto and Ermentrout, Blomquist et al. included inhibitory-to-inhibitory connections and did not assume that the firing rate function for the inhibitory population was linear. They did not include conduction delays, but obtained both stable and unstable breathers that were created in Hopf bifurcations, as seen here.
The common theme between our results and those of Coombes & Owen (2005), Bressloff et al. (2003),  and Blomquist et al. (2005) is the presence of a second variable describing either another population of neurons or a local field such as an adaptation current.

Travelling wave solutions
We now discuss travelling wave solutions, both fronts (which connect a region of high activity to one of low activity) and pulses (travelling bumps, before and after which the medium is quiescent). We introduce the coordinate ξ = x − ct and seek functions U(ξ, t) = u(x − ct, t) that satisfy the full integral model equations. In the (ξ, t) coordinates, these integral equations can be expressed as Different timings in neural field models

Fronts
We look for travelling front solutions such that q (ξ ) > h for ξ < 0 and q (ξ ) < h for ξ > 0. It is then a simple matter to show that (for the special case of the Heaviside firing rate function chosen earlier) The choice of origin, q (0) = h, gives an implicit equation for the speed of the wave as a function of system parameters: Both equations may be written as quadratic equations in c. To ensure that lim ξ →−∞ q (ξ ) > h, we must have R w(y)dy = 1 − > h. We note that a standing front (c = 0) occurs when 2h = 1 − . Also, it is not necessary to have different propagation velocities and synaptic processing time-constants for such fronts to exist; they are seen robustly in simpler systems where these parameters are equal (Pinto & Ermentrout 2001a). Once again the Evans function is easily obtained using the techniques described in Coombes & Owen (2004) and may be written in the form where H(λ) = H e (λ) − H i (λ) and For the exponential synaptic response we have simply that For a standing front it is a simple matter to check that E(λ) = 0 when Hence, there is a bifurcation of the standing front when α e σ e = α i σ i and 2h = 1 − . To determine the type of bifurcation, one examines various partial derivatives of the functions in Equations 33-34, evaluated at the bifurcation point (see, e.g., Sec. 3.1 of Wiggins 1990). It can be determined that the bifurcation is a simultaneous pair of transcritical bifurcations, each creating a branch with c = 0 (This bifurcation was misidentified as a pitchfork bifurcation in  and Coombes & Owen (2004)). The transcritical bifurcations are shown in Figure 8 (left) for specific values of the other parameters. When the condition 2h = 1 − does not hold, the transcritical bifurcations generically break into a single saddlenode bifurcation, as seen in Figure 8 (right). Similar results have been seen before in a one-layer model with a linear recovery variable Coombes & Owen 2004). Although we have not seen Hopf bifurcations of fronts, it is known that making the domain inhomogeneous can induce these (Bressloff et al. 2003). Figure 9 shows another plot of front speed, now as a function of h. We see that for these parameter values there is a minimum speed below which stable fronts do not exist. We will return to this figure later. Finally, we see in Figure 10 a result of choosing different parameters. For these parameters, only fronts with slow enough speeds are stable. This is in  strong contrast with results from networks with purely excitatory coupling (Bressloff 2000;Golomb & Ermentrout 2000) in which travelling structures cannot travel with arbitrarily slow speeds.

Travelling pulses
We now study travelling pulses, which are bumps of the form studied earlier, but which have speed c = 0. The travelling pulses have the form q (ξ ) ≥ h for ξ ∈ [0, ] and q (ξ ) < h otherwise. In this case the expression for ψ a (ξ ) is given by Coombes et al. (2003). where The dispersion relation c = c( ) is then implicitly defined by the simultaneous solution of q (0) = h and q ( ) = h ( > 0). For an exponential synapse these two conditions are Note that as c → 0, both 41 and 42 reduce to the equation governing the width of a stationary bump, Equation 20, as expected. The Evans function takes the same form as described previously, with Equation 21 replaced by Note that as c → 0, the matrix 43 reduces to the matrix associated with the stability of a stationary bump, 21, as expected. To evaluate the derivatives of q in Equations 46-49 we use q = q e − q i , where q a = α a (q a − ψ a )/c. Specifically, we have and We now discuss several examples. In Figure 11 we plot the speeds and widths of a pair of travelling pulses (one stable and one unstable) as v e is varied. The other parameter values are the same as those in Figure 3, and the instability seen in that figure can can now be understood with reference to Figure 11. For v e = 0.25 there is no stable travelling pulse, but for v e = 0.15 there is one, with speed c ≈ 0.05. (Of course, a similar stable pulse moving in the opposite direction also exists.) Interestingly, this pair of stable pulses seem to be created in a pair of transcritical bifurcations, in the same way that fronts were created in the previous Figure 11. Left: pulse speeds as a function of v e . Right: pulse widths as a function of v e . Solid lines indicate stable solutions, dashed unstable. Note that the speeds are always less than v e . Other parameters are section. This is in contrast with the pitchfork bifurcation in speed seen in other neural field models (Laing & Longtin 2001). Note that for these parameters a moving front does not exist, as the condition 1 − > h does not hold.
In Figure 12 we show the width and speed of a moving pulse as a function of h. For a range of values of h, there are two pulses, a fast, wide one and a slow, narrow one. By plotting the Evans function on the upper branch (not shown) we see that there is a Hopf bifurcation as h increases through h ≈ 0.07. We use this information to indicate the stability of the branches in Figure 12. This Hopf bifurcation appears to be subcritical. In Figure 13 we show a simulation that starts with h = 0.05, for which the moving pulse on the upper branch in Figure 12 is stable. At t = 100, h is switched to 0.07. The pulse starts to oscillate, but the oscillations grow in magnitude until the pulse is destroyed. Other families of travelling pulses are shown in Figure 8 (red curves).
Supercritical Hopf bifurcations of moving pulses, leading to travelling breathers or "lurching" waves, have been observed in several other neural systems (Golomb & Ermentrout 2000;Coombes & Owen 2005). However, we could not find parameters for the system under study for which this type of bifurcation occurred.

Other solutions
In this section, we discuss colliding fronts, "anti-pulses", and the patterns seen when inverted Mexican-hat connectivity is used. Returning to Figure 9, we see that for these parameters and h in some interval whose lower endpoint is h = 0.075, if we were to start a wave consisting of a large active interval the "front" of the wave would move more slowly than the "back".  The back would eventually catch up with the front. A simulation showing this can be seen in Figure 14, and we see that the result is a moving pulse from the stable family shown in Figure 8b.
If, however, we choose h slightly less than 0.075, the opposite will occur, i.e., the back moves more slowly than the front and the active region expands in width. On a periodic domain this occurs until the front meets the back, from behind, as seen in Figure 15 (top). We call the final solution a moving "anti-pulse", since all but a small part of the domain is active.
The formation of a moving pulse by the catching up of a back to a front was seen in Coombes & Owen (2004), but these authors did not mention the formation of anti-pulses, although they are most certainly expected even when v e = v i and α e = α i . We now analyse anti-pulses.

Anti-pulses
For anti-pulses, q (ξ ) < h for 0 < ξ < and q (ξ ) > h otherwise. Using Figure 14. A wide moving pulse for which the back travels faster than the front. The eventual solution is a moving pulse, of the type analysed previously. u e is plotted, with black being high and white low. Periodic boundary conditions are used. h = 0.115 and other parameters are as in Figure 9. where F a is given by Equation 40. The conditions q (0) = h = q ( ) are then easily determined using q (ξ ) = q e (ξ ) − q i (ξ ) and q a (ξ ) = ∞ 0 η a (s )ψ a (ξ + cs )ds (54) and the fact that these integrals have essentially been done in the determination of 41-42. They are and The width and speed of anti-pulses is determined by the simultaneous solution of 55-56. The Evans function for antipulses has the same form as that for pulses, since only the fact that q (0) = q ( ) = h is used, the sign of q (ξ ) − h for ξ ∈ (0, ) being irrelevant.
In fact, the analysis of anti-pulses does not bring any new results, as the system under study is symmetric about the "balanced" parameter point 2h = 1 − , at which there are stationary fronts. If we make the replacement h → 1 − − h in 55-56, we obtain 41-42, and vice versa. Thus, for a given value of h, say h * , at which there exists a pulse, there exists a corresponding antipulse when h = 1 − − h * (all other parameters being the same). It will have the same width (as determined by the zeros of q (ξ ) − h), speed and stability as the pulse.

Inverted Mexican hat connectivity
Although much work on pattern formation in neural field models has used one layer of neurons with Mexican-hat connectivity, for which inhibitory connections have a wider spatial extent than excitatory, there is evidence that the opposite is true, at least in some contexts (Swindale 1996). We now briefly analyse Equations 1-2 with coupling function 4, but with σ e > σ i , i.e., with the excitatory connections having a greater spatial extent than the inhibitory, which we refer to as inverted Mexican-hat connectivity.
As before, we can find families of moving pulses. In Figure 16 we show solutions of 41-42 as h is varied for inverted Mexican-hat connectivity. Note that not every point on the curves in Figure 16 corresponds to a one-pulse solution, as some of the solutions may have q (ξ ) > h for more than one interval. In Figure 17 we show a stable travelling one-pulse solution from the curve in Figure 16 pulses interact (bottom panel, same parameters). For this connectivity, we can also obtain families of travelling fronts (not shown).

Discussion
We have studied stationary and travelling bump and front solutions of a two-layer neural field model with different conduction velocities and synaptic processing time-constants for the two populations. By varying these parameters we have found bifurcations of stationary bumps to both travelling and breathing bumps. These bifurcations can be found by explicitly constructing an Evans function for these solutions and they cannot occur if the synaptic time-constants and conduction velocities are the same for both layers.
Our work has produced results similar to those of several other groups. For example, Curtu & Ermentrout (2004) recently studied an extension of the system first discussed by Hansel & Sompolinsky (1998). This model had one neural population, Mexican-hat type connectivity, an adaptation variable and no delays. The authors found travelling and standing waves, and stationary, spatially-periodic patterns. However, their results were derived by linearising about the spatially-uniform state, and are thus unable to say anything regarding spatially-localised patterns of the type studied here. Golomb & Ermentrout (2000) studied the effects of delays on propagating activity. They included a fixed delay and found that increasing this led to lurching waves. (We did not include a fixed delay, but see below.) However, they used a spiking neural network in which each neuron could only fire once, and only excitatory coupling. Because of this, they could not find stationary or arbitrarily slowly moving patterns, as we have. In later work Golomb & Ermentrout 2001, 2002, these authors studied a network with both excitatory and inhibitory populations but did not include conduction delays for most of their analysis, and still allowed neurons to fire at most once, thus precluding the existence of stationary patterns. One interesting result that they found was the coexistence of both fast and slow propagating pulses, which we have not found. However, these authors found that once neurons were allowed to fire multiple times this bistability disappeared, with only the slow pulses persisting. Blomquist et al. (2005) studied a two-layer neuronal network without delays (an extension of that studied by Pinto & Ermentrout (2001)), and found both subcritical and supercritical Hopf bifurcations of stationary bumps. Coombes & Owen (2005) studied a single neural layer with Mexican-hat connectivity and a variable representing spike frequency adaptation. They found both drifting and breathing bifurcations of stationary bumps, as we have, and also supercritical Hopf bifurcations of travelling bumps, which we have not found. We now discuss possible extensions of the work presented here.
One extension would be to break the homogeneity of the domain (reflected by the appearance of x and y in Equation 2 in only the combination x − y) as Jirsa & Kelso (2000) have done, but for the system studied here. These authors modelled heterogeneity by putting a direct connection from one part of the domain to another and studied the effect of varying the length of this connection. Roxin et al. (2005) recently studied a generalisation of a neural field model first presented by Hansel & Sompolinsky (1998) in which there is a fixed delay in the nonlinear term, as opposed to the space-dependent delays that we have. This delay is meant to mimic those due to synaptic and dendritic processing. These authors found a wide variety of spatiotemporal patterns. Including such a term in our equations would involve replacing u(x − y, t − |y|/v a ) in Equation 2 by u(x − y, t − |y|/v a − D), where D is a fixed positive delay. As mentioned above, this type of delay was studied by Golomb & Ermentrout (2000). As an example of the effect of including such a term, it is straightforward to show that the equations governing the speed of a front 33-34 would be modified to and it seems likely that all of the calculations performed here could also be done with such a term included. Hutt (2004) discussed a similar idea, but wrote the nonlinear term in 2 as a linear combination of a term whose delay depends on distance and one whose delay is fixed. Further extensions could include verifying some of the predictions here with a network of spiking neurons. This would be computationally intensive due to the inclusion of delays, but similar work has been performed (Golomb & Ermentrout 2000). An important extension that would make more valid the results here and elsewhere (Coombes & Owen 2004, would be to derive similar results for a continuous firing rate function, f (Curtu & Ermentrout 2004). Also important is the consideration of a two-dimensional domain since the cortex is best thought of as a two-dimensional sheet of interconnected neurons. There have been recent results on patterns (multi-bump solutions, breathers and spiral waves) in two-dimensional neural fields (Laing & Troy 2003;Laing 2005), but none of these models have included even one conduction velocity or delay. While it seems that including a finite conduction velocity does not destabilise travelling bumps or fronts in one spatial dimension (Golomb & Ermentrout 2000;Coombes 2005), it is not clear whether the same holds in two dimensions. A more general model, more clearly differentiating the two neural populations, would be u a = η a * ψ a , a ∈ {e, i } (59) Here we not only have different conduction velocities v e and v i and different synaptic filters η e and η i , but different firing rate functions f e and f i and four coupling functions, w ee , w ei , w ie and w ii instead of the two in Equations 1-2. Choosing f a (u) = (u − h a ), i.e., using the Heaviside function as the firing rate function, with two different thresholds, we should be able to analyse Equations 59-61 in much the same way as we have analysed 1-2 in this paper. Some of the analysis in Coombes & Owen (2005), in which the threshold is a dynamic variable, should be applicable to the analysis of 59-61 when h e = h i . The model 1-2 can be considered as intermediate between 59-61 and previous models in which there is only one population of neurons, one synaptic filter, and one conduction velocity (Coombes et al. 2003;Coombes & Owen 2004). Guo et al. (2005) have recently studied a pair of coupled delay-differential equations that are similar in structure to 59-61, as have Shayer & Campbell (2000), although those systems have no spatial structure. Note that after setting v e = v i = ∞ in 60-61 and choosing η e (t) = (t)e −t and η i (t) = (t)e −t/τ /τ we recover the model originally presented by Pinto & Ermentrout (2001) and later analysed by Blomquist et al. (2005).