Fixation in finite populations evolving in fluctuating environments.

The environment in which a population evolves can have a crucial impact on selection. We study evolutionary dynamics in finite populations of fixed size in a changing environment. The population dynamics are driven by birth and death events. The rates of these events may vary in time depending on the state of the environment, which follows an independent Markov process. We develop a general theory for the fixation probability of a mutant in a population of wild-types, and for mean unconditional and conditional fixation times. We apply our theory to evolutionary games for which the payoff structure varies in time. The mutant can exploit the environmental noise; a dynamic environment that switches between two states can lead to a probability of fixation that is higher than in any of the individual environmental states. We provide an intuitive interpretation of this surprising effect. We also investigate stationary distributions when mutations are present in the dynamics. In this regime, we find two approximations of the stationary measure. One works well for rapid switching, the other for slowly fluctuating environments.

Introduction

Evolutionary dynamics describes the change of populations over time subject to spontaneous mutation, selection and other random events [1,2]. Different phenotypes in the population can emerge spontaneously by mutation, i.e. through errors during reproduction of wild-types. In many cases, wild-type and mutant individuals are characterized by heritable differences in behavioural traits or strategies [2]. Selection acts on different phenotypes and changes the population composition. Changes in the state of the environment can alter these selective pressures over time.

Time-varying environments are relevant in the evolution of bacterial populations subject to environment modulations by a host [3,4], or varying antibiotic stress. An illustrative example is the evolution of normal (wild-type) cells and resistant ‘persister’ cells (mutant). This was examined by Kussell et al. [5], where periods of antibiosis were turned on and off. During times of antibiotic stress, the growth rate of normal cells was reduced, but the resistant cells sustain population levels. In addition, Acar et al. [6] provided further experimental evidence supporting the deterministic model used in [5]. More complicated studies of dynamics in switching environments rely on cells ‘sensing’ the environment [7] and on the history or information of the environment during a cell's life [8,9]. These examples illustrate that the assumption of an interaction structure independent of time is not always realistic. At the same time, it is largely an open question how complex interactions between phenotypes together with spontaneous changes in the environment influence the evolutionary dynamics.

The interactions of phenotypes in a population can be formalized in an evolutionary game [10,11]. Such games can be used to describe conflict over food or territory, cheating in resource allocation, as well as interactions between variants of a gene [12-16]. In an evolutionary normal-form game, each individual can be associated with one out of a finite set of strategies. A payoff matrix quantifies the reward received by a given individual when it interacts with another individual [11].

The dynamics of populations interacting in such a game are often described by deterministic replicator equations or similar differential equations [10,17,18]. While deterministic dynamics are useful to understand the action of selection per se, more interesting phenomena arise when stochastic effects are taken into account. A stochastic approach is appropriate—often even strictly required—to understand the impact of fluctuations in finite populations [19,20]. Deterministic approaches fail to capture effects such as fixation and extinction, or the convergence to a stationary distribution in systems with mutation [21-25].

There is an increasing body of literature on stochastic evolutionary games. For example, analytical results for the probabilities of a single mutant to reach fixation have been obtained [26-29]. However, most of this existing work focuses on games played in a fixed environment; the underlying payoff matrix itself remains unchanged in time. This assumption may not be appropriate in cases where external factors influence the environment.

External fluctuations in evolutionary games have been previously introduced by adding extrinsic noise to continuous model parameters [30], or by letting strategy space itself vary in time [31]. Environmental variability has also been the subject of investigation in predator–prey models [32,33].

In this article, we explore different theoretical approaches that allow calculations of fixation probabilities and mean fixation times of a rare mutation, under fluctuating environmental conditions. We use a generic birth–death framework, as described in §2, and our results thus apply to a wide class of population dynamics. In §3, the theory is developed for an environment that can transit between an arbitrary number of discrete states, and we expand on the two-environment scenarios in §4. To illustrate our theoretical results, we study the fixation properties in an evolutionary game that stochastically switches between a coexistence game and a coordination game in §5. We determine environmental conditions under which the success of a rare invading mutant is maximal. This is seen to occur at a non-trivial combination of switching rates. For the case in which mutations occur during the dynamics, as described in §6, we explore how the stationary distribution of the population changes in fluctuating environments. We derive approximations for the stationary distribution, valid for a large range of switching rates. We summarize our findings in §7 and put them into context.

Mathematical model

We seek to model evolutionary dynamics in finite populations of two species that are subject to environmental changes. The changes in the environment are such that at any given point in time, the system can be in one of a finite set of environmental states. The state that the environment is in determines the details of the birth and death dynamics. We focus on two cases: first, in the absence of mutations in the dynamics, we derive laws to predict the probability and mean time of the fixation of a mutant. Fixation describes the event in which mutants take over the population as opposed to going extinct. Fixation and extinction are the only two outcomes of dynamics of a rare mutant in a finite population [34]. In figure 1 we show a basic sketch of this scenario in which the two monomorphic states of the population are absorbing. Second, we study the case when mutations occur in the dynamics. There are then no absorbing states. Instead, the dynamics converges to a stationary distribution.

Figure 1.

A population undergoes a one-step birth–death process, such that given the population is in state i, in one time-step it may transit to i − 1 or i + 1, or remain at i. The states i = 0 and i = N are absorbing in both environments (no arrows out of these states). The transition probabilities (birth/death rates) are dependent on the state of the environment, indicated by solid versus dashed arrows in environments σ and σ′, respectively. The environment switches from state σ to σ′ with probability μ in any one time-step. The quantity ϕ represents the probability of fixation, as discussed in §3. (Online version in colour.)

Birth–death dynamics

We consider populations consisting of a fixed number of N individuals. Each individual can be of one of two types, A or B, which we refer to as ‘mutant’ and ‘wild-type’, respectively. The population is well mixed; every individual can interact with any other individual. The state of the population is fully characterized by the number, i, of individuals of type A. The remaining N − i individuals are of type B. We furthermore assume that at any one time the environment can be in one of Ω discrete states, labelled , where Λ is the space of states of the environment (|Λ| = Ω). Hence the state of the entire system at any time is given by the pair (i,σ).

The discrete-time birth–death dynamics of the population for a given environment, σ, is then specified by the transition probabilities and of a one-step process. Specifically, if the system is in state (i,σ) the population transitions to state i + 1 in the next time step with probability . Similarly, the state of the population in the next time step is i − 1 with a probability . These transitions are shown as black arrows in figure 1. With probability the population remains in state i. We always assume that and for all (i,σ).

With the exception of §6, we always assume that the states i = 0 (all-B) and i = N (all-A) are absorbing ( and for all ). In the absence of further mutation events, a type, once absent, can never be re-introduced. If mutations are present in the dynamics, then the states i = 0 and i = N are no longer absorbing and the system converges to a unique, non-trivial stationary state. We consider this case in §6.

Fluctuating environment

In our approach, the environment evolves from one state to another independently of the state of the population. This simplification still captures a wide array of natural scenarios. In the discrete-time set-up, we take the dynamics of the environment as a simple Markov chain, described by the transition matrix of size Ω × Ω. The entry μ′ represents the probability that the environment changes to state σ′ in the next time-step, if it is currently in state σ, as shown in figure 1. The matrix is a stochastic matrix, for all .

A switch of the environment effectively modifies the birth–death dynamics in the population. We do not specify the exact type of interaction at this point, but keep the rates general.

Fixation probability and time for a general birth–death process in a fluctuating environment

Fixation probability

Let us first consider a discrete-time evolutionary process. If the system is in state (i,σ) at a given time, it may transition to 3Ω possible states in any one time-step. These are given by (i,σ′), (i + 1,σ′) and (i − 1,σ′), where can be any of the Ω states of the environment. If we write R( for the probability of a transition from (i,σ) to (j,σ′), we have No transitions from (i,σ) to (j,σ′) can occur when |i − j| > 1. In this set-up, the birth–death probabilities are determined by the state of the environment at the beginning of the discrete-time-step.

The fixation probability, ϕi,, is the probability that the system ends up in the absorbing state with N individuals of type A, conditioned on initial state (i,σ). The probability of fixation of a single mutant, ϕ1,, is of particular interest [21]. It is briefly illustrated in figure 1. In our scenario with switching between Ω environmental states, following the lines of the earlier studies [2,28,35], the following balance equation for the fixation probabilities can be found:This is to be solved along with the boundary conditions ϕ0, = 0 and ϕ, = 1 for all .

To obtain a formal solution, we introduce , or in matrix form . The vectors and each have Ω components. The boundary conditions ϕ0, = 0 and ϕ, = 1 translate into ψ0, = 0 and ψ, = 1 for all . With this notation, we haveUsing , we obtain

where . We stress that the calculation of fixation probabilities and mean fixation times using this formalism requires the matrix to be invertible. We comment on this further in the context of a specific example in §4.

To keep the notation compact, we define the variable υ, = ψ, − ψ 1,. Using ψ0, = 0, we have . With this notation, we can write equation (3.4) in the following form:where is the Ω × Ω identity matrix. This relation expresses the vector + 1 in terms of the vectors 1, 2, … ,. We can therefore express all vectors (i = 2, … ,N) in terms of 1. The constraint then determines 1 self-consistently. We note that the resulting set of equations is linear in the set {υ1,σ}. Hence a solution can be obtained in a closed form, in principle. In practice, one inverts the linear system using one of the standard algebraic manipulation packages. Once 1 has been found, the other components , with i = 2, … ,N, can be computed via equation (3.5). One then uses to find the fixation probabilities starting with i individuals of type A in environment σ, {ϕ,}.

We note here that algebraically inverting the linear system (3.5) when N is large is difficult due to the very large number of terms in the corresponding expressions. Thus, at present, this theory is limited computationally to relatively small N.

In the case of a single environment, Ω = 1, the matrix is simply the 1 × 1 identity matrix, and equation (3.5) simplifies to the well-known result for discrete-time birth–death processes [2,28,35].

Unconditional fixation time

We write t, for the expected number of time-steps taken to reach any one of the two absorbing states, given that the system is started in state (i,σ). These fulfil the boundary conditions t0, = t, = 0. With these definitions we find the following relation:Introducing the variable , we haveand with the notation ν, = ξ, − ξ 1,, we arrive atThis relation allows one to express all vectors (i = 2, … ,N) in terms of ν1. The constraint then determines 1, and the mean unconditional fixation times are computed using .

Conditional fixation time

We write for the mean fixation time conditioned on absorption in the all-A state, given that the system is initially in state (i,σ). To find this conditional fixation time, we proceed along similar lines as before. Introducing the variable , which has boundary conditions θ0, = θ, = 0, the following balance equation can be found:We note that equations (3.6) and (3.9) appear to be very similar, but the difference is more than just a global pre-factor ϕ,; each term in the expression has different indices i and σ.

Introducing the variable , we haveand introducing η, = ζ, − ζ 1,, we arrive atThe set {θ,} can then be found using an approach similar to the one described above. Results for the mean conditional fixation time can then be obtained using .

Continuous-time model

In any of the elementary time steps of the above discrete-time model, both the state of the population (i) and the state of the environment (σ) can change. We next consider a continuous-time set-up. There are two types of discrete events that may occur at any time: (i) the state of the environment may change or (ii) a birth–death event may occur. The rate (per unit time) with which a transition from state σ to state σ′ occurs is denoted by m→. The occurrence of these events is independent of the state of the population. The rate with which a birth–death event of the type i → i + 1 occurs is , if the environment is in state σ. The rate with which i → i − 1 occurs is . We write P,(t) for the probability to find the system in state (i,σ) at time t, and find the master equation

Fixation probability

We write Q,(t) for the probability to find the system in state (j,σ′) a period of time t after it has been started in state (i,σ). The corresponding backward master equation [36,37] readsWe define as the probability that the system has reached fixation in the all-A state a period of time t after the dynamics has been started in state (i,σ). This includes fixation before time t. By setting j = N and summing over σ′ in equation (3.13), we obtainThe fixation probabilities are found as ϕ, = lim→φ,(t), and they can be obtained by setting the time derivative in equation (3.14) to zero. Introducing y, = ϕ, − ϕ 1,, one findswith and where we have used ϕ0, = 0 to write . For the case of a single environment, the second term on the right-hand side vanishes and one recovers again the well-known results in single environments [2,28,38]. Equation (3.15) has the same general structure as equation (3.5). Keeping in mind that ϕ, = 1, the fixation probabilities can hence be found by applying the approach outlined in §3.1.

Fixation times

Calculating the mean fixation time using a diffusion approximation [34,36,37,39] is not appropriate for our model. The environmental switching process has no continuum limit. Instead we work with the backward master equation (3.13) and adapt the calculation outlined by Antal & Scheuring [21].

We introduce ϑ,(t) = ∑[Q0,(t) + Q,(t)], the probability that the system has reached fixation in either of the two absorbing states a period of time t after being started in (i,σ). Again this includes fixation before t. We then have ρ,(t) = ∂ϑ,(t) for the probability density to reach fixation exactly at time t. From the backward master equation (3.13), we findThe mean unconditional fixation time is then found via , from which we findA similar iterative equation can be found for the mean fixation time conditioned on absorption in the all-A state. The only difference is the integral of is given by the fixation probability ϕ,, and that . The mean conditional fixation times, , therefore fulfil the relationStructurally, equations (3.17) and (3.18) are of the same form as the corresponding equations for the discrete-time model, and so they can be solved using an analogous procedure. In continuous time, however, the solution procedure no longer relies on the invertibility of the matrix of switching rates between the states of the environment. This is because we have split up the birth–death dynamics and the changes of the environment into separate events that occur successively.

Switching between two environments

We now focus on the case of environments which can be in one of two possible states, i.e. Ω = 2. We label the two states as σ = ±1 (Λ = {+1,−1}). We focus on the discrete-time scenario. The matrix can then be written aswhere the quantity is the probability that the state of the environment switches to −σ in a given time-step if it is in state σ at the beginning of this step. We recall that our theoretical results require the inversion of . Excluding the case when vanishes, this inversion can be carried out straightforwardlyFor the case Δ = 0, we have verified that there is no anomalous behaviour of simulation results.

Fixation probability and times

The general result of equation (3.5) now reduces to the recursionThe fixation probability is obtained from the set {ν,} viaSimilarly, equations (3.8) and (3.11) reduce toandrespectively. The mean unconditional and conditional fixation times are then found, respectively, asand

Effective description for fast switching

The environmental change is fast if the environmental states are short-lived, i.e. much shorter than the mean fixation time in either environment. Then we expect the population dynamics to be controlled by a set of effective transition probabilities, i.e. weighted averages of the original transition probabilities in the different environmental states. The weights are given by the fraction of time spent in each environmental state. As the dynamics of σ follows a simple telegraph process [36], the asymptotic fraction of time spent in the state σ is p/(p+p) for . Using this, the effective transition probabilities are given byWe note that p is the probability that in a given time-step the environment switches from state σ to −σ. Hence the time spent in state σ decreases with increasing p if p is held fixed.

We anticipate that expression (4.9) can formally be derived by introducing a relative scaling parameter between the switching probabilities and the birth–death probabilities, and then by taking a suitable limit in which the timescales of both processes are widely separated. We do not explore this route further here.

In this approximation, the dynamics of the population are mapped to a simple birth–death process on the set with absorbing states i = 0 and i = N. For such processes, explicit expressions for the fixation probabilities and mean fixation times are known [2,28,35]. In the fast-switching limit, we propose the following approximation for the fixation probability:We write here . The corresponding approximations for the mean unconditional and conditional fixation times of a single mutant, respectively, areandThese expressions exactly describe the fixation properties of a birth–death system with the effective transition probabilities; the nature of our approximation is to assume that the birth–death process in quickly changing environments can be described by the effective transition probabilities in equation (4.9).

Finally, we note that this theory is independent of the invertibility of the switching matrix .

Fixation in fluctuating two-player two-strategy games

Evolutionary games

As a direct application of the general theory we have developed, we now consider evolutionary game dynamics in well-mixed, finite populations. Any of the N individuals can be of one of two types, A or B. We limit the discussion to two-player games, but the extension to multi-player games (e.g. [25,40,41]) is straightforward.

At any point in time the environment is in one of two discrete states . This state fluctuates in time as specified above. The interaction between individuals is characterized by the payoff matrixThe subscript σ indicates the dependence on the environment. The matrix focuses on the column player: a type-A individual encountering another of its kind receives a, and it receives b when interacting with a type-B individual. In turn, an individual of type B interacting with an individual of type A obtains c, and d is the payoff for each individual if they are both of type B.

If the environment is in state σ, and if there are i individuals of type A in the population and N − i individuals of type B, the expected payoffs for each type of player are The reproductive fitness of any individual is a function of the individual's payoff in the evolutionary game. We use an exponential mapping [42,43] This is a common natural choice in which fitness is never negative and monotonically increasing with payoff. Any other functional form of the payoff-to-fitness mapping with these two properties would be equally appropriate. The constant parameter β > 0 is the so-called intensity of selection. Based on this definition of fitness, we model the evolutionary dynamics by the update rules of the Moran process [34,44], which has been widely used in evolutionary game theory [2,28,45]. The Moran process represents a simple birth–death process in which the population size remains constant, and by construction it has absorbing states at i = 0 and i = N. In a discrete-time setting, the frequency-dependent Moran process is specified by the transition probabilities [46]where is the average fitness in the population. We note that the framework of the previous section can be applied to microscopic evolutionary dynamics other than the Moran process. This includes, for example, pairwise comparison processes [47,48] or cases with constant selection in any one environment.

Switching between coexistence and coordination games

Rare mutations can introduce a previously absent strategy into the population. Typically, there is only one individual of this novel type initially. We say that B is the resident type, and that A is the invading mutant type. All results in this section are based on the initial condition i = 1. We chose a = d = 1 for the payoff matrix. The type of game is then determined by the off-diagonal terms. We chose b = 1 + σb and c = 1 + σc, where b and c are real-valued parameters. Thus, we have the payoff matrixOur parametrization does not span the entire space of all 2 × 2 games, but it covers some of the most common types (see below).

There exist three general types of two-player two-strategy evolutionary games. First, for the coexistence game b > 1 and c > 1, selection drives the population away from the absorbing boundaries. Second, for b < 1 and c < 1, the population dynamics exhibits bi-stability. This is also known as a coordination game; selection drives the population towards the monomorphic states. In both cases, there exists an internal point in frequency space for which the direction of selection changes its sign, i.e. at which the gradient of selection is zero. This point can be calculated by solving (or equivalently ) for i, and broadly speaking it is determined by the relative magnitudes of b and c. Third, for b > 1, c > 1 (or b < 1, c < 1) type A (or type B) always has the higher fitness irrespective of the composition of the population. This type is then always favoured by selection, which never changes direction.

For the remainder of this article, we focus on switching between coexistence and coordination games. More precisely, we choose b > 0 and c > 0 in (5.5). The coexistence game corresponds to σ = +1 and the coordination game to σ =−1.

Results

Sample trajectory of the dynamics

In figure 2a, we show a sample trajectory of a simulation in which a single mutant reaches fixation. The gradient of selection, , for the two games is shown in figure 2b and c. During periods when the environment is in the coexistence state (light background) the population fluctuates about the selection-balance point (dashed line), and during periods when the environment is in the coordination state (shaded background) the population is driven away from the selection-balance point. In the final period in the coordination state the mutant is driven to fixation.

Figure 2.

(a) A sample trajectory (time series) of the fraction of individuals of type A. White background corresponds to the environment being in the σ =+1 coexistence state, while the shaded background corresponds to the σ =−1 coordination state. Dashed line is the location of the point at which selection balances, which is the same in both states of the environment. (b) Gradient of selection in the σ =+1 coexistence state . Solid circle shows location of the point of selection balance, and arrows indicate the direction and magnitude of flow towards this point. (c) Gradient of selection in the σ =−1 coordination state . Empty circle shows location of the point of selection balance, and arrows indicate the direction and magnitude of flow away from this point. For the realization in panel (a) and the selection bias shown in (b) and (c), the payoff matrix elements are a = d = 1, b = 1 + 0.5σ and c = 1 + 0.9σ, the system size is N = 100, the selection intensity is β = 1, and the switching probabilities are p+ = 103 and p = 104. (Online version in colour.)

Fixation probability and conditional fixation time

In figure 3, we show the effect of switching the environment on the fixation dynamics. We choose c > b > 0. By equating the reaction probabilities (i.e. setting ), or equivalently equating the expected payoffs in equation (5.2), and looking at leading order terms in N, the gradient of selection is seen to change sign at i*/N ≈ b/(b + c) < 1/2. This point is closer to the extinction state of the mutant (i = 0) than to the fixation state (i = N). We next describe the key observations we make from these results, before we turn to their interpretation.

Figure 3.

(a) Fixation probability of a single mutant as a function of switching probabilities. The main panel shows simulation results (symbols; crosses correspond to σ(0) = +1 and circles to σ(0) =−1) for fixed p− = 0.01, along with the exact theoretical results (solid lines) from equation (4.4), and the effective theoretical result (dashed line) of equation (4.10). Inset panels show fixation probabilities from equation (4.4) over all combinations of p+ and p−. Left inset panel: initial condition σ(0) = +1. Right inset panel: σ(0) =−1. The horizontal lines correspond to the data shown in the main panel. (b) Mean conditional fixation time (in generations) of a single mutant as a function of switching probabilities. The main panel shows simulation results, as described above, for fixed p− = 0.01, along with the exact theoretical results (solid lines) of equation (4.8), and the effective theoretical result (dashed line) of equation (4.12). Inset panels show mean conditional fixation times from equation (4.8) over all combinations of p+ and p−. Left inset panel: initial condition σ(0) = +1. Right inset panel: initial condition σ(0) =−1. The horizontal lines correspond to the data shown in the main panel. The payoff matrix elements are a = d = 1, b = 1 + 0.5σ and c = 1 + 0.9σ, the system size is N = 50, and the selection intensity is β = 0.5. (Online version in colour.)

Fixation probability (figure 3a). The fixation probability in this example depends non-trivially on the rates with which the environment switches states; we find an optimal combination of switching rates, p+ ≃ p, for which fixation of a single mutant is most probable (figure 3a). The fixation probability is dependent on the initial state of the environment for .

Fixation time (figure 3b). Mean conditional fixation times show very little dependence on the initial state of the environment. The fixation time is small for p+ > p when the environment is found mostly in the coordination game, and large when the environment is mostly in the coexistence state.

Validity of the theoretical approach. As shown in both panels of figure 3, the theoretical predictions of equations (4.4) and (4.8), indicated by solid lines, are in convincing agreement with simulation data. Theoretical results from the model with effective transition rates (§4.2) reproduce the simulation data qualitatively. Quantitative agreement is obtained in the limit of large switching rates, but unsurprisingly, there are systematic deviations when switching is slow.

Interpretation

We now proceed and give an intuitive explanation for the observed effects.

Mean conditional fixation time is reduced as more time is spent in coordination environment. The behaviour of the fixation time can intuitively be understood from the deterministic gradient of selection of the two games (figure 2b,c). If fixation happens, it will generally be quicker in the coordination game than in the coexistence game [21,49]. This is due to the adverse selection bias in the coordination game at low mutant numbers (figure 2c). The more time the system spends in this region of adverse selection the less probable it is for the mutant to reach fixation. Thus, if fixation happens in a coordination game then it happens fast. In the coexistence game, on the other hand, the direction of selection is towards the balance point, as shown in figure 2b. The system can ‘afford’ to spend significant time in the region of small mutant numbers and still reach fixation eventually even after repeated excursions throughout frequency space. There is thus no need for fixation to occur quickly, and conditional fixation times can be long.

These observations make it plausible that the mean conditional fixation time will generally decrease when less time is spent in the coexistence game, which is exactly what we find in figure 3b. We have tested other choices of the parameters b and c for which the two games are a coexistence game and a coordination game, and we find that the behaviour of the mean conditional fixation times is robust under these changes.

Mean conditional fixation time is largely independent of the initial state of the environment. Systems started in the coordination environment will tend to reach extinction relatively quickly due to initial adverse selection, unless the environment switches to the coexistence state early on. Thus, the sample of runs that reach fixation started from the coordination-game environment will be dominated by runs in which the environment switches soon after the start of the run. Then we expect that the value of the mean conditional fixation time is close to the one obtained when starting in the coexistence game.

Dependence of fixation probability on the initial state of the environment. The data in figure 3a show that initiating the dynamics in the coexistence game favours fixation of the mutants for (when p = 0.01 is fixed), however, above this threshold the initial state of the environment has relatively little effect. The reason for this is as follows. When starting in the coordination game, selection pushes the mutant towards extinction. Hence, fixation is more probable if the initial state is the coexistence fitness landscape. Above p+ ≃ 0.1, the switching process of the environment is too fast for the initial condition to have any significant effect on the population dynamics. It is this regime in which we expect the effective description (§4.2) to approximate the system well. This is indeed confirmed in figure 3, the theoretical prediction of the effective theory, equation (4.10), agrees well with our simulation results in this fast-switching region.

Behaviour of fixation probability depends on location of the selection-balance point. If the environment is fixed to the coexistence-game state, fixation is more probable the closer the point of selection balance is to the fixated state (figure 4a). The location of this balance point is approximated by i*/N = 1/(1 + c/b), and so the fixation probability increases as c/b is decreased. In a fixed coordination-game environment, the reverse is the case. The range of adverse selection is to the left of the balance point, and so fixation is less probable the closer the point of selection balance is to the fixated state.

Figure 4.

(a) Illustration of selection bias in the two environments for different locations of the balance point; (b) value of p+ at which ϕ1, is maximal given p = 0.01 as a function of b and c. Remaining parameters are β = 0.5 and N = 50. (Online version in colour.)

For , i.e. a selection bias point close to i = 0, we therefore expect that the fixation probability will increase the more time that is spent in the coordination-game environment, i.e. ϕ1, is an increasing function of the probability p+ with which the system leaves the σ = +1 state (coexistence game). This is indeed what we find in simulations (data not shown). For , i.e. i* close to i = N, the reverse is the case. Fixation is more probable in the coexistence game (σ = +1), and the fixation probability is hence a decreasing function of p+ at fixed p−. Although we do not show the data here, this is again confirmed in simulations.

For b ≈ c the situation is less clear. The fixation probability will be comparable in both games if the environment is frozen. Two effects here conspire to produce a non-trivial outcome: These two effects taken together generate a maximum of the fixation probability at intermediate values of p+ ≈ p, which is exactly what we find in figure 3a. We would like to stress that the effect (i) is only present provided the selection-balance point is not too close to the extinct state. The phenomenon discussed under (ii) is present only if the selection-balance point is not too close to the fixated state. If the balance point is located too close to either boundary, the corresponding effect will be suppressed and the remaining effect dominates. One then finds monotonically increasing or decreasing dependences ϕ1,(p+).

(i) Consider the case in which the system is mostly in the coordination-game state, i.e. . It is plausible that an occasional switch to a coexistence game will make fixation more probable than in a constant coordination game. This is because the coexistence-game environment pushes the system away from extinction at low mutant numbers. In the regime of , we thus expect the fixation probability to increase as p+ is lowered. In other words, ϕ1,(p+) is a decreasing function at large p+.

(ii) Similarly, if the system is mostly in the coexistence-game environment , short periods of time in the coordination game can make fixation more probable. This is because selection at large mutant numbers is directed towards fixation in the coordination game. At , we expect ϕ1, to be an increasing function of p+.

To confirm our picture, we varied the payoff parameters b and c, and find the value of p+ that maximizes fixation probability for a given p = 0.01, as a function of b and c in figure 4b. The point of selection balance is approximately 1/(1 + c/b) (up to system-size corrections). The presence of diagonal structures in figure 4b shows that the behaviour of the fixation probability is dependent on the location of the selection-balance point. If this point is close to the fixation state i = N (, bottom-right in figure 4b), then the fixation probability is maximal for vanishing p+. If this point is close to the extinction state (, top-left in figure 4b), then the fixation probability is maximal for large p+. For intermediate locations of the selection-balance point (b ≈ c) fixation is maximized at a non-trivial combination of environment states.

The features observed in figure 3, i.e. the peak in the fixation probability and shape of the mean conditional fixation time as a function of p+, are found to be robust when the system size is increased. Fixation probabilities generally decrease with system size but the observed peak becomes sharper. The mean conditional fixation times in a frozen coexistence environment scale exponentially with N, whereas they scale as the logarithm of N in the coordination environment [50]. We observe these scalings in our system, with the mean conditional fixation time increasing exponentially with N for small p+, and increasing sub-linearly with N for large p+.

Mutation-selection equilibria under fluctuating environments

Mutations and stationary distributions

We now consider systems with mutations occurring during the dynamics. This removes the possibility of fixation and extinction. The combination of mutation, selection and noise can lead to non-trivial stationary states. We introduce mutation by modifying the discrete-time transition probabilities of equation (5.4) and now use where is the mutation rate. The transition probabilities are now non-zero, and so the states i = 0 and i = N are no longer absorbing.

The stationary probability ρ, of finding the system in state (i,σ) (i = 0,1, … ,N, ) is obtained as the solution of the balance equationThis equation is of the form , and it is solved by finding the eigenvector corresponding to the eigenvalue λ = 1 of the linear operator . In principle, this can be done analytically, but we use standard numerical packages to find the eigenvector. The stationary distribution for the state of the population is found by summing over all states of the environment, ρ = ∑ρ,. This solution is exact.

If the environment states are long-lived, the population will relax to the stationary state of the current environment before the next switching event. With this, one might expect that the overall stationary distribution is the weighted average of the stationary distributions one would obtain in the respective stationary environments. The stationary distribution in a single fixed environment, σ, can be found explicitly asThis can be derived for example from equation (6.2) assuming that the transition matrix of the environment is diagonal, μ = δ, where δ is the Kronecker delta. The average stationary distribution over many slow-switching environments can then be written aswhere ρ is the probability that the environment is in state σ.

Alternatively, if the switching probabilities per time-step are large one might expect the stationary distribution to be approximated by the distribution found in a system controlled by the effective transition rates, . These are obtained as described in §4.2, with suitable modifications to account for mutation. The resulting stationary distribution is found as

The distributions and ρ,eff are both approximations. To evaluate the validity of the assumptions leading to these approximations, we compute the distanceand similarly for ρ,eff, where is distribution of the population at time t obtained from simulations. To confirm our analytical approach, we also compute the distance of simulation data from the exact solution for ρ. We allow the system to run for a fixed time T, and we then use the time-averaged distance, to evaluate the accuracy of the approximations. We ignore the first half of the time series to remove remnants of the initial condition. We note that the time T, measured in generations, is equivalent to NT simulation time-steps and is chosen to be long enough such that the system relaxes to the stationary state before measurements start at time T/2.

We present results for the two-world scenario, where the environment switches between a coexistence game and a coordination game as described above. The stationary distribution of the environmental state is given by ρ =+1 = p/(p++p) and ρ =1 = p+/(p++p).

The stationary distributions of the population for the fixed environments (calculated using equation (6.3)) are shown in figure 5a. In a constant coexistence game (σ =+1), the stationary distribution is peaked about the point at which the gradient of selection changes sign, and in a fixed coordination-game environment (σ =−1) we find a distribution that is strongly peaked near the i = N state. The asymmetry is due to the imbalanced payoff matrix used, such that the basin of attraction for the i ≃ N state is much larger than for the i ≃ 0 state. For the parameters chosen in figure 5, the selection-balance point is at i* ≈ 18.

Figure 5.

(a) The stationary distributions in the single-environment coexistence game (dotted line) and coordination game (dashed lines) calculated from equation (6.3), along with the exact solution ρ (equation (6.2), evaluated numerically) and the ‘average’ (equation (6.4)) and ‘effective’ ρ,eff (equation (6.5)) approximate stationary solutions (solid lines). These distributions are for switching probabilities p+ = p = 103. (b) The time-averaged distance (equation (6.6)) between distributions obtained from an ensemble of 5 × 105 simulations, Q(t), and the analytic stationary distributions for symmetric switching, p+ = p−. Simulations run for T = 2 × 103 generations. Inset plots show the time-averaged distances over all switching-parameter space. Left inset panel shows the distance to the ‘averaged’ stationary solution (equation (6.4)), and right inset panel shows the distance to the ‘effective’ stationary solution (equation (6.5)). The payoff matrix elements are a = d = 1, b = 1 + 0.5σ and c = 1 + 0.9σ, the system size is N = 50, the selection intensity is β = 0.5, and the mutation probability is u = 0.02. (Online version in colour.)

For equal switching rates, p+ = p−, the averaged stationary distribution lies exactly in between the two single-environment distributions. The effective distribution is approximately uniform in the centre of the domain, with a lower probability of being found close to the domain boundaries. This reflects the fact that for equal switching probabilities the effective game is close to neutral, but frequent mutations push the population to the interior. The exact solution (equation (6.2)) matches the features of the single-environment distributions, with a peak at i ≃ N and at the coexistence point i*. Interestingly, this solution also predicts a peak at the i ≃ 0 state, a feature that is not seen in the single-environment distributions, or in the effective distribution.

As seen in the main panel of figure 5b, the exact solution is confirmed by simulations across many orders of magnitude of switching probabilities. Any deviations can be attributed to incomplete equilibration of the measure used in equation (6.6). For large switching probabilities, the effective stationary distribution, ρ,eff, matches the simulations. As expected the effective theory becomes inaccurate for slow switching, roughly below p ≃ 102, in our example. The weighted average stationary distribution, , shows the opposite behaviour. It is in reasonable agreement with simulations for slow switching, but shows systematic deviations when the switching process is too fast for the population to react adiabatically.

This picture is further corroborated by the data shown in the inset panels of figure 5b. The weighted average and the effective distributions accurately predict the stationary distribution obtained from simulations when the two switching rates are very disparate, i.e. or vice versa (top-left and bottom-right corners of the two insets). In these regions, the environment spends most of its time in one state, so that the model effectively reduces to the single-environment case. Simulation data, the exact solution, and both approximations then all collapse to the same result, the stationary distribution obtained in a single fixed environment.

The approach based on effective transition rates (right inset of figure 5b) is found to be accurate over a large range of switching probabilities away from the slow-switching scenario. Conversely, the weighted average distribution (left inset of figure 5b) becomes increasingly accurate if the dynamics of the environment is slow (p → 0).

Summary and conclusion

The dynamics of a population evolving under changing environmental conditions is an important concept in the study of bacterial populations. Previous work has focused on deterministic analyses [5], or on an environment following a continuous stochastic process [30]. Here, we have taken a different route and assumed that the environment switches between discrete states. We have developed the mathematical formalism to describe fixation properties in a general birth–death process in an environment fluctuating between an arbitrary number of states. The main results of this investigation are self-consistent expressions for the fixation probability of a mutant in a fixed-size population, as well as for the mean unconditional and conditional fixation times. For short-lived environments, we put forward an approximation based on effective transition probabilities.

As a specific application we discuss the fixation properties in the context of an evolutionary game in a two-world scenario. The two states of the environment then correspond to two different payoff matrices of the underlying games. Simulations confirm our exact solution over a wide range of switching probabilities. The approximation based on effective transition probabilities is seen to reproduce simulation data in the limit of fast switching.

Focusing on the case of switching between a coexistence game and a coordination game we find unexpected non-trivial behaviour of the fixation probability of a single mutant. We observe in our analytical results and in simulations that fixation can be more probable in a scenario in which the fitness landscape switches between the two games than in either of the two constant environments. We provide an intuitive explanation for this effect, and we have investigated in detail the circumstances under which this phenomenon can occur.

Adding mutations to the dynamics removes the possibility of fixation, but introduces non-trivial stationary states. We develop a method for the calculation of this distribution, along with approximations for long- and short-lived environmental states, respectively. These approximations are shown to agree well with simulations in their respective limits.

The general theory developed here now allows further investigation of evolutionary dynamics in time-varying environments. It provides a first mathematical characterization of the effects one may expect in such systems. The closed form self-consistent solutions will help to speed up future studies, and they may remove the need for extensive computer simulations.

While our work is mainly mathematical, we think that our theory can be used to interpret existing experimental studies such as those studied by Acar et al. [6]. For some biological systems, it may be more appropriate to use constant selection in each environment, as opposed to frequency-dependent selection. Our example of switching between coexistence and coordination games was chosen to illustrate the theory. We have also seen such cases lead to unexpected effects. We note that both types of game have been observed in systems of experimental evolution [14,51,52]. We hope the formalism we have developed will be useful to analyse models closer to other biological applications, and potentially to guide future experiments on evolutionary systems in time-dependent environments.

On a more general level constructing a mathematical theory of evolutionary dynamics is very much work in progress. An integral part of the evolution of microbes and higher organisms alike is frequency-dependent selection. At the same time, external factors determining the detailed mechanics of selection may vary in time. In this work, we have combined frequency-dependent selection, fluctuating environments and stochastic dynamics in finite populations into one model, and we have provided the analytical tools for its analysis. This, we hope, is a contribution towards a more complete understanding of evolutionary processes.