Heterogeneous network promotes species coexistence: metapopulation model for rock-paper-scissors game

Understanding mechanisms of biodiversity has been a central question in ecology. The coexistence of three species in rock-paper-scissors (RPS) systems are discussed by many authors; however, the relation between coexistence and network structure is rarely discussed. Here we present a metapopulation model for RPS game. The total population is assumed to consist of three subpopulations (nodes). Each individual migrates by random walk; the destination of migration is randomly determined. From reaction-migration equations, we obtain the population dynamics. It is found that the dynamic highly depends on network structures. When a network is homogeneous, the dynamics are neutrally stable: each node has a periodic solution, and the oscillations synchronize in all nodes. However, when a network is heterogeneous, the dynamics approach stable focus and all nodes reach equilibriums with different densities. Hence, the heterogeneity of the network promotes biodiversity.

The metapopulation model is popular in biology (ecology) [22][23][24] . The metapopulation consists of spatially separated habitats (patches or nodes); this is because the whole population of a species is usually separated into some nodes. Individuals can migrate between nodes. In most cases, the individuals move from higher-to lower-density nodes [43][44][45][46][47][48][49] . However, in the present model, we apply a "random migration": each agent randomly determines the destination of migration 50 . The RPS reactions only occur inside each node. By solving the reaction-migration equations analytically or numerically, we show that the RPS dynamics between homogeneous and heterogeneous graphs are significantly different. It is found that three species can stably coexist only on the heterogeneous graph. The heterogeneity can help to maintain the coexistence of species.

Models
In RPS game, each individual is either rock (R), scissors (S) or paper (P). Interactions (RPS games) take place inside each node as follows: where the parameters a, b, and c are victory rates 17,40,51 . If = = = a b c 1, then the system (3) becomes "standard" RPS system 5,6,13 . Three reactions represent a cyclic relation: species R beats S, S beats P, and P beats R.
First, we consider a case of no migration where the total population lives in a single habitat (node) 13,15,38 . When all individuals are assumed to be completely mixed (global interaction), the dynamics for system (1) are given by where ρ α t ( ) is the density of species α at time t (α = R, S, P). The total population size is assumed to be constant; thus, we put ρ ρ ρ In the equilibrium state, we set ρ = α d dt / 0. Therefore, the equilibrium densities (non-zero values) are given by If = = = a b c 1 (standard RPS game), then all densities take the same value (1/3). The equilibrium state is not stable 37 . The density of each species periodically oscillates around the equilibrium density. The time average of a species density over one period is given by equation (3).
Next, we consider metapopulation models which have N nodes ( ≥ N 2). Population dynamics are expected to be different whether the network is homogeneous or heterogeneous. When all nodes have the same degree (number of links), we call it homogeneous. We choose N = 3, because it is the simplest case to have both homogeneous and heterogeneous graphs. In Fig. 1, three circles mean habitats (nodes), and lines denote paths (links). The whole population consists of three nodes (N = 3). All individuals can move their nodes along a path. Let ρ αi be the density of species α in node i (α = R, S, P and = i 1, 2, 3), and ρ i be the total density in node i. From the definition, we have . Black, red and green curves indicate the densities of rock, scissors, and paper in node 1, respectively. (b) A single habitat case without migration. By the use of equation (2), the mean-field dynamics are displayed with initial values ρ ρ = = .
Provided that the migration rate for all individuals take the same value (unity), the density in node i is described by where k i is the degree of node i and N i represents the nearest neighbors of node i (sum over all paths) 50 . In equation (5), the terms on the right-hand side indicate the amount of incoming and outgoing individuals. According to mean-field theory, the densities of R, S, and P walkers in node i are described by In equations (6a-c), the first and second terms on the right-hand side represent the migration and reactions of RPS games, respectively.

Results
Case of homogeneous networks. Basic equation (6) is applied to two graphs in Fig. 1 (see Supplementary Materials). First, we report the results for the homogeneous graph [ Fig. 1(a)]. The victory rates are set to be equal ( = = = a b c 1: standard RPS game). Figure 2(a) shows the plots of R, S, P densities in node 1 against time t. Black, red and green curves indicate the densities of rock, scissors, and paper in node 1, respectively. The density of each species oscillates periodically around the equilibrium density (ρ = ). This value is the same for any species α and any node i (α = R, S, P and = i 1, 2, 3). For the sake of comparison, Fig. 2(b) displays the population dynamics obtained by equation (2) (no migration case). By the comparison between Fig. 2(a) and (b), we find the frequencies of oscillations in Fig. 2(b) are three times as large as those in Fig. 2(a). In the metapopulation model, the oscillation in each node becomes slower due to the random walk. Figure 3(a) shows an orbit for Fig. 2(a), where the scissors density ρ t ( ) in node 1. The dynamics finally approach a circle. The orbits rotate clockwise on the circle. The phase difference between rock and scissors densities is just π 3 /4. Figure 3 increases (decreases) simultaneously. Hence, the oscillations of rock densities are completely synchronized in both nodes 1 and 2. We can confirm that the densities of each species (R, S, or P) simultaneously oscillate in different nodes.
Case of heterogeneous networks. Next, we report the results for the heterogeneous graph [ Fig. 1(b)].
Individuals can migrate between a pair of nodes by random walks. However, they never migrate between nodes 2 and 3. Figure 4  where α = R, S, or P (see Supplementary Materials). From equation (5), the equilibrium density in each node is given by Namely, the density in node 1 is twice as large than that in node 2 (3). Individuals come together in node 1 (hub) by random walk.
In Fig. 4(b), two orbits are shown. i) the black orbit of (ρ t ( ) S,1 , ρ t ( ) R,1 ) in node 1, and ii) the red orbit of (ρ t ( ) S,2 , ρ t ( ) R,2 ) in node 2. Both orbits rotate clockwise and approach equilibrium points (stable focuses). Similarly, the densities in node 3 approach the same equilibrium as in node 2. The stable densities in node 1 and 2 (3) are given by equation (7). Thus, the dynamic behavior on the heterogeneous graph [see Fig. 3(a)] is significantly different from that on the homogeneous graph [see Fig. 4(b)].
Effect of the victory rates. In general, the victory rates take arbitrary values. When three rates a, b and c are all changed, the dynamics become very complicated. According to the previous works 17,38,51 we fix b = c = 1 and change the value of a which is the victory rate of rock. First, we report the results for the homogeneous graph [ Fig. 1(a)]. Figure 5 displays the effect of the victory rate, where (a) population dynamics and (b) orbits. The densities oscillate periodically around the equilibrium point expressed by R i e Si e P i e , , , , , , for any node ( = i 1, 2, 3) (see Supplementary Materials). Namely, these values are just one-third of equation (3). In Fig. 5(b), two orbits (ρ t ( ) R,1 , ρ t ( ) S,1 ) and (ρ t ( ) S,1 , ρ t ( ) P ,1 ) in node 1 are displayed. The final trajectories exhibit two circles where the upper and lower circles represent (ρ t ( ) R,1 , ρ t ( ) S,1 ) and (ρ t ( ) S,1 , ρ t ( ) P ,1 ), respectively. The phase difference between rock (scissors) and scissors (paper) is π 3 /4. The R, S or P densities synchronize among different nodes. Since we set = .
a 0 8, the average densities of all species over one period are not equal; the densities of paper are lower than those of rocks or scissors in all nodes. Next, we report the results for the heterogeneous graph [ Fig. 1(b)]. Figure 6 displays (a) population dynamics and (b) orbits. The upper three curves in Fig. 6(a) show the time dependence of three densities in node 1. The lower three curves represent the densities in node 2 (3). In node 1, three densities approach On the other hand, in node 2 (3), they approach  Figure 6(b) shows the four orbits which exhibit stable focuses. Here, the black curves are orbits of (ρ t ( ) R,1 , ρ t ( ) S,1 )) and (ρ t ( ) S,1 , ρ t ( ) P ,1 ) in node 1, and red curves denote the same orbits but in node 2 (3). Thus, the equilibrium densities expressed by equation (10) are stable on the heterogeneous graph.

Conclusion
We have developed the metapopulation model for RPS games with three subpopulations (nodes). In metapopulation models, the individuals usually migrate from higher-to lower-density nodes 25,26,[40][41][42][43][44][45][46] . However, in the present paper, we apply a random migration: each individual randomly determines the destination of migration 50 . The RPS reactions only occur inside each node. By solving the reaction-migration equations analytically or numerically, we show that the RPS dynamics between homogeneous [ Fig. 1(a)] and heterogeneous [ Fig. 1(b)] graphs are significantly different.
We first show the population dynamics for a two-node graph (N = 2). In this case, the dynamics are similar to the single-node case (N = 1): the system never approaches an equilibrium state (periodic oscillation). The two-node graph is homogeneous, because both nodes have the same degree (one link). Next, we choose N = 3, because it is the simplest case to have both homogeneous and heterogeneous graphs (see Fig. 1). In the case of the homogeneous graph [ Fig. 1(a)], the dynamics are similar to the single-node case (see Figs 2,3 and 5). The amplitudes of oscillations depend on the initial condition, and the frequency becomes one-third compared to the Upper three curves denote the three densities in node 1, while lower three curves mean those in node 2. (b) Four orbits. Black curves denote orbits of (ρ R,1 (t), ρ t ( ) S,1 )) and (ρ t ( ) S,1 , ρ t ( ) P ,1 ) in node 1, and red curves denote the orbits (ρ t ( ) R,2 , ρ t ( ) S,2 )) and (ρ t ( ) S,2 , ρ t ( ) P ,2 ) in node 2.
Scientific REPORTS | (2018) 8:7094 | DOI:10.1038/s41598-018-25353-4 single-node case (see Fig. 2). However, in the case of the heterogeneous graph [ Fig. 1(b)], the dynamics are represented by the stable focus (see Figs 4 and 6). Three species (R, P and S) stably coexist. Hence, we conclude that the heterogeneity of networks promotes stable coexistence in RPS systems. We can confirm this conclusion even for the star graph of N = 4; the star graph is also heterogeneous, because it has a hub at the center of the graph.
In nature, the coexistence of species with cyclic association is widely observed [31][32][33][34][35][36][37] . The well-mixed (global interaction) model never explains such a stable coexistence [see equation (2)]. To explain the biodiversity in RPS systems, plausible mechanisms for coexistence have been presented. A typical example is a local interaction: spatial distributions of species promote the coexistence of species 5,6,15,36,52 . In the present paper, we present another mechanism using the metapopulation model: the heterogeneity of the network promotes coexistence. In real ecosystems, the heterogeneous graph is more popular than the homogeneous one, since natural links are often disturbed by various obstructions.