Evolution's Update Order Problem: Why Where Selection Acts Mathematically Matters

Where Does Evolution Happen? The Hidden Mathematics of How Selection Acts

In 1958, a statistician named Patrick Moran introduced a deceptively simple model for how genes spread through a population. Pick one individual to die, let another reproduce, replace the dead with the offspring, repeat. The model became foundational to population genetics because it captures the interplay of random drift, mutation, and selection in a finite population—problems that had stumped earlier theorists working with infinite populations and smooth approximations.

Seven decades later, Braha and de Aguiar have uncovered something strange lurking inside this classical framework. They show that where you place selection within the update rule matters—not just a little, but in ways that change the entire mathematical structure of the model. In a population with three or more allele types, the same selection pressures acting through the same fitness values can produce reversible Markov chains in one scheme and irreversible, non-equilibrium dynamics in another. The math that describes each case looks fundamentally different. And these aren't peripheral edge cases: they're the default when you're studying real genetic loci with multiple variants.

The finding has practical consequences for anyone trying to model balancing selection, where multiple alleles persist in a population over evolutionary time. It raises questions about which mathematical framework to trust when analyzing self-incompatibility systems, immune response genes, or color-pattern loci in butterflies. And it suggests that the familiar toolbox of population genetics—stationary distributions, detailed balance, product forms—may need to be unpacked carefully before being applied.

The Science

The Moran process has been a workhorse of population genetics since its introduction. The basic setup involves a haploid population of fixed size: n individuals, each carrying one of m possible alleles at a single genetic locus. In each discrete time step, one individual is chosen for death and replaced by a new individual (an offspring). The population stays constant; only the composition changes.

This overlapping-generations structure makes the model mathematically tractable. The state of the system—the count vector (k₁, k₂, ..., kₘ) where kᵢ is the number of individuals with allele i—evolves as a finite Markov chain. With recurrent mutation and selection, the chain has a stationary distribution: a predictable long-run average composition around which the population fluctuates generation after generation.

The question Braha and de Aguiar ask is subtly different from the usual approach. Most studies of selection in Moran models focus on the strength of selection or the mutation rate. This paper fixes those parameters and asks: what happens if you change when selection acts within each update step?

To answer this, they define three distinct update kernels—three ways to structure the death-replacement event—that differ only in the stage at which fitness becomes relevant.

Scheme I (selection during reproduction) begins when an individual is chosen uniformly at random for death. Suppose this focal individual carries allele j. A second parent is then chosen uniformly from the remaining n−1 individuals. If the second parent carries allele r, the offspring copies allele j with probability Sⱼ/(Sⱼ+Sᵣ) and allele r with probability Sᵣ/(Sⱼ+Sᵣ). The copied allele may then mutate before replacing the dead individual. Here, fitness enters the copying decision itself: the offspring is more likely to inherit alleles from fitter parents.

Scheme II (fitness-biased mate choice) also begins with uniform death selection. But this time, when choosing the second parent, fitness enters differently: the probability of picking a type-r individual is proportional to their fitness Sᵣ, not uniform. However, once both parents are chosen, copying is neutral—the offspring picks either parental allele with equal probability ½, and mutation acts afterward. Selection here acts through partner selection, not through the inheritance probabilities.

Scheme III (selection at death) reverses the logic: the individual to be replaced is chosen with probability proportional to inverse fitness, so fitter individuals are less likely to die. Once the focal individual is removed, reproduction is fully neutral: the second parent is chosen uniformly, copying happens with probability ½, and mutation acts on the result.

All three mechanisms favor fitter alleles on average. All three reduce to the neutral model when all fitnesses are equal. But as Braha and de Aguiar demonstrate, they define genuinely different Markov chains—and for m ≥ 3, only one of them retains the mathematical structure that makes exact analysis tractable.

The analysis covers a well-mixed haploid population of constant size n with m labeled alleles, parent-independent mutation (meaning mutation from allele r to allele i occurs at rate uᵢ regardless of r), and allele-specific fitness coefficients S₁ through Sₘ. This is the standard framework for studying polymorphism maintenance under recurrent mutation.

What They Found

The first major result concerns what happens with just two alleles. When m = 2, each scheme reduces to a one-dimensional birth-death chain, which always admits an exact stationary law. The three stationary laws exist and are exact—but they are not the same. Even when both schemes favor allele 1 equally strongly, the probability distribution over how many copies of allele 1 exist in the population takes different mathematical forms in each scheme.

The neutral case (S₁ = S₂) serves as a benchmark. With parent-independent mutation and no selection, the stationary distribution is beta-binomial: the number of allele 1 copies follows a well-characterized distribution parameterized by the mutation rates and population size. This is classical result, derived originally by Wright and later refined by Kimura and others.

When selection is weak, Braha and de Aguiar derive expansions showing how each scheme perturbs this neutral benchmark. In Scheme I, for example, the stationary probability of finding k copies of allele 1 takes the form:

$${\pi }_k^{(I)}\approx \frac{\left(\genfrac{}{}{0ex}{}{n}{k}\right)({\alpha }_1{)}_k({\alpha }_2{)}_{n-k}}{({\alpha }_1+{\alpha }_2{)}_n}\left[1+\epsilon \left(\frac{k}{2}-\frac{n{\alpha }_1}{{\alpha }_1+{\alpha }_2}\right)\right]$$

where $(a)_k$ is the rising factorial (a)(a+1)...(a+k−1), α₁ and α₂ are scaled mutation parameters, and ε measures the relative strength of selection. The neutral beta-binomial appears as the leading term; the selection correction appears as an O(ε) perturbation that depends on k, the current count of allele 1.

Scheme II gives a different correction term—one that depends not just on k but on the full population structure through the fitness-weighted total $\xi(k) = \sum_ℓ S_ℓ k_ℓ$. The correction is no longer expressible as a simple polynomial in k.

The key insight emerges when the authors examine what happens as the number of allele types grows. For m ≥ 3, the three schemes diverge not just quantitatively but qualitatively.

Schemes I and II become non-reversible. A Markov chain is reversible if it satisfies detailed balance: the flow of probability from state A to state B equals the flow from B to A, weighted by their stationary probabilities. Reversibility is a powerful mathematical property because it implies the stationary distribution takes a product form—the probabilities for each allele type separate nicely, making them tractable to compute and interpret.

But Braha and de Aguiar prove that when fitnesses are unequal across three or more alleles, Schemes I and II generally fail to be reversible. There is no closed-form stationary distribution with a detailed-balance product structure. The Markov chains still have stationary distributions—the population doesn't explode or collapse—but those distributions don't factor into independent pieces, and they can't be written down in the familiar exponential-of-energy form that makes reversible chains so convenient.

Scheme III remains reversible for every m. Selection at death, they show, preserves reversibility regardless of how many allele types exist or how their fitnesses differ. The stationary distribution takes a clean form: a Dirichlet-multinomial core (the same structure that appears in the neutral case) multiplied by an explicit fitness-dependent factor:

$${\pi }^{(III)}(k)=\frac{n!}{k_1!\cdots k_m!}\frac{({\alpha }_1{)}_{k_1}\cdots ({\alpha }_m{)}_{k_m}}{({\alpha }_0{)}_n}\prod _{i=1}^m{S}_i^{k_i}$$

where α₀ = α₁ + ... + αₘ. The neutral Dirichlet-multinomial stationary structure survives intact; selection merely weights each configuration by the product of fitness factors raised to the power of its allele count. This is the Dirichlet-multinomial core with an explicit fitness factor multiplying it.

The contrast is stark: the same fitness pressures, the same mutation rates, the same population—yet only one of the three update schemes admits a closed-form stationary law that's both exact and reversible.

The authors push further, showing that for two alleles, all three mechanisms can operate simultaneously without losing exact solvability. When the death-target is biased against fitness, reproduction involves fitness-biased mating, and inheritance probabilities are fitness-dependent, the stationary distribution can still be written down in closed form. The three effects combine through a set of coupled equations that admit an explicit solution.

Selection Shifts Stationary Distribution Differently Across Schemes

Label	Value
Neutral (S₁=S₂=1)	0.02 probability
Scheme I - Reproduction	0.045 probability
Scheme II - Mate Choice	0.048 probability
Scheme III - Death	0.038 probability
All Three Combined	0.06 probability

Stationary distributions for allele 1 frequency in n=100 population with α₁=α₂=2 and S₁=1.10, S₂=1

This figure illustrates the practical consequences. In a population of n = 100 with symmetric mutation rates (α₁ = α₂ = 2, above the neutral critical value), the black line shows the neutral stationary distribution when S₁ = S₂ = 1: a bell-shaped curve centered around roughly k = 50 copies of allele 1. When selection favors allele 1 (S₁ = 1.10, S₂ = 1), the three schemes shift this distribution in different ways. Scheme I (orange) and Scheme II (red) produce similar but not identical shifts toward higher k values. Scheme III (blue) produces a distinctly different shift. The dashed line shows the combined effect of all three mechanisms acting together—remarkably, this combined distribution is also exact and computable.

Stronger Selection Produces Greater Frequency Shifts

Label	Value
ε = 0.02	0.15 shift
ε = 0.05	0.35 shift
ε = 0.10	0.6 shift

Shift in mean allele frequency under Scheme II selection for different ε values (n=100, α₁=α₂=2)

This figure zooms into Scheme II across a range of selection strengths. With α₁ = α₂ = 2 (panel a), above the neutral critical mutation rate, the neutral distribution (black) sits at intermediate frequencies. As selection strengthens from ε = 0.02 to ε = 0.10, the distribution shifts rightward, concentrating more probability mass on high-k states where allele 1 dominates. But at the critical mutation value α₁ = α₂ = 1 (panel b), the neutral distribution sits exactly at k = n/2 by symmetry. Here, selection can only break the symmetry—there is no "neutral benchmark" to shift toward. At lower mutation rates (panel c, α₁ = α₂ = 0.5), both the neutral and selected distributions concentrate near the boundaries k ≈ 0 and k ≈ n, reflecting the strong tendency toward fixation when mutation is weak.

Why This Changes Things

Population geneticists have long known that update rules matter. The Wright-Fisher model and the Moran model, though similar in spirit, can give different answers for fixation probabilities and timescales. Evolutionary graph theory has shown that changing the order of events in structured populations can alter whether a new mutant fixates or dies out. Work by Lieberman, Nowak, and others demonstrated that the same selective pressure can have different evolutionary consequences depending on whether it acts on birth, death, or replacement.

What Braha and de Aguiar add is a precise characterization of when these differences are merely quantitative and when they become qualitative. For two alleles, the differences are quantitative: all three schemes admit exact stationary laws, just with different numerical values. A researcher analyzing any of the three would arrive at correct predictions for allele frequency distributions, albeit with different formulas.

For three or more alleles with unequal fitnesses, the differences become qualitative. Schemes I and II are non-reversible; Scheme III is reversible. This isn't a technical subtlety—reversibility determines whether you can apply the powerful machinery of statistical mechanics to the problem, whether detailed balance holds, whether the stationary distribution has a product form. These are not optional features; they shape what calculations are possible and what interpretations are valid.

The implication for modelers is clear: the choice of update scheme is not just a matter of computational convenience or biological interpretation. It determines the mathematical structure of the entire analysis. If you're studying a locus with multiple allele types under selection—say, the MHC complex that governs immune recognition in vertebrates—using the wrong update scheme might lead you to a mathematically convenient but biologically inappropriate stationary distribution.

Consider the classic problem of balancing selection, where multiple alleles are maintained in a population over long evolutionary timescales. Self-incompatibility loci in plants, where dozens or hundreds of alleles may coexist, are a canonical example. The mathematical analysis of these systems has often relied on reversibility assumptions or diffusion approximations that may be justified for one update scheme but not another.

The paper also matters for understanding the relationship between forward-time and backward-time (genealogical) analyses. Much of the recent work on Moran models with selection uses ancestral processes—tracing lineages backward from present-day individuals to find their common ancestors. These methods have been extremely productive for two-type models, where the ancestral selection graph and related constructions give exact results. But for multiallelic models with allele-specific fitnesses, the non-reversibility of Schemes I and II may complicate the genealogical interpretation. If the forward-time chain is non-reversible, its time-reversal has a different structure—and the standard connections between forward-time equilibria and backward-time branching processes may not apply cleanly.

There's also a conceptual point worth dwelling on. The three schemes represent genuinely different biological mechanisms, not just mathematical variations. Scheme I, where selection acts during reproduction through inheritance probabilities, captures situations like gamete competition or fertilizing ability that depends on parental genotypes. Scheme II, where selection acts through mate choice, captures sexual selection or pollinator preference, where individuals with certain traits are more likely to be chosen as partners. Scheme III, where selection acts at death, captures differential survival—individuals with certain genotypes are more likely to survive long enough to reproduce.

These are evolutionarily distinct processes. A population where fitter genotypes are better at surviving (Scheme III) might have different temporal dynamics than one where fitter genotypes produce more competitive gametes (Scheme I), even if the net fitness differences are identical. The stationary distribution encodes these differences. A modeler who treats sexual selection and viability selection as mathematically equivalent—plugging both into the same update scheme—may be making an implicit assumption that the biology doesn't support.

What's Next

Several open questions emerge from this work.

First, the classification of reversibility suggests that Schemes I and II might be recoverable through alternative mathematical approaches. Non-reversible chains still have stationary distributions— they're just harder to compute and interpret. Perhaps there are other ways to express these distributions that don't rely on detailed balance. The authors mention that for two alleles, local detailed balance holds even when global reversibility fails; whether this extends to multiallelic cases is unclear.

Second, the biological interpretation of the three schemes deserves further exploration. The paper establishes the mathematical structures; the evolutionary interpretations remain to be fleshed out. When is selection during reproduction the appropriate model? When does fitness-biased mate choice better capture the biological mechanism? The comparative analysis in Section 10, which shows how the three distributions diverge for specific parameter values, is a start, but systematic exploration of the biological parameter space is needed.

Third, the results raise questions about continuous-time versus discrete-time formulations. The paper works in discrete time, with one death-replacement event per step. Many applications of the Moran model use continuous-time formulations, where each individual has a constant hazard of death and reproduction. The mapping between discrete-time and continuous-time versions is not always one-to-one, and the results here may or may not carry over.

Fourth, the authors' use of parent-independent mutation is standard but not universal. Parent-dependent mutation, where the mutation probability depends on the source allele, is biologically relevant in some contexts (notably, transitions between purines and between pyrimidines occur at different rates). Extending the analysis to parent-dependent mutation would broaden the applicability.

Finally, the weak-selection expansions reveal that the three schemes are asymptotically equivalent when selection is weak— they differ only in their O(ε) corrections. This suggests that for small fitness differences, any of the three schemes might serve as a reasonable approximation to the others. The practical question is whether real biological systems typically involve weak selection or strong selection at the loci where multiple alleles persist. If selection is typically weak, the choice of scheme may be less consequential for practical applications. If selection is strong, the scheme matters enormously.

The paper concludes by situating this work in the broader context of evolutionary dynamics in finite populations. The key contribution is not just the specific results but the methodological demonstration that the placement of selection within an update rule can be decisive—that the same selection pressures, the same mutation rates, and the same population size can yield qualitatively different mathematical structures depending on the order of events.

This is a reminder that mathematical models are not neutral containers. The choice of update rule is a modeling decision, and bad modeling decisions lead to bad conclusions even when the mathematics is internally consistent. For population genetics, where the Wright-Fisher and Moran frameworks have been used for decades to reason about everything from the maintenance of polymorphism to the evolution of sex, the lesson is clear: before applying a stationary distribution, make sure it corresponds to the right update scheme. The mathematics will be precise either way. Whether it describes your population is another question.

The paper opens up new territory at the intersection of mathematical population genetics and stochastic processes. It shows that even in one of the most thoroughly studied models in all of evolutionary biology, there are structural features that have been overlooked—features that matter not just for the specialist but for anyone trying to understand how genetic variation persists in natural populations.