Choose your preferred view mode

Please select whether you prefer to view the MDPI pages with a view tailored for mobile displays or to view the MDPI
pages in the normal scrollable desktop version. This selection will be stored into your cookies and used automatically
in next visits. You can also change the view style at any point from the main header when using the pages with your
mobile device.

Abstract

:
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 ℤ2 × ℤ2 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.

1. 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–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–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 ℤ4 × ℤ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 ℤ2 × ℤ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–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 Ju = 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
in a larger tetrahedron
, whose vertices are indicated on the figure: here K is a constant, defined in Theorem (23). The inner tetrahedron
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 ℤ2 × ℤ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 (ℤ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.

1.1. 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
that acts on the set of eigenvalues of the connection matrix in Section 11. Here we prove that
also acts on the set of regions defining the first local bifurcation. The group
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.

2. Gain Functions and Rate Models

In a rate model, the state of node j is a vector
xj=(xjE,xjH)∈ℝ2, where
xjE is an activity variable that that corresponds to the rate at which the node fires and
xjH is a fatigue variable. Nodes are coupled through a gain function. This has a sigmoid shape, and ideally we take
(x) = 0 for x < 0. In practice, it is common to assume slightly less:
(x) ≪ 1 for x < 0. A node j is then considered to be quiescent if it is in equilibrium at a small value of
xjE .

The original choice of
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

G(x)=a1+e−b(x−c)

(1)

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
(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
has an inverse
=
−1 given by

F(u)=c+1blog(ua−u)

(2)

Let
denote the derivative of
. There is an algebraic relation between
and
, which we exploit throughout this paper to obtain explicit formulas:

G′(F(u))=bau(a−u)

(3)

We also note that each member of this family of gain functions is symmetric about its inflection point. Specifically,

G(2c−x)=a−G(x)

(4)

Indeed,

G(x)=a2+a2tanh(b(x−c)2)

and tanh is odd.

We specify the rate model in terms of the connection matrix

A=(αij)1≤i,j≤n

where n is the number of nodes and αij ∈ ℝ 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 Ii. The rate equations are:

ɛx˙iE=−xiE+G(Ii−gxiH+∑jαijxjE)

(5)

x˙iH=xiE−xiH

(6)

for 1 ≤ i ≤ n. Since
∑jαijxjE=(AxE)i we can write Equations (5) and (6) in the abbreviated form

ɛx˙E=−xE+G(−gxH+AxE+I)x˙H=xE−xH

(7)

where
xE=[xiE]T ,
xH=[xiH]T ,
I=[Ii]T and the obvious convention is used to apply
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–37]. The methods developed here also apply to some Wilson networks.

3. 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 G = ℤ2 〈ρ〉 × ℤ2 〈τ〉 generated by

ρ=(12)(34)τ=(14)(23)

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
. 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,

∑jαij=r1≤i≤n

(8)

Proposition 1Suppose that a networkis homogeneous, with connection matrix A defining a rate model, so every row of A has the same sum r. Then the rateEquation (11)are invariant under the parameter symmetry η for which

η(x,I)=(a−x,2c+ga−ra−I)

(9)

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 yE = a − xE and yH = a − xH. 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 xE (there denoted u1) depends on I, is symmetric under rotation through 180°.

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 Ij to be equal, so that Ij = 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:

The equations inherit the symmetry of the network, so they are equivariant under the action of G = ℤ2 × ℤ2 defined by

(hx)j=xh(j)(h∈G)

(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
ℝE4 ,
ℝH4 are G-invariant subspaces with isomorphic G-actions.

Figure 3 shows the four phase patterns that occur generically in any ℤ2 × ℤ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 ℤ2 × ℤ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.

4. 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:xj=xj1≤i,j≤n}

is an invariant subspace for the dynamics. We call states in Δ fully synchronous.

The function
is monotonic strictly increasing and always positive. There is a unique inflection point at
(c,a2) and the slope at that point is
G′(c)=ab4 . The curve is asymptotic to 0 as x → −∞ and to a as x → +∞. See Figure 4 for the case a = 1, b = 8, c = 1.

We prove:

Theorem 1Letbe a homogeneous network whose connection matrix A has all row-sums equal to r. Let

σ=r−g

Then the fully synchronous equilibria u satisfy the equation

I+σu=F(u)

(12)

If σ ≤ 4/ab there is a unique solution for all I. If σ > 4/ab there exist IL, IU depending on σ for which there are three solutions when IL < I < Iu, and two solutions (one of multiplicity 2) when I = IL, IU.

The set of solutions toEquation (12)forms a cusp catastrophe surface in (σ, I, u)-space. Its bifurcation set in (σ, I)-space has the equation

IU=c+1blog(u−a−u−)−σu−IL=c+1blog(u+a−u+)−σu+

(13)

where

u±=a±a(a−4bσ)2

Consider a fully synchronous equilibrium at which all nodes have

xiE=uxiH=υ1≤i≤n

Equation (6) implies that
xiE−xiH=0 , so &upsilon; = u. Substituting in the equations for
x˙iE we obtain the single equation

The geometry of
, hence of
, implies that whenever σ ≤ c/2 there is a unique solution for all I. If σ > c/2 there exist IL, IU depending on σ for which there are three solutions when IL < I < IU , and two solutions (one of multiplicity 2) when I = IL, IU. See Figure 5 when a = 1, b = 8, c = 1.

In Figure 6a we plot the curves IL, IU 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.

If 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
σ>4ab : 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 IU 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

ξσ:ℝ→ℝξσ(u)=F(u)−σu

(18)

Geometrically, the image of ξσ is the constant-σ section of the cusp surface
. Define a projection

π:S→ℝ2π(σ,I,u)=(σ,I)

This maps
to the parameter space (σ, I), and its singularities determine the fold lines and cusp point. The bifurcation varieties for local bifurcation are the images under π ○ ξσ of the corresponding varieties in (σ, u)-space. Therefore, saddle-node (fold) bifurcations in the fully synchronous state u, which correspond to tangencies in Figure 5, are given by Equation (15), so that

I=ξσ(u)

along with Equation (16). The next lemma shows that we can invert ξσ on suitable domains to rewrite this as

Ifσ>4aband ξσ : (0,u−) → (−∞,IU) is a monotonic strictly increasing homeomorphism. It is a diffeomorphism when u < u−.

The proof is a restatement of the analysis in Theorem 1 and Equation (17).

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
σ>4ab and I reaches IU, the saddle-node bifurcation on the lower layer of the cusp surface
. We therefore have:

Proposition 2For any given α, β, γ the first local bifurcation as I increases quasistatically corresponds to the first local bifurcation as u increases.

By Lemma 1, if
σ≤4ab then all local bifurcations occur in the same order using either I or u as bifurcation parameter. If
σ>4ab then all local bifurcations occur in the same order when I ≤ IU , provided we consider quasistatic variation. But IU 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, 8, 9, 10, 11, 12, 13 and 14 below.

5. 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 ℝn. Denote the action of g ∈ G on x ∈ ℝn by gx. Consider a G-equivariant ODE

x˙=f(x,λ)

(19)

where f : ℝn × ℝ → ℝn is smooth. Here x represents the state of the system and λ is a bifurcation parameter. The equivariance condition is

f(gx,λ)=gf(x,λ)∀g∈G

It follows that if x(t) 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(λ), λ) : λ ∈ ℝ}. Local bifurcations (either steady-state or Hopf; occur at points (x0, λ0) for which the Jacobian

Df|(x0,λ0)

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.

5.1. 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 ∈isin;ℝn is

∑x={g∈G:gx=x}

and the isotropy subgroup of a solution x(t) is

∑x(t)={g∈G:gx(t)=x(t)∀t∈ℝ}

Conversely, if Σ is a subgroup of G then the fixed-point subspace of Σ is

Fix(∑)={x∈ℝn:gx=x∀g∈∑}

This space is invariant under the normaliser Ng(Σ). Since Σ acts trivially, the effective action is that of Ng(Σ)/Σ.

If G is a compact Lie group, and in particular if G is finite, which is the case for gaits, then G ×
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 ℝ, hence is isomorphic to one of ℝ, ℂ, and the quaternions (x0210D;. If this algebra is (x0211D; we call the representation absolutely irreducible. See Adams [42] 3.22.

Proposition 3Let f be a G-equivariant bifurcation problem(19), and suppose that along a branch of steady states (x(λ),λ) the Jacobian Df|(x0,λ0)has a zero eigenvalue, where x0 = x(λ0). Let E0be the corresponding real generalised eigenspace. Then generically E0is 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 2Consider a G-equivariant bifurcation problem(19)such that Df|(x0,x0)has a zero eigenvalue for which E0is absolutely irreducible and the eigenvalue passes through zero with nonzero speed. Let Σ be an isotropy subgroup of G acting on E0for which

dimℝFix(∑)=1

Then there exists a unique branch of steady-state solutions toEquation (19), bifurcating from (x0, λ0), with isotropy group Σ.

Note that it is the action of Σ on E0 that is required to have a 1-dimensional fixed-point space, not the action on ℝn.

5.2. 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
2π,
S2π1 be the loop spaces of continuous and once-differentiable 2π-periodic maps ℝ → ℝn. These are Banach spaces. Define the circle group

S1=ℝ/2πℤ

Then G ×
1 acts on
2π and
S2π1 by

(g,θ)υ(s)=gυ(s+θ)g∈G,θ∈S1

A map &upsilon; ∈
2π or
S2π1 is fixed by (g, θ) if and only if

gυ(s+θ)=υ(s)

(20)

or equivalently

gυ(s)=υ(s−θ)

So, transformation of &upsilon;(s) by g has the same effect as a phase shift θ. The spatio-temporal isotropy group of &upsilon; is the set of all (g, θ) ∈ G ×1 for which Equation (20) is valid. A representation W of G is said to be G-simple if either

It is clear that Eiω is always G-invariant. In the generic case it also supports a natural
1-action. Define the restricted Jacobian

J=Df|(x0,λ0)|Eiω

We can define an action of
1 on Eiω by

θυ=exp[θ|ω|J]υυ∈Eiω

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 ×1 acts on Eiω by

(γ,θ)υ=γexp[θ|ω|J]υ

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

ρ(0)=0σ(0)=ωσ′(0)≠0

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|(x0,λ0) has apair of purely imaginary eigenvalues ±iω for which Eiω is G-simple and the eigenvalues cross the imaginary axis with nonzero speed. Let Σ be an isotropy subgroup of G×
1 acting on Eiω for which

dimℝfix(∑)=2

Then there exists a unique branch of small-amplitude time-periodic solutions to Equation (19), bifurcating from (x0, λ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ℝ 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

The group H consists of all setwise symmetries of the trajectory through &upsilon;. 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 ℤ2 × ℤ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.

6. 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 ℤ2 × ℤ2 on ℝ4 is:

By equivariance or direct verification, these are the eigenvectors of A. The corresponding eigenvalues are:

1:hop.μH=α+β+γ2:jump.μJ=−α+β−γ3:run.μR=−α−β+γ4:walk.μW=α−β−γ

(21)

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

α=±βα=±γβ=±γ

(22)

As well as A having equal eigenvalues, there are other conditions that 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:

μP=μQP≠QμP=(1+1ɛ)g

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 ℤ2 × ℤ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 ℤ2 × ℤ2. Since the group is abelian, this is also the normaliser of Fix(Σ). Recall that
= ℝ/2πℤ so π is half the period.

The same four isotropy groups determine the possible actions on a critical eigenspace E0 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.

7. 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.

Theorem 3 Let

k=4(1+ɛ)abK=(1+1ɛ)g

(23)

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.

8. 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 Ju 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.

8.1. 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 Ju 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 Ju and those of A.

Write the variables in the order
x1E ,…,
xnE ;
x1H ,…,
xnH . From Equations (5) and (6) the Jacobian at a fully synchronous equilibrium u takes the block form

Ju=[1ɛ(−In+ΓA)−gɛΓInIn−In]

(24)

where

Γ=G′(I+σu)

(25)

Theorem 4 Let &upsilon; be an eigenvector of A with eigenvalue μ. Then [(λ + 1)&upsilon;, &upsilon;]T is an eigenvector of Ju with eigenvalue λ if and only if

ɛλ2+(1+ɛ−Γμ)λ+(1−Γμ+Γg)=0

(26)

Moreover, every eigenvector is of this form.

Let w, &upsilon; ∈ ℝn. Then [w, &upsilon;]T is an eigenvector of Ju with eigenvalue λ if and only if

Finally, we must show that every eigenvector is of this form. If an arbitrary [w, υ]T is an eigenvector of Ju with eigenvalue λ, conditions Equations (27) and (28) are valid. By Equation (28)w = (λ + 1)υ. Substituting for w in Equation (27) now proves that υ is an eigenvector of A, unless λ + 1 = 0; that is, λ = −1. We claim that −1 is not an eigenvalue of Ju. If it is, Equation (28) implies that w = 0, and Equation (27) reduces to gΓυ = 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

(1−ɛ−Γμ)2=4ɛΓg

(31)

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.

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, 8, 9, 10, 11, 12, 13 and 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 Ju is semisimple (diagonalisable) over ℂ.

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)υ, υ]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 Ju.

8.2. 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:

8.3. 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

Since A generically has four linearly independent eigenvectors (by equivariance under ℤ2 × ℤ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.

9. 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, 8, 9, 10, 11, 12, 13 and 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 IU is defined by u− and IL 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, 8, 9 and 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, 12, 13 and 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.)

10. Dominant Eigenvalues

Local bifurcation from a fully synchronous state u is determined by the eigenvalues of the Jacobian Ju. 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.

Theorem 8 Assume that the degeneracy conditions (22) do not hold. Then the largest eigenvalue of A is:

μH⇔α+γ>0α+β>0γ+β>0μJ⇔α+γ<0β>γβ>αμR⇔α+β<0γ>βγ>αμW⇔γ+β<0α>γα>β

(44)

This follows directly from Equation (21). For example, μH> μJ if and only if α+β+γ > −α+β−γ, that is, α + γ > 0. Five similar comparisons establish the stated conditions.

Remark 5 If any of the degeneracy conditions (22) is valid, two or more eigenvalues become equal.

The conditions (44) partition ℝ3 = {(α,β,γ)} into four open regions. Their closure is ℝ3 and the complement of their union comprises values of (α,β,γ) satisfying at least one of the degeneracy conditions (22).

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
with vertices at

(α,β,γ)=(−1,1,1)(1,−1,1)(1,1,−1)(−1,−1,−1)

For any constant c the four planes μP = c, for P = H, J, R, W, are parallel to the four faces of
. In fact, since the vertices of c are (±c, ±c, ±c) with 1 or 3 minus signs, they contain the faces of c. Define the face FP of
corresponding to gait P as follows:

(1)

FH is the face with vertices (−1, 1, 1), (1, −1, 1), (1, 1, −1)

(2)

FJ is the face with vertices (−1, 1, 1), (1, −1, 1), (−1, −1, −1)

(3)

FR is the face with vertices (−1, 1, 1), (1, 1, −1), (−1, −1, −1)

(4)

FW 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 RP of (α, β, γ) -space in which μP is larger than the other three eigenvalues of A is the positive cone

RΡ={rx:r>0,x∈int(FΡ)}

(45)

In more detail, RP is the interior of an infinite triangular pyramid obtained by extending to infinity the pyramid with base FP and vertex the origin. These regions are related by rigid motions, corresponding to symmetries of
. Figure 2b illustrates the three-dimensional geometry involved.

A direct derivation in three dimensions is sketched in Remark 7 below.

11. 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
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 ℤ2 × ℤ2 of the CPG network, and also richer than the parameter symmetries
3 acting on α, β, γ. It simplifies the proof of Theorem 3.

The CPG
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
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
4 on nodes induces a permutation of the three pairs of opposite edges, that is, of the symbols (α, β, γ). This arises via the standard homomorphism

ψ:𝕊4→𝕊3

with kernel the Klein four-group
≅ ℤ2 × ℤ2. See for example Rotman [43] page 42. Here
comprises the parameter symmetries that fix the parameters, that is, the symmetry group ℤ2 × ℤ2 of
.

There is an induced action of
4, also with kernel
, and
4 acts as rigid motions of

ℝ3≡{(α,β,γ):α,β,γ∈ℝ}

Since
4 acts by parameter symmetries, the partition of ℝ3 into regions in Figures 15, 16 and 17 are preserved. 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
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
of order 48, Hall [44]. Here the base group ℤ2 × ℤ2 × ℤ2 changes the signs of α, β, γ, and
3 permutes them. Geometrically,
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
has a subgroup
of order 24 that changes signs in pairs. Geometrically,
is the symmetry group of rigid motions of the tetrahedron
with vertices (±1, ±1, ±1) in which the number of minus signs is 0 or 2.

By definition,
acts on the set of eigenvalues {μP : P = H, J, R, W}. It is easy to prove that
≅
4, and it permutes these eigenvalues faithfully.

11.1. 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

X=⊕Xj

where the Xj 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 Xj invariant. (In particular, the connection matrix A has this property when G = ℤ2 × ℤ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 ℂ Therefore, any nonzero vector υj ∈ Xj is a common eigenvector for the actions of all g ∈ G, and also for A. That is,

Avj=μ(A)vjgvj=χ(g)vj∀j

with μ(A), χ(g) ∈ ℂ Suppose that a matrix B intertwines
. That is, whenever A ∈
there exists Ã ∈
for which

AB=BA˜

(46)

Proposition 6 Suppose that B intertwines
so there exists Ã satisfying Equation (46). Let υ be any eigenvector of Ã with eigenvalue μ̃. Then Bυ is an eigenvector of A with eigenvalue μ̃.

A(Bv)=(AB)v=(BA˜)v=B(A˜v)=B(μ˜v)=μ˜(Bv)

If B is invertible we can write the intertwining condition as

B−1AB=A˜

so B acts on
by conjugation. Proposition 6 implies that the group
= 〈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 = ℤ2 × ℤ2 and

A={A:α,β,γ∈ℝ}

We take
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
. 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
4 that acts on the eigenvalues and creates the tetrahedral symmetry in the space of connection strengths α, β, γ.

Remark 6 The tetrahedral group
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 NG(Σ)/Σ ≅ ℤ2. Thus
consists of symmetries of the bifurcation varieties, but not of the dynamics (even allowing changes of parameters).

12. 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

By Equation (33) a necessary and sufficient condition for a steady-state bifurcation from u of type P is

abu(a−u)=μΡ−g

(48)

By Equation (43) a necessary condition for a Hopf bifurcation from u of type P is

abu(a−u)=μΡ1+ɛ

(49)

By Theorem 6 Hopf bifurcation occurs if and only if Equation (43) is valid and

μΡ<(1+1ɛ)g

(50)

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 γ.

Which of these bifurcations occurs first (or none) is determined by the relative positions of the eight curves. Since they have explicit, simple equations, it is straightforward to derive necessary and sufficient conditions. The tetrahedral symmetry
of the eigenvalues of A simplifies this derivation and helps to explain the results, as we see below. First, we set up some notation.

Then u is a point of steady-state bifurcation for parameters (α, β, γ) if and only if
SΡαβγ . The other such point is a − u. Similarly, u is a point of Hopf bifurcation for parameters (α, β, γ) if and only if
ℋ¯Ραβγ. The other such point is a − u. We restrict u to the interval (0, a/2] because the first bifurcation as u increases must occur on this interval.

By Equations (48) and (50), each of the sets
SΡαβγ ,
ℋΡαβγ , and
ℋ¯Ραβγ is either a singleton or empty. Define
sΡαβγ ,
hΡαβγ by

SΡαβγ={sΡαβγ}ℋΡαβγ={hΡαβγ}

when the relevant set is non-empty, and leave them undefined if it is empty.

For given (α, β, γ), the ordering of the bifurcations as u increases from 0 is determined by which of the values
sΡαβγ ,
hΡαβγ (among those that exist) is smallest. If none exists, there is no bifurcation.

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 RH. Apply the group
to deduce the other cases RP 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.

12.1. 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
\
.

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
. We have
⊊
if and only if k < K, that is,

4(1+ɛ)ab<g(1+1ɛ)

which is equivalent to the condition
ɛ<abg4 of Equation (47). In the simulations of Section 14, ε = 0.67 and
abg4=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:

μH=[1,1,1]T⊥μJ=[−1,1,−1]T⊥μR=[−1,−1,1]T⊥μW=[1,−1,−1]T⊥

(52)

These vectors correspond to the vertices of
, 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
, extended to infinity, and when q = k they are the faces of
, 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 ℝ4, and are analogous to trilinear coordinates {(x1, x2, x3) : x1 + x2 + x3 = 0} in ℝ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 + δI4. The eigenvectors for the patterns H, J, R, W are mutually orthogonal (and become orthonormal if they are multiplied by
12). Now
,
are three-dimensional sections δ = 0 of the corresponding four-dimensional decomposition.

13. 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 Ju 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
than the family used here. This stems from the linearity of the coupling in the argument of
in Equation (5). This implies that if for given i, some linear combination of connection strengths αij is zero, and if moreover the corresponding cells j are in synchrony, then their combined input into cell i is zero. That is, cell i is decoupled from the cells j. We therefore have:

(1)

If α = − β and the gait is of types H or P, nodes 1, 2 decouple from 3, 4.

(2)

If α = −γ and the gait is of types H or J, nodes 1, 3 decouple from 2, 4.

(3)

If β = −γ and the gait is of types H or W, nodes 1, 4 decouple from 2, 3.

The decoupling implies that the dynamics depends on the initial conditions. In particular, phase shifts between the corresponding pairs of nodes can be arbitrary. This behaviour would be non-generic in a general dynamical system, even one with ℤ2 × ℤ2 symmetry. It follows that the behaviour of mixed modes (secondary bifurcations) in rate models might not resemble that of generic mixed modes in dynamical systems with given symmetry.

Cancellations involving more nodes can also occur, for example when α + β + 2γ = 0. Decoupling now requires at least three nodes to be synchronous, which generically occurs for type H only, and all four nodes being synchronous.

This type of decoupling is a general feature of all rate models.

There can be a similar “approximate decoupling” near a point of Hopf bifurcation when nodes are in antisynchrony. Here the linearised eigenfunctions are very close to p ± qsinωt for constants p, q and frequency ω. Therefore, their sum is close to 2p, a constant. For the exact waveform, only higher Fourier components will occur. For example, suppose α = β and the gait is a jump. Then nodes 2, 4 are in antisychrony with nodes 1, 3. So each pair of nodes transmits an almost constant signal to the other pair.

A transition of type (3) is associated with the other potential source of degeneracy in the associated eigenvalues λ of Ju, namely Equation (31). For convenience we repeat the equation here: the two solutions for λ coincide for given μ when

(1−ɛ−Γμ)2=4ɛΓg

(53)

Lemma 8If Equation (53) holds at a point of local bifurcation then

Γ=ɛg

Suppose the bifurcation is steady-state. By Equation (34) we then have

Theorem 10If Equation (53) holds at a point of local bifurcation then μ = K.

By Lemma 8, Γ = ε/g. At a steady-state bifurcation,

μ=g+1Γ=g+gɛ=(1+1ɛ)g=K

At a Hopf bifurcation

μ=1+ɛΓ=1+ɛɛg+K

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 Ju, and the trace of Ju 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 ℤ2 symmetry. So we conjecture that the above transition may be associated with such a bifurcation.

14. 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, 20, 21 and 22. Here the colours are node 1: red, node 2: blue, node 3: green, node 4: magenta.

Acknowledgments

I am grateful to Marty Golubitsky for introducing me to rate equations and for many helpful conversations.

Figure 2.
Tetrahedral partition of connection-space. (a) The two tetrahedra (img), (img); (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.

Figure 2.
Tetrahedral partition of connection-space. (a) The two tetrahedra (img), (img); (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.

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

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

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.

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.