Structure of principal eigenvectors and genetic diversity

Nonlinear Analysis 74 (2011) 7285–7295

Contents lists available at SciVerse ScienceDirect

Nonlinear Analysis journal homepage: www.elsevier.com/locate/na

Structure of principal eigenvectors and genetic diversity Peter W. Bates a , Fengxin Chen b,∗ a

Department of Mathematics, Michigan State University, East Lansing, MI 48824, United States

b

Department of Mathematics, University of Texas at San Antonio, San Antonio, TX 78249, United States

article

abstract

info

Article history: Received 16 April 2011 Accepted 19 July 2011 Communicated by S. Ahmad

The main concern of this paper is long-term genotypic diversity. Genotypes are represented as finite sequences (s1 , s2 , . . . , sn ), where the entries {si } are drawn from a finite alphabet. The mutation matrix is given in terms of Hamming distances. It is proved that the long time behavior of solutions for a class of genotype evolution models is governed by the principal eigenvectors of the sum of the mutation and fitness matrices. It is proved that the components of principal eigenvectors are symmetric and monotonely decreasing in terms of Hamming distances whenever the fitness matrix has those properties. Published by Elsevier Ltd

Keywords: Long time behavior Principal eigenvectors

1. Introduction Suppose that there are N possible genotypes labeled by the integers {1, 2, . . . , N } with reversible mutations between genotypes, differing fitnesses, and nonlinear interaction. We consider the following model for the evolution process dP(t ) dt

= (M + F )P(t ) + G(P(t )),

(1.1)

where P(t ) = (P1 (t ), P2 (t ), . . . , PN (t ))′ is the vector giving the population density of genotypes at time t , M = (Mij ) is an N × N mutation matrix, and F = diag(F1 , . . . , FN ) is the fitness matrix, with Fj the fitness rate of genotype j. The vector function G represents the interaction between genotypes. Recently, a great interest has been paid for the modeling of genetic evolution. We refer to [1–5] for references. In parallel models (see [6]), the mutation matrix is assumed to satisfy

 mji− Mij = − mik 

if i ̸= j, if i = j,

(1.2)

k̸=i

and mij (i ̸= j) is the probability that, in a given unit time, the genotype i mutates to genotype j. In various models, genotypes are represented as finite sequences (s1 , s2 , . . . , sn ), where the entries {si } are drawn from a finite alphabet. In the model for DNA evolution, the alphabet is that of {A, C , G, T } (or, {A, C , G, U } for RNA). One can also choose the alphabet as the list of known 20 amino acids to study protein evolution. For simplicity, we take a two letter alphabet {0, 1} for our mathematical model. Following the notation in [3], let Gn denote the set of possible genotypes of length n expressed in the given alphabet, {0, 1}; then there are N = |Gn | = 2n members in Gn and each member S can be represented as a sequence or string of 0s and 1s. If we interpret it as a binary representation of an integer then we obtain

∗

Corresponding author. E-mail address: [email protected] (F. Chen).

0362-546X/$ – see front matter. Published by Elsevier Ltd doi:10.1016/j.na.2011.07.045

7286

P.W. Bates, F. Chen / Nonlinear Analysis 74 (2011) 7285–7295

a labeling of genotypes. That is, if S = (s1 , s2 , . . . , sn ), then, S is the (1 + s1 + s2 2 + s3 22 + · · · + sn 2n−1 )th member. Consequently, S1 = (0, 0, . . . , 0), S2 = (1, 0, . . . , 0), S3 = (0, 1, . . . , 0), . . . , SN = (1, 1, . . . , 1).

(1.3)

If S = (s1 , s2 , . . . , sn ) and S¯ = (¯s1 , s¯2 , . . . , s¯n ) are two genotypes, then the Hamming distance d(S , S¯ ) = dH (S , S¯ ) between them is the∑ total number of entries where they differ. In the case of the alphabet with two members, {0, 1}, it is given by n ¯i |. It is attempting to label genotypes in terms of relative Hamming distance, but it is impossible to d(S , S¯ ) = i=1 |si − s do so since each genotype has n genotypes of distance 1. In Section 4, we show that many results do not depend on any particular order of labeling. Naturally, the mutation probabilities are often taken to be functions of the Hamming distance. For example, in the paramuse model [4], the mutation matrix is specified by

 µ mij = −nµ 0

if d(Si , Sj ) = 1, if i = j, otherwise,

where µ is a small parameter. In this model the probability-per-unit-time of changing one letter is µ and the probabilityper-unit-time of changing more than one letter is 0. In the Eigen model [7,2,8,9], the mutation rates are given by mij = µdH (Si ,Sj ) (1 − µ)n−dH (Si ,Sj ) ,

(1.4)

for i ̸= j. The reasoning behind this expression is that if µ is the probability of a mutation at one site, then the probability of mutating at exactly d specified sites (and not changing the other (n − d) sites) is obviously µd (1 −µ)n−d , since the mutations at different sites are independent. The long-time behavior of solutions to a certain class of nonlinear systems is essentially equivalent to the behavior of solutions to a related linear system. For instance (see [6]), if we use p(t ) = (p1 (t ), . . . , pN (t )) to denote the population sizes of all genotypes, and assume that these populations undergo exponential growth given by (F1 , . . . , FN ) and mutation from genotype i to genotype j given by mij , and let M = (mij ), then p(t ) changes according to dp dt

= (M + F )p.

(1.5)

If we now let P(t ) = p(t )/Σ pj (t ) then P satisfies (1.1) with G(P) = −F¯ (P)P,

(1.6)

and F¯ is given by F¯ (P) =

N −

Fj Pj

− N

j =1

Pj .

(1.7)

j =1

Furthermore, for nonnegative and nontrivial initial data, Σ Pi (t ) = 1 for all t ≥ 0. Of course, the denominator in the expression for F¯ is then unnecessary. However, we may consider (1.1) with G and F¯ as above and conclude that Σ Pi (t ) remains constant, whatever is its initial value. In Section 2, we exhibit a more general class of systems for which the long time behavior is equivalent to that of linear system (1.5). For linear system (1.5), if M has nonnegative off-diagonal entries and is irreducible, M + F + kI, for k large enough, is a nonnegative and irreducible matrix, by the Perron–Frobenius theorem, λ0 = ρ(M + F + kI ), the radius of spectrum of M + F + kI, is an algebraically simple eigenvalue with a nonnegative eigenvector. Moreover, for any eigenvalue µ of M + F + kI , |µ| ≤ λ0 . Note that eigenvalues of M + F is just those of M + F + kI minus k and the eigenvectors are the same as those of M + F + kI. The long time behavior of solutions p(t ) of (1.5) is governed by the principal eigenvectors associated with the principal eigenvalue λ0 − k of M + F . Therefore, the structure of the principal eigenvector is of interest. In Sections 3, 5 and 6 we assume that mij has the structure given by (1.2) and (1.4). Section 3 exhibits the complete set of eigenvalues and eigenvectors of M. Sections 4 and 5 are devoted to studying the symmetry and monotonicity of the principal eigenvector of M + F in terms of Hamming distance, when F has those properties. In Section 6, we study the fast and slow mutation model. Numerical simulation shows that the species with slower mutation prevail provided that the fitness environments remain the same. 2. Linear model Consider the following system of ordinary differential equations du dt

= Au − g (u)u

u(0) = u0

(2.1) (2.2)

P.W. Bates, F. Chen / Nonlinear Analysis 74 (2011) 7285–7295

7287

where A be an N × N real matrix and g is a Lipschitz continuous function from RN to R. We assume that A has the algebraically simple principal eigenvalue α0 with the eigenspace H + = span{n} for some n = (n1 , n2 , . . . , nN ). Note that the principal eigenvalue is characterized by re α0 > Re αj for j > 0. We denote the solution of (2.1) and (2.2) by u(t ) = u(t , u0 ). Theorem 2.1. Assume u(t , u0 ) is bounded, then limt →∞ u(t , u0 ) = θ n for some constant θ . Proof. Let α0 , α1 , . . . , αk be distinct eigenvalues of A and α0 be the principal eigenvalue. Let H − be the generalized eigenspace and H + = span{n}. Then RN = H + ⊕ H − . Write u(t ) as u(t ) = a(t )n + v(t ), where a(t ) = (u(t ), n∗ ) and v ∈ H − , and n∗ is the normalized principal eigenvector of A∗ corresponding to the simple eigenvalue α0 . Then (n∗ , v(t )) = 0 and a(t ) satisfies da dt

= α0 a − g (u)a

(2.3)

a(0) = (n∗ , u0 ).

(2.4) t

Therefore a(t ) = a(0)e 0 (α0 −g (u))ds . On the other hand, v(t ) satisfies dv dt

= PAv − g (u)v

(2.5)

where P is the projection from RN onto H − , that is, for any x ∈ Rn , Px = x2 where x2 is defined by the unique decomposition x = x1 + x2 with x1 ∈ H + and x2 ∈ H − . PA has eigenvalues α1 , . . . , αk and v has the form v(t ) =

k −

cj (t )eRe(αj )t −

t

0 g (u)ds

,

j =1

where cj (t ) are vectors whose components are bounded by polynomials of t. Since a(t ) is bounded,

t

t

t

α0 t − 0 g (u)ds is bounded from above. Note that Re(αj )t − 0 g (u)ds = −(α0 − Re(αj ))t + α0 t − that v(t ) → 0 exponentially as t → ∞. That is, there exist two positive constants C and γ0 such that

0

|v(t )| ≤ C e−γ0 t .

t

(α − g (u))ds = g (u)ds. We deduce 0

(2.6)

Therefore, if u¯ ∈ ω(u0 ), u¯ = limn→∞ u(tn ) = limn→∞ a(tn )n = θ n, for some tn → ∞, where

θ = (u0 , n∗ ) lim e

 tn 0

n→∞

(α0 −g (u))ds

.

Now we prove that the ω-limit set is a singleton. We first prove that the ω-limit set of a(t ) is included in the zero set of the function h(a) = f (a, 0), where f (a, v) = α0 a − g (an + v)a. Let a0 is in the omega limit set of a(t ). Assume h(a0 ) ̸= 0. Without loss of generality, we assume γ = h(a0 ) > 0. There exists a δ > 0 such that f (a, v) > γ /2 for (a, v) with a ∈ [a0 − δ, a0 + δ] and |v| ≤ δ . Since v(t ) → 0, there exists t1 such that|v(t )| ≤ δ for t > t1 . There is a t2 > t1 such that a(t2 ) ∈ [a0 −δ, a0 +δ]. t If a(t ) ∈ [a0 − δ, a0 + δ] and t > t2 , a(t ) = a(t2 ) + t f (a, v)dt ≥ a(t2 ) + γ /2(t − t2 ). Therefore, there exists t3 > t2 such 2 that a(t3 ) > a0 + δ . Then u(t ) ≥ a0 + δ for all t > t3 since f (a0 + δ, v) > γ /2 > 0. Therefore a0 is not a ω-limit point. We proved the claim. Since the ω-limit set is compact and connected, it is a single point or an interval [c , d] for some c < d. The second case cannot occur. If ω(u0 ) = [c , d] then f (a, 0) ≡ 0 for a ∈ [c , d]. Let a¯ = (c + d)/2 be the midpoint of the interval and δ0 > 0 be chosen such that [¯a − δ0 , a¯ + δ0 ] ⊂ (c , d). Since g is Lipschitz locally continuous there exist a constant M such that

|f (a, v)| = |f (a, v) − f (a, 0)| ≤ M |v|, for all a ∈ [c , d] and |v| ≤ δ1 . Choose t4 such that |v(t )| < δ1 for t > t4 and MC e−γ0 t4 ≤ γ0 δ0 . If a(t ) ∈ [c , d], then

∫ t  ∫ t ∫ t    f (a, v)ds ≤ M |v| ds ≤ MC e−γ0 t ds ≤ δ0 .   t0

t0

(2.7)

t0

Since a¯ is a ω-limit point of a(t ), there exists t5 > t4 such that a(t5 ) = a¯ . If a(t ) ∈ [c , d] for t > t5 , then a(t ) = a(t5 ) +

∫

t

f (a, v)ds ∈ [¯a − δ0 , a¯ + δ0 ]. t5

Therefore c is not a ω-limit point of a(t ), which contradicts the assumption.



7288

P.W. Bates, F. Chen / Nonlinear Analysis 74 (2011) 7285–7295

Corollary 2.2. Assume that A is irreducible with nonnegative off-diagonal entries and g (u) = (Au, 1) = any u0 = (u01 , u02 , . . . , u0N ) ≥ 0 satisfying

∑N

i=1

∑N

i =1

(Au)i . For

u0i = 1, then ω(u0 ) = {n}, that is, u(t , u0 ) → n as t → ∞, where

n = (n1 , n2 , . . . , nN ) is the normalized principal eigenvector of A satisfying

∑N

i=1

ni = 1.

Proof. Since A is irreducible with nonnegative off-diagonal entries, by the Perron–Frobenius theorem, A has a nonnegative eigenvector n = (n1 , n2 , . . . , nN ) corresponding to the simple principaleigenvalue α0 . It   is easily seen that ui (t ) ≥ 0 for each i, by the comparison principle. Since d

∑N

i =1

ui (t ) /dt = (Au, 1) 1 −

∑N

i=1

ui (t ) , we deduce that

∑N

i =1

ui (t ) = 1

for all t ≥ 0. Consequently, u(t ) is bounded and ω(u0 ) ̸= ∅. By the theorem, ω(u0 ) = {θ n} for some θ . Since ∑N ∑N  i=1 n0i = 1, we know that θ = 1. i=1 ui (t ) = Corollary 2.3. Assume that M satisfies (1.2) and is irreducible and assume G is given by (1.6). Let P0 = (P01 , P02 , . . . , P0N ) ≥ 0 ∑N with i=1 P0i = R > 0. Then for P(t , P0 ), the solution of (1.1) with initial value P(0) = P0 , one has ω(P0 ) = {Rn}, that is, P(t , P0 ) → Rn as t → ∞, where n = (n1 , n2 , . . . , nN ) is the normalized principal eigenvector of M + F . Proof. The conclusion is a direct consequence of Corollary 2.2 with u(t ) = P(t )/R and A = M + F .



Remark 2.1. From the proof of Theorem 2.1, we can deduce that ω(u0 ) is included in the line span{n} even without the assumption of boundedness of orbit u(t , u0 ). 3. Eigenvalues and eigenvectors of m In this section, we focus on the algebraic structure of mutation among genotypes of length n when they are expressed in the two member alphabet {0, 1} and the Eigen model of mutation is used. To begin with, let S1 , S2 , . . . , SN be the list of labels of the N = 2n genotypes according to the binary order (1.4) and let Mn be the corresponding mutation matrix given by (1.2) and (1.4). We will calculate all of the eigenvalues and eigenvectors of mutation matrix Mn through a simple recursion formula. Consider the Hamming  distance matrix Hn corresponding to the binary order of genotypes of length n; obviously, for n = 1 we have H1 =



0

1 H2 =  1 2

1 0 2 1

0 1

1 2 0 1

1 0

, and for n = 2



2 1 , 1 0

the vertices (or genome sequences) being S1 = (0, 0), S2 = (1, 0), S3 = (0, 1), and S4 = (1, 1). In general, Hn is given recursively by

[

Hn−1 ¯ n −1 H

Hn =

¯ n−1 H , Hn−1 ]

(3.1)

¯ n = Hn + En , and En is the matrix with all entries Eij = 1. Since mij = µdH (Si ,Sj ) (1 − µ)n−dH (Si ,Sj ) , we deduce that where H ˜ n−1 (µ) − Σn (µ)In−1 (1 − µ)M ˜ n−1 (µ) + µ(1 − µ)n−1 In−1 µM

[ Mn =

] ˜ n−1 (µ) + µ(1 − µ)n−1 In−1 µM ˜ n−1 (µ) − Σn (µ)In−1 , (1 − µ)M

(3.2)

where

Σn = −mii =

N − j=1,j̸=i

mij =

n   − n j =1

j

µj (1 − µ)n−j = 1 − (1 − µ)n ,

(3.3)

˜ n−1 (µ) = Mn−1 (µ)+ Σn−1 In−1 , i.e., M ˜ n−1 is obtained from Mn−1 by deleting the diagonal In is the n × n identity matrix, and M entries. ˜ (n−1) = Λ(n−1) + Σn−1 is an Let V(n−1) denote an eigenvector associated with an eigenvalue Λ(n−1) of Mn−1 . Then Λ (n)

˜ n−1 with the same eigenvector V(n−1) . Let V1 = eigenvalue of M (n)

Mn V1

 V(n−1) V(n−1)



(n)

and V2 =



V(n−1)

−V(n−1)



. Then

˜ n−1 V(n−1) − Σn V(n−1) + µM ˜ n−1 V(n−1) + µ(1 − µ)n−1 V(n−1) , = [(1 − µ)M ˜ n−1 V(n−1) − Σn V(n−1) + µM ˜ n−1 V(n−1) + µ(1 − µ)n−1 V(n−1) ]′ (1 − µ)M ˜ (n−1) − Σn + µΛ ˜ (n−1) + µ(1 − µ)n−1 ]V(1n) = [(1 − µ)Λ ˜ (n−1) ) + k(Λ ˜ (n−1) )]V(1n) , = [h(Λ

(3.4)

P.W. Bates, F. Chen / Nonlinear Analysis 74 (2011) 7285–7295

7289

where

˜ (n−1) ) = (1 − µ)Λ ˜ (n−1) − Σn (µ) h( Λ

(3.5)

˜ (n−1) ) = µΛ ˜ (n−1) + µ(1 − µ)n−1 . k(Λ

(3.6)

and

Also, (n)

˜ n−1 V(n−1) − Σn V(n−1) − µM ˜ n−1 V(n−1) − µ(1 − µ)n−1 V(n−1) , = [(1 − µ)M ˜ n−1 V(n−1) + Σn V(n−1) + µM ˜ n−1 V(n−1) + µ(1 − µ)n−1 V(n−1) ]′ − (1 − µ)M

Mn V2

˜ (n−1) − Σn − µΛ ˜ (n−1) − µ(1 − µ)n−1 ]V(1n) = [(1 − µ)Λ ˜ (n−1) ) − k(Λ ˜ (n−1) )]V(2n) . = [h(Λ

(3.7)

(n)

˜ (n−1) ) + k(Λ ˜ (n−1) ) and Λ(2n) = h(Λ ˜ (n−1) ) − k(Λ ˜ (n−1) ) are eigenvalues of Mn with eigenvectors V(1n) = h(Λ n −1 ( n − ˜ 1) = Λ(n−1) + Σn−1 , and V2 , respectively. Since Σn−1 − Σn + µ(1 − µ) = 0, and Λ Therefore, Λ1 (n)

(n) ˜ (n−1) ) + k(Λ ˜ (n−1) ) = Λ ˜ (n−1) − Σn + µ(1 − µ)n−1 = Λ(n−1) , Λ 1 = h( Λ

(3.8)

and (n)

(n)

(n−1)

˜ (n−1) ) = Λ1 Λ2 = Λ1 − 2k(Λ

˜ (n−1) + µ(1 − µ)n−1 ) − 2(µΛ

= (1 − 2µ)Λ(n−1) − 2µΣn−1 − 2µ(1 − µ)n−1 = (1 − 2µ)Λ(n−1) − 2µ.

(3.9)

Therefore we can find and eigenvectors of Mn by iteration.  eigenvalues  −µ µ

For n = 1, M1 = (1)

are v1

=

 1 1

(1)

, and v2

=

µ

(1)

. The eigenvalues of M1 are λ1

−µ  1

−1 , respectively. If (n)

(n−1)

respectively, then λj eigenvectors vj respectively, where (n)

λj

 (n−1) λj = (1 − 2µ)λ(j−n−2n1−) 1 − 2µ

= 0 and λ(21) = −2µ. The corresponding eigenvectors

λ(j n−1) for j = 1, 2, . . . , 2n−1 are eigenvalues of Mn−1 with orthogonal (n)

for j = 1, 2, . . . , 2n are eigenvalues of Mn with orthogonal eigenvectors vj

for 1 ≤ j ≤ 2n−1

(3.10)

for 2n−1 + 1 ≤ j ≤ 2n

and (n)

vj

 [vj(n−1) , v(j n−1) ]′ = [vj(n−1) , −v(j n−1) ]′

for 1 ≤ j ≤ 2n−1

(3.11)

for 2n−1 + 1 ≤ j ≤ 2n ,

where, for column vectors a and b, [a, b]′ is the column vector with the entries of a followed by those of b. Therefore we have the following theorem. Theorem 3.1. λj = −1 +(1 − 2µ)j are eigenvalues of Mn with geometric multiplicity

  n j

, for j = 0, 1, . . . , n. The eigenvectors

are the columns of the Hadamard matrix given by the recursion:

[ H1 = [1],

H2 =

1 1

]

1 , −1

and H2n = H2n−1 ⊗ H2 .

(3.12)

The mutation matrix M, of course, depends on the order in which we choose to list the genotypes. For any particular ordering S1 , S2 , . . . , SN of Gn , there is a permutation σ on Gn such that Sσ (1) , Sσ (2) , . . . , Sσ (N ) is the binary ordering given by (1.3). Since M is defined in terms of the relative distance between genotypes, it is similar to any other mutation matrix corresponding to a different ordering of Gn . In particular, if M0 is the mutation matrix corresponding to the binary ordering of Gn , and M is the permutation matrix corresponding to any other ordering of Gn , then M = Q T M0 Q , where Q is the product of permutation matrices. Therefore, the eigenvalues of M are the same as those for M0 and the eigenvectors are the images of the eigenvectors of M0 under the permutation matrix Q T .

7290

P.W. Bates, F. Chen / Nonlinear Analysis 74 (2011) 7285–7295

4. Symmetry in distance for the entries of the principal eigenvector of m + f Let S1 , S2 , . . . , SN be a list of elements of Gn , and suppose M = (mij ) with mij depending only on the Hamming distance between Si and Sj . Suppose that F = diag(F1 , . . . , FN ) and that Fi depends only on the Hamming distance between S1 and Si . We will show that the entries of the principal eigenvector of M + F also depends only on the Hamming distances from S1 . Suppose that Fi = f (d(Si , S1 )) for some function f on {0, 1, . . . , n} and mij = m(d(Si , Sj )) for i, j = 1, 2, . . . , N for some functions m satisfying m(x) ≥ 0 for x > 0. That is, we assume that Fi = Fj if d(Si , S1 ) = d(Sj , S1 ) and mij = mkl if d(Sk , Sl ) = d(Si , Sj ). Assume also that M is irreducible. Then, by the Perron–Frobenius theorem, M + F has a nonnegative principal eigenvector, which is denoted by V0 = (v1 , v2 , . . . , vN ). Then we have the following theorem. Theorem 4.1. There exists a function g on {0, 1, . . . , n} such that vi = g (d(Si , S1 )) for i = 1, 2, . . . , N. Remark 4.1. The theorem shows that the components of principal eigenvector only depends on the Hamming distances of the genotypes. To prove the theorem we first establish some lemmas. Lemma 4.2. Given t = (t1 , t2 , . . . , tn ) ∈ Gn , there exists a permutation σ on Gn such that σ (t ) = 0 = (0, 0, . . . , 0) and d(S , S¯ ) = d(σ (S ), σ (S¯ )), for any S , S¯ ∈ Gn . Proof. For given t = (t1 , t2 , . . . , tn ) ∈ Gn , we define n permutations on {0, 1} by

τi =



i0 , i1 ,

if ti = 0, if ti = 1,

where i0 is the identity permutation on {0, 1} and i1 is the permutation on {0, 1} defined by i1 (t ) =



1, 0,

if t = 0, if t = 1.

For any S = (s1 , s2 , . . . , sn ) ∈ Gn , define σ (S ) = (τ1 (s1 ), τ2 (s2 ), . . . , τn (sn )). Then σ is a permutation on Gn . Moreover,

σ (t ) = (τ1 (t1 ), τ2 (t2 ), . . . , τn (tn )) = (0, 0, . . . , 0).

(4.1)

Suppose S = (s1 , s2 , . . . , sn ) and S¯ = (¯s1 , s¯2 , . . . , s¯n ) are two arbitrary elements in Gn . Then, d(σ (S ), σ (S¯ )) = d((τ1 (s1 ), τ2 (s2 ), . . . , τn (sn )), (τ1 (¯s1 ), τ2 (¯s2 ), . . . , τn (¯sn )))

=

n −

|τi (si ) − τi (¯si )| =

i=1

n −

|si − s¯i | = d(S , S¯ ),

(4.2)

i=1

where we have used the fact |τi (si ) − τi (¯si )| = |si − s¯i |.



Lemma 4.3. Let S1 , S2 and S3 in Gn satisfying d(S1 , S2 ) = d(S1 , S3 ). Then there exists a permutation σ on Gn such that σ (S1 ) = S1 and σ (S2 ) = S3 and d(S , S¯ ) = d(σ (S ), σ (S¯ )), for any S , S¯ ∈ Gn . Proof. Applying Lemma 4.2 for the case t = S1 , we deduce that there exists a permutation σ1 on Gn such that σ1 (S1 ) = 0 and d(S , S¯ ) = d(σ1 (S ), σ1 (S¯ )), for any S , S¯ ∈ Gn . Let S˜i = σ1 (Si ) for Si ∈ Gn . Then d(S˜i , S˜1 ) = d(S˜i , 0) = d(Si , S1 ) for i = 2, 3. Since d(S1 , S2 ) = d(S1 , S3 ), it follows d(S˜2 , 0) = d(S˜3 , 0). That is, S˜2 and S˜3 contains the same number of 1s. Therefore, there exists a permutation τ on {1, 2, . . . , n} such that S˜3 = (˜sτ (1) , s˜τ (2) , . . . , s˜τ (n) ), where S˜2 = (˜s1 , s˜2 , . . . , s˜n ). τ induces a unique permutation σ2 on Gn defined by

σ2 (S ) = (sτ (1) , sτ (2) , . . . , sτ (n) ), for S = (s1 , s2 , . . . , sn ) ∈ Gn . Obviously σ2 (0) = 0, σ2 (S˜2 ) = S˜3 , and d(σ2 (S ), σ2 (S¯ )) = d(S , S¯ ) for all S , S¯ ∈ Gn . Let σ = σ1−1 σ2 σ1 . Then

σ (S1 ) = σ1−1 (σ2 (σ1 (S1 ))) = σ1−1 (σ2 (0)) = σ1−1 (0) = S1 , σ (S2 ) = σ1−1 (σ2 (σ1 (S2 ))) = σ1−1 (σ2 (S˜2 )) = σ1−1 (S˜3 ) = S3 ,

(4.3)

P.W. Bates, F. Chen / Nonlinear Analysis 74 (2011) 7285–7295

7291

and d(σ (S ), σ (S¯ )) = d(S , S¯ ).



Let K0N = {x = (x1 , x2 , . . . , xN ) ∈ RN | xi = xj whenever d(Si , S1 ) = d(Sj , S1 ) for i, j = 1, 2, . . . , N }.

(4.4)

That is, xi depends only on the Hamming distances of Si and S1 . Lemma 4.4. The cone K0N is invariant under the linear transformation T on RN defined by Tx = Mx. Proof. Assume x = (x1 , x2 , . . . , xN ) ∈ K0N and i, j are integers so that d(Si , S1 ) = d(Sj , S1 ) and i ̸= j. Let Si and Sj be the corresponding genotypes. We want to show that (Mx)i = (Mx)j , where

(Mx)i =

N −

mik xk ,

(4.5)

mjk xk ,

(4.6)

k=1

and

(Mx)j =

N − k=1

and mik = m(d(Si , Sk )), and mjk = m(d(Sj , Sk )). By Lemma 4.3, there is a permutation σ on Gn such that σ (S1 ) = S1 , σ (Si ) = Sj and d(σ (S ), σ (S¯ )) = d(S , S¯ ) for all S , S¯ ∈ Gn . Then d(Si , Sk ) = d(σ (Si ), σ (Sk )) = d(Sj , σ (Sk )), and d(σ (Sk ), S1 ) = d(σ (Sk ), σ (S1 )) = d(Sk , S1 ). Therefore, xσ (Sk ) = f (d(σ (Sk ), S1 )) = f (d(Sk , S1 )) = xk . Therefore, mik xk = mjσ (Sk ) xσ (Sk ) . From (4.5) and (4.6) we deduce that (Mx)i = (Mx)j .



Proof of Theorem 4.1. Since m(x) ≥ 0 for x > 0, mij = m(d(Si , Sj )) is nonnegative for i, j = 1, 2, . . . , N. Let A = M +F +k0 I, −j

where k0 is a constant large enough such that A is a nonnegative matrix. Let x0 > 0 ∈ K0N . Then λ0 Aj x0 approaches the principal eigenvectors of A which is also the principal eigenvector V0 of M + F as j → ∞. Since K0N is invariant under transformation A, we deduce that V0 ∈ K0N . That is, there exists a function g such that vj = g (d(Sj , S1 )) for j = 1, 2, . . . , N.  5. Monotonicity of the principal eigenvector In this section, we study the monotonicity of the principal eigenvectors of M + F . Let S1 , S2 , . . . , SN be a list of ordering of Gn . Suppose M = (mij ) is given by (1.2) and (1.4) and F = diag(F1 , . . . , FN ). Suppose that Fi = f (d(Si , S1 )) for some decreasing function f on {0, 1, . . . , n}, that is, Fi depends only on the Hamming distances. In the previous section, we show that the principal eigenvectors of M + F also depends only on the Hamming distances. Suppose the components of F decrease in terms of the Hamming distances. We show that the components of the principal eigenvectors of M + F also decrease in terms of the Hamming distances. Denote V0 = (v1 , v2 , . . . , vN ) the principal eigenvector of M + F . Then we have the following theorem. Theorem 5.1. Suppose µ < 12 . vj = f1 (dH (Sj , S1 )) for some positive decreasing function on {0, 1, 2, . . . , n}. f1 is strictly decreasing unless f is a constant function. To proof the theorem, we first consider the case S1 , S2 , . . . , SN is the binary ordering (1.3). From (3.1), we easily deduce the following lemma. Lemma 5.2. Hn is composed of 2n−1 × 2n−1 block matrices of the form each row of blocks is a permutation of the first row of blocks.



s s+1

s+1 s



, where s ∈ {0, 1, . . . , n − 1}. Moreover,

7292

P.W. Bates, F. Chen / Nonlinear Analysis 74 (2011) 7285–7295

Lemma 5.3. Let x ∈ RN . Then 2n−1

(Mx)2k−1 − (Mx)2k =

−

fjk (x2j−1 − x2j ),

j =1

for k = 1, 2, . . . , 2

n −1

fj 1 =



, where

−Σn − µ(1 − µ)n−1 m1(2j−1) − m2(2j)

if j = 1, if j = 2, . . . , 2n−1

and (f1k , . . . , f2kn−1 )′ = Pσ (f11 , . . . , f21n−1 )′ , and σ is a permutation on {1, 2, . . . , 2n−1 } and Pσ is the corresponding permutation matrix. Proof. 2n−1

(Mx)1 − (Mx)2 =

−

(m1(2j−1) x2j−1 − m2(2j−1) x2j−1 + m1(2j) x2j − m2(2j) x2j ).

(5.1)

j =1

Since m2(2j) = m1(2j−1) and m1(2j) = m2(2j−1) , by Lemma 5.2, we have m1(2j−1) x2j−1 − m2(2j−1) x2j−1 + m12j x2j − m2(2j) x2j = fj1 (x2j−1 − x2j ). By Lemma 5.2, for k > 1, the coefficients of (Mx)2k−1 − (Mx)2k are just a permutation of those of (Mx)1 − (Mx)2 .



Let K1N = {x = (x1 , x2 , . . . , xN ) ∈ RN | xi = xj whenever d(Si , S1 ) = d(Sj , S1 ), and xi ≥ xj whenever d(Si , S1 ) ≤ d(Sj , S1 ) for i, j = 1, 2, . . . , N }. Lemma 5.4. K1N is invariant under the transformation A = M + KI if K ≥ 2Σn and µ <

(5.2) 1 . 2

Proof. Let x ∈ K1N and y = Ax. Then by Lemma 5.3, 2n−1

y2k−1 − y2k =

−

f¯jk (x2j−1 − x2j ),

(5.3)

j =1

where (f¯1k , . . . , f¯2kn−1 )′ = Pσ (f¯11 , . . . , f¯21n−1 )′ , f¯11 = K + f11 and f¯j1 = fj1 for j = 2, 3, . . . , N /2. Since K ≥ 2Σn , f¯11 = K + f11 =

K − Σn − µ(1 − µ)n−1 > 0. And, for some positive integer s, fj1 = m1(2j−1) − m2(2j) = µs (1 − µ)n−s − µs+1 (1 − µ)n−s−1 =

µs (1 − µ)n−s−1 (1 − 2µ) > 0 for j = 2, 3, . . . , N /2. Therefore y2k−1 ≥ y2k , for k = 1, 2, . . . , N /2. For 0 < d < n, consider two integers j, k with d(Sj , S1 ) = d and d(Sk , S1 ) = d + 1. Any sequence of the form (0, s2 , . . . , sn ) has ∑n an odd index, by the binary indexing (1.3). Take S2i−1 = (0, s2 , . . . , sn ) with the Hamming distance from S1 being d, i.e., i=2 si = d. Note that S2i = (1, s2 , . . . , sn ) and its Hamming distance from S1 is d + 1. Lemma 4.4 shows that if d(Sj , S1 ) = d and d(Sk , S1 ) = d + 1, then yj = y2i−1 and y2i = yk . We have already proved that y2i−1 ≥ y2i for any i and it is noted that always d(S2i−1 , S1 ) + 1 = d(S2i , S1 ). We deduce yj = y2i−1 ≥ y2i = yk . Therefore, y ∈ K1N .  Proof of Theorem 5.1. First we consider the case of binary ordering. Let A = M + F + k1 I, where k1 is a constant large enough such that A is a nonnegative matrix. Since A is irreducible, A has a simple principal eigenvalue λ0 > 0 and a positive j principal eigenvector. Let x0 > 0 ∈ K1N . Then λ0 A−j x0 approaches a principal eigenvector of A which is also the principal eigenvector V0 of M + F as j → ∞. Apparently, F keeps K1N invariant. We deduce that V0 ∈ K1N since K1N is invariant under transformation M. That is, there exists a decreasing function f1 such that vj = f1 (d(Sj , S1 )) for j = 1, 2, . . . , N. Obviously, f1 is strictly decreasing if it is not a constant function. For general list of ordering, there exists a permutation on Gn that keeps Hamming distances invariant and that maps the order to the binary order.  6. Fast mutation In biological applications, the mutation rates of genotypes are often decided by a few nucleotides. Certain traits of the few nucleotides decide if the genotype has a fast or slow mutation rate. Consider a genotype S = (s1 , s2 , . . . , sn ) of geno-length n and assume that the nth node sn is such a node, that is, different sn of a genotype attributes to a different mutation rate. Let S1 , S2 , . . . , SN be the binary ordering of Gn given by (1.3). Assume that the genotypes have the faster mutation rate µ1

P.W. Bates, F. Chen / Nonlinear Analysis 74 (2011) 7285–7295

7293

when the nth node has the value 1 while the genotypes have the slower mutation rate µ ≪ µ1 < 1/2 when the nth node has the value 0. Assume M (µ, µ1 ) = (Mij (µ, µ1 )) is given by (1.2) and (1.4). Then we have Mij (µ, µ1 ) =



µdH (Si ,Sj ) (1 − µ)n−dH (Si ,Sj ) d (S ,S ) µ1H i j (1 − µ1 )n−dH (Si ,Sj )

if j ≤ N /2, if N /2 < j ≤ N

(6.1)

for i ̸= j. For the eigenvalues and eigenvectors of M (µ, µ1 ), we have the following theorem.

 Theorem 6.1.

(n−1) Vj (n−1) θj Vj



are eigenvectors of M (µ, µ1 ) with eigenvalues

(1 − µ)λ(j n−1) (µ) − µ + θj (µ1 λ(j n−1) (µ1 ) + µ1 ), (n−1) and θj = θj± , where (Vj , λ(j n−1) ), for j = 1, 2, . . . , N /2, are eigen-pairs of Mn−1 given by Theorem 3.1, and

K±

θj± =



K 2 + 4µµ1 (1 − 2µ)j (1 − 2µ1 )j 2 µ1 ( 1 − 2 µ1 ) j

,

(6.2)

and K = (1 − µ1 )(1 − 2µ1 )j − (1 − µ)(1 − 2µ)j . Proof. Note that Mn (µ, µ1 ) can be written as Mn (µ, µ1 ) =

[

˜ n−1 (µ) − Σn (µ)In−1 (1 − µ)M ˜ µMn−1 (µ) + µ(1 − µ)n−1 In−1

] ˜ n−1 (µ1 ) + µ1 (1 − µ1 )n−1 In−1 µ1 M ˜ n−1 (µ1 ) − Σn (µ1 )In−1 , (1 − µ1 )M

(6.3)

where Mn−1 (·) and Σn (·) are defined in Section 3. We seek eigenvectors of Mn (µ, µ1 ) of the form V (n) = [V , θ V ]′ , where ˜ n−1 (µ) associated with eigenvalue λ˜ (n−1) (µ), given by Section 3 and θ to be chosen. Note that V is V is an eigenvector of M independent of µ. Thus

˜ n−1 (s)V = λ˜ (n−1) (s)V , M for s = µ, µ1 . By simple computation, we have Mn (µ, µ1 )V (n) = [Λ1 V , Λ2 θ V ]′

˜ n−1 (µ) − Σn (µ) + [µ1 λ˜ n−1 (µ1 ) + µ1 (1 − µ1 )n−1 ]θ and Λ2 = [µλ˜ n−1 (µ) + µ(1 − µ)n−1 ]θ −1 + (1 − where Λ1 = (1 − µ)λ n − 1 µ1 )λ˜ (µ1 ) − Σn (µ1 ). Therefore, if θ is chosen such that Λ1 = Λ2 , V (n) is an eigenvector of Mn (µ, µ1 ) with eigenvalue Λ1 . Equation Λ1 = Λ2 is equivalent to the quadratic equation aθ 2 + bθ + c = 0 with a = µ1 λ(n−1) (µ1 ) + µ1 , b = (1 − µ)λ(n−1) (µ) − (1 − µ1 )λ(n−1) (µ1 )−µ+µ1 and c = µλ(n−1) (µ)+µ. Since λ(n−1) (µ) are of the form λj = −1 +(1 − 2µ)j , for j = 0, 1, . . . , n − 1, we deduce that the quadratic form has two distinct roots given by (6.2).  Now we will discuss the principal eigenvectors of M + F , where F = diag(F1 , F2 , . . . , FN /2 , F1 , F2 , . . . , FN /2 ). That is, the genotypes are divided into two groups with respect to fast and slow mutation rates respectively and assume that the fitness rates are same for each group. Let V = (v1 , v2 , . . . , vN ) denote the normalized principal eigenvector of M + F . If Fj = f (d(Sj , S1 )) for j = 1, 2, . . . , N /2 for some monotonely decreasing function f , then similar arguments to those in Sections 4 and 5 establish the symmetry and monotonicity results. That is, there exist two monotonely decreasing function g1 and g2 such that Vj = g1 (d(Sj , S1 )) for j = 1, 2, . . . , N /2 and Vj = g2 (d(Sj , S1 )) for j = N /2 + 1, N /2 + 2, . . . , N. The case for j = 1 in Theorem 6.1 gives the principal eigenvector of M (µ, µ1 ) + kI. It is a multiple of V = 1, . . . , 1, θ , . . . , θ



µ

′

,

(6.4)

where θ = µ . This indicates that when the fitness rates for every geno-type are same. The geno-types with slower mutation 1 rates prevail. Numerical simulation suggests that it is still true for general fitness matrix F . We use Matlab to compute the principal eigenvector of M (µ, µ1 )+ In + hF , where µ is chosen in [10−6 , 0.1] and µ1 = pµ with p a positive integer; F is the diagonal matrix whose diagonal entries are produced by the random function, sorting in descendent order in term of Hamming distances; h is a parameter varying from 0 to 1.

7294

P.W. Bates, F. Chen / Nonlinear Analysis 74 (2011) 7285–7295

Fig. 1. Eigenvector for n = 3.

Fig. 2. Eigenvector for n = 5.

Fig. 3. Eigenvector for n = 6.

Figs. 1–3 are the eigenvector of M (µ, µ1 ) + I + F corresponding to n = 3, 5, 6, respectively, where µ = 0.01 and

µ1 = 2µ1 . In Fig. 1,

F = diag[0.06650, 0.05887, 0.05887, 0.03352, 0.06650, 0.05887, 0.05887, 0.03352]. Fig. 4 is the eigenvector of M (µ, µ1 ) + I + F corresponding to n = 10, where µ = 0.1 and µ1 = 2µ1 and F is produced by the random function, sorting in descendent order in term of Hamming distances. For example, Fig. 1 shows that genotypes

P.W. Bates, F. Chen / Nonlinear Analysis 74 (2011) 7285–7295

7295

Fig. 4. Eigenvector for n = 10.

S2 = (1, 0, 0), S3 = (0, 1, 0) are 1 distance away from S1 = (0, 0, 0) and S6 = (1, 0, 1), S7 = (0, 1, 1) are 1 distance away from S1 = (0, 0, 1). The component values of the principal eigenvector for S2 and S3 are the same but are larger than those corresponding to S6 and S7 . Similar comparisons are true for other components respectively. In biological applications, µ varies from the scale of 10−3 to 10−6 . Our choices of µ is just to produce graphs for demonstration. For small µs, the graph shows the exact same qualitative property: the first half of vi s prevails the second half of vi s componentwise. In fact numerical results shows that vi > vi+N /2 for i = 1, 2, . . . , N /2 and

vN /2+1 1 ≤ . v1 p

(6.5)

This result can be proved by the perturbation theory for the cases (i) F = kI +ϵ F (1) for ϵ ≪ 1 (ii), F1 ≫ Fj for j = 2, . . . , N /2. Numerical simulation shows that the model presents the universal biological principle: the species with slower mutation prevail provided that the fitness environments remain the same. Acknowledgments The authors would like to thank Professor Richard Lenski for valuable discussions. The authors are also grateful to the referee for his useful suggestions and comments. References [1] [2] [3] [4] [5] [6]

E. Baake, M. Baake, H. Wagner, Ising quantum chain is equivalent to a model of biological evolution, Phys. Rev. Lett. 78 (1997) 559–562. M. Eigen, J. McCaskill, P. Schuster, Molecular quasi-species, J. Phys. Chem. 92 (1988) 6881–6891. C.L. Epstein, Anderson localization, non-linearity and stable genetic diversity, J. Stat. Phys. 124 (1) (2006) 25–46. K. Jain, J. Krug, Adaptation in simple and complex fitness landscapes, vol. 1, 2005, pp. 1–42. arXiv:q-bio,PE/0508008. L. Peliti, in: T. Riste, D. Sherrington (Eds.), Physics of Biomaterials: Fluctuations, Self-Assembly and Evolution, Kluwer, Dordrecht, 1996, p. 267. J.-M. Park, M.W. Deem, Schwinger Boson formulation and solution of the Crow–Kimura and Eigen models of quasispecies theory, J. Stat. Phys. (2006) 1–33. [7] M. Eigen, Molekulare selbstorganisation und evolution (self organization of matter and the evolution of biological macro molecules), Naturwissenschaften 58 (1971) 465–523. [8] M. Eigen, J. McCaskill, P. Schuster, The molecular quasi-species, Adv. Chem. Phys. 75 (1989) 149–263. [9] M. Eigen, P. Schuster, A principle of natural self-organization. Part A. Emergence of the hyper cycle, Naturwissenschaften 64 (1977) 541–565.

Further reading [1] [2] [3] [4] [5]

J.F. Crow, M. Kimura, An Introduction to Population Genetics Theory, Harper and Row, New York, 1970. S.S. Chow, C.O. Wilke, C. Ofria, R.E. Lenski, C. Adami, Adaptive radiation from resource competition in digital organisms, Science 305 (2004) 84–86. J. Hermisson, H. Wagner, M. Baake, Four-state quantum chain as a model for sequence evolution, J. Stat. Phys. 102 (2001) 315–343. P.B. Rainey, M. Travisano, Adaptive radiation in a heterogeneous environment, Nature 394 (1998) 69–72. D.B. Saakian, E. Muñoz, C.-K. Hu, M.W. Deem, Quasispecies theory for multiple-peak fitness landscapes, Phys. Rev. E 73 (2006) 041913-110.