Balancing selection with preferred frequency in nucleotide models

[Svitlana Braichenko]

Hi team,

I'm quite new to Slim and I work with version 3.7.1. Will appreciate any help.

I'm trying to study the balancing selection example in section 10.4 of the Slim manual and adapt it to the nucleotide-based models.

The adapted code is
initialize() {
  initializeSLiMOptions(nucleotideBased=T);
  defineConstant("L", 1e4);
  initializeAncestralNucleotides(randomNucleotides(L));
 
  // m1 mutation type: neutral
  initializeMutationTypeNuc("m1", 0.5, "f", 0.0);
  initializeMutationTypeNuc("m2", 0.5, "f", 0.1); // balanced
 
  // g1 genomic element type: uses m1 for all mutations
  initializeGenomicElementType("g1", c(m1,m2), c(999,1), mmKimura(1e-7, 2e-7));
 
  // uniform chromosome of length 10 kb with uniform recombination
  initializeGenomicElement(g1, 0, L-1);
  initializeRecombinationRate(1e-8);
}

// create a population of 500 individuals
1 {
  sim.addSubpop("p1", 100);
}
500 {
  sim.addSubpopSplit("p2", 50, p1);
}
1000 {
  sim.addSubpopSplit("p3", 25, p2);
}
100000 late() {
  // print mutant freqs
  cat("Frequencies m1:\n");
  print(sim.mutationFrequencies(p1, sim.mutationsOfType(m1)));
  cat("Frequencies m2:\n");
  print(sim.mutationFrequencies(p1, sim.mutationsOfType(m2)));
  }
Most of the output frequencies are around ~0.5 or ~1.0, indicating the possible preferred frequencies as far as I understand. So my question is if there is a way in Slim to study the balancing selection with the pre-set value of the preferred frequency in the population, for example, 0.1? Also, is there a way to set the preferred frequency of certain nucleotides, e.g. A?

Many thanks
Svitlana  

[Ben Haller]

Hi Svitlana!  OK, so when you say "section 10.4", I take it you mean section 10.4.1 ("Frequency-dependent selection").  I'm a bit puzzled by the adapted code that you post, because it does not, in fact, seem to implement balancing selection – it is missing the mutationEffect() callback that the original recipe in section 10.4 has:

mutationEffect(m2) {
   return 1.5 - sim.mutationFrequencies(p1, mut);
}

Without that callback, mutations of type m2 are just beneficial mutations with a selection coefficient of 0.1, as assigned to them by your initializeMutationTypeNuc() call.  (Note that I'm talking about a mutationEffect() callback because I am on SLiM 4.0; it looks like maybe you are on SLiM 3.7.1, so it would be called a fitness() callback, but it is the same thing.)

The frequencies reported by this model, in the test run I just did, are typically either 0.0 or 1.0.  This has nothing to do with balancing selection; it is because there is no migration at all between p1, p2, and p3, so mutations are either absent from p1 altogether (frequency 0.0) or fixed in p1 (frequency 1.0).  If you run the model under SLiMgui (which I always recommend, so that you can see model problems visually), you will see that mutations never fix (they can't fix across the whole population, since there is no migration); they just drift until they are either lost, or get locally fixed within the subpopulation they are in.

As for changing the "preferred frequency" of the mutation under balancing selection – the frequency at which it has the highest fitness – this is up to the code in the mutationEffect() callback.  It can define the fitness of the mutation to be whatever function of the mutation frequency is desired.  The version of that callback shown in section 10.4.1 makes 0.5 be the optimal frequency, because 1.5 - 0.5 is 1.0; at frequencies higher than 0.5 the callback will make the mutation deleterious, at frequencies lower than 0.5 it will make the mutation beneficial.  You can substitute whatever math you want there.  For example, if the script were:

mutationEffect(m2) {
   return 1.6 - sim.mutationFrequencies(p1, mut);
}

then the optimal frequency would then be 0.6, because then 1.6 - 0.6 = 1.0, so frequency 0.6 is the point at which the frequency-dependent selection crosses from deleterious to beneficial.  But you might also want to change the shape of this fitness function in other ways, to change the strength of the negative frequency-dependence, or make it non-linear, or whatever.  This is why scriptability is good; once you understand how the code works, you can modify it to do whatever you want it to do.

If all of this is less than clear, then I wonder whether you might want to do the SLiM Workshop.  It is available for free online at http://benhaller.com/workshops/workshops.html, and will introduce you to concepts like fitness() callbacks, working with Eidos scripting, migration between subpopulations, using SLiMgui, etc.  It is the best way to get started with SLiM.  If you do intend to do the workshop, then you should probably stay on SLiM 3.7.1 for now, since the workshop materials are not yet updated for SLiM 4.

Good luck, and happy modeling!

Cheers,
-B.

Benjamin C. Haller
Messer Lab
Cornell University

[Svitlana Braichenko]

Hi Benjamin

Thanks for your fast and detailed reply. Yes, you are right I mean section 10.4.1 and sorry I forgot to add the fitness() callback, thanks for the clarification of why it's necessary. I did look through the workshop materials and in the manual, apology if I missed something. Your suggestion to add

fitness(m2) {

  return 1.5 - sim.mutationFrequencies(p1, mut);
 }

would drive the frequencies in population p1 to the balancing selection at the preferred frequency of 0.5, and now I understand how to change this frequency. I'm curious, is there a way to introduce the same balancing selection in subpopulations p2 and p3 in the above code?

Many thanks,
Svitlana

[Ben Haller]

Hi Svitlana!

Yes, that's straightforward.  See section 25.2 (well, in the SLiM 4 manual, anyway) for discussion of mutationEffect() callbacks in more detail; I think it is probably the same section for fitness() callbacks in the SLiM 3.7.1 manual.  With that information under your belt, there are two typical strategies you could take:

- Declare a single fitness() callback that, instead of using `p1` as in the code below, uses the `subpop` pseudo-parameter, thus evaluating the frequency of the focal mutation within the subpopulation the individual is in:

fitness(m2) {
  return 1.5 - sim.mutationFrequencies(subpop, mut);
}

- Declare a fitness() callback for each subpopulation separately, such as:

fitness(m2, p1) {

  return 1.5 - sim.mutationFrequencies(p1, mut);
}

and similarly for p2, p3, etc.  The main difference between these approaches is that the first is simpler (just one callback handles all subpopulations), but the second allows you to vary the optimum frequency from one subpopulation to another since each subpopulation is handled by a separate callback that can have different script.  Of course, given scriptability, there are any number of other approaches you could take instead, such as giving each subpopulation a different `tag` value and then making the callback's formula depend upon the value of that `tag` in some way, etc.  Once you really wrap your head around Eidos scriptability, it becomes clear that all the recipes in the manual are just incredibly rudimentary starting points; the SLiM Workshop will help with getting to that point.  :->

I hope this clarifies things; happy modeling!

Cheers,
-B.

Benjamin C. Haller
Messer Lab
Cornell University

Svitlana Braichenko wrote on 8/23/22 4:18 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/ac00d54c-bfdb-47cc-89a1-f7557a210301n%40googlegroups.com.

[Svitlana Braichenko]

Hi Benjamin

Thanks, all seems to work!

Best

Svitlana

[Svitlana Braichenko]

Hi all

I’m posting here because it’s still relevant to this issue as I'm using nucleotide-based models. I’m still running SLiM 3.7.1 and now I would like to simulate GC-biased gene conversion (section 18.11), heterozygote advantage through overdominance and both of them combined together (code attached, please note that the rates in the script are not physiological to make simulations run faster). I run a similar setup with 4 subspecies evolving from one ancestor.

I can successfully simulate GC-bias and overdominance separately as shown in the site frequency spectra (SFS) plot attached. It is shown by the higher frequencies of G and C compared to A and T in the first case and by the polymorphic peak in the middle of every spectrum in the second case. However, when I try to combine two effects together, I get SFS very similar to the overdominance only. So it seems like the overdominance is overpowering the GC-bias effect. Could you please tell me if it’s possible to combine these two effects in SLiM or am I missing something?

Many thanks in advance,
Svitlana

[Ben Haller]

Hi Svitlana!

Sorry for the slow reply – you got me at a very complicated time, with a SLiM workshop and a conference and a vacation and a whole lot of driving.  :->  I was hoping that somebody else would have a good answer to this, and then it got buried in my inbox.  :->

I do not have a particularly good answer to this.  I think there's no reason why GC-biased gene conversion and overdominance wouldn't work together in the same model, from SLiM's perspective; that ought to be fine.  Perhaps what you are seeing is simply the expected result, from a popgen perspective, given the parameters you are using?  I guess I'd suggest working through this with the usual sort of investigatory methods:

- first, simplify the model down to a minimal model to make it easier to work with (get rid of all the demographic changes etc.)

- then, try making the overdominance extremely weak (not zero, but so close to zero that it shouldn't matter) and see if the GC bias is then detectable

- try other intermediate strengths for the overdominance to see whether the GC bias becomes harder and harder to see in the results

- think really hard about why the two phenomena might interact in the way you're observing; formulate hypotheses, and modify the model in different ways to test those hypotheses

I don't have any advice beyond that.  My guess is that what you're seeing is a popgen phenomenon, not any issue with SLiM, but if you gather evidence indicating otherwise, of course do let me know.  :->  Happy modeling!

Cheers,
-B.

Benjamin C. Haller
Messer Lab
Cornell University

[Svitlana Braichenko]

Hi Ben,

Thanks a lot for your response. Unfortunately, I couldn't solve the previous problem. It's okay though because my focus has shifted, and it's not as important now.

I'm currently working on simulating overdominance only (without GC bias) and adjusting mutation rates and recombination rates to be closer to physiological values (as indicated in the attached script). However, I'm struggling to achieve a strong enough overdominance effect. This means I'm not getting enough polymorphic sites in the simulations to produce a detectable effect. I tried increasing the dominance coefficient and selection coefficient in the balanced mutation through "initializeMutationTypeNuc("m2", 1.5, "f", 0.5)," but this resulted in the error "Simulation Runtime Error ERROR (Subpopulation::UpdateFitness): total fitness of subpopulation is not finite; numerical error will prevent accurate simulation."

Could you please advise if it's possible to overcome this error?

Best wishes, 

Svitlana

[Ben Haller]

Hi Svitlana,

No – if your simulation is using such extreme fitness values that the fitness of individuals ends up overflowing to infinity, there's no workaround for that.  :->  This is the sign of a bug in your script, or a problem with your modeling approach.  With a dominance coefficient of 1.5 and a selection coefficient of 1.5, 1+s is a fitness of 1.5 for homozygotes, and 1+hs is a fitness of 1.75 for heterozygotes.  Are those the values you want?  They are very large fitness effects, so for a reasonably biologically realistic simulation I'd expect that you'd have a relatively small number of such mutations segregating in the population.  For an individual with 1200 such heterozygous mutations, its fitness is about 4.4x10^291 times higher than the fitness of an individual with no mutations.  That's a VERY large fitness difference; it's not clear that that is biologically meaningful at all.  A bit north of 1200 such mutations, the fitness overflows to infinity due to the nature of floating-point calculations on computers; there is a limit to have large of a value can be represented as an IEEE double-precision float (which is what SLiM uses).

I would suggest that you probably need to think more about the biology of what you're trying to simulate and what it looks like, and move your model towards more realistic parameters so that it doesn't overflow to infinity in this manner.  Alternatively, if you really want to be able to represent arbitrarily large fitness values, you may need to shift to an analytical modeling approach.  Good luck, and happy modeling!

Cheers,
-B.

Benjamin C. Haller
Messer Lab
Cornell University

Svitlana Braichenko wrote on 1/31/24 2:24 PM:

To view this discussion on the web visit https://groups.google.com/d/msgid/slim-discuss/18721f8a-4933-4bed-bcaf-a3abe0d85b48n%40googlegroups.com.

[Svitlana Braichenko]

Hi Ben

Many thanks for your promt reply. Yes, you're right the fitness I used is indeed too high. Your reply made me thinking that instead of increasing heterozygosity and selection coefficient I probably need to wait longer for more mutations to accumaulate.

Many thanks 

Svitlana

[Ben Haller]

Great!  :->

A note for posterity: I should have written "With a dominance coefficient of 1.5 and a selection coefficient of 0.5, ...".

Cheers,
-B.

Benjamin C. Haller
Messer Lab
Cornell University

Svitlana Braichenko wrote on 2/2/24 7:10 AM:

To view this discussion on the web visit https://groups.google.com/d/msgid/slim-discuss/a454f89c-67d6-41be-8370-235d0417aca4n%40googlegroups.com.