- Email: [email protected]
Patterns for four-allele population genetics model
Theoretical Population Biology 81 (2012) 273–283
Contents lists available at SciVerse ScienceDirect
Theoretical Population Biology journal homepage: www.elsevier.com/locate/tpb
Patterns for four-allele population genetics model Linlin Su, Roger Lui ∗ Department of Mathematical Sciences, Worcester Polytechnic Institute, 100 Institute Road, Worcester, MA 01609, United States
article
info
Article history: Received 1 September 2011 Available online 26 February 2012 Keywords: k-allele model Asymptotically stable equilibrium Pattern Eigenvalue Computer simulation
abstract In this paper, we find and classify all existing patterns for a single-locus four-allele population genetics models in continuous time. An existing pattern for a k-allele model means a set of all coexisting asymptotically stable equilibria with respect to the flow defined by the system of equations p˙ i = pi (ri − r ), i = 1, . . . , k, where pi and ri are the frequency and marginal fitness of allele Ai , respectively, and r is the mean fitness of the population. It is well known that for the two-allele model there are only three existing patterns, depending on the relative fitness between the homozygotes and the heterozygote. For the three-allele model there are 14 existing patterns, and we shall show in this paper that for the fourallele model there are 117 existing patterns. We also describe the domains of attraction for coexisting asymptotically stable equilibria. The problem of finding existing patterns has been studied in the past, and it is an important problem because the results can be used to predict the long-term genetic makeup of a population. It should be pointed out that this continuous-time model is only an approximation to the corresponding discrete-time model. However, the set of equilibria and their stability properties are the same for the two models. © 2012 Elsevier Inc. All rights reserved.
1. Introduction Consider a diploid population whose members possess a gene that occurs in k different forms, called alleles, located at an autosomal locus. Then there are k(k + 1)/2 genotypes denoted by Ai Aj , i, j = 1, . . . , k. Suppose that Ai Aj has constant fitness denoted by rij . Let R = (rij ) be the (symmetric) fitness matrix, and let pi be the frequency of allele Ai . Then, assuming random mating and discrete-time non-overlapping generations, the frequency of allele Ai in the next generation is given by p′i = pi ri /r ,
i = 1, . . . , k,
(1)
where ri = j rij pj is the marginal fitness of allele Ai and r = i ri pi is the mean fitness of the population. For continuous-time overlapping generations, the above equations are replaced with p˙ i = pi (ri − r ),
i = 1, . . . , k.
(2)
In this paper, we focus our study on continuous-time models because we can apply the theory of dynamical systems. Although system (2) is only an approximation to (1), derived under the assumption of Hardy–Weinberg equilibrium (Bürger, 2000, Section 10.1), the set of equilibria and their stability properties are the same for the two models.
∗
Corresponding author. E-mail addresses: [email protected] (L. Su), [email protected] (R. Lui).
0040-5809/$ – see front matter © 2012 Elsevier Inc. All rights reserved. doi:10.1016/j.tpb.2012.02.002
Since the pi represent gene frequencies, they add up to 1. In fact, it can be shown that, if initially i pi (0) = 1, then the same holds for the solutions pi (t ), i = 1, . . . , k, of system (2) for t > 0. Therefore, (2) is well posed, and is a system of k − 1 equations with pℓ , 1 ≤ ℓ ≤ k, replaced with 1 − j̸=ℓ pj . By relabeling the alleles we may assume, without loss of generality, that the rii are nonincreasing. To avoid degeneracy, we assume throughout this paper that r11 > r22 > · · · > rkk
(3)
and that rij , where i ̸= j, is not equal to rii or rjj . An equilibrium of system (2) is aconstant vector p = (p1 , . . . , pk ) that satisfies the conditions i pi = 1 and pi (ri − r ) = 0, i = 1, . . . , k. An equilibrium is called an interior equilibrium if pi > 0 for all i and a boundary equilibrium if pi = 0 for at least one i. Note that Vi = (0, . . . , 1, . . . , 0), where the 1 is in the ith component, is always a boundary equilibrium. Let us consider the two-allele case. Since p1 + p2 = 1, system (2) may be reduced to p˙ 1 = p1 (1 − p1 )(∆12 p1 + r12 − r22 ),
(4)
where ∆12 = r11 + r22 − 2r12 . The following result is well known (Nagylaki, 1992, pp. 52–55). Proposition 1.1. For the two-allele model (4) with r11 > r22 , p1 = 0 and p1 = 1 are both equilibria, and a third interior equilibrium p∗ = (r22 − r12 )/∆12 exists if and only if r12 > r11 or r12 < r22 . There are three cases: (i) r22 < r12 < r11 : 1 is globally asymptotically stable, 0 is unstable and p∗ does not exist; (ii) r12 > r11 : 0, 1 are both
274
L. Su, R. Lui / Theoretical Population Biology 81 (2012) 273–283
unstable and p∗ exists and is globally asymptotically stable; (iii) r12 < r22 : 0, 1 are both asymptotically stable and p∗ exists and is unstable. We call these three cases the (heterozygote) intermediate, superior, and inferior case, respectively. The above proposition implies that for the two-allele model there are only three existing patterns: ⟨{1}⟩, ⟨{p∗ }⟩, and ⟨{0}, {1}⟩. An existing pattern is simply the list of all coexisting asymptotically stable equilibria. For example, the pattern ⟨{p∗ }⟩ means that the interior equilibrium p∗ is the only asymptotically stable equilibrium (and, consequently, is globally asymptotically stable) and, from Proposition 1.1, this happens in the superior case. The pattern ⟨{0}, {1}⟩ means there are two coexisting asymptotically stable equilibria, p1 = 0 and p1 = 1, and this happens in the inferior case. On the other hand, the pattern ⟨{p∗ }, {1}⟩ does not exist regardless of the choice of R. For the k-allele model there are 2k − 1 possible equilibria, but the number of coexisting asymptotically stable equilibria is considerably less. For example, for the four-allele model, there are 15 possible equilibria and potentially 215 − 1 patterns, but only 117 of them can exist, and none can contain more than four asymptotically stable equilibria. The problem of finding all patterns for a k-allele model has been studied extensively before. In a series of papers, (Broom et al., 1993; Cannings and Vickers, 1988, 1989, 1991; Vickers and Cannings, 1988), the authors studied maximal existing patterns (i.e., existing patterns which are not a subset of another existing pattern) for k = 3, 4, 5. Cannings and Vickers conjectured that any subset of an existing pattern is also an existing pattern (Cannings and Vickers, 1989; Vickers and Cannings, 1988). Thus, if the conjecture is true, one can specify the complete set of existing patterns by specifying the maximal existing patterns. This is the approach adopted in the literature. The main purpose of this paper is to give a complete classification of all existing patterns for the four-allele model without assuming the truth of the conjecture of Cannings and Vickers. The complete list of existing patterns for the case k = 3 is given by the authors in Su et al. (2011). We should also point out that condition (3) is not posited in the literature. Therefore, patterns that are rotations of each other, which can be obtained by multiplying the fitness matrix R on both sides by a permutation matrix, are considered as one pattern. However, under condition (3), rotation of an existing pattern is not necessarily an existing pattern. Therefore, we treat patterns that are rotations of each other as different patterns. For these reasons, the number of existing patterns we obtained is greater than those obtained in Cannings and Vickers (1989) and Vickers and Cannings (1988). The method we use in this paper, counting signs of eigenvalues of the Jacobian matrix evaluated at the existing equilibria, is also completely different from the methods used in Cannings and Vickers (1989) and Vickers and Cannings (1988). This paper is organized as follows. In Section 2, we review several important properties of system (2). Some of the results in this section may be new, since we cannot find them in the literature. In Section 3, we summarize our results for the threeallele model in Su et al. (2011). In Section 4, we present our results for the four-allele model. In Section 5, we describe the domains of attraction for some of the patterns. Section 6 is discussions. In order not to interrupt the flow of ideas, many tables and some of the proofs are given in the Appendix.
Proposition 2.1. Assume that det R ̸= 0. An interior equilibrium exists if and only if all the Ui are non-zero and have the same sign. Moreover, if an interior equilibrium exists, it is uniquely given by Ui pi = , Uj
i = 1, . . . , k.
j
The degenerate case det R = 0 was considered by an der Heiden (1975) and Hughes and Seneta (1975). In this paper, we always assume that the determinants of R and all its submatrices are nonzero. A submatrix of R is a matrix obtained by erasing any number of rows and the corresponding columns of R. Next, we consider the stability of existing equilibria of system (2). It is well known that the mean fitness function r is strictly increasing along non-constant solutions of (2); see Kingman (1961b), Mulholland and Smith (1959), and Scheuer and Mandel (1959). Therefore, r plays the role of a Lyapunov function, and asymptotically stable interior equilibria are identified with the local maxima of r. Quite a few necessary and sufficient conditions for the stability of an existing interior equilibrium can be obtained by studying the properties of local maxima of r; see Hughes and Seneta (1975), Mandel (1970) and the references therein. Another more direct approach is to investigate the eigenvalues of the Jacobian matrix J evaluated at an existing equilibrium of system (2). The difficulty with this approach is the complexity of finding the relationships between the fitness matrix R = (rij ), the existence of an equilibrium, and the signs of the eigenvalues of J evaluated at an equilibrium in a k-allele model. Besides finding and classifying all patterns for the four-allele model, another main contribution of this paper is to show readers how to overcome this difficulty with the help of computer simulations. Let p∗ = (p∗1 , . . . , p∗k ) be an equilibrium. We write K for the set {1, . . . , k}, and define L = {i ∈ K : p∗i > 0} and L¯ = {i ∈ K : p∗i = 0}. The following proposition may be known, but we are unable to find a reference for it, so we present a proof here.
Proposition 2.2. Let l ∈ K , and replace pl with 1 − i̸=l pi in (2). Then the ith row, jth column entry of the Jacobian matrix J of (2) (as a system of (k − 1) equations) evaluated at an equilibrium point p∗ is given by Jij = δij (rˆi − r )|p∗ + pˆ∗ [(rˆiˆj − rˆil ) − 2(rˆj − rl )|p∗ ], i i, j = 1, . . . , k − 1, where ˆi = i if 1 ≤ i < l and ˆi = i + 1 if i ≥ l; namely,
Jij =
pˆ∗ [(rˆiˆj − rˆil ) − 2(rˆj − rl )|p∗ ]
if ˆi ∈ L,
δij (rˆi − r )|p∗
if ˆi ∈ L¯ .
i
In particular, if p∗ is an interior equilibrium, then Jij = pˆ∗ (rˆiˆj − rˆil ), i
i, j = 1, . . . , k − 1.
Proof. Jij = ∂∂p [pˆi (rˆi − r )]|p∗ = δˆiˆj (rˆi − r )|p∗ + pˆ∗ ∂∂p (rˆi − r )|p∗ . i ˆj ˆj
On the other hand, r =
Now rˆi =
rpˆ = j
ℓ̸=l
ℓ̸=l rˆiℓ pℓ + rˆil
1−
=
ℓ̸=l
ℓ̸=l
pℓ . Therefore, (rˆi )pˆ = rˆiˆj − rˆil . j
(rℓˆj − rℓl )pℓ + rˆj +
k
ℓ=1
pl − rl
rlℓ pℓ pˆ j
= rˆj − rl − rˆjl pl + rll pl + rˆj + rlˆj pl − rll pl − rl = 2(rˆj − rl ).
ℓ̸=l pℓ . Therefore,
(rℓ )pˆj pℓ + rˆj + (rl )pˆj pl − rl
ℓ̸=l rℓ pℓ + rl 1 −
2. Theory Results in this section are valid for any number of alleles. We begin with the problem of the existence of interior equilibria. Let Rij be the cofactor of rij in the fitness matrix R, and let Ui = j R ij . The following result is due to Mandel (1959).
(5)
L. Su, R. Lui / Theoretical Population Biology 81 (2012) 273–283
Noting that ri = r for i ∈ L, and at an interior equilibrium that r1 = · · · = rk , the conclusion of the proposition follows. Remark 2.3. If we choose l to be in L, then
Jij =
pˆ∗ (rˆiˆj − rˆil )
if ˆi, ˆj ∈ L,
δij (rˆi − r )|p∗
if ˆi ∈ L¯ .
i
This implies that the eigenvalues of J consist of (ri − r )|p∗ (i ∈ L¯ ) and the eigenvalues of the matrix (p∗i (rij − ril ))i,j∈L−{l} . This formula greatly reduces the amount of work when we compute the eigenvalues of J. For example, for a five-allele model, if we want to find the eigenvalues of J at the boundary equilibrium on the side where p2 = 0, p5 = 0, all we have to do is find the interior equilibrium of the three-allele model with fitness matrix obtained by eliminating the second and fifth rows and columns of R. Let this interior equilibrium be denoted by (¯p1 , p¯ 2 , p¯ 3 ), and let p∗ = (¯p1 , 0, p¯ 2 , p¯ 3 , 0). Then, according to the above lemma with l = 4, the eigenvalues of J evaluated at p∗ consists of the eigenvalues of the matrix p¯ 1 (r11 − r14 ) p¯ 2 (r31 − r34 )
p¯ 1 (r13 − r14 ) p¯ 2 (r33 − r34 )
Remark 2.4. The eigenvalues of J evaluated at the boundary equilibrium p∗ = (q∗1 , . . . , q∗i−1 , 0, q∗i , . . . , q∗k−1 ) of a k-allele model with fitness matrix R consist of the k − 2 eigenvalues of the Jacobian matrix evaluated at q∗ = (q∗1 , . . . , q∗k−1 ) and ri · p∗ − R p∗ · p∗ . It turns out that this last eigenvalue is related to the interior equilibrium of the k-allele model by the formula ri · p∗ − R p∗ · p∗ = −
Ui D
,
(6)
where D is the sum of all the entries of the adjugate matrix of R1 , R1 being the matrix obtained from erasing the ith row and column of R. We do not use (6) in this paper. Let
= p∗i (rij − 2r )|p∗ . Let Γ = (Γij ) ≡ (rij − 2r ) p∗i p∗j , and let M = diag(1/ p∗1 , . . . , 1/ p∗k ). It is clear that J˜ = M −1 Γ M. Since Γ is symmetric, all the eigenvalues of J˜ are real. Finally, we observe that, if λ is an eigenvalue of J with eigenvector v = (v1 , . . . , vk−1 ), then λ is also an eigenvalue of J˜ with eigenvector v˜ = (v1 , . . . , vk−1 , −v1 −· · ·− vk−1 ). To see this, let Fi (p1 , . . . , pk ) = pi (ri − r ) for i = 1, . . . , k. Then Jij =
∂ Fi ∂ Fi − ∂ pj ∂ pk
for i, j = 1, . . . , k − 1
and k−1 ∂ Fi
∂ pj
−
∂ Fi ∂ pk
vj = λvi for i = 1, . . . , k − 1.
Therefore, for i = 1, . . . , k − 1, we have k−1
J˜ij vj + J˜ik (−v1 − · · · − vk−1 ) =
j =1
k−1 (J˜ij − J˜ik )vj j =1
=
k−1 ∂ Fi j =1
∂ pj
−
∂ Fi ∂ pk
vj = λvi .
Now, consider the function G(p1 , . . . , pk ) = i=1 Fi (p1 , . . . , pk ) defined in Rk . The set ∆ may be thought of as the plane determined by p1 + p2 + · · · + pk − 1 = 0 in the positive quadrant of Rk . Since G = 0 on ∆, ∇ G is parallel to the normal of ∆, which is the k-dimensional vector e = ⟨1, . . . , 1⟩. This implies that ∂ G/∂ p1 = · · · = ∂ G/∂ pk . From above, for i = k, we have
k
k−1
J˜kj vj + J˜kk (−v1 − · · · − vk−1 )
j =1
k−1 k−1 ∂ Fk ∂ Fk (J˜kj − J˜kk )vj = − vj ∂ pj ∂ pk j =1 j=1 k−1 k−1 k−1 ∂G ∂ Fl ∂G ∂ Fl = − − − vj ∂ pj ∂ pj ∂ pk ∂ pk j =1 l=1 l =1 k−1 k−1 ∂ Fl ∂ Fl =− − vj = λ(−v1 − · · · − vk−1 ). ∂ pj ∂ pk l=1 j=1
=
∆=
∂ ∂ [pi (ri − r )]|p∗ = δij (ri − r )|p∗ + p∗i (ri − r )|p∗ ∂ pj ∂ pj ∂ ri ∂ r ∗ = pi = p∗i (rij − 2rj |p∗ ) − ∂ pj p∗ ∂ p j p∗
J˜ij =
j =1
and ri · p∗ − R p∗ · p∗ , i = 2, 5, where ri is the ith row of the fitness matrix R. This decomposition was also used in Kingman (1961a) to study the stability of boundary equilibria for the discrete-time model.
275
(p1 , . . . , pk−1 ) : 0 ≤ pi ≤ 1 for i = 1, . . . , k − 1, k−1
pi ≤ 1 .
(7)
i=1
It can be shown that ∆ is invariant with respect to the flow defined by (2), and, given initial conditions in ∆, the solutions of (2) converge to an asymptotically stable equilibrium in ∆ as t → ∞ (Bürger, 2000, p. 38). The set ∆ contains a maximum possible 2k − 1 equilibria. We now prove that the eigenvalues of the Jacobian matrix J evaluated at any of these equilibria are real. The discrete-time version of this theorem may be found in Kingman (1961a, Section 3). Proposition 2.5. All the eigenvalues of the Jacobian matrix J at an equilibrium are real. Proof. Because of Proposition 2.2, it suffices to prove the result at an interior equilibrium. Let J be the Jacobian matrix of (2) with pk ˜ replaced with 1 − i̸=k pi , and let J be the Jacobian matrix of (2) as a system of k equations. At an interior equilibrium, all the rj are equal to r. Therefore,
˜ Thus, all the eigenvalues This shows that λ is also an eigenvalue of J. of J are real, and the proof of the proposition is complete. Proposition 2.6. The determinant of J evaluated at an interior equilibrium p∗ = (p∗1 , . . . , p∗k ) is given by det(J ) = p∗1 p∗2 · · · p∗k (U1 + U2 + · · · + Uk ).
(8)
Proof. Let ri , i = 1, . . . , k − 1, be the ith row of the matrix R without the last element rik , and let e be the row vector with (k − 1) 1s. From the formula for J at the interior equilibrium p∗ in Proposition 2.2, we have
det(J ) = p∗1 p∗2 · · · p∗k−1 det
r1 − r1k e r2 − r2k e
. ··· r(k−1) − r(k−1)k e
276
L. Su, R. Lui / Theoretical Population Biology 81 (2012) 273–283
Rule 1. For any k, if the interior equilibrium is asymptotically stable, then no other equilibrium is asymptotically stable.
To continue,
det
r1 − r1k e r2 − r2k e
Rule 2. The maximum number of coexisting asymptotically stable equilibria in a three-allele model is 3, which must be the pattern consisting of the three vertices.
··· r(k−1) − r(k−1)k e
r1 r2 − r2k e
e
Rule 3. In a three-allele model, if all the existing asymptotically stable equilibria are vertices, then they must include V1 .
r2 − r2k e − r1k det ··· ··· r(k−1) − r(k−1)k e r(k−1) − r(k−1)k e
= det
e
r1 r2 − r2k e
= det
r2 − r1k det · · ·
···
= det
e r2
r1 r2
4. Four-allele model
r(k−1)
r(k−1) − r(k−1)k e
Rule 4. In a three-allele model, if vertex V2 or V3 is asymptotically stable but not V1 , then P2 or P3 must also be asymptotically stable, respectively.
r1 r2
− · · · − r(k−1)k det · · · − r1k det ··· ··· r(k−1)
r(k−1) r1 r2
. = det .. r
(k−1)
e
rk−2 e
r1k r2k
.. .
r(k−1)k 1
= Uk = p∗ (U1 + U2 + · · · + Uk ). k
The second to last equality is due to the definition of Ui given before Proposition 2.1, and the last equality follows from (5). Thus, formula (8) follows. The following two results are helpful when we show the nonexistence of certain patterns. Theorem 2.7 (Kingman (1961a)). Suppose that ∆1 ⊂ ∆2 are two distinct invariant subsets of ∆ and that the corresponding interior equilibria exist; then they both cannot be asymptotically stable. In particular, if ∆ contains an asymptotically stable interior equilibrium, then no other equilibrium in ∆ can be asymptotically stable, and, therefore, the interior equilibrium is globally asymptotically stable. Theorem 2.8 (Kingman (1961a) and Mandel (1959)). If the interior equilibrium p∗ is globally asymptotically stable, then the mean fitness r has a global maximum there. Consequently, r |p∗ ≥ rii for i = 1, . . . , k, which is obtained by evaluating r at the vertex (0, . . . , 0, 1, 0, . . . , 0). 3. Three-allele model For the three-allele case, (2) can be written as p˙ i = pi (ri − r ),
i = 1, 2.
(9)
The invariant set ∆ defined by (7) is the triangle joining the vertices V1 = (1, 0), V2 = (0, 1) and V3 = (0, 0) in the (p1 , p2 )plane. There are seven possible equilibria in ∆, comprising three vertices, three boundary equilibria on the edges, and one interior equilibrium. The edge Si = {pi = 0} represents a two-allele case, and the interior equilibrium on that edge is denoted by Pi . In Su et al. (2011), we showed that there are 14 existing patterns, but none of them contains more than three asymptotically stable equilibria. We discovered four rules for the three-allele patterns, and they are summarized below. These rules will be used in the study of the four-allele case. Rule 1 follows from Theorem 2.7, and is valid for any number of alleles. Three-allele patterns are also studied in Mandel (1959) and Wright (1969, p. 46) without assuming condition (3).
For a higher number of alleles, there are too many possible patterns, and we cannot prove their existence as we did for the three-allele model (see Su et al. (2011)). We therefore resort to computer simulations to find examples of existing patterns and, with the help of the above four rules, prove that the patterns that do not show up in our computer simulations cannot exist. For our computer simulations, we wrote a Matlab program that randomly generates a fitness matrix R satisfying condition (3), and we use the method in Remark 2.3 to compute and count the number of negative eigenvalues of J evaluated at each existing equilibrium. We then write the results to a text file and run a sorting algorithm on it to display all the distinct patterns. We repeat this process several times with increasing number of simulations until no new patterns show up. Examples of patterns containing two edge equilibria and one face equilibrium are the most difficult to generate, and they require a more targeted approach which will not be described here. We introduce some notation. The invariant set ∆ defined by (7) is a tetrahedron. It contains 15 possible equilibria, comprising four vertices, six edge equilibria, four face equilibria, and one interior equilibrium. Each edge of ∆ is a two-allele model and each face of ∆ is a three-allele model. We denote by v a vertex equilibrium, e an edge equilibrium, f a face equilibrium, and i the interior equilibrium. A pattern containing α vertices, β edge equilibria, and γ face equilibria is denoted by αv + β e + γ f , where α, β, γ are nonnegative integers. We say that three edge equilibria are coplanar if the edges they lie on form the sides of a triangle. From Rule 2, coplanar edge equilibria cannot all be asymptotically stable. Also, we shall prove that the maximum number of coexisting asymptotically stable equilibria is 4, which we state as Rule 5 below. From Rule 1, if the interior equilibrium is asymptotically stable, then no other equilibrium can be asymptotically stable. Based on these properties, we list all the possible patterns for the four-allele model in Table 1. The existing ones are confirmed by examples from simulations (see Table 2 in Appendix A for details) and the rest are shown not to exist below. In each line of Table 1, Ns denotes the number of asymptotically stable equilibria and Np denotes the number of patterns of that type. If we add up the numbers in the column Np , we see that there are 117 existing patterns for the four-allele model. We also see from Table 1 that there are only two types of existing pattern with four asymptotically stable equilibria, one of which is the pattern containing all four vertices and the other is the pattern containing four asymptotically stable non-coplanar edge equilibria. Proof of Table 1. Appendix A gives examples of fitness matrices R obtained from simulations which produce the existing patterns in Table 1. Rule 1 is the most useful in ruling out the existence of a certain pattern. For example, the pattern 2v + f is not possible because, if two vertices are asymptotically stable, then any face in ∆ must contain at least one of the two asymptotically stable vertices, and this violates Rule 1 for three-allele case. Also, Rule 1 for
L. Su, R. Lui / Theoretical Population Biology 81 (2012) 273–283
277
Table 1 Four-allele model: Total 117 existing patterns. Ns
Pattern type
v
Np e
f
1 1
1
1 1 2
1 1
15 2
6
1 1
12 4 12
1 1
3
3 3 3
3
2 2 1
1 6 4 1 3
2 2
Explanations
i
1 1 2
must be (1, 0, 0), Rule 6 any of the six edges any of the four faces the interior equilibrium one of the two vertices must be (1, 0, 0), Rule 6
any two edges: any two faces:
1 1
1 1
1 2 2 1
4 4 4 3 3 1 4
1 1 3 3
1 2 2 2 1 1
1 2 2 1 2 1
1 3 3 2 2 1 1 2
the four-allele case states that, if the interior equilibrium of the tetrahedron is asymptotically stable, then no other equilibrium can be asymptotically stable. Therefore, we only list one pattern in Table 1 with the interior equilibrium i. The patterns 3f , e + 2f , 3e + f , and 4f were shown to be not possible in Vickers and Cannings (1988, Theorems 4 and 6), and we call this reason A in the explanation column of Table 1. We also use Rule 5 and Rule 6 below in establishing Table 1. We prove Rule 5 here and postpone the proof of Rule 6 to Appendix B.
4 2
= 15
=6
any edge with a vertex not in the edge: 6 × 2 = 12 pair any face with the vertex not in the face any face with an edge not in the face: 4 × 3 = 12 one of the three vertices must be (1, 0, 0), Rule 6 6 3
− 4 = 16
16
any three edges not coplanar:
0 6 0 12
does not exist because of reason A pair any edge with the two vertices not in the edge does not exist because of Rule 1 anyvertex with two edges not containing the vertex: 4 1
2
6
2
12 0 0 0 1 3 0 0 0 0 0 0 0 0 0 0 0 0 0
×
3 2
= 12
any face with two of the three edges not from the face does not exist because of Rule 1 does not exist because of reason A does not exist because of Rule 1 all four vertices any four edges without forming a triangle does not exist because of reason A does not exist because of Rule 1 does not exist because of Rule 1 does not exist because of Rule 1 and Rule 2 does not exist because of reason A does not exist because of Rule 1 does not exist because of Rule 1 does not exist because of Rule 1 does not exist because of Rule 1 does not exist because of Rule 1 does not exist because of Rule 1 does not exist because of Rule 1 does not exist because of Rule 1
we will violate the rule mentioned above that there cannot be coplanar asymptotically stable edge equilibria. The proof of Rule 5 is complete. Remark 4.1. From Table 1 we observe that, if ∆ contains three or more asymptotically stable equilibria, then none can be a face equilibrium except for the pattern 2e + f , where the two asymptotically stable edge equilibria do not lie on the face. This observation is useful in the study of the five-allele model.
Rule 5. The maximum number of coexisting asymptotically stable equilibria for the four-allele model is 4.
5. Domains of attraction
Rule 6. In a four-allele model, if all the asymptotically stable equilibria are vertices, then they must include the vertex (1, 0, 0). To prove Rule 5, first recall that ∆ has four vertices, six edges, and four faces. Suppose we have a pattern with five or more asymptotically stable equilibria. Then, by Rule 1, the pattern must not contain the interior equilibrium of ∆. Suppose that there is an asymptotically stable vertex; then, from Rule 1, all the other asymptotically stable equilibria are in the face opposite the vertex. From Rule 2, there can be at most three asymptotically stable equilibria in this face, so this is impossible because we need four more asymptotically stable equilibria. Suppose that a face equilibrium is asymptotically stable; then, according to Rule 1, the other asymptotically stable equilibria must be either be the vertex opposite the face, or the three edge equilibria not in the face. For the same reason, this is again impossible. Finally, if all the asymptotically stable equilibria are on the edges of ∆, then
Knowing all the patterns is not enough to predict the asymptotic behavior of solutions of (2). Given initial data, we also want to know which asymptotically stable equilibrium the solutions of (2) approach as t → ∞. The domain of attraction of an asymptotically stable equilibrium is the subset of ∆ such that solutions of (2) with initial data in this subset converge to the asymptotically stable equilibrium as t → ∞. For the three-allele model, the domains of attraction are separated by invariant curves, called separatrices, which are not analytically tractable. They are described in Theorem 5.5 of Su et al. (2011). For the four-allele model, the domains of attraction are formed by partitioning the tetrahedron ∆ by surfaces which are two-dimensional invariant manifolds. Each domain of attraction contains one asymptotically stable equilibrium. The intersection of the invariant manifolds with the four faces of ∆ are separatrices in the corresponding
278
L. Su, R. Lui / Theoretical Population Biology 81 (2012) 273–283
Fig. 1. Domains of attraction for the pattern 2v + e. M1 separates out the DA of v1 . On M1 , f2 is asymptotically stable. M2 separates out the DA of v2 . On M2 , f1 is asymptotically stable. The intersections of M1 and M2 are the separatrices connecting f3 , e4 and e4 , f4 , respectively. The remaining region in ∆ is the DA of e3 .
Fig. 2. Domains of attraction for the pattern v + 2e. M1 separates out the DA of v1 . On M1 , f3 is asymptotically stable. M2 separates out the DA of e6 . On M2 , f1 is asymptotically stable. The intersection of M1 and M2 are the separatrices connecting e4 , f4 , e5 and e5 , f2 , e3 , respectively. The remaining region in ∆ is the DA of e2 .
Fig. 3. Domains of attraction for the pattern 2e + f . M1 separates out the DA of f4 which is on the face {123}. On M1 , f3 is asymptotically stable. M2 separates out the DA of e3 . On M2 , f1 is asymptotically stable. The intersection of M1 and M2 is the separatrix connecting e1 , f2 , v3 . The remaining region in ∆ is the DA of e2 .
three-allele case. These separatrices may be identified by looking at Theorem 5.5 of Su et al. (2011) or by simulations. For each asymptotically stable equilibrium, we look at the separatrices that separate out the equilibrium on the faces where the equilibrium lies. The invariant manifolds are then put together by connecting these separatrices. In the following, we describe the domains of attraction for some of the existing patterns in the four-allele model with three or four asymptotically stable equilibria. If there is only one asymptotically stable equilibrium, the domain of attraction contains the interior of ∆, and, if there are two asymptotically stable equilibria, the domains of attraction are separated by an invariant manifold. These cases are easy to visualize, and they will not be described here. From Table 2, there are seven groups of patterns with three or four stable equilibria: 2v + e, v + 2e, 2e + f , 3v , 3e, 4v , and 4e. Because of space limitations, we only describe the domains of attraction for the first line of each group given in Table 2 (Figs. 1–7). For example, Fig. 1 below shows the domains of attraction for the pattern 2v + e with the fitness matrix r11 = 1, r12 = −4.544,
r13 = −0.714, r14 = −4.035, r22 = 0.065, r23 = −4.021, r24 = −2.943, r33 = 0.039, r34 = 4.420, and r44 = 0, which may be found in the first row of the group 2v + e. In the figures, the equilibria are denoted by v1 = {1}, v2 = {2}, v3 = {3}, v4 = {4}, e1 = {14}, e2 = {24}, e3 = {34}, e4 = {12}, e5 = {13}, e6 = {23}, f1 = {234}, f2 = {134}, f3 = {124}, f4 = {123}, and i = {1234}, which is the interior equilibrium. In Figs. 1–7, the asymptotically stable equilibria are indicated by dots and the unstable ones by open circles, except for unstable vertices. Solid curves represent separatrices on the face facing the reader and dotted curves represent separatrices on the other three faces of ∆. Note that an equilibrium which is unstable in ∆ may be stable on a face as an equilibrium in the three-allele case. This lowerdimensional stability is not reflected in the figures. The invariant manifolds together with the dynamics on them are also given. The figures are only schematic, because the invariant manifolds are in general analytically intractable. We use DA to denote a domain of attraction in the captions of the figures.
L. Su, R. Lui / Theoretical Population Biology 81 (2012) 273–283
279
Fig. 4. Domains of attraction for the pattern 3v . M1 separates out the DA of v3 . On M1 , e5 is asymptotically stable. M2 separates out the DA of v2 . On M2 , e2 is asymptotically stable. The remaining region in ∆ is the DA of v1 .
Fig. 5. Domains of attraction for the pattern 3e. M1 separates out the DA of e4 . On M1 , f4 is asymptotically stable. M2 separates out the DA of e1 . On M2 , f2 is asymptotically stable. The remaining region in ∆ is the DA of e5 .
Fig. 6. Domains of attraction for the pattern 4v . M1 separates out the DA of v3 . On M1 , e3 and e6 are asymptotically stable and their DA are separated by a separatrix connecting f2 and f1 . M2 separates out the DA of v1 . On M2 , e1 and e4 are asymptotically stable and their DA are separated by a separatrix connecting f2 and f3 .M3 separates out the DA of v2 and v4 from the remaining region of ∆. On M3 , e2 is asymptotically stable.
When we construct the domains of attraction, sometimes the choice of invariant manifolds is not obvious. This can be resolved by looking at the stable and unstable manifolds of the face equilibria. For example, in Fig. 5, suppose that M1 and M3 are the invariant manifolds rather than M1 and M2 shown in the figure, where M3 contains v1 , f4 , e6 , f1 , e3 , and f2 on its boundary. Then look at the face equilibrium f4 , which has two negative eigenvalues and one positive eigenvalue. The unstable direction corresponding to the positive eigenvalue lies on the face {123}. Furthermore, f4 is stable along the intersection of M1 and M3 , which is the separatrix on the face {123}, and f4 is stable on M1 , because no other equilibrium is stable on M1 . This implies that f4 must have an unstable direction on M3 , which is a contradiction, because it already has an unstable direction on the face {123}. One can also rule out the choice of M2 and M3 as the invariant manifolds by looking at the face equilibrium f2 . Therefore, in Fig. 5, M1 and M2 are the correct invariant manifolds which give the domains of attraction in ∆ (Figs. 6 and 7).
6. Discussions In this paper, we give all existing patterns for the four-allele selection model under conditions (3) and rij ̸= rii or rjj for i ̸= j. If, for some i, j, rij = rii or rjj , then the interior equilibrium on the edge {ij} disappears, and this is called complete dominance. Complete dominance reduces the total number of existing patterns in a k-allele model. Knowing all existing patterns is important, because it tells us the possible long-term genetic makeup of the population, that is, which alleles will survive and which will disappear under selection forces in the long run. The actual outcome will depend on the initial conditions of the population if there are multiple coexisting asymptotically stable equilibria. In such a case, the set ∆ is divided into domains of attraction, each containing one asymptotically stable equilibrium. The problem of finding all existing patterns for the k-allele model has been studied before, but it remains an open problem except when k = 2, 3, or 4. Our results are different from earlier work because we do not assume the truth of the conjecture of
280
L. Su, R. Lui / Theoretical Population Biology 81 (2012) 273–283
Fig. 7. Domains of attraction for the pattern 4e. M1 and M2 intersect at the invariant curve passing through e1 , i, and e6 . They divide ∆ into four subregions that correspond to the DA of e2 , e3 , e5 , and e4 , respectively. On M1 , f1 and f4 are asymptotically stable. On M2 , f2 and f3 are asymptotically stable.
Cannings and Vickers, that any subset of an existing pattern is also an existing pattern. In fact, this conjecture is not true in our case because, for example, the pattern ⟨{1}, {2}⟩ exists but not the pattern ⟨{2}⟩. Also, unlike earlier work, rotation of an existing pattern may not be an existing pattern. For example, ⟨{1}, {2}⟩ is an existing pattern but ⟨{2}, {3}⟩ is not. All these are consequences of condition (3), which may be satisfied by relabeling the alleles assuming that the rii are different from each other. In this paper, we study the continuous-time model instead of the discrete-time model (1), which allows us to adopt a completely different approach of counting signs of the eigenvalues of the Jacobian matrix evaluated at all existing equilibria. This approach motivates us to prove several interesting mathematical results, which are valid for any number of alleles. We believe that some of these results are new, because we could not find them in the literature. We have also made considerable progress towards solving the patterns problem for the five-allele model in Su et al. (2011). The number of existing patterns is around 2351. For the six-allele model, there will be more than 60,000 existing patterns, and we conjecture that the number of existing patterns grows exponentially as a function of k. It is worth pointing out that the two-locus two-allele selection model degenerates to the single-locus four-allele model in the absence of genetic recombination (Bürger, 2000, p. 47). Therefore, the four-allele model studied in this paper has a special role in population genetics as the reference model for the two-locus twoallele model. Acknowledgment We thank Allan E. Johannesen of the WPI Computing and Communications Center for writing the sorting algorithm. Appendix A. Examples of existing patterns In the following tables, we denote the vertices (1, 0, 0), (0, 1, 0), (0, 0, 1), (0, 0, 0) by {1}, {2}, {3}, {4}, respectively. {ij} denotes the interior equilibrium of the edge joining vertices {i} and {j}, and {ijk} denotes the interior equilibrium of the face with vertices {i}, {j}, and {k}. Note that, by scaling time, system (2) is unchanged under a linear transformation of the rij , so we may assume that r11 = 1 and r44 = 0. Appendix B. Proof of Rule 6 We first introduce the following notation: V1 = (1, 0, 0), V2 = (0, 1, 0), V3 = (0, 0, 1), V4 = (0, 0, 0), {ij} denotes the edge joining Vi and Vj , Pij denotes the edge equilibrium in {ij}, {ijk} denotes the face with vertices Vi , Vj , and Vk , and Pijk denotes the face equilibrium in {ijk}. It is easy to verify that the eigenvalues of
(j)
J evaluated at Vi are λi = rji − rii , j = 1, 2, 3, 4, j ̸= i. The matrix J evaluated at an edge equilibrium, say P13 , has three eigenvalues, one of which may be determined from applying Proposition 1.1, namely, P13 is asymptotically stable in {13} if r13 > r11 and unstable in {13} if r13 < r33 . Let the other two eigenvalues be denoted by λ213 and λ413 , where the superscript denotes that the vertex P13 is facing in the face where {13} is an edge. In this case, {13} lies on the faces {123} and {134}, so the vertices P13 is facing are V2 and V4 , respectively. We shall prove Rule 6 by showing that, if Vi , i ̸= 1, is asymptotically stable, then V1 must also be asymptotically stable. We only prove the case i = 2, since it is the most difficult. In the proof of the following lemma, we shall use Remark 2.3 repeatedly without mentioning it. Lemma B.1. Suppose that V2 is asymptotically stable in ∆, and that there is no non-vertex asymptotically stable equilibrium in ∆; then V1 must be asymptotically stable in ∆. Proof. Since V2 is asymptotically stable and r11 > r22 , {12} is in the inferior case, and hence V1 is asymptotically stable in {12}. We now focus on the invariant triangle {134}, which must contain at least one equilibrium asymptotically stable in {134}. V1 cannot be asymptotically stable in {134}, for otherwise, V1 is stable in ∆ and the lemma is proved. V3 or V4 (or both) cannot be the only asymptotically stable equilibrium (equilibria) in {134}, because Rule 3 of the three-allele case would force V1 to be asymptotically stable in {134}. Suppose that P13 exists. Then, since it is opposite to V2 in {123}, Rule 4 implies that P13 is asymptotically stable in {123}. Therefore, P13 cannot be asymptotically stable in {134}, because it will then be asymptotically stable in ∆, contradicting our assumption. Similarly, P14 cannot be asymptotically stable in {134}. (This argument does not work for P34 , since V2 is the most asymptotically stable vertex in {234}.) The above argument implies that either P134 or P34 is asymptotically stable in {134}. We now show that neither case can happen. From our assumption, V2 is asymptotically stable on the edges {12}, {23}, and {24}, so r22 > r12 , r32 , r42 . If P134 is asymptotically stable in {134}, then J evaluated at P134 = (¯p1 , 0, p¯ 3 ) has three eigenvalues, two of which are negative and the third equals
(r2 − r )|P134 < (r2 − r11 )|P134 = (r21 − r11 )¯p1 + (r23 − r11 )¯p3 + (r24 − r11 )(1 − p¯ 1 − p¯ 3 ) < (r22 − r11 )¯p1 + (r22 − r11 )¯p3 + (r22 − r11 )(1 − p¯ 1 − p¯ 3 ) = r22 − r11 < 0. The first inequality follows from Theorem 2.8. This implies that P134 must be asymptotically stable in ∆, contradicting our assumption.
L. Su, R. Lui / Theoretical Population Biology 81 (2012) 273–283
281
Table 2 Existing patterns of the four-allele model. Type
Pattern
v
⟨{1}⟩
e
⟨{12}⟩ ⟨{13}⟩ ⟨{14}⟩ ⟨{23}⟩ ⟨{24}⟩ ⟨{34}⟩
Values of fitness constants r12
r23
r24
0.645
0.866
0.230
0.572
−4.645
−3.889
0.080
−0.885
5.329 1.782 −3.167 3.279 1.810 −1.279
0.865 4.454 2.509 2.002 2.594 0.632
0.733 0.508 4.474 −1.299 1.870 3.383
0.944 0.256 0.074 0.712 0.010 0.606
3.858 1.835 0.588 5.881 −0.315 3.244
−2.076 −1.512
−3.868
0.835 4.022 5.732 2.483
0.716 0.230 0.063 0.689 0.005 0.114
1.856 0.873 −0.810 4.018 5.737
3.577 1.523 1.470 2.753
1.409 4.242 −2.998 2.020
−2.897 3.226 5.122 1.573
0.943 0.959 0.853 0.646
5.087
f
⟨{234}⟩ ⟨{134}⟩ ⟨{124}⟩ ⟨{123}⟩
1.264 1.555
4.180 2.970 4.589 −0.425
0.512 0.324 0.255 0.068
2.383 2.094 −2.331 −0.739
i
⟨{1234}⟩
4.057
5.051
3.359
0.950
2.632
3.038
0.930
3.735
2v
⟨{1}, {2}⟩ ⟨{1}, {3}⟩ ⟨{1}, {4}⟩
−0.205
0.793
0.608
−3.324
−0.032
−1.043
0.635
−1.045 −4.584
0.820 0.597 0.605
0.513
0.700 0.865
0.273 −0.647
0.156 0.432 0.451
−4.459 −0.770 −0.681
5.105 4.925 3.670 4.160 2.110 3.144 −1.632 −3.867 0.565 −2.815 −0.191 −2.305 2.414 −2.581 2.992
4.741 −1.871 −3.582 2.306 −0.275 3.972 4.761 1.418 4.954 −2.074 0.830 −4.559 1.039 1.685 −3.570
5.085 3.956 −0.182 −3.241 −3.338 3.528 −3.933 −1.383 −0.700 4.740 4.769 1.434 −1.627 3.234 1.022
0.925 0.688 0.737 0.511 0.910 0.221 0.314 0.406 0.068 0.932 0.195 0.150 0.782 0.742 0.508
−1.244
0.339
2.654 1.795 −0.146 −0.790 −0.186 5.533 −1.129 3.100 3.876 −3.193 1.200 2.959 4.508 −4.604
−3.474
0.326 0.397 0.288 0.260 0.835 0.058 0.230 0.063 0.051 0.097 0.162 0.127 0.625 0.558 0.408
−0.431 −4.996 −3.240 −4.437
2e
⟨{12}, {13}⟩ ⟨{12}, {14}⟩ ⟨{12}, {23}⟩ ⟨{12}, {24}⟩ ⟨{12}, {34}⟩ ⟨{13}, {14}⟩ ⟨{13}, {23}⟩ ⟨{13}, {24}⟩ ⟨{13}, {34}⟩ ⟨{14}, {23}⟩ ⟨{14}, {24}⟩ ⟨{14}, {34}⟩ ⟨{23}, {24}⟩ ⟨{23}, {34}⟩ ⟨{24}, {34}⟩
3.428 5.327 4.557 1.866 2.545 0.487
2.791 5.256 3.719 5.143 −2.247 2.829
3.783 5.264 2.460 −3.156 4.177 4.412
0.937 0.879 0.335 0.068 0.778 0.662
3.787 5.562 −3.421 4.382 4.242 5.165
0.691 0.186 0.245 0.052 0.457 0.360
−2.053
2f
⟨{123}, {124}⟩ ⟨{123}, {134}⟩ ⟨{124}, {134}⟩ ⟨{123}, {234}⟩ ⟨{124}, {234}⟩ ⟨{134}, {234}⟩
−0.402 −0.595 −0.044 −0.177 −4.643 −0.943
−0.709
−0.980 −4.050
v+e
⟨{1}, {23}⟩ ⟨{1}, {24}⟩ ⟨{1}, {34}⟩ ⟨{2}, {13}⟩ ⟨{2}, {14}⟩ ⟨{2}, {34}⟩ ⟨{3}, {12}⟩ ⟨{3}, {14}⟩ ⟨{3}, {24}⟩ ⟨{4}, {12}⟩ ⟨{4}, {13}⟩ ⟨{4}, {23}⟩
0.529 0.135 0.161 0.866 0.711 0.245 0.621 0.827 0.018 0.970 0.193 0.640
4.382 1.884 −0.945 −3.511 0.387 −3.275 −0.948 −3.813 −2.742 1.974 −1.388 3.392
⟨{1}, {234}⟩ ⟨{2}, {134}⟩ ⟨{3}, {124}⟩ ⟨{4}, {123}⟩
−3.360 −4.753
0.031 0.123 0.031 0.466
5.190
2.986
v+f
−2.672 −1.374
−0.857
1.259
0.415 0.299 0.595 0.919 0.392 0.830 0.019 0.130 0.819 0.892 0.181 0.132 0.065 0.690 0.891 0.787
e+f
2v + e
5.506 2.486 0.018 4.516 1.736 2.207
r13
0.592 −3.657 1.256 2.864 0.478 −0.194 −4.282 −4.984 1.478 4.036 1.238
r14
0.619 0.700 5.189 1.760 0.453 5.975 2.035 −3.825 −2.029 −0.844
−2.938
−2.309
1.268 3.208
3.338 −2.370 3.114
4.620 4.917 −3.204
5.807
−1.707
−1.412 −3.565
−0.828 −3.994
⟨{234}, {12}⟩ ⟨{234}, {13}⟩ ⟨{234}, {14}⟩ ⟨{134}, {12}⟩ ⟨{134}, {23}⟩ ⟨{134}, {24}⟩ ⟨{124}, {13}⟩ ⟨{124}, {23}⟩ ⟨{124}, {34}⟩ ⟨{123}, {14}⟩ ⟨{123}, {24}⟩ ⟨{123}, {34}⟩
5.417 −3.541 0.175 5.747 3.471 4.269 3.777 2.964 4.350
4.860 −3.540 4.824 1.816 3.855 5.669 −0.337 −0.371 2.546 3.140 1.870
⟨{1}, {2}, {34}⟩ ⟨{1}, {3}, {24}⟩ ⟨{1}, {4}, {23}⟩ ⟨{2}, {3}, {14}⟩
−4.544 −1.377 −4.799 −4.759
−0.714 −1.810 −4.711 −0.625
2.764 1.146 4.658 4.375 2.810 4.120 4.200 3.052 −3.268 −3.828
−4.035 −4.566 −1.911 1.630
r22
−0.423
0.701
1.059 5.583 −2.363 −3.856 4.697 5.490 −2.028 −4.505 5.197 3.455 2.954 −1.156 5.856 1.397
−4.202 5.428 0.693 5.176 4.116
−1.167 4.055 4.747 −3.083 −1.521 −3.205 −0.280 −1.136 5.523 −2.020 −1.072 −0.165
r33
0.142 0.076 0.084 0.644 0.097 0.105 0.049 0.177 0.017 0.442 0.157 0.182
r34
3.423
−4.934 0.964 3.130 4.518 2.463 3.137 4.797 −1.275 5.243 5.579 3.557 5.089 4.593 2.569 2.407 1.201
−0.423 5.914
−1.309 1.794 3.167 −1.029 −0.038 −3.242 −1.784 −0.136 −2.653
−0.765
0.013 0.040 0.013 0.168
4.376 4.580 −3.918 −3.012
1.719 3.181 5.357 −0.548 5.808 −0.039 −3.574 3.856 −3.336 3.773 2.270 3.590
3.836 3.559 5.038 −4.371 −4.131 4.912 5.682 0.924 2.356 −3.293 5.148 −4.576
0.034 0.016 0.553 0.287 0.301 0.147 0.017 0.074 0.069 0.404 0.117 0.118
3.853 5.052 1.344 5.504 3.626 1.051 1.335 −1.162 5.410 −1.993 2.336 1.429
−4.021 −0.114
−2.943
0.039 0.604 0.887 0.102
−4.368 −0.805 −4.345
3.840
−1.712
4.771
4.316 −3.124 −1.841
4.420
(continued on next page)
282
L. Su, R. Lui / Theoretical Population Biology 81 (2012) 273–283
Table 2 (continued) Type
v + 2e
Pattern
Values of fitness constants r12
r13
r14
r23
r24
⟨{2}, {4}, {13}⟩ ⟨{3}, {4}, {12}⟩
−3.766
5.156 −0.990
−3.753 −4.620
0.371 0.654
−4.471 −2.874
−4.570 −3.153
0.337 0.173
−2.315 −3.952
⟨{1}, {23}, {24}⟩ ⟨{1}, {23}, {34}⟩ ⟨{1}, {24}, {34}⟩ ⟨{2}, {13}, {14}⟩ ⟨{2}, {13}, {34}⟩ ⟨{2}, {14}, {34}⟩ ⟨{3}, {12}, {24}⟩ ⟨{3}, {12}, {14}⟩ ⟨{3}, {24}, {14}⟩ ⟨{4}, {23}, {12}⟩ ⟨{4}, {23}, {13}⟩ ⟨{4}, {12}, {13}⟩
−0.227 −1.856 −0.404
−3.813
0.979
2.773 4.171 −4.144 −4.437 −3.204 −2.479 −3.237 −3.274 −4.936 3.528 2.984 −4.099
3.606 −1.511 −0.495 0.296 5.655 −1.450 0.761 −4.017 −0.919 −3.797
0.300 0.660 0.009 0.270 0.143 0.049 0.685 0.042 0.052 0.375 0.161 0.329
−3.084
−3.799 −4.453
0.325 0.688 0.015 0.420 0.210 0.366 0.948 0.178 0.063 0.673 0.799 0.528
5.467
0.803 −0.344 4.287 4.494 −2.732 −3.793 −3.254 0.015 −1.499 3.847 5.188 1.500 1.895 −0.877 2.025 1.199 0.673 1.467 −0.041 5.628 0.901 1.392 1.111
−4.312
0.955 0.968 0.968 0.827 0.746 0.234 0.968 0.968 0.238 0.968 0.062 0.996
0.831 1.320 1.195 6.382 −8.892 6.965 −1.001 1.570 6.045 1.225 2.113 0.889
1.038 1.424 0.703 6.720 2.122 −8.108 1.415 0.704 2.724 −1.032 2.603 1.444
0.283 0.538 0.538 0.707 0.260 0.156 0.538 0.538 0.077 0.538 0.052 0.555
2.471 0.305 1.747 0.396 1.952 2.797 1.933 2.035 −8.365 1.859 −3.878 1.889
0.535 0.226 0.377
−3.852 −1.163 −4.284
−1.132 −4.188 −0.057
0.352 0.057 0.277
−2.018 −0.069 −1.120
0.395 0.237 0.130 0.484 0.205 0.541 0.548 0.682 0.607 0.766 0.551 0.860 0.311 0.111 0.926 0.519
−2.432
−2.823
2.979 5.705 −4.192 −3.617 4.652 −3.643 3.720 −3.254 −3.190 2.382 4.433 −2.043 1.999 −2.988 2.134
4.454 −2.352 4.009 2.338 −2.035 0.718 1.513 5.078 −2.231 −1.788 4.955 4.430 2.645 5.400 −2.193
0.141 0.049 0.009 0.454 0.008 0.186 0.161 0.605 0.038 0.238 0.097 0.026 0.192 0.033 0.762 0.277
−3.690 −1.943 4.197 2.014 −1.234 1.321 2.771 −3.446 −4.606 2.127 2.871 −1.110 3.784 −3.929 3.009 1.576
1.592
0.261
−1.023 −3.143 3.988 2.237 −2.002 1.574 −3.923 5.745
5.814
−4.302 2.570
−3.317 1.893 2.758 −4.840 −2.714 −1.322
r22
−3.201
r33
r34
3.402 3.986 −2.316 5.373 3.899 −3.816 −4.219 −2.662 −2.777 −2.127 −4.514
2e + f
⟨{123}, {24}, {34}⟩ ⟨{134}, {23}, {24}⟩ ⟨{124}, {23}, {34}⟩ ⟨{234}, {13}, {14}⟩ ⟨{124}, {13}, {34}⟩ ⟨{123}, {14}, {34}⟩ ⟨{134}, {12}, {24}⟩ ⟨{234}, {12}, {14}⟩ ⟨{123}, {24}, {14}⟩ ⟨{134}, {23}, {12}⟩ ⟨{124}, {23}, {13}⟩ ⟨{234}, {12}, {13}⟩
3v
⟨{1}, {2}, {3}⟩ ⟨{1}, {2}, {4}⟩ ⟨{1}, {3}, {4}⟩
−2.741 −2.121
3e
⟨{12}, {13}, {14}⟩ ⟨{12}, {23}, {24}⟩ ⟨{13}, {23}, {34}⟩ ⟨{14}, {24}, {34}⟩ ⟨{12}, {13}, {24}⟩ ⟨{12}, {14}, {23}⟩ ⟨{13}, {12}, {34}⟩ ⟨{13}, {14}, {23}⟩ ⟨{14}, {13}, {24}⟩ ⟨{14}, {12}, {34}⟩ ⟨{23}, {12}, {34}⟩ ⟨{23}, {13}, {24}⟩ ⟨{24}, {12}, {34}⟩ ⟨{24}, {23}, {14}⟩ ⟨{34}, {13}, {24}⟩ ⟨{34}, {14}, {23}⟩
2.912 2.128 −3.349 −4.071 4.656 1.096 3.267 −4.557 −0.518 2.617 1.367 −0.215 2.871 −4.427 0.158 −1.296
4.459 2.709 1.553 −4.102 −3.901 1.628 −3.588 −4.584 3.746 −2.321
4v
⟨{1}, {2}, {3}, {4}⟩
−2.828
−4.664
−3.314
0.890
−1.055
−0.512
0.463
−1.505
1.753 5.325 −3.728
2.081 −3.018 5.576
−4.961
0.128 0.776 0.573
−1.141
4e
⟨{12}, {13}, {24}, {34}⟩ ⟨{12}, {23}, {14}, {34}⟩ ⟨{13}, {23}, {14}, {24}⟩
2.237 −2.682 2.921
0.110 0.014 0.489
2.767 5.005 −4.481
1.332
−0.273 1.435
−5.116 1.849 6.488 1.020 1.063 0.680 1.027 0.724 1.022
0.524
2.049 1.715 2.252 0.815 2.249 0.835 1.169 2.217 1.504 1.959 0.550
−0.102
0.155
0.664 −1.133
−2.384 −1.795
4.833
3.753
−0.650
−0.636 −3.535
5.830
−3.325 2.506
−1.353
3.025 −4.977 1.398 −3.592 5.346 2.757 2.848 0.531 −4.676 −3.237 5.527 −4.351 1.368
5.653 3.137
Thus, P34 must be asymptotically stable in {134}. This implies that
{34} is in the superior case and λ234 > 0 > λ134 .
(10)
Since V1 is unstable in ∆ but asymptotically stable in {12}, it must be unstable in {13} or {14} or both. We have three cases to consider. (i) V1 is unstable on {13} and {14} (r13 > r11 , r14 > r11 ). Let P34 = (0, 0, p¯ 3 ). Then
λ234 − λ134 = (r2 − r1 )|P34 = (r23 − r13 )¯p3 + (r24 − r14 )¯p4 < (r22 − r11 )¯p3 + (r22 − r11 )¯p4 < 0.
(11)
Thus, λ234 < λ134 , which contradicts (10). In (11), we have used the fact that V2 is asymptotically stable in ∆. (ii) V1 is unstable in {13} but asymptotically stable in {14} (r13 > r11 , r14 < r11 ). Let P13 = (ˆp1 , 0, pˆ 3 ). Then
λ413 = (r4 − r )|P13 = (r4 − r1 )|P13
4.160 3.979
= (r41 − r11 )ˆp1 + (r43 − r13 )ˆp3 .
(12)
Since P13 is unstable in {134} but asymptotically stable in {13} and {123}, we have λ413 > 0. From (12) being positive and the fact that we are in case (ii), we have r43 > r13 > r11 . Altogether, we have established the inequalities r43 > r13 > r11 > r22 > r23 .
(13)
To continue, (10), (13), and the third term in (11) is positive imply that r24 > r14 . Now r = r3 = r4 at P34 . This and the fact that P34 is unstable in {234} imply that
λ234 = (r2 − r3 )|P34 = (r23 − r33 )¯p3 + (r24 − r34 )¯p4 > 0.
(14)
From (13) and V2 being asymptotically stable in ∆, we have r24 < r22 < r34 . Hence, (14) implies that r33 < r23 . From Proposition 1.1, we have p¯ 3 = (r44 − r34 )/∆34 and p¯ 4 = (r33 − r34 )/∆34 . Substituting these values into the third term of
L. Su, R. Lui / Theoretical Population Biology 81 (2012) 273–283
283
(11), then, from (10) and ∆34 < 0 (recall that {34} is in the superior case), we have
References
(r13 − r23 )(r34 − r44 ) < (r24 − r14 )(r34 − r33 ).
an der Heiden, U., 1975. On manifolds of equilibria in the selection model for multiple alleles. J. Math. Biol. 1, 321–330. Broom, M., Cannings, C., Vickers, G.T., 1993. On the number of local maxima of constrained quadratic form. Proc. R. Soc. Lond. A 443, 573–584. Bürger, R., 2000. The Mathematical Theory of Selection, Recombination, and Mutation. Wiley, Chichester. Cannings, C., Vickers, G.T., 1988. Patterns of ESS’s II. J. Theoret. Biol. 132, 409–420. Cannings, C., Vickers, G.T., 1989. Patterns and invasions of evolutionarily stable strategies. J. Appl. Math. Comput. 32, 227–253. Cannings, C., Vickers, G.T., 1991. The genealogy of patterns of ESS’s. Selected proceedings of the Sheffield symposium on applied probability. IMS Lecture Notes Monogr. Ser. 18, 193–204. Hughes, P.J., Seneta, E., 1975. Selection equilibria in a multi-allele single-locus setting. Heredity 35, 185–194. Kingman, J.F.C., 1961a. A mathematical problem in population genetics. Proc. Cambridge Philos. Soc. 57, 574–582. Kingman, J.F.C., 1961b. On an inequality in partial averages. Quart. J. Math. 12, 78–80. Mandel, S.P.H., 1959. The stability of a multiple allelic system. Heredity 13, 289–302. Mandel, S.P.H., 1970. The equivalence of different sets of stability conditions for multiple allelic systems. Biometrics 26, 840–845. Mulholland, H.P., Smith, C.A.B., 1959. An inequality arising in genetical theory. Amer. Math. Monthly 66, 673–683. Nagylaki, T., 1992. Introduction to Theoretical Population Genetics. Springer, New York. Scheuer, P.A.G., Mandel, S.P.H., 1959. An inequality in population genetics. Heredity 13, 519–524. Su, L., Sesanker, C., Lui, R., Coexisting stable equilibria in a multiple-allele population genetics model, arXiv:1108.5110v1 [q-bio.PE], http://arxiv.org/abs/1108.5110, 2011. Vickers, G.T., Cannings, C., 1988. Patterns of ESS’s I. J. Theoret. Biol. 132, 387–408. Wright, S., 1969. Evolution and the Genetics of Populations. In: The Theory of Gene Frequencies, vol. 2. University Press, Chicago.
(15)
Substituting the above formulas for p¯ 3 , p¯ 4 into (14) and using the fact that λ234 > 0, we have
(r34 − r24 )(r34 − r33 ) < (r23 − r33 )(r34 − r44 ).
(16)
It can be verified that the left-hand sides of inequalities (15) and (16) are positive. Multiplying these two inequalities, and upon cancellations, we have
(r13 − r23 )(r34 − r24 ) < (r24 − r14 )(r23 − r33 ). This, together with 0 < (r13 − r11 ) < (r13 − r23 ),
0 < (r34 − r13 ) < (r34 − r24 )
0 < (r24 − r14 ) < (r11 − r14 ),
0 < (r23 − r33 ) < (r13 − r33 ),
implies that
(r13 − r11 )(r43 − r13 ) < (r11 − r14 )(r13 − r33 ). From Proposition 1.1, we have pˆ 1 = (r33 − r13 )/∆13 and pˆ 3 = (r11 − r13 )/∆13 . If we substitute these values into (12), then the above inequality and ∆13 < 0 imply that λ413 < 0, which means that P13 is asymptotically stable in ∆, a contradiction. The proof of case (ii) is complete. (iii) V1 is asymptotically stable in {13} but unstable in {14} (r13 < r11 and r14 > r11 ). The proof of this case is similar to case (ii), and will be omitted. The proof of Rule 6 is complete.
Contents lists available at SciVerse ScienceDirect
Theoretical Population Biology journal homepage: www.elsevier.com/locate/tpb
Patterns for four-allele population genetics model Linlin Su, Roger Lui ∗ Department of Mathematical Sciences, Worcester Polytechnic Institute, 100 Institute Road, Worcester, MA 01609, United States
article
info
Article history: Received 1 September 2011 Available online 26 February 2012 Keywords: k-allele model Asymptotically stable equilibrium Pattern Eigenvalue Computer simulation
abstract In this paper, we find and classify all existing patterns for a single-locus four-allele population genetics models in continuous time. An existing pattern for a k-allele model means a set of all coexisting asymptotically stable equilibria with respect to the flow defined by the system of equations p˙ i = pi (ri − r ), i = 1, . . . , k, where pi and ri are the frequency and marginal fitness of allele Ai , respectively, and r is the mean fitness of the population. It is well known that for the two-allele model there are only three existing patterns, depending on the relative fitness between the homozygotes and the heterozygote. For the three-allele model there are 14 existing patterns, and we shall show in this paper that for the fourallele model there are 117 existing patterns. We also describe the domains of attraction for coexisting asymptotically stable equilibria. The problem of finding existing patterns has been studied in the past, and it is an important problem because the results can be used to predict the long-term genetic makeup of a population. It should be pointed out that this continuous-time model is only an approximation to the corresponding discrete-time model. However, the set of equilibria and their stability properties are the same for the two models. © 2012 Elsevier Inc. All rights reserved.
1. Introduction Consider a diploid population whose members possess a gene that occurs in k different forms, called alleles, located at an autosomal locus. Then there are k(k + 1)/2 genotypes denoted by Ai Aj , i, j = 1, . . . , k. Suppose that Ai Aj has constant fitness denoted by rij . Let R = (rij ) be the (symmetric) fitness matrix, and let pi be the frequency of allele Ai . Then, assuming random mating and discrete-time non-overlapping generations, the frequency of allele Ai in the next generation is given by p′i = pi ri /r ,
i = 1, . . . , k,
(1)
where ri = j rij pj is the marginal fitness of allele Ai and r = i ri pi is the mean fitness of the population. For continuous-time overlapping generations, the above equations are replaced with p˙ i = pi (ri − r ),
i = 1, . . . , k.
(2)
In this paper, we focus our study on continuous-time models because we can apply the theory of dynamical systems. Although system (2) is only an approximation to (1), derived under the assumption of Hardy–Weinberg equilibrium (Bürger, 2000, Section 10.1), the set of equilibria and their stability properties are the same for the two models.
∗
Corresponding author. E-mail addresses: [email protected] (L. Su), [email protected] (R. Lui).
0040-5809/$ – see front matter © 2012 Elsevier Inc. All rights reserved. doi:10.1016/j.tpb.2012.02.002
Since the pi represent gene frequencies, they add up to 1. In fact, it can be shown that, if initially i pi (0) = 1, then the same holds for the solutions pi (t ), i = 1, . . . , k, of system (2) for t > 0. Therefore, (2) is well posed, and is a system of k − 1 equations with pℓ , 1 ≤ ℓ ≤ k, replaced with 1 − j̸=ℓ pj . By relabeling the alleles we may assume, without loss of generality, that the rii are nonincreasing. To avoid degeneracy, we assume throughout this paper that r11 > r22 > · · · > rkk
(3)
and that rij , where i ̸= j, is not equal to rii or rjj . An equilibrium of system (2) is aconstant vector p = (p1 , . . . , pk ) that satisfies the conditions i pi = 1 and pi (ri − r ) = 0, i = 1, . . . , k. An equilibrium is called an interior equilibrium if pi > 0 for all i and a boundary equilibrium if pi = 0 for at least one i. Note that Vi = (0, . . . , 1, . . . , 0), where the 1 is in the ith component, is always a boundary equilibrium. Let us consider the two-allele case. Since p1 + p2 = 1, system (2) may be reduced to p˙ 1 = p1 (1 − p1 )(∆12 p1 + r12 − r22 ),
(4)
where ∆12 = r11 + r22 − 2r12 . The following result is well known (Nagylaki, 1992, pp. 52–55). Proposition 1.1. For the two-allele model (4) with r11 > r22 , p1 = 0 and p1 = 1 are both equilibria, and a third interior equilibrium p∗ = (r22 − r12 )/∆12 exists if and only if r12 > r11 or r12 < r22 . There are three cases: (i) r22 < r12 < r11 : 1 is globally asymptotically stable, 0 is unstable and p∗ does not exist; (ii) r12 > r11 : 0, 1 are both
274
L. Su, R. Lui / Theoretical Population Biology 81 (2012) 273–283
unstable and p∗ exists and is globally asymptotically stable; (iii) r12 < r22 : 0, 1 are both asymptotically stable and p∗ exists and is unstable. We call these three cases the (heterozygote) intermediate, superior, and inferior case, respectively. The above proposition implies that for the two-allele model there are only three existing patterns: ⟨{1}⟩, ⟨{p∗ }⟩, and ⟨{0}, {1}⟩. An existing pattern is simply the list of all coexisting asymptotically stable equilibria. For example, the pattern ⟨{p∗ }⟩ means that the interior equilibrium p∗ is the only asymptotically stable equilibrium (and, consequently, is globally asymptotically stable) and, from Proposition 1.1, this happens in the superior case. The pattern ⟨{0}, {1}⟩ means there are two coexisting asymptotically stable equilibria, p1 = 0 and p1 = 1, and this happens in the inferior case. On the other hand, the pattern ⟨{p∗ }, {1}⟩ does not exist regardless of the choice of R. For the k-allele model there are 2k − 1 possible equilibria, but the number of coexisting asymptotically stable equilibria is considerably less. For example, for the four-allele model, there are 15 possible equilibria and potentially 215 − 1 patterns, but only 117 of them can exist, and none can contain more than four asymptotically stable equilibria. The problem of finding all patterns for a k-allele model has been studied extensively before. In a series of papers, (Broom et al., 1993; Cannings and Vickers, 1988, 1989, 1991; Vickers and Cannings, 1988), the authors studied maximal existing patterns (i.e., existing patterns which are not a subset of another existing pattern) for k = 3, 4, 5. Cannings and Vickers conjectured that any subset of an existing pattern is also an existing pattern (Cannings and Vickers, 1989; Vickers and Cannings, 1988). Thus, if the conjecture is true, one can specify the complete set of existing patterns by specifying the maximal existing patterns. This is the approach adopted in the literature. The main purpose of this paper is to give a complete classification of all existing patterns for the four-allele model without assuming the truth of the conjecture of Cannings and Vickers. The complete list of existing patterns for the case k = 3 is given by the authors in Su et al. (2011). We should also point out that condition (3) is not posited in the literature. Therefore, patterns that are rotations of each other, which can be obtained by multiplying the fitness matrix R on both sides by a permutation matrix, are considered as one pattern. However, under condition (3), rotation of an existing pattern is not necessarily an existing pattern. Therefore, we treat patterns that are rotations of each other as different patterns. For these reasons, the number of existing patterns we obtained is greater than those obtained in Cannings and Vickers (1989) and Vickers and Cannings (1988). The method we use in this paper, counting signs of eigenvalues of the Jacobian matrix evaluated at the existing equilibria, is also completely different from the methods used in Cannings and Vickers (1989) and Vickers and Cannings (1988). This paper is organized as follows. In Section 2, we review several important properties of system (2). Some of the results in this section may be new, since we cannot find them in the literature. In Section 3, we summarize our results for the threeallele model in Su et al. (2011). In Section 4, we present our results for the four-allele model. In Section 5, we describe the domains of attraction for some of the patterns. Section 6 is discussions. In order not to interrupt the flow of ideas, many tables and some of the proofs are given in the Appendix.
Proposition 2.1. Assume that det R ̸= 0. An interior equilibrium exists if and only if all the Ui are non-zero and have the same sign. Moreover, if an interior equilibrium exists, it is uniquely given by Ui pi = , Uj
i = 1, . . . , k.
j
The degenerate case det R = 0 was considered by an der Heiden (1975) and Hughes and Seneta (1975). In this paper, we always assume that the determinants of R and all its submatrices are nonzero. A submatrix of R is a matrix obtained by erasing any number of rows and the corresponding columns of R. Next, we consider the stability of existing equilibria of system (2). It is well known that the mean fitness function r is strictly increasing along non-constant solutions of (2); see Kingman (1961b), Mulholland and Smith (1959), and Scheuer and Mandel (1959). Therefore, r plays the role of a Lyapunov function, and asymptotically stable interior equilibria are identified with the local maxima of r. Quite a few necessary and sufficient conditions for the stability of an existing interior equilibrium can be obtained by studying the properties of local maxima of r; see Hughes and Seneta (1975), Mandel (1970) and the references therein. Another more direct approach is to investigate the eigenvalues of the Jacobian matrix J evaluated at an existing equilibrium of system (2). The difficulty with this approach is the complexity of finding the relationships between the fitness matrix R = (rij ), the existence of an equilibrium, and the signs of the eigenvalues of J evaluated at an equilibrium in a k-allele model. Besides finding and classifying all patterns for the four-allele model, another main contribution of this paper is to show readers how to overcome this difficulty with the help of computer simulations. Let p∗ = (p∗1 , . . . , p∗k ) be an equilibrium. We write K for the set {1, . . . , k}, and define L = {i ∈ K : p∗i > 0} and L¯ = {i ∈ K : p∗i = 0}. The following proposition may be known, but we are unable to find a reference for it, so we present a proof here.
Proposition 2.2. Let l ∈ K , and replace pl with 1 − i̸=l pi in (2). Then the ith row, jth column entry of the Jacobian matrix J of (2) (as a system of (k − 1) equations) evaluated at an equilibrium point p∗ is given by Jij = δij (rˆi − r )|p∗ + pˆ∗ [(rˆiˆj − rˆil ) − 2(rˆj − rl )|p∗ ], i i, j = 1, . . . , k − 1, where ˆi = i if 1 ≤ i < l and ˆi = i + 1 if i ≥ l; namely,
Jij =
pˆ∗ [(rˆiˆj − rˆil ) − 2(rˆj − rl )|p∗ ]
if ˆi ∈ L,
δij (rˆi − r )|p∗
if ˆi ∈ L¯ .
i
In particular, if p∗ is an interior equilibrium, then Jij = pˆ∗ (rˆiˆj − rˆil ), i
i, j = 1, . . . , k − 1.
Proof. Jij = ∂∂p [pˆi (rˆi − r )]|p∗ = δˆiˆj (rˆi − r )|p∗ + pˆ∗ ∂∂p (rˆi − r )|p∗ . i ˆj ˆj
On the other hand, r =
Now rˆi =
rpˆ = j
ℓ̸=l
ℓ̸=l rˆiℓ pℓ + rˆil
1−
=
ℓ̸=l
ℓ̸=l
pℓ . Therefore, (rˆi )pˆ = rˆiˆj − rˆil . j
(rℓˆj − rℓl )pℓ + rˆj +
k
ℓ=1
pl − rl
rlℓ pℓ pˆ j
= rˆj − rl − rˆjl pl + rll pl + rˆj + rlˆj pl − rll pl − rl = 2(rˆj − rl ).
ℓ̸=l pℓ . Therefore,
(rℓ )pˆj pℓ + rˆj + (rl )pˆj pl − rl
ℓ̸=l rℓ pℓ + rl 1 −
2. Theory Results in this section are valid for any number of alleles. We begin with the problem of the existence of interior equilibria. Let Rij be the cofactor of rij in the fitness matrix R, and let Ui = j R ij . The following result is due to Mandel (1959).
(5)
L. Su, R. Lui / Theoretical Population Biology 81 (2012) 273–283
Noting that ri = r for i ∈ L, and at an interior equilibrium that r1 = · · · = rk , the conclusion of the proposition follows. Remark 2.3. If we choose l to be in L, then
Jij =
pˆ∗ (rˆiˆj − rˆil )
if ˆi, ˆj ∈ L,
δij (rˆi − r )|p∗
if ˆi ∈ L¯ .
i
This implies that the eigenvalues of J consist of (ri − r )|p∗ (i ∈ L¯ ) and the eigenvalues of the matrix (p∗i (rij − ril ))i,j∈L−{l} . This formula greatly reduces the amount of work when we compute the eigenvalues of J. For example, for a five-allele model, if we want to find the eigenvalues of J at the boundary equilibrium on the side where p2 = 0, p5 = 0, all we have to do is find the interior equilibrium of the three-allele model with fitness matrix obtained by eliminating the second and fifth rows and columns of R. Let this interior equilibrium be denoted by (¯p1 , p¯ 2 , p¯ 3 ), and let p∗ = (¯p1 , 0, p¯ 2 , p¯ 3 , 0). Then, according to the above lemma with l = 4, the eigenvalues of J evaluated at p∗ consists of the eigenvalues of the matrix p¯ 1 (r11 − r14 ) p¯ 2 (r31 − r34 )
p¯ 1 (r13 − r14 ) p¯ 2 (r33 − r34 )
Remark 2.4. The eigenvalues of J evaluated at the boundary equilibrium p∗ = (q∗1 , . . . , q∗i−1 , 0, q∗i , . . . , q∗k−1 ) of a k-allele model with fitness matrix R consist of the k − 2 eigenvalues of the Jacobian matrix evaluated at q∗ = (q∗1 , . . . , q∗k−1 ) and ri · p∗ − R p∗ · p∗ . It turns out that this last eigenvalue is related to the interior equilibrium of the k-allele model by the formula ri · p∗ − R p∗ · p∗ = −
Ui D
,
(6)
where D is the sum of all the entries of the adjugate matrix of R1 , R1 being the matrix obtained from erasing the ith row and column of R. We do not use (6) in this paper. Let
= p∗i (rij − 2r )|p∗ . Let Γ = (Γij ) ≡ (rij − 2r ) p∗i p∗j , and let M = diag(1/ p∗1 , . . . , 1/ p∗k ). It is clear that J˜ = M −1 Γ M. Since Γ is symmetric, all the eigenvalues of J˜ are real. Finally, we observe that, if λ is an eigenvalue of J with eigenvector v = (v1 , . . . , vk−1 ), then λ is also an eigenvalue of J˜ with eigenvector v˜ = (v1 , . . . , vk−1 , −v1 −· · ·− vk−1 ). To see this, let Fi (p1 , . . . , pk ) = pi (ri − r ) for i = 1, . . . , k. Then Jij =
∂ Fi ∂ Fi − ∂ pj ∂ pk
for i, j = 1, . . . , k − 1
and k−1 ∂ Fi
∂ pj
−
∂ Fi ∂ pk
vj = λvi for i = 1, . . . , k − 1.
Therefore, for i = 1, . . . , k − 1, we have k−1
J˜ij vj + J˜ik (−v1 − · · · − vk−1 ) =
j =1
k−1 (J˜ij − J˜ik )vj j =1
=
k−1 ∂ Fi j =1
∂ pj
−
∂ Fi ∂ pk
vj = λvi .
Now, consider the function G(p1 , . . . , pk ) = i=1 Fi (p1 , . . . , pk ) defined in Rk . The set ∆ may be thought of as the plane determined by p1 + p2 + · · · + pk − 1 = 0 in the positive quadrant of Rk . Since G = 0 on ∆, ∇ G is parallel to the normal of ∆, which is the k-dimensional vector e = ⟨1, . . . , 1⟩. This implies that ∂ G/∂ p1 = · · · = ∂ G/∂ pk . From above, for i = k, we have
k
k−1
J˜kj vj + J˜kk (−v1 − · · · − vk−1 )
j =1
k−1 k−1 ∂ Fk ∂ Fk (J˜kj − J˜kk )vj = − vj ∂ pj ∂ pk j =1 j=1 k−1 k−1 k−1 ∂G ∂ Fl ∂G ∂ Fl = − − − vj ∂ pj ∂ pj ∂ pk ∂ pk j =1 l=1 l =1 k−1 k−1 ∂ Fl ∂ Fl =− − vj = λ(−v1 − · · · − vk−1 ). ∂ pj ∂ pk l=1 j=1
=
∆=
∂ ∂ [pi (ri − r )]|p∗ = δij (ri − r )|p∗ + p∗i (ri − r )|p∗ ∂ pj ∂ pj ∂ ri ∂ r ∗ = pi = p∗i (rij − 2rj |p∗ ) − ∂ pj p∗ ∂ p j p∗
J˜ij =
j =1
and ri · p∗ − R p∗ · p∗ , i = 2, 5, where ri is the ith row of the fitness matrix R. This decomposition was also used in Kingman (1961a) to study the stability of boundary equilibria for the discrete-time model.
275
(p1 , . . . , pk−1 ) : 0 ≤ pi ≤ 1 for i = 1, . . . , k − 1, k−1
pi ≤ 1 .
(7)
i=1
It can be shown that ∆ is invariant with respect to the flow defined by (2), and, given initial conditions in ∆, the solutions of (2) converge to an asymptotically stable equilibrium in ∆ as t → ∞ (Bürger, 2000, p. 38). The set ∆ contains a maximum possible 2k − 1 equilibria. We now prove that the eigenvalues of the Jacobian matrix J evaluated at any of these equilibria are real. The discrete-time version of this theorem may be found in Kingman (1961a, Section 3). Proposition 2.5. All the eigenvalues of the Jacobian matrix J at an equilibrium are real. Proof. Because of Proposition 2.2, it suffices to prove the result at an interior equilibrium. Let J be the Jacobian matrix of (2) with pk ˜ replaced with 1 − i̸=k pi , and let J be the Jacobian matrix of (2) as a system of k equations. At an interior equilibrium, all the rj are equal to r. Therefore,
˜ Thus, all the eigenvalues This shows that λ is also an eigenvalue of J. of J are real, and the proof of the proposition is complete. Proposition 2.6. The determinant of J evaluated at an interior equilibrium p∗ = (p∗1 , . . . , p∗k ) is given by det(J ) = p∗1 p∗2 · · · p∗k (U1 + U2 + · · · + Uk ).
(8)
Proof. Let ri , i = 1, . . . , k − 1, be the ith row of the matrix R without the last element rik , and let e be the row vector with (k − 1) 1s. From the formula for J at the interior equilibrium p∗ in Proposition 2.2, we have
det(J ) = p∗1 p∗2 · · · p∗k−1 det
r1 − r1k e r2 − r2k e
. ··· r(k−1) − r(k−1)k e
276
L. Su, R. Lui / Theoretical Population Biology 81 (2012) 273–283
Rule 1. For any k, if the interior equilibrium is asymptotically stable, then no other equilibrium is asymptotically stable.
To continue,
det
r1 − r1k e r2 − r2k e
Rule 2. The maximum number of coexisting asymptotically stable equilibria in a three-allele model is 3, which must be the pattern consisting of the three vertices.
··· r(k−1) − r(k−1)k e
r1 r2 − r2k e
e
Rule 3. In a three-allele model, if all the existing asymptotically stable equilibria are vertices, then they must include V1 .
r2 − r2k e − r1k det ··· ··· r(k−1) − r(k−1)k e r(k−1) − r(k−1)k e
= det
e
r1 r2 − r2k e
= det
r2 − r1k det · · ·
···
= det
e r2
r1 r2
4. Four-allele model
r(k−1)
r(k−1) − r(k−1)k e
Rule 4. In a three-allele model, if vertex V2 or V3 is asymptotically stable but not V1 , then P2 or P3 must also be asymptotically stable, respectively.
r1 r2
− · · · − r(k−1)k det · · · − r1k det ··· ··· r(k−1)
r(k−1) r1 r2
. = det .. r
(k−1)
e
rk−2 e
r1k r2k
.. .
r(k−1)k 1
= Uk = p∗ (U1 + U2 + · · · + Uk ). k
The second to last equality is due to the definition of Ui given before Proposition 2.1, and the last equality follows from (5). Thus, formula (8) follows. The following two results are helpful when we show the nonexistence of certain patterns. Theorem 2.7 (Kingman (1961a)). Suppose that ∆1 ⊂ ∆2 are two distinct invariant subsets of ∆ and that the corresponding interior equilibria exist; then they both cannot be asymptotically stable. In particular, if ∆ contains an asymptotically stable interior equilibrium, then no other equilibrium in ∆ can be asymptotically stable, and, therefore, the interior equilibrium is globally asymptotically stable. Theorem 2.8 (Kingman (1961a) and Mandel (1959)). If the interior equilibrium p∗ is globally asymptotically stable, then the mean fitness r has a global maximum there. Consequently, r |p∗ ≥ rii for i = 1, . . . , k, which is obtained by evaluating r at the vertex (0, . . . , 0, 1, 0, . . . , 0). 3. Three-allele model For the three-allele case, (2) can be written as p˙ i = pi (ri − r ),
i = 1, 2.
(9)
The invariant set ∆ defined by (7) is the triangle joining the vertices V1 = (1, 0), V2 = (0, 1) and V3 = (0, 0) in the (p1 , p2 )plane. There are seven possible equilibria in ∆, comprising three vertices, three boundary equilibria on the edges, and one interior equilibrium. The edge Si = {pi = 0} represents a two-allele case, and the interior equilibrium on that edge is denoted by Pi . In Su et al. (2011), we showed that there are 14 existing patterns, but none of them contains more than three asymptotically stable equilibria. We discovered four rules for the three-allele patterns, and they are summarized below. These rules will be used in the study of the four-allele case. Rule 1 follows from Theorem 2.7, and is valid for any number of alleles. Three-allele patterns are also studied in Mandel (1959) and Wright (1969, p. 46) without assuming condition (3).
For a higher number of alleles, there are too many possible patterns, and we cannot prove their existence as we did for the three-allele model (see Su et al. (2011)). We therefore resort to computer simulations to find examples of existing patterns and, with the help of the above four rules, prove that the patterns that do not show up in our computer simulations cannot exist. For our computer simulations, we wrote a Matlab program that randomly generates a fitness matrix R satisfying condition (3), and we use the method in Remark 2.3 to compute and count the number of negative eigenvalues of J evaluated at each existing equilibrium. We then write the results to a text file and run a sorting algorithm on it to display all the distinct patterns. We repeat this process several times with increasing number of simulations until no new patterns show up. Examples of patterns containing two edge equilibria and one face equilibrium are the most difficult to generate, and they require a more targeted approach which will not be described here. We introduce some notation. The invariant set ∆ defined by (7) is a tetrahedron. It contains 15 possible equilibria, comprising four vertices, six edge equilibria, four face equilibria, and one interior equilibrium. Each edge of ∆ is a two-allele model and each face of ∆ is a three-allele model. We denote by v a vertex equilibrium, e an edge equilibrium, f a face equilibrium, and i the interior equilibrium. A pattern containing α vertices, β edge equilibria, and γ face equilibria is denoted by αv + β e + γ f , where α, β, γ are nonnegative integers. We say that three edge equilibria are coplanar if the edges they lie on form the sides of a triangle. From Rule 2, coplanar edge equilibria cannot all be asymptotically stable. Also, we shall prove that the maximum number of coexisting asymptotically stable equilibria is 4, which we state as Rule 5 below. From Rule 1, if the interior equilibrium is asymptotically stable, then no other equilibrium can be asymptotically stable. Based on these properties, we list all the possible patterns for the four-allele model in Table 1. The existing ones are confirmed by examples from simulations (see Table 2 in Appendix A for details) and the rest are shown not to exist below. In each line of Table 1, Ns denotes the number of asymptotically stable equilibria and Np denotes the number of patterns of that type. If we add up the numbers in the column Np , we see that there are 117 existing patterns for the four-allele model. We also see from Table 1 that there are only two types of existing pattern with four asymptotically stable equilibria, one of which is the pattern containing all four vertices and the other is the pattern containing four asymptotically stable non-coplanar edge equilibria. Proof of Table 1. Appendix A gives examples of fitness matrices R obtained from simulations which produce the existing patterns in Table 1. Rule 1 is the most useful in ruling out the existence of a certain pattern. For example, the pattern 2v + f is not possible because, if two vertices are asymptotically stable, then any face in ∆ must contain at least one of the two asymptotically stable vertices, and this violates Rule 1 for three-allele case. Also, Rule 1 for
L. Su, R. Lui / Theoretical Population Biology 81 (2012) 273–283
277
Table 1 Four-allele model: Total 117 existing patterns. Ns
Pattern type
v
Np e
f
1 1
1
1 1 2
1 1
15 2
6
1 1
12 4 12
1 1
3
3 3 3
3
2 2 1
1 6 4 1 3
2 2
Explanations
i
1 1 2
must be (1, 0, 0), Rule 6 any of the six edges any of the four faces the interior equilibrium one of the two vertices must be (1, 0, 0), Rule 6
any two edges: any two faces:
1 1
1 1
1 2 2 1
4 4 4 3 3 1 4
1 1 3 3
1 2 2 2 1 1
1 2 2 1 2 1
1 3 3 2 2 1 1 2
the four-allele case states that, if the interior equilibrium of the tetrahedron is asymptotically stable, then no other equilibrium can be asymptotically stable. Therefore, we only list one pattern in Table 1 with the interior equilibrium i. The patterns 3f , e + 2f , 3e + f , and 4f were shown to be not possible in Vickers and Cannings (1988, Theorems 4 and 6), and we call this reason A in the explanation column of Table 1. We also use Rule 5 and Rule 6 below in establishing Table 1. We prove Rule 5 here and postpone the proof of Rule 6 to Appendix B.
4 2
= 15
=6
any edge with a vertex not in the edge: 6 × 2 = 12 pair any face with the vertex not in the face any face with an edge not in the face: 4 × 3 = 12 one of the three vertices must be (1, 0, 0), Rule 6 6 3
− 4 = 16
16
any three edges not coplanar:
0 6 0 12
does not exist because of reason A pair any edge with the two vertices not in the edge does not exist because of Rule 1 anyvertex with two edges not containing the vertex: 4 1
2
6
2
12 0 0 0 1 3 0 0 0 0 0 0 0 0 0 0 0 0 0
×
3 2
= 12
any face with two of the three edges not from the face does not exist because of Rule 1 does not exist because of reason A does not exist because of Rule 1 all four vertices any four edges without forming a triangle does not exist because of reason A does not exist because of Rule 1 does not exist because of Rule 1 does not exist because of Rule 1 and Rule 2 does not exist because of reason A does not exist because of Rule 1 does not exist because of Rule 1 does not exist because of Rule 1 does not exist because of Rule 1 does not exist because of Rule 1 does not exist because of Rule 1 does not exist because of Rule 1 does not exist because of Rule 1
we will violate the rule mentioned above that there cannot be coplanar asymptotically stable edge equilibria. The proof of Rule 5 is complete. Remark 4.1. From Table 1 we observe that, if ∆ contains three or more asymptotically stable equilibria, then none can be a face equilibrium except for the pattern 2e + f , where the two asymptotically stable edge equilibria do not lie on the face. This observation is useful in the study of the five-allele model.
Rule 5. The maximum number of coexisting asymptotically stable equilibria for the four-allele model is 4.
5. Domains of attraction
Rule 6. In a four-allele model, if all the asymptotically stable equilibria are vertices, then they must include the vertex (1, 0, 0). To prove Rule 5, first recall that ∆ has four vertices, six edges, and four faces. Suppose we have a pattern with five or more asymptotically stable equilibria. Then, by Rule 1, the pattern must not contain the interior equilibrium of ∆. Suppose that there is an asymptotically stable vertex; then, from Rule 1, all the other asymptotically stable equilibria are in the face opposite the vertex. From Rule 2, there can be at most three asymptotically stable equilibria in this face, so this is impossible because we need four more asymptotically stable equilibria. Suppose that a face equilibrium is asymptotically stable; then, according to Rule 1, the other asymptotically stable equilibria must be either be the vertex opposite the face, or the three edge equilibria not in the face. For the same reason, this is again impossible. Finally, if all the asymptotically stable equilibria are on the edges of ∆, then
Knowing all the patterns is not enough to predict the asymptotic behavior of solutions of (2). Given initial data, we also want to know which asymptotically stable equilibrium the solutions of (2) approach as t → ∞. The domain of attraction of an asymptotically stable equilibrium is the subset of ∆ such that solutions of (2) with initial data in this subset converge to the asymptotically stable equilibrium as t → ∞. For the three-allele model, the domains of attraction are separated by invariant curves, called separatrices, which are not analytically tractable. They are described in Theorem 5.5 of Su et al. (2011). For the four-allele model, the domains of attraction are formed by partitioning the tetrahedron ∆ by surfaces which are two-dimensional invariant manifolds. Each domain of attraction contains one asymptotically stable equilibrium. The intersection of the invariant manifolds with the four faces of ∆ are separatrices in the corresponding
278
L. Su, R. Lui / Theoretical Population Biology 81 (2012) 273–283
Fig. 1. Domains of attraction for the pattern 2v + e. M1 separates out the DA of v1 . On M1 , f2 is asymptotically stable. M2 separates out the DA of v2 . On M2 , f1 is asymptotically stable. The intersections of M1 and M2 are the separatrices connecting f3 , e4 and e4 , f4 , respectively. The remaining region in ∆ is the DA of e3 .
Fig. 2. Domains of attraction for the pattern v + 2e. M1 separates out the DA of v1 . On M1 , f3 is asymptotically stable. M2 separates out the DA of e6 . On M2 , f1 is asymptotically stable. The intersection of M1 and M2 are the separatrices connecting e4 , f4 , e5 and e5 , f2 , e3 , respectively. The remaining region in ∆ is the DA of e2 .
Fig. 3. Domains of attraction for the pattern 2e + f . M1 separates out the DA of f4 which is on the face {123}. On M1 , f3 is asymptotically stable. M2 separates out the DA of e3 . On M2 , f1 is asymptotically stable. The intersection of M1 and M2 is the separatrix connecting e1 , f2 , v3 . The remaining region in ∆ is the DA of e2 .
three-allele case. These separatrices may be identified by looking at Theorem 5.5 of Su et al. (2011) or by simulations. For each asymptotically stable equilibrium, we look at the separatrices that separate out the equilibrium on the faces where the equilibrium lies. The invariant manifolds are then put together by connecting these separatrices. In the following, we describe the domains of attraction for some of the existing patterns in the four-allele model with three or four asymptotically stable equilibria. If there is only one asymptotically stable equilibrium, the domain of attraction contains the interior of ∆, and, if there are two asymptotically stable equilibria, the domains of attraction are separated by an invariant manifold. These cases are easy to visualize, and they will not be described here. From Table 2, there are seven groups of patterns with three or four stable equilibria: 2v + e, v + 2e, 2e + f , 3v , 3e, 4v , and 4e. Because of space limitations, we only describe the domains of attraction for the first line of each group given in Table 2 (Figs. 1–7). For example, Fig. 1 below shows the domains of attraction for the pattern 2v + e with the fitness matrix r11 = 1, r12 = −4.544,
r13 = −0.714, r14 = −4.035, r22 = 0.065, r23 = −4.021, r24 = −2.943, r33 = 0.039, r34 = 4.420, and r44 = 0, which may be found in the first row of the group 2v + e. In the figures, the equilibria are denoted by v1 = {1}, v2 = {2}, v3 = {3}, v4 = {4}, e1 = {14}, e2 = {24}, e3 = {34}, e4 = {12}, e5 = {13}, e6 = {23}, f1 = {234}, f2 = {134}, f3 = {124}, f4 = {123}, and i = {1234}, which is the interior equilibrium. In Figs. 1–7, the asymptotically stable equilibria are indicated by dots and the unstable ones by open circles, except for unstable vertices. Solid curves represent separatrices on the face facing the reader and dotted curves represent separatrices on the other three faces of ∆. Note that an equilibrium which is unstable in ∆ may be stable on a face as an equilibrium in the three-allele case. This lowerdimensional stability is not reflected in the figures. The invariant manifolds together with the dynamics on them are also given. The figures are only schematic, because the invariant manifolds are in general analytically intractable. We use DA to denote a domain of attraction in the captions of the figures.
L. Su, R. Lui / Theoretical Population Biology 81 (2012) 273–283
279
Fig. 4. Domains of attraction for the pattern 3v . M1 separates out the DA of v3 . On M1 , e5 is asymptotically stable. M2 separates out the DA of v2 . On M2 , e2 is asymptotically stable. The remaining region in ∆ is the DA of v1 .
Fig. 5. Domains of attraction for the pattern 3e. M1 separates out the DA of e4 . On M1 , f4 is asymptotically stable. M2 separates out the DA of e1 . On M2 , f2 is asymptotically stable. The remaining region in ∆ is the DA of e5 .
Fig. 6. Domains of attraction for the pattern 4v . M1 separates out the DA of v3 . On M1 , e3 and e6 are asymptotically stable and their DA are separated by a separatrix connecting f2 and f1 . M2 separates out the DA of v1 . On M2 , e1 and e4 are asymptotically stable and their DA are separated by a separatrix connecting f2 and f3 .M3 separates out the DA of v2 and v4 from the remaining region of ∆. On M3 , e2 is asymptotically stable.
When we construct the domains of attraction, sometimes the choice of invariant manifolds is not obvious. This can be resolved by looking at the stable and unstable manifolds of the face equilibria. For example, in Fig. 5, suppose that M1 and M3 are the invariant manifolds rather than M1 and M2 shown in the figure, where M3 contains v1 , f4 , e6 , f1 , e3 , and f2 on its boundary. Then look at the face equilibrium f4 , which has two negative eigenvalues and one positive eigenvalue. The unstable direction corresponding to the positive eigenvalue lies on the face {123}. Furthermore, f4 is stable along the intersection of M1 and M3 , which is the separatrix on the face {123}, and f4 is stable on M1 , because no other equilibrium is stable on M1 . This implies that f4 must have an unstable direction on M3 , which is a contradiction, because it already has an unstable direction on the face {123}. One can also rule out the choice of M2 and M3 as the invariant manifolds by looking at the face equilibrium f2 . Therefore, in Fig. 5, M1 and M2 are the correct invariant manifolds which give the domains of attraction in ∆ (Figs. 6 and 7).
6. Discussions In this paper, we give all existing patterns for the four-allele selection model under conditions (3) and rij ̸= rii or rjj for i ̸= j. If, for some i, j, rij = rii or rjj , then the interior equilibrium on the edge {ij} disappears, and this is called complete dominance. Complete dominance reduces the total number of existing patterns in a k-allele model. Knowing all existing patterns is important, because it tells us the possible long-term genetic makeup of the population, that is, which alleles will survive and which will disappear under selection forces in the long run. The actual outcome will depend on the initial conditions of the population if there are multiple coexisting asymptotically stable equilibria. In such a case, the set ∆ is divided into domains of attraction, each containing one asymptotically stable equilibrium. The problem of finding all existing patterns for the k-allele model has been studied before, but it remains an open problem except when k = 2, 3, or 4. Our results are different from earlier work because we do not assume the truth of the conjecture of
280
L. Su, R. Lui / Theoretical Population Biology 81 (2012) 273–283
Fig. 7. Domains of attraction for the pattern 4e. M1 and M2 intersect at the invariant curve passing through e1 , i, and e6 . They divide ∆ into four subregions that correspond to the DA of e2 , e3 , e5 , and e4 , respectively. On M1 , f1 and f4 are asymptotically stable. On M2 , f2 and f3 are asymptotically stable.
Cannings and Vickers, that any subset of an existing pattern is also an existing pattern. In fact, this conjecture is not true in our case because, for example, the pattern ⟨{1}, {2}⟩ exists but not the pattern ⟨{2}⟩. Also, unlike earlier work, rotation of an existing pattern may not be an existing pattern. For example, ⟨{1}, {2}⟩ is an existing pattern but ⟨{2}, {3}⟩ is not. All these are consequences of condition (3), which may be satisfied by relabeling the alleles assuming that the rii are different from each other. In this paper, we study the continuous-time model instead of the discrete-time model (1), which allows us to adopt a completely different approach of counting signs of the eigenvalues of the Jacobian matrix evaluated at all existing equilibria. This approach motivates us to prove several interesting mathematical results, which are valid for any number of alleles. We believe that some of these results are new, because we could not find them in the literature. We have also made considerable progress towards solving the patterns problem for the five-allele model in Su et al. (2011). The number of existing patterns is around 2351. For the six-allele model, there will be more than 60,000 existing patterns, and we conjecture that the number of existing patterns grows exponentially as a function of k. It is worth pointing out that the two-locus two-allele selection model degenerates to the single-locus four-allele model in the absence of genetic recombination (Bürger, 2000, p. 47). Therefore, the four-allele model studied in this paper has a special role in population genetics as the reference model for the two-locus twoallele model. Acknowledgment We thank Allan E. Johannesen of the WPI Computing and Communications Center for writing the sorting algorithm. Appendix A. Examples of existing patterns In the following tables, we denote the vertices (1, 0, 0), (0, 1, 0), (0, 0, 1), (0, 0, 0) by {1}, {2}, {3}, {4}, respectively. {ij} denotes the interior equilibrium of the edge joining vertices {i} and {j}, and {ijk} denotes the interior equilibrium of the face with vertices {i}, {j}, and {k}. Note that, by scaling time, system (2) is unchanged under a linear transformation of the rij , so we may assume that r11 = 1 and r44 = 0. Appendix B. Proof of Rule 6 We first introduce the following notation: V1 = (1, 0, 0), V2 = (0, 1, 0), V3 = (0, 0, 1), V4 = (0, 0, 0), {ij} denotes the edge joining Vi and Vj , Pij denotes the edge equilibrium in {ij}, {ijk} denotes the face with vertices Vi , Vj , and Vk , and Pijk denotes the face equilibrium in {ijk}. It is easy to verify that the eigenvalues of
(j)
J evaluated at Vi are λi = rji − rii , j = 1, 2, 3, 4, j ̸= i. The matrix J evaluated at an edge equilibrium, say P13 , has three eigenvalues, one of which may be determined from applying Proposition 1.1, namely, P13 is asymptotically stable in {13} if r13 > r11 and unstable in {13} if r13 < r33 . Let the other two eigenvalues be denoted by λ213 and λ413 , where the superscript denotes that the vertex P13 is facing in the face where {13} is an edge. In this case, {13} lies on the faces {123} and {134}, so the vertices P13 is facing are V2 and V4 , respectively. We shall prove Rule 6 by showing that, if Vi , i ̸= 1, is asymptotically stable, then V1 must also be asymptotically stable. We only prove the case i = 2, since it is the most difficult. In the proof of the following lemma, we shall use Remark 2.3 repeatedly without mentioning it. Lemma B.1. Suppose that V2 is asymptotically stable in ∆, and that there is no non-vertex asymptotically stable equilibrium in ∆; then V1 must be asymptotically stable in ∆. Proof. Since V2 is asymptotically stable and r11 > r22 , {12} is in the inferior case, and hence V1 is asymptotically stable in {12}. We now focus on the invariant triangle {134}, which must contain at least one equilibrium asymptotically stable in {134}. V1 cannot be asymptotically stable in {134}, for otherwise, V1 is stable in ∆ and the lemma is proved. V3 or V4 (or both) cannot be the only asymptotically stable equilibrium (equilibria) in {134}, because Rule 3 of the three-allele case would force V1 to be asymptotically stable in {134}. Suppose that P13 exists. Then, since it is opposite to V2 in {123}, Rule 4 implies that P13 is asymptotically stable in {123}. Therefore, P13 cannot be asymptotically stable in {134}, because it will then be asymptotically stable in ∆, contradicting our assumption. Similarly, P14 cannot be asymptotically stable in {134}. (This argument does not work for P34 , since V2 is the most asymptotically stable vertex in {234}.) The above argument implies that either P134 or P34 is asymptotically stable in {134}. We now show that neither case can happen. From our assumption, V2 is asymptotically stable on the edges {12}, {23}, and {24}, so r22 > r12 , r32 , r42 . If P134 is asymptotically stable in {134}, then J evaluated at P134 = (¯p1 , 0, p¯ 3 ) has three eigenvalues, two of which are negative and the third equals
(r2 − r )|P134 < (r2 − r11 )|P134 = (r21 − r11 )¯p1 + (r23 − r11 )¯p3 + (r24 − r11 )(1 − p¯ 1 − p¯ 3 ) < (r22 − r11 )¯p1 + (r22 − r11 )¯p3 + (r22 − r11 )(1 − p¯ 1 − p¯ 3 ) = r22 − r11 < 0. The first inequality follows from Theorem 2.8. This implies that P134 must be asymptotically stable in ∆, contradicting our assumption.
L. Su, R. Lui / Theoretical Population Biology 81 (2012) 273–283
281
Table 2 Existing patterns of the four-allele model. Type
Pattern
v
⟨{1}⟩
e
⟨{12}⟩ ⟨{13}⟩ ⟨{14}⟩ ⟨{23}⟩ ⟨{24}⟩ ⟨{34}⟩
Values of fitness constants r12
r23
r24
0.645
0.866
0.230
0.572
−4.645
−3.889
0.080
−0.885
5.329 1.782 −3.167 3.279 1.810 −1.279
0.865 4.454 2.509 2.002 2.594 0.632
0.733 0.508 4.474 −1.299 1.870 3.383
0.944 0.256 0.074 0.712 0.010 0.606
3.858 1.835 0.588 5.881 −0.315 3.244
−2.076 −1.512
−3.868
0.835 4.022 5.732 2.483
0.716 0.230 0.063 0.689 0.005 0.114
1.856 0.873 −0.810 4.018 5.737
3.577 1.523 1.470 2.753
1.409 4.242 −2.998 2.020
−2.897 3.226 5.122 1.573
0.943 0.959 0.853 0.646
5.087
f
⟨{234}⟩ ⟨{134}⟩ ⟨{124}⟩ ⟨{123}⟩
1.264 1.555
4.180 2.970 4.589 −0.425
0.512 0.324 0.255 0.068
2.383 2.094 −2.331 −0.739
i
⟨{1234}⟩
4.057
5.051
3.359
0.950
2.632
3.038
0.930
3.735
2v
⟨{1}, {2}⟩ ⟨{1}, {3}⟩ ⟨{1}, {4}⟩
−0.205
0.793
0.608
−3.324
−0.032
−1.043
0.635
−1.045 −4.584
0.820 0.597 0.605
0.513
0.700 0.865
0.273 −0.647
0.156 0.432 0.451
−4.459 −0.770 −0.681
5.105 4.925 3.670 4.160 2.110 3.144 −1.632 −3.867 0.565 −2.815 −0.191 −2.305 2.414 −2.581 2.992
4.741 −1.871 −3.582 2.306 −0.275 3.972 4.761 1.418 4.954 −2.074 0.830 −4.559 1.039 1.685 −3.570
5.085 3.956 −0.182 −3.241 −3.338 3.528 −3.933 −1.383 −0.700 4.740 4.769 1.434 −1.627 3.234 1.022
0.925 0.688 0.737 0.511 0.910 0.221 0.314 0.406 0.068 0.932 0.195 0.150 0.782 0.742 0.508
−1.244
0.339
2.654 1.795 −0.146 −0.790 −0.186 5.533 −1.129 3.100 3.876 −3.193 1.200 2.959 4.508 −4.604
−3.474
0.326 0.397 0.288 0.260 0.835 0.058 0.230 0.063 0.051 0.097 0.162 0.127 0.625 0.558 0.408
−0.431 −4.996 −3.240 −4.437
2e
⟨{12}, {13}⟩ ⟨{12}, {14}⟩ ⟨{12}, {23}⟩ ⟨{12}, {24}⟩ ⟨{12}, {34}⟩ ⟨{13}, {14}⟩ ⟨{13}, {23}⟩ ⟨{13}, {24}⟩ ⟨{13}, {34}⟩ ⟨{14}, {23}⟩ ⟨{14}, {24}⟩ ⟨{14}, {34}⟩ ⟨{23}, {24}⟩ ⟨{23}, {34}⟩ ⟨{24}, {34}⟩
3.428 5.327 4.557 1.866 2.545 0.487
2.791 5.256 3.719 5.143 −2.247 2.829
3.783 5.264 2.460 −3.156 4.177 4.412
0.937 0.879 0.335 0.068 0.778 0.662
3.787 5.562 −3.421 4.382 4.242 5.165
0.691 0.186 0.245 0.052 0.457 0.360
−2.053
2f
⟨{123}, {124}⟩ ⟨{123}, {134}⟩ ⟨{124}, {134}⟩ ⟨{123}, {234}⟩ ⟨{124}, {234}⟩ ⟨{134}, {234}⟩
−0.402 −0.595 −0.044 −0.177 −4.643 −0.943
−0.709
−0.980 −4.050
v+e
⟨{1}, {23}⟩ ⟨{1}, {24}⟩ ⟨{1}, {34}⟩ ⟨{2}, {13}⟩ ⟨{2}, {14}⟩ ⟨{2}, {34}⟩ ⟨{3}, {12}⟩ ⟨{3}, {14}⟩ ⟨{3}, {24}⟩ ⟨{4}, {12}⟩ ⟨{4}, {13}⟩ ⟨{4}, {23}⟩
0.529 0.135 0.161 0.866 0.711 0.245 0.621 0.827 0.018 0.970 0.193 0.640
4.382 1.884 −0.945 −3.511 0.387 −3.275 −0.948 −3.813 −2.742 1.974 −1.388 3.392
⟨{1}, {234}⟩ ⟨{2}, {134}⟩ ⟨{3}, {124}⟩ ⟨{4}, {123}⟩
−3.360 −4.753
0.031 0.123 0.031 0.466
5.190
2.986
v+f
−2.672 −1.374
−0.857
1.259
0.415 0.299 0.595 0.919 0.392 0.830 0.019 0.130 0.819 0.892 0.181 0.132 0.065 0.690 0.891 0.787
e+f
2v + e
5.506 2.486 0.018 4.516 1.736 2.207
r13
0.592 −3.657 1.256 2.864 0.478 −0.194 −4.282 −4.984 1.478 4.036 1.238
r14
0.619 0.700 5.189 1.760 0.453 5.975 2.035 −3.825 −2.029 −0.844
−2.938
−2.309
1.268 3.208
3.338 −2.370 3.114
4.620 4.917 −3.204
5.807
−1.707
−1.412 −3.565
−0.828 −3.994
⟨{234}, {12}⟩ ⟨{234}, {13}⟩ ⟨{234}, {14}⟩ ⟨{134}, {12}⟩ ⟨{134}, {23}⟩ ⟨{134}, {24}⟩ ⟨{124}, {13}⟩ ⟨{124}, {23}⟩ ⟨{124}, {34}⟩ ⟨{123}, {14}⟩ ⟨{123}, {24}⟩ ⟨{123}, {34}⟩
5.417 −3.541 0.175 5.747 3.471 4.269 3.777 2.964 4.350
4.860 −3.540 4.824 1.816 3.855 5.669 −0.337 −0.371 2.546 3.140 1.870
⟨{1}, {2}, {34}⟩ ⟨{1}, {3}, {24}⟩ ⟨{1}, {4}, {23}⟩ ⟨{2}, {3}, {14}⟩
−4.544 −1.377 −4.799 −4.759
−0.714 −1.810 −4.711 −0.625
2.764 1.146 4.658 4.375 2.810 4.120 4.200 3.052 −3.268 −3.828
−4.035 −4.566 −1.911 1.630
r22
−0.423
0.701
1.059 5.583 −2.363 −3.856 4.697 5.490 −2.028 −4.505 5.197 3.455 2.954 −1.156 5.856 1.397
−4.202 5.428 0.693 5.176 4.116
−1.167 4.055 4.747 −3.083 −1.521 −3.205 −0.280 −1.136 5.523 −2.020 −1.072 −0.165
r33
0.142 0.076 0.084 0.644 0.097 0.105 0.049 0.177 0.017 0.442 0.157 0.182
r34
3.423
−4.934 0.964 3.130 4.518 2.463 3.137 4.797 −1.275 5.243 5.579 3.557 5.089 4.593 2.569 2.407 1.201
−0.423 5.914
−1.309 1.794 3.167 −1.029 −0.038 −3.242 −1.784 −0.136 −2.653
−0.765
0.013 0.040 0.013 0.168
4.376 4.580 −3.918 −3.012
1.719 3.181 5.357 −0.548 5.808 −0.039 −3.574 3.856 −3.336 3.773 2.270 3.590
3.836 3.559 5.038 −4.371 −4.131 4.912 5.682 0.924 2.356 −3.293 5.148 −4.576
0.034 0.016 0.553 0.287 0.301 0.147 0.017 0.074 0.069 0.404 0.117 0.118
3.853 5.052 1.344 5.504 3.626 1.051 1.335 −1.162 5.410 −1.993 2.336 1.429
−4.021 −0.114
−2.943
0.039 0.604 0.887 0.102
−4.368 −0.805 −4.345
3.840
−1.712
4.771
4.316 −3.124 −1.841
4.420
(continued on next page)
282
L. Su, R. Lui / Theoretical Population Biology 81 (2012) 273–283
Table 2 (continued) Type
v + 2e
Pattern
Values of fitness constants r12
r13
r14
r23
r24
⟨{2}, {4}, {13}⟩ ⟨{3}, {4}, {12}⟩
−3.766
5.156 −0.990
−3.753 −4.620
0.371 0.654
−4.471 −2.874
−4.570 −3.153
0.337 0.173
−2.315 −3.952
⟨{1}, {23}, {24}⟩ ⟨{1}, {23}, {34}⟩ ⟨{1}, {24}, {34}⟩ ⟨{2}, {13}, {14}⟩ ⟨{2}, {13}, {34}⟩ ⟨{2}, {14}, {34}⟩ ⟨{3}, {12}, {24}⟩ ⟨{3}, {12}, {14}⟩ ⟨{3}, {24}, {14}⟩ ⟨{4}, {23}, {12}⟩ ⟨{4}, {23}, {13}⟩ ⟨{4}, {12}, {13}⟩
−0.227 −1.856 −0.404
−3.813
0.979
2.773 4.171 −4.144 −4.437 −3.204 −2.479 −3.237 −3.274 −4.936 3.528 2.984 −4.099
3.606 −1.511 −0.495 0.296 5.655 −1.450 0.761 −4.017 −0.919 −3.797
0.300 0.660 0.009 0.270 0.143 0.049 0.685 0.042 0.052 0.375 0.161 0.329
−3.084
−3.799 −4.453
0.325 0.688 0.015 0.420 0.210 0.366 0.948 0.178 0.063 0.673 0.799 0.528
5.467
0.803 −0.344 4.287 4.494 −2.732 −3.793 −3.254 0.015 −1.499 3.847 5.188 1.500 1.895 −0.877 2.025 1.199 0.673 1.467 −0.041 5.628 0.901 1.392 1.111
−4.312
0.955 0.968 0.968 0.827 0.746 0.234 0.968 0.968 0.238 0.968 0.062 0.996
0.831 1.320 1.195 6.382 −8.892 6.965 −1.001 1.570 6.045 1.225 2.113 0.889
1.038 1.424 0.703 6.720 2.122 −8.108 1.415 0.704 2.724 −1.032 2.603 1.444
0.283 0.538 0.538 0.707 0.260 0.156 0.538 0.538 0.077 0.538 0.052 0.555
2.471 0.305 1.747 0.396 1.952 2.797 1.933 2.035 −8.365 1.859 −3.878 1.889
0.535 0.226 0.377
−3.852 −1.163 −4.284
−1.132 −4.188 −0.057
0.352 0.057 0.277
−2.018 −0.069 −1.120
0.395 0.237 0.130 0.484 0.205 0.541 0.548 0.682 0.607 0.766 0.551 0.860 0.311 0.111 0.926 0.519
−2.432
−2.823
2.979 5.705 −4.192 −3.617 4.652 −3.643 3.720 −3.254 −3.190 2.382 4.433 −2.043 1.999 −2.988 2.134
4.454 −2.352 4.009 2.338 −2.035 0.718 1.513 5.078 −2.231 −1.788 4.955 4.430 2.645 5.400 −2.193
0.141 0.049 0.009 0.454 0.008 0.186 0.161 0.605 0.038 0.238 0.097 0.026 0.192 0.033 0.762 0.277
−3.690 −1.943 4.197 2.014 −1.234 1.321 2.771 −3.446 −4.606 2.127 2.871 −1.110 3.784 −3.929 3.009 1.576
1.592
0.261
−1.023 −3.143 3.988 2.237 −2.002 1.574 −3.923 5.745
5.814
−4.302 2.570
−3.317 1.893 2.758 −4.840 −2.714 −1.322
r22
−3.201
r33
r34
3.402 3.986 −2.316 5.373 3.899 −3.816 −4.219 −2.662 −2.777 −2.127 −4.514
2e + f
⟨{123}, {24}, {34}⟩ ⟨{134}, {23}, {24}⟩ ⟨{124}, {23}, {34}⟩ ⟨{234}, {13}, {14}⟩ ⟨{124}, {13}, {34}⟩ ⟨{123}, {14}, {34}⟩ ⟨{134}, {12}, {24}⟩ ⟨{234}, {12}, {14}⟩ ⟨{123}, {24}, {14}⟩ ⟨{134}, {23}, {12}⟩ ⟨{124}, {23}, {13}⟩ ⟨{234}, {12}, {13}⟩
3v
⟨{1}, {2}, {3}⟩ ⟨{1}, {2}, {4}⟩ ⟨{1}, {3}, {4}⟩
−2.741 −2.121
3e
⟨{12}, {13}, {14}⟩ ⟨{12}, {23}, {24}⟩ ⟨{13}, {23}, {34}⟩ ⟨{14}, {24}, {34}⟩ ⟨{12}, {13}, {24}⟩ ⟨{12}, {14}, {23}⟩ ⟨{13}, {12}, {34}⟩ ⟨{13}, {14}, {23}⟩ ⟨{14}, {13}, {24}⟩ ⟨{14}, {12}, {34}⟩ ⟨{23}, {12}, {34}⟩ ⟨{23}, {13}, {24}⟩ ⟨{24}, {12}, {34}⟩ ⟨{24}, {23}, {14}⟩ ⟨{34}, {13}, {24}⟩ ⟨{34}, {14}, {23}⟩
2.912 2.128 −3.349 −4.071 4.656 1.096 3.267 −4.557 −0.518 2.617 1.367 −0.215 2.871 −4.427 0.158 −1.296
4.459 2.709 1.553 −4.102 −3.901 1.628 −3.588 −4.584 3.746 −2.321
4v
⟨{1}, {2}, {3}, {4}⟩
−2.828
−4.664
−3.314
0.890
−1.055
−0.512
0.463
−1.505
1.753 5.325 −3.728
2.081 −3.018 5.576
−4.961
0.128 0.776 0.573
−1.141
4e
⟨{12}, {13}, {24}, {34}⟩ ⟨{12}, {23}, {14}, {34}⟩ ⟨{13}, {23}, {14}, {24}⟩
2.237 −2.682 2.921
0.110 0.014 0.489
2.767 5.005 −4.481
1.332
−0.273 1.435
−5.116 1.849 6.488 1.020 1.063 0.680 1.027 0.724 1.022
0.524
2.049 1.715 2.252 0.815 2.249 0.835 1.169 2.217 1.504 1.959 0.550
−0.102
0.155
0.664 −1.133
−2.384 −1.795
4.833
3.753
−0.650
−0.636 −3.535
5.830
−3.325 2.506
−1.353
3.025 −4.977 1.398 −3.592 5.346 2.757 2.848 0.531 −4.676 −3.237 5.527 −4.351 1.368
5.653 3.137
Thus, P34 must be asymptotically stable in {134}. This implies that
{34} is in the superior case and λ234 > 0 > λ134 .
(10)
Since V1 is unstable in ∆ but asymptotically stable in {12}, it must be unstable in {13} or {14} or both. We have three cases to consider. (i) V1 is unstable on {13} and {14} (r13 > r11 , r14 > r11 ). Let P34 = (0, 0, p¯ 3 ). Then
λ234 − λ134 = (r2 − r1 )|P34 = (r23 − r13 )¯p3 + (r24 − r14 )¯p4 < (r22 − r11 )¯p3 + (r22 − r11 )¯p4 < 0.
(11)
Thus, λ234 < λ134 , which contradicts (10). In (11), we have used the fact that V2 is asymptotically stable in ∆. (ii) V1 is unstable in {13} but asymptotically stable in {14} (r13 > r11 , r14 < r11 ). Let P13 = (ˆp1 , 0, pˆ 3 ). Then
λ413 = (r4 − r )|P13 = (r4 − r1 )|P13
4.160 3.979
= (r41 − r11 )ˆp1 + (r43 − r13 )ˆp3 .
(12)
Since P13 is unstable in {134} but asymptotically stable in {13} and {123}, we have λ413 > 0. From (12) being positive and the fact that we are in case (ii), we have r43 > r13 > r11 . Altogether, we have established the inequalities r43 > r13 > r11 > r22 > r23 .
(13)
To continue, (10), (13), and the third term in (11) is positive imply that r24 > r14 . Now r = r3 = r4 at P34 . This and the fact that P34 is unstable in {234} imply that
λ234 = (r2 − r3 )|P34 = (r23 − r33 )¯p3 + (r24 − r34 )¯p4 > 0.
(14)
From (13) and V2 being asymptotically stable in ∆, we have r24 < r22 < r34 . Hence, (14) implies that r33 < r23 . From Proposition 1.1, we have p¯ 3 = (r44 − r34 )/∆34 and p¯ 4 = (r33 − r34 )/∆34 . Substituting these values into the third term of
L. Su, R. Lui / Theoretical Population Biology 81 (2012) 273–283
283
(11), then, from (10) and ∆34 < 0 (recall that {34} is in the superior case), we have
References
(r13 − r23 )(r34 − r44 ) < (r24 − r14 )(r34 − r33 ).
an der Heiden, U., 1975. On manifolds of equilibria in the selection model for multiple alleles. J. Math. Biol. 1, 321–330. Broom, M., Cannings, C., Vickers, G.T., 1993. On the number of local maxima of constrained quadratic form. Proc. R. Soc. Lond. A 443, 573–584. Bürger, R., 2000. The Mathematical Theory of Selection, Recombination, and Mutation. Wiley, Chichester. Cannings, C., Vickers, G.T., 1988. Patterns of ESS’s II. J. Theoret. Biol. 132, 409–420. Cannings, C., Vickers, G.T., 1989. Patterns and invasions of evolutionarily stable strategies. J. Appl. Math. Comput. 32, 227–253. Cannings, C., Vickers, G.T., 1991. The genealogy of patterns of ESS’s. Selected proceedings of the Sheffield symposium on applied probability. IMS Lecture Notes Monogr. Ser. 18, 193–204. Hughes, P.J., Seneta, E., 1975. Selection equilibria in a multi-allele single-locus setting. Heredity 35, 185–194. Kingman, J.F.C., 1961a. A mathematical problem in population genetics. Proc. Cambridge Philos. Soc. 57, 574–582. Kingman, J.F.C., 1961b. On an inequality in partial averages. Quart. J. Math. 12, 78–80. Mandel, S.P.H., 1959. The stability of a multiple allelic system. Heredity 13, 289–302. Mandel, S.P.H., 1970. The equivalence of different sets of stability conditions for multiple allelic systems. Biometrics 26, 840–845. Mulholland, H.P., Smith, C.A.B., 1959. An inequality arising in genetical theory. Amer. Math. Monthly 66, 673–683. Nagylaki, T., 1992. Introduction to Theoretical Population Genetics. Springer, New York. Scheuer, P.A.G., Mandel, S.P.H., 1959. An inequality in population genetics. Heredity 13, 519–524. Su, L., Sesanker, C., Lui, R., Coexisting stable equilibria in a multiple-allele population genetics model, arXiv:1108.5110v1 [q-bio.PE], http://arxiv.org/abs/1108.5110, 2011. Vickers, G.T., Cannings, C., 1988. Patterns of ESS’s I. J. Theoret. Biol. 132, 387–408. Wright, S., 1969. Evolution and the Genetics of Populations. In: The Theory of Gene Frequencies, vol. 2. University Press, Chicago.
(15)
Substituting the above formulas for p¯ 3 , p¯ 4 into (14) and using the fact that λ234 > 0, we have
(r34 − r24 )(r34 − r33 ) < (r23 − r33 )(r34 − r44 ).
(16)
It can be verified that the left-hand sides of inequalities (15) and (16) are positive. Multiplying these two inequalities, and upon cancellations, we have
(r13 − r23 )(r34 − r24 ) < (r24 − r14 )(r23 − r33 ). This, together with 0 < (r13 − r11 ) < (r13 − r23 ),
0 < (r34 − r13 ) < (r34 − r24 )
0 < (r24 − r14 ) < (r11 − r14 ),
0 < (r23 − r33 ) < (r13 − r33 ),
implies that
(r13 − r11 )(r43 − r13 ) < (r11 − r14 )(r13 − r33 ). From Proposition 1.1, we have pˆ 1 = (r33 − r13 )/∆13 and pˆ 3 = (r11 − r13 )/∆13 . If we substitute these values into (12), then the above inequality and ∆13 < 0 imply that λ413 < 0, which means that P13 is asymptotically stable in ∆, a contradiction. The proof of case (ii) is complete. (iii) V1 is asymptotically stable in {13} but unstable in {14} (r13 < r11 and r14 > r11 ). The proof of this case is similar to case (ii), and will be omitted. The proof of Rule 6 is complete.