MesserLab / SLiM

SLiM is a genetically explicit forward simulation software package for population genetics and evolutionary biology. It is highly flexible, with a built-in scripting language, and has a cross-platform graphical modeling environment called SLiMgui.
https://messerlab.org/slim/
GNU General Public License v3.0
160 stars 33 forks source link

setConstraints() clashing with sexSegregation? #407

Closed jiseonmin closed 12 months ago

jiseonmin commented 12 months ago

Hi Ben, this is a follow-up to our conversation here: https://github.com/kr-colab/spatial_sims/pull/129 It seems like setConstraints() runs into an issue when it's used for mate choice, and there is "FM" distinction. Here is an example:

    initializeSLiMModelType("nonWF");
    initializeSLiMOptions(dimensionality="xy");
    initializeSex("A");
    // mate choice
    initializeInteractionType(2, "xy", reciprocal=T, maxDistance=SM * 3, sexSegregation="FM");
    i2.setInteractionFunction("n", 1.0/sqrt(2*PI*SM^2), SM);
    i2.setConstraints("both", minAge = 2);

...

reproduction(NULL, "F") {
    mate = i2.drawByStrength(individual, 1);
    can(mate.size());
    if (mate.size()) {
        nOff = rpois(1, FECUN);
        for (i in seqLen(nOff)) {
            pos = individual.spatialPosition + rnorm(2, 0, SD);
            if (p1.pointInBounds(pos)) {
                offspring = subpop.addCrossed(individual, mate);
                offspring.setSpatialPosition(pos);
            }
        }
    }
    return;
}

Running the script, it prints 0 and population goes extinct. (i.e. no one finds a mate). If I do the same for hermaphrodites, there is no issue. If I set constraints for a different interaction (e.g. i1, for competition), there is no issue. I hope this narrows down what to debug.

bhaller commented 12 months ago

Hi @jiseonmin! OK, so first of all, when you call setConstraints() the sexSegregation setting gets blown away. In the initializeInteractionType() doc we have this:

The sexSegregation parameter is shorthand for setting sex constraints on the interaction type using the setConstraints() method; see that method for a more extensive set of constraints that may be used.

and in the doc for setConstraints():

The constraints specified by a given call to setConstraints() override all previously set constraints for the category specified (receivers, exerter, or both).

So I'd recommend that you leave out the sexSegregation parameter (since it will get blown away anyhow), and do separate calls to setConstraints() for "exerter" and "receiver", setting the desired constraints on each, including the sex constraint for each. I take it that wasn't clear from the doc as written; if you have a recommendation for what would make it clearer, that would be great.

Second, I can't see the full model in the code you posted above – what is the initial age distribution? If you don't set any ages on any of the individuals, then by default they are age 0 – and so there will be no eligible mates. Is it possible that this is what's happening?

Post the full code (minimal reproducible example, if possible :->), so I can investigate further, if that is needed. Thanks!

jiseonmin commented 12 months ago

Oh sorry, I must have not read the doc carefully. That makes sense! For your reference, here is a full code. I removed the age constraint, but they still can't find mates. (printing 0s)

// Disperse only from a female, effective fecundity is halved due to FM distinction.
initialize() {

    initializeSLiMModelType("nonWF");
    initializeSLiMOptions(dimensionality="xy");
    initializeSex("A");

    if (!exists("SD")) defineConstant("SD", 0.3);  // sigma_D, the dispersal distance
    if (!exists("SI")) defineConstant("SI", SD);  // sigma_I, the spatial interaction distance
    if (!exists("SM")) defineConstant("SM", SI); // sigma_M, the mate choice distance
    if (!exists("K")) defineConstant("K", 5);  // carrying-capacity per unit square (roughly) 
    if (!exists("WIDTH")) defineConstant("WIDTH", 25.0);  // width of the simulated area
    if (!exists("HEIGHT")) defineConstant("HEIGHT", WIDTH); // height of the simulated area
    if (!exists("LIFETIME")) defineConstant("LIFETIME", 4);    // mean lifetime at stationarity ( = 1 / mean fecundity)
    if (!exists("L")) defineConstant("L", 1e8);  // genome length
    if (!exists("R")) defineConstant("R", 1e-8); // recombination rate
    if (!exists("MU")) defineConstant("MU", 1e-8); // mutation rate
    if (!exists("OUTPATH")) defineConstant("OUTPATH", "./data/");
    defineConstant("FECUN", 1/LIFETIME); // mean fecundity
    defineConstant("RHO", FECUN/2/((1+FECUN/2) * K)); // constant in spatial competition function
    initializeMutationType("m1", 0.5, "f", 0.0);
    initializeGenomicElementType("g1", m1, 1.0);
    initializeGenomicElement(g1, 0, L-1);
    initializeMutationRate(MU); 
    initializeRecombinationRate(R);

    // spatial competition
    initializeInteractionType(1, "xy", reciprocal=T, maxDistance=SI * 3);
    i1.setInteractionFunction("n", 1.0/sqrt(2*PI*SI^2), SI);
//  i1.setConstraints("both", minAge = 2);
    // mate choice
    initializeInteractionType(2, "xy", reciprocal=T, maxDistance=SM * 3);
    i2.setInteractionFunction("n", 1.0/sqrt(2*PI*SM^2), SM);
    i2.setConstraints("receiver", sex = "M");
    i2.setConstraints("exerter", sex = "F");

}

reproduction(NULL, "F") {
    mate = i2.drawByStrength(individual, 1);
    catn(mate.size());
    if (mate.size()) {
        nOff = rpois(1, FECUN);
        for (i in seqLen(nOff)) {
            pos = individual.spatialPosition + rnorm(2, 0, SD);
            if (p1.pointInBounds(pos)) {
                offspring = subpop.addCrossed(individual, mate);
                offspring.setSpatialPosition(pos);
            }
        }
    }
    return;
}

1 first() { 
    sim.addSubpop("p1", asInteger(K * WIDTH * HEIGHT));
    p1.setSpatialBounds(c(0, 0, WIDTH, HEIGHT));
    p1.individuals.setSpatialPosition(p1.pointUniform(p1.individualCount));
}

first() {
    // choose mates
    i2.evaluate(p1);
}

early() {
    i1.evaluate(p1);
    inds = p1.individuals;
    competition = i1.localPopulationDensity(inds);
    inds.fitnessScaling = 1/(1 + RHO * competition);
}

1000 late() {
    sim.simulationFinished();
}
bhaller commented 12 months ago

Looks like you have receivers constrained to males, and exerters constrained to females:

i2.setConstraints("receiver", sex = "M");
i2.setConstraints("exerter", sex = "F");

I think you want the opposite, since you call drawByStrength() to find males whose exerted interaction effect is received by the focal female. Interactions are always in this direction, from exerters to receivers; so the signature for drawByStrength() is:

– (object)drawByStrength(object<Individual> receiver, [integer$ count = 1], [No<Subpopulation>$ exerterSubpop = NULL], [logical$ returnDict = F])

The individual passed in is receiver. I reversed the constraints and the model then did find mates, although the population size gradually declines. Again, if there's a documentation angle here that could be improved, let me know. Thanks. :->

jiseonmin commented 12 months ago

Ohhhh okay that worked! Thanks, I will close this issue.