key: cord-0431860-chi7y297
authors: Postlethwaite, Claire M; Rucklidge, Alastair M
title: Stability of cycling behaviour near a heteroclinic network model of Rock-Paper-Scissors-Lizard-Spock
date: 2020-10-20
journal: nan
DOI: nan
sha: 00c1a6d3c08093df105aa91754fe52807a63cfbb
doc_id: 431860
cord_uid: chi7y297

The well-known game of Rock--Paper--Scissors can be used as a simple model of competition between three species. When modelled in continuous time using differential equations, the resulting system contains a heteroclinic cycle between the three equilibrium solutions representing the existence of only a single species. The game can be extended in a symmetric fashion by the addition of two further strategies (`Lizard' and `Spock'): now each strategy is dominant over two of the remaining four strategies, and is dominated by the remaining two. The differential equation model contains a set of coupled heteroclinic cycles forming a heteroclinic network. In this paper we carefully consider the dynamics near this heteroclinic network. We are able to identify regions of parameter space in which arbitrarily long periodic sequences of visits are made to the neighbourhoods of the equilibria, which form a complicated pattern in parameter space.

The well-known game of Rock-Paper-Scissors can be used as a simple model of competition between three species. This system has been studied extensively, in many different contexts, such as evolutionary game theory and biology [1] [2] [3] , stochastic models [2, 4] and spatially dependent models [4] [5] [6] [7] . For reviews see [8, 9] . When modelled in continuous time using ordinary differential equations, the resulting system contains a heteroclinic cycle between the three equilibrium solutions representing the existence of only a single species [1] . This same heteroclinic cycle has also been studied in the context of fluid dynamics [10] and equivariant dynamical systems [11, 12] .

The addition of two further strategies -where each of the now five strategies is dominant over two of the remaining four strategies, and is dominated by the remaining two -extends the Rock-Paper-Scissors system into a network of possible states. The game of Rock-Paper-Scissors-Lizard-Spock was popularised by the TV show 'The Big Bang Theory' in 2012 [13] (and credited to Sam Kass [14] ), but equivalent networks have been around for much longer: the Wuxing cycle of five phases from Chinese philosophy has been in existence since the second or first century BCE [15] . In modern literature, the earliest reference we can find to this network of interactions in an equivariant dynamical systems context is by Field and Richardson in 1992 [16] .

In this paper we present the first comprehensive description of the dynamics of this system when described by a system of ordinary differential equations (ODEs), namely the equations

where x j ∈ R, X = x 1 +x 2 +x 3 +x 4 +x 5 , and c A , c B , e A , e B > 0 are parameters. These ODEs contain a heteroclinic network between equilibrium solutions, shown schematically in figure 1. The stability of the entire network has been studied by Podvigina [17] , and Afraimovich and colleagues [18] , but they do not discuss the dynamics that occur as trajectories approach the network. Some results in a stochastic setting can be found in [19, 20] . The vector field generated by equations (1) (which have only linear and quadratic terms), is often thought of as being of 'Lokta-Volterra' (or 'May-Leonard') type. It is topologically equivalent, in the positive orthant, to one with Z 5 2 Z 5 symmetry which can be generated by the coordinate transformation x j → x 2 j . Vector fields of this latter type, with symmetry group Z k 2 or Z k 2 Z k (in both cases acting on R k ) can contain heteroclinic cycles, and some have been studied in the context of equivariant bifurcation theory -for example see [12] for k = 3, and [21] for k = 4. In all these cases the invariant sphere theorem applies [22] , and so all trajectories are attracted to an invariant (k − 1)-sphere. The results we present in this paper apply only to the dynamics of (1) in the positive orthant, and as such, are independent of whether we consider these equations, or the ones with the coordinates transformed to squared variables.

For the equivariant equations (i.e. those with linear and cubic terms), the existence and stability of heteroclinic cycles was studied in detail for k = 4 by Field and Swift [21] . The results of Field and Swift are complemented by studies of heteroclinic networks in four dimensions by Brannath [23] and Kirk and Silber [24] . Krupa and Melbourne generalised the stability results in the k = 3 and k = 4 cases to include cases without the Z k permutation symmetry [25] . In addition to the so-called edge cycles between equilibria with just one species present, the results of Field and Swift [21] show that face cycles between equilibria with two non-zero components can occur (see also [26] ). In the k = 5 case, Field and Richardson [16] show that these face cycles can co-exists with edge cycles, and there can also exist cycles between equilibria with three non-zero components. In this paper, we do not consider the face-cycles, because our choice of parameters (namely, that c A , c B , e A , e B > 0) does not allow it, but we do consider the cycle between equilibria with three non-zero components (3-face cycles in the terminology of [26] ).

Methods for determining the stability properties of an isolated robust heteroclinic cycle between equilibria are well-established [25, [27] [28] [29] [30] [31] [32] [33] [34] , and their implementation is generally straightforward, at least in principle, because there is only a single route around the cycle. Heteroclinic cycles which are proper subsets of a heteroclinic network cannot be asymptotically stable, because there must be some points on an unstable manifold of at least one of the equilibria in the cycle which are attracted to a different part of the network. This lead to the introduction of weaker notions of stability, firstly essential asymptotic stability, introduced by Ian Melbourne in 1991 [32] , and later fragmentary asymptotic stability, introduced by Olga Podvigina in 2012 [27] . Both of these types of stability require that the object only attract trajectories from a subset of its neighbourhood (a precise definition is given in section 2 below). In this paper, we are not only interested in the stability of subcycles of the heteroclinic network, but also of the existence of trajectories which approach the network following a particular sequence of equilibria: this sequence may include visiting the same equilibrium multiple times but, for instance, visiting two different equilibria afterwards.

The stability results of Podvigina [27] are for so-called Type Z heteroclinic cycles. The equivariant version of the vector field for (1) is in this class. Specifically, all of the heteroclinic connections lie in fixed-point subspaces, are all of the same dimension, and the appropriate isotypic decomposition is into one-dimensional components. For more details of the group theoretic details, see definition 8 in [27] .

There have been some recent results giving sufficient conditions on the asymptotic stability of classes of heteroclinic network [17, 18] , which include the network we consider in this paper. However, it turns out that the network can still have very strong attracting properties even when these conditions are not satisfied. There are several other results on the stability of heteroclinic networks, for examples see [23, 24, [35] [36] [37] [38] [39] [40] [41] , but these are, in general, partial results and confined to specific examples. One source of difficulty is that there may be many different routes by which a trajectory can traverse a heteroclinic network, and Figure 1 . The network of one-dimensional heteroclinic connectionsΣ, between equilibria ξ 1 , . . . ξ 5 . Eigenvalues are shown near ξ 1 , the remainder can be deduced by symmetry. The connection orbits coloured blue are of 'Type A', and those coloured red are of 'Type B'. Later, in section 4 we give a formal definition of these.

keeping track of all possibilities in the stability calculations can be challenging. Furthermore, there must be at least one equilibrium in a network for which the unstable manifold is twodimensional: again, keeping track of trajectories that travel close to all parts of this manifold can make the computations very involved: see [38, [42] [43] [44] [45] for examples.

The main contributions of this paper are as follows. Firstly, we give an explicit method for computing the stability, not just of sub-cycles of a heteroclinic network, but of arbitrarily complex (repeating) sequences of visits to the equilibria of the network. That is, we ask whether trajectories which visit neighbourhoods of the equilibria in the network in a particular order (the list of which can be arbitrarily long, but must eventually be repeating) are attracted to, or repelled from, the network. Using the terminology of a recent preprint from Podvigina [46] , each of these sequences is an omnicycle. Secondly, we are able to adapt the conditions for heteroclinic cycle stability given by Podvigina [27] for use with continuation software to compute boundaries of stability mentioned above. In order to do this, we need to write the stability conditions in such a way that they generate a scalar function, continuous in the parameters of interest. Finally, we find a region of parameter space in which there is a complicated set of stability regions for increasingly complex sequences in which trajectories visit equilibria; these regions are reminiscent of the sausage-shaped resonance tongues seen in piecewise continuous systems [47] [48] [49] [50] [51] . A complete explanation of the shape of these regions is beyond the scope of this paper, but leaves us with many open questions and future work.

We note that there have been many other papers which prove the existence of complicated cycling behaviour of trajectories close to heteroclinic networks, and switching between different sub-cycles of the network, (e.g. [38, 40, 45, [52] [53] [54] ). To the best of our knowledge this is the first paper to give regions of parameter space where different cycling behaviours are stable (and can be observed in numerical simulations), and to explain how that stability is lost.

The remainder of this paper is organised as follows. In section 2, we give the necessary terminology and definitions required throughout the paper. In section 3 we give an overview of the dynamics of equations (1), including the stability regions of various equilibrium solutions. We also show that in addition to the heteroclinic network between the singlepopulation equilibria, there may also exist a heteroclinic cycle between equilibria each of which have three species present; we find the existence region and boundaries of stability for this cycle. In section 4 we derive a Poincaré map which captures the behaviour of trajectories close to the heteroclinic network, regardless of which route between equilibria around the network is taken by a trajectory. The derivation of the Poincaré map is complicated because we have to use a mixture of Cartesian and polar coordinates in order to capture the whole of the two-dimensional unstable manifold.

We use transition matrices in section 4.4 to analyse the map, and apply results from Podvigina [27] to determine the existence of trajectories which approach the network while visiting equilibria in repeating patterns of arbitrary length. This allows us to find regions in parameter space where different patterns are stable. A summary of these results can be seen in figure 2 . The sausage-like tongues in the centre of the figure are shown in more detail in figures 7 and 8, and each one represents a different pattern of visiting the equilibria in the network. The grey shaded region is that for which the sufficient conditions for asymptotic stability of the heteroclinic network given in [17] and [18] apply. We conjecture that the grey (uncoloured) region is actually filled with infinitely many sausages for different patterns: obviously we could only compute a finite number of these. In addition, we note that there are regions of parameter space which are outside of the grey region in which the network is still strongly attracting. The numerical results and the algorithm used to generate them are discussed in detail in section 5. Finally, in section 6 we discuss regions of parameter space in which the network appears to be attracting, but no regular pattern of visiting equilibria can be found. This irregular behaviour appears to arise in at least two different ways. Section 7 concludes. Figure 2 . Stability boundaries of various types of behaviour, in c A -c B parameter space, with e A = 1 and e B = 0.8. The labels in each region correspond to stable objects: these are described in detail in later sections, but briefly, ξ T and ξ Q are equilibria, Σ T Q is a heteroclinic cycle and A, B and AAB are various patterns of approaching the heteroclinic network Σ. The green, purple and orange shaded regions are the regions in which ξ T , Σ T Q and ξ Q (as labelled) are asymptotically stable and the boundaries are the coloured curves: blue: δ T = 1; green: λ 4 = 0; purple: δ T Q = 1, orange: stability boundary of ξ Q . These are described in section 3. The set of tongues between the regions labelled AAB, A and B are regions of fragmentary asymptotic stability of more complicated ways of approaching the network Σ. These are described in more detail in sections 4 and 5. The grey shaded region is that for which the sufficient conditions for asymptotic stability of the network Σ given in [17] and [18] apply.

We briefly summarise some notions and give some definitions that are used throughout this paper. Consider a system of ordinary differential equationṡ

Definition 1 A heteroclinic cycle is a finite collection of equilibria {ξ 1 , . . . , ξ m } of (2), together with a set of heteroclinic connections {γ 1 (t), . . . , γ m (t)}, where γ j (t) is a solution of (2) such that lim t→−∞ γ j (t) = ξ j , lim t→∞ γ j (t) = ξ j+1 and ξ m+1 ≡ ξ 1 .

A heteroclinic network is a connected union of heteroclinic cycles. More generally, heteroclinic cycles and networks may connect invariant objects more complicated than equilibria, such as periodic orbits [55] or chaotic sets [56] , but we do not consider these possibilities here. In generic systems, heteroclinic cycles and networks are of high codimension, but when (2) contains invariant subspaces, then they may exist for open sets of parameters values, that is, they are robust.

We now give some notions of stability. The notion of fragmentary asymptotic stability was introduced by Podvigina in [27] . Let X ⊂ R n be a set which is invariant under (2), and let N be an -neighbourhood of X, that is:

denote the flow generated by (2), i.e., the solution x(t) to the initial value problem starting at x(0) = x 0 . We denote the δ-local basin of attraction of X as B δ (X):

Definition 2 An invariant set X is asymptotically stable if, for all δ > 0, there exists an > 0 such that

Definition 3 (from [27] ) An invariant set X is fragmentarily asymptotically stable (f.a.s.) if, for any δ > 0,

where µ is the Lebesgue measure of a set in R n .

If a set X is fragmentarily asymptotically stable, but not asymptotically stable, then not all the points which are arbitrarily close to X are attracted to X in forward time. In the context of heteroclinic networks and cycles, this usually arises because a cusp-(or anticusp-) shaped region of phase space which abuts the set X is excluded from the basin of attraction of X.

We now give an overview of the dynamics of the system of equations (1), including the calculation of the stability of some invariants sets. Some of these are objects are subsets of the heteroclinic network which is studied for the remainder of this paper, and some are not part of this heteroclinic network. Equations (1) are equivariant under the action of the group Γ = Z 5 , generated by the element ρ which has the action

In addition, each coordinate axis, coordinate plane and three-and four-dimensional hyperplane is invariant under the flow. We label the two-, three-and four-dimensional invariant subspaces respectively as

We label the equilibria of (1) with exactly one non-zero component as ξ 1 , . . . , ξ 5 , so ξ 1 = (1, 0, 0, 0, 0), ξ 2 = (0, 1, 0, 0, 0), etc, and note that ρξ 1 = ξ 2 . The eigenvalues of ξ 1 are −1, e A , −c B , e B and −c A , with eigenvectors in the x 1 , x 2 , x 3 , x 4 and x 5 directions respectively (see figure 1 ). Since we choose the parameters c A , c B , e A and e B to be positive, it can be shown that each of the ten two-dimensional invariant subspaces P jk contains either a one-dimensional heteroclinic connection from ξ j to ξ k , or a one-dimensional heteroclinic connection from ξ k to ξ j . The resulting network of these one-dimensional connections is shown schematically in figure 1 .

The positive orthant is invariant under the flow of (1), and throughout this paper we consider the dynamics restricted to R 5 + :

. . , 5}. Within R 5 + , each equilibrium ξ j has a two dimensional unstable manifold; we define the heteroclinic network Σ to be the union of the equilibria ξ j and their unstable manifolds, namely

We also give a name to the network of one-dimensional connections: Σ = Σ ∩ ∪ 5 j,k=1 P jk There are 10 three-dimensional invariant subspaces P jkl in which two of the co-ordinates are zero. The dynamics restricted to these subspaces falls into two classes, and is either symmetric (under a power of ρ) to the dynamics in P 123 or P 124 . The dynamics of (1) restricted to each of these subspaces is shown schematically in figure 3 . Figure 3 (a) shows the dynamics restricted to P 123 : it contains an equilibrium with three non-zero coordinates . Panel (a) shows a schematic of the dynamics of (1) restricted to the subspace P 123 ; the cycle Σ T is shown in bold. Panel (b) shows the same restricted to the subspace P 124 . Eigenvalues are shown at each equilibrium ξ j . The one-dimensional connections between the equilibria shown in bold (in both panels) are part of the networkΣ; the network Σ includes the two-dimensional manifold of connections betwen ξ j and ξ j+1 .

(labelled ξ T ), and a heteroclinic cycle labelled Σ T . We discuss the dynamics in this subspace further in section 3.2. Figure 3 (b) shows the dynamics restricted to P 124 : notice the twodimensional unstable manifold of ξ 1 . More specifically, notice that there is a one-dimensional connection ξ j → ξ j+3 , but a two-dimensional connection ξ j → ξ j+1 (throughout this paper, all subscripts on equilibria, subspaces, Poincaré sections and similar objects are taken mod 5). The non-cyclic triangle of heteroclinic connections in P 124 is called a ∆-clique in [57] .

In figures 4 and 5 we show some typical time series of trajectories of equations (1), for a variety of parameter values. Figure 4 (a) shows a trajectory approaching the heteroclinic network Σ, visiting the equilibria in the order ξ 1 , ξ 2 , ξ 3 , ξ 4 , ξ 5 . Note that the length of time spent near each equilibrium increases as time increases. Figure 4 (b) shows a trajectory approaching the heteroclinic cycle Σ T Q , in which each equilibrium has three nonzero components. Again, the time spent near each equilibrium increases. Figure 4 (c) and (d) show periodic and quasiperiodic solutions which have bifurcated from the equilibrium ξ Q . The objects ξ Q and Σ T Q are defined later in this section. In figure 5 we show trajectories plotted on a logarithmic scale, so that the closeness of the coordinates to zero as the equilibria are approached can be clearly seen. Figure 5 (a) shows the same trajectory as figure 4(a), but on a longer timescale. The geometric increase of time spent near each equilibrium solution can be clearly seen. Figure 5 (b) shows a trajectory which is approaching the heteroclinic network Σ, but the equilibria are visited in an irregular manner.

In [17] and [18] , it is shown that a sufficient condition for the asymptotic stability of Σ is that min(c A , c B ) > max(e A , e B ). This region is shown by the grey shading in figure 2. As we show in the remainder of this paper, the network can have very strong stability for many regions of parameter space when this condition does not hold.

In section 4, we compute a Poincaré map between the equilibria ξ j which includes the whole of the two-dimensional manifold, and allows us to analyse the dynamics near Σ + . First, for completeness, we review the remaining equilibria of (1), their stability, and other heteroclinic cycles which can exist in this system.

In this section we consider the stability of the equilibrium of (1) with five non-zero components. We label this equilibrium as ξ Q , and its coordinates are given by

The Jacobian of (1) evaluated at ξ Q is circulant, with first row equal to

The eigenvalues of circulant matrices are well known, and in this case are given by

are the fifth roots of unity. Since Σ 5 k=0 ω k j = 0, we can simplify to get

We thus find

and

Recall that x > 0, and

For ξ Q to be stable, we need to satisfy both of these conditions, which requires

This stability boundary is shown in parameter space in figure 2 by the orange lines, and the region in which ξ Q is stable is shaded orange. Both boundaries are Hopf bifurcations. A periodic solution resulting from one of these bifurcations is shown in figure 4(c). Figure 4 (d) shows a quasiperiodic solution to which this solution evolves after further bifurcations.

In this section we consider the dynamics within the subspace P 123 , as shown schematically in figure 3(a). This subspace contains a heteroclinic cycle between ξ 1 , ξ 2 and ξ 3 . We refer to this cycle as Σ T . There are four further symmetric copies of this sub-cycle contained in Σ, related by (powers of) the symmetry ρ. The subspace P 123 also contains an equilibrium with x 1 , x 2 and x 3 non-zero, which we label as ξ T (shown also in figure 3(a)). The full system (1) contains a further four symmetric copies of this equilibrium, also with three non-zero components, related to ξ T again by powers of the symmetry ρ.

The heteroclinic cycle Σ T is equivalent to the frequently-studied Guckenheimer-Holmes cycle [1, 12] , with the removal of the rotational symmetry relating the three equilibria. Regardless of this, the stability of Σ T is simple to compute using a Poincaré map, and stability results are similar. Namely, for the dynamics of (1) restricted to P 123 , the equilibrium ξ T and the cycle Σ T exchange stability at a degenerate Hopf bifurcation when δ T =

The cycle Σ T is asymptotically stable for the dynamics restricted to P 123 if δ T > 1 and the equilibrium ξ T is stable if δ T < 1. (Note that this degenerate Hopf bifurcation can be broken with the addition of higher order terms, resulting in a branch of periodic orbits.)

For the full dynamics in R 5 , Σ T cannot be asymptotically stable (since, for instance, any points arbitrarily close to ξ 1 that lie on the heteroclinic connection between ξ 1 and ξ 4 will asymptote to ξ 4 ), but it can still be fragmentarily asymptotically stable. We consider the stability of Σ T in detail in section 4.5.2, once we have computed the Poincaré map for the network.

The equilibrium ξ T can lose stability to perturbations transverse to the subspace P 123 . The Jacobian of ξ T is block diagonal, and so the eigenvalues of ξ T can be easily computed. We label the eigenvalues with eigenvectors in the x 4 and x 5 directions as λ 4 and λ 5 , respectively. When δ T < 1 and both λ 4 < 0 and λ 5 < 0, ξ T is asymptotically stable. In figure 2, we show the curve δ T = 1 in blue, and the curve λ 4 = 0 in green (λ 5 < 0 for all parameters shown in figure 2 ). The region shaded green shows where ξ T is stable. Now suppose that δ T < 1, λ 5 < 0, but λ 4 > 0, and consider the dynamics within the four-dimensional subspace P 1234 . This subspace contains the equilibria ξ T and ρξ T , and by the symmetry ρ, within this subspace, ρξ T will be a sink. ξ T will be a saddle with a onedimensional unstable manifold (the unstable eigenvector of ξ T points in the x 4 direction), and so there will exist a heteroclinic connection between ξ T and ρξ T . This connection will be robust to perturbations preserving the invariant subspaces in the system. Again, by the symmetry ρ, there are also robust connections between ρξ T and ρ 2 ξ T , ρ 2 ξ T and ρ 3 ξ T , etc, creating a heteroclinic cycle between the five equilibria ρ j ξ T , j = 0, . . . , 4; we label this heteroclinic cycle Σ T Q . Σ T Q is not part of a heteroclinic network: each equilibrium has a one-dimensional unstable manifold, and so we can compute its stability using Poincaré maps in the usual manner (see [58] for a very similar calculation). We find that it becomes unstable as the quantity δ T Q = −λ 4 /λ 5 decreases through 1, shown by a purple curve in figure 2. The stable region of Σ T Q is shaded purple in figure 2.

In this section, we construct a Poincaré map which captures the behaviour of trajectories close to the heteroclinic networkΣ (shown in figure 1 ) between the equilibria ξ 1 , ξ 2 , ξ 3 , ξ 4 and ξ 5 . We then use this map to compute the stability of any given periodic sequence of transitions around the heteroclinic network, using the results of Podvigina [27] .

We first define some terminology. We are interested in trajectories of (1) that remain close to Σ, and so we can describe the motion in terms of the itinerary around the network, that is, the sequence of visits to the equilibria ξ j . The following notation follows that used in [59] . For fixed 0 < H < 1 2 min i,j |ξ i − ξ j |, and x ∈ R 5 + , we define

Note that the choice of maximum allowed H means that M (x) is uniquely defined. When M (F t (x 0 ))) = k we say F t (x 0 ) is close to ξ k , and note that the set of points x for which M (x) = j forms a ball around ξ j . For a trajectory F t (x 0 ) of (1) we definẽ

If F t (x 0 ) is not close to any equilibria, that is, M (F t (x 0 )) = 0, then thet in the above expression will be the time when the trajectory was most recently close to an equilibrium, and henceM (t) gives the 'last visited equilibrium'. If a trajectory starts close to an equilibrium at t = 0 this will always be non-zero. The trajectory F t (x 0 ) can be thus characterised as an itinerary m(n) of epochs τ (n):

such thatM (t) = m(n) for the interval t ∈ [τ (n), τ (n + 1)), and m(n + 1) = m(n). We say that a trajectory transitions from ξ j to ξ k at time t if there exists an n ∈ N such that m(n) = j, m(n + 1) = k and τ (n + 1) = t. For the dynamics of (1), notice that if a trajectory stays close to the heteroclinic network Σ, it can make only two types of transitions between equilibria. The first type, which we call "Type A" is a transition from an equilibrium ξ j to the equilibrium ξ j+1 . The second type of transition is from an equilibrium ξ j to the equilibrium ξ j+3 . We label this type of transition as "Type B". Given any sequence of transitions between equilibria, we can thus translate the sequence m(n) into a word in the alphabet {A, B}. Specifically, for a trajectory F t (x 0 ) with itinerary m(n), we write

We say w(n) is eventually periodic with period p and root sequence w ∈ {A, B} p , if there exists an N ∈ N such that for all n > N , w(n + p) = w(n), and w(N + j) = w (j) for j = 1, . . . , p. The minimal periodp of an eventually periodic sequence w(n) is the smallest such p. The root sequence of an eventually periodic sequence is unique up to a cyclic permutation of letters.

We give a couple of examples. If a trajectory approaches the heteroclinic cycle in the subspace P 123 between the equilibria ξ 1 , ξ 2 and ξ 3 (called Σ T in section 3.2), then the resulting itinerary will be eventually periodic with root sequence AAB. If a trajectory approaches the heteroclinic cycle between all five equilibria ξ 1 , ξ 2 , ξ 3 , ξ 4 and ξ 5 , in that order, the resulting itinerary will be eventually periodic with root sequence A. If a trajectory approaches the same five equilibria but in the order ξ 1 , ξ 4 , ξ 2 , ξ 5 , ξ 3 , then the resulting itinerary will be eventually periodic with root sequence B.

In the calculations which follow, we determine for which parameter values particular sequences are f.a.s. When we talk about the stable region for a particular root sequence, we mean this region of parameter space. Note that if there exists at least one root sequence which is f.a.s., then the network Σ is also f.a.s. Figure 7 shows a more detailed summary of our results than figure 2, including a complicated region of strings of sausage-shaped stability regions for different sequences (shown in more detail in figure 8 ). Recall that the region of parameter space that the sufficient condition for stability of Σ given by Podvigina [17] and Afraimovich [18] is min(c A , c B ) > max(e A , e B ), which is min(c A , c B ) > 1 for our choice of parameters e A and e B : this is the region shaded gray in figure 2. There are clearly large regions of parameters space outside this region in which at least one root sequence is attracting, which means that the network is fragmentarily asymptotically stable.

We follow a standard procedure described, for instance, in [24, 45] . As in [45] , we have to use a combination of Cartesian and polar coordinates in order to capture the dynamics near the whole of the two-dimensional unstable manifold.

We begin our construction of a Poincaré map by defining Poincaré sections near ξ 1 , for some h 1, as follows:

. . , 5, |x 1 − 1| < h}. As is usual in these types of calculations, there is an attracting invariant sphere [22] , and thus we know that the radial direction will not play a part in the stability calculations. Discounting the radial direction then, the Poincaré sections are three-dimensional. We label points on each section as follows: and T can be found by solving

The computation of the global map between equilibria is more subtle: depending on the outgoing coordinates on H out 1 , the trajectory may hit H in 2 or H in 4 first. From figure 6 , we can see that if θ 24 is very close to π/2, we expect the trajectory to hit H in 2 first, and if θ 24 is very close to 0, we expect the trajectory to hit H in 4 first. More precisely, within the subspace P 124 , the invariant sphere theorem means that the dynamics are restricted to a two-dimensional manifold, and so there will be some θ , such that if θ 24 < θ the trajectory first hits H in 4 , and if θ 24 > θ the trajectory first hits H in 2 . As we will see in the calculations that follow, the invariant subspaces force the (lowest order) global map to be diagonal in the coordinates x 3 and x 5 , which means that the dynamics slightly outside of the subspace P 124 will be similar.

We thus define two global maps , where θ out 24 ∈ (θ , π/2). Due to the invariance of the subspaces P 1245 and P 1234 , it is clear that to lowest order, we will have

where A 3 and A 5 are order 1 'global' constants. Note also that θ 24 = π/2 defines an invariant subspace -which is equivalent to θ 14 = π/2 (both are the subspace P 1235 ). Thus, if we write θ in 14 = g 12 (x out 3 , x out 5 , θ out 24 ) then we know that g 12 (x out 3 , x out 5 , π/2) = π/2, for any x out 3 , x out 5 . If |θ out 24 − π/2| is small, we can Taylor expand g 12 , and get to lowest order

where A θ is another order 1 global constant. Next consider the global map ψ 14 , and write ψ 14 (x out 3 , x out 5 , θ out 24 ) = (x in 2 , x in 5 , θ in 13 ), and recall that θ out 24 ∈ (0, θ ). Again, by the invariance of the subspace with x 5 = 0, we can write

for an order 1 global constant B 5 . Since x in 2 is small on H in 4 (by the definition), then we can use the invariance of the subspace x 2 = 0, and Taylor expand about zero to get

Similarly, we find that to lowest order,

We can now write down a Poincaré map from a single Poincaré section to itself, using the symmetry ρ which can map H in j to H in j+1 . To do this, we write H in ≡ H in 1 , and introduce equivariant coordinates (see, e.g. [25] ) x e A , x e B and θ c on H in and, for consistency, we use x c A , x c B and θ e on H out 1 . When the trajectory is close to ξ 1 , we have the equivalencies 24 . The local map can be written in these coordinates as:

and T can be found by solving

To write the global map in these same coordinates, we can apply the symmetry ρ to the Poincaré sections H in 2 and H in 4 , namely ρ 4 H in 2 = H in 1 and ρ 2 H in 4 = H in 1 . In the equivariant coordinates, we then find that in a neighbourhood of ξ 2 , we have x e A ≡ x 3 , x e B ≡ x 5 , θ c ≡ θ 14 , and in a neighbourhood of ξ 4 , we have x e A ≡ x 5 , x e B ≡ x 2 , θ c ≡ θ 31 . This gives the following equivariant global map

x e B , θ c ) where, if 0 < θ e < θ , to lowest order we get:

and if θ < θ e < π/2, then

where we have limited information about the function g 12 . Specifically, we only know that if |θ e − π/2| is small, to lowest order

In the global map, B 2 , B 3 , B 5 , A 3 , A 5 and A θ are order 1 global constants.

In this section we consider the relevant approximations that can be made to the Poincaré map in the case where the trajectories remain always close to at least one two-dimensional invariant subspace. That is, the angles θ e and θ c are always either close to 0 or close to π/2 in a neighbourhood of any of the equilibria. There are four cases in total, which depend on the coordinates on H in . The coordinates on H in give information both about what type of transition between equilibria will be made once the trajectory leaves H out , and also about which type of transition was made by the trajectory as it approached H in . Namely, the relative size of x e A and x e B give you information about where the trajectory will go when it leaves H out , and whether θ c is close to 0 or π/2 tells you from which direction the section H in was approached.

More precisely, we consider the relative sizes of the coordinates x A and x B on H in , in order to approximate the time T . From equation (8), we can see that if x e A x e B e A e B , then to lowest order, T = − 1 e B log(x e B ). In this case, the trajectory leaving H out will make a Type B transition (as defined in at the start of section 4). We can then compose the local and global map to get a full return map on H in that we call Φ B : . In this case, the trajectory leaving H out will make a Type A transition. When we compose the local and global maps in this case, we find the second full return map on H in , which we call Φ A :

(again, dropping the superscripts on the coordinates).

We now consider the composition of these maps as different routes around the network are taken. Recall that x A and x B are small. If a Type B transition is made, then the resulting θ will be very small. If a Type A transition is made, the resulting θ will be close to π/2. We thus define ϕ = π/2 − θ, and are able to give four different Poincaré maps, which depend on both the previous equilibrium visited, and the one which will be visited next, that is, on both the previous and current transition types.

The final results give four complete Poincaré maps, H in → H in depending on the sequence of connections. Here, Φ X→Y (X, Y ∈ {A, B}) is the appropriate map to use when the previous transition is of type X, and the current transition is of type Y .

Let the sequence of letters Z = Z 1 Z 2 . . . Z m , where each Z i ∈ {A, B}, be a root sequence of transitions for an eventually periodic sequence corresponding to a trajectory close to the network. The map given by the composition

is the return map to H in which describs the dynamics after the sequence Z of transitions has been made. For any Z, this map has a fixed point at zero, and the stability of the zero fixed point under iterations of the map corresponds to the stability of the heteroclinic network, for trajectories following this particular sequence of transitions. The following section discusses how to compute the stability of the zero fixed point in maps of this type.

In order to analyse the maps given in the previous section, and apply the results of Podvigina [27] , we first need to define the transition matrix of a map [21, 25] . Let G be the set of mappings Ψ : R p → R p that have at lowest order the form Ψ(x 1 , . . . , x p ) = (C 1 x α 11 1 x α 12 2 · · · x α 1p p , . . . , C p x α p1 1 · · · x αpp p ) for real constants α ij ≥ 0 and C i non-zero, 1 ≤ i, j ≤ p. G is clearly closed under composition. We define the transition matrix of Ψ to be the p × p real matrix M (Ψ) with entries [M (Ψ)] ij = α ij . It is easily verified that if Ψ 1 , Ψ 2 ∈ G, then

Any Ψ ∈ G has a fixed point at x 1 = . . . = x p = 0. The zero fixed point of the map Ψ will be stable if all the row sums of M (Ψ) N diverge to +∞ as N → ∞. Conversely, if any of the row sums of M (Ψ) N tends to 0, then the fixed point is unstable.

Furthermore, if the vector v (with components v j ) is an eigenvector of M (Ψ), then the curve in R p with x j = q v j , q ∈ [0, ∞), j = 1, . . . , p, is invariant under the map Ψ. If the eigenvalue λ corresponding to v is real and greater than 1, then along this curve, points contract towards the origin. If, furthermore, λ is also the eigenvalue of M (Ψ) with largest magnitude, then this curve is attracting under iteration of Ψ. Podvigina [27] extends these results to show that if these conditions on the eigenvalue and eigenvector of M (Ψ) are satisfied, then the origin in Ψ is fragmentarily asymptotically stable. That is, the basin of attraction of the origin can be extended from this one-dimensional curve to a set with positive measure. The stability of the origin in Ψ is then related to the stability of the heteroclinic cycle for which the return map Ψ was derived. A summary of the results from [27] is given in definition 6 and lemma 1 below.

We now write down the transition matrices associated with the maps Φ A→A , Φ A→B , Φ B→B and Φ B→A :

B} be a root sequence. We write M Z 1 to be the following product of m transitions matrices:

and similarly also write M Z j , for j = 2, . . . , m to be the product of m transition matrices:

We say that the set of m matrices {M Z j , j = 1, . . . m} is a collection of transitions matrices describing the root sequence Z. When we compute the stability of a root sequence, we must consider the Poincaré maps associated with each of these m transitions matrices: this is equivalent to starting the Poincaré map near each of the m different equilibria in the sequence. The stability conditions (given below) depend on the eigenvalues and eigenvectors of these matrices. All of the matrices in the collection {M Z j } will have the same eigenvalues, but the eigenvectors will in general be different.

The following definition uses the conditions for stability in Lemma 5 of Podvigina's 2012 paper [27] .

Let M be a transition matrix for a map g. Let λ max be the eigenvalue with largest absolute value of the matrix M, and w max be the associated eigenvector. Suppose λ max = 1. Then M is of fragmentary asymptotic stability type if the following conditions hold:

Note that the last condition is equivalent to requiring all the entries of the eigenvector w max to be non-zero and of the same sign.

Using this definition, we then have the following result, which follows directly from Lemma 5 in Podvigina [27] .

. . m} be a collection of transition matrices for a root sequence Z of length m. If, for each j, the matrix M j is of fragmentary asymptotic stability type, then the root sequence Z is fragmentarily asymptotically stable.

A sequence which is f.a.s. can lose this stability by violating any of the three conditions listed in definition 6 for any of the transition matrices in its collection. As noted earlier, all the transition matrices in a collection will have the same eigenvalues, so if either condition (i) or (ii) are violated, this happens for all the transition matrices in the collection at the same parameter value. The violation of condition 1. corresponds to a resonance bifurcation, and would be expected to be associated with the appearance of a long-period periodic orbit. The violation of conditions (ii) or (iii) are not associated with the bifurcation of any other invariant objects (as noted by Podvigina [27] , and seen also in a particular example in [33] .) If condition (iii) is violated, generically it will only be violated for a single one of the matrices in a collection at one time. This gives information on how instability manifests, that is, at which point in the sequence a nearby trajectory exits. We give specific examples of this when discussing the losses of stability we observe in section 5.

4.5.1. Stability of five-cycles The network contains two cycles of length 5. As discussed earlier, these cycles correspond to root sequences A, and B respectively. The stability of these is computed in [28] , and the conditions are as follows. The root sequence A is f.a.s. if

When c A + c B = e A + e B , then the largest eigenvalue of M A→A becomes equal to one. This is equivalent to a resonance-type bifurcation (see, e.g. [34, 60] ). One might expect the appearance of a long-period periodic orbit to be associated with this bifurcation, but we have not been able to find one numerically. We suspect that this is due to the lack of higher order terms in equations (1) (as is the case with the resonance bifurcation of Σ T , discussed in section 3.2).

When 

The B sequence loses stability in the same way for the first and the third of these conditions, respectively, as the A sequence does, explained above. When c B e B = c A e A , the matrix M B→B has two real eigenvalues, with the same absolute value, one of which is positive and one of which is negative.

It is clear from the second condition in each list that it is impossible for both A and B to be stable for a single set of parameters. We show the regions of stability of these sequences in the c A -c B parameter plane in figure 7.

As discussed previously, the cycle Σ T has the root sequence AAB. To compute the stability of this cycle we need to compute the three matrix products 

The same conditions can also be computed using different methods, described in [24] and [40] . When the first of these conditions is broken, i.e. δ T < 1, a resonance bifurcation occurs, which, as described in section 3.2, is degenerate, and does not result in the appearance of a long-period periodic orbit (although we might expect it to with the addition of appropriate fifth order terms to equations (1)). The remaining two conditions are of a similar type, and correspond to one of the matrices M j having an eigenvector with a zero component. In [24, 40] , the same conditions are arrived at (for very similar systems) by considering perturbations to the cycle Σ T in the x 4 and x 5 directions respectively. For the parameters we consider in figure 7, ν 5 > 0. The boundary where ν 4 = 0 is shown. Along this curve, the matrix M 1 has an eigenvector with a zero in the third component. This means that there are no initial conditions for which a trajectory Figure 7 . Stability boundaries of various sequences, in c A -c B parameter space, with e A = 1 and e B = 0.8. We use the abbreviations D ≡ BB, T ≡ AAB and Q ≡ ABBB. The blue tongues are regions of the form (AAB) n1 (BB) n2 = T n1 D n2 , where from the middle region with n 1 = n 2 = 1 (marked (AAB)(BB)), n 1 increases to the left, and n 2 increases to the right. The deepness of the blue colour is proportional to the sequence length. As the parameter c B is increased through a pinch in each sausage string, an A in the sequence is replaced by a BB. A zoom showing more details of a portion of this figure is given in figure 8 .

which performs one A transition, then a B transition, and then an A transition. Or, in terms of the original coordinates, and the cycle Σ T , a trajectory which starts on H in 2 cannot visit ξ 3 and ξ 1 and then return to ξ 2 : what is observed is that after the visit to ξ 1 the trajectory performs a B transition to ξ 4 . In the absence of the x 5 coordinate (i.e. in the subspace P 1234 ), we would see an exchange of stability between the sequences AAB and ABBB.

Although it is in principal possible to generate analytic results about the stability of root sequences of longer length, it is not particularly enlightening (and the expressions are cumbersome at best), and so we proceed to compute the regions of stability for other root sequences numerically. Specifically, for any parameter set, and root sequence Z of length m, we can compute the collection of matrices M Z j and check the conditions in definition 6. In this section we explain how to extend this notion to numerically determine the stability boundaries in parameter space for the sequence Z.

Definition 6 gives a set of conditions which, if satisfied, means that a matrix is of fragmentary asymptotic stability type (abbreviated to simply 'stable' in what follows). This results in a dichotomy: either the matrix is stable or it is not. In order to compute stability boundaries using numerical continuation (we use the software MatCont [61] ), we require a continuous variable which determines the stability, and specifically we define a real-valued scalar, which is continuous across the stability boundary, and changes sign along the stability boundary.

From definition 6, we see that a matrix is stable if the eigenvalue with largest absolute value (λ max ) lies on the half-line in the complex plane extending from 1 to ∞ along the real axis (which we refer to as L), and if all the components of the corresponding eigenvector (w max ) are of the same sign. The matrix can change from stable to unstable in three different ways: (i) by λ max moving off the half-line L; (ii) by one of the components of w max changing sign; or (iii) an eigenvalue which does not satisfy the stability conditions becoming larger in magnitude than λ max .

In order to incorporate all these ways in which stability can change, for a matrix M, we define two quantities, s λ (M) in terms of the eigenvalues, and s w (M) in terms of the eigenvectors. We define the set V to be all vectors which satisfy condition (iii) of definition 6, namely, that all components of the vector are non-zero and of the same sign.

For a matrix M, let the eigenvalue with largest absolute value be λ max . If λ max ∈ R, we label the eigenvalue with second largest absolute value as λ 2 , otherwise (that is, if λ max is one of a complex conjugate pair of eigenvalues), then we label the eigenvalue with third largest absolute value as λ 2 . Further, let the eigenvector of λ max be w max , with components w max q . We first define, for any eigenvalue λ of M, with eigenvector v λ , the quantity

The first three lines of the definition of s d give the distance of λ from the half-line L. The fourth line says that if λ ∈ L but its eigenvector is not in V , then s d is set equal to 1. We consider those eigenvalues which have s d > 0 to be "unstable" eigenvalues, and those with s d = 0 to be "stable" eigenvalues. We then define

which measures how "stable" or "unstable" the largest eigenvalue of M is, giving a positive value when λ max is a "stable" eigenvalue, and a negative value when λ max is "unstable". Specifically, we have, in the first line, the difference in absolute values of the largest and second largest eigenvalues, in the case where λ max is stable. If λ max is unstable but λ 2 is stable, then the second line gives (minus) the smaller of s d (λ max ) and the difference in absolute values of the largest and second largest eigenvalues. When both λ max and λ 2 are unstable, the third line gives (minus) the smaller of the function s d evaluated at both eigenvalues. We then define

The quantity s w is only positive when all components of the eigenvector have the same sign, and λ ∈ L. It is zero if one of the components of the eigenvector is zero (and λ ∈ L). Otherwise it is negative. It will change sign in a continuous fashion for λ ∈ L when one of the components of the eigenvector changes sign. Similarly to s λ , those eigenvalues which have s w > 0 are "stable" eigenvalues, and those with s w < 0 are "unstable" eigenvalues. Finally, we define

This gives a positive number when both s λ and s w are positive, and a negative one otherwise, and is continuous across the boundary where s = 0. For a collection of matrices {M Z j , j = 1, . . . m} we compute s(M) for all matrices in the collection, and take the minimum, in order to compute the stability of the root sequence Z.

Note that it is possible for s(M) to be equal to zero when there is not a change in stability (specifically, in the case where |λ max | = |λ 2 | and both λ max and λ 2 lie on L), but it will not change sign.

Using this algorithm, we are able to find a series of 'tongues' of stability regions of different sequences, arranged into 'strings of sausages', thirty-four of which are shown in figures 7 and 8. These sausages are reminiscent of the shape of resonance tongues in piecewise smooth systems [47] [48] [49] [50] [51] . This is not overly surprising, because the map which we use to derive these stability boundaries, is of course, defined in a piecewise fashion. However, in those systems, in each string of sausages, the same object is stable, but this is not the case here.

We observe that the strings of sausages have the following properties:

• All but two of the strings occur entirely between the curves c A e A = c B e B (the upper stability boundary of the A sequence), the curve c 3 A e B = c B e 3 A (the left-most stability boundary of the B sequence), and the curve ν 4 = 0 (the right-hand stability boundary of the sequence AAB).

• The vast majority of the stable sequences are made up only of the components (BB), (AAB) and (ABBB). For clarity, we relabel these sequences as D, T and Q respectively.

• For tongues which abutt the curve c A e A = c B e B , between any two tongues with sequences S 1 and S 2 , there is another tongue which has the sequence S 1 S 2 , giving a Farey-like sequence. In particular, we predict the existence of stable regions for all sequences of the form (AAB) n 1 (BB) n 2 , n 1 , n 2 ∈ N. We found 12 of these numerically, and they are coloured blue in figure 7 . Between each of these blue-coloured tongues, there are additional tongues -see the zoomed figure 8 -here coloured purple. Notice the large purple tongue labelled T 2 DT D between the blue regions labelled T 2 D and T D. Similarly, there are smaller purple tongues which have the sequence of the adjacent tongues.

• As the parameter c B is increased through a 'pinch' in the sausage string, an A in the sequence transforms into a BB. Normally, this results in a T transforming into a Q. For instance, notice the string with sequences T 2 D, T QD and Q 2 D in figure 7 . The region of parameter space shown in figure 7 is mostly (but not all) within the region of parameter space for which the sufficient conditions for asymptotic stability of the network Σ computed in [17] and [18] apply (namely, that min(c A , c B ) > max(e A , e B ); recall the grey shaded region in figure 2 ). We conjecture that the remaining white space between the A, AAB and BB regions in figure 7 is actually completely filled with (overlapping) stability tongues.

As well as computing the boundaries of stability for each of the tongues shown in figure 7 , we can further identify how the stability is lost. In all cases we have observed, stability is lost through the breaking of condition (iii) of definition 6, that is, one of the matrices in the collection for that sequence has an eigenvector with a zero component. Furthermore, in all the cases we have checked, it is the third component of the eigenvector which is zero. As we noted at the end of section 4.4, identifying which of the matrices M Z j has the eigenvector with a zero eigenvalue can tell us what happens when stability is lost. We demonstrate this using the stability tongue for the sequence T 2 D ≡ AABAABBB as an example.

On the left hand side of the stability tongue for the sequence T 2 D, we find that the matrix

has an eigenvector with a zero in the third component. This means that there are no initial conditions which allow trajectories to complete this sequence in the order BAABBBAA. In particular, trajectories which perform the first seven of these transitions would then perform a B, rather than an A, as the next transition in the sequence. If we assume, for now, that the final A is actually replaced in the sequence by a BB, then we get the sequence BAABBBABB, which reordered is equivalent to ABBBAABBB ≡ QT D, which is the sequence observed in the next tongue up the sausage string.

On the right hand side of the stability tongue for the sequence T 2 D, we find that the matrix

has an eigenvector with a zero in the third component. This is different to the left hand boundary. Here, it means that there are no initial conditions which allow trajectories to complete the sequence in the order ABBBAABA. Taking a similar approach, and replacing the final A in this sequence with a BB, gives the sequence ABBBAABBB ≡ QT D, again.

We see a similar pattern arise on the boundaries of the other tongues: namely, that it is an A transition which is the first transition which cannot occur once the boundary has passed, and if this A is replaced in the sequence with a BB, then the resulting sequence is the one which is found up one level in the sausage string. A detailed analysis of why this occurs is beyond the scope of this paper, but will be investigated in future work.

6. Transitions to irregular cycling behaviour Figure 5 (b) showed an example of a trajectory which approached the network, but in a complicated manner, that is, the visits to the equilibria were a non-repeating sequence. Behaviour of this kind was also noted for a network between six equilibria in [40] , and there has been other work on the existence of such switching behaviour in other types of networks [38, [52] [53] [54] .

We identify a region of parameter space where this type of behaviour appears to be typical: we can find no stable regular sequences here. In figure 2 , this parameter region is the area coloured white in the lower right corner. Note that this region is not within the region of parameter space where the sufficient condition for asymptotic stability of Σ from [17, 18] applies. However, we note from numerical simulations that the network appears to be highly attracting here. We observe two different ways in which this irregular behaviour arises. Firstly, as a consequence of quasi-periodic behaviour (resulting from bifurcations from the equilibrium ξ Q ) which gets closer to the network. Secondly, as an instability of the heteroclinic cycle Σ T Q . We present numerical evidence of both of these mechanisms below, but leave a detailed investigation of this behaviour to later work.

In figure 4(d) we noted that following the Hopf bifurcation from ξ Q , we observe quasiperiodic behaviour at c A = 1.02, c B = 0.5. In figure 9 we show two further timeseries for increasing value of c A , both in logarithmic coordinates. As c A is increased, the quasiperiodic behaviour grows in amplitude, but due to the existence of the invariant sphere, it becomes essentially 'trapped' between the equilibria of the network Σ. As a result, the trajectory spends increasingly long times in neighbourhoods of the equilibria. In figure 9 (b) we can see clear visits of the trajectory to the five equilibria ξ j , during which all but one of the coordinates are very small. However, the value of the smallest of these coordinates, which can be thought of as an approximate distance of the trajectory from Σ, does not tend to decrease over time, as can be observed in both trajectories in figure 5 . As c A is increased further, we observe that the distance of the irregular trajectories from the network decreases, although the network does not appear to be 'attracting' until c A is greater than about 1.3, when we see a sustained approach of the trajectory towards the network.

In figure 4(b) we show a trajectory approaching the heteroclinic cycle Σ T Q (defined in section 3.2). Rather than approaching the equilibria on the coordinate axes, this cycle approaches equilibria which have three non-zero coordinates, and lie in the interior of the subspaces P j,j+1,j+2 . For the dynamics restricted to these three-dimensional subspaces, these equilibria loose stability when δ T is increased through 1, shown by a blue curve in figure 2. In figure 10 we show timeseries for two parameter sets, on each side of this line, with c A = 1.6. Specifically, in panel (a), we can see that within each visit to a subspace P j,j+1,j+2 (i.e. where two of the coordinates are very small), the oscillations are decaying: the trajectory is approaching the equilibrium in the interior. In panel (b), we see the opposite: during each visit to a subspace P j,j+1,j+2 , the amplitude of the oscillations increase, and in fact we see a sequence of visits to the three equilibria ξ j , ξ j+1 and ξ j+2 . However, the number of visits made to these equilibria before we switch to the next subspace is not the same each time. This behaviour was termed irregular cycling in [58] , and in that paper we proved its existence for an open region of parameter space. In this system, this particular type of irregular behaviour seems to be restricted to a fairly small region of parameter space: it rapidly breaks down into the much more irregular behaviour of the type seen in figure 5(b) . 

In this paper we have developed a method for determining regions of parameter space in which different patterns of approaching a heteroclinic network can be found, and have applied this to a model of five-species cyclic competition. We find some complicated and intriguing patterns of stability regions in parameter space, which are reminiscent of those found in other piecewise studies. The numerical technique we have developed here could easily be used to analyse the stability of sequences of visited equilibria in other heteroclinic networks with two-dimensional unstable manifolds. This work offers several potential avenues for further study. Firstly, we would like to be able to prove some of the observations made about the patterns of stability tongues found in section 5. This would likely involve further analysis of the Poincaré maps Φ B (9) and Φ A (10), before the approximations are made which assume we are close to the onedimensional connections inΣ. From the numerical results, it appears that the transitions between the tongues in each string of sausages is associated with a change from a A pattern to an BB pattern, which is exactly when the above approximation is not valid. Secondly, our calculations tell us only the regions where sequences have a basin of attraction of positive measure; we could further compute the shape of the δ-local basin of attraction (see equation (3)) using the 'stability index' [62, 63] . We have done these calculations numerically, and found that except for the A sequences, all sequences have a 'cusp' shaped local basin -that is, one in which the measure decreases to zero as the network is approached. These numerical calculations are supported by analytical calculations in [64] , for the A, B and AAB cycle. However, these sequences are not hard to find by randomly selecting initial conditions for numerical integrations, which seems somewhat counter-intuitive. In fact, what we observe is that many trajectories wander away from a particular sequence for some transient period before settling down, meaning that the entire basin of attraction is much larger than predicted using only a local analysis. This type of behaviour was noted in both [24] and [40] ; in the latter the behaviour was termed essential quasi-asympototic stability. Thirdly, the numerical simulations clearly indicate that the network has strong attractive properties in regions of parameter space where there are no stable root sequences, and where the sufficient conditions for stability derived in [17, 18] do not apply, namely the region where we observe irregular cycling. It would be of great interest to investigate this behaviour further. It is also natural to ask the question: are the sorts of dynamics we observe for the Rock-Paper-Scissors-Lizard-Spock network typical for larger networks, or for networks in which the Z 5 -symmetry is broken? To first address the issue of broken symmetry (i.e. there is a different set of eigenvalues at each equilibria, albeit with the same signs): we believe that much of the observed dynamics would remain similar. Numerical simulations (by integration) indicate that complicated dynamics are still possible. Although we have not computed any stability boundaries for the broken symmetry case, the technique would be exactly the same, only with more book-keeping: one would now need to keep track of equilibria visited in addition to whether A or B type connections were traversed.

We also expect to see similarly complicated dynamics in larger networks, with the added complication that there would now be more than one possible network topology, even if we consider only networks between equilibria with a single species present, and retain the restriction that there is a symmetry between the equilibria which preserves the network structure (that is, all the equilibria are the 'same'). It is possible to use group-theoretic methods to compute all the possible network topologies for 'small' number of equilibria k (as in, e.g. [65] [66] [67] ), and the number grows very rapidly with k. As for the non-symmetric case, the technique we have developed in this paper for analysing the behaviour could be applied to any of these networks. A full cataloguing of the behaviour for these networks will involve classifying the networks by some measures of their topology. Initial investigations we have made have shown that networks with k odd, which contain a ∆-clique (as the Rock-Paper-Scissors-Lizard-Spock network does), have a very similar pattern of tongues of stability regions arranged into strings of sausages. Understanding why this behaviour occurs would generalise the A → BB transition discussed above.

Nonlinear aspects of competition between three species

Local dispersal promotes biodiversity in a real-life game of rock-paper-scissors

The rock-paper-scissors game and the evolution of alternative male strategies

Evolutionary game theory: Theoretical concepts and applications to microbial communities

Mobility promotes and jeopardizes biodiversity in rock-paperscissors games

When does cyclic dominance lead to stable spiral waves?

Characterization of spiraling patterns in spatial rockpaper-scissors games

Cyclic dominance in evolutionary games: a review

Pattern formations driven by cyclic interactions: a brief review of recent developments

Convection in a rotating layer: A simple case of turbulence

Structural stability of equivariant vector fields on two-manifolds

Structurally stable heteroclinic cycles

The Big Bang Theory [television series

Heritage of China: Contemporary perspectives on Chinese civilization

Symmetry breaking and branching patterns in equivariant bifurcation theory II. Archive for rational mechanics and analysis

Asymptotic stability of robust heteroclinic networks

Two-dimensional heteroclinic attractor in the generalized lotka-volterra system

Diverging fluctuations in a spatial five-species cyclic dominance game

A golden point rule in rock-paper-scissors-lizard-spock game

Stationary bifurcation to limit cycles and heteroclinic cyles

Lectures on bifurcations, dynamics and symmetry. Longman Scientific and Technical

Heteroclinic networks on the tetrahedron

A competition between heteroclinic cycles

Asymptotic stability of heteroclinic cycles in systems with symmetry. ii

Patterns of desynchronization and resynchronization in heteroclinic networks

Stability and bifurcations of heteroclinic cycles of type Z

Classification and stability of simple homoclinic cycles in R 5

Asymptotic Stability of Pseudo-simple Heteroclinic Cycles in

Transverse bifurcations of homoclinic cycles

Asymptotic stability of heteroclinic cycles in systems with symmetry

An example of a non-asymptotically stable attractor

A new mechanism for stability loss from a heteroclinic cycle

Bifurcation d'orbites périodiquesà partir d'un cycle homocline symétrique. Comptes rendus de l'Académie des sciences

Stability in simple heteroclinic networks in R 4

A heteroclinic network in mode interaction with symmetry. Dynamical Systems: An international journal

Essentially asymptotically stable homoclinic networks

A mechanism for switching near a heteroclinic network

Nonasymptotically stable attractors in o(2) mode interactions

Regular and irregular cycling near a heteroclinic network

Stability of a heteroclinic network and its cycles: a case study from Boussinesq convection

Attractors for robust heteroclinic cycles with continua of connections

Cycling chaos: its creation, persistence and loss of stability in a model of nonlinear magnetoconvection

Cycling chaotic attractors in two models for dynamics with invariant subspaces

Resonance bifurcations of robust heteroclinic networks

Behaviour of trajectories near a two-cycle heteroclinic network

How the Arnold tongues become sausages in a piecewise linear circle map

Piecewise linear models for the quasiperiodic transition to chaos

Arnol'd tongues arising from a grazing-sliding bifurcation

Border-collision bifurcations in R n

The structure of mode-locking regions of piecewise-linear continuous maps: Ii. skew sawtooth maps

Dynamics near a heteroclinic network

Switching homoclinic networks

Switching in heteroclinic networks

The effect of symmetry breaking on the dynamics near a structurally stable heteroclinic cycle between equilibria and a periodic orbit

Cycles homoclinic to chaotic sets; robustness and resonance

Almost complete and equable heteroclinic networks

A codimension-two resonant bifurcation from a heteroclinic cycle with complex eigenvalues

Quantifying Noisy Attractors: From Heteroclinic to Excitable Networks

Resonance bifurcations from robust homoclinic cycles

Matcont: a MATLAB package for numerical bifurcation analysis of ODEs

On local attraction properties and a stability index for heteroclinic connections

Stability of quasi-simple heteroclinic cycles

Stability of cycles in a game of Rock-Scissors-Paper-Lizard-Spock

Groups of order at most 6,000 generated by two elements, one of which is an involution, and related structures

A census of 4-valent half-arc-transitive graphs and arc-transitive digraphs of valence two

A census of small transitive groups and vertex-transitive graphs

This project began during less traumatic times when international travel was still possible, and we acknowledge the London Mathematical Society for financial support through a Research in Pairs (Scheme 4) grant, and the hospitality of and financial support from the Departments of Mathematics at both the University of Auckland and the University of Leeds. The majority of this research was done during COVID lockdown in 2020, and both authors are grateful to their partners and children for giving them time and space to work. We are grateful to Bernd Krauskopf for helpful conversations about the stability 'sausages', to Gabriel Verret for discussions on how to compute the possible symmetric graphs with larger numbers of equilibria mentioned in the discussion, and for some constructive comments from anonymous referees. CMP is grateful for additional support from the Marsden Fund Council from New Zealand Government funding, managed by The Royal Society Te Apārangi, and from the London Mathematical Laboratory.