Higher Order Dynamics in the Replicator Equation Produce a Limit Cycle in Rock-Paper-Scissors

Recent work has shown that pairwise interactions may not be sufficient to fully model ecological dynamics in the wild. In this letter, we consider a replicator dynamic that takes both pairwise and triadic interactions into consideration using a rank-three tensor. We study {these} new nonlinear dynamics using a generalized rock-paper-scissors game whose dynamics are well understood in the {standard} replicator sense. We show that the addition of higher-order dynamics leads to the creation of a subcritical Hopf bifurcation and consequently an unstable limit cycle. It is known that this kind of behaviour cannot occur in the pairwise replicator in any three strategy games, showing the effect higher-order interactions can have on the resulting dynamics of the system. We numerically characterize parameter regimes in which limit cycles exist and discuss possible ways to generalize this approach to studying higher-order interactions.

In this letter, we show how to modify the standard game matrix replicator dynamic by defining fitness in terms of both a game matrix (for pairwise interactions) and a game tensor (for triadic interactions). Early work exhibiting limit cycles and higher-order interactions is by Hofbauer, Schuster and Sigmund [22]. Additionally, this work is related to the prior work of Gokhale and Traulsen [23] who studied evolutionary games with multiple (more than two) strategies and multiple players using a game tensor. In this work, they study the maximum number of mixed strategy equilibria that may emerge in the resulting replicator dynamic. More recently, Zhang et al. [24] study multiplayer evolutionary games in the context of asymmetric payoffs, which we do not consider. Additionally, Peixe and Rodrigues [25] study strange attractors and super-critical Hopf bifurcations in polymatrix replicators, modelling inter and intra group interactions in a * griffinch@psu.edu † ronglingwu@mail.tsinghua.edu.cn multi-group population. Polymatrix games are also considered in [26,27]. In contrast, we study the resulting dynamics on a generalized rock-paper-scissors game (RPS) [5] with higher-order interactions in a single population. We show that the resulting dynamics arising from triadic interactions are fundamentally different from those dynamics arising from RPS in the standard replicator equation with only pairwise interactions by showing that the HOIs lead to the emergence of a limit cycle.
Rock-paper-scissors (and its generalizations) has been studied in multiple different contexts  and there are at least two schools of partially compatible dynamics. Postlewaite and Rucklidge [30, 32-34, 37, 42] have extensively studied a dynamical systems model of RPS in both spatial and aspatial cases. Their dynamics are distinct [53] from the dynamics arising from the standard replicator equation, as given in [5][6][7]. We do not consider their dynamics, but instead focus on those arising from the replicator equation. In particular, Zeeman [54] showed that RPS dynamics under the standard replicator exhibit a degenerate Hopf bifurcation and cannot produce limit cycles. More generally, Zeeman showed that no three strategy game can produce a limit cycle under the standard replicator dynamics. However, since every dynamical system arising from the standard replicator equation is diffeomorphic to a generalized Lotka-Volterra system, limit cycles and chaos may emerge for games with more than three strategies.
The main results of this letter are: (i) We propose a method for modelling triadic interactions using a simple rank-three tensor. (ii) We show (numerically) that in generalized RPS with HOIs a subcritical Hopf bifurcation occurs, and a limit cycle emerges for appropriate parameter choices. This behaviour must be caused by the HOIs, since such dynamics cannot emerge with only pairwise interactions [54]. (iii) We use a statistical analysis to construct a two-dimensional bifurcation surface, showing parameter regions where the (unique) interior fixed point is stable and admits a limit cycle, is stable with no limit cycle, and is unstable.

II. GENERAL MODEL
Let ∆ n−1 be the n−1 dimensional unit simplex embedded in R n composed of vectors u = u 1 , . . . , u n so that u 1 + · · · + u n = 1 and u i ≥ 0 for all i ∈ {1, . . . , n}. In a biological context, u i is the proportion of species i when considered as part of the total biomass to be modelled.
Let A ∈ R n×n be a payoff matrix and assume that species i receives an expected payoff as a result of both inter and intra species interactions.
Here e i is the i th standard basis vector in Euclidean space. The standard replicator equation that arises from this fitness function is given bẏ wheref (u) = u T Au. Eq. (1) assumes simple binary interactions among species, with the payoff to species i resulting from an interaction between species i and species j given by A ij .
To model higher-order interactions, redefine f i (u) to be, where B i is a quadratic form (matrix) that takes two copies of the population proportion vector u and returns a payoff to species i that occurs when one member of species i meets two members of the population (at random). In general, we could think of B i as being a slice of a (0, 3) tensor B : ∆ n−1 × ∆ n−1 × ∆ n−1 → R. The replicator equation still has form, We now propose a biologically inspired approach to defining B i . In what follows, we will define for constants of proportionality p ijk and q ijk . That is, we assume that the three-way payoff is composed of payoffs from binary interactions that are scaled to model the effects of the more complex interactions. In our analysis of generalized rock-paper-scissors, we choose p ijk and q ijk so that the Nash equilibrium of RPS remains a fixed point. It is left as a question for future work whether there are sufficient conditions on the tensor B that ensure the Nash equilibria of the game matrix A are preserved as fixed points in three-way dynamics.
Before proceeding to the analysis of RPS, we note that Eq. (2) could be generalized to include n > 3-way interactions by using higher rank tensors [23]. However, it is unlikely that such interactions are biologically meaningful. Statistical tests for these kinds of interactions are discussed in [55].

III. HIGHER ORDER ROCK-PAPER-SCISSORS
The remainder of this paper is dedicated to showing that higher-order-interactions in a generalized rockpaper-scissors game produce dynamics not seen when only pairwise interactions are modelled. Fix the parameterized payoff matrix, where we assume a, b ≥ 0. When a = b = 0, this is the standard RPS matrix found in (e.g.) [51,53,56]. This matrix will govern the payoff from simple pairwise interactions. We now define the tensor B using its slices, so that Here we assume α ∈ (1, 2] and β ∈ (0, 1] for simplicity. The reasoning behind defining B in this way is justified by considering the payoff associated to rock (index 1). If a rock plays against two other rocks, it receives no payoff -as expected from Eq. (5). If it meets two papers (index 2), then it receives the same net negative payoff −b − 1 as if it met one paper. If it meets both a rock and paper, then its net negative payoff is cut in half, since the paper plays against two rocks. On the other hand, if it meets two scissors (index 3), then its payoff increases by a factor of α. If a rock meets a rock and scissors, then the rocks split the payoff and each receives a payoff decreased by a factor β, caused by competition between the two rocks. Finally, if a rock meets both a paper and scissors, then mutual destruction leads to no net payoff for any player. The same logic mutatis mutandis is used for all other players.
Eqs. (6) to (8) enforce symmetry in each slice and preserve a generalized cyclic property. Note that B 2 can be recovered from B 1 by rotating the rows and columns of B 1 down and to the right (just as row i + 1 can be obtained from row i in A by rotating to the right). Thus, B is a circulant (or Toeplitz) tensor.
It is straightforward to see that the general replicator dynamics given by Eqs. (2) to (4) using the matrix and tensor defined in Eqs. (5) to (8) hasū i = e i (i ∈ {1, 2, 3}) as fixed points. These correspond to single species populations or pure strategies. A fourth equilibrium is the interior pointū = 1 3 , 1 3 , 1 3 , corresponding to a perfectly mixed population or the Nash equilibrium of RPS. The fact that A is a circulant matrix and B is a circulant tensor ensures that these are the only equilibria in the system that occur in ∆ 2 andū is the only interior equilibrium.
Based on the assumptions on the values of a, b and α, these are hyperbolic saddles (see Theorem 3.3 of [57]). This is identical to the behaviour of the pure strategy equilibria in the case of ordinary rock-paper-scissors in the standard replicator dynamic [7].
When a = a * + , then If follows from our assumptions on α and β that if > 0, then a > a * and r 2 < 0 andū is attracting. Otherwise,ū is repelling. For = 0, the system has two pure imaginary eigenvalues, which satisfies the first requirement of Hopf's theorem (see [58], Page 152). Our assumptions that α ∈ (1, 2] and β ∈ (0, 1] ensures that Thus, the real part of the eigenvalues must cross the imaginary axis with non-zero speed, satisfying the second criterion of Hopf's theorem. Thus, the system exhibits a Hopf bifurcation.

B. Numerical Illustration of a Subcritical Limit Cycle
We can show numerically that an unstable limit cycle emerges for example parameters. For the remainder of this section, let b = 0. In our initial limit cycle construction, we assume α = 2 and β = 1 2 . In this case, the interaction of (e.g.) a rock with two scissors doubles the payoff to the rock, while two rocks interacting with one scissors will quarter the payoff associated to the interaction. We set = 1 100 , implying the interior fixed point will be stable. In Fig. 1 (top) we see an (approximated) limit cycle surrounding a stable interior fixed point. Outside the limit cycle, flow goes to the boundary. To complete the numeric proof, we apply the Poincaré-Bendixson theorem. In Fig. 1 (bottom), we compute the distance from the three trajectories to the interior fixed point after a ternary transform. This is a true representation of the distance seen in the trajectories in Fig. 1 (top). We can see that the distance from the (approximated) limit cycle to the interior fixed point oscillates around a constant mean. The trajectory outside the limit cycle increases its distance to the interior fixed point, while the trajectory inside the limit cycle decreases its distance to the interior fixed point as expected. Thus, the numerically identified limit cycle is unstable, implying a subcritical Hopf bifurcation. When < 0, the interior fixed point becomes unstable (as expected) and the limit cycle vanishes as shown in Fig. 2 (top). All trajectories approach the boundary, as confirmed numerically in Fig. 2 (bottom), which again shows the normalized distance from the trajectory to the interior fixed point.
When is increased beyond a certain value, the limit cycle disappears while the interior fixed point remains stable. This is illustrated in Fig. 3 (top). In Fig. 3 (bottom), we show that the distance from any trajectory to the interior fixed point collapses to zero as time goes to infinity. We can numerically approximate the value of where the limit cycle disappears for arbitrarily values of α ∈ (1, 2] and β ∈ (0, 1]. To do this, we use the following steps: 1. Input: α, β. 2. Initialize: = 1 100 . Compute a = a * + , where a * is given by Eq. (9).

Numerically integrate the time-inverted dynamics:
with u 0 = 0.99, 0.005, 0.005 . This accomplishes two things: (i) The modified dynamics invert the stability of the limit cycle, and (ii) the initial condition starts the trajectory close to the boundary. Let u(t) be the resulting solution. 5. Fit d(t) ∼ γ 0 + γ 1 t.
6. If γ 1 < 0 with p-value less than 0.001, then statistically the trajectory is decaying toward a limit cycle, and we set := + 1 1000 and go to Step 3. Otherwise, stop and return .
Using this technique, we can generate an interpolated surface showing the dependence of on the parameters α and β (see Fig. 4). Before continuing our analysis, we note that a similar procedure can be used to find limit cycles within this dynamical system. Mathematica code to generate all figures is provided in the SI.
The resolution of the algorithm prevents a complete characterization of max for all (α, β) input pairs. For β sufficiently large, the algorithm returns its smallest value, and we can only deduce that max < 0.01. Using this information, we fit the empirically determined (α, β, max ) points for which max > 0.01 using a generalized linear model. The resulting fit is given bŷ This statistical analysis shows that the maximum period of a limit cycle (correlated to max ) is negatively proportional to both α and β. However, since the interaction term is non-trivial, the interaction between α and β can increase the max before the limit cycle disappears. This is consistent with the interpolated surface shown in Fig. 4. Fig. 4 (top) illustrates an approximate threedimensional bifurcation diagram (for b = 0) for the triple (α, β, ). If < 0 (i.e., (α, β, ) lies below the α−β plane), then the interior fixed point is globally unstable. On the other hand, the inputs (α, β, ) generate a limit cycle if (α, β, ) falls between the surface and the α − β plane. If (α, β, ) lies above the surface, then the interior fixed point is globally attracting. This surface is accurate only to the resolution of the algorithm, as we have already noted. The contour plot in Fig. 4 (bottom) better illustrates the curvature of the bifurcation surface.

IV. CONCLUSION
In this letter, we used the replicator equation to model higher-order interactions between three (or more) species through the introduction of an interaction tensor. We showed that the dynamics that results in this case are fundamentally different from the ordinary binary interactions modelled by the standard replicator dynamics with a payoff matrix. In particular, we studied a generalized rock-paper-scissors model and showed the existence of a non-degenerate subcritical Hopf bifurcation that allows an unstable limit cycle to emerge when higher order interactions are allowed. This is in contrast to classical results on three strategy games in the standard replicator equation, in which limit cycles cannot emerge as a result of the degeneracy of a similar Hopf bifurcation [54].
While this letter provides a framework for modelling higher-order interactions within a replicator framework, there are several possible future directions. Generalizing the interaction rules used to construct the three-way interaction tensor could lead to a generalization of the Folk theorem in certain cases. It also would be interesting to introduce a spatial component as in [53,[59][60][61][62] and determine what spatial dynamics emerge as a result of these higher-order interactions.
V. ACKNOWLEDGEMENTS C.G. was supported in part by the National Science Foundation under grant DMS-1814876. R.W. was sup-ported by the Talents Grants at Beijing Forestry University.