Resource scaling

[Hanbin Lee]

Dear slim community

I'm simulating stabilizing selection with biallelic loci. I mimic a burn-in with a Beta distribution and uses the 2nd recipe from the slim manual. 

```

// Stabilizing Selection Simulation (Recipe 2 + Quantitative Trait)

initialize() {
    if (!exists("N")) defineConstant("N", 100);
    if (!exists("L")) defineConstant("L", 1000);
    if (!exists("MU")) defineConstant("MU", 2e-5);
   
    // Stabilizing Selection Parameters
    if (!exists("VS")) defineConstant("VS", 10000.0);
    if (!exists("VE")) defineConstant("VE", 0.0);
    if (!exists("OPTIMUM")) defineConstant("OPTIMUM", 0.0);

    // Duration and Output
    if (!exists("DURATION")) defineConstant("DURATION", 10);
    if (!exists("OUTPUT_FILE")) defineConstant("OUTPUT_FILE", "phenotypes.csv");
    if (!exists("INIT_FILE")) defineConstant("INIT_FILE", "init_cond.csv");
   
    // Beta parameters for burn-in
    defineConstant("THETA", 4.0 * N * MU);
    defineConstant("ALPHA", THETA);
    defineConstant("BETA", THETA);
   
    initializeTreeSeq();
   
    // m2: Derived Allele
    initializeMutationType("m2", 0.5, "n", 0.0, 1.0);
    m2.convertToSubstitution = F;
    m2.color = "red";
   
    // m3: Back mutation marker
    initializeMutationType("m3", 0.5, "f", 0.0);
    m3.convertToSubstitution = F;
    m3.color = "cornflowerblue";
   
    initializeGenomicElementType("g1", m3, 1.0);
    initializeGenomicElement(g1, 0, L-1);
    initializeMutationRate(MU);
    initializeRecombinationRate(0.5);
}

1 early() {
    sim.addSubpop("p1", N);
   
    // Generate Canonical Mutations
    target = p1.haplosomes[0];
    target.addNewDrawnMutation(m2, 0:(L-1));
    defineConstant("MUT", target.mutations);
    target.removeMutations();
   
    // Burn-in
    gen = p1.haplosomes;
    for (i in 0:(L-1)) {
        p = rbeta(1, ALPHA, BETA);
        k = rbinom(1, 2*N, p);
        if (k > 0) {
            target_genomes = sample(gen, k);
            target_genomes.addMutations(MUT[i]);
        }
    }
   
    sim.treeSeqOutput("burnin.trees");
    catn("Tick 1: Burn-in complete.");
   
    writeFile(OUTPUT_FILE, "tick,mean,var,v_genic", append=F);
   
    community.rescheduleScriptBlock(s1, start=2, end=DURATION);
    community.rescheduleScriptBlock(s2, start=DURATION, end=DURATION);
}

mutation(m3) {
    if (haplosome.containsMarkerMutation(m2, mut.position)) {
        return T;
    }
    return MUT[mut.position];
}

mutationEffect(m2) {
    return 1.0;
}

// Phenotype Fitness with Centering
fitnessEffect() {
    phenotype_raw = individual.sumOfMutationsOfType(m2);
    phenotype_centered = phenotype_raw - INITIAL_MEAN;
   
    env_noise = (VE > 0) ? rnorm(1, 0, sqrt(VE)) else 0.0;
    phenotype_final = phenotype_centered + env_noise;

    return exp( -1.0 * (phenotype_final^2) / (2.0 * VS) );
}

1 late() {
    // Initial Mean Centering
    inds = p1.individuals;
    phenotypes = inds.sumOfMutationsOfType(m2);
    mean_p = mean(phenotypes);
   
    defineConstant("INITIAL_MEAN", mean_p);
   
    catn("Initial Mean Phenotype: " + mean_p);
    catn("Centering logic initialized.");
   
    // Export Initial Conditions (Effects and Frequencies)
    output_lines = c("locus,effect,freq");
   
    all_freqs = sim.mutationFrequencies(p1, MUT);
    effects = MUT.selectionCoeff;
   
    for (i in 0:(size(MUT)-1)) {
        line = i + "," + effects[i] + "," + all_freqs[i];
        output_lines = c(output_lines, line);
    }
   
    writeFile(INIT_FILE, output_lines);
}

s1 2: late() {
    // Vectorized Cleanup
    m3_muts = sim.mutationsOfType(m3);
   
    if (size(m3_muts) > 0) {
        all_genomes = p1.haplosomes;
        counts = all_genomes.countOfMutationsOfType(m3);
        indices = which(counts > 0);
       
        if (size(indices) > 0) {
            dirty_genomes = all_genomes[indices];
            for (g in dirty_genomes) {
                local_m3 = g.mutationsOfType(m3);
                local_m2 = g.mutationsOfType(m2);
                pos_m3 = local_m3.position;
                pos_m2 = local_m2.position;
                s = (match(pos_m2, pos_m3) >= 0);
                stacked_m2 = local_m2[s];
                g.removeMutations(c(local_m3, stacked_m2));
            }
        }
    }
   
    // Logging
    if (community.tick % 1 == 0) {
        genetic_values = p1.individuals.sumOfMutationsOfType(m2) - INITIAL_MEAN;
        mean_p = mean(genetic_values);
        var_p = var(genetic_values);
       
        // Calculate Genic Variance (Vg)
        muts = sim.mutationsOfType(m2);
        if (size(muts) > 0) {
            freqs = sim.mutationFrequencies(p1, muts);
            effects = muts.selectionCoeff;
            v_genic = sum(2.0 * freqs * (1.0 - freqs) * (effects^2));
        } else {
            v_genic = 0.0;
        }
       
        log_line = paste(c(community.tick, mean_p, var_p, v_genic), sep=",");
        writeFile(OUTPUT_FILE, log_line, append=T);
       
        if (community.tick % 10 == 0) {
           catn("Tick " + community.tick + ": Var=" + var_p + " Vg=" + v_genic);
        }
    }
}

s2 2000 late() {
    sim.treeSeqOutput("stabilizing.trees");
    catn("Tick " + community.tick + ": Finished.");
    community.simulationFinished();
}

```

I noticed that with biallelic loci, the program demands substantially more memory than without it. With the following setting

config.yaml

```

N: 5e3

L: 1000

MU: 2e-6

VS: 1000.0

VE: 0.0

```

It consumes ~20GB of memory. When I increase L (number of loci) to 2000, it exceeds 80GB. I have mainly two questions to this

(1) Why is forcing biallelic so memory heavy?

(2) Does slim generally scale quadratically as a number of loci?

[Ben Haller]

Hi Hanbin,

Without getting too deep into the weeds, my first guess would be that with biallelic loci the underlying "mutation run" structure of the model is no longer leveraging shared IBD runs effectively.  You can read about mutation runs in section 22.4.  Profiling, whether in SLiMgui (section 22.5) or at the command like (section 22.7), should provide a breakdown of memory usage that might (or might not) be illuminating as to where the memory usage is going.  As that discussion of mutation runs perhaps illuminates, no, memory usage in SLiM is much more complex than just a scaling – whether linear or quadratic – with the number of loci, and depends upon factors such as the mutation density, the amount of IBD sharing, etc.  Feel free to post if you have further questions; if you do, perhaps providing the profiles for the two models in question would be useful.  I hope this helps; 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/3af4d895-d986-42b9-bf89-142c724d8913n%40googlegroups.com.

[Hanbin Lee]

I read through and it makes a lot of sense. I think I should think harder about what I'm really trying to do.