3.13 Code from Week 3

Here is the final simulation developed in class:
// set up a simple neutral simulation
initialize() {
    initializeMutationRate(1e-7);
    
    // m1 mutation type: neutral
    initializeMutationType("m1", 0.5, "f", 0.0);
    
    // g1 genomic element type: uses m1 for all mutations
    initializeGenomicElementType("g1", m1, 1.0);
    
    // uniform chromosome of length 100 kb with uniform recombination
    initializeGenomicElement(g1, 0, 99999);
    initializeRecombinationRate(1e-8);
}

// create a population of 500 individuals
1 early() {
    sim.addSubpop("p1", 500);
    sim.addSubpop("p2", 3000);
}

// allow migration from p1 to p2 
100 early() {
    p2.setMigrationRates(p1, 0.2);
}

// allow migration from p2 to p1 
200 early() {
    p1.setMigrationRates(p2, 0.5);
}

// stop allowing migration from p2 to p1
300 early() {
    p1.setMigrationRates(p2, 0);
}
// view these dynamics by clicking the bar chart button on the right and then "Graph Population Visualization"

// create a new population 
// allow to migrate into both p1 and p2 
// don't receive any migration
400 early() {
    sim.addSubpop("p3", 3000);
    p2.setMigrationRates(c(p1, p3), c(0.2,0.1)); //allow migration from p1 and p3 into p2
    p1.setMigrationRates(p3, 0.25);
}


// output samples of 10 genomes periodically, all fixed mutations at end
1000 late() { p1.outputSample(10); }
2000 late() { p1.outputSample(10); }
2000 late() { sim.outputFixedMutations(); }