Evolutionary dynamics on strongly correlated fitness landscapes
aa r X i v : . [ q - b i o . P E ] S e p Evolutionary dynamics on strongly correlated fitness landscapes
Sarada Seetharaman and Kavita Jain ∗ Theoretical Sciences Unit, Jawaharlal Nehru Centre for AdvancedScientific Research, Jakkur P.O., Bangalore 560064, India and Theoretical Sciences Unit and Evolutionary and Organismal Biology Unit,Jawaharlal Nehru Centre for Advanced Scientific Research, Jakkur P.O., Bangalore 560064,India (Dated: November 8, 2018)We study the evolutionary dynamics of a maladapted population of self-replicating sequenceson strongly correlated fitness landscapes. Each sequence is assumed to be composed of blocks ofequal length and its fitness is given by a linear combination of four independent block fitnesses. Amutation affects the fitness contribution of a single block leaving the other blocks unchanged andhence inducing correlations between the parent and mutant fitness. On such strongly correlatedfitness landscapes, we calculate the dynamical properties like the number of jumps in the mostpopulated sequence and the temporal distribution of the last jump which is shown to exhibit ainverse square dependence as in evolution on uncorrelated fitness landscapes. We also obtain exactresults for the distribution of records and extremes for correlated random variables.
PACS numbers: 87.23.Kg, 02.50.Cw, 02.50.Ey
I. INTRODUCTION
Fitness is a measure of an organism’s ability to survive and reproduce [1, 2]. Fit organisms produce more offspringand can dominate the population while the less fit ones can be lost. Mathematically, fitness is a non-negative realnumber associated with a sequence which is a string of L letters whose meaning is context dependent. For example,fitness represents the stability of a sequence of amino acids in case of proteins, activity for an enzyme or replicationrate for a genetic sequence of nucleotides. On plotting the fitness as a function of the sequence, a fitness landscapeis obtained. Empirical measurement of fitness landscapes is very hard since the number of sequences increasesexponentially with the sequence length L . However several qualitative features particularly the topography of thefitness landscapes has been deduced in experiments on proteins and microbes either by an explicit construction of thefitness landscapes for small L ( < ∼
5) or indirect measurements of relevant quantities. These experiments show that thefitness landscapes can have smooth hills as evidenced by fast adaptation in some proteins [3] or multiple peaks as seenin microbial populations that evolve towards different fitness maxima [4–6] and enzymes with short uphill paths tothe global fitness peak [7]. Detailed studies in which all or a set of mutants from wild type to an optimum are createdand their fitness measured [8] have also indicated the smooth [9] and rugged [10, 11] nature of the fitness landscapes.The topography of the fitness landscapes is related to the correlations between the fitness of the sequences. If thefitness of the mutants of a sequence is correlated to that of the sequence so that the mutant fitness does not differappreciably from the parent sequence, a smooth fitness landscape is generated whereas if the mutant fitnesses areindependent random variables so that the fitness of one sequence has no influence over the fitness of other sequencesdiffering from it by even a single mutation, a highly rugged fitness landscape with multiple optima is obtained. Severaltheoretical models such as NK model [12], block model [7] and rough Mt. Fuji-type model [13] in which correlationscan be tuned via a parameter have been proposed. Although realistic fitness landscapes exhibit intermediate degreeof correlations [14], much of the theoretical work has focused on the limiting cases of fitness landscapes with highdegree of correlation but single fitness peak [15] and no correlations but a large number of local optima [16–19].In this article, we study the evolutionary dynamics on the fitness landscapes generated by the block model [7] inwhich a sequence of length L is assumed to be composed of independent units or blocks of length ℓ . As explained later(see Sec. II), the correlations can be varied by changing the block length ℓ from maximally correlated case with ℓ = 1to maximally uncorrelated one with ℓ = L . Here we focus on the block model with ℓ = 2 which generates fitnessesthat are strongly correlated but to a lesser degree than the maximally correlated case and the fitness landscape ismoderately rugged i.e. exhibits several peaks.The evolution model that we work with here describes the deterministic evolution of an infinitely large populationof asexually replicating sequences. In this model, the population is initially distributed in such a manner that the ∗ Electronic address: [email protected], [email protected] high fitness sequences have lower initial population and vice versa but the population of all the sequences increaseslinearly with time [16]. As time goes on, a highly fit subpopulation is able to overcome the poor initial condition anddominate the population until an even fitter population overtakes it. This process goes on until the globally fittestsequence becomes the most populated one. The stepwise dynamics of such leadership changes termed jumps havebeen studied when the fitness variables are completely uncorrelated [16–18]; here we are interested in this problemwhen the fitnesses are strongly correlated. As explained in the next section, in the context of this problem, it is alsorelevant to consider the sequence with largest fitness amongst sequences carrying D mutations relative to a referencesequence and whose fitness is a record in that its fitness exceeds the fitness of all the sequences with less than or equalto D mutations. Thus we are led to study the statistics of maximum [20] and records [21] when the random variablesare not independent, both of which have been much less studied unlike the problem when the random variables areindependently distributed.Our detailed analysis presented in Sec. IV shows that the statistical properties studied depend only on whetherthe number of mutations D are odd or even and whether D lies below or above L/
2. This simplification allowsus to tackle the problem analytically and to find exact expressions for various quantities. On uncorrelated fitnesslandscapes, it has been shown that the average number of leadership changes increases as √ L and the timing of thelast jump exhibits a 1 /t dependence [17, 18]. For evolution on the class of strongly correlated fitness landscapesstudied here, we find that the average number of jumps is a constant independent of L but the time dependence ofthe distribution of the last jump remains unaffected. The average number of records is found to increase linearly with L as in maximally rugged case albeit with a larger prefactor. II. SHELL MODEL ON CORRELATED FITNESS LANDSCAPES
Consider a microbial population evolving in a complex environment that can be modeled by rugged fitness land-scapes. At large times, most of the population resides at the globally fittest sequence of the fitness landscape and dueto mutations, a suite of mutants is also present. If the population size is infinite, a nonzero population is present at all the sequences whereas a finite population produces only a small number of mutants around the fittest sequence [19].Now if the environment is changed by changing (say) the nutrient medium of the microbial population, the fittestsubpopulation before the environment change will be typically maladapted to the new environment and dependingon the total population size, a small population may be present at the new fittest sequence. We are interested infinding how the new global maximum is reached starting with an initial condition in which all the population is atthe sequence that was globally fittest before the environmental change. The exact evolutionary dynamics of averageHamming distance and overlap function has been studied on permutationally invariant [22] and uncorrelated [23]fitness landscapes. Here we will be tracking the evolution of the most populated sequence in time on strongly corre-lated fitness landscapes. The dynamics of the adaptation process is studied in the setting when the population size isinfinite so that the fluctuations in the population frequency of a sequence can be neglected and one can work with theaverages. In the following, we begin with the quasispecies model of biological evolution [24, 25] and proceed to relateit to the shell model [16]. We then define and explain some properties of the block model [7] of correlated fitnesslandscapes that we shall use in the paper.We consider an infinitely large population of binary sequences where a sequence σ ≡ { σ , ..., σ L } , σ i = 0 , L letters. The population evolves by the elementary processes of replication and mutation. If the fitness A ( σ )of the sequence σ is defined as the average number of copies produced per generation and p σ ← σ ′ is the probability thata sequence σ ′ mutates to the sequence σ at a Hamming distance D ( σ, σ ′ ) = P Li =1 ( σ i − σ ′ i ) from it, the populationfraction X ( σ, t ) of sequence σ at time t evolves according to the following quasispecies equation [17, 24]: X ( σ, t + 1) = P σ ′ p σ ← σ ′ A ( σ ′ ) X ( σ ′ , t ) P σ ′ A ( σ ′ ) X ( σ ′ , t ) (1)where the denominator on the right hand side ensures the normalisation condition P σ X ( σ, t ) = 1 is satisfied atall times. Assuming that the mutations occur independently at each locus with a probability µ , the mutationalprobability p σ ← σ ′ = µ D ( σ,σ ′ ) (1 − µ ) L − D ( σ,σ ′ ) . In the following discussion, we will use the unnormalised population Z ( σ, t ) defined through the relation X ( σ, t ) = Z ( σ, t ) / P σ ′ Z ( σ ′ , t ) as it obeys a linear equation given by Z ( σ, t + 1) = X σ ′ µ D ( σ,σ ′ ) (1 − µ ) L − D ( σ,σ ′ ) A ( σ ′ ) Z ( σ ′ , t ) (2)As discussed at the beginning of this section, we consider the evolution of the dominant population starting with theinitial condition X ( σ,
0) = Z ( σ,
0) = δ σ,σ (0) where σ (0) is the fittest sequence before the change in the environment.Earlier work [16–18] has shown that the statistical properties of the most populated sequence in the quasispeciesmodel are accurately described by a simplified shell model which approximates the solution of (2) by Z ( σ, t ) ∼ µ D ( σ,σ (0) ) A t ( σ ) (3)The above equation can be heuristically obtained as follows: on iterating (2) with the given initial condition, thepopulation Z ( σ, ∼ µ D ( σ,σ (0) ) A ( σ (0) ) for small µ . Thus all the mutants become available in one generation for aninfinitely large population even after starting with a highly localised population. If the mutations are neglected forfurther evolution i.e. Z ( σ, t + 1) = A ( σ ) Z ( σ, t ), the solution (3) is immediately obtained. A detailed analysis hasshown that the behavior of Z ( σ, t ) in shell model matches the quasispecies dynamics (2) only for highly fit sequencesand at short times. However it captures the behavior of the most populated genotype exactly at all times [18] andtherefore we will work with the shell model in the rest of the article.Taking the logarithm of both sides of (3) and rescaling the time by | ln µ | , the logarithmic population E ( σ, t ) ∼ ln Z ( σ, t ) is seen to increase linearly in time with a slope ln A ( σ ) = w ( σ ), E ( σ, t ) = − D ( σ, σ (0) ) + w ( σ ) t (4)According to the above equation, there are (cid:0) LD (cid:1) populations in a shell of radius D from the initial sequence which havethe same initial condition but different growth rates. As the fittest population in each shell grows the fastest, one canwork with the largest fitness w (max) ( D ) in each shell. Labeling the fittest sequence in a shell by its shell number, (4)can be rewritten as E ( D, t ) = − D + w (max) ( D ) t (5)Thus we arrive at a model in which the fitness variables w (max) ( D ) are independent but non-identically distributed.We mention in passing that the above linear dynamics when the slope variables are independent and identicallydistributed (i.i.d.) have appeared in a shell model with one-dimensional fitness [16], a gas of elastically colliding hardcore particles [26] and a spin glass model with random entropy [27].As mentioned earlier, we are mainly interested in the dynamics of the most populated sequence whose fitnesschanges abruptly or jumps in time. Due to (5), the leader in shell D ′ is overtaken by a fitter population in shell D > D ′ at time T ( D, D ′ ) given by T ( D, D ′ ) = D − D ′ w (max) ( D ) − w (max) ( D ′ ) (6)Initially the sequence σ (0) is the leader. As the overtaking time must be positive, the population in shell D = 1can be a leader provided w (max) (1) > w (max) (0). Similarly, the fittest sequence in shell 2 can be the most populatedsequence if w (max) (2) = max { w (max) (0) , w (max) (1) , w (max) (2) } . In general, a population at Hamming distance D > D ′ < D or in other words, the fitness in shell D is a record. As noted in earlier works, it is not sufficientto be a record fitness in order that the corresponding sequence can be the dominant sequence [16, 17] and a jumpoccurs only when the current leader is overtaken in minimum time. Due to this constraint, not all record sequencesparticipate in the jump process and thus the number of records is an upper bound on the number of jumps.We next define the block model introduced by Perelson and Macken who were motivated by the observation thatmany biomolecules such as proteins and antibodies are composed of domains or partitions [7]. As shown in Fig. 1,a sequence of length L is divided into B independent blocks of equal length ℓ = L/B , ≤ ℓ ≤ L . Each blockconfiguration is assigned a fitness value which may also depend on the position of the block (locus-dependent blockfitness model). In this article, we assume that a block configuration at any location in the sequence carries the sameblock fitness (see Sec. V also). These 2 ℓ block fitnesses are chosen independently from a common distribution withsupport on the interval [ l, u ] where l and u are respectively the lower and upper limits of the block fitness distribution.The sequence fitness is given by the average of the corresponding block fitnesses.The topographical features such as the number of local maxima depends on ℓ . For a sequence to be a local maximum,each of its B blocks must also be a local maximum. Since a sequence is composed of independent blocks and theaverage number of local optima of a sequence of length ℓ with i.i.d. fitness is 2 ℓ / ( ℓ + 1), it follows that the averagenumber n opt of local maxima of a sequence of length L and block length ℓ is given by (2 / ( ℓ + 1) /ℓ ) L [7]. Exceptfor ℓ = 1 for which there is a single local (same as global) fitness peak, n opt increases with increasing ℓ and L (seeFig. 1). For ℓ = 2 with which we work in this article, there are ≈ . L local optima on an average. Arguing as abovefor local maximum, it can be seen that the globally fittest sequence is composed of identical blocks with the largest FIG. 1: (Color online) Left panel: Block model for ℓ = 2. The initial sequence and its mutant have correlated fitnesses asthey have several blocks in common. Right panel: Average number n opt of local maxima in block model as a function of ℓ forvarious L . block fitness and has the average fitness given by 2 ℓ R ul df f p ( f )[ R fl df ′ p ( f ′ )] ℓ − . Thus the initial sequence σ (0) can bechosen to be any one of the 2 ℓ sequences with same blocks. For convenience, we choose it to be the one with all 0s.An attractive feature of the block model is that the correlations can be tuned with the block length ℓ . As illustratedin Fig. 1, when two sequences have at least one common block, their respective fitnesses are correlated. For ℓ = 1, thesequence fitnesses are maximally correlated while for ℓ = L , we obtain the model with maximally uncorrelated fitnesses.This statement can be quantified by considering the correlation function C ,j between the fitness w = w ( σ (0) ) of theinitial sequence and the fitness w j of a sequence at Hamming distance D = 1 from it given by w j = ( L − ℓ ) f + ℓf j L , j = 0 , ..., ℓ (7)where f j is the fitness of the block of length ℓ with 1 in the j th position. Using the fact that f j ’s are i.i.d. variables,we can write the correlation function as [7] C ,j = h w w j i − h w ih w j i = L − ℓL σ (8)where σ is the variance of the block fitness distribution p ( f ). The above correlation function is largest at ℓ = 1 andvanishes at ℓ = L . Similarly the correlation function C i,j amongst the one mutant neighbors given by [28] C i,j = h w i w j i − h w i ih w j i = "(cid:18) L − ℓL (cid:19) + δ i,j (cid:18) ℓL (cid:19) σ (9)is a monotonically decreasing function of ℓ for i = j .In the following, we will study the shell model defined by (4) where the fitness w ( σ ) is chosen from the block model.In the next section, we briefly discuss the dynamics of the shell model for the two limiting cases namely ℓ = 1 and L . Section IV which forms the major part of the paper discusses the evolutionary dynamics when the block length ℓ = 2. Finally we conclude with a discussion of our results in Sec. V. In the rest of the article, we will assume thatthe sequence length L is an even integer. III. SHELL MODEL DYNAMICS WHEN BLOCK LENGTH ℓ = 1 AND ℓ = L In this section, we briefly discuss the evolutionary dynamics on the fitness landscapes for the two limits of the blockmodel namely block length ℓ = 1 and L . When block length is equal to one, the sequence fitnesses are maximallycorrelated. Let f and f denote the block fitness of the two blocks { } and { } respectively. Then the fitness w ( D )of a sequence at Hamming distance D from the initial sequence is given by w (max) ( D ) = ( L − D ) f + Df L (10)The fitness landscape thus generated is permutationally invariant since there is a single distinct fitness at each D fromthe initial sequence. It is easy to see that the average number of jumps on fitness landscapes with ℓ = 1 is half. Thisis because if f > f , a jump cannot occur after D = 0. If f > f , as the time taken by the population at D > D = 0 given by T ( D,
0) = Dw (max) ( D ) − w (max) (0) = 1 f − f , D > D , all the populations overtake E (0 , t ) at the same time and hence one jump occurs with probability1 /
2. Thus the average number of jumps is 1 / L . The average number of records from the aboveconsiderations is given by 1 + ( L/ ℓ = L has been studied earlier [16, 17, 25]. It hasbeen shown that the average number of records is given by (1 − ln 2) L for any underlying block fitness distribution[17] and the average number of jumps by √ Lπ/
IV. SHELL MODEL DYNAMICS WHEN BLOCK LENGTH ℓ = 2 For the rest of the article, we will consider the case when the sequences are built by blocks of length ℓ = 2. Theblock fitness is given by f , f , f and f corresponding to the blocks { , } , { , } , { , } and { , } respectively.Let n i denote the number of blocks with fitness f i , i = 0 , ...,
3. Then the fitness of a sequence of length L with D mutations obtained by averaging over B = L/ w n ,n ( D ) = ( L − D − n − n ) f + 2 n f + 2 n f + ( D − n − n ) f L (12)In the above expression, since the total number of blocks equals L/ D of a sequence fromthe initial sequence is given by n + n + 2 n = D , we get n = ( D − n − n ) / n = ( L − D − n − n ) /
2. As n and n must be integers, for even D , both n , n must be either even or odd whereas for odd D , either n shouldbe odd and n even or vice versa. Besides, for D ≤ L/
2, the conditions n + n ≤ D, n ≤ D must be satisfied as n ≥ D > L/ n + n ≤ L − D, n ≤ L − D are required to ensure the non-negativity of n .As mentioned in Sec. II, in order to be the globally fittest sequence, a sequence must be composed of blocks of thesame type. For ℓ = 2, the global maximum can thus occur at D = 0 , L/ L corresponding to f , either f or f and f being the largest of the block fitnesses respectively. Starting with all the population at the initial sequence,we wish to find the properties of the jumps by which the most populated sequence reaches the global maximum. Inthe following subsections, we discuss the statistics of extremes (Sec. IV A), records (Sec. IV B) and jumps (Sec. IV C)on correlated fitness landscapes. A. Distribution of the largest fitness at constant Hamming distance
It was shown in [29] that the total number of distinct fitnesses at a fixed D increases as D . However for questionsof interest, we need to consider only the sequence with the largest fitness. To identify such sequences, we first considerfitnesses with fixed n + n = n , n > n , n satisfy the conditions described above. As the coefficient of f and f in fitness w n ,n depends on n + n , assuming f > f and comparing w n ,n − n and w n ′ ,n − n ′ for all n ′ = n , we find that w n ,n − n > w n ′ ,n − n ′ for n ′ < n . The fitness w n ,n − n can be the largest w (max) n ,n ( D ) of allthe fitnesses at fixed D and n provided n = n . We next compare w n, and w n + k, , k = 0 for D ≥
2. Since w n + k, ( D ) = w n, ( D ) − ( k/L )( f + f − f ), it follows that for D ≤ L/ w (max) n ,n ( D ) = w D, ( D ) if f > f and f − f + f < w , ( D ) if f > f and f − f + f > D is odd (13b) w , ( D ) if f > f and f − f + f > D is even (13c)The above conditions are independent of D (except for the parity) and as we shall see, this property simplifies theproblem considerably. For D > L/
2, the largest possible fitness is obtained on replacing D by L − D in the abovediscussion. The corresponding conditions for the case when f < f are obtained by interchanging fitnesses f and f and the indices n and n in the preceding equation. w P E ( w , D ) δ =1 (a) w P E ( w , D ) r=0.1(b) FIG. 2: (Color online) Distribution P E ( w, D ) of maximum fitness for (a) r = 0 . r = 0 . δ = 1 and(b) δ = 1 (solid) and δ = 2 (broken) with r = 0 . We consider the cumulative distribution P E ( w, D ) that all the fitnesses at constant D are smaller than w . Asargued above, for even D ≤ L/
2, only one of the three fitnesses w D, , w ,D and w , can be the largest. For unboundedunderlying distribution p ( f ) with f >
0, we can thus write P E ( w, D ) = Z ∞ df p ( f ) Z ∞ df p ( f ) Z ul df p ( f ) Z ∞ df p ( f )Θ( w − w D, )Θ( w − w ,D )Θ( w − w , ) (14)= Z w − r df p ( f ) Z w − (1 − r ) f r df p ( f ) "Z w − (1 − r ) f r df p ( f ) (15)where Θ( ... ) is the Heaviside step function and r = D/L < /
2. Specifically, for p ( f ) = δf δ − e − f δ , δ >
0, we have P E ( w, D ) = δ Z w − r df f δ − e − f δ (cid:20) − e − ( w − (1 − r ) fr ) δ (cid:21) (cid:20) − e − ( w − (1 − r ) f r ) δ (cid:21) (16)= a Z df e − af " − e − a (cid:18) (1 − r )(1 − f /δ ) r (cid:19) δ − e − a (1 − r )(1 − − r − r f /δ )2 r ! δ (17)where a = ( w/ (1 − r )) δ . The probability P E ( w, D ) = d P E /dw that the largest sequence fitness with D mutations hasa value w can be easily computed for δ = 1 and is given by P E ( w, D ) = e − w − r − e − wr + 2 e − w r − r − e − w r − r − (cid:16) e − w − r − e − wr (cid:17) r − r + 6 r + 4 e − w − r ) r − r + 8 r (18)The above distribution shown in Fig. 2a for two values of r shifts towards right with increasing r as the average h w i E = 1 + (55 r/
36) [29]. Figure 2b shows that the extreme value distribution at fixed r peaks at larger w as δ increases. This is contrary to the general expectation that if the tail of the underlying distribution decays fast, theprobability of finding a large maximum value of a set of S random variables should also decrease when S ≫
1. Here asthe number of independent random variables is merely four, the tail of the block fitness distribution is not adequatelysampled and the block fitnesses lie closer to the average value which increases with increasing δ thus resulting in thebehavior seen for P E ( w, D ).When D is odd, one can write down an expression for the extreme distribution but for large L , it reduces to thatobtained for even D . The results for extreme statistics when D > L/ D by L − D inthe above discussion [29]. B. Statistics of record fitnesses
In this subsection, we are interested in finding the probability that a fitness w n ,n ( D ) is a record i.e. it exceedsall the fitnesses in the shell D ′ ≤ D . As only the largest fitness at constant D can possibly be a record, we need toconsider only such fitnesses. Unless otherwise mentioned, we assume f > f so that the largest fitness at constant D can be one of the following: w , ( D ) if D is even and w , ( D ) otherwise, w D, ( D ) for D ≤ L/ w L − D, ( D ) for D > L/ D ≤ L/
2, the fitness w D, ( D ) can be a record if it exceeds all the fitnesses at constant D as well as the oneswith number of mutations D ′ < D . The first condition is met if (13a) is satisfied. As the conditions in (13a) areindependent of D (barring parity), the largest fitness in a shell with D ′ mutations is also w D ′ , ( D ′ ) , < D ′ < D .Then w D, ( D ) > w D ′ , ( D ′ ) for all D ′ ≥ f > f . Thus the probability of w D, ( D ) being a record can be writtenas P ( w D, is a record) = Z ul Y i =0 df i p ( f i )Θ( f − f )Θ( f − f )Θ(2 f − f − f ) (19)= Z ul df p ( f ) Z uf df p ( f ) Z f l df p ( f ) Z f − f l df p ( f ) (20)For D > L/
2, the fitness w L − D, ( D ) can be record if w L − D, ( D ) > w L − D ′ , ( D ′ ) for D ′ ≥ L/ w L − D, ( D ) >w D ′ , ( D ′ ) for D ′ < L/ f > f and f − f + f < f > f and f < f . Thus we can write P ( w L − D, is a record) = Z ul Y i =0 df i p ( f i )Θ( f − f )Θ( f − f )Θ( f − f )Θ(2 f − f − f ) (21)= Z ul df p ( f ) Z uf df p ( f ) Z f l df p ( f ) Z f − f f df p ( f ) (22)For even D , the fitness w , ( D ) can be a record if w , ( D ) > w , ( D ′ ) for even D ′ and w , ( D ) > w , ( D ′ ) for odd D ′ besides satisfying (13c). If f < f , the fitness w , ( D ) can be a record if f > f and f > f − f . The lasttwo conditions can be split into two cases, namely f > f if f > f and f > f − f if f < f . Similarly, for f > f , the conditions for w , ( D ) to be a record are obtained by interchanging f and f . Combining all the aboveconditions, we get P ( w , is a record) = 2 Z ul Y i =0 df i p ( f i )Θ( f − f )Θ( f − f )Θ( f + f − f ) (23)= 2[ Z ul df p ( f ) Z uf df p ( f ) Z f l df p ( f ) Z u f − f df p ( f )+ Z ul df p ( f ) Z f l df p ( f ) Z f l df p ( f ) Z f l df p ( f )] (24)For odd D , the fitness w , ( D ) , D > w , ( D ) > w , ( D ′ ) for odd D ′ < D and w , ( D ) > w , ( D ′ ) for even D ′ < D . The last two conditions are satisfied if f < f and f < f respectively. Thenthe probability of w , ( D ) , D > P ( w , is a record) = Z ul Y i =0 df i p ( f i )Θ( f − f )Θ( f − f )Θ( f − f )Θ( f + f − f ) (25)= Z ul df p ( f ) Z uf df p ( f ) Z f l df p ( f ) Z u f − f df p ( f ) (26)The above expression holds for D = 1 also as w , (1) is a record if w , (1) > w , (0) which implies f < f besides f < f . D P R ( D ) FIG. 3: (Color online) Variation of record occurrence probability P R ( D ) with Hamming distance D for L = 64.
1. Record occurrence distribution
Using the results derived above, we now calculate the probability P R ( D ) that a record occurs in the shell with D > P R (0) = 1. Figure 3 shows that P R ( D ) is not a smooth function - the value of P R ( D ) dependson whether D is odd or even and whether it is below or above L/
2. Thus four distinct cases arise due to this characterof P R ( D ) which we will discuss below. We shall find that the distribution P R ( D ) is universal i.e. does not dependon the choice of the underlying distribution of the block fitness. As the global maximum is the last record and theonly global maximum for D > L/ /
4, we may expect the record occurrence probability for
D > L/ D ≤ L/ Even D : When D is even, either w D, ( D ) or w ,D ( D ) can be a record for D ≤ L/ w L − D, ( D ) or w ,L − D ( D ) for D > L/ w , ( D ) for any even D .Thus the probability of even D for D ≤ L/ P R ( D ) = 2 P ( w D, is a record) + P ( w , is a record) (27)= 2 Z ul df p ( f ) Z uf df p ( f ) Z f l df p ( f ) + 2 Z ul df p ( f ) Z f l df p ( f ) Z f l df p ( f ) Z f l df p ( f ) (28)= 23 + 112 = 34 , D ≤ L/ D > L/
2, the record occurrence probability is given by P R ( D ) = 2 P ( w L − D, is a record) + P ( w , is a record) (30)= 2 Z ul df p ( f ) Z uf df p ( f ) Z f l df p ( f ) Z uf df p ( f ) + 112 = 14 , D > L/ Odd D : For w D, ( D ) , D > D is odd, the same conditions as for even D are required so that(20) holds. Thus the probability of a shell with odd D, < D ≤ L/ P R ( D ) = 2 [ P ( w D, is a record) + P ( w , is a record)] (32)= 2 Z ul df p ( f ) Z uf df p ( f ) Z f l df p ( f ) = 23 , D ≤ L/ D > L/
2, the probability that w L − D, ( D ) is a record is given by (22) and w , ( D ) is a record by (26). Thusthe probability of a record occurring for odd D > L/ P R ( D ) = 2 [ P ( w L − D, is a record) + P ( w , is a record)] (34)= 2 Z ul df p ( f ) Z uf df p ( f ) Z f l df p ( f ) Z uf df p ( f ) = 16 , D > L/
2. Record value distribution
In this subsection, we calculate the probability P R ( w, D ) that the record value in shell D is smaller than or equal to w . For this purpose, we will need the probability P R ( w ( D ) ≤ w ) that the fitness w ( D ) in shell D does not exceed w . Asthe record value distribution is not expected to be universal, we will restrict ourselves to distributions with support onthe interval [0 , ∞ ). It can be checked that the cumulative distribution P R ( w, D ) gives the probability P R ( w ) obtainedin the last subsection when w → ∞ . Below we present the expressions for D ≤ L/ D > L/
Even D : As seen for the distribution of extreme values in Sec. IV A, the distribution for the record value is a functionof the ratio r = D/L for even D . Since either w D, ( D ) or w , ( D ) can be a record for even D ≤ L/
2, the cumulativeprobability P R ( w, D ) = 2 P ( w D, ≤ w ) + P ( w , ≤ w ) where P ( w D, ≤ w ) = Z ∞ Y i =0 df i p ( f i ) Θ( w − w D, )Θ( f − f )Θ(2 f − f − f )Θ( f − f )= Z w df p ( f ) Z w − f r + f f df p ( f ) Z f df p ( f ) Z f − f df p ( f ) (36)and P ( w , ≤ w ) = 2 Z ∞ Y i =0 df i p ( f i ) Θ( w − w , )Θ( f − f )Θ( f − f )Θ( f + f − f )= 2 Z w df p ( f ) Z f df p ( f ) Z f df p ( f ) Z w − f r + f f df p ( f )+ 2 Z w df p ( f ) Z w − f r + f f df p ( f ) Z f df p ( f ) Z w − f r + f f − f df p ( f ) (37)Using these expressions, it is straightforward to see that P R ( w, D ) = 2 Z w df p ( f ) Z f df p ( f ) Z f df p ( f ) Z w − f r + f f df p ( f )+ 2 Z w df p ( f ) Z w − f r + f f df p ( f ) Z f df p ( f ) Z w − f r + f df p ( f ) (38)Taking the derivative of the last expression with respect to w , we obtain the distribution P R ( w, D ) that the recordvalue equals w . For p ( f ) = e − f , the distribution P R ( w, D ) is given by P R ( w, D ) = e − w + 2 e − w r − e − wr − r − e − w r − r + e − w (3 − r )1 − r + 8 r + e − wr r − e − w (3 − r )1 − r (5 − r ) (39)The above result for the record value distribution is compared with the extreme value distribution P E ( w, D ) givenby (18) in Fig. 4 for two values of r . Though the record fitness is also the extreme fitness in shell D , the converseis not true and the distribution P R ( w, D ) < P E ( w, D ) for all w at a given D . We also note that the most probablerecord value in shell D is smaller than the corresponding extreme value - this behavior is unlike that for uncorrelatedfitnesses for which record is a maximum of a larger set of independent variables.0 w FIG. 4: (Color online) The probability distribution of the extreme value (solid lines) given by (18) and record value (dashedlines) by (39) for r = 0 . r = 0 . p ( f ) = e − f . Odd D : To find the record value distribution for odd D , besides P ( w D, ≤ w ), we require the cumulative probability P ( w , ≤ w ) that the fitness w , ( D ) in shell D does not exceed w . The latter can written as P ( w , ≤ w ) = Z ∞ Y i =0 df i p ( f i ) Θ( w − w , )Θ( f − f )Θ( f + f − f )Θ( f − f )Θ( f − f )= Z w df p ( f ) Z L ( w − f D + f f df p ( f ) Z f df p ( f ) Z Lw − ( L − D − f − f D − f − f df p ( f ) (40)which reduces to the second integral in (37) for L ≫
1. Thus for large L , the cumulative distribution P R ( D, w ) forodd D is also a function of r . However unlike extreme value distribution for odd D , the distributions for even and odd D do not match for L ≫ P ( w , ≤ w ) and P ( w , ≤ w )do not coincide.
3. Distribution of the number of records
To find the probability N R ( n ) that the total number of records equals n , we first calculate the record configurationprobability Q ( { w n ,n ( D ) } ) defined as the probability that all the elements in the set { w n ,n ( D ) } are records. Thisdistribution depends on the location of the global maximum. If f is the largest block fitness, the global maximumoccurs at D = 0 and obviously there are no records beyond D = 0 in this case.When f is not a global maximum and f > f , only four record configurations occur with a nonzero probability.When the fittest block has a fitness f , a record cannot occur beyond D = L/ f − f − f must be positive. Thus the fitness w D, ( D ) for all D ≤ L/ Q ( w , (1) , ..., w L/ , ( L/ f is the largest, the records occur until D = L at a spacing of one or two depending on thesign of f − f as explained below:(i) From the discussion at the beginning of Sec. IV B, it is evident that when f < f , the only set of fitnesses thatcan be a record are w , ( D ) for all D . Using the conditions in (24), it follows that when f < f < f < f , a record1occurs only in even D shells. As f i ’s are independent and identically distributed (i.i.d.) random variables, all 4! blockfitness configurations are equally likely and therefore we get Q ( w , (2) , w , (4) , ..., w , ( L )) = 124 (42)(ii) If f > f (and f ), the fitness w , (1) is a record. The next record depends on the sign of 2 f − f − f . From(24) and (26), it follows that if 2 f − f − f <
0, the fitness w , ( D ) is a record for all even D and w , ( D ) for allodd D with probability Q ( w , (1) , w , (2) , ..., w , ( L − , w , ( L )) = Z ul df p ( f ) Z uf df p ( f ) Z f l df p ( f ) Z u f − f df p ( f ) (43)If 2 f − f − f >
0, due to (20) and (22), the fitnesses w D, ( D ) for all D ≤ L/ w L − D, ( D ) for all D > L/ Q ( w , (1) , ..., w L/ , ( L/ , w L/ − , ( L/ , ..., w , ( L )) = Z ul df p ( f ) Z uf df p ( f ) Z f l df p ( f ) Z f − f f df p ( f )(44)From the above discussion, it is evident that the total number of records (ignoring the one at D = 0) can beeither L/ L (see (43) and (44)). The probability N R ( n ) of total number n of records isindependent of underlying block fitness distribution and is given by N R ( L/
2) = 2 (cid:18)
14 + 124 (cid:19) = 712 , N R ( L ) = 212 = 16 (45)where we have used that twice the sum of (43) and (44) equals (35). The average number R of records can be foundusing N R ( n ) or P R ( D ) and is given by R = L X n =1 nN R ( n ) = L X D =1 P R ( D ) = 11 L ≈ . L (46)for any even L . C. Reaching the global maximum
As discussed in Sec. II, all records are contenders for being a leader; however only those records for which theovertaking time is minimised qualifies to be a jump [16–18]. Like records, the statistics of jumps depends on thelocation of the global maximum. If f is the fittest block, the unmutated sequence with fitness w , (0) = f is theleader throughout.If f ( > f ) is the global maximum, the last record and hence the last jump occurs at D = L/
2. Since the time ofintersection T (0 , D ) of the population E ( D, t ) , D ≤ L/ E (0 , t ) given by T = T (0 , D ) = Dw D, ( D ) − w , (0) = L f − f ) , D ≤ L/ D , all the populations overtake the population of the initial sequence at the same point. Thus allthe record populations participate in the evolutionary race. But as the population E ( L/ , t ) has the largest fitness,it becomes the final leader thus leading to a single jump when f (or f ) is the largest fitness.If the global maximum is f which occurs at D = L , the following cases as discussed in Sec IV B 3 arise:(i) If f < f , the population with the record fitness w , ( D ) , D ≤ L overtakes that with the initial fitness w , (0)at a time given by T , = T (0 , D ) = Dw , ( D ) − w , (0) = Lf − f , D ≤ L (48)so that all the populations with record fitness w , ( D ) intersect at the same time and the population of the globalmaximum at D = L takes over in a single jump.2(ii) If f > f and 2 f − f − f <
0, the population with fitness w , ( D ) for all even D and w , ( D ) for all odd D intersects E (0 , t ) at the following intersection time: T (0 , D ) = Dw , ( D ) − w , (0) = LD ( D − f + 2 f − ( D + 1) f , for odd D (49) T (0 , D ) = Dw , ( D ) − w , (0) = Lf − f , for even D (50)By virtue of the condition 2 f − f − f <
0, the intersection time for odd D is greater than that for even D . Thereforethe current leader at D = 0 is overtaken by D = L resulting in a single jump at time T , = L/ ( f − f ).If 2 f − f − f >
0, the record fitnesses are w D, ( D ) for D ≤ L/ w L − D, ( D ) for D > L/
2. The populationscorresponding to these fitnesses overtake the leader at D = 0 at time T (0 , D ) = L f − f ) , D ≤ L/ T (0 , D ) = DL (2 D − l ) f + 2( L − D ) f − Lf , D > L/ D ≤ L/ w L/ , ( L/
2) is the largest fitness, the firstjump occurs when the population of the sequence with fitness w L/ , ( L/
2) overtakes E (0 , t ). The next change inleader occurs at the point of intersection of populations involving the fitness w L − D, ( D ) , D > L/ T , = T ( L/ , D ) = L f − f ) , D > L/ D independent. Thus the population E ( L, t ) is the leader after E ( L/ , t ) and the global maximum isreached in two jumps.
1. Distribution of the number of jumps
It is obvious that when any block fitness other than f is the globally largest fitness, there will be at least one jump(corresponding to globally fittest being the final leader) so that the probability of at least one jump equals 3 /
4. Inaddition, there can be one more jump when f is the global maximum and 2 f − f − f > p of the second jump is given by p = 2 Z ul df p ( f ) Z uf df p ( f ) Z f l df p ( f ) Z f − f f df p ( f )Θ( u − f + f ) (54)Thus the average number J of jumps is given by (3 /
4) + p . As p is independent of L , the average number of jumpsis of order unity for any underlying distribution but the constant p is not universal. For instance, when the blockfitnesses are chosen from an exponential probability distribution, p = 5 / ≈ .
069 while for uniform distribution, itequals 5 / ≈ .
2. Temporal jump distribution
We are interested in the probability P ( t ) that the last jump occurs at time t > p ( f ) = e − f .This distribution is a sum of the probability P A ( t ) that the last jump occurs at t when f or f is a global maximumand P B ( t ) when f is a global maximum. We first consider the cumulative probability P A ( t ) = R t dt ′ P A ( t ′ ) which onusing that f (or f ) is a global maximum and (47) gives P A ( t ) = 2 Z ul Y i =0 df i p ( f i )Θ( t − T )Θ( f − f )Θ( f − f )Θ( f − f ) (55)= 2 Z ul + L t df p ( f ) Z f − L t l df p ( f ) Z f l df p ( f ) Z f l df p ( f ) (56)3
10 100 1000 t P ( t ) FIG. 5: (Color online) Log-log plot of the distribution P ( t ) of the last jump for p ( f ) = e − f and L = 100. The broken line hasa slope − Differentiating P A ( t ) with respect to time t yields P A ( t ) = − L t d P A dǫ = Lt Z ul + ǫ df p ( f ) p ( f − ǫ ) Z fl dgp ( g ) ! (57)where we have defined ǫ = L/ (2 t ). For large times t ≫ L/
2, the integral on the right hand side of the above equationreduces to the probability G (0) that the gap between the globally largest and the second largest in a set of i.i.d.random variables is zero [17]. Thus the probability P A ( t ) decays as ∼ LG (0) /t at large times.When f is the largest fitness (and f > f ), the last jump can occur at times given by (48), (50) and (53). As T , = T , , the corresponding conditions (discussed in Sec. IV B 3) on the block fitnesses can be combined to give thefollowing cumulative probability P ( t ) = Z ul +2 ǫ df p ( f ) Z f − ǫl df p ( f ) Z f f l df p ( f ) Z f l df p ( f ) (58)and the probability distribution P ( t ) = Lt Z ul +2 ǫ df p ( f ) p ( f − ǫ ) Z f − ǫl df p ( f ) Z f l df p ( f ) (59)which also decays as 1 /t at large times. An expression for the distribution for the last jump time T , can also bewritten down in an analogous manner and reads as P ( t ) = L t Z ul + ǫ df p ( f ) p ( f + ǫ ) Z f − ǫl df p ( f ) Z f l df p ( f ) ǫ → → L t G (0) (60)Clearly the distribution P B ( t ) = 2( P ( t ) + P ( t )) ∼ t − . Thus the probability distribution P ( t ) = P A ( t ) + P B ( t ) obeysthe inverse square law for any block fitness distribution. V. CONCLUSIONS
In this article, we studied a deterministic model [16] describing the evolution of a population of self-replicatingsequences on a class of strongly correlated fitness landscapes with several fitness peaks [7]. The broad questions4addressed in this paper have been studied on completely uncorrelated fitness landscapes in previous works [16–18].Here we are interested in finding how the various evolutionary properties are affected when the sequence fitnesses arecorrelated.We are primarily interested in the evolutionary dynamics and in particular, the properties of jumps that occur inthe population fitness when the most populated sequence changes. As discussed in Sec. II, the largest fitness at aconstant Hamming distance from the initial sequence only need to be considered for this purpose. This led us toconsider the problem of the extreme statistics of correlated random variables [20, 29] which has been much less studiedthan its uncorrelated counterpart. We found that the extreme value distribution is not of the Gumbel form which isobtained when the random variables are i.i.d. and their distribution decays faster than a power law. In fact, we expectthat the universal scaling distributions which depend only on the nature of the tail of the underlying distribution donot exist for such correlated random variables as the number of independent variables namely the block fitnesses istoo small.As the minimum requirement of a sequence to qualify as a leader is that it must be a record, we also studied severalrecord properties of correlated variables. Recently the statistics of record events when the number of observationsadded at each time step increases either deterministically [30] or stochastically [31] have been studied. The recordsdefined in the shell model are an example of the former category as the number of observations changes as (cid:0) LD (cid:1) with D . It was shown that the probability for a record to occur in a shell with D mutations is not a continuous functionunlike the record distributions for independent random variables [17]; however the universality property that thedistribution is independent of block fitness distribution continues to hold. The average number of records was foundto increase linearly with L as in the maximally uncorrelated case but with the prefactor given by (1 − ln 2) ≈ . L dependence of the average number of jumps was seen to depend on the classof the fitness distribution p ( f ). For p ( f ) decaying faster than a power law, the average number of jumps increasedas √ L [16, 17]. In contrast, here the average number of jumps was shown to be independent of L for any choice ofblock fitness distribution p ( f ) although the value of the constant was found to be nonuniversal. These results suggestthat for block fitness distributions decaying faster than a power law, the average number of records increases but theaverage number of jumps decreases with increasing correlations. It is also interesting to see how the average numberof jumps change when the block fitness depends not only on the block configuration but also on its position in thesequence [7]. The result of our numerical simulations for this general model shows that the average number of jumpsincreases linearly with the number of blocks. However the prefactor is given by the average number of jumps obtainedin the locus-independent block fitness model namely (3 /
4) + p (see Fig. 6). This suggests that the different blocksbehave independently in the locus-dependent block fitness model; a detailed understanding of this model is beyondthe scope of this article and will be presented elsewhere.The temporal distribution for the last jump to occur at time t obeys t − law for infinite (and finite) populationsevolving on uncorrelated fitness landscapes [16–18]. Here we have shown that on a class of strongly correlated fitnesslandscapes, the same law is obeyed. The origin of this power law can be understood using a simple scaling argumentwhen the fitness variables are independent variables [16] but it is not obvious at the outset that such an argumentcan be used here since the sequence fitnesses are correlated. But it turns out that the jump time involves the i.i.d.block fitnesses and therefore t − law is obtained here as well.We close this article by a discussion of the deterministically evolving populations of infinite size studied here vis-a-vis finite populations that are subject to stochastic fluctuations on multi-peaked fitness landscapes. As discussedin Sec. II, the basic difference between a finite and an infinite population is that while the former has a finitemutational spread in the sequence space, all the mutants are available at all times in the deterministic case. As aconsequence, on rugged fitness landscapes, a finite population can get trapped at a local optimum from which it canescape by tunneling through a fitness valley [19]. In fact at late times, most of the population passes exclusivelythrough the local fitness peaks and thus such sequences are the most populated ones when the population size isfinite. In contrast, as the entire sequence space is occupied for infinite population, a transition to a higher fitnesspeak takes place by overtaking the less fitter populations as explained in Sec. II. Thus the underlying mechanismfor the punctuated increase of fitness on fitness landscapes with multiple peaks is different in the two situations[18]. Moreover the most populated sequence involved in the jump event is not necessarily a local maximum (for anycorrelation) for infinite populations. To see this, consider the fittest sequence with fitness w ( max ) ( D ) at Hammingdistance D from the initial sequence σ (0) . Barring the initial sequence, all the one-mutant neighbors of sequence withfitness w ( max ) (1) are at Hamming distance two from the initial sequence. Consider the scenario when the sequencewith fitness w ( max ) (2) is a nearest neighbor of sequence with fitness w ( max ) (1). Then the fittest sequence at distanceunity from the initial sequence can be a jump if at least w ( max ) (1) > w ( σ (0) ) and the minimum intersection timecondition ( w ( max ) (1) − w ( σ (0) )) − < w ( max ) (2) − w ( σ (0) )) − is obeyed. Clearly the latter condition rewritten as w ( max ) (2) − w ( max ) (1) < w ( max ) (1) − w ( σ (0) ) can be satisfied even when w ( max ) (1) is not a local maximum. Thus thenumber of jump events are not related to the number of local optima for an infinite population.5 FIG. 6: (Color online) Average number of jumps as a function of B for the block model with locus-dependent block fitnesschosen from exponential distribution. The line has a slope given by 3 / p = 0 . Acknowledgments: The authors thank Gayatri Das, Joachim Krug and Sanjib Sabhapandit for useful discussions. [1] S. Gavrilets,
Fitness Landscapes and the Origin of Species (Princeton University Press, New Jersey, 2004).[2] H. Orr, Nat. Rev. Genet. , 531 (2009).[3] P. Romero and F. Arnold, Nat. Rev. Mol. Cell Biol. , 866 (2009).[4] R. Korona, C. H. Nakatsu, L. J. Forney, and R. E. Lenski, Proc. Natl. Acad. Sci. USA , 9037 (1994).[5] C. L. Burch and L. Chao, Nature , 625 (2000).[6] G. Fernandez, B. Clotet, and M. Martinez, J. Virol. , 2485 (2007).[7] A. Perelson and C. Macken, Proc. Natl. Acad. Sci. USA , 9657 (1995).[8] F. Poelwijk, D. Kivet, D. Weinreich, and S. Tans, Nature , 383 (2007).[9] M. Lunzer, S. P. Miller, R. Felsheim, and A. M. Dean, Science , 499 (2005).[10] D. M. Weinreich, N. F. Delaney, M. A. DePristo, and D. L. Hartl, Science , 111 (2006).[11] J. de Visser, S.-C. Park, and J. Krug, Am. Nat. , S15 (2009).[12] S. A. Kauffman, The Origins of Order (Oxford University Press, New York, 1993).[13] T. Aita, H. Uchiyama, T. Inaoka, M. Nakajima, T. Kokubo, and Y. Husimi, Biopolymers , 64 (2000).[14] C. Carneiro and D. Hartl, Proc. Natl. Acad. Sci. USA , 1747 (2010).[15] G. Woodcock and P. G. Higgs, J. theor. Biol. , 61 (1996).[16] J. Krug and C. Karl, Physica A , 137 (2003).[17] K. Jain and J. Krug, J. Stat. Mech.: Theor. Exp., P04008(2005).[18] K. Jain, Phys. Rev. E , 031922 (2007).[19] K. Jain and J. Krug, Genetics , 1275 (2007).[20] H. David and H. Nagaraja, Order Statistics (Wiley, New York, 2003).[21] B. C. Arnold, N. Balakrishnan, and H. N. Nagaraja,
Records (Wiley-Interscience, 1998).[22] D. B. Saakian, O. Rozanova and A. Akmetzhanov, Phys. Rev. E , 041908 (2008).[23] D. B. Saakian and J. F. Fontanari, Phys. Rev. E , 041903 (2009).[24] M. Eigen, Naturwissenchaften , 465 (1971).[25] K. Jain and J. Krug, in Structural Approaches to Sequence Evolution: Molecules, Networks and Populations , edited byU. Bastolla, M. Porto, H. Roman, and M. Vendruscolo (Springer, Berlin, 2007) pp. 299–340, arXiv:q-bio.PE/0508008.[26] I. Bena and S. Majumdar, Phys. Rev. E , 051103 (2007).[27] F. Krzakala and O. Martin, Eur. Phys. J. B , 199 (2002).[28] G. Das, Dynamical properties of a quasispecies model on correlated fitness landscapes (M.S. thesis, JNCASR, Bangalore,2010).[29] K. Jain, A. Dasgupta, and G. Das, J. Stat. Mech., L10001(2009).[30] J. Krug, J. Stat. Mech.:Theor. Exp., P07001(2007).[31] I. Eliazar and J. Klafter, Phys. Rev. E80