Speed of evolution in large asexual populations with diminishing returns

Speed of evolution in large asexual populations with diminishing returns

Maria R. Fumagalli Genomic Physics Group, CNRS (UMR 7238) “Microorganism Genomics”, 15 Rue de l’École de Médecine, 75006 Paris, France Dipartimento di Fisica, Università degli Studi di Milano, Via G. Celoria 16, Milano, Italy Dipartimento di Fisica, Università degli Studi di Torino, Via P. Giuria 1, Torino, Italy    Matteo Osella Genomic Physics Group, CNRS (UMR 7238) “Microorganism Genomics”, 15 Rue de l’École de Médecine, 75006 Paris, France Université Pierre et Marie Curie, 4 Place Jussieu, 75005 Paris, France.    Philippe Thomen Université Pierre et Marie Curie, 4 Place Jussieu, 75005 Paris, France. Laboratoire Pierre Aigrain, Ecole Normale Supérieure, CNRS (UMR 8551), Université P. et M. Curie, Université D. Diderot, 24 rue Lhomond, 75005 Paris, France    Francois Heslot Université Pierre et Marie Curie, 4 Place Jussieu, 75005 Paris, France. Laboratoire Pierre Aigrain, Ecole Normale Supérieure, CNRS (UMR 8551), Université P. et M. Curie, Université D. Diderot, 24 rue Lhomond, 75005 Paris, France    Marco Cosentino Lagomarsino marco.cosentino-lagomarsino@upmc.fr Genomic Physics Group, CNRS (UMR 7238) “Microorganism Genomics”, 15 Rue de l’École de Médecine, 75006 Paris, France Dipartimento di Fisica, Università degli Studi di Torino, Via P. Giuria 1, Torino, Italy Université Pierre et Marie Curie, 4 Place Jussieu, 75005 Paris, France.
July 5, 2019

The adaptive evolution of large asexual populations is generally characterized by competition between clones carrying different beneficial mutations. This interference phenomenon slows down the adaptation speed and makes the theoretical description of the dynamics more complex with respect to the successional occurrence and fixation of beneficial mutations typical of small populations. A simplified modeling framework considering multiple beneficial mutations with equal and constant fitness advantage captures some of the essential features of the actual complex dynamics, and some key predictions from this model are verified in laboratory evolution experiments. However, in these experiments the relative advantage of a beneficial mutation is generally dependent on the genetic background. In particular, the general pattern is that, as mutations in different loci accumulate, the relative advantage of new mutations decreases, trend often referred to as “diminishing return” epistasis. In this paper, we propose a phenomenological model that generalizes the fixed-advantage framework to include in a simple way this feature. To evaluate the quantitative consequences of diminishing returns on the evolutionary dynamics, we approach the model analytically as well as with direct simulation. Finally, we show how the model parameters can be matched with data from evolutionary experiments in order to infer the mean effect of epistasis and derive order-of-magnitude estimates of the rate of beneficial mutations. Applying this procedure to two experimental data sets gives values of the beneficial mutation rate within the range of previous measurements.

I Introduction

Although the links between the statistical theory of evolution and statistical physics can hardly be considered novel Park et al. (2010), in the recent years they are attracting the attention of researchers in a renewed way. What changed is the increased availability of genomic and experimental data, in amounts and levels of precision that (despite there is still a considerable room for improvement), were difficult to envisage in the past century. These innovations are consolidating the biology and genomics of evolution into a more stable field of application of modeling ideas from statistical physics.

In particular, contemporary technology allows the realization of controlled laboratory evolution experiments, which can guide theoretical investigations and in principle makes the validation and falsification of phenomenological theories feasible Hindré et al. (2012). These experiments are often performed with large populations of microorganisms, and allow to explore adaptation under well-defined sources of natural selection. Moreover, phenotypic characterization and high-throughput sequencing give a quantitative insight into the genomic adaptation of these microbial populations Barrick et al. (2009); Tenaillon et al. (2012), with notable consequences in a wide range of bio-technological and ecological contexts. For these large asexual (or rarely mating) populations, a high number of beneficial mutations emerge in different clones, and cannot be mixed because of slow or absent recombination. These beneficial mutations appearing in parallel coexist and compete to drive adaptation. This phenomenon of concurrent beneficial mutations (sometimes generically termed “clonal interference”), is related to the Fisher-Muller hypothesis (or Hill-Robertson effect) for the advantage of recombination Felsenstein (1974). In general, beneficial mutations also arise with a distribution of fitness advantage, which is generally believed to be exponential Orr (2003).

Recent models have generally dealt with the competition between mutations of different strengths and the competition between mutations that arise on different fitness backgrounds separately. The first effect, the role of a distribution of fitness effects, is analyzed by so-called “clonal interference” models 111Note that the term has a stricter sense in this case. In the following we will reserve the term clonal interference to this stricter meaning of competition between mutations of different strength, and talk of interference between beneficial mutations in the generic case. Gerrish and Lenski (1998); Wilke (2004); Park et al. (2010). In these models any individual is either the wild type or a mutant derived directly from the wild type. Thus, multiple mutations arising from the extant clones are neglected. Conversely, models that explicitly deal with multiple mutations typically assume that all mutations have the same effect Tsimring et al. (1996); Desai and Fisher (2007); Brunet et al. (2008); Park et al. (2010). The latter kind of model has the advantage of being simpler to treat and accessible analytically, and is characterized by a Gaussian-like traveling wave for the histogram of log-fitness throughout the population. In absence of epistasis, this wave moves towards higher log-fitness with a constant speed and shape. New mutations are fixed in the population if they occur in the high-fitness edge (or “nose”) of the distribution. Consequently, a quantitative law relates the width of the log-fitness histogram and the adaptation speed. Recent work incorporating both effects Good et al. (2012) has shown that if the distribution of fitness advantages given by beneficial mutations is sufficiently peaked around a single characteristic value, the evolutionary dynamics can be described by an effective theory with a single typical selection coefficient and a rescaled beneficial mutation rate.

While the models described so far have been matched successfully with the diversity and adaptation speed of short-time laboratory evolution experiments Desai et al. (2007), there appears to be one important discrepancy between the models and the behavior of bacteria evolved in the laboratory for longer times (roughly, generations). Furthermore, two recent experimental studies Khan et al. (2011); Chou et al. (2011) have shown a general pattern in the advantage of combined beneficial mutations occurring in different genes. This combined advantage is lower than the sum of that of individual mutations. In other words, when mutations of loci in different genes accumulate, the effective advantage of each of them is lower. This was shown by combinatorial genetics techniques, by constructing all the possible configurations of a small set of mutations, and evaluating their advantage through competition experiments. This trend, referred to as “diminishing returns” epistasis, had been previously suggested theoretically on the basis of the general pattern of adaptation observed in long-term microbial experiments Kryazhimskiy et al. (2009), using a modeling framework that neglected concurrent or multiple mutations. Another study predicts the same principle on the basis of a simple fitness landscape model combined with the distribution of single mutation effects measured experimentally Martin et al. (2007). The actual pattern in the fitness associated to the same mutation in different backgrounds observed by the two studies is complex, as, on top of the diminishing return effect, the advantage appears to depend on the mutation identity. Based on their data, Chou and coworkers Chou et al. (2011) defined an expression for the fitness made of two additive components. They consider the growth rate as the sum of two “constituent phenotypes”: metabolic rate and protein expression burden. When mutations are combined, the metabolic and expression contributions to the fitness combine multiplicatively in an independent way. They show that this model generates an accurate prediction for the fitness of all their multi-allele strains. Even more recent systematic experiments Tenaillon et al. (2012) are unveiling a complex scenario where different mechanisms coexist for the interactions of mutations between and within functional “blocks”, which can span multiple genes along the genome. One interesting modeling approach that was developed recently Schiffels et al. (2011) incorporates genetically linked multiple mutations explicitly, and describes interference interactions between multiple beneficial and deleterious mutations. This description has enough degrees of freedom to possibly account for many of the phenomena observed experimentally. However, the full experimental complexity is difficult to incorporate in a treatable model, and experimental data on linked mutations and interference between them are not easy to be obtained. A simplified descriptions in the spirit of the multiple mutations model can be useful in order to model evolving populations using a minimal quantity of information on mutations and fitness advantage.

Here, we take a simplified approach to investigate the diversity and speed of adaptation in presence of diminishing returns. We define a framework that can account for multiple mutations, and incorporates the effect of diminishing return epistasis in a simplified way. Namely, the fitness of a mutation depends only on its order of appearance in a clone, and decreases with it. This generalizes the standard multiple mutation model, which we recover in the case in which the fitness decrease with mutation index is zero. We preserve the model assumption that evolution is driven by beneficial mutations which appear with constant rate (see the Discussion for an evaluation of these assumptions in light of the results). We study the mean-field behavior of this model using standard techniques Park et al. (2010), and derive a relation that is analogous to Fisher’s theorem. Moreover, we discuss how the model assumptions cease to be valid when the rate of fixed beneficial mutations drops to the level where deleterious mutations should be accounted for. We also explore the behavior at finite population sizes, and obtain simple analytical quantitative estimates of the speed of adaptation that generalize that of the standard multiple mutation model. Finally, we show how the model parameters can be matched with data from two different long-time laboratory evolution experiments and how, assuming the model, one can derive order-of-magnitude estimates of the beneficial mutation rate and the mean effect of the diminishing returns.

Ii Basic features of the model

ii.1 Model definition

Figure 1: (Color online) Basic features of the model. (A) Because of competition between beneficial mutations, the population is divided into sub-populations with different frequencies (left panel), defined by the number of mutations . is the difference, in number of mutations, between the maximum number of mutations found in a clone and the mean. This induces a distribution for the log-fitness (right panel). Both distributions travel in time, driven by established beneficial mutations, with instantaneous velocities and . (B) The dynamics is driven by beneficial mutations and the rate of establishment of new classes. The dashed rectangles represent the class frequency histogram at a subsequent time. (C) Sketch of the algorithm used in the simulations (from ref. Park et al. (2010), see text). The instantaneous frequencies of mutation classes define, through Eq. 8 (which incorporates selection and mutation) the probability that an individual with mutations is found at the subsequent generation. The frequencies at the subsequent generations are then sampled from a multinomial distribution with parameters . This procedure allows to access large population sizes.

We model an asexual population of haploid individuals, or sequences, in which each individual of type produces a random number of offspring with average equal to its fitness . Inheritance is introduced by assigning the fitness of their parents to the offspring. The evolutionary dynamics of the population is based on the well-known Wright-Fisher model Wright (1931); Fisher (1930). We used the following reproduction scheme Park and Krug (2007). Each individual at generation is chosen as the offspring of an individual in class present at generation with probability , where is the relative fitness of class . This prescription for the evolutionary dynamics assumes non-overlapping generations since at each generation there is a complete replacement of parents with progeny. The population size is kept constant, and there is no recombination move. The presence of fitness differences in the population leads to natural selection since classes of individuals with higher fitness will generate increasingly larger fractions of the population as the dynamics proceeds, while classes with low fitness will progressively disappear.

Each offspring has a constant probability per generation (the “beneficial mutation rate”) of acquiring a beneficial mutation. Beneficial mutations increase the parental fitness to a higher value following the relation , where the positive parameter is the “selection coefficient” whose typical values in laboratory evolution experiments with bacteria are in the range  Hegreness et al. (2006); Perfeito et al. (2007). We assume that each new beneficial mutations hits a new site on the genome. While the fitness of new mutations is a complex issue Eyre-Walker and Keightley (2007), in presence of abundant beneficial mutations, deleterious mutations (negative effect on fitness, ) do not typically contribute to the adaptation of large populations and are customarily neglected Park et al. (2010); Desai and Fisher (2007); Brunet et al. (2008).

The conditions for the emergence of the interference phenomenon between mutations for a large population can be understood with simple scaling arguments as a competition between processes occurring at different time scales Park et al. (2010). When a new mutant with advantage arises, the probability that its lineage grows sufficiently in size to overcome genetic drift (stochastic fluctuations in the reproductive process) and start to deterministically expand in the population (i.e. it is “established”) is , where c is a constant factor that depends in the specific model used Haldane (1927). For the algorithm used here, (see Appendix, Fig. A.2222 see also ref. Park et al. (2010) for a self-contained motivation of these formulas.. Furthermore, since the extinction probability is , one can also estimate the population size conditioned to certainty of survival of the lineage as . In absence of additional beneficial mutations, once a lineage is established (i.e. has survived genetic drift and has roughly size ) it will take over the population (go to “fixation”) logistically at initial rate . Consequently, the scale of the fixation time can be estimated by imposing that , giving a characteristic time  Desai and Fisher (2007). On the other hand, the time scale for appearance and establishment of a new beneficial mutation is . Therefore, when () a beneficial mutation can fix before any new mutation can establish, making the evolutionary dynamics driven by successive sweeps of new lineages arising on an essentially clonal population. This regime is called “selective sweeps” or “periodic selection”. Instead, a sufficiently large population with high beneficial mutation rate evolves in the opposite regime, in which multiple mutations can establish before fixation of any of them and interfere with each other. This is the regime considered here, which is believed to be relevant for laboratory evolution experiments with microorganisms.

As mentioned in the introduction, the standard approach of multiple mutation models is to consider constant advantage for all mutations. This entails that, in absence of epistasis, the fitness can be assigned to individuals (sequences) with mutations, which (roughly) propagate with rate  Brunet et al. (2008). is the average number of accumulated mutations in a single realization of the process, i.e. the average over the distribution of shown in Fig. 1. For a range of parameter values, this constant advantage framework can be seen as an effective theory for a model with a distribution of fitness for beneficial mutations Good et al. (2012).

Following the spirit of this modeling framework, we build a minimal phenomenological model including diminishing return epistasis in presence of competition between beneficial mutations. The model assumes that successive mutations do not lead to the same fitness gain, but the fitness gain is dependent on the number of mutations already occurred in an individual Kryazhimskiy et al. (2009). A simple way to implement this feature is to extend the standard model for multiple mutations and consider selection coefficients dependent on the number of mutations, i.e. , where is a decreasing function of . As in the standard multiple-mutation model Desai and Fisher (2007); Brunet et al. (2008), the population can be divided into classes of individuals with the same number of mutations that are in one-to-one correspondence to fitness advantage classes, as represented in Fig. 1A. An individual with beneficial mutations has fitness


Since the rules defining the dynamics contain the relative fitness , the model is unaffected by multiplication of all the by a common factor. Therefore, the fitness value can be arbitrarily assigned to the genotype with no mutation (“ancestral” or “wild-type”).

We considered typical values of the parameter between and  Perfeito et al. (2007); Park et al. (2010). For the parameter , we explored the range  Perfeito et al. (2007); Hegreness et al. (2006), and we considered population sizes between and  Hindré et al. (2012). These ranges, which can be considered plausible with respect to what is known empirically, impose a hierarchy of scales in the effective parameters. For example, typically and is large. As for , we are interested in exploring sufficiently large values to ensure the regime of concurrent mutations. However, since the advantage given by a beneficial mutation decreases during the time-evolution of the model, this hierarchy can in part change over time. When the advantage becomes too small, i.e. for large times, some of the assumptions of the model become unrealistic. Namely, for sufficiently large times, deleterious mutations cannot be neglected, and the establishment size can eventually become comparable to . The inequality , which ensures concurrent mutations, can also break down for very large times, when decreases. This is not a problem - as any experiment spans a finite time. However, in contrast to the standard multiple mutation model, which describes a stationary state where the regime is entirely defined by the parameters, in this model the time-range of validity of the assumptions has to be kept in mind. We will comment on these issues in the following sections, when reviewing the results.

In order to fully define the model, one has to choose a specific functional form for the function , describing the strength of the negative epistasis between mutations. In general, every function leading to a sum that is sub-linear in defines a model with diminishing return. A simple example is given by the choice of a fitness gain that depends on the number of the extant mutations as a power law. In this case, with , where the epistasis grows in strength with decreasing from the non-epistatic case . In the case , the fitness of an individual with beneficial mutations (see Appendix A for the derivation) has the form


Appendix A describes two additional examples of the diminishing returns function that we have considered, leading to logarithmic (i.e. the case ) or geometric increase of the log-fitness. The comparison between the three variants will be useful in the second part of this work, when the model parameters are matched with experimental data.

ii.2 Background on the multiple mutation model with no epistasis

For , the power law return model reduces to the particular case of absence of epistasis (i.e. , hence ). This case coincides with the standard multiple mutation model for the evolutionary dynamics of a large asexual population Rouzine et al. (2003); Desai and Fisher (2007), in which every mutation carries a constant fitness advantage . The phenomenology of the multiple mutation model with constant advantage has been explored in detail theoretically Desai and Fisher (2007); Brunet et al. (2008) and some of the resulting predictions have been experimentally tested Desai et al. (2007).

As mentioned in the Introduction, one of the main results of these studies is that a steady state exists, where the population is organized in fitness classes (and the corresponding mutation classes, as illustrated in Fig. 1A), forming an approximately Gaussian traveling wave encompassing a constant number of mutation-classes and propagating to higher average fitness with constant velocity . Technically, this wave has all the characteristics of a soliton Rouzine et al. (2003); Park et al. (2010). Importantly, the wave is driven by the mutations that are fixed at the edge of the distribution of the advantage. Fig. 1B depicts schematically how the fitness and mutation class waves propagate. Mutation class is fed by beneficial mutations from the previous class, . The fittest class produces a new class of mutants that will be subjected to genetic drift until its size reaches the establishment size that depends on its relative advantage. Assuming new class senses a constant background from the class with mean fitness, the advantage can be estimated as , where indicates its distance from (see Fig. 1). Under this argument, the establishment size is . During the time in which a new fittest class establishes, the peak of the distribution will advance and the less-fit clones die out, thus keeping constant the distribution width (i.e. constant ).

The speed of adaptation can be estimated by a scaling argument Desai and Fisher (2007); Brunet et al. (2008); Park et al. (2010) imposing that the time during which a new class establishes (and the class histogram moves by one bin) has to correspond to the larger-scale movement of the histogram by deterministic growth.

The first part of the argument estimates from the condition


which imposes that order one mutations are produced (and establish with probability ) from the foremost bin while it grows exponentially according to its advantage.

The foremost bin is treated as it were growing in a constant background with mean fitness (as above) and is taken as the size of the foremost bin at its birth. Integration of the above equation yields the expression


which is valid for sufficiently large and advantage of the fittest class, since the contribution of the integration boundary has been neglected.

The second part of the argument can be formulated as a normalization condition of the fitness class histogram, joint with the condition that each bin grows exponentially according to its advantage,


note that this condition neglects the contribution of beneficial mutations from the neighboring bins. Approximating the above formula as an integral computed with saddle-point (i.e. assuming that the largest term in the sum, the class with mean fitness, dominates Desai and Fisher (2007); Brunet et al. (2008)) yields the expression


Equating the two expressions Eq. (4) (or some variants that can be derived more rigorously) and (6) for or generates implicit estimates for the speed of the mutation classes histogram, , which in general work very well, despite of the approximations taken. In particular, neglecting the addends in Eq. (4) and (6), i.e. up to logarithmic corrections in , one obtains the closed formula Desai and Fisher (2007); Brunet et al. (2008)


Finally, can be estimated by .

In the general case of diminishing return epistasis, the variables of interest become a function of the genetic background. Specifically, the speed of adaptation in log-fitness space and the speed of the mutation classes, are two non-trivially distinct quantities, and the same holds for the width of the mutation and the advantage histograms ( and ) relative to the population (Fig. 1).

ii.3 Simulation algorithm and effective parameters

The model was simulated with the algorithm of Park and Krug, as described in refs Park and Krug (2007); Park et al. (2010). The simulation scheme is sketched in Fig. 1C. In a typical initial configuration, all clones have mutations. At each subsequent time step the progeny of the individuals of all fitness classes are sampled from a multinomial distribution of parameters . The parameters take into account the relative fitness computed at time and the contribution of beneficial mutations. Specifically


(see also Eq. (9)). Together, these definitions are equivalent to a Wright-Fisher model with separate selection and mutation steps. In particular, the microscopic step described above reduces to the standard Wright-Fisher model for genetic drift when . The multinomial random numbers are generated by iteratively drawing binomial random numbers with parameters starting from .

The model defined above is invariant by suitable rescaling of time, provided the other model parameters are also rescaled correctly. This feature is useful in comparison with experiments (see following) in order to understand the correspondence between time steps in the model and a generation in the experiment. It can also be useful in increasing the efficiency of simulations.

Let us suppose one desires to map a reference empirical time into the time steps of the model. This requires the mapping , in the expression describing the growth of a sub-colony from a clone in terms of the experimental generation time. The coefficient is the number of empirical generations corresponding to a time step, and represents the simulation time scale. A simple choice of time scale is the case , implying a one-to-one correspondence between empirical generations and time steps in the model. Since an established clone with advantage grows as , the advantage function is proportional to the time scale as follows, . This means that the value of the fitness in the model depends exponentially on . Similarly, since the beneficial mutation rate is defined as the number of expected mutations per genome per generation, the map between time steps and empirical generations implies . Finally, the correct rescaled population size is . In practice can be rather large in a typical experiment (e.g. ), so that the rescaling of population size does not affect much the dynamics provided is not too large. The rescaling of the parameters described above can easily be rationalized keeping in mind that the basic time scales of the model are set by the products and . In order to verify that the invariance discussed above is valid for our model, where the advantage varies with the number of mutations as , we ran some simulations choosing the model parameters (, , ) using different maps between the time units (see Appendix, Fig. A.1). The results indicate that for any practical purpose, the invariance is effective.

Iii Results

iii.1 Adaptation slows down until it reaches an effective arrest.

Figure 2: (Color online) The mean number of mutations and the mean advantage grow sub-linearly with time. The top panels show the increase in time of (A) and of (B) obtained by direct simulation of the diminishing return model. The plot in panel B is in log-log scale, and the data are compared with a reference line (dashed blue line, with slope ) to highlight the sublinear growth of . The continuous red line shows the the asymptotic long-time linear behavior with slope corresponding to . (C) Long-time behavior of the mean speed of fixed mutations (green symbols), averaged over different realizations. For long times, this quantity decreases (as a power law) towards the limit value , where the assumptions of the model break down and deleterious mutations need to be accounted for Rouzine et al. (2008). This limit also corresponds to the limit value of obtained by a mean-field estimate (see text). Simulations are carried out using the parameters , , , . Averages are computed over realizations (these averages are implied in the notations for the y-axis labels).

Direct simulation of the model (Fig. 2A), shows that the mean advantage grows sub-linearly with time, as expected from the diminishing returns pattern of the advantage. The average number of mutations also grows sub-linearly with time, for intermediate to long times (Fig. 2B). This trend is independent from (or from the specific model of the of decreasing and is due to the fact that decreasing advantage and the consequent rise of the establishment threshold for clones together slow down adaptation. The time derivatives of and estimate the adaptation speed and the mutation-accumulation speed of a typical realization. Figure 2C shows an average of over one hundred realizations, plotted as a function of time. The simulations indicate that relaxes to a plateau which is close to the beneficial mutation rate . Equivalently, for long times, the mean number of fixed mutations shows a linear behavior in time with a rate close to (red line in Fig. 2B). In the same long-time limit, the advantage of a mutation drops asymptotically to zero.

The long-time regime where follows a linear trend is outside of the limit of validity of the model, and has to be regarded as unphysical, since when , deleterious mutations cannot be neglected Rouzine et al. (2008). Thus, for any finite the asymptotic trend of has to be interpreted as an effective signature of an evolutionary arrest for both and , where beneficial mutations should be in equilibrium with deleterious one, which could possibly be captured by a variant of the model including deleterious mutations Rouzine et al. (2008). Note also that realistically the beneficial mutation rate itself could decrease in the later stages of evolution Park et al. (2010); Park and Krug (2008).

Despite the fact that the assumptions of the model break down asymptotically, the long-time phenomenology appears to be theoretically consistent. Indeed, the long-time regime can be easily rationalized as a transition to an “effectively neutral” evolution regime. In fact, the selection coefficient becomes irrelevant at long times, so that the fixation dynamics is driven solely by genetic drift. Given that deleterious mutations are not accounted for, and the benefit carried by the allowed mutations becomes negligible, in this regime effectively becomes a neutral mutation rate. The probability of fixation of an essentially neutral mutation is while the rate of appearance of new mutations is . Therefore, the pace at which new mutations are accumulated is approximately . As we will see, this asymptotic velocity is also recovered from a mean-field estimate of , valid in the infinite population limit.

iii.2 A mean-field analysis gives a relation between the variance of the advantage distribution and the adaptation speed.

Figure 3: (Color online) The mean field approximation captures a relation between and the width of the fitness class distribution, which is valid at intermediate times for moderate and until longer times for large population sizes. (A) Simulated variance of the mutation classes distribution, shown as a function of the expected variance from the mean-field estimate (, see Eq (11)). To avoid ambiguities, averages over realizations are indicated by a suffix . The continuous red line represents the theoretical prediction . The error bars (standard deviations over realizations, ) become larger with . (B) While for increasing times diverges, the relative variability over realizations decreases with increasing population size , for any fixed time. This suggests that the infinite mean-field estimate is well-defined. Simulations are carried out using the parameters , , . Population size in panel A is .

In the limit of infinite population, the dynamics of the model can be described using the following mean-field equation Park et al. (2010),


where is the frequency of individuals with beneficial mutations at generation , , and the mean fitness.

Multiplying Eq. (9) by and summing over gives the following expression for the dynamics of the mean number of mutations ,


This expression can be further simplified assuming that the frequency distribution is narrowly peaked around the mean (which travels in time) and that it can be expressed as . In this case it can be easily verified that .

A different, more instructive, relation, which keeps into account the width of can be obtained starting from Eq. 10, and expanding under the assumption that , for every index of non-empty classes. This assumption is verified by simulations (see Fig. A.3 in Appendix) and by further considerations on the finite- width of the distribution given in the following sections.

Computing averages in Eq. (10) to first order in , and noticing that and , it is possible to obtain an expression for the speed of accumulated mutations as a function of the variance of the fitness class distribution. Estimating as gives, for the case


According to this equation, is driven by two terms, the increase of mutations due to the beneficial mutation rate and the selection of individuals with larger fitness. This result is somewhat reminiscent of Fisher’s fundamental theorem and the Guess relation, which relate the speed of adaptation to the variance of the fitness. In this case, the speed of accumulation of successive mutations is related to the width of the mutation class histogram, but rescaled by the factor , which decreases with time. For the non-epistatic case () one recovers the more usual linear proportionality, since the advantage is linear in the mutation class index Park et al. (2010).

The mean-field limit of the model with has been previously addressed by Park et al Park et al. (2010) with a moment generating function approach. In particular, they estimated the distribution variance as , which substituted in Eq. 11 gives precisely their expression for the speed (in the limit of small ) suggesting that our result is a consistent generalization to the epistatic case.

Therefore, for the speed of adaptation does not depend on the mutation rate for infinite populations. On the other hand, for diminishing returns, the increase in the width of the distribution of , does not compensate for the term (which tends to 0), and, for long times, Eq. (11) predicts the limit velocity , as observed in simulations (Fig. 2C). Since the increase in the width of the distribution depends on the size of the population and on the advantage , the unphysical limit velocity is approached with different laws depending on the specification of .

It is necessary to discuss under which conditions this mean-field limit can be considered a valid estimate of the typical behavior of a realization. At finite , simulated data are in good accordance for intermediate time with the mean field estimate for the variance of the mutation class distribution (Fig. 3). Additionally, decreases quickly with time (see Fig. A.3 in Appendix), justifying the assumption of small (since it is expected that for every ). However, the variability of over different realizations increases quickly. In order to verify whether these fluctuations are well-behaved, we have evaluated , where indicates the variance of the mutation class distribution in a single realization (roughly analogous to ), and the suffix indicates averages over realizations. Specifically, indicates the average of the quantity over different realizations, while is its standard deviation. Therefore, , plotted in Fig. 3B as a function of , represents the relative variability over the realizations of the variance of the distribution. For any fixed time, this quantity decreases with , suggesting that the mean-field limit is well-defined for infinite populations. Conversely, fixing and increasing , appears to reach finite values, hence (and hence ) seems to be non-self averaging in time. This indicates that a mean-field description of the population dynamics might be appropriate for longer time-scales at increasing population size.

Note that these questions are superfluous for any empirical purpose of the model, since the infinite mean-field limit does not describe correctly the finite dynamics already for  Park et al. (2010), and the time scale of any experiment is orders of magnitude smaller than those considered here. “Real” evolutionary time scales are longer, but presumably subject to changing conditions, and therefore beyond the scopes of this model. However, the infinite limit is instructive on purely theoretical grounds, for understanding the main mechanisms driving the model.

To sum up, the mean-field regime predicts that the infinite-population/infinite-time speed for the mutation classes is . At finite , the long-time averages for the number of mutations and the associated speed agree with the predictions of the mean-field regime, but are characterized by strong fluctuations. Finally, an extended expression accounting for the finite width of the class relates to the width itself, and is verified by simulations. In this long-time regime, a number of assumptions of the model break down, and more complex models need to be employed. For practical purposes, neither the infinite population size nor the infinite time limit are relevant, as even in long-term experiments the number of generations of interest is typically sufficiently low. This will be discussed in a later section, when the model parameters are matched with empirical data from two laboratory evolution experiments. The next sections discuss in greater detail the behavior of the model at finite .

iii.3 At finite population size, the advantage and fitness-class distributions remain Gaussian, but their width is unsteady, with opposed trends.

The large- long-time phenomenology described so far and captured in terms of mean-field estimates is not very relevant empirically, given the experimental range of times and population sizes. In the following we offer a more detailed summary of the behavior of finite-size populations at intermediate times (i.e. on the relevant experimental time scale generations, with ).

As already illustrated, the frequency distributions of the mutation and advantage classes are approximatively Gaussian for a constant advantage model. Direct simulations of the diminishing return model show that this still holds (Fig. 4A). Given the discrete nature of these distributions, their widths are well represented by the distances and of the foremost bin from the average (shown in Fig. 1A). While in a fixed advantage model these distances, in terms of both mutation and advantage classes, between the edge and the mean are constant, for a diminishing return model is an increasing function of the mean number of mutation classes , while decreases with (Fig. 4B).

The model without epistasis (and the mean-field estimates) make us expect that the widths of these histograms are related to the speed of adaptation and consequently of mutation accumulation. The instantaneous velocities of the two distributions and both decrease with .

Qualitatively, one expects that, with changing mean number of mutations, the two speeds are connected by the decay of the advantage between and . However, the relation between them is also affected by the change in the width of the distributions with increasing . Indeed decreases more rapidly and, in the long-time limit, even if , it vanishes. This implies the existence of a large number of sub-populations with different but essentially the same fitness, which leads to the effectively neutral behavior discussed in the previous section. Since the advantage decreases with time, and the distributions change width, it is interesting to ask under which conditions the biologically relevant regimes are preserved. Specifically, the condition of multiple concurrent mutations could be lost over time.

As discussed in the Introduction, the clonal interference condition can be roughly defined (neglecting irrelevant numerical factors) by the scaling relation , where is the advantage of the fittest class. This condition is satisfied when the population size and the edge advantage are large enough, i.e. when and . Considering a fixed advantage model, for typical values of experimental population size () and advantage (), these inequalities are always satisfied, for estimated in the range mutations per genome per generation. However, when the contribution of a newly arising mutation becomes too small, and the advantage in fitness of the edge () is smaller than , the clonal interference hypothesis breaks down. Note that this implies also an establishment size larger than . These considerations confirm that, while the long-time limit of this model is biologically implausible, the experimentally relevant parameter regime should be captured correctly by our theoretical framework.

We will now attempt to generalize the analytical estimates for the finite- adaptation speed available for the multiple mutation model to the case of diminishing returns. These “stochastic edge” estimates are based on the hypothesis that the only class subjected to substantial stochastic effects is the fittest one, i.e. that  Rouzine et al. (2008); Desai and Fisher (2007). When the single stochastic class approximation fails, and the dynamics is more complex. Supposing that , (this is a much stronger epistatic effect than the one we estimate from experimental data, see the following) and that (of the order of the values obtained from simulations with this parameter set and population size ), the advantage becomes close to a beneficial mutation rate of (which can be considered very large Hegreness et al. (2006); Perfeito et al. (2007)), for . This exceeds the interesting experimental range of (), suggesting that for the relevant range we can always suppose that the stochastic edge approximation is valid.

iii.4 Analytical estimate for the adaptation speed at finite .

Figure 4: (Color online) The histograms of fitness advantage and mutation classes have nearly Gaussian forms. While adaptation slows down, the latter histogram expands while the former becomes increasingly peaked. (A) Histograms of the mutation classes (top) and advantage classes (bottom) obtained from simulations averaging over 200 realizations at different generations (different symbols, see legend). The parabolic form in the semi-log plot indicates that they are approximately Gaussian. The establishment size is represented as a dashed line in the top panel. Solid lines connecting the symbols are guides to the eye. (B) Simulated data for the widths of the mutation class histogram (green squares, top) and of the fitness advantage histogram, (blue circles, bottom), plotted as a function of the mean number of mutations . The continuous line represents the theoretical estimates of the width (see Eq. (21)). (C) Plots of the speed of mutation accumulation (top, green squares) and of adaptation (bottom, blue circles), as a function of the mean number of mutations . Continuous lines are the corresponding theoretical estimates (see Eq.(16) and (17)). Note that since depends logarithmically on , the estimates for both speeds are in satisfactory agreement with the simulated data even if is approximated more roughly. The parameters used in the simulations are , compatible with those estimated (see Sec.III.5). Averages are performed over 200 realizations.

Since we have verified that the advantage and mutation class histograms are both nearly Gaussian (but not stationary in width), it is possible to attempt a generalization of the estimates applied for the standard multiple mutation model, Eqs. (3) and (5). Supposing a slow increase of the width of the distribution with and assuming that the width of the histogram is stable while the new class is being established, the mean of the distribution moves from to during the establishment time and . For (and not too high due to the condition of sufficiently large ) the advantage of the edge respect to the mean class is .

Thus, imposing the condition that one new mutation class is established at the edge of the histogram gives


As in Eq. (3), the last term in this integral is the establishment probability of a new fittest class, while the first is the rate of beneficial mutations from the previous fitness class, which is born with size and grows exponentially. Note that imposing , Eq. (12) becomes Eq. (3), recovering the standard multiple-mutation model.

An estimate of can be obtained integrating the above expression


where we used the approximation . Under the assumption that is not too large (i.e. the advantage of the fittest class is sufficiently high) the contribution of the integration boundary can be neglected and the expression of is the equivalent of Eq. (4) for the non-epistatic case,


The above expression can be further simplified assuming that both and are large enough so that and and neglecting the logarithmic term in leading to the following expression


This indicates that for intermediate , the estimate of the standard multiple mutation model is valid provided the diminishing return advantage function is substituted to the constant advantage. For sufficiently large , expansion of the exponential in Eq (13) gives , compatible with the mean-field result.

The time is related to the instantaneous speed of the mutation class histogram . The speed , can be obtained knowing that during the time the mutation class histogram travels by one class, the advantage histogram has to move by the relative fitness between the newly added fittest class and the previous one, . Using the simplified Eq. (15) (which neglects logarithmic terms in ), one obtains the speed




The second part of the estimate involves the normalization condition (5). Assuming (as in the standard estimate) that the largest term of the histogram dominates, we need to evaluate the time necessary for the fittest class to become the class with mean advantage, whose size is order . If the fittest class has mutations, its establishment size is . Its growth will be roughly exponential, with a rate that decreases while it gets closer to the mean. We estimate its growth by its mean growth rate during the time . Immediately after establishment, its relative growth rate will be , while its rate will tend to zero when it gets close to the mean. Thus, on average, we can assume that it grows exponentially with rate .

This argument leads to the equation


which corresponds to Eq. (5), and implies the equivalent of Eq. (6),


In order to estimate , we need to determine how much the histogram of fitness advantage has progressed during the time from the establishment of the -th mutation class. We assume that during this time is roughly constant, so that after time , mutations are established, and the advantage of the edge has reached .

This allows to estimate as the advantage gained divided by the time , i.e.


Eq. (17) and (20) together allow to determine , which can subsequently be used to obtain the speed of adaptation , or of fixed mutations , using e.g. Eq. (17) and (16). Assuming that and neglecting the logarithmic corrections in , as made to obtain (20), we can obtain the following closed expression for ,


Comparison with simulated data shows that the above expression for , despite the rough approximations made, is a reasonably good estimate of the width of the distribution (Fig. 4B). In particular, the speed of adaptation and mutation accumulation (Fig. 4C) are well captured by our analytical description. We should stress again that these expressions can be considered valid for intermediate values of , as explained above. For large , the integration boundary in Eq. (13) cannot be neglected, and it can be verified that tends to a constant, with the approximations taken, restoring the correct mean-field limit.

iii.5 Parameter matching of the model with data from laboratory evolution experiments.

Having explored some of the main features of the diminishing return model using theoretical arguments and simulations, we now proceed to compare it to experimental data. We have made very clear that this model is a very crude description of any realistic situation. However, its advantage is that it is simple and based on few parameters, so that estimating some of the model’s fundamental parameters from data is a relatively simple task, which could be instructive Desai et al. (2007).

We developed a parameter matching procedure that produces an estimate of the beneficial mutation rate. This estimate depends on the specific advantage model chosen, but also allows to select models according to how well they resemble the data. Comparing the simulated dynamics with the experimental one, it is possible, for each specific data set, to determine an optimal functional form of the advantage , thus narrowing down the estimate of .

We analyzed fitness/mutations data from two laboratory evolution experiments using the three different variants of the diminishing return model described in Appendix A. The power law model, where the advantage is described by , is the main case presented in the former sections. The advantage functions of other two models considered have a logarithmic () and an exponential () dependence on . These functions can be derived, with some approximations, as partial sum of the harmonic and geometric series, which allows to produce simple expressions of the antagonistic effect of each added mutation. These two model variants share most of the qualitative features of the main formulation.

iii.5.1 Experiments analyzed

We applied our model to two experiments of controlled evolution. The first are the first 20000 generations of the well-known “E. coli long-term evolution” experiment Barrick et al. (2009); LTE (); Lenski et al. (1991), while the second is a chemostat experiment performed by some of the authors Jezequel et al. (Article accepted).

The two experiments concern two bacteria (respectively Escherichia coli and Acinetobacter baylyi), were performed using distinct propagation techniques (serial dilution in batch and chemostat respectively) and also their duration is quite different both in terms of generations ( and ) and experimental time ( years and months). Another remarkable difference between the two experiments is the population size, which varies every day by two orders of magnitude (between ) for the serial dilution experiment while it is very large and approximately fixed in the A. baylyi experiment (). The two experiments share interesting features suggesting that a diminishing return model might be applicable to describe the evolutionary dynamics of the populations. Firstly, their duration in terms of generations is long enough to observe a deceleration of fitness increase Kryazhimskiy et al. (2009). Secondly, the large (effective) population size suggests that the clonal interference regime might be relevant in these experiments. The simultaneous presence of different genotypes within the population has been verified in both experiments Jezequel et al. (Article accepted); Elena and Lenski (1997). Moreover, the decrease of the beneficial effect of the first five fixed mutations in the E. coli serial dilution experiment has been recently demonstrated Khan et al. (2011). The effect is more complex than the description of diminishing returns given here. However, as suggested by the authors, one can surmise that a simplified model including epistatic interactions might be useful to roughly describe this phenomenon.

A detailed description of the data and the choices made for the analysis procedure is given in Appendix B. Here we give a brief description of the main features of the two experiments.

The A. baylyi experiment studied the population dynamics in a chemostat using a minimal medium supply for about four months, at a dilution rate . The use of chemostat allowed to grow a large population () under controlled conditions for a fairly long time. Since the number of individuals is large, it is expected that different sub-populations will grow in parallel in clonal interference regime. This has been confirmed by population sequencing data Jezequel et al. (Article accepted). Additionally, several single clones where isolated from samples collected and frozen at different time.

The maximum growth rate () of 21 isolated clones and of the the original strain introduced in the chemostat (wild type) has been measured in batch, fitting the growth curve during the exponential phase. We considered the maximum growth rate measured in batch as indicative of the fitness of the population into the chemostat, defining the normalized fitness of the clone as . A more detailed description of this procedure is given in App. B. Note that, since the reference fitness value is given by the ancestral growth rate, . A whole-genome sequencing was performed on two single clones isolated at the end of the experiment (AB2800b and AB2800a) and on the ancestral strain that served as a reference for identifying mutations in evolved clones. A total of 11 mutations were detected into the evolved clones, eight of them in common between the two. Additional sequencing has been performed in the remaining clones selected at different generations on PCR fragments encompassing the mutated loci identified in the two end-point clones in order to reconstruct a sketch of the history (or the genealogy) of the mutation appearances.

The E. coli long-term evolution experiment concerns twelve E. coli populations evolved in parallel for about generations in batch. Serial dilution was performed daily, allowing generations each day and an effective population size of individuals. We used the mutations and fitness data of the population designated Ara-1 referred to the first generations, as given in ref. Barrick et al. (2009). The relative log-fitness of the ancestor and evolved populations(that is essentially the ratio between the Malthusian parameters, ) was measured through competition experiments (see ref. de Visser and Lenski (2002)) every generations. We define in this case the mean fitness of the sub-population as , where is the approximate mean growth rate (see ref. Lenski et al. (1991) and Appendix B). The fitness of the reference strain is again . Genome sequencing was performed on samples from generation , , and as well as on the ancestor. A total of mutations were found in the most evolved strain, most of which were stable in later clones Barrick et al. (2009).

For all the analyzed clones in both experiments, it is possible to associate fitness values with numbers of mutations, and thus bridge with the parameter in our model. For the A. baylyi experiment the estimate of the number of accumulated mutations are inferred from the genealogy, and thus have to be considered as a lower bound for all the clones, except for the ancestor and the two fully sequenced clones where they are measured directly. However, the indications of from the population sequencing data are compatible with the inferred values of  Jezequel et al. (Article accepted).

iii.5.2 Matching procedure and estimate of

Figure 5: (Color online) Sketch of the parameter matching procedure leading to an estimate of the beneficial mutation rate from experimental fitness data. Top: the experimental data for fitness as a function of the number of mutation (blue symbols) allow to estimate the advantage fitness function from a fit (continuous green line). Middle: simulations are run using the estimate advantage function for different values of , which remains undetermined. This leads to different predicted dynamics for the number of acquired mutations (dotted,dashed, and dashed-dotted lines) and the fitness increase in time. These predictions can be matched with the experiment to estimate . In particular, the predicted number of time steps necessary to reach the final experimental number of mutations varies with the beneficial mutation rate (filled circles). Bottom: The estimate of is obtained matching the predicted final time with the experimental one. The best value of (blue star) is obtained when the predicted number of time steps necessary to reach the final number of mutations (determined in the experiment) corresponds to the number of experimental generations (horizontal continuous red line). We applied this procedure to the three model variants for the diminishing return described in the main text.
Figure 6: (Color online) Estimate of the beneficial mutation rate from different experiments. Estimate of the beneficial mutation rate using data from Jezequel et al. Jezequel et al. (Article accepted) (left panels) and Barrick et al. Barrick et al. (2009) (right panels). Top: estimate of the fitness advantage function from data of advantage as a function of mutation number (red symbols), as described in the top panel of Fig. 5. The different lines correspond to the three model variants for the diminishing-return advantage (power-law dot-dashed purple lines, logarithmic long-dashed green lines, geometric blue dashed lines, as in legend). For sake of simplicity for experiment A. baylyi (panel A1) for each only the mean value of the advantage is shown. Complete data are reported in ref. Jezequel et al. (Article accepted). Middle: estimates of , from different model variants, obtained matching the time for which predicted by simulations with a given to the experimental number of generations, as described in the bottom panel of Fig. 5 (the different line styles refer to the model variants as above). For A. baylyi experiment we considered a total of generations. Bottom: comparison of the performance of the three model variants with the estimated values of . For the Jezequel et al. data (left panel), the power law and geometric variants are quite close, but power law model for the diminishing return gives the best agreement, especially considering the data relative to the number of acquired mutations as a function of time (not shown). This choice leads to the estimate . For the Barrick et al. data, the logarithm and power-law models perform better, and give equivalent estimates of .

Figure 5 summarizes the parameter matching procedure. The first step finds the best-fitting parameters for each functional form of the fitness advantage function . This uses data on fitness values and number of acquired mutations, using the definitions given in the previous section. Since the number of experimental points is low ( for the E. coli and for the A. baylyi experiment- corresponding to different values of ), it is possible to obtain good fits with different functional forms of the advantage . The results of the fit using the three different functional forms considered here (power law, logarithm and geometric) is shown in the top panel of Fig. 6. In a second step of the procedure, simulations are repeated for a wide range of values of . An estimate of this parameter can be obtained using the mean number of mutations at the end of the two experiments (Fig. 5). In our case, and for the A. baylyi experiment and E. coli experiment respectively. In the model, for a fixed interval of time steps, the number of accumulated mutations decreases significantly and monotonically with , and the number of time steps necessary to reach and depends on . Thus, for each model there is a single value of the beneficial mutation rate that verifies (Fig. 5 and middle panel of Fig. 6). Note that we are referring to the estimated values of as mean values assuming that the sequenced clones are representative of the populations. For the A. baylyi experiment, the uncertainty on the number of mutations present in the clones is a relevant source of error. On the other side, even if for the E. coli experiment the number of mutations is exact, and referred to single clones, the fitness data we associated are population mean fitness data, that can introduce an error in the estimate.

The procedure applied so far gives some estimated values of the beneficial mutation rate, which however,depend on the advantage model, varying at this stage up to three orders of magnitude. Nevertheless, all the estimates are roughly in the expected biological range (, see Ref. Perfeito et al. (2007); Hegreness et al. (2006)). The third step of the estimate procedure allows to select between advantage model, and is based on the comparison of simulated and experimentally measured dynamics for the fitness and the number of mutations as a function of time. Note that this is not trivially equivalent to the fitness as a function of , but contains the effects of the population dynamics in presence of clonal interference, for which the model provides a description. Since different functional forms of the advantage imply a different behavior of the population, one can check which of the advantage functions gives the closest description to the experiments. Thus, the experimental data of the advantage as a function of the time can be used to discern between power law, logarithmic and geometric model, ultimately selecting an estimate of .

A qualitative comparison between data and simulations shows that for the A. baylyi experiment the power law and geometric advantage model (with parameters , and , respectively) best describe the increase of the fitness with time (bottom panel of Fig. 6). Comparing the increase of the number of mutations as a function of time of the two models suggests that the power law model better resembles the data. For this reason, we speculate that, despite of the large errors in the fits, the power law model could perform better in this parameter-estimation procedure and the preferred estimated value of beneficial mutation rate is around the value . In the E. coli experiment, the power-law and logarithm model result describe the data best, and are roughly equivalent. The estimate parameters are , for the power law and for the logarithm model. The estimates for the beneficial mutation rate are also essentially equivalent (considering the errors connected to the procedure) for the two preferred models for the power law and for the logarithmic advantage model. Note that for this step of the procedure the data used are much more abundant, since we used all the fitness data measured every generations. A logarithmic increase of the fitness advantage for E. coli long term evolution experiment is also suggested by a parallel work Wielgoss et al. (In press).

Iv Discussion and conclusions

Different laboratory evolution experiments show a decrease of the fitness advantage due to newly acquired mutations and a decrease of the speed of evolution Elena and Lenski (2003); Kryazhimskiy et al. (2009); Hindré et al. (2012). This effect can be explained in different ways, and accordingly different models have been formulated in this context. For example, the speed of evolution could decrease because beneficial mutations with larger advantage fix sooner in the population Schiffels et al. (2011) and because the mutation rate or the number of possible beneficial mutations decreases with time Rouzine et al. (2008); Park and Krug (2008). Another explanation, suggested by recent experimental observations Khan et al. (2011); Chou et al. (2011), could be that increasing the number of accumulated mutations, epistatic interactions cause a decrease of their effect on the fitness. These different explanations are not necessary mutually exclusive, and could be stratified in actual laboratory evolution experiments Tenaillon et al. (2012).

We considered a simplified model, using a minimal number of parameters, which is a direct generalization of the multiple mutations model at constant advantage, but describes diminishing returns. Specifically, it is assumed that the selective advantage of all individuals having beneficial mutations is identical, but decreases with . While more complex and realistic descriptions applicable to laboratory evolution exist an advantage of the approach taken here is that, as the multiple mutation model, it depends on few parameters, and as we have shown, in line of principle allows matching of these parameters with data. On theoretical grounds, this class of models is interesting because its dynamics is driven by rare events Desai and Fisher (2007); Brunet et al. (2008), and because its behavior at finite is qualitatively distinct from its mean-field behavior Park et al. (2010).

We have shown that basic phenomenology of the model entails a sublinear decrease of the mean number of fixed mutations and a steeper sublinear decrease of the mean advantage. This is in qualitative agreement with previous results using a similar model applicable in a regime where concurrent mutations do not occur Kryazhimskiy et al. (2009).

The two evolutionary speeds are related to the width of the distribution of coexisting advantage classes, and thus of coexisting mutation classes. We showed how a theoretical mean-field argument produces a relation between the speed of fixed mutations and the second moment of the histogram of mutation classes, confirmed by simulations. In the limit case of constant advantage, this relation reduces to the one previously derived in ref. Park et al. (2010), related to Fisher’s fundamental theorem. Interestingly, simulations indicate that for any finite , different model realizations behave increasingly differently with time in terms of both and width of the mutation class histogram. This non-self-averaging property implies that even at intermediate times, the behavior of a realization can be quite different from the average.

Within a mean-field framework, we also derived arguments showing that the mutation speed approaches the beneficial mutation rate for very large times, which corresponds to the mean behavior observed in simulations at large times. In fact, as the original multiple mutations model, the present model does not include deleterious mutations. However, in the case of diminishing returns the assumption is more delicate, because the advantage of fixed mutations will decrease indefinitely, until no beneficial mutation is able to fix. Empirically, before this happens, beneficial mutations will be balanced by deleterious ones Rouzine et al. (2003). This makes the (effectively neutral) long-time limit of the model unphysical. We have discussed how this limitation should not affect time-scales and parameter values that are relevant for current laboratory evolution experiments focused on adaptation Hindré et al. (2012). A generalization of the model that includes deleterious mutations is possible and could be the subject of future investigations, since it allows to formulate evolutionary questions related to the balance of mutations of different kinds Rouzine et al. (2003, 2008).

Finally, for finite population size , we were able to define through analytical arguments the regime where the stochastic edge estimate of the adaptation and mutation speeds can be extended to the case of diminishing returns. In the simplest case it is possible to obtain closed expressions for the mutation class and adaptation speed and for the supports and of the mutation class and advantage histogram respectively. These expressions are direct generalizations of those obtained in the case of constant advantage Desai and Fisher (2007); Brunet et al. (2008), and compare well with numerical simulations. Once again, in presence of diminishing returns these expressions, besides the usual limitations in terms of parameter values, have a limited applicability in time, but are typically valid in the experimental range.

The constant-advantage multiple mutation model has previously been applied to short-term laboratory evolution experiments Desai et al. (2007). In those early stages, adaptation does not slow down, and the assumption of constant advantage is justified, allowing to use the model to estimate empirical parameters. However, on longer but experimentally observable time scales, the assumption breaks down, and arguably the closest extension of that model that can be used to estimate parameters is the one described here.

In order to show a proof-of-principle parameter-matching procedure for the diminishing return model, we considered two different experimental data sets from long-term evolution experiments, and defined a procedure that can be used to infer the order of magnitude of the beneficial mutation rate from genomic and fitness data, assuming the model. This procedure allows to estimate the beneficial mutation rate and also compare different functional forms of the decreasing advantage functions against data, through simulations of the model. The values obtained for the beneficial mutation rate fall within the range of the available measurements Perfeito et al. (2007); Hegreness et al. (2006), on the order of mutations per genome per generation for the E. coli long-term evolution experiment, and between and in the case of A. baylyi.

MO and MCL acknowledge support from the International Human Frontier Science Program Organization (Grant RGY0069/2009-C). The work of MRF, MCL, FH, and PT was supported by a “Convergence” grant from the University Pierre and Marie Curie, Paris. We thank G. Malaguti, L. Peliti and S. Wielgoss for the useful comments and discussions on this work.

Appendix A Definition of the three model variants

This Appendix describes in additional detail the three model variants used in the main text.

All variants start from the main assumption that the selection coefficient is dependent on the number of mutations , where is a decreasing function of , and corresponds to different specifications of . Every individual with beneficial mutation has fitness


and we can always arbitrarily define as the “wild-type” fitness. The main example of we considered is given by the choice of a fitness gain that depends on the number of mutations occurred as a power law. In this case, with , where the epistasis grows in strength decreasing from the non-epistatic case of . In this case, the fitness is


with for diminishing return, and for no epistasis.

In the general case , one can write


where indicates the Riemann zeta function. Thus, the relative fitness can be expressed as


For the particular case ,


which is obtained truncating the harmonic number expansion, and where is the Euler-Mascheroni constant. This approximation neglects the terms . Under this assumption, the relative fitness can be expressed as


These expressions show that it is essentially equivalent to directly assume fitness functions of the form


that can be expressed, neglecting terms of order , as sum of powers of the number of accumulated mutations. In the case the factor results a multiplicative factor that can be canceled since it does not affect the relative fitness, that is the relevant quantity for the dynamics. On the other hand, this notation requires an appropriate rescaling of in the case of . Under these assumptions, the relative fitness, which ultimately defines the dynamics, is identical to the one obtained considering the explicit sum of the contributions of single mutations, neglecting corrections for small . Moreover, the fitness functions 28, using instead of , allow to automatically include the case with the correct normalization condition . We refer to the case as the power law model, while the case corresponds to the logarithmic model.

In the case of an approximately Gaussian distribution of mutation classes with distance between the fittest class and the class with mean log-fitness, the relative fitness of the fittest class can be expressed as


where the expansion is performed for which should be the case in our system. In the main text, we use this approximation to perform the estimate of the adaptation speed. In the non-epistatic case () the expression gives the relative advantage expression of the fittest class used in the standard multiple mutation model. Similarly, in the limit


The third model we use is characterized by a “geometric dependence” of the fitness advantage on the number of acquired mutations. In this case the advantage function is with and the fitness is given by


In other words, the advantage accumulates following a geometric sum. As for the former models, considering only the final form of the advantage function , allows to directly satisfy the condition . The constant factor could be adsorbed in , as for the power law model. In the same quasi-Gaussian approximation described for the power law and logarithmic model, the relative fitness is given by


and, in the case of the fittest class it becomes


Note that, while in the previous models the fitness is not bounded from above, the fitness for the geometric model is finite for infinite (i.e. ).

Figure A.1: (Color online)The model is effectively invariant for rescaling of the parameters. The figure shows the mean number of mutations (, panel A) and the fitness advantage (panel B) obtained from simulations run using three different rescaling factors (=1,10,100 different symbols as in legend). The simulated dynamics is almost unaffected by the rescaling procedure. Simulations are made using the parameters , , , . The data are averaged over iterations, and error bars (standard error) are smaller than symbols.
Figure A.2: The fixation probability is proportional to the fitness advantage. The figure shows the fixation probability of a single clone that grows in a uniform background having advantage (symbols), measured from our simulations. The fixation probability is obtained as the ratio between the number of realizations where the beneficial mutator fixes and the total number of realizations (). The results give (continuous red line), in accordance with Park and Krug (2007). Simulations are performed using .
Figure A.3: The relative variance of the distribution decreases in time, while the relative variance of the velocity follows an increasing trend. The figure shows the ratio between the variance and the mean number of mutations (left panel) and the ratio between (calculated over different realizations) and the velocity (right panel) as a function of time. Simulations results confirm that, despite increases in time, the mean number of mutations increases more rapidly. Since , simulations results confirm the hypothesis used in the main text (see also Fig. 4B). The increase in variance of the velocity reflects both the increase of the variance of the distribution and the decrease of . However the error is sufficiently small not to affect the results on experimentally relevant time scales. Simulations are performed using the parameters , , , . Data are averaged over realizations of the process.

Appendix B Experimental data and definition of fitness

This section briefly summarizes the main features of the evolutionary experiments with A. baylyi and E. coli considered here, in order to describe the assumptions made to link genomic and growth rate measurements to fitness in the model. More detailed information about the long-term evolution experiment with E. coli, and in particular about fitness measurements can be found in refs. Lenski et al. (1991); LTE (). Fitness data used during the parameter-matching procedure are reported in ref. de Visser and Lenski (2002); Barrick et al. (2009). A comprehensive description of the experimental methods used for the experiments with A. baylyi can be found in ref. Jezequel et al. (Article accepted). Note that in our modeling framework, as in most standard evolutionary models, the population size is kept constant and all individuals are substituted by newly generated offsprings at every generation, assuming a constant time interval between generations, although none of these assumptions are completely verified in the two experiments we considered.

b.1 Acinetobacter baylyi evolution experiment

In the A. baylyi experiment the population size is almost constant and, even if generations are not synchronous, the mean growth rate is held fixed by the dilution rate in the chemostat. Maximum growth-rate measurements were performed in batch on single-clone colonies. These values differ from the effective growth rate that each population can reach in chemostat because of competition between different sub-populations for a limited amount of nutrients. However, the measured values of are supposed to be indicative of the effective fitness of the different sub-populations inside the chemostat.

When a monoclonal population grows in the chemostat, its growth rate equals the dilution rate , while in presence of different sub-populations the dynamics becomes more complex Dykhuizen and Hartl (1983). For the sake of simplicity, and since the detailed dynamics is not experimentally accessible, we assume that the growth rate of all clones is fixed by dilution rate in the chemostat over the whole experiment and that the generation time is defined as . Thus, the different maximum growth rates of clones measured in batch are interpreted as different survival probabilities of their offspring inside the chemostat.

Defining the fitness as the mean expected number of (surviving) offspring per generation inside the chemostat we assume that it is estimated by the maximum growth rate of mutant, i.e. , where the index indicates the sub-population, is the growth rate expressed in , experimental data of have the units of , and (). In simple words the fitness of an individual is defined as the amount of offspring assuming it can grow by the maximum growth rate over an average generation defined by the chemostat dilution (and thus fitness is a dimensionless quantity). Since the relevant quantity of the dynamics is the relative fitness, during the parameter-matching procedure we used the normalized fitness () to infer the functional form of the advantage.

b.2 Escherichia coli evolution experiment

The long-term E. coli evolution experiment is performed in batch, and serial dilution is performed daily. The maximum population size () is fixed by the total amount of nutrient in the medium. The number of doubling for each day is as derived from the 100-fold daily increase Lenski et al. (1991). The effective population size can be estimated around individuals. Since the speed of evolution depends logarithmically on the population size, the error on the estimate of given by this approximation should be very low. Competition experiments are performed between samples at different times and a spontaneous mutant of the ancestor, which is easier to track visually and has been verified to have almost the same fitness de Visser and Lenski (2002). Relative log-fitness of population respect to the wild type is defined as


where and indicate the frequencies at the beginning and the end of the competition experiment de Visser and Lenski (2002); Barrick et al. (2009). can be seen as the ratio between the growth rate of the two populations and corresponds to log-fitness in the model. Note again that the experimental data are dimensionless. The difference between the Malthusian parameters can be approximatively deduced as where is the mean Malthusian parameter Lenski et al. (1991). Since we are interested in expressing the fitness as the mean number of offspring per generation, we use as mean Malthusian parameter . Thus we define the normalized fitness as .


  • Park et al. (2010) S.-C. Park, D. Simon, and J. Krug, Journal of Statistical Physics 138, 381 (2010), ISSN 0022-4715, 10.1007/s10955-009-9915-x, URL http://dx.doi.org/10.1007/s10955-009-9915-x.
  • Hindré et al. (2012) T. Hindré, C. Knibbe, G. Beslon, and D. Schneider, Nat Rev Microbiol 10, 352 (2012), URL http://dx.doi.org/10.1038/nrmicro2750.
  • Barrick et al. (2009) J. E. Barrick, D. S. Yu, S. H. Yoon, H. Jeong, T. K. Oh, D. Schneider, R. E. Lenski, and J. F. Kim, Nature 461, 1243 (2009), URL http://dx.doi.org/10.1038/nature08480.
  • Tenaillon et al. (2012) O. Tenaillon, A. Rodríguez-Verdugo, R. L. Gaut, P. McDonald, A. F. Bennett, A. D. Long, and B. S. Gaut, Science 335, 457 (2012), URL http://dx.doi.org/10.1126/science.1212986.
  • Felsenstein (1974) J. Felsenstein, Genetics 78, 737 (1974).
  • Orr (2003) H. A. Orr, Genetics 163, 1519 (2003).
  • Gerrish and Lenski (1998) P. J. Gerrish and R. E. Lenski, Genetica 102-103, 127 (1998).
  • Wilke (2004) C. O. Wilke, Genetics 167, 2045 (2004), URL http://dx.doi.org/10.1534/genetics.104.027136.
  • Tsimring et al. (1996) L. S. Tsimring, H. Levine, and D. A. Kessler, Phys Rev Lett 76, 4440 (1996).
  • Desai and Fisher (2007) M. M. Desai and D. S. Fisher, Genetics 176, 1759 (2007), URL http://dx.doi.org/10.1534/genetics.106.067678.
  • Brunet et al. (2008) E. Brunet, I. M. Rouzine, and C. O. Wilke, Genetics 179, 603 (2008), URL http://dx.doi.org/10.1534/genetics.107.079319.
  • Good et al. (2012) B. H. Good, I. M. Rouzine, D. J. Balick, O. Hallatschek, and M. M. Desai, Proc Natl Acad Sci U S A 109, 4950 (2012), URL http://dx.doi.org/10.1073/pnas.1119910109.
  • Desai et al. (2007) M. M. Desai, D. S. Fisher, and A. W. Murray, Curr Biol 17, 385 (2007), URL http://dx.doi.org/10.1016/j.cub.2007.01.072.
  • Khan et al. (2011) A. I. Khan, D. M. Dinh, D. Schneider, R. E. Lenski, and T. F. Cooper, Sci 332, 1193 (2011).
  • Chou et al. (2011) H.-H. Chou, H.-C. Chiu, N. F. Delaney, D. Segrè, and C. J. Marx, Science 332, 1190 (2011), URL http://dx.doi.org/10.1126/science.1203799.
  • Kryazhimskiy et al. (2009) S. Kryazhimskiy, G. Tkacik, and J. B. Plotkin, Proc. Natl. Acad. Sci. 106, 18638 (2009), URL http://dx.doi.org/10.1073/pnas.0905497106.
  • Martin et al. (2007) G. Martin, S. F. Elena, and T. Lenormand, Nat Genet 39, 555 (2007).
  • Schiffels et al. (2011) S. Schiffels, G. Szöllösi, V. Mustonen, and M. Lässig, Genetics 189, 1361 (2011).
  • Wright (1931) S. Wright, Genetics 16, 97–159 (1931).
  • Fisher (1930) R. Fisher, The genetical theory of natural selection (Clarendon Press, Oxford, 1930).
  • Park and Krug (2007) S.-C. Park and J. Krug, Proc. Natl. Acad. Sci. 104, 18135 (2007), URL http://dx.doi.org/10.1073/pnas.0705778104.
  • Hegreness et al. (2006) M. Hegreness, N. Shoresh, D. Hartl, and R. Kishony, Science 311, 1615 (2006), URL http://dx.doi.org/10.1126/science.1122469.
  • Perfeito et al. (2007) L. Perfeito, L. Fernandes, C. Mota, and I. Gordo, Science 317, 813 (2007), URL http://dx.doi.org/10.1126/science.1142284.
  • Eyre-Walker and Keightley (2007) A. Eyre-Walker and P. D. Keightley, Nat Rev Genet 8, 610 (2007), URL http://dx.doi.org/10.1038/nrg2146.
  • Haldane (1927) J. B. S. Haldane, Proc. Camb. Philos. Soc. 28, 838 (1927).
  • Rouzine et al. (2003) I. M. Rouzine, J. Wakeley, and J. M. Coffin, Proc Natl Acad Sci U S A 100, 587 (2003), URL http://dx.doi.org/10.1073/pnas.242719299.
  • Rouzine et al. (2008) I. Rouzine, E. Brunet, and C. O. Wilke, Theor. Popul. Biol. 72, 24–46 (2008).
  • Park and Krug (2008) S.-C. Park and J. Krug, Journal of Statistical Mechanics: Theory and Experiment 2008, P04014 (2008).
  • (29) E. coli long-term experimental evolution project site, URL http://myxo.css.msu.edu/index.html.
  • Lenski et al. (1991) R. E. Lenski, M. R. Rose, S. C. Simpson, and S. C. Tadler, Am Nat 138, 1315 (1991).
  • Jezequel et al. (Article accepted) N. Jezequel, M. C. Lagomarsino, F. Heslot, and P. Thomen, Genome Biology and Evolution (Article accepted).
  • Elena and Lenski (1997) S. F. Elena and R. E. Lenski, Evolution 51, 1058 (1997).
  • de Visser and Lenski (2002) J. A. G. M. de Visser and R. E. Lenski, BMC Evolutionary Biology 2, 19 (2002).
  • Wielgoss et al. (In press) S. Wielgoss, J. Barrick, O. Tenaillon, M. Wiser, W. Dittmar, S. Cruveiller, B. Chane-Woon-Ming, C. Médigue, and R. Lenski, Proc Natl Acad Sci USA (In press).
  • Elena and Lenski (2003) S. F. Elena and R. E. Lenski, Nat Rev Genet 4, 457 (2003), URL http://dx.doi.org/10.1038/nrg1088.
  • Dykhuizen and Hartl (1983) D. E. Dykhuizen and D. L. Hartl, Microbiol. Mol. Biol. Rev. 47, 150 (1983).
Comments 0
Request Comment
You are adding the first comment!
How to quickly get a good reply:
  • Give credit where it’s due by listing out the positive aspects of a paper before getting into which changes should be made.
  • Be specific in your critique, and provide supporting evidence with appropriate references to substantiate general statements.
  • Your comment should inspire ideas to flow and help the author improves the paper.

The better we are at sharing our knowledge with each other, the faster we move forward.
The feedback must be of minimum 40 characters and the title a minimum of 5 characters
Add comment
Loading ...
This is a comment super asjknd jkasnjk adsnkj
The feedback must be of minumum 40 characters
The feedback must be of minumum 40 characters