|Home | About | Journals | Submit | Contact Us | Français|
We consider a fixed size population that undergoes an evolutionary adaptation in the weak mutation rate limit, which we model as a biased Langevin process in the genotype space. We show analytically and numerically that, if the fitness landscape has a small highly epistatic (rough) and time-varying component, then the population genotype exhibits a high effective diffusion in the genotype space and is able to escape local fitness minima with a large probability. We argue that our principal finding that even very small time-dependent fluctuations of fitness can substantially speed up evolution is valid for a wide class of models.
Organisms adapt to their environment by sequential fixation of beneficial mutations. This process is often visualized as motion of a population, specified by multi-dimensional genomic variables corresponding to the dominant genotype in the population, in the fitness landscape, where the height of the landscape corresponds to the reproductive fitness of an average individual in the population . Fitness landscapes are believed to be rough with many local maxima, and a population may get stuck in one, so that every plausible mutation is deleterious [2–5]. In such cases adaptation to a global fitness maximum requires the rare fixation of deleterious mutations. Even when there is a path of only neutral or weakly selective mutations to the global optimum, it may be difficult to find it, and navigating such paths can be slow due to the low fixation probability (≈ 1/N for a neutral mutation in a population of N individuals ). In view of this, one of the central questions of evolutionary biology is how the current diversity of the living world has emerged in a short few billion years since the life on Earth has started, especially since thousands of generations of laboratory evolutionary experiments lead to only a few dozen fixated mutations .
It has been recognized that temporal fluctuations in the fitness landscape can drive the population out of a local fitness maximum, thereby accelerating evolutionary processes. For example, the maximum of fitness at one time may be on a fitness slope at another time, allowing the population to leave the area. The effects of fluctuating selective pressure on mutation accumulation dynamics have been studied extensively [8, 9], starting with the introduction of the concept of adaptive topography by Wright . More recently the evolutionary dynamics of density regulated populations in fluctuating environments have been elucidated in more ecologically realistic models [10–12], bridging the gap between the classical population dynamics [13, 14] and population genetics models.
In a recent pioneering numerical evolution experiment , these ideas were further developed to show that certain types of fluctuating environmental pressures may speed up evolution many times. However, it remains unknown to what extent these results generalize. Is the speedup a general property? How does it depend on the spatiotemporal structure of the fluctuating environment? Can a population escape any local maximum? How does the motion in the genotype space depend on time? Do fitness fluctuations have to be dramatic, as in Ref. , or can small fluctuations still speed the evolution up?
In this article, we answer some of these questions in the context of a model of evolutionary dynamics that is simple enough to allow a thorough analytical and numerical treatment, but is at the same time general enough so that at least some of our predictions hold for a wide class of evolutionary models. We consider the limit of a weak mutation rate, when the time scales are well separated. The time between successive mutations is longer that the typical fixation time, and the characteristic time scale of the fitness landscape changes is the longest. Such a situation is relevant for microbial populations under seasonal environmental changes, or for host-pathogen interactions, where environmental changes may correspond to phases of transmission, unhindered growth in a new host, and activation of the host immune response. Further, we assume that the total population size is independent of the genotype (though not necessarily fixed), so that the evolutionary dynamics depends only on the relative fitness differences between the genotypes. We consider adaptation in a highly epistatic genotypic space, such that the evolution takes place on a one-dimensional pathway with large, local fitness differences. Under these assumptions, we show that the evolutionary search can be sped up substantially when only a small component of the fitness landscape undergoes temporal variations.
Our model of a fluctuating environment is based on an overdamped particle in a potential. The position x is some generalized coordinate that describes the dominant genotype in the population, and hence a change in x is a fixation event. This genotype changes according to an equivalent Langevin dynamics given by
where U is the potential, γ is the “friction”-like scale factor, η is a white Gaussian noise with variance 2D, where D is the intrinsic diffusivity, presumably related in the population biology context to the population size . Notice that in the usual physics language, the process will minimize the potential so that U is the negative fitness. Motion to minimize U represents fixation of beneficial mutations, while Langevin noise allows low probability fixation of neutral and deleterious mutations. The first phenomenon is called natural selection/drift in evolutionary/physics languages, and the second is unfortunately referred to as drift/diffusion, respectively. To avoid confusion, in the remainder of the article we use the physics terminology.
We write the potential as
and we focus on the following range of parameters:
This models the emergence of novel functions in a population. Namely, the fitness is largely independent of time, as described by U0. However, small temporal changes in fitness are allowed. For example, acquiring a new enzyme is generically advantageous if its substrate is present, but deleterious if it is absent due to generic costs associated with protein overproduction [16, 17]. We model this by adding a small fitness component Φ(x) that fluctuates as S(t), representing, for example, changes in the availability of the metabolite due to seasonal or geological variations. Finally, we choose to separate the global, almost non-epistatic, fitness from the local, possibly highly-epistatic (but small) effects by making the gradient of U0 smaller than that of Φ, even though the scale of Φ itself is smaller than that of U0. For example, this agrees with the observation that the ability of proteins to bind to DNA or to metabolic substrates is highly sensitive to the details of the protein sequence [2–5].
With the conditions above, we can redefine U0, Φ, and S without much loss of generality, so that St = 0. We then consider the simplest form of Φ(x) and S(t) that satisfies these conditions, and we will discuss how our results generalize to some other forms of the functions in Sect. 5. Namely, we choose Φ to be a zero-mean periodic saw-tooth potential, and S to be a zero-mean periodic telegraph signal. These considerations allow us to write near a particular point x in the genotype space
where v is the intrinsic drift (in physics terms) or bias, defined as positive for the drift to the right, see Fig. 1. We always assume that h/L > v, so that the fluctuating component of the potential can actually create local maxima and minima on top of the global landscape U0(x). In what follows, we denote by T the time between subsequent potential flips (the half-period of the fluctuations), and L is half of the spatial period.
This model is similar to various stochastic ratchets considered in the literature [18–21]. Thus, the question of whether the fitness fluctuations can speed up the evolutionary search is a question similar to whether a rectified or a high-variance motion can appear due to ratcheting. We know from prior analysis  that any unbiased spatially variable but temporally constant potential cannot give rise to rectified motion, and it will always slow down diffusion. Hence temporal fluctuations are an essential component of the model.
Using the choices above, we can rewrite the equation of motion, (1)as
where η is a Gaussian white noise of unit variance. In (8), the dynamics explicitly depends on five different parameters L, T, v, h, and D. Nevertheless, by rescaling the time, the space, and the potential as x/L → x, t/T → t, , , we can reduce the number of parameters to only three: the ratio of the typical diffusion time over half the spatial period to half of the temporal period, , the height of the fluctuating barriers in diffusivity (temperature) units, , and the ratio between the slope of the average, large scale potential to the slope of the fluctuating perturbation, v. In physical terms, if ω is large, the particle has time to explore the entire valley of ϕ before the potential flips. Further, β measures the difficulty of crossing the peaks by diffusion. Finally, the condition that the perturbation induces local optima is v < 1.
Using the rescaled variables, the dynamics becomes
We will use these rescaled variables in the rest of the article, unless noted otherwise. From this equation, it is easy to recover the dynamics in the original, non-scaled units by simple multiplications. In what follows we present simulation results obtained using first order Euler integration scheme of the dynamics defined in rescaled variables, (9).
Numerical simulations suggest that the behavior of x(t) at large times is diffusive, and anomalous scaling is not seen [23–26]. This is consistent with general analytical results obtained by multiscale techniques . We can characterize the genotype coordinate motion by effective drift and diffusion constants, which depend on the spatial and the temporal periods of the fluctuations and the barrier height. To quantify the enhancement or the suppression of the motion at long spatial and temporal scales compared to the intrinsic diffusivity and drift, we define
where the time-dependent mean and variance of the trajectories x(t) are obtained numerically. As seen in the Fig. 2, both the effective drift and the effective diffusion can be enhanced with respect to the intrinsic values, this enhancement having a maximum for fluctuation periods comparable to the average time to travel between two extrema.
When β 1, the sawtooth peaks are very small, and the diffusion has no trouble crossing them. When β is large, the behavior is more interesting. For ω → 0, the particle has ample time to fall into a minimum of ϕ before s(t) flips. When the potential flips, the particle can now go either left or right, with unequal but comparable probabilities, which creates a biased random walk behavior with the effective diffusion coefficient ≈ L2/(2T), or, equivalently, rD ≈ ω. The fluctuating potential allows the particle to diffuse against the drift, so that the speed of the evolutionary search is strongly enhanced when the environment oscillates.
For low ω, rD ≈ ω is also small, so the diffusion is suppressed compared to the internal value. However, for h 1 and without the changing sign of s(t), the particle would get stuck at a minimum of ϕ almost immediately, and the overall diffusion would be essentially zero. It is hard to exit a deep potential well only with the help of diffusion. Whether rD is greater or less than 1, the fact that oscillations make it non-zero in the long term is our most important finding, suggesting that temporal fluctuations can make fixation of rare-to-fixate mutations a much more common process. Figure 2 demonstrates these findings for different values of ω and for β = 10 and v = 1/10.
When the flipping is fast compared to the diffusion time (large ω in the figure), the fluctuating potential averages out to zero. The only motion is due to the internal diffusion, and D and rD go to one.
In the limit β 1, the fluctuating peaks are very high, the particle almost never crosses them due to noise, and analytical progress can be made. First consider a particle that starts close to a local minimum x0 of the sawtooth. After some time ~ D/v2, which goes to zero as β → ∞, the particle equilibrates near x0 with a probability density
Then the ratio between the probability, k>, that a particle is located to the right of x0 and will move to the right after the potential flips to the probability, k<, that it is to the left of x0 and will move to the left is
When s(t) changes sign, for a small D the particle then glides down with a constant velocity in its chosen direction, reaching the next minimum to the right (>) or to the left (<) in time
If τ < 1, then the particle has the time to reach the minimum on either the left or the right hand side (assuming, as always, that the sawtooth actually forms the local minima, i.e., 0 < v < 1). Then when the potential flips the next time, the process repeats. This results in a discrete random walk between the extrema of the sawtooth, and (in unscaled variables)
In dimensionless units,
Using (13) results in:
Notably, when D → 0, but rD/β = Deff/h is finite.
In Fig. 3, we compare the analytical results to the numerically estimated rD and rv (here rD is normalized by β/2). As β → ∞, the agreement is clearly seen for small ω, that is, for an infrequently flipping potential.
When the potential changes faster, and τ> > 1 > τ<, the particle fails to make it to the minimum to the right and, after a subsequent flip, always comes back to where it started from. However, it always reaches the left minimum before the flip. When τ> > τ< > 1, it does not reach the left minimum either, but goes the distance to the left, then reverses and travels to the right, reverses again and repeats until it reaches the left minimum. After one flip, it moves , after three flips it moves , and so on. Eventually, when β/(2ω)[(2n + 1)v + 1] > 1, the particle reaches the next minimum. Thus every time 2ω/(βv) crosses an odd integer, more periods are needed to travel between the nearby extrema, and the diffusive behavior changes, resulting in the discontinuities in Fig. 3. These dynamics can be described by a master equation
where Pi(t) stands for being at an extremum i at the end of the flip t, and exactly 2n + 1 flips are needed to travel between the i + 1’th and the i’th extrema. Solving this for the drift and the diffusion (see Appendix)gives
The analytical results and the numerical simulations verifying them are shown in Fig. 3 with ω normalized by β/2.
We have shown that the typical diffusive/drift behavior of the system is enhanced by the fluctuating component of the potential. However, what kind of an effect does this enhancement have on the probability of rare, atypical events, such as escape from a suboptimal fitness maximum? To model a barrier in fitness, we place the overdamped particle in a flipping sawtooth potential and a constant drift v, and we observe the mean first passage time for the particle to reach an unscaled distance , against the drift, with a reflecting boundary at x = 0, see Fig. 1.
If during one half-period the particle is able to travel between the two extrema (independent of the direction and without the help of the noise) then τ < 1, and the probability that the particle travels over multiple periods against the (effective) drift is finite even at a very low noise. As seen in Fig. 4, the escape time over the average barrier in U0 depends on the length of the barrier . In order to understand this dependence, we consider our model in the limit of zero fluctuating potential, which allows us to use an analytical expression for the escape time 
where Pe is the Péclet number
In Fig. 4, we show the simulation results of the exit time for slow fluctuations. We compare the results to what we would have expected under a continuous diffusion with the drift and the diffusion constant given by the values of the effective parameters as obtained in the previous section. We observe that the escape times are well approximated with the “effective” continuous diffusion model. The conclusion is valid even for a very small intrinsic noise D → 0, which implies that the average time to escape over the global barrier is fluctuation activated (as opposed to noise activated), and is much faster.
We emphasize again that the fact that even rare escape events in the fluctuating potential are modeled well with the drift/diffusion approximation should not necessarily be expected from our earlier observation that the typical behavior in this system is diffusive.
The qualitative behavior of particle trajectories changes if the fluctuating potential is fast enough such that the particle can only move less than one period in one direction in the absence of the noise. Even though, on average, the variance of the particle position grows linearly, and one can define a proper effective diffusion coefficient, the particle never crosses a barrier against drift in the absence of the intrinsic noise, β → ∞. Hence the probability of rare excursions against the effective drift cannot be described using the same effective drift/diffusion model.
Figure 5 reports results of numerical simulations at different values of the additive noise. The escape time as a function of still can be fitted well with a drift/diffusion model, (24), with a noise dependent effective Péclet number Pe(β). We conclude that, for small noises, the effective Péclet number is approximatively inversely proportional to the noise strength (the plot seems to reach a constant for β → ∞). This implies that the escape times will become infinite at zero noise. The dependence is consistent with a model without any fluctuating potential, but with some effective parameters. The parameters are such that, with the fluctuations, the escapes are significantly faster. Indeed as shown in Fig. 5, the effective Péclet number, defined as , is always smaller than the equivalent quantity in (24), such that
Using a one-dimensional model of diffusion in the presence of a constant force perturbed by small periodic fluctuations, we have shown that time dependent potentials can significantly speed up the large scale drift, diffusion, and barrier escape times. This conclusion is valid even for a very small intrinsic diffusivity, the long time statistics of the particle trajectories being mainly determined by the properties of the potential fluctuations, but not of the Langevin noise.
Our model is a caricature of evolutionary dynamics in the limit of low mutation rates and constant population sizes. Even in this limit, there are several simplifying assumptions in our model that can be relaxed.
First, the periodicity of the potential time dependence is not crucial. Based on the similarity with Brownian motor models [18, 19], we expect that our conclusions will still be valid for a nonperiodic s(t): the nonperiodic flipping will mix together and average behaviors from the different regions in Fig. 3. Further, Dubkov et al.  studied a randomly flip-ping sawtooth with no drift. Their potential flipped with dichotomous Markovian noise with rate v. They found an analytic expression for the diffusion enhancement rD which grows slower than in our model for small ω, and the peak of rD is lower and at larger ω. This is consistent with the averaging over different regions in Fig. 3. Similarly, if the potential is not spatially periodic, this will introduce a quenched noise and is likely to result in emergence of regions in x that are very hard to cross, similar to .
Our model is based on a piecewise linear periodic potential with discontinuous first derivatives (sharp minima and maxima). We expect that our conclusions stay valid for smoother perturbations as long as, upon a flip, a particle can leave the vicinity of a maximum in a time much faster than a typical travel time between the extrema.
The piecewise-constant perturbation in our model is symmetric with respect to reflection, thus no mean rectified motion can be created [18, 19, 22]. Simply put, the system is not a ratchet, and crossing of large scale barriers instead is achieved by creating a large effective diffusion. Asymmetry of the perturbation may create additional rectification, but the results in Fig. 3 suggest that the diffusive effects will not change qualitatively. This is crucial for evolutionary modeling since crossing barriers by ratcheting requires tuning the ratchet in a specific direction in the genotype space, while diffusion will explore the space isotropically.
The genomic space is expected to be high dimensional. For simplicity we have considered a one-dimensional model of mutation accumulation, which could still be relevant for highly epistatic genomic landscapes with a large proportion of highly deleterious or non-viable genotypes. Moreover, we expect that systematic application of multiscale homogenization methods , would allow one to extend the results to these more general formulations.
If the assumption of well separated time scales and constant population size are not satisfied, the Langevin microscopic dynamics used here is not valid any more. Then different fitnesses can give rise to different population sizes, and hence to varying fixation rates and to a space dependent D in our language, which would require more complex models with position and time dependent noise strengths. Such models would modify the probability of individual trajectories [26, 28], and more work is needed in order to identify the regimes in which such diffusion models and their predictions apply to population genetics dynamics.
We are grateful to S. Boettcher, D. Cutler, F. Family, H.G.E. Hentschel, B. Levin, and A. Singh for stimulating discussions, and to the anonymous referee for valuable comments. We also thank J. Lebowitz for organization of the 102nd Statistical Mechanics Conference in December 2009, which initiated our interest in the subject.
In this Appendix we derive the effective drift and diffusion coefficients for the model of discrete time jumps among discrete sites described by the master equation (21). In this model, a particle at site i at time t can either start moving to the right (with probability b < 1), in which case it returns to the site i after two potential flips (at time t + 2), or it can start moving to the left (with probability 1 − b), in which case it reaches the site i − 1 after 2n + 1 flips. The corresponding discrete time master equation is:
Multiplying by i and summing over it, we get
Assuming that the average can be written as
Now multiplying (27) by i2 and summing over i, we get
This allows to write for the variance of i at moment t + 1
Open Access This article is distributed under the terms of the Creative Commons Attribution Noncommercial License which permits any noncommercial use, distribution, and reproduction in any medium, provided the original author(s) and source are credited.
Jakub Otwinowski, Department of Physics, Emory University, Atlanta, GA 30322, USA ; Email: jotwinowski/at/physics.emory.edu.
Sorin Tanase-Nicola, Department of Physics, Emory University, Atlanta, GA 30322, USA ; Email: sorintan/at/physics.emory.edu.
Ilya Nemenman, Departments of Physics and Biology and Computational and Life Sciences Initiative, Emory University, Atlanta, GA 30322, USA.