Drift, selection and mutation

Three forces in a fixed order

The Combined module applies all three forces to a population of N diploid individuals, in this order every generation. Selection first: pₛ = (p²wAA + pq·wAa)/w̄. Mutation second: p′ = pₛ(1 − μ) + (1 − pₛ)ν. Sampling third: the count of A copies is a draw from Binomial(2N, p′). The first two steps are deterministic and the third is random.

The case in this tutorial is a harmful recessive allele. A is the normal allele with wAA = wAa = 1, and a is harmful with waa = 0.8, so the selection coefficient is s = 0.2. A mutates to a at μ = 0.001, and there is no back mutation (ν = 0). Selection removes a, and mutation re-creates it, so a population can settle where the two flows are equal. That is mutation-selection balance.

The default values of the module are different: wAA = wAa = 1, waa = 0.8, μ = 0.0001 and ν = 0.00001. The deterministic path from p₀ = 0.5 reaches 0.946 at generation 100 and 0.977 at generation 500, which heads for a balance of about 0.978.

Where the balance sits

For a recessive harmful allele the balance is q̂ ≈ √(μ/s). Here q̂ = √(0.001/0.2) = 0.071, so p̂ = 0.929. The square root has a clear cause: selection only sees aa individuals, and they are rare, q² = 0.005 of the population at balance. The a copies hidden in Aa heterozygotes are invisible to selection.

If the heterozygote also pays a cost, selection reaches every carrier. With wAa = 0.9, the heterozygote cost is hs = 0.1, and the balance becomes q̂ ≈ μ/(hs) = 0.001/0.1 = 0.01, seven times lower than 0.071. In the same N = 5000 setting the mean frequency of A is 0.9907 after 300 generations, against a no-drift value of 0.9900.

A recessive harmful allele persists at a much higher frequency than a harmful allele that acts in heterozygotes. The same mutation rate therefore leaves a recessive allele at a higher frequency, 0.071 against 0.010 in this example.

Which force wins

The balance is a prediction for large populations only. Drift is negligible when 4Ns is far above 1, but for a rare recessive a the relevant strength is closer to 4Nsq, because selection only acts on the q² fraction that is aa. At N = 5000 and q = 0.07 that is 4 × 5000 × 0.2 × 0.07 = 280. At N = 50 it is 2.8, and drift begins to dominate.

Selection also slows as it weakens. With waa = 0.99 (s = 0.01) the balance is q̂ = √(0.001/0.01) = 0.316, so p̂ is about 0.684. The path takes much longer: after 1,000 generations at N = 5000 the mean is 0.690 with a standard deviation of 0.049, wider than the 0.011 for s = 0.2. A weaker restoring force leaves more room for drift.