Standing and travelling waves in a spherical brain model: The Nunez model revisited.

The Nunez model for the generation of electroencephalogram (EEG) signals is naturally described as a neural field model on a sphere with space-dependent delays. For simplicity, dynamical realisations of this model either as a damped wave equation or an integro-differential equation, have typically been studied in idealised one dimensional or planar settings. Here we revisit the original Nunez model to specifically address the role of spherical topology on spatio-temporal pattern generation. We do this using a mixture of Turing instability analysis, symmetric bifurcation theory, centre manifold reduction and direct simulations with a bespoke numerical scheme. In particular we examine standing and travelling wave solutions using normal form computation of primary and secondary bifurcations from a steady state. Interestingly, we observe spatio-temporal patterns which have counterparts seen in the EEG patterns of both epileptic and schizophrenic brain conditions.

Introduction

Modern neuroimaging methodologies give us a window on the activity of the brain that may reveal both structure and function. Despite the recent advances in technologies for magnetic resonance imaging (MRI) for assessing anatomy, and functional MRI for assessing functional changes over seconds or minutes, the historical predecessor of electroencephalography (EEG), and its more recent magnetic counterpart magnetoencephalography (MEG), is still a hugely practical non-invasive tool for studying brain activity on the millisecond time-scale. The compromise being the relatively poor centimetre scale spatial resolution. However, from a modelling perspective this can actually be beneficial as it favours a more coarse grained description of neural tissue without recourse to detailed neuronal modelling. Indeed this is the current mode of thinking in cognitive neuroscience studies where single scalp electrodes (used in arrays across the head) are used to record the activity of neurons. Models that capture the large scale dynamics of neural tissue are often referred to as neural field models, and see  [1] for a recent discussion.

The model is, in its most general setting, described as a dynamical system with space-dependent delays that invariably is thought of as an integro-differential equation—which reduces to a damped inhomogeneous wave equation for a particular choice of exponentially decaying spatial interactions. Perhaps, Paul Nunez was one of the first to realise the importance of modelling the long range cortico–cortico connections for generating the all important -rhythm of EEG (an 8–13 Hz frequency)  [2]. Moreover, he recognised that because the cortical white matter system is topologically close to a sphere that a model that respected this (with periodic boundaries) should naturally produce standing waves (via interference)  [3], [4]. For a more recent perspective on this work see  [5].

Given the importance assigned by Nunez to the boundary conditions  [6], the model has, surprisingly, been studied more often than not in scenarios that have different topologies to that of the sphere, e.g.  [7], [8], [9], [10]. Understandably this facilitates both mathematical and numerical analyses, though the results have less relevance to the application of standing waves seen in EEG. One exception to this is the numerical study of Jirsa et al.  [11], though even here analysis and simulations are performed by using the less general partial differential equation (PDE) formulation of the model. Despite the significance of geometry in the Nunez model for understanding EEG, a thorough exploration of its pattern forming properties has not been performed since its inception roughly forty years ago. In this paper we undertake a first step along this path.

Neural fields and symmetry

Analysis of spatio-temporal patterns in dynamical systems goes hand in hand with identifying the various symmetries in the model. Both the internal structure of the model, lattice-structure for example, and the domain under consideration, e.g. a disk, will impact the system’s symmetries. Since neural fields are primarily studied on either infinite or periodic domains, the group of translations and rotations (Euclidean group) arises naturally in this setting. Ermentrout and Cowan used this to good effect in developing their original model for visual hallucinations in primary visual cortex (V1), arising from a Turing instability  [12]. Apart from describing the time-evolution of activity in a strictly anatomical space, neural fields have been extended to ‘feature spaces’, which allow one to represent abstract attributes of neural activity. An outstanding example of this are the models of V1 of Cowan and collaborators, who included the cortical columns’ orientation preference for visual input in the framework of neural fields. A detailed analysis of the shift-twist symmetry group, which is at the heart of these models, has yielded an understanding on the origin of visual hallucinations characterised as lattices of locally oriented contours or edges  [13], [14]. Extending the model to account for spatial frequency preference of the visual stimuli, a feature distinct from the orientation, has also resulted in the formulation of a neural field on a sphere  [15]. Yet the differences with Nunez’ interpretation are marked: Nunez focuses on direct anatomical connectivity rather than interactions in an abstract feature space. In this work, we follow Nunez and take the sphere as the physical domain of the neural field.

Spherical symmetries have a long history of application, e.g. they play a role in morphogenesis (how an initial spherical ball of cells develops into a mature shape)  [16], as well many other biological and physical systems including the understanding of tumour growth  [17], sphere buckling under pressure  [18], and Rayleigh–Bénard convection in a spherical shell  [19] to give but a few examples. Typical models for these systems take the form of PDEs, such as reaction–diffusion or Swift–Hohenberg, and the techniques for understanding bifurcations from spherically symmetric states have included group theory  [20], scientific computation  [21] and Turing instability analysis  [22]. For a further discussion we refer the reader to the article by Matthews  [23].

Role of time delays

In contrast to other studies of pattern formation on a sphere we are concerned not with PDEs, but rather with non-local models. The Nunez model can be parsimoniously expressed as an integro-differential equation with time delays. Although the very first formulations of neural fields already include transmission delays  [24], [25], they have often been disregarded in rigorous analysis, due to lack of a proper mathematical setting for these problems. Recently, Faugeras and collaborators made their first steps in formulating a rigorous framework for these models  [26], [27], [28]. Subsequently, van Gils et al. proved that this class of dynamical systems can be cast as abstract delay differential equations  [29]. As a result of this, many of the mathematical techniques developed for the analysis of ODEs and PDEs, such as Turing analysis, symmetric bifurcation theory, and centre manifold reduction can be adapted for use in the delayed integro-differential equation we study in this article.

In  [30] a delayed neural field is studied on a one-dimensional interval and symmetries are used to simplify the computations of spectral values and normal form coefficients for a pitchfork–Hopf bifurcation. The inhomogeneities in their model (i.e. the boundaries) complicate the analytical computation of eigenfunctions and critical normal form coefficients—numerics have to be used instead. As a contrast, we will employ a Turing analysis in this paper, which enables us to express the eigenfunctions in closed form using spherical harmonics. Consequently, we are able to identify closed expressions for the critical normal form coefficients, where numerics are only required for computing the eigenvalues as a solution of a transcendental equation (which is common practice for delayed systems).

Although we are able to perform the centre manifold reduction with minimal numerical effort, forward-time simulations of the model are a whole different challenge. Indeed, the toolbox of numerical schemes is as yet relatively underdeveloped and so here we apply a novel scheme for the simulation of (discretised) integro-differential equations on large meshes  [31]: linear features of Cubic-Hermite spline interpolation and numerical integration are exploited to express the majority of operations in sparse matrix–vector products.

Outline

In Section  2 we give a brief review of the relevant neocortical anatomy and physiology to set the scene for the mathematical formulation of the large-scale Nunez model of EEG. We consider the case that the anatomical connectivity function is invariant with respect to the symmetry group of the sphere and show this can naturally be constructed using a spherical harmonic basis set. The non-instantaneous interaction between cortical regions is described with the use of a space-dependent delay determined by the speed of an action potential along an axonal fibre. Next in Section  3 we perform a linear Turing analysis of the steady state to show that both spatial and spatio-temporal neural activity patterns can occur (as linear combinations of spherical harmonics), depending on the precise shape of the connectivity function and the speed of the action potential. The techniques from equivariant bifurcation theory, which enable us to identify the possible planforms that can arise at the bifurcation, are reviewed in Section  4. Similarly, in Section  5 we offer a comprehensive overview of the framework of sun–star calculus that is required to perform the centre manifold reduction and critical normal form computation in neural fields with time delays. We apply these techniques to our model in Section  6 to obtain explicit expressions for the critical normal form coefficients, which enable us to classify the system’s bifurcation. In particular, we determine the first Lyapunov coefficient of a Hopf bifurcation and, by continuation, we subsequently find two different codimension two bifurcations: a generalised Hopf bifurcation and a double Hopf bifurcation. Both bifurcations give rise to bistability in the system, which we investigate analytically as well as numerically. Finally in Section  7 we discuss natural extensions of the work in this paper.

A model of cortex with axonal delays

The columnar organisation of the neocortex has been appreciated for some time, and for a review see  [32]. These (internally connected) macrocolumns consist of neurons with similar response properties and tend to be vertically aligned into columnar arrangements of roughly 1–3 mm in diameter. Columns in cortical areas located far from one another, but with some common properties, may be linked by long-range, intracortical connections (1–15 cm). Thus, to a first approximation the cortex is often viewed as being built from a dense reciprocally interconnected network of corticocortical axonal pathways, of which there are typically 1010 in a human brain  [33]. These fibres make connections within the roughly 3 mm outer layer of the cerebrum, and this wrinkled and folded cortical structure contains about 1010 neurons. Approximately 80% of these connections are excitatory and the remainder inhibitory. Excitatory pyramidal cells generally send their myelinated axons to other parts of the cortex (forming the white matter), so that most long-range synaptic interactions are excitatory. Roughly 95% of these connections target the same cerebral hemisphere, whilst the remaining ones either cross the corpus callosum to the other hemisphere or connect to the thalamus. In contrast inhibitory interactions tend to be much more short-ranged. It is the combination of local synaptic activity and non-local interactions within the cortex that is believed to be the major source of large-scale EEG and MEG signals recorded at (or near) the scalp.

Perhaps the most definitive model of EEG generation to date is that of Paul Nunez (reviewed in  [4]), which has culminated in the brain-wave equation for EEG generation. Indeed this and more general neural field models (reviewed in  [1]) are the major frameworks for the forward generation of EEG signals. At heart these modern biophysical theories assert that EEG signals from a single scalp electrode arise from the coordinated activity of pyramidal cells in cortex  [34]. EEG resolution (from the scalp) is typically in the 6 cm range for unprocessed EEG and 2–3 cm for high resolution EEG  [35]. Thus the number of neurons contributing to each scalp electrode is expected to be roughly 109 for unprocessed EEG and 108 for high resolution EEG. These are arranged with their dendrites in parallel and perpendicular to the cortical surface. When synchronously activated by synapses at the distal dendrites extracellular current flows (parallel to the dendrites), with a net membrane current at the synapse. For excitatory (inhibitory) synapses this creates a sink (source) with a negative (positive) extracellular potential. Because there is no accumulation of charge in the tissue this distal synaptic current is compensated by other currents flowing in the medium causing a distributed source in the case of a sink and vice versa for a synapse that acts as a source. Hence, at the population level the potential field generated by a synchronously activated population of cortical pyramidal cells behaves like that of a dipole layer. The interneurons’ contribution to the electrical field, on the other hand, is negligible due to both the small cell volume and the lack of a clear dipolar (or other orientation-dominant) morphology.

Nunez has convincingly argued that the dynamics of neural membrane alone cannot credibly account for the robust human EEG rhythms seen in the 1–15 Hz range, primarily because there is no such thing as a fixed membrane time constant in vivo (since for voltage gated membrane ion channels this is a time and state dependent attribute). However, local delays arising from synaptic processing (seen in the rise and decay of post synaptic potentials) as well as global delays arising from action potential propagation along corticocortical fibres are believed to be far more important. The former typically have time-scales from 1 to 100 ms and the latter of up to 30 ms in humans. The Nunez model of EEG respects the physiology and anatomy described above and has been particularly successful for describing standing EEG waves. Indeed these motivate one form of the model as a damped inhomogeneous wave equation whereby standing waves arise naturally by interference in a system with periodic boundary conditions. Nunez considered each cortical hemisphere (together with its white matter connections) to be topologically equivalent to a sphere, with the speed of an action potential fixed for all fibres—thus ignoring known anisotropy in the form of a preferred anterior–posterior orientation  [3], [4]. The radius of each cortical hemisphere was calculated from the known surface area of as . Taking a value for the corticocortical action potential speed in the range Nunez used simple interference arguments (using an analogy with vibrations on a string) to predict that fundamental cortical frequencies (for standing waves) would lie in the range , using the relationship . Another version of the model takes the form of an integro-differential equation and it is this formulation of the model that we shall consider in this paper.

The model

Here we give a modern perspective on the Nunez model in its integral form. Using some artistic licence we also slightly generalise it to include a simple model of synaptic processing, to bring it more in line with popular formulations of neural field theories. For a review of these we refer the reader to the recent book  [36].

We represent synaptic activity by where is a point on the surface of a sphere and . We shall consider simple neural field models that, after rescaling time and space, can be written in the form with either or We also have to specify the initial condition: where , the state space defined below.

Here is the surface of the unit sphere in , the integration measure and denotes function composition, with a firing rate function. The weight distribution specifies the anatomical connectivity between points and , whilst specifies the corresponding delay arising from the finite speed of the action potential travelling along the fibre connecting the two points. The model defined by (2a) is often referred to as a voltage based model, whereas (2b) is referred to as an activity based model  [37]. In either case the models are qualitatively similar in their behaviour, and the analytical and numerical techniques we develop throughout this paper can be adapted to either case. For concreteness we shall work with (2a).

Functional analytic setting

Throughout this paper we fix the following assumptions: These assumptions are nearly identical to those in the seminal work  [29], apart from one: where the original work sets the spatial space as the continuous functions, we have chosen the Hölder continuous functions instead. This additional constraint on the spatial domain is in our case required to guarantee the convergence of the spherical harmonics expansion (see next section). Moreover, we point out this change of spatial function space has no apparent impact on the outcomes of the original work and a discussion of these minute adjustments is therefore omitted.

the firing rate function is smooth and bounded on ,

the domain is the unit sphere in and the corresponding metric is the great circle distance,

the connectivity kernel , the Banach space of Hölder continuous functions with exponent , on .

the delay function is non-negative and not identically zero,

the maximal delay ,

the spatial space the Banach space of Hölder continuous functions with exponent , , with the standard norm:

the state space the Banach space of Lipschitz continuous functions with the standard norm:

a function that satisfies (1) with initial condition is a global solution.

We include for completeness the proofs of two lemmas that guarantee that the equations (1)  +  (2a)  +  (2c) are well-posed. Using the notations of  [29], define the nonlinear operator as We have the following lemma.

is well-defined by   (3).

Given , we consider two points and of and write  □

Using the definition of the operator , the system (1)  +  (2a)  +  (2c) can be rewritten as the following initial value problem Then (4) is of the form of a Delayed Differential Equation when we define by It is then well-known that (4) has a unique solution or equivalently that (1)  +  (2a)  +  (2c) has a unique global solution if the operator defined by (5) is Lipschitz continuous. We prove the following lemma.

The operator defined by(5)   is Lipschitz continuous.

If and are elements of and , consider for some positive constant . Consider now We also have so that for some positive constant . Combining (6), (7) we obtain for some positive constant , i.e.  is Lipschitz continuous. □

Spherical geometry

For the remainder of this paper we will assume that is a point on the sphere with polar angle , azimuthal angle and radius 1 and similar for . Furthermore, we make use of the complex-valued spherical harmonics of degree and order , for which a representation is given in Appendix A. As much as Fourier series form an orthonormal basis on the circle, spherical harmonics form an orthonormal basis on the sphere.

Spherical Harmonics Expansion  [38, Thm. 5]

Let , then the spherical harmonics expansion of converges uniformly to , that is for The coefficients are given by projections on the basis functions:where the overline denotes complex conjugation.

We give a very short proof for completeness. Let be the linear operator defined by and be the infinite norm of the difference : Ref.  [39, Cor. 2.2] implies that , while Ref.  [38, Thm. 4] implies that the operator norm of the linear operator w.r.t. the uniform norm on satisfies . Finally  [38, Thm. 1] shows that and this yields hence the uniform convergence on of to when . □

We shall consider a homogeneous neural field, where both the weight kernel and transmission delays depend on the relative position of and . Naturally, and are chosen as functions of distance along the surface, but more generally we set We remark that, on a unit sphere, the inner product is equal to the cosine of the angular separation (and therefore great circle distance) between two points; i.e.  . This leads to the following theorem.

Let be a Hölder continuous function on with exponent , , and . Then the series with coefficientswhere is the Legendre polynomial of degree , converges to , uniformly on , when .

It is easily checked that the hypothesis on the function implies that is Hölder continuous on with exponent : An easy extension of Theorem 1 shows that the series converges uniformly to when , where the coefficients are given by the projections: Since is continuous on hence in (integrable), we can apply the Funk–Hecke theorem  [40, Vol. 2, p. 247] to (12) and obtain due to orthogonality. □

We note that a similar mathematical representation has previously been used in  [14], [15] for a neural field model describing orientation and spatial frequency tuning in a cortical hypercolumn. However, the physical differences between our studies are marked, with ours focusing on direct anatomical connectivity rather than interactions in a neural feature space.

Concrete choices

While the majority of the following results are independent of the specific choices for , and , we will illustrate our results with concrete choices in computations and simulations. Where required, we choose the following forms.

The firing rate function is sigmoidal and increasing: with steepness parameter and threshold .

A natural choice for the connectivity kernel is with and . For we have a wizard-hat shape, whilst for we have an inverted wizard-hat shape.

We have the following lemma.

The function defined by   (16)   is Hölder continuous with exponent , .

The space of Hölder continuous functions with the same exponent being a vector space, it is sufficient to check that is Hölder continuous with exponent . We have Since , we have Since is Lipschitz continuous, , for some positive constant , and since we have , . It follows that for some positive constant and for all . □

Furthermore, we call the synaptic kernel balanced if Other choices than (16), such as a difference of Gaussians, are also natural.

A common choice for the transmission delay is which incorporates both a constant offset delay and a constant propagation speed of action potentials. An onset delay has been shown to lead to dynamics reminiscent of those seen in simulations of large-scale spiking networks  [41], and its physiological interpretation can be connected to the relaxation time-scale for which spiking networks can reasonably allow for a firing rate description. We have the following lemma.

The function defined by   (17)   is Hölder continuous with exponent , .

It follows from the proof of Lemma 3 that is Hölder continuous with exponent , . □

Linear stability analysis

It is sensible to begin the investigation of the spherical Nunez model with a standard Turing analysis, treating the instability of a homogeneous steady state. This will allow us to determine the conditions for the onset of spatial patterned states or more dynamic spatio-temporal patterns in the form of standing or travelling waves. This approach has a long tradition in neural field modelling, as exemplified in  [12], [42], [13] and reviewed in  [37], [43]. A one-dimensional system with space-dependent delays and periodic domain is studied in  [28] using the same methodology.

Spectral problem

Stability and pattern formation in the model are dictated by the spectral values, or eigenvalues, of the semigroup generator underlying the dynamical system. We continue with a heuristic derivation of a set of characteristic equations whose roots are the eigenvalues of (1)  +  (2a). A detailed treatise on the validity of our result, as well as special properties of the spectral problem, is available in  [29].

A homogeneous steady state of (1)  +  (2a) satisfies For our choice (15), up to three homogeneous equilibria might be present. Note that in the special case that the kernel is balanced, i.e.  , is the only homogeneous equilibrium.

Linearising (1)  +  (2a) around gives where . Solutions of this linear equation are separable and we set . In this case satisfies the linear equation where is the characteristic function with . Note that the particular structure of allows application of Theorem 2, which yields coefficients —see Appendix B. Indeed, Lemma 3, Lemma 4 show that the functions and , are Hölder continuous with exponent , . It is then easily verified that the function , is also Hölder continuous with exponent , and the beginning of the proof of Theorem 2 shows that the function , is Hölder continuous with exponent . Now, is an eigenvalue of (19) if has non-trivial solutions. This occurs when for some . In this case, (19) has solutions (i.e. eigenfunctions) of the form , . Hence, the algebraic and geometric multiplicity of the eigenvalues are the same.

The spectrum corresponding to (19) consists of both a point spectrum and a (Browder) essential spectrum , see also  [29]. The point spectrum consists of all complex numbers which solve (21), We call a regular value.

Resolvent

Recall that the point spectrum corresponds to values of for which the operator is not invertible. For all regular values the following theorem holds.

Resolvent

For given and , there exists a unique which solvesIf both and are expanded in spherical harmonics, see   Theorem  1, then coefficients are given by

Existence and uniqueness of the solution are shown in  [29, Prop. 14]. To identify the solution, we start with , substitute (20) and expand using Theorem 2 where integration and summation are interchanged. Next, we multiply both sides by , integrate over the domain w.r.t. and change the order of integration: where only one term remains in the summation due to orthonormality of the spherical harmonics. Application of Theorem 1 yields Since , the first factor is non-zero and, hence, we obtain the final identity via division. □

Although we will not use the resolvent in the remainder of this section, it will play a prominent role in the evaluation of the normal form coefficients in Section  5 and Appendix C.

Stability region

The homogeneous steady state is stable if for all . If real eigenvalues vanish, due to a parameter variation, a fold or transcritical bifurcation will occur. If, instead, eigenvalues have a vanishing real part but a nonzero complex part, then a Hopf bifurcation is expected. In the former case one would expect the formation of time-independent patterns, and in the latter the emergence of travelling or standing waves. Note that in the absence of delays (), all eigenvalues are real and given explicitly by . In this case Hopf bifurcations are not possible. The focus of the remainder of this paper will be on the emergence of spatio-temporal patterns as expected from the general theory of Hopf bifurcations with symmetry  [44], [45].

For the chosen connectivity and delay functions (16), (17) we are able to find explicit expressions for the coefficients , of the form with according to Appendix B. Using (21), one can identify the parameters at which the homogeneous steady state undergoes a fold or transcritical bifurcation—in this case . Indeed, defines a line in the -plane which corresponds to candidate fold/transcritical bifurcations corresponding to the mode number . In Fig. 1, these lines are plotted with different colours according to their mode number. Similarly, we trace out the candidate Hopf bifurcations by solving (21) for , whose solutions are expressed parametrically in terms of by equating real and imaginary parts: These parametric curves are plotted for increasing in Fig. 1 as solid lines.

Fig. 1

Stability and bifurcations of the homogeneous steady state in the absence of the offset delay . Parametric curves in the -plane mark the boundary of the stability region, which is coloured grey. Solid (dashed) coloured lines represent parameters at which the steady state undergoes a Hopf (fold/transcritical) bifurcation with respect to the spherical harmonics of degree . Grey parallel lines in the background mark lines along which is constant. In particular, the line passing through the origin is dashed, which corresponds with a balanced kernel for which explicit calculations can be made. Parameters: , , , . (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

Furthermore, for small and , the linear equation (19) is dominated by the term such that solutions near the origin of the -plane are asymptotically stable. Since instabilities only occur at the curves (23), (24), the stability region can be extended from the origin to the first bifurcation. This region is coloured grey in Fig. 1.

Recall that the fixed points of the system (18) depend on and, hence, has an implicit dependency on and as well. As such, the coordinates shown in Fig. 1 are conditional on solutions of the transcendental equation (18). For parameter variations, however, along a line in the set , the fixed point structure–and therefore –remains unaltered. This collection of parallel lines is illustrated in Fig. 1 for various . One line is of particular importance: the line corresponds to a balanced kernel, such that is the only (homogeneous) fixed point in the system and can be determined explicitly.

Remarks

From Fig. 1 we conclude that predominantly spatially homogeneous instabilities are expected to occur, since the stability region is largely bounded by fold/transcritical and Hopf bifurcations corresponding to . Only a small part of the stability region, as shown in the inset, is bounded by curves relating to bifurcations of higher degree in .

Fig. 2 depicts a similar diagram to Fig. 1 for different parameters values. In particular, the offset delay is non-zero. First of all, it is noted that this increment results in a shift of the stability region: the stability region is now positioned more symmetrically around the origin. Furthermore, the Hopf bifurcation curves have a richer structure than in the case where , giving rise to many intersections (corresponding to double Hopf bifurcations). As a consequence, the stability region is bounded by curves relating to Hopf bifurcations of spherical harmonics of degrees . The close-up shows the part of the stability region which is bounded by the Hopf curve for degree .

Fig. 2

Stability and bifurcations of the homogeneous steady state for non-zero . Similar to Fig. 1, but for different parameters. The parallel lines along which is constant are not plotted for clarity. The inset shows a marker at parameter values for which the spectrum is shown in Fig. 3. Parameters: , , , . (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

For parameters indicated with the marker in the inset of Fig. 2, we compute the spectrum to verify the foregoing analysis. The result is depicted in Fig. 3 where we show eigenvalues in the plane for , as determined by solving (20) numerically. A pair of complex eigenvalues, corresponding to can be seen in the positive half-plane, which matches with leaving the stability region by crossing the Hopf curve of the same degree (magenta in Fig. 2).

Fig. 3

Spectrum of the spatially homogeneous steady state after a Hopf bifurcation. The eigenvalues are determined by solving (20) numerically for . The pair in the right half-plane corresponds with . Parameters as in Fig. 2 and additionally and .

In Fig. 4 we show a direct simulation of an instability to a pattern state with and , highlighting the emergence of a standing wave. Predictions of the bifurcation point as predicted by the linear stability analysis are found to be in excellent agreement with the results from direct numerical simulations in all cases. These were performed using a bespoke numerical scheme, which combines standard techniques for tessellating the sphere with a new approach for solving integro-differential equations with delays on large meshes  [31].

Fig. 4

A direct simulation of the spherical Nunez model just beyond the point of an instability with showing the onset of a standing wave. Left to right, top to bottom shows eight (equally spaced in time) snapshots of the standing wave for one period of oscillation. Warm (cold) colours correspond to high (low) values of . Parameters as in Fig. 3 and additionally and . Simulated with a mesh of 5120 triangles. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

Intermezzo: planforms

Near a Hopf bifurcation, as identified by the foregoing analysis, the solutions destabilise tangent to the critical eigenspace. Therefore, we expect a dynamical pattern of the form where stands for complex conjugate, with and determined from the spectral equation (21) such that while for all other . Sufficiently close to the bifurcation point we expect certain classes of solutions to emerge that break the symmetry of the homogeneous steady state. The tools of equivariant bifurcation theory help to identify solution candidates based on symmetry arguments alone  [44, Chapter XVIII Section 5],  [46]. Typically this is done by developing a system of ordinary differential equations for the evolution of the amplitudes . If we denote the space spanned by the spherical harmonics of degree by then the action of on is determined by its action on and we will consider , where is a smooth function that commutes with the action of on . By this we mean that is equivariant with respect to the action of on :

In order to determine the structure of the equivariant normal form on to a given order we use the action of elements of on the spherical harmonics of order and the raising and lowering operators  [47]. Recall that and note that the inversion element acts on spherical harmonics of degree as multiplication by (for the natural action, which is the action which we shall be considering throughout). The group contains a maximal torus corresponding to rotations in the direction (-axis) through an arbitrary angle which act on as multiplication by . The raising and lowering operators act on the amplitudes as  [47]

Symmetry can also be used to determine branches of periodic solutions of . Let, without loss of generality, be a periodic solution of with period . Then is a spatiotemporal symmetry of if The set of all spatiotemporal symmetries of a solution is a subgroup of called the isotropy subgroup of and denoted . An isotropy subgroup is -axial if (i.e. the subspace of which is invariant under the action of is two dimensional). The equivariant Hopf theorem (see  [44]) states that is guaranteed to have a branch of periodic solutions with symmetry (isotropy) bifurcating from the point of dynamic instability if is a -axial isotropy subgroup of for the action of on . This theorem also requires that the eigenvalues cross the imaginary axis with non-zero speed. Thus the symmetries of the branches of periodic solutions bifurcating at a dynamic instability where the modes of degree become unstable correspond to the subgroups which are -axial under the action of on . Which subgroups are -axial isotropy subgroups depends on the value of and have been determined for all values of   [44], [45]. The isotropy subgroups of are twisted subgroups where is a subgroup of and is a group homomorphism. Elements of can be thought of as spatial symmetries whilst elements of are temporal symmetries acting on periodic solutions by a phase shift. An element is a spatial symmetry if and a spatiotemporal symmetry if . The spatial symmetries form a normal subgroup of and the quotient group is isomorphic to a closed subgroup of (i.e.  , or ). The -axial isotropy for values of between 1 and 6 is given in Table 1 which is reproduced from  [45]. All notations for the subgroups of here and throughout are consistent with that in  [44], [45].

Table 1

The -axial subgroups of for the natural representations on for . Here .

n_c	J	K	θ(H)	Number of branches given by equivariant Hopf theorem
1	O(2)	O(2)^−	Z_2	2
	SO(2)	Z2n−	S^1[n=1]
2	O(2)	O(2)×Z2c	1	5
	SO(2)	Z_n×Z2c	S^1[n=1,2]
	T	D_2×Z2c	Z_3
	D_4	D_2×Z2c	Z_2
3	O(2)	O(2)^−	Z_2	6
	SO(2)	Z2n−	S^1[1≤n≤3]
	O	O^−	Z_2
	D_6	D6d	Z_2
4	O(2)	O(2)×Z2c	1	10
	SO(2)	Z_n×Z2c	S^1[1≤n≤4]
	O	O×Z2c	1
	T	D_2×Z2c	Z_3
	D_8	D_4×Z2c	Z_2
	D_6	D_3×Z2c	Z_2
	D_4	D_2×Z2c	Z_2
5	O(2)	O(2)^−	Z_2	11
	SO(2)	Z2n−	S^1[1≤n≤5]
	T	D_2	Z_6
	D_10	D10d	Z_2
	D_8	D8d	Z_2
	D_6	D6d	Z_2
	D_4	D4d	Z_2
6	O(2)	O(2)×Z2c	1	15
	SO(2)	Z_n×Z2c	S^1[1≤n≤6]
	I	I×Z2c	1
	O	O×Z2c	1
	O	T×Z2c	Z_2
	T	D_2×Z2c	Z_3
	D_12	D_6×Z2c	Z_2
	D_10	D_5×Z2c	Z_2
	D_8	D_4×Z2c	Z_2
	D_6	D_3×Z2c	Z_2

Observe for example that if the spherical harmonics of degree become unstable at the dynamic instability, then branches of periodic solutions with at least ten distinct symmetry types bifurcate. These solutions are listed in Table 2 and illustrated in Fig. 5. There are six standing wave solutions (where the image of subgroup under the homomorphism is 1 or ) and four travelling wave solutions (which have since the spatial symmetry of a travelling wave does not change over time). Note from Table 2 that the travelling wave solutions consist of a single spherical harmonic rotating in one direction whereas standing wave solutions result from the sum of spherical harmonics and . The resulting standing wave can be thought of as arising due to interference between the waves travelling in opposite directions around the sphere.

Table 2

The -axial subgroups of for the natural representations on . Here .

Σ	J	K	θ(H)	Fix(Σ)
O(2)˜	O(2)	O(2)×Z2c	1	{(0,0,0,0,z,0,0,0,0)}
O˜	O	O×Z2c	1	{(5z,0,0,0,14z,0,0,0,5z)}
T˜	T	D_2×Z2c	Z_3	{(7z,0,12iz,0,−10z,0,12iz,0,7z)}
D_8˜	D_8	D_4×Z2c	Z_2	{(z,0,0,0,0,0,0,0,z)}
D_6˜	D_6	D_3×Z2c	Z_2	{(0,z,0,0,0,0,0,z,0)}
D_4˜	D_4	D_2×Z2c	Z_2	{(0,0,z,0,0,0,z,0,0)}
SO(2)_4˜	SO(2)	Z_4×Z2c	S^1	{(z,0,0,0,0,0,0,0,0)}
SO(2)_3˜	SO(2)	Z_3×Z2c	S^1	{(0,z,0,0,0,0,0,0,0)}
SO(2)_2˜	SO(2)	Z_2×Z2c	S^1	{(0,0,z,0,0,0,0,0,0)}
SO(2)_1˜	SO(2)	Z2c	S^1	{(0,0,0,z,0,0,0,0,0)}

Fig. 5

The ten periodic solutions guaranteed to exist at a dynamic instability with corresponding to axial isotropy subgroups as listed in Table 2. (i)–(vi) illustrate the evolution of the six standing wave solutions over one period and (vii)–(x) illustrate the travelling wave solutions indicating the apparent axis and direction of rotation.

Intermezzo: centre manifold reduction

In order to determine which patterns are stable near the bifurcation, we apply the method developed in  [29], which provides a generic method for normal form computation in delayed neural fields. The functional analytic setting of this work is based on sun–star calculus, which is described in  [48] for traditional delay differential equations. Since the forementioned works are particularly technical, we aim to offer the reader a rudimentary overview of the sun–star framework leaving out many details; these can be found in the original works. In particular, we will introduce and discuss all components of Eqs. (35a)  +  (35b), which are required to compute the critical normal form coefficients of bifurcations.

Sun-star calculus

The space is the state space of the delayed equation (1) and elements relate to the system’s history via . Next, consider the linear system with initial condition and linear operator . For our particular system, is given by1 Note that, in general, the evolution of the state of a time-delayed system involves two actions. The first equation of (28) is a rule for the extension to the future. Secondly, the present state shifts through the history as time progresses:

Formally, consider the strongly continuous semigroup , which solves (28), that is . Associated with is the abstract differential equation where is the generator defined as and corresponding domain such that this limit exists. Hence, With the generator in this form, we can indeed see that the solution is generated by shifting and extending. Namely, the action of is differentiation, which is the generator of the shift semigroup . The extension component, on the other hand, is incorporated in the domain of , suggesting that the solution of (28) is generated by shifting only those functions that satisfy the differential equation. Here, we stress the fact that the appearance of the differential equation as a condition on the domain is cumbersome. If we, at this point, were to proceed with a linear parameter-dependent perturbation of the differential equation, the domains of definition of the solution spaces would change, obstructing standard bifurcation analysis.

This particular problem is overcome if we study the problem in a new and ‘larger’ space . The space (pronounced: sun–star) is best thought of as the double dual space of with additional canonical restrictions geared towards maintaining strong continuity of the semigroup. In our case , which cannot be represented in terms of known functions or measures. Yet, is canonically embedded in via given by —here we exploited the fact that , which suffices for our purpose.

The flow on is generated by . If , then and . Comparing with as in (30), we see that the condition on the domain, which deals with the right hand side of the differential equation, is transformed into an action of the operator (i.e. the first component of ). In this setting, we can readily perturb the linear system (28) without altering the space or the domain of .

Centre manifold and homological equation

We proceed to the non-linear model which generalises the neural field (1)  +  (2a): for some Lipschitz continuous . Note that the differential equation is an equation in and, upon assuming , we can extend (31) to an equation in The last equation defines an abstract differential equation on , where is a true non-linearity in its first component.

Whenever the system (31) is at a critical point, i.e.  has eigenvalues with vanishing real part, then there exists a locally invariant centre manifold , which is tangent to the critical eigenspace at the origin. Moreover, if there are also critical eigenfunctions , then there exists a smooth such that . We expand where is a multi-index of length and . (Note that, in the case of complex eigenvalues, the expansion is usually chosen differently; see Appendix C.) On the centre manifold, the system satisfies some ODE in that is equivalent to the normal form with unknown critical normal form coefficients . Due to invariance of the centre manifold, we can restrict the dynamics to the centre manifold, i.e.  , and substitute (32) to obtain the homological equation Next, admits the expansion where and are the bi- and tri-linear operators corresponding respectively to the second and third derivatives of , cf. Appendix C.1. Note that, since is zero in its second component, also and will be zero in their second components.

Finally, (34) allows us to find expressions for the critical normal coefficients. Equating coefficients of powers of on both sides, one recursively obtains expressions for and . In particular, one encounters linear systems of two forms which can be solved explicitly: We note that the two foregoing statements are not explicitly stated in the original works  [29], [30], but are discussed and implemented in §4.4 and §4.3 respectively.

For a regular value (i.e.  ), and with the resolvent as in Theorem 3.

For an eigenvalue of (i.e.  ), with corresponding eigenfunction , and , Fredholm solvability yields a condition for where is a punctured disk around , with sufficiently small radius such that .

Bifurcations and normal form computation

In this section we combine the various results from the foregoing sections to identify and classify several bifurcations in the model (1)  +  (2a). Using the stability results from Section  3, we identify the precise parameter values of the bifurcation points. Since we know the mode number(s) involved in the bifurcation, the planforms in Section  4 inform us about the relevant symmetries and the corresponding normal form equation. This normal form is, in turn, at the heart of the centre manifold reduction reviewed in Section  5, which enables us to computate the relevant coefficients in the normal form. All that remains at this point, is to analyse the behaviour of the normal form’s low-dimensional dynamics. While these have been studied and documented in great detail for bifurcations of codimension 1 and 2 in systems without symmetry, e.g.  [49], the normal forms of symmetric bifurcations are numerous and, therefore, lacking a clear overview. As a consequence, we will devote part of this section to the analysis of the low-dimensional system that arises from a double Hopf bifurcation with a special symmetry.

More concretely, we study the following two bifurcations: A third case, the Hopf bifurcation due to only, is treated in Appendix C.3.

Hopf bifurcation where , and

double Hopf bifurcation for mixed interactions, where and .

Hopf bifurcation,

As a starting point we will analyse the simplest Hopf bifurcation in the system, namely the one without symmetry. The reasons for this are twofold: Firstly, Fig. 1 reveals that for the stability region is largely bounded by the Hopf bifurcation corresponding to . Secondly, the centre manifold is more accessible because we can apply the results for the generic Hopf bifurcation given in  [29]. We are, however, in a better position than the authors of the original work: since our model allows its eigenfunctions to be expressed analytically, we are able to find a closed expression of the first Lyapunov coefficient .

If for , is a simple root of and for all other and , then the system undergoes a Hopf bifurcation w.r.t. the mode . In this case, the eigenvalue has both algebraic and geometric multiplicity equal to 1, and its eigenfunction is . Since the eigenfunction is constant across space, the oscillations originating from this bifurcation will also be homogeneous across space. Therefore, we denote this pattern as a bulk oscillation.

Because the eigenvalue has multiplicity one, the normal form for this Hopf bifurcation is given by: with the coefficient determining the criticality. Here denotes the real part of the critical eigenvalue in the neighbourhood (in parameter space) of the bifurcation; clearly at the bifurcation. Moreover, the first Lyapunov coefficient is defined as . Using the techniques described in Section  5, we are able to find an explicit expression for this coefficient (cf. Appendix C.2): where Now, a negative (positive) sign of the first Lyapunov coefficient corresponds with a supercritical (subcritical) bifurcation.

Using the parametric expression for the Hopf bifurcations in parameter space, as formulated in Section  3, and the choice of as in (15), we are able to determine the first Lyapunov coefficient along these curves. In particular, we identify parameters for which the first Lyapunov coefficient vanishes, corresponding to a generalised Hopf bifurcation. Although we do not compute the second Lyapunov coefficient, we investigate this bifurcation numerically; see Fig. 6. The bistability, as observed in Fig. 6 (B3), between a focus (red) and a limit cycle (blue) suggests that the second Lyapunov coefficient is negative. Note that, since the pattern is homogeneous across space, it suffices to plot time series at only one point on the sphere.

Fig. 6

Generalised Hopf bifurcation occurs for and . (A1) shows that for the Hopf bifurcation is supercritical. Time series for two different initial conditions (red and blue) reveal one stable focus for (B1) while for a stable limit cycle is observed (B2). For , the criticality of the Hopf bifurcation has changed, resulting in a multistable regime (A2). Indeed, in (B3) simulations show that for there exist two stable solutions: a stable focus (red) and a stable limit cycle (blue). For , the stable focus has destabilised and only the stable limit cycle remains. Parameter values: , , , and , , and . (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

Double Hopf bifurcation, and

If we, starting at the supercritical Hopf bifurcation for , vary a different parameter, we will not observe a change of sign of the first Lyapunov exponent. Instead, we find that another pair of eigenvalues, corresponding to , passes through the imaginary axis. We analyse this double Hopf bifurcation using symmetry techniques.

For , let be such that is a simple root of and, moreover, has no other roots for and . There is one eigenfunction related to , while there are three corresponding to , namely for .

In the non-resonant case, i.e.  for all with , symmetry arguments yield the truncated normal form up to cubic order: where Expressions for the critical coefficients, , , , , of the model (2b) are given in Appendix C.4. Moreover, setting gives the centre manifold: It is important to stress the fact that, since we are neglecting terms of order 4 and higher, we can only study the ‘simple’ case  [49, Section 8.6.2]. Therefore, our analysis is only valid whenever .

We continue to identify possible solutions of (39). By computing all isotropy subgroups of for the representation on we find that there are three -axial subgroups which are numbered 1–3 in Table 3. This table gives all isotropy subgroups, , for this representation, along with a representative fixed point subspace, , for one particular choice of set of generators of the subgroup. Branches of solutions with these symmetry types are guaranteed to exist by the equivariant Hopf theorem. The patterns corresponding to solution types 1–3 are respectively bulk oscillations, travelling waves and standing waves.

Table 3

Isotropy types for possible solutions of (39) together with their representation in the reduced system (42).

	Σ	H	K	Fix(Σ)	|w|^2	|z02−2z_−1z_1|^2	|z|^2
1	O(3)˜	O(3)	O(3)	(w;0,0,0)	> 0	0	0
2	SO(2)˜	SO(2)×Z2c	Z2−	(0;z_−1,0,0)	0	0	> 0
3	O(2)˜	O(2)×Z2c	O(2)^−	(0;0,z_0,0)	0	> 0	> 0
4	Z_2˜	Z_2×Z2c	Z2−	(0;z_−1,0,z_1)	0	> 0	> 0
5	1˜	Z2c	1	(0;z_−1,z_0,z_1)	0	≥ 0	> 0
6	O(2)^−	O(2)^−	O(2)^−	(w;0,z_0,0)	> 0	> 0	> 0
7	Z2−	Z2−	Z2−	(w;z_−1,0,z_1)	> 0	> 0	> 0
8	1	1	1	(w;z_−1,z_0,z_1)	> 0	≥ 0	> 0

Other solutions to (39) may exist depending on the values of parameters. These correspond to isotropy types which have fixed point subspaces of dimension greater than 2 and are types 4–8 in Table 3. We make the following observations about the solutions with such symmetries. We have two types of solutions with only modes (using the numbering in Table 3):

This solution has contributions from two different modes. If a solution with this symmetry were to exist generically (i.e. when ) it would have and for some fixed phase shift .

This solution has contributions from two different modes. Thus . In this case the equations for the amplitudes and phases , do not decouple. However the equations have effective dimension 4 since they only depend on the combination of phases .

There could also exist three types of mixed mode solutions. These can only have spatial symmetries since we assume no resonance between the and modes.

Removing the time shift symmetries of the standing wave solution with symmetry, we are left only with the symmetry which allows for a nonzero amplitude of the mode.

Solutions with just a reflection symmetry include (for the case of reflection in the plane) . In this case the equations for the amplitudes again decouple but a solution with and only exists generically in the case that and for a constant phase shift .

For solutions with no symmetry we have Here the amplitude and phase equations are coupled but as for solution type 5 they depend only on and so the effective dimension is 6 rather than 8.

Although it is not possible to decouple the system in its original form, a different set of coordinates does. Indeed, the transformation yields a decoupled system of ODEs: where and denote the real parts of and respectively. Fixed points of the system can be related to all of the eight cases described above; Table 3 shows this relation.

It is straightforward to show that, generically, the system (42) has at most 6 equilibria in the first octant . (Note that solutions outside the first octant are irrelevant, for are non-negative amplitudes.) For parameters corresponding with the double Hopf bifurcation in the full non-linear model (2a), we find that the real parts of all normal form coefficients in (39) are negative. The phase plane of the corresponding reduced system (42) is shown in Fig. 7, which reveals the occurrence of all 6 steady states. Continuing with a linear stability analysis of these six equilibria, we find that only points ii and iii are asymptotically stable (i.e.  and respectively). Therefore, we conclude from Table 3 the existence of a multistable regime in which we observe either bulk oscillations, travelling waves, or a mixture of different 1-modes (cases 1, 2 and 5 respectively).

Fig. 7

Nullplanes of (42) in the plane, depicting , and in red, green and blue respectively. Note that the plane is a -nullplane, such that shown intersections of - and -nullplanes correspond with steady states; i.e. points i–iv. For , the -nullplane is given by the dashed green line, yielding two more equilibria in the unseen dimension: v and vi, for which . Parameters as follows: and for all . (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

This result is confirmed using direct simulations of the model near the bifurcation. We find two stable solutions, namely bulk oscillations (case 1) and rotating waves (case 2); the latter is depicted in Fig. 8. We have not identified a stable solution consisting of mixed 1-modes (case 5).

Fig. 8

Emergence of a rotating wave in the presence of multistability. Direct simulation of the initial condition for evolves into the rotating wave . The time-course of the sphere’s centre of mass is shown in , projected in the -plane. The initial condition, marked with the red cross, is near the (asymptotically) unstable standing wave pattern, causing the system to remain close to this pattern—as is shown by the centre of mass moving oscillating in an almost straight line. In the course of the simulation, the centre of mass approaches its limit cycle (red circle), which corresponds to the rotating wave solution. This stable solution co-exists with stable bulk oscillations (not shown). Parameter values: , , , , , and , , , and . Simulated with a mesh of 5120 triangles. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

Discussion

In this paper we have analysed pattern formation in a spherical model of brain activity and shown how a rich repertoire of spatio-temporal states can be supported in neural models of Nunez type that contain only simple representations for anatomical connectivity, axonal delays and population firing rates. Even though these models are naturally formulated as nonlinear evolution equations of integro-differential form with delays, which have been little studied, they are amenable to similar analyses used for studying pattern-formation in PDEs. Perhaps surprisingly this is the first time that a combination of linear stability analysis, symmetric bifurcation theory, centre manifold reduction and direct numerical simulations have been combined to explore such a popular model for EEG generation. Most certainly this is because it has proved easier to study versions of the model on the line  [50] or the plane  [51], despite the obvious motivation to study the model on more brain like topologies.

Importantly, our study is the first one to carry out a detailed centre manifold reduction on the full integral formulation of the model, including transmission delays, in a two-dimensional setting—the one-dimensional case is discussed in  [28], [29], [30]. As a result, this new work has shed light on the importance of delays in generating patterns with a high degree of spatial structure, as well as developed the bifurcation theory that can be used to ascertain the emergence of a given symmetric structure via the destabilisation of a homogeneous steady state. Our approach has also allowed us to explicitly pin-point the conditions for codimension-2 bifurcations, where two or more distinct patterns can be excited and then subsequently interact.

Secondary bifurcations, especially those related to multi-stability, have received considerable attention in the modelling of EEG; especially in the context of epileptic seizures. The generalised Hopf bifurcation, with bistability between rest and oscillation, has a pivotal position in explaining ictal transitions in models of cortical columns  [52], [53]. Since these models have been studied primarily in the absence of a spatial component, it has remained unclear what the impact of both spatial structure and transmission delays would be on oscillations and synchrony in the model. In this article we have shown that the generalised Hopf bifurcation can still occur in the extended setting, enabling the model to switch between rest and a fully synchronous periodic solution. The double Hopf bifurcation, which yields a type of multistability where multiple stable periodic solutions can exist, was studied in a neural field context as early as 1980  [54]. Therein, the authors emphasised that the appearance of quasi-periodic behaviour, which occurs in a special case of the bifurcation as a result of mixing between two stable oscillations, could explain the transition between the tonic and clonic seizure states. Being largely theoretical, their work does not provide a methodology for computing or classifying these transitions. Although we have been able to identify these quasi-periodic oscillations in our work (points IV and VI in Fig. 7), we have been unable to find parameters for which they would be stable—it remains an open question whether such parameters exist for our model. Other work on secondary bifurcations in one-dimensional neural field models can be found in  [55].

Of course, the Nunez model is also able to generate a whole host of more exotic behaviour in regimes where our analytical approaches have less sway. For example, Fig. 8 shows the emergence of a large scale rotating wave which is reminiscent of those reported in EEG studies of schizophrenic patients  [56]. In these instances the development of our bespoke numerical scheme pays further dividends since it is not restricted to perfectly spherical models. The scheme is sufficiently general to handle real folded cortical structures with more detailed white matter fibre tractography data, of the type that is increasingly available in public repositories such as the Human Connectome Project  [57]. Indeed one natural extension of the work presented in this paper is the development of a primarily computational model of brain activity that can incorporate more biological detail, such as folded cortical hemispheres  [58], heterogeneous connection topologies  [59] and a distribution of axonal speeds  [60]. Such a model is relevant to interpreting modern whole brain neuroimaging signals, and its exploration would set the scene for formulating the relevant mathematical questions about how best to understand the behaviour of a complex brain model. For example, it would be interesting to explore how wave dispersion and interaction on a folded cortex affects the frequency of emergent rhythms and the dependence of these on underlying activity that is primarily either in the form of travelling waves or standing waves. Even before developing such a programme the model explored here has other solutions that are of interest to the neuroscience community in the form of localised states (such as spots), often invoked in the context of working memory and easily established for a steep sigmoidal firing rate and Mexican-hat connectivity, see for example  [61]. Once again progress in this arena might take as a starting point ideas recently developed for the study of spots for PDEs on a sphere  [62]. These and other topics, including the response of the model to forcing, will be reported upon elsewhere.