## Abstract

We analyze the dynamics of two alternative alleles in a simple model of a population that allows for large family sizes in the distribution of offspring number. This population model was first introduced by Eldon and Wakeley, who described the backward-time genealogical relationships among sampled individuals, assuming neutrality. We study the corresponding forward-time dynamics of allele frequencies, with or without selection. We derive a continuum approximation, analogous to Kimura’s diffusion approximation, and we describe three distinct regimes of behavior that correspond to distinct regimes in the coalescent processes of Eldon and Wakeley. We demonstrate that the effect of selection is strongly amplified in the Eldon–Wakeley model, compared to the Wright–Fisher model with the same variance effective population size. Remarkably, an advantageous allele can even be guaranteed to fix in the Eldon–Wakeley model, despite the presence of genetic drift. We compute the selection coefficient required for such behavior in populations of Pacific oysters, based on estimates of their family sizes. Our analysis underscores that populations with the same effective population size may nevertheless experience radically different forms of genetic drift, depending on the reproductive mechanism, with significant consequences for the resulting allele dynamics.

WHEREAS natural selection undoubtedly plays an important role in adaptation, genetic drift is recognized as an equally or even more important force in shaping the patterns of heritable variation in a population. At a qualitative level, genetic drift implies that variation will tend to be lost from a population over time, even in the absence of selection, and that, in the presence of selection, the fitter type in a population is not guaranteed to fix. Our understanding of genetic drift and its interplay with natural selection has been shaped primarily by the population models introduced by Wright and Fisher (Wright 1931; Fisher 1958) and Moran (Moran 1958). Most work in population genetics rests implicitly on the Wright–Fisher framework, including Kimura’s work on fixation probabilities (Kimura 1955), Ewens’ sampling formula (Ewens 1972; Lessard 2007), Kingman’s coalescent (Kingman 1982), tests of neutrality (Hudson *et al.* 1987; Tajima 1989; McDonald and Kreitman 1991; Fu and Li 1992; Fay and Wu 2000), and techniques for inferring mutation rates and selection pressures (Sawyer and Hartl 1992; Yang and Bielawski 2000; Bustamante *et al.* 2001; Desai and Plotkin 2008).

The behavior of the Wright–Fisher model is remarkably robust to modifications of the model’s underlying assumptions. Indeed, many other population genetic models, including the Moran process (Moran 1958), and a number of Karlin–Taylor and Cannings processes (Karlin and McGregor 1964; Cannings 1974; Ewens 2004), share the same diffusion limit as the Wright–Fisher model (Möhle 2001) and therefore have very similar behavior in sufficiently large populations. Likewise, the Kingman coalescent—that is, the backward-time description of genealogical relationships among a sample of individuals drawn from a Wright–Fisher population—is also robust to deviations in the underlying model assumptions (Kingman 1982).

Nevertheless, the Wright–Fisher model and its diffusion limit are not appropriate in all circumstances. As a result, in some cases the standard model of genetic drift does not provide an accurate description of variation in a population or accurately predict how drift disrupts the force of natural selection. The Wright–Fisher model is founded on a set of assumptions about how organisms reproduce and the way in which the constraint of finite population size enforces dependencies among individuals’ offspring. In particular, the Wright–Fisher model assumes that individuals each produce a Poisson-distributed number of offspring each generation, subject to the constraint of a constant population size (Karlin and McGregor 1964). This formulation excludes the possibility of a highly skewed distribution of offspring numbers, which has been observed empirically in some species. Marine species in particular, as well as some plants and fungi, sometimes produce a very large number of offspring when faced with high mortality early in life (Hedgecock 1994).

To accommodate species with the possibility of large family sizes, researchers have developed a form of coalescent theory that describes the genealogical relationships among samples from such a population, assuming neutrality (Pitman 1999; Sagitov 1999; Schweinsberg 2003; Eldon and Wakeley 2006, 2008, 2009; Sargsyan and Wakeley 2008). These generalized coalescents can differ dramatically from the standard Kingman coalescent, and they include the possibility of multiple mergers in the genealogy. As a result, for species with skewed offspring numbers the inferences reached by assuming the standard coalescent are likely to be incorrect (Eldon and Wakeley 2006). The forward-time dynamics of alleles in such populations have been defined in an abstract setting (Donnelly and Kurtz 1999; Birkner and Blath 2009), but their population-genetic consequences have yet to be worked out, especially those concerning the interaction between selection and genetic drift.

In this article, we study the forward-time population genetics of the population model introduced in Eldon and Wakeley (2006), and we describe their relationship to coalescent theory. Even though the backward-time perspective of such a population facilitates inferences from sampled genetic data, there are many merits to developing a corresponding forward-time perspective for such populations. Most importantly, a forward-time analysis allows us to incorporate the effects of natural selection, whose inclusion is difficult in a backward-time theory. Remarkably, we demonstrate that selection operates very differently in the Eldon–Wakeley model than it does in the standard Wright–Fisher or Moran model. In particular, the form of genetic drift that arises in populations with reproductive skew tends to amplify the effects of selection, relative to the standard form of Wright–Fisherian drift—even when both models are normalized to have the same variance-effective population size. The extent of selection amplification in the Eldon–Wakeley model can be so dramatic that advantageous alleles may even be guaranteed fixation in a population, despite the presence of genetic drift.

The Eldon–Wakeley model is an example of a generalized Wright–Fisher (GWF) process, a class of population processes first introduced and studied in Der *et al.* (2011). While the general methods and bounds of that earlier work apply to the Eldon–Wakeley process, the presentation in this article is considerably more refined and specific to the case of a skewed offspring distribution. In particular, we prove the convergence of the discrete chain to a particular Lambda population process, derive the appropriate scaling conditions on selection necessary for the existence of the continuum limit, and present a new analytical formula to quantify the amount of selection amplification exhibited by the Eldon–Wakeley model, which may be generalized to the entire class of Lambda processes (see *Discussion*).

Our presentation is organized as follows. We start by describing the population model introduced by Eldon and Wakeley, and we specify it precisely in terms of a Markov transition matrix for the change in allele frequency between subsequent time steps. We next describe how to introduce selection into the model, by modifying the Markov transition matrix in a standard way. To analyze the forward-time behavior of the Eldon–Wakeley model, either with or without selection, we develop a continuum approximation, analogous to Kimura’s diffusion approximation of the Wright–Fisher model. We apply this continuum approximation to study the probability that a beneficial mutant will fix in the Eldon–Wakeley model, in three distinct regimes. We end by discussing how our results relate to empirical populations and to coalescent theory.

## The Eldon–Wakeley Model

Eldon and Wakeley (2006) introduced a simple model of a population, based on the Moran process, that incorporates the possibility of large family sizes in the offspring distribution. They considered a population of *N* individuals, each of two types. At each time step, a single individual is chosen uniformly from the population and produces a random number *U* of offspring, drawn from a distribution of offspring numbers, *P _{U}*. The subsequent generation is then composed of the

*U*offspring from the chosen individual supplemented by

*N*−

*U*other individuals, randomly selected without replacement from the remainder of the population. To be clear, only a single individual contributes offspring in each reproduction event—the remaining individuals who neither contribute offspring nor die simply persist to the next time step. A schematic of this reproduction mechanism is shown in Figure 1.

This model coincides with the classical Moran process when the distribution *P _{U}* is concentrated at two individuals. Eldon and Wakeley (2006) considered the following distribution of offspring numbers,

*λ*≤ 1 and

*γ*> 0 are parameters of the model. The Eldon–Wakeley process is thus a mixture of two mechanisms: a typical Moran step in which one individual replicates and displaces another, which occurs with probability 1 −

*N*

^{−}

*, and a rarer event in which a fixed fraction*

^{γ}*λ*of the population is replaced by the offspring from a single individual, which occurs with probability

*N*

^{−}

*. The rare large-family-size events can alternatively be thought of as rapid bottlenecks, in which a fraction*

^{γ}*λ*of the population is removed and immediately replaced by the offspring of a single individual. The parameter

*γ*controls the frequency of such “bottleneck” events.

Our analysis in this article focuses on the Eldon–Wakeley offspring distribution given by Equation 1, but many of our results can be generalized to a large class of offspring distributions *P _{U}* (see

*Discussion*).

Eldon and Wakeley (2006) analyzed the genealogical relationships in a sample from the population, assuming neutrality. We introduce selection into the Eldon–Wakeley model and study the forward-time behavior of the population. We are especially concerned with understanding the fixation probability of an adaptive mutant and how this probability differs from that of the standard Moran model.

### The forward-time transition matrix

Here we study the forward-time evolution of a population under the Eldon–Wakeley model with selection, in the absence of mutation. We keep track of the number of individuals in generation *k* of type 1, denoted by *X _{k}*. The process

*X*is a Markov chain on the states {0, … ,

_{k}*N*}, with some transition matrix

**P**

*, 0 ≤*

_{ij}*i*,

*j*≤

*N*. We decompose the transition matrix

**P**into a product of two matrices: one transition matrix due to selection alone and another due to the (neutral) reproduction process alone. In other words, we introduce selection by first altering the relative frequencies of the two types and then proceeding as in the neutral Eldon–Wakeley case. This form of viability selection coincides with the standard form of selection in the Moran model. (Other forms of selection, such as fecundity selection, generally coincide as well in the limit of large population size; see

*Discussion*.)

More specifically, we define the transition matrix **P** as**S** can be thought as the transition distribution for the standard logistic growth process, which, originating in state *i*, moves to state [(1 + *s*/*N*)*i* ⋅ *N*/((1 + *s*/*N*)*i* + (*N* − *i*))] with probability 1. We use the notation [*x*] to denote the integer part of the real number *x*. Under neutrality, **S** = **I**, the identity matrix. The matrix **Q** describes the neutral genetic drift of the Eldon–Wakeley model; it may be written as the mixture**M** is the standard transition matrix for the two-type neutral Moran process, and **J** is a “jump” matrix modeling the rare bottleneck events. The rows of **J** are in general mixtures of hypergeometric distributions (see *Appendix*); to a high degree of approximation, **J** is a matrix of zeroes except for the index pairs (*i*, *i* + [*λ*(*N* −*i*)]) and (*i*, *i* − [*λi*]), whose entries are respectively *i*/*N* and 1 − *i*/*N*, denoting jumps in the allele frequency caused by the bottleneck events.

The “offspring distribution” of the process is the distribution represented by the row **Q**_{1,}* _{j}*, for

*j*= 0, 1, 2, … . Its variance, given by (4), is an important quantity that will determine the appropriate time scaling of the continuum limit. For large

*N*the offspring variance is

From the above, *N* arising from the Moran part of the process and a second term arising from the bottleneck component. Note that when *γ* < 2, *i.e.*, when the bottleneck events are sufficiently common, the offspring variance is dominated by the contribution from the jump process, whereas for *γ* > 2, it is dominated by the Moran component. This transition point in the offspring variance has important consequences for the nature of the continuum limit, discussed below.

## A Continuum Approximation of the Eldon–Wakeley Process

Analysis of the Moran (or Wright–Fisher) model is typically facilitated by a continuum limit, which becomes accurate in the limit of large population size *N* → ∞ (Kimura 1955; Ewens 2004). For the Moran model, the continuum limit is a diffusion—*i.e.*, it has continuous sample paths and it can be described by the second-order partial differential equations introduced by Kimura and others. By contrast, the continuum approximation of the Eldon–Wakeley model is complicated by the fact that no diffusion process (with continuous sample paths) is an appropriate limit point, due to the occasional large family sizes. These extreme reproduction events produce jumps in the allele frequencies, and hence discontinuous sample paths appear in the limit. Nonetheless, the techniques described in Der (2010; Der *et al.* 2011) can be applied to the Eldon–Wakeley model, with or without selection, to derive a continuous-time/space-limiting process that accurately approximates the discrete Eldon–Wakeley process.

### State, time, and selection scalings

To derive a continuum approximation of a discrete process *X _{k}* we must choose how to scale time and space, as

*N*→ ∞. Such scalings involve replacing the number

*i*of individuals with the frequency

*x*=

*i*/

*b*and the generation number

_{N}*k*by the time

*t*=

*kc*, for some choice of sequences {

_{N}*b*}, {

_{N}*c*}. The continuum limit is then the process

_{N}*b*=

_{N}*N*and

*c*=

_{N}*N*

^{−2}and the resulting continuum limit

*et al.*2011) that the choice of state-space scaling dictates the choice of time scaling—up to a constant factor—in the sense that no other scaling leads to a nontrivial limiting process. Since we are interested in studying fixation events, we choose to normalize the state space by the population size, so that

*b*=

_{N}*N*; the general theory then indicates that the time scaling must be proportional to

*t*units of continuous time correspond to

*t*/

*c*discrete generations. Using (4), we see that in the Eldon–Wakeley model, when 0 <

_{N}*γ*< 2, time must be scaled proportional to

*N*

^{−}

*and, when*

^{γ}*γ*≥ 2, proportional to

*N*

^{−2}, which coincides with the usual Moran scaling. In the

*Appendix*, we give an independent derivation of these constraints on

*c*.

_{N}The choice of time scaling then dictates the appropriate scaling for the selection coefficient *s*, to produce a nontrivial balance between selection and drift. In the *Appendix*, we show that this scaling must satisfy

### The limiting generator

In the *Appendix* we derive the continuum limit of the Eldon–Wakeley process. Such limits are characterized by an operator *G*, the infinitesimal generator for the continuum process, and an associated Kolmogorov backward equation, analogous to the diffusion equation of Kimura. As discussed above, we scale time by *c _{N}* =

*N*

^{−}

*for*

^{γ}*γ*< 2 and by

*c*=

_{N}*N*

^{−2}for

*γ*≥ 2. In accordance with (7), we define the corresponding scaled selection coefficients as

*β*= lim

_{N}_{→∞}

*sN*for

*γ*≥ 2 and

*β*= lim

_{N}_{→∞}

*sN*

^{γ}^{−1}when

*γ*< 2. Let

*t*. Then

*u*(

_{N}*x*,

*t*) converges to

As with the original discrete process, the limiting continuum process splits into three regimes, depending on the jump frequency *γ*. In all three regimes the operator on the right-hand side of (8) decomposes into two terms: a portion *G*_{D} that describes genetic drift. Note that the term describing selection is the standard first-order advection term in Kimura’s diffusion equation. The portion describing genetic drift, however, may differ from the standard form of Wright–Fisher drift.

We first study the behavior of the Eldon–Wakeley model without selection, meaning the continuum process with *G*_{D} alone.

For *γ* > 2, the discrete Eldon–Wakeley model is well approximated by the standard Kimura diffusion, as *N* → ∞. In other words, when large family sizes are sufficiently rare, then the continuum limit corresponds to the standard diffusion used in population genetics. In this regime, the Eldon–Wakeley process behaves much like the standard Moran model.

When *γ* < 2, the generator *G*_{D} corresponds to a neutral two-type model known as the Λ-process (Pitman 1999; Möhle 2001; Birkner and Blath 2009) with Λ a Dirac measure at the point *λ*. In this process, allele frequencies exhibit a “jump and hold” behavior: jumps that increase the allele frequency occur with probability *x* and have size (1 − *x*)*λ*, and jumps that reduce the allele frequency occur with probability (1 − *x*) and have size *λx*. In between these jumps the allele frequency remains fixed, for an exponentially distributed amount of time whose mean does not depend on the state *x*.

At the critical point *γ* = 2, the continuum limit is a mixture of the Λ-process and Kimura’s diffusion: jumps still occur as in the case *γ* < 2, but the allele frequency follows classical genetic drift between the jump events.

The splitting into three regimes that we observe in the forward-time process is mirrored by an analogous trichotomy in the associated coalescent theory of these processes, as described in Eldon and Wakeley (2006). For *γ* > 2, the retrospective genealogy is the classical Kingman coalescent, whereas *γ* < 2 is described by a Λ-coalescent with Dirac measure at *λ*. The critical case *γ* = 2 again mixes the two processes. This correspondence between the forward-time and backward-time limits is of course not accidental, but rather reflects the duality theory of Möhle (2001).

### An intuitive explanation of forward-time dynamics

There is an intuitive explanation for the trichotomy of regimes discussed above. Recall that a process moving only by Moran steps (*i.e.*, when the offspring distribution *P _{U}* is concentrated at 2) requires at least

*O*(

*N*

^{2}) time steps before the allele frequency changes a fixed fraction of the population size

*N*. More generally, in the discrete Eldon–Wakeley model (1), jump events occur every

*O*(

*N*) steps, and Moran steps occur in between the jump events. As a result, when

^{γ}*γ*< 2, the jumps occur too rapidly for the intervening Moran steps to change the allele frequency substantially—so that, in the continuum limit, allele frequencies are constant in between the jumps. In the case

*γ*> 2, by contrast, strictly more than

*O*(

*N*

^{2}) Moran steps occur between jump events, and the process is likely to have fixed at 0 or 1 during the intervening time as a result of the Moran steps. Thus the process with

*γ*> 2 is completely governed, in the limit, by the Moran process. Finally, at

*γ*= 2 a balance is achieved: jumps occur every

*O*(

*N*

^{2}) steps, but this time period is of precisely the correct order for the intervening Moran steps to move allele frequencies some fraction of the population size. Thus allele frequencies follow Moran-type drift in between jumps, when

*γ*= 2.

The preceding discussion applies in the absence of selection (*i.e.*, *β* = 0). Adding selection does not change the frequency or the size of jumps. However, when selection is present and *γ* < 2, then the periods of constant allele frequency between jumps are replaced with deterministic logistic growth. Similarly, when selection is present and *γ* - 2, the periods of Kimuran drift between jumps are replaced by Kimuran drift with directional selection. Sample paths illustrating these behaviors are shown in Figure 2, along with comparisons to the Moran process.

## Selection Amplification in the Eldon–Wakeley Process

The forward-time continuum limit (8) allows us to study the Eldon–Wakeley model in the presence of selection, which was not possible from a continuum limit of the (neutral) genealogies alone (Eldon and Wakeley 2006). In this section we focus on the case when large family sizes occur sufficiently often that the Eldon–Wakeley model diverges from the standard Moran model; *i.e.*, *γ* < 2. We show that such populations exhibit a remarkable behavior: an advantageous allele has markedly higher probability of fixation relative to the classical Moran process with the same selection pressure and the same variance effective population size.

In the following, we restrict ourselves to a discussion of the case of positive selection; the techniques presented can be easily adapted to investigate the effects of deleterious mutations. As with the standard Moran model, the Eldon–Wakeley process *X _{t}* will eventually absorb at either 0 or 1, and the fixation probability, as a function of initial allele frequency

*x*, can be expressed as

*π*(

*x*) = ℙ(

*X*

_{∞}= 1 |

*X*

_{0}=

*x*). From the general theory of Markov processes, this quantity is the unique stationary solution of the backward equation (8),

*π*(0) = 0,

*π*(1) = 1. Equation 10 is simply the analog of the diffusion equation Kimura used to compute the fixation probability of the Moran model (Kimura 1955). However, since

*G*

_{D}is now a complicated operator, (10) does not possess a simple closed-form solution. Nonetheless, without explicitly solving (10), we can still use it to study properties of the fixation probability

*π*. To begin, we can derive a lower bound on the fixation probability

*π*(

*x*) for a mutant initially at frequency

*x*(see

*Appendix*), which takes the form, for positive selective values

*β*≥ 0,

*p*=

*p*

_{r}is the positive root of the equation

The bound (11) provides significant information about the fixation probability of an adaptive allele. For fixed *λ*, the root *p*_{r} is a decreasing function of selection pressure *β*, and so our bound on the fixation probability increases with the strength of selection, as expected. When *β* > *λ*, then *p*_{r} < 1 and *π* has infinite derivative at *x* = 0, whereas the analogous fixation curve in the standard Moran model has finite derivative at zero, for all fitness values. In other words, for sufficiently strong selection (*β* > *λ*) the fixation probability of a rare mutant is much higher under the Eldon–Wakeley model than under the Moran model with the same selection pressure and the same variance effective population size. Thus, the Eldon–Wakeley process amplifies the impact of selection, relative to the standard Moran model, even though the two processes have the same effective population size.

The extent of selection amplification can be quite dramatic. In fact, the fixation of the advantageous allele can even be ensured, in some regimes. For the Eldon–Wakeley model with *γ* < 2 there exists a *finite* value of the selection pressure *β*, which we call *β**, for which *p*_{r} decreases all the way to zero, and hence *π*(*x*) = 1 for all *x* > 0. This value occurs at precisely*Appendix*). In other words, alleles introduced at any positive fraction of the population size will *assuredly* fix when the selection pressure exceeds the finite threshold, *β**, despite the presence of genetic drift. Figure 2B shows a typical sample path illustrating this behavior: there, *β* = 0.5 exceeds the threshold *β** = −log(0.75). Moreover, we demonstrate in the *Appendix* that the value *β** is indeed critical: if *β* < *β**, then fixation is not ensured for some starting frequency *x*. Thus *β** gives the precise value of selection pressure under which assured fixation occurs for all positive initial frequencies *x*, in the Eldon–Wakeley process with *γ* < 2. By contrast, assured fixation never occurs in the classical Moran model or in the Eldon–Wakeley model with *γ* > 2.

Assured fixation of an advantageous mutant is a property of the continuum limit; in the discrete model, assured fixation is achieved only as the population size increases to infinity. Figure 3 displays an example of this phenomenon, in which *π*(*x*) is plotted as a function of initial frequency *x* with fixed selection pressure *β* > *β**, for both the Eldon–Wakeley model and a variance-equivalent Moran model (in the sense that the reference Moran model is chosen to have the exact same variance effective population size as the corresponding Eldon–Wakeley model, for each population size *N*). As the population size gets large, the Eldon–Wakeley fixation curves approach *π*(*x*) ≡ 1, whereas the equivalent Moran fixation probability approaches Kimura’s standard equation:

### Fixation of new mutants

The analysis above focused on the fate of beneficial mutants initiated at a fixed fraction of the population size. We are naturally also interested in the fate of a new mutant, initially at frequency 1/*N*. While the continuum theory cannot be used formally to obtain results on the fixation probability of new mutants, the extremely large derivative of *π*(*x*) near *x* = 0 suggests that selection amplification also occurs when the Eldon–Wakeley process is initiated at a single mutant. We have investigated this behavior numerically.

Figure 4 shows the fixation probability of a new mutant with fixed selective advantage *s*, as a function of population size, for both the Eldon–Wakeley model and the variance-equivalent Moran model. Note that since the offspring variance

Figure 5 shows the fixation probability of new mutants as a function of their selective advantage *s*, again for an Eldon–Wakeley model and for a variance-equivalent Moran process. For all values of *s*, the Eldon–Wakeley model amplifies the efficiency of selection relative to the Moran model and especially so for larger values of *s*. By contrast, the Moran fixation probability is relatively small and approximately linear with *s* in the regime explored.

Figure 6 shows the fixation probability of a new mutant as a function of the model parameters *λ* and *γ*, with the population size and selection coefficient fixed. In Figure 6A, when *γ* < 2, the fixation probability is a decreasing function of the jump size *λ*. Here, offspring variance increases with *λ*, which decreases the fixation probability. A similar situation occurs in Figure 6B, where the jump frequency *γ* is varied, and all other parameters are fixed. For *γ* > 2, the process is similar to the classic Moran model and it has fixation probability ≈*s*. As the rate of bottleneck events increases (smaller *γ*), the offspring variance decreases and so the fixation probability declines.

### A fitness estimate for Pacific oysters

Eldon and Wakeley (2006) applied their model of large family sizes to sampled genetic data from Pacific oysters (Hedgecock 1994). They estimated the model parameters *θ* = 2*μN ^{γ}*

^{−1}and

*λ*as

Equation 13 can be used to estimate a critical selection coefficient that ensures a high probability that an advantageous mutant will fix, assuming *γ* < 2. Substituting *β* = *sN ^{γ}*

^{−1}, we find that the critical selective value is

*μ*- 10

^{−5}, we find that

*N*

^{γ}^{−1}= 1.54 × 10

^{3}, and so

*s** = 5.24 × 10

^{−5}. This critical value of the selection coefficient is quite small, in the sense that such a mutant would have very little chance of fixing under the standard Moran model. However, under an Eldon–Wakeley model for

*γ*< 2, the same mutant has a reasonable probability of fixation, provided the population size is large enough. Thus, the offspring mechanism of oysters not only changes the expected lengths of gene genealogies (Eldon and Wakeley 2006), but also radically alters the behavior of advantageous mutants, compared to the standard Moran model.

## Discussion

Eldon and Wakeley observed that a skew in an organism’s offspring distribution can drastically change the predictions for neutral genetic diversity in a population (Eldon and Wakeley 2006). Their analysis was based on the genealogical properties of a sample from the population. Here we have analyzed the forward-time properties of the same model introduced by Eldon and Wakeley. To do so, we developed a continuum approximation, which generalized the diffusion approximation of Kimura. Our analysis of the Eldon–Wakeley process in forward time has uncovered three distinct regimes of behavior, analogous to the distinct coalescent limits discovered by Eldon and Wakeley.

The virtue of a forward-time analysis is that it allows us to study the effects of natural selection, in a population with reproductive skew. We have found that the skewed offspring distribution greatly amplifies the effects of selection. In particular, when the frequency of large family sizes is sufficiently large (*γ* < 2), the fixation probability of an advantageous allele is far greater in the Eldon–Wakeley model than it is in a corresponding Moran or Wright–Fisher model with the same selection coefficient and variance-effective population size. Moreover, in this regime the fixation of a beneficial allele can even be ensured in the (continuum limit of the) Eldon–Wakeley model—a phenomenon that never occurs in the standard population-genetic model or its diffusion limit.

Our results on the Eldon–Wakeley model can be seen as a special case of general studies initiated in Der *et al.* (2011), where it was shown that, under minimal assumptions on the structure of drift, the Moran model possesses the *strongest* type of genetic drift and minimizes the fixation probability of advantageous mutants. The intuition developed there, in terms of the strength of higher-order moments of the offspring distribution and their tendency to lengthen the holding time of states, also holds for the Eldon–Wakeley model. As Equation 3 shows, the Eldon–Wakeley process can be viewed a mixture between two processes: (1) a Moran component, which is effectively a low-variance (high *N*_{e}) process, and (2) a rarer jump component, which is a high-variance (low *N*_{e}) process. These two components are mixed in such a way that, as an aggregate, they possess the same variance as the standard Moran process. Nonetheless, selection amplification occurs because new mutants can rise to an appreciable frequency through the low-variance Moran component, which becomes logistic growth in the limit of large populations. The high-variance bottleneck events disrupt this trajectory and reduce the fixation probability, but not sufficiently often enough for a commensurate decrease in the fixation probability of the advantageous mutant.

### Accuracy of the continuum limit

Our analysis has relied in parts on a continuum approximation of the discrete process. For a given model, these approximations are guaranteed to be accurate for large population size, *N*. The accuracy of this representation, however, is not uniform over the space of model parameters *λ* and *γ*.

Some idea of the accuracy of the limit can be obtained by recalling the offspring variance (Equation 4) of the process:*N*^{−1}, corresponding to the Moran component, is small relative to the jump term *λ*^{2}*N*^{−}^{γ}^{+1}, we can expect the *γ* < 2 continuum limit to be accurate, and vice versa for *γ* > 2. In other words, Kimura’s diffusion is a good approximation as *γ* becomes large, and the *λ*-process is a good approximation as *γ* becomes small. Near *γ* = 2, however, very large values of the population size are required for the continuum limit on either side to be an accurate representation of the discrete process. Large values of population size are also required when *γ* < 2 and the jump-size *λ* is small or, conversely, when *γ* > 2 and the jump-size *λ* is large.

Regardless of the accuracy of the continuum approximation, the Eldon–Wakeley process exhibits marked selection amplification and a qualitatively different behavior from the Moran model, even at modest population sizes. These results are demonstrated in Figures 2–5, which were produced by simulations of the discrete Eldon–Wakeley Markov chains, as opposed to using their continuum approximation.

### Alternative mechanisms of selection

We introduced fitness differentials into the Eldon–Wakeley model, using a form of viability selection. This involved biasing the allele frequencies logistically before applying the action of genetic drift and has the interpretation of introducing a preparatory stage in which not all individuals survive to the moment of reproduction. This form of selection has the virtue that it changes the continuum limit simply by adding the usual advective term

Another way to incorporate fitness differences is fertility selection, where the sampling of gametes is biased itself. In this model, individuals of type 1 would have (1 + *s*) times greater likelihood of being selected to reproduce than type 2 individuals, per capita. It can be shown that in this version of selection, the continuum limit still involves adding only the selective term *γ* < 2. For *γ* ≤ 1, however, fertility selection alters the neutral drift component of the continuum limit as well. It can be demonstrated, however, that all of our results concerning amplification still hold: indeed the selection amplification is as strong or even *stronger* under fertility selection, compared to viability selection.

### The absorption time and coming down from infinity

There is an important connection between the absorption time of the forward-time process and the concept of “coming down from infinity” (Pitman 1999; Schweinsberg 2000) in the associated coalescent. A coalescent initialized with a partition possessing an infinite number of blocks comes down from infinity when there are only finitely many blocks in the coalescent partition for all positive times. Let *τ* be the time to absorption for the continuum Eldon–Wakeley model, for *γ* < 2, started at some interior state. When *λ* < 1, the process cannot jump from an interior state to the boundary, but rather only to another interior state. Since the time between jumps is independent and identically exponentially distributed, and since the process is driven by the strictly interior logistic curve between jumps, it follows that the time to absorption is almost surely infinite. Note that this is in contrast to the Moran model, whose diffusion limit has a finite absorption time. This pathology occurs only in the continuum model, however; indeed it can be shown that the discrete Eldon–Wakeley model has finite absorption time with expectation at most *O*(*N ^{γ}* log

*N*) steps (Der

^{γ}*et al.*2011).

Likewise, it can be demonstrated that the Eldon–Wakeley model for *γ* < 2 corresponds to a generalized coalescent that “does not come down from infinity,” so that the time to coalescence for the entire population tends to be longer than that of a Moran process (Birkner and Blath 2009). This lengthening of fixation and absorption times does not appear to be merely accidental, and we speculate that for a general population process, there is a close relationship between the characteristics of (1) selection amplification, (2) forward-time absorption, and (3) the property of coming down from infinity.

### Model generalizations

Our analysis has focused on the case where the offspring distribution *P _{U}* is a mixture of two components: a Moran component and a bottleneck component where the reproducing individual contributes a fixed fraction

*λ*of the population size to the succeeding generation (Eldon and Wakeley 2006). More generally, the offspring distribution

*P*can be conceived as a mixture of bottleneck components over a range of jump sizes

_{U}*λ*, where

*λ*is drawn from some probability measure Λ on [0, 1]. Under this scenario, the continuum limit for the neutral Eldon–Wakeley processes takes the full form of a Λ-process:

Each of the continuum limits described above is a special case of this general formula: when *γ* > 2, Λ is a Dirac function concentrated at 0 and *σ*^{2} = 2, giving the standard limit for the Moran model; when *γ* < 2, Λ is a Dirac function concentrated at the point *λ* with *σ*^{2} = *λ*^{2}; and finally when *γ* - 2, Λ is mixture of atoms at 0 and at *λ*.

In general, the full Λ-process also exhibits selection amplification relative to the Kimura diffusion model. The strength of this amplification can be characterized by the number *β** we have defined above—the selective value required for assured fixation of a mutation initiated at positive frequency. One can demonstrate that the following formula holds for *β**,*λ*. Note that when Λ has any mass at zero, then *β** = ∞—that is, any process for which the continuum limit possesses a nontrivial Moran component cannot exhibit assured fixation for finite values of the selective pressure *β*.

## Acknowledgments

The authors thank two anonymous reviewers for their feedback, which improved the manuscript. J.B.P. acknowledges funding from the Burroughs Wellcome Fund, the David and Lucile Packard Foundation, the James S. McDonnell Foundation, the Alfred P. Sloan Foundation, and grant D12AP00025 from the U.S. Department of the Interior and Defense Advanced Research Projects Agency.

## Appendix

In the sections that follow we give an outline of the arguments involved, without specifying all the mathematical details.

### Description of the Jump Matrix

The exact specification of the forward-time jump matrix **J** appearing in (3) is as follows. Let *Y* be a random variable with distribution corresponding to the *i*th row of **J**. Conditioned on the event *E*_{1} that a type 1 individual is sampled for extreme reproduction,*H*_{1} is a random variable with the hypergeometric distribution of sample size *N* − 1, initial configuration *i* − 1, and *N* − *λN* draws. Conditioned on the event *E*_{2} that a type 2 individual is sampled for extreme reproduction,*H*_{2} is a hypergeometric random variable with sample size *N* − 1, initial configuration *i*, and *N* − *λN* draws. The full variable *Y* is then a mixture of *Y* | *E*_{1} and *Y* | *E*_{2}, with respective weights *i*/*N* and 1 − *i*/*N*.

### Derivation of the Limiting Generator

Let **P*** _{N}* =

**S**

_{N}**Q**

*be the transition matrix of the Eldon–Wakeley process with selection, where we have made explicit the dependence on population size*

_{N}*N*. We choose to study a continuum limit whose state space corresponds to allele frequencies lying in [0, 1];

*i.e.*, we scale the discrete state space by the population size

*N*. With this choice, the continuum representation for the sequence

**P**

*under time-scaling*

_{N}*c*involves finding an operator

_{N}*G*acting on smooth functions

*u*for which

*u*are samples from the smooth test function

_{N}*u*on the unit interval, so that

More precisely, it is necessary to prove that for any smooth function *u* on [0, 1], one has*N*}. From here on we use the expression (A3) as a short-hand for (A4).

We derive the conditions on *c _{N}* that result in a nontrivial continuum limit. Let us first consider the process in the absence of selection, so that

**P**

*=*

_{N}**Q**

*. The neutral drift component*

_{N}**Q**

*is a convex combination of the Moran process*

_{N}**M**

*and a jump process*

_{N}**J**

*:*

_{N}**Q**

*= (1 −*

_{N}*N*

^{−}

*)*

^{γ}**M**

*+*

_{N}*N*

^{−γ}

**J**

*. Thus*

_{N}#### Case 1: *γ* > 2:

From the classical theory of the Moran model, we know that if the state space is scaled by *N*, then *N*^{2}, then the first term goes to zero, as must the second term, since ‖(**J*** _{N}* −

**I**)

*u*‖ ≤ 2‖

_{N}*u*‖ < ∞, and the second term will be

*o*(1) as

*γ*> 2. Thus, the only scaling that leads to a nontrivial limit must be

#### Case 2: *γ* < 2:

Let us analyze the jump term. Fix a number 0 < *x* < 1 and consider the row *i* = [*xN*]. From the earlier description of the jump matrix, we can decompose*Y*_{1}= *Y* | *E*_{1} and *Y*_{2}= *Y* | *E*_{2} are the hypergeometric distributions described earlier in the *Appendix*. Computing the mean and the variance of the hypergeometric distributions, one finds that var(*Y*_{1}/*N*) = *O*(*N*^{−1}), var(*Y*_{2}/*N*) = *O*(*N*^{−1}), and hence go to zero; thus*γ* < 2, and because of the result (A10), any scaling *N ^{γ}* must make the second term of (A5) unbounded. On the other hand, using a Taylor expansion, one finds ‖(

**M**

*−*

_{N}**I**)

*u*‖ ≤

_{N}*C*‖

*u*″‖/

*N*

^{2}, and so if we take

*N*, then all terms must converge to zero. Thus, the only scaling leading to a nontrivial continuum limit is

^{γ}#### Case 3: *γ* = 2:

The mixture continuum limit follows by combining cases 1 and 2. Here the only scaling leading to a nontrivial limit is

To summarize the time scalings of the preceding cases, *c _{N}* must have the asymptotic behavior

### Adding selection

First we address the component of the process due to selection alone. We seek an operator *G*_{S} acting on smooth functions such that*x* < 1. As discussed in the main text, **S*** _{N}* is the process that biases allele frequencies according to fitness differences. For a population process for which each transition corresponds to a single reproductive event, it can be described as moving from state

*i*= [

*xN*] to state

*j*with probability 1, where

*j*is given by

We have*γ* > 0, *s*/*N* → 0. Using a Taylor expansion in *u*,*c _{N}* be expressed as the selection scaling limits:

Finally, we can address the product **P*** _{N}* =

**S**

_{N}**Q**

*. From the algebraic identity*

_{N}**S**

*→*

_{N}**I**, one sees that the total generator is just a sum of the parts due to selection and pure drift,

*G*

_{S}is given by (A20), and

*G*

_{D}is given by the operators in (A6), (A10), or their mixture (A12), depending on

*γ*.

### Selection Amplification

We consider here the continuum process when *γ* < 2, so that the generator for the Eldon–Wakeley process *X _{t}* has the form

*β*≥ 0.

The fixation probability *π*(*x*) = ℙ(*X*_{∞} = 1 | *X*_{0} = *x*) as a function of initial frequency *x* is a harmonic function relative to *G* and satisfies the equation *Gπ* = 0 with *π*(0) = 0, *π*(1) = 1. To examine its properties we make use of the following maximum principle.

### Theorem 1

*Suppose that a function u*, *continuous on* [0, 1], *and continuously differentiable on* (0, 1), *satisfies Gu*(*x*) ≥ 0 *for x* ε (0, 1). *Then* max_{x}_{ε[0,1]}*u*(*x*) = max{*u*(0), *u*(1)}.

#### Proof:

This result follows from general results on Markov processes, but we provide a direct verification here. Say *u*(*c*) = max_{x}_{ε[0,1]}*u*(*x*) = *M*, and *c* is an interior point. We prove *u*(1) = *M*. For if *u*(1) < *M*, define the function *z*(*x*) = *x*^{2} − *c*^{2}. A computation reveals that *Gz*(*x*) = *x*(1 − *x*)(*λ*^{2} + 2*βx*), which is strictly positive on the open interval. Now look at *w*(*x*) = *u*(*x*) + *εz*(*x*), for *ε* > 0. We have the strict inequality *Gw*(*x*) > 0 throughout the open interval, and *w*(0) < *M*, *w*(*c*) = *M*, and *w*(1) = *u*(1) + *εz*(1) < *M* whenever *ε* is sufficiently small. It follows now that *w* must have an interior global maximum. Let *x*_{0} be the location of that interior maximum. But then*Gw*(*x*) > 0 on the open interval. ▪

Now as a corollary we have the following comparison principle.

### Corollary 1

*Let* π *be the fixation probability as above*, *and assume that it is continuous on* [0, 1] *and continuously differentiable on* (0, 1). *Suppose that u also has the same regularity*, *Gu* ≥ 0, *and u*(0) = 0, *u*(1) = 1. *Then u*(*x*) ≤ *π*(*x*)*. Similarly if Gu* ≤ 0 *with the same boundary conditions*, *then u*(*x*) ≥ *π*(*x*).

#### Proof:

Since *Gπ* = 0, *G*(*u* − *π*) ≥ 0, and *u* − *π* has value zero at the endpoints. Applying the maximum principle above, we find that *u* − *π* ≤ 0. If *Gu* ≤ 0, look at *G*(*π* − *u*) to get the result. ▪

We note that while the proof above had to assume *π* satisfied a regularity condition, it is not in fact necessary for the conclusion to hold.

Corollary 1 is the result we use to establish lower and upper bounds on the fixation probability, by finding functions *u* super- and subharmonic relative to *G*, with the boundary conditions *u*(0) = 0 and *u*(1) = 1.

To do this, we decompose the operator *G* as*G*_{0}*u*(*x*) = *βxu*′(*x*) − *u*(*x*) + *u*(*x* − *λx*). Now observe that *x ^{p}* is an eigenfunction of

*G*

_{0}for every

*p*≥ 0,

*S*has generally two roots, one at

*p*= 0 and the other a nontrivial root

*p*

_{r}(

*β*) that is positive for small

*β*. As

*β*increases,

*p*

_{r}decreases, and there is a critical point

*β**, for which the nontrivial root passes through zero. This occurs at precisely

*S*′(0) = 0, which can be solved to obtain

*β** = −log(1 −

*λ*).

For *β* < *β**, *p*_{r} is strictly positive. Now it is easy to show the following.

### Theorem 2

*If* *then Gu*(*x*) ≥ 0, *and so u is a lower bound on* π.

#### Proof:

Using the decomposition (A28), we find*u* is now a lower bound on *π* follows from Corollary 1. ▪

Theorem 2 establishes the result (11) and also shows, since *p*_{r} → 0 as *β* ↑ *β**, that the selection parameter *β* = *β** must produce assured fixation. To show that *β** is critical, *i.e.*, that assured fixation does not occur when the selection pressure is *strictly* smaller than *β**, one can demonstrate the following.

### Theorem 3

*Let β* < *β*. Then constants ρ* > 0 *and C* > 0 *exist such that u*(*x*) *= Cx ^{ρ}* + (1 −

*C*)

*x satisfies Gu*≤ 0.

We defer the proof of this theorem here. It follows from Corollary 1 that the *u* constructed in Theorem 3 is an upper bound on *π*, and hence *π* must be strictly <1 in a neighborhood of 0. Thus *β* = *β** is indeed the critical selection pressure for assured fixation.

## Footnotes

*Communicating editor: Y. S. Song*

- Received February 24, 2012.
- Accepted May 28, 2012.

- Copyright © 2012 by the Genetics Society of America