ation data to test our approach. All simulation scripts are available at 
https://github.com/leospeidel/twigstats_paper.
f4-ratio admixture simulation. Our simulation in Fig. 1b and Extended 
Data Fig. 3b simulates five populations named PI, PO, P1, P2 and PX, 
where PO splits from all other populations 10,000 generations ago, P1 
and P2 represent two proxy source groups that split from each other at 
250 generations or 500 generations ago, PI splits from P1 100 genera-
tions ago and PX emerges from a pulse admixture between P1 and P2 
50 generations ago. All populations have a constant diploid population 
size of 5,000, a variable human-like recombination map, in which our 
simulation only covers chromosome 1, and a human-like mutation rate 
of 1.25 × 10−8 mutations per base per generation. We additionally have 
a modified simulation with a lower mutation rate of 4 × 10−9 mutation 
per base per generation, emulating a transversions-only dataset, and a 
simulation in which P2 has a diploid population size of 1,000 in the last 
50 generations, emulating a recent bottleneck in this population. We 
sample 20 haploid sequences from all populations. The ‘large sample 
size’ simulation samples 100 haploid sequences from all populations.
f4-ratio admixture simulation with genotype and phasing errors. 
We emulate the data quality we expect in imputed ancient genomes 
(Extended Data Fig. 3b). We implement a simple error model in which 
every haploid genotype at any segregating site can switch with a certain 
error probability. We can theoretically compute the predicted squared 
correlation coefficient (r2) between the true simulated genotypes and 
the genotypes that include error, stratified by minor allele frequency, 
to generate a plot similar to those used for evaluating imputation accu-
racy using downsampled high-coverage ancient genomes72 (Extended 
Data Fig. 3a). As imputation accuracy varies for each individual in real 
settings, we randomly sample the error probability for each individual 
uniformly between 1 × 10−4 and 1 × 10−3 (errors per SNP per haplotype). 
This yields r2 curves that are comparable to those observed in real data. 
We additionally simulate a high error case, for which we sample error 
probabilities between 1 × 10−3 and 1 × 10−2.
In real settings, we are additionally required to computationally 
phase genomes. We emulate this by combining two haploid sequences 
to construct a diploid individual. We then computationally rephase 
these diploid individuals without a reference panel. This approach is 
expected to result in suboptimal phasing and should therefore be well 
suited to test robustness to phase-switch errors.
qpAdm simulation. Our simulation in Extended Data Fig. 3c uses the 
simulation model and script provided with ref. 23, although we changed 
this script to use the human hotspot recombination map. We simulate 
only chromosome 1. In the original simulation model, admixing sources 
split 1,200 generations ago, with admixture occurring 40 generations 
ago. We additionally simulate a version in which all population split 
times and admixture times are reduced by a factor of 5. We sample 20 
haploid sequences per population.
Stepping-stone separation by distance simulation. We adapt the 
simulation model provided previously23 to simulate a stepping-stone 
model of nine populations organized on a 1D grid, in which individu-
als are able to migrate between adjacent populations (Extended Data 
Fig. 3d). We changed this script to use the human hotspot recombina-
tion map and simulate only chromosome 1. We simulate under migra-
tion rates of 0.001 and 0.005, corresponding to average FST values of 
0.01 and 0.002, respectively23. We sample 20 haploid sequences per 
population. We then fitted population 4 using pairs of other popu -
lations as sources in a rotational qpAdm scheme such that unused 
populations are assigned to the reference set.
We expect that this simulation model violates qpAdm assumptions