We consider that at the outset (i.e. before the shift in optimal phenotype), the population has attained mutation-selection-drift balance. We follow previous work modeling this balance in making several plausible assumptions about parameter ranges (e.g. Simons et al., 2018), which ensure that genetic variation in the trait is highly polygenic, and subject to effective but not catastrophically strong selection. First, we assume that the per generation, population scaled mutational input is sufficiently large to guarantee high polygenicity (specifically, that $\sqrt{2NU}\gg 1$ ). Second, we assume that the expected number of mutations affecting the trait per generation, per gamete, is small (specifically, that $U=Lu\le 0.02$ ), such that the loss in mean population fitness (i.e. the genetic load) is not too large. As an example, for this assumption to be violated in humans, the mutational target size, $L$ , would have to exceed ∼1.5 Mb (assuming that $u\approx 1.25\cdot {10}^{-8}$ per bp per generation; Kong et al., 2012; Besenbacher et al., 2016). We note that we expect our results to hold for substantially greater values of $U$ in extensions of our model in which genetic variation in the trait under consideration has pleiotropic effects on other selected traits. Third, we make the standard assumption that the selection coefficients of all alleles satisfy $s\ll 1$ , which implies that $s_e=a^2/V_S\ll 1$ (subscript $e$ for equilibrium; see below and Wright, 1931; Wright, 1935a; Turelli, 1984). Fourth, we assume that a substantial proportion of mutations are not effectively neutral, i.e., have $S=2N{s}_e⪆1$ (by " ${\displaystyle ⪆}$ "/" ${\displaystyle ⪅}$ " we mean greater/smaller than or on the same order as). This last assumption is supported by empirical estimates of persistence time of mutations contributing to quantitative genetic variation for a variety of traits and taxa (Walsh and Lynch, 2018; Sella and Barton, 2019) and by inferences based on human GWASs (Simons et al., 2018; Zeng et al., 2018; O’Connor et al., 2019; Zeng et al., 2021), which indicate that quantitative genetic variation is not predominantly neutral. Under these assumptions, the phenotypic distribution at mutation-selection-drift balance is symmetric and tightly centered on the optimal phenotype (Figure 1). Specifically, the mean phenotype exhibits tiny, rapid fluctuations around the optimal phenotype with variance ${\delta }^2=V_S/(2N)$ (Simons et al., 2018); the phenotypic standard deviation is considerably greater than these fluctuations, i.e., ${\displaystyle \sqrt{V_A}\gg \delta }$ (Section 3.2 of Appendix 3), but it is substantially smaller than the width of the fitness function, that is, ${\displaystyle V_S\gg V_A}$ (Simons et al., 2018 and Section 3.2 of Appendix 3).

Figure 1

Download asset Open asset

We then consider the response to an instantaneous shift of $\mathrm{\Lambda }$ in optimal phenotype at time $t=0$ (Figure 1), ensuring that the shift is substantial yet not immensely large when compared to the genetic variance in the trait and in terms of the reduction in mean population fitness. To these ends, we first assume that the shift in optimum is greater than the equilibrium fluctuations in mean phenotype, i.e., that ${\displaystyle \mathrm{\Lambda }>\delta }$ . Second, we assume that the shift is smaller than, or on the order of, the width of the fitness function ( ${\displaystyle \mathrm{\Lambda }⪅\sqrt{V_S}}$ ) and that the vast majority of mutations move the phenotype by much less than the distance over which fitness declines, $\sqrt{V_S}$ (i.e. that ${\displaystyle a\ll \sqrt{V_S}}$ ). These two assumptions guarantee that the maximal directional selection coefficient of alleles, attained immediately after the shift, satisfies ${\displaystyle s_d=2\cdot (\mathrm{\Lambda }\cdot a)/V_S\ll 1}$ (see below and Wright, 1931; Barton and Turelli, 1986b). Our current condition on the magnitudes of mutations (i.e. that ${\displaystyle a\ll \sqrt{V_S}}$ ) is stronger than our earlier one at equilibrium (i.e., that ${\displaystyle s_e=a^2/V_S\ll 1}$ ), but is still not particularly restrictive, as it allows for the selection coefficients of alleles at equilibrium, ${\displaystyle s_e}$ , to be as large as 1%. As a concrete example, with a population size of $N={10}^4$ , both this condition and the condition that a substantial proportion of mutations not be effectively neutral ( $S=2N{s}_e⪆1$ ) are satisfied if selection coefficients of alleles at equilibrium are exponentially distributed with mean between ${10}^{-5}$ and ${10}^{-3}$ . Third, we assume that adaptation to the new optimum requires only a small average frequency change per segregating site, which translates into the requirement that ${\displaystyle \mathrm{\Lambda }/\sqrt{V_A(0)}⪅1/2\cdot \sqrt{2NU}}$ (see Section 3.2 of Appendix 3). Both this and the previous bound on the shift size (that ${\displaystyle \mathrm{\Lambda }⪅\sqrt{V_S}}$ ) are again not particularly restrictive, in that they allow for shifts of several equilibrium phenotypic standard deviations (given that $\sqrt{2NU}\gg 1$ and $V_S\gg V_A(0)$ for this and the former bounds, respectively). Our assumptions on parameter values are summarized in Appendix 1—table 2.

When we study the allelic response, we use units based on the dynamics at mutation-selection-drift balance (i.e. before the shift in optimum). Using arbitrary units, and denoting corresponding parameter values with tildes, the population-scaled selection coefficient at mutation-selection-drift balance is ${\displaystyle S=2N{s}_e=2N{\tilde{a}}^2/{\tilde{V}}_S={\tilde{a}}^2/{\tilde{\delta }}^2}$ . We measure the trait in units of $\delta$ , the typical deviation of the population mean from the optimum. In these units, the magnitude of effect size $a=\sqrt{S}$ , the stabilizing selection parameter $V_S=2N$ , and the contribution of a segregating allele to variance is $v^{*}(a,x)=2a^{2}x(1-x)$ . As shown below, stating our results in these terms makes their form invariant with respect to the population size, $N$ , and the strength of stabilizing selection, $V_S^{-1}$ .

We compare our analytical results to three kinds of simulations, which differ in their simplifying assumptions and computational tractability (see Section 2 of Appendix 3 for further detail). The first realizes the full model described above; it is run for a burn-in period of $10N$ generations before the shift and for a period of $12N$ generations after, to attain equilibrium both before and after. The second traces all alleles rather than individuals, assuming linkage equilibrium (rather than free recombination). Changes to allele frequencies every generation are modeled according to the Wright-Fisher process. These second simulations are also run with a burn-in period of ${\displaystyle 10N}$ generations before the shift, then for $12N$ generations after. The third kind of simulation traces the trajectory of a single allele segregating at the time of the shift. To that end: (i) given the effect size of the allele, we sample initial minor allele frequencies from the closed form, equilibrium distributions (Equation A3-27 in Appendix 3), using importance sampling based on the density of variance contributed by different minor allele frequencies (Section 3.1 of Appendix 3); and (ii) the trajectory of the mean phenotype of the population over time, on which the allele dynamics depend, is given as input, based on either an analytical approximation (see below) or on an average over simulations of the second kind. The single allele simulation is run until the focal allele fixes or goes extinct. Documented code for simulations can be found at https://github.com/sellalab/PolygenicAdaptation1D (copy archived at swh:1:rev:35d0857272a3929bad9fad0856e90c24e032b5ff; Hayward and Sella, 2022).

In the main text, we use the simulations that afford the highest resolution in comparisons with analytical predictions, whereas in Section 2.2 of Appendix 3 we validate our main results against simulations that realize the full model (at a lower resolution). Specifically, we compare most of the predictions about the allele dynamics with the results of single allele simulations, running 250,000 replicas for any given allele effect size and optimum shift size (see parameter choices below). The single allele simulations do not describe phenotypic change or the trajectories of mutations that arise after the shift in optimum, however. We therefore compare the predictions for these processes with the results of the all alleles simulations; in these simulations, we run 2500 replicas with any given set of parameters.

The simulations used in the main text correspond to the two qualitative phenotypic responses described below, which we refer to as the Lande and non-Lande (see Results). Specifically, we use the following simulation parameter values:

• In all simulations, we take a populations size of $N={10}^4$ and a shift size of ${\displaystyle \mathrm{\Lambda }=2\sqrt{V_A(0)}}$ or $4\sqrt{V_A(0)}$ . Since we work in units of $\delta$ , the typical deviation of the population mean from the optimum at equilibrium, we take $V_S=2N$ .

• The single allele simulations always assume the Lande phenotypic response, which is determined by the initial genetic variance; we take an initial variance such that $\sqrt{V_A(0)}=17\cdot \delta$ .

• The all alleles simulations are specified by the mutation rate, $U$ , and the distribution of allele effect sizes squared, for which we use an exponential distribution (measuring the trait in units of $\delta$ ). We use the following parameter values:

  Lande case: ${\displaystyle U=0.03}$ and ${\displaystyle E(a^2)=E(S)=1}$ .

  Non-Lande case: ${\displaystyle U=0.01}$ and ${\displaystyle E(a^2)=E(S)=16}$ .

  These parameter choices yield the same genetic variance at equilibrium (before the shift) in both cases; specifically, $\sqrt{V_A(0)}=29\cdot \delta$ .

We first consider how the population’s mean phenotype approaches the new optimum. In Section 1.2 of Appendix 3, we express the mean distance from the new optimum, $D(t)$ , as a sum over alleles’ contributions. We show that under our assumptions, the expected, per generation change in this distance is well approximated by

where $V_A(t)$ and ${\mu }_3(t)$ denote the 2^nd and 3^rd central moments of the phenotypic distribution. The 1^st term on the right-hand side reflects selection to reduce the distance between the mean phenotype and the new optimum, which is proportional to this distance and to the additive genetic variance (Lande, 1976). The 2^nd term reflects the effect of stabilizing selection on an asymmetric (skewed) phenotypic distribution. In particular, when the mean phenotype is near the optimum (and this 2^nd term is approximately ${\mu }_3(t)/(2V_S)$ ), stabilizing selection pushes the mean phenotype in the direction opposite to the thicker tail of the phenotype distribution (see Section 1.2 of Appendix 3 for further discussion of Equation 3). Similar expressions were derived by Barton and Turelli, 1986b under the rare-alleles approximation and by Bürger, 1991 under the assumption of a parabolic fitness function.

We rely on Equation 3 to describe the phenotypic response to selection. This response takes a simple form in the infinitesimal limit (Fisher, 1918) in which genetic variation at equilibrium arises from infinitely many segregating alleles with infinitesimally small effect sizes (see Section 8 of Appendix 3 for details). In this limit, the equilibrium phenotypic distribution is Normal and it remains Normal with the same variance after the shift, because the change in mean phenotype is achieved by infinitesimally small changes to allele frequencies at infinitely many loci, with no change to the frequency distribution (Lande, 1976, and Section 8 of Appendix 3). Under these assumptions, Equation 3 reduces to

which (in continuous time) is solved by

This solution was first derived by Lande, 1976, and in what follows, we refer to it as Lande’s solution or approximation. When genetic variance is dominated by loci with small and intermediate effect sizes (as defined below), the trait is highly polygenic, and the shift in optimum is not too large relative to the phenotypic standard deviation, changes to the 2^nd and 3^rd central moments of the phenotypic distribution are small and the expected phenotypic response is well approximated by Lande’s solution (Figure 2A and Appendix 3—figures 26 and 27).

Figure 2

Download asset Open asset

More generally, given our assumptions that polygenicity is high and that the shift is not too large (see Section 6 of Appendix 3), the deviations from Lande’s approximation are usually small and their magnitude is determined by the distribution of allele effect sizes (see section on The allelic response in the equilibration phase and Section 6 of Appendix 3). Specifically, changes to the 2^nd and 3^rd central moments of the phenotypic distribution are greater when alleles with large effects contribute markedly to genetic variance (Figure 2C, Appendix 3—figures 26 and 27, and Barton and Turelli, 1986b). For some intuition, consider a pair of minor alleles with the same initial frequency and magnitude of effect, where the effect of one is aligned with the shift and the effect of the other opposes it. After the shift, directional selection increases the frequency of the aligned allele relative to that of the opposing one. The frequency increase of the aligned allele increases variance more than the frequency decrease of the opposing allele decreases it, resulting in a net increase to variance (Figure 2C; Barton and Turelli, 1986b; de Vladar and Barton, 2014; Jain and Stephan, 2017a). The relative changes in frequency and thus the net increase in variance are greater for alleles with larger effects. Next consider the 3^rd central moment. At equilibrium, the contributions of alleles with opposing effects to the 3^rd central moment cancel out. After the shift, the frequency increase of aligned alleles relative to opposing ones introduces a non-zero 3^rd central moment (Figure 2C). Large effect alleles contribute substantially more to this 3^rd central moment, plausibly because their individual contribution to the 3^rd central moment at equilibrium is greater (see Section 3.3 of Appendix 3 and Appendix 3—figure 7B) and because they exhibit larger relative changes in frequency after the shift (see section on The allelic response in the rapid phase and Section 4 of Appendix 3). The same reasoning suggests that very large shifts in optima or low polygenicity could also lead to substantial changes to the 2^nd and 3^rd central moments of the phenotypic distribution (Section 6 of Appendix 3 and Barton and Turelli, 1986b), but these cases violate our assumptions and are beyond the scope of this manuscript.

The increase in 2^nd and 3^rd central moments after the shift result in a phenotypic dynamic with two distinct phases. First, immediately after the shift, the mean phenotype rapidly approaches the new optimum, akin to the exponential approach in Lande’s approximation. In this case, however, genetic variance increases and thus the exponential rate of approach may increase, making the expected approach even faster (Equation 3). Shortly thereafter, when the mean phenotype nears the optimum, the decreasing distance and increasing 3^rd central moment reach the point at which

The two terms on the right-hand side of Equation 3 then approximately cancel out, and the dynamic enters a second, prolonged phase, in which the approach to the optimum nearly grinds to a halt (Figure 2D). During this phase, the expected change in mean phenotype can be described in terms of a quasi-static approximation given by Equation 6 (Figure 2D and Appendix 3—figure 25). The rate of approaching the optimum is then largely determined by the rate at which the 3^rd central moment decays. This roughly corresponds to the rate at which the allele frequency distribution equilibrates and mutation-selection-drift balance is restored around the new optimum (see section on Other properties of the equilibration process).

We now turn to the allele dynamics that underlie the phenotypic response. These dynamics can be described in terms of the first two moments of change in frequency in a single generation (Ewens, 2004, Chapter 4). For an allele with effect size ${\displaystyle \pm a}$ and frequency ${\displaystyle x}$ , we calculate the moments by averaging the fitness of the three genotypes over genetic backgrounds (Section 1 of Appendix 3). Under our assumptions, the moments are well approximated by

and

which is the standard drift term. Similar expressions for the first moment trace back to Wright, 1935b and have been used previously to study the response to selection on quantitative traits (Barton, 1986a; Bürger, 1991; Charlesworth, 2013; de Vladar and Barton, 2014).

The two terms in the first moment reflect different modes of selection: directional and stabilizing, respectively. The first term arises from directional selection on the trait and takes a semi-dominant (additive) form with selection coefficient $s_d=\pm 2a\cdot D(t)/V_S$ . Its effect is to increase the frequency of alleles whose effects are aligned with the shift (and vice versa) and its strength weakens as the distance to the new optimum, $D$ , decreases. The second term arises from stabilizing selection on the trait and takes an under-dominant form with selection coefficient $s_e=a^2/V_S\cdot (1-D^2(t)/V_S)$ . Its effect is to decrease an allele’s contribution to phenotypic variance, $2a^{2}x(1-x)$ , by reducing minor allele frequency (MAF); it becomes weaker as the MAF approaches 1/2.

The relative importance of the two modes of selection varies as the mean distance to the new optimum, $D$ , decreases. We therefore divide the allelic response into two phases: a rapid phase, immediately after the shift, in which the mean distance to the new optimum is substantial and changes rapidly, and a subsequent, prolonged equilibration phase, in which the mean distance is small and changes slowly (Jain and Stephan, 2017a). We define the end of the rapid phase as the time, ${\displaystyle t_1}$ , at which Lande’s approximation for the distance to the optimum $D_L\left(t_1\right)$ equals the typical deviation of the population mean from the optimum at equilibrium ${\displaystyle \delta =\sqrt{V_S/2N}}$ , i.e.,

(in Section 3.2 of Appendix 3 we show that $V_S/V_A(0)\sim 1/U$ ). This definition is somewhat arbitrary, as the transition between phases is gradual, but it roughly captures the change in allele dynamics (Appendix 3—figure 25). Moreover, our analysis is insensitive to this particular choice (we only use it in comparing analytic and simulation results for the rapid phase).

The change in mean phenotype during the rapid phase is driven by the differential effect of directional selection on minor alleles whose effects are aligned and opposed to the shift in optimum (Figure 3). Considering a pair of minor alleles with opposing effects of the same magnitude and the same initial frequency, selection increases the frequency of the aligned allele relative to the opposing one. By the end of the rapid phase, the frequency differences across all aligned and opposing alleles drive the mean phenotype close to the new optimum (Figure 2A). Deviations from Lande’s approximation manifest as prolonged, weak directional selection during the equilibration phase, which further increases the expected frequency difference between aligned and opposing alleles. However, given that we are considering a highly polygenic trait, the expected frequency difference between a pair of opposing alleles will be small. This small difference causes aligned alleles to have a slightly greater probability of eventually fixing during the equilibration phase (Figure 3). Over a period on the order of $2N$ generations (see below), the frequency differences between aligned and opposing alleles are replaced by a slight excess of fixed differences between them, and the equilibrium genetic architecture is restored around the new optimum. In the following sections, we describe these processes quantitatively: Specifically, we ask how the relative contribution of alleles to phenotypic change during the two phases depends on their effect size and initial frequency.

Figure 3

Download asset Open asset

We can describe changes to allele frequencies during the rapid phase with a simple deterministic approximation. The duration of the rapid phase is much shorter than the time scale over which genetic drift has a substantial effect ( $t_1\sim 1/U\ll 2N$ generations; see Equation 9), allowing us to rely only on the first moment of change in allele frequency (Equation 7). Additionally, deviations of the distance $D(t)$ from Lande’s approximation during this phase have negligible effects (Figure 2B and Appendix 3—figure 8D–F), allowing us to assume that ${\displaystyle D(t)=D_L(t)}$ (Equation 5). Lastly, when relative frequency changes are small, we can substitute the frequency in the first moment by its initial value. With these simplifications, we can integrate the first moment over time to obtain an explicit linear approximation for frequency changes.

Consider a pair of minor alleles with opposing effects of size ${\displaystyle \pm a}$ and initial frequency ${\displaystyle x_0}$ before the shift in optimum. Using our linear approximation, we find that the frequency difference between them at the end of the rapid phase is

The contribution of the pair to the change in mean phenotype is

where $v^{*}(a,x_0)=2a^2{x}_0(1-x_0)$ is the contribution to variance of an allele with magnitude $a$ and frequency ${\displaystyle x_0}$ . Thus, the pair’s contribution to phenotypic change is proportional to its contribution to phenotypic variance before the shift in optimum.

The expected contribution of alleles with a given magnitude and initial frequency is therefore proportional to their expected contribution to phenotypic variance at equilibrium, before the shift occurs. We focus on the contribution of alleles divided by the mutation rate at which they are introduced into the population, that is the ‘contribution per unit mutational input’. To this end, we measure the trait value in units of ${\displaystyle \delta =\sqrt{V_S/(2N)}}$ and express allele magnitudes in terms of the scaled selection coefficients at equilibrium (when $D=0$ ); in these units $S=2N{s}_e=a^2$ (see Choice of units). Expressing our results in this form makes them invariant with respect to changing the population size, $N$ , stabilizing selection parameter, $V_S$ , mutational input per generation, $2NU$ , and distribution of magnitudes, $g(a)$ . In these terms, the expected contribution of alleles with given magnitude and initial MAF to phenotypic change is

and the marginal contribution of alleles with a given magnitude is

where ${\displaystyle v(a,x_0)\approx 4a^2\cdot \text{Exp}\left[-a^2{x}_0(1-x_0)\right]}$ and ${\displaystyle v(a)\approx 4a^2\cdot {\int }_0^{1/2}\text{Exp}\left[-a^{2}x(1-x)\right]dx=4a\cdot D_{+}\left(a/2\right)}$ are the corresponding densities of variance per unit mutational input at equilibrium, and $D_{+}$ is the Dawson function (Section 3.2 of Appendix 3). The expected absolute contributions follow from multiplying these expressions by the mutational input per generation, $2NU\cdot g(a)$ . Specifically, as we would expect, the total change in mean phenotype during the rapid phase is ${\displaystyle 2NU\cdot \mathrm{\Delta }{z}_{t_1}=2NU\cdot {\int }_0^{\mathrm{\infty }}\mathrm{\Delta }{z}_{t_1}(a)\cdot g(a)da\approx \mathrm{\Lambda }-D_L(t_1)}$ , as ${\displaystyle V_A(0)=2NU\cdot {\int }_0^{\mathrm{\infty }}v(a)\cdot g(a)da}$ .

The relative contribution of alleles with a given magnitude and initial MAF to phenotypic change follows from their expected contribution to variance at equilibrium (Equations 12 and 13, and Figure 4). The properties of $v(a)$ imply that (Figure 4A): (i) the relative contribution of alleles with small effect sizes ( $a^2\ll 1$ ) scale linearly with $S=a^2$ ( ${\displaystyle v(a)\approx 2a^2}$ , measured in units of ${\delta }^2$ ); (ii) the contribution of alleles with moderate and large effect sizes (roughly ${\displaystyle S=a^2>3}$ ) are much greater, and fairly insensitive to the effect size (with $v(a)\approx 4$ ); and (iii) the contribution is maximized for $a^2\approx 10$ ( ${\displaystyle v(\sqrt{10})\approx 5.2}$ ) (see Simons et al., 2018, for intuition about these properties). While large and moderate effect alleles make similar contributions to phenotypic change, MAFs of large effect alleles before the shift are much lower than MAFs of moderate ones (Figure 4B), because they are subject to stronger stabilizing selection. The expected frequency difference between pairs of opposing alleles is greatest for moderate effect sizes (Figure 4C), because it is proportional to $E\left(2a{x}_0(1-x_0)\right)\propto v(a)/a$ (Equation 10), and $v(a)$ is similar for moderate and large effect sizes. Additional properties of the allelic response during the rapid phase are presented in Section 4 of Appendix 3.

Figure 4

Download asset Open asset

When the polygenicity is low and/or the shift in optimum or effect sizes are large our linear approximation becomes less accurate (Figure 4A). Specifically, minor alleles exhibit large relative changes in frequency such that substituting the initial MAF for the frequency in Equation 7 for the 1^st moment is inaccurate. In Section 4.1.2 of Appendix 3 we derive a nonlinear approximation that is more accurate in these cases (Figure 4A and Appendix 3—figures 9 and 10). Nonetheless, the qualitative behaviors we outlined remain intact.

Over the long run, the small frequency differences between opposite alleles that accrue during the rapid phase translate into small differences in their fixation probabilities (Figure 3). In the non-Lande case, prolonged weak directional selection during the equilibration phase amplifies these differences in fixation probabilities. We approximate fixation probabilities in two steps. First, we model the effect of directional selection on frequency as an instantaneous, deterministic pulse. Second, we apply the diffusion approximation for the fixation probability (Ewens, 2004, Chapter 4), assuming stationary stabilizing selection ( $D=0$ ), genetic drift, and the initial frequency after the pulse. We further assume that the relative changes in allele frequencies due to directional selection are small, such that we can use approximations that are linear in this change; but in Section 5 of Appendix 3, we derive nonlinear approximations that relax this assumption.

When Lande’s approximation is accurate, directional selection is non-negligible only briefly after the shift. This justifies approximating its effects as if they were caused by an instantaneous pulse. It also suggests that mutations that arise after the shift in optimum contribute negligibly to phenotypic change, because when directional selection is non-negligible, few of them arise and their fixation probabilities are tiny (given that they start from an initial frequency of $1/2N$ ).

Consider a pair of opposite minor alleles, with magnitude $a$ and initial frequency ${\displaystyle x_0}$ . Analogously to our derivations for the rapid phase (Equations 10, 11 and 12) by modeling the effects of directional selection on their frequencies as an instantaneous pulse, and assuming that these effects are small, we find that the resulting frequency differences between them is approximated by

Consequently, the pair’s expected contribution to phenotypic change is approximated by

and the contribution of such pairs per unit mutational input is approximated by

where, as before, $v^{*}(a,x)=2a^{2}x(1-x)$ is an allele’s contribution to genetic variance and $v(a,x)$ is the density of variance per unit mutational input at equilibrium, and we use the superscript $L$ to denote that this applies to the Lande case.

We approximate a pair’s expected long-term, fixed contribution to phenotypic change by calculating the difference in fixation probabilities of the opposite alleles given their frequency after the pulse, again assuming that the effects of the pulse are small. Namely,

where $\pi (a,x)$ denotes the fixation probability of an allele with magnitude $a$ and initial frequency $x$ under stationary stabilizing selection and drift. In Section 3 of Appendix 3 we derive the diffusion approximation for $\pi (a,x)$ and show that

where ${\displaystyle f(a)\equiv 2a^3\cdot \text{Exp}\left[-a^2/4\right]/\left(\sqrt{\pi }\cdot \mathrm{E}\mathrm{r}\mathrm{f}\left[a/2\right]\right)}$ . From Equation 14-18, we find that the expected fixed contribution per unit mutational input of pairs of alleles is

Note that this expression does not depend on the initial frequency! The expected marginal contribution of alleles with a given magnitude follows and is

Hence, the function $f$ approximates how the relative long-term contribution of alleles depends on their magnitudes (Figure 5A).

Figure 5

Download asset Open asset

We expect the total long-term allelic contribution to equal the shift in optimum, $\mathrm{\Lambda }$ . In our linear Lande approximation, the total contribution is

where we use the fact that ${\displaystyle V_A(0)=2NU\cdot {\int }_0^{\mathrm{\infty }}v(a)\cdot g(a)da}$ and define

$C$ measures the extent to which our approximation underestimates the long-term phenotypic change. We note that ${\displaystyle C>0}$ for any distribution of allele magnitudes, because ${\displaystyle v(a)>f(a)}$ for any magnitude $a$ (Figure 5A), and further that $C\ll 1$ is a necessary condition for our approximation to be accurate. Given that $v(a)$ is appreciably greater than $f(a)$ only for $a^2⪆4$ (Figure 5A), this condition implies that for our approximation to be accurate, the vast majority of the genetic variance at equilibrium must arise from alleles with ${\displaystyle a^2<4}$ .

When alleles with larger effects contribute substantially to the variance at equilibrium (and $C$ is appreciable), Lande’s approximation becomes inaccurate. The prevalence of large effect alleles leads to a quasi-static decay of the mean distance from the new optimum, $D$ , during the equilibration phase (Figure 2D and the section on the Phenotypic response). For the distance during the equilibration phase, and therefore for the deviations from Lande’s phenotypic approximation to be substantial, would require that $C\gg 1$ (see Section 6 of Appendix 3), which implies that the vast majority (e.g. > 90%) of the genetic variance at equilibrium arises from alleles with large effect sizes, say with ${\displaystyle a^2\gg 4}$ (see Figure 5A and Equation 22). Even a small distance $D$ during the equilibration phase, however, would result in prolonged, weak directional selection that could markedly amplify the difference in fixation probabilities between opposite alleles. Our linear Lande approximation does not account for this amplification, and it could therefore greatly underestimate the total long-term allelic contribution.

We can, however, extend our approximation to account for the amplification in the non-Lande case. To this end, we modify our instantaneous pulse approximation for a pair of opposite alleles (Equation 14) to

where the factor ${\displaystyle A>0}$ accounts for the greater effect of directional selection relative to the Lande case (with the caveat that the justification for the instantaneous pulse approximation is less obvious in the non-Lande case, given prolonged, weak directional selection; see Section 5.2 of Appendix 3). Following the same steps as taken in the Lande case, we then find that the expected long-term (fixed) contribution per unit mutational input of pairs of alleles with a given magnitude and initial MAF is

and that their expected marginal contribution for a given magnitude is

where the superscript $v$ denotes that these contributions originate from variation that segregated before the shift in optimum.

In the non-Lande case, the fixation of mutations that arise after the shift in optimum can also contribute substantially (Appendix 3—figure 21B), because prolonged, weak directional selection can produce a substantial difference in the numbers of fixations of mutations with opposite effects. In Section 5.2.2 of Appendix 3, we follow the same approach that applied for standing variation to show that the relative long-term contribution of new mutations with a given magnitude can be approximated by

where the factor ${\displaystyle B>0}$ does not dependent on the magnitude (Section 5.2.2 of Appendix 3). Thus, similar to what we found in the Lande case, the function $f$ approximates how the relative long-term contribution of alleles depends on their magnitudes, but here it applies to both standing variation and new mutations (Figure 5B).

To gain further understanding of the non-Lande case, we consider the joint contribution of standing variation and new mutations. Equating the total contribution with the shift in optimum we find that

with $C$ defined in Equation 22. This implies that $A+B=C$ and that the proportional contributions of standing variation and new mutations are $(1+A)/(1+C)$ and $B/(1+C)$ , respectively. It also implies that the contribution per unit mutational input of alleles with a given magnitude $a$ is

Thus, in the linear non-Lande approximation, prolonged weak directional selection amplifies the relative contribution of alleles of any given magnitude by the same factor of $(1+C)$ (Figure 5B and Appendix 3—figure 19). $C$ is therefore an allelic measure of the deviation from Lande’s approximation (see Appendix 3—figures 20 and 21), and, intriguingly, it depends only on the distribution of mutation magnitudes (Equation 22).

In Section 5 of Appendix 3, we show that the linear approximations are accurate when ${\displaystyle a\cdot (1+A)\cdot \mathrm{\Lambda }/V_A(0)\ll 1}$ (with $A=0$ in the Lande case) and we derive nonlinear approximations that are more accurate when this condition is violated (Figure 4A and Appendix 3—figures 13 and 16). When polygenicity is low, the shift in optimum is large, or effect sizes are large, directional selection causes large relative changes in MAFs, such that the use of the initial MAF in the instantaneous pulse approximation (i.e. in Equations 14 and 23) and the Taylor approximation of fixation probabilities (Equation 17) become inaccurate. Even in these cases, however, the linear approximations capture the salient features of the long-term allelic contribution to phenotypic adaptation (Figure 5 and Appendix 3—figures 13–19).

Notably, our linear approximations capture the dramatic turnover in the genetic basis of adaptation during the equilibration phase (Figure 6). In the long run, the short-term contribution of large effect alleles ( $S=a^2⪆30$ ) is almost entirely wiped out, and is supplanted by the contribution of moderate effect alleles ( $a^2\approx 5$ ) (Figure 6A and Appendix 3—figures 15 and 19). Moreover, for any given magnitude, the proportional long-term contribution of minor alleles that segregated at low frequencies before the shift is diminished relative to their short-term contribution, all the more so for large effect sizes (Figure 6B and Appendix 3—figure 14). For instance, for an effect size $a^2=35$ , minor alleles with initial frequencies below 0.05 account for more than 99% of the short-term contribution but for only ∼10% of the (much smaller) long-term contribution (Figure 6B and Appendix 3—figure 14).

Figure 6

Download asset Open asset

We can understand this turnover by considering the effects of stabilizing selection during the equilibration phase (Appendix 2—figure 2). As noted, stabilizing selection on the trait induces selection against minor alleles, which weakens as MAF increases and vanishes at MAF $=1/2$ . Now consider how it affects a pair of alleles with opposite, moderate or large effect. If their initial frequencies are very low, both alleles will have low MAFs at the end of the rapid phase. Consequently, they will both be strongly selected against during the equilibration phase and will almost certainly go extinct (Appendix 2—figure 2). In the long run, their expected contribution to phenotypic adaptation is therefore diminished. In contrast, if the alleles’ initial MAF is sufficiently high, the relative increase in the aligned allele’s frequency by the end of the rapid phase causes it to be subject to substantially weaker selection than is the opposing allele. In the extreme in which the aligned allele has exceeded frequency $1/2$ , the direction of selection on it is even reversed (Appendix 2—figure 2). In such cases, the pair’s expected contribution to phenotypic adaptation will be amplified. This reasoning suggests that, for a given magnitude, there is a critical initial MAF such that the long-term contribution of alleles that start above it is amplified and the contribution of those that start below it is diminished (Figure 6B and Appendix 2—figure 3). This critical frequency is lower in the non-Lande case, because prolonged, weak directional selection amplifies the long-term contribution from alleles with any initial MAF.

The turnover among alleles with different magnitudes can be explained in similar terms. Alleles with large magnitudes almost always start from low MAF, because they are subject to strong stabilizing selection before the shift (Figure 4B). Consequently, they are highly unlikely to exceed the critical initial frequency and their expected long-term contribution to phenotypic adaptation is diminished (Figures 5A and 6A). However, changes to the frequencies of these alleles skew the phenotypic distribution, leading to prolonged, weak directional selection that amplifies the long-term contribution of small and moderate effect alleles (Figure 6B). In our linear approximation, this amplification will occur for allele magnitudes that satisfy $(1+C)\cdot f(a)⪆v(a)$ (Figure 5B). These considerations explain why the contributions of alleles with small and moderate effects supplant those of alleles with large effects (Figures 5A and 6A). They also highlight that deviations from Lande’s approximation are critical to understanding the allelic response, even when they have small phenotypic effects.

While long-term phenotypic adaptation arises from an excess in fixations of aligned relative to opposing alleles, this excess and the increase in the total number of fixations are typically small relative to the number of fixations at equilibrium (Figure 7). The relative excess of aligned fixations (defined as $(\pi (a)-\pi (-a))/(\pi (a)+\pi (-a))$ ) and the relative increase in total fixations (compared to the expected number at equilibrium) both decrease with increased polygenicity and increase with the shift in optimum and allele magnitude (Figure 7B and C). For sufficiently large effect sizes, practically all fixations are caused by the shift and are of aligned alleles (Figure 7B and C). However, with the exception of extreme cases in which the contribution of alleles of small and moderate effects to genetic variance is negligible, the number of fixations of such large effect alleles and their contribution to phenotypic change will be small (Figures 5 and 7). Typically, most fixations and contribution to phenotypic change will arise from alleles with small and moderate effects (Figures 5 and 7), for which the proportional excess of aligned and total fixations is modest (Figure 7B and C).

Figure 7

Download asset Open asset

In the long run, these fixations move the mean phenotype all the way to the new optimum, and genetic variation around the new optimum returns to equilibrium. A proxy for the approach to equilibrium is the ‘fixed distance from the optimum’, defined as the phenotypic distance of an individual that is homozygous for the ancestral allele at every segregating site; at equilibrium, we expect the fixed distance to be $0$ . Our simulations suggest that, under a broad range of parameter values, the change in fixed distance after the shift is well approximated by an exponential decay with a rate of $1/(2N)$ per generation (Figure 8 and Appendix 3—figure 28). This approximation is remarkably accurate in the Lande case. In the non-Lande case, the decay is initially slower than the approximation suggests, possibly because the long-term contribution of new mutations (as opposed to standing variation) takes longer to amass. In both cases, the return to equilibrium occurs on a time scale of $2N$ generations after the shift (Figure 8 and Appendix 3—figure 28).

Figure 8

Download asset Open asset

1. The basic equations

   • 1.1 Derivation of the allelic equation

   • 1.2 Derivation of the phenotypic equation

2. Simulations

   • 2.1 Choice of simulation parameters

   • 2.2 Validation of our main results with full model simulations

3. Allele dynamics under stabilizing selection

   • 3.1 Summaries of architecture

   • 3.2 Phenotypic variance at equilibrium

   • 3.3 3^rd phenotypic moment

4. Allele dynamics during the rapid phase

   • 4.1 Approximate solutions for allele trajectories

   • 4.2 Contribution to phenotypic change

5. Allele dynamics during equilibration

   • 5.1 The linear Lande approximation

   • 5.2 The linear non-Lande approximation

   • 5.3 The nonlinear Lande approximation

   • 5.4 The nonlinear non-Lande approximation

6. Phenotypic deviations from Lande’s approximation are determined by C

7. Variation over realizations of the adaptive response

   • 7.1 Variation in the phenotypic response

   • 7.2 Variation in the allelic response

8. The infinitesimal limit and Lande’s approximation

   • 8.1 The phenotypic equation in the infinitesimal limit

   • 8.2 The infinitesimal limit in our model

9. Additional figures

Our analysis relies on two basic equations:

• The allelic equation: describing the expected change in frequency per generation of an allele with effect size $\pm a$ and frequency $x$ (Equation 7 in the main text):

  (A3-1) $E\left(\mathrm{\Delta }x\right)\approx \left(\pm a\cdot D(t)/V_S\right)\cdot x(1-x)-\left(a^2/V_S\right)\cdot (1-D^2(t)/V_S)\cdot x(1-x)(1/2-x).$

• The phenotypic equation: describing the expected change in the mean distance from the optimum, $D(t)$ , per generation (Equation 3 in the main text):

  (A3-2) $E\left(\mathrm{\Delta }D(t)\right)\approx -\left(V_A(t)/V_S\right)\cdot D(t)+(1-D^2(t)/V_S)\cdot {\mu }_3(t)/(2V_S),$

  where $V_A(t)$ and ${\mu }_3(t)$ denote the 2^nd and 3^rd central moments of the phenotypic distribution.

Previous work relied on similar allelic and phenotypic equations. Assuming a fitness landscape with a single optimum, the expected change in allele frequency per generation has commonly been approximated by

(Barton, 1986a; Charlesworth, 2013; de Vladar and Barton, 2014). This equation was used to derive phenotypic equations similar to ours (see below and Barton and Turelli, 1986b; Bürger, 1991). The allelic equation has been justified for a Gaussian fitness function, i.e., $W(z)=\text{Exp}\left[-z^2/(2V_S)\right]$ (our Equation 2 in the main text), further assuming that: (i) phenotypes are normally distributed in the population; (ii) the width of the distribution is far smaller than the width of the fitness function (i.e. $V_A\ll V_S$ ); and (iii) the mean distance from the optimum is very small (i.e., $D^2\ll V_S$ ) (Barton and Turelli, 1986b; Simons et al., 2018). As we show in the main text, the assumption that the phenotype distribution is normally distributed can be violated in ways that affect the course of adaptation qualitatively, giving rise to the non-Lande case. Also, requiring that the mean distance from the optimum is very small can be restrictive.

We derive our basic equations without making the assumptions of a Normal phenotypic distribution or that the distance to the optimum is very small. Similar to previous work, we do not provide a rigorous treatment of the errors introduced by our approximations (see Hayward, 2020, for such treatment). We therefore validate our main results against simulations of the full model (Section 2.2).

We compared our analytic results with results from three kinds of simulations. The first kind implements the full model. The second traces all alleles rather than individuals. Specifically, the frequency of an allele in the next generation follows a binomial distribution $x^{\prime }\sim B(2N,x+E\left(\mathrm{\Delta }x\right))$ with $E\left(\mathrm{\Delta }x\right)$ approximated by Equation A3-9, and the number of new mutations introduced per generation is Poisson distributed with mean $2NU$ . Tracing alleles rather than individuals in this way entails two potential sources of error. First, our approximation for $E\left(\mathrm{\Delta }x\right)$ replaces the distribution of background phenotypic contributions by its mean (Section 1.1), thus neglecting the effects of variation in this distribution. Second, we assume linkage equilibrium rather than free recombination, thus neglecting the effects of short-term linkage disequilibrium (Robertson, 1956a).

The 3^rd kind of simulations traces a single allele that was segregating at the time of the shift. We sample the initial MAFs from the approximate closed form, equilibrium distribution (Equation A3-27), using importance sampling based on the density of variance at different MAFs (Equation A3-28). Changes in frequency are modeled as we described for the all alleles simulations, but here we use the mean distance to the optimum, $\overline{D}(t)$ , which is given from the outset. In the Lande case, we take the mean distance to equal the expectation under Lande’s approximation (Equation 5 in the main text), which depends only on the initial phenotypic variance ${\displaystyle V_A(0);V_A(0)}$ follows from the distribution of magnitudes and extent of polygenicity as described in Equation A3-30. In the non-Lande case, we estimate the mean distance by averaging over 2500 all alleles simulations with the same model parameters. During the equilibration phase, we average over the quasi-static approximation for the distance (Equation 6 in the main text) rather than the distance itself, because the former is considerably less variable among simulations than the latter (Section 7.1 and Appendix 3—figure 22A and B). Single allele simulations share sources of error with the all alleles simulations, in addition to potential errors arising from relying on the analytic approximation for the initial MAF and on the mean rather than on the actual (fluctuating) distance $D(t)$ .

We use the all alleles and single allele simulations throughout because they are considerably more computationally tractable, allowing us to obtain larger sample sizes and thus more precise comparisons with our analytic predictions. For example, running a single simulation, with the parameter values of the non-Lande case used in the main text (see section on Simulations and resources), takes ∼340 hours with the full model simulations and ∼10 with the all allele simulations (for a burn-in period of $10N$ generations before the shift and a period of $12N$ generations after). Unfortunately, precise estimates of fixation probabilities of alleles with a particular effect size may require hundreds or thousands (or even orders of magnitude more, if the effect size is extremely rare) of full or all alleles simulations. This problem is overcome by the single allele simulations. For example, with the same parameters, running a single allele simulation for 1000 aligned and 1000 opposing alleles segregating at the time of the shift, with each of 13 different effect sizes (those used in many of the figures, e.g. Appendix 3—figure 8) takes ∼10 hours on the same computational platform.

2.1. Choice of simulation parameters

The simulation results in the appendix are based on a wider range of model parameter values than those shown in the main text (see section on Simulations and resources). Specifically:

All simulation types: In all simulations, we use the same population size and shifts as in the main text, i.e. a populations size of $N={10}^4$ and a shift size of $\mathrm{\Lambda }=2\sqrt{V_A(0)}$ or $4\sqrt{V_A(0)}$ . Since we work in units of $\delta$ , the typical deviation of the population mean from the optimum at equilibrium, we take $V_S=2N$ (see section on Choice of units).

All alleles simulations: We run 2500 simulations for each of the following combinations of parameters:

• An exponential distribution of effect sizes squared (with the trait measured in units of $\delta$ ), with mean $E(a^2)=E(S)=1$ in the Lande case , and $E(a^2)=E(S)=16$ in the non-Lande case .

• Three mutation rates (per gamete per generation), which correspond to different degrees of polygenicity:

  Low polygenicity: ${\displaystyle U=0.005}$ ,

  Medium polygenicity: ${\displaystyle U=0.01}$ , and

  High polygenicity: ${\displaystyle U=0.02}$ .

Single allele simulations: We run 250,000 simulations for each effect size used, where:

• In the Lande case , we assume the Lande phenotypic response, which is determined by the initial genetic variance (Equation 5 in the main text). The variances are calculated from Equation A3-30 assuming an exponential distribution of effect sizes squared (with ${\displaystyle E(a^2)=1}$ ) and the same three mutation rates used in the all alleles simulations yielding:

  Low polygenicity: ${\displaystyle \sqrt{V_A(0)}=12\cdot \delta }$

  Medium polygenicity: ${\displaystyle \sqrt{V_A(0)}=17\cdot \delta }$ , and

  High polygenicity: ${\displaystyle \sqrt{V_A(0)}=24\cdot \delta }$ .

• In the non-Lande case , $\overline{D}(t)$ is the average over the corresponding all alleles simulations.

These parameter values were primarily chosen to span as wide a range as possible given our conditions on them (Appendix 1—table 2). The one exception is our choice of a population size of $N=10000$ , which was motivated by being on the order of the effective size in human populations (Schiffels and Durbin, 2014), though on the lower side in order to allow for greater computational tractability. The choice of the effect size distribution of new mutations is constrained by the requirement that a substantial proportion of them are not effectively neutral, i.e., $a^2⪆1$ (Appendix 1—table 2). In simulations representative of the Lande case, we further require that most new mutations satisfy $a^2⪅4$ (see Equation A3-66). Our choice of an exponential distribution of effect size squared with $E(a^2)=1$ satisfies both requirements, with proportion $0.63$ of mutations with $a^2\le 1$ , $0.35$ with ${\displaystyle 1<a^2\le 4}$ , and $0.02$ with ${\displaystyle a^2>4}$ (Appendix 3—figure 1A). In simulations representative of the non-Lande case, we require that a substantial proportion of mutations have effect sizes ${\displaystyle a^2>4}$ while being much smaller than the width of the fitness function, i.e., $a/\sqrt{V_S}\ll 1$ (Appendix 1—table 2). Our choice of an exponential distribution of effect size squared with $E(a^2)=16$ satisfies these requirements, with proportion $0.78$ of mutations with ${\displaystyle a^2>4}$ and only $4\times {10}^{-6}$ with ${\displaystyle a^2>V_S/100}$ .

Appendix 3—figure 1

Download asset Open asset

Given the population size, the choice of mutation rate per haploid genome per generation, $U$ , is bounded from below by the assumption that the trait is highly polygenic, i.e., $\sqrt{2NU}\gg 1$ (Appendix 1—table 2). Consequently, the lowest mutation rate that we used was $U=0.005$ , corresponding to our low end of polygenicity, i.e., $\sqrt{2NU}=10\gg 1$ . The choice of mutation rate is also bounded from above by the assumption that $U\le 0.02$ (Appendix 1—table 2). The highest mutation rate that we used was therefore $U=0.02$ , corresponding to our high end of polygenicity, i.e., $\sqrt{2NU}=20$ . As an intermediate value, we used $U=0.01$ , corresponding to $\sqrt{2NU}\approx 15$ (Appendix 3—figure 1B).

Lastly, we require the shift in optimum to be smaller than or on the order of the width of the fitness function, i.e., $\mathrm{\Lambda }⪅\sqrt{V_S}$ (Appendix 1—table 2 and Appendix 3—figure 1C). As we detail in Section 3.2, this requirement constrains how large the shift can be relative to the phenotypic standard deviation, $\sqrt{V_A(0)}$ , in a manner that depends on the mutation rate and distribution of allele magnitudes: ${\displaystyle \mathrm{\Lambda }/\sqrt{V_A(0)}\le (1/\sqrt{U})/\sqrt{{\int }_0^{\mathrm{\infty }}v\left(a\right)g(a)da}}$ (Equation A3-31 and Appendix 3—figure 5). As our largest shift, we chose $\mathrm{\Lambda }=4\sqrt{V_A(0)}$ , which, for our highest mutation rate ( $U=0.02$ ), is near the upper bound given our choices of distribution of allele magnitudes for the Lande and non-Lande cases (with $\mathrm{\Lambda }/\sqrt{V_S}=0.7$ and $1.2$ , respectively). As a smaller shift, we used $\mathrm{\Lambda }=2\sqrt{V_A(0)}$ .

2.2. Validation of our main results with full model simulations

In Appendix 3—figure 2, we compare our main results about the phenotypic dynamics with those from the full model simulation. Because of the computational intensity of the full model simulation, we use only one set of parameters, the non-Lande case used in the main text (see section on Simulations and resources) with a shift of $\mathrm{\Lambda }=2\sqrt{V_A(0)}$ , and run a sample of only 500 simulations, which is why our standard errors are sizable. The phenotypic variance before the shift is greater in the full model simulation than in our analytic approximation and all alleles simulations (Appendix 3—figure 2A), plausibly because only the full model simulation incorporates the effect of the variance of the phenotypic background of a focal site (see Section 1.1). Specifically, the variance of the phenotypic background reduces the efficacy of stabilizing selection—akin to how environmental contributions to phenotypic variance effectively increase $V_S$ (Turelli, 1984; Bürger, 2000)—which increases alleles’ initial MAFs and contributions to phenotypic variance (a similar slight underestimation of [equilibrium] genetic variance is seen in Fig. A6(a) of Simons et al., 2018).

The same effect plausibly accounts for most of the other small phenotypic differences between the full model and all alleles simulations. The short-term term phenotypic response is slightly faster in the full than in the all alleles simulations (Appendix 3—figure 2E), because of the greater initial phenotypic variance in the full model. The increase in the 3^rd phenotypic moment is smaller in the full than in the all alleles simulations (Appendix 3—figure 2B), possibly because the greater initial phenotypic variance allows the same shift in mean phenotype to be achieved by smaller relative changes in allele frequency (e.g. Equation 10 in the main text). Lastly, the greater variance and smaller 3^rd central moment in the full model simulations plausibly explain the smaller distance from the optimum during the equilibration phase in the full model relative to the all alleles simulations, seen both directly and in the quasi-static approximation (Appendix 3—figure 2D and C). Small differences notwithstanding, the results of the full model simulations support all of our main results about the phenotypic dynamic. The same is true of our main analytic results about the allele dynamics (Appendix 3—figure 3).

Appendix 3—figure 2

Download asset Open asset

Appendix 3—figure 3

Download asset Open asset

Understanding how stabilizing selection on the trait affects allele dynamics is crucial to characterizing the allelic response after a shift in the optimal trait value. Notably, the equilibrium allele dynamic under stabilizing selection shapes the genetic architecture of the trait (i.e., the joint distribution of phenotypic effect sizes and frequencies of alleles affecting the trait) prior to the shift and this genetic architecture shapes the short term allelic and phenotypic response to the shift. Also, after the rapid, initial response, when the mean phenotype is close to the new optimum but the genetic architecture is out of equilibrium, the longer-term allele dynamics are largely determined by the effects of stabilizing selection. Here we therefore describe the allele dynamics under stabilizing selection, both at and out of equilibrium.

The analysis of the effects of stabilizing selection is greatly aided by it being stationary. In particular, the first two moments of change in the frequency of an allele with effect size $\pm a$ and frequency $x$ in a single generation do not depend on time:

Consequently, we can use the diffusion approximation to calculate the sojourn time: the density of time a derived allele spends at any given frequency $x^{\prime }$ until fixation or loss, having been introduced at some initial frequency $p$ (Ewens, 2004, page 141). Namely,

where

and $\mathrm{Erf}$ is the Gaussian error function. We can also use the diffusion approximation to calculate the probability that an allele fixes, $\pi$ , or goes extinct, $1-\pi$ , and the sojourn time conditional on these outcomes. Namely,