Migration Networks from Neutral Allele Frequencies

During the Delta wave, genomic surveillance in England was producing tens of thousands of SARS-CoV-2 genomes a week. Follow a common mutation in two regions, tracking the fraction of genomes in each that carry it: while travel between them is suppressed the two frequencies tend to wander apart, and once travel resumes they converge again. Nobody recorded the journeys that caused the convergence, yet their imprint was left for us to read.

Whether that imprint is there at all depends on which mutations we follow, and the ones tracked here are the ones that do not matter. A neutral mutation is neither favored nor purged, so on average the only thing that moves its frequency in a region is the arrival and departure of the individuals carrying it. Left alone it still wanders (that is genetic drift, and it is the noise the rest of this note has to survive), but its expected change is zero, and an estimator can be built on an expectation. Migration is the only thing left pushing in a direction, and that is what makes a neutral mutation readable.

What follows is the single linear equation that turns the convergence of two such frequencies into an estimate of who moved where, and its application to two datasets with almost nothing in common beyond neutrality: weekly SARS-CoV-2 genomes across regions, and ancient human genomes sampled centuries apart.

A Metapopulation of Neutral Alleles

We consider a population divided into \(n\) subpopulations (regions, districts, age groups, etc.) and follow the frequency \(X_i(t)\) of one neutral allele in each. Under neutrality the frequencies at a later time depend linearly on the frequencies at an earlier one:

\[ X_i(t+\Delta t) = \sum_{j=1}^n A_{ij}\, X_j(t) + \eta_i \]

where \(A\) is a right-stochastic matrix, its elements non-negative and summing to one within each row, and \(\eta_i\) is genetic drift, about which the least-squares estimators below need nothing except that its expectation vanishes. The coefficient \(A_{ij}\) is the proportion of individuals in population \(i\) that originated from population \(j\) during the interval, the fraction of \(i\) replaced by migrants from \(j\). Read epidemiologically it is the proportion of infections that \(i\) imports from \(j\), which is why the COVID work¹ calls \(A\) the importation-rate matrix and the ancient DNA work² the migration matrix.

The linear form is not an extra assumption. Under neutrality the frequency of any union of lineages must obey the same evolution equation, which makes the map additive over disjoint groups of lineages, and the conditional expectation of a neutral frequency is extensive: starting with half as many copies of the allele is expected to lead to half as many copies later. Additive and extensive leaves nothing but a linear map. Neutrality also fixes \(A\) to be the same matrix for every allele we track, since it describes who moved and not what they carried, and that is what allows many loci to be pooled into one estimate.

The same matrix can be read in the opposite direction of time, and that reading is worth having. Take a genome sampled in population \(i\) and follow its lineage of ancestors into the past. Each row of \(A\) sums to one over non-negative entries, which makes it a probability distribution over source populations, so \(A_{ij}\) is the probability that the lineage sat in \(j\) one interval earlier and \(A_{ii}\) the probability that it stayed put. Forward in time \(A\) is a table of migrant proportions; backward in time it is the transition matrix of a Markov chain on demes, and the history of a single lineage is one walk of that chain.

If we subtract \(X_i(t)\) from both sides, the mechanism is in plain view:

\[ X_i(t+\Delta t) - X_i(t) = \sum_{j=1}^n A_{ij}\,\bigl[X_j(t) - X_i(t)\bigr] + \eta_i \]

where the bracket shows that \(j\) influences \(i\) only when the two frequencies differ, and that a larger \(A_{ij}\) makes \(X_i\) converge on \(X_j\) faster. Migration is a spring and drift is the noise that keeps it from ever settling.

Time Points Versus Loci

The model is the easy part. What separates one application from the next is which statistic of the data survives the noise, and the two datasets we are comparing fail in opposite directions.

SARS-CoV-2 surveillance is dense in time and thin in loci. England was sequenced weekly across nine regions, so \(\Delta t\) is one week and there are dozens of consecutive time points, but a 30-kilobase genome under a strong clonal background offers only a few tens of approximately independent segregating alleles. Ancient DNA is the mirror image. The Allen Ancient DNA Resource³ genotypes individuals at up to 1.23 million positions, which is an enormous number of loci, and delivers them as a handful of skeletons per 300-year window, pseudohaploid, with the great majority of positions missing in any one sample and the dates themselves uncertain by decades to centuries.

Frequencies We Can Infer

Suppose for a moment that the true frequencies are known. Then the estimator writes itself: minimizing the squared difference between predicted and observed frequencies over all right-stochastic matrices gives

\[ A^{\rm (LS)} = \operatorname*{arg\,min}_{A} \sum_{i,\,\mu,\,t} \Bigl[ X^\mu_i(t+\Delta t) - \sum_j A_{ij} X^\mu_j(t) \Bigr]^2 \]

where \(\mu\) labels the independent loci and the minimization runs subject to \(A_{ij}\ge 0\) and \(\sum_j A_{ij}=1\). Summing over a window of time points regularizes the fit and sets the resolution at which \(A\) is allowed to vary. This convex problem is the target in both applications, and the rest of this section and the next are the two ways of reaching it when the frequencies cannot be taken at face value.

The estimator wants the true frequencies and we only ever see a noisy sample of them. The remedy in the COVID work is to stop treating what we measured as the state of the system. A hidden Markov model keeps the true frequencies as hidden states and gives genetic drift and sampling their own noise terms, with drift variance set by an effective population size and sampling variance set by the number of sequences and a per-region overdispersion. Gaussian noise makes the model a Kalman filter with an analytic likelihood, and the parameters follow from MCMC or, more cheaply, from expectation-maximization. Figure 1 shows the structure of this model.

Fig. 1 State-space model: hidden frequencies, noisy observations

true frequencies X(t) in every region, unobserved sampled frequencies Y(t), observed

Standard graphical model notation: open nodes are hidden, shaded nodes are observed. The hidden state X(t) is the whole vector of true frequencies, one per region. Every transition is the same step of Eq. 1, migration by A plus drift η, labeled once on the diagram; A is what couples the regions and it is what we are after. The drift η has variance x(1−x)/N_e per step, and observing the frequency in S sequences adds a further c x(1−x)/S, with c a per-region overdispersion. With Gaussian approximations in both terms the model becomes a Kalman filter, and its analytic likelihood makes parameter estimation cheap.

Distances We Can Trust

None of that survives contact with ancient DNA. A few genomes in a 300-year bin do not determine a frequency, the noise at that depth is not remotely Gaussian, and each sample covers a different subset of sites. We give up on individual frequencies, and write the dynamics directly in terms of quantities that many loci can estimate together.

The \(F_2\) statistic⁴ is the simplest of them: the squared difference in allele frequency between two populations, averaged over the whole genome. No single site is measurable at this sample size, but the average over a million of them is, and because every site contributes its own unbiased term, samples covering different subsets of positions still combine into one number. Under neutrality its expectation grows with the drift that separates two populations, which is why the field reads it as a genetic distance and builds admixture graphs out of it.

Write \(F_{ik}=\langle (X_i-X_k)^2\rangle\) for that distance and \(F'_{ik}=\langle (X'_i-X_k)^2\rangle\) for its lagged version, where the prime denotes time \(t+\Delta t\) and the brackets are the average over sites. The lag is the addition here, and it is what tells a donor from a recipient.

Figure 2 runs both statistics in the browser for two demes, \(i\) and \(k\), where the two sliders are the only off-diagonal entries \(A\) has. The top panel is the relaxation of the model itself, the neutral frequencies in \(i\) and \(k\) pulled together by migration, and the bottom panel is the distance between them with its two lagged versions, \(F'_{ik}\) and \(F'_{ki}\). Symmetric migration pulls all three curves down together. Set the sliders so that flow runs in one direction only and the two separate, because only the receiving population has moved toward where the donor used to be.

Fig. 2 Relaxation, and the asymmetry that dates it

i ← k 0.005 k ← i 0.040 Resample

population i population k F_ik F′_ik F′_ki

240 neutral loci in two demes of effective size 500, six of them drawn in the top panel. The sliders set the fraction of each population replaced by migrants per generation, which for two demes are the two off-diagonal entries of A; the lag of the primed statistics is ℓ = 4 generations. With k importing from i and not the reverse, k moves toward where i used to be, so F′_ki falls below F_ik while F′_ik does not. Set the two sliders equal and the asymmetry disappears.

Squaring the lagged difference and averaging over loci turns the relaxation form of the model into a relation between distances alone²:

\[ F'_{ik} - F'_{ii} = \sum_{j=1}^n A_{ij}\,\bigl(F_{jk} - F_{ji}\bigr) \]

where \(F'_{ii}\) is the mean squared displacement of population \(i\) over the interval, every quantity on both sides is a measurable genetic distance, and the dependence on \(A\) is still linear. Genetic drift and sampling noise both drop out of the expectation, because they shift the two sides of the difference by the same amount, once the rows of \(A\) sum to one. This is what makes the statistic usable at a few genomes per bin: the relation holds for the noisy estimates unchanged, with no correction to apply and nothing about the sampling to model.

From there \(A\) is fit as it was before, by least squares on the residuals of this relation over all \(i\), \(k\) and \(t\), subject to the same non-negativity and row-sum constraints: the same convex program, run on distances instead of frequencies.

Wherever Something Is Sampled Twice

The two analyses are one model and two estimators, run at opposite ends of a single axis. SARS-CoV-2 gives many time points and few independent alleles, so the noise has to be modeled, and a Kalman filter pulls the true frequencies out from under drift and sampling. Ancient DNA gives a million loci and almost no time points, so the noise has to be cancelled, and the difference of two lagged distances does it without a noise model at all. What comes back is the same object either way. In England the networks track geography but carry long-range links stronger than mobility data predicts, and they change between variant waves¹. In Europe the Steppe ancestry that reaches Britain in the third millennium arrives by way of Central Europe rather than directly: a route with a direction and a date².

Nothing in the derivation knows that the populations are places. It asks only for subpopulations, alleles nobody is selecting on, and samples at more than one time. Communities split by age or socioeconomic status rather than by geography would serve as well¹, as would any pathogen under routine surveillance, or the microbial populations that metagenomics reads repeatedly from the same hosts and sites, where dispersal between them does exactly what migration does here.

One limit travels with the method wherever it goes. It sees movement only through the divergence that movement destroys, so exchange between populations that already look alike leaves no imprint, however much of it there is. How much there is to read is settled upstream, by how the populations are carved up and how far apart in time they are sampled, long before any matrix is fit.

References

1. Okada et al., PNAS 122, e2500663122 (2025), the SARS-CoV-2 transmission networks, with a commentary in the same journal.
2. Isacchini et al., bioRxiv (2026), the ancient DNA analysis and the \(F_2\) formulation used here.
3. Mallick et al., Scientific Data 11, 182 (2024), the Allen Ancient DNA Resource, released through Harvard Dataverse.
4. Patterson et al., Genetics 192, 1065 (2012), the \(f\)-statistics.