Phenotypic evolution implies sequential rise in frequency of new genomic sequences. The speed of the rise depends, in part, on the relative fitness (selection coefficient) of the mutant versus the ancestor. Using a simple population dynamics model, we show that the relative fitness in dynamical environments is not equal to the geometric average of the fitness over individual environments. Instead, it includes a term that explicitly depends on the sequence of the environments. For slowly varying environments, this term depends only on the oriented area enclosed by the trajectory taken by the system in the environment state space. It is closely related to the well-studied geometric phases in classical and quantum physical systems. We discuss possible biological implications of these observations, focusing on evolution of novel metabolic or stress-resistant functions.
Organisms react to long-term changes in environmental conditions by sequential rises in frequency of new genome sequences, mostly corresponding to increasingly more adapted phenotypes. However, often environmental changes are faster than the characteristic time for mutation-selection cycles needed to evolve an optimal phenotype. In such cases, depending on the structure and timescales of the fluctuations, a dynamic environment creates dynamic fitness landscapes , promotes sensing , modularity [3,4], switching , and can change the speed of adaptation [6,7].
The effect of fluctuating selection and/or population size on the population-genetics dynamics has been extensively studied over the years [7,8], starting with the introduction of the concept of adaptive topography by Wright . Dempster was probably the first to relate a population growth rate to the geometric mean of its fitness across the experiences environments . More recently, the evolutionary dynamics of populations with density-dependent dynamics in fluctuating environments has been elucidated in ecologically realistic models [11–13]. These bridge the gap between the classical population dynamics exhibiting very diverse responses to fluctuating environments [14,15] and classical population genetics models. However, a complete understanding of the effect of fluctuations on population and evolutionary dynamics has not been achieved yet.
Some of the relevant parameters describing evolutionary response of a population to a changing environment are the rate at which new genotypes are created (mutation rate), the relative fitness of new phenotypes and the total population size. Each of these parameters depend on the magnitude and speed of the changing environment. We concentrate on the case of environments changing on scales longer than an individual's lifetimes. This is relevant, in particular, for bacterial populations confronted with daily environmental changes (natural or artificial) , for longer living organism affected by seasonal variations or for pathogens experiencing transmission, uncontrolled growth in a new host and then effects of the host immune system. For example, in the now classic long-term Escherichia coli evolution experiment , bacterial cultures are diluted daily, and the environment (i.e. cell growth and death rates) changes during dilution events and between them owing to depletion of resources, cell density growth and cell-to-cell interactions. These experiments are a great model to study clonal competition . Interestingly, the number of accumulated beneficial mutations is relatively small, considering that every single point and many possible double mutations have happened thousands of times in the 25-year history of the experiment. This discrepancy is likely largely accounted for by strong bottlenecks at dilution times, when most new mutations disappear by chance. However, all clones, even beneficial ones, experience additional huge fluctuations in their reproductive rates during the course of the experiment. It remains to be seen if such fluctuations can contribute to the slowing down of the evolutionary adaptation as well.
In this article, we make a step in this direction by studying effects of fluctuating environments (represented by birth and death rates) on the effective selection coefficient. Using analytic and computational tools, we investigate a model of a heterogeneous population (a background strain and a newly emergent mutant) under the assumption that the timescales of the clonal frequency dynamics on the one hand and the environment fluctuation on the other hand are both much larger than the division time, but not necessary well separated from each other. We start by showing that the selection coefficient in infinitely slowly changing environments (we call this the adiabatic limit) is given by a time-average of static selection coefficients corresponding to each environment. This time-average is equivalent to the geometric mean result of Dempster . It is independent of the order or the speed with which different environmental states are visited. However, for environments varying at a slow but finite rate, this time-average is not the whole story. A new contribution emerges. For example, in a cyclically oscillating environment, this new contribution to the selection coefficient now depends on the sequence of environments visited during each cycle, while still remaining independent of the speed of variation. The contribution is non-zero only for non-trivial couplings between the environment and the population dynamics, represented as a multi-dimensional trajectory in the space of birth and death rates. The contribution changes sign when the sequence of the visited environments is reversed. It is largely independent on the speed of the dynamics. Finally, it scales quadratically with the amplitude of the environmental fluctuations. In other words, the contribution is geometric in nature and parallels the discussion of geometric phases in mathematical physics . We believe that this has not been noticed before in the context of population dynamics. In particular, the concepts of geometric mean and geometric phases are totally unrelated.
We will focus on the deterministic approximation to population dynamics. This focus allows us to state beyond doubt that the observed geometric phases are unrelated to the stochastic treatment of the temporally variable selection problem . Geometric effects are well-known for slowly changing deterministic dynamical systems [19,21]. While evolutionary dynamics of a population driven together by forces of mutation, drift and selection cannot be accurately described deterministically, we believe that our model is meaningful even for a stochastic case for a large population, low mutation rate and strong selection. Indeed, the recent observation that stochastic dynamical systems are also subject to geometric corrections suggests that deterministic versus stochastic treatment of population dynamics is not crucial for the phenomenon [22,23].
In what follows, we develop our results in a relatively simple two species population model with bilinear, symmetric competition, which we believe is general enough to capture the main effects of fluctuations for a large class of related models. We first solve the system in the limit of small differences between the birth and the death rates of the competing species. We derive expressions for the selection coefficient in the limit of stationary, very slowly continuously and infrequently discontinuously varying environments. The selection coefficient for an arbitrary timescale of the environment fluctuations can be derived then using a perturbative approach.
Let xi be the number of individuals of genotypes i, i = 1,2, in a large asexual population. We assume that xi ≫ 1, so that demographic (phenotypic) fluctuations and random genetic drift can be neglected. We refer to x1 as an ancestral phenotype, and to x2 as a mutant. The competition between the two is described by a driven two-dimensional Lotka–Volterra (logistic) model [14,15] 2.1Here bi(t) represents the birth rates, and di(t) parametrize the death rates for each of the genotypes. Generally, all parameters are time-dependent. This equation describes the ecological competition between populations of ancestral and mutant phenotypes, which may eventually result in displacement of the former by the latter.
Following classical models of ecological population genetics, we view our model as a particular form of the more general dynamics. Defining the total population size, x(t) = x1(t) + x2(t), we write 2.2Here, g is the generalized growth rate. For this system of equations to represent the dynamics of a realistic self-sustaining population, gi(x) must be negative for large x, and it must have at least one zero. Our approach applies to a very general subset of such growth rate functions provided that the system, equation (2.2), has exactly one fixed point on each of the axes xi = 0 in addition to the trivial unstable extinction point (0,0).
One traditionally takes  2.3where ri's are the intrinsic maximum growth rates of each genotype, if unconstrained by limited resources. The terms rif(x)/f(Ki(t)) represent the reduction of these rates due to competition for resources. This reduction depends only on the total population size x(t) and on Ki, which are stable total populations of the isolated phenotypes i supported by stationary resource-limited environments. K's are referred to as the carrying capacities. Our approach applies for any non-negative, monotonously increasing f(x), as explained earlier. However, for simplicity, we now concentrate on f(x) = x. In this case, the competition is linear and symmetric, and the simple Lotka–Volterra model (2.1) is recovered with bi(t) = ri(t) and di(t) = ri(t)/Ki(t).
We are interested in modelling competition of the ancestral genotype with the mutant one. The two are very close in the genotype space, essentially one mutation away. Because mutation effects are, in general, small , we assume that the differences between g1 and g2 are also small, 2.4This corresponds to small differences in the parameters bi, di, ri and Ki. We assume this from now on. In particular, it is possible that differences between the mutant and the ancestor parameters at any particular time are much smaller than the variations of each of the parameters over time.
In order to determine the conditions under which the mutant, initially present in small numbers relative to the ancestor, invades the population, we explicitly integrate the model, equation (2.1). We write the dynamics of the total population size x = x1 + x2: 2.5To the zeroth order in ε≪1, this does not depend on the individual values x1 and x2: 2.6where we have defined 2.7We also define 2.8the fraction of the mutant in the whole population. This obeys The model then reduces to 2.9where we have used the notation δ(b,d) for small (order ε) time-dependent differences between the corresponding mutant and ancestral rates. To simplify the notation, for any pair of parameters (𝒫1, 𝒫2) describing the ancestor and the mutant, we write 𝒫 = (𝒫1 + 𝒫2)/2, and δ𝒫 = 𝒫2 − 𝒫1. In addition, we always assume |δ𝒫/𝒫| = O(ε) ≪ 1. The second of equations (2.9) is one of the simplest models of evolutionary game dynamics , whose outcome entirely depends on the time average of the function s(t) = [δb(t) − δd(t)x]. In what follows, we describe how the environment-determined fluctuations of birth and death rates determine this time-dependent selection coefficient s(t), and hence the outcome of the competition between the mutant and ancestor populations.
To the zeroth order in ε, the dynamics of the total population size defined by equations (2.1) is now uncoupled from the dynamics of the mutant fraction 2.10Owing to the small variation assumption, equation (2.4), |δb(t) − δd(t)x| ≪ |b(t) − d(t)x| and p(t) changes on timescales much longer than x(t). On these timescales, x(t) converges to a unique (up to the first order in ε) attractor xa(t), independent of the initial conditions, 2.11Then the slower dynamics of p is 2.12where logit p = logp − log(1 − p). The obvious first lesson from this equation is that the clone with the largest average growth rate, for some large 𝒯, will have an advantage.
2.3. Selection coefficient
For coefficients varying periodically with a period T, we write for the logarithmic change of the mutant-to-ancestor ratio, logit p, over time 𝒯 ≫T, 2.13where the last equality defines the selection coefficient, s. It is the sign of s that decides the stability of the fixed points p = 1 and p = 0. For example, for s > 0, p = 0 is unstable, and the mutant phenotype invades the population towards a stable fixed point p = 1.
In a constant environment, and for ε ≪ 1, the selection coefficient s can be rewritten in terms of the ecological parameters defined in equation (2.3): 2.14We have a classical result that selection favours phenotypes with larger carrying capacities (larger Ki) independent of the magnitude of the intrinsic growth rates ri [11,12]. To derive this, we rely on the fact that the total population is given at all times by K, and it is independent of the frequency of the mutants in the population.
In this paper, we are interested in the values of the selection coefficient for temporally varying environments. As a consequence, the selection coefficient is now given by the interaction between several varying quantities. To simplify the discussion, we focus on limiting cases of large time-scale separation between the environment fluctuations and individual lifetimes.
In the regime of infinitely fast environmental fluctuations, for 𝒯 → ∞, we approximate the general driven model, equation (2.2), as 2.15We assume here that the environment variation attains a well-defined, constant average for every state (x1, x2). We denote this by 〈 … 〉T, where the subscript T stands for averaging over a period. We assume that x does not change appreciably over this time. For the specific case of the Lotka–Volterra model, the selection coefficient for fast fluctuations, sf, can be computed using the formula for the constant case, equation (2.14), keeping in mind that one has to use the average values of the relevant coefficients: 2.16
In the opposite limit of an infinitely slow parameter variation, the total population is equal to the carrying capacity at all times, x(t) = K(t). In this case, the quasi-stationary (qst) selection coefficient s is 2.17where the period T is much longer than the individual's lifetime. This allows for a proper average to be attained.
In both limits, the sign of the selection coefficient does not depend on the average carrying capacity [12,13]. Indeed, it is possible to have a slowly varying environment, in which the mutant has, on average, a larger carrying capacity but a lower fitness. In both limits, the selection coefficient becomes independent of the speed of environmental variations, and it is symmetric with respect to time reversal for the driving parameters.
3.1. Continuous, deterministic, oscillatory environments
We now proceed to a more realistic case of an environment fluctuating slowly, but not infinitely slowly, compared with an individual's lifetime. This condition allows us to derive a perturbative approximation for the selection coefficient valid when b(t), r(t) ≫ 1/T are satisfied at every t. Our approximation is based on a simplified solution for the dynamics of the total population size xa(t), equation (2.11). By making a variable change , we write 3.1In the limit of slow environmental changes, the carrying capacity K(y) varies slowly, and the integral in the denominator is dominated by the value of 1/K(z) around z = y. In this regime, 3.2
Using equation (3.2), we now derive an approximation for the total population trajectory xa valid in the qst regime. We denote it as xqa, 3.3This solution represents the correction to the quasi-stationary result xqst(t) = K(t) as a first-order perturbation in the small ratio between the rate of change of the environment and the typical rate of change of the total population. Note that the approximation is consistent with the intuition that the instantaneous total population falls behind the instantaneous carrying capacity.
The selection coefficient can be expressed now as 3.4where 3.5is a geometric contribution to the selection rate. The geometric nature of this term can be better understood if we express the change in the mutant-to-ancestor ratio as 3.6We note that, for any reparametrization of time, λ = λ(t); Δgeom can be written in a very similar form: 3.7This expression shows that Λgeom explicitly depends on the trajectory, [r1,2(λ), K1,2(λ)], but not on how this trajectory is traversed, as also can be seen in figures 1 and 2. Just as any other closed contour integral, this expression can be transformed into a surface integral over any two-dimensional domain bounded by the trajectory [r1(t), r2(t), K1(t), K2(t)] in the parameter space. In particular, using variables 3.8and the Stokes theorem, we can equate Δgeom(𝒯 ) = ∫𝒳d𝒴 with the oriented area bounded by the trajectory for times t ∈ [0,𝒯) in the plane (𝒳, 𝒴).
In other words, Δgeom is a truly geometric term in the spirit of geometric phases in quantum or classical mechanics [19,21]. The geometric nature of the change in the population composition over long times, equation (3.7), is the main result of the paper. It allows us to make important macroscopic predictions about the population dynamics that will hold generally irrespective of the microscopic details of the model. First, the geometric changes in the relative fraction of the mutant depend on the sequence of the environmental states in addition to their identity: same environmental states may have very different effects depending on the order in which the states are visited. At an extreme, a reversal of the order (time-reversal) would change the sign of the geometric contribution, which may make a deleterious mutation advantageous, and vice versa. To our knowledge, such dependence of the effective selection coefficient on the sequence of the environmental states has not been noticed before in simple Lotka–Volterra models. Second, the contribution to Δgeom depends only on the oriented area covered in the parameter space (and thus, in particular, on the number of periodic oscillations), but not on the speed of traversal of the trajectory. Figure 1 illustrates these features: even when the environmental dynamics involves backtracking, the overall contribution per period still does not change. The dependence on the area in the parameter space also suggests that the geometric contribution scales as the square of the fluctuation amplitudes. Finally, to achieve a non-zero area, more than one parameter must be changing, and they must change incoherently. We illustrate some of these features in figure 2.
3.2. Switching among discrete environment states
The approach can be extended to a more common model of piecewise constant environments [2,7]. Consider the case of parameters abruptly changing between m sets indexed by μ = 1, … , m, (r1μ, r2μ, K1μ, K2μ), at possibly random times ta. The state occupied between ta and ta+1 will be denoted by μa. We assume that the interval (ta+1 − ta) is long enough so that the total population x(t) reaches the carrying capacity long before the environment switches again—that is, 1/rμi ≪ (ta+1 − ta). In this case, one can derive the qst contribution as a sum over all of the environment states 3.9At each switch, there is an extra contribution because x(t > ta) reaches the value Kμa with a delay. That is, from equation (2.10), we derive 3.10Integrating equation (3.10) results in a geometric contribution after M environment state changes 3.11The fact that equation (3.11) is independent of the actual time spent in each state and depends only on the sequence of environmental states is the signature of its geometric nature, illustrated in figure 3. Importantly, unlike in the continuous variation case, equation (3.7), Δgeom in equation (3.11) can have a finite value even if parameters change only between two states. Hence, it is unclear whether the contribution can be interpreted as an oriented area enclosed by the trajectory in the parameter space.
3.3. Continuous stochastic environments
Often environments change in a continuous but unpredictable way, such that the typical rate of change is still small. This scenario is modelled with Gaussian fluctuations of the parameters [13,25]. Denoting all parameters with a single symbol γα, α = 1, … , A, we generalize our result, equation (3.7), and represent the geometric contribution for randomly driven equations (2.2) as a line integral [22,23,26] 3.12Here 𝒯 is a long time that allows for averaging, and fα are some model-dependent functions. Because fluctuations are small, we expand fα to the first order in the fluctuations around the average parameters 3.13Now using suitable continuity properties of the parameters' trajectory, we transform the geometric contribution to 3.14The geometric properties of Δgeom are clear from equation (3.14): Δgeom(𝒯) depends only on the length of the parameters' trajectory, is antisymmetric with respect to time reversals and is non-zero only if multiple parameter vary simultaneously and incoherently. Note that equation (3.14) is valid only for parameter variations with small (bounded) speeds. Therefore, if the parameter dynamics, γα(t), are modelled as multi-dimensional Wiener processes, care must be taken to regularize and properly define the stochastic integrals in equations (3.12) and (3.14) .
Equations (3.12) and (3.14) represent a natural extension of the geometric correction to acyclic trajectories . While now the geometric term Δgeom(𝒯) is aperiodic, for parameters dynamics with a stationary distribution of γα and , Δgeom(𝒯) still has a mean linear dependence on 𝒯 for large times: 3.15where Cαβ(t) = 〈γα(0)γβ(t)〉 are time-dependent correlations of the environment. Note that the derivatives dCαβ(t)/dt|t=0 are inversely proportional to the correlation times of the process. Moreover, one can identify the terms in the right-hand side of equation (3.15) as products of the Berry curvature, κβα − καβ, previously introduced in the classical and quantum geometric phases literature [19,23], and, for α ≠ β, the rates of growth of the oriented areas bounded by the process dCαβ(t)/dt|t=0.
3.4. Possible experimental effects
The existence of geometric corrections to fitness in a time-dependent environment requires that changes in the environment are felt by the population on multiple timescales. In the model, equation (2.1), the immediate change in the growth rates and the delayed effect of the population reaching the carrying capacity provide these scales, but other mechanism would work as well. Similar effects will be encountered in almost any situation when a population responds to asynchronous changes in multiple external stresses or nutrient supplies. Therefore, the geometric effects must be considered when modelling emergence or fixation of new metabolic or stress-resistance functions in the presence of environmental changes. We suggest that the relative timing of fluctuations of extracellular nutrient/stressor concentrations will affect the relative fitness advantage of these functions.
Of particular interest is the emergence of antibiotic resistance in bacteria. Mutations conferring antibiotic resistance often decrease the ability of cells to grow in the absence of antibiotics, but provide a growth advantage in their presence . At the same time, delivery of antibiotics is hardly ever uniform, and nutrient supplies also fluctuate. Focusing for simplicity on periodic nutrient and antibiotics concentration changes, we see that the time delay, or the phase lag, between the changing concentrations will join their amplitudes and the period in selecting whether a resistant strain will fix or not. We illustrate this in figure 4: depending on the phase difference between the nutrient and the antibiotic influx, either the resistant or the faster growing bacterium will be selected for. A robust prediction of our theory is that the difference in the logarithmic fractional population changes between an environmental trajectory and its time reversed counterpart will grow almost linearly in time with the number of periods. We emphasize that the effect is different from episodic selection , where only frequencies and magnitudes of antibiotic selection episodes determine fixation of the resistant strain.
Another experimental system where our predictions can be important is the evolution of a metabolic pathway corresponding to a new metabolite, when both the old and the new metabolite concentrations change in time. In such a case, one would need to take into account possible effects of catabolite repression and di-auxic growth in addition to instantaneous effects of metabolite concentrations on the birth/death rates. Nevertheless, we expect that careful modelling of these effects will also uncover the fitness sensitivity to the timing of pathway activation.
An important characteristic of the geometric effect is that it is much harder to be observed in typical serial dilution experiments, especially when the environment changes are only imposed at the dilution points. Such experimental protocols will miss important effects that may be relevant for wild-type conditions.
Fixation dynamics of mutants in a large class of mathematical models is governed by a single effective parameter, the selection rate, obtained as a time average of the instantaneous growth rate difference between the mutant and the ancestral population. In population dynamics with symmetric competition, and in the limit of small differences between the mutant and the ancestor, the total population size is decoupled from changes in the population composition. Instead, the total population enters the fixation dynamics only as a time-dependent parameter. Then the population growth rates and the selection coefficient depend on the interplay between the timescales of the population dynamics and the environmental fluctuations. For infinite separation between the timescales, the selection depends only on values of environmental parameters. More specifically, here the fitness difference can be expressed as a function of growth rates and carrying capacities averaged over all of the environmental states and independent of the period of the fluctuations. Nonetheless, owing to the nonlinear dependence of the growth rates on the environmental parameters, the average fitness difference is not necessarily the same as the fitness difference for the average environment.
This quasi-steady-state approximation breaks down for faster environmental changes. The mutant fraction dynamics is now dependent not only on the period of environmental changes but also on the sequence of successive environmental states. In particular, the first non-adiabatic correction is always anti-symmetric with respect to time reversals, and it is geometric in nature. As long as the fluctuations in the parameters are large, this non-adiabatic correction can be of the same order of magnitude in the birth and death rate variations as the qst contribution to the fitness difference. The geometric nature of this term constrains the effect that environment fluctuations can have on fitness differences. Indeed, like other geometric contributions [19,21], this effect is independent of the instantaneous speed of variation of parameters. In ecological terms, this implies that the geometric contribution to the mutant ratio drift does not depend on how fast the environment changes, but only on the sequence of environmental states. We illustrate this in figures 1 and 3. Further, we note that the mutant fraction drift, Δ, can be seen as a line integral in the parameter space, cf. equation (3.7). This implies that only multidimensional and off-phase parameter variations can give non-zero long-term contributions to the population dynamics.
For the results derived in this work, the assumption of an oscillatory environment is not essential. Our conclusions, and the concept of geometric phase in general, are valid for non-cyclic environment dynamics . Typically such dynamics is represented with a Gaussian and, in general, uncorrelated noise [8,12,13]. While a detailed extension of the present results to random trajectories is beyond the scope of this paper, we have shown here that the geometric contribution to the selection coefficient is present generically if and only if the population dynamics contains multiple correlated parameters driven by a coloured noise, cf. equation (3.15).
In this article, we have focused on deterministic population dynamics with small parameter differences among the competing species, which is equivalent to frequency-independent selection. We expect that a similar geometric phase contribution to the fixation dynamics is present in stochastic Fisher–Wright type models, as well as models that exhibit phenotypic switching or frequency and density-dependent selection effects. Indeed, it has been shown recently that evolutionary processes, including mutations, genetic drift and time-dependent selection satisfy a fluctuation relation  that is similar to the Jarzynski equality . Building upon the deep mathematical connections between fluctuation theorems and geometric phases , we expect that the geometric contributions also shape the evolutionary dynamics in more complex, stochastic models. In particular, the measures of selection for populations of switching phenotypes [2,31] likely contain a geometric term that is awaiting an explicit identification.
Our detailed analysis of multiple specific models allows us to make a general conclusion that is independent of the exact variation of the parameters and the exact details of the model. Namely, for clones with the same mean fitness, the clone that has a higher growth rate when the environment is abundant (increasing carrying capacity) will have a selective advantage over the clone that performs well when the carrying capacity decreases. This is important during acquisition of new metabolic or stress-response functions, as discussed earlier. Further, in the case of the long-term E. coli evolution experiment , we point out that unless mutations manifest themselves in a positive way during the exponential growth phase following a serial dilution, daily variability of the environment would make it harder for mutations to fixate even without stochastic effects associated with the dilution bottlenecks.
We conclude with an observation that species with the fitness advantage in the average environment, with the average fitness advantage over all environments, and with the average fitness advantage for a particular time course of the environment are not necessarily the same species. In particular, a naively deleterious mutation can fixate in a population owing to these temporal effects. We believe this to hold true independent of many of the simplifying assumptions of our toy model.
We thank R. Austin, B. Levin, J. Otwinowski, M. Tchernookov and N. Sinitsyn for important discussions that have shaped this work. We are particularly grateful to B. Levin for his insightful critique of the manuscript and the approach.
- Received October 10, 2011.
- Accepted October 31, 2011.
- This journal is © 2011 The Royal Society