A General Extreme Value Theory Model for the Adaptation of DNA Sequences Under Strong Selection and Weak Mutation

Abstract

Recent theoretical studies of the adaptation of DNA sequences assume that the distribution of fitness effects among new beneficial mutations is exponential. This has been justified by using extreme value theory and, in particular, by assuming that the distribution of fitnesses belongs to the Gumbel domain of attraction. However, extreme value theory shows that two other domains of attraction are also possible: the Fréchet and Weibull domains. Distributions in the Fréchet domain have right tails that are heavier than exponential, while distributions in the Weibull domain have right tails that are truncated. To explore the consequences of relaxing the Gumbel assumption, we generalize previous adaptation theory to allow all three domains. We find that many of the previously derived Gumbel-based predictions about the first step of adaptation are fairly robust for some moderate forms of right tails in the Weibull and Fréchet domains, but significant departures are possible, especially for predictions concerning multiple steps in adaptation.

ADAPTATION occurs at the level of DNA sequences. Recent efforts in the theory of adaptation have, therefore, focused on patterns that might characterize the movement of a population through DNA sequence space when evolution is driven by natural selection. Building on Gillespie's seminal work (Gillespie 1983, 1984, 1991), Orr (2002, 2003a, 2005) and Rokytaet al. (2006) derived a number of predictions about the adaptation of DNA sequences. Their models consider the evolution of individual genes and yield predictions that are generally independent of most biological details. In fact, predictions typically depend only on the number of beneficial mutations available to a starting wild-type sequence.

The key assumption underlying this theory is that the distribution of fitness, while unknown, belongs to a large class of probability distributions known as the Gumbel domain of attraction. The Gumbel domain is broad and includes most familiar probability distributions, including the normal, exponential, logistic, and gamma. Extreme value theory shows that the right tails of such distributions have similar behavior. In particular, values in excess of a high threshold are approximately exponentially distributed. In the case of adaptation, then, the distribution of effects among beneficial mutations should be nearly exponential so long as the (unknown) fitness distribution belongs to the Gumbel domain and the wild-type sequence is highly fit (i.e., represents a high threshold).

Although the Gumbel assumption represents a natural starting point for the theory of adaptation, it also represents a possible limitation: extreme value theory shows that non-Gumbel tail behavior can also occur. In particular, extreme value theory shows that a distribution, so long as it meets minimal criteria, can belong to one of three domains: Gumbel, Fréchet, or Weibull. The Fréchet domain loosely corresponds to distributions with heavy tails (heavier than exponential), while the Weibull domain loosely corresponds to distributions with right truncated tails (although the Gumbel domain can also include some right truncated distributions). Although previous workers (e.g., Orr 2006) suggested that the Fréchet and Weibull domains are less natural biologically than the Gumbel, these arguments are far from conclusive. It is therefore important to consider adaptation through DNA sequence space when the Gumbel assumption is relaxed and fitness distributions belong to the Fréchet and Weibull domains. Here we study this problem.

We begin by considering adaptation under extreme forms of the Fréchet and Weibull domains. Though clearly biologically unrealistic, these forms provide valuable intuitions about the more realistic cases that follow. We then derive generalized forms of the results of Orr (2002, 2005) and Rokytaet al. (2006). Most of these general results require introducing only one new parameter, denoted by κ, into the models. By altering the value of κ, we can tune whether the right tail of a fitness distribution behaves like the tail of a distribution belonging to the Gumbel (κ = 0), the Fréchet (κ > 0), or the Weibull (κ < 0) domain.

We find that some of the previous results derived for the Gumbel domain are robust to modest departures into the other domains. Indeed if κ ranges between

and

\(\frac{1}{2}\)

⁠, our results are qualitatively similar to those of Orr (2002), who considered κ = 0. We do, though, find qualitatively different patterns of adaptation in the Weibull domain when using parameter values suggested by the empirical work of Rokytaet al. (2008). For some types of Fréchet tails (κ ≥

\(\frac{1}{2}\)

⁠), we obtain unstable results that are difficult to interpret. However, for κ <

\(\frac{1}{2},\)

we derive a simple set of results in which findings from Orr (2002, 2005) and Rokytaet al. (2006) represent special cases.

THE MODEL

General assumptions:

We consider a scenario that is identical to that considered by Orr (2002)—a haploid population of DNA sequences of length L adapting under strong-selection weak-mutation (SSWM) conditions; i.e.,

and Nμ < 1, where N is the population size, s is a typical selection coefficient, and μ is the per site per generation mutation rate. Under these conditions, a population adapts through a series of selective sweeps, and the rate of adaptation is limited by the appearance of new beneficial mutations. Also, double mutants occur at too low a rate to influence adaptation; thus only the 3L single-mutant neighboring sequences to the currently fixed sequence need be considered. Deleterious and neutral mutations as well as recombination events are ignored.

The process begins with a population consisting of a wild-type sequence displaced from its fitness optimum, perhaps by an environmental change. If some number of the 3L single-mutant neighbors of this wild type are more fit, then adaptation will occur. If all of the 3L possible mutations and the wild type are ranked by their fitnesses such that the fittest has rank 1, the second fittest has rank 2, etc., then the wild type will have some rank i. Thus i − 1 beneficial mutations are available. One of these i − 1 beneficial mutations will eventually fix in the population, increasing fitness. At this point, the process begins anew.

The generalized Pareto distribution:

To build a model of adaptation, one must posit a distribution for the effect sizes of beneficial mutations. It is safe to assume that the vast majority of mutations available to a wild-type sequence decrease fitness; consequently, beneficial mutations are rare. The classical approach to studying rare events is to consider the upper-order statistics, for example, the maxima of large samples. Extreme value theory shows that the maxima from very large samples converge on one of three possible extreme value distributions (Pickands 1975). Each distribution corresponds to one of the domains of attraction discussed earlier.

Importantly, the distributions for the maximum under all three domains of attraction can be described by a single generalized extreme value distribution (see Embrechtset al. 1997). Results under the special case of the Gumbel domain were used by Gillespie (1991, 1983, 1984) and Orr (2002, 2003a) in their implementations of the model we consider. Specifically, they relied on the distributions of the spacings between the order statistics that have certain convenient properties under the Gumbel domain (in particular, the spacings are independent exponential random variables). Unfortunately, the spacings for distributions in the other two domains lack such simple forms, necessitating an alternative approach.

Our approach is based on the “peaks over threshold” formulation of extreme value theory. This approach considers the distribution of values greater than some high threshold. This approach is natural in the context of adaptive evolution, as we are interested in those mutations that have fitnesses greater than the wild type. If the fitness of the wild type is used as the threshold, the tail of the “excess” distribution describes the distribution of beneficial fitness effects.

RESULTS

Adaptation as κ →−∞:

We first consider the limiting form of the transition probabilities, based on Equation 7, for the Weibull domain as κ → −∞. If S ∼ GPD(κ, τ), then as κ decreases, S becomes less variable. Ultimately, Var(S) → 0 as κ → −∞, as can be seen from Equation 6. Despite this, we will see that P_ij(S) converges to the discrete uniform distribution on the integers 1, 2, … , i − 1, which is the most variable of all of the possible distributions for P_ij(S).

In words, as κ decreases adaptation becomes more and more random in terms of the identity of the beneficial mutation fixed, although it becomes less variable in terms of the fitness effects [since Var(S) → 0]. In the limit, each beneficial mutation has the same probability of being “grabbed” by natural selection. Interestingly, this “equally often” or “random” move rule has been studied extensively in various models of adaptation, including NK models (Kauffman and Levin 1987; Kauffman 1993) and the block model (Perelson and Macken 1995).

Adaptation as κ → ∞:

We have thus found that the two extreme (limiting) forms of adaptation correspond to two well-characterized move rules in adaptation: when κ → −∞ (Weibull), evolution gives no preference to any of the more fit alleles and uses a random (i.e., equally often) move rule, whereas when κ → ∞ (Fréchet), evolution uses a perfect move rule. Assuming that Equation 7 holds, adaptation can be viewed as a continuum between random and perfect adaptation with the location along the continuum specified by the shape parameter κ of the GPD (Figure 2). As pointed out by Orr (2002), adaptation under the Gumbel domain (κ = 0) falls exactly between these two extremes (see Figure 2).

Mean transition probabilities in the general case:

The jth largest fitness effect S_j will have finite mean only if κ < 1 and finite variance for κ <

However,

\(P_{ij}(\mathbf{\mathrm{S}}){=}S_{j}/({\sum}_{k{=}1}^{i{-}1}S_{k})\)

is a uniformly bounded random variable, since 0 ≤ P_ij(S) ≤ 1. Therefore the expected value of P_ij(S) and in fact all of its moments will exist regardless of κ. However, Equation 16 is a valid probability distribution only for κ < 1. For

\(\frac{1}{2}\)

≤ κ < 1, the above approximation is extremely poor (Figure 3), as Var(S) is infinite and the law of large numbers approximation fails. For 0 < κ <

\(\frac{1}{2}\)

⁠, the approximation is better, though still not good. Equation 16 is plotted in Figure 2 for i = 10 and −3 ≤ κ ≤

\(1/2\)

⁠.

While Equation 16 provides a simple formula for the mean transition probabilities, Equation 14 provides the key insight into their behavior. Equation 14 shows, for example, that when κ = −1, there is a decreasing linear relationship between E(P_ij(S)) and U. If κ = −2, there is a quadratic relationship, and when κ =

⁠, there is a one minus a square-root relationship, and when κ =

\(\frac{1}{2}\)

⁠, there is a reciprocal of a square root.

Expected rank after a single step:

One complete step of adaptation:

After a population substitutes a beneficial mutation through natural selection, it acquires a new set of neighboring sequences. Thus the fitness rank of the new wild-type sequence can change depending on the number of new neighbors that have higher fitnesses. The process through which a new wild-type sequence acquires new neighbors is described in detail by Rokytaet al. (2006) in the case in which the fitness distribution belongs to the Gumbel domain of attraction. To illustrate the idea, consider a sequence of length L = 1000. Initially, the wild-type sequence has some rank, say i = 10; i.e., if the fitnesses of the wild-type sequence and its 3000 single-mutation neighbors are listed in increasing size, the wild type's fitness is 10th from the top. Now assume natural selection moves the population up the list to the sequence having rank Y = 5. The new sequence had rank 5 among the 3001 sequences including the original wild type and its single-mutation neighbors, but now those sequences, except for the original wild type, differ from the new wild type by two mutations and are thus no longer accessible to natural selection. The new wild-type sequence now has a new set of 2999 accessible sequences plus the original wild type and thus, potentially, a new fitness rank. We can imagine that these 2999 sequences are all assigned fitnesses from the same distribution as for the original set and are again listed in increasing order.

Probability of adaptation—exact vs. approximate:

In the above sections we assumed that the transition probabilities are given by

⁠. However, this formula is valid only when s is small and the probability of fixation, denoted by Π(s), is proportional to s. In this section we investigate some of the implications on the theory when the small s approximation does not hold. We consider Π(s) = 1 − e^−2s (Kimura 1983) to be a more valid approximation.

Under the Weibull domain,

⁠, which is only approximately GPD when ln(1 − x) ≈ −x. This distribution is still in the Weibull domain of attraction. However, Martin and Lenormand (2008) showed that to have extreme observations under a right-bounded distribution, the wild type must have fitness fairly close the boundary, which dictates that Π(s) ≈ 2s. Under the Weibull domain it appears that extreme value theory is valid whenever the approximation Π(s) ≈ 2s is valid.

Average fitness improvement for a single step:

Why do the results in earlier sections match Orr (2002) exactly when κ = 0, while the results of this section only approximate those of Orr (2002) at κ = 0? The distribution of Y is derived using the approximation

⁠. The validity of this approximation (used in Equation 14) is justified by the law of large numbers and is rigorously established by Equation C2 in  appendix c with r = 1. However, it is a remarkable fact that when κ = 0, this is not an approximation, but rather an exact result. This implies that the mean and the variance of the rank distribution E(Y) and Var(Y) agree with the results derived in Orr (2002) when κ = 0. However, in this section we use the approximations

\((1/(i{-}1))E(S_{j}^{2}/{\bar{S}}){\approx}(1/(i{-}1))E(S_{j}^{2})/E({\bar{S}})\)

and

\((1/(i{-}1))E(S_{j}^{3}/{\bar{S}}){\approx}(1/(i{-}1))E(S_{j}^{3})/E({\bar{S}})\)

⁠, which also follow from Equation C2 in  appendix c with r = 2 and r = 3, but are no longer exact results when κ = 0. However, as noted after Equation C2 in  appendix c, the above approximations are quite accurate even for small sample sizes. As stated previously, the accuracy of the above approximations is determined by Var(S), and since Var(S) decreases as κ decreases (see Equation 6), the approximation improves as κ becomes small.

The probability of parallel evolution:

Mean fitness improvement after multiple steps:

We now compare the average fitness improvement over multiple steps for the Gumbel and Weibull domains. Similar calculations assuming the Fréchet domain are generally not possible, as most of the expectations do not exist (see discussion). In general, for single-step results, the model yields results that change continuously with κ (Figures 2–4). The differences between the domains become much more pronounced when considering multiple steps in adaptation.

If we denote the absolute fitness of allele j by W_j and view it as a random variable from some unknown distribution, then we define the fitness effect of allele j by W_j − W_i, where W_i is the fitness of the wild-type sequence. Selection coefficients are given by (W_j − W_i)/W_i. Recall that in all of the previous results the fitness of the wild type was assumed to be fixed and nonrandom (W_i = w_i). Therefore, any result previously derived for selection coefficients involving the GPD(κ, τ) is valid for fitness effects, provided τ is replaced by w_iτ. Results that do not depend on the scale parameter τ will be the same under both conventions. However, in this section we view the fitness of the wild type as a random variable, which changes over multiple steps. It is inconvenient to rescale by W_i at each step, and therefore it is more appropriate to deal with fitness improvement. In what follows we admit to a slight abuse of notation. Throughout we use the same notation for fitness effects that was previously used for selection coefficients. We define the fitness effect of the jth fittest allele by S_j ≡ W_j − W_i. If we define δ = E(W − W_i | W_i, W > W_i), then the fitness improvement results above (Equations 33 and 35) still hold.

DISCUSSION

Previous adaptation theory has lived with a certain tension. On the one hand, the theory's most important conclusion is that adaptation is characterized by robust patterns: we can say a good deal about the evolution of DNA sequences independent of knowledge of many biological details (e.g., fitness distribution, gene size, etc.). On the other hand, these conclusions have followed from one major assumption: that the distribution of fitness, while unknown in detail, belongs to the Gumbel domain of attraction. As noted earlier, the Gumbel domain is broad and includes most familiar distributions (normal, exponential, gamma, logistic, etc.). The Gumbel domain also represents the focus of most classical extreme value theory and allows many elegant, closed-form solutions, e.g., to the distribution of spacings between upper-order statistics. For these reasons and others, the Gumbel domain of attraction provides a natural starting point for a theory of adaptation based on extreme value theory. Not surprisingly, Gillespie (1983, 1984, 1991), Orr (2002, 2003b, 2005), and Rokytaet al. (2006) all assumed that the distribution of fitness, while unknown, belongs to the Gumbel domain.

Extreme value theory shows, however, that three domains of attraction are possible for distributions meeting minimal conditions: Gumbel, Fréchet, and Weibull. Although the Fréchet and Weibull domains might first seem less biological than the Gumbel—distributions belonging to the Fréchet domain, in particular, are somewhat exotic, lacking all or higher moments—this impression is obviously far from convincing. Moreover, history suggests caution about casual inferences about the naturalness of various distributions: while early work in mathematical finance often employed Gumbel statistics, data from equity markets now reveal heavy tail distributions of stock returns; i.e., these distributions likely belong to the Fréchet domain of attraction. Analogously, early theoretical work from Gillespie (1984) suggests Gumbel statistics, but empirical data in Rokytaet al. (2008) reject the Gumbel domain in favor of a uniform distribution (κ ≈ −1) in the Weibull domain. It is important, therefore, to study the consequences for the theory of adaptation of departures into the Fréchet and Weibull domains of attraction.

We have studied this problem here. Our approach took advantage of a general version of extreme value theory that encompasses all three domains of attraction. In particular, we analyzed adaptation via the GPD, which provides the distribution of values above a high threshold. (In the case of adaptation, the high threshold represents the fitness of the wild-type allele, while “excesses” above this threshold represent the fitnesses of rare beneficial mutations.) By varying a single parameter, κ, one can “tune” the GPD for the cases of Gumbel, Fréchet, or Weibull domains of attraction.

Our findings show that Gumbel-based results in the genetics of adaptation are often robust to modest departures into the Weibull and, especially, the Fréchet domains of attraction. But in some cases (e.g., κ <

or extended adaptive walks through sequence space), real departures from Gumbel results arise. Because our findings for the Fréchet and Weibull domains also differ from each other, it is worth briefly summarizing our results for these domains separately, emphasizing when they are and are not similar to those for the Gumbel domain.

The Fréchet domain:

Finally we consider extremely heavy-tailed distributions, i.e., κ > 1. A sample drawn from such a distribution will depend heavily on the scale parameter τ. If the scale parameter is moderate, then the selection coefficient for the fittest allele will be unreasonably large. If the scale parameter is small, then the fitness for the higher ranked alleles will be undetectably small. Therefore, it seems biologically impossible to produce a scenario in which one could draw reasonable inferences under the Fréchet domain that depart significantly from the Gumbel domain. This is one of our most important conclusions.

The Weibull domain:

As pointed out in Rokytaet al. (2008), a right-truncated distribution is not at all unreasonable. For example, even if the phenotypic effects of mutations show no apparent upper bound, the translation of phenotype into fitness might yield a truncated distribution. A pattern of diminishing returns (i.e., a concave mapping of phenotype onto fitness), such as is commonly seen for biochemical reactions (Hartlet al. 1985), could produce an apparent right truncation point for fitness effects, regardless of the underlying distribution of phenotypic effects. Recent theoretical work by Martin and Lenormand (2008) shows that selection to a phenotypic optimum in a “geometric” model of evolution can lead to a Weibull-type distribution of mutational effects, although a Gumbel-type distribution remains a good approximation in the case of many independent characters. There is also empirical evidence from two different viruses to suggest that the Weibull domain may be appropriate in some circumstances (Rokytaet al. 2008).

While the pattern of variability in samples may differ substantially between the Gumbel and Weibull domains, predictions for the first step in an adaptive walk developed by Orr (2002) for the Gumbel domain easily generalize to the Weibull domain (Figure 2). In general, the results describing an average gain after one step are smaller under the Weibull domain than under the Gumbel domain (Figure 3). For example, κ = −1 (as suggested by Rokytaet al. 2008) represents a uniform distribution and κ = 0 represents the exponential distribution. By comparing Equation 21 to Equation 22, we see that the first step in adaptation moves the population on average two-thirds toward the one-step optimum under the uniform model vs. three-fourths under the exponential. Also, using Equation 37 we see that parallel evolution is 1.5 times more likely under the exponential than under the uniform, although both models predict that the probability of parallel evolution is inversely proportional to the rank of the wild type. The mean fitness improvement after one step in adaptation under the exponential distribution is again 1.5 times larger than under the uniform (see Equation 41).

Predictions under the two models differ even more when considering multiple steps. Under the Weibull domain, the mean of the fitness distribution decreases geometrically in the number of steps; i.e., after n steps the mean of the fitness distribution is δ_n = δ/(1 − 2κ)ⁿ. The memoryless property of the exponential distribution, on the other hand, guarantees that the mean of the fitness distribution remains the same after each step.

Much of the present theory of adaptation rests on the assumption of SSWM conditions. However, the geometrically decreasing change in mean fitness suggests that even if SSWM conditions hold during the first step of adaptation, they may fail to do so during later steps. This represents a marked change from results seen under the Gumbel domain.

APPENDIX A: ORDER STATISTICS OF THE UNIFORM DISTRIBUTION

APPENDIX B: MOMENTS OF THE GENERALIZED PARETO DISTRIBUTION

APPENDIX C: VALIDITY OF THE APPROXIMATION USED TO CALCULATE THE MEAN TRANSITION PROBABILITIES

APPENDIX D: DERIVATIONS OF EQUATIONS 22 AND 24

APPENDIX E: THE ORDER STATISTICS FROM THE FRÉCHET DOMAIN

Footnotes

Footnotes

Acknowledgement

We thank Holly A. Wichman for many useful discussions and support on this project. This work was supported by grants from the National Institutes of Health (NIH) to P.J. (R01GM076040) and H.A.O. (GM51932). C.J.B. was supported in part by NIH grant P20 RR16448 and a grant from the National Science Foundation to P.J. (DEB-0515738). D.R.R. was supported in part by NIH grant P20 RR16454.

References