Effect of Spatial Dispersion on Evolutionary Stability: A Two-Phenotype and Two-Patch Model

In this paper, we investigate a simple two-phenotype and two-patch model that incorporates both spatial dispersion and density effects in the evolutionary game dynamics. The migration rates from one patch to another are considered to be patch-dependent but independent of individual’s phenotype. Our main goal is to reveal the dynamical properties of the evolutionary game in a heterogeneous patchy environment. By analyzing the equilibria and their stabilities, we find that the dynamical behavior of the evolutionary game dynamics could be very complicated. Numerical analysis shows that the simple model can have twelve equilibria where four of them are stable. This implies that spatial dispersion can significantly complicate the evolutionary game, and the evolutionary outcome in a patchy environment should depend sensitively on the initial state of the patches.


Introduction
In order to explain the evolution of animal behavior, Maynard Smith and Price [1] developed the concept of evolutionarily stable strategy (ESS) (see also [2][3][4][5]). Prior et al. [6] investigated an evolutionary game model that incorporates both spatial dispersion and density effects in the evolutionary dynamics. In this model, the population is considered to be dispersed in a patchy environment, where the background fitness and payoff matrix in each patch can be different. Migration from region to region is considered as an incidental aspect of the population, i.e., the migration is a chance event unrelated to an individual's phenotype (strategy) or the fitness of the patch. As pointed out by Prior et al. [6], their assumptions differ from that of Ludwig and Levin [7] who treat the tendency to migrate as an individual characteristic subject to selection (see also [8][9][10][11][12][13]), and also differ from that of Hines and Maynard Smith [14] who interpret the effect of spatial dispersion as an increased tendency to interact with opponents sharing one's own characteristics (see also [15]). Recently, Cressman and Krivan [16] investigated the migration dynamics for the ideal free distribution (IFD) in a patchy environment. They showed that IFD is evolutionarily stable under the assumptions that individuals never migrate from patches with a higher payoff to patches with a lower payoff and some individuals always migrate to the best patch. But migration does not necessarily lead to IFD if migration rates are independent of the payoffs of the patches.
For the evolutionary game dynamics in a patchy environment, Prior et al. [6] mainly focused their analysis on the stability of the homogeneous states, where they assumed that all patches have the same payoff matrix and density-dependent background fitness. Their main results showed that a stable equilibrium (e.g. an evolutionarily stable strategy) of the non-dispersed frequency dynamics becomes a stable equilibrium of the large system if population density stabilizes at these fixed frequencies.
In this paper, following Prior et al. [6], a simple two-patch and two-phenotype model is investigated. Three basic assumptions for this model are: (i) The environment consists of two patches, called patch 1 and patch 2, respectively. Individuals can move from one patch to the other at any time t. The migration rates are patchdependent but independent of individual's phenotype [6]. Let c 1 denote the probability that an individual moves from patch 1 to patch 2, and, similarly, c 2 the probability that an individual moves from patch 2 to patch 1.
(ii) In each of two patches, individuals display two possible phenotypes (strategies), denoted by R 1 and R 2 , and individuals interact in random pairwise contests. The payoff matrix is off of an individual displaying phenotype R i when it plays against an individual displaying phenotype R j in patch 1 (or in patch 2) for all i, j = 1, 2. Without loss of generality, we also assume that a ij ! 0 and b ij ! 0 for all i, j = 1, 2.
(iii) In each of the two patches, the background fitness is density-dependent [3][4], which is defined as α 1 − β 1 n 1 in patch 1, and α 2 − β 2 n 2 in patch 2, where n 1 is the total population size in patch 1 and n 2 the total population size in patch 2. We also assume that α i > c i for all i = 1, 2. That is, the migration rates are small enough to ensure that population size in a patch will increase when there are few individuals in the patch.
Let x i denote the number of individuals with phenotype R i in patch 1, and y i the number of individuals with phenotype R i in patch 2 (i = 1, 2). Clearly, n 1 = x 1 + x 2 and n 2 = y 1 + y 2 . Similarly, let p denote the frequency of phenotype R 1 in patch 1, and q the frequency of phenotype R 1 in patch 2, i.e., p = x 1 /n 1 and q = y 1 /n 2 . According to the basic assumption (ii), the expected payoff of an individual displaying phenotype R i is f i = pa i1 + (1 − p)a i2 in patch 1, and g i = qb i1 + (1 − q)b i2 in patch 2 for i = 1, 2. Similarly, according to the basic assumption (iii), the (total) fitness of an individual displaying phenotype R i is defined as F i = (α 1 − β 1 n 1 ) + f i in patch 1, and G i = (α 2 − β 2 n 2 ) + g i in patch 2 for i = 1, 2. Thus, the time evolution of x i and y i can be given by respectively, for i = 1, 2. Dynamics (1) also equivalent to the following system expressed in terms of phenotypic frequency and population size in each patch. dp dt ¼ pð1 À pÞðf 1 À f 2 Þ þ c 2 ðq À pÞ n 2 n 1 ; where f ¼ pf 1 þ ð1 À pÞf 2 and g ¼ qg 1 þ ð1 À qÞg 2 are the average payoffs in patch 1 and patch 2, respectively.
In this paper, the equilibria of dynamics (2) and their stabilities are analyzed. Different from Prior et al. [6] who focused on the homogeneous states, we are primarily interested in analyzing the heterogeneous states, where two patches have different ESSs. Our main goal is to reveal the dynamical properties of the evolutionary game in a heterogeneous patchy environment.

Symmetric equilibria of dynamics (2)
For given p and q with 0 p, q 1, an equilibrium of dynamics denoted by (n 1 (p, q), n 2 (p, q)), satisfies It is clear that (n 1 , n 2 ) = (0, 0) is always a solution of Eq (4) for any given p and q. Furthermore, (p, q, 0, 0) must be unstable under dynamics (2) since α 1 > c 1 and α 2 > c 2 . Notice that n 2 is a parabolic function of n 1 and vice versa, Eq (4) also has a unique positive solution, denoted by ðn 1 ;n 2 Þ ¼ ðn 1 ðp; qÞ; n 2 ðp; qÞÞ withn 1 > 0 andn 2 > 0 (see Fig 1). When this solution corresponds to an equilibrium ðp; q;n 1 ;n 2 Þ of dynamics (2), we call it a positive equilibrium. In the rest of this paper, we only focus on the number and stabilities of these positive equilibria. A positive equilibrium, ðp; q;n 1 ;n 2 Þ, is called a symmetric equilibrium of dynamics (2) if p = q. That is, at a symmetric equilibrium, population compositions in the two patches are the same.
It is easy to see that dynamics (2) always have two symmetric boundary equilibria, ð1; 1;n 1 ;n 2 Þ and ð0; 0;n 1 ;n 2 Þ, where at these equilibria, all individuals in the two patches display the same phenotype. Stabilities of the boundary equilibria can be characterized by analyzing the Jacobian matrix of dynamics (2) (see Method section). The main result is that if R 1 (or R 2 ) is an ESS for both payoff matrices A and B, then the boundary equilibrium ð1; 1;n 1 ;n 2 Þ (or ð0; 0;n 1 ;n 2 Þ) must be asymptotically stable. This result is consistent with that of Prior et al. [6], where the stable equilibrium of the non-dispersed system becomes a stable equilibrium of dynamics (2). However, under the influence of migration, a symmetric boundary equilibria could be stable even if the corresponding phenotype is not an ESS in either patch. For instance, the boundary equilibrium ð1; 1;n 1 ;n 2 Þ is asymptotically stable if a 11 < a 21 as long as a 21 − a 11 is small enough.

General cases
For more general situations (i.e., p Ã 6 ¼ q Ã ), it is very tedious to determine the equilibria of dynamics (2). In fact, numerical simulations show that dynamics (2) may have twelve equilibria. To investigate the properties of the equilibria of dynamics (2), two cases are considered below. The first case is special in that there is no migration in one direction (i.e., one of c 1 and c 2 is 0) and we analyze the number of stable equilibria for all possible payoff structures. The second case is more general (i.e., c 1 > 0 and c 2 > 0) and we show the equilibria of dynamics (2) and their stabilities for 0 < p Ã , q Ã < 1.
We first look at the stability of dynamics (5). It is easy to see that: (i) the boundary equilibrium ð1;n 1 Þ (or ð0;n 1 Þ) is locally asymptotically stable if and only if p = 1 (or p = 0) is an ESS for the payoff matrix A, i.e. a 11 > a 12 (or a 22 > a 12 ), From dynamics (6), the frequencyq in a stable equilibrium of dynamics (2), ðp;q;n 1 ;n 2 Þ, should obey the equation wherep 2 f0; 1; p Ã g corresponds to the stable equilibrium, ðp;n 1 Þ, of dynamics (5). In the Method section, we analyze the solutions of Eq (7) and the stabilities of the corresponding equilibria under dynamics (6). According to the stability conditions of dynamics (5) and (6), the equilibria of dynamics (2) and their properties can be summarized as follows: 1. If R 1 (or R 2 ) is the only ESS for A and B, then the symmetric boundary equilibrium ð0; 0;n 1 ;n 2 Þ (or ð1; 1;n 1 ;n 2 Þ) is unstable and the other ð1; 1;n 1 ;n 2 Þ (or ð0; 0;n 1 ;n 2 Þ) is stable. Furthermore, Eq (7) has no interior solution. This implies that ð1; 1;n 1 ;n 2 Þ (or ð0; 0;n 1 ;n 2 Þ) is also globally asymptotically stable, i.e., all individuals in the two patches will eventually display R 1 (or R 2 ) under evolutionary dynamics (2) (see Fig 2A and 2B).
3. If p Ã is the only ESS for A (0 < p Ã < 1) and R 1 (or R 2 ) is the only ESS for B, then both the symmetric boundary equilibria ð0; 0;n 1 ;n 2 Þ and ð1; 1;n 1 ;n 2 Þ are unstable. Furthermore, Eq (7) has a unique (interior) solutionq, which corresponds to an asymptotically stable equilibrium ðp Ã ;q;n 1 ;n 2 Þ of dynamics (2). Numerical simulation shows that this equilibrium is also globally stable. This implies that the two phenotypes will stably coexist in the system (see Fig 2E and 2F).
5. If both R 1 and R 2 are ESSs for A and B, then both the symmetric boundary equilibria ð0; 0;n 1 ;n 2 Þ and ð1; 1;n 1 ;n 2 Þ are asymptotically stable. Furthermore, Eq (7) has at most four (interior) solutions, where two are stable equilibria of dynamics (2) and the other two are unstable. This implies that the evolutionary outcome in this situation is very difficult to predict since the system can have four stable states (see Fig 2I).
6. If both R 1 and R 2 are ESSs for A and q Ã is an ESS for B, then both the symmetric boundary equilibria ð0; 0;n 1 ;n 2 Þ and ð1; 1;n 1 ;n 2 Þ are unstable. Furthermore, Eq (7) has at most two (interior) solutions, where both of them are stable equilibria of dynamics (2) (see Fig 2J). 7. If p Ã is an ESS for A and q Ã is an ESS for B, then both the symmetric boundary equilibria ð0; 0;n 1 ;n 2 Þ and ð1; 1;n 1 ;n 2 Þ are unstable. Furthermore, Eq (7) has a unique (interior) solutionq, which corresponds to an asymptotically stable equilibrium ðp Ã ;q;n 1 ;n 2 Þ of dynamics (2). Similarly as (3), this equilibrium is also globally stable and two phenotypes will stably coexist in the system (see Fig 2K).
10. If R 2 (or R 1 ) is the only ESS for A and both R 1 and R 2 are ESSs for B, then the symmetric boundary equilibrium ð0; 0;n 1 ;n 2 Þ (or ð1; 1;n 1 ;n 2 Þ) is stable and the other ð1; 1;n 1 ;n 2 Þ (or ð0; 0;n 1 ;n 2 Þ) is unstable. Furthermore, Eq (7) has at most two (interior) solutions, where one corresponds to a stable equilibria of dynamics (2) and the other is unstable. Similarly as (2), either all individuals in the system display R 2 (or R 1 ), or individuals in patch 1 display R 2 (or R 1 ) and two phenotypes coexist in patch 2 (see Fig 2O and 2P).
Case 2. c 1 > 0 and c 2 > 0. We now consider the case with c 1 > 0 and c 2 > 0. It is easy to check that dynamics (2) only have two boundary equilibria, ð0; 0;n 1 ;n 2 Þ and ð1; 1;n 1 ;n 2 Þ, and the existence of asymmetric boundary equilibrium is impossible, for instance, if p = 0 and q 6 ¼ 0, then dp dt > 0. We then focus on the number and stability of interior equilibria. Notice that an equilibrium of dynamics (2) should be the solution of equation So if both n 1 and n 2 are positive, then from the first two equations of Eq (8), an interior equilibrium of dynamics (2) should obey the equations Furthermore, from the third and the forth equations of Eq (8) where we assume that both p Ã and q Ã are in the interval 0 < p Ã , q Ã < 1 (i.e., we consider the most complicated payoff structures). From Eq (9), it is easy to see that for the situation with p Ã 6 ¼ q Ã , if the interior equilibrium exists, then it should be in the region (0, p Ã ) × (q Ã , 1), or (p Ã , 1) × (0, q Ã ) if Δ 1 Δ 2 > 0, and in the region (0, p Ã ) × (0, q Ã ), or (p Ã , 1) × (q Ã , 1) if Δ 1 Δ 2 < 0. Of course, it is very difficult to get the exactly analytic solutions of Eqs (9) and (10) in general. The numerical analysis suggests that ten interior equilibria can exist (see Fig 3H). To show this, some examples are plotted in Fig 3. All of these examples show clearly that the equilibrium structure of dynamics (2) could be very complicated.
We further look at the bifurcation behaviors of system (2) for the case that both R 1 and R 2 are ESSs for A and B (i.e., the most complicated case), and assume equal migration rates between regions, i.e., c 1 = c 2 = c. In this case, both the symmetric boundary equilibria ð0; 0;n 1 ;n 2 Þ and ð1; 1;n 1 ;n 2 Þ are asymptotically stable, and the system can have ten interior equilibria. When c = 0, it is easy to see that the system has nine equilibria in total, including four stable (boundary) equilibria and five unstable equilibria. The number of equilibria jumps from nine to twelve as soon as the migration rates becomes positive although the number of stable equilibria keeps unchanged (see Fig 4A and 4B). Furthermore, numerical simulation shows that both the numbers of stable equilibria and unstable equilibria decrease as c increases. In particular, when c > 0.019, the system has only two stable equilibria (i.e., the two symmetric boundary equilibria), and all interior equilibria are unstable (see Fig 4). These results suggest that small migration rates make the dynamical behavior of the system more complex.

Discussion
A vast amount of research has been devoted to analyze the influence of spatial diffusion on the evolutionary stability of ecology systems. One well known mathematical approach is the reaction-diffusion equation, where in this framework, individuals are dispersed in a continues space [10][11][17][18][19][20][21]. For instance, Hofbauer et al. [11,21] considered a population of two types of individuals distributed in an one-dimensional space, and assumed that the migration (or diffusion) rate is both individual-independent and location-independent. They showed that in two-strategy coordination games, if the reaction term of the reaction-diffusion equation is taken as replicator dynamics, then one strategy will drive out the other strategy in form of a traveling wave front, although there is no simple rule to decide which strategy can survive. In this paper, we assume that individuals are distributed in a (discrete) patchy environment. Following Prior et al. [6], we investigate a simple two-phenotype and two-patch model, where individuals compete only with their immediate neighbors and the migration rates between patches are individual-independent but patch-dependent. Different from Prior et al. [6] who focused on the homogeneous patchy environment, we here are more interested in the dynamical stability in heterogeneous environment. Our main results show that: (i) if the pure strategy R 1 (or R 2 ) is an ESS for both two patches, then the boundary equilibrium corresponding to p = 1 and q = 1 (or p = 0 and q = 0) must be asymptotically stable; (ii) if the payoff matrices A and B satisfy p Ã ¼ a 12 Àa 22 a 12 Àa 22 þa 21 Àa 11 ¼ q Ã ¼ b 12 Àb 22 b 12 Àb 22 þb 21 Àb 11 2 ð0; 1Þ, then the interior equilibrium corresponding to (p Ã , q Ã ) is asymptotically stable if p Ã is an ESS for A and q Ã an ESS for B; (iii) as a special case with c 1 > 0 and c 2 = 0 (or c 1 = 0 and c 2 > 0), i.e. individuals can only move from patch 1 to patch 2 (or from patch 2 to patch 1), all possible situations for the existence and stability of boundary and interior equilibria are considered, and we find that dynamics (2) can have six equilibria where four of them are stable; (iv) for c 1 > 0 and c 2 > 0, the numerical analysis shows that the equilibrium structure and dynamical behavior of the system could be very complicated in general. In particular, dynamics (2) can have twelve equilibria where four of them are stable.
Our analysis provides an insight for understanding the effect of spatial dispersion on the evolutionary stability of patchy environment. Both the analytical analysis and the numerical simulation indicate that the original ESS formulations which ignore the dispersion process cannot be applied to predict the evolutionary outcome of the dispersion system even for small migration rates. For instance, in the case that both patches have multiple ESS's and no dispersal between patches, the system has four (boundary) stable equilibria and five unstable equilibria. However, if one of c 1 and c 2 becomes positive, the system can have two to four stable equilibria and four to eight unstable equilibria. Furthermore, we found that both the numbers of stable equilibria and unstable equilibria decrease in the migration rates. This observation has an intuitive biological interpretation [6]. In a heterogenous patchy environment, the effect of selection is to make the overall population more heterogeneous in the sense of different patches have different population compositions, while the effect of migration is to move the population composition in each patch towards the mean of the overall population, i.e., migration promotes homogeneity. Thus, when the migration rates are small (i.e., the effect of selection is strong), similarly as the case of no dispersal, the system has two stable symmetric boundary equilibria and two stable asymmetric equilibria; and when the migration rates are large (i.e., the effect of migration is strong), the existence of stable asymmetric equilibrium is impossible, and the system has only two stable symmetric boundary equilibria, where at these equilibria all individuals display the same phenotype.
In this paper, we focus on the effect of spatial dispersion on two-patch system only. A natural extension would be to consider the three-patch system. However, analyzing the dynamical behavior of the three-patch system may be an even more difficult issue because the equilibrium structure of the two-patch system is already very complex. Another possible development would be to compare the evolutionary stability of the patchy environment under different migration rules. One commonly used migration rule is that individuals know perfectly the payoff in all patches and they always move to the patch with the highest payoff (i.e., ideal animals) [15]. In contrast, a more realistic model is that individuals do not migrate to patches with lower payoff [22]. Recent studies have shown that these two migration rules can lead to the IFD [16,23]. Since that the IFD corresponds to a stable equilibrium of the non-dispersed evolutionary dynamics, we can then expect that these migration rules may also lead to the ESS of the non-dispersed evolutionary dynamics [23].

Stability of the symmetric equilibria
The Jacobian matrix of the dynamics (2) about the symmetric boundary equilibrium ð1; 1;n 1 ;n 2 Þ, denoted by J (1,1) , is Àða 11 À a 21 Þ À c 2n and similarly, the Jacobian matrix about ð0; 0;n 1 ;n 2 Þ, denoted by J (0,0) , is A must be negative. So, if the pure strategy R 1 is an ESS for both payoff matrices A and B, then the eigenvalues of J (1,1) must have negative real parts [6]. Similar to the matrix J (0,0) , if the pure strategy R 2 is an ESS for both payoff matrices A and B, then the eigenvalues of J (0,0) have negative real parts. The Jacobian matrix about the symmetric interior equilibrium ðp Ã ; q Ã ;n 1 ;n 2 Þ, denoted by J (p Ã , q Ã ) , is p Ã ð1 À p Ã ÞD 1 À c 2n where Δ 1 = a 11 − a 12 − a 21 + a 22 and Δ 2 = b 11 − b 12 − b 21 + b 22 . Also similar to the matrix J (1,1) (or the matrix J (0,0) ), the eigenvalues of J (p Ã , q Ã ) have the negative real parts if p Ã (= q Ã ) is an ESS for both payoff matrices A and B, i.e., the equilibrium ðp Ã ; q Ã ;n 1 ;n 2 Þ is asymptotically stable if p Ã (= q Ã ) is an ESS for both A and B.
asymptotically stable under dynamics (6). On the hand, let ðq;n 2 Þ be an interior equilibrium of dynamics (6), and the Jacobian matrix about ðq;n 2 Þ, denoted by J ðq;n 2 Þ , is given by Thus, for given parameter values, stabilities of the interior equilibria of dynamics (6) (i.e., interior solutions of Eq (7)) can be analyzed numerically according to the above conditions (see the figure caption of Fig 2 for detailed parameters, note that the following results may not be true for all parameter values).