## Abstract

A population that adapts to gradual environmental change will typically experience temporal variation in its population size and the selection pressure. On the basis of the mathematical theory of inhomogeneous branching processes, we present a framework to describe the fixation process of a single beneficial allele under these conditions. The approach allows for arbitrary time-dependence of the selection coefficient *s*(*t*) and the population size *N*(*t*), as may result from an underlying ecological model. We derive compact analytical approximations for the fixation probability and the distribution of passage times for the beneficial allele to reach a given intermediate frequency. We apply the formalism to several biologically relevant scenarios, such as linear or cyclic changes in the selection coefficient, and logistic population growth. Comparison with computer simulations shows that the analytical results are accurate for a large parameter range, as long as selection is not very weak.

FOR adaptive evolution to proceed, it is not enough that new beneficial mutations enter a population. To complete an adaptive step, these mutations also need to escape stochastic loss due to genetic drift, get established, and finally rise to fixation. The fixation process of beneficial (or neutral or deleterious) alleles is one of the building blocks of population genetic theory and many of the key results on fixation probabilities and times date back to its early days. Two alternative mathematical frameworks have been developed to derive analytical expressions for these quantities: branching processes (Fisher 1922, 1930; Haldane 1927) and diffusion theory (Kimura 1962; Kimura and Ohta 1969). Today, a large body of literature exists to study fixation under various ecological scenarios and genetic conditions (reviewed in Patwa and Wahl 2008), such as the effects of population structure (Whitlock 2003) and spatial heterogeneity (Whitlock and Gomulkiewicz 2005) and interference due to selection on linked loci (Barton 1995) or due to epistatic interaction (Takahasi and Tajima 2005).

In this article, we consider the fixation process in a variable environment, leading to time-dependent selection coefficients and population sizes. Aspects of this problem have already been studied in previous work: In particular, the impact of various scenarios of demographic change (growth, decline, cycles) on the fixation probability has been treated in a series of articles (Ewens 1967; Chia 1968; Kimura and Ohta 1974; Otto and Whitlock 1997; Pollak 2000; Parsons and Quince 2007a; Orr and Unckless 2008). Studies on time-dependent selection mostly concentrate on stochastic fluctuations of the selection coefficient (Jensen 1973; Karlin and Levikson 1974; Takahata and Ishii 1975; Huillet 2011). Since the distribution of the selection coefficients is constant across generations, these models are still time-homogeneous in a probabilistic sense. In contrast, surprisingly little is known when the changes of the selection coefficient *s* = *s*(*t*) follow an explicit trend. Ohta and Kojima (1968) derive an expression for the fixation probability of a mutation with time-dependent selective advantage in the context of the evolution of chromosomal inversions. Apart from that, only particular functions have been considered: Kimura and Ohta (1969) discuss the case where selection decreases exponentially in time, and Pollak (1966) derives expressions for the fixation probability under two alternating selection pressures. No previous work seems to exist where both population size and selection strength are variable, although this is a generic case under realistic ecological conditions. Also, there does not seem to be an investigation of fixation or passage times in variable environments.

In the following, we present a formalism to describe the fixation process under a wide range of scenarios of environmental variation. We use branching processes in continuous time to derive analytical approximations for the fixation probability and the passage time needed for the mutant allele to reach some intermediate frequency *x*_{c}. After the introduction of our model, we describe how a general approximation for the fixation probability can be obtained from known mathematical results on inhomogeneous branching processes. Afterward, we discuss applications of this result to several biologically relevant scenarios. In the second part of the article, we introduce and apply a method to calculate the distribution of the passage time needed for a beneficial mutation to reach an intermediate frequency. The method works by combining the stochastic fluctuations from the branching approximation with the deterministic growth of the full model. This technique has been used before (for constant selection and population size) by Desai and Fisher (2007) in a model of clonal evolution. All analytical results are complemented by computer simulations, which are briefly described in a separate section. We close with a short discussion. In the *Appendix*, we discuss a generalized version of the model to include allele-frequency–dependent population demographies. We illustrate how the formalism can be used by applying it to an illustrative ecological scenario: the fixation probability of a “rescue mutation” in a population that is threatened by extinction (*cf*. Orr and Unckless 2008). Additional material is devoted to Supporting Information, File S1, Figure S1, Figure S2, and Figure S3. We discuss in some detail the scope and limits of the approach and the accuracy of the approximation. In section S2 of File S1, we present an alternative treatment to derive the fixation probability in a variable environment from a diffusion approach.

## The Model

We consider a large population of haploid individuals with time-dependent population size *N _{t}*. The population dynamics are modeled as a time-inhomogeneous birth–death process with birth and death rates

*b*(

*t*,

*N*) and

_{t}*d*(

*t*,

*N*):

_{t}The impact of the changes in the external environment on the population size is reflected in the explicit time-dependence of the rates on *t*. The dependence on *N _{t}* accounts for density-dependence [

*e.g.*, logistic: ]. We call the growth parameter. Obviously, the expected change of

*N*over a small time interval d

_{t}*t*reads

Consider now two alleles, a beneficial mutant allele *A* and the ancestral (resident) allele *a*, that segregate in the population at a single locus. Recurrent mutations in both directions are ignored. In general, birth and death rates might be different for residents and mutants. These rates can depend on time and on the (absolute) frequencies of both allelic types, allowing for general frequency-dependent selection. As a consequence, also the population dynamics depend on the allelic composition and cannot be described by Equation 1 anymore. We discuss this model in the *Appendix*. For the main part of the article, however, we assume that the rates are the same for mutants and residents and that all model parameters are independent of allele frequencies. This means in particular that selection is soft; *i.e.*, changes in the allelic composition due to selection or drift do not interfere with the population dynamics. Population growth and decline of the polymorphic population are then correctly described by Equation 1.

In this setting, selection is modeled as competitive replacement between individuals, which does not change the population size, and is implemented as follows: At per capita rate ξ(*t*, *N _{t}*) +

*s*(

*t*,

*N*), a mutant additionally reproduces and succeeds in replacing a randomly chosen individual from the population by its offspring. Residents do the same at rate ξ(

_{t}*t*,

*N*). Again, the selective advantage

_{t}*s*(

*t*,

*N*) of the mutant may thus depend on the external environment (modeled by the dependence of

_{t}*s*(

*t*,

*N*) on

_{t}*t*) and the population size (modeled by the dependence on

*N*). Changes in the number of mutants then occur at rates

_{t}The model corresponds to a continuous-time Moran model, but with a population size that may change in time. Putting *b*(*t*, *N _{t}*) =

*d*(

*t*,

*N*) = 0, ξ(

_{t}*t*,

*N*) = 1, and

_{t}*s*(

*t*,

*N*) =

_{t}*s*= const. reproduces the standard Moran model (Moran 1958a,b; Novozhilov

*et al.*2006). The free parameter ξ(

*t*,

*N*) has been introduced to our model to allow for easy interpolation to other models (see below) and additionally to make the analysis of density-dependent competition possible.

_{t}To further clarify the relation to other models, we calculate how the frequency of mutants *x _{t}* :=

*n*/

_{t}*N*changes over time. Let Δ

_{t}*x*be its change in an infinitesimal time interval d

*t*. The expectation and the variance of Δ

*x*are calculated to be(4a)(4b)with the time-dependent variance effective population size(5)In the last step we approximated

*N*+ 1 ≈

_{t}*N*and

_{t}*N*− 1 ≈

_{t}*N*(see section S3 of File S1 for the derivation of Equations 4a and 4b).

_{t}We see that the strength of drift, measured as , is proportional to the total rate of events in the model. The choice coincides with the strength of drift in the standard Moran model, while is consistent with the scaling in the Wright–Fisher model. In contrast to many diffusion or coalescent approaches, we do not rescale time with the effective population size (which would be impractical since itself depends on *t*). Generation time in the continuous-time Moran model is defined as the inverse of the total death rate of an individual, , and may again depend on time in our model.

## Fixation Probability

### Analytical theory

Following pioneering work by Haldane (1927) and Fisher (1930), there has been a long tradition in population genetics to calculate fixation probabilities by branching process methods (reviewed in Haccou *et al.* 2005 and Patwa and Wahl 2008). For the general time-dependent case, the relevant results have long been known in the mathematical literature (*e.g.*, Kendall 1948; Allen 2011). However, only specific cases (usually in the context of changing population sizes) have been discussed in the population genetics context (Ewens 1967; Otto and Whitlock 1997; Pollak 2000; Wahl and Gerrish 2001). We therefore give a brief outline of the general theory below and show how it applies to the biological problem at hand. Previous results are recovered as special cases.

The branching process approximation is based on the following reasoning: Initially, the fate of a new beneficial mutation arising in a population will be strongly determined by genetic drift. In most cases it will actually get lost again. Once the mutation has survived this early phase, it is, however, almost sure to get fixed given its selective advantage is large enough. To calculate the fixation probability it is therefore often sufficient to consider the stage at which the mutant population size *n _{t}* is still small relative to the total population size

*N*. In this early phase, the mutant individuals suffer nearly independent fates, as do the individuals in a Galton–Watson branching process (this assumption is precisely met in an infinite population). The extinction probability of the latter can therefore be used as an approximation for the probability that the mutation gets lost. Because a mutation in a finite population is in the long term either fixed or lost, the fixation probability is the complementary probability. (We exclude the unbiological case of a population that increases without bounds, where possibly neither wild types nor mutants become extinct.)

_{t}Ignoring terms proportional to *x _{t}* =

*n*/

_{t}*N*in the birth–death model (Equation 3), which corresponds to the limit

_{t}*N*→∞, leads to transition rates that are proportional to

_{t}*n*. Following standard practice in ecological modeling, we further ignore stochastic fluctuations in the population dynamics. This is done by replacing the stochastic variable

_{t}*N*by its deterministic approximation denoted as

_{t}*N*(

*t*), with dynamics(6)Inserting the deterministic solution

*N*(

*t*) into the rates for birth, death, and selection reduces the dependence of these rates on

*t*and

*N*to a dependence on

_{t}*t*only [, etc.]. We arrive at a branching process with time-dependent per capita birth and death rates:

As explained above for the birth–death process, the total rate of events determines the strength of genetic drift, while λ(*t*) − μ(*t*) = *s*(*t*) + *r*(*t*) corresponds to the absolute expected rate of increase of a small mutant population (its Malthusian fitness parameter).

To derive the extinction probability of this process, we follow Allen (2011, p. 278ff). Let *p _{i}*(

*n*

_{0},

*t*) be the probability that there are

*i*individuals at time

*t*when the process started with

*n*

_{0}individuals at time

*t*= 0. Using the Kolmogorov forward equation,(8)a differential equation for the probability generating function of the branching process can be derived (see Allen 2011, p. 279):

The solution is known from the mathematical literature (Kendall 1948; Allen 2011, p. 280) and given as(10)with

(11)To keep notation short, we introduce the abbreviations

(12a)(12b)The extinction probability *p*_{0}(*n*_{0},*t*) is immediately obtained from the generating function via

The probability that the mutation will eventually fix in the population is therefore given by(14)where we introduced(15)The last equality is valid for . This condition is met in all examples considered below. For a single mutant (*n*_{0} = 1), it thus holds:(16a)(16b)where we have used for the last equality and *N*_{e}(*t*) is defined analogously to Equation 5. A similar expression was also derived by Ohta and Kojima (1968). Restricted to a constant population size and a Poisson offspring distribution, their result is, however, less general.

The result depends on two independent parameters, which are compositions of three biologically relevant factors: the strength of selection given by *s*(*t*), the combined rate of birth and death events (λ + μ)(*t*), and the changes in the total population size [modeled by *r*(*t*) or *N*(*t*)]. In the above equation we formulated the result via two different combinations of these three variables. In the first version (Equation 16a) it is expressed in terms of the absolute rate of increase of mutants in the population (*s* + *r*)(*t*) and the total rate (λ + μ)(*t*) at which events happen, which defines the time scale of the problem and also quantifies the influence of drift. In the second formulation (Equation 16b) of the result we combined the birth and death rates and the changing population size with the time-dependent variance effective population size. The second decisive parameter is the selection coefficient of the mutation. Depending on the question to be answered, one or the other version is more favorable. In the first version, the correspondence between a mutation with time-dependent selective advantage and a mutation in a population of changing size can be easily seen: a mutation with time-dependent selective advantage *s*(*t*) in a population of (on average) constant size *b*(*t*) = *d*(*t*) = const. has the same chance to reach fixation as a mutation with constant selective advantage *s*_{0} in a population with time-dependent death rate *d*(*t*) and growth parameter *r*(*t*) = *s*(*t*) − *s*_{0}. [Note, however, that *s*_{0} must be >0, such that fixation is (almost) certain once the mutation has survived genetic drift.] The second version is closer to the traditional view in population genetics. It is advantageous if the variance effective population size is directly given and allows, in particular, also for the treatment of discontinuous changes in the population size.

We see that the fixation probability is independent of many details of the individual-level dynamics of the original process (Equation 3), which depends on four rates [*b*(*t*, *N _{t}*),

*d*(

*t*,

*N*),

_{t}*s*(

*t*,

*N*), and ξ(

_{t}*t*,

*N*)]. Further aspects of the stochastic model are ignored by our deterministic approximation for the population dynamics. In particular, the analytical results become independent of the particular form of density regulation. As an example, consider three scenarios: (1) a population with inherently constant size

_{t}*N*

_{0}with

*b*(

*t*,

*N*) =

_{t}*d*(

*t*,

*N*) = 0 and ξ(

_{t}*t*,

*N*) =

_{t}*c*; (2) a population with density regulation according to and

*d*(

*t*,

*N*) =

_{t}*c*, with initial size

*N*

_{0}, and ξ(

*t*,

*N*) = 0; and (3) a population without density regulation with , initial size

_{t}*N*

_{0}, and ξ(

*t*,

*N*) = 0. In all three cases, the deterministic dynamics of the population size are the same (

_{t}*N*(

*t*) =

*N*

_{0}) and the analytical predictions for the fixation probability coincide for arbitrary

*s*(

*t*,

*N*). Simulation results of all three scenarios indeed showed no significant difference, justifying the approximation (see section S1 of File S1). This observation agrees with findings by Parsons and Quince (2007a) that demographic stochasticity does not significantly influence the fixation probability of advantageous alleles [Parsons and Quince (2007a) discuss this issue for a population that starts in the vicinity of the dynamic equilibrium]. For concreteness, we use the notion of “constant population size” in the following to refer to the case of a strictly constant population size . We further set ξ(

_{t}*t*) = ξ = const. in all applications. If not stated otherwise, we use ξ

*=*1 (corresponding to the Moran model scaling) for the results shown in the figures.

A necessary condition for the branching approximation to yield meaningful results is that the fate of the mutation is decided as long as the mutant frequency *x _{t}* =

*n*/

_{t}*N*is still small. As

_{t}*x*increases, mutants are no longer independent in the original birth–death process (Equation 3), and both processes differ significantly from each other. In particular, the birth–death process has a second absorbing boundary at

_{t}*x*= 1. As a consequence, neutral or even deleterious mutations, too, can become fixed by genetic drift. In the corresponding branching model, an upper absorbing boundary does not exist. Consequently, neutral or deleterious alleles must go extinct in the long term. For an allele with a general time-dependent selection coefficient, the branching process approximation is thereby valid only if the allele is “sufficiently beneficial on average.” We can formalize this condition as follows: Since divergence of the (first) integral in Equation 15 leads to the wrong prediction of a zero fixation probability, we need to require that the integral converges. If we assume that μ(

_{t}*t*) is bounded below by a constant

*C*

_{μ}> 0, a necessary condition for the convergence of the integral is divergence of the integral in Equation 11. More precisely, if we assume that μ(

*t*) has bounds

*C*

_{μ}> 0,

*C*

^{μ}such that

*C*

_{μ}< μ(

*t*) <

*C*

^{μ}, it must hold that for some constant

*C*> 1. A prominent example where this condition is not fulfilled is a mutation with an exponentially decreasing selection coefficient in a population of constant size [Kimura and Ohta (1969), cf. also Pollak (1966) and Ohta and Kojima (1968) for mutations that lose their advantage over time]. The condition is, however, not sufficient to obtain good results: It is also necessary that the contribution of times beyond the initial phase be negligible. In particular, changes in the environment at times larger than the fixation time must not have an impact on the results. As shown in section S1 of File S1, the approximation will usually be excellent if the estimate of

*p*

_{fix}from Equation 16 fulfills . By a similar reasoning it is clear that the requirement that the mutation be beneficial all the time is at this point unnecessarily strict. It is sufficient that the extinction probability is in both processes—the original birth–death process and the branching approximation—negligible once the number of mutants has reached a certain size.

### Applications

#### Constant selection and constant population size:

Setting λ(*t*) = ξ + *s* and μ(*t*) = ξ yields for the fixation probability

The well-known result found by Haldane (1927) is therefore reproduced for *N*_{e} = *N*. In the special case of a constant environment, it is possible to calculate the exact fixation probability from the transition matrix defined via Equation 3 (see Ewens 2004, p. 90). One obtains(18)which shows that the branching approximation is very accurate for large values of *sN*/ξ. Furthermore, it is possible to calculate the fixation probability of a deleterious mutation with selective disadvantage −*s*. To do so, we switch the roles of the transition rates defined in Equation 3. Every path leading to fixation *n* = *N* has now a chance to be realized, which is by a factor of (ξ + *s*)^{−}* ^{N}* lower than the corresponding path for a beneficial mutation. It immediately follows that

#### Constant selection and changing population size:

Otto and Whitlock (1997) analyzed the fixation probability for several important scenarios of demographic change. Our general result (Equation 16b) contains all those scenarios, among arbitrary others. In addition, it provides some insight into the influence of the individual-level dynamics on the fixation probability. As an example, we consider a population that follows logistic growth (or decline) until it has reached its new carrying capacity *K*. There are different ways to describe these global dynamics at the individual level. We discuss two possibilities, which arise naturally in a biological context: The first one assumes that a decreasing availability of resources per individual leads to a lower birth rate while the death rate stays constant (*e.g.*, fertility is reduced). The second one assumes that the same circumstances lead to a higher death rate while the birth rate stays constant. The selection coefficient *s* of the beneficial mutant is constant in both cases. In the first scenario, the birth and death rates are given by(20)The total population size thus changes according to(21)and using Equation 16b yields(22)with .

In the second scenario, birth and death rates are then given by

(23)While the deterministic dynamics at the population level are the same for both scenarios (given by Equation 21), this is not true for the fixation probability. From Equation 16b, we obtain for the fixation probability in the second scenario(24)*i.e.*, a reduced probability relative to the first scenario. This is explained by the fact that genetic drift is stronger in the second scenario. For small values of s and *r* both fixation probabilities are approximately the same, and reproduce the result found by Otto and Whitlock (1997) if we choose *b* + ξ = 0.5 such that the variance of the increase in the mutant number Var[Δ*n*|*n _{t}*] is the same as in their model.

Note that for a sudden jump in population size, the result of our theory coincides with the result derived by Otto and Whitlock (1997) (where selection is effective during the change in population size), while Wahl and Gerrish (2001) consider a slightly different situation (where selection is switched off during the bottleneck; *cf*. Patwa and Wahl 2008).

#### Linearly increasing selection:

When environmental conditions develop continuously in a given direction, the selective advantage of a mutation may gradually increase during the fixation process. Let us assume that the total population size stays constant and that the selection coefficient of the mutation linearly increases in time, thus *s*(*t*) = *s*_{0} + *s*_{1}*t*. We obtain the following expression for *p*_{fix} from Equation 16a,(25)where is the complementary error function. For the special case *s*_{0} = 0, the result simplifies,

In Figure 1 the analytical results are compared to simulation results, showing very good agreement. The fixation probability increases significantly with *s*_{1}. Generally,(27)with . The function *f*(*y*) is shown in the inset of Figure 1. For *y* = 1 it is evaluated to be *f*(1) ≈ 1.52. This means that if *y* ≈ 1, *i.e.*, , we observe an increase in the fixation probability of ∼ in comparison to a mutation with constant selective advantage. *E.g*., for initially moderately strong selection the fixation probability is still increased by ≈50% if selection increases as slightly as *s*(*t*) = 0.01 + 0.0001*t*.

#### Periodically changing selection:

Cyclic environmental changes, such as seasonal changes or cyclic climate fluctuations (like the El Niño phenomenon), are frequent in nature. In the following, we consider strictly periodic changes, although our theory does not rely on this condition. Let us assume that again the total population size stays unaffected, but that the selective advantage of a mutation changes periodically with time; *e.g.*,

Depending on the parameter values, it is thus possible that the mutation is disadvantageous at certain periods of time. We obtain for the fixation probability:(29)The integral can be evaluated only numerically. Comparison to simulated data (see Figure S4) shows that the theory provides an accurate prediction of the fixation probability also for scenarios in which the mutation temporarily becomes disadvantageous. For *s*_{0} ≤ 0, however, it predicts a fixation probability of zero and therefore underestimates the true value.

To slightly reduce the parameter space, we concentrate now on the special case *s*_{max} = 2*s*_{0}, *i.e.*,

Figure 2 shows how the fixation probability changes with ω for various values of φ. For small and intermediate values of ω, the value of φ has a strong impact on the result. If ω increases, the fixation probability converges to the value for a mutation with selective advantage *s*(t) = *s*_{0} for all values of φ. This behavior is further illustrated in Figure 3, which shows the fixation probability and the initial selection strength in dependence of φ for various values of ω. In the limit ω → 0, the fixation probability equals *s*(0)/(1 + *s*(0)) ≈ *s*(0). It therefore follows approximately the curve of *s*(0) for small values of ω. For large values it converges to *s*_{0} for all values of φ. In an intermediate regime , more complex behavior is found: Here, not only the initial value *s*(0) but also the following time development of the selection strength [*i.e.*, not only *s*(0) = *s*_{0}(1 + cos(φ)), but φ itself] becomes important. The extrema of *p*_{fix}(φ) are attained for smaller values of φ than the extrema of *s*(0). A straightforward calculation shows furthermore that(31)where φ_{e} is the value where the extremum is attained. We thus see that the fixation probability at the extrema φ_{e} is the same for cyclic selection and constant selection with *s* = *s*(0) = *s*_{0}(1 + cos(φ_{e})).

## Time to reach an intermediate frequency x_{c}

### Analytical theory

In models of the adaptive process, it is often necessary to know the time that it takes for a successful mutation to become established and to reach a certain threshold frequency (*e.g.*, Desai and Fisher 2007; Kopp and Hermisson 2009a). It is well known that the deterministic model of the allele frequency increase yields a poor approximation even for the expected value of this time. The reason is that a mutation that has survived the initial phase will on average have grown faster during this phase than predicted by the deterministic model. When the frequency is finally large enough that stochasticity can be neglected and the path be modeled deterministically, the frequency will thus be larger than if it had *always* grown following the deterministic path. In that phase, it is well described by the deterministic path that is not started from a single individual, but from an (on average) larger “effective initial population size.” The method presented here builds on work by Cohn and Jagers (1994) and Desai and Fisher (2007). It consists of subsuming all the stochasticity of the path under this effective initial population size as a single random variable and then modeling the path deterministically. The procedure consists of two steps: In a first step the distribution of the effective initial population size is estimated via a branching process. In a second step the deterministic approximation of the full birth–death model is used to describe the allele frequency path starting from the effective initial population size.

Let us again consider the phase in which the mutation is rare and in which the dynamics can be described by the time-inhomogeneous branching process (Equation 7). The key to the method is that it is possible to separate stochastic fluctuations and deterministic growth for this process. Here, deterministic growth coincides with the time-development of the expected number of mutants:

(32)In the following, we restrict to the case *n*_{0} = 1. Define now a new random variable:(33)ν* _{t}* describes all the stochasticity that has accumulated in the branching process until time

*t*. Crucially, a theorem by Cohn and Jagers (1994) guarantees that for

*t*→ ∞, ν

*converges (almost surely) to a positive random variable ν that summarizes the entire stochasticity of the process. Since*

_{t}*n*=ν exp(ρ(

_{t}*t*)) in the limit

*t*→ ∞, we can interprete ν as the random initial population size of an ensemble of deterministically growing paths that approximate the original process

*n*for large

_{t}*t*.

For the fixation process, in particular, we are interested in the distribution of ν conditioned on non-extinction. We proceed as follows: From the probability generating function *P*(*z*, *t*) ≡ *P*_{1}(*z*, *t*) (see Equation 10), we obtain for the probability *p _{n}*(

*t*) ≡

*p*(1,

_{n}*t*) to have

*n*individuals at time

*t*(proof by induction, see section S3 of File S1):

Conditioning on non-extinction [which requires the biologically meaningful condition *p*_{0}(*t*) ≠ 1] leads to

By induction, we find(36)*i.e.*, for large times the *k*th moment grows in the conditioned as well as in the unconditioned process as *e ^{kρ}*

^{(}

^{t}^{)}. In particular, we obtain for the expected value(37)[Since , a comparison with Equation 32 confirms that on average the conditioned path grows faster than the unconditioned path.] The distribution function of is immediately obtained from the distribution function of :

As , we obtain in the limit the stationary distribution

(40)This effective initial mutant population size may now be used as the starting value for the deterministic solution of the full model. If the birth and death rates are the same for mutants and residents, the latter can be obtained from Equation 4a and is given by a generalized logistic,(41)where *x*(*t*) is the mutant frequency and

We want to calculate the time needed to reach an intermediate frequency *x*_{c}. Here, intermediate frequency means a frequency at which the dynamics are well described by Equation 41, *i.e.*, not too close to either 0 or 1. Using Equation 41, we can express the random variable defined via in terms of ν:

For the definition of to be unique for all *x*_{c,} the condition *s*(*t*) > 0 (up to single points) is necessary. It is possible that the process grows so quickly that the effective initial mutant population size is already larger than *x*_{c}*N* [*cf*. the exponential tail in the distribution function for the effective initial population size ν (Equation 40)]. In this case, we have set so that the distribution of will have a mass on . However, if *x*_{c} is not too small, this is very unlikely and the mass will be small. As the deterministic path (Equation 41) ignores stochastic fluctuations, there is only a single passage time in the deterministic approximation in contrast to the stochastic path in which the frequency *x*_{c} can be hit several times. as defined in Equation 43 is best interpreted as the average over all times at which the path crosses the frequency *x*_{c} from lower to higher values.

By a parameter transformation ( using Equation 43), the distribution of is found from Equation 40:(44)Equation 44 constitutes the main result of this section.

β(*t*) is a function of *s*(*t*), and *p*_{fix} can be expressed in terms of *s*(*t*) and *N*_{e}(*t*) (Equation 16b). Like the fixation probability, the distribution of therefore depends on two parameters that summarize the time-dependence of the model. For the fixation probability, we have seen in Equation 16a that there is a second way to summarize the dynamics in terms of the absolute rate of increase (*s* + *r*)(*t*) and the total rate of birth and death events. This is not possible for the distribution of , linked to the fact that variable selection and demographic change are not equivalent in this case. While the selection function *s*(*t*) has a strong influence on the deterministic part of the frequency path *x*(*t*), changes in the total population size (with equal effect on mutants and residents) have no influence on the latter.

To calculate moments of the distribution, we perform the parameter transformation(45)Solving for yields(46)and we obtain for the moments(47)For constant selection and therefore β(*t*) = *st* we obtain for the mean of :

### Applications

#### Constant selection and constant population size:

For the distribution of the passage time , we obtain with (see Equation 17)

(49)Plots of the probability density for *x*_{c} = 0.5 and various values of *s* are depicted in Figure 4A. With increasing selection strength the distribution gets narrower and is shifted to the left.

In the particular situation of constant selection and population size, additional results can be obtained. First, we note that can also be interpreted as the age of a derived allele that is currently found at frequency *x*_{c} in the population. This is a consequence of the time-homogeneity of the model. Second, the distribution of the time to reach fixation at *x* = 1 starting from a frequency (1−*x*_{c}) equals the distribution of (see section S3 of File S1 for an explanation). In particular, the times needed from 0 to 0.5 and from 0.5 to 1 follow the same distribution. Since both random variables are independent, we obtain the probability density of the whole fixation time from the density *p*(*T*_{1/2}) as

An alternative expression for the distribution of *t*_{fix} has been derived before by Wang and Rannala (2004), using diffusion techniques. We note that our approximation (Equation 50) is simpler than the series expansion in terms of Gegenbauer polynomials, using the eigenvalues of the oblate spheroidal angular function and the intermediate coefficients of the spheroidal harmonics obtained by these authors. As discussed in section S1 of File S1, it is nevertheless highly accurate if selection is not very weak or the population size small. The cumulants of the fixation time are just two times the cumulants of the time to reach frequency 0.5. For the mean fixation time, in particular, we obtain(51)where γ ≈ 0.577 is the Euler constant.

For the higher cumulants, approximations of leading order in *sN* can be obtained by the approximation exp(*sT*) + 1 ≈ exp (*sT*) in the distribution function (Equation 49) and by extending the integral (Equation 47) to infinity. We state the results for the variance, the skewness, and the kurtosis,(52a)(52b)(52c)where ζ(*z*) is the Riemann Zeta function and ζ(3) ≈ 1.202 is the Apéry constant. It can be shown that all cumulants of order ≥2 are to leading order of the form(53)with some numerical constant *c _{j}*.

Since *s* > 0 is required for the branching approximation, deviations from the exact values must be expected for weak selection. To estimate the precision of our approximation in this case, we compare the result for the mean fixation time to the result obtained from the diffusion approximation. For this purpose we scale time with *N*(*t*_{fix} → τ_{fix} = *t*_{fix}/*N*) and take the diffusion limit (*i.e.*, *s* → 0, *N* → ∞ with α := *sN*/ξ remaining constant) of Equation 51. We obtain(54)This is in perfect agreement with the approximation for the mean fixation time found in Hermisson and Pennings (2005), using the diffusion approximation. As shown in section S2 of File S1, the deviation of our branching process approximation from the exact diffusion result is exponentially small in α/2.

Figure 4, B and C, shows the mean time to reach a frequency *x*_{c} in dependence of *x*_{c}. Significant deviations between the theoretical predictions of Equation 49 and simulation results are found only for values of *x*_{c} that are extremely close to 1 (or 0). Generally, it holds that the approximation becomes less good for *x*_{c} ≳ 1−1/α (or *x*_{c} ≲ 1/α), as in this regime the deterministic frequency path (Equation 41), which we used in the derivation, is not a good approximation to the real frequency path of the mutation (see Ewens 2004, p. 167f).

It is interesting to note that the distribution of is the same for beneficial mutations with selective advantage *s* and deleterious mutations with selective disadvantage −*s* (*cf*. Ewens 2004): For every point of a particular path leading from *n* = 1 to *n* = *n*_{c}, the sum of the birth and the death rate is the same for both scenarios. This implies that the distribution of the run time of this particular path is the same. The chance to be realized differs for each such path by a constant factor If we therefore condition on *n* = *n*_{c} being reached, the distribution of the time to reach *n*_{c} is the same for beneficial and deleterious mutations.

#### Constant selection and changing population size:

Using the expression for *p*_{fix} calculated for a logistically growing population, we immediately get the distribution of the passage time . Figure S5 shows a comparison to simulation results for the second scenario. With Equation 48 we obtain the result for the mean of :(55)

#### Linearly increasing selection:

Putting the expression obtained for *p*_{fix} (Equation 25) into Equation 44, one obtains the distribution of the time at which a fraction *x*_{c} of the population consists of mutant individuals. The distribution can be used to calculate the mean of numerically. A comparison to simulation results for the distribution and the mean of in dependence of *s*_{1} can be found in Figure S6. From Figure 5, it can be seen that for small values of s_{0} even a very slight increase of the selection strength over time leads to a drastically smaller mean value of in comparison to an environment with constant selection *s* = *s*_{0}.

#### Periodically changing selection:

Figure 6 shows the mean time to reach frequency *x*_{c} = 0.5 for a periodic selection coefficient . As for the fixation probability, the initial selection strength *s*(0) is decisive if the environment changes very slowly; φ itself gains importance with increasing (but still small) ω. If the environment changes fast, becomes independent of φ and converges to its value for constant selection s(t) = s_{0}. Convergence is much faster than for the fixation probability.

## Simulations

To test our analytical results, we performed individual based computer simulations for which we use a Gillespie algorithm (Gillespie, 1977). Events happen at rate where ξ(*t, N _{t}*)

*s*(

*t, N*),

_{t}*b*(

*t, N*) and

_{t}*d*(

*t, N*) are assumed to be constant between events. Once the time of an event is fixed, which kind of event takes place is determined using the respective probabilities. Subsequently, the values of

_{t}*N*and

_{t}*n*are updated and the rates are set to their new values (additional update steps between events did not change the results).

_{t}For most of the simulation runs, we used the first passage time of a given frequency *x _{c}* to determine , i.e., we neglected fluctuations around

*x*, which for large values of α and not too small or large values of

_{c}*x*has no significant effect. For the data shown in panel C of Figure 4 and in Figure S2 and Figure S3, where we pushed the boundaries of the theory, we took the average over all passage times, more precisely the average over all times at which the path crossed the frequency

_{c}*x*from lower to higher values (cf. section on the analytical theory).

_{c}For all simulation results, the number of runs was chosen sufficiently large that the standard error bars vanish in the symbols. All programs were written in C, making use of the Gnu Scientific Library (Galassi *et al*., 2009). We used Mathematica (Wolfram Research, Champaign, USA) for all numerical evaluation of integrals.

## Discussion

Adaptation is the evolutionary response of a population to an environmental challenge: variation in the environmental conditions leads to altered selection pressures and changing population sizes. In nature, these changes occur on all time scales, from rapid shifts within a single generation to long-term geological trends. It therefore seems natural that models on the genetics of adaptation, too, should account for the ecological dynamics that drive the process. In contrast to this expectation, however, the vast majority of published studies with a genetic focus assume a constant environment (reviewed, *e.g.*, in Orr 2010). They rely on the idea that a fast—almost instantaneous—change in the environment is followed by a period of environmental stasis. If this period is long compared to the total time it takes for one or several beneficial mutations to appear and to rise to fixation, ecology and evolution are effectively decoupled. In many cases of ecological interest, however, this separation of time scales is not appropriate. In this case, the parameters of the evolutionary model (such as selection coefficients or population sizes) are turned into time-dependent variables, which must be determined from an underlying ecological model.

In a series of recent publications by several authors, the impact of the ecological dynamics on the genetics of adaptation has been studied for the so-called moving optimum model (Bello and Waxman 2006; Collins *et al.* 2007; Kopp and Hermisson 2007, 2009a,b). The consensus, at least for this model, is that the adaptive process is strongly affected by the dynamics of the selective environment. In this article, we present a more detailed treatment of the most basic aspect of the genetics of adaptation: the fixation process of a single beneficial mutation in a variable environment. All relevant parameters of the process, *i.e.*, selection coefficient, population size, and the variance in offspring number (or, equivalently, the birth and death rates of wild types and mutants) may depend explicitly on time. An example of how this can result from an explicit ecological model is given in the *Appendix*, where we discuss the fixation probability of a “rescue mutation” in a population that is otherwise doomed to extinction. Again, the results are strongly influenced by the ecological dynamics, leading to qualitative differences relative to previous studies that assume constant selection (Orr and Unckless 2008).

Even if the external environment is constant, individual alleles in a population can experience variable selection pressures if multiple selected alleles segregate in the population and if these alleles interfere due to either epistasis or linkage (*cf*. Hartfield and Otto 2011). Examples include the evolution of compensatory mutations or classical problems of clonal interference, where the beneficial mutation rate is high relative to the recombination rate. Another potential application is adaptive gene flow across a genetic barrier, where adaptations need to be purged from linked deleterious alleles by recombination.

Using branching process techniques, we obtain analytical approximations for the fixation probability and the distribution of the time for the mutant to reach a given frequency *x*_{c} in the population.

### Fixation probability

The derivation of fixation probabilities of rare beneficial mutations from a supercritical branching process is a standard approach in population genetics. In particular, results have previously been obtained for the most important modes of demographic changes (*e.g.*, Otto and Whitlock 1997; Pollak 2000). A look into the mathematical literature reveals, however, that a general formalism with arbitrary time-dependent birth and death rates has been available since the work of Kendall (1948). In contrast to most of the previous studies in population genetics, which use a branching process in discrete time (following Haldane 1927), this approach is based on a continuous-time framework, which simplifies some technical aspects. Our adaptation of this formalism to the genetic context leads to a compact, yet general formula for the fixation probability *p*_{fix} (Equations 16a and 16b) that covers previous results as special cases. In section S2 of File S1, we show how an analogous result can also be obtained using the diffusion approach.

It turns out that the ecological dynamics affect *p*_{fix} through two independent variables. In two alternative formulations of the result, these can be either the time-dependent birth and death rates of rare mutants (Equation 16a) or the selection coefficient *s*(*t*) and the variance effective population size *N*_{e}(*t*) (Equation 16b). Our applications to various scenarios show that relative to the case with constant selection *s*_{0}, consistent changes in the selection coefficient Δ*s* per generation of the order of have a strong effect on *p*_{fix.} This is to be expected since, for constant selection, the fate of a new mutation (fixation or loss) is decided once the mutant allele reaches a frequency on the order of (*Ns*_{0})^{−1}. According to Equation 48, this will take, on average, in the order of generations. The observation is of practical importance since it shows that predictions about fixation probabilities cannot be based on short-term fitness assays. Unmeasurably small fitness changes across generations may have a large effect. The limitations of the branching approach are the usual ones: The fate of the mutation must be decided while the mutant frequency is small and the independence assumption of the branching process applies. Deviations from the simulation results are found in the case of mutations that are almost neutral on average, with fixation probabilities *p*_{fix} ≲ 10/*N* (*cf*. File S1).

### Time *T*_{xc} to reach frequency *x*_{c}

The usual method in population genetics to derive fixation times or (more generally) first-passage times for mutant alleles to reach certain frequencies is diffusion theory. Since this proves difficult for the time-inhomogeneous case, we again turn to a branching process approach. We can use the fact that (almost) all stochasticity of a mutant trajectory *x*(*t*) is due to its early phase while the mutant number is still small. We can therefore (approximately) describe this stochasticity in the branching framework and combine it with the deterministic growth of the full process. For constant selection and population sizes, this idea has previously been used by Desai and Fisher (2007) in the context of a model for clonal evolution.

For the general case, the crucial step of the method has again been anticipated in the mathematical literature (Cohn and Jagers 1994). There, it has been shown that a clean separation exists for the random variable *n _{t}* of an inhomogeneous supercritical branching process into a growth function that describes the asymptotic growth and a time-constant random variable ν that describes the stochastic fluctuations in the large-time limit. We can interpret ν as the effective initial size of the mutant population. It turns out that the distribution function of ν,

*P*(ν ≤ ν

_{0}|not extinct) = 1 − exp(−

*p*

_{fix}ν

_{0}) (Equation 40), is pleasingly simple even in a variable environment. In particular, the impact of the ecological dynamics is conveniently summarized in the fixation probability

*p*

_{fix.}When this initial size is combined with the deterministic growth of the full model, an approximation for the distribution function for is obtained by a simple transformation of the probability density (Equation 44).

Computer simulations show that the results from the branching approach are usually highly accurate as long as selection is not very weak. This is also confirmed by a comparison with previous results for the expected value of the fixation time *T*_{1} for constant *s* and *N* (Kimura and Ohta 1969; Ewens 2004) derived from the diffusion approximation: The order of the error term for the average fixation time is *o*(exp(−*Ns*/2)) (see File S1). Previous results also exist for the entire probability density of *T*_{1} in a constant environment (Wang and Rannala 2004). Here, the branching approximation leads to a much simpler analytical result.

As a major caveat of the method, we note that for variable *s* or *N*, we do not obtain a closed expression for the time to reach full fixation *x*_{c} = 1. More generally, the approximation is accurate only as long as 1/(*sN*) ≲ *x*_{c} ≲ 1 − 1/(*sN*). For *x*_{c} close to 0, estimates can easily be found directly from the conditioned branching process, if needed. For mutant frequencies near 1, the fixation process enters another stochastic phase when the wild-type alleles become rare. Although this phase can be described, in principle, as a subcritical branching process, its details depend on the random time when this stochastic phase is entered, which complicates the treatment. In natural systems, an allele frequency of *x*_{c} = 1 is not an absorbing state and may never be reached in the face of back mutation, migration, or ongoing adaptation. For many practical applications, other thresholds (such as 5%, 50%, or 95%) are therefore more relevant and have previously been used. For example, Desai and Fisher (2007) use a threshold of *x*_{c} = 1/(*Ns*) to characterize their establishment time τ_{est}. Kopp and Hermisson (2007, 2009a) use *x*_{c} = 0.5 as the critical value where the mutant comes to dominate the population in an analysis of the order and step sizes of adaptive substitutions.

Although we have formulated our results for haploids, they also apply to diploids as long as the mutant is not fully recessive (*i.e.*, *hs* must be sufficiently large on average). Only the selection coefficients (or birth and death rates) of heterozygotes enter into the stochastic part of the model described by the branching process. In contrast, the deterministic frequency path used in the derivation of depends on the fitness values of both heterozygous and homozygous mutants. In the purely recessive case, the dynamics of rare mutants can no longer be described by a supercritical branching process. Diffusion methods (for constant selection) show that also the scaling of the fixation time with the selection parameters is altered in this case (Ewing *et al.* 2011).

To summarize, our results show that inhomogeneous branching processes provide a powerful framework to describe the fixation process under a wide range of ecological scenarios and genetic conditions. We therefore hope that the methods and results can provide a step to combine ecological and genetic modeling traditions to study the genetics of adaptation under realistic ecological conditions.

*Note*: After acceptance of this article we became aware of a recent preprint by Waxman on the fixation process under variable selection and demography [see D. Waxman (pp. 907–913) in this issue]. Both articles are complementary: While our approach focuses on analytical approximations for the fixation process of a definitely beneficial mutant, Waxman (2011) presents numerical methods to derive fixation probabilities for alleles with arbitrary selection coefficients.

## Acknowledgments

We thank Nick Barton, Ulf Dieckmann, Peter Pfaffelhuber, and Matthew Hartfield for helpful discussions and Lindi Wahl and two anonymous referees for many useful comments. This work was made possible with financial support from the Vienna Science and Technology Fund and from the Deutsche Forschungsgemeinschaft, Research Unit 1078, *Natural selection in structured populations*.

## Appendix

#### Fixation in General Ecological Models

In this *Appendix*, we describe how the branching process approach can be applied to a general ecological scenario where the genetic composition of the population and the population dynamics are mutually dependent. Assume that at time *t* a population consists of *w _{t}* wild-type “residents” and

*n*mutants. Both types reproduce and die at different per capita rates, which may depend on the number of residents and mutants,

_{t}*w*and

_{t}*n*, and on external factors that are independent of the population, but also may change over time. A general framework, which takes the full dynamics with all kinds of interactions between types into account, is given by the following scheme of transition rates:

_{t}Higher-order interaction terms with transition rates proportional to any polynomial in *n _{t}* and

*w*can be added as needed. The model corresponds to a general two-dimensional birth–death process with time-dependent coefficients. For such a model, the corresponding deterministic evolution of

_{t}*n*and

_{t}*w*is generically given by a system of two coupled differential equations. In most cases, it is impossible to solve this system analytically. Since the determinstic frequency path of the mutant is needed for the theory developed to calculate , derivations for must resort to numerical solutions to apply the methods outlined in the main text. In contrast, an analytical approximation for the fixation probability (and thus also the distribution of the effective initial population size ν,

_{t}*cf*. Equation 40) is often still possible.

In the branching limit of rare mutant individuals, all interactions of either residents or mutants with (other) mutant individuals can be neglected. Mathematically, this corresponds to the approximation *n _{t}* ≈ 0 in the birth and death rates of the residents and negligence of all terms of order (

*n*) > 1 in the transition rates of the mutants. As a consequence, the evolution of the residents is decoupled from the mutants. In the deterministic limit, the time development of the resident population size

_{t}*w*(

*t*) ≈

*N*(

*t*) is then given by an ordinary differential equation that can frequently be solved explicitly. We can then describe the dynamics of rare mutants by a one-dimensional branching process. Inserting the deterministic solution for

*w*(

*t*) into all mutant transition rates of linear order in the mutant number

*n*, the theory developed in this article can be applied by choosing λ(

*t*) and μ(

*t*) appropriately. Note, however, that the first expression for

*p*

_{fix}(Equation 16a) must be used, since the step leading from Equation 16a to Equation 16b is not possible if the growth rates of mutants and wild types are not the same. Previous results on fixation probabilities in two-dimensional birth–death chains (with constant coefficients) have been obtained by Parsons and Quince (2007a). Their approach is based on singular pertubation methods and has the advantage that it can be applied to neutral and deleterious mutation, too (Parsons and Quince 2007a,b), which is beyond the reach of the branching model. However, it cannot easily be extended to the time-inhomogeneous case that is our focus here.

#### Rescue mutation

Consider, as an example, the following scenario: A population is subject to a new and severe selection pressure (*e.g.*, new insecticide, drug, parasite, competitor, …). As a consequence, the population size decreases, and selection will drive the population to extinction unless a rescue mutation can be established that confers immunity to its carriers. A similar scenario has been discussed by Orr and Unckless (2008).

A plausible ecological model for this situation could be as follows: The original resident population evolves under logistic population dynamics with growth parameter *r* and time-dependent carrying capacity *K*(*t*). Due to the environmental challenge, the carrying capacity decreases exponentially like *K*(*t*) = *K*_{0}exp(−*at*). At the individual level, the model can be described in the scheme of Equation A1,(A2)with all other coefficients and equal to zero. As long as the number of beneficial mutants is small, the wild-type population size changes according to

Assume now that a beneficial mutation—if it survives stochastic loss—restores population size to its original carrying capacity *K*_{0} = *K*(0); *i.e.*,(A4)(and again all other coefficients are equal to zero). The fixation probability of a single mutation arising τ time units after the decline in the carrying capacity started is then given by

In Figure A1 the fixation probability of a single mutation and the ratio γ_{0}(τ) = *w*(τ)/*K*_{0} = *N*(τ)/*K*_{0} are shown in dependence of τ. Note that the fixation probability increases with τ: While the mutant does not have a higher intrinsic growth rate than the resident (in a variant of the model, it could even be lower), it thrives due to relaxed competition as the resident population declines. However, for a given mutation probability per time unit *u* from resident to mutant, the new total rate of beneficial mutations is proportional to the resident population size and thus decreases with τ. The rate at which new successful mutations arise is given by θ_{fix}(τ) = *uw*(τ)*p*_{fix}(τ), which reaches a maximum at some intermediate time τ* (see Figure A1). The probability that a successful mutation arises before time *T* is

The probability that a successful mutation arises at all, *i.e.*, that the population is rescued from extinction, is therefore given by(A7)This probability is equal to 1 if and only if the integral diverges, *i.e.*, if the rate θ_{fix}(τ) declines sufficiently slowly. Otherwise the probability that the population is saved from extinction by a rescue mutation is <1, as was also found by Orr and Unckless (2008). For the parameter values used in Figure A1 the survival probability of the population is calculated to be *P*_{rescue} ≈ 0.39. Given the population is rescued, the probability that the rescuing mutation arises in the infinitesimal time interval [τ, τ + *d*τ] is given by *p*(τ)*d*τ with

For not too large mutation probabilities (*u* ≲ 0.5/*K*_{0} for the chosen parameter values), this function has a maximum at an intermediate value of τ = τ**, as also shown in Figure A1. This result is in contrast to Orr and Unckless (2008), who find a monotonic decrease in the (conditioned) probability that a successful mutation arises in a given generation independently of the mutation rate. This is due to the assumption in Orr and Unckless (2008) that the selection coefficient of mutants (and thus the fixation probability *p*_{fix}) of a single mutant is constant over time. In our model, selection pressure on mutants is not simply included by assumption, but results from the underlying ecological dynamics. As such, it explicitly depends on time. The scenario by Orr and Unckless (2008) is approximately met only for large τ > τ** when competition from the resident population gets very small.

- Received February 17, 2011.
- Accepted May 17, 2011.

- Copyright © 2011 by the Genetics Society of America