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 ∼10 8 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 http 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, latticestructure 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 integrodifferential 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 ∼10 6 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 10 10 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 10 10 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 10 9 for unprocessed EEG and 10 8 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 ∼800-1500 cm 2 as R = √ A/(4π ) ∼ 8-10 cm. Taking a value for the corticocortical action potential speed in the range v ∼ 6-9 m/s 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 f ∼ 13-25 Hz, using the rela- 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 u = u(t, r) ∈ R where r is a point on the surface of a sphere and t ∈ R + . We shall consider simple neural field models that, after rescaling time and space, can be written in the form with either We also have to specify the initial condition: where η ∈ X , the state space defined below.
Here Ω := S 2 is the surface of the unit sphere in R 3 , dr ′ the integration measure and • denotes function composition, with f a firing rate function. The weight distribution w(r, r ′ ) specifies the anatomical connectivity between points r and r ′ , whilst τ (r, r ′ ) 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: • the firing rate function f is smooth and bounded on R, • the domain Ω := S 2 is the unit sphere in R 3 and the corresponding metric d is the great circle distance, • the connectivity kernel w ∈ C 0,α (Ω × Ω), the Banach space of Hölder continuous functions with exponent α, 1/2 < α ≤ 1 on Ω × Ω.
• the delay function τ ∈ C 0,α (Ω × Ω) is non-negative and not identically zero, • the maximal delay h := sup{τ (r, r ′ ) : r, r ′ ∈ Ω}, • the spatial space V := C 0,α (Ω, C) the Banach space of Hölder continuous functions with exponent α, 1/2 < α ≤ 1, with the standard norm: • the state space X : These assumptions are nearly identical to those in the seminal work [29], apart from one: where the original work sets the spatial space V 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. 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 G : X → V as We have the following lemma.
Proof. Given u ∈ X , we consider two points r andr of Ω and write Using the definition of the operator G, 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 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 F defined by (5) is Lipschitz continuous. We prove the following lemma. Proof. If Ψ andΨ are elements of X and r,r ∈ Ω, consider for some positive constant C 1 . Consider now We also have for some positive constant C 2 . Combining (6) and (7) we obtain for some positive constant C , i.e. F is Lipschitz continuous.

Spherical geometry
For the remainder of this paper we will assume that r = r(θ , φ) ∈ S 2 is a point on the sphere with polar angle θ ∈ [0, π], azimuthal angle φ ∈ [0, 2π ) and radius 1 and similar for r ′ . Furthermore, we make use of the complex-valued spherical harmonics Y m n (r) of degree n ≥ 0 and order |m| ≤ n, 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.

Theorem 1 (Spherical Harmonics Expansion [38, Thm. 5]). Let v ∈ V , then the spherical harmonics expansion of v converges uniformly to
The coefficients v m n are given by projections on the basis functions: where the overline denotes complex conjugation.
Proof. We give a very short proof for completeness. Let P N be the We shall consider a homogeneous neural field, where both the weight kernel w(r, r ′ ) and transmission delays τ (r, r ′ ) depend on the relative position of r and r ′ . Naturally, w and τ are chosen as functions of distance along the surface, but more generally we set w(r, r ′ ) := w(r · r ′ ) and τ (r, r ′ ) := τ (r · r ′ ). (10) 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. cos(α) = r · r ′ . This leads to the following theorem.

Theorem 2. Let g be a Hölder continuous function on
where P n is the Legendre polynomial of degree n, converges to G(r, r ′ ), uniformly on Ω × Ω, when N → ∞.
Proof. It is easily checked that the hypothesis on the function g implies that G is Hölder continuous on Ω × Ω with exponent α: An easy extension of Theorem 1 shows that the series when N → ∞, where the coefficients G mm ′ nn ′ are given by the projections: Since g is continuous on [−1, 1] hence in L 1 (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 f , w 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 β > 0 and threshold δ.
A natural choice for the connectivity kernel is with σ 1 > σ 2 > σ min > 0 and J 1 J 2 < 0. For J 1 + J 2 > 0 we have a wizard-hat shape, whilst for J 1 + J 2 < 0 we have an inverted wizard-hat shape.
We have the following lemma.
Proof. The space of Hölder continuous functions with the same exponent being a vector space, it is sufficient to check that ww : for some positive constant C and for all s, s ′ 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 τ 0 and a constant propagation speed c 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. (17) is Hölder continuous with exponent α, 1/2 < α ≤ 1.

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 E n (λ) 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 u(t, r) =û of (1) + (2a) satisfieŝ 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. w 0 = 0,û = 0 is the only homogeneous equilibrium.
Linearising (1) where κ = f ′ (û). Solutions of this linear equation are separable and we set ξ (t, r) = e µt q(r). In this case q(r) satisfies the linear is the characteristic function with G(r, r ′ ; µ) = w(r · r ′ ) exp (−µτ (r · r ′ )). Note that the particular structure of G allows application of Theorem 2, which yields coefficients G n (µ)-see Appendix B. Indeed, Lemmas 3 and 4 show that the functions w and τ , [−1, 1] → R are Hölder continuous with exponent α, 1/2 < α ≤ 1. It is then easily verified that the function [−1, 1] → R, s → w(s)e −µτ (s) is also Hölder continuous with exponent α, and the beginning of the proof of Theorem 2 shows that the function solutions. This occurs when for some n ≥ 0. In this case, (19) has 2n + 1 solutions (i.e. eigenfunctions) of the form ξ m (t, r) = e λt Y m n (r), m = −n, . . . , n.
Hence, the algebraic and geometric multiplicity of the eigenvalues are the same.

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.
If both q and y are expanded in spherical harmonics, see Theorem 1, then coefficients q m n are given by Proof. Existence and uniqueness of the solution are shown in [29,Prop. 14]. To identify the solution, we start with ∆(µ)q = y, substitute (20) and expand G using Theorem 2 where integration and summation are interchanged. Next, we multiply both sides by Y m n (r), integrate over the domain w.r.t. r 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 Re λ < 0 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 (τ = 0), all eigenvalues are real and given explicitly by λ n = −1 + κw n . 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) and (17) we are able to find explicit expressions for the coefficients G n (λ), of the form with h n (λ; σ ) according to Appendix B. Using (21), one can identify the parameters (J 1 , J 2 ) at which the homogeneous steady state undergoes a fold or transcritical bifurcation-in this case λ = 0. Indeed, defines a line in the (κJ 1 , κJ 2 )-plane which corresponds to candidate fold/transcritical bifurcations corresponding to the mode number n. In Fig. 1 whose solutions are expressed parametrically in terms of ω by equating real and imaginary parts: These parametric curves are plotted for increasing n in Fig. 1 as solid lines.
Furthermore, for small κJ 1 and κJ 2 , the linear equation (19) is dominated by the term −v(t, r) such that solutions near the origin of the (κJ 1 , κJ 2 )-plane are asymptotically stable. Since instabilities only occur at the curves (23) and (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 w 0 = J 1 h 0 (0; σ 1 ) + J 2 h 0 (0; σ 2 ) and, hence, κ has an implicit dependency on J 1 and J 2 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 {J 1 h 0 (0; σ 1 )+ J 2 h 0 (0; σ 2 ) = C , C ∈ R}, the fixed point structure -and therefore κ -remains unaltered. This collection of parallel lines is illustrated in Fig. 1 for various C . One line is of particular importance: the line w 0 = 0 corresponds to a balanced kernel, such thatû = 0 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 n = 0. Only a small part of the stability region, as shown in the inset, is bounded by curves relating to bifurcations of higher degree in n. 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 τ 0 = 0, 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 n ≤ 5. The close-up shows the part of the stability region which is bounded by the Hopf curve for degree n = 4.
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 λ = ρ + iω in the (ρ, ω) plane for n ≤ 5, as determined by solving (20) numerically. A pair of complex eigenvalues, corresponding to n = 4 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).
In Fig. 4 we show a direct simulation of an instability to a pattern state with n = 4 and ω ̸ = 0, 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].

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 cc stands for complex conjugate, with n c and ω c determined from the spectral equation (21) such that E n c (iω c ) = 0 while Re λ < 0 for all other λ ∈ σ . Sufficiently close to the bifurcation point we expect certain classes of solutions to emerge that break the O(3) 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 In order to determine the structure of the equivariant normal form on V n c ⊕ V n c to a given order we use the action of elements of Symmetry can also be used to determine branches of periodic solutions ofż = f(z). Let, without loss of generality, z(t) be a periodic solution ofż = f(z) with period 2π . The where H is a subgroup of O(3) and θ : H → S 1 is a group homomorphism. Elements of O(3) can be thought of as spatial symmetries whilst elements of S 1 are temporal symmetries acting on periodic solutions by a phase shift. An element (h, θ (h)) ∈ O(3) × S 1 is a spatial symmetry if θ (h) = 0 and a spatiotemporal symmetry if θ (h) ̸ = 0. The spatial symmetries form a normal subgroup K = ker(θ ) of H and the quotient group H/K is isomorphic to a closed subgroup of S 1 (i.e. 1, Z n (n ≥ 2) or S 1 ). The C-axial isotropy for values of n c between 1 and 6 is given in Table 1 which is reproduced from [45]. All notations for the subgroups of O(3) here and throughout are consistent with that in [44,45]. Observe for example that if the spherical harmonics of degree n c = 4 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 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 Y m n c and Y −m n c . The resulting standing wave can be thought of as arising due to interference between the waves travelling in opposite directions around the sphere.

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 X is the state space of the delayed equation (1) and elements x t ∈ X relate to the system's history via x t (θ ) := x(t + θ ) ∀t ≥ 0, θ ∈ [−h, 0]. Next, consider the linear system with initial condition η ∈ X and linear operator L : X  → V . For our particular system, L is given by 1 (Lη)(r) = −η(0, r) + κ  Ω w(r · r ′ )η(−τ (r · r ′ ), r ′ )dr ′ ∀η ∈ X , ∀r ∈ Ω. (29) Note that, in general, the evolution of the state x t ∈ X of a timedelayed 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:  Table 2 The C-axial subgroups of O(3) × S 1 for the natural representations on Formally, consider the strongly continuous semigroup T (·), which solves (28), that is x t = T (t)η, ∀t ≥ 0. Associated with T is the abstract differential equation With the generator in this form, we can indeed see that the solution is generated by shifting and extending. Namely, the action of A is differentiation, which is the generator of the shift semigroup T 0 . The extension component, on the other hand, is incorporated in the domain of A, 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 X ⊙ * . The space X ⊙ * (pronounced: sun-star) is best thought of as the double dual space of X with additional canonical restrictions geared towards maintaining strong continuity of the semigroup. In our case X ⊙ *  j : X  → X ⊙ * given by jη = (η(0), η)-here we exploited the fact that V × L ∞ ([−h, 0]; V ) ⊂ X ⊙ * , which suffices for our purpose.
The flow on X ⊙ * is generated by A ⊙ * . If η ∈ C 1 ([−h, 0]; V ) ⊂ X , then jη ∈ D(A ⊙ * ) and A ⊙ * jη = (Lη,η). Comparing A ⊙ * with A 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 A ⊙ * jη). In this setting, we can readily perturb the linear system (28) without altering the space X ⊙ * or the domain of A ⊙ * .

Centre manifold and homological equation
We proceed to the non-linear model which generalises the neural field (1) + (2a): for some Lipschitz continuous F : X  → V . Note that the differential equation is an equation in V ⊂ V * * and, upon assuming The last equation defines an abstract differential equation on X ⊙ * , where R : X  → X ⊙ * is a true non-linearity in its first component.
where ν is a multi-index of length N and h ν ∈ X . (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 R N that is equivalent to the normal forṁ with unknown critical normal form coefficients g ν ∈ R N . Due to invariance of the centre manifold, we can restrict the dynamics to the centre manifold, i.e. x t = H(z(t)), and substitute (32) to obtain the homological equation Next, R admits the expansion where B and C are the bi-and tri-linear operators corresponding respectively to the second and third derivatives of F , cf. Appendix C.1. Note that, since R is zero in its second component, also B and C will be zero in their second components. Finally, (34) allows us to find expressions for the critical normal coefficients. Equating coefficients of powers of z on both sides, one recursively obtains expressions for h ν and g ν . In particular, one encounters linear systems of two forms which can be solved explicitly: • For µ a regular value (i.e. µ ̸ ∈ σ (A)), η ∈ X and v ∈ V with ∆ −1 (µ) the resolvent as in Theorem 3.
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.

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.
A third case, the Hopf bifurcation due to E 1 only, is treated in Appendix C.3.

Hopf bifurcation, E 0
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 τ 0 = 0 the stability region is largely bounded by the Hopf bifurcation corresponding to n = 0. 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 l 1 .
If for ω 0 > 0, iω 0 is a simple root of E 0 and E n (iω) ̸ = 0 for all other n and ω ≥ 0, then the system undergoes a Hopf bifurcation w.r.t. the mode n = 0. In this case, the eigenvalue λ = iω 0 has both algebraic and geometric multiplicity equal to 1, and its eigenfunction is ξ (t, r) = e iω 0 t Y 0 0 (r) = e iω 0 t / √ 4π . 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 g 21 the coefficient determining the criticality. Here ρ 0 denotes the real part of the critical eigenvalue in the neighbourhood (in parameter space) of the bifurcation; clearly ρ 0 = 0 at the bifurcation. Moreover, the first Lyapunov coefficient is defined as Re g 21 . 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 f 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.

Double Hopf bifurcation, E 0 and E 1
If we, starting at the supercritical Hopf bifurcation for n = 0, 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 n = 1, passes through the imaginary axis. We analyse this double Hopf bifurcation using symmetry techniques.
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 Re a 1 Re b 1 > 0.
We continue to identify possible solutions of (39). By computing all isotropy subgroups of O(3) × S 1 for the representation on we find that there are three C-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, Fix(Σ), 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. 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 n = 1 modes (using the numbering in Table 3): 4. This solution has contributions from two different n = 1 modes.
If a solution with this symmetry were to exist generically (i.e. when Re b 2 ̸ = 0) it would have z −1 = R(t)e iφ(t) and z 1 = R(t)e i(φ(t)+ψ ) for some fixed phase shift ψ. 5. This solution has contributions from two different n = 1 modes.
Thus w = 0, z = (R −1 (t)e iφ −1 (t) , R 0 (t)e iφ 0 (t) , R 1 (t)e iφ 1 (t) ). There could also exist three types of mixed mode solutions. These can only have spatial symmetries since we assume no resonance between the n = 0 and n = 1 modes.  1 (t) ). In this case the equations for the amplitudes again decouple but a solution with R −1 ̸ = 0, R 1 ̸ = 0 and r ̸ = 0 only exists generically in the case that R −1 = R 1 and φ −1 = φ 1 + ψ 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 2φ 0 − φ 1 − φ −1 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 (x 1 , x 2 , x 3 ) = (|w| 2 , |z 2 0 −2z −1 z 1 |, |z| 2 ) ∈ R 3 yields a decoupled system of ODEs: whereã i andb i denote the real parts of a i and b i respectively. Fixed points of the system x 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 x 1 , x 2 , x 3 ≥ 0. (Note that solutions outside the first octant are irrelevant, for x i are nonnegative 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. x = (−µ 1 /a 1 , 0, 0) and x = (0, 0, −µ 2 /b 1 ) 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). 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).

Discussion
In this paper we have analysed pattern formation in a spherical model of brain activity and shown how a rich repertoire of spatiotemporal 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 multistability, 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 quasiperiodic 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 Fig. 8. Emergence of a rotating wave in the presence of multistability. Direct simulation of the initial condition A1 for t ≤ 0 evolves into the rotating wave A2. The time-course of the sphere's centre of mass is shown in B, projected in the x, y-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: α = 1.08, β = 4, δ = 0.1, τ 0 ≈ 3.483, c = 1, σ 1 = 1 and σ 2 = 1/2, J 1 ≈ 1. 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 Mexicanhat 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.
The spherical harmonics are represented by the functions: where P m is the associated Legendre function. The spherical harmonics satisfy They also satisfy the orthogonality condition

Appendix B. Basis coefficients
For given G(r, r ′ ; µ) = w(r · r ′ ) exp(−µτ (r · r ′ )), the expansion into spherical harmonics, see Theorem 2, yields coefficients: For the evaluation of the integrals, we apply the following result. , This yields, moreover, the two-term recurrence relation: Proof. First, consider for n ≥ 0 H n (a) :=  π 0 e at P n (cos t)dt.

(B.2)
We will use mathematical induction to show that H n satisfies the following recurrence relation H n+1 (a)H n (a) = 1 − e 2aπ a n + (n Starting with a recurrence relation for Legendre polynomials, we find (n + 1)P n+1 (cos t) + nP n−1 (cos t) = (2n + 1) cos t P n (cos t), where we have expressed cos t as a complex exponential. Another identity for Legendre polynomials gives P ′ n+1 (cos t) − P ′ n−1 (cos t) = (2n + 1)P n (cos t),  π 0 e at sin t where the left hand side is obtained via integration by parts and the right hand side via expressing sin t as a complex exponential. Adding (B.4) and i times (B.5) yields This identity will be the starting point for the mathematical induction.
Assume that (B.3) holds for all n ≤ N for some N > 0; we will prove that it also holds for n = N + 1. Combining (B.3) at n = N and at n = N − 1, yields

C.1. Multilinear operators B and C
In the remainder of Appendix, we denote the second and third derivatives of the full non-linear system by B and C respectively.
This eigenvalue has multiplicity 1 and the corresponding eigenfunction is ξ (t, r) = e iω 0 t Y 0 0 (r) = e iω 0 t / √ 4π . This case corresponds with the ordinary Hopf bifurcation formulated in [29] and, hence, we can apply their result directly.
The quadratic terms in the normal form yield equations The last step results from properties of the spherical harmonics. Hence: where we have defined Q n (µ) := G n (µ) µ + 1 − κG n (µ) .
In a similar manner, we solve (C.3b) The cubic terms of the Hopf bifurcation give rise to the equation because B and C are zero in the second component. We apply (35b) to obtain Expansion of both sides into spherical harmonics results in an equation for each coefficient From the characteristic equation it follows that iω 0 + 1 − κG 0 (iω 0 ) = 0, and we also have that µ + 1 − κG n (µ) ̸ = 0 for all n > 0. Therefore, in the second equation the integrand is analytic for all n > 0 and these contour integrals vanish. The first equation is evaluated using Cauchy's integral formula: where G ′ 0 denotes the derivative of G 0 . Next, (C.7) yields v(r) = 1 16π 3/2 Therefore, we find the critical normal form coefficient to be
Following the same procedure as in Appendix C.3, we construct the homological equation. At quadratic order, we obtain: From the solvability conditions, we find expressions for the critical normal form coefficients: