The Moran process and fixation probability

evolutionary ecology

population dynamics

natural selection

ecology tutorial

Model selection in a finite population in base R: the Moran birth-death process, the exact fixation probability of a single mutant, and how long fixation takes.

Author

Tidy Ecology

Published

2026-08-10

Every model in the last few posts assumed an infinite population. The replicator equation tracks frequencies, frequencies change smoothly, and a strategy with higher payoff always increases. Real populations are finite, and in a finite population the fitter type sometimes disappears anyway, because the individual carrying it happened to die before reproducing.

That is not a small correction. In a population of 100, a mutant with a ten per cent fitness advantage fixes about nine times in a hundred and vanishes the other ninety-one. Selection is still doing its work: nine per cent is nine times the neutral rate. But the modal outcome for a beneficial mutation is loss, and no deterministic model can say that.

The Moran process is the smallest model that gets this right. It has one parameter beyond the population size, it can be written in ten lines, and its fixation probability has a closed form that you can check against simulation.

The model

Fix the population size at N individuals, each of one of two types: A with fitness r, and B with fitness 1. One elementary step does two things at once. An individual is chosen to reproduce, with probability proportional to fitness, and its offspring replaces an individual chosen uniformly at random to die. The population size never changes, and generations overlap completely.

If there are i copies of A, the total fitness is ri + (N - i). The count goes up when an A reproduces and a B dies, and down when a B reproduces and an A dies:

\[p_{i,i+1} = \frac{ri}{ri + N - i} \cdot \frac{N-i}{N}, \qquad p_{i,i-1} = \frac{N-i}{ri + N - i} \cdot \frac{i}{N}\]

Everything else is the probability of staying put. States 0 and N are absorbing: once a type is gone it cannot come back, because there is no mutation.

That structure makes the fixation probability a standard birth-death calculation. Writing $\gamma_i = p_{i,i-1}/p_{i,i+1} = 1/r$, which does not depend on i here, the probability that a single A takes over is

\[\rho_1 = \frac{1 - r^{-1}}{1 - r^{-N}}\]

and the general starting point i gives $(1 - r^{-i})/(1 - r^{-N})$. At r exactly 1 the expression is 0/0 and the limit is i/N, which is the neutral result: a mutant fixes with the frequency it starts at.

rho <- function(r, N, i = 1) if (r == 1) i / N else (1 - r^-i) / (1 - r^-N)

for (r in c(0.90, 0.99, 1.00, 1.01, 1.10, 1.50, 2.00)) {
  cat(sprintf("r = %.2f   rho_1 = %-11.6g   1 - 1/r = %s\n",
              r, rho(r, 100),
              if (r > 1) sprintf("%.6g", 1 - 1 / r) else "not defined"))
}

r = 0.90   rho_1 = 2.95134e-06   1 - 1/r = not defined
r = 0.99   rho_1 = 0.005832      1 - 1/r = not defined
r = 1.00   rho_1 = 0.01          1 - 1/r = not defined
r = 1.01   rho_1 = 0.0157087     1 - 1/r = 0.00990099
r = 1.10   rho_1 = 0.0909157     1 - 1/r = 0.0909091
r = 1.50   rho_1 = 0.333333      1 - 1/r = 0.333333
r = 2.00   rho_1 = 0.5           1 - 1/r = 0.5

Three things in that table are worth stopping on.

A ten per cent advantage buys a fixation probability of 0.0909157 in a population of 100, against a neutral 0.01. Selection multiplied the mutant’s chances by nine, and still lost nine times out of ten.

A ten per cent disadvantage does not give zero. It gives 2.95134e-06, which is small but not impossible, and in a large population of such mutants it happens. Deleterious variants do fix.

For r above about 1.1 the exact answer and the large-population limit 1 - 1/r agree to four decimal places at N = 100. That limit is the number worth memorising: a mutant with a ten per cent advantage fixes about a tenth of the time, whatever the population size, once the population is not tiny. The agreement breaks down for weak selection, where r = 1.01 gives 0.0157087 against a limit of 0.00990099, because the relevant quantity is N(r - 1) rather than the selection coefficient on its own, and here that product is only 1.

A simulation check

The formula is short enough to be wrong in a way that looks right, so run the process itself. Each replicate walks until it hits 0 or N, and the fraction ending at N estimates $\rho_1$. The binomial standard error over 2000 replicates is about 0.0064 near a probability of 0.09, so agreement to within a couple of standard errors is all that can be asked.

sim_moran <- function(r, N, i0 = 1, reps = 2000, seed = 1) {
  set.seed(seed)
  hits <- 0L
  for (m in seq_len(reps)) {
    i <- i0
    while (i > 0 && i < N) {
      fa <- r * i; fb <- N - i; tot <- fa + fb
      pup <- (fa / tot) * ((N - i) / N)
      pdn <- (fb / tot) * (i / N)
      u <- runif(1)
      if (u < pup) i <- i + 1L else if (u < pup + pdn) i <- i - 1L
    }
    if (i == N) hits <- hits + 1L
  }
  hits / reps
}

cfg <- data.frame(r = c(1.1, 1.0), seed = c(1, 2))
for (k in seq_len(nrow(cfg))) {
  rr <- cfg$r[k]
  p  <- sim_moran(rr, 100, reps = 2000, seed = cfg$seed[k])
  ex <- rho(rr, 100)
  err <- sqrt(ex * (1 - ex) / 2000)
  cat(sprintf("r = %.2f   simulated %.4f   exact %.7f   se %.5f   gap %.2f se\n",
              rr, p, ex, err, abs(p - ex) / err))
}

r = 1.10   simulated 0.0850   exact 0.0909157   se 0.00643   gap 0.92 se
r = 1.00   simulated 0.0120   exact 0.0100000   se 0.00222   gap 0.90 se

Both land inside one standard error. That is the check that matters, and it also shows why simulating this model to answer a fine question is a poor use of a processor: distinguishing r = 1.01 from neutrality means separating 0.0157 from 0.0100, and at 2000 replicates the standard error is larger than the difference.

The shape of the curve

library(ggplot2)
rs <- seq(0.8, 2, by = 0.01)
curves <- do.call(rbind, lapply(c(10, 100, 1000), function(N)
  data.frame(r = rs, N = factor(N), p = sapply(rs, rho, N = N))))
lim <- data.frame(r = rs[rs > 1], p = 1 - 1 / rs[rs > 1])

ggplot(curves, aes(r, p, colour = N)) +
  geom_line(linewidth = 1.1) +
  geom_line(data = lim, aes(r, p), inherit.aes = FALSE,
            colour = "#3a3a33", linetype = "22", linewidth = 0.8) +
  scale_colour_manual(values = c(`10` = "#c98a3a", `100` = "#275139",
                                 `1000` = "#b5534e")) +
  labs(x = "Relative fitness of the mutant (r)", y = "Fixation probability",
       colour = "N") +
  theme_minimal(base_size = 12) +
  theme(legend.position = "top", panel.grid.minor = element_blank())

Three rising curves of fixation probability against relative fitness; the two larger population sizes coincide with a dashed reference curve while the smallest sits well above it — Figure 1: Fixation probability of a single mutant against its relative fitness, for three population sizes. The dashed curve is the large-population limit 1 - 1/r, which the N = 100 and N = 1000 curves sit on top of above about r = 1.1.

Below r = 1 the curves separate sharply and the large population is the harshest: a deleterious mutant in a population of 1000 is essentially certain to be lost, while in a population of 10 it has a real chance. Above r = 1 the population size stops mattering almost immediately. This asymmetry is the whole reason effective population size appears in every discussion of molecular evolution. Big populations are efficient at removing bad variants and no better than small ones at spreading good variants.

The N = 10 curve staying above the others on the left is drift protecting deleterious mutations, and it is the same effect that lets nearly neutral theory work: whether a mutation behaves as neutral depends on Ns, not on s.

Trajectories

The averages hide what individual runs look like. Twelve replicates starting from ten copies, with a ten per cent advantage:

traj <- function(r, N, i0, seed, maxstep = 60000) {
  set.seed(seed)
  out <- integer(maxstep); i <- i0; k <- 0L
  while (k < maxstep) {
    k <- k + 1L; out[k] <- i
    if (i == 0 || i == N) break
    fa <- r * i; fb <- N - i; tot <- fa + fb
    pup <- (fa / tot) * ((N - i) / N)
    pdn <- (fb / tot) * (i / N)
    u <- runif(1)
    if (u < pup) i <- i + 1L else if (u < pup + pdn) i <- i - 1L
  }
  out[seq_len(k)]
}

runs <- do.call(rbind, lapply(1:12, function(s) {
  v <- traj(1.1, 100, 10, seed = s)
  data.frame(run = factor(s), step = seq_along(v), i = v,
             fate = if (v[length(v)] == 100) "fixed" else "lost")
}))
cat("fixed:", length(unique(runs$run[runs$fate == "fixed"])), "of 12",
    " exact probability from i = 10:", sprintf("%.4f", rho(1.1, 100, 10)), "\n")

fixed: 8 of 12  exact probability from i = 10: 0.6145

ggplot(runs, aes(step, i, group = run, colour = fate)) +
  geom_step(linewidth = 0.5, alpha = 0.85) +
  scale_colour_manual(values = c(fixed = "#275139", lost = "#b5534e")) +
  labs(x = "Elementary steps", y = "Copies of the mutant", colour = NULL) +
  theme_minimal(base_size = 12) +
  theme(legend.position = "top", panel.grid.minor = element_blank())

Twelve jagged random walk paths of mutant count against time; four fall to zero early while eight climb and reach the population size — Figure 2: Twelve independent runs of the Moran process at r = 1.1 in a population of 100, each starting from ten copies of the mutant. Runs are coloured by their outcome.

Eight of the twelve fixed, against an exact probability of 0.6145 from ten copies. Starting from ten rather than one has changed the problem: a mutant that reaches ten copies is already past the dangerous phase, and the fixation probability jumped from 0.09 to 0.61.

The lost runs are short and the surviving ones are long, which is the usual pattern and the reason a plot of “typical” trajectories from a fixed time budget is misleading. Failure is quick.

How long does it take?

Fixation probability answers whether, not when. The expected number of steps to absorption from state i satisfies

\[t_i = 1 + p_{i,i-1}\,t_{i-1} + p_{i,i+1}\,t_{i+1} + (1 - p_{i,i-1} - p_{i,i+1})\,t_i\]

with $t_0 = t_N = 0$. Rearranged, that is a tridiagonal linear system in the N - 1 interior states, which solve handles directly. No simulation, no approximation.

abs_time <- function(r, N) {
  pup <- function(i) (r * i) / (r * i + N - i) * (N - i) / N
  pdn <- function(i) (N - i) / (r * i + N - i) * i / N
  A <- matrix(0, N - 1, N - 1); b <- rep(-1, N - 1)
  for (i in 1:(N - 1)) {
    A[i, i] <- -(pup(i) + pdn(i))
    if (i > 1)     A[i, i - 1] <- pdn(i)
    if (i < N - 1) A[i, i + 1] <- pup(i)
  }
  solve(A, b)
}

t_neu <- abs_time(1.0, 100)
t_sel <- abs_time(1.1, 100)
cat(sprintf("neutral   from i = 1: %9.4f steps    from i = 50: %9.4f\n",
            t_neu[1], t_neu[50]))

neutral   from i = 1:  517.7378 steps    from i = 50: 6881.7218

cat(sprintf("r = 1.1   from i = 1: %9.4f steps    from i = 50: %9.4f\n",
            t_sel[1], t_sel[50]))

r = 1.1   from i = 1:  738.2261 steps    from i = 50: 3307.7862

cat(sprintf("diffusion approximation for the neutral midpoint, N^2 log 2 = %.4f\n",
            100^2 * log(2)))

diffusion approximation for the neutral midpoint, N^2 log 2 = 6931.4718

times <- rbind(
  data.frame(i = 1:99, steps = t_neu, model = "neutral"),
  data.frame(i = 1:99, steps = t_sel, model = "r = 1.1"))
ggplot(times, aes(i, steps, colour = model)) +
  geom_line(linewidth = 1.1) +
  scale_colour_manual(values = c(neutral = "#8d8d80", "r = 1.1" = "#275139")) +
  labs(x = "Starting number of mutants", y = "Expected steps to absorption",
       colour = NULL) +
  theme_minimal(base_size = 12) +
  theme(legend.position = "top", panel.grid.minor = element_blank())

Two humped curves of expected absorption time against starting mutant count; the neutral curve peaks near the middle and is roughly twice the height of the selected curve — Figure 3: Expected number of elementary steps to absorption against the starting number of mutants, in a population of 100, for a neutral mutant and one with a ten per cent advantage.

The neutral midpoint takes 6881.7218 steps against a diffusion prediction of 6931.4718, an agreement of better than one per cent from a formula that assumes N is large and 100 is not especially large. Divide by N to get generations in the usual sense and the neutral midpoint is about 69 generations, which for most organisms is nothing.

Selection cuts the midpoint time roughly in half, to 3307.7862. Starting from a single copy the ordering reverses: the selected mutant takes 738.2261 steps against the neutral 517.7378, because most neutral runs are lost within a few steps while the selected ones that survive have further to travel. Conditioning changes the answer, and unconditional expected times mix two very different populations of runs.

An honest limit

The Moran process buys its closed form with assumptions that are all worth naming.

Population size is fixed. Nothing about a real population holds N constant to the individual, and the fixation probability depends on N only through r^-N, so the result is not sensitive to the exact value; but a population that crashes and recovers behaves like the smaller size for drift purposes, which is what effective population size is trying to capture. Feeding a census count into these formulas will overstate the efficiency of selection, often by an order of magnitude.

Generations overlap completely, one birth and one death at a time. The Wright-Fisher model makes the opposite choice, replacing the whole population at once, and gives $\rho_1 \approx (1 - e^{-2s})/(1 - e^{-2Ns})$ with s = r - 1, differing from the Moran answer by a factor of two in the exponent. Neither is more correct in general. The factor of two is a modelling choice about life history that is easy to lose track of when comparing published numbers.

Fitness is constant. The whole apparatus assumes r does not depend on how common the type is, which rules out every frequency-dependent case, and those are the interesting ones in behaviour and in competition. That extension is the subject of the next post.

There is no mutation. Absorbing states are absorbing only because nothing regenerates the lost type. With recurrent mutation the process has a stationary distribution instead, and “fixation probability” becomes a statement about the rate of transitions between mostly-A and mostly-B states rather than a final outcome.

The population is well mixed, every individual competing with every other. The lattice posts showed how much structure can change a conclusion, and structure changes fixation probabilities too.

--- title: "The Moran process and fixation probability" description: "Model selection in a finite population in base R: the Moran birth-death process, the exact fixation probability of a single mutant, and how long fixation takes." date: "2026-08-10 09:00" categories: [evolutionary ecology, population dynamics, natural selection, R, ecology tutorial] image: thumbnail.png image-alt: "Fixation probability rising with mutant fitness for four population sizes, with the large-population limit drawn as a dashed curve" --- Every model in the last few posts assumed an infinite population. The replicator equation tracks frequencies, frequencies change smoothly, and a strategy with higher payoff always increases. Real populations are finite, and in a finite population the fitter type sometimes disappears anyway, because the individual carrying it happened to die before reproducing. That is not a small correction. In a population of 100, a mutant with a ten per cent fitness advantage fixes about nine times in a hundred and vanishes the other ninety-one. Selection is still doing its work: nine per cent is nine times the neutral rate. But the modal outcome for a beneficial mutation is loss, and no deterministic model can say that. The Moran process is the smallest model that gets this right. It has one parameter beyond the population size, it can be written in ten lines, and its fixation probability has a closed form that you can check against simulation. ## The model Fix the population size at `N` individuals, each of one of two types: `A` with fitness `r`, and `B` with fitness 1. One elementary step does two things at once. An individual is chosen to reproduce, with probability proportional to fitness, and its offspring replaces an individual chosen uniformly at random to die. The population size never changes, and generations overlap completely. If there are `i` copies of `A`, the total fitness is `ri + (N - i)`. The count goes up when an `A` reproduces and a `B` dies, and down when a `B` reproduces and an `A` dies: $$p_{i,i+1} = \frac{ri}{ri + N - i} \cdot \frac{N-i}{N}, \qquad p_{i,i-1} = \frac{N-i}{ri + N - i} \cdot \frac{i}{N}$$ Everything else is the probability of staying put. States 0 and `N` are absorbing: once a type is gone it cannot come back, because there is no mutation. That structure makes the fixation probability a standard birth-death calculation. Writing $\gamma_i = p_{i,i-1}/p_{i,i+1} = 1/r$, which does not depend on `i` here, the probability that a single `A` takes over is $$\rho_1 = \frac{1 - r^{-1}}{1 - r^{-N}}$$ and the general starting point `i` gives $(1 - r^{-i})/(1 - r^{-N})$. At `r` exactly 1 the expression is 0/0 and the limit is `i/N`, which is the neutral result: a mutant fixes with the frequency it starts at. ```{r} #| label: exact rho <- function(r, N, i = 1) if (r == 1) i / N else (1 - r^-i) / (1 - r^-N) for (r in c(0.90, 0.99, 1.00, 1.01, 1.10, 1.50, 2.00)) { cat(sprintf("r = %.2f rho_1 = %-11.6g 1 - 1/r = %s\n", r, rho(r, 100), if (r > 1) sprintf("%.6g", 1 - 1 / r) else "not defined")) } ``` Three things in that table are worth stopping on. A ten per cent advantage buys a fixation probability of 0.0909157 in a population of 100, against a neutral 0.01. Selection multiplied the mutant's chances by nine, and still lost nine times out of ten. A ten per cent disadvantage does not give zero. It gives 2.95134e-06, which is small but not impossible, and in a large population of such mutants it happens. Deleterious variants do fix. For `r` above about 1.1 the exact answer and the large-population limit `1 - 1/r` agree to four decimal places at `N` = 100. That limit is the number worth memorising: a mutant with a ten per cent advantage fixes about a tenth of the time, whatever the population size, once the population is not tiny. The agreement breaks down for weak selection, where `r` = 1.01 gives 0.0157087 against a limit of 0.00990099, because the relevant quantity is `N(r - 1)` rather than the selection coefficient on its own, and here that product is only 1. ## A simulation check The formula is short enough to be wrong in a way that looks right, so run the process itself. Each replicate walks until it hits 0 or `N`, and the fraction ending at `N` estimates $\rho_1$. The binomial standard error over 2000 replicates is about 0.0064 near a probability of 0.09, so agreement to within a couple of standard errors is all that can be asked. ```{r} #| label: simulate sim_moran <- function(r, N, i0 = 1, reps = 2000, seed = 1) { set.seed(seed) hits <- 0L for (m in seq_len(reps)) { i <- i0 while (i > 0 && i < N) { fa <- r * i; fb <- N - i; tot <- fa + fb pup <- (fa / tot) * ((N - i) / N) pdn <- (fb / tot) * (i / N) u <- runif(1) if (u < pup) i <- i + 1L else if (u < pup + pdn) i <- i - 1L } if (i == N) hits <- hits + 1L } hits / reps } cfg <- data.frame(r = c(1.1, 1.0), seed = c(1, 2)) for (k in seq_len(nrow(cfg))) { rr <- cfg$r[k] p <- sim_moran(rr, 100, reps = 2000, seed = cfg$seed[k]) ex <- rho(rr, 100) err <- sqrt(ex * (1 - ex) / 2000) cat(sprintf("r = %.2f simulated %.4f exact %.7f se %.5f gap %.2f se\n", rr, p, ex, err, abs(p - ex) / err)) } ``` Both land inside one standard error. That is the check that matters, and it also shows why simulating this model to answer a fine question is a poor use of a processor: distinguishing `r` = 1.01 from neutrality means separating 0.0157 from 0.0100, and at 2000 replicates the standard error is larger than the difference. ## The shape of the curve ```{r} #| label: fig-rho #| fig-cap: "Fixation probability of a single mutant against its relative fitness, for three population sizes. The dashed curve is the large-population limit 1 - 1/r, which the N = 100 and N = 1000 curves sit on top of above about r = 1.1." #| fig-alt: "Three rising curves of fixation probability against relative fitness; the two larger population sizes coincide with a dashed reference curve while the smallest sits well above it" library(ggplot2) rs <- seq(0.8, 2, by = 0.01) curves <- do.call(rbind, lapply(c(10, 100, 1000), function(N) data.frame(r = rs, N = factor(N), p = sapply(rs, rho, N = N)))) lim <- data.frame(r = rs[rs > 1], p = 1 - 1 / rs[rs > 1]) ggplot(curves, aes(r, p, colour = N)) + geom_line(linewidth = 1.1) + geom_line(data = lim, aes(r, p), inherit.aes = FALSE, colour = "#3a3a33", linetype = "22", linewidth = 0.8) + scale_colour_manual(values = c(`10` = "#c98a3a", `100` = "#275139", `1000` = "#b5534e")) + labs(x = "Relative fitness of the mutant (r)", y = "Fixation probability", colour = "N") + theme_minimal(base_size = 12) + theme(legend.position = "top", panel.grid.minor = element_blank()) ``` Below `r` = 1 the curves separate sharply and the large population is the harshest: a deleterious mutant in a population of 1000 is essentially certain to be lost, while in a population of 10 it has a real chance. Above `r` = 1 the population size stops mattering almost immediately. This asymmetry is the whole reason effective population size appears in every discussion of molecular evolution. Big populations are efficient at removing bad variants and no better than small ones at spreading good variants. The `N` = 10 curve staying above the others on the left is drift protecting deleterious mutations, and it is the same effect that lets nearly neutral theory work: whether a mutation behaves as neutral depends on `Ns`, not on `s`. ## Trajectories The averages hide what individual runs look like. Twelve replicates starting from ten copies, with a ten per cent advantage: ```{r} #| label: fig-paths #| fig-cap: "Twelve independent runs of the Moran process at r = 1.1 in a population of 100, each starting from ten copies of the mutant. Runs are coloured by their outcome." #| fig-alt: "Twelve jagged random walk paths of mutant count against time; four fall to zero early while eight climb and reach the population size" traj <- function(r, N, i0, seed, maxstep = 60000) { set.seed(seed) out <- integer(maxstep); i <- i0; k <- 0L while (k < maxstep) { k <- k + 1L; out[k] <- i if (i == 0 || i == N) break fa <- r * i; fb <- N - i; tot <- fa + fb pup <- (fa / tot) * ((N - i) / N) pdn <- (fb / tot) * (i / N) u <- runif(1) if (u < pup) i <- i + 1L else if (u < pup + pdn) i <- i - 1L } out[seq_len(k)] } runs <- do.call(rbind, lapply(1:12, function(s) { v <- traj(1.1, 100, 10, seed = s) data.frame(run = factor(s), step = seq_along(v), i = v, fate = if (v[length(v)] == 100) "fixed" else "lost") })) cat("fixed:", length(unique(runs$run[runs$fate == "fixed"])), "of 12", " exact probability from i = 10:", sprintf("%.4f", rho(1.1, 100, 10)), "\n") ggplot(runs, aes(step, i, group = run, colour = fate)) + geom_step(linewidth = 0.5, alpha = 0.85) + scale_colour_manual(values = c(fixed = "#275139", lost = "#b5534e")) + labs(x = "Elementary steps", y = "Copies of the mutant", colour = NULL) + theme_minimal(base_size = 12) + theme(legend.position = "top", panel.grid.minor = element_blank()) ``` Eight of the twelve fixed, against an exact probability of 0.6145 from ten copies. Starting from ten rather than one has changed the problem: a mutant that reaches ten copies is already past the dangerous phase, and the fixation probability jumped from 0.09 to 0.61. The lost runs are short and the surviving ones are long, which is the usual pattern and the reason a plot of "typical" trajectories from a fixed time budget is misleading. Failure is quick. ## How long does it take? Fixation probability answers whether, not when. The expected number of steps to absorption from state `i` satisfies $$t_i = 1 + p_{i,i-1}\,t_{i-1} + p_{i,i+1}\,t_{i+1} + (1 - p_{i,i-1} - p_{i,i+1})\,t_i$$ with $t_0 = t_N = 0$. Rearranged, that is a tridiagonal linear system in the `N - 1` interior states, which `solve` handles directly. No simulation, no approximation. ```{r} #| label: fig-time #| fig-cap: "Expected number of elementary steps to absorption against the starting number of mutants, in a population of 100, for a neutral mutant and one with a ten per cent advantage." #| fig-alt: "Two humped curves of expected absorption time against starting mutant count; the neutral curve peaks near the middle and is roughly twice the height of the selected curve" abs_time <- function(r, N) { pup <- function(i) (r * i) / (r * i + N - i) * (N - i) / N pdn <- function(i) (N - i) / (r * i + N - i) * i / N A <- matrix(0, N - 1, N - 1); b <- rep(-1, N - 1) for (i in 1:(N - 1)) { A[i, i] <- -(pup(i) + pdn(i)) if (i > 1) A[i, i - 1] <- pdn(i) if (i < N - 1) A[i, i + 1] <- pup(i) } solve(A, b) } t_neu <- abs_time(1.0, 100) t_sel <- abs_time(1.1, 100) cat(sprintf("neutral from i = 1: %9.4f steps from i = 50: %9.4f\n", t_neu[1], t_neu[50])) cat(sprintf("r = 1.1 from i = 1: %9.4f steps from i = 50: %9.4f\n", t_sel[1], t_sel[50])) cat(sprintf("diffusion approximation for the neutral midpoint, N^2 log 2 = %.4f\n", 100^2 * log(2))) times <- rbind( data.frame(i = 1:99, steps = t_neu, model = "neutral"), data.frame(i = 1:99, steps = t_sel, model = "r = 1.1")) ggplot(times, aes(i, steps, colour = model)) + geom_line(linewidth = 1.1) + scale_colour_manual(values = c(neutral = "#8d8d80", "r = 1.1" = "#275139")) + labs(x = "Starting number of mutants", y = "Expected steps to absorption", colour = NULL) + theme_minimal(base_size = 12) + theme(legend.position = "top", panel.grid.minor = element_blank()) ``` The neutral midpoint takes 6881.7218 steps against a diffusion prediction of 6931.4718, an agreement of better than one per cent from a formula that assumes `N` is large and 100 is not especially large. Divide by `N` to get generations in the usual sense and the neutral midpoint is about 69 generations, which for most organisms is nothing. Selection cuts the midpoint time roughly in half, to 3307.7862. Starting from a single copy the ordering reverses: the selected mutant takes 738.2261 steps against the neutral 517.7378, because most neutral runs are lost within a few steps while the selected ones that survive have further to travel. Conditioning changes the answer, and unconditional expected times mix two very different populations of runs. ## An honest limit The Moran process buys its closed form with assumptions that are all worth naming. Population size is fixed. Nothing about a real population holds `N` constant to the individual, and the fixation probability depends on `N` only through `r^-N`, so the result is not sensitive to the exact value; but a population that crashes and recovers behaves like the smaller size for drift purposes, which is what effective population size is trying to capture. Feeding a census count into these formulas will overstate the efficiency of selection, often by an order of magnitude. Generations overlap completely, one birth and one death at a time. The Wright-Fisher model makes the opposite choice, replacing the whole population at once, and gives $\rho_1 \approx (1 - e^{-2s})/(1 - e^{-2Ns})$ with `s` = `r` - 1, differing from the Moran answer by a factor of two in the exponent. Neither is more correct in general. The factor of two is a modelling choice about life history that is easy to lose track of when comparing published numbers. Fitness is constant. The whole apparatus assumes `r` does not depend on how common the type is, which rules out every frequency-dependent case, and those are the interesting ones in behaviour and in competition. That extension is the subject of the next post. There is no mutation. Absorbing states are absorbing only because nothing regenerates the lost type. With recurrent mutation the process has a stationary distribution instead, and "fixation probability" becomes a statement about the rate of transitions between mostly-`A` and mostly-`B` states rather than a final outcome. The population is well mixed, every individual competing with every other. The lattice posts showed how much structure can change a conclusion, and structure changes fixation probabilities too. ## Related tutorials - [Evolutionary games in finite populations](../games-in-finite-populations/) - [The replicator equation](../the-replicator-equation/) - [Checking a finite population model](../checking-a-finite-population-model/) - [Adaptive dynamics and evolutionary branching](../adaptive-dynamics-and-branching/) ## References Moran 1958 Math Proc Camb Phil Soc 54(1):60-71 (10.1017/S0305004100033193) Kimura 1962 Genetics 47(6):713-719 (10.1093/genetics/47.6.713) Nowak 2006 "Evolutionary Dynamics: Exploring the Equations of Life" ISBN 978-0674023383 Ewens 2004 "Mathematical Population Genetics 1: Theoretical Introduction" ISBN 978-0387201917