Pop Gen Calculator

genotype   aa      Aa        AA
fitness    1      1 + h s   1 + s

deno run -A html/tools/popgen/sim/tests.ts
deno run -A html/tools/popgen/sim/worked-chain.ts
deno run -A html/tools/popgen/make-tables.ts

What each calculator does

Twenty-six calculators in nine groups.  Every one names its primary source with a page or an equation number, because "standard result" is not a citation.

Fixation

Probability a new mutation fixes.  Kimura's diffusion result: a single new copy, in a population of effective size Ne, eventually taken up by everybody.

u(p0) = (1 - e^(-2 Ne s p0)) / (1 - e^(-2 Ne s)),   p0 = 1/(2N)

Kimura 1962, equation (8), page 715, written on the homozygote s.  His own equation (8) reads (1 - e^(-4Nsp))/(1 - e^(-4Ns)) because his equation (7) defines s per allele copy.  Setting p0 = 1/(2N) gives his equation (10).  Arbitrary dominance goes through his equation (13), page 716, by quadrature.

Haldane's 2s, with the offspring-variance correction.  A beneficial mutation survives with probability about 2s no matter how large the population is -- but only if reproduction is Poisson.

P(fixation) ~ 2 s_carrier / sigma^2

Haldane 1927, Proc. Camb. Phil. Soc. 23:838-844.  Kimura 1962 page 715 rederives it as the large-4Ns limit of his equation (11).  The s here is the carrier advantage, which is h times the homozygote coefficient, and that is the one place on this page where s is not on the homozygote.

How long fixation takes.  Two answers that could hardly be further apart:

neutral, given it fixes:  t1 = 4 Ne neutral, given it is lost: t0 ~ 2 (Ne/N) ln(2N)

Kimura & Ohta 1969, equations (14) and (15), page 765, and equation (16), page 766.  In a population of ten thousand, a neutral mutation that fixes took forty thousand generations; the 99.99% that are lost were gone in about twenty.  Selected fixation times come from the numerical evaluation of their equation (12), page 765.

Drift

The drift barrier.  How small a fitness effect has to be before drift overwhelms selection.  There is no single number, and this is where a great deal of published confusion lives:

master parameter:  4 Ne s_het = 2 Ne s_hom the strict cut:    |s_hom| = 1/(2 Ne),  which is |s_het| = 1/(4 Ne) Kimura's own band: twice as wide again

Ohta 1973, Nature 246:96-98, for the argument.  Kimura 1979, PNAS 76:3441, states his convention outright -- "the mutant is semidominant in fitness so that selection coefficients against the mutant homo- and heterozygotes at a site are 2s' and s', respectively" -- and then integrates his effectively-neutral class out to 1/(2 Ne) in that s', which is twice the strict cut.  Both are his.  The way out is to print the master parameter, which is unambiguous once you say which s built it.  The barrier is not a wall either: a mutation exactly at 2 Ne s = 1 still fixes at 58% of the neutral rate.

Where mutation, selection and drift balance.

fraction of sites sitting wrong = 3 / (3 + e^(2 Ne s))

The Wright 1931 stationary distribution with three wrong letters for every right one.  A site selection cannot see is wrong three times in four.

Harmful fixations per generation, over a measured distribution of fitness effects.

rate = N U E[u(s)],   s ~ Gamma(shape, mean)

Kimura's integral averaged over a gamma DFE with shape near 0.2 (Eyre-Walker & Keightley 2007; Boyko et al. 2008; Kim et al. 2017).  Change the shape parameter and watch it move: that is the largest real uncertainty in the whole genetic-entropy argument.

Bottlenecks

Heterozygosity and inbreeding after a bottleneck.

one pair:  F(t) = 1/4 (1 + 2 F(t-1) + F(t-2));   general Ne:  H(t) = H0 (1 - 1/(2 Ne))^t

Two formulas, and the card picks by what you type.  Enter Ne = 2 and it is read as ONE BREEDING PAIR, which is not a Wright-Fisher population of size two: the pair's children each take one allele from each parent, and those are two independent draws from the source, so nothing is lost at generation one.  That case uses the exact dioecious recursion, Crow & Kimura 1970 Table 3.1.2, page 64 -- the series 1, 3/4, 5/8, 1/2, 13/32, 21/64 that the population-genetics article publishes, and the Tables tab prints in full.  Any other Ne uses the geometric form, Crow & Kimura chapter 3 problem 15, page 114.  The two are far apart:  at generation five the geometric shortcut says 23.7% where the exact answer is 40.6%, and at generation ten it is out by a factor of 2.5.  Use the exact power rather than e^(-t/2Ne) in the general case too;  at Ne = 2 the exponential shortcut is 46% wrong after ten generations, and it errs in the direction that flatters the bottleneck.

Effective size across generations of unequal size.

Ne = T / (1/N1 + 1/N2 + ... + 1/NT)

Crow & Kimura 1970, equation (7.6.3.34), page 360, after Wright 1938a.  Their own worked example, quoted exactly: a population that "increases in five generations in geometric series from 10 to 10^6 and then returns to 10 and repeats the cycle" has a representative effective size of 54, "much nearer to the minimum number than to the maximum".  This is the single most useful fact about Ne in any bottleneck argument, in either direction: one generation at Ne = 2 inside ten thousand generations at Ne = 10,000 gives 6,667.

Alleles surviving a founder event.

E[k] = k (1 - (1 - 1/k)^(ploidy x founders))

Nei, Maruyama & Chakraborty 1975, Evolution 29:1-10.  Their point, which the calculator makes visible: allelic richness collapses in a bottleneck far faster than heterozygosity does, because heterozygosity is carried by common alleles and richness by rare ones.  Four gene copies drawn from a locus with four equally common alleles keep 2.73 of them, while heterozygosity falls only from 0.750 to 0.5625.

Lethal equivalents.  Inbreeding depression, and the factor of two that catches nearly everybody:

excess mortality = 1 - e^(-B F),   B per GAMETE

Morton, Crow & Muller 1956, PNAS 42:859, define B per gamete in their own words: "the total genetic damage here measured is 1.5-2.5 lethal equivalents per gamete, or 3-5 per zygote."  Conservation genetics quotes the diploid count -- 3.14, 12, 16 -- so a published figure has to be halved before it goes into the exponent.  The halving is what reproduces the published answers.  Frankham et al. 2014, Biol. Conserv. 170:58, print a 26% fitness decline for 12 diploid lethal equivalents at F = 0.05, and that comes out at 25.9% halved against 45.1% unhalved.

Reading Ne off data

Effective size from nucleotide diversity.

theta = 2 x ploidy x Ne x mu    ->    Ne = theta / (2 ploidy mu)

Tajima 1983 for E[pi] = theta; Hartl & Clark 1997 page 302.  pi is per site.  Ploidy 2 gives the familiar 4 Ne mu; ploidy 1 is mtDNA and the Y, where using 4 halves the answer, and that is the commonest single error in mitochondrial-Eve arithmetic.

Watterson's theta from segregating sites.

theta = S / (a(n) L),   a(n) = 1 + 1/2 + ... + 1/(n-1)

Watterson 1975; the form used is Hartl & Clark 1997 equation (7.34), page 302.  n is the number of sequences, two per diploid individual.  The harmonic correction is the whole estimator: drop it and the answer is out by a factor of nineteen at n = 100.

Recombination and linkage

Map distance to recombination fraction, Haldane and Kosambi side by side.

Haldane  r = (1 - e^(-2d)) / 2          Kosambi  r = tanh(2d) / 2

Haldane 1919, J. Genetics 8:302, and Kosambi 1944, Ann. Eugenics 12:172-175.  A Morgan is a count: one expected crossover on that interval per gamete per meiosis.  So a genome's map length in Morgans is the average number of crossovers in a gamete, and "34 Morgans" and "34 crossovers per meiosis" are the same sentence.  The defaults are the whole human autosome as deCODE measured it -- Halldorsson et al. 2019, Science 363:eaau1043, built from 4,531,535 directly counted crossovers: male 2,602.3 cM, female 4,180.4 cM, sex-averaged 3,391.4 cM over 2,753.9 Mb, which is 1.23 cM per Mb.  At that distance both functions return exactly 1/2, which is the honest answer: two loci a whole genome apart are simply unlinked.  The two functions disagree by tens of percent over long intervals because Kosambi allows crossover interference and Haldane forbids it, so a Kosambi centimorgan must never be summed into a Haldane map.  The Verify tab treats that mixture as a negative control.

LD decay, and the effective size it implies.

E[r^2] = 1 / (1 + 4 Ne c),    Ne = (1/r^2 - 1) / (4 c)

Sved & Feldman 1973, restated and rederived in Sved & Hill 2018, Genetics 209:632-633.  Linkage disequilibrium decays with distance, so an r^2 measured between two sites reports the effective size of roughly 1/(2c) generations ago (Hayes et al. 2003, quoted by Sved & Hill on page 633).  Tightly linked pairs speak about ancient Ne and loose pairs about recent Ne, which is why "Ne from LD" is a history rather than a number, and why an Ne quoted without its c is meaningless.  A finite sample of n chromosomes inflates r^2 by about 1/n even at complete equilibrium, which is why Tenesa et al. 2007 fit E(r^2) = 1/(1+4Nec) + 1/n.  Ohta & Kimura 1971, Genetics 68:574, equation (9), is printed beside the answer as an independent route through a different paper; it is sigma_d^2, the ratio of expectations, not E[r^2], and the two converge at large 4Nec.

Background selection: how much linked selection costs.

B = e^(-U / M)

Charlesworth 2012, Genetics 190:10, equation (5b), in his own words: "the effect of BGS depends only on the density of deleterious mutations per unit map length, as first pointed out by Hudson and Kaplan (1994) and Barton (1995)."  It does not depend on how strong the selection is.  Across the whole human autosome -- U about 1.5 per diploid genome, M = 33.914 Morgans -- the answer is a reduction of only 4%.  Wind M down and the cost appears: the same U over one Morgan leaves 22%.  Measured genome-average reductions are 19-26%, larger than 4%, because deleterious mutations are concentrated near genes where the local U/M is far higher than the genome average.  So Ne is not one number for a genome.  M is in Morgans; typing centimorgans divides the exponent by a hundred, and that too is a negative control.

Structure and Hardy-Weinberg

Hardy-Weinberg proportions, with inbreeding.

AA = p^2(1-F) + pF,   Aa = 2pq(1-F),   aa = q^2(1-F) + qF

Wright 1951, Ann. Eugenics 15:323, Table 1, in the panmictic-and-fixed-components column.  The F = 0 case is Hardy 1908, Science 28(706):49, and Weinberg 1908, page 378 -- and a point of credit worth getting right is that Hardy never writes p^2 : 2pq : q^2;  he works in counts, as p : 2q : r, and the modern algebraic form is Weinberg's.  The proportions arrive in one generation of random mating, which was Hardy's actual point.  Inbreeding does not touch the allele frequency at all;  it moves copies out of heterozygotes into homozygotes, pqF into each homozygote class.  The same algebra with F_ST in place of F gives the average genotype frequencies across subpopulations (Hartl & Clark 1997, equation 4.7, page 128), and those are two different quantities sharing one formula.

Testing Hardy-Weinberg: chi-square and the exact test.

chi2 = SUM (obs - exp)^2 / exp,  df = 1

Hartl & Clark 1997, equation (3.3) and the df rule, pages 81-82; the exact test is Wigginton, Cutler & Abecasis 2005, Am. J. Hum. Genet. 76:887, equation (1), computed by their recurrence, equation (2), page 888.  One degree of freedom, not two: three classes, minus one for the total, minus one more for the allele frequency estimated from the same data.  The defaults are the MN blood group in 1,279 Britons (Race & Sanger 1962, via Crow & Kimura 1970 Table 2.2.2, page 36), which gives chi-square 0.0272 and P = 0.869.  The identity chi2 = n F-hat^2 is exact for two alleles, which is why "chi-square test of Hardy-Weinberg" and "test that F = 0" are the same test.  Use the exact test when an allele is rare: Wigginton et al. page 892 measured the chi-square test's true error rate at 0.0208 for a nominal 0.01.  And a significant result names no cause -- the classic failure, 18,121 Oslo schoolchildren scored for colour blindness, violates the two-allele assumption, not random mating, because "colour blindness" is four conditions at two loci.

F_ST from subpopulation frequencies, and the migration inversion.

F_ST = Var(p) / [pBar (1 - pBar)] = 1 - H_S/H_T Nm = (1/F_ST - 1)/4 diploid,   (1/F_ST - 1)/2 haploid

Wahlund 1928, Hereditas 11:65-106, as Hartl & Clark 1997 equation (4.6), page 126 and equation (4.3), page 118.  Pooling groups that differ in allele frequency makes an excess of homozygotes even though every group is internally in Hardy-Weinberg proportions, and that excess is F_ST.  The inversion to migrants per generation is Wright 1951, page 331, for diploids and Hudson, Slatkin & Maddison 1992, equation (2), page 583, for haploids, who state the whole relation in one sentence: "For diploids the 2 in Equation 1 is replaced by a 4."  A factor of two, so an Nm printed without its ploidy is unreadable.

The inversion is also fragile, and the warning is part of the answer.  Whitlock & McCauley 1999, Heredity 82:122, took five loci, ten demes and fifty diploids per deme from a metapopulation whose true Nm was 50, and found 95% of the estimates between 11 and 283.  Their verdict, on page 117: "it is rare that F_ST can be translated into an accurate estimate of Nm."  The calculator also prints Wright's exact island-model recursion beside his own 1/(4Nm+1) approximation, which needs m small rather than Nm small and runs 17% high at Nm = 10.

Structure read as a bottleneck.

apparent decline = n + (n-1)^2 / (n M)

Mazet, Rodriguez, Grusea, Boitard & Chikhi 2016, Heredity 116:364, equation (19).  A population that never changed size, but is divided into demes, is read by any method assuming panmixia as a decline by this factor.  The defaults are Mazet's own human case, page 367: ten demes of 530 diploids, 5,300 diploids in total, constant forever, at M = 0.55 migrant genes per deme per generation.  That reports an apparent ancestral size of 13,105, which is the published human PSMC peak of about 13,000.  The apparent current size is the deme size, 530, not the metapopulation.  Sampling two lineages from the same deme fakes a decline;  from different demes it fakes an expansion.  Their abstract states the general result: "any demographic model with structure will necessarily be interpreted as a series of changes in population size by inference methods ignoring structure."  It cuts both ways, and this is the calculator the Babel and ark arguments need.

The site frequency spectrum

The neutral spectrum, raw and scaled.

E[xi_i] = theta / i,   and therefore E[i xi_i] = theta for every i

Achaz 2009, Genetics 183:250, equations (2)-(4), attributing the expectation to Fu 1995; the folded form is his equation (10), page 251.  The spectrum is a histogram of polymorphisms by allele count, and its whole neutral content is one fact about tree shape: the coalescent devotes branch length 4Ne/i to branches with i descendants.

Plot i x xi_i, not xi_i.  Under the standard neutral model the scaled version is flat, and on the raw plot every spectrum in existence looks like the same decaying hyperbola, so nothing is visible.  Achaz endorses exactly this on page 250, after Nawa & Tajima.  The Verify tab treats reading the raw spectrum as flat as a negative control.  n is chromosomes: 108 people give 216.  When n is even, folding puts derived count n/2 onto itself and it must not be doubled -- that is the Kronecker delta in Achaz's equation (10), and leaving it out inflates exactly the bin a bottleneck would be read off.

Tajima's D.

D = (pi - S/a_n) / sqrt(Var)

Tajima 1989, Genetics 123:589, equation (38).  Two estimators of the same theta, subtracted:  pi weights common variants, Watterson's S/a_n counts them all equally.  The defaults are the measured chr22 spectrum of the 108 Yoruba in 1000 Genomes phase 3 -- 278,607 segregating sites in 216 chromosomes -- computed from the stored array rather than typed in, and they give D = -0.512 with 22.0% singletons against a neutral 16.8%.  Negative means an excess of rare variants: growth, a sweep, or purifying selection, and the statistic does not say which.  A negative genome-wide D is not by itself evidence of growth.

Ancestral misorientation, which is a one-sided error.

xi'_i = (1 - pM) xii + pM xi(n-i)

Baudry & Depaulis 2003, Genetics 165:1621, who simulate exactly this "by exchanging the frequency of the derived and ancestral variant ... with probability given by the misorientation parameter".  The unfolded spectrum needs an outgroup to say which allele is ancestral, and the outgroup can be wrong.  Because xi_1 is theta and xi_215 is theta/215 at n = 216, the sites flow almost entirely one way: a 5% misorientation rate barely dents the singletons and multiplies the top bin 11.7-fold.  Their measured rates were 3.3% mean in Drosophila simulans and 5.1% in Arabidopsis thaliana.  What it does to the tests, in their words: Fu & Li's D and F are barely affected, but Fay & Wu's H "if PM = 0.15 with a sample size of 50 ... will be rejected ~25% of the time" -- five times the nominal 5%, from homoplasy alone, with no selection anywhere.  The folded spectrum is exactly immune.

Autosome, X, Y and mtDNA: the 4 : 3 : 1 : 1 table.

theta = 4 Ne mu autosome,  3 Ne mu X,  1 Ne mu Y and mtDNA -- all on ONE Ne

Wilson Sayres et al. 2014, PLoS Genet. 10(1):e1004064, page 2, state the ratio in words; the sex-ratio effective sizes are Wright 1931, Genetics 16:109.  Count gene copies at an even sex ratio: 2N autosomal, 1.5N X, 0.5N Y, 0.5N mitochondrial.  The Y and the mitochondrion come out equal, because there are exactly as many mitochondrial genomes as Y chromosomes.  Move the sex ratio and polygyny raises mtDNA and lowers the Y.  Using the autosomal theta = 4 Ne mu on mitochondrial data is four times too high, not eight;  eight is what you get by applying the ploidy correction (4 to 2) and the Ne correction (Ne to Ne/4) as if they were two different corrections, when they are one correction stated twice.  In the real world the Y sits about ten times below its 0.25 expectation, from purifying selection on a chromosome that never recombines, while mtDNA sits at expectation -- an asymmetry any mitochondrial-Eve or Y-Adam argument has to carry.

Load

Haldane mutation load.

mean fitness = e^(-U),   load = 1 - e^(-U)

Haldane 1937, Am. Nat. 71:341, in his own words: "the loss of fitness to the species depends entirely on the mutation rate and not at all on the effect of the gene upon the fitness of the individual carrying it."  U is per diploid genome per generation.  This is the equilibrium load -- the standing cost of mutations selection is still removing.  It says nothing about mutations that have already fixed, and confusing the two has wasted decades of this argument.

Muller's ratchet: the least-loaded class.

n0 = N e^(-U/s)

Haigh 1978, Theor. Popul. Biol. 14:251-267.  This is the one haploid formula on the page, and it must never be fed the diploid U every other box takes.  Asexual populations and non-recombining chromosomes only: recombination rebuilds the least-loaded class, which is why the classical ratchet is not by itself an argument about the human autosomes.

The whole chain: from a measured diversity to a date

The last calculator is the one worth the most, and the only one that is not a single formula.  It carries one measured number through four steps and reports how much of the final answer is data and how much is convention.

measured pi  ->  Ne = pi/(4 mu)  ->  TMRCA = 4 Ne (1 - 1/n)  ->  years

Every input is read out of the source itself:

Under the conventions this project uses, that chain lands at 2.80 million years.

Three things about that number matter more than the number.

Only the first row is a count of anything.  Yu et al. sequenced DNA and Porubsky et al. counted mutations in a pedigree.  Everything below those two lines is arithmetic, and arithmetic inherits whatever went into it.

The error bar is dominated by a term most published dates leave out.  The mutation rate's confidence interval and the diversity's sampling error together give about +/- 16%.  The standard deviation of the time to the most recent common ancestor, from the coalescent alone, is 55% -- and it does not shrink with more sequencing, because there is one human genealogy and its TMRCA is a single draw from a very wide distribution.  A date quoted with only the measurement interval understates its own uncertainty by roughly a factor of four.

The conventions move the answer further than the measurements do.  Switch one convention at a time, holding every measured input fixed:

Together those span a factor of 4.3, and every one of the answers in that range is arithmetically correct.  What separates them is entirely which convention each was computed under.  That is the argument for printing the convention beside every number, and it is why this whole toolkit has one.

The phylogenetic branch carries an extra floor nobody can sequence their way out of.  dos Reis & Yang 2013, J. Syst. Evol. 51:31, fitted 8,708,584 sites of six primate genomes and found that "for every one million years of divergence, 0.322 million years are added to the CI width of time estimates."  That is +/- 16% of the age, permanently, on any date whose clock was calibrated on fossils.

What the Verify tab does

Ninety-one checks in ten groups, run in the browser and identically under Deno.  Four of the groups exist to make the other six worth something.

The dice.  A Monte Carlo check is only as good as its generator.  Uniform, binomial and gamma draws are graded against their own known means and variances, and the exact binomial is compared term by term with the exact probability mass function over 200,000 draws.

Paper against paper.  Identities and limits that must hold if the formulas were transcribed correctly: u(s,p) + u(-s,1-p) = 1 across three dominance values, the general-dominance quadrature collapsing onto the closed form at h = 1/2, the selected fixation-time integral collapsing to 4 Ne as s goes to zero, Crow & Kimura's own worked harmonic-mean example, and their Table 3.1.2 founding-pair series at generations five and ten, 13/32 and 9/64 exactly.

The conventions, asserted.  Every line of the notation block is a claim about how the code behaves, and a convention stated only in a comment is a convention that will be broken.  s_hom = 2 s_het is checked at eighteen parameter combinations to exact equality; the theta ploidy factors are checked as the ratio 1 : 4/3 : 2; the lethal-equivalent halving is checked against Frankham's two published percentages; the generation multiplier is checked to be carried in the return value rather than applied silently.

One question, every path.  Four independent routes to "a neutral mutation takes 4 Ne generations"; the harmonic number against the explicit sum; Ewens' expected allele count by direct summation against the digamma identity; the worked chain rebuilt step by step from the calculator.  This group exists because five defects shipped in this directory in a single week, and every one of them was a quantity computed twice with nothing anywhere comparing the two answers.

Paper against dice.  Every formula a one-locus Wright-Fisher run can reach is measured by running it, with the Monte Carlo standard error printed beside the measurement.  A validation without an error bar is not a validation.

Recombination and linkage.  Eight checks, and the map anchors are counts rather than model outputs: deCODE's sex-averaged map must be the mean of its paternal and maternal maps, the human genome must average 1.23 cM per Mb, and the mapping functions must invert themselves exactly.  Ohta & Kimura's own worked value -- R = 1/2 gives sigma_d^2 = 3/13, printed on the page they derive it -- pins every coefficient in their equation (9), and their formula and Sved's must converge in the large-4Nec limit, which is two papers agreeing rather than one formula agreeing with itself.

Structure and Hardy-Weinberg.  Thirteen checks against published worked examples: the MN blood group's chi-square and P, Wigginton's own exact-test table (including his point that the chi-square test rejects that configuration and the exact test does not), Wright's Linanthus parryae F-statistics from Hartl & Clark's Table 4.1 by two different routes, Wright's island model exact against his own approximation, and Mazet's ten demes of 530 diploids reproducing the human PSMC peak at 13,105.

The site frequency spectrum.  Nine checks, four of them against a real measured spectrum -- chr22 in the 108 Yoruba of 1000 Genomes phase 3, 278,607 segregating sites in 216 chromosomes, 22.0% singletons against a neutral 16.8%, Tajima's D of -0.512.  Those are the numbers the population-genetics article prints, so a change here is a change to published text.

Already published here.  Six numbers this site has printed elsewhere, recomputed from independent code: Drift Lab's four validation targets, and the genetic-entropy article's 1.11 harmful fixations per generation and its 60,780-against-56,255 check on Masel's parameters.

Negative controls -- required to fail.  A checker that cannot fail proves nothing.  Twelve deliberately wrong inputs are fed in and the comparison is required to break: 4 Ne s where 2 Ne s belongs, a cost passed as a benefit, the exponential heterozygosity shortcut at Ne = 2, the geometric shortcut applied to a dioecious founding pair, Watterson without his harmonic correction, theta = 4 Ne mu applied to mitochondrial DNA, the Monte Carlo graded against the wrong neutral expectation, the haploid Nm inversion applied to diploid autosomal data, the raw spectrum read as though it were the flat one, a Kosambi map distance summed into a Haldane map, B = e^(-U/M) with M in centimorgans, and the folded balanced class doubled.  Each prints how far off it is, so the page can say which mistakes matter and by how much.

Companions

• Drift Lab -- the WebGPU Wright-Fisher simulator. This page is its closed-form companion, and the two are checked against each other in the Verify tab.
• Genetic Entropy -- the article whose central table this calculator regenerates.
• A general population-genetics reference article is in preparation; its research dossiers and the headless research modules that produced their numbers live beside it in html/unpublished-ai/population-genetics/. Those modules import from sim/popgen.ts above rather than reimplementing it, and a cross-module consistency test in research/sim/consistency.ts asks all of them the same thirty-one questions and fails if any two disagree.