A next generation neural field model: The evolution of synchrony within patterns and waves

Neural field models are commonly used to describe wave propagation and bump attractors at a tissue level in the brain. Although motivated by biology, these models are phenomenological in nature. They are built on the assumption that the neural tissue operates in a near synchronous regime, and hence, cannot account for changes in the underlying synchrony of patterns. It is customary to use spiking neural network models when examining within population synchronisation. Unfortunately, these high dimensional models are notoriously hard to obtain insight from. In this paper, we consider a network of $\theta$-neurons, which has recently been shown to admit an exact mean-field description in the absence of a spatial component. We show that the inclusion of space and a realistic synapse model leads to a reduced model that has many of the features of a standard neural field model coupled to a further dynamical equation that describes the evolution of network synchrony. Both Turing instability analysis and numerical continuation software are used to explore the existence and stability of spatio-temporal patterns in the system. In particular, we show that this new model can support states above and beyond those seen in a standard neural field model. These states are typified by structures within bumps and waves showing the dynamic evolution of population synchrony.


I. INTRODUCTION
The act of passing information between brain regions produces waves of neural activity. These waves are readily observed using noninvasive techniques such as electroencephalography (EEG) and magnetoencephalography (MEG) [1], as well as in brain slices [2]. Both experimental and theoretical work has shown that EEG and MEG recordings and evoked potentials can exhibit traveling and standing waves [3]. In particular, traveling waves are seen in EEG sleep recordings propagating across the cortex at a speed of about 1.2-7.0 m/s [4]. Standing waves are often associated with idle brain states. For example, standing waves at α frequency (8)(9)(10)(11)(12)(13) are observed in the vicinity of the visual cortex when the subject has their eyes closed [5]. Another commonly observed spatial pattern is the so-called bump attractor. This spatially localized increase in population firing is produced in working memory tasks and the location of the bump can be linked to memory location [6]. * aine.byrne@nyu.edu Traditionally, neural field models are used to describe wave and bump states in the brain. Although inspired by biology, these models are entirely phenomenological in nature. Even so, they have been particularly successful in describing neurophysiological phenomena, such as EEG and MEG rhythms [7], working memory [8], binocular rivalry [9], and orientation tuning in the visual cortex [10]. They are typically cast as a system of nonlocal differential equations which describe the spatiotemporal evolution of coarse-grained population variables, such as the firing rate of a neuronal population, the average synaptic current, or the mean membrane potential [11].
The first attempt at a neural field mode is attributed to Beurle [12], who built a model to describe the propagation of activation in a given volume of neural tissue. This model was purely excitatory, but even so allowed him to examine the propagation of large-scale brain activity. In the 1970s Wilson and Cowan [13,14] extended this model to include a second inhibitory layer. Unlike Beurle, they were interested in spatially localized bump solutions. In his seminal paper, Amari [15,16] created what is now known as the standard neural field equation. By introducing a Mexican hat-type coupling function (local excitation and long-range inhibition), he reduced the model to a single equation with a mixture of excitatory and inhibitory connections. This allowed him to construct explicit solutions for spatially localized patterns and assess their stability (at least for a Heaviside firing rate function). For a review of the Amari model and bumps in one spatial dimension, see Ref. [17], and for a discussion of bubbles in two spatial dimensions, see Ref. [18].
One of the main assumptions in many types of neural field models, and especially those for describing EEG, is that pointwise they describe a density of neurons operating in a near synchronous regime [19]. This wholly reasonable assumption can be traced back to the observation that an EEG scalp electrode, which typically experiences the activity of roughly 10 9 cortical pyramidal cells, can only detect an electric field if all the individual cell dipoles add coherently [20]. However, a near synchrony assumption in any neural mass model precludes its use in describing the increase and decrease of power commonly seen in given EEG and MEG frequency bands. These temporal variations are believed to be the result of changes in synchrony within the neural tissue. The former phenomenon is called event-related synchronization (ERS), and the latter event-related desynchronization (ERD) [21]. Consequently, there is a pressing need to develop the next generation of neural field models, which include this notion of a dynamic within-population synchrony (not fixed to be near synchronous), to more accurately describe the evolution of large-scale spatiotemporal brain rhythms.
When looking at within-population synchronization one typically uses a spiking neural network model. However, these high-dimensional models are almost impossible to gain insight from. In an ideal world there would be a mathematical procedure for linking microscopic dynamics to macroscopic dynamics. This link has proved elusive for the majority of spiking models. However, Luke et al. [22] showed that the θ -neuron model is amenable to such a reduction for pulsatile coupling. Montbrió et al. [23] used a similar approach to reduce a network of quadratic integrate-and-fire neurons. Laing [24] has also shown that the same approach can be applied to a network of spatially extended θ -neurons, in the presence of gap junction coupling. In previous papers [25,26], we have shown that the approach of Luke et al. can be extended to incorporate a biologically realistic form of synaptic coupling. Here we build on this work to construct a neural field model that incorporates within population synchrony and a realistic form of synaptic coupling. We shall refer to this model as a next generation neural field model. Esnaola-Acebes et al. recently showed that a similar model of pulsecoupled excitatory and inhibitory quadratic integrate-and-fire neurons displayed damped temporal oscillations, which modulate a spatially nonlocalized profile [27]. Such oscillations are, however, not self-sustained. Here, we show, through the inclusion of synaptic dynamics, that the model can support self-sustained oscillations.
Unlike more traditional neural field models the one studied here has a population firing rate at a point in the tissue that depends on the degree of synchrony at that point. This firing rate is a real valued function of the complex Kuramoto order parameter, and as such different values of the order parameter can give rise to the same population rate. Thus, in contrast to phenomenological neural field models that are closed under the choice of a firing rate that is a nonlinear (typically sigmoidal) function of mean-membrane potential or synaptic activity, the one presented here is expected to have a richer set of responses, due to the more sophisticated representation of the response of the population rate to synchrony levels. This allows the model to represent underlying patterns of spike-train synchrony that cannot be accounted for in standard neural field models. In particular, we show, using a mixture of analysis and simulation, that this new neural field model can support exotic patterned states more reminiscent of highdimensional spiking networks, with spatiotemporal patterns showing the evolution of synchrony.
We begin with an overview of the model formulation in Sec. II and outline the necessary steps for reduction to a neural field. A Turing instability analysis of the model is covered in Sec. III. Here, we show that the system can be unstable to both static and dynamic Turing patterns for a wide window of parameter space. More interestingly, we show that when the Turing bifurcation collides with a Hopf bifurcation patterned states emerge in which there exists an oscillating structure within a spatially localized bump. The Turing analysis is complimented with a numerical bifurcation analysis for the full nonlinear model in Sec. IV. Here we further examine the emergent patterns away from bifurcation, as well as consider localized traveling waves. Finally, in Sec. V we discuss our main results as well as natural extensions of the work presented.

II. THE MODEL
We first consider a network of N coupled quadratic integrate-and-fire (QIF) neurons, uniformly distributed along a line of length such that the j th neuron is at position x j = − /2 + (j − 1) x, where j = 1, . . . , N and x = /(N − 1) is the spacing between neurons. The coupling between neuron i and neuron j depends only upon the distance between the two neurons, w m ij = w m (|x i − x j |), where m is a label used to keep track of neural subpopulations. We write the network dynamics for the voltage, subject to reset, v i → v r , whenever a firing threshold v th is reached by v i . The time at which the ith cell reaches threshold from below for the sth time will be denoted by T s i , s ∈ N. Here the background drives η i will be assumed to be heterogeneous and chosen from a Lorentzian distribution, where η 0 is the center of the distribution and is the half width. We assume a synaptic input current of the form for a global conductance g i m , synaptic reversal potential v m syn and local voltage v i . In what follows we will assume m = {1, 2} or m = 1, but one should note that this framework can be extended to include multiple types of synapses. The interplay of excitation and inhibition, on different spatial 012313-2 scales, is known to play an important role in the generation of global spatially patterned states [28][29][30]. Hence, we separate our synaptic conductance into two parts, such that we have one excitatory synaptic current (v syn > 0) and one inhibitory (v syn < 0), each with different spatial ranges. Each of the synaptic conductances will be taken to mimic that of a synapse with a finite rise and fall time that evolves according to for some coupling strength κ m , where Q is the linear secondorder differential operator Here τ m is the synaptic time scale, and Q m has a response (Green's function) s(t ) = τ −2 m te −t/τ m for t 0 (and is zero otherwise), which is a popular choice for many synapse models in computational neuroscience [31].
It is well known that the QIF model is formally equivalent to the θ -neuron model [32] under the transformation v i = tan(θ i /2), for θ i ∈ [0, 2π ) (when the threshold v th and reset v r are set to +∞ and −∞, respectively). This relationship allows us to construct a θ -neuron network dynamics aṡ To obtain Eq. (7) we have made use of the fact that δ(t − T s j ) = δ(θ j − π )|θ j (T m j )|. As such, we say that neuron j "spikes" whenever θ j increases through π . The network formulation in terms of dynamics on a circle is particularly useful since we no longer have to worry about handling the discontinuous reset process as we would have to do for a QIF network.

Mean field limit
We take the large N limit, N → ∞, which allows us to describe the system in terms of a continuous probability distribution function ρ(x, η, θ, t ), with x, η ∈ R, θ ∈ [0, 2π ), and t ∈ R + , which satisfies the continuity equation: where v θ is the following realization of Eq. (6), The mean-field representation of the synaptic inputs g m are written as follows: where we have used the result 2πδ(θ − π ) = l∈Z e il(θ−π ) to write the right hand side in terms of exponentials. The formula for v θ Eq. (9) may be conveniently written in terms of e ±iθ as where φ = (η + I − 1)/2 and χ = η + I + 1. Note that a similar approach has previously been considered by Laing [24]. However, his focus was on smooth (nonpulsatile) interactions with a first-order model of the synapse, and he did not consider reversal potentials. To reduce the system we make use of the Ott-Antonsen (OA) ansatz [33]. This decomposes ρ in the form ρ(x, η, θ, t ) = L(η)F (x, η, θ, t )/(2π ), where F is 2π -periodic in θ with a Fourier series representation F (x, η, θ, t ) = n F n (x, η, t )e inθ . The OA ansatz restricts the choice of F n such that F n (x, η, t ) = α(x, η, t ) n , with |α(x, η, t )| < 1 to ensure convergence. Hence, ρ can be written as where cc denotes the complex conjugate of the previous term. Substituting Eq. (11) into the continuity Eq. (8) and balancing terms in e iθ gives an evolution equation for α: We define the Kuramoto order parameter as follows: where |z| 1. The Kuramoto order parameter is a complex number z = Re i , whose magnitude R represents the degree of within population synchrony and angle represents the average phase of the population. Substituting Eq. (12) into Eq. (14) we find that wherez denotes the complex conjugate of z. As the Lorentzian has two simple poles η ± = η 0 ± i , the above integral may be performed by choosing a large semicircular contour in the lower half η-plane and using the residue theorem, to yield We use Eq. (14) to write Eq. (10) as where ⊗ represents a spatial convolution and f is the population firing rate Note that Eq. (16) takes the form of a generalised neural field equation, where the firing rate function f is a derived quantity that depends on the within population synchrony. Also noteworthy is the fact that this firing rate function is not a sigmoid. It is a highly nonlinear function that depends on the intrinsic population dynamics. The firing rate can be plotted as a function of the Kuramoto order parameter (Fig. 1). As expected, the firing rate is highest at z = e iπ (where a single neuron fires as θ increases through π ). The dynamics of z are obtained by evaluating Eq. (13) at η − = η 0 − i , and taking the complex conjugate, which gives As in Ref. [26], we interpret Eq. (17) as describing the intrinsic population dynamics and Eq. (18) as the dynamics generated by synaptic coupling.
In summary, we have the following evolution equations: where G m (z, g m ) ≡ G(z, g m ; v m syn ) and we have omitted the dependence on control parameters. The dynamical system Eqs. (19) and (20) have therefore one complex (z) and |m| real (g m ) state variables. In what follows, we choose normalized exponentially decaying synaptic kernels of the form w m (x) = exp(−β m |x|)/(2/β m ), with Fourier transform w m (k) = 1/[1 + (k/β m ) 2 ]. The spatial scale is set by the choice of 1/β m in the exponentially decaying kernels. We shall choose to be several times larger than max m β −1 m . As the Fourier transform is a rational function, the evolution Eqs. (19) and (20) are equivalent to the PDEs [34]: In what follows, we will show that the model can support Turing patterns, traveling waves, and other complex global spatiotemporal patterns, for m = {1, 2} and traveling fronts for m = 1. In the former case, this is achieved by choosing β 1 > β 2 such that w 1 has a shorter spatial scale than w 2 . For v 1 syn > 0 and v 2 syn < 0 this has the overall effect of short-range excitation and long-range inhibition, and for v 1 syn < 0 and v 2 syn > 0 inhibition dominates at short distances and excitation at longer distances.

III. INSTABILITY OF THE HOMOGENEOUS STEADY STATE
We require excitation and inhibition to generate Turing patterns, hence we let m = {1, 2} in this section. We will also assume β 1 = 1 (without loss of generality) and define β = β 2 < 1 as the parameter which measures the difference in the spatial scales of w 1 and w 2 . We first recast the system Eqs. (21) and (22) as a first-order evolution equation in the state variables u = (a, b, K 1 , g 1 , K 2 , g 2 ), where a = Re(z), b = Im(z), and K m = (1 + τ m ∂/∂t )g m , and seek stationary and spatially homogeneous states u(x, t ) = u * for all x and t.
We apply a small perturbation of the form u(x, t ) = u * + u(x, t ), where u(x, t ) = Ae λt e ikx , λ ∈ C, k ∈ R, and A ∈ C 6 . Using the identity we obtain, to leading order, ∂ t u(x, t ) = J (k) u(x, t ), where J is the following 6 × 6 (k-dependent) Jacobian, written as a block matrix. We highlight the k-dependence as follows: and we show the other block elements in Appendix A. The complex eigenvalues λ = ν + iω satisfy the characteristic equation where I 6 is the 6 × 6 identity matrix. A homogeneous steady state u * is linearly stable to perturbations e ikx if ν(k) < 0 for all k. By the implicit function theorem, a branch of solutions λ(k) to Eq. (23) touches the imaginary axis when where M = Re(E ) and N = Im(E ). We observe a Hopf bifurcation of the spatially uniform state u * if Eq. (23) holds for ν(0) = 0 and k = 0, i.e., if there is a nonzero solution of Here, p i are scalars which depend on u * and the control parameters of the problem (see Appendix B). Solving Eqs. (25) and (26), for fixed width of the Lorentzian , synaptic coupling strengths κ 1 , κ 2 , synaptic time constants τ 1 , τ 2 , and relative width of the synaptic kernel β −1 , gives the locus of Hopf bifurcations as a function of the synaptic reversal potential v syn and background drive η 0 [ Fig. 2(a), green curve]. We fix the synaptic reversal potentials to be equal and opposite, and define v syn = v 1 syn = −v 2 syn . For simplicity, we set κ 1 = κ 2 and τ 1 = τ 2 . For this choice of parameters, excitation dominates for short-range interactions for v syn > 0 and inhibition dominates for short-range interactions for v syn < 0. Note, also, that the Hopf bifurcation occurs for the same value of η 0 for all v syn . This is a consequence of the choice of equal coupling strengths and time constants. If this balance is disrupted, the Hopf bifurcation will depend upon v syn and for some parameter choices we see several Hopf curves in the (v syn , η 0 ) plane.
The homogeneous steady state u * undergoes a static Turing bifurcation if there exists a nonzero critical wave number k c such that Eqs. (23) and (24) hold for ν(k c ) = 0, ω(k c ) = 0. This leads to the conditions where w i and its derivative depend on k c and β, and the scalars p i , q j i depend on u * and the control parameters of the problem (see Appendix B).
As with the Hopf bifurcation, the locus of Turing bifurcations can be plotted as a function of the synaptic reversal potential v syn and background drive η 0 (Fig. 2, red curve). The curve exhibits a turning point, therefore the system has two Turing bifurcations for sufficiently large values of v syn . If we fix v syn at a value in this region, and ascend the bifurcation diagram by increasing η 0 , the spectrum λ(k) touches the imaginary axis of the (ν, ω) plane with wave number k c (lower Turing bifurcation) [ Fig. 2(c)] and continues to move to the right, implying that u * is unstable to a range of perturbations with wave numbers k ∈ (k 1 , k 2 ). As η 0 is further increased, the spectrum returns to the left-hand plane (upper Turing bifurcation), which restores the stability of u * . Note the value of k c is not necessarily the same on the upper and lower branches. The scenario described above is robust to perturbations in other parameters: we found that changes in τ 1 , τ 1 and do not significantly affect the location of the Turing bifurcations. Increasing κ 2 results in an increase of the gap between the two critical values of η 0 , and hence a larger window of instability, while increasing κ 1 decreases this gap and at a certain point the shape is inverted, such that the unstable region lies predominantly in the η 0 < 0 region of the plot. As κ 1 corresponds to the synaptic strength for the excitatory current for v syn > 0, when it becomes significantly larger than κ 2 excitation dominates even for long-range interactions. As patterns arise from the interplay of excitation and inhibition, an inhibitory external drive (η 0 < 0) is needed to observe Turing patterns in this regime.
A Turing-Hopf instability of the homogeneous steady state occurs if there exists a nonzero wave number k c such that Eqs. (23) and (24) hold with ν(k c ) = 0 and ω(k c ) = ±ω c = 0. This instability, which we will also refer to as dynamic Turing bifurcation, elicits wave trains with wave number k c and phase velocity ω c /k c (near bifurcation).
From Eqs. (23) and (24) we obtain a system of the form where i−2 ω i−2 , and p i , q j i are scalars which depend on the control parameters (see Appendix B). In passing, we note that the characteristic equation now has an imaginary part, resulting in the two conditions, Eqs. (29) and (30). Solving Eqs. (29)-(31) allows us to plot the locus of the Turing-Hopf bifurcations in the (v syn , η 0 ) plane (Fig. 2, blue curve), together with the Hopf and Turing bifurcation curves. Note that the Turing-Hopf curve intersects with the Hopf curve at v syn ≈ 8. The value of k c deceases along the Turing-Hopf curve as v syn is increased, such that at the point where the two curves collide k c = 0 on the Turing-Hopf curve.
The curves in Fig. 2(a) were computed using XPPAUT [35], by continuing a suitable algebraic problem in one of the parameters of the system (v syn ). The equations defining an equilibrium are solved simultaneously with Eqs. (25) and (26), (27) and (28), or (29)-(31) and continued in parameter space. The curves partition the (v syn , η 0 ) space into five sectors (labelled I-V). In addition to sectors where u * is stable (I), u * is unstable to bulk oscillations (II), Turing instabilities (III), and Turing-Hopf instabilities (IV), a fifth sector (V) is generated by the crossing of the Turing and Hopf curves. At the intersection (codimension-2) point, the spectrum of the linearized operator around u * has one zero eigenvalue and two complex conjugate eigenvalues, where the critical wave number k c is nonzero for the zero eigenvalue and equal to zero for the complex conjugate pair [ Fig. 2(e)]. Direct numerical simulations close to the instability confirm the predictions of the linear stability analysis: bulk oscillations are observed in in Sec. I [ Fig. 2(b)] stationary Turing patterns are observed in sector III [ Fig. 2(c)] and both wave trains and standing waves are seen in region IV [ Fig. 2(d)]. In sector V, the simultaneous Turing and Hopf unstable modes compete, resulting in a characteristic complex spatiotemporal pattern, where temporal oscillations develop within each bump of a Turing pattern [ Fig. 2(e)]. The existence of the intersection point of the Turing and Hopf curves and the observation of the exotic spatiotemporal patterns are robust to changes in parameters. In passing we note that, close to onset, we could find only large amplitude patterns of the type shown in Figs. 2(c)-2(e), indicating that the Turing and Turing-Hopf bifurcations are subcritical (as confirmed below by numerical continuation). Interestingly, we also observe dynamic Turing patterns in region II, implying that the Hopf bifurcation does not stabilize the system, instead it creates bistability in this region, where the system supports both bulk oscillations and dynamic global patterns.
We also found other complex spatiotemporal patterns away from bifurcation onset, using direct numerical simulation (Fig. 3). The system supports time-periodic patterns containing structures within bumps [ Fig. 3(a)]. This pattern is observed as we move away from the Turing bifurcation but stay close to the Hopf curve. Structures of this type, in which the bumps of a Turing pattern are periodically modulated in time, were observed initially in sector IV of Fig. 2(a). Spatiotemporal patterns of this form, with modulation at the core, were found only in parameter sets where the Turing and Hopf codimension-2 point is present, as expected. Structures within bump solutions are not seen in standard neural field models. They are however commonly observed in spiking neuron models [36,37], which emphasises that this next generation neural mass model retains information about the underlying spiking model.
In region IV, we also observe an number of interesting patterns, such as the wandering bump [ Fig. 3(b)] and a form of standing wave, where both the width and the height of the bumps is periodically modulated [ Fig. 3(c)]. Numerical simulations indicate that the wandering bump is stable for a large area of parameter space and can coexist with a periodic wave train. The standing wave, however, persists for a long time but ultimately is unstable and transitions to a periodic traveling wave (not shown).

IV. NUMERICAL BIFURCATION ANALYSIS
To study patterns away from bifurcations, we employ numerical bifurcation techniques, which allows us to compute coherent structures, determine their stability, and track their dependence on control parameters. Here, we employ the numerical tool kit developed by Avitabile [38] and employed in the context of standard neural fields in Refs. [39,40]. This tool kit can be used to compute waves and patterns and their stability. We refer the reader to recent reviews on numerical bifurcation analysis for coherent structures [41,42]. First, we rescaled space such that x ∈ [−0.5, 0.5], to highlight the dependence on the scale of the domain size . The equivalent PDE formulation Eqs. (21) and (22) We perform numerical bifurcation analysis for heterogeneous spatially patterned steady states and periodic traveling waves of the system above. We construct the stationary patterns by solving for (z, g 1 , g 2 ) the boundary value problem F (z) + G 1 (z, g 1 ) + +G 2 (z, g 2 ) = 0, with Neumann boundary conditions. Wave trains are solutions to Eq. (32) are found by defining ξ = x − ct for c ∈ R, and setting z( where Z(ξ ), G 1 (ξ ) and G 2 (ξ ) are -periodic. To compute wave trains, we solve for (Z, G 1 , G 2 , c) the boundary value problem c∂ ξ z + (F (Z) + G 1 (Z, G 1 ) + G 2 (Z, G 2 )) = 0, posed on ξ ∈ [−1/2, 1/2] with periodic boundary conditions. The last equation in Eqs. (34) is a standard phase condition [43], where ( Z, G 1 , G 2 ) is a reference template solution, such as one of the solutions obtained via direct simulation. We discretized the differential operators Eqs. (34) using standard differentiation matrices, which are also used to compute linear stability of the coherent structures [38].

A. Turing patterns
We first analyze the stationary patterns seen in Sec. III, using numerical continuation to verify the analytical results, determine the criticality of the Turing bifurcation, and examine the behavior of these solutions away from the onset of the instability. We continued solutions to the boundary value problem Eq. (33) in the parameter η 0 with = 12π ,  Fig. 2(a). v syn = 15, and all other parameters as described in the caption of Fig. 2(a). This corresponds to making a vertical excursion through the two-parameter bifurcation diagram Fig. 2(a). Solving Eqs. (27) and (28) we found two Turing bifurcations at η 0 = −0.648, k c = 0.738 and η 0 = 12.67, k c = 0.969. Hence, we numerically continued patterns states with k = 0.738 (blue) and k = 0.969 (red) (Fig. 4). As expected, the homogeneous steady state bifurcates to patterns at η 0 = −0.648 (blue) and η 0 = 12.67 (red), corresponding to the Turing bifurcations found in Sec. III. The stability of the patterned states, as well as the homogeneous state, was numerically calculated along the continuation branch. Along with the Turing bifurcations, the homogeneous steady state undergoes a Hopf bifurcation at η 0 = 3.298, which matches the value found analytically. The patterned Turing solutions were found to go unstable to a globally oscillating periodic pattern, through a Hopf bifurcation at η 0 = 2.8380 (blue) and η 0 = 5.8546 (red). The patterned solutions are also unstable (through a symmetry breaking bifurcation) to broader patterns at low values of η 0 . As anticipated, the bifurcations from the homogeneous steady state are subcritical, hence we have bistability between a spatial pattern and the homogeneous state and the bistability region occurs in a wide region of parameter space. Further numerical continuation results (not shown) indicate that the region of bi-stability increases as the reversal potential v syn is increased. As v syn decreases toward 0, the static Turing bifurcation points collide, and the patterned states cease to exist, as predicted from the Turing analysis in Sec. III. We found the scenario presented above to be robust to changes in the other parameters.

B. Wave trains
We now shift our focus to wave train solutions originating at a Turing-Hopf bifurcation of the homogeneous steady state. We recall that we find these states as solutions to Eq. (34) with periodic boundary conditions, hence corresponds to the spatial period of the wave train profile. In addition, the phase velocity c of a wave train is accessible from the boundaryvalue problem solution.
As in Sec. IV A, the spatial frequency k c at bifurcation was found by solving Eqs. (29)-(31), and we continued patterns with this wave number by setting = 2π/k c , where k c = 0.739 [Figs. 5(a) and 5(b), green curve). We also continued spatial patterns with k = 2k c (red curve) and k = k c /2 (blue curve), this was achieved by posing the system on domains = π/k c and = 4π/k c , respectively. The wave trains bifurcate from the homogeneous steady state with a nonzero phase velocity, at different values of η 0 [ Fig. 5(b)]. The value for which periodic waves emerge for = 2π/k c corresponds to the Turing-Hopf bifurcation (green dot) found in Sec. III, and the speed is equal to ω c /k c . The wave trains are unstable for low values of η 0 . For the smaller domain size they are unstable to the stable steady state, whereas for the larger domain size they transition to finer patterns, with a smaller spatial wavelength k.
To fully explore the relationship between the wave speed and the spatial period, we fix the value of η 0 = 0 and use as a continuation parameter. This opens up the possibility to trace branches of solutions in the (c, ) plane, and hence, approximate the dispersion curve of the waves [Fig. 5(c)]. We found that wave trains occur in isolas, and therefore only exist for a finite range of spatial periods. For 5.2 < < 13.1, a fast (stable) wave train coexist with slower (unstable) one. Whereas, when 3.2 < < 5.2 and 13.1 < < 28.0 the two waves are unstable. For all other values of we do not see wave train solutions. If η 0 is increased, the dispersion curve is no longer an isola, but rather a monotonically increasing function. Thus, in this regime, the system supports periodic traveling waves for all values of the spatial period, above a threshold value.

C. Fronts
Standard neural field models are known to support traveling fronts (for exponentially decaying kernels), which are travelling waves whose profile connects a uniform highactivity state to another low-activity state [44,45]. It is therefore natural to search for these coherent structures in our new neural field model Eqs. (19) and (20). As in any other nonlocal neural field, the existence of the high-and lowactivity states depends on the choice of the synaptic kernel. As inhibition is not necessary for the generation of localized patterns, such as traveling fronts, in this section we consider only a single synaptic conductance, m = 1 and suppress the m label, leading to the equivalent PDE formulation which we pose on ξ ∈ R. We now set z( solutions U (ξ ) = [U 1 (ξ ), . . . , U 6 (ξ )] to the boundary-value problem where N : R 6 → R 6 is the real-valued nonlinear function In this spatial-dynamical system formulation of the problem, the first three components of U have a direct interpretation in terms of the state variables (z, g) of Eq. (36), whereas U 4 , U 5 , and U 6 are auxiliary variables, necessary to cast the problem as a system of first-order differential equations in ξ . The equilibria U ± = (U ± 1 , U ± 2 , U ± 3 , U ± 3 , 0, 0) of the spatial-dynamical system Eq. (37) correspond to highand low-activity homogeneous steady states of Eq. (36) and are completely determined by solving for (U 1 , U 2 , U 3 ) the algebraic problem We define U + as the high activity state, which displays high synaptic conductance and U − as the low activity state, displaying lower synaptic conductance. Note that we also seek the symmetric counterpart solution where lim ξ →±∞ U (ξ ) = U ± . We have continued solutions to Eq. (38) in η 0 and v syn using XPPAUT [35] [ Fig. 6(a)]. The system has three fixed points in the region enclosed by the saddle-node curves, two of which are stable (U ± ). Therefore, we look for traveling waves in this shaded region of parameter space. We note that, as τ is decreased, a Hopf bifurcation of the U − state arises, opening up the possibility of creating heteroclinic connections to periodic orbits, which we will not consider further here.
The state at ξ → −∞ (ξ → ∞) displays a highconductance (low-conductance) G, hence a high (low) firing, and it is therefore referred to as the high (low) activity state [ Fig. 6(b)]. We computed traveling waves using the routines from Ref. [38], which allows us to study the spectral stability of the waves. Computations are performed on a large truncated domain = 60, and plotted on [−10, 10] for convenience. Interestingly, we see ripples in the wake of the front, which indicates that the solution connects a node to a focus. The solution in the phase space (A, B, G) illustrates that the high-activity state (U − ) is a focus, and the low activity one (U + ) a node. We continued the traveling wave shown in Fig. 6(b) (blue), and its symmetric counterpart (red), in the mean background drive η 0 , using c as a solution measure [ Fig. 6(d)]. Fronts live on an isolated branch and destabilise at saddle-node bifurcations. The bifurcation diagram is symmetric with respect to the axis c = 0, as solutions on the branch with c > 0 and c < 0 are related via the transformation U (ξ ) → U (−ξ ). Up to six coexisting waves exist in a wide region of η 0 parameter space, albeit only two of them are stable. In addition, we found stable waves in which the high activity state U + is moving across the tissue, invading the low activity state U − , and vice versa. The inversion of velocity occurs where the blue and red stable branches overlap.
As the isola of traveling fronts is traced, these solutions gain or lose oscillations in the wake of the wave. To investigate further this aspect, we monitor the spectrum of U ± [as equilibria of the spatial-dynamical system Eq. (37)] as we move along the branch. Taking a vertical excursion in Fig. 6(d), we find three independent solutions, one of which is stable and two of which are unstable. We examine the stability of the fixed points U ± in the traveling wave frame for each of these solutions. As expected, oscillations in the wake of the wave are absent where the unstable spectrum of U − is purely real [ Fig. 7(a)], and they develop when a complex-conjugate pair of eigenvalues crosses the imaginary axis [ Fig. 7(b)]. We also find that the amplitude of the oscillations is small when the real part of the complex eigenpair is small [ Fig. 7(c)].

V. DISCUSSION
We have presented the derivation of an atypical neural field model from a network of spatially distributed θ -neurons. In the reduced model, which we dub a next generation neural field model, within population synchrony drives the population firing rate. This model supports a range of patterns, such as bumps, waves and breathers. Noteworthy is the state characterized by structures within bumps, as these states are not seen in standard neural mass models. These structures instead typify patterns seen in networks of spiking neurons [36,37], which signifies that by maintaining the notion of within population synchrony this neural field model can retain information about the underlying spiking network. Exotic states, with within-bump oscillations, were found in the region where the Hopf and Turing bifurcations collided.
A Turing instability analysis provided us with an understanding of how the system behaved close to bifurcation points, allowing us to determine when the system transitioned from the homogeneous steady state. However, unlike the Amari model we cannot use the Heaviside approximation to make further analytical progress since the firing rate is now a fixed real valued function of the Kuramoto order parameter. As such, we have moved to numerical continuation techniques to analyze the behavior of the system away from these bifurcation points. Numerical techniques were also used to examine the existence and stability of travelling fronts.
Previous work [25] illustrated that the point version of this model could support β-rebound, an event-related modulation of the β rhythm, as seen in MEG. This work showed that an external time-dependent drive desynchronized the population and when switched off the population synchrony peaked before returning back to its steady-state level. These changes in population synchrony resulted in modulation of the β rhythm, similar to modulations observed in neuroimaging data. An interesting extension of this work would be to include a similar spatially dependent drive in the neural field model presented here and examine how the inclusion of space affects β-band modulations, and perhaps explain why β rebound is seen in both hemispheres during movement of a single side of the body. More generally, the model parameters can be altered so that the population oscillates at other frequencies, and hence, used to explain other event-related desynchronization and synchronization phenomena in the brain.
In Sec. IV, we pointed to the existence of a front which connects periodic orbits to nodes and focuses, for the model with an exponential coupling kernel. We observed such fronts by decreasing the synaptic time constant τ (Fig. 8). The numerical machinery used here does not allow for the continuation of these solutions, known as defects. The analysis of defects is still an open problem. Close examination of Fig. 8 reveals that there are two fronts in the connection between the node and the limit cycle, which appear to be moving at different speeds. Even more interesting, would be the analysis of fronts which connect two periodic orbits of different amplitudes. The spreading of such a wave across the cortex could be viewed as the spreading of an epileptic seizure. Another numerical challenge would be to continue the exotic patterned states seen in Sec. III when the Turing and Hopf bifurcations collide. These patterns have both a spatial FIG. 8. Travelling wave connecting a periodic orbit to a node: Surface plots showing the evolution of (a) the synaptic conductance g and (b) the synchrony R for a front which connects an oscillatory state to a fixed point state. Simulations for the system defined by Eq. (36), with η 0 = −5, v syn = 10, = 0.5, κ = 5, τ = 0.2 and a temporal period, which would require extending the numerical machineries [38] to continue both a spatial and a temporal pattern. This has been achieved in Ref. [46] for the Brusselator model.
A natural extension to the work presented in Secs. III and IV would be to include a second spatial dimension. It is more natural to view the cortex as a two dimensional sheet and examine the propagation of waves across it. We would expect that the 2D system supports the two dimensional versions of the patterned states presented here, but also potentially some more exotic states. Extending both the Turing analysis and the numerical machinery to include a second spatial extension is worthy of further exploration. To calculate the Jacobian we first write the system Eqs. (19) and (20) in full six-dimensional form,