Population Genetics
Allele frequency dynamics — drift, selection, mutation, and the coalescent
A Brief History
Deterministic Selection at a Single Locus
Consider a single locus with two alleles, $A$ and $B$, at frequencies $p$ and $q = 1 - p$ in a large, randomly mating population. Assign fitnesses $w_{AA}$, $w_{AB}$, and $w_{BB}$ to the three genotypes. Under Hardy–Weinberg mating, the mean fitness of the population is:
$$\bar{w} = p^2 w_{AA} + 2pq \, w_{AB} + q^2 w_{BB}$$The change in allele frequency per generation is given by the selection equation:
$$p' = \frac{p^2 w_{AA} + pq \, w_{AB}}{\bar{w}}$$
or equivalently, the change per generation:
$$\Delta p = \frac{pq \left[ p(w_{AA} - w_{AB}) + q(w_{AB} - w_{BB}) \right]}{2\bar{w}}$$
This is zero when $p = 0$, $p = 1$, or when there is an interior equilibrium satisfying $\hat{p}(w_{AA} - w_{AB}) + \hat{q}(w_{AB} - w_{BB}) = 0$.
Three classical cases arise depending on the dominance and fitness relationships:
- Directional selection ($w_{AA} > w_{AB} > w_{BB}$): allele $A$ sweeps to fixation. The trajectory is sigmoidal — slow at low frequency (few $AA$ homozygotes), fastest at intermediate frequency, and slow again near fixation.
- Overdominance ($w_{AB} > w_{AA}, w_{BB}$): the heterozygote is fittest, producing a stable polymorphic equilibrium at $\hat{p} = (w_{AB} - w_{BB}) / (2w_{AB} - w_{AA} - w_{BB})$. Selection maintains both alleles indefinitely.
- Underdominance ($w_{AB} < w_{AA}, w_{BB}$): the heterozygote is least fit, creating an unstable equilibrium. The population is driven to fixation of whichever allele starts above the threshold — a bistable system.
Note the connection to dynamical systems: the selection equation is a one-dimensional discrete map on $[0,1]$, and its fixed points and stability can be analysed exactly as in Flows on the Line.
Genetic Drift: The Wright–Fisher Model
In a finite population of $N$ diploid individuals (so $2N$ gene copies), allele frequencies fluctuate from generation to generation even without selection. The Wright–Fisher model captures this by treating each new generation as a binomial sample from the parental gene pool.
If the current frequency of allele $A$ is $p$, the number of $A$ copies in the next generation, $k$, is drawn from:
$$k \sim \text{Binomial}(2N, p)$$The new frequency is $p' = k / (2N)$. Because the variance of a binomial is $2Np(1-p)$, the variance in allele frequency change per generation is:
$$\text{Var}(\Delta p) = \frac{p(1-p)}{2N}$$This reveals the central tension of population genetics: drift is inversely proportional to population size. In small populations, allele frequencies wander erratically; in large populations, drift becomes negligible and deterministic selection dominates.
Absorption: $p = 0$ and $p = 1$ are absorbing states.
Every neutral allele is eventually lost or fixed.
Fixation probability (neutral): An allele present in
$i$ copies out of $2N$ fixes with probability $i/(2N)$. A single new
mutant fixes with probability $1/(2N)$.
Time to fixation: The expected time to fixation of a
new neutral mutant, conditional on fixation, is approximately $4N$
generations.
Heterozygosity decay: Expected heterozygosity decays
as $H_t = H_0 \left(1 - \frac{1}{2N}\right)^t$.
The Moran Process
The Moran process provides an alternative finite-population model with overlapping generations. At each time step, one individual is chosen to reproduce (proportional to fitness) and one is chosen to die (uniformly at random). The offspring replaces the dead individual, keeping population size constant at $N$.
For two types $A$ and $B$ with fitnesses $f_A$ and $f_B$, if there are currently $i$ copies of $A$, the transition probabilities are:
$$T^+_i = \frac{i \, f_A}{i \, f_A + (N-i) f_B} \cdot \frac{N-i}{N}, \qquad T^-_i = \frac{(N-i) f_B}{i \, f_A + (N-i) f_B} \cdot \frac{i}{N}$$where $T^+_i$ is the probability of gaining one $A$ and $T^-_i$ of losing one. This is a birth–death chain on $\{0, 1, \ldots, N\}$ with absorbing barriers at $0$ and $N$.
For a single mutant $A$ with relative fitness $r = f_A / f_B$:
$$\rho_A = \frac{1 - 1/r}{1 - 1/r^N}$$
For neutral drift ($r = 1$): $\rho = 1/N$. For advantageous mutants ($r > 1$) in large populations: $\rho \approx 1 - 1/r$. For deleterious mutants ($r < 1$): $\rho \to 0$ exponentially in $N$.
The Moran process connects directly to the evolutionary game theory of later chapters: when fitness depends on the composition of the population (frequency-dependent selection), the Moran process becomes the basis for stochastic evolutionary dynamics on finite populations.
The Neutral Theory
Kimura’s neutral theory of molecular evolution (1968) proposed that the vast majority of mutations observed at the molecular level are selectively neutral — their fate is governed by drift alone. This was controversial because it challenged the prevailing view that most evolutionary change is adaptive.
The neutral theory makes a striking quantitative prediction. If neutral mutations arise at rate $\mu$ per gene per generation in a diploid population of size $N$, then $2N\mu$ new mutations enter per generation. Each fixes with probability $1/(2N)$, so the rate of neutral substitution is:
$$k = 2N\mu \cdot \frac{1}{2N} = \mu$$The substitution rate equals the mutation rate, independent of population size. This elegant result explains the approximate constancy of the molecular clock: proteins evolve at roughly constant rates across lineages, regardless of differences in population size, generation time, or ecology — as long as the fraction of neutral sites remains similar.
The neutral theory does not claim that all mutations are neutral. It claims that among mutations that reach detectable frequency or fix, the overwhelming majority are neutral. Strongly deleterious mutations are removed by purifying selection before they can drift to appreciable frequency; strongly advantageous mutations are rare. What remains — and what shapes most molecular variation — is drift acting on neutral (or nearly neutral) variants.
The Diffusion Approximation
For large populations, tracking individual allele counts becomes unwieldy. The diffusion approximation replaces the discrete Wright–Fisher or Moran dynamics with a continuous stochastic differential equation. If $p$ is the allele frequency, the change per generation is approximated by:
$$dp = M(p)\,dt + \sqrt{V(p)}\,dW$$where $M(p)$ is the deterministic drift (selection, mutation), $V(p)$ is the variance due to genetic drift, and $dW$ is a Wiener process increment. For the Wright–Fisher model with selection coefficient $s$ favouring allele $A$:
$$M(p) = sp(1-p), \qquad V(p) = \frac{p(1-p)}{2N_e}$$
The ratio $\theta = 4N_e s$ determines the relative strength of selection versus drift. When $|\theta| \gg 1$, selection dominates; when $|\theta| \ll 1$, drift dominates; the transition occurs around $|s| \approx 1/(2N_e)$.
From the diffusion equation, Kimura (1962) derived the fixation probability of a new mutant with selective advantage $s$:
$$u(p) = \frac{1 - e^{-4N_e s p}}{1 - e^{-4N_e s}}$$For a single new mutant ($p = 1/(2N)$) with $N_e = N$: $u \approx 2s$ when $s \ll 1$ and $N$ is large — Haldane’s classical result from 1927. For neutral alleles ($s = 0$): $u = p = 1/(2N)$, recovering the Wright–Fisher result.
The Coalescent
The coalescent (Kingman, 1982) reverses the perspective: instead of following allele frequencies forward in time, it traces a sample of $n$ gene copies backward to their most recent common ancestor (MRCA). This backward view is enormously powerful because it conditions on the observed sample, ignoring the vast number of lineages that left no descendants.
In a Wright–Fisher population of $2N$ gene copies, consider $n$ sampled lineages. Looking one generation back, any two lineages coalesce (share a parent) with probability $1/(2N)$. When $n \ll 2N$, the probability that exactly one coalescent event occurs is approximately:
$$\binom{n}{2} \cdot \frac{1}{2N} = \frac{n(n-1)}{4N}$$The waiting time until the next coalescence, when there are $k$ lineages, is approximately exponential with rate $\binom{k}{2}/(2N)$. The total time from $n$ lineages to the MRCA (2 lineages coalescing to 1) is:
$$\mathbb{E}[T_{\text{MRCA}}] = 4N\left(1 - \frac{1}{n}\right)$$which approaches $4N$ for large samples. Half of the total tree height is contributed by the final coalescence of two lineages — a striking and counterintuitive result.
The genealogy of a sample is independent of the mutational process. First, generate the coalescent tree (which depends only on $N$); then, sprinkle mutations onto the branches (which depends only on $\mu$). This separation makes the coalescent computationally and analytically tractable, and is the basis of modern tools like BEAST, ms, and msprime.
Explore: Drift vs. Selection
Exercises
Conceptual Questions
- Explain why the fixation probability of a neutral allele equals its initial frequency. What assumption about reproduction makes this result exact?
- The neutral theory predicts that the rate of molecular substitution equals the mutation rate $\mu$, independent of population size. Derive this result and explain its biological significance for the molecular clock.
- Under overdominance ($w_{AB} > w_{AA}, w_{BB}$), selection maintains a stable polymorphism in an infinite population. What happens in a finite population? Is the polymorphism truly stable?
- In the coalescent, the expected time to the most recent common ancestor of a sample of $n$ lineages is approximately $4N(1 - 1/n)$. Explain why this is dominated by the final coalescence of two lineages and what this implies about the shape of gene genealogies.
- Discuss the relationship between the effective population size $N_e$ and the census population size $N$. Give three biological factors that cause $N_e < N$ and explain their effects on the rate of genetic drift.
Computer Problems
- Wright–Fisher Simulation. Implement the Wright–Fisher model with binomial sampling. Simulate 1000 replicate populations with $N = 50$, $p_0 = 0.5$, $s = 0$ for 200 generations. Plot the mean and variance of $p$ over time and compare to the theoretical predictions $\mathbb{E}[p] = p_0$ and $\text{Var}(p_t) = p_0(1-p_0)\left[1 - (1-1/(2N))^t\right]$.
- Fixation Probability. For the Moran process with $N = 50$ and a single mutant ($i=1$), estimate the fixation probability by running 10,000 simulations for $r = 0.5, 0.8, 1.0, 1.2, 1.5, 2.0$. Plot the empirical fixation probability against $\rho = (1 - 1/r)/(1 - 1/r^N)$ and assess the agreement.
- Drift vs. Selection. Simulate Wright–Fisher trajectories with $p_0 = 0.1$, $N = 100$, and $s = 0, 0.01, 0.05, 0.1$. For each, run 500 replicates and compute the fraction that fix. Plot fixation probability versus $s$ and compare to the diffusion approximation $u \approx (1 - e^{-4Nsp_0})/(1 - e^{-4Ns})$.
- Coalescent Simulation. Implement Kingman’s coalescent for a sample of $n = 10$ from a population of $N = 1000$. Draw 5000 coalescent trees, record the time to MRCA, and compare the distribution to the theoretical expectation $\mathbb{E}[T_{\text{MRCA}}] = 4N(1 - 1/n)$.
- Molecular Clock Test. Using the neutral substitution rate $k = \mu$, simulate sequence evolution on a star phylogeny with 5 taxa diverging $T$ generations ago. Vary $N$ (50, 500, 5000) while keeping $\mu = 10^{-4}$. Show that the number of substitutions per lineage is approximately $\mu T$ regardless of $N$.
References
- Hardy, G. H. (1908). Mendelian proportions in a mixed population. Science, 28, 49–50.
- Fisher, R. A. (1930). The Genetical Theory of Natural Selection. Clarendon Press, Oxford.
- Wright, S. (1931). Evolution in Mendelian populations. Genetics, 16, 97–159.
- Haldane, J. B. S. (1932). The Causes of Evolution. Longmans, Green & Co., London.
- Moran, P. A. P. (1962). The Statistical Processes of Evolutionary Theory. Clarendon Press, Oxford.
- Kimura, M. (1968). Evolutionary rate at the molecular level. Nature, 217, 624–626.
- Kimura, M. (1962). On the probability of fixation of mutant genes in a population. Genetics, 47, 713–719.
- Kingman, J. F. C. (1982). The coalescent. Stochastic Processes and their Applications, 13, 235–248.
- Ewens, W. J. (2004). Mathematical Population Genetics. 2nd ed. Springer.
- Nowak, M. A. (2006). Evolutionary Dynamics: Exploring the Equations of Life. Harvard University Press.