// Model based on Tennessen et al. 2012
initialize() {
initializeMutationRate(10e-8);
initializeMutationType("m1", 0.0, "f", -0.1);
initializeGenomicElementType("g1", m1, 1.0);
for (index in 1:10)
initializeGenomicElement(g1, index*1000, index*1000 + 499);
initializeRecombinationRate(1e-8);
}
// INITIALIZE the ancestral African population of size 7310
1 early() { sim.addSubpop("p1", asInteger(round(7310.370867595234))); } // paper rounds to 7310
// END BURN-IN period of 10*N=73104 generations (specific to SLiM recipe); EXPAND the African population
// This occurs (5919.131117 generations)*(25 years)=147978 yr ago; paper rounds to 5920 gens (148000 yr)
// Thus, simulation should end at generation 1+73104+5919.131117=79024
73105 early() { p1.setSubpopulationSize(asInteger(round(14474.54608753566))); } // paper rounds to 14474
// SPLIT Eurasians (p2) from Africans (p1) and SET UP MIGRATION between them
// This occurs 2056.396652 generations (51409.9163 years) ago; paper rounds to 2040 gens (51000 yr)
// Relative to beginning, this is generation 79024-2056.396652=76968
76968 early() {
sim.addSubpopSplit("p2", asInteger(round(1861.288190027689)), p1); // paper rounds to 1861
p1.setMigrationRates(c(p2), c(15.24422112e-5)); // paper rounds to 15e-5
p2.setMigrationRates(c(p1), c(15.24422112e-5)); // paper rounds to 15e-5
}
// SPLIT p2 into European (p2) and East Asian (p3) subpopulations; RESIZE; SET UP MIGRATION between them
// This occurs 939.8072428 generations (23495.18107 years) ago; paper rounds to 920 gens (23000 yr)
// Relative to beginning, this is generation 79024-939.8072428=78084
78084 early() {
sim.addSubpopSplit("p3", asInteger(round(553.8181989)), p2); // paper rounds to 554
p2.setSubpopulationSize(asInteger(round(1032.1046957333444))); // reduce European size; paper rounds to 1032
// Set migration rates for the rest of the simulation
p1.setMigrationRates(c(p2, p3), c(2.54332678e-5, 0.7770583877e-5)); // paper rounds to c(2.5e-5, 0.78e-5)
p2.setMigrationRates(c(p1, p3), c(2.54332678e-5, 3.115817913e-5)); // paper rounds to c(2.5e-5, 3.11e-5)
p3.setMigrationRates(c(p1, p2), c(0.7770583877e-5, 3.115817913e-5)); // paper rounds to c(0.78e-5, 3.11e-5)
}
// SET UP (mild) EXPONENTIAL GROWTH in Europe (p2) and East Asia (p3)
// Where N(0) is the base subpopulation size and t = gen - 78084:
// N(Europe) should be int(round(N(0) * (1 + 0.003784324268)^t)), i.e., growth is r=0.307% per generation
// N(East Asia) should be int(round(N(0) * (1 + 0.004780219543)^t)), i.e., growth is r=0.48% per generation
78084:78819 early() {
t = sim.cycle - 78084;
p2_size = round(1032.1046957333444 * (1 + 0.00307)^t); // change from Gravel, reduced initial exponential growth
p3_size = round(553.8181989 * (1 + 0.004780219543)^t); // paper rounds to N(0)=554 and r=0.0048
p2.setSubpopulationSize(asInteger(p2_size));
p3.setSubpopulationSize(asInteger(p3_size));
}
// set up explosive exponential growth in Europe (p2) and Africa (p1)
// in Tennessen model, r=1.95% per generation in Europe and r=1.66% per generation in Africa (same as last period for p3)
78819:79024 early() {
t = sim.cycle - 78819;
p2_size_new = round(1032.1046957333444 * (1 + 0.0195)^t); //increase growth in Europe pop
p1_size = round(14474.54608753566 * (1 + 0.0166)^t); //add growth in African pop
p3_size_new = round(553.8181989 * (1 + 0.004780219543)^t); // East Asia same growth as before
p2.setSubpopulationSize(asInteger(p2_size_new));
p1.setSubpopulationSize(asInteger(p1_size));
p3.setSubpopulationSize(asInteger(p3_size_new));
}
// OUTPUT AND TERMINATE
// Generation 79024 is the present, i.e., 1 initialize + 73104 burn-in + 5919 evolution
79024 late() {
sample1 = sample(p1.genomes, 1000);
sample2 = sample(p2.genomes, 1000);
sample3 = sample(p3.genomes, 1000);
c(sample1, sample2, sample3).outputVCF([filepath]);
}