Understanding evolutionary and ecological dynamics using a continuum limit.

Continuum limits in the form of stochastic differential equations are typically used in theoretical population genetics to account for genetic drift or more generally, inherent randomness of the model. In evolutionary game theory and theoretical ecology, however, this method is used less frequently to study demographic stochasticity. Here, we review the use of continuum limits in ecology and evolution. Starting with an individual-based model, we derive a large population size limit, a (stochastic) differential equation which is called continuum limit. By example of the Wright-Fisher diffusion, we outline how to compute the stationary distribution, the fixation probability of a certain type, and the mean extinction time using the continuum limit. In the context of the logistic growth equation, we approximate the quasi-stationary distribution in a finite population.

INTRODUCTION

The huge computational power available today allows researchers to develop individual‐based models of high complexity to explore dynamical processes in ecology and evolution. Here, we aim to make a link between these individual‐based descriptions and continuous models like (stochastic) differential equations that remain amenable to analysis. We review these techniques and apply them to some frequently used models in ecology and evolution.

In ecology, probably the most common description of population dynamics is the logistic growth equation (Verhulst, 1838). Its attractiveness draws from its simplicity. It has a globally attractive fixed point (when started with any nonzero population size), the carrying capacity of the population. This simplicity comes at the cost that biological observations like population size fluctuations or even extinction events are not captured by this deterministic model. To account for these stochastic effects, one needs to change to a stochastic differential equation which can be derived from individual‐based reactions (Champagnat et al., 2006). Stochastic differential equations consist of a deterministic and a stochastic term, the latter also often referred to as noise. In theoretical ecology, the effect of random environmental fluctuations on deterministic dynamics like the logistic equation has been studied in great detail (e.g., Allen et al., 1993; Schaffer et al., 1986). Demographic stochasticity, which arises from the inherent randomness of birth and death events of individuals, is often modeled in discrete‐time models (e.g., Melbourne & Hastings, 2008; Schreiber et al., 2018). Models including demographic stochasticity in continuous time are more scarce in ecology. Typically, these models are studied in the context of evolutionary ecology, for example to model co‐evolution in ecological communities (Dieckmann & Law, 1996). More generally, questions related to extinction and coexistence, transitions between different stable deterministic equilibria or the maintenance of quasi‐cycles can be studied within the framework of stochastic differential equations (Boettiger, 2018; Jeltsch et al., 2019). Here, we outline how to derive stochastic differential equations from an individual‐based model. We do this along similar lines as reviewed in Black and McKane (2012) but with a larger emphasis on the technical details.

In evolutionary game theory, the Moran process has become a popular model for stochastic dynamics in finite populations (Nowak et al., 2004). It is a model describing the dynamics of different alleles in a population of fixed size and overlapping generations. As this is a birth–death process, quantities like fixation probabilities, fixation or extinction times, and the stationary distribution can be calculated based on recursions (Allen, 2011; Goel & Richter‐Dyn, 1974; Karlin & Taylor, 1975; Traulsen & Hauert, 2009). A continuum approximation for quantities that are known exactly may thus make limited sense at first sight—but it can provide a very useful new perspective. Another important process in population genetics is the Wright–Fisher process—a model for allele frequency dynamics in a population of fixed size and nonoverlapping generations (Wright, 1931). It is more popular in population genetics but is also used in evolutionary game theory (e.g., Imhof & Nowak, 2006; Taylor & Maciejewski, 2012; Traulsen et al., 2006; Wakano & Lehmann, 2014). The Wright–Fisher process is mathematically more challenging to analyze exactly than the Moran process. Therefore, continuum approximations resulting in stochastic differential equations are used to compute typical quantities of interest such as the probability of fixation of a certain genotype or the mean time until this fixation event occurs (Crow & Kimura, 1970; Ewens, 2004). Additional to the derivation of the continuum limit, we also present how to compute these quantities. In evolutionary game theory, demographic stochasticity is particularly important for behaviors that only evolve in small populations, such as spite (Gardner & West, 2004; Nowak et al., 2004). Initially, stochastic replicator dynamics have been developed to address the problem of equilibrium selection in evolutionary games (Cabrales, 2000; Foster & Young, 1990; Fudenberg & Harris, 1992). With the more recent development to look at games with a larger number of possible strategies, an approach based on weak mutation which allows to reduce the analysis to pairwise comparison between strategies became popular (Fudenberg & Imhof, 2006; Wu et al., 2011). These methods are particularly useful when many strategies have to be taken into account (García & Traulsen, 2019; Rand & Nowak, 2011; Sigmund et al., 2010; Traulsen et al., 2012; Vasconcelos et al., 2017). However, while the focus on weak mutation allows an analysis that takes demographic noise into account, additional approaches are necessary when, for example, internal equilibria are present (Black et al., 2012; Vasconcelos et al., 2017). The stochastic tools discussed here may be a further step into that direction.

Even though similar in the questions they try to answer, evolutionary game theory and population genetics are developing in parallel, sometimes with little interaction between them. As this is partly arising from the different methods applied, here we aim to provide an introduction to the continuum limit for those less comfortable with these methods and hesitant to go into the extensive, more mathematical, literature. In summary, for the theoretical population geneticist with a probabilistic background, we provide a summary of some key results on stochastic differential equations; for the evolutionary game theorist, we give a new perspective on the derivations of results obtained when using discrete birth–death processes; and lastly, for the theoretical ecologist familiar with deterministic modeling, we outline how to derive and work with stochastic versions of classical ecological and evolutionary processes.

Since our goal is to illustrate how to apply a continuum limit to individual‐based descriptions of a biological process, the calculations and derivations below may remain vague where more mathematical theory is necessary. For a mathematically rigorous presentation of this topic, we refer to the lecture notes by Etheridge (2012) or the book by Ewens (2004). A more application‐oriented treatment of stochastic processes in biology can be found in the books by Lande et al. (2003), Otto and Day (2007), and Allen (2011).

EVOLUTIONARY AND ECOLOGICAL PROTO‐TYPE PROCESSES

We outline the derivation of continuum limits by application to exemplary processes from evolution and ecology, the Wright–Fisher and Moran process, and the logistic equation. By showing the explicit derivation in these examples, we provide the necessary tool set to derive continuum limits of more complex processes motivated by individual‐based models. In this section, we define the models by their microscopic descriptions, that is, we describe the model dynamics as viewed from an individual's perspective.

Wright–Fisher and Moran process

The two most popular processes to model (stochastic) evolutionary dynamics are the Wright–Fisher and the Moran process. While in the Wright–Fisher process generations are nonoverlapping and time is measured in discrete steps, generations in the Moran model are overlapping and measured in either discrete or continuous time. Originally, both processes describe the stochastic variation of allele frequencies due to finite population size effects referred to as genetic drift.

Wright–Fisher model

One of the oldest population genetics model is the finite size Wright–Fisher process (Fisher, 1930; Wright, 1931). Given a population of constant size, it describes the change in frequencies of alleles in nonoverlapping generations over time, measured in (discrete) generations.

Classically, one considers a population of N individuals where each individual is of type A or B. The population is considered to be in its ecological equilibrium. The population size N is therefore constant over time. One interpretation of the dynamics is that every generation each individual chooses, independently of all other individuals, an ancestor from the previous generation and inherits its type. Under selection, the likelihood of drawing type A individuals increases (or decreases) which introduces a sampling bias. The probability for an offspring to have a parent of type A, conditional on k individuals being of type A in the parental generation, iswhere is the selective advantage of type A. The number of type A individuals in the next generation is then given by a binomial distribution with sample size N and success probability p. Denoting the number of type A individuals in generation n by X we have

Unfortunately, the Wright–Fisher model, even though very illustrative, is difficult to study analytically. Through the developments in stochastic modeling in the last century, a lot of this new theory could be adopted to overcome this problem (e.g., Ewens, 2004; Kimura, 1983).

Moran model

Another way to resolve the difficulties associated with the Wright–Fisher process is provided by the Moran process (Moran, 1958). The setup is the same as for the Wright–Fisher process (constant population size N with two types or – in population genetics – alleles A and B) with one exception: Time is not measured in generations but each change in the population configuration affects only one individual, the one that dies and gets replaced by an offspring of another randomly selected individual. Therefore, generations are overlapping, and time can be measured in discrete steps or continuously.

Discrete time

The Moran process in discrete time progresses as follows. Every time step, one individual is randomly chosen to reproduce and the offspring replaces a randomly chosen individual among the remaining N − 1 individuals (sometimes, the replacement mechanism is not restricted to the remaining individuals but also includes the parent). Therefore, in a population with k type‐A individuals, the probability that one of these replaces a type‐B individual is given bywith p as defined in Equation 1. Analogously, the probability for the number of type‐A individuals to decrease from k to k − 1 reads

We have implemented selection on the reproduction step, but it could also affect the replacement step. In a nonspatial setting, as considered here, this leads to the same transition probabilities. However, the Moran model can also be studied on a graph which aims to model spatial structure. In that case, the order of reproduction and replacement, and which of these steps is affected by selection, matters and can potentially give rise to different evolutionary dynamics (Kaveh et al., 2015; Lieberman et al., 2005). We note further that without selection (s = 0) we have p = k/N, that is, the increase and decrease probabilities are equal for any choice of k. Dynamics with this property are called neutral.

Continuous time

The same dynamics (albeit on a different time scale) are obtained by assuming that each pair of individuals is associated with a random exponentially distributed time (also described as exponential clocks). The next pair to update their types is determined by the smallest random time (or the clock that rings first). At these updating times, one of the two individuals is chosen to reproduce, the offspring replacing the other individual of the pair. There is no standard choice when it comes to choosing the rate of these exponential times.

Both formulations of the Moran process are Markov chains, either in discrete or continuous time, with the special property of having jumps of ±1 only. These processes are called birth–death processes. The theory of these is well developed, see for example the books of Karlin and Taylor (1975, 1981), Gardiner (2004), or Allen (2011), so that the dynamics of Moran processes are often amenable to analysis (typically by solutions of recursion equations).

To sum up the introduction to these two evolutionary models, the difference between the Wright–Fisher model and the Moran model is the progression of populations in time. In the Wright–Fisher process, generations are nonoverlapping, that is, all individuals update their type at the same time. Therefore, the distribution of types in the offspring generation is binomial. In contrast, generations are overlapping in the Moran model and the dynamics are described by a birth–death process since only one individual is updated at a time.

Logistic growth

In ecology, one is typically interested in population sizes or densities rather than allele frequencies. The simplest population growth model is that of exponential growth. Obviously, a population cannot grow exponentially forever. Its growth will be limited at some point, for example due to spatial constraints or resource depletion. This form of density regulation suffices to stabilize a population around its carrying capacity, the positive population size at which in the deterministic process the growth rate equals zero.

Here, we give the mechanistic basis that could potentially describe such a process. We denote a single individual of the population by Y. The birth and death processes can then be written as

The parameters β and δ correspond to the rates at which the two reactions happen, that is, each reaction corresponds to an exponential clock with rate either β or δ. For β > δ, the population grows to infinity (exponential growth), whereas for β <δ, it goes extinct.

Population regulation is achieved through a nonlinear term that is interpreted as an interaction between two individuals, for example, competition for space. The corresponding microscopic process is given by

The parameter γ is referred to as the intraspecific competition coefficient, and K is a measure of the number of individuals at carrying capacity. The division by K in the competition rate is accounting for the probability of interaction of two individuals in a well‐mixed population where space is measured by the parameter K so that Y/K becomes a density (or rate of encountering an individual when randomly moving in space). For a more detailed derivation of these type of interaction rates, we refer to Anderson and Kurtz (2015).

The logistic process is, like the Moran process, a birth–death process. We will see in the next section that in the infinite population size limit (we let K tend to infinity), the mechanistic description above yields the logistic equationwhere r = β − δ is the per‐capita growth rate, c = (β − δ)/γ is the rescaled carrying capacity, and y = Y/K is the density of the population.

INFINITE POPULATION SIZE LIMIT

The microscopic descriptions can be implemented by a stochastic simulation algorithm. Yet, the theoretical analysis of finite size populations can be challenging. A common technique to overcome this challenge is to consider a continuum approximation, that is, studying the limiting model for N (or K) to infinity. The limit is a (stochastic) differential equation of the formwhere (W)  ≥ 0 is a standard Brownian motion. This equation describes the population dynamics, that is, the macroscopic evolution of a certain model. For a general introduction to stochastic differential equations, see for example the books by Karlin and Taylor (1981) and Gard (1988).

The term μ(x) is called the infinitesimal mean, that is, the expected change of the stochastic process (x) ≥ 0 in a very short (infinitesimal) time interval. It represents the deterministic dynamics of the process. The term σ 2 (x) is called the infinitesimal variance, that is, the expected variation of the continuum limit in very small time steps. It quantifies the random fluctuations. In the case where σ 2 (x) is zero, the limit is deterministic and Equation 8 reduces to an ordinary differential equation.

One can formally show that the terms μ(x) and σ 2 (x) indeed correspond to the changes of the mean and the variance in infinitesimally small time steps if they are derived formally (Appendix 1). This allows us to compute them by the following identities:where and denote the expectation and variance of the process x.

A solution of a stochastic differential equation of the form in Equation 8 is called a diffusion. Another common representation of Equation 8 is the following integral equationwhere the stochastic integral is interpreted in the sense of Itô. For a discussion of the different choices of stochastic integrals and their consequences in terms of modeling, we refer for example to Turelli (1977).

We now present how to derive Equation 8 for the three introduced models. The strategy is rather simple: Compute the infinitesimal mean and variance as given in Equation 9 for the individual‐based model.

Discrete‐time derivation: Wright–Fisher model

For the reason of illustration, we assume a Wright–Fisher model as outlined in Section 2.1 without selection, s = 0. We need to compute the expectations in Equation 9 using the probability distribution given in Equation 2 (with s = 0). Let us ignore the time step Δt for the moment and simply compute the change in expectation and variance from one generation to the other. Writing X for the number of individuals of type A and setting ΔX = X − X −1, we findwhere we used that the number of individuals of a certain type in the next generation is binomially distributed. Analogously, the infinitesimal variance is

It remains to account for the transition from discrete to continuous time. In this case, the natural choice is Δt = N −1. This can be seen by examining the infinitesimal variance. To obtain a limit different from zero or infinity in Equation 12 for N to infinity, we need to divide that equation by N. Comparing Equation 12 to the corresponding line in Equation 9, we set Δt = N −1. Replacing k/N by x and taking the limit N → ∞, we find the neutral Wright–Fisher diffusion for the allele frequency dynamics

This equation is called Wright–Fisher diffusion and describes the dynamics of a neutral allele due to genetic drift. In other words, the allele frequency behaves like a random walk in continuous time and space. A similar derivation as above can be done by including selection and mutation. The more lengthy computation is relegated to Appendix 2.

Conclusion 1: To derive a continuum limit of a finite population size model in discrete time, one computes the infinitesimal mean and variance as given in Equation 9 and rescales time so that the two quantities converge in a meaningful way, that is, do not tend to infinity.

Continuous‐time derivation: general case

In principle, the same methodology as above is applicable for the derivation of the continuum limit of a process measured in continuous time. However, in view of our subsequent analysis of the limit Equation 8, we will introduce a new tool, the infinitesimal generator.

The change in infinitesimal time of any continuous‐time Markov process (x) ≥ 0 can be described by the infinitesimal generator, denoted (Ethier & Kurtz, 1986, chapter 1, eq. (1.10)). Intuitively, one can think of it as the derivative of the expectation (of an arbitrary function) of a stochastic process. Formally, it is defined bywhere denotes the conditional expectation of the stochastic process f(x) at time Δt given the initial value x 0 = x. Here, f is an arbitrary function so that the limit is well‐defined. For example, applying to f(x) = x describes the dynamics of the mean of x, and for f(x) = x 2, we obtain the dynamics of the second moment of x. From the first two moments, we can recover the variance, so that from Equation 14, one can derive the infinitesimal mean and variance.

The infinitesimal generator is useful in our context since it can be related to a diffusion process. More precisely, the infinitesimal generator associated with the stochastic diffusionis given by (we refer to Appendix 1 for a derivation by the Itô formula)

Our strategy is to find a limit of the infinitesimal generator associated with a finite population size process, which corresponds to the form given in Equation 16. We consider a continuous‐time birth–death process with transition rates T + and T − for 0 ≤ k ≤ N. Due to the exponentially distributed waiting times, the probability for a single update until time t is λt exp(−λt), where λ is the rate of the corresponding exponential clock. Setting x = X/N, the frequency of type‐A individuals, we find the infinitesimal generator for the model with finite population size N, , to be of the form

We have used the Landau notation o(Δt) to summarize processes that scale with order (Δt)1 +  for any ε > 0. We will also use the Big‐O notation O(Δt) for processes that scale with order Δt or higher. Doing a Taylor expansion for large N and neglecting the terms of order higher than 1/N 2, we find

Translating this equation to a stochastic differential equation, we identify the single components as

Note that we have made no assumption on the dependence of the transition probabilities on the frequency x, such that this approach is applicable for constant selection, linear frequency dependence arising in two player games (Traulsen et al., 2005) or multiplayer games with polynomial frequency dependence (Gokhale & Traulsen, 2010; Peña et al., 2014).

Conclusion 2: For time‐continuous finite population size models with jumps of ±1, that is, a birth–death process, the terms of the continuum limit can be computed by Equation 19.

Example: Moran process with selection and mutation

Returning to our proto‐type processes, we explicitly derive the stochastic differential equation corresponding to the Moran model with selection and mutation. We decouple reproduction and mutation processes, but similar derivations can be made if we assume a coupling of mutations to reproduction events. The selection coefficient is denoted by s, and the mutation rates from type A to B and type B to A are given by u → and u → , respectively. Then, the transition rates are

Inserting these into Equation 19 yields

Depending on the choice of selection and mutation rates, these equations result in different limits. Typically, one is interested in nontrivial limits for these equations, that is, a limit so that not both components equal zero. Often this can be achieved by rescaling time (Section 3.3) and/or defining the strength of selection and mutation in terms of the population size N. As an example, we will focus on two specific limits: (a) strong selection and strong mutation and (b) weak selection and weak mutation.

Strong selection and mutation

We consider large selection and mutation rates. We assume that s and u do not depend on N but are constant. To obtain a limit equation for the frequency of individuals of type A, x = X/N, we rescale time by N, that is, t→Nt, which transforms Equation 21 to

The first equation is independent of N and the vanishing variance in the second equation implies that the limit process is deterministic. We find the ordinary differential equationwhich describes the change in allele frequency in a population under strong selection and mutation over time.

Weak selection and mutation

In contrast to the previous scenario, we now assume that both selection and mutation are weak. We let s and u scale inversely with N and define the constants α = sN and ν = u. Inserting these into Equation 21 (and here without rescaling time) yieldswhich gives the diffusion limit

Note that compared to the Wright–Fisher diffusion in Equation 13 this limit has twice as much variance. This difference is explained by the different sampling schemes in the individual‐based description of the model. To see this, we assume no selection, s = 0, and no mutation, u →  = u →  = 0.

In the Wright–Fisher process, individuals are updated by binomial sampling. The variance of this sampling procedure is Nx(1 − x), where the factor N vanishes by rescaling the time. This gives σ 2 (x) = x(1 − x).

In the Moran model, or more generally for a birth–death process, the variance is computed by the sum of the transition rates, cf. Equation 19. In our example, both transitions happen at rate 1, which explains the additional factor 2. The difference between the variances is therefore a result of the different sampling schemes of the individual‐based models.

As a consequence of this difference in the variance σ 2 (x) between the two models, the Moran diffusion limit progresses twice as fast as the Wright–Fisher diffusion limit which can be seen by the scaling property of the Brownian motion. In terms of the original discrete processes, this means that N individual jumps, like in the Moran process, accumulate more variance than one update of the whole population, like in the Wright–Fisher process. The sampling therefore determines the variance and consequently the speed of the continuum limit.

Conclusion 3: The Moran process, by definition, has the same mean behavior as the Wright–Fisher model. However, its variance in the diffusion limit is twice the variance of the corresponding Wright–Fisher diffusion. This difference arises from the different sampling schemes of the individual‐based models.

Change of time scales in the derivation of a continuum limit

In the derivation of the continuum limit, we have repeatedly rescaled time to obtain a nontrivial limit, for example, right before Equations (13), (22) and 22. Rescaling the time speeds up (or slows down) the original process so that the dynamics of interest, for example, allele frequency changes, become observable. For example, if the dynamics were to be very fast in the original process, we would need to slow down time appropriately to observe the changes of the quantity of interest more gradually. In general, we are free to chose any scaling of time. However, it is important to keep in mind the scaling when interpreting results obtained in the limiting process in terms of the original process. Especially so, if one is interested in quantities involving time, for example, fixation or extinction times.

Conclusion 4: Different assumptions on the model dynamics, for example, on selection and mutation, can lead to different continuum limits on the population level. To identify parameter combinations that result in a reasonable continuum limit, one needs to study Equation 19 to match the orders of the scaling parameter. Rescaling time gives an additional degree of freedom when trying to match these orders to obtain a reasonable limit.

DIFFUSION APPROXIMATION

We have seen that if we let the population size N tend to infinity, we can derive a (stochastic) differential equation describing the studied evolutionary or ecological process. A natural question that arises is how these results relate to finite population size models. To study this difference between the finite population process and the continuum limit, we consider the logistic growth equation. The transition rates are given by

Repeating the steps from the previous section with y = j/K, we find the following expressions for the infinitesimal mean and variance:

Thus, the classical deterministic logistic equation is obtained in the infinite population size limit, K → ∞:

How well does the finite population size description approximate this deterministic limit? One way to approach this question is to simply not take the limit of K to infinity. The finite population size logistic equation, derived from Equation 27, is then approximated bywhere the superscript K in indicates the order of magnitude of the carrying capacity.

The approximation in Equation 29 is called Diffusion approximation. One can prove formally that this approximation, under the assumptions that the function μ(x) and σ 2(x) are (twice) continuously differentiable, performs equally well as a more accurate analysis based on the central limit theorem (Ethier & Kurtz, 1986, Theorem 11.3.2). For a rigorous discussion of diffusion approximations and their relation to the central limit theorem, we refer to Ethier and Kurtz (1986, Chapter 11).

In terms of performance of the diffusion approximation, the population size measure K does not need to be very large for the individual‐based model to approach the deterministic limit (K ≈ 1,000 is enough in this example, Figure 1).

FIGURE 1

Individual‐based simulations of the logistic growth model. (a) For low population sizes, the individual‐based simulation (solid lines) fluctuates strongly around the deterministic solution of the population (dashed lines) given by Equation 28. (b) Increasing the scaling parameter K, the stochastic fluctuations around the deterministic prediction decrease, until eventually the individual‐based simulation is indistinguishable from the deterministic curve. The parameter values are chosen as follows: β = 2, δ = 1, γ = 1, and (a) K=100, (b) K=1,000. The initial population sizes are stated in subfigure (b)

Conclusion 5: Not taking the limit in Equation 19 yields the diffusion approximation of the studied model. This approximation is a stochastic differential equation (Equation 29) where the variance (typically) scales inversely with the square root of the scaling parameter.

STATIONARY DISTRIBUTIONS

For the Moran model, we have derived two different limits that differ in their assumptions on selection and mutation. If both selection and mutation are strong, the infinite population size limit is an ordinary differential equation. For weak selection and weak mutation, we derived a stochastic differential equation. One qualitative difference between these two limits is that trajectories of the deterministic limit will always converge to a fixed point (other limits are possible in general, e.g., limit cycles) while the stochastic differential equation fluctuates indefinitely for positive mutation rates. The deterministic fixed point of the Moran model is given by the solution of Equation 23 equal to zero. In our example, a single fixed point x* lies within the interval between 0 and 1 and is stable. Therefore, all trajectories will approach this value, for example, Figure 2a.

FIGURE 2

Allele frequency dynamics with selection and mutation. (a) The deterministic system given by Equation 23 converges to the fixed point (dashed line) and remains there. (b) The stochastic process given by Equation 25 fluctuates strongly in frequency space and for the chosen parameter values spends most time close to the monotypic states x = 0 and x = 1

In contrast, Equation 25 is a stochastic equation. Thus, even if the trajectories approach or hit the deterministic fixed point they will not stay there due to the stochasticity of the Brownian motion, cf. Figure 2b. Still, we can make predictions about the time a trajectory spends in certain allele configurations. This information is summarized in the stationary distribution, the stochastic equivalent of a deterministic fixed point. If the initial state of the population is given by the stationary distribution, then the distribution of all future time points will not change. For birth–death processes, the stationary distribution can be calculated based on detailed balance, that is, the incoming and outgoing rates need to be equal for every state of the process (Antal et al., 2009; Claussen & Traulsen, 2005; Gardiner, 2004). Formally, the stationary distribution, denoted ψ, is defined as the solution ofwhere x 0∼ψ denotes that x 0 is distributed according to ψ and f is an arbitrary function. Importantly, this condition means that the distribution of allele frequencies does not change over time because its derivative in time is zero (for any choice of f). The above equation can be solved with the infinitesimal generator (Etheridge (2012, Chapter 3.6)). The solution is expressed in terms of the speed measure density m(x), which we introduce in Box 1, and is given by

Intuitively, the speed measure at a point x, m(x), quantifies the time which the process spends in this state. Therefore, ψ(x) is nothing but the average time spent in state x.

Conclusion 6: The stationary distribution of a one‐dimensional diffusion can be expressed in terms of the density of the speed measure m(x) through Equation 31. The density of the speed measure is given by the scale function corresponding to the stochastic diffusion process, Equations (32), (34) and (22), (34) in Box 1.

A one‐dimensional stochastic diffusion can be transformed into a standard Brownian motion. Since the Brownian motion is well‐studied, a lot of results can then be translated to the stochastic diffusion by the transformation functions, the scale function, and the speed measure.

First, we rescale the space of the original process by the scale function. It is defined by where the lower boundaries of the integrals can be chosen arbitrarily. The name of this function derives from the fact that for a one‐dimensional diffusion x satisfying

the scaled process S(x) becomes a time‐changed Brownian motion on the interval [S(0), S(1)], that is, there is no deterministic contribution in the scaled process. The process S(x) is a Brownian motion with a “non‐standard” time scale. To map this time‐changed Brownian motion to the time scale of a standard Brownian motion, one needs to rescale time by the speed measure M. It defines how much faster (or slower) the process S(x) is evolving compared to a standard Brownian motion. The speed measure is defined by the density of the speed measure. The time is then rescaled by

Compactly written, we have changed the stochastic diffusion x to the standard Brownian motion by the following steps:

Stationary distribution of the Wright–Fisher diffusion

As an example, let us consider the Wright–Fisher diffusion with selection and mutation as derived in Equation 13 (and Equation 25 when derived from the Moran model), that is,

Computing Equation 31 with help of the quantities defined in Box 1, one obtainswhere is the Gamma function and 1 F 1 (a, b, z) is the generalized hypergeometric function.

The equation itself does not provide much insight. To illustrate the possible shapes of stationary distributions, we plot several choices of mutation rates and selection coefficients in Figure 3. We see that for higher mutation rates, more probability mass is allocated to intermediate allele frequencies (compare the solid and dashed lines). In this case, the Wright–Fisher diffusion spends more time in states of coexistence than in monotypic states (the boundaries of the allele frequency space in Figure 3) because temporary extinction events are prevented by recurrent mutations. If the mutation rates are asymmetric (dotted line), the stationary distribution is skewed toward the type with the lower mutation rate. If one type is favored selectively (dash‐dotted line), the stationary distribution is skewed toward the favored type.

FIGURE 3

Stationary distribution of the Wright–Fisher diffusion with selection and mutation. The lines are given by Equation 36. Larger mutation rates accumulate more probability on intermediate allele frequencies (compare solid and dashed lines). Selection (or asymmetric mutation) skews the stationary distribution toward the selectively favored type (or type with the lower mutation rate), dash‐dotted (dotted) line

Stationary distributions are the stochastic equivalents of deterministic fixed points and as such provide a basic description and a starting point for further analysis of the qualitative behavior of a stochastic model, especially in situations where polymorphisms of alleles, coexistence of species, or spatial population distributions are to be expected (e.g., Czuppon & Rogers, 2019; Gaston & He, 2002; Lehmann, 2012; Polansky, 1979; Turelli, 1981).

Quasi‐stationary distribution of the logistic process

Next, we consider the diffusion approximation of the logistic growth model, that is,

The logistic process for finite K has a (unique) absorbing state, y = 0 because there is no transition from this state to positive population densities. Once the population went extinct, it remains so. Since the extinction state is accessible from all values y > 0 (and because the population size remains finite for all times), the population will go extinct with probability 1. The only stationary distribution is the point measure on 0, that is, ψ ∼ δ 0 (δ is the Dirac measure at point x).

In contrast, the positive deterministic population equilibrium, y* = (β − δ)/γ, is a stable fixed point of the deterministic system. Considering large values of the deterministic equilibrium (), we expect the finite population size process from Equation 37 to remain close to this value for long times. In fact, the expected extinction time of the logistic growth process when started in the positive population equilibrium is of order exp(K) (Champagnat, 2006). This suggests that the process will be in a quasi‐stationary state, that is, before its extinction the population is described by the stationary distribution of the corresponding logistic process conditioned on nonextinction.

Formally, the quasi‐stationary distribution is computed by conditioning the original process on its survival. This means that the transition rates change and the novel process can be analyzed by the techniques described above. However, this method goes beyond the scope of this manuscript. For a theoretical treatment of this topic in the context of the logistic equation, we refer to Cattiaux et al. (2009), Assaf et al. (2010), Méléard and Villemonais (2012). For a general review on methods related to quasi‐stationary distributions, see Ovaskainen and Meerson (2010).

Another way to approximate the quasi‐stationary distribution when extinction is very unlikely for long times (which is the case for large K) is provided by the central limit theorem (sometimes also called linear noise approximation in this context). Here, the distribution of the process is derived from its local dynamics around the deterministic fixed point y* (Ethier & Kurtz, 1986; van Kampen, 2007). The underlying assumption is that the population stays close to its positive steady state and just slightly fluctuates around this value. This is only a valid assumption when the probability of extinction within the studied time frame is essentially zero. These small fluctuations are described by a Gaussian distribution. Formally this translates towhere U is a Gaussian random variable and y the deterministic trajectory. Writing μ(y) = (β − δ − γy)y and σ 2 (y) = (β + δ + γy)y, the dynamics of U can be rewritten as

Evaluating the process U at , we obtain a description of the variance in the deterministic fixed point. For this fixed choice of y and , U becomes an Ornstein–Uhlenbeck process. Its stationary distribution is then given by

This distribution describes the fluctuations of the process around the deterministic steady state y*. Therefore, when plugging the distribution ψ into the original process from Equation 38, we find the quasi‐stationary distribution of around the deterministic equilibrium y* which yields

For increasing population sizes K, the variance is decreasing and vanishes in the limit K→∞, as is to be expected by the deterministic limit.

Conclusion 7: If the deterministic process has a stable steady state but is almost surely going extinct for finite population sizes, a quasi‐stationary distribution can be computed to describe the behavior of the process conditioned on survival. If the extinction probability is very low, an approximation of this distribution is given by the linear noise approximation where the variance around the deterministic steady state is modeled by the Ornstein–Uhlenbeck process derived from Equation 39.

FIXATION PROBABILITIES

We have seen that stochastic descriptions of processes can lead to outcomes that are different from their deterministic counterparts. Here, we study one of these phenomena: the probability for a certain type to become fixed in a population. For one‐dimensional stochastic differential equations, fixation probabilities can be computed explicitly. As before, we denote by x the frequency of type A individuals at time t ≥ 0 in the population. We use the fact that the mean of the scaled process S(x) in Equation 32 does not change over time. Stochastic processes with this property are called martingales. With we findwhere we have used the martingale property of the scaled process S(x) in the first equality. The second equality is explained by the process being absorbed at one of the two boundaries x = 0 or x = 1 at large times. For a formal derivation, we refer to Otto and Day (2007, Chapter 15.3.3) or Etheridge (2012, Lemma 3.14).

As an example, let us consider the Wright–Fisher diffusion with selection and without mutations (ν →  = ν →  = 0) given in Equation 35. We have μ(x) = αx(1 − x) and σ 2(x) = x(1 − x) such that the scale function simplifies to

Recalling the definition of α = sN for a finite population of size N and plugging this into Equation 42 yieldswhich for x 0 = 1/N becomes , the result of Haldane for the fixation of a single mutant copy in a population of size N (Haldane, 1927). The first line of Equation 44, the classical result of fixation probabilities when derived from diffusion theory, and its applicability has been the subject of extensive research, for example, Bürger and Ewens (1995) and references therein; for a more general review on fixation probabilities, we refer to Patwa and Wahl (2008).

Of course, the fixation probability can also be calculated for more complicated stochastic differential equations where the sign of the deterministic dynamics μ(x) depends on the population configuration. Most classically, these frequency‐dependent problems were studied in deterministic evolutionary game theory introduced by Maynard Smith and Price (1973) (see also Hofbauer and Sigmund (1998) for an introduction to evolutionary game dynamics). In Appendix 3, we (re‐)derive the fixation probability in case of frequency‐dependent selection.

Conclusion 8: The fixation probability of a one‐dimensional diffusion is given by the scale function as stated in Equation 42.

MEAN TIME TO FIXATION

A related quantity of interest is the expected time to fixation (or extinction from the other type's point of view), that is, the average time of coexistence of two types. Again, the calculation relies on a special function, this time Green's function G(x, y), which can be interpreted as the average time that a diffusion started in x spends in the interval [y, y + dy) before reaching one of the boundaries (Etheridge, 2012, Chapter 3.5). It is therefore also called sojourn time density (Ewens, 2004). It is defined aswhere S(x) is the previously defined scale function and m(x) denotes the speed measure density (Box 1).

The expected time to fixation for a process started at frequency x, denoted , is then given by (Ewens, 2004, section 4.4).

This integral corresponds to the summation of the sojourn times in the discrete case, for example, we refer to Ohtsuki et al. (2007) for an application in finite populations. In some cases, the result of this equation yields an analytically tractable result, for example, for the neutral Wright–Fisher diffusion

In this case, the scale function and speed measure density are given by

Then, the expected time to fixation of one of the two alleles can be expressed as

In Appendix 3, we consider the more involved example of frequency‐dependent selection (Altrock & Traulsen, 2009; Pfaffelhuber & Wakolbinger, 2018).

Similarly to fixation probabilities, the mean time to fixation has been studied extensively through stochastic diffusions, for example, Kimura and Ohta (1969). It is especially important in population genetics where one is interested in the time to extinction or fixation of newly arising alleles (van Herwaarden & van der Wal, 2002). In ecology, the mean time to extinction or fixation is for example applied in the context of population extinction (Lande, 1994) and speciation events (Yamaguchi & Iwasa, 2013).

Conclusion 9: Expected unconditional fixation times, that is, the expected time of coexistence of two types in a population, can be calculated by integrating over Green's function (the mean occupation time of a certain frequency until extinction), as shown in Equation 46.

DISCUSSION AND CONCLUSION

We have outlined how to derive a stochastic differential equation from an individual‐based description of two classical models in evolutionary theory and theoretical ecology, the Wright–Fisher diffusion and the logistic growth equation. The resulting stochastic differential equations in one dimension describe the dynamics of the allele frequency and population density, respectively. Using probabilistic properties of this equation, that is, transforming it to a standard Brownian motion (Box 1), it is possible to analytically derive the (quasi‐) stationary distribution, fixation probability, and the mean time to fixation. As an example, we derived these quantities for the Wright–Fisher diffusion.

The diffusion process emerges as the infinite population size limit. However, as we have shown in Section 4, one can also derive a finite population size approximation of the dynamics, the diffusion approximation. The fixation probability, mean extinction time, and stationary distribution are accessible by the same means as for the continuum limit. Applications of diffusion approximations are abundant and cover diverse topics (e.g., Assaf & Mobilia, 2011; Constable et al., 2016; Czuppon & Gokhale, 2018; Czuppon & Traulsen, 2018; Débarre & Otto, 2016; Houchmandzadeh, 2015; Kang & Park, 2017; Koopmann et al., 2017; McLeod & Day, 2019; Parsons et al., 2018; Reichenbach et al., 2007; Schenk et al., 2020; Serrao & Täuber, 2017; Traulsen et al., 2005).

Apart from the fixation probability and the mean time to fixation, the (quasi‐)stationary distribution is a commonly used measure to describe stochastic processes. Its calculation through the speed measure of the associated scaled process (Box 1) is (in many cases) numerically straightforward. If the process has an absorbing state, for example, an extinction boundary of the population, the stationary distribution is not meaningful. Here, the quasi‐stationary distribution describes the stationary distribution conditioned on the survival of the population. For negligible extinction probabilities, that is, very large survival probabilities of the population, the functional central limit theorem (or linear noise approximation) can be used to approximate this quasi‐stationary distribution. In the theoretical biology literature, this method is frequently used in models of gene regulatory networks (see Anderson and Kurtz (2015) for a mathematical introduction), and less so in the context of ecology or evolution (e.g., Boettiger et al. (2010); Kopp et al. (2018); Wienand et al. (2018); Czuppon and Constable (2019); and Assaf and Meerson (2017) for a review of the physics literature related to this topic).

Lastly, we did not cover multi‐dimensional or spatially explicit stochastic differential equations in this methods review. These processes are often much more complicated to analyze. Here, we aimed to give a brief introduction into the derivation of a continuum limit from an individual‐based model. We hope, that with our basic comparisons between different approaches used in different subfields of theoretical biology, we help to clarify the common foundation, the individual‐based model, on which these different methods are based.

CONFLICT OF INTEREST

The authors declare no conflict of interest.

AUTHOR CONTRIBUTION

Peter Czuppon: Conceptualization (equal); methodology (lead); writing–original draft (lead); writing–review and editing (equal). Arne Traulsen: Conceptualization (equal); writing–review and editing (equal).