Adapting Gravel et al Recipe to Include Increased Growth Rate In Recent Epoch

[Grace Rhodes]

Hi everyone!

I'm working on a simulation where I would like to adapt the Gravel et al. (2011) recipe based on a paper from Tennessen et al. (2012) (link) which adapts the Gravel model to "allowing for a reduced initial European expansion that is compensated for by accelerated growth starting after the split of European and Asian populations. Similarly, we introduced a phase of exponential growth in the African population starting at the same time." Because I want the growth rate to be precise, I use the example on page 121 of the SLiM manual. However, I think perhaps I am missing something in adding the second period of exponential growth because the PCA plots I am getting from running PLINK on the resulting VCF file is not what I expect to see. Has anyone done something similar and have any advice? Thank you in advance!

Current script:

// 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]);

}

PCA:

(done in PLINK with linkage pruning)

[Chase W. Nelson]

Hi Grace! I'm too ignorant to comment on the downstream analyses with PLINK, but a potentially useful resource is the paper by Jouganous et al. 2017 (Gravel is last author) which estimates many of the demographic parameters for similarly complex demographic models — could be useful as a comparison. The Gravel and Tennessen supplementary info (which you've likely already consulted) also have potentially relevant info.

Regarding the model you've presented, I only had a chance to look briefly, but I do notice that the mutation rate of 10e-8 is approximately an order of magnitude higher than in the original model; that could be right, but perhaps worth double-checking, particularly if you're keeping the recombination rate down at 1e-8. Additionally, a fixed deleterious selection coefficient of -0.1 (-10%!) seems really extreme, even if fully recessive (d=0). It might be worth trying a neutral model first as a baseline for expected behavior? If I'm interpreting the code right, your first 1000 nucleotides (indices 0 through 999) will be totally mutation-free, after which chunks of 500-mutable, 500-immutable bases are evenly interspersed, for a total of 5,000 mutable positions, implying 0.0005 mutations per genome per generation. Because the mutations are fully recessive, they'll first segregate as if totally neutral, i.e., most will be lost to drift, but any which luckily increased in frequency would face a prohibitive barrier from homozygosity, i.e., selection might not allow them to reach even intermediate frequencies because their s = -0.1. Perhaps this is all obvious, but just wanted to make some observations in case they help thinking through the implications of these parameters for the downstream methods. Sorry I don't have any more specific insights!

Chase

[Peter Ralph]

I agree with Chase (thanks!).

But, let me add one more thing - "what's expected" from PCA is pretty hard to pin down. For instance, what you get out of PCA depends on sample sizes. What you're getting isn't totally unexpected to me: it's saying (roughly) that most variation is within AFR, which matches what we know to be true (and, in the model). If your sample sizes were less balanced (e.g., only 20 from AFR) then probably some of the top PCs would distinguish AFR from the others (which is maybe what you were expecting?). So - it's better to look first at closer-to-the-source summaries: e.g., are diversity and divergence what you expect?

--peter

---

From: slim-d...@googlegroups.com <slim-d...@googlegroups.com> on behalf of Chase W. Nelson <cwnel...@gmail.com>
Sent: Saturday, June 15, 2024 5:56 AM
To: slim-discuss <slim-d...@googlegroups.com>
Subject: Re: Adapting Gravel et al Recipe to Include Increased Growth Rate In Recent Epoch

 

--
SLiM forward genetic simulation: http://messerlab.org/slim/
---
You received this message because you are subscribed to the Google Groups "slim-discuss" group.
To unsubscribe from this group and stop receiving emails from it, send an email to slim-discuss...@googlegroups.com.
To view this discussion on the web visit https://groups.google.com/d/msgid/slim-discuss/e8e70699-c55b-4488-a079-ca7e5244802bn%40googlegroups.com.