Prevention of outcrossing ability depending on the genotype of the individual

[Yujean Kim]

Hi!

Currently I am trying to encode a haploinsufficient population suppression gene drive scenario. I am trying to introduce a mutant allele which affects the ability of snails to produce sperms. I want to prevent outcrossing between two different snails when parent 1 has one or more mutant allele and when parent 2 has one or more mutant allele. I have tried to encode this through using the following code in my modifyChild() callback:

if (parent1.countOfMutationsOfType(m2) >= 1) {

     if (parent2.countOfMutationsOfType(m2) >=1)

          return F;

}

Theoretically, the snail population should decrease as the total mutant allele frequency in the population increases. However, my snail population remains constant. I was wondering if I did something wrong and if there is an alternative approach to take.

Kind regards,

Yujean Kim

[Ben Haller]

I don't see anything obviously wrong with that code, but then I'm not seeing the rest of the script for context.  There are lots of debugging tools available; you could insert print() statements to check what's going on, use the Eidos console to look at the state of objects, and use debug points in SLiMgui to get debugging logs for key lines of code as the model runs.  You could also add assert() calls to assert things that should be true if the model is working correctly, and see if any of those get hit; and of course graphical visualization in SLiMgui can be useful (what does a frequency plot show that your drive mutation is doing?).  Check your assumptions and do detective work.  Maybe your drive allele is getting lost from the population, and thus no longer affects dynamics?  Maybe it is at such low frequency that drawing two parents that both carry it is so rare that it doesn't matter – perhaps because the drive is not copying itself properly?  Maybe the modifyChild() callback isn't even being called, because it is not active or does not apply to the subpopulation in question?  Maybe other logic in the modifyChild() callback is interfering with the correct operation of this code?  And then of course you have not said whether this is a WF or a nonWF model; in a WF model the population size is fixed, and returning F from a modifyChild() callback will simply draw new parents and try again, so you will only see population suppression with a nonWF model.  In other words, there are a million reasons why this problem could be happening; with a little investigation in SLiMgui, and some reading of the documentation, you can probably figure it out.  If not, feel free to post the full model here – but please try to figure it out for yourself first.  :->  Good luck, and happy modeling!

Cheers,
-B.

Benjamin C. Haller
Messer Lab
Cornell University

[Yujean Kim]

Hi Ben,

Thank you for your email! I have tried multiple ways to figure this out during the past 12 hours but have not been able to understand why. Could you please help me out? I have highlighted the specific code where I affect outcrossing ability in red!

initialize() {

   defineConstant("cut_rate", 0.98); //rate of double strand break occuring

   defineConstant("resistance_rate", 0.078); //rate of repair with NHEJ

   defineConstant("mutation_rate", 0); //drive specific mutation rate

   defineConstant("nonres", 10000); //non-resistant allele (target population)

   defineConstant("capacity" , 10000);

   defineConstant("drive_size", 1000); //gene drive release

   defineConstant("growth_rate", 53);

   defineConstant("max_offspring", 14000);

   defineConstant("mortality", c(1,0));

   defineConstant("selfing_rate", 0.8);

   defineConstant("inbreeding_depression", 0.063);

   

   initializeSLiMModelType("nonWF");

   

   //Assume no fitness cost

   

   awt = initializeMutationType("m1", 0.5, "f", 0.0); //wildtype fertility allele

   adr = initializeMutationType("m2", 0.5, "f", 0.0); //mutant fertility allele (drive allele)

   ar2 = initializeMutationType("m3", 0.5, "f", 0.0); //resistant allele

   

   avariations = c(awt, adr, ar2);

   

   //Defining different alleles into one locus    

   

   initializeGenomicElementType("g1", avariations, c(1, 1, 1)); //locus a

   

   //Defining each locus into a chromosome position

   

   initializeGenomicElement(g1, 0, 0); //chromsosome position 0

   

   avariations.mutationStackPolicy = "l"; //makes sure last mutation added is the most recent

   avariations.mutationStackGroup = 1; //makes sure there is only one mutation/allele in one locus

   

   initializeMutationRate(0);

   initializeRecombinationRate(0.5);

}

//Reproduction rules - ones with no wild type alleles are infertile (doublesex genes are haplosufficient)

reproduction() {

   

   //Calculation for carrying capacity

   

   capacity_fitness_scaling = growth_rate / (((growth_rate-1) * (p1.individualCount / capacity)) + 1);

   p = capacity_fitness_scaling * 2 / max_offspring ;

   num_offspring = rbinom(1, max_offspring, p); //number of offspring depends on population and carrying capacity

   

   for (i in seqLen(num_offspring)) {

       if (subpop.size())

           if (runif(1) < selfing_rate)

               offspring = subpop.addSelfed(individual); //Add selfing rate

           else

               offspring = subpop.addCrossed(individual, subpop.sampleIndividuals(1)); //Add outcrossing rate

   }

}

//Gene drive implementation in offspring in the germline - this is a double drive design

//convert to single gene drive

1: modifyChild() {

   parent = c(parent1, parent2);

   child_chromosome = c(childGenome1, childGenome2);

   for (i in c(0,1))

       if (sum(parent[i].genomes.countOfMutationsOfType(m2)) & sum(child_chromosome[i].countOfMutationsOfType(m1))) {

           if ((runif(1) < cut_rate) & (runif(1) > resistance_rate)){

               if (runif(1) < mutation_rate)

                   child_chromosome[i].addNewDrawnMutation(m3, 0);

               else

                   child_chromosome[i].addNewDrawnMutation(m2, 0);

           }

       }

       if (sum(parent[i].genomes.countOfMutationsOfType(m2)) & sum(child_chromosome[i].countOfMutationsOfType(m1))) {

           if ((runif(1) < cut_rate) & (runif(1) < resistance_rate))

               childGenome2.addNewDrawnMutation(m3, 0);

   }

   

   return T;

   

   if (individual.countOfMutationsOfType(m2) >= 1) {

       if (parent1 == parent2)

           return F;

       return T;

   }

   if (individual.countOfMutationsOfType(m3) >= 1) {

       if (parent1 == parent2)

           return F;

       return T;

   }

   if (individual.countOfMutationsOfType(m1) == 0) {

       if (parent1 == parent2)

           return F;

       return T;

   }

   if (parent1.countOfMutationsOfType(m2) >= 1) {

       if (parent2.countOfMutationsOfType(m2) >= 1)

           return F;

   }

   if (parent1.countOfMutationsOfType(m2) >= 1) {

       if (parent2.countOfMutationsOfType(m3) >= 1)

           return F;

   }

   if (parent1.countOfMutationsOfType(m3) >= 1) {

       if (parent2.countOfMutationsOfType(m3) >= 1)

           return F;

   }

   if (parent1.countOfMutationsOfType(m3) >= 1) {

       if (parent2.countOfMutationsOfType(m2) >= 1)

           return F;

   }

   if (parent1.countOfMutationsOfType(m1) < 2) {

       if (parent2.countOfMutationsOfType(m1) < 2)

           return F;

   }

   

   all = sim.subpopulations.individuals;

   for (ind in all){

       if (parent1 == parent2){

           if (runif(1) < inbreeding_depression)

               ind.fitnessScaling = 1 - inbreeding_depression;

       }

   }

}

//fitness scaling - discrete generations (non-overlapping generations)

late() {

   all = sim.subpopulations.individuals;

   for (ind in all) {

       age_mortality_rate = mortality[ind.age];

       ind.fitnessScaling = 1 - age_mortality_rate;

   }

}

//introduce of drive allele

1 {

   sim.addSubpop("p1", nonres);

   p1.individuals.genomes.addNewDrawnMutation(m1, 0);

   

   sim.addSubpop("p2", drive_size);

   p2.individuals.genomes.addNewDrawnMutation(m2, 0);

   

   drop = p2.individuals;

   p1.takeMigrants(drop);

}

late() {

   

   //Sum of each mutation type in population

   

   all = sim.subpopulations.individuals;

   num_awt = sum(all.genomes.countOfMutationsOfType(m1) > 0);

   num_adr = sum(all.genomes.countOfMutationsOfType(m2) > 0);

   num_ar2 = sum(all.genomes.countOfMutationsOfType(m3) > 0);

   

   //Propotion of each mutation type in population

   

   rate_awt = num_awt / (2 * size(all));

   rate_adr = num_adr / (2 * size(all));

   rate_ar2 = num_ar2 / (2 * size(all));

   

   //output

   cat(asInteger(sim.generation) + "\t " + rate_awt + "\t " + rate_adr + "\t " + rate_ar2 + "\t " + size(all) + "\n");

   

   

   if (rate_awt == 0.0) {

       catn("OUTCOME: DRIVE SUCCESS");

       sim.simulationFinished();

   }

   if (rate_adr == 0.0) {

       catn("OUTCOME: DRIVE LOST");

       sim.simulationFinished();

   }}

//// MAIN END CONDITION.

111 late() {

   catn("OUTCOME: DRIVE PERSISTS");

   sim.simulationFinished();

}

[Ben Haller]

Hi Yujean.  It looks like your script has a stray "return T;" in the modifyChild() callback before the point where the code checks for the condition that would return F.  I think you might have been confused by the indentation of your code, which appears to be wrong; the close brace that is indented as if it closes the for loop actually closes the compound statement on an if statement, and that if statement is indented as if it is inside the for loop, but in fact it is not.

So, just so you and others know about the tools that could have caught this mistake for you:

- SLiMgui can prettyprint your code for you, which would have immediately revealed the indentation problems that were perhaps misleading you; click the prettyprint button (third from the left on the "Input Commands" row of buttons) or choose the corresponding menu item.  If you hold down the "option"/"alt" key it will completely reformat your code for you, which is often even more useful.

- SLiMgui's "debug points" facility would have immediately indicated to you that the code you were interested in was never being hit.  Just click in the left gutter of the editor to set a debug point on the relevant line of code, and then run the model.  If you open the debug output window, you will then see that the debug point is never hit – no debug output is generated.  This would lead you directly to that "return T;" line.

- A print statement inserted in the relevant code would also not have been hit, if you prefer that to using debug points.

I don't mean to be unkind about this; debugging is a skill that takes time and effort to develop.  But please do work to develop it; it will pay off handsomely.  When I was first starting out, I hit bugs more than once that took me several *months* to find; perhaps you can imagine how frustrating that was!  :->  But it has been a couple of decades since that happened to me, because I developed the necessary skills to find my bugs, and you can too.  I would suggest that you try out prettyprinting and debug points now, with the version of your model that has the bug, to see how to use them to find the problem.  And next time, don't resort to asking the list until you've spent at least a full week trying first; you learn the skills, sadly, by beating your head against the wall until you figure it out.  Such is life.  :->

Good luck, and happy modeling!

Cheers,
-B.

Benjamin C. Haller
Messer Lab
Cornell University

Yujean Kim wrote on 4/27/22 9:14 PM:

--
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/81a5529e-d3d4-49e8-a369-77198ac385cfn%40googlegroups.com.