Continuum of alleles: support as new StackPolicy, or simply a mention in the manual/extra?

[vsudbrack]

Hello,
following a chat to Ben, I'm writing here a bit about the continuum of alleles (CA) models, and opening the space for ideas of optimization in terms of implementation, the possibility of native support if it could be further improved under the hood (say a new stackPolicy 'a' for accumulate or 'ar' accumulate and replace).

---

The continuum of alleles (CA) model in population genetics describes a scenario in which a quantitative trait is influenced by a very large number of alleles, each contributing with an extremely small (ideally infinitesimal) effect to the quantitative trait. In CA models, this set of many alleles is assumed to be inherited together without recombination among them, meaning that their individual identities are less important than their cumulative effect. Because many alleles are thought to underlie this trait, the mutation rate in CA models is typically quite large (say 0.1 or 0.01), and this model is often summarized as "a trait evolving by many mutations of small effect".

This model is particularly useful when the interest is the distribution of a trait in a population (like its mean, variance, or possible polymorphisms) rather than the specific genetic architecture (i.e., which genes or loci are involved; although there are studies of genetic architecture that use CA models or modifications of these models). Furthermore, it allows for simplifying assumptions that make theoretical modeling and prediction more tractable, so people using SLiM to verify and expand mathematical results might also be very interested in efficient implementations of CA models in SLiM.

At this point, you're probably thinking that the CA model would be equivalent to our (default)StackPolicy as "stack", where within a linkage group there can be many non-recombining mutations whose overall effect is the sum of the effect of each mutation. And indeed it is. But I don't think that this implementation is the most efficient way to implement the CA model on SLiM, both in terms of time as well as memory.

Instead, because we are not interested in the individual identities of each mutation (and because mutation rate is high there will be many many many of these), I find that the most efficient way is to have a single "representative allele" at each position. In practice, this means that the StackPolicy is changed to "last", and whenever any new mutation appear at a position that already contains a representative allele, it's selectionCoeff accumulates the previous effect together with its own new effect. The mutation callback looks like this

mutation(m1) /* Modeling continuum of alleles: a mix between stack (accumulating selCoeff) and last (having a single mutation per position) policies */
{
	// Get position of new mutation
	mutpos = mut.position;
	
	// Gets which previous mutation was in that position
	allMuts = haplosome.mutationsOfType(m1);
	parentalMut = allMuts[allMuts.position == mutpos];
	
	// If there was no previous mutation, just accept
	if (size(parentalMut) == 0) return T;
		
	// If there was a previous mutation, accumulate selection coefficients
	prev = parentalMut.selectionCoeff;
	mut.setSelectionCoeff(prev + mut.selectionCoeff);
	
	// Always accept 
	return T;
}

Testing on my PC the code I'm attaching to the end of this message (a neutral simulation of a quantitative trait affected by 50 unlinked CA loci), I find a 7% reduction in running time (here measured between tick 1'000 and 5'000, with mutation rate 10^-3):

• Using last policy + accumulating selectionCoeff: 216.0s
• Using stack policy: 233.4s
  and I imagine as the number of mutations per position increase in the "stack policy" version of the code, the benefit of this implementation can also increases. I also imagine there is a saving in terms of memory because there are less mutation objects coexisting within the simulation.
  Also, I think converting mutations to substitutions here wouldn't save us from having many mutations segregating at a given time because as mutation rate is high and there is full linkage between mutations, things like compensation and redundancy among mutations would decrease the chances that a mutation is present across all individuals (I imagine it would be similar to the neutral probabilities of fixation in fact, but it probably depends on selective forces. Under disruptive selection, where multiple morphs can be maintained, mutations are in fact likely to never fix).

I'd be happy to hear any thoughts on how to improve it, and if you think this could be a native supported Stack Policy (especially if you think this could be further accelerated under the hood).

@bhaller , you mentioned this could be a recipe in the manual -- let's see if we gather more feedback on this Issue first and I'll compile some paragraphs similar to those written here describing the CA model and an example for you.

---

// set up a simple quantitative trait simulation
initialize() {
	setSeed(0);
	initializeMutationRate(1e-3);
	initializeMutationType("m1", 0.5, "n", 0.0, 0.01);
	initializeGenomicElementType("g1", m1, 1.0);
	m1.mutationStackPolicy = 'l'; //comment this line and the mutation callback for the stacking version
	m1.convertToSubstitution = F;
	initializeGenomicElement(g1, 0, 49);
	initializeRecombinationRate(0.5);
}

// create a population of 500 individuals
1 early() { sim.addSubpop("p1", 500); }

// calculates the phenotype at every generation (despite not using it, just to test)
1:5000 late() { p1.individuals.tagF = p1.individuals.sumOfMutationsOfType(m1); }

// turning off the QTL default fitness effect
mutationEffect(m1) { return 1.0; }


mutation(m1) /* Modeling continiuum of alleles: a mix between stack and last */
{
	// Get position of new mutation
	mutpos = mut.position;
	
	// Gets which previous mutation was in that position
	allMuts = haplosome.mutationsOfType(m1);
	parentalMut = allMuts[allMuts.position == mutpos];
	
	// If there was no previous mutation, just accept
	if (size(parentalMut) == 0) return T;
		
	// If there was a previous mutation, accumulate selection coefficients
	prev = parentalMut.selectionCoeff;
	mut.setSelectionCoeff(prev + mut.selectionCoeff);
	
	// Always accept 
	return T;
}

[bhaller]

I think that I like this idea, and that it might merit a new stacking policy. I'll circle back to this after the multitrait work is done. There will be some tricky things to ponder, however: what to do when the dominance coefficient is not compatible with combining the effects of the two mutations – i.e., the two mutations combined are not equivalent to the two mutations stacked – and what to do for multiplicative traits, since this seems more compatible with additive traits. I'll have to do a bit of head-scratching.

There is also a potential problem in that each time this occurred for A and B you'd get a new mutational lineage for (A+B) that would be distinct from other possible (A+B) lineages. This is a problem anyway, for stacked mutations; on a background of A you'd get new, distinct mutational lineages each time B arose. This works unintuitively in SLiM at present; if you get two mutations, B1 and B2, an individual with a B1/B2 genotype gets a heterozygous fitness effect for B1 and a heterozygous fitness effect for B2, rather than the homozygous fitness effect for B/B that might be expected. With complete additivity (h=0.5) this is all moot, but in other cases, and for multiplicative traits, it gets confusing. I'm not sure whether it's more of a problem with the proposed stacking policy than it would be otherwise. I'll need to think through whether this will work as it ought to, whether there are cases where an error ought to be raised rather than allowing the stacking policy to function, etc.

Any thoughts you have on this would be welcome, @vsudbrack!

[petrelharp]

In my view this is not a new stacking policy, but it is a suggestion for how to make sumOfMutationsOfType more efficient? The suggestion being to perhaps cache sums of effects of mutation above, or something. At least, the argument above was not "this is how it works biologically" but rather "this makes the code go faster".

Well, okay, thinking more, implementing it as a stacking policy ("cumulative") would be fine - better, probably. Just proposing another way of looking at it.

[bhaller]

In my view this is not a new stacking policy, but it is a suggestion for how to make sumOfMutationsOfType more efficient? ...

I agree that it isn't about biology, but about performance. And at least for the simple purely-additive-effects case – which is seems like this optimization is probably limited to – the multitrait work I'm doing now will achieve the same optimization anyway, in a different way, by caching the additive effects of entire mutation runs. So I suspect once the multitrait work is finished this optimization will no longer really be worthwhile. So right now I'm waiting to see how everything plays out.

[bhaller]

A note in favor of this proposal: in the continuum-of-alleles example that Vitor has worked out, mutation run effect caching would not achieve the same optimization as this proposal, because he is modeling 160 unlinked single-site loci, each of which effectively "stacks" quite intensively. Those single-site loci would not be represented with one mutation run each (in the present design, at least), so the mutation run effect caching would be entirely ineffective.

[bhaller]

OK, now that multitrait has pretty much settled down, I'm returning to this issue. I've updated the test model for SLiM 6, and made it run somewhat longer to get cleaner data:

// set up a simple quantitative trait simulation
initialize() {
	initializeTrait("trait1", "additive");
	initializeMutationRate(1e-3);
	initializeMutationType("m1", NAN, "n", 0.0, 0.01);
	initializeGenomicElementType("g1", m1, 1.0);
	m1.mutationStackPolicy = 'l'; //comment this line and the mutation callback for the stacking version
	initializeGenomicElement(g1, 0, 49);
	initializeRecombinationRate(0.5);
}

// create a population of 500 individuals
1 early() { sim.addSubpop("p1", 500); }

// calculates the phenotype at every generation (despite not using it, just to test)
1:50000 late() {
	sim.demandPhenotype(NULL, NULL);
}

mutation(m1) /* Modeling continuum of alleles: a mix between stack and last */
{
	// Get position of new mutation
	mutpos = mut.position;
	
	// Gets which previous mutation was in that position
	allMuts = haplosome.mutationsOfType(m1);
	parentalMut = allMuts[allMuts.position == mutpos];
	
	// If there was no previous mutation, just accept
	if (size(parentalMut) == 0) return T;
		
	// If there was a previous mutation, accumulate selection coefficients
	mut.trait1EffectSize = parentalMut.trait1EffectSize + mut.trait1EffectSize;
	
	// Always accept 
	return T;
}

I took out the setSeed() since the two versions of the model do not produce the same results (presumably due to roundoff). And I took out convertToSubstitution=F since that is no longer necessary for this type of model in SLiM 6, due to the new "baseline accumulation" feature. Averaged across three runs using different seeds:

With stacking: 75.3 seconds, 765 mutations segregating, 2377 substitutions

With accumulation: 62.9 seconds, 654 mutations segregating, 3 substitutions

So that's a 16.5% speedup. This is not affected by the length of the run, since substitutions now happen; the speed difference is not due to an ever-increasing buildup of fixed but un-substituted mutations. It is purely from the combining of the mutations themselves. There are fewer mutations segregating overall, and that presumably makes the phenotype calculation a bit faster as well as making offspring generation with recombination faster.

So, this seems quite worthwhile. The next question is whether it is worth building in as a stacking policy, or if it is fine just being a recipe in the manual to demonstrate the technique. So, I threw together a preliminary implementation; not ready for prime time yet, but it let me run timings (again averaged across three runs):

With stacking (to see if my changes slowed down the default case): 75.5 seconds, 764 mutations segregating, 2407 substitutions

With stacking policy 'a': 59.6 seconds, 658 mutations segregating, 2 substitutions.

So, the default case with stacking is unchanged, within statistical significance. The implementation of stacking policy a is ~3 seconds faster, or about 5% faster, than the implementation of accumulation with a mutation() callback as shown above. That seems worthwhile to me; it's faster, and the C implementation handles accumulation across multiple traits, etc., which the mutation() callback does not. The implementation was quite straightforward, so there's no downside there. So I'm going to polish this up and add it! Good suggestion @vsudbrack! I do think I'll make a recipe demonstrating it; I will try to adapt what you wrote above, Vitor, to discuss why this is useful and what it means. :->

[petrelharp]

Reconsidering here: I'm in favor, this would simplify some things greatly. But: what will you do about dominance?

[bhaller]

Reconsidering here: I'm in favor, this would simplify some things greatly. But: what will you do about dominance?

This mechanism only works for the independent-dominance case, so that will be required. But given that, it seems to work for both additive (and thus also logistic) and multiplicative traits, so that's nice. It appears to work for the hemizygous case only if the hemizygous dominance coefficient is either 0.0 or 1.0, values in between do not work; so that will also be required. But those requirements are still broad enough that it will be useful for quite a few models, I think.

[vsudbrack]

Hi Ben!

I'm happy to hear that your optimizations led to a fast CoA simulation. I think this can be really helpful for people doing simulations to extend theoretical predictions, and also for people wanting to predict polygenic trait dynamics.

Let me know if I can help with the description and/or a recipe to illustrate the new policy (the CoA example from QGS could work here to illustrate both the new policy as well as the calcVX functions).

[bhaller]

Hi Ben!

I'm happy to hear that your optimizations led to a fast CoA simulation. I think this can be really helpful for people doing simulations to extend theoretical predictions, and also for people wanting to predict polygenic trait dynamics.

Let me know if I can help with the description and/or a recipe to illustrate the new policy (the CoA example from QGS could work here to illustrate both the new policy as well as the calcVX functions).

Yes, a combo recipe that demonstrates both of those things would probably be great! :-> I'll get the features in first and then you can draft a recipe, maybe? Thanks!

[bhaller]

OK, I'm making progress on this in a bit of free time today. Here is my test model for the final feature:

initialize() {
	setSeed(10);
	initializeTrait("trait1", "additive");
	initializeMutationRate(1e-3);
	initializeMutationType("m1", NAN, "n", 0.0, 0.01);  // note independent dominance
	initializeGenomicElementType("g1", m1, 1.0);
	m1.mutationStackPolicy = 'a'; //comment this line for the stacking version
	initializeGenomicElement(g1, 0, 49);
	initializeRecombinationRate(0.5);
}
1 early() { sim.addSubpop("p1", 500); }
1:50000 late() {
	// calculate the phenotype at every generation (despite not using it, just to test)
	sim.demandPhenotype(NULL, NULL);
}

I'm back to using setSeed(), since it turns out to be useful. The results with stacking policy s versus a are very close; a trait value from one might be -0.784249 where the trait value from the other is -0.784248. The reason for that is now well-understood, and has to do with stacking policy a rounding down to single-precision float internally after each new mutation that comes in, whereas with stacking policy s each whole stack of mutations has its trait effect calculated in double-precision float without rounding, since all the mutations still exist. Anyhow, it's understood and it's fine. :-> It does mean that if test code checks for equivalence between the a and s approaches, that code should compare results with some error tolerance.

I've got it implemented for multiplicative traits too, now, and it does indeed work there as well; that has now been checked empirically and the results do match (with the same error tolerance issue). The runtimes for multiplicative traits are pretty much identical to those I posted above for additive traits; you get the same speedup. :->

I'll just fix the doc to cover this new stacking policy and then I'll commit and push.

[bhaller]

OK, this is in! I made a draft recipe 13.15 as part of this commit. I'll email it to you, @vsudbrack (and to you if you want to see it in this early draft form, @petrelharp). It'd be great if you could polish it up, correct any errors of terminology/math that I might have committed, and add whatever use of your calcX() functions and/or live plotting you feel like doing. That should allow it to become a more fully-fleshed demonstration of how to write a Continuum of Alleles model in SLiM, which would be great. :-> But I think it is probably already an ok recipe as it is, so no worries. I'll close this issue now; any additional work on the recipe can be done through email. Thanks for this very useful issue; I think this is a nice addition!