Symmetry-Breaking in a Rate Model for a Biped Locomotion Central Pattern Generator

The timing patterns of animal gaits are produced by a network of spinal neurons called a Central Pattern Generator (CPG). Pinto and Golubitsky studied a four-node CPG for biped dynamics in which each leg is associated with one flexor node and one extensor node, with Z2 × Z2 symmetry. They used symmetric bifurcation theory to predict the existence of four primary gaits and seven secondary gaits. We use methods from symmetric bifurcation theory to investigate local bifurcation, both steady-state and Hopf, for their network architecture in a rate model. Rate models incorporate parameters corresponding to the strengths of connections in the CPG: positive for excitatory connections and negative for inhibitory ones. The three-dimensional space of connection strengths is partitioned into regions that correspond to the first local bifurcation from a fully symmetric equilibrium. The partition is polyhedral, and its symmetry group is that of a tetrahedron. It comprises two concentric tetrahedra, subdivided by various symmetry planes. The tetrahedral symmetry arises from the structure of the eigenvalues of the connection matrix, which is involved in, but not equal to, the Jacobian of the rate model at bifurcation points. Some of the results apply to rate equations on more general networks.


Introduction
Animals use many different patterns of locomotion, known as gaits, Gambaryan [1], Muybridge [2].Mathematically, gaits are primarily characterised by the pattern of phase shifts among distinct legs, McGhee [3], although other features such as the duty factor (fraction of the gait cycle for which the foot remains on the ground, Hildebrand [4]) are also relevant.
It is generally believed, with much circumstantial evidence and a growing amount of more specific evidence, that gait patterns in mammals, reptiles, and indeed many other animals are determined by a relatively simple network of neurons in the spinal cord, known as a Central Pattern Generator (CPG).In a series of papers, Kopell and Ermentrout [5][6][7][8] studied models of lamprey and fish CPGs based on chains of phase oscillators, paying particular attention to differences between diffusive coupling, which vanishes for oscillators that are in phase with each other, and synaptic coupling, where it does not.They concluded that synaptic coupling is required to produce locomotion.
Collins and Stewart [9][10][11][12] analysed gaits as time-periodic states created by Hopf bifurcation (Hassard et al. [13]) from a rest state.Such gaits are said to be primary, and in quadrupeds they include the main "symmetric" gaits such as walk, bound, pace, trot, and pronk.The pronk is symmetry-preserving; the other gaits are symmetry-breaking.More complex gaits such as canter and transverse and rotary gallop are said to be secondary, Buono [14,15] and Buono and Golubitsky [16].Collins and Stewart used methods from symmetric bifurcation theory [17,18] to classify gait patterns for a variety of hypothetical CPG architectures.They discussed the gaits of quadrupeds [9,11] and hexapods [10], and applied the same methods to more general biological oscillators [12].
Golubitsky et al. [19,20] observed that none of the architectures that Collins and Stewart considered is satisfactory for quadruped gaits, because any four-node network supporting walk and either trot or pace must support all three gaits.Trot and pace are then dynamically conjugate, and therefore coexist for the same parameter values in any model with those architectures.This conflicts with observations.By stating a list of properties that a suitable CPG should possess, they deduced the simplest symmetry class for a CPG consistent with observations.This network has eight nodes and Z 4 × Z 2 symmetry, so that two nodes are associated with each leg.A physiological interpretation was suggested: the two nodes associated with a given leg control the timing of flexor and extensor muscle groups.
These results are model-independent, in the sense that no specific equations are employed to model the CPG dynamics.The analysis applies to all such equations, provided a Hopf bifurcation of the appropriate kind occurs and appropriate nondegeneracy conditions are satisfied, as they are generically.Buono [14,15] and Buono and Golubitsky [16] also performed numerical simulations with specific equations to show that this network architecture produces all eight primary gaits with either synaptic or diffusive coupling, for suitable node dynamics and coupling terms.
Golubitsky et al. [19,20] deduced several predictions about gaits, valid for all such models, which appear to be consistent with observations.They briefly discussed analogous CPG architectures for bipeds, hexapods, and 2n-legged organisms for any n.They remarked that the analogous CPG architecture predicts the half-integer wave numbers observed in centipede locomotion.Subsequently, Pinto and Golubitsky [21] studied an analogous CPG network for bipeds, which has Z 2 × Z 2 symmetry.See Figure 1.They used Hopf bifurcation and the H/K Theorem of Buono and Golubitsky [16] to classify the spatio-temporal symmetries of periodic states for this symmetry class, and obtained a list of four primary states arising from Hopf bifurcation, plus seven secondary states predicted by the H/K Theorem and expected to arise through mode interactions between primary states.Here we revisit this four-node CPG, interpreting it as a coupled cell network in the sense of [22][23][24].We investigate a specific system of ODEs consistent with this network structure, which takes the form of a rate model.Such equations model neural activity in terms of the rate at which a neuron fires.They have two main ingredients: a matrix of connection strengths between pairs of nodes, which in this case has three adjustable parameters because of symmetry, and a "gain function" which couples the nodes together.
We consider local bifurcations from a fully synchronous equilibrium u, which arise when the Jacobian J u = Df | u has either a zero eigenvalue (steady-state bifurcation) or a complex conjugate pair of purely imaginary eigenvalues (Hopf bifurcation).The symmetries of the network lead to four distinct symmetry classes of bifurcating branches in each case.The Hopf branches are the primary gaits.
We focus on one approach to these bifurcations which seems a reasonable way to model gaits, and reveals an elegant underlying structure having a richer group of symmetries than the rate equations themselves.Specifically, we consider the first local bifurcation from a fully synchronous equilibrium as an input parameter I increases, assuming that the system reacts quasistatically-that is, it remains in a continuously varying equilibrium state provided that state remains stable, even if other states exist globally.
In this context, we derive explicit conditions on connection strengths for the first local bifurcation to be one of four symmetry types of steady-state bifurcation, or one of the four symmetry types of Hopf bifurcation identified by Pinto and Golubitsky [21] as primary gaits.We also identify the region of parameter space (whose coordinates are the three connection strengths between nodes) in which no local bifurcations occur: the synchronous steady state remains stable and changes continuously whatever value I takes.This region corresponds to small values of the connection strengths.
The specific form of rate equations, in which the argument of the gain function is a linear combination of state variables, implies that the four regions corresponding to primary gaits are connected polyhedra.They all have the same shape and size: each is a truncated pyramid with an equilateral triangle as base.They fit together to form a hollow tetrahedron centred at the origin, as in Figure 2a.This is the complement of a small tetrahedron T * in a larger tetrahedron T * , whose vertices are indicated on the figure: here K is a constant, defined in Theorem (23).The inner tetrahedron T * is similar, with K replaced by k, also defined in Theorem (23).The four primary gait regions are relatively large when K is significantly greater than k.This is the case for reasonable values of the connection strengths and other parameters in the model, so the structure provides a robust way to select specific gait patterns.
Despite the nonlinear form of the equations, it turns out to be possible to derive much of the bifurcation behaviour analytically.The possibility of doing this stems from a pleasant property of a gain function commonly used in rate equations, combined with linearity of the argument of the gain function in the equations.The rate equations themselves have a symmetry group Z 2 × Z 2 of order 4, but the bifurcation behaviour is controlled by a larger group, the symmetries of a tetrahedron, of order 24.We describe how this group arises and why it determines many key aspects of the bifurcation behaviour.Ultimately, it is related to the eigenvalues of the connection matrix, which are closely related to the eigenvalues of the Jacobian for the rate equations, evaluated at a synchronous equilibrium.
Rate models in a network with the same symmetries are analysed in Diekman et al. [25] in a model of binocular rivalry using Wilson networks (the "monkey/text" experiment).They assume a more general form of the gain function, but take β = 0, γ < 0 in Figure 1 because some connections in Wilson networks are always inhibitory, and the others are determined by a learning process.
We do not discuss the stability of the bifurcating states here, to avoid complicating an already lengthy paper.However, simulations suggest that all of the symmetry types of states that we describe can exist stably for some ranges of parameters.
Other types of dynamics can occur in this model, including apparent chaos, but the parameter regions concerned seem to be small.We do not discuss more exotic possibilities here.Curtu [26,27] and Curtu et al. [28] have analysed the bifurcations of a rate model for a network with two identical nodes (Z 2 symmetry) in considerable detail.Even in this case there is a richer range of dynamic behaviour than equilibria and periodic states arising from local bifurcation.Note, however, that they assume that certain connections are inhibitory, so the corresponding connection strength is negative.This restriction implies in particular that there is a unique fully synchronous steady state for all parameter values under consideration.This is not the case when the signs are not restricted, as we show in Section 4. This is the first in a projected series of papers.Others will discuss stability issues and investigate quadruped and hexapod gaits from a similar point of view.New features enter into the analysis, in part because the eigenvalues of the connection matrix can now be complex.

Structure of the Paper
In Section 2 we introduce gain functions, the connection matrix, and rate models.We establish some key properties of the family of gain functions used here.Section 3 describes the CPG network used to model biped gaits, sets up the corresponding rate equations, and discusses parameter symmetries used later to simplify the calculations.We also describe the four primary gaits hop, jump, run, and walk in the context of symmetric Hopf bifurcation, following Pinto and Golubitsky [21].
Section 4 analyses the fully synchronous steady states u of the rate equations, a preliminary step towards studying local bifurcations from such states.We show that these are described by a cusp catastrophe surface in the sense of Poston and Stewart [29] and Zeeman [30].The geometry of this surface is used to show that the first local bifurcation as an input parameter I increases can be obtained from the corresponding bifurcation as u increases.This simplifies later calculations considerably and is central to the main theorem, a formal statement of Figure 2 that is stated in Section 7 and proved in Section 12.
Section 5 reviews basic existence results for local bifurcation (steady-state and Hopf) in symmetric ODEs.Section 6 summarises basic data for primary gaits: eigenvalues and eigenvectors of the connection matrix, symmetries of bifurcating branches for steady-state and Hopf bifurcation.We state the main theorem in Section 7.
Section 8 examines how the eigenvalues of the Jacobian evaluated at a fully synchronous equilibrium u relate to the eigenvalues of the connection matrix.It derives necessary and sufficient conditions for steady-state and Hopf bifurcations of given symmetry type to occur.These conditions are used to plot typical bifurcation curves for the primary gaits in Section 9.
Section 10 analyses which eigenvalue of the connection matrix is largest, as the connection strengths vary.The largest eigenvalue is closely associated with the first bifurcation to a state with the corresponding symmetry type.This principle is exploited in the proof of the main theorem, but first we discuss a symmetry group T that acts on the set of eigenvalues of the connection matrix in Section 11.
Here we prove that T also acts on the set of regions defining the first local bifurcation.The group T is not a symmetry group or parameter symmetry group of the rate equations in the sense of equivariance [18].We also describe some general representation-theoretic reasons that explain this structure.We then prove the main theorem in Section 12.
We briefly discuss degeneracies in the connection strengths in Section 13.Most of these correspond to transitions between different first bifurcations, that is, to mode interactions.We observe that in some cases disjoint sets of nodes decouple, a feature of rate models that is non-generic in general symmetric dynamical systems.Finally, Section 14 shows simulations of the rate equations giving the four primary gaits.

Gain Functions and Rate Models
In a rate model, the state of node j is a vector x j = (x E j , x H j ) ∈ R 2 , where x E j is an activity variable that that corresponds to the rate at which the node fires and x H j is a fatigue variable.Nodes are coupled through a gain function G.This has a sigmoid shape, and ideally we take G(x) = 0 for x < 0. In practice, it is common to assume slightly less: G(x) 1 for x < 0. A node j is then considered to be quiescent if it is in equilibrium at a small value of x E j .The original choice of G in Wilson [31] was a ratio of two quadratic terms.However, throughout this paper we follow later practice and assume a gain function of sigmoid form where a, b, c > 0. See Laing and Chow [32] and Shpiro et al. [33], who introduced such models in the context of binocular rivalry.With reasonable choices of a, b, c the value of G(x) is very small when x ≤ 0. A standard choice in the literature is a = 0.8, b = 7.2, c = 0.9.For simplicity, we use a = 1, b = 8, c = 1 in simulations.In our analysis, a, b, c are arbitrary positive constants.We expect similar results to apply for any reasonable sigmoidal gain function, but proving this could be difficult, and some fine detail would be likely to change.
The function G has an inverse F = G −1 given by Let G denote the derivative of G.There is an algebraic relation between G and F, which we exploit throughout this paper to obtain explicit formulas: We also note that each member of this family of gain functions is symmetric about its inflection point.Specifically, Indeed, We specify the rate model in terms of the connection matrix where n is the number of nodes and α ij ∈ R is the strength of the connection from node j to node i (positive for excitation, negative for inhibition).Let g > 0 be a parameter that determines the strength of reduction of the activity variable by the fatigue variable.Assume that nodes have some kind of input, with the input to node i being I i .The rate equations are: for 1 ≤ i ≤ n.Since j α ij x E j = (Ax E ) i we can write Equations ( 5) and ( 6) in the abbreviated form where , and the obvious convention is used to apply G to a vector.The model Equations ( 5) and ( 6) is very similar to, and was partly motivated by, the concept of a Wilson network, modelling the brain's decision-making process when presented with conflicting, or potentially conflicting, inputs.See Wilson [34].Wilson networks and generalisations have been used recently to model binocular rivalry and visual illusions [25,[35][36][37].The methods developed here also apply to some Wilson networks.

The CPG Network
Following Pinto and Golubitsky [21] we consider the 4-node network of Figure 1.Here, each leg corresponds to two nodes: one flexor and one extensor.
This network has symmetry group There are also parameter symmetries, of two kinds.First: any permutation of the four nodes produces the same network topology, but the connection strengths are permuted.For example (12) interchanges α and β, leaving γ fixed.Second: there is a parameter symmetry in which α, β, and γ are fixed but I changes, induced by the symmetry (4) of G.This leads to the following proposition, which is valid for rate equations on any homogeneous network, that is, one on which every node has the same number of input arrows of each type (connection strength).For example, in Figure 1 every node has one input arrow of type α, one of type β, and one of type γ.For rate models, it is convenient to employ a slightly more general condition: a network is homogeneous if and only if each row of the connection matrix A has the same sum r.That is, Proposition 1 Suppose that a network N is homogeneous, with connection matrix A defining a rate model, so every row of A has the same sum r.Then the rate Equation (11) are invariant under the parameter symmetry η for which and α, β, γ are fixed.
We employ the abbreviated notation of Equation (7), and scalars such as a are interpreted as [a, a, . . ., a] T .Let y = a − x, so that y E = a − x E and y H = a − x H . Then These are the same equations as Equation ( 7) with the stated parameter changes.
Here the only parameter that changes is I.The connection strengths α, β, γ are unchanged; the waveform of the solution is turned upside down and shifted by a constant.This parameter symmetry explains why the numerical bifurcation diagram in Curtu [27] Figure 2A, showing how the activity variable x E (there denoted u 1 ) depends on I, is symmetric under rotation through 180 • .
For the network of Figure 1, the connection matrix is This network is bidirectional (identical arrows point both ways) so A is symmetric.Therefore, A has real eigenvalues.We state them in Section 6.
To preserve symmetry, we take all inputs I j to be equal, so that I j = I for 1 ≤ j ≤ 4. We use I as bifurcation parameter and consider the connection strengths as subsidiary parameters; these are arbitrary, but held fixed as I varies.Explicit rate equations (written in a convenient order for later calculations, in which the activity variables and fatigue variables are collected together) are: Here the parameters are: Equation ( 11) have phase space The equations inherit the symmetry of the network, so they are equivariant under the action of (It is usual to use the inverse h −1 to define a left action, but in this case every element of G is equal to its inverse.)The subspaces R 4 E , R 4 H are G-invariant subspaces with isomorphic G-actions.Figure 3 shows the four phase patterns that occur generically in any Z 2 × Z 2 -symmetric system at Hopf bifurcation, as classified in Pinto and Golubitsky [21].We have abbreviated their "two-legged hop" and "two-legged jump" to "hop" and "jump".Which pattern occurs is determined by the representation of Z 2 × Z 2 on the critical (purely imaginary) eigenspace [17,18,38].See Section 5.This system can also undergo steady-state bifurcation, either symmetry-preserving or symmetry-breaking.This possibility is examined in Section 8.

Synchronous Steady States
As a preliminary step, we analyse fully synchronous steady states in the rate Equation (11) for the CPG network in Figure 1.That is, we determine conditions under which all nodes are in the same equilibrium state.Later we will consider local bifurcation from such states.
The main result generalises to any homogeneous network.When the network is homogeneous, the diagonal ∆ = {x : is an invariant subspace for the dynamics.We call states in ∆ fully synchronous.
Let the gain function be as in Equation ( 1), namely: The function G is monotonic strictly increasing and always positive.There is a unique inflection point at (c, a 2 ) and the slope at that point is G (c) = ab 4 .The curve is asymptotic to 0 as x → −∞ and to a as x → +∞.See Figure 4

We prove:
Theorem 1 Let N be a homogeneous network whose connection matrix A has all row-sums equal to r.
Then the fully synchronous equilibria u satisfy the equation If σ ≤ 4/ab there is a unique solution for all I.If σ > 4/ab there exist I L , I U depending on σ for which there are three solutions when I L < I < I U , and two solutions (one of multiplicity 2) when I = I L , I U .The set of solutions to Equation (12) forms a cusp catastrophe surface in (σ, I, u)-space.Its bifurcation set in (σ, I)-space has the equation where Consider a fully synchronous equilibrium at which all nodes have Substituting in the equations for ẋE i we obtain the single equation by Equation (8).Since σ = r − g we can rewrite this as Apply F to deduce Equation ( 12).The geometry of G, hence of F, implies that whenever σ ≤ c/2 there is a unique solution for all I.If σ > c/2 there exist I L , I U depending on σ for which there are three solutions when I L < I < I U , and two solutions (one of multiplicity 2) when I = I L , I U .See Figure 5   The derivative of F is The tangencies with lines of slope σ occur when Solve Equation ( 16) for x to get From Figure 5 we want solutions where 0 < u ± < a.
Solutions are real when σ < 0 or σ ≥ 4 ab .When σ < 0 both solutions are outside [0, 1].When σ ≥ 4 ab both solutions are in [0, a], which we require.Now substitute in Equation ( 15) to get In Figure 6a we plot the curves I L , I U against σ.Not surprisingly, we recognise the bifurcation curve of a cusp catastrophe surface, Zeeman [30], Poston and Stewart [29].Figure 6b shows this surface.There are three solutions to Equation ( 12), that is, three fully symmetric equilibria, when (σ, I) lies inside the cusp; a unique solution when (σ, I) lies outside the cusp; two solutions (one with multiplicity 2) when (σ, I) lies on the curves except at the cusp point; and one solution (with multiplicity 3) when (σ, I) lies at the cusp point.ab : this can be seen in constant-I cross-sections of Figure 6b, especially at the right of the figure.With quasistatic variation of the bifucation parameter I, the first steady state bifurcation occurs on crossing the upper branch I U of the cusp curve, which corresponds to u − .As we show later, there is always a Hopf bifurcation before that point.However, that does not rule out a synchrony-breaking steady-state bifurcation occurring before the Hopf bifurcation.We consider this possibility in Section 12 and show that it does not occur.
Remark 1 Curtu [26,27] study a rate model on a network with two identical nodes, and Theorem 1 applies to this network.However, in those papers there is a unique fully synchronous equilibrium for all I.The reason for this apparent discrepancy is that those papers assume an inhibitory connection, so σ < 0. Therefore, σ < c/2 since c > 0.
The Cusp Catastrophe.The geometry of the cusp catastrophe surface makes it relatively straightforward to understand how the first bifurcation from a synchronous equilibrium (as I varies quasistatically) depends on the connection strengths.See Zeeman [30] and Poston and Stewart [29] for further discussion in the context of the canonical cusp catastrophe.More precisely, the cusp geometry motivates the following procedure.
For each (α, β, γ) recall that σ = α + β + γ − g.Define a map Geometrically, the image of ξ σ is the constant-σ section of the cusp surface S. Define a projection This maps S to the parameter space (σ, I), and its singularities determine the fold lines and cusp point.
We can therefore use u as bifurcation parameter in place of I when considering the first bifurcation, provided we consider quasistatic changes as I increases from stable equilibrium.That is, I remains in equilibrium and varies continuously while this is possible.By Lemma 1 this first ceases to be possible only when σ > 4  ab and I reaches I U , the saddle-node bifurcation on the lower layer of the cusp surface S. We therefore have: Proposition 2 For any given α, β, γ the first local bifurcation as I increases quasistatically corresponds to the first local bifurcation as u increases.

By Lemma 1, if σ ≤ 4
ab then all local bifurcations occur in the same order using either I or u as bifurcation parameter.If σ > 4  ab then all local bifurcations occur in the same order when I ≤ I U , provided we consider quasistatic variation.But I U is itself a bifurcation point-here a saddle-node.
This proposition is decisive because the geometry of the bifurcation curves is much simpler in terms of u than it is using I.
Local bifurcations that are not the first can be converted from u to I using Equation (17) no matter what the value of σ may be.But when ξ σ is not invertible, the order gets mixed up.See the crossings in some of the Hopf bifurcation curves in Figures 7-14 below.

Review of Symmetric Steady-State and Hopf Bifurcation
Our main goal is to analyse the time-periodic gait patterns generated by the four-cell CPG of Figure 1.The symmetries of the primary gaits are spatio-temporal patterns, combining permutations of the cells with phase shifts.We assume that these patterns arise through symmetry-breaking Hopf bifurcation from a fully synchronous state, and use the Equivariant Hopf Theorem [17,18,38] to determine the corresponding bifurcation varieties: the subsets of parameter space at which bifurcations occur.The relevant parameters are I,which we take to be the main bifurcation parameter and the connection strengths α, β, γ.We assume that a, b, c, g, ε are fixed and positive.
However, steady-state bifurcations also occur in the rate equations.These correspond to various "stand" gaits, where the animal is motionless, which may not be of great interest.(However, there is a difference between passive standing and "getting ready" to initiate a periodic gait, which is one way to interpret these symmetry-breaking equilibria.)At any rate, the analysis has to take them into account.We discussed the bifurcation of synchronous (that is, symmetry-preserving) steady states in Section 5.
Here we review general theorems about symmetry-breaking local bifurcations.We state the Equivariant Branching Lemma for the existence of steady-state bifurcations with specified symmetry types.We also summarise the Equivariant Hopf Theorem and specify our notational conventions.This helps to avoid potential confusion in the later analysis, because several different conventions exist in the literature.We follow [17,18,38], with some minor changes.We apply these results in Section 8.
Let G be a group (technically, a compact Lie group, which includes all finite groups) acting linearly on R n .Denote the action of g where f : R n × R → R n is smooth.Here x represents the state of the system and λ is a bifurcation parameter.The equivariance condition is is a solution of a G-equivariant ODE and g ∈ G, then gx(t) is also a solution.
We do not consider such solutions to be significantly different from x(t), and by abuse of terminology call them conjugate solutions to x(t).The name arises because their symmetry groups are conjugate in G.
Assume that as λ varies near λ 0 there is a path of equilibria {(x(λ), λ) : λ ∈ R}.Local bifurcations (either steady-state or Hopf) occur at points (x 0 , λ 0 ) for which the Jacobian has at least one eigenvalue on the imaginary axis.Here the steady state loses stability.A zero eigenvalue corresponds to steady-state bifurcation, and a complex conjugate pair of purely imaginary eigenvalues corresponds to Hopf bifurcation to a time-periodic state.

Symmetric Steady-State Bifurcation
Steady-state bifurcation in the presence of symmetry is governed by the Equivariant Branching Lemma of Cicogna [39] and Vanderbauwhede [40,41].To state it we need some group-theoretic concepts.
The isotropy subgroup of x ∈ R n is and the isotropy subgroup of a solution x(t) is This space is invariant under the normaliser N G (Σ).Since Σ acts trivially, the effective action is that of If G is a compact Lie group, and in particular if G is finite, which is the case for gaits, then G × S 1 is a compact Lie group.Finite-dimensional real linear representations of compact Lie groups occur in three classes, distinguished by the space of commuting linear maps.By Schur's Lemma this is a division algebra over R, hence is isomorphic to one of R, C, and the quaternions H.If this algebra is R we call the representation absolutely irreducible.See Adams [42] 3.22.
Proposition 3.2 of Chapter XIII of [18] states: Proposition 3 Let f be a G-equivariant bifurcation problem (19), and suppose that along a branch of steady states (x(λ), λ) the Jacobian Df | (x 0 ,λ 0 ) has a zero eigenvalue, where x 0 = x(λ 0 ).Let E 0 be the corresponding real generalised eigenspace.Then generically E 0 is absolutely irreducible under the action of G.
A proof of the Equivariant Branching Lemma of Cicogna [39] and Vanderbauwhede [40,41] is given in Proposition 3.3 of Chapter XIII of [18], which states: Lemma 2 Consider a G-equivariant bifurcation problem (19) such that Df | (x 0 ,λ 0 ) has a zero eigenvalue for which E 0 is absolutely irreducible and the eigenvalue passes through zero with nonzero speed.Let Σ be an isotropy subgroup of G acting on E 0 for which Then there exists a unique branch of steady-state solutions to Equation (19), bifurcating from (x 0 , λ 0 ), with isotropy group Σ.
Note that it is the action of Σ on E 0 that is required to have a 1-dimensional fixed-point space, not the action on R n .

Symmetric Hopf Bifurcation
The Equivariant Hopf Theorem provides information about the spatio-temporal symmetries of bifurcating branches of periodic states.It also applies in circumstances where the classical Hopf bifurcation theorem does not, namely multiple eigenvalues caused by symmetry.
Let S 2π , S 1 2π be the loop spaces of continuous and once-differentiable 2π-periodic maps R → R n .These are Banach spaces.Define the circle group Then G × S 1 acts on S 2π and S 1 2π by A map v ∈ S 2π or S 1 2π is fixed by (g, θ) if and only if or equivalently So, transformation of v(s) by g has the same effect as a phase shift θ.The spatio-temporal isotropy group of v is the set of all (g, θ) ∈ G × S 1 for which Equation ( 20) is valid.
A representation W of G is said to be G-simple if either (1) W ∼ = V ⊕ V where V is absolutely irreducible.
(2) W is irreducible of type C or H.
It is clear that E iω is always G-invariant.In the generic case it also supports a natural S 1 -action.Define the restricted Jacobian We can define an action of S 1 on E iω by where we use |ω| to fix the orientation, an issue that becomes important for quadruped gaits to distinguish gaits like buck and walk from their reverses.By equivariance, J commutes with any g ∈ G.So G × S 1 acts on E iω by We say that the eigenvalues ±iω cross the imaginary axis with nonzero speed if there is a family of eigenvalues ρ(λ) ± σ(λ)i of D|f (x,λ(x)) for λ near λ 0 , such that The Equivariant Hopf Theorem of [17,38] is proved in Theorem 4.1 of Chapter XVI of [18], and states: Theorem 2 Consider a G-equivariant bifurcation problem (19) such that Df | (x 0 ,λ 0 ) has a pair of purely imaginary eigenvalues ±iω for which E iω is G-simple and the eigenvalues cross the imaginary axis with nonzero speed.Let Σ be an isotropy subgroup of G × S 1 acting on E iω for which Then there exists a unique branch of small-amplitude time-periodic solutions to Equation (19), bifurcating from (x 0 , λ 0 ), with period near 2π |ω| , having spatio-temporal symmetry group Σ in loop space.
Notice that we use the linearised action of J to define the circle group action on the critical eigenspace, for which dim R Fix(Σ) is computed, and infer the existence of a periodic solution of the full nonlinear equation with the same group Σ as its spatio-temporal symmetry group in loop space.So the symmetries are exact, not linearised approximations.
In the generic case where the critical eigenspace is G-simple, the spatio-temporal symmetry groups can be characterised as twisted subgroups is a group homomorphism.See Proposition 7.2 page 300 of [18].
The group H consists of all setwise symmetries of the trajectory through v.The kernel K = ker φ is the subgroup of all pointwise symmetries of that trajectory.The pair (H, K) provides a useful classification of the time-periodic solutions, and is characterised by the H/K Theorem of Buono and Golubitsky [16].However, there can exist different (that is, non-conjugate) solutions with the same pair (H, K).Examples in gaits are the walk and reverse walk.This distinction arises in particular for analogous models of quadruped and hexapod locomotor CPGs.
When G is abelian (as it is for Z 2 × Z 2 ) a G-simple critical eigenspace generically corresponds to a simple eigenvalue, so the classical Hopf bifurcation theorem applies.In particular, the bifurcating branch is locally unique.However, the Equivariant Hopf Theorem is a stronger result because the spatio-temporal isotropy subgroup provides additional information about the symmetries along the bifurcating branch.Similar remarks apply to steady-state bifurcation.

Data for Primary Gaits
Bifurcation depends on the eigenvalues of the Jacobian, and we will see that these are closely related to the eigenvalues of the connection matrix A. The eigenvectors and eigenvalues of A are given in Buono [14], Buono and Golubitsky [16], and Golubitsky et al. [19,20].Using T to denote column vectors, a basis adapted to the four irreducible representations of Z 2 × Z 2 on R 4 is: By equivariance or direct verification, these are the eigenvectors of A.
The corresponding eigenvalues are: where we use H, J, R, W to denote the four gait patterns hop, jump, run, walk respectively.The next proposition is trivial but important.Equalities between eigenvalues of A correspond to degeneracies in the bifurcations or the dynamics: see Section 13.
Proposition 5 The eigenvalues of A are distinct unless As well as A having equal eigenvalues, there are other conditions cause degeneracies in the bifurcations.We list them here: Definition 1 A triple of connection strengths (α, β, γ) is degenerate if any of the following conditions is satisfied for some patterns P, Q: The first of these conditions is equivalent to Equation (22).We consider the second in Section 13.
Pinto and Golubitsky [21] list the spatio-temporal isotropy subgroups of the primary gaits obtained by Hopf bifurcation in this Z 2 × Z 2 -symmetric network.We reproduce that table in the notation of this paper as Table 1, using Σ in place of their K to denote the spatial (that is, pointwise) isotropy subgroup.For all primary gaits, the setwise isotropy subgroup H is Z 2 × Z 2 .Since the group is abelian, this is also the normaliser of Fix(Σ).Recall that S 1 = R/2πZ so π is half the period.
Table 1.Spatio-temporal isotropy subgroups of the primary gaits.

Gait
Twist φ Σ = ker φ Fix(Σ) Type The same four isotropy groups determine the possible actions on a critical eigenspace E 0 for steady-state bifurcation, in the generic case when this space has dimension 1.We say that such a steady-state bifurcation is of type P if the isotropy subgroup is that of the primary gait P (= H, J, R, or W).Steady-state bifurcation of type H is symmetry-preserving; the other three are symmetry-breaking.

Main Theorem
We can now state the main theorem of this paper.It provides a complete characterisation of the regions of parameter space for which the first bifurcation is of a given type, as illustrated in Figure 2. Most of the remainder of the paper will provide a proof.
Assume that k < K, that is, ε < abg/4.Suppose that µ P = k, K and none of the conditions (22) holds.
Then the first local bifurcation (or its absence) is given by Table 2.
We develop the proof of this theorem in stages in the following sections.

Remark 2
The assumption that ε < abg/4 is natural because it is usual to assume that ε 1, leading to fast/slow dynamics.If ε > abg/4, no parameter values lead to Hopf bifurcation.

Local Bifurcation Analysis
We now derive conditions for local bifurcations of given symmetry type (that is, spatial isotropy subgroup) H, J, R, or W. Subject to technical genericity conditions, a local bifurcation occurs when the Jacobian J u at a fully synchronous equilibrium u has an eigenvalue on the imaginary axis.If it is zero, we get a steady-state bifurcation and the bifurcating branch consists of equilibria; if it is purely imaginary, we get a Hopf bifurcation and the branch consists of periodic states.

Eigenvalues of the Jacobian
It turns out that for any rate model with n nodes, there is a strong connection between the eigenvalues and eigenvectors of the Jacobian J u at a fully synchronous equilibrium u, and those of the connection matrix A. We therefore begin with a general theorem characterising the relation between eigenvalues and eigenvectors of J u and those of A.
Write the variables in the order x E 1 , . . ., x E n ; x H 1 , . . ., x H n .From Equations ( 5) and ( 6) the Jacobian at a fully synchronous equilibrium u takes the block form where Theorem 4 Let v be an eigenvector of A with eigenvalue µ.Then [(λ + 1)v, v] T is an eigenvector of J u with eigenvalue λ if and only if Moreover, every eigenvector is of this form.
Let w, v ∈ R n .Then [w, v] T is an eigenvector of J u with eigenvalue λ if and only if Then Equation ( 28) implies w = v + λv = (λ + 1)v so that Equation (27) becomes which simplifies to This suggests taking v to be an eigenvector of A, say Av = µv.We then get Therefore, if Equation ( 26) holds we obtain an eigenvector of J u .Solving Equation ( 26) yields: which we can rewrite as Finally, we must show that every eigenvector is of this form.If an arbitrary [w, v] T is an eigenvector of J u with eigenvalue λ, conditions Equations ( 27) and ( 28) are valid.By Equation ( 28), w = (λ + 1)v.Substituting for w in Equation ( 27) now proves that v is an eigenvector of A, unless λ + 1 = 0; that is, λ = −1.We claim that −1 is not an eigenvalue of J u .If it is, Equation (28) implies that w = 0, and Equation ( 27) reduces to gΓv = 0, implying that gΓ = 0.But the parameter g in Equation ( 5) is greater than 0, and Equation (25) implies that Γ > 0, a contradiction.
Remark 3 Generically, the quadratic Equation ( 26) has two distinct solutions.The solutions coincide if We discuss this condition in Theorem 10, and prove that at any local bifurcation it holds only when µ = K.This is a standard degeneracy condition for µ and corresponds to the transition from Hopf to steady-state for a given symmetry type P, where µ = µ P .Generically in the connection strengths it does not occur.
Corollary 1 For any α, β, γ, and sufficiently small I, the unique synchronous steady state u is linearly stable.
Let I → −∞.Then Γ → 0, so (26) becomes so all eigenvalues λ at u are negative for sufficiently small I.
In fact, the subsequent analysis of local bifurcations shows that u is linearly stable when I ≤ 0, indeed for larger I, provided we consider quasistatic variation of I from −∞ and ignore non-local bifurcations (such as jumps to the upper sheet of the cusp catastrophe surface).

Remark 4
The Hopf bifurcation curves in Figures 7-14 below, corresponding to hop and run, have segments with negative I.This does not contradict stability of the steady state for sufficiently small I because these segments are not relevant when I varies quasistatically.
Corollary 2 If (α, β, γ) is not degenerate in the sense of Definition 1, the Jacobian J u is semisimple (diagonalisable) over C.
Generically in connection-strengths, the eigenvalues µ P are distinct as P ranges over the four patterns.Generically in u, the above construction yields two distinct λ for each µ.The form [(λ + 1)v, v] T of the eigenvectors implies that distinct patterns P, Q lead to distinct eigenvectors µ P , µ Q , even if the associated eigenvalues are equal.This is the case even when there are degeneracies in the connection strengths (α, β, γ).
Corollary 2 shows that in the local bifurcation analysis for nondegenerate connection strengths, Theorem 4 determines all eigenvalues and eigenvectors of J u .

Conditions for Steady-State Bifurcation
We seek conditions for the corresponding λ eigenvalues to be purely imaginary.Here we work with a general network, making the assumption that µ is real.In particular, this is the case for all eigenvalues of A when the network is the one in Figure 1.
Theorem 5 Suppose that µ is real.Then the eigenvalue λ (evaluated at the synchronous equilibrium u) is zero if and only if the following equivalent conditions hold: Set λ = 0 in Equation ( 26) to obtain This leads to Equation (32), and also implies that Then use Equations ( 3) and ( 25):

Conditions for Hopf Bifurcation
We seek conditions for the corresponding eigenvalues λ to be purely imaginary.Again, we work with a general network, but now we assume that µ is real.A similar but more complicated result for complex µ can be proved, but we omit it here because it is not required.Theorem 6 Suppose that µ is real.Then the eigenvalue λ (evaluated at the synchronous equilibrium u) is purely imaginary if and only if By Equation (30) λ is purely imaginary if and only if: Equation ( 38) is equivalent to Substitute in Equation (39) to get bearing in mind that ε > 0. Therefore, Γ > ε/g, bearing in mind that g > 0. Substitute in Equation ( 40) to get In the simulations of Section 14, (1 + 1 ε )g is approximately 4.5.We can now derive a necessary condition for Hopf bifurcation from a synchronous equilibrium, for the primary gaits.
Theorem 7 Let µ P be the eigenvalue for the primary gait P, where P = H, J, R, or W. Then the necessary condition (38) for Hopf bifurcation to gait P is valid if and only if where Let u be the synchronous equilibrium.By Equation ( 12) (the first by definition, the second by Theorem 6).Therefore by Equation (3).That is, This has solutions and then Since A generically has four linearly independent eigenvectors (by equivariance under Z 2 × Z 2 ) we obtain 8 eigenvectors by this method, except when there is a solution of multiplicity 2 to Equation (26), which is non-generic as discussed in Theorem 10 below.

Plots of Hopf Bifurcation Curves
We now use Equations ( 41) and ( 42) to plot the Hopf bifurcation curves in (γ, I)-space for selected choices of α, β.More precisely, we plot curves given by the necessary condition (38) for Hopf bifurcation.For the moment we ignore the inequality (39), which leads to sufficient conditions by guaranteeing that the eigenvalue µ is imaginary.
These plots illustrate the main possibilities and motivate the subsequent analysis.Note that µ is real for all four gaits: it is the corresponding eigenvalue λ of J that is imaginary.In simulations, we choose to fix α, β and plot the curve of I against γ.Figures 7-14 show such plots for representative choices of α, β.
These plots share some general features.The curves for jump and walk extend towards γ = −∞, while those for hop, run, and the fold curve of the cusp extend towards γ = +∞.This feature is a consequence of the sign of γ in µ P and in σ.The hop and run curves fold over with a self-intersection.This is a consequence of the map ξ σ in Equation (18).The upper branch of the curve corresponds to the lower sheet of the cusp surface, and the lower branch of the curve corresponds to the upper sheet of the cusp surface: this follows since I U is defined by u − and I L is defined by u + in Equation (13).
In some of these plots, there is no Hopf bifurcation for some parameter values.For example, this is the case for Figures 7-10, where there is a gap (corresponding to some interval of values of γ) between the bifurcation curves for jump/walk and hop/run.However, other plots, such as Figures 11-14, do not have this gap.In Section 12 we show that the gaps correspond to connection strengths for which there is no local bifurcation from the branch of symmetric equilibria as I increases.(This is not yet obvious because there might be a symmetry-breaking steady-state bifurcation in this range.)

Dominant Eigenvalues
Local bifurcation from a fully synchronous state u is determined by the eigenvalues of the Jacobian J u .However, Equation (26) implies that the overall "skeleton" for the partition of (α, β, γ)-space, according to the type of the first bifurcation (relative to I, or equivalently to u), is determined by which of the four eigenvalues µ P of the connection matrix A is largest.We discuss this issue first.
Remark 5 If any of the degeneracy conditions ( 22) is valid, two or more eigenvalues become equal.
It is convenient to draw these regions by considering constant-γ sections of (α, β, γ)-space, since these are planes.For fixed γ the inequalities (44) let us plot the regions of (α, β) values corresponding to the four gaits.Figure 15 shows these regions for γ > 0, γ < 0. In three dimensions the structure of this partition is tetrahedral, and we can describe it in terms of the tetrahedron T with vertices at For any constant c the four planes µ P = c, for P = H, J, R, W, are parallel to the four faces of T .In fact, since the vertices of cT are (±c, ±c, ±c) with 1 or 3 minus signs, they contain the faces of cT .Define the face F P of T corresponding to gait P as follows: (1) F H is the face with vertices (−1, 1, 1), (1, −1, 1), (1, 1, −1) (2) F J is the face with vertices (−1, 1, 1), (1, −1, 1), (−1, −1, −1) Routine calculations, or an examination of Figure 15, show that the inequalities (44) partition (α, β, γ)-space into four connected regions: Theorem 9 For each gait pattern P, the region R P of (α, β, γ)-space in which µ P is larger than the other three eigenvalues of A is the positive cone In more detail, R P is the interior of an infinite triangular pyramid obtained by extending to infinity the pyramid with base F P and vertex the origin.These regions are related by rigid motions, corresponding to symmetries of T .Figure 2b illustrates the three-dimensional geometry involved.
A direct derivation in three dimensions is sketched in Remark 7 below.

Tetrahedral Structure
As Figure 15 suggests, the entire problem has an elegant tetrahedral structure, which arises from the expressions for the eigenvalues of A and the linear form of the argument of G in Equation (5).Indeed, there is a general representation-theoretic reason for it, which we explain at the end of this section.The tetrahedral symmetry involved is richer than the symmetry group Z 2 × Z 2 of the CPG network, and also richer than the parameter symmetries S 3 acting on α, β, γ.It simplifies the proof of Theorem 3.
The CPG N of Figure 1 can be represented as a tetrahedron, with the nodes corresponding to vertices and the arrows to edges.In fact, Figure 1 can be seen as the projection of a tetrahedron in which nodes 1 and 4 lie in a different horizontal plane from, nodes 2 and 3. Arrows of a given type (α, β, or γ) correspond to pairs of opposite edges-not having a common vertex.
The symmetry group S 4 of the tetrahedron acts on the four vertices by permuting them.We will see that this action of the tetrahedral group does not fully explain the tetrahedral structure of the first bifurcations, although it goes part way.Instead, we require a subtler (though related) action.We first discuss the above action.
The action of S 4 on nodes induces a permutation of the three pairs of opposite edges, that is, of the symbols (α, β, γ).This arises via the standard homomorphism ψ : S 4 → S 3 with kernel the Klein four-group V ∼ = Z 2 ×Z 2 .See for example Rotman [43] page 42.Here V comprises the parameter symmetries that fix the parameters, that is, the symmetry group Z 2 × Z 2 of N .
There is an induced action of S 4 , also with kernel V, and S 4 acts as rigid motions of This action preserves the face of the large tetrahedron that forms the base of the "hop" region, and permutes the other three regions.These parameter symmetries explain the S 3 symmetry of Figure 15, but not the tetrahedral symmetry.To see how this arises, we consider a different group acting on the set of linear forms ±α ± β ± γ that include the four eigenvalues of A. This is the wreath product Z 2 S 3 of order 48, Hall [44].Here the base group Z 2 × Z 2 × Z 2 changes the signs of α, β, γ, and S 3 permutes them.Geometrically, Z 2 S 3 is the symmetry group of rigid motions of the cube with vertices (±1, ±1, ±1).
Four of these forms are the eigenvalues of A, namely µ 1 = α + β + γ, µ 2 = −α + β − γ, µ 3 = −α − β + γ, µ 4 = α − β − γ.These are the forms for which the number of minus signs is even.The group Z 2 S 3 has a subgroup T of order 24 that changes signs in pairs.Geometrically, T is the symmetry group of rigid motions of the tetrahedron T with vertices (±1, ±1, ±1) in which the number of minus signs is 0 or 2.
By definition, T acts on the set of eigenvalues {µ P : P = H, J, R, W}.It is easy to prove that T ∼ = S 4 , and it permutes these eigenvalues faithfully.

Representation-Theoretic Generalities
We now describe a general representation-theoretic context that explains this structure, and can be used for other networks with an abelian symmetry group.
Suppose that G is any finite abelian group acting linearly on a complex vector space X.Then where the X j are irreducible representations.We further assume that the representations on each summand are non-isomorphic, so that these irreducibles are also the isotypic components of X, that is, the sums of all isomorphic irreducible subspaces.Theorem 3.5 of Chapter XII of [18] proves that the isotypic components are invariant under any matrix A that commutes with G.So A leaves each irreducible component X j invariant.(In particular, the connection matrix A has this property when G = Z 2 × Z 2 .)Let A denote the algebra of all matrices that commute with the G-action.
When G is abelian, the complex irreducibles are 1-dimensional over C. Therefore, any nonzero vector v j ∈ X j is a common eigenvector for the actions of all g ∈ G, and also for A. That is, with µ(A), χ(g) ∈ C. Suppose that a matrix B intertwines A. That is, whenever A ∈ A there exists Ã ∈ A for which Proposition 6 Suppose that B intertwines A, so there exists Ã satisfying Equation (46).Let v be any eigenvector of Ã with eigenvalue μ.Then Bv is an eigenvector of A with eigenvalue μ.
If B is invertible we can write the intertwining condition as so B acts on A by conjugation.Proposition 6 implies that the group B = B generated by all such B acts on the set of common eigenvectors of G, and also on the functions µ(A) expressing the eigenvalues of A.
For the four-cell network, G = Z 2 × Z 2 and We take B to consist of all 4 × 4 matrices that permute α, β, γ in the matrices A and multiply them by ±1.This is isomorphic to the wreath product Z 2 S 3 .However, minus the identity (−id ) fixes all eigenvalues.Therefore, the effective action is by The extension of −id splits, with complement the subgroup of signed permutations having an even number of minus signs.This is the group isomorphic to S 4 that acts on the eigenvalues and creates the tetrahedral symmetry in the space of connection strengths α, β, γ.

Remark 6
The tetrahedral group T is not a symmetry group of the rate equations in the sense of equivariance.More generally, it is not induced by a parameter symmetry.This follows since the first bifurcation to a hop gait is a saddle-node, whereas the first bifurcation to a jump, run, or walk gait can be shown to be a pitchfork, as suggested by the normalizer symmetry N G (Σ)/Σ ∼ = Z 2 .Thus T consists of symmetries of the bifurcation varieties, but not of the dynamics (even allowing changes of parameters).

Proof of the Main Theorem
We now complete the proof of Theorem 3. By Proposition 2 we can derive the geometry of bifurcation varieties for first bifurcation by employing u, the synchronous steady state, as the bifurcation parameter, as a proxy for I.
First, we set up some notation for the various bifurcation varieties and curves that arise.Let P be a gait pattern (P = H, J, R, W for hop, jump, run, walk respectively).We assume that a, b, c, g, ε > 0 are fixed and ε is sufficiently small; to be precise, that The remaining parameters define the parameter space R 3 = {(α, β, γ)}.We also consider R 4 = {(α, β, γ, u)}.By Equation (33) a necessary and sufficient condition for a steady-state bifurcation from u of type P is By Equation ( 43) a necessary condition for a Hopf bifurcation from u of type P is By Theorem 6 Hopf bifurcation occurs if and only if Equation ( 43) is valid and For fixed (α, β) Equations ( 48) and (49) define curves of u against γ. Figure 16 shows a typical set of these curves for fixed values of (α, β).Here u runs horizontally and γ vertically, so the first bifurcation, for a given value of γ, occurs as u crosses the first curve when u increases horizontally from 0 towards a at the level given by γ.Lemma 5 Suppose that ε < abg/4.Let P be a gait pattern. (1) since s + h ≤ a, with equality only when s = h = a/2.Therefore h < s when µ P < K, and h > s when µ P > K, as claimed.

Lemma 6
The group T permutes the sets S P , H P , and H P .The permutations act on the patterns {H, J, R, P} according to the S 4 -action on the corresponding eigenvalues of A.
For all four primary gaits P, the sets S P , H P , and H P are defined by equations and inequalities (48)-(50) in u that depend only on µ P , and do so in exactly the same manner for each gait type.

Corollary 3
The partition of (α, β, γ)-space into regions that determine the first local bifurcation as I increases from 0 is preserved by the action of T. It therefore has tetrahedral symmetry.
Lemma 7 Let P be a gait pattern.Let R P be as in Equation (51).Then for (α, β, γ) ∈ R P the first local bifurcation is: (1) None if µ P < k.
(3) Steady-state of type P if µ P > K.
Let P = H and specialise to the region R H . Equation (51) implies that µ Q < µ H when (α, β, γ) ∈ R H . Therefore, µ Q < k for all four patterns Q when (α, β, γ) ∈ R H .By Lemma 5 (2), bearing Lemma 4 in mind, there are no local bifurcations in R H when µ H < k.
If µ H > k there exists a Hopf bifurcation of type H.By Lemma 5 (1,3,4) this is the first bifurcation provided that k < µ P < K.
Finally, the first bifurcation is steady-state of type H when µ H > K by Lemma 5 (2,4).This proves the lemma for region R H . Apply the group T to deduce the other cases R P when P = J, R, W. Theorem 3 is an explicit statement of this lemma in terms of the connection strengths, so the proof is complete.

Plot of Gait Regions
We now complete the geometry of Figure 15 by incorporating the bifurcation varieties, both steady-state and Hopf, for the different symmetry types of pattern.Analysis of the figure shows that the geometry of the first bifurcation regions in (α, β, γ)-space is determined by two concentric and similar tetrahedra: For fixed γ, the regions of (α, β)-space for which the first bifurcation is a Hopf bifurcation are shown in Figure 17.These are constant-γ slices of T * \ T * .
Figure 18 shows how the four steady-state regions are related to Figure 17, using the same two-dimensional plots.The "no bifurcation" region is T * .We have T * T * if and only if k < K, that is, which is equivalent to the condition ε < abg 4 of Equation (47).In the simulations of Section 14, ε = 0.67 and abg 4 = 3.6 so this condition holds.In fact, k ∼ 0.835 and K ∼ 4.5.
Remark 7 An alternative approach is possible.We have derived the geometry of first bifurcations using two-dimensional sections in order to make it easier to find parameter values with given behaviour.There is, however, a more direct approach in three dimensions.For each pattern P, the linear form µ P defines a family of planes µ P = q for constant q.These planes are the orthogonal complements of four vectors: These vectors correspond to the vertices of T , so the planes µ P = q are parallel to the opposite faces.
Regions in which µ P < q are open half-spaces bounded by such planes.It is then easy to check that when q = K the planes are the faces of T * , extended to infinity, and when q = k they are the faces of T * , extended to infinity.For distinct patterns P and Q the planes µ P = µ Q determine the transitions between first bifurcation to P and first bifurcation to Q, as in Section 13 below.The rest of the geometry then follows.
The four vectors Equation (52) sum to zero in R 4 , and are analogous to trilinear coordinates {(x 1 , x 2 , x 3 ) : x 1 + x 2 + x 3 = 0} in R 3 , see Loney [45] chapter III.The analogy can be developed further by introducing connections of strength δ from each node to itself in Figure 1.Now the connection matrix A becomes A + δI 4 .The eigenvectors for the patterns H, J, R, W are mutually orthogonal (and become orthonormal if they are multiplied by 1  2 ).Now T * , T * are three-dimensional sections δ = 0 of the corresponding four-dimensional decomposition.

Degeneracies
We now discuss degenerate bifurcations, which occur when one or more of the degeneracy conditions (22) holds.By Corollary 2, degeneracies in the eigenvalues µ of A do not cause nilpotent Jordan blocks in J u even if they occur.However, the condition (53) below might signal such a Jordan block.
From the geometry of Figure 2, or from the detailed equations, degeneracies correspond to transitions between different local bifurcations.That is, they represent mode interactions.
Three main types of mode interaction occur along a generic path in (α, β, γ): (1) Transition from one primary Hopf mode to a different primary Hopf mode: change of gait.
(2) Transition from one primary steady-state mode to a different primary steady-state mode: change of equilibrium type.
(3) When µ P = K there is a transition between a primary Hopf mode of symmetry-type P and a steady-state mode of the same symmetry-type.Also, when µ P = k there is a transition between a primary Hopf mode of symmetry-type P and a fully synchronous steady-state mode.These transitions correspond to the onset or cessation of gait P.
We discuss transitions (1, 2) first.Rate equations have a special form, even for more general gain functions G than the family used here.This stems from the linearity of the coupling in the argument of G in Equation ( 5).This implies that if for given i, some linear combination of connection strengths α ij By Lemma 8, Γ = ε/g.At a steady-state bifurcation, In other words, the degeneracy in λ occurs when the first bifurcation is in transition from Hopf of type P to steady-state of type P. Such transitions occur when µ P = K.By Lemma 5(3) this occurs when the frequency ω of a Hopf bifurcation tends to zero.That is, we find an infinite-period or "blue sky" bifurcation.Typically, this is associated with a homoclinic orbit and arises via a Takens-Bogdanov bifurcation [46].The restriction of the Jacobian to the critical eigenspace at such a point is not semisimple but nilpotent, which is consistent with the linear degeneracy implied by Equation (53).
Theorem 4 characterises all eigenvectors of J u , and the trace of J u varies continuously with parameters, so coincidence of the values of λ implies that there exists a nontrivial 2 × 2 Jordan block for that value of λ.Diekman et al. [36] find a Takens-Bogdanov bifurcation in the "reduced network" for a rate model with two nodes and Z 2 symmetry.So we conjecture that the above transition may be associated with such a bifurcation.

Simulations
From Figure 18 it is now straightforward to select suitable connection strengths to obtain all four primary gaits in simulations of Equation (11).See Figures 19-22.Here the colours are node 1: red, node 2: blue, node 3: green, node 4: magenta.

Figure 1 .
Figure 1.Central Pattern Generator (CPG) network for bipedal gaits, with explicit connections corresponding to a rate model.

Figure 2 .
Figure 2. Tetrahedral partition of connection-space.(a) The two tetrahedra T * , T * ; (b) one of the four primary Hopf regions, in this case the hop gait.Base triangle left unshaded to show interior; (c) Extending the pyramid to define the corresponding steady-state regions, which extend to infinity in the direction indicated.
of medial connection β = strength of lateral connection ε = fast/slow dynamic timescale g = strength of reduction of activity variable by fatigue variable

Figure 3 .
Figure 3. Phase shift patterns for the four primary bipedal gaits.The numbers 0 and 1/2 indicate relative phases of the nodes.

Figure 4 .
Figure 4. Graph of gain function in the typical case a = 1, b = 8, c = 1.Here x runs horizontally and the vertical axis represents G(x).

Figure 6 .
Figure 6.(a) The curves I U , I L plotted against σ.Parameter values in the gain function are a = 1, b = 8, c = 1; (b) the cusp catastrophe surface in (σ, I, u)-space.
we fix a value of σ and use I as bifurcation parameter (one reasonable choice) then the bifurcation diagram is the usual S-shaped curve when σ >4

Figure 18 .
Figure 18.Regions of (α, β) values corresponding to first bifurcation to the four types of Hopf and four types of steady-state bifurcation, or to no bifurcation.(a) γ > 0; (b) γ < 0. Inner rectangle is empty when γ > k; outer rectangle is empty when γ > K.

Table 2 .
Necessary and sufficient conditions for each type of first bifurcation when no degeneracy condition holds.