SLiM's Wright-Fisher model and formation of next generation

[Hanbin Lee]

Dear SLiM community,

I have a question on the next generation is formed in SLiM's WF model.

Based on the experiments described below, I think what I assumed about SLiM's mechanism is wrong.

A common modeling assumption for the life-cycle of a population is

Parent generation (finite) -> Selection (finite) -> Gamete generation (infinite via sampling with replacement) -> Offspring generation (finite). The gamete pool in which the offspring generation is formed is considered infinite.

One such example is Waples (1989, see Figure 1, https://onlinelibrary.wiley.com/doi/abs/10.1111/j.1558-5646.1989.tb02571.x).

This serves as the basis of quantitative genetics modeling. We project the parent population's moments (e.g. mean, variance) to a corresponding Gaussian distribution. Then, we obtain the new population's mean and variance by sampling from the Gaussian distribution. For the variance side, this recursion yields a geometric brownian motion (when selection is absent). At stationarity, this variance follows an inverse-Gamma distribution with a polynomially decaying tail (much slower than normal's exponential tail).

In SLiM, I found that the population variance's distribution over replicate simulation is much narrower than what is predicted by the previous paragraphs. I simulated the following stabilizing selection scenario with very weak selection and high recombination to mimic an effectively neutral evolution. 

Legend: Each point is a simulation replicate. The distribution of V_t (population additive variance) over replicates is more spread in theory quantile compared to actual slim outcome.

The slim code:

```

// Simulate quantitative trait with infinite sites model
initialize() {
    if (!exists("N")) defineConstant("N", 1000);
    if (!exists("L")) defineConstant("L", 1000000);
    if (!exists("MUT_RATE")) defineConstant("MUT_RATE", 1e-8);
    if (!exists("REC_RATE")) defineConstant("REC_RATE", 1e-8);
    if (!exists("SIGMA_A")) defineConstant("SIGMA_A", 1.0);
    if (!exists("VS")) defineConstant("VS", 10000.0);
    if (!exists("MU")) defineConstant("MU", 10.0);
    if (!exists("GENERATIONS")) defineConstant("GENERATIONS", 10000);
    if (!exists("OUTPUT_FILE")) defineConstant("OUTPUT_FILE", "sim_output.csv");
    if (!exists("TRAIT_OUTPUT_FILE")) defineConstant("TRAIT_OUTPUT_FILE", "sim_traits.csv");
    if (!exists("TRAIT_SNAPSHOT_INTERVAL")) defineConstant("TRAIT_SNAPSHOT_INTERVAL", 100);

    initializeMutationType("m2", 0.5, "n", 0.0, SIGMA_A);
    initializeGenomicElementType("g1", m2, 1.0);
    initializeGenomicElement(g1, 0, L-1);
    initializeMutationRate(MUT_RATE);
    initializeRecombinationRate(REC_RATE);
}

1 early() {
    sim.addSubpop("p1", N);
   
    // Initial centering to 0, optimum is MU
    defineConstant("INITIAL_MEAN", 0.0);
    writeFile(OUTPUT_FILE, "tick,mean,var", append=F);
    writeFile(TRAIT_OUTPUT_FILE, "tick,individual_id,genetic_value", append=F);
   
    community.rescheduleScriptBlock(s1, start=1, end=GENERATIONS);
    community.rescheduleScriptBlock(s2, start=GENERATIONS, end=GENERATIONS);
}

mutationEffect(m2) {
    return 1.0;
}

s1 1: late() {
    genetic_values = p1.individuals.sumOfMutationsOfType(m2);
    mean_p = mean(genetic_values);
    var_p = var(genetic_values);
   
    log_line = paste(c(community.tick, mean_p, var_p), sep=",");
    writeFile(OUTPUT_FILE, log_line, append=T);

    if (community.tick % TRAIT_SNAPSHOT_INTERVAL == 0) {
        for (i in seqAlong(genetic_values)) {
            trait_line = paste(c(community.tick, i, genetic_values[i]), sep=",");
            writeFile(TRAIT_OUTPUT_FILE, trait_line, append=T);
        }
    }
}

fitnessEffect() {
    phenotype_raw = individual.sumOfMutationsOfType(m2);
    return exp( -1.0 * ((phenotype_raw - MU)^2) / (2.0 * VS) ); # stabilizing selection
}

s2 10000 late() {
    catn("Simulation finished at tick " + community.tick);
    community.simulationFinished();
}

```

Does SLiM follow a different offspring generation mechanism other than forming an infinite gamete pool? If so, can I adjust the option to follow this line of idea?

Best,

Hanbin

[Ben Haller]

Hi Hanbin!

The Wright–Fisher model is explicitly defined (by Wright and Fisher, not by me :->) as involving generating offspring independently, up to the specified population size, by choosing parents randomly from the parental population with probability proportional to fitness.  In other words, there is no separate selection step prior to offspring generation; selection is expressed in the relative probabilities that individuals will be chosen as parents, independently for each offspring.  This does lead to different results than the scheme you describe, I think.  With mortality due to selection prior to offspring generation, a given individual will either (a) not be a parent at all, because they are dead, or (b) will be a parent with equal probability to everybody else who survived the selection stage.  That is not true of the WF model; instead, the number of offspring produced by an individual will tend to be proportional to their fitness rather than being roughly bimodal (0 or x, where x is the same expectation for all individuals that survive selection).  I'm guessing this is the point of confusion you're hitting, but apologies if I've misunderstood.

The SLiM manual discusses this in various places; the reference doc for the WF model's mechanics is in chapter 24, but it is also covered in various "recipes" (e.g., section 15.12).

The solution is to use a nonWF model instead, which provides the flexibility to do whatever you want in terms of selection and reproduction.  Section 1.6 of the manual discusses differences between the WF and nonWF model types, chapter 15 provides an introduction to nonWF recipes, and chapter 25 is the reference doc for the nonWF model's mechanics.  If you're new to all this, though I would strongly suggest that you do the SLiM Workshop, which is available for free online with recorded lectures and worksheets.  Good luck, and happy modeling!

Cheers,
-B.

Benjamin C. Haller
Messer Lab
Cornell University

--
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 visit https://groups.google.com/d/msgid/slim-discuss/e7ffe2e0-6880-4523-8072-ff345f51d7b8n%40googlegroups.com.

[Gregor Gorjanc]

Hanbin,

The SLiM quantity you are calculating is additive genetic variance which includes genic component and LD component (due to linkage). Stabilising selection is inducing specific LD (I can’t recall if it is positive or negative). Is your analytic? SDE quantity considering LD? Your SDE quantity is often much higher than SLiM quantity indicating LD component would be negative, but you also sometimes underestimate…

gg

With regards,
  Gregor

[Hanbin Lee]

To Ben,

Thank you for pointing out the manual's chapter. That was exactly what I was looking for. 

To Gregor,

I now think this is indeed recombination related, though I need a bit more in-depth investigation. 

[Hanbin Lee]

Yes. It was indeed due to recombination, and more precisely due to finite-population LD.

Finite population LD is an extremely transient process (I won't put the maths here precisely but I plan to a preprint out pretty soon).

The genic variance evolves relatively slowly (the order of Brownian noise in std is ~1/N_e^2) while the finite same LD has noise magnitude of ~1/N_e. This drastic difference in noise level makes the variance distribution dominated by the massive finite-population LD fluctuation. 

Reading the manual carefully really helped organizing my thoughts. Thx!