Genetic drift, effective population size and neutral theory

Natural selection is not the only force that changes allele frequencies. In every finite population, the alleles that happen to be passed on are a random sample of those in the parents, and chance alone shifts their frequencies from one generation to the next. This random change is genetic drift. BIOL-108 simulated drift with the Wright–Fisher model and showed that it is stronger in small populations. This lesson makes the theory quantitative: how fast drift removes variation, what the effective size of a real population is, how likely a new mutation is to become fixed, and how Motoo Kimura turned drift into a theory of molecular evolution that remains the starting point for interpreting genome data.

Drift removes heterozygosity

In the Wright–Fisher model, a diploid population of N individuals carries 2N copies of each gene, and each generation is formed by drawing 2N copies at random, with replacement, from the previous one. The allele frequency then follows a binomial process: its expected value stays the same, but its variance after one generation is p(1 − p)/2N. Over time, frequencies wander until one allele is fixed and the other lost. Because drift is random, it does not favor either allele; the probability that an allele eventually becomes fixed equals its current frequency.

Drift steadily erodes genetic variation. The expected heterozygosity, H, the probability that two randomly chosen gene copies differ, declines by a factor of 1 − 1/2N each generation, so after t generations H_t = H_0 (1 − 1/2N)^t ≈ H_0 e^(−t/2N). A population of 50 individuals loses about 1 percent of its heterozygosity per generation; a population of 5,000 loses about 0.01 percent. The time scale of drift is therefore measured in units of 2N generations.

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

def wright_fisher(N, p0=0.5, gens=200, reps=2000):
    p = np.full(reps, p0)
    H = []
    for t in range(gens):
        p = rng.binomial(2*N, p) / (2*N)
        H.append(np.mean(2*p*(1 - p)))
    return np.array(H), p

for N in (25, 100, 1000):
    H, p = wright_fisher(N)
    theory = 0.5*(1 - 1/(2*N))**200
    fixed = np.mean((p == 0) | (p == 1))
    print(f"N = {N:4d}: mean H after 200 gens = {H[-1]:.3f} (theory {theory:.3f}); "
          f"{100*fixed:5.1f}% of replicates fixed or lost")

N =   25: mean H after 200 gens = 0.011 (theory 0.009);  96.9% of replicates fixed or lost
N =  100: mean H after 200 gens = 0.184 (theory 0.183);  45.5% of replicates fixed or lost
N = 1000: mean H after 200 gens = 0.452 (theory 0.452);   0.0% of replicates fixed or lost

The simulations match the theory closely. In populations of 25, most replicates have lost one allele within 200 generations; in populations of 1,000, almost none have. Drift is the reason small, isolated populations, such as those of many endangered species, lose genetic variation even without selection, a theme lesson 12 returns to.

Effective population size

Real populations do not behave like the idealized Wright–Fisher model. Not every individual breeds, the sexes may be unequal in number, family sizes vary and population size changes over time. Population geneticists summarize all these departures in a single number, the effective population size, N_e: the size of an ideal population that would lose variation by drift at the same rate. N_e is almost always smaller than the census size, N, often by a factor of ten or more.

These examples explain why population bottlenecks and founder events have such strong genetic effects. Northern elephant seals were hunted to perhaps 20 to 30 animals in the 1890s; they have recovered to well over 100,000, but their genetic diversity remains extremely low. The cheetah shows a similar signature of ancient bottlenecks. In humans, founder effects explain why certain genetic diseases are unusually common in groups descended from small founding populations, such as several conditions in parts of Quebec, including Saguenay–Lac-Saint-Jean, and among Old Order Amish communities in Pennsylvania.

The fate of new mutations

Every new mutation starts as a single copy, at frequency 1/2N. If it is neutral, its probability of eventual fixation is simply that starting frequency, 1/2N, so nearly all new neutral mutations are lost, but a steady trickle is fixed. A beneficial mutation does better, but not as much better as intuition suggests. For a new advantageous allele with selection coefficient s in a large population, the probability of fixation is approximately 2s (for additive effects, Haldane’s result). A mutation with a 1 percent advantage is lost by chance about 98 percent of the time, usually in the first few generations when it exists in only a few copies.

Kimura derived a general formula for the fixation probability of an allele at initial frequency p: u(p) = (1 − e^(−4N_e s p)) / (1 − e^(−4N_e s)). The key quantity is the product N_e s. When N_e s is much smaller than 1, the allele behaves almost as if neutral; when it is much larger than 1, selection dominates. The same mutation can therefore be effectively neutral in a small population and strongly selected in a large one.

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

def kimura(Ne, s, p):
    if s == 0: return p
    return (1 - np.exp(-4*Ne*s*p)) / (1 - np.exp(-4*Ne*s))

def simulate(N, s, reps=20000):
    p = np.full(reps, 1/(2*N)); alive = np.ones(reps, bool)
    while alive.any():
        q = p[alive]
        q = q*(1 + s) / (q*(1 + s) + (1 - q))          # selection (haploid-style, additive)
        q = rng.binomial(2*N, q) / (2*N)                # drift
        p[alive] = q
        alive = (p > 0) & (p < 1)
    return np.mean(p == 1)

N = 200
for s in (0.0, 0.001, 0.01, 0.05):
    print(f"s = {s:5.3f}: N s = {N*s:5.1f}; fixation, Kimura {kimura(N, s, 1/(2*N)):.4f}, "
          f"simulated {simulate(N, s):.4f}, 2s approximation {2*s:.4f}")

s = 0.000: N s =   0.0; fixation, Kimura 0.0025, simulated 0.0022, 2s approximation 0.0000
s = 0.001: N s =   0.2; fixation, Kimura 0.0036, simulated 0.0028, 2s approximation 0.0020
s = 0.010: N s =   2.0; fixation, Kimura 0.0198, simulated 0.0202, 2s approximation 0.0200
s = 0.050: N s =  10.0; fixation, Kimura 0.0952, simulated 0.0918, 2s approximation 0.1000

With N = 200, a neutral mutation is fixed with probability 1/400 = 0.0025. A mutation with s = 0.001 has N s = 0.2 and is barely more likely to be fixed than a neutral one: drift dominates. With s = 0.05, fixation probability is close to 2s = 0.1, about 40 times the neutral value, but even this strongly favored mutation is lost nine times out of ten.

The neutral theory of molecular evolution

When amino acid sequences of proteins from different species became available in the 1960s, they showed that each protein changes at a roughly constant rate over time, the basis of the molecular clock (BIOL-108). Kimura argued in 1968 that this regularity was hard to explain if most changes were driven by selection, which would depend on changing environments. He proposed instead that most substitutions are neutral or nearly so, and are fixed by drift.

Neutral theory makes a remarkable prediction. The rate of substitution, k, is the number of new mutations per generation, 2Nμ, times the probability that each is fixed, 1/2N. These cancel: k = μ. For neutral mutations, the rate of evolution equals the mutation rate, independent of population size. This explains why the molecular clock ticks steadily in very different species. It also predicts that parts of the genome under less functional constraint should evolve faster, and they do: synonymous sites, which do not change amino acids, evolve faster than non-synonymous sites; pseudogenes evolve fastest of all; and functionally critical proteins such as histones barely change.

Neutral theory does not deny that adaptive changes happen; it claims that they are a minority of molecular changes and that most variation within species is neutral. In 1973, Tomoko Ohta extended the idea to the nearly neutral theory, which recognizes that many mutations are slightly harmful. In small populations, where N_e s is small, these slightly harmful mutations behave as if neutral and can be fixed; in large populations, selection removes them. The nearly neutral theory explains why some molecular clocks run at different rates in species with different population sizes and generation times.

Testing for selection with neutral theory

Because neutral theory makes precise predictions, it serves as a null hypothesis: departures from it are evidence of selection. One widely used comparison is the ratio of non-synonymous to synonymous substitution rates, dN/dS, often written ω. If amino acid changes are neutral, ω ≈ 1. Purifying selection, which removes harmful amino acid changes, makes ω less than 1, and most genes have ω well below 1. Positive selection, which repeatedly favors new amino acids, can make ω greater than 1. Genes with ω above 1 include those in arms races with pathogens, such as antiviral genes and HLA binding sites, and those involved in reproduction.

# The McDonald–Kreitman test compares polymorphism (within species) with divergence (between species)
# Invented counts for a gene compared between two fruit fly species
from scipy.stats import fisher_exact
Pn, Ps = 2, 42     # non-synonymous and synonymous polymorphisms within species
Dn, Ds = 7, 17     # non-synonymous and synonymous fixed differences between species
odds, pval = fisher_exact([[Dn, Ds], [Pn, Ps]])
alpha = 1 - (Ds*Pn)/(Dn*Ps)    # estimated proportion of substitutions fixed by positive selection
print(f"Dn/Ds = {Dn/Ds:.2f}, Pn/Ps = {Pn/Ps:.2f}")
print(f"Fisher exact test p = {pval:.4f}; estimated alpha = {alpha:.2f}")

Dn/Ds = 0.41, Pn/Ps = 0.05
Fisher exact test p = 0.0073; estimated alpha = 0.88

The McDonald–Kreitman test, introduced in 1991 using the Adh gene of fruit flies, compares the ratio of non-synonymous to synonymous changes within a species with the same ratio for fixed differences between species. Under neutrality the two ratios should be equal. In the invented example, amino acid changes make up a much larger share of fixed differences than of polymorphisms, which suggests that many fixed amino acid changes were driven by positive selection; the estimate alpha is the proportion of substitutions that were adaptive. Genome-wide studies using this approach suggest that in fruit flies roughly half of amino acid substitutions were adaptive, while in humans, with a smaller long-term effective size, the proportion is much lower. Methods that use the full distribution of allele frequencies and patterns of linked variation, covered in BIOL-322 and BIOL-410, extend these tests.