Percolation and cooperation with mobile agents: Geometric and strategy clusters
aa r X i v : . [ phy s i c s . s o c - ph ] A ug Percolation and cooperation with mobile agents: geometric and strategy clusters
Mendeli H. Vainstein, ∗ Carolina Brito, † and Jeferson J. Arenzon ‡ Instituto de F´ısica, Universidade Federal do Rio Grande do Sul CP 15051, 91501-970 Porto Alegre RS, Brazil (Dated: August 16, 2018)We study the conditions for persistent cooperation in an off-lattice model of mobile agents playingthe Prisoner’s Dilemma game with pure, unconditional strategies. Each agent has an exclusionradius r P that accounts for the population viscosity, and an interaction radius r int that defines theinstantaneous contact network for the game dynamics. We show that, differently from the r P = 0case, the model with finite sized agents presents a coexistence phase with both cooperators anddefectors, besides the two absorbing phases in which either cooperators or defectors dominate. Weprovide, in addition, a geometric interpretation of the transitions between phases. In analogy withlattice models, the geometric percolation of the contact network (i.e., irrespective of the strategy)enhances cooperation. More importantly, we show that the percolation of defectors is an essentialcondition for their survival. Differently from compact clusters of cooperators, isolated groups ofdefectors will eventually become extinct if not percolating, independently of their size. I. INTRODUCTION
Network reciprocity [1, 2] is a general mechanismresponsible for the development of spatial correlationswithin a viscous population, opening the possibility ofpersistent cooperation. Several specific models have beenproposed showing how these correlations are related tostable groups of cooperating individuals, whose bulk ben-efits of self-defense and mutual support outcompete thesurface exploitation by defectors [3–7]. Although actualexperiments have been performed [8–11], most of ourknowledge comes from these simple models. In partic-ular, a prevailing characteristic in real systems and animportant ingredient for cooperation is the heterogeneouscontact in systems whose interactions are given by com-plex [12, 13] or diluted networks [14]. When we considerthe Prisoner’s Dilemma (PD) dynamics [1] on a dilutedlattice that, albeit heterogeneous, has only short rangeinteractions, intermediate densities present an enhance-ment of cooperation [14–16] and, in the presence of asmall amount of noise, the optimal dilution is closely re-lated to the (random site) percolation threshold for thatlattice [16].Whatever the level of heterogeneity, the contact net-work topology may evolve in time. Although severalrewiring mechanisms can be devised (see Ref. [7] andreferences therein), this may also be accomplished whenthe high viscosity restriction is relaxed and the agentsbecome mobile. Mobility patterns on different scales ofhuman activity, and their far fetched consequences, havebeen studied in recent decades. For example, airplanedisplacement and its connection with disease spread [17],on a global level, can be contrasted with the more localdynamics of pedestrians, crowds or traffic [18, 19]. Ofparticular interest is how the observed patterns can af- ∗ Electronic address: [email protected] † Electronic address: [email protected] ‡ Electronic address: [email protected] fect the outcome of the competition between agents and,in turn, be influenced by it as well. Within the evo-lutionary game theory framework, after several sparse,early attempts to include mobility [20–27], it was onlyrecently that the interest in the combined effects of mo-bility and cooperation in the PD game had a significantincrease. Some level of information processing capabilityis required, for example, when the movement is strat-egy dependent [28, 29] or driven by payoff [30–33], suc-cess [34–37] or the neighborhood composition [38–45].However, the simplest scenario is when mobility is dif-fusive [29, 46–53]. Indeed, as hypothesized in Ref. [54],random mobility may have evolved prior to contingentmobility, allowing bacteria to move away from each otherwhile exploring new resources. Our previous results ona lattice [46, 48, 55] show that even in the framework ofrandom, noncontingent mobility of unconditional agents,diffusion is favorable to cooperation, under rather broadconditions, if velocities are not too high. Analogous con-clusions, attesting the robustness of the results, were alsofound in off-lattice models [45, 49, 56].When diffusion occurs on a lattice and the one agentper site constraint holds, this area exclusion couples thediffusivity of the agents with the free area. This de-pendence on density, on the other hand, is not imme-diately present in off-lattice systems with point-like par-ticles [33, 42, 45, 49, 56]. Moreover, while on a latticethe number of simultaneous interactions is limited by itscoordination number, there is no such restriction on thenumber of point-like particles within the range of interac-tion in off-lattice systems (unless it is explicitly includedas in Ref. [45]). A relevant question concerns the uni-versal effects of such geometrical hindrance on the emer-gence and persistence of cooperation. For example, let-ting the average body size be a coevolving trait, theremay be some evolutive pressure for not too small co-operators because, assuming random diffusion, a groupof small individuals will more easily evaporate from thecluster surface. They should not be too large either and,consequently, not be able to evade defectors and avoidexploitation. Analogously, for defectors, they should nei-ther be too small in order to stay closer to their prey forlonger periods, nor too large so new, more promising re-gions will not be explored. Therefore, one intuitively ex-pects that intermediate, optimal sizes may be beneficialto cooperation and thus be selected for. Another possibleinterpretation for an exclusion zone around each agent isits protected region, and the resources within. Whateverthe interpretation, it is important to better understandthe relevance of area exclusion in these games. As a firstapproximation, we consider an effective radius of exclu-sion, modeled as a hard disk.Here we study, by explicitly taking into account theexcluded area of the agents, the interplay between geom-etry, density and mobility on the capability of a sim-ple model to sustain cooperation and the question ofwhether the transitions in this class of model have a geo-metric interpretation. Although the connection betweenthe threshold of geometric percolation, that is indepen-dent of the game dynamics, and cooperation has alreadybeen reported in Refs. [16, 57, 58], we also explore thegeometry of clusters of cooperators and defectors, andthe connection between their critical properties and thetransition between regions with and without cooperation,thus providing a geometric interpretation of these tran-sitions.
II. THE MODEL
We study an off-lattice model [45, 49, 59] in whichthe N agents living in a square of side L (with periodicboundary conditions) are characterized by an uncondi-tional strategy (cooperate, C, or defect, D) and two in-dependent geometric parameters: an interaction radius r int and a hard disk radius r P to account for excludedarea. The radius r int determines the neighborhood ofeach agent and, as a consequence, its instantaneous con-tact network. The area fraction occupied by the harddisk particles is φ = N πr P /L . We use d ≡ L/ √ N asour length scale. These geometric parameters are illus-trated in Fig. 1. The particular case studied by Meloni et al [49] is recovered in the limit of point-like particles, r P = 0.Initially, N individuals with probability 1/2 of beingeither C or D are randomly placed in such a way thatthere is no overlap between any two individuals i and j ,i.e., their center-to-center distance r ij satisfies r ij ≥ r P .Moreover, they are allowed to randomly diffuse whileplaying the PD game with their neighbors. Two agentsare considered neighbors if r ij < r int . A time step isdefined as a sequence of N attempts of diffusion and acomplete, synchronous round of the PD in which eachof the N agents plays with all its neighbors. During thediffusive part, the position ( x i , y i ) of the center of parti-cle i at time t is updated if there is no overlap between PSfrag replacements r P r int d FIG. 1: (Color online) Geometric parameters of the model,the hard core radius r P and the interaction radius r int . Whilethe former sets an exclusion area around each agent, the latterdefines its contact network. The characteristic length d , asindicated on the right, is defined by dividing the availablearea, L , by the area of the square box around each of the N particles of diameter d , that is, L /d = N . particles in the final position: x i ( t + 1) = x i ( t ) + µr int cos θ i ( t ) y i ( t + 1) = y i ( t ) + µr int sin θ i ( t ) . Each step has a constant size, µr int , and a random ori-entation θ i ( t ) drawn from a uniform distribution in theinterval [ − π, π ]. When µ = 0 there is no mobility andwe consider here that µ is small enough so that jumpsover other agents do not occur. Under mutual coopera-tion (defection), both receive payoff R ( P ) as a reward(punishment); if one cooperates and the other defects,then the latter receives T (temptation) and the former, S .To characterize the PD, the following inequalities shouldhold T > R > P ≥ S and 2 R > T + S . In particu-lar, we use R = 1, P = S = 0 and T >
1, a commonparametrization known as the weak form of the PD game.The evolution follows the finite population analog of thereplicator dynamics [13]. Each individual i , after accu-mulating the payoff from all combats, randomly choosesa neighbor j with whom to compare their respective pay-offs P i and P j . If P i ≥ P j , then i maintains its strategy.On the other hand, if P j > P i , i will adopt the strategyof j with probability proportional to the payoff differenceΠ ij = P j − P i max { k j , k i } T , (1)where k i and k j are the number of neighbors of i and j ,respectively. Under this update rule, the total number ofindividuals is kept constant.Most of our results are for N = 32 , T = 1 . µ = 0 . L = 1. We then check the robustness of the model bytesting finite size effects with up to N = 128 particles, aswell as the dependence on T and µ . Averages are takenover 100 or more samples. III. COOPERATION AND PERCOLATION
Two macroscopic asymptotic quantities, once aver-aged, are used to characterize the system: the fraction ofcooperators ρ C (those, among the N agents, that cooper-ate) and the fraction f C ≤ ρ C of initial conditions whoseevolution ends in the absorbing state ρ C = 1. Their dif-ference, ρ C − f C , is a measure of the coexistence of bothstrategies. Four regimes are present in the time evolution,as seen by the behavior of ρ C ( t ) in Fig. 2. As is often thecase for this class of model, there is an initial drop in thefraction of cooperators from ρ C (0) = 1 /
2, since small co-operator clusters are easily preyed on in the beginning ofthe simulation. As ρ C ( t ) approaches its minimum valueat t ∼ , fluctuations may lead to extinctions in finitesize systems. Away from the minimum, the survivingclusters of cooperators resume growth. These two initialregimes are quite independent of the occupied area frac-tion, as indicated in Fig. 2 by the close proximity of allcurves up to t ∼ . In the third regime, ρ C ( t ) attains aplateau where it stays, indicating persistent coexistenceof both strategies, if φ is large enough (in the case ofFig. 2, the threshold is at φ ≃ . φ is belowthe threshold, after the stasis period on the plateau, thesystem enters the last regime in which cooperators takeover the system, that is, ρ C ( ∞ ) = 1. The time to reachthis asymptotic state seems to diverge as the thresholdarea fraction is approached from below. φ ≅ PSfrag replacements t ρ C ( t ) FIG. 2: (Color online) Average fraction of cooperators as afunction of time (in Monte Carlo steps). N = 32 , T = 1 . r int = √ . d , and µ = 0 .
01 for various area fractions φ . Below φ ≃ .
37, cooperators eventually invade the whole system.As φ approaches this threshold, the time spent close to thecritical plateau at ρ C ≃ .
85 also increases.
Fig. 3 summarizes our most important results, show-ing, in the stationary regime, the average fraction of co-operators ρ C as a function of both r int /d and φ . The linesare independent measures of percolative properties thatwill be explained below. Several regimes may be identi-fied: two absorbing phases in which all agents eventuallybecome either defectors ( ρ C = f C = 0, labeled ‘D’) orcooperators ( ρ C = f C = 1, ‘C’) and two coexistence ones(0 < ρ C < f C = 0, ‘C and D’). Finite systems may also present a bi-stable phase, in which all initialconditions lead to one of the absorbing states. In thiscase, although the average cooperativeness still obeys0 < ρ C <
1, it differs from the coexistence state since f C = ρ C . However, by increasing the system size, theprobability of becoming dominated by defectors goes tozero inside the ’C’ region in Fig. 3 and, taking this intoaccount (see Fig. 5), we already properly labeled it. Wenow concentrate on the results for a finite system with N = 32 particles, and discuss the finite size effects atthe end of this section.There are two limits of the diagram that have trivialresults. For very large values of r int , very distant particlesinteract, increasing the number of contacts and decreas-ing the effects of spatial correlation. Thus, we recoverthe mean field result in which all agents become defec-tors. This trivial region of the phase diagram was notexplored. In the other limit, when r int < r P , the hardcore prevents any interaction between the agents and thefraction of cooperators remains equal to the initial one, ρ C = 1 / f C = 0. This is the trivial coexistenceregion, labeled ‘C and D’ on the left of Fig. 3, abovethe line φ = ( π/ r int /d ) , corresponding to 2 r P = r int .Immediately to the right of this line, interaction, albeitweak, is possible and clusters are formed. However, theyare small and do not favor cooperation, therefore ρ C = 0,and this all defector phase is labeled ‘D’. As we discussbelow, the transitions between the other phases have ge-ometric origins and are closely related to the percolatingproperties of the contact networks. FIG. 3: (Color online) Phase diagram in the plane r int /d and φ for N = 32 , µ = 0 .
01 and T = 1 .
1. The color code isthe average fraction of cooperators, ρ C . Besides the phasescorresponding to the absorbing states, ‘C’ ( ρ C = f C = 1) and‘D’ ( ρ C = f C = 0), there are two coexistence phases, ‘C andD’. Geometric percolation.
For the area fractions consid-ered here, non-trivial cooperation first appears around r int ≈ d . For small φ , the transition between phases Dto C corresponds to a change in stability of the absorb-ing state, from the defector to the cooperator dominatedphase, while for larger φ , roughly φ > .
25, the emergentphase is one in which cooperators and defectors coex-ist, ’C and D’. This onset of cooperation is strongly cor-related with the appearance of a geometric percolatingcluster, indicated by the steep line in Fig. 3 at the pointwhere the probability of finding a percolating cluster is50%. This is a purely geometric problem of disks withboth an inner hard core in addition to a soft, penetrableregion, and thus independent of the game dynamics. Inother words, the network of contacts of the percolatingcluster spans the whole length of the system. For φ = 0,our result is consistent both with the percolation thresh-old obtained numerically in Refs. [60, 61], and the exactbounds of Ref. [62]. For finite φ , the threshold is slightlysmaller than for φ = 0, since as r P increases, the overlapbetween the disks decreases and percolation is attainedwith a smaller r int . Percolating clusters of defectors.
Besides the clustersof particles irrespective of their strategies, we also con-sider the geometry of clusters composed only by defec-tors (that are, in turn, intimately connected with thegeometry of cooperator clusters), which depends on theparticular strategy evolving dynamics. Fig. 4 shows theprobability of percolation of D clusters as a function oftime, P D ( t ), for r int = √ . d . Notice that in this regionthere always is a percolating geometric cluster. Initially,as the fraction of cooperators decreases towards the min-imum, there is a sea of defectors that obviously perco-lates and all curves overlap at P D = 1. It is only whenthe curves of ρ C ( t ), for different values of φ , start to sep-arate, around t ∼ , that P D ( t ) starts decreasing. Theasymptotic probability of there being a percolating clus-ter of defectors attains a limiting value, as shown in theinset of Fig. 4, from which the threshold can be obtained.For the particular value of r int shown in this figure, whenthe area fraction is below the threshold at (roughly) 0.37, P D ( ∞ ) = 0 and, as can be seen in Fig. 2, ρ C ( ∞ ) = 1. Thedotted transition line in the phase diagram (Fig. 3) is ob-tained in the same manner: it is the asymptotic value of φ at which P D ( ∞ ) = 0 as a function of r int . In the C re-gion, finite size fluctuations sometimes lead to the all-Dstate ( r int /d ' φ . This transition line suggests an important ingredientto understand the coexistence between cooperators anddefectors: only under the presence of a percolating seaof defectors is that a stable coexistence between coopera-tors and defectors is possible. In other words, differentlyfrom compact clusters of cooperators, isolated groups ofdefectors either grow and percolate or eventually becomeextinct. Finite size effects.
We now discuss how the systemsize can affect the phase diagram, Fig. 3. For finite sizedsystems, a small region inside the ’C’ phase has a bi-stable equilibrium, in which all initial conditions leadto an absorbing state, either ρ C = 0 or 1. The all-Dstate is due to the fact that the population of coopera-tors becomes quite small during the initial drop in thefirst generations and, therefore, sensitive to fluctuations PSfrag replacements tφ P D ( t ) P D ( ∞ ) FIG. 4: (Color online) Probability of percolation of D clustersas a function of time, P D ( t ), for several area fractions φ andthe same parameters of Fig. 2. Inset: Asymptotic limit of P D ( t ) as a function of φ . The linear fit indicates that thepercolation probability goes to zero at φ ≃ .
37. The smalldeviation seen close to this point is due to the long time ofconvergence. which occasionally cause extinctions. Fig. 5 shows f C fortwo different values of φ and several system sizes. Forpoint particles ( φ = 0), the absorbing all-C state regiongrows as the system size increases. On the other hand,for ( φ ≈ . N = 32 , which explains why this phase is nothomogeneously painted in Fig. 3. Notice that when thesystem is bi-stable, the average fraction of cooperators isnot a good measure since it represents neither one of thefinal states [63]. On the other hand, in the coexistencestate, both strategies are present in the asymptotic stateand, while 0 < ρ C < f C = 0. For the r P = 0 casestudied by Meloni et al [49], the φ = 0 line in the phasediagram, there is no coexistence phase and the systemeventually enters an absorbing state. Robustness against mobility and temptation.
We fi-nally consider the robustness of our results when the mo-bility µ and the temptation T are varied, Fig. 6. Thetop panel shows several values of µ , with the temptationfixed at T = 1 .
1. Whatever the velocity, no cooperationis possible for r int ≤ d due to the absence of a percolatinggeometric cluster. Notwithstanding, for small mobilities,existence of cooperators is possible in a wide range of r int .When comparing the low mobility case with the one withimmobile agents ( µ = 0), it can be seen that there is animprovement only for low values of the r int /d ratio. Forhigh values of this ratio, the curves overlap, which is ex-pected, since the agents have a large neighborhood thatis minimally perturbed by the small random movements. PSfrag replacements r int /d f C φ = 0 φ ≃ . FIG. 5: (Color online) Fraction of initial states that go tothe ρ C = 1 absorbing state as a function of r int /d for severalsystem sizes and two area fractions φ = 0 and 0.28. For largevalues of r int /d , as f C increases for larger system sizes, theexistence of the ρ C = 0 state for point particles ( φ = 0) is dueto finite size fluctuations. For φ ≃ .
28 the behavior is theopposite and the transition becomes sharper when the systemsize increases.
On the other hand, for low values of the ratio, the contactrange is small and is more affected by the diffusion. Ran-dom movements lead to the evaporation of cooperativeclusters and, as the velocity increases, cooperation lev-els decrease until a threshold above which it is no longerpossible. This effect of cooperation enhancement drivenby a low mobility is in accordance with previous simula-tions on a lattice [46, 48, 49, 55]. The bottom panel ofFig. 6 shows the effect of the parameter T . As expected,as the temptation to defect increases, the fraction of co-operators decreases. PSfrag replacements r int /d ρ C µT FIG. 6: (Color online) Average asymptotic fraction of coop-erators as a function of r int /d for φ ≃ .
28 and several valuesof µ (top), including µ = 0 (empty symbols), and T (bottom).In both cases, N = 32 ; for the top panel, T = 1 . µ = 0 . IV. DISCUSSION AND CONCLUSION
We presented numerical results for an off-lattice modelof mobile agents playing the PD game while randomlymoving across a closed region. This model combinesingredients found in two distinct models for such sys-tems: the excluded area found in lattice simulations andan available continuous space. We arrive at two impor-tant results: first, if the agents are not point-like, a non-trivial coexistence phase with cooperators and defectorsbecomes possible (besides the trivial one when the exclu-sion range is larger than the interacting radius); second,it is possible to geometrically interpret, in terms of per-colation, the observed transitions. An important role isplayed not only by geometric percolation, irrespective ofthe game dynamics, but also by the percolative proper-ties of defector clusters, whose threshold depends on thedetails of the game.When the agents interaction is prevented by the hardcore, trivial coexistence between cooperators and defec-tors is possible. Beyond that region of the phase diagram,no cooperation is possible in the absence of a percolatinggeometric cluster. Random mobility provides an evapo-rating mechanism for groups of cooperators what is detri-mental to cooperation as isolated cooperators are easilypreyed on. However, in the presence of a percolatingclusters, it is easier for a detaching cooperator to be incontact with another cluster and be protected. Whenhard cores are included, movements are hindered andthe agents spend some time rattling around the sameregion, while large displacements become less probable,what increases the correlation among agents and benefitscooperation even further. This localization is probablythe mechanism responsible for the coexistence betweencooperators and defectors that is not present for φ = 0.Interestingly, the presence of defectors is only possible ifthey form a percolating cluster and no finite cluster ofdefectors is stable: they either grow, merge with othersand span the whole lattice or the isolated cluster becomesextinct. This is our main result: finite size cooperatorsand defectors, whose hard core is an effective, averagedinteraction restraining their movements, are able to co-exist over a broad region of the phase diagram only ifdefectors are organized in a interconnected cluster, a per-colating sea of defectors. A similar effect was found forthe public goods game played on a lattice with emptysites and no mobility [16]. It would be interesting to in-vestigate whether this condition for coexistence betweencooperators and defectors also occurs in lattice models,where excluded area is inherent to the formulation of theproblem.In this manuscript we focused on the particular homo-geneous case of equal sizes and equal velocities for bothcooperators and defectors. Following Refs. [45, 55], itis essential to explore the whole ( S, T ) parameter spaceand the dependence on the chosen dynamic rule in orderto check the robustness of cooperation. Furthermore,several extensions are possible. For example, velocitiesmay not be constant [29] and depend on the neighbor-hood [45] or strategy. The hard core radius may alsocorrelate with strategy, r C and r D for cooperators anddefectors, respectively. In particular, if individuals co-evolve with mutations: is there an optimal equilibriumradii ratio to which the system converges to or, instead,a permanent arms race? What are the effects of havingsize dispersion? If velocity and size coevolve along withstrategies, defectors may become small and fast while co-operators become large and slow. Finally, what happensif there is a fraction of fundamentalists (both coopera-tors and defectors or just defectors) whose strategies orpositions never change? In all these cases, it is importantto study also the geometric properties of the interfacesbetween cooperator and defector clusters since these arethe places where all strategy flips occur. From a more physical perspective, it would be interesting to find out,for each transition line, which is the dynamical universal-ity class that it belongs to [58]. These and other relevantquestions are being considered. Acknowledgments
Research partially supported by the Brazilian agenciesCNPq, CAPES and Fapergs. JJA also thanks the INCT-Sistemas Complexos (CNPq) for partial support. Wethank the supercomputing laboratory at UniversidadeFederal da Integra¸c˜ao Latino-Americana (LCAD/Unila),where the simulations were run, for computer time. [1] R. Axelrod,
The Evolution of Cooperation (BasicBooks,New York, 1984).[2] M. A. Nowak and R. M. May, Nature , 826 (1992).[3] M. Doebeli and C. Hauert, Ecology Letters , 748 (2005).[4] M. A. Nowak, Science , 1560 (2006).[5] G. Szab´o and G. F´ath, Phys. Rep. , 97 (2007).[6] C. P. Roca, J. A. Cuesta, and A. S´anchez, Phys. LifeRev. , 208 (2009).[7] M. Perc and A. Szolnoki, BioSystems , 109 (2010).[8] A. Traulsen, D. Semmann, R. D. Sommerfeld, H.-J.Krambeck, and M. Milinski, Proc. Natl. Acad. Sci. ,2962 (2010).[9] D. G. Rand, S. Arbesman, and N. A. Christakis, Proc.Nat. Acad. Sci. , 19193 (2011).[10] C. Gracia-L´azaro, A. Ferrer, G. Ruiz, A. Tarancon, J. A.Cuesta, A. S´anchez, and Y. Moreno, Proc. Nat. Acad.Sci. , 12922 (2012).[11] J. Grujic, T. Rohl, D. Semmann, M. Milinski, andA. Traulsen, PLOS One , e47718 (2012).[12] F. C. Santos and J. M. Pacheco, Phys. Rev. Lett. ,098104 (2005).[13] F. C. Santos, J. M. Pacheco, and T. Lenaerts, Proc.Natl. Acad. Sci. , 3490 (2006).[14] M. H. Vainstein and J. J. Arenzon, Phys. Rev. E ,051905 (2001).[15] M. Lin, N. Li, L. Tian, and D.-N. Shi, Physica A ,1753 (2010).[16] Z. Wang, A. Szolnoki, and M. Perc, Nature Sci. Rep. ,369 (2012).[17] V. Colizza, A. Barrat, M. Barth´elemy, and A. Vespig-nani, Proc. Nat. Acad. Sci. , 2015 (2006).[18] D. Helbing, I. Farkas, and T. Vicsek, Nature , 487(2000).[19] S. Bouzat and M. N. Kuperman, Phys. Rev. E , 032806(2014).[20] L. A. Dugatkin and D. S. Wilson, Am. Nat. , 687(1991).[21] M. Enquist and O. Leimar, Anim. Behav. , 747 (1993).[22] R. Ferri`ere and R. E. Michod, Am. Nat. , 692 (1996).[23] S. J. Majeski, G. Linden, C. Linden, and A. Spitzer,Complexity , 16 (1999).[24] I. M. Hamilton and M. Taborsky, Proc. R. Soc. B , 2259 (2005).[25] J. C. Koella, Proc. R. Soc. B , 1979 (2000).[26] C. A. Aktipis, J. Theor. Biol. , 249 (2004).[27] J.-F. Le Galliard, F. Ferri`ere, and U. Dieckmann, Am.Nat. , 206 (2005).[28] L.-L. Jiang, W.-X. Wang, Y.-C. Lai, and B.-H. Wang,Phys. Rev. E , 036108 (2010).[29] H. Cheng, H. Li, Q. Dai, Y. Zhu, and J. Yang, New J.Phys. , 123014 (2010).[30] H.-X. Yang, Z.-X. Wu, and B.-H. Wang, Phys. Rev. E , 065101(R) (2010).[31] H. Cheng, Q. Dai, H. Li, Y. Zhu, M. Zhang, and J. Yang,New J. Phys. , 043032 (2011).[32] H. Lin, D.-P. Yang, and J. Shuai, Chaos, Solitons &Fractals , 153 (2011).[33] Y.-T. Lin, H.-X. Yang, Z.-X. Wu, and B.-H. Wang, Phys-ica A , 77 (2011).[34] D. Helbing and W. Yu, Proc. Nat. Acad. Sci. , 3680(2009).[35] W. Yu, Phys. Rev. E , 026105 (2011).[36] Y. Liu, X. Chen, L. Zhang, F. Tao, and L. Wang, Chaos,Solitons & Fractals , 1301 (2012).[37] P. Buesser, M. Tomassini, and A. Antonioni, Phys. Rev.E , 042806 (2013).[38] Z. Chen, J. Gao, Y. Cai, and X. Xu, Physica A ,1615 (2011).[39] J. Zhang, W.-Y. Wang, W.-B. Du, and X.-B. Cao, Phys-ica A , 2251 (2011).[40] C. Zhang, J. Zhang, F. J. Weissing, M. Perc, G. Xie, andL. Wang, PloS one , e35183 (2012).[41] R. Cong, B. Wu, Y. Qiu, and L. Wang, PLOS One ,e35776 (2012).[42] H. Cheng, Q. Dai, H. Li, X. Qian, M. Zhang, andJ. Yang, Eur. Phys. J. B , 127 (2013).[43] G. Ichinose, M. Saito, H. Sayama, and D. S. Wilson, Sci.Rep. , 1 (2013).[44] C. Zhang, J. Zhang, and G. Xie, Chaos, Solitons & Frac-tals , 103 (2014).[45] A. Antonioni, M. Tomassini, and P. Buesser, J. Theor.Biol. , 40 (2014).[46] M. H. Vainstein, A. T. C. Silva, and J. J. Arenzon, J.Theor. Biol. , 722 (2007). [47] M. Droz, J. Szwabi´nski, and G. Szab´o, Eur. Phys. J. B , 579 (2009).[48] E. A. Sicardi, H. Fort, M. H. Vainstein, and J. J. Aren-zon, J. Theor. Biol. , 240 (2009).[49] S. Meloni, A. Buscarino, L. Fortuna, M. Frasca,J. G´omez-Garde˜nes, V. Latora, and Y. Moreno, Phys.Rev. E , 067101 (2009).[50] S. Suzuki and H. Kimura, J. Theor. Biol. , 42 (2011).[51] P. E. Smaldino and J. C. Schank, Theor. Pop. Biol. ,48 (2012).[52] H. Yang and B. Wang, Chin. Sci. Bull. , 3693 (2011).[53] A. Gelimson, J. Cremer, and E. Frey, Phys. Rev. E ,042711 (2013).[54] Y. Wei, X. Wang, J. Liu, I. Nememan, A. H. Singh,H. Weiss, and B. R. Levin, Proc. Natl. Acad. Sci. ,4047 (2011).[55] M. H. Vainstein and J. J. Arenzon, Physica A , 145 (2014).[56] D. Amor and J. Fort, Phys. Rev. E , 066115 (2011).[57] Z. Wang, A. Szolnoki, and M. Perc, Phys. Rev. E ,037101 (2012).[58] H.-X. Yang, Z. Rong, and W.-X. Wang, New J. Phys. , 013010 (2014).[59] M. A. Nowak, S. Bonhoeffer, and R. M. May, Int. J. ofBifurcation and Chaos , 33 (1994).[60] J. Quintanilla, S. Torquato, and R. M. Ziff, J. Phys. A:Math. Gen. , L399 (2000).[61] S. Mertens and C. Moore, Phys. Rev. E , 061109(2012).[62] P. Balister, B. Bollob´as, and M. Walters, Rand. Struct.Algor. , 392 (2005).[63] M. N. Kuperman and S. Risau-Gusman, Phys. Rev. E86