Thoughts about gamete selection

[bguo068]

SLiM is so powerful for its flexibility and realisticity. I have been thinking about using SLiM+pyslim+msprime to simulate different kinds of selection after skimming the comprehensive SLiM manual and the excellent tutorial in pyslim/msprime documentation.

When it comes to the species I am currently working with - plasmodium - which presents as a haploid for the majority of its complicated life cycle but as a diploid only briefly during its sexual stage, I ran into a simulation case very different from the typical cases of simulating recombining diploid (for which SLiM is designed for) and simulating nonrecombining haploid (for which SLiM can still deal with by removing the mutations from one genome of each diploid individual).

The main features of Plasmodium that cause issues for simulation are

• Plasmodium experience sexual reproduction during which recombination can occur. So for this stage, we would like to simulate as diploids (with recombination), and thus it is different from simulating the nonrecombining haploid;
• In the majority part of its life cycle, Plasmodium presents as haploid, which means the selection stress is acting on its haploid form (in my opinion, it is similar to gamete selection). Thus the mating process is affected not by the diploid fitness but instead by the haploid fitness.

Based on my superficial understanding of SLiM, it could be challenging to simulate plasmodium selection with these features considered. I was imagining some naive ways of doing gametes selection so that simulating selection in Plasmodium could be realistic:

1. Simulating gamete selection directly by rejecting offspring proposals with a probability based on gamete fitness​.
2. Simulating recombination between a pair of (fake) diploid individual​s. We can manually maintain fake diploids (each has an effective genome and an empty genome. Fitness is calculated based on the non-empty genome (gamete fitness)​. When generating offspring, based on the breakpoints on one parent, manually recombine the non-empty genomes of the two parents.​

Both methods could work in theory but could be very slow. I am wondering if there is a better way of implementing this type of simulation - "selection in recombing haploid" just with Eidos scripting or it is necessary to dive deep into the C++ code.

Thanks,
Bing

[bhaller]

Hi! Actually, this sort of thing is now pretty straightforward with facilities that will be released in SLiM 3.7 (which should be out before Christmas). I've attached two excerpts from the draft manual for 3.7. Together, they contain three sections; the sections on alternation of generations and haplodiploidy are probably the most central for you.

SLiM_Manual_excerpt.pdf
SLiM_Manual_excerpt2.pdf

So, have a look through those and then let me know if you have questions remaining. It doesn't show exactly your system, of course, but I think things would be very similar for you – using addRecombinant() to produce whatever type of offspring you want from whatever type of individuals are involved, using haploidDominanceCoeff to control how selection works in the haploids, etc.

[bguo068]

That is awesome!! I think the addRecombinant() method and haploidDominanceCoeff attribute could be the exact functionality I am looking for. I am more familiar with the WF model but I guess a nonWF model is required as reproduce() callback only allowed in nonWF model. I am able to compile the SLiM GitHub repo following manual 3.6 and will try to learn nonWF model and test addRecombinant() method and haploidDominanceCoeff for my simulation cases. Thank you so much!

[bhaller]

Oh, and also re: recombination between haploids, maybe this is relevant:

SLiM_Manual_excerpt3.pdf

This recipe is also in the 3.6 manual, but I think it maybe had some changes in it due to the changes in SLiM 3.7. I definitely very strongly urge you to develop this model in SLiM 3.7; a lot of improvements have been made in the modeling of unusual mating systems since 3.6. You can build the current GitHub head and it ought to work, although there are some issues left to shake out surrounding tree-sequence recording and tskit and pyslim.

[bguo068]

Wow, addRecombinant was already in SLiM 3.6. Thanks for the updated manual section!

[bhaller]

Yes, addRecombinant() was in SLiM 3.6, but it worked fairly differently. The way that the empty genomes of haploids are managed, the way that allele frequencies are computed, and other such things have been greatly improved, and that entailed some breaks in backward compatibility. So again, I would urge you not to use SLiM 3.6 for this work; it will not work nearly as well, and if you then try to migrate to 3.7, the changes in behavior will be a confusing source of bugs. Just start working on 3.7 to begin with. If you're interested in the history here, see #205, especially the long comment near the end.

[bguo068]

Thanks for linking the history. The changes about handling haploids look great!! I was tripped by simulating recombining haploids (wrongly following manual section 14.9) and experienced some strange behaviors. I will definitely start to work with v3.7 and follow up here once I have more experience with it. :)

[bguo068]

I think it works with at least the following balancing selection on recombining haploids.
This allele freqency trajectory is like this:

One thing I still don't quite understand is that the balancing selection will not work as I would imagine when I try to kill all parents in early events (try to mick manual v3.6 section 16.15 implementing a Wright-Fisher model with a nonWF model). In this case, I got allele frequency like this:

Any comments on this or suggested readings related to this behavior? Thanks.

// SLiM v3.7
initialize(){
	initializeSLiMModelType("nonWF");
	initializeTreeSeq();
	initializeMutationRate(0.0);
	initializeMutationType("m1", 0.5, "f", 0.0);
	initializeMutationType("m2", 1.0, "f", 0.2);
	m2.haploidDominanceCoeff = 1.0;
	initializeGenomicElementType("g1", m1, 1.0);
	initializeGenomicElement(g1, 0, 100000 - 1);
	initializeRecombinationRate(1e-8);
}

reproduction(){
	breaks = sim.chromosome.drawBreakpoints(individual);
	mate = p1.sampleIndividuals(1, exclude=individual);
	p1.addRecombinant(individual.genome1, mate.genome1, breaks, NULL, NULL, NULL);
}

1 early(){
	// add initial population size
	sim.addSubpop("p1", 1000, haploid=T);
	p1.individuals.age = 0;
	// add mutations
	pos = 50000;
	target = sample(p1.individuals.genome1, 5);
	target.addNewDrawnMutation(m2, pos);
}

early(){
	// kill older generations
	// inds = p1.individuals;
	// inds[inds.age > 0].fitnessScaling = 0.0;
	// adjust subpop fitness scaling factor to keep population size around N.
	p1.fitnessScaling = 1000 / p1.individualCount;
}

fitness(m2){
	// simple negative frequncy-dependent selection
	return 1.5 - sim.mutationFrequencies(p1, mut);
}

2000 late(){
	sim.simulationFinished();
	sim.treeSeqOutput("slim_out.trees");
}

[bhaller]

Ha, so, this is a bit subtle. I added a line to the early() event that prints:

catn("generation " + sim.generation + ": killing " + sum(inds.age > 0) + " of " + size(inds) + " individuals.");

It shows:

generation 1: killing 0 of 1000 individuals.
generation 2: killing 1000 of 2000 individuals.
generation 3: killing 496 of 992 individuals.
generation 4: killing 496 of 992 individuals.
generation 5: killing 496 of 992 individuals.
generation 6: killing 496 of 992 individuals.
generation 7: killing 496 of 992 individuals.
generation 8: killing 496 of 992 individuals.
...

So the model quickly equilibrates into a mode where there are, in this specific run (it will vary slightly), exactly 496 parental individuals. The reproduction() callback generates exactly one offspring per parental individual, so that gives you 992 individuals. Then the age limit kills off all the parents, so you're back to 496. And then the density-dependence doesn't kill anyone at all, since the population size (still measured at 992, at this point in the code, since fitnessScaling has not yet taken effect) is less than 1000.

Why does this mean that the frequency doesn't change substantially, and doesn't balance out to ~0.5? Well, it's basically because (a) fecundity does not depend upon fitness – every individual is having exactly one offspring as a first parent, and an average of one offspring as a mate, and (b) survival also does not depend upon fitness – every individual is surviving, regardless of the fitness effect returned by the m2 fitness() callback. So the mutations get passed on from parent to offspring without selection, and thus change very little in frequency; it's just a model of neutral drift, really.

Point (b) needs further explanation: why does survival not depend upon fitness? So, first of all, note that the selection coefficient provided for m2 mutations, 0.2, does not matter at all; it is completely overridden by the fitness() callback. In such cases I suggest using 0.0 in the code to make it clear that the value does not matter. So, what does the fitness() callback do? It says "If the frequency of an m2 mutation is 0.0, make the fitness of that m2 mutation 1.5; if its frequency is 1.0, make its fitness 0.5; and shade continuously in-between." So when the m2 mutation is at low frequency, it is supposed to rise to higher frequency by having a higher-than-1.0 fitness. But fitness values greater than 1.0 don't matter in nonWF models; you can't have a survival probability greater than 1.0, so all final fitness values >= 1.0 have an identical effect of guaranteed survival. Section 16.4 of the manual has a long discussion of this; it is perhaps the single most important thing to really internalize when doing nonWF models.

So, how to fix this? The simplest way is to increase the fecundity:

reproduction(){
	for (i in seqLen(rpois(1, 5.0)))
	{
		breaks = sim.chromosome.drawBreakpoints(individual);
		mate = p1.sampleIndividuals(1, exclude=individual);
		p1.addRecombinant(individual.genome1, mate.genome1, breaks, NULL, NULL, NULL);
	}
}

Now there are lots of excess individuals, which means that density-dependence sets the population fitnessScaling to be well below 1.0. As a result, the beneficial effect of m2 mutations at low frequency now does matter. And so the balancing selection works.

I'd urge you to work through this very carefully to make sure you understand the problem and the solution completely. It is not necessarily the right solution, it is just a solution. In practice, you want to think about the biology of your system very carefully and figure out how your system really works, in terms of fitness, selection, density-dependence, fecundity, etc., and model that. If you do that properly, such problems should not occur.

Good luck!

[bguo068]

Thank you for your patient and the excellent explanation.

I think I got it by explicitly caculating the absoluate fitness ( applying the subpopulation fitnessScaling and individual fitnessScaling to relativeFintess returned by fitness callback). Now I can clearly see when density-dependent scaling is not effective (subpopulation fitnessScaling >= 1,), individuals without beneficial mutation got a default absFitness >= 1.0 (1.0 * p1.fitnessScaling) and individuals with beneficial mutation is > 1.0 (>1.0 * p1.fitnessScaling) , all of which have of 100% survival probabilty. No offspring will die and there will be no selection.

However, when we increase fecundity, the density-dependent scaling will work – subpopulation fitnessScaling < 1. Individuals without the beneficial mutation will get a default relFitness 1.0 , but their absFintess will be scaled down by subpopulation fitnessScaling to less than 1.0 (have some probabilty to die). Individuals with the beneficial mutation will have a relFitness > 1.0 and, when scaled down by subpopulation fitnessScaling, it might get an absFitness > or < 1.0 but will always be > the absFitness of indiviual without the beneficial mutation. Thus, the "fitness" disadvantage of the individual without beneficial allele will make them die faster and less likely to survive to produce offspring in the next generation; the beneficial mutaton allle frequency will go up as individuals with the mutation are more likely survive and produce.

Thank you again for the help!!

early(){
	// kill older generations
	inds = p1.individuals;
	inds[inds.age > 0].fitnessScaling = 0.0;
	// adjust subpop fitness scaling factor to keep population size around N.
	p1.fitnessScaling = 1000 / length(inds[inds.age==0]);
	
	// default offspring absFitness
	defaultAbsFitness = 1.0 * 1.0 * p1.fitnessScaling;  / <===========================================
	catn('Generation: ' + sim.generation + 
		'\t id(default) AbsFitness: ' +  defaultAbsFitness +  
		'\t 100% Survival: ' +  (defaultAbsFitness>=1) );
}

fitness(m2){
	// simple negative frequncy-dependent selection
	relFit = 1.5 - sim.mutationFrequencies(p1, mut);
	
	updatedAbsFintess = relFit * individual.fitnessScaling * p1.fitnessScaling;  // <=====================
	if (individual.age == 0)
		catn('Generation: ' + sim.generation + 
		   '\t id(' + individual.index+') AbsFintess: ' +  updatedAbsFintess +
		   '\t 100% Survival: ' + (updatedAbsFintess>1));

	return relFit;
}

[bhaller]

Yes, it sounds like you've got it. Happy modeling!