Networks of interacting nodes connected by edges arise in almost every branch of scientific inquiry. The connectivity structure of the network can force the existence of invariant subspaces, which would not arise in generic dynamical systems. These invariant subspaces can result in the appearance of robust heteroclinic cycles, which would otherwise be structurally unstable. Typically, the dynamics near a stable heteroclinic cycle is non-ergodic: mean residence times near the fixed points in the cycle are undefined, and there is a persistent slowing down. In this paper, we examine ring graphs with nearest-neighbor or nearest-$m$-neighbor coupling and show that there exist classes of heteroclinic cycles in the phase space of the dynamics. We show that there is always at least one heteroclinic cycle that can be asymptotically stable, and, thus, the attracting dynamics of the network are expected to be non-ergodic. We conjecture that much of this behavior persists in less structured networks and as such, non-ergodic behavior is somehow typical.

From the coupled map lattices proposed in the 1980s to the modern discipline of complex networks, the study of simple systems connected in some way forms a fundamental paradigm in dynamical systems. Applications are plentiful and diverse and include spatially extended systems, chemical reactions, biological, and ecological networks and span many length scales.^{1} In many cases, a goal of the study is *spatiotemporal chaos*,^{2} and this might typically mean computing a long-term average (for example, a Lyapunov exponent). However, it is well known that a class of networks—those with invariant subspaces forced by symmetries in the system—permit *heteroclinic cycles*,^{3} that is, trajectories along which time averages do not converge, instead slowing down as they repeatedly and systematically get closer and closer to fixed points. We investigate a family of coupled map lattices defined on ring networks and establish stability properties of the many possible families of heteroclinic cycles.

## I. INTRODUCTION

We consider how the structure of the architecture, or topology, of a network of physical nodes determines the architecture of a heteroclinic network in phase space between fixed points. We note that this question was asked in the reverse by Ashwin and Postlethwaite,^{4} who showed how to construct a system of ordinary differential equations into which was embedded a heteroclinic network of any specified topology.

In particular, we consider systems of the form

where $\gamma \u22650$ is a parameter. Each of the $xi(k)\u2208[0,1]$ measures the activity at time $i$ of the $k$th “node” in a network, and $Ajk$ is an adjacency matrix, which determines the connectivity between the different nodes. We take the function $f$ to be the logistic map,

where $r\u2208(0,4]$ is a parameter. The dynamics of the uncoupled system with $\gamma =0$ is well known,^{5} but briefly, for $r\u2208(0,1]$, the origin is an asymptotically stable fixed point. At $r=1$, there is a transcritical bifurcation creating a second fixed point at $x=r\u22121r$. This fixed point is asymptotically stable for $r\u2208(1,3]$, and at $r=3$ undergoes a period-doubling bifurcation, which leads to a period-doubling cascade followed by the onset of chaos at $r\u22483.56995$.

When $x(j)$ is non-zero for some $j$, then the term $e\u2212\gamma \u2211jAjkx(j)$ in (1), with $\gamma >0$, has an inhibitory effect upon any node connected to $j$, i.e., those for which $Ajk=1$. Specifically, if $\gamma x(j)$ is large enough, then this term can have the same effect as reducing the value of $r$ in the $x(k)$ equation to less than 1 and, hence, causing the values of those $x(k)$ to decrease toward zero. Heuristically, it is then clear that oscillatory behavior is possible, as nodes can alternately be active (have a non-zero value) and, hence, inhibit those nodes they are connected to; decay, when other nodes in turn inhibit them; and finally grow again to an active state as the nodes inhibiting them decay in turn.

In Fig. 1, we show a time series from a cycle of three such coupled nodes, specifically, the set of equations,

In panel (a), we use $r=2$, so the dynamics of the uncoupled system contains a non-zero stable fixed point. The time series clearly shows the trajectory cycling between three fixed points, in a manner essentially identical to that seen in the well-known Guckenheimer–Holmes heteroclinic cycle.^{3} In panel (b), we use $r=3.5$, so the uncoupled system is in the chaotic regime, and we see cycling between three chaotic attractors. This phenomenon was previously described by Ashwin *et al*.^{6,7}

In this paper, we extend the work of Ashwin *et al.*^{6,7} and consider larger networks of coupled systems in the form of (1). We refer to these equations as describing the network of connections between nodes in *physical space*, and for the remainder of the paper, refer to this network as instead a *directed graph* with *directed edges* between *nodes*. We begin in Sec. II, by considering another example: a five-node ring graph with one-way nearest-neighbor coupling. We determine the fixed points and the heteroclinic connections which exist between them. We refer to this network of connections as the *phase space network*, or *heteroclinic network*, which has *heteroclinic connections* (or sometimes simply *connections*) between *fixed points*. In Sec. III, we consider general systems of the form of (1) and describe how to find the fixed points and heteroclinic connections for such a system. In general, this procedure results in a very complex heteroclinic network that is difficult to analyze, so in Sec. IV, we look in detail at $N$-node directed graphs with one-way nearest-neighbor coupling in the physical space. Here as well as determining the structure of the heteroclinic network in phase space, we are able to analyze the dynamic stability of subcycles within the network. We use results from Podvigina^{8} and some classical results on solutions to polynomials^{9,10} to prove Theorem 1, which shows that only one of the subcycles can ever be stable, and then, only if $\gamma $ is large enough. In Sec. V, we make some conjectures about larger networks. Section VI concludes.

## II. EXAMPLE: FIVE-NODE RING GRAPH WITH NEAREST-NEIGHBOR COUPLING

We give, in this section, an example of a five node directed graph with one-way nearest-neighbor coupling. We determine the possible fixed points in phase space and the heteroclinic connections between them. We show time series of typical trajectories close to the resulting heteroclinic network but defer the computation of the stability of each of the sub-cycles to Sec. IV.

We begin with a few formal definitions. Consider the map

with fixed points $\zeta 1,\u2026,\zeta M$. A *heteroclinic connection* between $\zeta j$ and $\zeta j\u2032$ is a solution to (3) for which $xi\u2192\zeta j$ as $i\u2192\u2212\u221e$ and $xi\u2192\zeta j\u2032$ as $i\u2192\u221e$. Suppose that there exist heteroclinic connections between $\zeta jk$ and $\zeta jk+1$ for $k=1,\u2026,M\u22121$ and also one between $\zeta jM$ and $\zeta j1$. Then, the set $H$ consisting of the fixed points $\zeta j$ and the connecting orbits is a *heteroclinic cycle*. A *heteroclinic network*^{11} is usually defined to be a union of heteroclinic cycles. In this paper, we relax the definition somewhat: we allow a heteroclinic network to consist of a set of fixed point solutions and heteroclinic connections between them, which contains at least one heteroclinic cycle. Note that this means that not every heteroclinic connection in the network need be part of a cycle.

For the example, we consider in this section, the network of nodes in physical space is shown in Fig. 2(a). The equations governing this system are

which are equivariant with respect to a rotation symmetry of the coordinates.

Equation (4) has two different types of fixed points solution in phase space, which are of interest to us, namely, those with one node active (that is, with a single component that is $O(1)$) or those with two nodes active. There may be other fixed points in this system, but they are not part of the heteroclinic network of interest in this paper. More precisely, assume for now that $r\u2208[1,3]$ and let $x^=r\u22121r$. Then, using coordinates $(x(1),x(2),x(3),x(4),x(5))$, we label the fixed points with only $x(1)$ active as

and similarly, we have $\xi 2,\u2026,\xi 5$, where for each $\xi j$, the $j$th component is equal to $x^\u22600$, and the remainder are equal to zero. Next, we label

and similarly define $\xi j,m$, which has the $j$th and $m$th components equal to $x^$, where $|j\u2212m|\u22601$ (i.e., $j$ and $m$ are not adjacent nodes in the ring graph). Note that $\xi j,m\u2261\xi m,j$ but we typically list $j$ and $m$ in increasing numerical order.

Note that in this example, there cannot be any fixed points with more than two components equal to $x^$ because of the connectivity of the graph: two nodes that are connected by an edge cannot both be active at a fixed point. In larger, more general graphs, we would expect to see fixed points with more active components (see Sec. III for the general setup).

We next consider the dynamics of (4) in two two-dimensional subspaces and show that there exist heteroclinic connections from $\xi j$ to $\xi j\u22121$ and from $\xi j$ to $\xi j,j\xb12$ (where indices are taken $mod5$).

For the first, consider the dynamics in the subspace where $x(2)=x(3)=x(4)=0$, namely, the system

System (5) has fixed points at $(x(1),x(5))=(x^,0)\u2261\xi 1$ and $(x(1),x(5))=(0,x^)\u2261\xi 5$. Consider an initial condition close to $\xi 1$ but with $x0(5)\u22600$. Since the $x(5)$ equation is decoupled from $x(1)$, then it behaves as it would in the uncoupled logistic map, specifically, the $x(5)$ component initially grows and approaches the value $x^$. As $x(5)$ grows, the coupling term in the $x(1)$ equation has the effect of essentially reducing the $r$ value of the logistic map in the $x(1)$ equation to, eventually, $re\u2212\gamma x^$. Thus, if $re\u2212\gamma x^<1$, then $x(1)$ will eventually decay to zero, and the trajectory approaches $\xi 5$. There is, thus, a heteroclinic connection from $\xi 1$ to $\xi 5$, and by symmetry, heteroclinic cycles from $\xi j$ to $\xi j\u22121$.

For the second type of connection, consider the dynamics in the subspace where $x(2)=x(4)=x(5)=0$, namely, the system

In this subspace, both $x(1)$ and $x(3)$ are decoupled from each other. There are three fixed points, $(x(1),x(3))=(x^,0)\u2261\xi 1$, $(x(1),x(3))=(0,x^)\u2261\xi 3$, and $(x(1),x(3))=(x^,x^)\u2261\xi 1,3$. Both $\xi 1$ and $\xi 3$ are saddle points, and perturbations close to these fixed points will result in trajectories, which approach $\xi 1,3$. There are, thus, heteroclinic connections between $\xi 1$ and $\xi 1,3$, $\xi 3$ and $\xi 1,3$, and by analogy, heteroclinic connections between any $\xi j$ and $\xi j,m$ or $\xi m,j$.

Finally, we consider the dynamics of (4) in a three-dimensional subspace and show that there is a heteroclinic connection from $\xi j,j+2$ to $\xi j\u22121,j+2$. Consider the dynamics in the subspace where $x(2)=x(4)=0$, namely, the system

There are, as before, fixed points in this system with one component non-zero, but of interest right now are the two fixed points with $(x(1),x(3),x(5))=(x^,x^,0)\u2261\xi 1,3$ and $(x(1),x(3),x(5))=(0,x^,x^)\u2261\xi 3,5$. Note first that the $x(3)$ and $x(5)$ components are decoupled. Both have stable fixed points at $x(i)=x^$. Thus, perturbations close to $\xi 1,3$ will have an $x(3)$ component, which remains close to $x^$, but, as in the case of the connection between $\xi 1$ and $\xi 5$, the $x(5)$ component will grow, and again, so long as $re\u2212\gamma x^<1$, the $x(1)$ component will decay to zero.

Due to the rotational symmetry of the system (4), we, thus, have three families of heteroclinic connections, namely,

Note that here, indices are taken $mod5$. In later sections, indexes are taken $modn$, where $n$ is the number of nodes in the graph. We refer to the families of heteroclinic connections as $1\u21921$ connections, $1\u21922$ connections, and $2\u21922$ connections, respectively, where a $p\u2192q$ connection is a transition from a fixed point with $p$ nodes active to a fixed point with $q$ nodes active.

The complete set of connections between fixed points is shown in panel (b) of Fig. 2. Notice that there are two heteroclinic cycles, one between fixed points of type $\xi j$, with $1\u21921$ connections, and the other between fixed points of type $\xi j,k$, with $2\u21922$ connections. The $1\u21922$ connections are not part of any cycles. In Fig. 3, we show numerical simulations showing trajectories close to each of these cycles. We will compute the stability of cycles of these types in general in Sec. IV that follows. We shall see that the cycle between fixed points of type $\xi j,k$ can have some form of stability if parameters are chosen correctly, specifically, if $\gamma >3log\u2061r2x^$, but the cycle between fixed points of type $\xi j$ can never be stable. However, if initial conditions are chosen carefully (in a manner described in that section), we can, as seen in panel (b) of Fig. 3, observe this cycle for a reasonably long period of time.

## III. ENUMERATION OF FIXED POINTS AND HETEROCLINIC CONNECTIONS FOR A GENERAL DIRECTED GRAPH

In this section, we describe how to find the fixed points and heteroclinic connections in phase space for any directed graph with inhibitory coupling. In Sec. IV that follows, we apply this to an $n$-node graph with nearest-neighbor coupling.

### A. Enumeration of fixed points

For convenience and readability, we restate the general system (1) with $N$ nodes,

where $Ajk$ is an adjacency matrix. We enumerate the fixed points that can occur in this system, specifically, those with one or more non-zero coordinates. Fixed points can have any number of non-zero coordinates, so long as the corresponding nodes are not adjacent in the physical network. More formally, consider a partition of the first $N$ natural numbers into two sets,

with $J<N$, $\alpha k,\beta k\u2208{1,\u2026,N}$. Then, the point with

is a fixed point of (1) if $A\alpha ^\alpha =0$ for all pairs $(\alpha ^,\alpha )\u2208Z+\xd7Z+$. We label this fixed point $\xi Z+$. In the language of graph theory, $Z+$ is called an *independent set*.

### B. Existence of heteroclinic connections in phase space

Consider a fixed point $\xi Z+$ in system (1), with $|Z+|=J$ (i.e., there are $J$ nodes active). We label the set of *suppressed* nodes at $\xi Z+$ to be $Zs(Z+)$, where

where for each $al$, there exists at least one $\alpha k\u2208Z+$ with $Aal\alpha k=1$. We further define the *suppression number* of a node $al$ to be the number of different $\alpha k\u2208Z+$ with $Aal\alpha k=1$. The remaining nodes are the *growing* nodes, and we define $Zg(Z+)={1,\u2026,N}\u2216(Z+\u222aZs(Z+))$.

It is simple to check that the linearization of system (1) about $\xi Z+$ has the following eigenvalues:

$J$ eigenvalues equal to $2\u2212r$, with eigenvectors in each of the directions corresponding to the active nodes.

$s(Z+)$ eigenvalues in the suppressed directions. Each of these will be equal to $re\u2212ns\gamma x^$, where $ns$ is the suppression number of that node.

$N\u2212J\u2212s(Z+)$ eigenvalues equal to $r$. These are the

*growing*nodes.

We assume that $1<r<3$ and $re\u2212\gamma x^<1$, so each fixed point is a saddle. There are two ways in which heteroclinic connections between fixed points can arise.

Consider a fixed point $\xi Z+$ and let $b\u2208{1,\u2026,N}\u2216Z+$. Consider the subspace in all components in $Z+\u222a{b}$ are fixed at zero. There are then three possible cases:

$b\u2208Zs(Z+)$ and so $\xi Z+$ is a sink in this subspace. There are no heteroclinic connections from $\xi Z+$ in this subspace.

$b\u2208Zg(Z+)$ and $Z+\u2229Zs({b})=\u2205$. Then, there is a heteroclinic connection from $\xi Z+$ to $\xi Z+\u222a{b}$. This is a heteroclinic connection of type $|Z+|\u2192|Z+|+1$.

$b\u2208Zg(Z+)$, and $Z+\u2229Zs({b})=Z\u2212$ is non-empty. Then, initial conditions near $\xi Z+$ will have an increasing $x(b)$ component. All $x(a)$ (for $a\u2208Z\u2212$) will eventually decay, and the trajectory will asymptote toward the fixed point $\xi Z+\u222a{b}\u2216{Z\u2212}$ (which has node $b$ active but nodes in $Z\u2212$ inactive). This is a heteroclinic connection of type $|Z+|\u2192|Z+|\u2212|Z\u2212|+1$.

Note that for the directed graph with nearest-neighbor coupling, each node only inhibits one single other node, and so in case 3 above, $|Z\u2212|=1$ always. Thus, for those examples, the number of active nodes can increase via way of heteroclinic connections but it can never decrease.

## IV. DIRECTED GRAPH WITH NEAREST-NEIGHBOR COUPLING

In this section, we consider the case of a general $n$-node ring graph, with one-way nearest-neighbor coupling. That is, system (1) with $A$ a cyclic permutation matrix, given by equation

We refer to these graphs as $(N,1)$-graphs.

Note that system (9) is equivariant with respect to the group $ZN$, generated by the element $\sigma $, which has the action

As for the five-node system, Eq. (9) has $n$ fixed point solutions $\xi j$ each with the $j$th component equal to $x^$, and all other components zero. Other fixed points are found using the method in Sec. III A.^{12} Note that for any set of nodes $Z+$, $Zs(Z+)=|Z+|$, that is, the number of suppressed nodes is the same as the number of active nodes, and each suppressed node has a suppression number of one.

### A. Examples of heteroclinic networks

Using the methods described in Sec. III, we give some further examples of the heteroclinic networks, which occur in $(N,1)$-graphs.

#### 1. (6, 1)-graph

The $(6,1)$-graph has six fixed points with one non-zero component, nine with two non-zero components, and two with three non-zero components. The latter are stable and the remainder are saddles. Each of the $\xi j$ fixed points with one non-zero component has heteroclinic connections to $\xi j\u22121$, and to the three fixed points $\xi j,j+2$, $\xi j,j+3$, and $\xi j,j+4$ with two non-zero components. There is, thus, a heteroclinic cycle between the nodes of type $\xi j$. Note that the fixed points with two non-zero components can be divided into two classes depending on whether the spaces between the active components is $2$ and $2$ (the fixed points $\xi j,j+3$) or $1$ and $3$ (the fixed points $\xi j,j+2$). The subset of the heteroclinic network between just these fixed points is shown in Fig. 4 .

#### 2. (7, 1)-graph

The $(7,1)$-graph has 7 fixed points with 1 non-zero component, 14 with 2 non-zero components, and 7 with 3 non-zero components. The fixed points with two non-zero components can further be divided into two classes of types $\xi j,j+2$ and $\xi j,j+3$. The heteroclinic network between the fixed points with two or three non-zero components is shown in Fig. 5. Each of the $\xi j$ fixed points, which is not shown here, will have heteroclinic connections to $\xi j\u22121$, $\xi j,j+2$, $\xi j,j+3$, and $\xi j,j+4$, again, creating a heteroclinic cycle between the nodes of type $\xi j$.

### B. Symmetric subcycles

As can be seen from the examples given so far of the $(5,1)$-, $(6,1)$-, and $(7,1)$-graphs, there can exist many different heteroclinic cycles within the heteroclinic network in phase space. As the number of nodes in the ring graph increases, so too will the number of heteroclinic cycles. In the following, we establish stability results for some of these subcycles. Although in theory our method can be used for any subcycle, in practice, it is much easier to compute explicit stability conditions if the cycle is *symmetric*. That is, if the cycle is between $n$ fixed points $\zeta 1,\u2026,\zeta n$, then there exists a symmetry $\rho =\sigma M$ (for some $M$) such that $\zeta j=\rho \zeta j\u22121$ for all $j\u22082,\u2026,n$ (and $\zeta 1=\rho \zeta n$).

Note that for both the $(6,1)$- and the $(7,1)$-graph, such cycles exist between the $\xi j$ fixed points. For the $(7,1)$-graph, symmetric cycles also exist between the $\xi j,j+3$ fixed points and the $\xi j,j+2,j+4$ fixed points. However, for the $(6,1)$-graph, there is no symmetric cycle between fixed points with two non-zero components (see again Fig. 4).

In this section, we enumerate the possible symmetric cycles in $(N,1)$-graphs. Clearly, the symmetry requires the number of active coordinates at each fixed point to be the same, and additionally, that the spacing of the active coordinates around the ring is the same at each fixed point. This gives restrictions on the allowed spacing between the active nodes, as follows.

In Fig. 6(a), we show an $n$-node ring graph and suppose that there are three nodes which are active (colored blue). The nodes colored red are those which are being suppressed by the blue nodes, and all others are colored green. Suppose that there are $n1$, $n2$, and $n3$ green nodes between each pair of blue and red nodes (as marked in the figure) (where $n1\u22651$, $n2,n3\u22650$). In order for the next fixed point to have the same number of active components, the next node to reach $O(1)$ much be adjacent to a blue node: the heteroclinic connection must be of type $3\u21923$. Without loss of generality, we suppose that the next node to reach $O(1)$ is the node marked with a yellow star. The next fixed point in the cycle, thus, has nodes colored as in Fig. 6(b) (the yellow star marks the same node). There are then two options for how we could rotate the arrangement at the second fixed point to match the arrangement at the first. Either (i) we rotate panel (b) anticlockwise by $n1+1$ nodes, and must, therefore, have

which implies $n2=n3=n1\u22121$; or, (ii), we rotate panel (b) clockwise by $n1+1$ nodes, and then we must have

which implies $n1=n2=n3+1$.

If there are instead $J>3$ nodes active, rather than just three, similar arguments can be made, and the results give the same two possible cases. In case (i), we write $n1=p\u22121$, $n2,\u2026,nJ=p\u22122$, and the total number of nodes is $N=pJ+1$, for some $p\u22652$. In case (ii), we write $nJ=s\u22123$, $n1,\u2026,nJ\u22121=s\u22122$ and the total number of nodes is $N=sJ\u22121$, for some $s\u22653$.

In the case where there is only a single active node ($J=1$), then clearly all fixed points are symmetric. In the case where $J=2$, the same arguments apply as for $J=3$ or more, except that there is no distinction between cases (i) and (ii) above.

Note that in both cases (i) and (ii), the spacing between the active nodes is such that all gaps between active nodes are of equal length except one that is one greater or fewer than the others. In Fig. 7, we show four possible types of fixed points in the graph with $N=11$, with two, three, four, and five active nodes, respectively. The fixed points with three [Fig. 7(b)] and four [Fig. 7(c)] active nodes are in case (ii): one gap is smaller than the others. The fixed point with five active nodes [Fig. 7(d)] is in case (i): one gap is bigger than the others.

The maximum value that $J$ can take (the maximum number of active nodes at a fixed point) is equal to $N/2$ if $N$ is even, and $(N\u22121)/2$ if $N$ is odd. We refer to these fixed points as *maximally active* fixed points.

### C. Construction and analysis of heteroclinic cycles

We now specifically construct a heteroclinic connection between two fixed points and use this construction to show how the stability of a heteroclinic cycle between fixed points can be computed. We first consider the dynamics within an *epoch*, which we define to be the length of time a trajectory spends in a neighborhood of a single fixed point. We then discuss how the trajectory transitions between epochs. This method echoes the construction of Poincaré maps, which is typical in the analysis of heteroclinic cycles in continuous-time systems.^{8,11,13,14}

#### 1. Dynamics within one epoch

Consider the dynamics of (9), with $N$ nodes in the graph, near a fixed point at which $J$ nodes are active, that is, an equilibrium at which $J$ components are equal to $x^$. Then, by the arguments given in Sec. III B, there will be $J$ negative eigenvalues and, hence, $J$ components, which are decaying. Similarly, there will be $N\u22122J$ components that are growing. We give a schematic sketch of this in Fig. 8 and note that the only initial condition we have specified is that one of the decaying components starts at $O(1)$. We have further specified that one of the growing components reaches $O(1)$ at the end of the period of time shown in the figure. We label this period of time one *epoch* and note that after this epoch, the trajectory will be in the neighborhood of a different fixed point, which may have the same, or one more, nodes that are active.

We label each of the growing components, in order of largest to smallest initial conditions,

where $L=N\u22122J\u22121$, and we label each of the decaying components, again in order of largest to smallest initial conditions,

We use a superscript “in” and “out” to indicate the initial and final conditions for each of the components. Recall that each of the growing components grows at a rate $r>1$, and each of the decaying components decays at a rate $rexp\u2061(\u2212\gamma x^)<1$. Let the number of iterations in the epoch shown in Fig. 8 be $T$, and then we have

or, assuming that $xein\u226a1$ and, hence, $T$ is large

We can then use this expression for $T$ to compute the “out” coordinates of the other components in terms of the “in” components. Specifically, we find

We write $\delta =(\gamma x^log\u2061r\u22121)$, $Xc=log\u2061(xc)$, etc. and then have the following linear map from the “in” variables to the “out” variables,

where $1m$ is a length $m$ column vector of $1$’s, and $Im$ is the $m\xd7m$ identity matrix.

#### 2. Transitions between epochs

In this section, we discuss how to map the “out” variables given in Eq. (11) onto a new set of “in” variables for the next epoch.

In Fig. 9, we show how the nodes are labeled (as in the time series schematic shown in Fig. 8) for case (i). The node with the yellow star is $xe$, as described above. The remaining growing nodes are $xs1,\u2026,xsL$, where $L=N\u22122J\u22121$. As per the labeling scheme in Fig. 8, $xs1$ is the component that will become $O(1)$ next in the sequence (following $xe$), and by applying the necessary rotations between fixed points, it can be seen how this node is selected. Similar arguments explain how the other growing nodes are labeled. The labeling of the contracting nodes is done in a similar fashion: $xc$ is the component that was $O(1)$ at the start of the epoch and, hence, was blue in the previous fixed point: application of the rotation between panels (a) and (b) in Fig. 6 gives us this label. The $xti$ labels are found in a similar way.

This results, for case (i), in the following transformation between the “out” coordinates of the last fixed point, and the “in” coordinates of the next one,

Combining this with the linear map in (11), we get that the logarithmic coordinates in each epoch are transformed under the map,

where $M$ is called a *transition matrix* and is given by

The matrix $M$ is a $q\xd7q$-square matrix, where $q=N\u2212J\u22121=J(p\u22121)$. There are $L=N\u22122J\u22121=J(p\u22122)$ rows starting with a $\u22121$ and $J$ rows starting with a $\delta $, and $1$’s on the upper diagonal. Here, $\delta =\gamma x^log\u2061r\u22121$ as before.

In case (ii), the labeling of the coordinates can be computed in the same way (although the labeling turns out to be different), but the resulting map is exactly the same. That is, $M$ is a $q\xd7q$-square matrix, with $q=N\u2212J\u22121=J(s\u22121)\u22122$, the number of rows of $M$ starting with a $\u22121$ is $L=N\u22122J\u22121=J(s\u22122)\u22122$, the number of rows starting with a $\delta $ is still $J$.

In the original $x$ coordinates, the map (12) has a fixed point at $x=0$, which corresponds to the heteroclinic cycle in the original system. Podvigina^{8} gives results on the stability of this heteroclinic cycle, dependent on the properties of the eigenvalues and eigenvectors of the transition matrix. In Sec. IV D, we give a brief heuristic argument explaining Podvigina’s results and then state the precise requirements for stability.

### D. Transition matrices and fragmentary asymptotic stability

We begin this section with some formal definitions, referring back to a generic dynamical system of form (3).

We define the $\epsilon $-local basin of attraction of a set $H$, invariant under $f$, as $B\epsilon (H)$,

From Podvigina,^{8} we also have the following:

*fragmentarily asymptotically stable*if, for any $\epsilon >0$,

Now, suppose that for a heteroclinic cycle $H$, we have derived a Poincaré map in logarithmic coordinates, as in Sec. IV C, of the form

Here, the subscript index $i$ now counts *epochs*, rather than individual iterations of the original map. Let $Xi$ have dimension $q$, so $M$ is an $q\xd7q$ matrix and then we can write the initial condition $X0$, in the basis of eigenvectors $vj$ of $M$, i.e.,

where the $cj$ are scalars. Then, we find

where $\lambda j$ is the eigenvalue corresponding to the eigenvector $vj$. Let $\lambda max$ be the eigenvalue with the largest absolute value; then, the leading order term of $Xi$ is

Recall that $Xi$ are logarithmic variables, so $xi\u21920$ if $Xi\u2192\u2212\u221e$, that is, in order for the heteroclinic cycle to be stable, we require at least that $|\lambda max|>1$. However, in addition, the $Xi$ values are required to stay real and negative, so additional conditions are required, namely, that $\lambda max$ is real, and that all the entries in the eigenvector $vmax$ are of the same sign. Podvigina shows that if these conditions are satisfied, then there exists an open set of initial conditions that remain close to the heteroclinic cycle for all time, more specifically, the heteroclinic cycle is *fragmentarily asymptotically stable*.

#### (Adapted from Podvigina8)

^{8})

Let $M$ be a transition matrix for a heteroclinic cycle $H$. Let $\lambda max$ be the eigenvalue with the largest absolute value of the matrix $M$, and $vmax$ be the associated eigenvector. Suppose $\lambda max\u22601$. Then, $H$ is fragmentarily asymptotically stable if the following conditions hold*:*

$\lambda max$ is real,

$\lambda max>1$, and

$vmaxlvmaxj>0$ for all $l,j$.

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

### E. Stability calculations

In this section, we will prove the following:

Let $\delta \u2217=q+1\u2212JJ$. If $q=J$, then the corresponding heteroclinic cycle is f.a.s. if $\delta >\delta \u2217$ and is unstable otherwise. All heteroclinic cycles with $q>J$ are unstable.

Note that $\delta >\delta \u2217$ is equivalent to $\gamma >log\u2061rx^N\u2212JJ$.

We prove Theorem 1 by presenting results about the eigenvalues of the matrix $M$. First, note that the characteristic polynomial of $M$ is

(this follows, e.g., from Claim 1 in Postlethwaite and Dawes,^{15} p. 629), and recall that $q=J(p\u22121)$, $p\u22652$, $J\u22651$. To establish the properties of the roots of the polynomial $P(\lambda )$, we will appeal to three classical results.

#### (Descartes rule of signs)

**(Descartes rule of signs)**

For a polynomial with real coefficients, ordered by descending variable exponent, the number of positive roots of the polynomial is either equal to the number of sign changes between consecutive $($nonzero$)$ coefficients, or is less than it by an even number. A root of multiplicity $k$ is counted as $k$ roots.

Of the $q$ complex roots of $P(\lambda )$, exactly one, $r+$, is real and positive, by Theorem 2. Note that since $P(1)=(q\u2212J)\u2212\delta J$, we have $r+>1$ if and only if $\delta >\delta \u2217$. Next, given a polynomial

with $an\u22600$, define

and note that Theorem 2 shows that $f+(x)$ has exactly one real positive root, $f+^>0$, and $f\u2212(x)$ has exactly one real positive root $f\u2212^>0$.

#### (Cauchy9)

**(Cauchy**

^{9})We will also use the following:

#### (Rouché10)

**(Rouché**

^{10})Let $f$ and $g$ be functions analytic inside and on a simple closed contour $C$, and suppose $|g(z)|<|f(z)|$ on $C$. Then, both $f$ and $f+g$ have the same number of zeros inside $C$ $($with each zero counted as many times as its multiplicity$)$.

We will prove the following:

The polynomial $P(\lambda )$ (with $\delta >0$) has roots satisfying the following*:*

When $q=J$, $r+>1$ is the root of the largest magnitude.

When $q>J$ we have the following cases

*:*If $J$ is odd and $q$ is even, the root of largest magnitude is real and negative. Call this root $r\u2212$; then, there are $q\u22121$ roots inside $|\lambda |=r\u2212$, one of which is $r+$, and $q\u22122$ of which are complex.

If $J$ and $q$ are both even, or if $q$ is odd, then there are $J\u22121$ roots inside $|\lambda |=max{r+,1}$, of which exactly one, $r\u2212$, is real if and only if $J$ is even, and there are $q\u2212J$ complex roots outside $|\lambda |=r+$.

First, in the case $q=J$, $P(\lambda )=P\u2212(\lambda )$, and so by Theorem 3, all roots of $P(\lambda )$ are bounded in magnitude by $r+$. Recall that $r+>1$ when $\delta >\delta \u2217$.

- When $q>J$, Theorem 3 is no longer of use. It is convenient to study the related polynomialwhich has the same roots as $P(\lambda )$, plus a root at $\lambda =1$. Considering $Q\u2032(\lambda )=(q+1)\lambda q\u2212J(1+\delta )\lambda J\u22121$, and recalling that a double root of a polynomial is also a root of the polynomial’s derivative, we see that $r+=1$ when $\delta =\delta \u2217=(q+1\u2212J)/J$, and that$Q(\lambda )=(\lambda \u22121)P(\lambda )=\lambda q+1\u2212(1+\delta )\lambda J+\delta ,$when $\delta >\delta \u2217$ (and $Q\u2032(r+)<0$ when $\delta <\delta \u2217$).(16)$Q\u2032(r+)=(q+1)r+q\u2212J(1+\delta )r+J\u22121>0$
- In the case $J$ odd and $q$ even, we haveso $Q(\u2212r+)=\u2212Q(r+)+2\delta =2\delta >0$. Since $Q(\lambda )\u2192\u2212\u221e$ as $\lambda \u2192\u2212\u221e$, there must exist a real root of $Q(\lambda )$ between $\u2212\u221e$ and $\u2212r+$. Thus, $Q(\lambda )$, and hence also $P(\lambda )$, has a real negative root, $r\u2212$, of greater magnitude than $r+$. The same argument applies to $Q(\u22121)$, so $r\u2212$ is also of greater magnitude than the root at $\lambda =1$. To show that (in this case) $r\u2212$ is the root of largest magnitude, we apply Theorem 4, setting $f(\lambda )=\lambda q+1$ and $g(\lambda )=\u2212(1+\delta )\lambda J+\delta $. We have, on the circle of radius $|\lambda |=(1+\u03f5)|r\u2212|$,$Q(\u2212\lambda )=\u2212\lambda q+1+(1+\delta )\lambda J+\delta =\u2212Q(\lambda )+2\delta ,$where the first inequality follows since $|r\u2212q+1|=|(1+\delta )r\u2212J|+\delta $ and the second by the triangle inequality. Then, by Theorem 4, $Q(\lambda )=f(\lambda )+g(\lambda )$ has the same number of complex zeros inside that circle as $f(\lambda )$, namely, $q+1$. Hence, $r\u2212$ is the root of largest modulus.$|f(\lambda )|=|(r\u2212(1+\u03f5))q+1|>|(1+\delta )(r\u2212(1+\u03f5))J|+\delta \u2265|(1+\delta )(r\u2212(1+\u03f5))J\u2212\delta |=|g(\lambda )|,$
- In the other cases, we will again apply Theorem 4. First, consider $\delta <\delta \u2217$, so $r+<1$ and take $f(\lambda )=\u2212(1+\delta )\lambda J$ and $g(\lambda )=\lambda q+1+\delta $. For $\u03f5<2(1+q\u2212J(1+\delta ))/(q(q+1))$ (observing that $\u03f5>0$ since $\delta <\delta \u2217$), we have, on the circle $|\lambda |=1\u2212\u03f5$,$|f(\lambda )|=|\u2212(1+\delta )(1\u2212\u03f5)J|>(1+\delta )(1\u2212J\u03f5)=1+\delta +\u03f5(\u2212J(1+\delta ))>1+\delta +\u03f5(\u2212(q+1)+\u03f5q(q+1)/2)=1\u2212\u03f5(q+1)+\u03f52q(q+1)/2+\delta >(1\u2212\u03f5)q+1+\delta >|(1\u2212\u03f5)q+1+\delta |=|g(\lambda )|.$Similarly, if $\delta >\delta \u2217$, so $r+>1$, we inspect $f(\lambda )$ and $g(\lambda )$ on the circle $|\lambda |=(1\u2212\u03f5)r+$, with $\u03f5<2(1+q\u2212J(1+\delta )/rq+1\u2212J)/q(q+1)$ [which is again positive by (16)]. We have, similarly,$|f(\lambda )|=|\u2212(1+\delta )(1\u2212\u03f5)Jr+J|>(1+\delta )(1\u2212J\u03f5)r+J=r+q+1+\delta +\u03f5r+J(\u2212J(1+\delta ))>(1\u2212\u03f5)q+1r+q+1+\delta >|(1\u2212\u03f5)q+1rq+1+\delta |=|g(z)|.$

Theorem 4 implies that $Q(\lambda )$ has $J$ roots inside the circle $|\lambda |=1\u2212\u03f5$ (respectively, $|\lambda |=(1\u2212\u03f5)r+)$ if $\delta <\delta \u2217$ (respectively, $\delta >\delta \u2217)$, and so $P(\lambda )$ has $J\u22121$ roots inside those circles.

#### Proof of Theorem 1

When $q>J$, the conditions of Lemma 1 are not met, as $\lambda max$ is no longer real and positive.

### F. Appearance of instabilities

As described in the section above, many of the heteroclinic cycles we find are unstable. However, if initial conditions are carefully chosen, then the cycles can be observed for reasonably long times in numerical simulations. Specifically, suppose that a heteroclinic cycle $H$ has transition matrix $M$, with eigenvalues and corresponding eigenvalues $\lambda i$ and $vi$. Suppose that the heteroclinic cycle is unstable so that the eigenvalue with largest magnitude, which without loss of generality, we assume to be $\lambda 1$, does not satisfy the conditions of Lemma 1. Assume further that $\lambda 2$ *does* satisfy the conditions of Lemma 1. Then, if we choose initial conditions $X0=c2v2$, then the forward trajectory will remain close to $H$.

In numerical simulations, of course, errors accumulate, and the trajectory can only remain close to $H$ for a finite time. In Fig. 10(a), we show an example of a trajectory, which remains close to an unstable cycle for a long time. Here, the cycle in question is the one between fixed points with one node active in the $(5,1)$-graph. In Fig. 10(b), we show the coordinates from (a) at the bottom of each “valley” in the time series: this corresponds to the coordinates at the transition between epochs and, hence, the coordinates $Xj$. For the particular transition matrix for this cycle with the noted parameters, $\lambda 1$ and $\lambda 2$ satisfy the assumptions given above, and $\lambda 1$ is complex. Using these values for $\lambda 1$ and $\lambda 2$, we use least squares to estimate the values of $c1$, $c2$, and $c3$ to fit the curve $Xj=c1\lambda 2j+c2|\lambda 1|jcos\u2061(jarg\u2061(\lambda 1)+c3)$ to the obtained data. The dashed line in (b) shows the curve $Xj=c1\lambda 2j$, that is, the data that would be expected if there were no numerical errors and we were able to start exactly on the required eigenvector. The solid curve includes the second term and is clearly an excellent fit to the data points.

## V. ANALYSIS OF RING GRAPH WITH *m*-NEAREST-NEIGHBOR COUPLING

In this section, we expand on our results from Sec. IV to discuss ring graphs with $m$-nearest-neighbor coupling ($m<N/2$). We find that, depending on the number of nodes, $N$, in the graph, and the number $m$ of neighbors coupled, different types of heteroclinic networks can arise in the dynamics. Some of these have dynamics, which can be described using the same methods as in Sec. IV, and some of these are more complex. We refer to the $N$-node graph, with $m$-nearest-neighbor coupling as the *$(N,m)$-graph*.

The smallest graph, which falls into this category (with $m\u22601$), is the a five-node graph, with $m=2$, shown in Fig. 11(a). In this example, it is not possible to have any fixed points, which have more than one component non-zero, and the network of heteroclinic connections between these fixed points is shown in Fig. 11(b). There are two-subcycles in this network between five fixed points, and the transition matrices for these cycles can be found using the methods in Sec. IV. We find, for the cycle $\xi 1\u2192\xi 5\u2192\xi 4\u2192\xi 3\u2192\xi 2$, that the transition matrix is

The eigenvalues can be found explicitly as $\xb1\delta $, $\u22121$, and, thus, this cycle can never be fragmentarily asymptotically stable. For the cycle $\xi 1\u2192\xi 4\u2192\xi 2\u2192\xi 5\u2192\xi 3$, the transition matrix is

which has eigenvalues $\delta $, $\xb1i$, but the eigenvector for the eigenvalue $\delta $ has a zero in the second component and so this cycle can also never be fragmentarily asymptotically stable. There are other routes trajectories can take while still approaching the network: in fact, this network is equivalent to the Rock–Paper–Scissors–Lizard–Spock network investigated by Postlethwaite and Rucklidge (for ODEs),^{16} which has some very complicated dynamics: see Fig. 12 for a typical time series.

As a second example, consider the seven-node graph with two-nearest-neighbor coupling, shown in Fig. 13(a). In this example, there are seven fixed points, which have exactly one non-zero component, and seven with exactly two non-zero components. The heteroclinic network between these fixed points is shown in Fig. 13(b). Note the similarity in structure between this network and the network shown in Fig. 2(b). The stability of the cycles between the fixed points with either one or two non-zero components can be computed in exactly the same way as shown previously, and we find that the cycle between the fixed points with two non-zero components can be stable if $\delta $ is large enough but the other cycle cannot.

As in the examples of the $(N,1)$ graphs, when transitioning from one fixed point to another along a heteroclinic connection, the number of active nodes may increase, but it can never decrease. We refer to those fixed points with the largest number of active nodes as the *maximally active* fixed points, and in the following, discuss the possible network of heteroclinic connections between these nodes for a general $(N,m)$-graph. We have the following:

If $N=0mod(m+1)$, then all maximally active fixed points are asymptotically stable, and there are no heteroclinic connections.

If $N=1mod(m+1)$, then all maximally active fixed points have an unstable manifold of dimension one, and there exists a heteroclinic cycle between the fixed points. The transition matrix for this heteroclinic cycle takes the form of $M$ in Eq. (13), with all rows starting with a $\delta $. It can, thus, be asymptotically stable if $\delta $ is sufficiently large.

If $N=pmod(m+1)$, $p\u22600,1$, then all maximally active fixed points have an unstable manifold of dimension $p$, and there exists a heteroclinic network between the fixed points. We conjecture that this network can be asymptotically stable for large enough $\delta $ but may have complex dynamics.

## VI. DISCUSSION

In this paper, we have shown that heteroclinic networks can typically arise in the phase space dynamics of certain types of symmetric (physical space) graphs with inhibitory coupling. We further showed that at most, one of the subcycles can be stable for an open set of parameters and, hence, observable in simulations. Many studies of coupled map lattices and complex networks seek asymptotic behavior described by a Sinai–Ruelle–Bowen (SRB) invariant measure.^{17} However, the dynamics associated with a stable heteroclinic cycle preclude this behavior—the dynamics is not ergodic, and long-term averages do not converge. In particular, averaged observed quantities, such as Lyapunov exponents, are ill-defined and will oscillate at a progressively slower rate.

From this work arises the more general question of whether or not a stable heteroclinic cycle is likely to be found in the corresponding phase space network of a randomly generated physical network of nodes. We performed some preliminary investigations on this question numerically, for randomly generated Erdos–Rényi graphs (where each edge exists with some fixed probability). We find that the probability of the existence of heteroclinic cycles in the phase space network increases both as the number of nodes in the physical network increases and also as the density of edges in the physical network decreases. However, even in cases where the probability of the existence of heteroclinic cycles is very high, there is also a very high chance of the existence of a stable fixed point in the phase space. Thus, the question of the stability of the heteroclinic cycle is important in determining whether or not the heteroclinic cycle, and associated slowing down of trajectories, will be observed in the phase space associated with a randomly generated graph. The methods we describe in this paper can be used to determine the stability of any specific heteroclinic cycle but as yet, it is not clear how one would determine the likelihood of a heteroclinic cycle to be stable in such a randomly generated network.

In this work, we consider specific dynamics for each single node in the directed graph in physical space; specifically, we suppose that there is only an “on” state and an “off” state. If more general dynamics are allowed, then other types of heteroclinic cycles can be found in ring graphs.^{18}

An obvious extension of this work would include different types of coupling and/or different dynamics in the uncoupled nodes. For example, neural networks may also be modeled with dynamics, such as phase-coupling and pulse-coupling.^{19} Specifically, a situation that might better exemplify neuronal dynamics could include both inhibitory and excitatory types of connections, and nodes could require a “kick” from an “on” excitatory connection in order to leave a stable zero state. Networks of this type were investigated by Ashwin and Postlethwaite,^{20} although they made no attempt to classify the possible heteroclinic networks that could occur. Further work on this avenue of investigation is ongoing.

## ACKNOWLEDGMENTS

The authors thank two anonymous referees for helpful comments on the original version of the manuscript. C.M.P. acknowledges support from the Marsden Fund Council from New Zealand Government funding, managed by The Royal Society Te Apārangi (Grant No. 17-UOA-096).

## AUTHOR DECLARATIONS

### Conflict of Interest

The authors have no conflicts to disclose.

## DATA AVAILABILITY

Data sharing is not applicable to this article as no new data were created or analyzed in this study.

## REFERENCES

*Mathematical Proceedings of the Cambridge Philosophical Society*(Cambridge University Press, 1988), Vol. 103, pp. 189–192.