Von-Neumann’s and related scaling laws in Rock-Paper-Scissors type models
We introduce a family of Rock-Paper-Scissors type models with symmetry ( is the number of species) and we show that it has a very rich structure with many completely different phases. We study realizations which lead to the formation of domains, where individuals of one or more species coexist, separated by interfaces whose (average) dynamics is curvature driven. This type of behavior, which might be relevant for the development of biological complexity, leads to an interface network evolution and pattern formation similar to the ones of several other nonlinear systems in condensed matter and cosmology.
The mechanisms leading to the enormous biodiversity observed in nature are still not fully understood. Rock-Paper-Scissors (RPS) type models incorporate some of the crucial ingredients associated with the dynamics of a network of competing species and they have been used as a powerful tool in the study of complex biological systems. In its simplest form, the RPS game describes the evolution of 3 species which cyclically dominate each other Kerr et al. (2002); Reichenbach et al. (2007); Shi et al. (2010) (see also Volterra ((1931); Lotka (1920) for the pioneer work by Lotka and Volterra). The basic interactions are Motion, Reproduction and Predation but generalizations, incorporating more than 3 species and/or new interactions between them, have been proposed in the literature Szabó et al. (2007); Peltomaki and Alava (2008); Szabo et al. (2008); Wang et al. (2010, 2011); Hawick (2011a, b); Zia (2010); Dobrinevski and Frey (2012). Particularly, in Szabó et al. (2007) segregation processes and phase transitions in predator-prey models with an even number of species have been investigated.
RPS type models naturally lead to the formation of complex spatial patterns observed in some biological systems Kerr et al. (2002). Complex spatial structures also arise in many other systems. For example, interfaces in ideal soap froths and grain growth have a velocity proportional to the mean curvature at each point, which is at the core of the Von Neumann’s law von Neumann (1952) for the evolution of the area of individual domains. The evolution with time of the characteristic scale of such interface networks obeys the scaling law , leading to the formation of cellular patterns of increasing size Glazier (1989); Stavans and Glazier (1989); Glazier and Weaire (1992); Flyvbjerg (1993); Monnereau and Vignes-Adler (1998); Weaire and Hutzler (2000); Kim et al. (2006); Arenzon et al. (2007). Despite the different dynamical scaling laws, the evolution of interface networks in a cosmological context may also generate similar spatial patterns to the ones observed in soap froths Avelino et al. (2008, 2011).
Consider a family of models where individuals of various species are distributed on a square lattice of size at some initial time. The different species are labelled by the number (or ) with , and we make the cyclic identification where is an integer. The number of individuals of the species will be denoted by . In addition to individuals, there are also empty spaces which shall be denoted by . In this paper, except if stated otherwise, we shall assume that the initial distribution is random and that the number densities of the various species are all identical at the initial time. The number density of empty spaces is given by and its initial value is equal to in all the network simulations described in this paper. At each time step a random individual (active) interacts with one of its four nearest neighbors (passive). The unit of time is defined as the time necessary for interactions to occur (one generation time). The possible interactions are classified as Motion
where with if is even or if is odd. We shall denote the corresponding probabilities by (Motion), (Reproduction), (left-handed Predation) and (right-handed Predation). This family of models has a symmetry if , , and for all . Figure 1 represents the different predation interactions in a model with species having symmetry. Throughout this paper, we shall assume that and that the symmetry is preserved with the following interacting probabilities: , (or zero, if the passive is not ), (or zero, if the passive is or if the passive is not a prey of the active). The symmetry is, in general, not associated with curvature driven dynamics. For example, the standard RPS model has a symmetry but the dynamics of the spatial patterns is not curvature driven in this model. In fact, we shall show that the dynamics is curvature driven only in realizations which result in the formation of domains, where individuals of one or more species coexist, separated by interfaces whose dynamics is controlled by interactions of equal strength between competing species.
Let us start by introducing a simple model with having (model I). This model has 2 competing species which tend to distribute themselves into separate domains bounded by interfaces where most of the action occurs. Predation happens mainly at the interfaces whose thickness is directly related to the mobility. This is clearly shown on Fig. 2 where a snapshot of the network evolution of model I is presented. On the left panel each species is represented by a different color while on the right panel the black dots represent the distribution of empty spaces. Fig. 2 clearly shows that the empty spaces, which are a result of Predation between individuals of the competing species, are located at the domain borders. The average thickness of the interfaces does not change with time.
For an interface with curvature radius and thickness , the average number of attacks per unit time from individuals outside the border is proportional to the outer interface length (proportional to ) while the average number of attacks from individuals inside the border is proportional to the inner interface length (proportional to ). The difference between the average number of attacks per unit of time from outside and inside the border is proportional to . On the other hand, the number of attacks necessary to modify the domain radius by is proportional to the interface length (). This implies that the value of the velocity of the interface is on average proportional to its curvature (, where is a positive constant), which is typical of non-relativistic interfaces in condensed matter Avelino et al. (2011).
The average evolution of the area of a compact simply connected domain with no vertices is then given by , where is the infinitesimal interface arc length and a dot denotes a derivative with respect to . Note that the domain area can be calculated at any given time by counting the number of individuals inside the domain. If the domain is compact, but not necessarily simply connected, then
where represents the area of a compact domain with genus and no vertices. Eq. (1) implies that the area decreases (increases) with time depending on whether (). The genus dependency accounts for the contribution of the decrease of the area of each hole to the growth of the area of the domain. If then, on average, the evolution of the area with time is given by
where is the time of collapse. The average time evolution of the area of an initial circular domain with is illustrated in Fig. 3. The domain area was calculated by counting the number of individuals of the species inside the domain. We verified that a nearly identical result is obtained using the relation . The solid red line represents the average result in an ensemble of 43 simulations using model I and the dashed blue line shows the theoretical evolution given by Eq. (2), with calculated as the median collapse time. Fig. 3 shows that the agreement between the analytical and numerical results is very good.
In the case of an interface network without junctions the evolution of the total number of domains, , with time can be obtained using Eq. (1). If the fraction of the total number of domains with genus is a constant () then Glazier (1989)
Here is the average area of a compact domain with genus , is the average domain area, is the area of the entire system and is the total area occupied by domains with genus . If both and are time independent, then . The characteristic length of the network defined by evolves as
Fig. 4 shows the evolution of the number of empty spaces with time for model I. The cyan dots (light grey in black and white) represent the results of 20 interface network simulations while the solid line represents the average. The scaling exponent defined by is at one-sigma. The number of domains is equal to the ratio between the total area and the average domain area . It is also proportional to the ratio between the total interface length and the average domain perimeter which is proportional to . Hence, implying that . On the other hand, the average thickness of the interfaces does not change with time and consequently the total interface length is proportional to . Hence, is inversely proportional to the number of empty spaces () and consequently one would expect an average scaling solution with , with . The fact that the numerical value of is very close to the theoretical one demonstrates that the interface network evolution is already attaining the expected scaling regime. The increase with time of the dispersion between the values of for the various simulations is associated with the growth of the characteristic length scale of the network (see Avelino et al. (2005) for a detailed analytical discussion).
If then the evolution of the area of a single domain with genus and vertices is given by
where represent each of the discontinuous angle changes at the vertices. If the interface network has only Y-type junctions which meet at an angle of , then one obtains the Von Neumann’s law von Neumann (1952)
implying that the area domains with () decreases (increases) proportionally to time. The evolution of interface networks in RPS type models is not deterministic and consequently the Von Neumann’s law can only be valid on average. Fig. 5 shows the evolution of a Y-type junction in model IV where for all (this model will be described later in more detail), starting from the initial configuration shown in the left panel. The right panel represents the most frequent species at each lattice point from to , showing that the average angles at the vertex all tend to an average value of .
The evolution of the total number of domains with time, in the case of an interface network with Y-type junctions, can be obtained using Eq. (6). If the fraction of the total number of domains with edges is a constant () then
Here is the average area of a domain with edges, is the average domain area, is the area of the entire system, is the total area occupied by domains with edges and the function accounts for the fact that the collapse of an individual domain might lead to the merger of some of the surrounding domains ( if all domains are of different types and otherwise). If the interface network is in a scaling regime with time-independent , and , then again . Consequently, the scaling law given in Eq. (4) also applies to interface networks with Y-type junctions.
Now we shall demonstrate that the curvature driven interface dynamics of model I is common in RPS type models. Let us consider a model with . By taking and (model II) one ensures that the domains which appear in the simulations are populated with one of the following non-interacting species partnerships: and . The interactions (of equal strength) between the 2 different partnership domains occur mainly at the border where a large number of empty spaces are continuously being generated (see video vid (a) and Fig. 6 (upper left)). If one considers the case with , (model III) then there are 3 possible partnership domains containing non-interacting species: , and thus leading to an interface network with Y-type junctions (see video vid (b) and Fig. 6 (upper right)). In the case with for (model IV), there are no partnerships since every species is linked in a bidirectional way to all other species thus generating an interface between any 2 given domains. This gives rise to an interface network with Y-type junctions and 6 different domain types (see video vid (c) and Fig. 6 (lower left)). If one now takes and (model V) then there are 2 possible domains defined by and but due to the non-zero unidirectional probability of interaction between species in the same domain, spiral patterns do form (see video vid (d) and Fig. 6 (lower right)). Inside the 2 different domains the interactions are those of the standard RPS model with 3 species. In the boundary there is an (average) equilibrium between the predation probabilities from individuals on either side of the wall (for example, individuals of the species 1 predate and are predated by individuals of the species 2 and 6). In models IV and V the colors light blue, dark blue, red, maroon, green and yellow represent species 1 to 6, respectively. This study represents a significant extension with respect to previous work and, to the best of our knowledge, it is the first time that models III, IV and V have been studied in the literature. Note that model V is very different from the 6 species extension of the 5 species Rock-Paper-Scissors-Lizard-Spock model proposed in Hawick (2011b), defined by and (see video vid (e)), whose dynamics is not curvature driven.
The scaling parameters calculated numerically for all models are: (model I), (model II), (model III), (model IV) and (model V). In the case of model V, the empty spaces also appear associated to the spiral patterns. Hence, the parameter for this model was obtained using only the empty spaces which have as some of the four immediate neighbors individuals from the 2 groups: and . In all cases the scaling parameter is reasonably close to the expected value . The deviations can be attributed to the finite size and dynamical range of the simulations. We have verified that the results obtained for models I, II, III, and IV are weakly dependent on the value of , as long as the thickness of the interfaces is much smaller than the box size. In the case of model V we have found that the network evolution is curvature dominated only if . For smaller values of there are other effects which have a significant impact on the network dynamics, such as local partnerships along the borders of the domains Szabó et al. (2007) (for example, between species 1 and 4 which do not predate each other). These effects are outside the scope of the present paper and will studied in detail in future work. In summary, we have shown that curvature driven interface dynamics, analogous to that observed in other physical systems in condensed matter and cosmology, is common in the family of RPS type models investigated in this paper and may be crucial to the understanding of biological complexity.
We thank Breno Oliveira for computational support and useful discussions. This work is partially supported by FCT-Portugal, CAPES and CNPq-Brazil.
- Kerr et al. (2002) B. Kerr, M. A. Riley, M. W. Feldman, and B. J. M. Bohannan, Nature 418, 171 (2002).
- Reichenbach et al. (2007) T. Reichenbach, M. Mobilia, and E. Frey, Nature 448, 1046 (2007).
- Shi et al. (2010) H. Shi, W.-X. Wang, R. Yang, and Y.-C. Lai, Phys. Rev. E 81 (2010).
- Volterra ((1931) V. Volterra, Lecons dur la Theorie Mathematique de la Lutte pour la Vie (Gauthier-Villars, Paris, (1931)), ed.
- Lotka (1920) A. J. Lotka, Journal of the American Chemical Society 42, 1595 (1920).
- Szabó et al. (2007) G. Szabó, A. Szolnoki, and G. A. Sznaider, Phys. Rev. E 76, 051921 (2007).
- Peltomaki and Alava (2008) M. Peltomaki and M. Alava, Phys. Rev. E 78 (2008).
- Szabo et al. (2008) G. Szabo, A. Szolnoki, and I. Borsos, Phys. Rev. E 77 (2008).
- Wang et al. (2010) W.-X. Wang, Y.-C. Lai, and C. Grebogi, Phys. Rev. E 81 (2010).
- Wang et al. (2011) W.-X. Wang, X. Ni, and C. Lai, Y.-C. an Grebogi, Phys. Rev. E 83 (2011).
- Hawick (2011a) K. Hawick, Proceedings of the IASTED International Conference on Modelling and Simulation pp. 129–136 (2011a).
- Hawick (2011b) K. A. Hawick, CSTN Computational Science Technical Note Series (2011b), URL http://www.massey.ac.nz/~kahawick/cstn/129/cstn-129.pdf.
- Zia (2010) R. K. P. Zia, General properties of a system of species competing pairwise (2010), eprint arXiv:1101.0018.
- Dobrinevski and Frey (2012) A. Dobrinevski and E. Frey, Phys. Rev. E 85, 051903 (2012).
- von Neumann (1952) J. von Neumann, in Metal Interfaces (1952), pp. 108–110.
- Glazier (1989) J. A. Glazier, Phd thesis, University of Chicago (1989), URL http://biocomplexity.indiana.edu/jglazier/cv.php?pg=2.
- Stavans and Glazier (1989) J. Stavans and J. A. Glazier, Phys. Rev. Lett. 62, 1318 (1989).
- Glazier and Weaire (1992) J. A. Glazier and D. Weaire, J. Phys. Cond. Mat. 4, 1867 (1992).
- Flyvbjerg (1993) H. Flyvbjerg, Phys. Rev. E 47, 4037 (1993).
- Monnereau and Vignes-Adler (1998) C. Monnereau and M. Vignes-Adler, Phys. Rev. Lett. 80, 5228 (1998).
- Weaire and Hutzler (2000) D. Weaire and R. Hutzler, The physics of foams (Oxford University Press, Oxford, 2000).
- Kim et al. (2006) S. G. Kim, D. I. Kim, W. T. Kim, and Y. B. Park, Phys. Rev. E 74, 061605 (2006).
- Arenzon et al. (2007) J. J. Arenzon, A. J. Bray, L. F. Cugliandolo, and A. Sicilia, Phys. Rev. Lett. 98, 145701 (2007).
- Avelino et al. (2008) P. P. Avelino, C. J. A. P. Martins, J. Menezes, R. Menezes, and J. C. R. E. Oliveira, Phys. Rev. D 78, 103508 (2008).
- Avelino et al. (2011) P. P. Avelino, R. Menezes, and J. C. R. E. Oliveira, Phys. Rev. E 83, 011602 (2011).
- Avelino et al. (2005) P. P. Avelino, J. C. R. E. Oliveira, and C. J. A. P. Martins, Phys. Lett. B 610, 1 (2005).
- vid (a) URL http://www.youtube.com/watch?v=Q3Y43kKHC_8.
- vid (b) URL http://www.youtube.com/watch?v=LeAz2i2FUu8.
- vid (c) URL http://www.youtube.com/watch?v=kZRWQwDNq70.
- vid (d) URL http://www.youtube.com/watch?v=__jSZOyMtxE.
- vid (e) URL http://www.youtube.com/watch?v=jFEmm7W3b2Q.