Coevolutionary dynamics of a variant of the cyclic Lotka-Volterra model with three-agent interactions

TL;DR

This paper explores a variant of the cyclic Lotka-Volterra model with three-agent interactions, analyzing its bifurcation behavior, stochastic effects, and spatial dynamics, including the impact of mobility and nonlinear diffusion in ecological systems.

Contribution

It introduces a new model with three-agent interactions, analyzes bifurcations and stochastic effects, and studies spatial and mobility influences on cyclic dominance dynamics.

Findings

Degenerate Hopf bifurcation identified in well-mixed environments.

Stochastic noise influences extinction probabilities near bifurcation.

Spatial locality eliminates bifurcations, affecting stability and dynamics.

Abstract

We study a variant of the cyclic Lotka-Volterra model with three-agent interactions. Inspired by a multiplayer variation of the Rock-Paper-Scissors game, the model describes an ideal ecosystem in which cyclic competition among three species develops through cooperative predation. Its rate equations in a well-mixed environment display a degenerate Hopf bifurcation, occurring as reactions involving two predators plus one prey have the same rate as reactions involving two preys plus one predator. We estimate the magnitude of the stochastic noise at the bifurcation point, where finite size effects turn neutrally stable orbits into erratically diverging trajectories. In particular, we compare analytic predictions for the extinction probability, derived in the Fokker-Planck approximation, with numerical simulations based on the Gillespie stochastic algorithm. We then extend the analysis of the phase portrait to heterogeneous rates. In a well-mixed environment, we observe a continuum of degenerate Hopf bifurcations, generalizing the above one. Neutral stability ensues from a complex equilibrium between different reactions. Remarkably, on a two-dimensional lattice, all bifurcations disappear as a consequence of the spatial locality of the interactions. In the second part of the paper, we investigate the effects of mobility in a lattice metapopulation model with patches hosting several agents. We find that strategies propagate along the arms of rotating spirals, as they usually do in models of cyclic dominance. We observe propagation instabilities in the regime of large wavelengths. We also examine three-agent interactions inducing nonlinear diffusion.

Tables

Table 1. Table 1: A three-player variant of RPS .

outcome	payoffs
R R R, P P P, S S S	all players receive 0 points
R R P, P P S, S S R	the dominant player
	receives 1 point
R R S, P P R, S S P	each dominant player
	receives 1/2 point
R P S	all players receive 0 points

Table 2. Table 2: Transition rates and density variations.

reaction	$w$	$N\delta r$	$N\delta p$
R R P $\to$ R P P	$d_{e}N{r}^{2}p$	$-1$	$1$
P P S $\to$ P S S	$d_{e}N{p}^{2}s$	$0$	$-1$
S S R $\to$ S R R	$d_{e}N{s}^{2}r$	$1$	$0$
R R S $\to$ R R R	$d_{p}N{r}^{2}s$	$1$	$0$
P P R $\to$ P P P	$d_{p}N{p}^{2}r$	$-1$	$1$
S S P $\to$ S S S	$d_{p}N{s}^{2}p$	$0$	$-1$

Equations

$$\left\{\begin{array}[]{l}X_{1}X_{2}\dots X_{k}\to Y_{1}Y_{2}\dots Y_{\ell},\text{ for }k,\ell\geq 3\,,\\[4.30554pt] X_{i},Y_{j}\in\{{{{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}\tt R}},{{\color[rgb]{0,0.7,0.3}\definecolor[named]{pgfstrokecolor}{rgb}{0,0.7,0.3}\tt P}},{{\color[rgb]{0.15,0.15,0.85}\definecolor[named]{pgfstrokecolor}{rgb}{0.15,0.15,0.85}\tt S}},\text{{{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}$\emptyset$}}}}\}\end{array}\right.$$

$$\begin{array}[]{llcl}\text{{{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}\tt R}}\,{{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}\tt R}}\,{{\color[rgb]{0,0.7,0.3}\definecolor[named]{pgfstrokecolor}{rgb}{0,0.7,0.3}\tt P}}}&\to\ \ \text{{{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}\tt R}}\,{{\color[rgb]{0,0.7,0.3}\definecolor[named]{pgfstrokecolor}{rgb}{0,0.7,0.3}\tt P}}\,{{\color[rgb]{0,0.7,0.3}\definecolor[named]{pgfstrokecolor}{rgb}{0,0.7,0.3}\tt P}}}&\text{\ with rate}&d_{\scriptscriptstyle\rm RRP}\,,\\[3.01385pt] \text{{{\color[rgb]{0,0.7,0.3}\definecolor[named]{pgfstrokecolor}{rgb}{0,0.7,0.3}\tt P}}\,{{\color[rgb]{0,0.7,0.3}\definecolor[named]{pgfstrokecolor}{rgb}{0,0.7,0.3}\tt P}}\,{{\color[rgb]{0.15,0.15,0.85}\definecolor[named]{pgfstrokecolor}{rgb}{0.15,0.15,0.85}\tt S}}}&\to\ \ \text{{{\color[rgb]{0,0.7,0.3}\definecolor[named]{pgfstrokecolor}{rgb}{0,0.7,0.3}\tt P}}\,{{\color[rgb]{0.15,0.15,0.85}\definecolor[named]{pgfstrokecolor}{rgb}{0.15,0.15,0.85}\tt S}}\,{{\color[rgb]{0.15,0.15,0.85}\definecolor[named]{pgfstrokecolor}{rgb}{0.15,0.15,0.85}\tt S}}}&\text{\ \, "}&d_{\scriptscriptstyle\rm PPS}\,,\\[3.01385pt] \text{{{\color[rgb]{0.15,0.15,0.85}\definecolor[named]{pgfstrokecolor}{rgb}{0.15,0.15,0.85}\tt S}}\,{{\color[rgb]{0.15,0.15,0.85}\definecolor[named]{pgfstrokecolor}{rgb}{0.15,0.15,0.85}\tt S}}\,{{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}\tt R}}}&\to\ \ \text{{{\color[rgb]{0.15,0.15,0.85}\definecolor[named]{pgfstrokecolor}{rgb}{0.15,0.15,0.85}\tt S}}\,{{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}\tt R}}\,{{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}\tt R}}}&\text{\ \, "}&d_{\scriptscriptstyle\rm SSR}\,,\\[3.01385pt] \text{{{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}\tt R}}\,{{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}\tt R}}\,{{\color[rgb]{0.15,0.15,0.85}\definecolor[named]{pgfstrokecolor}{rgb}{0.15,0.15,0.85}\tt S}}}&\to\ \ \text{{{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}\tt R}}\,{{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}\tt R}}\,{{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}\tt R}}}&\text{\ \, "}&d_{\scriptscriptstyle\rm RRS}\,,\\[3.01385pt] \text{{{\color[rgb]{0,0.7,0.3}\definecolor[named]{pgfstrokecolor}{rgb}{0,0.7,0.3}\tt P}}\,{{\color[rgb]{0,0.7,0.3}\definecolor[named]{pgfstrokecolor}{rgb}{0,0.7,0.3}\tt P}}\,{{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}\tt R}}}&\to\ \ \text{{{\color[rgb]{0,0.7,0.3}\definecolor[named]{pgfstrokecolor}{rgb}{0,0.7,0.3}\tt P}}\,{{\color[rgb]{0,0.7,0.3}\definecolor[named]{pgfstrokecolor}{rgb}{0,0.7,0.3}\tt P}}\,{{\color[rgb]{0,0.7,0.3}\definecolor[named]{pgfstrokecolor}{rgb}{0,0.7,0.3}\tt P}}}&\text{\ \, "}&d_{\scriptscriptstyle\rm PPR}\,,\\[3.01385pt] \text{{{\color[rgb]{0.15,0.15,0.85}\definecolor[named]{pgfstrokecolor}{rgb}{0.15,0.15,0.85}\tt S}}\,{{\color[rgb]{0.15,0.15,0.85}\definecolor[named]{pgfstrokecolor}{rgb}{0.15,0.15,0.85}\tt S}}\,{{\color[rgb]{0,0.7,0.3}\definecolor[named]{pgfstrokecolor}{rgb}{0,0.7,0.3}\tt P}}}&\to\ \ \text{{{\color[rgb]{0.15,0.15,0.85}\definecolor[named]{pgfstrokecolor}{rgb}{0.15,0.15,0.85}\tt S}}\,{{\color[rgb]{0.15,0.15,0.85}\definecolor[named]{pgfstrokecolor}{rgb}{0.15,0.15,0.85}\tt S}}\,{{\color[rgb]{0.15,0.15,0.85}\definecolor[named]{pgfstrokecolor}{rgb}{0.15,0.15,0.85}\tt S}}}&\text{\ \, "}&d_{\scriptscriptstyle\rm SSP}\,,\end{array}$$

$$\displaystyle\hskip-19.63246pt\begin{array}[]{l}\dot{r}={\cal F}_{\rm r}(r,p,s)\,,\ \dot{p}={\cal F}_{\rm p}(r,p,s)\,,\ \dot{s}={\cal F}_{\rm s}(r,p,s)\,;\\[8.61108pt] {\cal F}_{\rm r}=r^{2}\left(d_{\scriptscriptstyle\rm RRS}s-d_{\scriptscriptstyle\rm RRP}p\right)+r\left(d_{\scriptscriptstyle\rm SSR}s^{2}-d_{\scriptscriptstyle\rm PPR}p^{2}\right),\\[4.30554pt] {\cal F}_{\rm p}=p^{2}\left(d_{\scriptscriptstyle\rm PPR}r-d_{\scriptscriptstyle\rm PPS}s\right)+p\left(d_{\scriptscriptstyle\rm RRP}r^{2}-d_{\scriptscriptstyle\rm SSP}s^{2}\right),\\[4.30554pt] {\cal F}_{\rm s}=s^{2}\left(d_{\scriptscriptstyle\rm SSP}p-d_{\scriptscriptstyle\rm SSR}r\right)+s\left(d_{\scriptscriptstyle\rm PPS}p^{2}-d_{\scriptscriptstyle\rm RRS}r^{2}\right).\end{array}$$

$$\left.\begin{array}[]{l}d_{\scriptscriptstyle\rm RRP}=d_{\scriptscriptstyle\rm PPS}=d_{\scriptscriptstyle\rm SSR}\equiv d_{\rm e}\,,\\[4.30554pt] d_{\scriptscriptstyle\rm RRS}=d_{\scriptscriptstyle\rm PPR}=d_{\scriptscriptstyle\rm SSP}\equiv d_{\rm p}\,.\end{array}\right.$$

$$\displaystyle\left.\begin{array}[]{l}\displaystyle{\dot{x}_{\text{r}}=-\frac{d_{\rm e}}{3}x_{\text{r}}-\frac{d_{\rm e}+d_{\rm p}}{3}x_{\text{p}}}\,,\\[8.61108pt] \displaystyle{\dot{x}_{\text{p}}=\ \,\frac{d_{\rm p}}{3}x_{\text{p}}+\frac{d_{\rm e}+d_{\rm p}}{3}x_{\text{r}}}\,.\end{array}\right.$$

$$\left.\begin{array}[]{l}\displaystyle{\lambda=\frac{1}{6}(d_{\rm p}-d_{\rm e})+\frac{\text{i}}{2\sqrt{3}}(d_{\rm p}+d_{\rm e})}\,,\\[8.61108pt] \displaystyle{\bar{\lambda}=\frac{1}{6}(d_{\rm p}-d_{\rm e})-\frac{\text{i}}{2\sqrt{3}}(d_{\rm p}+d_{\rm e})}\,.\end{array}\right.$$

$$S=\left(\begin{array}[]{cc}1/\sqrt{3}&-1/\sqrt{3}\\[4.30554pt] 1&1\end{array}\right)\,.$$

$$\left.\begin{array}[]{l}\dot{y}_{\text{r}}=\text{Re}(\lambda)y_{\text{r}}-\text{Im}(\lambda)y_{\text{p}}+f_{\text{r}}(y_{\text{r}},y_{\text{p}},\epsilon)\,,\\[8.61108pt] \dot{y}_{\text{p}}=\text{Im}(\lambda)y_{\text{r}}+\text{Re}(\lambda)y_{\text{p}}+f_{\text{p}}(y_{\text{r}},y_{\text{p}},\epsilon)\,,\end{array}\right.$$

$$f_{\text{r}}(y_{\text{r}},y_{\text{p}},ϵ)$$

$$-\frac{\sqrt{3}}{4}(1+ϵ)y_{\text{p}}^2-\frac{3ϵ}{4}{y}_{\text{r}}^3-\frac{3ϵ}{4}{y}_{\text{r}}{y}_{\text{p}}^2,$$

$$f_{\text{p}}(y_{\text{r}},y_{\text{p}},ϵ)$$

$$-\frac{ϵ}{4}{y}_{\text{p}}^2-\frac{3ϵ}{4}{y}_{\text{p}}^3-\frac{3ϵ}{4}{y}_{\text{p}}{y}_{\text{r}}^2.$$

$$\dot{z}=\lambda z+f(z,\bar{z}),$$

$$f(z,\bar{z})$$

$$\equiv R{\bar{z}}^2-S\bar{z}{z}^2.$$

$$\dot{z}=\lambda z-c(ϵ)z|z{|}^2,$$

$$a(ϵ)=\frac{3}{2}\frac{4{ϵ}^2+15ϵ+15}{7{ϵ}^2+27ϵ+27}ϵ,$$

$$b(ϵ)=\frac{9\sqrt{3}}{4}\frac{(ϵ+2)({ϵ}^2+3ϵ+3)}{7{ϵ}^2+27ϵ+27}.$$

$$\dot{\rho }=\frac{ϵ}{6}\rho -a(ϵ){\rho }^3,\dot{\theta }=\frac{2+ϵ}{2\sqrt{3}}-b(ϵ){\rho }^2.$$

$${\partial }_t\mathcal{P}(\varphi ,t)$$

$$-w(\varphi \to \varphi +\delta \varphi )\mathcal{P}(\varphi ,t)\},$$

$${\partial }_t\mathcal{P}(\varphi ,t)$$

$$+\frac{1}{2}\sum _{i,j=\text{r},\text{p}}{\partial }_i{\partial }_j\left[B_{ij}(\varphi )\mathcal{P}(\varphi ,t)\right].$$

$${\alpha }_i(\varphi )$$

$$B_{ij}(\varphi )$$

$${\alpha }_{\text{r}}(\varphi )=r(s-p),$$

$${\alpha }_{\text{p}}(\varphi )=p(r-s),$$

$$B_{\text{rr}}(\varphi )=\frac{1}{N}r(p+s-2ps),$$

$$B_{\text{rp}}(\varphi )=B_{\text{pr}}(\varphi )=-\frac{1}{N}rp(r+p),$$

$$B_{\text{pp}}(\varphi )=\frac{1}{N}p(r+s-2rs).$$

$${\partial }_t\mathcal{P}(x,t)$$

$$+\frac{1}{2}{\mathcal{B}}_{ij}{\partial }_i{\partial }_j\mathcal{P}(x,t),$$

$${\cal A}=-\frac{1}{3}\left(\begin{array}[]{rr}\!\!1&2\\ \!\!-2&1\end{array}\right)\,,\quad{\cal B}=\frac{2}{27N}\left(\begin{array}[]{rr}\!\!2&-1\\ -1&2\end{array}\right)\,.$$

$$\displaystyle{\cal A}\to\tilde{\cal A}={\cal S}{\cal A}{\cal S}^{-1}=\omega_{0}\left(\begin{array}[]{rr}0&-1\\ -1&0\end{array}\right)\,,$$

$$\displaystyle{\cal B}\to\tilde{\cal B}={\cal S}{\cal B}{\cal S}^{\text{T}}=\frac{1}{9N}\left(\begin{array}[]{rr}1&0\\ 0&1\end{array}\right)\,.$$

$${\partial }_t\mathcal{P}(y,t)$$

$$+\frac{1}{18N}[{\partial }_{\text{r}}^2+{\partial }_{\text{p}}^2]\mathcal{P}(y,t).$$

Full text

“Three at play. That’ll be the day!”

(a child in Wings of desire [W. Wenders, 1987])

I Introduction

Cyclic competition is distinctively associated with closed relational chains, describing aspects of the struggle for life such as feeding, hunting, mating. The prototype is a system made of three different species, interacting with one another like children playing Rock-Paper-Scissors (RPS), the famous game where paper (P) wraps rock, scissors (S) cut paper and rock (R) crushes scissors, see Fig. 1. Such schemes are at heart of biological systems spanning a wide range of length scales and complexity. Examples include the repressilator Elowitz1 , the E. coli colicin E2 system Kerr1 , the common side-blotched lizard Sinervo1 , several plant systems Taylor1 ; Cameron1 ; Lankau1 and so forth.

In principle, the absence of apex predators and bottom prey in cyclic chains allows species to dominate in turn. As soon as one of them outperforms the others, it becomes itself a source of nourishment for the next one along the chain. This feature suggests that cyclic competition may serve as a fundamental mechanism facilitating coexistence and biodiversity. Experiments performed in Ref. Kerr1 on three cyclically interacting strains of E. coli confirmed this thesis and also made it more precise: species can coexist provided ecological processes (interaction and dispersal) develop locally. In practice, to ensure coexistence in the experiments, colonies of different strains had to first grow in separate spatial domains and then be left free to invade neighboring colonies.

Theoretical insight into the role of locality and individual mobility for coexistence was achieved in a ground-breaking paper Reichenbach1 thanks to an evolutionary model based on a three-state RPS game. Results were then generalized to a four-state variant of the model Peltomaki1 . Models unifying pairwise cyclic dominance with removal (as in Ref. Reichenbach1 ) and replacement (as in Ref. Peltomaki1 ) were later studied in Refs. Szczesny1 ; Szczesny2 ; Rulands1 . These studies revealed that:

• $\circ$

individual mobility promotes coexistence by inducing self-organization of the strategies into spiral waves, traveling across the environment;

• $\circ$

stochastic noise, arising in finite populations, produces local inhomogeneities; nevertheless, it cannot prevent the creation of spiral waves;

• $\circ$

beyond a critical threshold, individual mobility leads to species extinction.

Subsequent research dealt with a plethora of either induced or independent issues. They include the effects of competition on pattern formation Jiang1 , the observation of multi-armed spirals Jiang2 , the emergence of convective instability Reichenbach2 , the role of stochastic noise Reichenbach3 and mutations Szczesny1 ; Szczesny2 , the analysis of coexistence and extinction basins Rulands1 ; Shi1 ; Ni1 , the effects of uniform and nonuniform intra-specific competition Yang1 ; Park1 , the influence of directional mobility on coexistence Avelino1 and so forth. We omit to mention other relevant developments only for the sake of conciseness, while we refer the reader to Ref. Szolnoki1 for a systematic review of concepts and results.

A common assumption in all the above-mentioned studies is that predation is a pairwise interaction, involving a single predator and a single prey. While this is the case for many living organisms, either animals or microbes, it is not so for others. An alternative strategy, favored by natural selection, is cooperative predation. Mammals such as wolves Mech1 , chimpanzees Stanford1 , dolphins Samuels1 and lions Schaller1 cooperate in hunting. Some insects, including ants Holldobler1 , behave analogously. Even bacteria can practice group hunting. Among them we mention Saprospira Perez1 , Myxococcus Xanthus Berleman1 ; Munoz1 and Lysobacter Jurkevitch1 . They essentially “require a quorum of predators to degrade the prey, using excreted hydrolytic enzymes” Jurkevitch1 .

Another recently esablished mechanism for species coexistence involving multiagent interactions in microbial communities is the protection spillover. Kelsic et al have shown in Ref. Kelsic that bacteria protecting themselves from potential predators by excreting inhibitory enzyms, can protect bacteria of other species in the nearby area from the same predators. If we rephrase the the RPS chain as “P survives an attack by R, which survives an attack by S, which in turn survives an attack by P”, then protection spillover means that R helps P survive an attack by S. This protection mechanism stabilizes the species in a well-mixed environment Bergstrom and gives rise, in structured populations, to spatial patterns which cannot be observed in the classical RPS game Szolnoki0 .

Other spatially structured models of cyclic competition featuring group predation were considered in Refs. Szolnoki2 ; Cheng1 ; Cazaubiel1 ; Lutz1 . Specifically, the authors of Ref. Szolnoki2 studied two variants of a three-state lattice model in which predation entails simultaneous pairwise interactions. One of them assumes that agents interact with their four von-Neumann neighbors, the other with their eight Moore neighbors Toffoli1 . They found that increasing the interaction range can decelerate the propagation of predators and even revert the direction of species invasion contrary to its natural definition. The authors of Ref. Cheng1 observed the emergence of mesoscopic subgroups of coexisting species in a five-state lattice model based on the Rock-Paper-Scissors-Lizard-Spock game. They developed a mean-field theory to show that group interactions at the mesoscopic scale must be taken into account to justify the observed states of coexistence. The authors of Ref. Cazaubiel1 studied a stochastic lattice version of a model introduced by Lett, Auger, Gaillard Lett1 . In this model, the abundances of prey and predators are constant, while the fractions of each population using either an individual or a collective strategy coevolve. Their results include a complex phase diagram in which four different strategies (corresponding to prey/predators behaving collectively/individually in all possible combinations) dominate or coexist. Remarkably, in the pure coexistence phase, they found cyclic dominance of the four strategies. Finally, the same authors further explored the spatial version of the model to quantify some geometrical and percolative properties of the clusters formed by the four strategies Lutz1 .

The reader will recognize that none of the abovementioned approaches considers ab initio multiagent microscopic reactions with independent rates, such as

[TABLE]

(in the literature $\emptyset$ conventionally denotes an empty site)111Reactions considered in Ref. Szolnoki2 can be in fact recast in the form of Eq. (I.1) for $k=\ell=5,9$. Yet, their rates are somewhat intertwined in that they are functions of the pairwise rates $\delta_{0},\delta_{1}$.. In this paper, we study a simple variant of the cyclic Lotka-Volterra model Reichenbach5 ; Frachebourg1 ; Frachebourg2 ; Frachebourg3 ; Provata1 ; Tsekouras1 ; Szabo1 ; Szabo2 with three-agent interactions, like Eq. (I.1) for $k=\ell=3$. To explore the model, we use mathematical techniques developed in similar contexts, including nonlinear analysis of bifurcations, stochastic partial differential equations and numerical simulations.

We find that the underlying rate equations in a well-mixed environment exhibit a reactive fixed point falling in the universality class of the Hopf bifurcations. As such, they induce a macroscopic phenomenology qualitatively similar to models already studied in the literature. Nonetheless, the internal dynamics of the model is original and interesting. Two opposing forces drive the system. Interactions involving two prey plus one predator pull it towards the reactive fixed point, thus playing an equilibrating role. By contrast, interactions involving two predators plus one prey push the system away from the reactive fixed point, thus producing a polarizing effect. The relative strength of the two forces controls the evolution of the system in a spatially structured version of the model, with patches hosting several motile agents (metapopulation model). When polarizing interactions dominate over equilibrating ones, coexistence is achieved through the development of spatiotemporal patterns, taking the usual form of rotating spiral waves. Similar to Refs. Szolnoki10 ; Szolnoki11 ; Szolnoki12 ; Szolnoki13 , the dynamics of three-agent interactions on a spatially structured population shows a largely different behavior from predictions based on well-mixed calculations.

We also find that spatial topology is critical for the evolution of species. Indeed, the phase portrait of the model changes drastically on a two-dimensional lattice with single agent per site and nearest-neighbor interactions. Locality makes the reactive fixed point stable for every choice of reaction rates. The disappearance of Hopf bifurcations indicates that local patches hosting several agents are essential for the development of patterns.

The paper is organized as follows. First of all, in section II, we define the model and study its equations in a well-mixed environment for homogeneous rates. In section III, we estimate the magnitude of the stochastic noise affecting finite agent populations. As known Reichenbach5 , fluctuations turn neutrally stable orbits into diverging trajectories, eventually resulting in the extinction of two species. We compare the extinction probability, derived in the Fokker-Planck approximation, with numerical simulations. In section IV, we discuss the phase portrait of the model for heterogeneous rates. We also compare results in a well-mixed environment with those on a two-dimensional lattice. In section V, we introduce mobility reactions in a lattice metapopulation model with patches hosting several agents. In section VI, we examine three-agent interactions inducing nonlinear diffusion. Finally, in section VII, we draw conclusions.

II Rate equations

Some time ago we happened to observe three children playing a three-player variant of RPS in the hall of their school. Intrigued by the game, we asked them about it. In order not to leave anyone out—they proudly explained—they were playing all at once. The rules of the game were as follows. On each round, children had to deliver simultaneously one of the usual hand signals, representing R, P and S. Round by round they received payoffs, based on cyclic dominance and depending on the combination of delivered signals, as we report in Table 1222While writing the paper we realized that essentially the same three-player variant of RPS is described in Ref. Walker1 . Similar variants can be found on various websites. In a scene of Sonatine (ソナチネ), a 1993 film by T. Kitano, three yakuza gangsters play a three-player variant of RPS on the beach..

To translate the game into the language of evolutionary game theory, we assume a population of $N\gg 1$ agents, each adopting one of the competing strategies. We let $r$, $p$, $s$ denote the relative abundances of R, P, S respectively. As such they fulfill $r+p+s=1$. In our model, densities evolve in time as a consequence of microscopic interactions inspired by Table 1. More precisely, in place of payoffs we consider stochastic reactions mediated by a dominance-replacement mechanism, namely

[TABLE]

where rates (transition probabilities per unit time) are independent of one another. For the time being, we leave out transitions in which agents carry three different strategies before interacting. We also leave out transitions in which they all carry the same strategy, since such configurations yield no payoff in the classic formulation of the game. Rate equations (RE) including all contributions listed in Eqs. (II.1) read as

[TABLE]

Altogether, these cubic equations depend on six parameters. We can absorb one of them into a redefinition of time, thus obtaining five effective parameters. Moreover, we can drop the equation for $\dot{s}$ provided we insert $s=1-r-p$ into those for $\dot{r}$ and $\dot{p}$. Eqs. (II.6) describe correctly the evolution of strategies, induced by Eqs. (II.1), in a well-mixed environment as $N\to\infty$. They represent a variant of the cyclic Lotka-Volterra equations Frachebourg1 ; Frachebourg2 ; Frachebourg3 ; Provata1 ; Tsekouras1 ; Szabo1 ; Szabo2 ; Reichenbach1 encompassing three-agent interactions.

To study the model, we first consider a simplified version in which rates are homogeneous. More precisely, we let

[TABLE]

As can be seen, homogeneity is not complete. Eq. (II.7) follows from separating reactions into two disjoint sets. The first three reactions in Eq. (II.1) correspond to the first line of Eq. (II.7). They are functionally homogeneous in that they involve two prey and one predator in the initial state. Their occurrence produces a local change of majority. The last three reactions in Eq. (II.1) correspond to the second line of Eq. (II.7). They start with two predators and one prey. Their occurrence results in a further increase of the local majority. The two groups of reactions play antagonistic roles. When the system is in a macroscopic state in which a strategy dominates, the former contributes to equilibrating it, whereas the latter contributes to further polarizing it.

Resulting RE have three absorbing fixed points, $(r_{1},p_{1})=(1,0)$, $(r_{2},p_{2})=(0,1)$, $(r_{3},p_{3})=(0,0)$ and a reactive one, $(r_{*},p_{*})=(1/3,1/3)$. Since we are interested in the behavior of the system near the latter point, we let $r=1/3+x_{\text{r}}$ and $p=1/3+x_{\text{p}}$. Upon expanding the RE in Taylor series to first order around $x_{\text{r}}=x_{\text{p}}=0$, we obtain

[TABLE]

The eigenvalues of the Jacobian matrix are

[TABLE]

The real part vanishes for $d_{\rm e}=d_{\rm p}$, i.e. when equilibrating and polarizing forces compensate exactly. For this choice of rates we have a Hopf bifurcation. Starting near $(r_{*},p_{*})$, the system spirals inwards for $d_{\rm e}>d_{\rm p}$, whereas it spirals outwards for $d_{\rm e}<d_{\rm p}$. For $d_{\rm e}=d_{\rm p}(=1)$ it travels on neutrally stable orbits enclosing the fixed point. They have angular frequency $\omega_{0}=1/\sqrt{3}$ and first integral $rps=\text{constant}$.

To determine the type of the Hopf bifurcation, we need to bring the RE to normal form according to a standard procedure, detailed, e.g., in Refs. Wiggins1 ; Strogatz1 . To this aim, we let $d_{\rm e}=1$ and $d_{\rm p}=1+\epsilon$. We perform a linear change of variables, namely we let $y=Sx$ with

[TABLE]

It yields equations

[TABLE]

where the nonlinear functions $f_{\text{r}}$ and $f_{\text{p}}$ read as

[TABLE]

Next, we perform a change to complex variables $z=y_{\text{r}}+\text{i}y_{\text{p}}$ and ${\bar{z}}=y_{\text{r}}-\text{i}y_{\text{p}}$. Inverse relations are given by $y_{r}(z,\bar{z})=(z+\bar{z})/2$ and $y_{p}(z,\bar{z})=(z-\bar{z})/2\text{i}$. Inserting them into the above formulas yields

[TABLE]

with

[TABLE]

As can be seen, $f$ has quadratic and cubic terms. We can remove the former by performing an additional change of variable, viz. $z\to z+A\bar{z}^{2}$. We choose $A$ such that $(\lambda-2\bar{\lambda})A=R$. After some algebra333Recall that $y=x+Ax^{2}$ cannot be inverted exactly. By iteration we have $x=y-Ax^{2}=y-A(y-Ax^{2})^{2}=y-Ay^{2}+2A^{2}y^{3}+\text{O}(y^{4})\,$. we arrive at

[TABLE]

with $c(\epsilon)=a(\epsilon)+\text{i}b(\epsilon)$ and

[TABLE]

Since $a(\epsilon)\approx(5/6)\,\epsilon$ as $\epsilon\to 0$, we conclude that the Hopf bifurcation is degenerate. In polar coordinates $z=\rho\text{e}^{\text{i}\theta}$ Eq. (II.18) turns into

[TABLE]

In Fig. 2, we compare by numerical integration Eq. (II.21) with the full RE for $\epsilon=0.1$. In both plots, dashed and dotted lines correspond respectively to the leading order approximation (LO) and the full next-to-leading order (NLO) version of Eq. (II.21). We always choose initial conditions with densities lying near the reactive fixed point. We observe no limit cycle. The asymptotic saturation of $\rho$ in Fig. 2 (a) corresponds to heteroclinic cycles approaching the absorbing fixed points.

A Hopf bifurcation of the same type as Eqs. (II.18)-(II.20) characterizes the dynamics of the May-Leonard model May1 . In that model, the bifurcation point corresponds to the vanishing of dominance-removal reactions. Therefore, the May-Leonard model is well-defined only on one side of the bifurcation. Ours makes sense on both sides. More importantly, on the unstable side, both models feature heteroclinic cycles. Owing to this similarity, their spatially structured versions, discussed respectively in Ref. Reichenbach2 and section V, are phenomenologically equivalent.

III Stochastic noise

Intuition suggests that multiagent interactions should produce larger stochastic fluctuations than pairwise ones since they involve more fluctuating degrees of freedom. A non-trivial and surprising consequence is that strategies have to fight longer before one of them prevails on the others. As a result, the probability of coexistence increases. Using the Fokker-Planck equation, we derive an accurate estimate of the magnitude of the stochastic noise in our model. We follow Ref. Reichenbach5 , where calculations are fully detailed. The idea is to consider a specific setting in which the RE predict neutrally stable orbits, namely $d_{\rm e}=d_{\rm p}=1$ in our case. For $N<\infty$ the conservation law $rps=\text{const}.$ is broken by $\text{O}(N^{-1})$ terms. Hence, the system follows an erratic trajectory, interpolating between different neutrally stable orbits. Eventually, it ends up on one of the absorbing fixed points. We can derive the intrinsic magnitude of the stochastic noise from the Master Equation

[TABLE]

where $\phi=(r,p)$ and $w(\phi\to\phi^{\prime})$ denotes the transition probability per unit time (rate) from $\phi$ to $\phi^{\prime}$. The sum over $\delta\phi$ includes all possible microscopic transitions characterizing the model. We choose conventionally the macroscopic time unit as the interval including $N$ reactions. In Table 2, we list rates and density variations corresponding to this choice (one has to replace $s=1-r-p$ in all expressions).

Eq. (III.1) yields an exact description for $N<\infty$. Unfortunately, we cannot solve it analytically. We can obtain an effective approximation by means of the Kramers–Moyal expansion, which, upon truncation to second order, yields the Fokker-Planck equation (FPE)

[TABLE]

The functions $\alpha_{i}$ and $B_{ij}$, respectively known as drift and diffusion functions, are given by

[TABLE]

From Table 2 it follows that

[TABLE]

[TABLE]

Just like in section II, we let $\phi=(1/3,1/3)+x$, then we expand the whole FPE around $x=0$ (Van Kampen’s linear noise approximation VanKampen1 ). Accordingly, we obtain

[TABLE]

with

[TABLE]

As a final step, we perform another change of variables, namely we let $x\to y={\cal S}x$, with ${\cal S}=\frac{\sqrt{3}}{2}\left(\begin{smallmatrix}2\omega_{0}&\omega_{0}\\ 0&1\end{smallmatrix}\right)$. Accordingly, ${\cal A}$ and ${\cal B}$ turn into

[TABLE]

The diffusion matrix $\tilde{\cal B}$ is now diagonal. Hence, the FPE takes the simplified form

[TABLE]

The diffusion constant $D=1/(18N)$ is smaller than the analogous constant in the original cyclic Lotka-Volterra model by a factor of 2/3, whereas $\omega_{0}$ is the same, see Eq. (25) of Ref. Reichenbach5 444The drift term in Eq. (III.18) has an overall minus sign with respect to Eq. (25) of Ref. Reichenbach5 . It is due to the opposite invasion order: Ref. Reichenbach5 assumes $a\to b\to c\to a$, while here we have {{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}\tt R}}\leftarrow{{\color[rgb]{0,0.7,0.3}\definecolor[named]{pgfstrokecolor}{rgb}{0,0.7,0.3}\tt P}}\leftarrow{{\color[rgb]{0.15,0.15,0.85}\definecolor[named]{pgfstrokecolor}{rgb}{0.15,0.15,0.85}\tt S}}\leftarrow{{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}\tt R}}.. The reduction factor is rather suggestive, as it equals the ratio of simultaneous players in the two models. The value of 2/3 can be easily explained. The original model has three microscopic reactions, ours has twice this number. Each reaction contributes positively to the stochastic noise. This yields a multiplicative factor of 2 in the diffusion matrix. Moreover, the contribution from each reaction is quadratic in the original model while cubic in ours. In the linear noise approximation, we calculate the diffusion matrix at the reactive fixed point. This yields an additional multiplicative factor of 1/3 and that is all. In general, the more agents partake in the interactions, the larger the number of possible microscopic reactions. The magnitude of the induced stochastic noise depends eventually on combinatorial factors, including the number and degree of reactions. Of course, it depends as well on the value of the strategy densities at equilibrium.

We can calculate the probability ${\cal P}_{\rm ext}$ that two species go extinct after a certain time in the Fokker-Planck approximation as the authors of Ref. Reichenbach5 do. There is no need to repeat the derivation here, since it applies identically to our model. The final formula is

[TABLE]

where LT stands for Laplace Transform, $u=t/N$ is a scaling variable measuring time in units of $N$ and $I_{0}$ denotes the Bessel function of first kind and 0th order. $R$ represents the distance traveled by the system on its random walk from the reactive fixed point to one of the absorbing fixed points. We can regard it as the radius of an absorbing sphere. The authors of Ref. Reichenbach5 adopt three possible definitions of $R$, namely $R_{0}=1/3$, $R_{1}=1/\sqrt{3}$ and $R_{2}=(R_{0}+R_{1})/2$. They yield three different probability functions. Fig. 3 shows them together with the results of numerical simulations based on Gillespie’s algorithm Gillespie1 ; Gillespie2 . To compute the extinction probability from Eq. (III.19) we use a numerical implementation of the inverse LT. Moreover, we expand $I_{0}$ asymptotically to the 10th order, as also the authors of Ref. Reichenbach5 do. To simulate extinction times correctly, we need to rescale all rates by appropriate volume factors. Notice that all reactions in our model involve two reactants of the same species and one of another. For such reactions, the right definition of the reaction parameters $c_{\text{XXY}}$ in Gillespie’s algorithm is $c_{\text{XXY}}=2d_{\text{XXY}}/N^{2}$, for \text{X,Y}\in\{{{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}\tt R}},{{\color[rgb]{0,0.7,0.3}\definecolor[named]{pgfstrokecolor}{rgb}{0,0.7,0.3}\tt P}},{{\color[rgb]{0.15,0.15,0.85}\definecolor[named]{pgfstrokecolor}{rgb}{0.15,0.15,0.85}\tt S}}\} Gillespie1 ; Gillespie2 . Similar to the original Lotka-Volterra model, the analytic prediction that best fits numerical data is the one corresponding to $R_{2}$.

As anticipated, the extinction probability is uniformly lower in the cyclic Lotka-Volterra model with three-agent interactions than with two-agent ones, although the former are intrinsically noisier than the latter. This result is only apparently counterintuitive. In fact, it has a simple interpretation. When three agents interact, strategies fluctuate longer around the reactive fixed point before one of them prevails on the others. Group interactions help the system stay in equilibrium. Hence, we conclude, they promote species coexistence.

IV Heterogeneous rates

The reactive fixed point is symmetric for homogeneous rates, independently of whether $d_{\rm e}>d_{\rm p}$ or $d_{\rm e}<d_{\rm p}$. Things become more interesting as soon as we break homogeneity. Unfortunately, studying the model in full generality is not simple, due to the high dimensionality of the parameter space. A reasonable compromise is to let $d_{\scriptscriptstyle\rm RRP}=d_{\scriptscriptstyle\rm PPS}=d_{\scriptscriptstyle\rm SSR}\equiv d_{\rm e}=1$ and to leave all other rates unconstrained. Under this assumption, the RE still exhibit four fixed points, three absorbing plus one reactive. The former coincide with the vertices of the ternary diagram, as can be easily checked from Eqs. (II.6). The latter has a complex algebraic structure. Indeed, it reads

[TABLE]

with $X$ fulfilling the cubic equation

[TABLE]

We wish to examine how strategies rank at equilibrium depending on the polarizing rates. To this aim, we proceed like the authors of Ref. Szolnoki2 . Specifically, we denote by ($\tt A$,$\tt B$,$\tt C$) a domain in parameter space for which $a_{*}\leq b_{*}\leq c_{*}$, for $a_{*},b_{*},c_{*}$ a permutation of $r_{*},p_{*},s_{*}\equiv 1-r_{*}-p_{*}$ and $\tt A$, $\tt B$, $\tt C$ the corresponding permutation of R, P, S. In accordance with Ref. Szolnoki2 , we adopt the name of phase for all such domains, although this word is somewhat misleading (at least in our case). Indeed, both the equilibrium densities and their derivatives are smooth functions. In particular, they do not show discontinuities at phase transitions. In Fig. 4, we report the phase structure for a sequence of values of $d_{\scriptscriptstyle\rm RRS}$, ranging from 0.2 to 2.1 in steps of 0.1. In each plot, we notice a white region and a variously colored one. The former corresponds to unstable spirals, diverging from the reactive fixed point (they eventually turn into heteroclinic cycles), whereas the latter corresponds to stable spirals, converging to it. A continuum of Hopf bifurcations lies along the lines separating colored and white regions, as we discuss below. Colored regions split into six phases or less, depending on $d_{\scriptscriptstyle\rm RRS}$. Each phase makes contact with all others for $d_{\scriptscriptstyle\rm PPR}=d_{\scriptscriptstyle\rm SSP}=d_{\scriptscriptstyle\rm RRS}$. The contact point lies in the stable region only for $d_{\scriptscriptstyle\rm RRS}\leq 1$. The overall surface occupied by stable phases reduces as $d_{\scriptscriptstyle\rm RRS}$ increases (the larger $d_{\scriptscriptstyle\rm RRS}$ the stronger the polarizing force corresponding to fixed values of $d_{\scriptscriptstyle\rm PPR}$ and $d_{\scriptscriptstyle\rm SSP}$).

A glance to plot (i), corresponding to $d_{\scriptscriptstyle\rm RRS}=1$, suggests a recipe for predicting the phase of the system for given values of $d_{\scriptscriptstyle\rm RRS},d_{\scriptscriptstyle\rm PPR},d_{\scriptscriptstyle\rm SSP}$. It consists of three steps:

1. $1.$

rank the polarizing rates in ascending order (e.g. $d_{\scriptscriptstyle\rm PPR}<d_{\scriptscriptstyle\rm RRS}<d_{\scriptscriptstyle\rm SSP}$); 2. $2.$

extract the losing strategy from each rate label (it yields ({{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}\tt R}},{{\color[rgb]{0.15,0.15,0.85}\definecolor[named]{pgfstrokecolor}{rgb}{0.15,0.15,0.85}\tt S}},{{\color[rgb]{0,0.7,0.3}\definecolor[named]{pgfstrokecolor}{rgb}{0,0.7,0.3}\tt P}}) in the above example); 3. $3.$

turn each strategy into the immediately inferior one (it finally yields ({{\color[rgb]{0.15,0.15,0.85}\definecolor[named]{pgfstrokecolor}{rgb}{0.15,0.15,0.85}\tt S}},{{\color[rgb]{0,0.7,0.3}\definecolor[named]{pgfstrokecolor}{rgb}{0,0.7,0.3}\tt P}},{{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}\tt R}})).

The recipe has a straightforward interpretation: the more a species is preyed on by its predator, the more its prey has room to develop. Its main drawback is that it is only approximately exact due to nonlinear effects. For instance, the transition line separating ({{\color[rgb]{0.15,0.15,0.85}\definecolor[named]{pgfstrokecolor}{rgb}{0.15,0.15,0.85}\tt S}},{{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}\tt R}},{{\color[rgb]{0,0.7,0.3}\definecolor[named]{pgfstrokecolor}{rgb}{0,0.7,0.3}\tt P}}) from ({{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}\tt R}},{{\color[rgb]{0.15,0.15,0.85}\definecolor[named]{pgfstrokecolor}{rgb}{0.15,0.15,0.85}\tt S}},{{\color[rgb]{0,0.7,0.3}\definecolor[named]{pgfstrokecolor}{rgb}{0,0.7,0.3}\tt P}}) should bisect the phase plane for $d_{\scriptscriptstyle\rm RRS}=1$, while it does not. Moreover, all other transition lines develop a slope for $d_{\scriptscriptstyle\rm RRS}\neq 1$. We notice that by no reason the phases of Fig. 4 should be symmetric for $d_{\scriptscriptstyle\rm SSP}\leftrightarrow d_{\scriptscriptstyle\rm PPR}$. The only symmetry of the system is the cyclic one, namely {{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}\tt R}}\leftarrow{{\color[rgb]{0,0.7,0.3}\definecolor[named]{pgfstrokecolor}{rgb}{0,0.7,0.3}\tt P}}\leftarrow{{\color[rgb]{0.15,0.15,0.85}\definecolor[named]{pgfstrokecolor}{rgb}{0.15,0.15,0.85}\tt S}}\leftarrow{{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}\tt R}}. If we fix $d_{\scriptscriptstyle\rm RRS}$, we lose it.

Altogether, the phase structure looks more complex than observed in Ref. Szolnoki2 , although our model features simpler multiagent interactions. As far as we understand, the rationale behind this is that we consider fully independent rates, while the authors of Ref. Szolnoki2 build multiagent rates as accumulated payoffs, depending on the pairwise rates. As a consequence, they explore a two-dimensional manifold embedded in a larger and more complex phase space.

We can provide a better characterization of the transition lines separating white and colored regions. The Jacobian of the linearized RE has two complex conjugate eigenvalues $\lambda_{\text{r}}\pm\text{i}\lambda_{\text{i}}$. We can express both the real and the imaginary part as functions of the fixed-point densities, namely

[TABLE]

[TABLE]

[TABLE]

with coefficients reading

[TABLE]

Both $(r_{*},p_{*})$ and $\lambda_{\text{r}}$ depend nonlinearly upon $d_{\scriptscriptstyle\rm RRS}\,,d_{\scriptscriptstyle\rm PPR}\,,d_{\scriptscriptstyle\rm SSP}$. Any point, belonging to a transition line separating white and colored regions in Fig. 4, corresponds by definition to a reactive fixed point yielding $\lambda_{\text{r}}=0$. Keeping $(r_{*},p_{*})$ fixed means imposing two constraints. Since we have three degrees of freedom, we remain with one. We conclude that there is a one-dimensional (non-planar) manifold, characterized by $(r_{*},p_{*})=\text{const.}$, crossing the transition line at the chosen point. All other points belonging to this manifold correspond to reactive fixed points yielding $\lambda_{\text{r}}\neq 0$. Apart from possible exceptions (that we never observed in our numerical experiments), the manifold splits into two parts, one having $\lambda_{\text{r}}>0$, the other $\lambda_{\text{r}}<0$. Therefore, the system undergoes a Hopf bifurcation along the manifold for $\lambda_{\text{r}}=0$. Hence, the RE exhibit a continuum of Hopf bifurcations at all transition lines of Fig. 4. Unfortunately, we have been unable to clarify whether RE at the bifurcation points have first integrals, like $rps$ for homogeneous rates and, if affirmative, how they look analytically.

So far, we have assumed that $d_{\scriptscriptstyle\rm RRP}=d_{\scriptscriptstyle\rm PPS}=d_{\scriptscriptstyle\rm SSR}=1$. We now relax this constraint to examine another feature of Eqs. (II.6). We know that the cyclic Lotka-Volterra model with pairwise interactions has only neutrally stable orbits Reichenbach5 . Moreover, the ensemble of fixed points fills the ternary diagram, as can be seen from Eq. (4) of Ref. Reichenbach5 : for each pair $(r_{*},p_{*})$, there exists a choice of rates for which $(r_{*},p_{*})$ is a reactive fixed point. The reader may ask whether this feature holds similarly in our model. The problem is non-trivial because $(r_{*},p_{*})$ has a complex dependency on three-agent rates, as Eq. (IV.3) shows. To answer the question, we let ${\cal H}$ denote the ensemble of neutrally stable fixed points corresponding to a given choice of equilibrating rates, namely

[TABLE]

where integration over polarizing rates is understood. In Fig. 5, we reconstruct ${\cal H}$ numerically for several values of $d_{\scriptscriptstyle\rm RRP},\,d_{\scriptscriptstyle\rm PPS},\,d_{\scriptscriptstyle\rm SSR}$. Continuous lines correspond to sequences of polarizing rates. Points lying between neighboring lines belong to the ensembles as well. They just correspond to polarizing rates we did not consider numerically. Plot (i) shows that ${\cal H}(1,1,1)$ is a proper subset of the ternary diagram, symmetric under cyclic permutations. Changing one or two rates distorts its shape: the smaller the rates, the closer ${\cal H}$ shifts towards the boundaries. By extrapolation, Fig. 5 suggests that $\bigcup_{d_{\scriptscriptstyle\rm RRP},d_{\scriptscriptstyle\rm PPS},d_{\scriptscriptstyle\rm SSR}}{\cal H}(d_{\scriptscriptstyle\rm RRP},d_{\scriptscriptstyle\rm PPS},d_{\scriptscriptstyle\rm SSR})$ covers the whole diagram. However, one-to-one correspondence between fixed points and reaction rates is lost: several distinct sets of rates yield the same reactive fixed point. Therefore, we conclude, our model has a more complex algebraic structure than the cyclic Lotka-Volterra model with pairwise interactions.

Besides, Fig. 5 highlights that the the law of the weakest, first described in Ref. Berr1 , holds here as well. Take for instance plots (g)-(h)-(i): the lower $d_{\scriptscriptstyle\rm PPS}$, the closer ${\cal H}$ approaches the $(p,s)$ boundary. Recall that $d_{\scriptscriptstyle\rm PPS}$ mediates the equilibrating transition {{\color[rgb]{0,0.7,0.3}\definecolor[named]{pgfstrokecolor}{rgb}{0,0.7,0.3}\tt P}}\,{{\color[rgb]{0,0.7,0.3}\definecolor[named]{pgfstrokecolor}{rgb}{0,0.7,0.3}\tt P}}\,{{\color[rgb]{0.15,0.15,0.85}\definecolor[named]{pgfstrokecolor}{rgb}{0.15,0.15,0.85}\tt S}}\to{{\color[rgb]{0,0.7,0.3}\definecolor[named]{pgfstrokecolor}{rgb}{0,0.7,0.3}\tt P}}\,{{\color[rgb]{0.15,0.15,0.85}\definecolor[named]{pgfstrokecolor}{rgb}{0.15,0.15,0.85}\tt S}}\,{{\color[rgb]{0.15,0.15,0.85}\definecolor[named]{pgfstrokecolor}{rgb}{0.15,0.15,0.85}\tt S}}. As such, it yields a relative measure of the strength of {{\color[rgb]{0.15,0.15,0.85}\definecolor[named]{pgfstrokecolor}{rgb}{0.15,0.15,0.85}\tt S}} versus {{\color[rgb]{0,0.7,0.3}\definecolor[named]{pgfstrokecolor}{rgb}{0,0.7,0.3}\tt P}}. Any neutrally stable orbit surrounds the corresponding fixed point and flows counterclockwise on the ternary diagram. Therefore, the lower $d_{\scriptscriptstyle\rm PPS}$ the higher the probability that the system leaves its orbit by stochastic noise and falls eventually on the $(p,s)$ boundary, where $r=0$. When this happens, the dynamics terminates with $s=1$. Hence, S, the weakest strategy, survives, while R and P go extinct.

We finally investigate how the phase structure of the model changes off a well-mixed environment. To this aim, we make use of the Gillespie algorithm to simulate the dynamics on a two-dimensional lattice. We consider a square grid with $N=L\times L$ sites and periodic boundary conditions. We assume that each lattice site is occupied by a single agent carrying one of the strategies. Local interactions involve one agent and two of its nearest neighbors, placed in one of six possible configurations, as shown in Fig. 6.

Fig. 7 (a)-(b) shows phases corresponding to $d_{\scriptscriptstyle\rm RRS}=0.2$ and $d_{\scriptscriptstyle\rm RRS}=1.0$ for $N=256\times 256$. The most remarkable difference with respect to the well-mixed environment is the absence of bifurcations. In both plots the phase plane is entirely filled by stable fixed points with global densities (here meaning strategy fractions over the whole lattice) spiralling inwards. A comparison with plots (a) and (i) of Fig. 4 indicates that the relative position of phases is essentially the same in the region of stable equilibrium, even if their shapes are deformed. For instance, for $d_{\scriptscriptstyle\rm RRS}=1.0$, $d_{\scriptscriptstyle\rm PPR}<1$, $d_{\scriptscriptstyle\rm SSP}\gtrsim 1$, {{\color[rgb]{0,0.7,0.3}\definecolor[named]{pgfstrokecolor}{rgb}{0,0.7,0.3}\tt P}} has a relative majority on the lattice, whereas {{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}\tt R}} has a relative majority in a well-mixed environment. We conclude that changing spatial topology induces distortive effects analogous to those observed in Ref. Szolnoki2 .

As we discussed in section III, stochastic noise perturbates the evolution of strategies for $N<\infty$. On a two-dimensional lattice, stable spirals do not converge exactly. Global densities end up fluctuating erratically around the fixed point. Fig. 7 (c)-(d) shows this effect for two sets of rates, both corresponding to unstable equilibria in a well-mixed environment. In both plots, continuous lines are representative trajectories of the global densities, while dashed/dotted lines represent averages over 100 sample trajectories. As can be seen, stochastic noise averages to zero. The amplitude of similar erratic oscillations was studied with full detail in a model featuring species mutations Mobilia2 . In that context, a resonance amplification, occurring at a specific frequency, influenced fluctuations. The resonant frequency was estimated by the power spectrum method in the Fokker-Planck approximation. Here, we have less analytic information concerning the reactive fixed point. Hence, we do not attempt such an analysis. We only report that we observed a very light dependence of the fluctuation amplitude upon the polarizing rates in our numerical experiments. However, this result might depend on the excessive coarseness of the grid of values we chose.

V Agent mobility

In the previous sections, we have allowed strategies to propagate as a result of predation. Now, we let them diffuse explicitly via additional pair exchange reactions, namely

[TABLE]

all occurring with rate $\gamma_{2}$. We set up a lattice metapopulation model along the lines of Refs. McKane1 ; Lugo1 ; Szczesny1 ; Lamouroux1 ; Rulands1 . Specifically, we consider a two-dimensional square lattice with $N=L\times L$ sites and periodic boundary conditions. We interpret lattice sites as patches having a carrying capacity of $M\leq\infty$ agents. We consider all reactions listed in Eq. (II.1) as local processes, meaning that they always involve agents lying on the same patch. For simplicity, we assume homogeneous rates, as specified by Eqs. (II.7). By contrast, we interpret exchange reactions as bilocal processes, involving agents lying on two neighboring patches. Finally, we choose the lattice spacing to be $h=1/L$ so that the whole lattice has length one. As $N\to\infty$, the density field $\phi\equiv\left(r,p\right)(\mathbf{x},t)$ is governed by stochastic partial differential equations (SPDE), reading as

[TABLE]

where $D_{2}=\gamma_{2}/N$ is a scaling diffusion constant (in order to keep it finite we have to scale $\gamma_{2}\propto N$ as $N\to\infty$), $\Delta=\partial^{2}_{x}+\partial^{2}_{y}$ is the two-dimensional Laplace operator and $\xi_{i}$ denotes uncorrelated Gaussian noise. The matrix $C$ fulfills $CC^{\rm\scriptscriptstyle T}=B$, where $B$ is the diffusion matrix of the FPE. In section III, we obtained a simplified formula for $B$, holding for $d_{\rm e}=d_{\rm p}=1$. The most general expression to be used in Eqs. (V.4) reads

[TABLE]

We compute $C$ from $B$ via Cholesky decomposition whenever $B$ is a positive definite matrix. If $\phi$ falls on a vertex of the ternary diagram, we let $C_{ij}=0$ by continuity.

The dynamics of the system is trivial for $d_{\rm e}>d_{\rm p}$. In this case, $\phi$ converges uniformly to the reactive fixed point $(1/3,1/3)$ on all patches, up to small fluctuations, independently of the initial conditions. A convenient way to studying the dynamics for $d_{\rm p}\geq d_{\rm e}$ is to set up initial conditions such that each strategy occupies exclusively a finite portion of the lattice, as originally proposed and implemented in Refs. Szczesny2 ; Szczesny4 . We proceed identically, namely for $\mathbf{x}=(x,y)$ we let

[TABLE]

The advantage of such initial conditions is related to the special role of the four lattice points $\mathbf{x}=(0,0),(L/2,0),(0,L/2),(L/2,L/2)$. Here, all strategies meet and give rise to spiral waves for $t>0$. Just for the sake of completeness, we briefly review why this happens. To this aim, we introduce the topological current

[TABLE]

We assume that Greek indices take values $\{0,1,2\}$, while Latin indices take values $\{1,2\}$. Repeated indices are conventionally understood to be summed over their respective domains. The symbols $\epsilon_{\mu\nu\rho}$ and $\epsilon_{ab}$ denote totally antisymmetric tensors with three and two components respectively. We let $\partial_{0}\equiv\partial_{t}$. It takes no effort to show that $J_{\mu}$ fulfills the local conservation law $0=\partial_{\mu}J_{\mu}=\partial_{0}J_{0}+\partial_{k}J_{k}$. We define the topological charge in the thermodynamic limit as

[TABLE]

From the conservation of $J_{\mu}$ and the periodicity of the boundary conditions, it follows that $\text{d}Q/\text{d}t=\int\text{d}\mathbf{x}\,\partial_{0}J_{0}=-\int\text{d}\mathbf{x}\,\partial_{k}J_{k}=0$, hence $Q$ is invariant in time. We then let $J_{0}^{\pm}(\mathbf{x},t)=\max\{\pm J_{0}(\mathbf{x},t),0\}$. We notice that $J_{0}^{\pm}$ is strictly positive only near points around which the three strategies follow cyclically in counterclockwise/clockwise order. In particular, $J_{0}^{+}(\mathbf{x},0)$ is positive (infinite) for $\mathbf{x}=(0,L/2)$, $(L/2,0)$, while $J_{0}^{-}(\mathbf{x},0)$ is positive (infinite) for $\mathbf{x}=(0,0)$, $(L/2,L/2)$ and $J_{0}^{\pm}(\mathbf{x},0)=0$ elsewhere. As a result, we have four topological charges, two positive and two negative, localized on the four mentioned points, yielding $0=Q(0)=Q(t)$ for all $t>0$. While the topological density is sharply peaked for $t=0$, it becomes somewhat smooth for $t>0$. This corresponds to the appearance of two spirals plus two anti-spirals originating from the four points. Whether these objects last forever or disappear sooner or later is not a matter of topology; it depends on Eqs. (V.4).

In Fig. 8 we report a sequence of snapshots of $\phi$ in RGB representation. The sequence refers to the following choice of parameters:

• $\circ$

$N=2048\times 2048$;

• $\circ$

$\gamma_{2}=16$ ($D_{2}\simeq 3.81\times 10^{-6}$);

• $\circ$

$d_{\rm e}=1\,,\ d_{\rm p}=1.2^{k}$, for $k=0,1,\ldots,14$;

• $\circ$

$M=\infty$.

All pictures represent $\phi$ in the central area of the lattice, namely for $L/4<x_{1},x_{2}<3L/4$. Plot (a) corresponds to $t=0$, all others to $t=400$. The latter time is sufficiently large to ensure that transient effects have disappeared and the system has evolved to a steady state. To obtain Fig. 8, we integrated Eqs. (V.4) numerically via the ETD2RK scheme, introduced in Ref. Cox1 . We summarize below the main features of the plots:

• $\circ$

chaotic patterns with blurred and stretched shapes are predominant for $\epsilon<1$;

• $\circ$

blurring reduces as $\epsilon$ increases;

• $\circ$

small spirals emerge from chaos for $\epsilon\gtrsim 1$;

• $\circ$

a central spiral arises in a background of smaller spirals for $\epsilon\gtrsim 2$;

• $\circ$

the propagation radius of the central spiral increases progressively, until it becomes largest for $\epsilon\approx 10$;

• $\circ$

for even larger $\epsilon$ the central spiral keeps definitively stable;

• $\circ$

the wavelength of spatial patterns decreases monotonically as $\epsilon$ increases.

Patterns characterized by very similar behavior emerge in the spatial version of the May-Leonard model, featuring different particle interactions but analogous Hopf bifurcation Reichenbach2 . This similarity provides a confirmation that macroscopic phenomena induced by cyclic competition are robust. They depend only on the type of the Hopf bifurcation, whereas the details of the interactions have no qualitative (and little quantitative) influence Reichenbach1 ; Frey1 . We follow Ref. Reichenbach2 for the analysis of spiral waves in terms of the complex Ginzburg-Landau equation (CGLE). In our case, the CGLE reads (after rescaling the complex amplitude and shifting its phase)

[TABLE]

with $\lambda_{\text{r}}=\epsilon/6$ being the real part of the eigenvalues of the linearized RE, see Eq. (II.11), and with $\alpha(\epsilon)$ a function parameter given by

[TABLE]

It is important to recall that the CGLE is accurate only in the vicinity of a supercritical Hopf bifurcation. In that case, the role of the Laplacian is to synchronize limit cycles at neighboring lattice sites, thus giving rise to coherent spatiotemporal patterns. In our model, we have heteroclinic cycles, just like in Ref. Reichenbach2 . Hence, the CGLE is not guaranteed to describe correctly the dynamics of the system. Nonetheless, we can compare predictions of the CGLE with numerical observations. In particular, in the approach of Ref. VanSaarloos the wavelength $\ell(\epsilon)$ of spiral waves and their propagation speed $v(\epsilon)$ read

[TABLE]

Fig. 9 shows a comparison of Eq. (V.13) with numerical data corresponding to $N=512\times 512$, $M=\infty$ and $\gamma_{2}=1,2,4$. With regard to the wavelength, the agreement is perfect at all tested scales up to a multiplicative constant $\mu=1.55\div 1.60$, depending on $\gamma_{2}$ but not on $\epsilon$. This constant includes nonlinear effects which are not properly captured by the CGLE. Surprisingly, $\mu$ is very similar in size to the analogous constant in Ref. Reichenbach2 . As for the propagation speed, the agreement becomes very good only in the asymptotic regime $\epsilon\to\infty$. Moreover, the multiplicative constant $\mu$ is slightly smaller. This yields evidence that nonlinear effects may change from one observable to another555A different theoretical approach can be found in Ref. Szczesny2 . Here, $\ell$ and $v$ are expressed in the plane wave approximation as functions of $|z|^{2}$. The latter quantity is computed from numerical solutions of the CGLE via global averaging over the whole space..

An intuitive explanation of why the RGB representation of $\phi$ blurs as $\epsilon\to 0$ follows from the observation that $\dot{\rho}$ and $\dot{\theta}$ in Eq. (II.21) are monotonic functions of $\epsilon$ at LO. This feature has straightforward consequences for the dynamics in the metapopulation model. Indeed, the smaller $\epsilon$, the closer $\phi$ keeps wandering chaotically around the reactive fixed point before walking over heteroclinic cycles. As a result, RGB colors representing $\phi$ shift progressively to gray (the color corresponding to equal densities), while spatial patterns fade out.

Closer inspection suggests that two different physical mechanisms interfere with the spatiotemporal coherence of rotating spirals in regimes of intermediate and small $\epsilon$, as discussed in Refs. Reichenbach2 ; Szczesny2 ; Szczesny4 . For $2\lesssim\epsilon\lesssim 5$ small disturbances propagate radially together with the wavefronts. They intensify while traveling away from the core of spirals until they result in a far-field break-up of spatial coherence (convective instability). As $\epsilon\to 0$, the blurring of $\phi$ becomes very strong. In this limit, coherent propagation is fully compromised. Small disturbances grow locally. Their overall effect is to twist and stretch spatial patterns (absolute instability).

In Fig. 10, three snapshots of $\phi$, corresponding respectively to $M=256$ (a), $128$ (b) and $64$ (c), illustrate how stochastic noise induced by finite $M$ perturbates the stability of the central spiral wave. We obtained all plots from the numerical integration of Eqs. (V.4) for $N=2048\times 2048$, $\gamma_{2}=16$, $\epsilon=6.43$, $t=400$ and initial conditions given by Eq. (V.8). The increase rate of the spiral radius with $M$ indicates relatively slow convergence to the asymptotic limit, represented by Fig. 8 (m).

VI Three-agent chase reactions

Spatially structured games, admitting empty sites, usually assume that agents move via site hopping and pair exchange. These two processes correspond to distinct mobility reactions, namely

[TABLE]

with {\tt X},{\tt Y}\in\{{{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}\tt R}},{{\color[rgb]{0,0.7,0.3}\definecolor[named]{pgfstrokecolor}{rgb}{0,0.7,0.3}\tt P}},{{\color[rgb]{0.15,0.15,0.85}\definecolor[named]{pgfstrokecolor}{rgb}{0.15,0.15,0.85}\tt S}}\}. As $N\to\infty$, they induce differential variations of the strategy densities, quantified by

[TABLE]

for $(a,b,c)$ a permutation of $(r,p,s)$. Eq. (VI.3) reduces to Gaussian diffusion for $\gamma_{\text{h}}=\gamma_{\text{e}}$. Nonlinear terms, arising for $\gamma_{\text{h}}\neq\gamma_{\text{e}}$, can disrupt the stability of spiral waves by producing perturbations that result in the far-field break-up of spatial coherence Szczesny1 . Hence, they could be important to clarify in which circumstances growing bacterial colonies develop coherent patterns and in which they do not Szczesny1 . Anyway, hopping and pair exchange yield homogeneous and isotropic diffusion, regardless of how $\gamma_{\text{h}}$ and $\gamma_{\text{e}}$ are chosen.

In this section, we address the issue of whether we can introduce nonlinear diffusion in our model. Since we admit no empty sites, Eq. (VI.1) is ruled out. Hence, Eq. (VI.3) makes sense no more. The mobility operator induced by pair exchange is the Laplacian, $\delta a_{\text{ex}}=(\gamma_{2}/N)\Delta a$ for $a=r,p,s$. We have three possibilities: either we give up species homogeneity in pair exchange, or we break isotropy, or, less trivially, we introduce more complex mobility reactions. Going back to Table 1 and Eq. (II.1), we realize that we can tentatively use for our purpose transitions where agents carry initially three different strategies. Given a lattice site $\mathbf{x}$ and two nearest neighbors $\mathbf{y}_{1},\,\mathbf{y}_{2}$, we introduce the chase reactions

[TABLE]

for $({\tt X},{\tt Y},{\tt Z})$ an even permutation of ({{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}\tt R}},{{\color[rgb]{0,0.7,0.3}\definecolor[named]{pgfstrokecolor}{rgb}{0,0.7,0.3}\tt P}},{{\color[rgb]{0.15,0.15,0.85}\definecolor[named]{pgfstrokecolor}{rgb}{0.15,0.15,0.85}\tt S}}) and

[TABLE]

for $({\tt X},{\tt Y},{\tt Z})$ an odd permutation of ({{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}\tt R}},{{\color[rgb]{0,0.7,0.3}\definecolor[named]{pgfstrokecolor}{rgb}{0,0.7,0.3}\tt P}},{{\color[rgb]{0.15,0.15,0.85}\definecolor[named]{pgfstrokecolor}{rgb}{0.15,0.15,0.85}\tt S}}). Incidentally, the interpretation of Eqs. (VI.4)-(VI.5) as reactions where species chase one another cyclically is not the only possible one. We could, equivalently and perhaps more imaginatively, regard them as an evolutionary version of the playground singing game Ring a ring of roses. Anyway, we assume that all such reactions occur with rate $\gamma_{3}$. Fig. 11 shows the initial state corresponding to all possible permutations.

We wish to derive the mobility operator arising from Eqs. (VI.4)-(VI.5), as $N\to\infty$. To this aim, it is sufficient to focus on one of the strategies, e.g. R. The average variation of $r(\mathbf{x})$, due to the above reactions, receives two negative and two positive contributions, corresponding respectively to initial configurations 1.-2. and 5.-6., namely

[TABLE]

with

[TABLE]

In principle, there are six possible ways of choosing $\mathbf{y}_{1},\,\mathbf{y}_{2}$ given $\mathbf{x}$. An elegant expression for the mobility operator follows provided we discard configurations (a)-(b) of Fig. 6 and keep all others. The variation of $r(\mathbf{x})$ corresponding to configuration (c) of Fig. 6 reads as

[TABLE]

Expanding it in Taylor series around $h=0$ yields

[TABLE]

with

[TABLE]

We find analogous expressions $\delta r_{\text{k}}(\mathbf{x})$ for k = d, e, f, corresponding respectively to initial configurations (d), (e), (f) of Fig. 6. For the sake of conciseness, we leave the reader with the exercise of deriving them. Upon adding $\delta r_{\text{c}}$ and $\delta r_{\text{d}}$, the mixed-derivative terms at $\text{O}(h^{2})$ cancel, hence we obtain

[TABLE]

Analogously, adding $\delta r_{\text{e}}$ and $\delta r_{\text{f}}$ yields

[TABLE]

As can be seen, the $\text{O}(h)$ contributions, Eqs. (VI.17), (VI.20), are equal and opposite, while the $\text{O}(h^{2})$ ones, Eqs. (VI.18), (VI.21), are just equal. Consequently, the former cancel whereas the latter add up. At the end, we get

[TABLE]

with scaling diffusion constant $D_{3}=2\gamma_{3}/N$.

Eq. (VI.22) yields a continuous chase operator. It is just one of several possible definitions. We could produce others, for instance, by differently choosing configurations from Fig. 6. Unsurprisingly, $\delta r_{\text{chase}}$ is nonlinear (cubic) in the strategy densities. At first sight, the expression in square brackets looks smaller than $\Delta r(\mathbf{x})$ because two densities multiply the Laplace operators. In view of this, we should reasonably expect that Eq. (VI.22) contributes as a small perturbation to Eq. (V.4) for $\gamma_{3}\approx\gamma_{2}$. To make $\delta r_{\text{chase}}$ comparable in strength to $\delta r_{\text{ex}}$, we should let $\gamma_{3}$ larger that $\gamma_{2}$ by a factor of about $3\div 5$. Generally, reaction-diffusion equations with nonlinear diffusion operators cannot be mapped onto a CGLE. Hence, the formation of stable spirals is ruled out.

Fig. 12 shows six snapshots of $\phi$ corresponding to $N=2048\times 2048$, $\gamma_{2}=0$, $\gamma_{3}=64,128,256$, $\epsilon=9.70$ and $M=\infty$. We obtained them by numerically integrating the SPDE with random initial conditions. To this aim, we used a semi-implicit Runge-Kutta scheme of order two, with time step $\text{d}t=0.005$ and mobility reactions represented in configuration space (exponential time differencing is not possible in this case).

Plots (a)-(b)-(c) show $\phi$ at an early stage. They correspond respectively to $t=5,17,36$. As can be seen, strategy densities try to arrange initially into spatial patterns with size and shapes depending on $\gamma_{3}$. In particular, as $\gamma_{3}$ increases, patterns grow. Soon, chaos replaces all regular patterns. Plots (d)-(e)-(f) represent $\phi$ for $t=400$. No sign of spatial coherence is left. Patterns are chaotic for all values of $\gamma_{3}$. They keep evolving rapidly in time. Their size increases with $\gamma_{3}$ as expected, while their boundaries, as we observed in unreported snapshots, can be equally sharp or blurred.

Fig. 13 shows four snapshots of $\phi$ corresponding to $N=3072\times 3072$, $\gamma_{2}=36$, $\gamma_{3}=18,36,54,72$, $\epsilon=9.70$ and $M=\infty$. We obtained them by numerically integrating the SPDE with initial conditions given by Eq. (V.8). To this aim, we used the same integration scheme as explained above. We chose $\epsilon=9.70$ since we know from Fig. 8 (o) that the evolution of the strategy densities for $\gamma_{3}=0$ yields a fully extended and perfectly stable central spiral. Hence, this choice allows us to assess the effects of the chase operator in a controlled set up.

Similar to stochastic noise, nonlinear diffusion is responsible for the far-field break-up of the central spiral. The break-up mechanism seems to be essentially the same: nonlinear disturbances propagate from the core outwards; they grow in size while propagating; spatial coherence breaks down as soon as disturbances overcome the carrier signal. However, a comparison with Fig. 10 highlights two differences regarding the breaking pattern:

• $\circ$

stochastic noise deforms the profile of wavefronts near the break-up radius. This effect is absent in Fig. 13, where all wavefronts up to the break-up radius are nearly perfect. If our interpretation is correct, the effects of nonlinear diffusion reveal suddenly, whereas stochastic noise is somewhat progressive;

• $\circ$

in Fig. 10, spatial coherence is essentially lost beyond the break-up radius. Small spirals fill the environment without large emerging patterns. This behavior is at odds with Fig. 13, where spatial coherence reappears beyond the break-up radius for $\gamma_{3}\lesssim\gamma_{2}$. External patterns consist of almost-linear wavefronts. Nearly oriented towards the main lattice directions, they form cusps along the main diagonals. Moreover, they resemble the boxy-shaped patterns of Fig. 12 (c). The effect of coherence recovery reduces progressively as $\gamma_{3}$ increases. Paradoxically, nonlinear diffusion has the power to rectify spiral waves.

In this paper, we do not go beyond the above qualitative observations. Our goal was to show that mobility reactions, different from hopping and pair exchange, can have non-trivial and interesting consequences on the propagation of spiral waves. Reaction-diffusion equations with nonlinear diffusion mechanisms have attracted much attention in the past years, thanks to their potential applications to many natural phenomena. We refer the reader to the Refs. Gilding1 ; Gilding2 ; Gilding3 ; Gilding4 ; Vazquez1 for in-depth studies.

VII Conclusions

To summarize, we have investigated several aspects of a variant of the cyclic Lotka-Volterra model, featuring three-agent interactions. We aimed at studying cyclic dominance, mediated by cooperative predation, in a simple theoretical setting.

Nonlinear analysis of the underlying rate equations in a well-mixed environment has revealed the existence of degenerate Hopf bifurcations. They occur for specific values of the rate constants. More precisely, in our model, reactions involving two prey and one predator equilibrate the system, while reactions involving two predators and one prey polarize it. Bifurcations correspond to non-trivial equilibria between the rates of the former and the latter reactions. If equilibrating and polarizing reactions have homogeneous rates, respectively $d_{\rm e}$ and $d_{\rm p}$, the rate equations bifurcate for $d_{\rm e}=d_{\rm p}$. This condition describes predators hunting in a group or alone with equal propensity. In a metapopulation model with patches hosting several agents, rotating spiral waves appear only for $d_{\rm e}<d_{\rm p}$, i.e., when the propensity for cooperative predation is stronger than for individual hunting.

Theoretical methodologies used in this paper have been developed elsewhere in the literature. We have just applied them to our model to compare features of the underlying dynamics with other existing models.

In particular, we have derived the magnitude of the stochastic noise at the bifurcation point for homogeneous rates (where the rate equations predict neutrally stable orbits), to make a comparison with the original cyclic Lotka-Volterra model Reichenbach5 . Three-agent interactions are intrinsically noisier than two-agent ones since they involve more fluctuating degrees of freedom. Nevertheless, the extinction probability is uniformly lower in our model than in Ref. Reichenbach5 . This apparent paradox has a simple solution. When three agents interact, strategies fluctuate longer around the reactive fixed point before one of them prevails on the others. Stochastic noise has no preferred direction. Hence, it acts as an equilibrating force. Group interactions help the system stay in equilibrium. Doing so, they promote species coexistence.

Similarly, we have studied the phase portrait for heterogeneous rates, to make a comparison with a model in which group interactions involve an agent and its four von-Neumann neighbors or its eight Moore neighbors Szolnoki2 . Although we consider only three interacting agents, the phase portrait of our model shows a richer structure. As far as we understand, the reason for this result is that our reaction rates are fully independent. Besides, we have shown that spatial topology plays a critical role in shaping the phase portrait. Indeed, Hopf bifurcations disappear on a two-dimensional lattice (with one agent per lattice site) as a consequence of the locality of the interactions.

Then, we have studied the effects of individual mobility in a lattice metapopulation model, with patches hosting several agents. It turns out that our rotating spirals are qualitatively similar to those arising in the spatially structured version of the May-Leonard model Reichenbach2 . Although our model differs from Ref. Reichenbach2 in that we assume no empty sites, no birth and no selection-removal interactions, the observed similarity is not surprising, for three reasons. First of all, models of cyclic dominance without mutations usually undergo degenerate Hopf bifurcations. Secondly, it is usually possible to map a system of reaction-diffusion equations, with supercritical Hopf bifurcation, onto a complex Landau-Ginzburg equation Reichenbach1 ; Frey1 . The map is possible provided diffusion is linear. Finally, the latter equation yields a good description, even when the bifurcation is degenerate. As a result, the regimes of species coexistence in our model and the spatially structured May-Leonard model fall in the same universality class. As such, they are in one-to-one correspondence.

To conclude, we have shown that one can build nonlinear continuous mobility operators starting from three-agent chase reactions on a lattice. We have focused on one of several possible definitions. In particular, we have studied the effects of chase reactions on the propagation of spiral waves. Similar to nonlinear diffusion, arising in models where hopping and pair exchange occur with different rates, chase reactions produce far-field break-up of spiral waves. However, the breaking pattern is peculiar. Nonlinear diffusion could play a role in explaining structural patterns in growing bacterial colonies. Our study shows that group interactions provide viable mechanisms for both predation and dispersal.

Acknowledgments

F.P. acknowledges his little daughter Giulia for (re)inventing RPS with three players and for suggesting it as a research topic. The computing resources used for our numerical simulations and the related technical support have been provided by the CRESCO/ENEAGRID High Performance Computing infrastructure and its staff Ponti . CRESCO (Computational RESearch centre on COmplex systems) is funded by ENEA and by Italian and European research programmes. Author contribution statement F.P. performed theoretical calculations, numerical simulations, and data analysis with assistance from S.F. and S.T. All authors contributed to the discussion and the writing of the manuscript.