Effects of cyclic allele dominance rules and spatial structure on the dynamics of cyclic competition models.

Barreto et al. (2017) showed that the genotypic cyclic competition model with three phenotypes appearing as three morphs of male lizards' throats had the same equilibrium but a wider stability region as the corresponding phenotypic model. In this paper we re-examine stability conditions under the symmetric choice of parameters for the phenotypic model so we can set the same internal equilibrium densities for all three phenotypes. In this setting we compare the stability regions of cyclic allele dominance rule. Next we consider the dynamics on a two-dimensional square lattice space and then show the effect of this spatial structure on the stability of phenotypic model. We obtain the following results: (i) Cyclic allele dominance rule in a genotypic model gives a wider stable region of internal equilibrium than the allele dominance rule observed in lizards; and (ii) spatial structure drastically changes dynamical behavior, especially when all three phenotypes coexist in almost all the parameter spaces when both competition and dispersal occur locally.


Introduction
One of interactions between three types of populations is known as non-transitive, cyclic competition or a rock-paper-scissors relationship, that is, type A is stronger than type B, and type B is superior to type C, which in turn is better than A. We can find several examples of this kind of interaction in nature; the most well-known are bacteria (e.g., [1][2][3][4]) and lizards (e.g., [5][6][7]).
Barreto et al. [8] proposed phenotypic and genotypic cyclic competition models to clarify the fact that three throat morphs in male lizards, known as orange, blue, and yellow, are maintained with a rock-paper-scissors relationship. Both models have the same equilibrium, but the latter gives a wider  stability region than the former. Barreto et al. [8] analyzed their model by fixing all payoffs except two elements, which are given by variable parameters. However, such a case is difficult to manage because an internal equilibrium depends on these parameters. Then we consider an alternative setting as an internal equilibrium with equal frequencies for all phenotypes independent of parameters, namely, we use symmetric conditions of parameters to obtain global insights into the dynamics in this simplified system for the convenience of analysis.
In this paper, we re-examine stability conditions under the following situations: symmetric parameters for the phenotypic model, cyclic allele dominance rule, and spatial structure. We obtain the following results: (i) Cyclic allele dominance rule in a genotypic model gives a wider stable region of internal equilibrium than the allele dominance rule observed in lizards. (ii) Spatial structure drastically changes the dynamical behavior, especially when all three phenotypes coexist in almost all the parameter spaces when both competition and dispersal occur locally.
In the next section, we review the phenotypic model and its corresponding genotypic model proposed by [8] and then restrict these models with two parameters to generate an internal equilibrium independent of the parameters. In section 3, we derive local stability conditions for each model. We consider cyclic allele dominance rule in section 4. In section 5, we investigate the effects of spatial structure on the stability of the internal equilibrium.

The phenotypic model and its corresponding genotypic model
First, we consider a phenotypic model of cyclic competition, or a rock-paper-scissors game, following [8]. There are three phenotypes O, B, and Y. O has a better strategy than B, and B is superior to Y, which is, in turn, surpasses O. When we use a payoff matrix in Table 1 Using these payoffs obtained by the competition between two phenotypes, the dynamics of phenotypic model are defined as This is a standard model of the rock-paper-scissors game (e.g., [9,10] (2.14) We have another three equilibria (1, 0, 0), (0, 1, 0), and (0, 0, 1) as the boundary equilibria (p O ,p B ,p Y ).
Next we move to a genotypic model corresponding to the above phenotypic model by the genotypic allele dominance rule in which an allele o dominates over two others, y and b, and an allele y dominates b. In other words, genotypes oo, oy, and ob are phenotype O, genotypes yy and yb are phenotype Y, and genotype bb is phenotype B ( Table 2).
Then phenotypic frequencies can be expressed by the genotypic frequencies g i j,n (i, j ∈ {o, y, b}) at generation n as follows: p O,n = g oo,n + g oy,n + g ob,n , (2.15) p Y,n = g yy,n + g yb,n , Third, assuming random mating gives the genotypic frequencies at generation n + 1, g i j,n+1 (i, j ∈ {o, y, b}) as where W O = M OO (g oo,n + g oy,n + g ob,n ) + M OB g bb,n + M OY (g yy,n + g yb,n ), (2.33) W B = M BO (g oo,n + g oy,n + g ob,n ) + M BB g bb,n + M BY (g yy,n + g yb,n ), (2.34) Using the above equations, we can obtain the following relations: Therefore, because we also have the same relationships (2.8) for genetypic model, we have the same phenotypic equilibrium (2.9)-(2.11) with eqs.(2.12)-(2.14) as the corresponding phenotypic model.

Stability condition for equilibrium (1/3, 1/3, 1/3)
The payoff matrix includes nine payoffs and is a little complicated. Following [8], we reduce it to two values. We fix the values of the payoffs at 1 for the competitions between the same phenotypes, but we also adopt the symmetric cases -namely, we give a larger payoff α > 1 to a stronger competitor Table 3. Payoff matrix with only two parameters. We obtain the same equilibrium fractions of 1/3 for all three phenotypes. Parameter ranges are given by 0 < β < 1 < α.

Opponent Individual
and a smaller one β < 1 to a weaker competitor (Table 3), which results in the same equilibrium fractions of 1/3 for all three phenotypes [9]. By using this reduced payoff matrix, eqs.(2.1)-(2.7) become By eqs.(2.12)-(2.14), so that an internal equilibrium (2.9)-(2.11) becomes The linearized system of eqs.(3.1)-(3.3) about an internal equilibrium (3.4) gives the following Jacobian matrix: Then the characteristic equation of the linearized system becomes Jury conditions or Schur-Cohn criteria (e.g., [11]) for second-degree characteristic equations λ 2 + a 1 λ + a 2 = 0 are known as Because a 1 < 0 from eq.(3.5), the left inequality of eq.(3.7) becomes −a 1 < 1+a 2 , that is, a 1 +a 2 +1 > 0, which holds by the calculation using eqs.(3.5)-(3.6): The right inequality of eq.(3.7) 1 + a 2 < 2, that is, Notice that it is symmetric with respect to α and β; namely, eq.(3.8) does not change when α and β are exchanged. Indeed, it is known that eq.(3.8) becomes a globally asymptotic stable condition, which can be checked using the Lyapunov function [9]. If eq.(3.8) does not hold, then we can observe the trajectory of the heteroclinic cycle between three vertices in the triangular space for three phenotypic frequencies [9] similar to the three-species Lotka-Volterra cyclic competition model [12] (See Figure 1  On the other hand, for a genotypic model, we get the following unique internal equilibrium by eqs.(2.27)-(2.36): 10) (3.14) Of course, We can say that an internal equilibrium does not depend on the parameters α nor β on either the phenotyic or genotypic models. The characteristic equation of the linearized system of eqs.(2.27)-(2.32) with eqs.(2.33)-(2.36) about an internal equilibrium (3.9)-(3.14) becomes Noticing that a 1 < 0 from eq.(3.15), the left inequality of eq.(3.7) holds because, using eqs.(3.15)-(3.16), The right inequality of eq.(3.7) a 2 − 1 < 0 becomes √ 2 which is a quadratic inequation on α or β and is also symmetric with respect to α and β. Figure 2(a) shows the local stability region obtained by eq.(3.8) and eq.(3.17). We can observe their dynamics by numerical simulations (Figure 1). . α and β indicate payoffs for strong and weak competitors, respectively. Green, red and yellow regions indicate stable for both phenotypic and genotypic models, stable for genotypic model but unstable for phenotypic one and unstable for both phenotypic and genotypic models, respectively. Allele dominance rule (II) has larger stability region than (I).

Cyclic allele dominance rule expands the stability region of the internal equilibrium
Common side-blotched lizards have an allele dominance rule (I), that is, an allele o for an orange throat is the most dominant over others, an allele y for a yellow one is intermediate, and an allele b for a blue one is the most recessive (Table 2). In other words, this allele dominance rule has three kinds of alleles: most dominant, intermediate, and most recessive. We refer to these as o, y, and b, respectively.
Although perhaps we have not yet discovered it in realistic genetic systems, we can theoretically consider another allele dominance rule (II), that is, all alleles having intermediate dominance: o dominant over y, y over b, and b over o (Table 4). This rule could be called as a "cyclic allele dominance rule." Similarly, in this rule, we name these alleles o, y, and b, and any allele can be replaced by another one.
A calculation for allele dominance rule (II) similar to that for allele dominance rule (I) reveals that the boundary of local stability obtained from the Jury condition (3.7) is determined by which only differs from the first calculation in the smaller coefficient of αβ − 1 compared to eq.(3.17). Figure 2(b) shows the local stability region by allele dominance rule (II) using eq.(4.1), which clearly has a larger stability region than rule (I).
From a theoretical point of view, we should also consider the cases with the same allele dominance rules but with different cyclic competitive strengths. Then another possible combination exists, that is, Table 2 with α < 1 and β > 1. However, this combination gives the same stability regions as Figure  2(a) by the symmetry of the boundary equation (3.17) on the parameters α and β; eq.(3.17) does not change when α and β are exchanged. On the other hand, it is trivial that no qualitative difference exists between the two cases: Table 4 with α > 1 and β < 1 and Table 4 with α < 1 and β > 1.

Spatial simulation
Barreto et al. [8] showed the same stabilizing effect in a genetic system of lizards as the previous section of this paper. Here we introduce a spatial structure into a phenotypic model and consider its effects on population dynamics.
The phenotypic model includes two processes: competition and reproduction. When we add a spatial structure into this model, we restrict the spatial range, both for the opponents against a focal individual and for the dispersion of offspring by reproduction. When we set such a spatial restriction as a whole or as a neighborhood, we consider four distinct cases:  (Table 3). All the payoffs are transferred to their offspring, and they are gathered in one offspring pool in the cases of (a) and (c) or, in the cases of (b) and (d), in local offspring pools. Here, we use lower cases for the sites and corresponding upper cases for the states of the individuals on those sites.
(iii) The second process is dispersion, which gives positions of offspring selected from the offspring pool. The model is discrete-generation: all individuals die, and the states at the next time steps for all the sites i (= 1, . . . , N) in a whole population with size N are replaced in order. For each site, an individual is randomly chosen from a whole pool for (a) and (c) or from each local pool for (b) and (d); this choice is made in proportion to the relative payoff against total payoffs in a whole population in the cases of (a) and (c) or in the nearest neighboring four individuals in the cases of (b) and (d). In other words, this random selection depends on W O W , W Y W , and W B W . After the previous process (ii), all the individuals have their own payoffs. In the cases of (a) and (c), we can then calculate W, W O , W Y , and W B as the total payoffs for, respectively, a whole population, phenotype O, phenotype Y, and phenotype B. In the cases of (b) and (d), we can obtain W, W O , W Y , and W B as the total payoffs for, respectively, a local population, phenotype O in the local population, phenotype Y in the local population, and phenotype B in the local population. Here a local population is restricted to individuals on five sites: the focal site and its nearest neighboring four sites.
We show the results of the Monte Carlo simulation by the above algorithm in Figure 3. Case (a) completely coincides with Figure 2(a). Case (b) gives a wider coexistence region due to local dispersion. Case (c) has the same result as case (a) because the global dispersal produces the same effect with random choices of opponents from the population. Case (d) shows coexistence in almost all the simulations for the entire parameter space. Therefore, we can conclude that locally limited interactions between individuals strongly promote the coexistence of all phenotypes.
In case (d), individuals with the same phenotypes tend to cluster and adopt collective behavior, so changes of state occur only on the boundaries between these clusters, causing a slower change of total population and stabilization (See Figure 4). Indeed [13] already reached such a conclusion from a continuous-time cyclic competition model that also calculated the average domain size or boundary length at a steady state. This spatial structure in the population is gradually produced by both local competition and local dispersion, and, in turn, also gives different results in competition and dispersion, either globally or locally.

Discussion
Similar to [8], we show the difference between the phenotypic model and its corresponding genotypic model, but we can also give an explicit condition using a simplified system. In addition, we investigate other allele dominance rules to clarify their effects on the stability of internal equilibrium compared to the allele dominance rule observed in lizards. Unfortunately, however, we cannot specify the mechanism by which the genotypic model stabilizes the system more than the phenotypic model. To discover it, we need to find models other than cyclic competition in which the genotypic models give different dynamic behaviors than the corresponding phenotypic models.
A recent proposal by [14] may have important implications for future work in this area. As they point out, the present three-strategy payoff matrix can be built up as the sum of nine independent Fourier components, g(1), . . . , g (9), with orthogonality and normalization: where we have g(1), g (8), and g(9) as irrelevant constant terms, a coordination (Potts) component that equally favors the formation of one of the homogeneous states, and a cyclic (rock-paper-scissors) component, respectively. This formulation may offer more detailed insights into dynamical behaviors if we examine cases with various combinations of coefficients corresponding to the strengths of these components. The intuitive reason local competition and local dispersion both stabilize an internal equilibrium is that several dozens of identical phenotypic individuals gather and adopt collective behavior, restricting phenotypic changes on the boundary. To investigate their stabilization mechanisms in detail, we should further rely on other kinds of analysis, especially spatial pattern formation, which may play an important role in stabilizing the dynamics. So far, the concept of "vortex" has been proposed to characterize the spatial pattern [15][16][17], and it has also been shown that clockwise or anti-clockwise rotating spiral patterns with characteristic sizes emerge that yield finite size effects by cyclic competition models or rock-paper-scissors games in physics [18][19][20][21] and evolutionary games [22,23]. More precisely, the coexistence of three strategies transforms into one of the homogeneous states after a suitable relaxation time if the characteristic length of patterns exceeds a threshold value comparable to the system size.
Here we use a small lattice space size (100 × 100) for Figure 3, but we should further investigate these finite size effects using a larger lattice space, such as 1000 × 1000, as in Figure 4.
In this paper, we only show whether three phenotypes can coexist or not in the long run using a Monte Carlo simulation. However, we can expect an internal equilibrium to be asymptotically stable for the case of coexistence in larger Monte Carlo systems than in smaller ones ( Figure 5). Here we can only show the condition of the locally asymptotic stability of a genotypic model. We expect that condition to be replaced by global stability, but we leave it as a future problem.
In addition, we would like to examine a genotypic model in a lattice space and clarify whether it becomes more stable than a phenotypic model and to what degree its stability can be evaluated by, for example, the return time of perturbation to an internal equilibrium.
Throughout this paper, we use only symmetric conditions of parameters to obtain global insights into the dynamics, then we should study more complicated models with symmetry-breaking interactions, including the parameter setting of [8], in future research. I sincerely appreciate Prof. Akira Sasaki for giving me useful comments and encouraging me. I would also like to thank three anonymous reviewers whose comments were helpful in revising the manuscript.