Predominance of the weakest species in Lotka-Volterra and May-Leonard implementations of the rock-paper-scissors model

TL;DR

This study compares Lotka-Volterra and May-Leonard models of rock-paper-scissors dynamics, showing that the 'weakest' species often prevails in large systems, with size-dependent effects observed in smaller simulations.

Contribution

It demonstrates that both models yield similar long-term species abundance patterns despite different dynamics, emphasizing the role of system size in species survival.

Findings

Weakest species tends to dominate in large systems.

Initial oscillations influence survival probabilities in small systems.

Both models produce qualitatively similar results for large lattices.

Abstract

We revisit the problem of the predominance of the 'weakest' species in the context of Lotka-Volterra and May-Leonard implementations of a spatial stochastic rock-paper-scissors model in which one of the species has its predation probability reduced by $0 < \mathcal{P}_w < 1$. We show that,despite the different population dynamics and spatial patterns, these two implementations lead to qualitatively similar results for the late time values of the relative abundances of the three species (as a function of $\mathcal{P}_w$), as long as the simulation lattices are sufficiently large for coexistence to prevail --- the 'weakest' species generally having an advantage over the others (specially over its predator). However, for smaller simulation lattices, we find that the relatively large oscillations at the initial stages of simulations with random initial conditions may result in a significant dependence of the probability of species survival on the lattice size and total simulation time.

Equations

$$i(i+1)\to ii,$$

$${\dot{\rho }}_i=p_i{\rho }_i{\rho }_{i+1}-p_{i-1}{\rho }_{i-1}{\rho }_i,$$

$${\rho }_{i+1}=\frac{p_{i-1}}{p_i}{\rho }_{i-1},\sum _{i=1}^3{\rho }_i=1.$$

$${\rho }_1={\rho }_2=\frac{1}{2+{\mathcal{P}}_w}{\rho }_3=\frac{{\mathcal{P}}_w}{2+{\mathcal{P}}_w},$$

$$i(i+1)\to i0,$$

$$i0\to ii,$$

$${\dot{\rho }}_i=r{\rho }_i{\rho }_0-p_{i-1}{\rho }_{i-1}{\rho }_i,$$

$${\dot{\rho }}_0=-r{\rho }_0\sum _{i=1}^3{\rho }_i+\sum _{i=1}^3{p}_{i-1}{\rho }_{i-1}{\rho }_i.$$

$$p_{i-1}{\rho }_{i-1}=r{\rho }_0,{\rho }_0+\sum _{i=1}^3{\rho }_i=1.$$

$${\rho }_0$$

$${\rho }_1$$

$${\rho }_2={\rho }_3$$

$$i\odot \to \odot i,$$

Full text

I Introduction

Non-hierarchical competition interactions have been shown to play a crucial role in the preservation of coexistence. The classical rock-paper-scissors (RPS) model Kerr et al. (2002); Reichenbach et al. (2007), the simplest cyclic predator-prey model, describes the dynamics of three different species subject to interspecific competition (see Lotka (1920); Volterra (1926); May and Leonard (1975) for the pioneer work by Lotka and Volterra, and May and Leonard). It allows for the stable coexistence of all three species, and successfully reproduces some of the main dynamical features observed in simple biological systems with cyclic selection interactions Sinervo and Lively (1996); Kerr et al. (2002); Kirkup and Riley (2004).

Simulations of the spatial RPS model are usually performed on a square lattice (see Szabó et al. (2004); Zhang et al. (2009); Laird (2014); Rulquin and Arenzon (2014), however, for other lattice configurations) and consider nearest-neighbor cyclic predator-prey interactions. In a Lotka-Volterra implementation of the RPS model, each site is occupied by a single individual of one of the three species, and there is a conservation law for the total number of individuals. On the other hand, in a May-Leonard implementation each site may either be occupied by a single individual or left empty, and the total number of individuals is, in general, not conserved Peltomäki and Alava (2008).

Generalizations of the RPS model involving additional species and interactions, have also been investigated in recent years Szabó et al. (2008); Allesina and Levine (2011); Avelino et al. (2012a, b); Li et al. (2012); Roman et al. (2012); Lütz et al. (2013); Roman et al. (2013); Cheng et al. (2014); Szolnoki et al. (2014); Kang et al. (2016); Roman et al. (2016); Brown and Pleimling (2017); Park et al. (2017); Bazeia et al. (2017); Souza-Filho et al. (2017); Esmaeili et al. (2018); Avelino et al. (2019). Complex dynamical spatial structures (such as spirals with an arbitrary number of arms Avelino et al. (2012b); Bazeia et al. (2017, 2019), domain interfaces, with or without non-trivial internal dynamics Avelino et al. (2014a), and string networks, with or without junctions Avelino et al. (2014b); Avelino et al. (2017)), diverse scaling laws Avelino et al. (2012b); Brown and Pleimling (2017), and phase transitions Szabó and Czárán (2001); Szabó and Arial Sznaider (2004); Szolnoki and Szabó (2004); Perc et al. (2007); Szabó et al. (2007); Szabó and Szolnoki (2008); Szolnoki et al. (2011); Vukov et al. (2013); Bazeia et al. (2018) have been shown to naturally emerge in some of these scenarios.

In most of these models the species may be characterized as having equal strength, with the survival probability being mainly dependent on initial conditions. However, there are other situations in which there is a competitive difference between species, such as in the case of a RPS model in which one of the species (often termed the ‘weakest’) has a reduced predation probability ${\mathcal{P}}_{w}$. It has been shown in Refs. Frean and Abraham (2001); Berr et al. (2009) that in a Lotka-Volterra implementation of this model, the ‘weakest’ species tends to be the most abundant. These results have recently been challenged in Ref. Menezes et al. (2019), with the authors claiming that some of the model parameters have a significant impact on which species survives in a May-Leonard implementation.

In this paper we shall address the question of whether the predominance of the ‘weakest’ species generally occurs in both Lotka-Volterra and May-Leonard implementations of the RPS model. The outline of this paper is as follows. In Sec. II we start by considering a non-spatial RPS model in which one of the species has a reduced predation probability, discussing the properties of its stationary solutions in both Lotka-Volterra and May-Leonard implementations. In Sec. III we describe these two implementations of the spatial stochastic RPS model and present the corresponding results. Special emphasis is given to the dependence of the survival probability on the size of the simulation lattice and total simulation time, and to the way the average densities of the three species depend on the reduced predation probability, parameterized by ${\mathcal{P}}_{w}$, for sufficiently large simulation lattices. Finally, we conclude in Sec. IV.

II Non-spatial RPS model

Let us start by considering Lotka-Volterra and May-Leonard implementations of the non-spatial RPS model.

II.1 Lotka-Volterra

A Lotka-Volterra implementation of the RPS model considers three species with densities $\rho_{1}$, $\rho_{2}$ and $\rho_{3}$, such that $\rho_{1}+\rho_{2}+\rho_{3}=1$ (the total density is normalized to unity). At each timestep an individual of one of the species $i$ is selected at random and the predation interaction

[TABLE]

with $i=1,...,3$, is performed with probability $p_{i}$. In this paper, modular arithmetic, where integers wrap around upon reaching $1$ or $3$, is assumed (the integers $i$ and $j$ represent the same species whenever $i=j\,{\rm mod}\,3$, where $\rm mod$ denotes the modulo operation).

With an appropriate choice of time unit, the equations for the evolution of the densities of the different species may be written as

[TABLE]

where a dot represents a derivative with respect to time. Stationary solutions to Eq. (1) satisfy the condition ${\dot{\rho}_{i}}=0$, and are, therefore, characterized by

[TABLE]

Here, we shall be interested in the case where $p_{1}={\mathcal{P}}_{w}\,p$ and $p_{2}=p_{3}=p$, with $0<{\mathcal{P}}_{w}<1$, so that $\rho_{1}=\rho_{2}=\rho_{3}/{\mathcal{P}}_{w}$. Hence, Eq. (2) implies that the stationary solutions of Eq. (1) are characterized by

[TABLE]

with $\rho_{1}=\rho_{2}>\rho_{3}$.

II.2 May-Leonard

In a May-Leonard implementation of the RPS model the total density of individuals is no longer conserved. In this case, $\rho_{0}+\rho_{1}+\rho_{2}+\rho_{3}=1$ where, for uniformity of notation, $\rho_{0}$ shall be referred to as the density of empty sites — denoted by a ‘0’ — even when considering a non-spatial RPS model. At each timestep an individual of one of the species $i$ is selected at random and an interaction is performed: either predation

[TABLE]

with probability $p_{i}$, or reproduction

[TABLE]

with probability $r$ (assumed to be the same for all the species) — notice that predation has a different meaning in Lotka-Volterra and May-Leonard implementations of the RPS model. Again, with an appropriate choice of time unit, the equations for the evolution of the densities of the different species may be written as

[TABLE]

while the evolution of the density of empty sites is given by

[TABLE]

Stationary solutions to Eqs. (4) and (5) satisfy the conditions ${\dot{\rho}_{i}}=0$ and ${\dot{\rho}_{0}}=0$, and are, therefore, characterized by

[TABLE]

Again, we shall be interested in the case where $p_{1}={\mathcal{P}}_{w}p$, with $p_{2}=p_{3}=p$, with $0<{\mathcal{P}}_{w}<1$. Equation (6) implies that the stationary solutions to Eqs. (4) and (5) are characterized by

[TABLE]

with $\rho_{1}=\rho_{2}/{\mathcal{P}}_{w}=\rho_{3}/{\mathcal{P}}_{w}$, so that $\rho_{1}>\rho_{2}=\rho_{3}$.

Hence, we may conclude that in both Lotka-Volterra and May-Leonard implementations of the non-spatial RPS model the ‘weakest’ species ($1$) has a competitive advantage. In a May-Leonard implementation the stationary density of individuals of the ‘weakest’ species is larger than that of the other two species. On the other hand, in a Lotka-Volterra implementation the competitive advantage is less pronounced, since the stationary density of the ‘weakest’ species is only larger than that of its predator (its prey having an equal density).

III Spatial rock-paper-scissors model

In this section we shall describe Lotka-Volterrra and May-Leonard implementations of the spatial RPS model which we shall investigate in the present paper. To this end, we shall consider a $N^{2}$ square lattice — $N$ shall be referred to as its linear size — with ${\mathcal{N}}$ sites and periodic boundary conditions. In a Lotka-Volterra implementation every site is occupied by a single individual of one of the three-species, while in a May-Leonard implementation there is also the possibility of a site being empty. The number of individuals of the species $i$ and the number of empty sites will be denoted by $I_{i}$ and $I_{0}$, respectively — the density of individuals of the species $i$ and the density of empty sites shall be defined by $\rho_{i}=I_{i}/{\mathcal{N}}$ and $\rho_{0}=I_{0}/{\mathcal{N}}$, respectively (note that $\rho_{0}=0$ in a Lotka-Volterra implementation). The possible interactions are the ones described in the case of the non-spatial RPS model, plus mobility

[TABLE]

where $\odot$ represents either an individual of any species or an empty site.

At every simulation step, the algorithm randomly picks an occupied site to be the active one, randomly selects one of its adjacent neighbour sites to be the passive one, and randomly chooses an interaction to be executed by the individual at the active position: predation, mobility or reproduction with probabilities $p$, $m$ and $r$, respectively ($r=0$ in a Lotka-Volterra implementation) — except if stated otherwise, in this paper we use the von Neumann neighbourhood (or 4-neighbourhood) composed of a central cell (the active one) and its four non-diagonal adjacent cells. These three actions are repeated until a possible interaction is selected — note that in both implementations of the RPS model the interaction cannot be carried out whenever predation is selected and the passive is not a prey of the active, while in a May-Leonard implementation the interaction is not completed also if reproduction is selected and the passive is not an empty site.

A generation time (our time unit) is defined as the time necessary for $\mathcal{N}$ successive interactions to be completed.

III.1 Results

Figures 1 and 2 show, respectively, the population network evolution in Lotka-Volterra and May-Leonard implementations of the spatial stochastic RPS model with random initial conditions — each site being initially occupied by a randomly chosen single individual of any of the three species with a uniform discrete probability of $1/3$.

In Fig. 1 (Lotka-Volterra implementation) the upper and lower panels display snapshots of the spatial distribution of the different species on a $1000^{2}$ lattice at $t_{0}=0$, $t_{1}=16$, $t_{2}=23$, $t_{3}=30$, $t_{4}=40$, $t_{5}=48$, $t_{6}=98$, and $t_{7}=500$, for a single Lotka-Volterra realization of the spatial stochastic RPS model with $m=0.5$, $p=0.5$, and ${\mathcal{P}}_{w}=0.5$ — species $1$, $2$, and $3$ are represented in red, blue and green, respectively. Notice the change in the overall color tone which takes place in the early stages of the simulation, associated to changes in the densities of the three species. Such oscillations are clearly visible in the central panel of Fig. 1 which shows the evolution of the density $\rho_{i}$ of the different species for the entire timespan of the simulation — the red, blue and green lines (from top to bottom, respectively) representing the densities of species 1, 2 and 3, respectively. Fig. 1 shows that after short transient initial stage, with relatively large coherent oscillations, the densities of the three species quickly approach nearly constant values, with $\rho_{1}>\rho_{2}>\rho_{3}$. It reveals the predominance of the ‘weakest’ species ($1$), specially over its predator ($3$).

In Fig. 2 (May-Leonard implementation) the upper and lower panels show snapshots of the spatial distribution of the different species on a $1000^{2}$ lattice at $t_{0}=0$, $t_{1}=17$, $t_{2}=51$, $t_{3}=118$, $t_{4}=264$, $t_{5}=454$, $t_{6}=740$, and $t_{7}=5000$, for a single May-Leonard realization of the spatial stochastic RPS model with $m=0.5$, $p=0.25$, $r=0.25$, and ${\mathcal{P}}_{w}=0.5$ — species $1$, $2$, and $3$, and empty sites are represented in red, blue, green and white, respectively. The most prominent feature in the snapshots shown in Fig. 2 is the presence of distinctive compact spatial domains of increasing characteristic size in a variable one-species background up to $t\sim 10^{3}$. At larger $t$ the percolation between three-species spatial domains eventually leads to a population network of spiral patterns. The central panel of Fig. 2 depicts the evolution of the density of the different species and empty sites ($\rho_{i}$ and $\rho_{0}$, respectively). As in Fig. 1, the red, blue and green lines represent the densities of species 1, 2 and 3, respectively, but in this case there are also empty sites whose density is given by the grey bottom line. Fig. 2 shows that in a May-Leonard implementation there is also a transient initial stage prior to an asymptotic regime in which the densities of the three species quickly approach nearly constant values, with \rho_{1}\mathrel{\hbox to0.0pt{\lower 4.0pt\hbox{\hskip 1.0pt\sim}\hss}\raise 1.0pt\hbox{>}}\rho_{2}>\rho_{3}. However, the evolution is considerably slower and the fluctuations are considerably larger compared to a Lotka-Volterra implementation.

In the case of a May-Leonard implementation, the large coherent oscillations of the abundances of the various species in the early stages of simulations of the spatial RPS model with random initial conditions may result in a significant dependence of the surviving/most abundant species on the linear size of the lattice bellow a given linear size threshold $N_{th}$, even in the case of simulations with a large total simulation time. This is shown in Fig. 3 which depicts the probability of single species survival and coexistence as a function of the linear size of the simulation lattice for a May-Leonard implementation of the spatial stochastic RPS model with $m=0.3$, $p=0.35$, $r=0.35$ and ${\mathcal{P}}_{w}=0.5$. Each point was estimated from $10^{4}$ simulations with a total simulation time equal to $2\times 10^{4}$ generations, starting from random initial conditions with $\rho_{1}=\rho_{2}=\rho_{3}=1/3$. The error bars are always much smaller than the size of the symbols: the one-sigma uncertainty in the value of $P$, at each point, may be estimated as $(P(1-P)/10^{4})^{1/2}$, with a maximum of $5\times 10^{-3}$ for $P=0.5$. Fig. 3 shows that for linear sizes $N>N_{th}\sim 30$, the ‘weakest’ species has the largest probability to survive, but this no longer holds for $N<N_{th}$.

Figure 4 is analogous to Fig. 3 but considers a different choice of model parameters: $m=0.5$, $p=0.25$, $r=0.25$, and ${\mathcal{P}}_{w}=0.5$. The larger mobility leads to an increase of the lattice linear size above which the ‘weakest’ is the most likely to survive (in this case, $N_{th}\sim 110$), thus showing that this threshold is strongly dependent on the choice of models parameters. We also verified that the use of a Moore neighbourhood — composed of a central cell (the active one) and the eight cells that surround it — leads to similar qualitative results to the ones presented in Figs. 3 and 4 for a von Neumann neighbourhood, albeit with significantly larger linear thresholds ($N_{th}\sim 70$ and $N_{th}\sim 370$, respectively). Hence, the small linear size ($N=50$) associated to the limited total simulation time ($t=250$) of the simulations performed in Ref. Menezes et al. (2019) using a Moore neighbourhood explains the reported impact of some of the model parameters on the determination of the surviving species in a May-Leonard implementation of the RPS model.

Figure 5 shows the value of the average density of the three species as a function of $\mathcal{P}_{w}$ for a Lotka-Volterra implementation of the spatial stochastic RPS model. The data points result from an average over the last $10^{4}$ generations of simulations with a timespan equal to $1.5\times 10^{3}$ generations performed on a $1000^{2}$ lattice. The results for ${\mathcal{P}}_{w}=1$ were computed first, starting from random initial conditions (as in Fig. 1). The final conditions of each simulation with ${\mathcal{P}}_{w}=1$ were used as initial conditions for a new simulation with ${\mathcal{P}}_{w}=1-0.01$. This procedure was repeated until ${\mathcal{P}}_{w}=0.01$ was reached. Such an approach was used in order to allow for a fast convergence (we verified that, with such conditions, $5\times 10^{3}$ simulations are sufficient for $\langle\rho_{i}\rangle$ to attain its asymptotic value). In this way the large oscillations at the initial stages of simulations with random initial conditions shown in Figs. 1 and 2 — which, depending on the value of ${\mathcal{P}}_{w}$ could be responsible for the loss of coexistence on a relatively short timescale — are avoided. Hence, this choice of initial conditions allowed us to obtain results which, in the case random initial conditions, would require larger simulation lattices. Figure 5 shows that the ‘weakest’ species is always the most abundant, thus having a competitive advantage over the others, specially over its predator. Figure 5 also shows that competitive advantage over the other species increases as ${\mathcal{P}}_{w}$ decreases — the ‘weakest’ species and its prey having similar abundances for ${\mathcal{P}}_{w}>0.6$. The lines in Fig. 5 represent the stationary solution, given in Eq. (3), for the density of the species $1$ and $2$ (solid line) and $3$ (dashed line) obtained in the context of a Lotka-Volterra implementation of the non-spatial RPS model. Notice the remarkable agreement between the spatial and non-spatial results in a Lotka-Volterra implementation of the RPS model.

Figure 6 is analogous to Fig. 5, except that, in this case, a May-Leonard realization of the RPS model with $m=0.5$, $p=r=0.25$ is considered. Notice that, despite the considerably different population dynamics and spatial patterns, the late time asymptotic values of the relative abundances of the three species (as a function of ${\mathcal{P}}_{w}$) obtained for a May-Leonard implementation are qualitatively similar to the ones shown in Figure 5 for a Lotka-Volterra implementation. In both cases the ‘weakest’ species generally has a competitive advantage over the others — specially over its predator. Again, this is particularly true at low values of ${\mathcal{P}}_{w}$ in both implementations. However, in a May-Leonard implementation there is a regime, for $0.6<{\mathcal{P}}_{w}<1$, in which the prey of the ‘weakest’ species (species 2) is the dominant one, albeit only by a small margin. The lines in Fig. 6 represent the stationary solution, given in Eqs. (7-9), for the density of the species $1$ (solid line) and of species $2$ and $3$, and empty sites (dashed line) obtained in the context of a May-Leonard implementation of the non-spatial RPS model with $r=p$. In the case of a May-Leonard implementation of the RPS model the differences between the spatial and non-spatial results are significant. This is a result of the distinct spatial structure and of the associated dynamics generated in a May-Leonard implementation of the spatial RPS model.

We verified that the results shown in Figs. 5 and 6 would remain essentially unchanged if a Moore neighbourhood had been used instead of a von Neumann one. We also checked that for other values of $p$, $m$, and $r$ (with $r=0$ in the case of a Lotka-Volterra implementation) the results obtained for the dependence of the values of the average densities on $\mathcal{P}_{w}$ are qualitatively similar to the ones shown in Figs. 5 and 6. This is partially explained by the fact that different choices of the parameters may, to some extent, be absorbed by a redefinition of time and spatial units. In particular, the stationary solution for the values of the average densities in the non-spatial Lotka-Volterra implementation given in Eq. (3) only depends on ${\mathcal{P}}_{w}$ — any dependence on $p$ may be absorbed in the choice of a different time unit, which does not affect the stationary solutions. In the case of a May-Leonard implementation the stationary average densities given in Eqs. (7-9) depend both on ${\mathcal{P}}_{w}$ and $r/p$, but the dependence on $r/p$ has no impact on which species is the most abundant one. In a spatial version of the PRS model mobility also plays an important role. However, in a mean field description of a May-Leonard implementation of the RPS model changes of $m$ may be absorbed by an appropriate redefinition of spatial units Reichenbach et al. (2007).

IV Conclusions

In this paper we revisited the problem of the predominance of the ‘weakest’ species in the context of Lotka-Volterra and May-Leonard implementations of a spatial stochastic RPS model in which one of the species has a reduced predation probability. We have shown that, despite the significant dynamical differences between Lotka-Volterra and May-Leonard implementations of the RPS model, for sufficiently large lattices the late time values of the relative abundances of the various species display similar qualitative dependencies on the reduced predation probability (parameterized by ${\mathcal{P}}_{w}$) — with the ‘weakest’ species being the most abundant or having an average density extremely close to that of the most abundant species. We have also found that if the linear size of the lattice or the total simulation time are not sufficiently large, then the probability of species survival is strongly dependent on initial conditions, in which case a higher probability of survival of the ‘weakest’ species does not generally happen.

Acknowledgements.

P.P.A. acknowledges the support by FCT/MCTES through Portuguese funds (PIDDAC) by the grant UID/FIS/04434/2019. B.F.O. and R.S.T. thank CAPES, Fundação Araucária, and INCT-FCx (CNPq/FAPESP) for financial and computational support.