Matching recombination rates for recapitation of SLiM with msprime

[freshspaceoctopus]

Hello, I am currently trying to simulate neutral and positive selection on a complex demography with multiple chromosomes.

My understanding is that for msprime, the recombination rate is the mean probability of a split event using a poisson distribution, so to simulate breakpoints for multiple chromosomes, the recombination rate at the breakpoint should be given as log(2).

For SLiM, my understanding is that it is the exact binomial probability, so the value should be given as 0.5.

I want to give a recombination rate in msprime 1e-08 within a chromosome, and this is where I am confused.

I wonder for which simulation should I give the value 1e-08 and for which simulation should I modify the value to match 1e-08 so that the two simulations agree in terms of recombination. My guess is that 1e-08 is the expected number of recombinations, so the msprime values should be 1e-08 and SLiM should be 1 - exp(-1e-08), which is anyway almost 1e-08 where the probability and the rate values are not so different in the first place, so it isn't much of an issue, but I just wanted to be sure.

Thank you for your time.

[petrelharp]

Good question. I'll transfer this to a Discussion.

First, for between chromosomes, as you say, for msprime the value should be log(2), and for SLiM the value should be 1/2 (see bits of the SLiM manual below).

So, I think your question is for along chromosomes, when rates are very small (e.g., 1e-8). Both software uses a Poisson model of recombination; the difference is that SLiM takes the user-provided per-bp probability p and transforms it to a rate by λ = - log(1 - p). So, if msprime has a given value of λ, what value of p should we use? As you say, 1 - exp(-λ). And, there is virtually no difference at λ=1e-8, since |exp(x) - 1 - x| is O(x^2):

  > 1e-8 - expm1(1e-8)
[1] -5e-17

There are certainly bigger things to worry about in the simulation (e.g., discrete vs continuous time, bigger but still not worth worrying about).

Two caveats. First, note that since msprime's default model uses continuous time, we don't need to worry about "what if two breakpoints fall between two adjacent positions at once" like SLiM does. Second, I don't know why the dtwf model doesn't need to worry about this (i.e., why is log(2) still the right answer; see the sorts of considerations in the SLiM manual). Probably whoever wrote that bit of the docs was considering this (@apragsdale?).