Setting birth/death rate

[Daiki Tagami]

To whom it may concern,

It seems that the SLiM simulates the discrete generations based on the fitness values. I'm wondering, but would it be possible for us to set the birth and death rate to the model? Thank you very much for your time.

[Ben Haller]

Hi Daiki!

If you write a non-Wright-Fisher model in SLiM, you have complete control over births and deaths in the model, and can do whatever you wish.  I would suggest that you do the SLiM Workshop to learn more.  It is free, available online, and you can download that materials from http://benhaller.com/workshops/workshops.html.

Cheers,
-B.

Benjamin C. Haller
Messer Lab
Cornell University


Daiki Tagami wrote on 7/6/22 5:16 PM:

To whom it may concern,

It seems that the SLiM simulates the discrete generations based on the fitness values. I'm wondering, but would it be possible for us to set the birth and death rate to the model? Thank you very much for your time. --
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/1f1fd000-13ee-485e-9f92-aba2b24d1131n%40googlegroups.com.

[Daiki Tagami]

Dear Ben,

Thank you very much for your prompt response. It seems that the offspring are generated by using discrete generations, but would it be possible for us to generate an offspring through continuous time frame as well? I looked through the workshop materials and the manual as well, but I'm having some issues with implementing the birth/death rate inside the model.

Thank you,

Daiki

[Ben Haller]

Hi Daiki,

SLiM is a discrete-time framework; time proceeds in steps ("ticks").  However, generations can overlap and individuals can live for as many ticks as you wish, so you can approximate continuous time by having individuals live for a relatively large number of ticks on average.  If individuals have some probability of generating an offspring in each tick, and some probability of dying in each tick, then that leads directly to birth and death rates, I think.  Scale those probabilities down by a constant factor as much as you want to, to take the limit as the timestep approaches zero; that gives you continuous time.  But since SLiM is still a discrete-time framework, it will become increasingly inefficient as the time scale gets increasingly short.  This is not really a standard mode of operation for SLiM, although it should work fine.

Have you actually done the workshop, or just looked through the materials?  It is designed to teach you all the fundamentals, but you need to actually do it – learning is most effective when it is hands-on.  I'm happy to help, but at present your questions are so general that it is hard to know how to do so, beyond pointing you to the workshop.

Cheers,
-B.

Benjamin C. Haller
Messer Lab
Cornell University

Daiki Tagami wrote on 7/6/22 5:37 PM:

To view this discussion on the web visit https://groups.google.com/d/msgid/slim-discuss/8cf9e981-6f32-4743-9507-86c793d9d318n%40googlegroups.com.

[Daiki Tagami]

Dear Ben,

Thank you very much for your prompt reply. I did all the workshop materials and worked on some examples in the manual as well. I specifically want to implement the allee effects through setting

Birth rate = r(1-s)^k (A+1)

Death rate = r (n/N) (A+1)

where n is the number of individuals, N is the carrying capacity, s is the fitness effect, and k is the number of mutations (I only have 1 deleterious mutation).

Here, r and A are constants which I can set.

I tried imposing the carrying capacity by

subpop.fitnessScaling = N / subpop.individualCount;

but I got stuck afterwards. I tried setting the fitnessScaling through the birth/death rate that I have specified, but it didn't work out.

Do you have any suggestions? Thank you very much for your help.

Thank you,

Daiki

[Ben Haller]

Here is a simple model that implements a given birth and death rate:

initialize() {
    defineConstant("BIRTH_RATE", 0.001000);
    defineConstant("DEATH_RATE", 0.000990);
    initializeSLiMModelType("nonWF");
    initializeMutationType("m1", 0.5, "f", 0.0);
    m1.convertToSubstitution = T;
    initializeGenomicElementType("g1", m1, 1.0);
    initializeGenomicElement(g1, 0, 99999);
    initializeMutationRate(1e-7);
    initializeRecombinationRate(1e-8);
}
reproduction() {
    // assign birth randomly, at the birth rate
    if (runif(1) < BIRTH_RATE)
        subpop.addCrossed(individual, subpop.sampleIndividuals(1));
}
1 early() {
    sim.addSubpop("p1", 10000);
}
early() {
    // assign death randomly, at the death rate
    inds = p1.individuals;
    inds.fitnessScaling = 1 - DEATH_RATE;
}
20000 late() { }

I didn't understand your description of the equations you're trying to implement.  If you want further help you need to be much more specific.  (A) Provide your actual script: complete, minimal, and without any required command-line arguments (I debug in SLiMgui).  (B) Don't say "constants which I can set", specify the values you are using for all constants.  (C) Don't say "got stuck afterwards" or "it didn't work out", say what you expected to happen, what happened instead, what you tried to do to fix that problem, what that resulted in, etc.  If you got an error message, tell us the error message.  Completely specify the problem that you are trying to solve – we cannot guess what is in your head.

In any case, I imagine if you work on it you'll figure it out – please try to work on the problem for several days before asking more questions.  Learning to solve such problems is an essential part of coding; you cannot expect others to just give you a complete solution to the model you're trying to write.  But if, after that, you remain stuck, please feel free to ask again – I'm not trying to be mean, or discourage you from asking for help, I'm just trying to encourage you to do your best to solve your problem on your own, because that is a skill that must be developed if one is to become proficient.  :->

Happy modeling!

Cheers,
-B.

Benjamin C. Haller
Messer Lab
Cornell University

Daiki Tagami wrote on 7/6/22 5:58 PM:

To view this discussion on the web visit https://groups.google.com/d/msgid/slim-discuss/b428e3ef-f4df-4274-892b-4eb97f3e0c3an%40googlegroups.com.

[Daiki Tagami]

Dear Ben,

I'm very sorry about the late reply. I was surprised with how you have modeled the rate by using the uniform distribution. I looked up multiple papers, and I decided to use a Poisson distribution to create the model and it went well. Thank you very much for your help, and my issues are resolved now.

I just have a question, but would it be possible for SLiM to implement recombination among haploids? I ideally would like to create a haploid offspring which is a mixture of their parents. I would like to know if this is possible.

Thank you again for your time and help.

[Ben Haller]

Hi Daiki,

Let's talk about this uniform versus Poisson distribution issue, because I'm not sure what I think about it.  In the model I posted, I used a uniform distribution in this context:

reproduction() {
    // assign birth randomly, at the birth rate
    if (runif(1) < BIRTH_RATE)
        subpop.addCrossed(individual, subpop.sampleIndividuals(1));
}

I guess this is what you're referring to.  Here, my intention (perhaps incorrect, as I will discuss below) was to use BIRTH_RATE as the probability that a given individual will generate a new offspring; given that intention, using a uniform distribution is correct.  If BIRTH_RATE is 0.001, for example, as in the model I posted, then I wanted the focal individual for the reproduction() callback to generate an offspring 0.1% of the time.  But I think you are technically correct that since the model is intended to be an approximation of a continuous-time model, conceptually each reproduction event is independent and the distribution of the number of offspring for a focal individual, across one timestep, ought to be Poisson; for example, in a given timestep, it ought to be possible for an individual to give birth to more than one offspring (even if the probability of that, for a rate like 0.0001, is extremely small).

In my defense, though, I don't think you would want to scale the timesteps up to be large enough that this matters.  The whole point is to approximate continuous time by using very small timesteps, and if you scale the model up too far you are effectively out of continuous time and back in the world of discrete time.  As long as you keep the timescale of the model fine enough, my understanding (perhaps someone better versed than I in math can confirm?) is that using runif() is going to be essentially the same as using rpois(), I guess because of https://en.wikipedia.org/wiki/Poisson_distribution#Law_of_rare_events.  And drawing from a uniform distribution is substantially faster than drawing from a Poisson distribution, so I tend to take advantage of that approximation when I can.

Note that the way death is handled is also, in effect, making this approximation assumption:

early() {
    // assign death randomly, at the death rate
    inds = p1.individuals;
    inds.fitnessScaling = 1 - DEATH_RATE;
}

For each individual, SLiM will do a uniform draw from [0,1] and compare it to the fitnessScaling value, and if it is greater, the focal individual will die; that is, conceptually, the same as what I did for the birth code, using runif().  That is also technically wrong, for a continuous-time model, I suppose; and it will also malfunction if the timescale of the model is made too coarse.

I guess overall, my thinking is that since SLiM is not, in fact, a continuous-time simulator and makes various assumptions in its design that are based upon the discrete-time model, if you are going to try to approximate a continuous-time model in SLiM by making timesteps represent a very small time interval, you need to do that – keep your timesteps small.  In that regime, using runif() is a good approximation of the Poisson (I think?) and will run faster (perhaps negligibly faster, in this case, I don't know; I know that it makes a difference in some models).  Outside of that regime, probably a lot more will break than just that approximation assumption; so don't do that.  :->

But this is all debatable (assuming I am not simply wrong about the math), and if you want to use rpois() instead I think that is fine.  I just don't think it should matter; if it does matter, then you might be making your timesteps too large.  But please, if I am mistaken about all this, somebody correct me; I certainly don't want to be spreading misinformation!

-----------------

Regarding recombination among haploids, that has a simpler answer: see section 16.14 for an example.  :->

Cheers,
-B.

Benjamin C. Haller
Messer Lab
Cornell University

Daiki Tagami wrote on 7/19/22 9:07 AM:

To view this discussion on the web visit https://groups.google.com/d/msgid/slim-discuss/dd8a29ae-aba8-4654-afe9-11fd5996617fn%40googlegroups.com.

[Daiki Tagami]

Dear Ben,

Thank you very much for your detailed reply on uniform vs poisson distribution. It was very helpful for me, and I learned many new things from you. I really appreciate your help.

Also, thank you for pointing out to the haploid recombination section. I did not know that before, and I will read it in detail.

I really appreciate all your kind support, and it has been very helpful for me.

Thank you,

Daiki