Coombes, Stephen (2009) Large-scale neural dynamics: simple and complex. NeuroImage . ISSN 1053-8119 (In Press)

,


Introduction
Systems of equations that can model the brain at the very large scale are becoming increasingly important for underpinning experimental techniques including those of EEG, fMRI, MEG and optical imaging using voltage sensitive dyes.All of these can reveal patterns of spatio-temporal activity that span centimetres of tissue in the cerebral cortex.One brute-force approach to building models that can help make sense of these rich dynamical patterns is to consider the detailed properties of synapses, neurons, microcircuits and cortical columnar organisation and build vast computational compartmental models [1].Although such an approach may bring one close to the neurobiology, the complexity of such models is a severe barrier to gaining insight into emergent neurodynamics.
An alternative approach is to consider coarse-grained cortical units, called neural masses, and use these as building blocks for cortical tissue models.Such an approach has two main advantages.Firstly, as opposed to the brute-force approach there is a drastic reduction in the dimensionality of both the parameter and variable space.Secondly, since detailed knowledge of single neuron properties (ionic currents and dendritic structure) is not required, much of the experimental data (on connectivity, neurotransmitter types, distribution of axonal speeds) needed to constrain such models is already available.The current mathematical approach for understanding coarse-grained activity of large ensembles of neurons in cortex is based around the work of Wilson and Cowan [2,3], Amari [4,5] and Nunez [6] in the 1970s.Because the number of neurons and synapses in even a small piece of cortex is so vast a natural first approximation is to take a continuum limit and treat cortical space as continuous, giving rise to the notion of a neural field model.These models typically take the form of integro-differential equations.Their non-local nature (arising from long range axonal connections) has led to the development of a set of analytical and numerical tools for the study of waves and patterns based around natural extensions of those used for local partial differential equa-tion (PDE) models.Indeed they have been used in a number of neural contexts including understanding mechanisms for short term working memory [7], motion perception [8], representations in the headdirection system [9], and feature selectivity in the visual cortex [10].For a review of such models in both one and two dimensions that can incorporate realistic forms of axo-dendritic interactions we refer the reader to [11,12].Albeit a framework that is limited to coarse-grained or mean-field activity, neural field models provide a direct connection from neural activity to EEG and fMRI data [13,14] (unifying data from different imaging modalities) as well as a providing a bridge to cognitive theories of brain function [15].Moreover, dynamic causal modelling (DCM) (see [16] for a review), which is frequently invoked for the interpretation of fMRI data, is now being extended from a data-driven perspective to incorporate activity models based upon neural field equations [17].In light of the recent and rapid advances in the imaging of large scale cortical dynamics it is thus timely to review how neural field models can underpin empirical research that emphasises brain structure, dynamics and function.
We begin this overview of the practical uses of neural field models by first describing their basic building block, namely a neural mass model of a homogeneous neuronal population (which we may loosely think of as being a part of a cortical column).Next we describe the dynamics of an interacting set of neural masses (now at the scale of a whole cortical column) and explore their dynamics using bifurcation theory to uncover the natural time-scales for emergent rhythms.[19].In two spatial dimensions we show how the full non-local dynamics of a neural field model can be approximated with a local PDE model and consider its extension to treat patchy connections of the type that arise when isotropic connectivity is periodically modulated.We also discuss the effects that more general inhomogeneous connectivities can have on wave propagation through cortex, using mathematical analysis to emphasise the conditions for wave propagation failure.
As an exemplar of the gains to be made with coarsegrained modelling we report on recent work of Bojak et al. [13] that utilises neural field models with realistic anatomical and physiological parameters for a folded cortex and a realistic head model to predict EEG and fMRI responses.Finally we discuss future directions for the mathematical descriptions of neural tissue relevant to neuroimaging.

Population modules
It is common practice to define a neural mass as a collection of thousands of near identical interconnected neurons with a preference to operate in synchrony.
The spatial extent of this population is taken to be on the order of a few hundred micrometers.The state variable describing the activity of the population is the average membrane potential.Perhaps the most well known neural mass model is that of Jansen and Rit [20] based on the original work of Lopes Da Silva et al. [21].Such lumped parameter models are capable not only of producing EEG style alpha rhythms but can generate more complex signals ranging from delta to gamma seen in EEG and MEG recordings (by appropriately modifying the population kinetics within a physiologically plausible range) [22].Moreover they are amenable to mathematical analysis using techniques from dynamical systems theory and notably bifurcation theory [23].However, it is well to mention here that because of the "near to synchrony" assumption, neural mass models are unable to say anything about more complex behaviours within a single population, such as phase-locked states (away from synchrony) or clustering.Recent approaches that improve upon this situation have been developed by Stefanescu and Jirsa (using mode decomposition techniques) [24] and Laing and Kevrekidis (using "equation free modelling" and generalised polynomial chaos expansions) [25,26].
One of the more successful population models for generating rhythms consistent with those found in the human EEG spectrum is that of Liley et al. [18].In this mesoscopic model cortical activity is locally described by the mean soma membrane potentials of an interacting excitatory and an inhibitory population.The interaction is through a model of the synapse that treats both shunting currents and a realistic time course for post-synaptic conductance changes.Referring to the diagram in Fig. 1 the model can be written in a succinct form as with a, b ∈ {E, I}, where E (I) is the mean membrane potential in the excitatory (inhibitory) population.The relaxation time constants for the populations are given by τ a , whilst h a describes a reversal potential such that h E (h I ) is positive (negative) with respect to the resting state.The weights W ab are the product of a static strength factor and a dynamic con-ductance W ab = W ab g ab , where Here Q ab represents a linear differential operator: and the conductances are considered to be driven by a combination of firing from populations to which they are connected and some external drive.The former is modelled using a sigmoidal function: where H is a Heaviside step function.We recognise η ab (t) as a delayed α-function, commonly used in computational neuroscience to mimic the rise and fall of a post-synaptic conductance change.
To determine the properties of such a population model it is natural to first consider the stability of the steady state, which we shall denote by a ss (defined by ȧ = 0).Linearising about this fixed point and seeking solutions of the form e λt gives a characteristic equation det E(λ) = 0 where the entries of the 2 × 2 matrix Here ) and η ab (λ) is a Laplace transform given by: It is instructive to consider the limit τ a → 0, α ab = α, κ a = κ and ∆ ab = ∆ in which case the characteristic equation takes the simple form where Γ ± are the two eigenvalues of W .If a pair of complex conjugate eigenvalues λ = ν ± iω crosses the imaginary axis ν = 0 from left to right in the complex plane, then a Hopf bifurcation can occur, leading to the formation of periodic oscillations.This is almost generic in the presence of delays, see for example [27].For an interesting discussion on the role of delays in models for generalised epileptic seizures we refer the reader to Breakspear et al. [28].For zero-delay (∆ = 0) let us suppose that W has a pair of complex conjugate eigenvalues re ±iθ with 0 < θ < π.In this case a Hopf bifurcation occurs when 1 = r/κ cos(θ/2) independent of α with a non-zero frequency ω = α r/κ sin(θ/2).However, in general a Hopf bifurcation will depend on the relative timescales in the model and should be determined as the solution of det E(iω) = 0.In practice it is much easier to numerically determine the stability of fixed points of the full ten dimensional model, defined by equations (1-5), using software such as XPPAUT [29].
Moreover, this allows us to construct bifurcation diagrams like those in Fig. 2, which shows the coexistence of a large and small amplitude periodic orbit, with time-course shown in Fig. 3.The large amplitude, (∼ 5 Hz) orbit has been suggested to correspond to a form of epileptic dynamics, whilst the smaller amplitude (∼ 10 Hz) oscillation is more consistent with the alpha band of the EEG spectrum.Moreover, a period doubling cascade is supported and beyond this the model is known to support chaos.Indeed the chaotic behaviour has been extensively investigated in a series of papers [31,32,30], highlighting that the route to chaos is actually via a novel Shilnikov saddle-node bifurcation.Note that in contrast the chaotic EEG patterns in the Freeman model of the olfactory system are generated via the Ruelle-Takens-Newhouse route [33].
For further discussion on the use of neural mass models in epileptic modelling and human brain rhythms we refer the reader to [34,35] and [36]  Stable (unstable) branches are solid (dashed).Chaotic solutions are found after the period doubling cascade.
Parameters are modified from [30] as P IE = 0.005763Here parameters are as in Fig. 2 with P EE = 0.006216 ductances and we may write E = E(g EE , g EI ) and Symbolically the model now takes the closed form where we suppress indices and use the notation {g} to emphasise that the firing rate depends on the set of network conductances.A whole host of models can be described with such a system of equations, ranging from networks with many more than two modules and a complicated dependence of the firing rate on the set of network conductances down to a single population with self-excitation.A very simple example of this latter case for a single conductance g = u can be obtained for the choice of an exponential synapse η(t) = αe −αt H(t) and assuming f = f (u) with P = 0, so that In some sense we may regard this model as one of the most basic to arise in mathematical neuroscience.
Interestingly it has been subjected to modifications that improve its ability to model neuronal dynamics without recourse to abandoning the fast relaxation assumption.One of the main extensions of this type has been the inclusion of a term to describe refractoriness [2].A particular example is that of Curtu and Ermentrout [37]: where R is the absolute refractory period of the neurons in the population.The fixed point u ss satisfies the equation −u ss + (1 − u ss )f (u ss ) = 0.For a sigmoid with 0 < f < 1 then there is at least one solution of this equation for u ss ∈ (0, 1/2).The characteristic equation determining the linear stability of the steady state is calculated as E(λ) = 0, with where a dynamic bifurcation defined by E(iω) = 0 can occur [37].Interestingly periodic solutions that emerge beyond such a bifurcation can be analysed explicitly in some singular limit.To see this we write (11) as a twodimensional delay differential equation: where we have re-scaled time according to t → t/R and set = (αR) −1 .Formally setting = 0 gives the graph z = 1 − u/f (u), which we show an example of in Fig. 4. Note the 'cubic' shape of this curve, reminiscent of the nullcline for the fast variable seen in many excitable systems, and particularly those of FitzHugh-Nagumo type [38].In the limit → 0 dynamics consists of slow evolution along the left and right branches of the 'cubic' with fast transitions from one branch to another.In this case it can be shown that the period of oscillation satisfies T < 2R (and numerical experiments show that the actual oscillation period scales linearly with R).Hence, the inclusion of refractoriness in activity based models is a natural way to induce oscillations with an emergent period that is largely set by the refractory time-scale.Another extension of the basic model ( 10) is to recognise that in cortical tissue there are an abundance of metabolic processes whose combined effect is to modulate neuronal response.It is convenient to think of these processes in terms of local feedback mechanisms.For example, spike frequency adaptation (SFA) is a property of many single neurons and has been linked to the presence of a Ca 2+ gated K + current, I AHP [39].
The generation of an action potential leads to a small calcium influx that increments I AHP , with the end result being a decrease in the firing rate response to persistent stimuli.A simple phenomenological model of this process is to add a current, that activates in the presence of high activity, I − gs, to the right hand side of (10) [40], where and I is a constant bias.This form of SFA leads naturally to the emergence of periodic behaviour.Indeed in the limit of large τ (slow adaptation) and a for a steep sigmoid (large β) this period can be calculated explicitly as τ log((g − I)(1 + I)/(I(g − 1 − I)) for g > 1 + I and shows a classic relaxation oscillation between up and down states like those seen in biophysical models of slow (< 1 Hz) oscillations [41].An alternative approach, consistent with observations first made by Hill in 1936 [42], is to treat the threshold in the firing function to be state-dependent.This is explored in detail in [43,44] and shown to lead to exotic dynamics at the network level, including the emergence of dissipative solitons.However, rather than pursue these extensions in more detail we shall instead next show how to build tissue level models taking as a starting point models of the form (9).

Tissue models
Here we will view a macroscopic part of the neocortex as being adequately modelled as a spatial assembly of population models -a neural field model.Such a viewpoint has already proven useful in understanding spatial aspects of the alpha rhythm and in particular cortical travelling waves [6,45].For a recent perspective on the use of neural field models in interpreting extrinsic optical imaging data (from in vitro exper-iments on pharmacologically treated brain slices) we refer the reader to [46].
To develop the extension of ( 9) to treat spatially continuous neural sheets (such as a two-dimensional cortex) we will adopt the continuum assumption and treat a density of neurons at a point with inputs that arise from the delayed and weighted contribution of activity at other points in the tissue.Because these interactions are mediated by long-range axonal fibres the resulting tissue-level model is inherently non-local and is often cast in the form of an integral equation.We represent this symbolically in the form where the operator ⊗ captures information about both anatomical connectivity patterns and the distribution of axonal delays.As a concrete example consider a set of two-dimensional interacting layered sheets (with both self and layer to layer interactions), each containing only one cell type (either excitatory or inhibitory).
The activity in layer a induced by that in layer b (generalising equation ( 9)) then takes the form where ψ ab = ψ ab (r, t) is given by Cowan model we refer the reader to [47].Suffice to say that, apart from their relevance to neuroimaging, neural field models have found many applications in neuroscience, including to understanding the generation of visual hallucinations [48,49], modelling orientation tuning in visual cortex area v1 [10], describing travelling waves of activity in v1 during binocular rivalry [50,51], models of working memory [7] and encoding of continuous stimuli [52], motion perception [8], somatosensory illusions [53], and developing a theory of cognitive robotics [54] for example.
Neural field models of the type ( 16) are nonlinear spatially extended systems and thus have all the necessary ingredients to support pattern formation.The analysis of such behaviour is typically performed with a mixture of linear Turing instability theory, weakly nonlinear perturbative analysis and numerical simulations (see [55] for a review).In the absence of detailed anatomical data it is common practice to consider cortico-cortical connectivity functions to be homogeneous and isotropic so that w ab (r, r ) = w ab (|r − r |).In this case a homogeneous steady state is expected and can be defined by h ss a = b W ab f b (h ss b ), where W ab = R 2 drw ab (r).For concreteness we shall take w ab (r) = w 0 ab e −r/σ ab /(2π), where r = |r|.Linearising around the steady state and considering perturbations of the form h a (r, t) ∼ e λt e ik•r , gives an equation for the continuous spectrum λ = λ(k), for k = |k|, in the form E(k, λ) = 0 [56], where E(k, λ) = det(D(k, λ) − I), and Here γ a = f a (h ss a ) and G ab (k, ω) is the Fourier transform of G ab (r, t) = w ab (r)δ(t − r/v ab ) defined by where A ab (ω) = 1/σ ab + iω/v ab .Here we use a notation that distinguishes functions and their transforms simply by their arguments, namely (r, t) for the original space and (k, ω) for the Fourier space.An instability occurs when for the first time there are values of k at which the real part of λ is non-negative.In neocortex the extent of excitatory connections W aE is broader than that of inhibitory connections W aI , and so we take σ aI = 1 and σ aE = 2.In Fig. 5 we show a plot of the critical curves in the (v, γ) plane above which the homogeneous steady state, h ss E,I = h ss , is unstable to dynamic instabilities with |k c | = 0 (bulk oscillations) and |k c | = 0 (travelling waves).Figure 6 shows a pattern of parallel moving stripes seen beyond the Turing-Hopf  Despite the natural framework for neural field models being non-local integral equations with spacedependent delays the techniques for analysing them are nowhere near as developed as they are for local PDE models.Thus for this reason it is sometimes worth constructing equivalent PDE models.To date progress in this area has been made by Nunez [6] and Jirsa and Haken [19] for neural field models in one spatial dimension with axonal delays, links to the theory of damped inhomogeneous wave equations.For example in a one-dimensional setting with where Equation ( 22) is the Jirsa-Haken-Nunez brain wave equation [19,6,57].Turing instabilities in such models (built from two populations) have been exhaustively analysed in [58,59] and lead to bifurcation scenarios consistent with those in Fig. 5 dental equations (as it would for an ordinary differential equations with a fixed delay), though can be analysed making use of Lambert functions [60].
For a single self-coupled population with σ aa = σ, v aa = v and a simple exponential synaptic filter also supports travelling fronts [61].
Moreover, for steep sigmoids (β → ∞ and θ a = θ) the speed can be calculated in closed form as This strong dependence of the wave speed on the threshold θ has now been indirectly established in rat cortical slices (bathed in the GABA A blocker picrotoxin) [62].An applied positive (negative) electric field across the slice increased (decreased) the speed of wave propagation, consistent with (24) assuming that a positive (negative) electric field reduces (increases) the threshold θ.Of course travelling waves are also possible solutions of the more general non-local formalism described by (16), and have been explored mathematically in a number of papers, reviewed in [11], with a particular regard to studying waves in cortical slices [63], and notably epileptiform activity [64].
As well as supporting periodic travelling waves the brain wave equation ( 22) can support propagating pulses.A detailed mathematical analysis of waves of this type has been performed in [65] for the case of a Heaviside firing rate function (a sigmoid with infinite gain).The treatment of smooth firing rates can be pur-   sued for waves in one space-dimension using a standard travelling wave analysis (in which one looks for stationary profiles in a co-moving frame ξ = x − ct of speed c < v).This can often only be done numerically, say using the techniques of spatial dynamics and constructing homoclinic connections, as in [61].Interestingly propagating pulses of (22) can scatter in novel ways, leading to an elevated firing across the whole tissue [65].An example of this is shown in Fig. 7.
Finding equivalent brain wave equations in twospatial dimensions has proved far more challenging [18], with various approximations being made to obtain a local model.The most common of these is the so-called long-wavelength approximation, although this has recently been improved upon in [56].In the former case this amounts to expanding (21) around k = 0 for small k, yielding a "nice" rational polynomial structure to give (A 2 ab (ω) + 3k 2 /2)ψ ab (k, ω) = f b (k, ω) which may then be inverse transformed to give the telegraph PDE: This model has been intensively studied by a number of authors in the context of EEG modelling, see for example [66,67,68].
Undoubtedly the assumption of isotropic connectivity is a strong one for the modelling of cortical tissue.That it has been pursued so aggressively to date is more a reflection of the mathematical tractability of such models as opposed to their relation to real tissue.Indeed questions about the existence, uniqueness and absolute stability of solutions for truly inhomogeneous models are only just beginning to be addressed using tools from functional analysis [69].
However, one symmetry breaking effect that can be tackled without too much effort is that of the loss of continuous rotation symmetry.This is an important issue to treat in light of the fact that it is now known that (visual) cortex has a crystalline micro-structure at the millimeter length scale (reviewed in [70]).This has given rise to the notion of patchy connections that break continuous rotation symmetry (but not necessarily continuous translation symmetry).The natural way to model this is to introduce a modified connectivity kernel w P ab (r, r ) as where J ab (r) varies periodically with respect to a regular planar lattice L. Note that the patchy kernel w P ab is homogeneous, but not isotropic.The generalisation of brain wave equations to treat patchiness that arise via periodic modulation of an isotropic kernel has recently been developed by Robinson [71] and analysed further in [56].In essence the brain wave equation (25) is replaced by an infinite set of PDEs -indexed by the reciprocal lattice vectors q of the underlying lattice L, that arise in the Fourier series representation In the long-wavelength approximation this set of PDEs is obtained from (25) under the replacement ∇ → ∇ − iq, ψ ab → ψ q ab and u ab → η ab * q J q ab ψ q ab .Assuming that there is a natural cut-off in q, then we need only evolve a finite subset of these PDEs to see the effects of patchy connections on solution behaviour.For example with a square lattice we would need only to specify two reciprocal lattice vectors.
In this case a Turing instability analysis shows that, compared to the unmodulated case shown in Fig. 5, the Hopf bifurcation is transformed to a Turing-Hopf bifurcation with critical wavevectors coinciding with those of the lattice.With increasing v the dominant bifurcation is also of Turing-Hopf type.However, in this case it is a ring of wavevectors surrounding the reciprocal lattice vectors that go unstable first.In both cases this suggests the emergence of periodic travelling waves aligned to the lattice size and direction, which are indeed observed in direct numerical simulations [56].For a more general inhomogeneous kernel of the form w I ab (r, r ) = w ab (|r − r |)J ab (r ), then the PDE formulation goes over with the source term f b (r, t) in ( 25) replaced by J ab (r)f b (r, t).When the modulating kernel J ab (r) varies both weakly and rapidly in space then one may use techniques from homogenisation theory to study wave propagation and its failure [72].Indeed recent work avoiding the assumption of rapid spatial variation shows that for a one dimensional neural field model with J(x) = 1 + sin(2πx/σ) and f (u) = H(u − θ) then pulsating waves can occur (where say a front edge is modulated in space) with a wave speed c 0 1 − ( / c ) 2 so that propagation failure occurs for some > c = c (σ, θ) [73].Here c 0 is the speed of the wave in the homogeneous case (when = 0 and see equation ( 24) for example).Further treatment of heterogeneous connection topologies can be found in [74,75,76,77].
One of the more recent and compelling uses of neural field modelling is by Bojak and colleagues [13,14].
They relate different (co-registered) imaging modalities to one another, namely EEG (with its excellent temporal resolution) and fMRI (with its superior spatial resolution) by modelling an underlying neural generator that is based upon the Liley model [18] discussed earlier.Importantly regional connectivity data is incorporated using the CoCoMac database [78], which contains information on structural connectivity in the macaque brain (from tract-tracing experiments).In order to be compatible with true fMRI images their computational model uses triangular spatial grids matched to cortical geometries extracted from structural MR images.From the activity of the  [13,14].Associated animations may be found at http://www.mbfys.ru.nl/neuropi/cns09/.model cortical sheet, the simulated EEG at the scalp is found using a realistic volume conductor model.The response for fMRI BOLD is constructed using an established model for neurovascular coupling and convolving with "Balloon-Windkessel" haemodynamics [79].This is one of the first attempts to seriously combine neuronal dynamics and brain connectivity along the lines advocated in [80].An example of output from their modelling approach is shown in Fig. 8.For a recent review on the challenges in combining EEG and fMRI imaging modalities using neural models we refer the reader to Valdes-Sosa et al. [81], and for work on folded three-dimensional cortical sheets and their use in forward EEG and MEG prediction see Jirsa et al. [82].Other work on using realistic anatomical connectivities in large scale neural models (especially as it relates to the resting brain state) can be found in [83,84,85].
There are a number of natural extensions to neural field models to incorporate ever more biological realism, including the incorporation of slow intrinsic currents that underlie bursting behaviour in single neurons [86,11], synaptic depression and adaptation [87], dendritic processing [61], neuromodulation, say by endo-cannabinoids [88], and anaesthetic drug action [66,89,68,90].Moreover, it is possible to gener-alise such models even further to treat feature selectivity such as that observed in visual cortex for orientation [10], spatial frequency [91] and texture [92].
In this case neural activity is now also considered to be a function of some set of features χ and interactions must be specified by a feature dependent kernel w ab (r , t , χ |r, t, χ).Indeed with the increasing march of experimental progress in recording from populations of neurons we are now ideally poised to push forward with the development of multi-scale models that can integrate macroscopic models at large spatial scales with models at the microscopic scale [93].
One outstanding challenge for the development of tissue level firing rate models is how best to include gap junction coupling [94].

Discussion
A known limitation of neural field models is that they only try to track mean activity levels and cannot, by definition, track the higher order correlations of any underlying spiking model.Of course one approach would be to abandon them altogether in favour of more biophysically detailed models.However, as we have discussed here they can go a remarkably long way to providing a framework for interpreting neuro-imaging data whilst maintaining contact with known brain structure, dynamics and function.Rather it is preferable to work with spiking models and establish more concretely the link to mean activity models, and when more sophisticated kinetic models of brain activity are required [95,96].
Moreover, in some instances it is possible to analyse spiking networks directly (usually under the assumption of global coupling and fast synaptic interactions) as in the spike-density approach [97,98,99], which makes heavy use of the numerical solution of coupled partial differential-integral equations.In other situations equations going beyond the mean-field approach have been proposed that govern second-order correlations [100,101,102,103].Indeed there has been a recent upsurge of interest in this area adapting methods from non-equilibrium statistical physics to determine corrections to mean-field theory involving equations for two-point and higher-order cumulants [104,105].One immediate, yet potentially tractable, challenge would be to develop a framework for understanding networks of synaptically interacting nonlinear integrate-and-fire networks.At the single neuron level such models are already known to be able to capture many different physiological firing patterns and at the network level are computationally far cheaper to implement than their conductance based cousins yet still able to generate the rich repertoire of behaviour seen in a real nervous system [106].One step in this direction has already been undertaken for the absolute-integrate-and-fire model [107] and this may well be a natural spiking model to pursue for the development of a specific soluble spiking neurodynamics.

Focusing
on the Liley et al. model [18], we review the success of such descriptions in generating oscillations consistent with the alpha band of the human EEG spectrum.Extensions to this population model to treat refractoriness and spike-frequency adaptation are also discussed.Next we show to model a cortical area as a two-dimensional continuous network of such cortical columns, defining a neural field.After reviewing the basic instability mechanism that can lead to the formation of travelling patterns of activ-ity we show how to formulate the model in terms of a PDE, recovering the Jirsa-Haken-Nunez brain wave equation in one spatial dimension

Figure 1 :
Figure 1: A diagram of a local cortical module represented by the interaction of two neuronal populations, one excitatory (E) and the other inhibitory (I).
and the latter, P ab , is considered constant.Note the inclusion of delays ∆ ba in (2) that represent the fixed axonal communication lag for action potentials propagating from population b to a. Exploiting the linearity of Q ab the model for the conductance can be integrated to give g ab = η ab * [f b + P ab ], where * denotes a temporal convolution:

Figure 2 :
Figure 2: Bifurcation diagram for the Liley population model showing the absolute maximum of E in terms of P EE for steady state (red), small amplitude periodic (blue) and large amplitude periodic (green) orbits.

Figure 3 :
Figure 3: An example of multistability in the Liley population model.A small and large amplitude stable rhythm co-exist over a range of parameter values.
This transcendental equation allows for the possibility of complex roots, and not surprisingly it is possible to choose values for the slope and threshold of the sigmoid such that

R 2 dr
w ab (r, r )f b (r , t − |r − r |/v ab ) (18) and f b is the firing rate in layer b.Here r ∈ R 2 and w ab (r, r ) prescribes the coupling strength between position r in layer a and position r in layer b.The velocity of an action potential travelling along a fibre connecting layer b to layer a is denoted v ab and underlies the space-dependent delay |r − r |/v ab for signals propagating over a distance |r − r |.Note that (18) provides meaning for the operator on the right hand side of the expression ψ = w ⊗ f .To close the system of equations one could choose the firing rate to depend on some dynamic mean-membrane potential as in the Liley model.Alternatively, to recover a purely activity based model in the spirit of Wilson-Cowan and Amari one could set f a = f a (h a ), where h a = b u ab .For simplicity we shall restrict further discussion to this case.For a recent historical perspective of the Wilson-

Figure 5 :
Figure 5: Critical curves showing the instability borders for dynamic instabilities in the (v, γ) plane, where γ = f (h ss ).

Figure 6 :
Figure 6: Snapshots of a periodic travelling wave in u EE , each 1/4 of a period later than the previous one (ordered top left, top right, bottom left, bottom tight), for the 2D model of Fig. 5 at v = 12 and γ = 15 (on a 30 × 30 domain).
, when considering analogous architectures.Moreover, a weakly non-linear analysis of the travelling and standing waves that develop beyond the point of instability has been developed.The appropriate amplitude equations are found to be the coupled mean-field Ginzburg-Landau equations describing a Turing-Hopf bifurcation with modulation group velocity of O(1).In particular this has allowed an investigation of Benjamin-Feir modulational instabilities in which a periodic travelling wave (of moderate amplitude) loses energy to a small perturbation of other waves with nearly the same frequency and direction.For asymmetric kernels (w(x) = w(−x)) or symmetric kernels with a peak away from the origin then the stability analysis of the homogeneous steady state involves the solution of transcen-

Figure 7 :
Figure 7: Numerical simulations of a coupled neural field model in one spatial dimension showing the interaction of two travelling pulses and the emergence of transient complex behaviour ultimately leading to an elevated firing rate across the whole tissue.The model is defined by [A 2 a − ∂ xx ]ψ a = Γ a A a H(u e − u i − θ)/σ a , A a = (σ −1 a + v −1 a ∂ t ), (1 + α −1 a ∂ t )u a = ψ a, with a ∈ {e, i}, describing a model with short range inhi- with a ∈ {e, i}, describing a model with short range inhibition and long range excitation.The plot shows the evolution of u e with v e = 0.2, v i = 1, α e = α i = 1, Γ e = 1, Γ i = 0.7, σ e = 2, σ i = 1 and θ = 0.07.

Figure 8 :
Figure 8: An EEG scalp isopotential (left) generated from the mean membrane excitatory potential of a neural field model (middle) and the corresponding fMRI BOLD signal construction on the cortical surface (right), illustrating the capabilities of the recent computational framework developed by Bojak and colleagues[13,14].Associated animations may be found