Mean-Field Models for EEG/MEG: From Oscillations to Waves

Neural mass models have been used since the 1970s to model the coarse-grained activity of large populations of neurons. They have proven especially fruitful for understanding brain rhythms. However, although motivated by neurobiological considerations they are phenomenological in nature, and cannot hope to recreate some of the rich repertoire of responses seen in real neuronal tissue. Here we consider a simple spiking neuron network model that has recently been shown to admit an exact mean-field description for both synaptic and gap-junction interactions. The mean-field model takes a similar form to a standard neural mass model, with an additional dynamical equation to describe the evolution of within-population synchrony. As well as reviewing the origins of this next generation mass model we discuss its extension to describe an idealised spatially extended planar cortex. To emphasise the usefulness of this model for EEG/MEG modelling we show how it can be used to uncover the role of local gap-junction coupling in shaping large scale synaptic waves. Supplementary Information The online version contains supplementary material available at 10.1007/s10548-021-00842-4.


Introduction
The use of mathematics has many historical successes, especially in the fields of physics and engineering, where mathematical concepts have been put to good use to address challenges far beyond the context in which they were originally developed. Physicists in particular are well aware of the "The Unreasonable Effectiveness of Mathematics in the Natural Sciences" (Wigner 1960). One recent breakthrough in the field of large-scale brain modelling has come about because of advances in obtaining exact mean-field reductions of certain classes of coupled oscillator networks via the so-called Ott-Antonsen (OA) ansatz (Ott and Antonsen 2008). This is especially important because the mathematical step from microscopic to macroscopic dynamics has proved elusive in all but a few special cases. Indeed, many of the current models used to describe coarse-grained neural activity, such as the Wilson-Cowan (Wilson and Cowan 1972), Jansen-Rit (Jansen and Rit 1995), or Liley (Liley et al. 2002) model are phenomenological in nature. Nonetheless they have been used extensively to study and explore the potential mechanisms that coordinate brain rhythms underlying cognitive processing and large scale neuronal communication (Fries 2005). For example, such neural mass models have recently been used to understand cross-frequency coupling between brain areas (Jedynak et al. 2015), understand how patterns of functional connectivity may arise in brain imaging studies (Forrester et al. 2020), and are a key ingredient of the Virtual Brain project that aims to deliver the first open simulation of the human brain based on individual large-scale connectivity (Sanz-Leon et al. 2015). Moreover they have been used to uncover how hyper-and hypo-synchrony of neuronal network firing may underpin brain dysfunction including epilepsy (Wendling et al. 2016).
Making use of the OA reduction Luke and colleagues (Luke et al. 2013;So et al. 2014) were able to obtain exact asymptotic dynamics for networks of pulse-coupled theta-neurons (Ermentrout and Kopell 1986). Although the theta-neuron model is simplistic, it is able to capture Handling Editor: Katharina Glomb. This is one of several papers published together in Brain Topography on the "Special Issue: Computational Modeling and M/EEG". some of the essential features of cortical firing pattern, such as low firing rates. As such, this mean-field reduction is a candidate for a new type of cortical neural mass model that makes a stronger connection to biological reality than the phenomenological models mentioned above. The theta-neuron is formally equivalent to the quadratic integrate-and-fire (QIF) model (Latham et al. 2000), a mainstay of many studies in computational neuroscience, e.g. Dipoppa and Gutkin (2013). Interestingly an alternative to the OA approach has been developed by Montbrió et al. (2015) that allows for an equivalent reduction of networks of pulse-coupled QIF neurons, and establishes an interesting duality between the two approaches. In the OA approach the complex Kuramoto order parameter is a fundamental macroscopic variable and the population firing rate is a function of the degree of dynamically evolving within-population synchrony. Alternatively in the approach of Montbrió et al. average voltage and firing rate couple dynamically to describe emergent population behaviour. Given that both approaches describe the same overall system exactly (at least in the thermodynamic limit of an infinite number of neurons) there must be an equivalence between the two macroscopic descriptions. Montbrió et al. have further shown that this relationship takes the form of a conformal map between the two physical perspectives. This correspondence is very useful when dealing with different types of neuroimaging modality. For example, when looking at power spectrograms from electro-or magneto-encephalograms (EEG/MEG), it is useful to contemplate the Kuramoto order parameter since changes in coherence (synchrony) of spike trains are likely to manifest as changes in power. On the other hand the local field potential recorded by an extracellular electrode may more accurately reflect the average population voltage. A model with a perspective on both, simply by a mathematical change of viewpoint, is not only useful for describing experimental data, it may also help the brain imaging community develop new approaches that can exploit a non-intuitive link between seemingly disparate macroscopic variables. Importantly, for this to be relevant to the real world some further features of neurobiology need to be incorporated, as purely pulsatile coupling is not expected to capture all of the rich behaviour seen in brain oscillations and waves. In particular synaptic processing and gap-junction coupling at the level of localised populations of neurons, and axonal delays at the larger tissue scale are all well known to make a major contribution to brain rhythms, both temporal and spatio-temporal (Nunez and Srinivasan 2005;Buzsáki 2011). Fortunately, these biological extensions, that generalise the initial theta-neuron and QIF network models with pulsatile coupling, are natural and easily accommodated. Work in this area has already progressed, e.g. with theoretical work by Laing (2015) and Pietras et al. (2019) on how to treat gap junctions, and by Coombes and Byrne (2019) on the inclusion of realistic synaptic currents (governed by reversal potentials and dynamic conductance changes). Recent work by Byrne et al. (2020) has also considered the inclusion of finite action potential speeds. In this paper we consider a synthesis of modelling work to date on developing a new class of mean-field models fit for use in complementing neuroimaging studies, and present some new results emphasising the important role of local gap-junction coupling in shaping brain rhythms and waves.
Even without the inclusion of gap junctions a first major success of this so-called next generation neural mass and field modelling approach has been in explaining the phenomenon of beta-rebound. Here a sharp decrease in neural oscillatory power in the 15 Hz EEG/MEG beta band is observed during movement followed by an increase above baseline on movement cessation. Standard neural mass models cannot readily reproduce this phenomenon, as they cannot track changes of synchrony within a population. On the other hand the next-generation models treat population coherence as fundamental, and are able to track and describe changes in synchrony in a way consistent with movementrelated beta decrease, followed by an increase above baseline upon movement termination (post-movement beta rebound) (Byrne et al. 2017). Moreover, these models are capable of explaining the abnormal beta-rebound seen in patients with schizophrenia . Beta decrease and rebound are special cases of event related synchrony/desynchrony (ERS/ERD), as measured by changes in power at given frequencies in EEG/MEG recordings (Pfurtscheller and da Silva 1999), and as such this class of model clearly has wider applicability than standard neural mass models that cannot describe ERD/ERS because their level of coarsegraining does not allow one to interrogate the degree of within-population synchrony. By merging this new dynamical model of neural tissue with anatomical connectome data it has also been possible to gain a perspective on whole brain dynamics, and preliminary work in Byrne et al. (2020) has given insight into how patterns of resting state functionalconnectivity can emerge and how they might be disrupted by transcranial magnetic stimulation.
Despite the success of the next generation models that include synaptic processing it is well to recognise the importance of direct electrical communication between neurons that can arise via gap junctions. Without the need for receptors to recognise chemical messengers gap junctions are much faster than chemical synapses at relaying signals. The communication delay for a chemical synapse is typically in the range 1-100 ms, while that for an electrical synapse may be only about 0.2 ms. Gap junctions have long been thought to be involved in the synchronisation of neurons (Alvarez et al. 2002;Bennet and Zukin 2004) and are believed to contribute to both normal (Hormuzdi et al. 2004) and abnormal physiological brain rhythms, including epilepsy (Velazquez and Carlen 2000;Martinet et al. 2017).
In the "Neural Mass Model" section we introduce the mathematical description for the microscopic spiking cell dynamics as a network of QIF neurons with both synaptic and gap-junction coupling. We present the corresponding mean-field ordinary differential equation model with a focus on the bifurcation properties of the model under variation of key parameters, including the level of population excitability and the strength of gap-junction coupling. A simple cortical model built from two sub-populations, one excitatory and the other inhibitory, is shown to produce robust oscillations via a Hopf bifurcation. The derivation of the macroscopic equations of motion is deferred to a technical appendix. This new class of neural mass model is used as a building block in "Neural Field Model" section to construct a continuum model of cortical tissue in the form of an integro-differential neural field model. Here, long-range connections are mediated by action potentials giving rise to space-dependent axonal delays. For computational ease we reformulate the neural field as a brain-wave partial differential equation, and pose it on idealised one-and two-dimensional spatial domains. A Turing analysis is performed to determine the onset of instabilities that lead to novel patterned states, including bulk oscillations and periodic travelling waves. These theoretical predictions, again with details deferred to a technical appendix, are confirmed against direct numerical simulations. Moreover, beyond bifurcation we show that the tissue model can support rich rotating structures, as well as localised states with dynamic cores. Finally, in the "Discussion" section we outline further applications and extensions of the work presented in this paper.

Neural Mass Model
Here we describe a new class of neural mass model that can be derived from a network of spiking neurons. The microscopic dynamics of choice is the QIF neuron model, which is able to replicate many of the properties of cortical cells, including a low firing rate. In contrast to the perhaps more well studied linear or leaky IF model it is also able to represent the shape of an action potential. This is important when considering electrical synapses, whereby neurons directly "feel" the shape of action potentials from other neurons to which they are connected. An electrical synapse is an electrically conductive link between two adjacent nerve cells that is formed at a fine gap between the pre-and post-synaptic cells known as a gap junction and permits a direct electrical connection between them. They are now known to be ubiquitous throughout the human brain, being found in the neocortex (Galarreta and Hestrin 1999), hippocampus (Fukuda and Kosaka 2000), the inferior olivary nucleus in the brain stem (Sotelo et al. 1974), the spinal cord (Rash et al. 1996), the thalamus (Hughes and Crunelli 2007) and have recently been shown to form axo-axonic connections between excitatory cells in the hippocampus (on mossy fibers) (Hamzei-Sichani et al. 2007). It is common to view the gap junction as nothing more than a channel that conducts current according to a simple ohmic model. For two neurons with voltages v i and v j the current flowing into cell i from cell j is proportional to v j − v i . This gives rise to a state-dependent interaction. In contrast, chemical synaptic currents are better modelled with event-driven interactions. If we denote the mth firing time of neuron j by T m j then the current received by neuron i if connected to neuron j would be proportional to ∑ m∈ℤ s(t − T m j ) , where s is a temporal shape that describes the typical rise and fall of a post synaptic response. This is often taken to be the Green's function of a linear differential operator Q, so that Qs = where is a delta-Dirac spike. Throughout the rest of this paper we shall take where H is a Heaviside step function. In this case the operator Q is second order in time and given by where −1 is the time-to-peak of the synapse.
We are now in a position to consider a heterogeneous network of N quadratic integrate-and-fire neurons with voltage v i and both gap-junction and synaptic coupling: Here, firing times are defined implicitly by v j (T m j ) = v th . The network nodes are subject to reset: v i → v r at times T m i . The parameter is the membrane time constant. The strengths of gap-junction and synaptic coupling are v and s respectively. The background inputs i are random variables drawn from a Lorentzian distribution with median 0 and half width . The value of 0 can be thought of as setting the level of excitability, and as the degree of heterogeneity in the network. The larger 0 is, the more neurons would fire if uncoupled, and the larger is, the more dissimilar the inputs are. A schematic of a QIF network and its reduction to a neural field model is shown in Fig. 1, with details of the neural field formulation described in "Neural Field Model" section.
The mean-field reduction of (2) can be achieved by using the approach of Montbrió et al. (2015). This is described in detail in Appendix 1, and is valid for globally coupled cells in the thermodynamic limit N → ∞ . The network behaviour can be summarised by the instantaneous mean firing rate R(t) (the fraction of neurons firing at time t), the average membrane potential V(t) , and the synaptic activity U(t). The synaptic activity U is driven by mean firing rate according to QU = R , with the mean-field dynamical equations for (R, V): Interestingly this (R, V) perspective on population dynamics can be mapped to one that tracks the degree of within-population synchrony described by the complex Kuramoto order parameter Z according to the conformal map (Montbrió et al. 2015): where W * is the complex conjugate of W. The corresponding dynamics for Z is given by equation (23) in Appendix 1. Alternatively, one can evolve the model for (R, V, U) and then obtain results about synchrony |Z| by the use of (5).
The mean field model accurately describes the underlying spiking network (Fig. 2). A network of 1000 synaptically and electrically coupled QIF neurons (blue), as described by (2), was simulated and compared to the mean field dynamics (red), as described by (3) and (4). The finite size fluctuations are most apparent for the membrane potential V. However, the overall behaviour is similar. As expected, increasing the population size reduces the finite size fluctuations.
A previous instance of this model, without gap-junction coupling and with synaptic reversal potentials, was applied to describe beta-rebound, as seen in real MEG data (Byrne et al. 2017). Beta-rebound is a special case of event-related desynchronisation and synchronisation, whereby power in the beta band decreases at movement initiation and rebounds above baseline after movement termination. For our model, which does not incorporate synaptic reversal potentials, we find that gap-junction coupling is important for betarebound. In particular, our results suggest that there is a delicate balance between too little gap-junction coupling and too much gap-junction coupling (Fig. 3). A temporally filtered square pulse of length 400 ms and magnitude 3 μA was applied to the model at t = 0 ms to mirror a movement.
For intermediate values of gap-junction coupling, there is a reduction in beta power at movement onset (0 ms), followed by a sharp increase in power shortly after movement termination (500 ms). The transient increase in beta band power presents with high within population synchrony, confirming the link between rebound and synchronisation. With weak gap-junction coupling, the system does not oscillate, and as such, beta rebound is not possible. With strong gapjunction coupling, the system oscillates, but after movement termination the system returns, almost immediately, to its original behaviour. The population is overly synchronised at steady state, and as such, a transient of high synchrony is not possible.

Single Population: Bifurcation Analysis
We first consider a single excitatory population ( s > 0 ). In contrast to a scalar rate model with self-feedback (as an exemplar of a single population model), the next generation model has at least two variables (to describe either synchrony Z, or the pair (R, V)) and thus is of high enough dimension to support oscillations in time (Fig. 4). Fig. 1 Model schematic. At each point in a two-dimensional spatial continuum there resides a density of QIF neurons whose mean-field dynamics are described by the triple (R, V, U), where R represents population firing rate, V the average membrane potential, and U the synaptic activity. The non-local interactions are described by a kernel w, taken to be a function of the distance between two points. The space-dependent delays arising from signal propagation along axonal fibres are determined in terms of the speed of the action potential c Examining the profile of these oscillations, we observe that the peaks and troughs of the firing rate R and the synchrony |Z| roughly coincide. This indicates, rather unsurprisingly, that when a population is highly synchronised the population firing rate will be high.
As the strength of gap-junction coupling v is decreased the system undergoes a Hopf bifurcation and oscillations disappear (Fig. 5). Note that to the right of the Hopf bifurcation the amplitude and frequency of the oscillations increases with v . Increasing the level of excitability 0 also leads to oscillatory behaviour. A continuation of the Hopf bifurcation in v and 0 is shown for different values of (Fig. 5). The system oscillates for parameter values to the right of these curves. Remembering that sets the level of heterogeneity, we note the window for oscillations gets smaller as the heterogeneity of network is increased.

Excitatory-Inhibitory Network: Bifurcation Analysis
The single population model can be easily extended to a two population network, consisting of an excitatory and an inhibitory population, labelled by E and I respectively. Synaptic coupling is present both within and between populations, while gap-junction coupling only exists between neurons in the same population. The augmented system of equations, describing the mean firing rate R and the average membrane potential V of each population, as well as 4 distinct synaptic variables U for each of the synaptic connections, is presented in Appendix 2.
The excitatory-inhibitory network possesses a rich repertoire of dynamics. For example, it is possible to generate bursts of high frequency and high amplitude activity at a slow burst rate (Fig. 6). This pattern of activity is typical in epileptic seizures, e.g. (Krishnan et al. 2013). Decreasing the gap-junction coupling strengths E v and I v results in smoother lower amplitude oscillations, more in line with healthy brain oscillations. We note that E v and I v are not the only parameters that can change the profile of the oscillations; reducing E 0 (the median background drive to the excitatory population) can also eradicate the seizure-like oscillations.
Next we examine the bifurcation structure of the excitatory-inhibitory network for different combinations of gap-junction coupling strengths E v and I v (Fig. 7). With no gap-junction coupling in either population ((a) , the amplitude of oscillation increases significantly for this branch of oscillatory solutions. An additional branch of oscillatory solutions emerges for low I 0 , with moderate amplitude oscillations for the firing rate of the inhibitory population R I and low amplitude oscillations for the excitatory population R E . Interestingly, the two oscillatory solutions co-exist for I 0 ≈ −7.5 to 2.5. Jansen and Rit (Jansen and Rit 1995) demonstrated that transitions between seizures and healthy brain activity could be  With a good understanding of the behaviour of the spatially clamped system, we move on to consider the spatially extended neural field model.

Neural Field Model
Brain waves are inherently dynamical phenomena and come in a vast variety of forms that can be observed with a wide range of neuroimaging modalities. For example, at the mesoscopic scale it is possible to observe a rich repertoire of wave patterns, as seen in voltage-sensitive dye imaging data from the primary visual cortex of the awake monkey (Muller et al. 2014), and local field potential signals across the primary motor cortex of monkeys (Rubino et al. 2006). At the whole brain scale they can manifest as EEG alpha oscillations propagating over the scalp (Hindriks et al. 2014), and as rotating waves (defined as a significant increase in phase offset with rotation about a wave center) seen during human sleep spindles with intracranial electrocorticogram recordings (Muller et al. 2016). Waves are known to subserve important functions, including visual processing (Sato et al. 2012), saccades (Zanos et al. 2015), and movement preparation (Rubino et al. 2006). They can also be associated with dysfunction and in particular epileptic seizures (Martinet et al. 2017). Computational modelling is a very natural way to investigate the mechanisms for their occurrence in brain tissue, as well as how they may evolve and disperse (Heitmann et al. 2013(Heitmann et al. , 2017Liou et al. 2020).
The study of cortical waves (at the scale of the whole brain) is best advanced using a continuum description of neural tissue. The most common of these are referred to as neural fields, and are natural extensions of neural mass models to incorporate anatomical connectivity and the associated delays that arise through wiring up distant regions using axonal fibres. The study of waves, their initiation, and their interactions is especially pertinent to the study of epileptic brain seizures and it is known that gap junctions are especially important in this context (Martinet et al. 2017). Phenomenological neural field models with gap-junction coupling have previously been developed and analysed by Steyn-Ross et al. (2007, and more principled ones derived from theta-neuron models by Laing (2014Laing ( , 2015. In the latter approach it was necessary to overcome a technical difficulty by regularising the shape of the action potential. However, with the approach used in "Neural mass model" section this is not necessary and the neural field version of (3)-(4) is constructed by replacing full temporal derivatives by partial temporal derivatives and replacing the temporal dynamics for U with the dynamics QU = Ψ , where Ψ denotes a spatio-temporal drive. For example, in the plane we might consider where R( , t) is the population firing rate at position ∈ ℝ 2 at time t and c represents the speed of an action potential, as illustrated in Fig. 1. Typical values for cortico-cortical axonal speeds in humans are distributed, and appear to peak in the 5 − 10 m/s range (Nunez and Srinivasan 2014). Here, w represents structural connectivity as determined by anatomy. For example, long-range corticocortical interactions are predominantly excitatory whilst inhibitory interactions tend to be much more short-ranged, suggesting a natural choice for the shape of w as an inverted Mexican hat (Stepanyants et al. 2009). A similar equation would hold in one spatial dimension. We emphasise that in the continuum model presented here, the gap-junction coupling has no spatial extent. The mass model (defined locally) incorporates gap junctions, while the only coupling between masses is via synaptic currents. The model in Laing (2014) treats a linear array with nearest neighbour electrical interactions (representing cells that touch) as well as allowing for interactions beyond nearest neighbour. For large scale brain modelling it is more natural to view the brain as a network of synaptically interacting neural masses, each with its own local synaptic and gap-junction currents, with longer range interactions mediated only by synaptic currents.
In this section we shall work with the explicit choices of structural connectivity w(x) = (|x| − 1)e −|x| in 1D and w(r) = (r∕2 − 1)e −r ∕(2 ) in 2D (where x represents distance in 1D, and r represents radial distance in 2D). For convenience we have chosen spatial units so that the scale of exponential delay is unity, though note that typical values for the decay of excitatory connections between cortical areas (at least in macaque monkeys) is ∼ 10 mm (Markov et al. 2010). Both of the above kernel shapes have an inverted wizard hat shape and are balanced in the sense that the integral over the whole domain is zero. They also allow for a reformulation of the neural field model as a partial differential equation, as detailed in Appendix 3. The resulting brain-wave equation is very amenable to numerical simulation using standard (e.g. finite difference) techniques. Before we do this, it is first informative to determine some of the patterning properties of the neural field model using a Turing instability analysis. Below we outline the results of the analysis and discuss the ensuing patterns for the neural field model in both 1D and 2D.

One Spatial Dimension
Turing instability analysis, originally proposed by Turing in 1952(Turing 1952, is a mechanism for exploring the emergence of patterns in spatio-temporal system, including neural fields. Similar to the bifurcation analysis for the neural mass model, it allows us to determine the parameter values for which oscillations and patterns occur. Bulk oscillations, whereby synchronous activity across the spatial domain varies uniformly at the same rate, emerge at a Hopf bifurcation. Static patterns, which do not change with time, emerge at a Turing bifurcation. Dynamic patterns, that oscillate in time and space, emerge at a Turing-Hopf bifurcation. The 1D neural field model, given in Appendix 3 by (38), supports both bulk oscillations and spatio-temporal patterns. Using the inverted wizard hat connectivity kernel (longrange excitation and short-range inhibition), we find Hopf and Turing-Hopf bifurcations (Fig. 8 left). See Appendix 4 for details of the analysis. For the chosen parameter values and weak gap-junction coupling ( v ≲ 0.8 ), the spatiallyuniform steady state is always stable and neither patterns nor oscillations exist. Increasing the median background drive 0 moves the Hopf and Turing-Hopf curves down in the c-V plane, allowing for oscillations and patterns in the absence of gap junctions ( v = 0 ). For slow action potential speeds ( c ≲ 0.2 ), the system first undergoes a Hopf bifurcation as v is increased and bulk oscillations emerge (Fig. 8I). As v is increased further, the system undergoes a Turing-Hopf bifurcation and standing waves emerge (Fig. 8II). For faster action potential speeds ( c ≳ 0.2 ), the Turing-Hopf bifurcation occurs before the Hopf, and we see periodic travelling and standing waves between the two bifurcations (Fig. 8III).
To assess the role of gap junctions, we fixed the action potential speed c = 0.11 and explored the dynamics of the synchrony variable |Z| for different gap-junction coupling strengths v (Fig. 9). For weak gap-junction coupling (I), there is a regular standing wave and the level of synchronisation is low. As v is increased bulk oscillations emerge and the level of synchronisation increases (II). Increasing v further leads to the emergence of mixed dynamics, with both spatial and temporal patterning. The tissue is now highly synchronised, confirming the belief that gap-junction coupling increases the level of synchronisation.
For a standard wizard hat coupling kernel (long-range inhibition and short-range excitation) the neural field model can undergo a Turing bifurcation, as well as Hopf and Turing-Hopf bifurcations (see Supplementary material 1 panel (a)). Changing the sign of the synaptic coupling strength s changes the coupling to long-range inhibition and short-range excitation. When Turing and Hopf instabilities occur simultaneously, interesting patterns emerge. In particular, we see stationary bumps where the activity at the centre of the bump oscillates in both space and time (see Supplementary material 1 panel (b)). We will discuss the two dimensional version of such patterns in more detail below.

Two Spatial Dimensions
A Turing analysis was also performed for the 2D neural field equation, given in Appendix 3 by (37), and a very similar bifurcation structure was found (see Supplementary material 2 panel (a)). As expected, close to the Hopf bifurcation the activity of the tissue oscillates in time, but no spatial pattern emerges (see Supplementary material 3). Near the Turing-Hopf bifurcation we see both planar waves (Supplementary material 4) and radial waves (Supplementary material 5) depending on initial conditions. Close to the intersection of the Turing-Hopf and Hopf bifurcation we see mixed spatio-temporal dynamics (Supplementary material 6). Away from bifurcation, more interesting patterns emerge.
We fix the action potential speed c = 1 and vary the gapjunction coupling strength v to assess how gap-junction coupling affects patterning. For weak gap-junction coupling, we observe rotating waves with source and sink dynamics where the waves collide with each other (Fig. 10). The domain shown contains 12 rotating cores. Periodic boundary conditions were used. Hence, the cores at the edge of the domain wrap around to those on the other side. Supplementary material 7 shows the temporal evolution of the synchrony variable |Z|, from which the cores and rotations are readily observed. The direction of rotation alternates, such that every second core rotates clockwise/anti-clockwise.
As the gap-junction coupling strength v is increased robust spirals emerge at the centre of the rotating cores. The spiral is tightly wound with a diffused tail of high amplitude activity that propagates into the rest of the domain and interacts with the other rotating waves (Fig. 11). The time course of a point close to the centre of a rotating core (green dot) depicts higher amplitude oscillations for the firing rate R, mean membrane potential V and level of synchronisation |Z| when compared to the simulations for lower gap-junction coupling strength v (Fig. 10). In addition, the peaks in R are sharper and the minimum level of synchrony |Z| is substantially higher. The temporal evolution for the full tissue can be seen in Supplementary material 8.
We again note that increasing the gap-junction coupling strength increases the level of synchronisation across the tissue. For v = 0.695 , the synchrony variable oscillates between 0.02 and 0.36. For v = 0.8 , it oscillates between .120 and 0.56. Increasing v further does not change the overall dynamics, but does result in higher levels of synchronisation. This supports the hypothesis that gap-junction coupling lends to more synchronous activity.
As mentioned in the "One Spatial Dimension" section, for a regular wizard hat connectivity kernel (short-range excitation and long-range inhibition) the neural field supports static Turing patterns, periodic bumps of high activity in 1D and a periodic lattice of high-activity spots in 2D. When a Hopf instability coincides with this Turing instability, patterns form at the centre of these localised states. In 2D, patterns of concentric circles can appear within spots when the two bifurcations coincide (Fig. 12). Activity within a localised state can oscillate in time, while the activity in the surround is constant with a low firing rate. These patterns are reminiscent of chimeras (Abrams and Strogatz 2004;Kuramoto and Battogtokh 2002;Laing 2009a, b;Omelchenko et al. 2011), as seen in networks of coupled oscillators, where a fraction of the oscillators are phase-locked or silent while the others oscillate incoherently. Note how the peaks in firing rate coincide with peaks in synchrony. However, in the surround synchrony is high, but the firing rate is minuscule. This indicates that the neurons are also synchronised at rest. A video illustrating how these exotic patterns evolve on the entire spatial domain is provided in Supplementary material 9 and the bifurcation diagram is given in Supplementary material 2 panel (b). The patterns presented here persist with the refining of the numerical mesh so we are confident that the relatively sharp changes between spatial points are not just a numerical artefact.
When not perfectly synchronised, the relative timing of oscillations across single areas or distant regions in the cortex can give rise to a range of flexible phase offsets which can manifest as travelling waves of various shapes, including plane, radial and spiral waves (Muller et al. 2018). The numerical simulations presented above highlight the ease with which these can be generated within the next generation neural field model with local gap-junction currents, and in particular spiral waves. The latter are thought to be highly relevant to status epilepticus characterised by the formation of spiral waves that emerge after wavefront annihilation and exhibit complex interactions (Liou et al. 2020).

Discussion
Mean-field models have proven invaluable in understanding neural dynamics. Although phenomenological in nature, coarse-grained neural mass/field models have proven particularly useful in describing neurophysiological phenomena, such as EEG/MEG rhythms (Zhang 1996), cortical waves (Wilson et al. 2001;Roberts et al. 2019), binocular rivalry (Laing and Chow 2002;Bressloff and Webber 2012), working memory  and visual hallucinations (Ermentrout and Cowan 1979;Bressloff et al. 2001). The exclusion of synchrony in standard neural mass/field models prohibits them from describing event-related synchronisation and desynchronisation; the increase and decrease of oscillatory EEG/MEG power due to changes in synchrony within the neural tissue. Here we presented and analysed a recently developed neural mass/ field model that incorporates within population synchrony. In contrast to other reductive approaches for describing the behaviour of populations of spiking neurons the one described here is exact (in the thermodynamic limit) for realistic event-driven models of (non-instantaneous) synaptic process. For example, the spike-density formalism for reducing networks of linear integrate-and-fire neurons requires a moment closure approximation (Ly and Tranchina 2007), whilst Fokker-Planck approaches for describing renewal-type spiking neurons often only reduce after the truncation of some eigenfunction expansion (Pietras et al. 2020). The mean-field model presented here has previously been applied to real MEG data in Byrne et al. (2017) (at the neural mass level) and used with MRI-derived structural connectivity in Byrne et al. (2020) (for a network of neural masses), though not with the inclusion of gap junctions. The main benefit of such a model is that it is derived from a population of interacting spiking neurons, with the QIF model incorporating a reasonable representation of the action potential shape. This further allows for the inclusion of realistic gap junctions at the cellular level. Gap junctions are known to promote synchrony within neural tissue (Watanabe 1958;Bennett 1977) and the strength of these connections has been linked to the excessive synchronisation driving epileptic seizures (Mylvaganam et al. 2014;Volman et al. 2011). Nonetheless, it is also important to recognise the important effects that the extracellular space has on seizure dynamics, as discussed in Wei et al. (2014). Recent work by Martinet et al. (2017) has emphasised the usefulness of bringing models to bear on this problem, and coupled the Steyn-Ross neural field model (Steyn-Ross et al. 2013) to a simple dynamics for local extracellular potassium concentration. Here, gap junctions are modelled by appending a diffusive term to a standard neural field and increases in the local extracellular potassium concentration act to decrease the inhibitory-to-inhibitory gap-junction diffusion coefficient (to model the closing of gap junctions caused by the slow acidification of the extracellular environment late in seizure). A more refined version of this phenomenological approach would be to replace the Steyn-Ross model with the neural field described here. This would allow a more principled study of how slow changes in the extracellular environment could initiate wave propagation, leading to waves that travel, collide, and annihilate. Indeed, simulations of the next-generation neural field model (without coupling to the extracellular space) have already shown such rich transient dynamics including seizure-like oscillations (and their dependence on the strength of gap-junction coupling). It would be interesting to explore this further, and in particular the transitions whereby spatio-temporal wave patterns are visited in sequence. This has already been the topic of a major modelling study by Roberts et al. (2019) who considered a variety of more traditional neural mass models in a connectome inspired network using the 998-node Hagmann et al. dataset (Hagmann et al. 2008) with a single fixed axonal delay. A similar computational study, with a focus on spiral waves and sinks/sources from which activity emanates/converges, could also be undertaken using the alternative neural mass model presented here, and with the further inclusion of space-dependent axonal delays. Moreover, electrical stimulation can easily be integrated into the model, by returning to the microscopic voltage dynamics model given by (2) (which ensure current balance) and including a time-dependent drive, say A(t), which could represent a pattern of applied transcranial direct current. This modifies the background drive in the mean-field model according to 0 → 0 + A(t) . In Byrne et al. (2020) this approach was used to determine the effects of transcranial magnetic stimulation (with an induced electrical form for A(t)) on patterns of network functional connectivity. Finally, it is well to note the assumption throughout our modelling study that chemical and electrical synapses operate independently. However, there is now accumulating evidence to suggest that this might not be the case (Pereda 2014). For example, neurotransmitter modulators released by nearby synaptic terminals can regulate the synaptic strength of co-localised chemical and electrical synapses through the activation of G protein-coupled metabotropic receptors. All of the above are topics of ongoing investigation and will be reported upon elsewhere.

Appendix 1: Mean-Field Reduction
Consider a heterogeneous network of N quadratic integrate-and-fire neurons with voltage v i and both gap-junction and synaptic coupling: Here, the mth firing time of the jth neuron is defined implicitly by v j (T m j ) = v th . The network nodes are subject to reset: v i → v r at times T m i . The strengths of gap-junction and synaptic coupling are v and s respectively. The function s(t) represents the shape of a post synaptic response (to a delta-Dirac spike) and will be taken to be the Green's function of a linear differential operator Q. For an alpha-function s(t) where H is a Heaviside function, Q = (1 + −1 d∕dt) 2 , whilst for an exponential response s(t) = exp(− t)H(t) , Q = (1 + −1 d∕dt) . In (7) the i are random variables drawn from a Lorentzian distribution: with median 0 and half-width . The threshold v th and reset v r values are taken to be ∞ and −∞ , respectively.
To derive the mean-field equations we follow closely the exposition by Montbrió et al. (2015). Consider the 1 ( − 0 ) 2 + 2 , thermodynamic limit N → ∞ with a distribution of voltage values ( | , t) . The continuity equation for is where and which represent the average voltage and population firing rate respectively. We now assume a solution (v| , t) of the form For a fixed the firing rate r( , t) can be calculated as (v → ∞| , t)̇v(v → ∞| , t) , from which we may establish that By exploiting the structure of (14), with poles at v ± = y ± ix , a contour integration shows that where PV denotes the Cauchy principal value. After averaging over the distribution of single neuron drives given by (8) we obtain For fixed , substitution of (14) into the continuity equation and balancing powers of v shows that x and y obey two coupled differential equations that can be written as .
where ( , t) = x( , t) + iy( , t) . After evaluating the integrals in (17) and (18) using contour integration (and using the fact that has poles at ± = 0 ± i ) the coupled equations for (R, V) can be found as The complex quantity W = R + iV is known to be related to the Kuramoto order parameter Z by the conformal map (Montbrió et al. 2015): where W * is the complex conjugate of W. The evolution equation for Z is given by the complex differential equation where QU = R(Z) and

Appendix 2: Interacting Sub-populations
Consider an excitatory population labelled by E coupled to an inhibitory one labelled by I. In this case there are four distinct synaptic inputs with connection strengths ab s , a, b ∈ {E, I} , with aE s > 0 and aI s < 0 . Each population has a background drive drawn from a Lorentzian with median a 0 and half-width a , a ∈ {E, I} . Generalising the mean-field model derived in section Appendix 1, for gap-junction coupling only within a given sub-population, we have that For a second order synapse with time-scale −1 ab we would set Note that in a slightly more general setting that would allow for electrical connections between excitatory and inhibitory sub-populations then we would replace (26) and (27)

Appendix 3: Brain Wave Equation
A simple continuum model for an effective single population dynamics can be written in the form where Ψ = w ⊗ R . The symbol ⊗ is used to describe spatial interaction within the neural field model, while w represents structural connectivity. For example, in the plane we might consider (R, V, U) = (R( , t), V( , t), U( , t)) , with ∈ ℝ 2 and t ≥ 0 with where c represents the speed of an action potential. We note that (35) can be written as a convolution: where G(r, t) = w(r) (t − r∕c) . For certain choices of w it is possible to exploit this convolution structure to obtain a PDE model, often referred to as a brain-wave equation (Nunez 1974;Jirsa and Haken 1997