## Abstract

Much of the current theory of adaptation is based on Gillespie’s mutational landscape model (MLM), which assumes that the fitness values of genotypes linked by single mutational steps are independent random variables. On the other hand, a growing body of empirical evidence shows that real fitness landscapes, while possessing a considerable amount of ruggedness, are smoother than predicted by the MLM. In the present article we propose and analyze a simple fitness landscape model with tunable ruggedness based on the rough Mount Fuji (RMF) model originally introduced by Aita *et al.* in the context of protein evolution. We provide a comprehensive collection of results pertaining to the topographical structure of RMF landscapes, including explicit formulas for the expected number of local fitness maxima, the location of the global peak, and the fitness correlation function. The statistics of single and multiple adaptive steps on the RMF landscape are explored mainly through simulations, and the results are compared to the known behavior in the MLM model. Finally, we show that the RMF model can explain the large number of second-step mutations observed on a highly fit first-step background in a recent evolution experiment with a microvirid bacteriophage.

THE genetic adaptation of an asexual population to a novel environment is governed by the number and fitness effects of available beneficial mutations, their epistatic interactions, and the rate at which they are supplied (Sniegowski and Gerrish 2010). Despite the inherent complexity of this process, recent theoretical work has identified several robust statistical patterns of adaptive evolution (Orr 2005a,b). Most of these predictions were derived in the framework of Gillespie’s mutational landscape model (MLM), which is based on three key assumptions (Gillespie 1983, 1984, 1991; Orr 2002). First, selection is strong enough to prevent the fixation of deleterious mutations and mutation is sufficiently weak such that mutations emerge and fix one at a time [the strong selection/weak mutation (SSWM) regime]. Second, wild-type fitness is high, which allows one to describe the statistics of beneficial mutations using extreme value theory (EVT). Third, the fitness values of new mutants are uncorrelated with the fitness of the ancestor from which they arise. This last assumption implies that the fitness landscape underlying the adaptive process is maximally rugged with many local maxima and minima (Kauffman and Levin 1987; Kauffman 1993; Jain and Krug 2007), a limiting situation that is often referred to as the house of cards (HoC) landscape (Kingman 1978). Thus, the MLM is concerned with a population evolving in a HoC landscape under SSWM dynamics, starting from an initial state of high fitness.

The validity of the SSWM assumption depends primarily on the population size *N*. Denoting the mutation rate by *u* and the typical selection strength by *s*, the criterion for the SSWM regime reads *Nu* ≪ 1 ≪ *Ns*, which can always be satisfied by a suitable choice of *N* provided *u* ≪ *s*, as is usually the case. On the other hand, whether the other two assumptions underlying the MLM are realistic is an empirical question that has been addressed in a number of experimental studies of microbial evolution. Investigations aimed at determining the distribution of effect sizes of beneficial mutations have generally found support for the EVT hypothesis (Orr 2003a, 2010; Joyce *et al.* 2008), and examples for all three EVT universality classes have been reported in the literature (Rokyta *et al.* 2005, 2008; Kassen and Bataillon 2006; MacLean and Buckling 2009; Schenk *et al.* 2012; Bank *et al.* 2014; Foll *et al.* 2014). At the same time, however, it has become increasingly clear that the extreme assumption of uncorrelated fitness values between genotypes connected by single mutational steps cannot be upheld in the face of empirical evidence.

Indications for the presence of correlations in real fitness landscapes derive from two types of experimental studies. In one approach, a subset of the fitness landscape is explicitly generated by constructing genotypes containing all combinations of a small group of mutations chosen for either individual or collective effects and measuring their fitness or some proxy thereof (Weinreich *et al.* 2006, 2013; Lozovsky *et al.* 2009; Franke *et al.* 2011; Schenk *et al.* 2013; Szendro *et al.* 2013b; De Visser and Krug 2014). Although the topographic properties of the resulting landscapes vary over a broad range, in most cases they display a degree of ruggedness that is intermediate between a smooth, additive landscape and the maximally rugged landscape assumed by the MLM (Szendro *et al.* 2013b; De Visser and Krug 2014). In a second approach, properties of the underlying landscape are inferred from the observed dynamics of adaptation as manifested, for example, in the trajectories of fitness increase (Kryazhimskiy *et al.* 2009) or the number of substitutions in an adaptive walk (Gifford *et al.* 2011; Schoustra *et al.* 2012). Of particular interest in the present context is the recent study of Miller *et al.* (2011) on the microvirid bacteriophage ID11, where the MLM was tested by comparing the distribution of beneficial mutations from the wild type to the corresponding distribution after one step of adaptation. According to the MLM, the two distributions should be identical up to a rescaling, but this hypothesis was clearly refuted by the experiment.

The observation that most empirical fitness landscapes display an intermediate degree of ruggedness implies that there is a need for simple, mathematically tractable landscape models in which the ruggedness can be tuned. A frequently used model with tunable ruggedness is Kauffman’s “NK” landscape, where each of *L* binary loci interacts randomly with *K* other loci, and the interaction degree *K* serves to interpolate between the additive limit *K* = 0 and the HoC limit *K* = *L* − 1 (Kauffman and Weinberger 1989; Kauffman 1993; Ohta 1997; Welch and Waxman 2005; Aita 2008; Franke and Krug 2012; Østman *et al.* 2012). In the original definition of the NK model the letter “N” stands for the sequence length, which we denote by *L* in the present article. While this model has been shown to be capable of describing various features of empirical fitness landscapes (Hayashi *et al.* 2006; Rowe *et al.* 2010; Franke *et al.* 2011), its mathematical complexity is such that even rather elementary properties—for example, the mean number of local fitness maxima (Evans and Steinsaltz 2002; Durrett and Limic 2003; Limic and Pemantle 2004)—are difficult to derive in closed form (but see Perelson and Macken 1995, Orr 2006a, and Schmiegelt and Krug (2014) for a variant of the NK model that is simpler to analyze).

In the present article we therefore propose the rough Mount Fuji (RMF) model as an alternative description of fitness landscapes with tunable ruggedness. The model is a simplified version of the RMF fitness landscape originally introduced by Aita *et al.* (2000) in the context of protein evolution. In essence, the RMF model superimposes an additive fitness landscape and an uncorrelated random HoC landscape, and the ruggedness is tuned by changing the ratio of the additive selection coefficient to the standard deviation of the random fitness component (Aita and Husimi 2000; Franke *et al.* 2011; Szendro *et al.* 2013b). Below we derive simple, explicit formulas for various quantitative measures of the RMF topography such as the number and location of local fitness maxima and fitness correlations. Moreover, assuming SSWM conditions, we show how the adaptation of a population on the RMF landscape can be efficiently simulated for realistic numbers of loci by locally generating the mutational neighborhood of the current genotype along the adaptive trajectory. Finally, as an example for the application of the RMF model to empirical fitness landscapes, we estimate the parameters of the fitness landscape of the microvirid bacteriophage ID11 studied by Miller *et al.* (2011) by matching the number of secondary beneficial mutations available after one adaptive step (the *number of exceedances*) predicted by the model to the experimentally observed value.

## Model

### Definition

Following a common practice in the description of empirical fitness landscapes, we represent genotypes by binary sequences *σ* = (*σ*_{1}, *σ*_{2}, … , *σ _{L}*) of fixed length

*L*, composed of elements taken from the set {0, 1} with

*σ*= 1 (

_{i}*σ*= 0) if a mutation is present (absent) at the

_{i}*i*th locus. The set of all binary sequences of length

*L*is known as the Hamming space. It is endowed with a natural distance measure, the

*Hamming distance*, defined as (1)which simply counts number of loci at which

*σ*and

*σ*′ differ. It is convenient to introduce the antipodal (or reversal) sequence of

*σ*through A sequence and its antipode are maximally distant from each other,

To introduce the rough Mount Fuji model we first choose a reference sequence *σ**, which represents the state of maximal fitness of the additive part of the fitness landscape. The fitness *F*(*σ*) of genotype *σ* is then defined through (2)where *c* > 0 is a constant parameter and the *η*’s are 2* ^{L}* independent and identically distributed (i.i.d.) random variables. Equation 2 describes an average decrease of fitness with increasing distance from

*σ** by an amount of

*c*per mutational step, superimposed by a random fitness variation. For

*c*= 0 the RMF model reduces to an uncorrelated HoC landscape, while for large

*c*it becomes essentially additive, as the random fitness component

*η*is then negligible compared to the mean fitness gradient. The competition between the additive and random contributions is governed by the parameter (3)defined as the ratio between the additive selection coefficient

*c*and the standard deviation of the random fitness component

*η*(Franke

*et al.*2011). With increasing

*θ*the landscape becomes less rugged.

It is important to note that the RMF landscape is anisotropic, in the sense that the mutational neighborhood of a sequence *σ* depends on its distance from *σ**. To be specific, we define the neighborhood *ν*(*σ*) of *σ* as the set *ν*(*σ*) = {*σ*′|*D*(*σ*, *σ*′) = 1} ∪ {*σ*}. Denoting the distance of *σ* from the reference sequence by *d* ≡ *D*(*σ*, *σ*^{∗}) > 0, the set *ν*(*σ*) is split into three parts: (1) the *downhill neighborhood* that consists of the *L* − *d* sequences at distance *d* + 1 from *σ**, which have an expected fitness disadvantage of *c* compared to *σ*; (2) *σ* itself; and (3) the *uphill neighborhood* that consists of the *d* sequences at distance *d* − 1 from *σ*, which have an expected fitness advantage of *c* compared to *σ*.

This decomposition implies that, in contrast to the MLM, the fitness values of the mutational neighbors of *σ* are not i.i.d. random variables. We will see in the following how this leads to new properties of the fitness landscape and of the adaptive process on that landscape.

### Fitness distribution and extreme value theory

To complete the definition of the RMF model we need to specify the statistics of the random fitness component *η* in terms of its probability distribution function *P*(*x*) ≡ ℙ(*η* < *x*) and the corresponding probability density Following the approach developed previously for the MLM, we invoke EVT to classify the fitness distribution according to its tail behavior. The underlying reasoning is that viable organisms must have high fitness in absolute terms, which implies that the beneficial mutations that drive adaptation reside in the tail of the (usually unknown) distribution of fitness effects of all possible mutations (Gillespie 1983, 1984; Orr 2002, 2003a, 2010). The theorems of EVT then show that, irrespective of the detailed form of the full distribution, the tail shape has to fall into one of three classes (De Haan and Ferreira 2006): the *Gumbel class* containing all distributions with unbounded support and a density vanishing faster than a power law, for example, the exponential, normal and gamma distributions; the *Fréchet class* containing all distributions with unbounded support and a density vanishing as a power law; and the *Weibull class* containing all distributions with a truncated upper tail.

The three classes are conveniently represented through the generalized Pareto distribution (GPD) defined by the distribution function (4)(Pickands 1975; Beisel *et al.* 2007; Joyce *et al.* 2008) with the extreme value index *κ*. For *κ* > 0 the support of *P _{κ}* is [0, ∞) and the distribution belongs to the Fréchet class, while for

*κ*< 0 the support is and the distribution belongs to the Weibull class. For

*κ*→ 0 the GPD reduces to an exponential, which is a representative of the Gumbel class.

In the initial formulation of the EVT approach it was argued that the Gumbel class is the most likely case to be realized biologically, and empirical support for this hypothesis was found in several studies (Rokyta *et al.* 2005; Kassen and Bataillon 2006; MacLean and Buckling 2009; Orr 2010). However, subsequently systems showing truncated (Weibull-type) fitness distributions were discovered (Rokyta *et al.* 2008), and recently several examples for heavy-tailed (Fréchet-type) fitness distributions have been reported (Schenk *et al.* 2012; Bank *et al.* 2014; Foll *et al.* 2014). A truncated fitness distribution arises naturally in the adaptation toward a single fitness peak, as assumed in Fisher’s geometric model, and analysis of this model shows that the distribution becomes Gumbel-like as the dimensionality of the underlying phenotype increases (Orr 2006b; Martin and Lenormand 2008). A similar mechanistic interpretation is not known for distributions falling into the Fréchet class, but the available empirical evidence suggests that such heavy-tailed distributions may generally be associated with situations of strong selection pressure, *e.g.*, in the adaptation of pathogens to new drugs (Schenk *et al.* 2012; Bank *et al.* 2014; Foll *et al.* 2014).

We will see below that several properties of the RMF fitness landscape take on a particularly simple form when the random fitness component is chosen from a particular representative of the Gumbel class, the *Gumbel distribution*, defined by (5)This distribution arises in EVT as the limit law of the maximum of i.i.d. random variables drawn from a distribution in the Gumbel class (De Haan and Ferreira 2006). It approaches an exponential for large positive values and vanishes very rapidly (as the exponential of an exponential) for negative values. The key property of interest here is the behavior of *P*_{G} under shifts, (6)For completeness we note that the mean of the Gumbel density *p*_{G}(*x*) is the Euler–Mascheroni constant *γ* ≈ 0.5772156649 … and its variance is We emphasize that our use of this distribution in the following sections is motivated primarily by its mathematical convenience.

### Comparison to empirical fitness landscapes

Several recent studies using the RMF model to quantify properties of empirical fitness landscapes provide some guidance as to what range of model parameters can be expected in applications to biological data. Franke *et al.* (2011) and Szendro *et al.* (2013a) analyzed an eight-locus fitness landscape composed of individually deleterious mutations in the filamentous fungus *Aspergillus niger*. Based on the statistical properties of selectively accessible mutational pathways and a direct fit to the data, respectively, and assuming a normal distribution for the random fitness component *η*, they obtained consistent estimates *θ* ≈ 0.25 and *θ* ≈ 0.21 for the parameter defined in Equation 3. In a metaanalysis of 10 different fitness data sets the *A. niger* landscape was found to be among the more rugged empirical landscapes, indicating that these estimates for *θ* are at the lower end of the range of biologically relevant values (Szendro *et al.* 2013b). However, we will see below that the properties of the RMF landscape are strongly affected by the EVT index *κ*, which implies that both *θ* and *κ* are generally required to characterize an empirical data set (see *The Number of Exceedances* for an example).

## Structure of The Fitness Landscape

In this section we present results concerning the main structural features of the RMF fitness landscape, in particular its local maxima and fitness correlations.

### Fitness maxima

Local fitness maxima play a key role in adaptation, as they present obstacles to an evolving population, and the number of maxima is a commonly used measure of landscape ruggedness. In the HoC model all genotype fitness values are independent and statistically equivalent. The probability that a given sequence is a local fitness maximum is then simply since each of the *L* + 1 sequences in the neighborhood is equally likely to have the largest fitness, and the expected number of maxima is (Kauffman and Levin 1987; Kauffman 1993). These expressions for the maximally rugged HoC model serve as a benchmark for the corresponding results for the RMF model that are presented in the following section. Detailed derivations are given in *Appendix A*.

#### Density of local maxima:

In the RMF model the probability that a given genotype is a local fitness maximum depends on its distance *d* to the reference sequence. When the random component is Gumbel distributed, this quantity can be computed exactly, with the result (7)The limits of this expression for small and large *c*, (8)correspond to the HoC model and the additive landscape with a single maximum at *d* = 0, respectively. Here the Kronecker symbol *δ _{x}*

_{,}

*is defined by (9)The behavior of (7) for intermediate values of*

_{y}*c*is illustrated in Figure 1A.

It is also of interest to consider the probability that the neighboring genotype of largest fitness for a sequence at distance *d* from the reference sequence is located in the uphill or downhill part of the neighborhood, denoted by and respectively. These probabilities determine the fate of a “greedy” adaptive walk that chooses the neighboring sequence of highest fitness at each step (Orr 2003b). They can be explicitly evaluated for the Gumbel distribution, with the result (10)Note that and In the absence of a fitness gradient (*c* = 0) the greedy walker is more likely to go uphill (downhill) for *d* > *L*/2 (*d* < *L*/2) because of the greater availability of neighboring sequences in that direction. For *c* > 0 the crossing point where moves toward the reference sequence and is generally located at see Figure 1B.

#### Total number of maxima:

To determine the expected number of local fitness maxima ℳ in the entire landscape, the distance-dependent probability has to be averaged over *d* with the appropriate weights, giving the number of genotypes at distance *d*, (11)Using the exact expression in Equation 7 for the Gumbel distribution, the sum in Equation 11 can be expressed in terms of a hypergeometric function; see Equation A12. To obtain a more tractable expression we note that for large *L* the binomial weights in Equation 11 become sharply peaked around such that (12)For fixed *L* and large *c* the expression in Equation 12 violates the obvious bound ℳ ≥ 1. A simple modification that cures this deficiency is (13)which provides a very good approximation to the exact number of maxima over the entire range of parameters; see Figure 2A. Interestingly, apart from a multiplicative constant the large *L* asymptotics of the number of maxima are the same as for the HoC model, This is in contrast to the corresponding results for Kauffman’s NK model, where (depending on how the number of interacting loci *K* is scaled with the sequence length *L*) different exponential and algebraic dependencies of the number of maxima on *L* can be found (Perelson and Macken 1995; Evans and Steinsaltz 2002; Durrett and Limic 2003; Limic and Pemantle 2004; Schmiegelt and Krug 2014).

An exact expression for the expected number of maxima is derived in *Appendix A* for exponentially distributed randomness. In that case the asymptotic behavior for large *L* is seen to be identical to that in Equation 12, and in fact this behavior arises whenever the tail of the distribution is exponential (see Figure A1A). For tails heavier than exponential it is shown that the asymptotics are exactly those of the HoC model, independent of *c*, which implies, remarkably, that the mean fitness gradient has no effect on the number of maxima when the genotype dimension is sufficiently large. This result applies in particular to distributions in the Fréchet class of EVT (Figure A1B). For Gumbel-class distributions with tails thinner than exponential the number of maxima is significantly reduced for any value of *c*, and the ratio of ℳ to the HoC value 2* ^{L}*/

*L*vanishes subexponentially in

*L*(see Equation A21).

The effect of the mean fitness gradient on the number of local maxima is most pronounced for distributions with bounded support belonging to the Weibull class of EVT. An exact expression for the total number of maxima is derived in *Appendix A* for the case of a uniform fitness distribution, a simple representative of this class, and the result is illustrated in Figure 2B. To leading order the number of maxima is proportional to (2 − *c*)* ^{L}*/

*L*, which interpolates smoothly between the HoC result for

*c*= 0 and the additive limit ℳ = 1 attained at

*c*= 1; for

*c*> 1 the increase in fitness gained in one step toward the reference sequence exceeds the support of the distribution of the random component, and the landscape becomes strictly monotonic in

*d*. For other distributions in the Weibull class with support on the unit interval the behavior is similar, to leading order (Equation A26). The behavior of the number of maxima for distributions with bounded support is thus reminiscent of the NK model at fixed

*K*, where it has been shown that with a

*K*-dependent constant with 1 <

*λ*< 2 for 0 <

_{K}*K*<

*L*− 1 (Durrett and Limic 2003; Limic and Pemantle 2004; Schmiegelt and Krug 2014).

#### Location of the global maximum:

We next ask where the global maximum is located. For large *c* it will be found close to the reference sequence *σ**, while in the HoC limit *c* = 0 it is equally probable to be located anywhere. Since most sequences lie near the distance *L*/2 from the reference sequence, also the global maximum is then most likely at In general, the probability for the global maximum to lie at Hamming distance *d* from the reference sequence is given by (14)where is the probability for some specific genotype at distance *d* to be the globally fittest state and the binomial coefficient accounts for the multiplicity of states at distance *d*. In *Appendix A* it is shown that for the Gumbel distribution (15)With this quantity at hand we proceed to calculate the mean distance of the global maximum from *σ**, (16)which interpolates smoothly between the two limiting cases discussed above, and the corresponding variance (17)These calculations could be extended to include the positions of suboptimal local maxima and thus to address the possible clustering of maxima discussed previously in the context of the NK model (Kauffman 1993), but we do not pursue this question here. For general distributions an expansion for small *c* shows that the shift in the position of the maximum from the HoC value *L*/2 is of order *cL*2^{−}* ^{κL}*, where

*κ*is the extreme value index defined in Equation 4. For distributions in the Weibull class (

*κ*< 0) this implies that minute values of

*c*∼ 2

*suffice to bring the global maximum close to the reference sequence with high probability.*

^{κL}### Fitness correlations

In addition to the number of local maxima, a commonly used measure for fitness landscape ruggedness is the decay of fitness correlations (Weinberger 1990; Stadler and Happel 1999). Here we consider the correlation function defined by (18)where angular brackets denote an average over sequence space as well as over the realizations of the random fitness component *η* in Equation 2, and 〈⋅〉* _{r}* denotes an average over pairs of sequences

*σ*,

*σ*′ with

*D*(

*σ*,

*σ*′) =

*r*. The normalization of the expression in Equation 18 ensures that

*C*(0) = 1. The derivation in

*Appendix B*shows that the correlation function depends on the underlying random fitness distribution only through its variance, and it is given by the expression (19)The correlation function is a superposition of a peak at

*r*= 0, which originates from the uncorrelated HoC component in Equation 2, and a linearly decaying piece that reflects the global fitness gradient. It is instructive to compare this result to the correlation function for the NK model, which reads (20)(Campos

*et al.*2002, 2003). This displays a linear decay, in the nonepistatic limit

*K*= 0. However, in contrast to Equation 20, which is nonnegative for any

*r*, Equation 19 becomes negative for (see also Neidhart

*et al.*2013 for further discussion of the relation between the two models).

## Adaptation on the RMF Landscape

In the previous section we studied properties that depend purely on the topography of the landscape. Now we shift the focus to the implications of the landscape structure for the evolutionary dynamics. More specifically, we consider the dynamics in the SSWM limit. Here, only one genotype is populated at any time. If a beneficial mutation occurs, either it is fixed in the entire population or the mutant goes extinct before another mutation arises. Hence, the population behaves as a single entity that performs an “adaptive walk” (AW) on the fitness landscape (Gillespie 1983).

The adaptive walk is a sequence of single adaptive steps. In each step, the population moves from the currently populated sequence to a neighboring one with a transition probability given by the fixation probability normalized by the fixation probabilities of all other available beneficial mutants. Following Gillespie (1983) and Orr (2002) we consider the rank-ordered fitness values *F _{j}* in the current mutational neighborhood, where the rank of the resident genotype is

*i*and beneficial mutants have ranks

*j*<

*i*. Since the fixation probability of a beneficial mutation in the SSWM regime is proportional to its selection coefficient, the probability for a transition from the current

*i*th fittest genotype to the

*j*th fittest mutant is (21)Based on this expression a number of results have been obtained for the adaptation dynamics on the uncorrelated HoC landscape (Orr 2002; Rokyta

*et al.*2006; Joyce

*et al.*2008). In the following we ask how these results are modified by the fitness gradient in the RMF model.

### A single step of adaptation

Before we turn to the full adaptive walk we consider a single step of adaptation, specifically the change in the rank of the resident genotype during an adaptive step. Using the transition probability in Equation 21, Orr (2002) calculated the expectation value and the variance of the rank *j* of the next populated sequence after an adaptive step starting from the genotype with rank *i*. For HoC landscapes with fitness values drawn from a Gumbel class distribution he obtained (22)These results were subsequently generalized to the other extreme value classes by Joyce *et al.* (2008), who found, for EVT index (23)For the approximation used to derive the expression for Var(*j*) breaks down because the fitness distribution ceases to have a second moment; similarly the expression for (*j*) breaks down for *κ* → 1. For further discussion of this case of extremely heavy-tailed fitness distributions we refer to Schenk *et al.* (2012).

To see to what extent these results are modified in the RMF landscape we refer to *Appendix C*, where it is shown that the distribution of fitness values in a mutational neighborhood of the RMF model with exponential randomness approximately remains a simple exponential that is shifted by a constant amount, depending on the distance *d* to the reference sequence as well as on the model parameters *c* and *L*. This implies that the statistics of the fitness *spacings F _{k}*

_{−1}−

*F*are approximately independent of

_{k}*c*and

*d*and take the same form as in the HoC model with exponential fitness distribution. In fact, the exponential nature of the fitness spacings in this case is guaranteed by a general theorem of order statistics (Smid and Stam 1975). Since the transition probability in Equation 21 can be written in terms of these spacings, it seems reasonable to conjecture that the properties of single adaptive steps in the RMF model are well approximated by the HoC results at least for moderate values of the fitness gradient

*c*.

To test this conjecture, simulations consisting of the following steps were carried out: (i) create a random RMF neighborhood, (ii) determine the current rank of the initial genotype, (iii) carry out an adaptive step according to Equation 21, and (iv) determine the new rank in the old neighborhood. A comparison of the RMF simulation results to Orr’s analytical formulas in Equation 22 is shown in Figure 3. Obviously, the HoC expressions are in good agreement with the simulation data even in cases where the slope *c* of the RMF landscape is comparable to the variance of the random fitness contributions (for the exponential distribution considered here the variance equals unity).

For other distributions the approximation in *Appendix C* is not directly applicable, but the results shown in Figure 4 indicate that the generalized HoC formulas in Equation 23 provide a good approximation to the RMF model for a range of EVT indexes around *κ* = 0. Thus we conclude, somewhat unexpectedly, that the statistical properties of single adaptive steps are not strongly affected by the fitness gradient in the RMF model.

### Adaptive walks

In the context of AWs, a property of interest is the *mean walk length* ℓ, which is the average number of steps performed until the process reaches a local fitness maximum and terminates. The mean walk length in the HoC landscape has been analyzed using various approaches (Gillespie 1983; Flyvbjerg and Lautrup 1992; Orr 2002; Jain and Seetharaman 2011; Jain 2011; Neidhart and Krug 2011; Seetharaman and Jain 2014). Because of the lack of isotropy in the RMF landscapes, analytical results for the walk length are much harder to derive for this model, and the results described in the following paragraphs were therefore obtained from simulations.

The simulation algorithm is analogous to that described above for the first step of adaptation, but now a new neighborhood is created after each adaptive step and the procedure is repeated until a local maximum has been found, that is, until the current genotype has rank 1 in its new neighborhood. Since the distribution of the fitness values in the new neighborhood depends on the distance to the reference genotype, the direction of the adaptive steps (uphill or downhill) has to be kept track of. Creating the neighborhoods “on the fly” during the adaptive walk implies that the memory of previously visited genotypes is lost beyond the second step. However, the error associated with this approximation is expected to be negligible for large *L* (Flyvbjerg and Lautrup 1992; Neidhart and Krug 2011) and it is the only feasible approach for simulating walks on landscapes with thousands of loci.

On HoC landscapes the mean walk length is determined primarily by the starting rank *r*, that is, the rank that the first populated state has in its initial neighborhood, and was shown in previous work to be proportional to log(*r*) for *r* ≪ *L*; see below for further discussion. On the other hand for a smooth landscape with only one maximum, the walk length equals of course the distance *d* of the starting point to this maximum. Due to its anisotropic structure, on the RMF landscape one expects the walk length to depend on both initial rank and initial distance to the reference sequence *σ**. Specifically, for small *c* (in the sense *θ* ≪ 1), the mean adaptive walk length should increase logarithmically in the starting rank and be approximately independent of the initial Hamming distance *d* from *σ**, while for *θ* ≫ 1, it should increase linearly in *d* and be approximately independent of the starting rank, since with high probability only a single maximum exists in the fitness landscape at *σ** or close to it. Analysis of a simplified version of the problem where the walks are assumed to start from the antipode of the reference sequence shows that the two regimes are separated by a sharp transition under certain conditions (Park *et al.* 2014).

Before discussing the process governed by the transition probability in Equation 21, we briefly consider the simpler case of a *greedy* adaptive walk in which the transition occurs deterministically to the available genotype of maximal fitness in each step. On HoC landscapes, the mean walk length for this process is known to be asymptotically constant for large *r* and *L*, approaching the universal limit ℓ = *e* − 1 ≈ 1.72 (Orr 2003b). For the RMF model, simulations displayed in Figure 5 show that for very small *c*, the walk length is still on average equal to ℓ = *e* − 1. For larger *c*, the mean walk length first decreases slightly (see curve corresponding to *c* = 0.3) and then increases rapidly, until ℓ = *d*, the limit expected for a smooth landscape. For all values of *c*, the average walk length remains independent of the starting rank.

Numerical results obtained from simulations of the full fitness-dependent AW with transition probability in Equation 21 are shown in Figure 6. These simulations were carried out for *L* = 2000 loci and the random component of the fitness was drawn from a normal distribution with unit variance. While in Figure 6A the starting rank was kept constant, the initial Hamming distance to *σ**, *d*, was varied. For *c* = 0.01 (Figure 6A, inset), the behavior seems to be independent of *d*, while a *d* dependence starts to emerge for *c* = 0.3. For *c* = 1 the relation between initial distance *d* and the walk length is roughly linear with a slope smaller than unity.

For constant and various choices of the starting rank, one can observe the following behavior (Figure 6B). While for *c* = 1 the walk length seems to be independent of the starting rank, the data for *c* = 0.01 can be fitted by the relation (24)which was first obtained by Orr (2002) for the HoC landscape with Gumbel-distributed fitness values. For *c* = 0.3 the dependence on the starting rank is similar but the constant in Equation 24 is significantly larger.

Equation 24 is a special case of the general relation (25)derived by Neidhart and Krug (2011) and Jain (2011), which holds for HoC landscapes with fitness distributions characterized by an EVT index *κ* ≤ 1; for *κ* > 1 the mean walk length is asymptotically independent of starting rank. To see how this behavior is modified in the RMF model, we carried out simulations with the random fitness component drawn from the GPD, keeping *c* = 0.5 fixed and varying the EVT index *κ*; see Figure 7. The *d* dependence of the data seems to be well fitted by functions linear in *d*, with a slope that increases rapidly with decreasing *κ* when *κ* < 0; as was pointed out previously, the effect of the mean fitness gradient is particularly pronounced for distributions in the Weibull class. Somewhat surprisingly, the dependence on the starting rank at fixed *d* can be approximately described by the functional form in Equation 25 obtained for the HoC model but with a constant depending on *c* and *κ*; see Figure 7B.

Summarizing, inspection of the numerical data suggests the following dependencies of the mean adaptive walk length: a linear dependence on *d* with a slope that increases with increasing *c* and decreasing *κ* and a logarithmic dependence on the starting rank, similar to that known from the HoC model, with a constant offset depending on *c*, *κ*, and *d*. These findings are captured in the following conjectured expression for the adaptive walk length on RMF landscapes, (26)with so far unknown, nonlinear functions *α*, *β* with *α*(0, *κ*) = 0 and *β*(0, *κ*) > 0 (see Neidhart and Krug 2011 for a discussion of the *κ* dependence of the constant term in Equation 25).

### Crossing probability

While the adaptive walk length is a measure of the length of *typical* adaptive trajectories, it is also of interest to ask how likely it is for the population to traverse the entire landscape. To quantify this feature we introduce the *crossing probability*, which is the probability that an AW starting at the maximal distance *L* from the reference sequence *σ** reaches it and terminates there. For such an event to happen, three conditions must be fulfilled: the reference sequence *σ** must be a local maximum; there must exist at least one fitness-monotonic path connecting the antipodal sequence to *σ**; and finally, such a path must be chosen by the AW. The probability for the first condition was evaluated above and is given by (27)for Gumbel-distributed random fitness components, and this obviously constitutes an upper bound on the crossing probability. The probability for the existence of fitness-monotonic pathways in the RMF model has been investigated previously for the case when the paths end at the global fitness maximum (Franke *et al.* 2011), and it has been shown that such paths exist with unit probability for large *L* and any *c* > 0 (Hegarty and Martinsson 2014). It is not clear whether this result applies in the present setting, however, because the probability that the global maximum coincides with the reference sequence vanishes for large *L* (see Equation A29). Figure 8 shows numerical data for the crossing probability in comparison with the upper bound in Equation 27. Both quantities follow a sigmoidal behavior with a fairly sharp transition from zero to unity around a characteristic value of *c*, which increases roughly logarithmically with *L*, as would be expected on the basis of Equation 27.

## The Number of Exceedances

In this section we consider a feature of the fitness landscape that provides a distribution-free statistical test of the assumption of the MLM that fitness values of different genotypes are i.i.d. random variables (Miller *et al.* 2011). To define the quantity of interest, suppose that an adaptive step is taken from a starting genotype *σ* with rank *i* in its neighborhood *ν*(*σ*) to a genotype *σ*′ with rank *j* < *i* in *ν*(*σ*). The number of exceedances (NoE) is then equal to the number of neighboring genotypes in *ν*(*σ*′) that are fitter than *σ*′, that is, the rank of *σ*′ in its own neighborhood minus one. Since the only sequences present in both neighborhoods *ν*(*σ*) and *ν*(*σ*′) are *σ* and *σ*′, the NoE can principally vary between 0 and *L* − 1. Under the HoC assumption that fitness values are i.i.d. random variables, a classic result due to Gumbel and Von Schelling (1950) states that the distribution of the NoE is independent of the fitness distribution, and for large *L* the mean NoE is equal to the rank of *σ*′ in the initial neighborhood *ν*(*σ*) (Rokyta *et al.* 2006). In other words, if the first adaptive step goes to a genotype of rank *j*, the expected number of beneficial mutations available for the second step is *j*. In contrast, in a purely additive landscape such as the RMF model for very large *θ*, the NoE of a genotype is equal to its distance *d* from the global fitness maximum and independent of its rank in the initial neighborhood.

In their evolution experiments with the microvirid bacteriophage ID11, Miller *et al.* (2011) identified 9 beneficial second-step mutations on the background of a mutation, named g2534t, that had been found to have the largest effect among 16 beneficial first-step mutants. To apply the result of Gumbel and Von Schelling (1950), the rank of this mutation among *all* possible first-step mutations (not only the 16 observed in the experiment) has to be estimated. Assuming conservatively that the rank of g2534t among all beneficial first-step mutations is at most 3, at most 3 beneficial second-step mutations would then have been expected if fitness values were identically and independently distributed. Thus, the observation of 9 beneficial second-step mutations allowed Miller *et al.* (2011) to reject the HoC hypothesis with high confidence (*P* < 0.02).

Here we ask whether the RMF model is capable of yielding predictions for the NoE that are compatible with the experimental findings of Miller *et al.* (2011). Note that, unlike HoC landscapes, RMF landscapes are not isotropic and the NoE will depend on the position of genotypes *σ* and *σ*′ on the landscape, *i.e.*, their distance to the reference state, and on whether the adaptive step was taken in the uphill or the downhill direction. In contrast to the universal result of Gumbel and Von Schelling (1950), there will also be a dependence on the probability distribution of fitness values.

In *Appendix C* we present an approximate analytic calculation of the expected NoE for the RMF model, assuming an exponential distribution for the random fitness component. While the complete expressions for the expected NoE displayed in Equations C8–C13 are fairly complex, for small *c* they reduce to the simple form (28)Here ^{up} () is the expected number of exceedances after an adaptive step in the uphill (downhill) direction, and *r* is the rank of the mutated genotype *σ*′ in the initial neighborhood. For the HoC landscape (*c* = 0) Equations 28 yield = = *r* + 1, which differs slightly from the exact result = *r* as a consequence of the approximations involved in the derivation. Figure 9 compares the full expressions derived in *Appendix C* to numerical simulations, showing good agreement. Interestingly, for the case of an uphill step, the expected number of exceedances is maximal for landscapes of intermediate ruggedness.

Equation 28 shows that a considerable enhancement of the NoE is possible for moderate *c*, provided the adaptive step is in the uphill direction and the random fitness components are exponentially distributed. However, as we do not know whether the beneficial mutations that were observed in the experiments of Miller *et al.* (2011) correspond to uphill or downhill steps in a presumed RMF landscape, we need to average the predictions for the two cases, weighted with the probabilities for each of the two types of transitions to have happened. Furthermore, for an unbiased comparison it is appropriate to consider more general distributions for the random component than the exponential distribution discussed above. Here we choose the GPD, which allows us to cover distributions corresponding to all extreme value classes by varying a single parameter. Extending the approximate derivation in *Appendix C* to this general case is complicated and not very enlightening. Therefore, in the following we use numerical simulations to estimate the NoE.

To find parameter combinations that match the experimentally observed NoE, we constructed RMF landscapes, fixing the GPD index in the interval *κ* ∈ [−1, 0.5] and sampling *κ* with a resolution of Δ*κ* = 0.01. For each value of *κ*, the NoE was calculated for a range of *c* values. Following Miller *et al.* (2011) we assume that the first-step mutation has rank *r* = 1 in the initial neighborhood and determine for each choice of *κ* the smallest value of *c* for which (1) ≥ 9 (Figure 10A). For comparison, results assuming initial rank *r* = 3 are shown in Figure 10B. Note that, similar to the case of the exponential distribution displayed in Figure 9, for some *κ* there exists a second, larger *c* yielding the same value for . The results displayed in Figure 10 show that the strength of the fitness gradient *c* required to reproduce the experimentally observed NoE depends sensitively on the EVT index and becomes very small for negative *κ* deep inside the Weibull domain.

This is in accordance with the generally stronger effect of the additive fitness contribution for negative *κ* discussed above in the context of the statistics of fitness maxima and reflects the fact that for strongly negative *κ* the random variables drawn from the GPD probability density tend to crowd near the upper boundary of its support. For the EVT index *κ* = −0.29 estimated by Miller *et al.* (2011) from a maximum-likelihood analysis of the fitness values of 16 first-step beneficial mutations, the *c* value required to produce at least nine exceedances varies between 0.86 and 1.04 for initial rank 1 and between 0.44 and 0.76 for initial rank 3, depending on the assumed distance *d* to the reference sequence. Normalizing *c* by the standard deviation of the random fitness component, which equals 0.617 for this value of *κ* (compare to Equation C14), this translates into the intervals 1.39 ≤ *θ* ≤ 1.69 and 0.71 ≤ *θ* ≤ 1.23 for the parameter *θ*, respectively, which in comparison to estimates for other empirical data sets is indicative of a relatively smooth landscape. To further narrow down the range of RMF model parameters that can provide a consistent description of the fitness landscape of the ID11 microvirid system would require including other experimental observables into the analysis, which is beyond the scope of this article.

## Discussion

Motivated by the increasing availability of empirical information about the structure of adaptive landscapes, we have presented a detailed analysis of a one-parameter family of tunably rugged fitness landscapes. The model is a variant of the rough Mount Fuji landscape, in which the locus-specific additive effects considered in the original version (Aita *et al.* 2000; Aita and Husimi 2000) are replaced by a single parameter, and ruggedness is governed by the ratio *θ* of the additive fitness effect to the standard deviation of the random fitness component.

### Landscape topography

The mathematical simplicity of the model allowed us to derive several explicit results for the number and positions of local and global maxima; such results are much harder to come by for other generic fitness landscape models, notably Kauffman’s NK model. In particular, we have arrived at a complete classification of the expected number of fitness maxima in the limit of large genotype dimensionality *L*, which highlights the importance of the tail behavior of the distribution of the random fitness component in determining the ruggedness of the landscape. For distributions with tails heavier than exponential, the additive fitness component becomes asymptotically irrelevant, in the sense that the number of maxima is equal to the value expected for an uncorrelated HoC landscape. In contrast, for distributions with bounded support the number of maxima is reduced compared to the HoC expectation by a factor that varies exponentially in *L*. The interplay between additive effects and tail behavior has been a recurrent theme of this article, which manifests itself in various quantities of interest. Our results show that both the parameter *θ* characterizing the additive effects and the EVT index *κ* must be specified for a comprehensive description of the landscape topography in the RMF model.

An exception to this rule is provided by the fitness correlation function, which, similar to the NK model, is independent of the type of randomness. In contrast to the NK model, however, the correlations become negative at large distances, which reflects the inherent anisotropy of the RMF landscape and the long-range effect of the fitness gradient.

Another important measure of fitness landscape ruggedness not discussed so far in this article is the existence and abundance of selectively accessible mutational pathways, defined as paths composed of single mutational steps along which fitness increases monotonically (Weinreich *et al.* 2006; Franke *et al.* 2011). Using an approach similar to that of the present work, an explicit expression for the expected number of accessible pathways in the RMF model with Gumbel-distributed randomness can be derived (Franke *et al.* 2010, 2011). Subsequently Hegarty and Martinsson (2014) presented a rigorous proof that accessible pathways exist with unit probability for large *L* in the RMF model for any *c* > 0, independent of the distribution of the random fitness component. This is in stark contrast to the behavior in the HoC model (*c* = 0), where the probability for existence of accessible paths tends to zero for large *L*. Analyses in which the genotypes are placed on a regular tree show that this strong dichotomy between the HoC and RMF models is specific to the hypercube topology of sequence space (Nowak and Krug 2013; Roberts and Zhao 2013).

### Dynamics of adaptation

Apart from being amenable to rigorous analysis, the RMF model is useful for exploring various aspects of evolutionary dynamics in rugged fitness landscapes through simulations. Recent applications in this context include studies of evolutionary predictability (Lobkovsky *et al.* 2011; Szendro *et al.* 2013a), epistatic interactions between mutations occurring along an adaptive walk (Greene and Crona 2014), and the effect of recombination on rugged fitness landscapes (Nowak *et al.* 2014). Here we have focused on adaptation in the SSWM regime and considered both single adaptive steps and adaptive walks to local fitness maxima. Interestingly, while the statistics of single adaptive steps largely conform to the classic results obtained for the MLM (Orr 2002; Joyce *et al.* 2008), adaptive walks in the RMF are much longer than in the MLM already for small values of *c*. Specifically, the heuristic expression in Equation 26 that summarizes the simulation results suggests a linear dependence of the walk length on the initial distance to the reference sequence.

The qualitatively different effects that the fitness correlations in the RMF have on single *vs.* multiple adaptive steps highlight the fact that a step in an adaptive walk involves two distinct random processes (Rokyta *et al.* 2006; Neidhart and Krug 2011). The first process is the selection of a fitter neighbor according to the transition probability in Equation 21, and the second process is the change of the mutational neighborhood after the fixation of the mutated genotype. In the MLM the effect of the second process is relatively weak, and as a consequence adaptive walks are well described by an approximation that ignores the neighborhood change (Orr 2002; Neidhart and Krug 2011).

To understand the role of the neighborhood change in the RMF model we refer to the analysis in *Appendix C*, where it is shown (for exponentially distributed randomness) that the effect of the fitness gradient can be approximately subsumed into an overall shift of all the fitness values constituting a neighborhood by the same amount. This implies that the transition probability in Equation 21, which depends only on fitness differences within the neighborhood, is approximately independent of *c*. However, since the shift is a function of the distance *d* to the reference sequence that changes in the adaptive step, the rank of the mutant genotype in the new neighborhood is strongly dependent on whether the step occurred in the uphill (*d* → *d* − 1) or downhill (*d* → *d* + 1) direction. Further investigations are required to elucidate how this effect gives rise to the observed dependence of the walk length on the landscape parameters. A promising approach is to consider walks starting from the antipode of the reference sequence, which can be assumed to move almost exclusively in the uphill direction provided the walk length remains short compared to *L* (Park *et al.* 2014).

For large populations and/or high mutation rates the SSWM approximation underlying the simple adaptive walk picture breaks down, and additional processes such as the competition between beneficial clones and the crossing of fitness valleys have to be taken into account. The resulting complex population dynamics are governed by the tension between a tendency toward greater determinism induced by clonal competition (Jain and Krug 2007; Jain *et al.* 2011) and the increasing role of nonmonotonic pathways that become accessible through valley crossing (Szendro *et al.* 2013a; De Visser and Krug 2014). In contrast to the SSWM regime, the trapping of large populations at local fitness maxima is only a transient occurrence, and the rate of escape from such peaks plays an important role in the comparison between recombining and nonrecombining populations on rugged fitness landscapes (Nowak *et al.* 2014).

### Application to experiments

The usefulness of the RMF model for the quantitative description of empirical fitness landscapes has been documented in several recent studies. Franke *et al.* (2011) applied the model to an eight-locus fitness landscape for the fungus *A. niger* and extracted an estimate of *c* from a subgraph analysis of pathway accessibility. In a study of amplitude spectra of fitness landscapes Neidhart *et al.* (2013) showed that the correlation function of a six-locus fitness landscape obtained by Hall *et al.* (2010) for the yeast *Saccharomyces cerevisiae* is well described by the RMF model. Finally, in a metaanalysis of 10 empirical fitness landscapes Szendro *et al.* (2013b) used the RMF model to interpolate the behavior of various ruggedness measures between the limits of a completely random (HoC) and an additive landscape and found good agreement with the trends in the empirical data.

In the present article we have complemented these analyses by considering the effect that the fitness gradient in the RMF model has on the number of secondary beneficial mutations that are available after an adaptive step. We have identified model parameters for which the RMF prediction matches the large number of fitness exceedances observed in the experiment of Miller *et al.* (2011). A small fitness gradient suffices to explain the experiments when the distribution of the random fitness component is assumed to belong to the Weibull class of EVT, as is suggested by the analysis of the distribution of first-step mutational effects. We believe that more work along these lines, focusing on the changes in the statistics of mutational neighborhoods along an adaptive trajectory, will provide important insights into the role of epistatic interactions during adaptation and the viability of schematic models like the one considered here.

## Acknowledgments

We thank Su-Chan Park for discussions and two anonymous reviewers for helpful comments. J.K. acknowledges the kind hospitality of the Simons Institute for the Theory of Computing, Berkeley, CA, and the Galileo Galilei Institute for Theoretical Physics, Florence, Italy, during the completion of this article. This work has been supported by the Deutsche Forschungsgemeinschaft within SFB 680, SFB TR12, and SPP 1590 and the Bonn–Cologne Graduate School of Physics and Astronomy.

## Appendix A

### Properties of Fitness Maxima

#### Density of local maxima

The probability that a genotype at distance *d* from the reference sequence is a local fitness maximum is given by the integral (A1)This is simply the probability that the genotype’s fitness *x* exceeds that of its uphill and downhill neighbors, averaged with respect to the probability density *p*(*x*). Unless specified otherwise, here and in the following the domain of integration is equal to the support of the probability distribution. For the Gumbel distribution the integral (A1) can be evaluated exactly using the shift property (6), yielding the result in Equation 7 of the main text.

Similarly, the probabilities to find the neighboring genotype of largest fitness in the uphill or downhill direction, respectively, are given by (A2) (A3)which can be evaluated for the Gumbel distribution to yield the expressions in Equation 10.

#### Total number of maxima

Inserting Equation A1 into the sum in Equation 11 and exchanging the order of integration and summation, one arrives at the compact expression (A4)which is evaluated in the following sections for various special cases. Equation A4 makes it evident that ℳ is an even function of *c*: changing *c* to −*c* produces a fitness landscape with the antipodal reference sequence that is statistically equivalent to the original landscape.

##### Exponential distribution:

For the exponential distribution defined by (A5)the expression in Equation A4 becomes (A6)Substituting *z* = *e*^{−}* ^{x}*, this yields (A7)It is straightforward to check that ℳ → 1 for

*c*→ ∞. Moreover, the asymptotics for large

*L*at fixed

*c*are identical to those derived in Equation 12 for the Gumbel distribution, . A graphical illustration of Equation A7 is shown in Figure A1A.

##### Gumbel distribution:

Inserting the exact expression in Equation 7 into Equation 11, one obtains (A8)We first show that the sum in Equation A8 can be expressed in terms of a hypergeometric function defined by(Graham *et al.* 1994) where the Pochhammer symbol is defined byThe defining feature of the hypergeometric function is that the terms *t _{n}* satisfy

*t*

_{0}= 1 and (A9)To bring Equation A8 into this form we write (A10)ensuring that

*t*

_{0}= 1, and compute the fractions

*t*

_{d}_{+1}/

*t*according to (A11)By comparison with Equation A9 the arguments

_{d}*a*,

*b*,

*c*, and

*z*can be identified and it follows that (A12)The asymptotic behavior for large

*L*is most conveniently extracted from the general expression in Equation A4, which in this case can be brought into the form (A13)Recognizing that the dominant contribution to the integral comes from the region near

*P*= 1 and expanding the integrand around this point, one readily finds that for large

*L*, in agreement with Equation 12. Along similar lines it can be shown that the same asymptotics obtain for any distribution with a strictly exponential tail.

##### Fréchet class:

Here we show that for distributions in the Fréchet class the expected number of maxima is asymptotically equal to that in the HoC model, for large *L* and any *c* > 0. To simplify the notation we use a standard Pareto distribution defined by (A14)where the exponent *α* is related to the EVT index through Similar to the case of the exponential distribution, the domain of integration in Equation A4 has to be subdivided into the interval [1, 1 + *c*), where *P*(*x* − *c*) = 0, and the remainder [1 + *c*, ∞). The dominant contribution comes from the second part, which is given by (A15)for large *L*. Substituting *y* = *Lx*^{−}* ^{α}*, this yields (A16)as claimed; see Figure A1B for illustration.

##### Gumbel class:

We are now in the position to generalize the results obtained for the exponential and Gumbel distributions to the entire Gumbel class of EVT. For calculational convenience we choose the one-parameter family of Weibull distributions defined by (A17)to represent Gumbel-class distributions with tails that are heavier (for *β* < 1) or lighter (for *β* > 1) than exponential. As in the cases considered above, the dominant contribution to the number of maxima comes from the part of the integral in Equation A4 that extends from *c* to infinity, which now reads (A18)The crucial difference between the cases *β* < 1 and *β* > 1 lies in the behavior of the ratio of the two exponential terms inside the inner square brackets. For *β* < 1 (A19)which implies that the shift by ±*c* becomes irrelevant asymptotically and independent of *c*, as in the Fréchet class. On the other hand, for *β* > 1 one finds that for large *x* and Equation A18 simplifies to (A20)For large *L* the integrand is effectively zero below a cutoff scale *x** determined by or Thus (A21)to leading order in *L*.

##### Weibull class:

We next consider distributions with bounded support, which we take to be the unit interval [0, 1]. In the evaluation of the general expression in Equation A4 it has to be taken into account that *P*(*x* − *c*) = 0 for *x* < *c* and *P*(*x* + *c*) = 1 for *x* > 1 − *c*, which implies in particular that the cases (where *c* < 1 − *c*) and (where *c* > 1 − *c*) have to distinguished. We begin with the uniform distribution and assume which yields (A22)Analogously for we obtain (A23)Comparing Equations A22 and A23, the two expressions are seen to coincide at Moreover, Equation A22 reduces to for *c* = 0 and Equation A23 confirms that ℳ = 1 for *c* = 1, as expected.

As representatives of Weibull-class distributions with a general EVT index *κ* < 0 we consider the Kumaraswamy distributions defined by (A24)with We focus on the leading-order behavior of the number of maxima for large *L*, which is given by the part of the integral in Equation A4 that contains the upper boundary of the support at *x* = 1 (compare to Equations A22 and A23). We assume and obtain (A25)which is dominated by the region near *y* = 0 for large *L*. Expanding (*y* + *c*)* ^{ν}* for small

*y*, it follows that (A26)for large

*L*, where Γ denotes the Gamma function. The result for is asymptotically the same. Thus to leading order the number of maxima is in the Weibull class.

##### Expansion for small *c*:

A unified view of the effect of the fitness gradient on the expected number of maxima can be obtained by expanding the general expression in Equation A4 for small values of *c*. To leading order in *c* the integrand is *P*(*x* − *c*) + *P*(*x* + *c*) = 2*P*(*x*) + *c*^{2}*P*″(*x*) = 2*P*(*x*) + *c*^{2}*p*′(*x*), where primes indicate derivatives with respect to *x*. Integrating and keeping terms of order *c*^{2}, we thus obtain (A27)Evaluating the integral on the right-hand side for the GPD in Equation 4, one finds that it is of the order of *L*^{−(3+2}^{κ}^{)}, and hence the entire correction term is of order 2^{L}L^{−(1+2}^{κ}^{)}*c*^{2}. Thus for *κ* > 0 (Fréchet class) the correction becomes negligible compared to the HoC term for large *L*. In contrast, for *κ* < 0 (Weibull class) the correction term eventually dominates the HoC term, showing that the leading-order behavior of ℳ is modified.

#### Location of the global maximum

To compute the probability that a sequence at distance *d* is the *global* fitness maximum one has to demand that its fitness exceeds the fitnesses of all other sequences, which leads to the expression (A28)Inserting the Gumbel distribution in Equation 5 and using its shift property in Equation 6, this can be evaluated according to (A29)For general fitness distributions one has to resort to an expansion in *c*. Starting from Equation A28 and collecting terms linear in *c*, one obtains (A30)where (A31)Expressions for *I _{L}* for representatives of the three extreme value classes have been derived by Franke

*et al.*(2010). For large

*L*they behave as (A32)(Wergen

*et al.*2011). The weighted average with respect to

*d*then yields the approximate expression (A33)for the mean distance of the global maximum to the reference sequence. Using the asymptotic behavior of

*I*given in Equation A32, it follows that the shift in the position of the maximum from the HoC value

_{L}*L*/2 is of order

*cL*2

^{−}

*.*

^{κL}## Appendix B

### Calculation of the Fitness Correlation Function

For notational convenience we substract the mean value (*η*) of the random fitness component from Equation 2 and define (B1)such that (*ξ*) = 0. It is then easy to see that (B2)where *v* = Var(*η*) is the variance of *η* (or *ξ*), which is assumed in the following to exist, and *d* = *D*(*σ*, *σ**), *d*′ = *D*(*σ*′, *σ**). While the sequence space averages of *d* and *d*′ are obviously equal to *L*/2, the evaluation of 〈*dd*′〉* _{r}* requires a double summation, first over sequences

*σ*at distance

*d*from

*σ** and then over sequences

*σ*′ at distance

*r*from

*σ*. The latter sequences are grouped according to a number

*k*that counts how many of the

*r*point mutations distinguishing

*σ*′ from

*σ*fall on alleles that are different in

*σ*and

*σ**. Obviously, each such mutation decreases the distance

*d*′ by 1, while each of the

*r*−

*k*mutations acting on previously unaltered alleles increases

*d*′ by 1, such that

*d*′ =

*d*−

*k*+ (

*r*−

*k*) =

*d*+

*r*− 2

*k*. The number of sequences

*σ*′ with a given value of

*k*is equal to . Thus the sum to be evaluated is (B3)where the combinatorial identitieshave been used (Graham

*et al.*1994). Putting everything together finally yields Equation 19.

## Appendix C

### Calculation of the Number of Exceedances

We have seen above in *Model* that the full neighborhood *ν*(*σ*) of a sequence *σ* is divided into the uphill neighborhood with the corresponding distribution function *P*^{↑}(*x*) = *P*(*x* + *c*(*d* − 1)) of fitness values, *σ* itself with distribution function *P*^{•}(*x*) = *P*(*x* + *cd*), and the downhill neighborhood with distribution function *P*^{↓}(*x*) = *P*(*x* + *c*(*d* + 1)). The full distribution function of fitness values is then given by (C1)and the expectation of the *k*th largest fitness value is obtained as (C2)(David and Nagaraja 2003). In general, the evaluation of this expression is complicated, because the different components of Π do not have the same support. Here we show how this problem can be circumvented in an approximate way for the special case of an exponential distribution fitness distribution *P*(*x*) = 1 − *e*^{−}* ^{x}*. Naively inserting this into Equation C1, we obtain (C3)with (C4)Our approximation consists in ignoring the fact that Equation C3 holds only on the intersection of the supports of

*P*

^{↑},

*P*

^{•}, and

*P*

^{↓}and instead defining the

*common support*of Π(

*x*) by [log(

*ξ*(

*c*,

*d*,

*L*)), ∞), such that Π(log(

*ξ*)) = 0. Thus the full distribution of fitness values defined in Equation C1 is replaced by a simple exponential that is shifted in a

*d*-dependent way.

To proceed we recall that the expected value of the *k*th largest of *n* exponentially distributed random variables is given by (C5)(David and Nagaraja 2003), where are the harmonic numbers, and we use the convention that *H*_{0} = 0. In the second part of Equation C5 the logarithmic approximation *H _{n}* ≈ log(

*n*) +

*γ*valid for large arguments has been applied, with

*γ*≈ 0.5772156649 … . It follows that the mean of the

*k*th largest fitness value in the neighborhood, defined in Equation C2, is approximately given by

#### After a step up

To obtain the mean NoE after a transition from a sequence *σ* at distance *d* to *σ*′ at distance *d* − 1, where *σ*′ has rank *r* in the old neighborhood, we need to count how many times *F*(*σ*′) is exceeded in the new neighborhood. As an estimation, *F*(*σ*′) ≈ *μ _{r}*(

*L*,

*c*,

*d*) is compared to the mean rank-ordered fitness values from the uphill and the downhill parts of the new neighborhood, which are sets of

*d*− 1 and

*L*−

*d*+ 1 independent exponential variables, respectively. In the uphill neighborhood at distance

*d*− 2 from the reference sequence the exponential random variables are shifted by

*c*(

*d*− 2), and therefore the number of exceedances

*k*

_{>}derived from this part of the neighborhood is obtained by solving the relation (C7)for

*k*

_{>}. Using the approximation in Equation C6, this yields the expression (C8)Similarly, the contribution

*k*

_{<}to the exceedances from the downhill neighborhood is obtained from the relation which yields (C9)To complete the calculation we have to take into account the fact that, by construction, and which is not always satisfied by the approximate expressions in Equations C8 and C9. Incorporating these constraints, we arrive at our final result (C10)A simpler and more transparent expression can be obtained by assuming that

*d*≫ 1 and

*e*is not too large. Under the first assumption the combination of Equations C8 and C9 reduces to while the second assumption ensures that the min-constraints in Equation C10 can be ignored, such that see Equation 28 in the main text.

^{c}#### After a step down

The calculation of the NoE after a step is taken in the downhill direction is analogous to the previous one. In this case the contribution from the upper part of the new neighborhood is obtained from the relation which yields (C11)Correspondingly, the contribution from the downhill part is obtained from evaluating with the result (C12)The final estimate for the number of exceedances reads (C13)and the simplified expression in Equation 28 arises from the same approximations employed previously for .

#### Analysis of the experiment of Miller *et al.* (2011)

In Figure 10 we report the RMF model parameter combinations (*κ*, *c*) required to explain the nine beneficial second-step mutations observed on the background of a highly fit first-step mutation in the experiment of Miller *et al.* (2011). Figure C1 displays the same data with the additive selection coefficient *c* scaled by the standard deviation of the random fitness component, which in the case of the GPD is given by

## Footnotes

*Communicating editor: J. Hermisson*

- Received June 23, 2014.
- Accepted August 8, 2014.

- Copyright © 2014 by the Genetics Society of America