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.

Introduction

Coexistence of multiple species is observed in nature. The mechanism of coexistence has been discussed by many authors[1-4]. Plausible mechanisms for the coexistence have been presented, such as spatial and temporal segregations[5-8], cooperative interactions[4,9-12] and so on. In the present paper, we focus on the biodiversity in rock-paper-scissors (RPS) systems[13-21]. We apply metapopulation dynamics[22-24] to RPS systems, and report that the heterogeneity of a metapopulation network promotes the coexistence of species.

The cyclic balance in ecosystems maintains biodiversity. An example is a relationship among plant, herbivore and carnivore. Herbivores eat plants, but they are eaten by carnivores; when the carnivores die, they become nutrition for plants. Such a cyclic association is very common in real ecosystems[25-30]. More concrete examples of RPS games are the mating strategies of side-blotched lizards[31], marine sessile organisms[32,33], mutant strains of yeast[34], grass-tree systems[35], three strains of Escherichia coli[36,37] and fish in fresh water[38]. These species in cyclic relation can coexist in nature.

Much literature exists for theoretical works on the population dynamics in RPS systems. In 1973, Itoh has presented a “well-mixed” (“global interaction”) model for rock (R), scissors (S) and paper (P)[13]. He mathematically proved that the population dynamics are represented by classical Lotka-Volterra equation: the densities of three species (R, S, P) oscillate periodically (“neutrally stable”). However, when the total population size is finite, three species cannot coexist. In contrast, Tainaka has presented a lattice model (stochastic cellular automaton) for RPS games[5,6]. The collision occurs between adjacent sites (“local interaction”). He showed that three species can stably coexist. Such a difference between global and local interactions has been verified by the experiments with Escherichia coli[36]. In the case of local interaction, the population dynamics are largely affected by spatial pattern formation of species[14,20,39,40].

The spatial RPS game has been extended to network models[41,42]. In these cases, each node means an individual (agent), and a link means the interaction between agents. It is, however, very rare to discuss the relationship between network structure and biodiversity. Masuda and Konno have studied RPS games on complex networks, and discussed the relation between network structure and species coexistence[42]. In the present paper, we discuss the same relation, applying a metapopulation model[22].

The metapopulation model is popular in biology (ecology)[22-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-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 , 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 bywhere 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 . Therefore, the equilibrium densities (non-zero values) are given by

If (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 (). 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 be the density of species α in node i (α = R, S, P and ), and be the total density in node i. From the definition, we haveand . Provided that the migration rate for all individuals take the same value (unity), the density in node i is described bywhere k is the degree of node i and N 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

Figure 1

Schematic illustration of the metapopulation model (N = 3). Each graph is composed of three subpopulations (nodes) 1, 2, and 3. Each line connecting a pair of nodes denotes the path (link). Two configurations are possible for N = 3. (a) Homogeneous graph: all nodes have the same degree (number of links). (b) Heterogeneous graph.

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 (: 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 ). 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 2

Results of population dynamics at rate . (a) Metapopulation model for homogeneous graph in Fig. 1(a). The densities of three species (R, S, P) in node 1 are plotted against time t. The initial values are for species α in node i (α = R, S, P and ), but the rock density in node 3 is only set as . 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 and .

Figure 3(a) shows an orbit for Fig. 2(a), where the scissors density is plotted against rock density 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 . Figure 3(b) shows the plot of in node 2 against in node 1. When increases (decreases), then 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.

Figure 3

The orbits for homogeneous graph. (a) The scissors density is plotted against the rock density in node 1. (b) Plots of rock density in node 2 against the rock density in node 1.

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(a) shows the population dynamics for the heterogeneous graph. The densities of three species exhibit damping oscillations. The amplitudes are gradually decreased. The densities approach different equilibrium values depending on different nodes:where  = R, S, or P (see Supplementary Materials). From equation (5), the equilibrium density in each node is given by

Figure 4

Results of population dynamics for heterogeneous graph [Fig. 1(b)] at rate . (a) Three (R, S, P) densities in nodes 1 and 2 are depicted against time t. (b) Two orbits are displayed. (i) the black orbit: (, ) in node 1, and (ii) the red orbit: (, ) in node 2. Initial values are the same as used in Fig. 2(a).

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 (, ) in node 1, and ii) the red orbit of (, ) 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 byfor any node () (see Supplementary Materials). Namely, these values are just one-third of equation (3). In Fig. 5(b), two orbits (, ) and (, ) in node 1 are displayed. The final trajectories exhibit two circles where the upper and lower circles represent (, ) and (, ), respectively. The phase difference between rock (scissors) and scissors (paper) is . The R, S or P densities synchronize among different nodes. Since we set 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.

Figure 5

Population dynamics at rate and for homogeneous graph. (a) Three densities in node 1 are plotted against time t. (b) Two orbits (, ) and (, ) in node 1 are displayed.

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

Figure 6

Same as Fig. 5, but for heterogeneous graph. (a) Population dynamics. 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 (ρ (t), )) and (, ) in node 1, and red curves denote the orbits (, )) and (, ) in node 2.

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 (, )) and (, ) 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-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 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-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.