Genetic drift and effective population size

Genetic drift is the change in allele frequency caused by the random sampling of gene copies in a finite population. BIOL-321 introduced the Wright–Fisher model, the decay of heterozygosity at rate 1/2N per generation, the standard formulas for effective population size and the result that a neutral allele is fixed with probability equal to its frequency. This lesson takes the mathematics further. We treat drift as a stochastic process with a variance that grows over time, compute exact absorption probabilities and times with a transition matrix, compare the Wright–Fisher and Moran models, distinguish different kinds of effective size, and learn how N_e is estimated from genetic data in real populations, the step that turns theory into a conservation tool.

Drift as a stochastic process

In the Wright–Fisher model, a diploid population of N individuals has 2N gene copies. Each generation, the number of copies of allele A is a binomial draw with 2N trials and success probability equal to the current frequency p. The expected frequency next generation is p, so drift has no direction, but the variance of the change is p(1 − p)/2N. Over many generations these variances accumulate. Starting from p_0, the variance of the frequency among replicate populations after t generations is:

At t = 0 the variance is zero; as t grows it approaches p_0(1 − p_0), the variance of a population that is either fixed for A (with probability p_0) or for a (with probability 1 − p_0). This formula has a second interpretation. If a species is divided into many isolated subpopulations of size N, all founded at frequency p_0, the variance among them grows exactly this way. Lesson 7 will define F_ST as this variance scaled by p_0(1 − p_0), so under pure drift F_ST = 1 − (1 − 1/2N)^t, the same quantity as the loss of heterozygosity within each subpopulation. Drift within populations and differentiation among them are two faces of one process.

Exact results from the transition matrix

Because the number of copies next generation depends only on the number now, the Wright–Fisher model is a Markov chain with 2N + 1 states, from 0 copies to 2N. The entry in row i and column j of its transition matrix is the binomial probability of moving from i to j copies. The states 0 and 2N are absorbing: once an allele is lost or fixed, drift alone cannot bring variation back. With the matrix in hand, we can compute exact probabilities of fixation and expected times to absorption, without simulation and without approximations. For small populations, such as captive breeding groups, these exact results are directly useful.

import numpy as np
from scipy.stats import binom

N = 10; M = 2 * N                                # 20 gene copies
states = np.arange(M + 1)
P = np.array([binom.pmf(states, M, i / M) for i in states])   # transition matrix

# fixation probabilities: iterate the chain for many generations from each start
Pt = np.linalg.matrix_power(P, 2000)
fix = Pt[:, M]
print("start copies  P(fix)  equals p?")
for i in (1, 5, 10, 15):
    print(f"{i:12d}  {fix[i]:.4f}  {abs(fix[i] - i / M) < 1e-6}")

# expected time to absorption from transient states: t = (I - Q)^-1 1
Q = P[1:M, 1:M]
t = np.linalg.solve(np.eye(M - 1) - Q, np.ones(M - 1))
for i in (1, 10):
    p = i / M
    approx = -4 * N * (p * np.log(p) + (1 - p) * np.log(1 - p))   # Kimura-Ohta
    print(f"start at p = {p:.2f}: exact mean time to absorption = {t[i-1]:.1f}, diffusion approx = {approx:.1f}")

start copies  P(fix)  equals p?
           1  0.0500  True
           5  0.2500  True
          10  0.5000  True
          15  0.7500  True
start at p = 0.05: exact mean time to absorption = 7.2, diffusion approx = 7.9
start at p = 0.50: exact mean time to absorption = 26.2, diffusion approx = 27.7

The exact calculation confirms that the fixation probability of a neutral allele equals its starting frequency. It also checks a classic result of Motoo Kimura and Tomoko Ohta, derived from the diffusion approximation, which treats allele frequency as varying continuously. The expected time until an allele at frequency p is either lost or fixed is about −4N[p ln p + (1 − p) ln(1 − p)] generations. Even in a population of only ten individuals the approximation is close.

Times to fixation and loss

Kimura and Ohta also derived the expected time to fixation for alleles that do become fixed. For a new neutral mutation, starting as a single copy, the conditional time to fixation is about 4N generations, while the expected time to loss, for the much larger number of mutations that are lost, is only about 2 ln(2N) generations. The contrast is striking: in a population of 10,000, a new neutral mutation that is lost survives on average about 20 generations, but one that fixes takes on the order of 40,000 generations. Most new mutations flicker briefly and vanish; the few that succeed take a very long time. This explains why most polymorphisms in a sample are young and rare, while the rare old ones that reach intermediate frequency carry information about deep history.

Wright–Fisher and Moran models

The Wright–Fisher model assumes discrete, non-overlapping generations. In the Moran model, often used for organisms with overlapping generations and in the mathematics of the coalescent, one individual is chosen at random to reproduce and one to die at each time step. A Moran population of N haploid individuals experiences drift twice as fast per generation (N time steps) as a Wright–Fisher population of the same size, because the variance in offspring number is twice as large. In other words, a Moran population of size N behaves like a Wright–Fisher population of size N/2. This is a first glimpse of the most important idea in the lesson: what matters for drift is not the census size but the variance in reproductive success.

What effective size really measures

The effective population size, N_e, is the size of an idealized Wright–Fisher population that would show the same amount of drift as the real one. ‘The same amount of drift’ can be measured in different ways, and they do not always agree. The inbreeding effective size matches the rate of increase in inbreeding, or equivalently the loss of heterozygosity, and depends on the number of parents in the previous generation. The variance effective size matches the variance in allele frequency change from one generation to the next and depends on the number of offspring produced. In a stable population the two coincide; in a growing or shrinking population they can differ substantially. The coalescent effective size, used in lesson 6, matches the rate at which gene lineages find common ancestors and is the quantity estimated from DNA sequence diversity.

A general expression for a population of N parents, where V_k is the variance in the number of gene copies each contributes to the next generation, shows the central role of reproductive variance:

When family sizes follow a Poisson distribution with mean 2 (V_k = 2), N_e ≈ N, which is the ideal case. When some individuals have many offspring and most have none, as in species with high fecundity and high early mortality, V_k can be enormous and N_e tiny relative to N. Marine fishes and oysters, with millions of eggs per female and survival governed by chance currents, sometimes show N_e/N ratios below 10⁻³, a phenomenon called sweepstakes reproduction. In the other direction, captive breeding programs that equalize family sizes, so that every pair contributes exactly two offspring (V_k = 0), achieve N_e ≈ 2N, twice the census size. This is why zoo and hatchery managers try to equalize contributions among founders.

Estimating effective size from genetic data

Census counts are not enough, so conservation geneticists estimate N_e directly from genetic data. Three families of methods are widely used.

The temporal method compares allele frequencies in samples taken several generations apart. Because drift causes frequencies to change with variance proportional to t/2N_e, the standardized variance in frequency change, F, estimates t/(2N_e) after correcting for the noise added by sampling S_0 and S_t individuals at each time. The estimate is then N_e = t / [2(F − 1/(2S_0) − 1/(2S_t))]. Museum specimens and archived scales or tissue make this method possible for many species.

The linkage disequilibrium method uses a single sample. In a finite population, drift creates random associations between alleles at different loci, and for unlinked loci the expected squared correlation is about 1/(3N_e) plus a sampling term of about 1/S. Subtracting the sampling term and inverting gives an estimate of N_e for the parental generation. This approach, implemented in the program NeEstimator, is now the most widely used single-sample method. Lesson 9 explains linkage disequilibrium in detail.

Sibship and pedigree methods reconstruct family relationships from multilocus genotypes; the more often random pairs of individuals turn out to be siblings, the smaller the number of breeders that produced them. Finally, coalescent methods use sequence diversity to estimate long-term N_e and, through changes in the site frequency spectrum or in the density of coalescent events across the genome, to reconstruct how N_e has changed over thousands of generations.

import numpy as np
rng = np.random.default_rng(322)

Ne_true, t, L = 150, 6, 400                 # effective size, generations apart, loci
S0 = St = 60                                # individuals sampled at each time
p0 = rng.uniform(0.1, 0.9, L)
p = p0.copy()
for _ in range(t):
    p = rng.binomial(2 * Ne_true, p) / (2 * Ne_true)

x = rng.binomial(2 * S0, p0) / (2 * S0)     # sample frequencies at time 0
y = rng.binomial(2 * St, p) / (2 * St)      # sample frequencies at time t
z = (x + y) / 2
keep = (z > 0.05) & (z < 0.95)
Fc = np.mean((x[keep] - y[keep])**2 / (z[keep] * (1 - z[keep])))
Ne_hat = t / (2 * (Fc - 1 / (2 * S0) - 1 / (2 * St)))
print(f"standardized variance F = {Fc:.4f}")
print(f"temporal estimate of Ne = {Ne_hat:.0f} (true value {Ne_true})")

standardized variance F = 0.0345
temporal estimate of Ne = 168 (true value 150)

The estimate lands near the true value, but the confidence intervals on N_e estimates are typically wide, especially when N_e is large relative to the number of generations and individuals sampled, because the drift signal is then small compared with sampling noise. A general rule of genetic estimation of N_e is that small effective sizes, the ones of greatest conservation concern, are estimated most precisely.

Ratios of effective to census size and conservation rules of thumb

Reviews of hundreds of estimates, beginning with Richard Frankham’s in 1995, found that N_e for wildlife populations is typically around one tenth of census size, with very wide variation among species. Fluctuations in population size, unequal sex ratios and variance in family size each contribute, and they compound one another. A population of 1,000 animals may therefore lose variation like a population of 100.

These numbers feed into the ‘50/500’ rule proposed in 1980: an effective size of at least 50 to avoid inbreeding depression in the short term, and 500 to retain enough additive genetic variance for long-term adaptation. Frankham and colleagues argued in 2014 that these figures are too low and proposed at least 100 and 1,000, respectively, a revision since adopted in some assessments of genetic diversity in conservation policy. Because N_e is typically a fraction of N, an effective size of 1,000 can require several thousand adults. The Kunming–Montreal Global Biodiversity Framework of 2022 includes a headline indicator based on the proportion of populations within species that have N_e above 500, which has made effective-size estimation part of national biodiversity reporting.