Spatial dynamics of higher order rock-paper-scissors and generalisations

We introduce and study the spatial replicator equation with higher order interactions and both infinite (spatially homogeneous) populations and finite (spatially inhomogeneous) populations. We show that in the special case of three strategies (rock–paper–scissors) higher order interaction terms allow travelling waves to emerge in non-declining finite populations. We show that these travelling waves arise from diffusion stabilisation of an unstable interior equilibrium point that is present in the aspatial dynamics. Based on these observations and prior results, we offer two conjectures whose proofs would fully generalise our results to all odd cyclic games, both with and without higher order interactions, assuming a spatial replicator dynamic. Intriguingly, these generalisations for N⩾5 strategies seem to require declining populations, as we show in our discussion.


I. INTRODUCTION
Replicator dynamics have been used extensively in theoretical ecology to model ecosystem interactions at a high level [1][2][3].Surprisingly, these models intersect those from theoretical physics, with tournament dynamics in ecology [4,5] also occurring in the analysis of Schrödinger operator [6] and in the discrete KdV equation [7].
Most biological and ecological models assume pairwise interactions [8][9][10][11], leading naturally to generalized Lotka-Volterra equations or the (equivalent) replicator dynamics in which the interaction matrix and the payoff matrix become synonymous.In this case, the payoff from interactions is used to define species fitness, as discussed in Section II.This simple assumption is invalidated by strong evidence for the existence of higher order interactions [2,[12][13][14][15][16][17][18][19][20][21].Higher order interactions occur when three or more species (not necessarily distinct) interact with each other simultaneously to produce an additional payoff, which may increase or decrease fitness in the constituent species [2,9,[12][13][14][15][22][23][24].In particular, higher order interactions have the potential to alter the established relationship between diversity and stability [13].
While the replicator equation has been studied extensively [25][26][27][28], the replicator dynamic with higher order interactions has recently been considered by Griffin and Wu [29].In this work, they show that the presence of higher order interactions in rock-paper-scissors can change the well-known dynamics of this game to allow the emergence of a sub-critical Hopf bifurcation as compared to the known degenerate Hopf bifurcation that characterizes the dynamics of rock-paper-scissors under the ordinary (pairwise) replicator dynamics [30].Before this, Gokhale and Traulsen [31] studied evolutionary games with multiple (more than two) strategies and multiple players, while Zhang et al. [32] study multiplayer evolutionary games with asymmetric payoffs.In related but distinct work, Peixe and Rodrigues [33] study strange attractors and super-critical Hopf bifurcations in polymatrix replicators.Polymatrix games are also discussed in [34,35].However, to our knowledge, no one has yet studied a spatial replicator with higher order interactions, which is the goal of this paper.
Spatial evolutionary dynamics using partial differential equations have been studied by several authors, with [36][37][38][39][40][41][42][43][44][45][46][47] providing a small example of the body of work.Most of these models assume an infinite, spatially homogeneous, population in so far as the state variables of the model are the proportions of the population playing a given strategy at a given location and time.Durrett and Levin were the first to point out the fundamental differences between discrete and continuous evolutionary game models and finite and infinite population assumptions [48].These distinctions have been further by Griffin et al. [49,50], where it is shown that finite populations can destroy travelling wave solutions (in rock-paper-scissors) or even reverse the direction of travelling waves (in prisoner's dilemma).
Alternative approaches to studying finite populations frequently use discrete (grid) based methods and are based on the early work of Nowak and Martin [39], with extensions by several authors [51][52][53][54][55][56][57][58][59].These models often focus on the interplay between concepts from statistical mechanics and evolutionary games via updating rules that use (among other mechanisms) the Boltzmann distribution.We will not consider these models in this paper.
Instead, we will use the models of Vickers [42] for infinite (or spatially homogeneous) populations and Griffin, Mummah and Deforest for finite (or spatially inhomogeneous) populations.In this paper, we study spatial replicator equations with higher order interactions for both infinite (spatially homogeneous) and finite (spatially inhomogeneous) populations.Formal definitions for these cases are provided in Section II.While Griffin et al. [49] show that rockpaper-scissors under the ordinary spatial replicator dynamic can only admit travelling waves if the net population is decreasing, we show that the introduction of higher order interactions allows travelling waves to emerge in spatially homogeneous and inhomogeneous populations with no decline.Interestingly, when we generalise to cyclic games with more strategies (e.g., rock-paper-scissors-Spock-lizard), we see that this property of both the existence of travelling waves and a non-declining population seems to be a property of the three strategy case only.Nevertheless, we use observations made in this paper to pose two general conjectures on travelling waves and cyclic games under both the ordinary spatial replicator and the spatial replicator with higher order dynamics.
The remainder of this paper is organized as follows.In Section II, we introduce notation needed in the remainder of the paper.We formulate the higher order spatial replicator in Section III.Our analysis on rock-paper-scissors is carried out in Section IV.We generalise this analysis in Section V, proposing two conjectures on odd cyclic games.Conclusions and future directions are discussed in Section VI.There is also an appendix (A) that contains a derivation of the first Lyapunov coefficient for the Hopf bifurcation identified in the Section IV.

II. BACKGROUND
Let ∆ n−1 be the n − 1 dimensional unit simplex embedded in R n composed of vectors u = ⟨u 1 , . . ., u n ⟩ where u 1 + • • • + u n = 1 and u i ≥ 0 for all i = 1, . . ., n.We assume that an ecosystem supports a population size of M total species.Then u i M is the size of the population of species i, where we allow fractional species counts for simplicity.
Suppose the fitness of species i is given by the function f i (u).The replicator equation with fitness f is then, where, is the mean fitness of the population.Assuming a finite population, the dynamics of the whole population are given by, If A ∈ R n×n is a payoff (or interaction) matrix, then, produces the classic replicator from evolutionary game theory, When u is a function of space x and time t.Vickers [42,43,[60][61][62][63][64] (and many others) study the spatial replicator with form, where D is a diffusion constant.Without loss of generality, we assume that all species share a diffusion constant.Griffin, Mummah and DeForest [49] generalised the work of Durrett and Levin [48] to show that when the total population M (x, t) is neither homogeneous nor infinite, the species and total population are governed by the system of equations, where n is the spatial dimension.It is straightforward to see that when M (x, t) is homogeneous (or infinite), then the Eq. ( 8) is recovered.As we will discuss these two cases in the remainder of the paper, we will refer to equations of the form given in Eq. (3) as the finite population spatial replicator and equations of the form given in Eq. (2) as the infinite population spatial replicator, even though we may really be considering finite populations that are spatially inhomogeneous vs. homogeneous.A biased rock-paper-scissors (RPS) payoff matrix is given by, where we assume that a > −2 to maintain a rock-paper-scissors dynamics.It is well known that with this payoff matrix, the aspatial replicator, Eq. ( 1), has a single interior fixed point at u 1 = u 2 = u 3 = 1 3 and this fixed point is stable when a > 0, unstable when a < 0 and elliptic when a = 0 [25].
Griffin, Mummah and DeForest [49] showed that a travelling wave solution exists for Eq. ( 2) using the biased RPS matrix just in case a < 0. In a finite population case, this is biologically unrealistic since, which is negative just in case a < 0. Since, M ̇= f ¯M < 0 the population will collapse in this case.Moreover, [49] shows numerically that the travelling waves can be destroyed in the finite declining population case.However, we know that spatial travelling waves exist in real, non-declining populations [65][66][67][68].Our objective is to show that higher order interactions lead to the existence of travelling waves in cyclic competition (rock-paper-scissors) under both the finite and non-finite spatial replicator dynamics.

III. HIGHER ORDER INTERACTIONS IN ROCK-PAPER-SCISSORS IN SPACE
In [29], Griffin and Wu introduce a higher order interaction dynamic modelled by, 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 randomly meets two members of the population.We think of B i as being a slice of a (0, 3) tensor The mean fitness is then given by, Following Vickers, [43] we can construct a spatial model for higher order interactions that assumes a homogeneous population by appending a diffusion term to the replicator to obtain, Let A be the standard rock-paper-scissors matrix, obtained by setting a = 0 in Eq. ( 4).Now, u T Au = 0. Generalising from Griffin and Wu [29], we assume that the quadratic forms B i (i = 1, 2, 3) can be written as, where we assume, α, β, γ, δ > 0. As in [29], the tensor B, composed of slices B 1 , B 2 and B 3 , has cyclic structure.When we assume that γ = 2β and δ = 2α, then, and consequently, f ¯= 0. Thus, any finite population would be stable assuming these dynamics.The resulting spatial dynamics for a homogeneous (infinite) population are, The corresponding finite population model is then, Notice that our assumption on A and B implies that f ¯= 0 and so the bulk population is governed by the diffusion equation.

IV. TRAVELLING WAVE SOLUTIONS EXIST IN ONE DIMENSION
We begin by analysing the aspatial dynamics.Under the assumption that γ = 2β and δ = 2α, the aspatial dynamics are given by, Just as with ordinary rock paper scissors, the three extreme points of ∆ 2 are fixed points as is the interior point . First order analysis of the Jacobian matrix at the interior fixed point gives eigenvalues, Thus the interior fixed point is unstable when β > α and stable if β < α.When β = α, the Hartman-Grobman theorem cannot be used.In this case, the dynamics simplify to, which is just an ordinary rock-paper-scissors dynamic with payoff matrix (1 + 2α)A.Therefore, the interior fixed point is elliptic in this case.Moreover, we have shown that the higher order dynamics we consider have analogous dynamics to the ordinary RPS system, except that by construction f ¯= 0. Now consider the spatial replicator with infinite (homogeneous) population.Let z = x + ct, where c is a wave speed to be determined.If we have u(x, t) = u(z), then the resulting system becomes, where u ′ i is the derivative in terms of z.Let v i = u ′ i .Then we have the modified system of differential equations, This system has a fixed point at v i = 0, u i = 1 3 for i = 1, 2, 3. Computing the eigenvalues of the Jacobian at this point gives, The zero eigenvalue arises because we necessarily have u r (z) + u p (z) + u s (z) = 1 and v r (z) + v p (z) + v s (z) = 0 and thus the dynamics play out on a 4 dimensional manifold and λ 1 and λ 2 can be ignored.Focusing on the term under the outer radical, assume there is some r so that, Then we obtain the equations, We can compute the wave speed and the parameter r as, We conclude that the wave speed is non-imaginary, just in case β > α.That is, a travelling wave can emerge when the interior fixed point of the aspatial dynamics is unstable and hence stabilised by the diffusion term.This is similar to the condition found by Griffin, Mummah and Deforest for the ordinary spatial replicator with rock-paper-scissors [49].
We can simplify the eigenvalues λ 3,4 and λ 5,6 using the negative branch of (r ⋆ , c ⋆ ) to obtain, Thus we have three eigenvalues with negative real part indicating a three-dimensional stable manifold with two additional eigenvalues that are pure imaginary.The presence of a stable manifold with imaginary eigenvalues satisfies the first criterion of Hopf's theorem [69] (Page 152).We use the negative branch because that will ensure that solutions to the PDE are (locally) attracted to the limit cycle and hence the travelling wave solution.Now consider the specific eigenvalues, Differentiating with respect to c and evaluating at the identified wave speed yields, Then, To see that this is always non-zero, note that the equation, is quadratic in α and β.Solving for α in terms of β leads to a quadratic equation with discriminant, Thus, there are no real values of α and β that make this expression zero.As such, the eigenvalues must cross the imaginary axis with non-zero speed, satisfying the second criterion of Hopf's theorem.Thus, we have proved the existence of a Hopf bifurcation at the fixed point, which implies the existence of an isolated attracting period orbit (stable limit cycle) just in case the first Lyapnuov coefficient of the system's normal form is non-zero and negative [69].The first Lyapunov coefficient can be constructed using techniques in [70,71], as, , where, The details of this construction are provided in A and the SI, where it is also shown that this quantity is always negative.Thus, by Hopf's theorem, we have proved that the travelling wave system Eq.( 16) has an attracting periodic solution (because ℓ 1 < 0) and consequently a travelling wave solution must exist for Eq.(8).
We now consider the one-dimensional finite population model from Eq. ( 9).In one dimension we have, Following work by Griffin [50], we have a travelling wave solution for the diffusion equation, where A and B are arbitrary constants and kc ∈ R is the population wave speed, and c ∈ R will be the species wave speed.Assume B = 0. Then Then in the finite dimensional case, the resulting travelling wave equation for the species is, leading to the system of equations, This is identical to Eq. ( 16) but with a modified wave speed and consequently, our proof of the existence of a travelling wave solution applies mutatis mutandis.Thus, for small diffusion (D < 1 2 ), the finite and infinite populations will share solutions but travelling at different speeds.
We now illustrate this for = 1  10 , β = 3 2 and α = 1, and simultaneously show the existence of the predicted limit cycle solution for Eq. ( 16).Consider the following initial conditions for the PDE's Eq. ( 8) and Eq. ( 9), )︃]︃ , and assume periodic boundary conditions u i (−π, t) = u i (π, t).Then four snapshots of the resulting travelling wave solution are shown in Fig. 1.We see a perturbation of the initial condition that quickly settles into the travelling  FIG. 1.Both the finite and infinite population spatial replicator with higher order interactions in RPS converge to a travelling wave solution with distinct wave speeds depending on whether a finite and spatially inhomogeneous or spatially homogeneous (or infinite) population is modelled.
wave solution in both the finite and infinite population cases.We can numerically investigate solutions for Eq. ( 16).Let u i (x, t) be the (numerical) travelling wave solutions to the infinite (finite) population spatial replicator.For u i (z) and v i (z) in Eq. ( 16), we set, where T = 50 is sufficiently large to ensure that the resulting numeric solution is (effectively) on the limit cycle.When we plot u i (0, t) for i = 1, 2, 3 (in an appropriate projection) we see that the solutions to the partial differential equation(s) approach the limit cycle, as expected.This is shown in Fig. 2 Recall that when β > α, the interior fixed point u i = 1 3 (i = 1, 2, 3) in the aspatial dynamics is unstable.We conclude that the travelling wave solution arises because the diffusion is stabilising the growing oscillations that would arise at all points in space and (under certain initial conditions), allowing the stabilised oscillations to synchronize.
We can prove this stabilisation occurs by first order analysis of the infinite population system.Let J 0 be the Jacobian of the equation system given in Eqs. ( 10) to ( 12) evaluated at the interior fixed point.Then we have, Let υ i = u i − 1 3 with υ = ⟨u 1 , u 2 , u 3 ⟩ and let D = DI, where I is the identify matrix.Following [72], we analyse the linearised stationary problem with Neumann boundary conditions, by computing the roots of the characteristic polynomial, The PDE solutions converge to a travelling wave solution, which corresponds to an attracting limit cycle in the equivalent system of ordinary differential equations.The (left) shows the spatially homogeneous or infinite population spatial replicator with higher order interactions in RPS.The (right) shows the finite spatially inhomogeneous population spatial replicator with higher order interactions in RPS.
Here k is a wave number in a Fourier basis of a proposed solution ansatz and λ is an eigenvalue.We find three eigenvalues, The fact that −Dk 2 appears in the real parts of all three eigenvalues is sufficient to show that the diffusion exerts only a stabilising effect.Moreover, is positive only if β > α + 3Dk 2 .That is β > α, which we already knew.Thus, we have not only shown that the diffusion exerts a stabilising effect on the system, but also that Turing patterns cannot emerge in this system as a result of diffusion induced instability.Interestingly, this seems also explain the occurrence of travelling waves when no higher order interactions are present but a < 0 in the interaction matrix in Eq. (4) using the infinite population spatial replicator as shown in [49].We discuss this as a possible framework for generalising these results in future directions.
While it is generally difficult to construct the amplitude of a limit cycle, and thus a travelling wave, we can show numerically that the amplitude of the travelling wave (limit cycle) varies inversely with (a function of) the diffusion constant.Thus, as D increases, we expect to see travelling wave solutions that approach the fixed point solution u i (x, t) = 1  3 , further demonstrating the stabilising effect of the diffusion.This is illustrated in Fig. 3 for in the infinite population case.these solutions lead to a globally oscillating solution that asymptotically approaches the boundary of the simplex at all spatial positions.For the finite population case, we are using the travelling wave solution used in our prior numerical illustration.Interestingly, the finite population case takes longer to approach the globally oscillating solution than the infinite population case, most likely as a result of the bulk movement of the finite population.This is illustrated in Fig. 4(bottom).This phenomenon may warrant investigation in future work.

V. GENERALISATION
To generalise the work in this paper, recall that a circulant matrix [73] has structure, That is, the entire matrix structure is determined by the first row.The set of all circulant matrices forms an algebra under addition and (commutative) matrix multiplication.
Let N = 2n + 1 with n ≥ 1.Consider the N dimensional row vector, where a is the biasing term.Then the N × N circulant matrix A N defined by A N1 is the payoff matrix for the N strategy generalisation of rock-paper-scissors.Griffin and Fan [74] showed that the replicator dynamics Eq. ( 1) have a unique interior fixed point at u = ⟨︁ 1 N , . . ., 1 N ⟩︁ that is stable when a > 0 and unstable when a < 0. Then we have the following conjecture, which is proved for the case N = 3.
Conjecture 1.For all odd N , the one-dimensional spatial replicator equation, admits a travelling wave solution when A N is defined as above and a < 0.
We hypothesize that As in the three strategy case (rock-paper-scissors), when a < 0, uA N u < 0, which implies a globally decreasing population.
Generalising our result to the higher order interaction case produces a surprising result.To generalise the interaction tensor to the case of N strategies, let Σ = {1, . . ., N } be the strategy set and let w(i) denote the set of strategies that are beaten by strategy i and let l(i) be the set of strategies that beat strategy i.Then, if j = i and k ∈ w(i) or k = i and j ∈ w(i) −α if j = i and k ∈ l(i) or k = i and j ∈ l(i) 0 j ∈ w(i) and k ∈ w(j) and i ∈ w(k) 0 j ∈ l(i) and k ∈ l(j) and i ∈ l(k) The last three cases produce the 0 diagonals and the case when the three strategies form a winning/losing cycle (like rock, paper and scissors).As before, α, β, γ, δ > 0. When N = 3, we recover the higher order interaction matrices we have already studied.
Consider N = 5.Then evaluating at γ = 2β and δ = 2α gives, This value is 0 if and only if α = β and otherwise its sign is equivalent to sgn(α − β).Notice, the triples are composed of entries of the form u i u j u k where j, k ∈ l(i) and therefore cannot occur in the case when N = 3, which is why f ¯= 0 in our prior analysis.Simple computation shows that the only way for f ¯to be zero is in the case when α = β.Further analysis shows that the eigenvalues of the Jacobian at the (unique) interior fixed point u i = 1 5 (i = 1, . . ., 5) are, Thus the interior fixed point is unstable just in case β < α, as in the three strategy case, which we conjecture will lead to a stable travelling wave solution in the infinite population spatial case.We summarize this in the following conjecture.
Conjecture 2. For all odd N , the one-dimensional higher order spatial replicator equation, admits a travelling wave solution when A N and B Ni (i = 1, . . ., N ) are defined as above and a = 0, γ = 2β and δ = 2α and β > α.
We note that the computation of the first Lyapunov coefficient will most likely be the most difficult component of any proofs of the generalised conjectures.
What is perhaps the most interesting aspect of this is the fact that higher order interactions as defined by Eq. ( 20) seem to be able to simultaneously produce travelling wave solutions and maintain a constant population size only for the three strategy rock-paper-scissors game.While we do not rule out the possibility that a more complex interaction mechanism may be able to simultaneously accomplish this, it is surprising that this property seems to hold for only the three strategy case and thus may warrant additional study.

VI. CONCLUSIONS AND FUTURE DIRECTIONS
this paper, we merged the higher order interaction model first discussed by Griffin and Wu [29] with the spatial replicator equation model of Vickers [43] and the finite population spatial model of Griffin, Mummah and Deforest [49].For higher order interactions in rock-paper-scissors, we showed that travelling wave solutions exist in both the finite and infinite population cases, with the important model feature that the net population was stable (as opposed to declining).This suggests that if replicators are models of real-world cyclically interacting populations, then travelling waves in such populations can be explained by either a declining population count or the presence of higher order (i.e., non-pairwise) interactions.In discussing a generalisation of this approach, we provided two conjectures on the existence of travelling wave solutions in spatial replicator dynamics with an arbitrary odd number of strategies.Most interestingly, we found that the property of population size preservation and the existence of travelling wave solutions appears to be present only in the rock-paper-scissors game (three strategy case), with higher order interactions.Games with more than three strategies (e.g., rock-paper-scissors-Spock-lizard) seem to admit travelling wave solutions only when the total population is decreasing and higher order interactions cannot remediate this.
Proving the conjectures raised in this paper is clearly an area of future work.We argue that stable Turing patterns will not be admitted by the infinite population spatial replicator with higher order interactions as defined in this paper.However, Griffin and Wu [29] show that a (subcritical) Hopf bifurcation can emerge in the aspatial higher order dynamics using a related but distinct payoff matrix and higher order interaction matrices.If parameter regimes exist where a supercritical Hopf bifurcation exists in the aspatial case, then it may be possible that a diffusion mediated transition maybe possible from periodic solutions to asymptotically stable solutions as in the work of Ginzburg and Landau equation [75] or in the work of Dilão [76].Investigating this possibility of significant interest for future work.
To prove this value is always negative, and thus that the limit cycle is always attracting, it suffices to show that, Solving the inequality for α yields the requirement that, We know by our assumptions that 0 < α < β.Therefore, it suffices to show that the left-hand-side of the inequality is always less than zero.To prove this, note that for all β > 0 we have, Thus multiplying the left and right-hand sides of Eq. (A1) yields, )︃ < 0.
Therefore, it follows that and for all 0 < α < β, ℓ 1 < 0 and thus the limit cycle is always attracting.

FIG. 3 .
FIG.3.The amplitude of the travelling wave solution decreases, approaching the fixed point solution ui(x, t) =1  3 as k increases.