Rock-paper-scissors models with a preferred mobility direction
RRock-paper-scissors models with a preferred mobility direction
P.P. Avelino,
1, 2, 3
B.F. de Oliveira, and J.V.O. Silva Departamento de Física e Astronomia, Faculdade de Ciências,Universidade do Porto, Rua do Campo Alegre s/n, 4169-007 Porto, Portugal Instituto de Astrofísica e Ciências do Espaço, Universidade do Porto,CAUP, Rua das Estrelas, 4150-762 Porto, Portugal School of Physics and Astronomy, University of Birmingham, Birmingham B15 2TT, United Kingdom Departamento de Física, Universidade Estadual de Maringá,Av. Colombo 5790, 87020-900 Maringá, PR, Brazil
We investigate a modified spatial stochastic Lotka-Volterra formulation of the rock-paper-scissorsmodel using off-lattice stochastic simulations. In this model one of the species moves preferentiallyin a specific direction — the level of preference being controlled by a noise strength parameter η ∈ [0 ,
1] ( η = 0 and η = 1 corresponding to total preference and no preference, respectively) —while the other two species have no preferred direction of motion. We study the behaviour of thesystem starting from random initial conditions, showing that the species with asymmetric mobilityhas always an advantage over its predator. We also determine the optimal value of the noise strengthparameter which gives the maximum advantage to that species. Finally, we find that the criticalnumber of individuals, below which the probability of extinction becomes significant, decreases asthe noise level increases, thus showing that the addition of a preferred mobility direction studied inthe present paper does not favour coexistence. I. INTRODUCTION
Cyclic predator-prey models, also known as rock-paper-scissors (RPS) models [1, 2], have provided insightinto the role of non-hierarchical competition interactionsin the preservation of coexistence, successfully reproduc-ing some of the main properties of simple biological sys-tems with cyclic selection [1, 3, 4]. The simplest mod-els in this class describe the dynamics of a three speciespopulation subject to cyclic interspecific competition (seealso [5–23] for generalizations of the standard RPS modelinvolving additional species and interactions and [24, 25]for recent reviews).The role of mobility has been the subject of many stud-ies which have shown that it may promote or jeopardizebiodiversity (see, e.g., [2]). Although most of these stud-ies considered a uniform isotropic mobility, it has beenshown in [21] that a non-uniform anisotropic mobilitymay affect coexistence in a (four state) May-Leonardformulation of the RPS model using lattice based sim-ulations. In this model the direction of motion for eachindividual was assumed to be the one with a larger den-sity of preys in the surrounding neighborhood. In thiscontext, anisotropic mobility has been shown to have aprofound impact on the dynamics of the population andon the emerging spatial patterns.In the present paper rather than attempting to simu-late the ability of individuals to detect surrounding preyusing their senses, as done in [21], we investigate the po-tential impact of strong correlations between the motionof individuals in a given (large) neighbourhood. This isgenerally observed in nature, specially among the mostdeveloped species, as a result of predator-prey interac-tions. To this end, we shall investigate modified spatialstochastic RPS models in which individuals of one of thespecies move preferentially in a specific direction while the other two species have an isotropic mobility. In thepresent paper we shall consider a (three state) Lotka-Volterra formulation of the RPS model and perform off-lattice stochastic simulations, which, unlike lattice basedones, allow individuals to move in a continuous spatialarea (see [26–30] for numerical studies using simulationsof this type).The outline of this paper is as follows. In Sec. II westart by describing the RPS models with a preferred mo-bility direction that shall be investigated here. In Sec. IIIwe present the results, characterizing the main featuresassociated to the time and spatial dynamics of the popu-lations and showing how the existence of a preferred mo-bility direction of one species investigated in the presentpaper can both have a positive impact on that speciesabundance and a negative impact on the preservation ofcoexistence. Finally, we conclude in Sec. IV.
II. THE MODEL
Here, we consider a modified spatial stochastic Lotka-Volterra formulation of the RPS model. To this end, weshall perform off-lattice simulations in which the individ-uals of the various species (labelled by the letters A , B and C ) are initially randomly distributed in a square-shaped cell of linear size, L , with periodic boundary con-ditions. At the start of the simulations the number ofindividuals of any of the three different species is thesame ( N A = N B = N C = N/
3, where N representsthe total number of individuals). The fractional abun-dance of individuals shall be defined by ρ i = N i /N , with i ∈ A, B, C .At each simulation time step, a single individual (ac-tive) is chosen at random and one action — either mo-bility or predation — is selected, with any of these two a r X i v : . [ q - b i o . P E ] J un AB C
Figure 1: Illustration of non-hierarchical predator-prey inter-actions in RPS models. possible actions carrying the same probability. Figure 1shows the standard scheme of cyclic predation employedin our model: A predates B , B predates C , and C pre-dates A . Whenever predation is selected a circle withradius ‘ is drawn around the active individual and thenearest prey, if it exists, is replaced by an individualof the same species of the active individual — if thereis no prey inside the circle then nothing happens (notethat in a Lotka-Volterra formulation of the RPS model,predation and reproduction take place simultaneously).On the other hand, if mobility is selected and the ac-tive individual belongs to species B (blue) or C (yel-low), then it moves in a random direction with a stepsize ‘ = 2 × − . If the active individual belongs tospecies A (red) then we shall consider two possible mo-bility implementations, which shall be referred as modelI and model II. In model I the random direction asso-ciated to the mobility of species A is restricted to anangle ∆ θ = ξ ( t ) × η around the ˆ x direction, where ξ isa random variable uniformly distributed on [ − π : π ) and η ∈ [0 ,
1] is the noise strength. Note that for η = 1there is no preferred direction of motion while for η = 0the individuals of the species A move always in the sameˆ x direction. Model II, inspired in Vicsek’s model [26], issimilar to model I, except that the ˆ x direction is replacedby the average direction of motion of the individuals in aneighborhood of radius r of the selected individual of thespecies A — for the sake of definiteness we shall consider r = 0 . N actions to take place.Except for the special case discussed in Figure 3, allsimulations performed in the present paper consider ran-dom initial conditions and have a total time span of1 . × generations — the first 10 generations beingdiscarded in the derivation of our main results. The to-tal number of individuals considered in the simulations is N = 3 × , except in the study of extinction probabilitywhere different values in the interval [99 , III. RESULTS
Figure 2 presents three off-lattice simulation snapshotstaken after 5 × generations, considering: (a) the RPS (a)(b) (c) Figure 2: Snapshots of spatial stochastic Lotka-Volterra nu-merical off-lattice simulations of RPS models with a preferredmobility direction. The snapshots were taken after 5 × generations of simulations with 3 × individuals and ran-dom initial conditions, considering: (a) the RPS model withisotropic mobility — η = 1 . η = 0 . η = 0 . model with isotropic mobility ( η = 1 .
0) [30] (b) model Iwith η = 0 . η = 0 .
1. Note that thedistinct spiral patterns present in the top panel of Fig. 2(case (a)) are absent in the two bottom panels (cases (b)and (c)). On the other hand, the characteristic size of thespatial structures seems to decrease with η , being largerin cases (b) and (c) than in case (a). Also, in cases (b)and (c) there are regions with a large density of emptypatches that do not seem to occur in case (a).In order to get a better understanding of the processresponsible for the formation of regions with a high den-sity of empty patches we consider a simulation of model Iwith η = 0 . x direction isresponsible for an initial fast decrease of the blue speciespopulation and for the high density of empty patches inthe boundary region separating the yellow and red speciespopulations. When the blue region gets sufficiently thinit becomes permeable to the passage of individuals of theyellow species which end up engulfing the individuals ofthe red species until it finally becomes extinct.The evolution of the fractional abundances ρ i of thevarious species as a function of time is displayed in Figure4 for a single realization of model I considering η = 1 (toppanel) and η = 0 . η = 1 (top panel) t = 0 t = 10 t = 20 t = 30 Figure 3: Snapshots of the evolution of a single Lotka-Volterraoff-lattice simulation of a RPS model with a preferred mobilitydirection (model I, with η = 0 . N = 3 × ), consideringthe initial configuration shown in the top left panel. The sub-sequent snapshots were taken after 10, 20, and 30 generations(top right, bottom left and bottom right panels, respectively).After 30 generations the red species is extinguished. all the abundances oscillate around the common averagevalue of 1 /
3, while for η = 0 . A , both over its prey and its predator species.Moreover, the characteristic time and amplitude of theoscillations is larger in the case with η = 0 . η = 1. This is associated to the largercharacteristic size of the spatial structures in the formercase. The oscillations are also less sinusoidal for η = 0 . η = 1.The dependence of the characteristic frequency on η may be further quantified using the power spectrum. Thetemporal discrete Fourier transform is defined as ρ ( f ) = N G − X t =0 ρ ( t ) e − πitf/N G , (1)where ρ ( t ) is the fractional abundance of a species, N G =10 generations and f is the frequency. Figure 5 displaysthe power spectrum h| ρ A ( f ) | i for the time evolution ofthe fractional abundance ρ A of the red species for thecases with η = 0 . η = 1 (solid ma-genta line). It shows that maximum of the power occursat a smaller frequency — referred to as the fundamentalfrequency or first harmonic — in the presence of a pre-ferred directional mobility for species A ( η = 0 .
7) thanin the case where the mobility of all species is isotropic( η = 1). On the other hand, the width of the powerspectrum is larger in the former than in the later case. . . . . ρ i tA B C . . . . ρ i tA B C Figure 4: Evolution of the fractional abundances ρ i of thevarious species as a function of time for a single realizationof our model considering η = 1 . η = 0 . . . . . − h | ρ A ( f ) | i f η = 1 . η = 0 . . . . . . . − h | ρ A ( f ) | i f Figure 5: Power spectrum h| ρ A ( f ) | i for the time evolution ofthe fractional abundance ρ A of the red species for the caseswith η = 0 . η = 1 (solid magenta line).The inset highlights the power around the second harmonic. These properties of the power spectrum are a result ofthe larger characteristic time, size and amplitude of theoscillations observed for η = 0 .
7. The inset in Figure 5highlights the power around the second harmonic. Notethat the peak ratio between the first and second harmon-ics is larger for η = 1 than for η = 0 .
7, which is associatedto the less sinusoidal nature of the later case comparedto the former one (see [31, 32] for an application of thepeak ratios in an astrophysical context).We also consider the dependency of the average frac-tional abundances h ρ i i of the various species on the noiseparameter η . To this end, we perform simulations with a . . . . . .
45 0 . N = 30000 h ρ i i η A B C . . . . . .
45 0 . N = 30000 h ρ i i η A B C Figure 6: The average fractional abundances h ρ i i of thespecies i = A , B and C obtained for model I (top panel)and model II (bottom panel) as a function of the noise pa-rameter η . All simulations have a total number of individuals N equal to 3 × and the average was performed over thelast 5 × generations of simulations with a time span equalto 5 . × generations. Notice that in both cases for η < A ) alwayshas an advantage over its predator (species C ). total number of individuals N equal to 3 × — the aver-age abundances are computed considering a time averageover the last 5 × generations of simulations with a timespan equal to 5 . × generations. Figure 6 shows theaverage fractional abundances h ρ i i of the species i = A , B and C obtained for model I (top panel) and model II(bottom panel) as a function of the noise parameter η .Although for maximum noise strength ( η = 1) all specieshave the same average fractional abundance, for η < A ) al-ways has an advantage over its predator (species C ) inboth model I and model II (in model I the maximumadvantage of the red species occurs for a noise strength η ∼ .
7, while for model II it occurs for a value of η closer to unity). The results for model I and model II arequalitatively similar, except for values of η close to unity.We also verified that for r ∼ > . . ∼ < η ∼ < η ∼ < .
35) the preying efficiency of the redspecies is negatively affected by the anisotropic mobil-ity and the blue species becomes the most abundant (atthe expense of the yellow species). Species C and species B are the least abundant for η ∼ < . . ∼ < η ∼ < . . . .
81 10 P N η = 1 . η = 0 . Figure 7: The extinction probability P is depicted as a func-tion of the total number of individuals N . Each point wasobtained from 10 simulations. Notice, that the critical num-ber of individuals, below which the probability of extinctionbecomes significant, decreases as η increases. This result im-plies that the anisotropic mobility considered in the presentpaper does not favour coexistence. respectively.The impact of the existence of a preferred mobility di-rection (for species A ) on coexistence may be studied byestimating the extinction probability P as a function ofthe number of individuals N . This is shown in Fig. 7.Each point was obtained from 10 simulations and, conse-quently, the one-sigma uncertainty in the value of P maybe estimated as [ P (1 − P ) / ] / , with a maximum ofapproximately 0 .
016 for P = 0 .
5. Figure 7 shows that thecritical number of individuals below which the probabil-ity of extinction becomes significant decreases as the levelof anisotropy decreases (or, equivalently, as η increases)— note that the results obtained for η = 1 . η = 0 . N c above which the prob-ability of extinction becomes significant ( P ( N c ) ≡ .
5) isapproximately equal to N c = 9 . × for η = 1 .
0, and N c = 1 . × for η = 0 .
7. This behaviour implies thatthe preferred mobility direction of species A consideredin the present paper does not favour coexistence. IV. CONCLUSIONS
In this work we investigated the dynamics of a pop-ulation of individuals from three different species in thecontext of a spatial stochastic Lotka-Volterra formula-tion of the RPS model where one of the species hasa preferred mobility direction. This has been accom-plished using off-lattice stochastic simulations, startingfrom random initial conditions. We have shown that theanisotropic mobility has a significant impact on the spa-tial patterns which form as a result of the population dy-namics, with the distinct spiral patterns, common in theisotropic case, becoming unrecognizable in the present ofa significant asymmetric mobility. We characterized therelative abundance of the three species as a function ofthe noise level, showing, in particular, that the specieswith asymmetric mobility has always an advantage overits predator. We have also determined the optimal valueof the noise strength parameter which is associated to themaximum advantage of that species relative to the othertwo. Finally, we have found that the threshold numberof individuals, below which the probability of extinctionbecomes significant, decreases as the noise level increases,thus showing that the preferred mobility direction stud-ied in the present paper does not favour coexistence.
ACKNOWLEDGMENTSAcknowledgments
P.P.A. acknowledges the support from FCT throughthe Sabbatical Grant No. SFRH/BSAB/150322/2019 and through the Research Grants No.UID/FIS/04434/2019, UIDB/04434/2020 andUIDP/04434/2020. B.F.O. and J.V.O.S. thank CAPES- Finance Code 001, Fundação Araucária, and INCT-FCx (CNPq/FAPESP) for financial and computationalsupport. [1] B. Kerr, M. A. Riley, M. W. Feldman, and B. J. M.Bohannan, Nature , 171 (2002).[2] T. Reichenbach, M. Mobilia, and E. Frey, Nature ,1046 (2007).[3] B. Sinervo and C. M. Lively, Nature , 240 (1996).[4] B. C. Kirkup and M. A. Riley, Nature , 412 (2004).[5] G. Szabó, A. Szolnoki, and I. Borsos, Phys. Rev. E ,041919 (2008).[6] H. Shi, W.-X. Wang, R. Yang, and Y.-C. Lai, Phys. Rev.E , 030901 (2010).[7] S. Allesina and J. M. Levine, PNAS , 5638 (2011).[8] P. P. Avelino, D. Bazeia, L. Losano, and J. Menezes,Phys. Rev. E , 031119 (2012).[9] P. P. Avelino, D. Bazeia, L. Losano, J. Menezes, and B. F.Oliveira, Phys. Rev. E , 036112 (2012).[10] Y. Li, L. Dong, and G. Yang, Physica A: Statistical Me-chanics and its Applications , 125 (2012).[11] A. Roman, D. Konrad, and M. Pleimling, Journal ofStatistical Mechanics: Theory and Experiment ,P07014 (2012).[12] A. F. Lütz, S. Risau-Gusman, and J. J. Arenzon, Journalof Theoretical Biology , 286 (2013).[13] A. Roman, D. Dasgupta, and M. Pleimling, Phys. Rev.E , 032148 (2013).[14] H. Cheng, N. Yao, Z.-G. Huang, J. Park, Y. Do, andY.-C. Lai, Scientific Reports , 7486 (2014).[15] Y. Kang, Q. Pan, X. Wang, and M. He, Entropy , 284(2016).[16] A. Roman, D. Dasgupta, and M. Pleimling, Journal ofTheoretical Biology , 10 (2016).[17] B. L. Brown and M. Pleimling, Phys. Rev. E , 012147(2017).[18] J. Park, Y. Do, B. Jang, and Y.-C. Lai, Scientific Reports , 7465 (2017).[19] Bazeia, D., Menezes, J., de Oliveira, B. F., and Ramos, J. G. G. S., EPL , 58003 (2017).[20] C. A. Souza-Filho, D. Bazeia, and J. G. G. S. Ramos,Phys. Rev. E , 062411 (2017).[21] P. P. Avelino, D. Bazeia, L. Losano, J. Menezes, B. F.de Oliveira, and M. A. Santos, Phys. Rev. E , 032415(2018).[22] S. Esmaeili, B. L. Brown, and M. Pleimling, Phys. Rev.E , 062105 (2018).[23] J. Park and B. Jang, Chaos: An Interdisciplinary Journalof Nonlinear Science , 051105 (2019).[24] A. Szolnoki, M. Mobilia, L.-L. Jiang, B. Szczesny, A. M.Rucklidge, and M. Perc, Journal of The Royal SocietyInterface , 20140735 (2014).[25] U. Dobramysl, M. Mobilia, M. Pleimling, and U. C. Täu-ber, Journal of Physics A: Mathematical and Theoretical , 063001 (2018).[26] T. Vicsek, A. Czirók, E. Ben-Jacob, I. Cohen, andO. Shochet, Phys. Rev. Lett. , 1226 (1995).[27] X. Ni, R. Yang, W.-X. Wang, Y.-C. Lai, and C. Grebogi,Chaos: An Interdisciplinary Journal of Nonlinear Science , 045116 (2010).[28] X. Ni, W.-X. Wang, Y.-C. Lai, and C. Grebogi, Phys.Rev. E , 066211 (2010).[29] L. You, J. S. Brown, F. Thuijsman, J. J. Cunningham,R. A. Gatenby, J. Zhang, and K. Staňková, Journal ofTheoretical Biology , 78 (2017), ISSN 0022-5193.[30] P. P. Avelino, D. Bazeia, L. Losano, J. Menezes, andB. F. de Oliveira, EPL (Europhysics Letters) , 48003(2018).[31] T. Reinhold and R. Arlt, Astronomy & Astrophysics ,A15 (2015).[32] A. R. G. Santos, M. S. Cunha, P. P. Avelino, R. A. Gar-cía, and S. Mathur, Astronomy & Astrophysics599