Forcing neutral dynamics

[Manas G. A.]

Hi Ben,

I am attempting to do the following:

1. Generate neutral diversity and assign selection coefficients from a DFE using msprime and Pyslim

2. Read the tree-sequences generated by msprime

3. Run the population forward using SLiM (with no new mutations), such that the dynamics are forced to be neutral for the first five genrations, but are subsequently governmed by selection coefficients assigned using msprime.

I tried using the fittnessEffect() callback for the first five generations to force neutral dynamics as follows:

#########################################

1:5 fitnessEffect() { return exp(1.0); }

1 late() {

sim.readFromPopulationFile("neutral_burnin.trees"); // This is the.trees file generated by msprime

sim.treeSeqRememberIndividuals(p0.individuals);

}

1:51 late() {

sim.outputFull( "Output_" + sim.cycle + ".txt");

}

##########################################

I am not sure if this is the correct way to implement what I need. Note that the simulation does run along expected lines if I chose not to force neutral dynamics in the first 5 generations (by deleting the fitnessEffect() callback).

Kindly help!

Cheers,

Manas

 

[Ben Haller]

Hi Manas!

There are a few issues here.

- First of all, I'm not sure why you're returning exp(1.0) from your callback.  Neutrality is a fitness effect of 1.0, so that is what you would want to use.  (I'm a bit curious where the idea to use exp(1.0) came from, actually!)

- Second, a fitnessEffect() callback is not quite the right tool.  That callback provides an *additional* fitness effect for individuals, which gets multiplied together with all of the other existing fitness effects – notably, the fitness effect from the non-neutral mutations.  It does not cancel out those other fitness effects.  So returning 1.0 from a fitnessEffect() callback has no effect on overall fitness values.  (Returning exp(1.0) multiplies the individual fitness values that would otherwise exist by 2.718..., which presumably also is not what you want.)

- Instead, you want to use a mutationEffect() callback (which used to be called fitness() callbacks, for those still living in the Before Times :->).  This callback type does override the effect of mutations, replacing their default effect with a different effect.  So you could define a mutationEffect(m1) callback (or whatever mutation type is relevant, if it is not m1), and return 1.0, and that should do what you want.

- Assigning the fitness effects on the Python side and then temporarily overriding them to be neutral in SLiM is fine.  It is a little bit inefficient, perhaps, since that mutationEffect() callback will take a little time to run, but that won't be noticeable unless you have a large population size and many such mutations.  Possibly a cleaner/faster design would be to assign the fitness effects in SLiM, at the end of tick 5, using setSelectionCoefficient().  But unless you see a performance problem with your current approach, that probably isn't worth worrying about.

I hope this clears things up.  Happy modeling!

Cheers,
-B.

Benjamin C. Haller
Messer Lab
Cornell University


Manas G. A. wrote on 8/7/23 6:57 AM:

--
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/30b3d13e-89c0-4510-8d03-803590a76c73n%40googlegroups.com.