Unexpected gene drive behaviour in simple panmictic model

[Jack Sissons]

Hi Ben and drive modellers,

I have a model of hermaphroditic broadcast spawner to test gene drive dynamics and I'm experiencing unexpected model behaviour. I've found that the drive, with neutral fitness, can reach very high frequencines (80%+) in the population and still be lost in subsequent generations. I've also observed that implementing serial releases of drive carriers actually suppresses (!) the proliferation of the drive, and causes the frequency to remain low. To get to the bottom of this, I've pared back the model to what is pasted below. As a result, I've narrowed down the activity to being caused by the reproduction function. If I execute random mating (red text below, gene drive has no effect), then the drive fixes immediately. However, when I implement an effect of the drive, i.e., carriers cannot produce eggs, with the script block below, the behaviour gets weird. 

reproduction() {

if (individual.tagL0==F) {

mates = sample(p1.individuals, 2);

for (mate in mates);

subpop.addCrossed(individual, mate);

}

}

The intention here is that only WT individuals produce eggs (are parent1) and they choose random mates (parent2) from the entire population (WT or GD+). The drive copies itself in the child and they become homozygotic for the drive. What I observe now is high stochasticity in the mutation frequency trajectory, sometimes it fixes, sometimes it is lost, and sometimes it fluctuates all over the place for 5000 generations. I would expect it to be a rare event for the drive to then extinguish from frequencies >50%. I've checked several things:

- carriers are indeed homozygotes

- carriers aren't preferentially dying (i've tried excluding them from age-related moratlity)

- tagL0 is written for all individuals

- varying number of offspring & mates doesn't drastically change the profile

I'm starting to wonder if this is a feature, not a bug, but to me it seems rather buggy. It seems that there is a bias for parent1 to choose WT individuals to be parent2 over drive carriers. If you have any insight into the implentation or behaviour here that would be greatly appreicated!

I'm running SLiM 5.1 and I get the same behaviour on both mac and windows.

Thanks, 

Jack

Full script:

initialize() {

initializeSLiMModelType("nonWF");

defineConstant("K", 1000);

initializeMutationRate(1.3e-6); // scaled

initializeMutationType("m1", 0.5, "f", 0.0); // neutral

initializeMutationType("m2", 0.5, "f", 0.0); // Gene drive

initializeMutationType("m3", 0.5, "f", 0.0); // S-locus mutation

m2.convertToSubstitution = T;

initializeGenomicElementType("g1", m1, 1.0);

initializeGenomicElementType("g2", m3, 1.0); // S-locus mutations

initializeGenomicElement(g1, 0, 999);

initializeGenomicElement(g1, 1001, 20000);

initializeGenomicElement(g2, 20001, 21000); // S-locus

initializeGenomicElement(g1, 21001, 99999);

initializeRecombinationRate(5e-6); //scaled

}

1 early() {

sim.addSubpop("p1", 100);

}

1: early() {

sim.killIndividuals(p1.subsetIndividuals(minAge=9));

}

early() {

p1.fitnessScaling = K / p1.individualCount;

}

2: early () {

inds = p1.individuals;

inds.color = ifelse(inds.haploidGenome1.containsMarkerMutation(m2, 1000), "red", "green");

}

reproduction() {

//subpop.addCrossed(individual, sample(p1.individuals, 1));

///// gene drive carriers produce no eggs

if (individual.tagL0==F) {

mates = sample(p1.individuals, 1);

for (mate in mates);

subpop.addCrossed(individual, mate);

}

}

120 late() {

driveSeeders = p1.individuals[0:499];

driveSeeders.haplosomes.addNewDrawnMutation(m2, 1000);

catn(size(driveSeeders) + " carriers released");

}

121: late() {

catn(sim.cycle + ": " + 100*(sim.mutationFrequencies(NULL, sim.mutationsOfType(m2))) + "% carriers");

}

120:10000 late() {

if (sim.countOfMutationsOfType(m2) == 0)

{

fixed = any(sim.substitutions.mutationType == m2);

cat(ifelse(fixed, "FIXED\n", "LOST\n"));

catn("in " + (sim.cycle - 120) + " generations");

sim.simulationFinished();

}

}

modifyChild() {

mut = sim.mutationsOfType(m2);

if (size(mut) == 1)

{

hasMutOnChromosome1 = child.haploidGenome1.containsMutations(mut);

hasMutOnChromosome2 = child.haploidGenome2.containsMutations(mut);

if (hasMutOnChromosome1 & !hasMutOnChromosome2)

child.haploidGenome2.addMutations(mut);

else if (hasMutOnChromosome2 & !hasMutOnChromosome1)

child.haploidGenome1.addMutations(mut);

}

return T;

}

1: late() {

p1.individuals.tagL0 = p1.individuals.haploidGenome1.containsMarkerMutation(m2, 1000) & p1.individuals.haploidGenome2.containsMarkerMutation(m2, 1000);

}

[Ben Haller]

Hi Jack!

Before delving into this too carefully, I notice you have a typo in your script:

for (mate in mates);

subpop.addCrossed(individual, mate);

You've got an extra semicolon at the end of the first line.  So the loop does not generate two mates, as intended; it loops first, doing nothing (the semicolon), and then it calls addCrossed() just once, making one offspring with the second mate in mates.  Is that the problem?

Cheers,
-B.

Benjamin C. Haller
Messer Lab
Cornell University

--
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 visit https://groups.google.com/d/msgid/slim-discuss/fefbbde5-dacc-44d9-b1b8-9a6bc24d8bb3n%40googlegroups.com.

[Jack Sissons]

Hey Ben, 

My mistake, good spotting. I'm away from my machine for a while but I don't think this is the issue. Drawing a single mate from a majority of drive carriers ought to still lead to fixation. Thanks for taking a look!

[Isabel Kim]

Hey Jack,

I am a recent PhD graduate from the Messer lab who focused on gene drive modeling. I took at a look at your code and don't see any bugs besides the semicolon that Ben pointed out. I think the issue is that you're simulating embryo conversion (a child of a w/w x d/d cross is converted from d/w to d/d immediately). With a major fitness cost in drive homozygotes (only reproducing at half the rate of wild-types), you need heterozygotes (with equal fitness to wild-types) to be able to build up in the population and spread the drive allele. The way you have it now, the drive is initialized such that 50% of the population is d/d and 50% is w/w. In the next generation, 50% of offspring will be from w/w x w/w crosses and be w/w and 50% of offspring will be from w/w x d/d crosses and, due to perfect embryo conversion, be d/d. Thus, you're stuck at this gridlock of 50/50 drive homozygotes to wild-types in the population, and the population can only break out of it through random fluctuations, many of which will cause the drive allele to be lost.

I added a boolean parameter, "GERMLINE_CONVERSION" below. When this is T, an offspring of a w/w x d/d cross is a heterozygote and reproduces normally, but due to perfect germline conversion, they only pass down the drive allele. In all the simulation runs I've done with this turned on, the drive always fixes.

Best,

Isabel

initialize() {

initializeSLiMModelType("nonWF");

defineConstant("K", 1000);

defineConstant("GERMLINE_CONVERSION", T); // if F, default embryo conversion

if (!GERMLINE_CONVERSION){

hasMutOnChromosome1 = child.haploidGenome1.containsMutations(mut);

hasMutOnChromosome2 = child.haploidGenome2.containsMutations(mut);

if (hasMutOnChromosome1 & !hasMutOnChromosome2)

child.haploidGenome2.addMutations(mut);

else if (hasMutOnChromosome2 & !hasMutOnChromosome1)

child.haploidGenome1.addMutations(mut);

} else {

// handle germline conversion in mother and father separately

mother_passed_wt = !(child.haploidGenome1.containsMutations(mut));

mother_drive_on_1 = parent1.haploidGenome1.containsMutations(mut);

mother_drive_on_2 = parent1.haploidGenome2.containsMutations(mut);

if (mother_passed_wt & (mother_drive_on_1 | mother_drive_on_2))

// wt chrom was converted to drive in germline

child.haploidGenome1.addMutations(mut);

father_passed_wt = !(child.haploidGenome2.containsMutations(mut));

father_drive_on_1 = parent2.haploidGenome1.containsMutations(mut);

father_drive_on_2 = parent2.haploidGenome2.containsMutations(mut);

if (father_passed_wt & (father_drive_on_1 | father_drive_on_2))

// wt chrom was converted to drive in germline

child.haploidGenome2.addMutations(mut);