initialize() {
initializeSLiMModelType("nonWF"); // Non-WF model
initializeSLiMOptions(dimensionality="xy");
initializeSex("A"); // Initialise sex
initializeMutationRate(1e-7);
initializeMutationType("m2", 0.5, "f", 0.0);
initializeGenomicElementType("g1", m2, 1.0);
initializeGenomicElement(g1, 0, 99999);
initializeRecombinationRate(1e-8);
// QTL parameters
defineConstant("QTL_sigma", 0.1);
defineConstant("SIGMA_K", 0.5);
// logging locality parameters
defineConstant("outx", 0.2);
defineConstant("outy", 0.2);
defineConstant("logbuffer", 0.1); // Buffer distance around locality of interest
// Other parameters
defineConstant("FECUN", 2.0);
defineConstant("N", 1000); // Initial population size
defineConstant("K_density", 20); // Number of individuals per gridcell
defineConstant("logfilename", "ind_n.txt");
// InteractionType used extract mean individual phenotype at outx, outy
initializeInteractionType(3, "xy", reciprocal=T, maxDistance=logbuffer);
}
1 first() {
sim.addSubpop("p1", N);
// initial positions are random in ([0,1], [0,1])
p1.individuals.x = runif(p1.individualCount);
p1.individuals.y = runif(p1.individualCount);
// Create map of random environmental variable
mapVar1 = c();
for (i in seq(1,12))
{
line = c();
for (j in seq(1,12)){
line = c(line, runif(1, (1/12)*(j-1), (1/12)*j));
}
mapVar1 = rbind(mapVar1, line);
}
// MAP CREATION and MAP VARIATION
defineConstant("var1", mapVar1);
map1 = p1.defineSpatialMap("map1", "xy", var1, interpolate=F,valueRange=c(0.0, 1.0), colors=c("red", "yellow"));
defineConstant("OPTIMUM", map1);
}
reproduction() {
// parents are chosen proportional to fitness
female=which(p1.individuals.sex=="F");
male=which(p1.individuals.sex=="M");
parents1 = sample(p1.individuals[female], N, replace=T);
parents2 = sample(p1.individuals[male], N, replace=T);
for (i in seqLen(N),
c in rpois(N, FECUN)){
p1.addCrossed(parents1[i], parents2[i], count = c);
}
self.active = 0; // call the reproduction callback only once
}
mutation(m2) {
// draw mutational effects for the new m2 mutation
effects = rnorm(1, 0, QTL_sigma);
mut.setValue("e0", effects);
return T;
}
early(){
inds = p1.individuals;
location = inds.spatialPosition;
optimum = OPTIMUM.mapValue(location);
fitness_norm = dnorm(0.0, 0.0, SIGMA_K);
for (ind in inds,
opt0 in optimum){
muts = ind.haplosomes.mutationsOfType(m2);
ind.setValue("phenotype0", size(muts) ? sum(muts.getValue("e0")) else 0.0);
ind.fitnessScaling = dnorm(ind.getValue("phenotype0"), opt0, SIGMA_K)/fitness_norm;
}
res = c(12, 12);
bounds = p1.spatialBounds;
density = summarizeIndividuals(inds, res, bounds, operation="individuals.size();", empty=0.0, perUnitArea=F);
d_map = p1.defineSpatialMap("density", "xy", density, F);
density_inds = d_map.mapValue(location);
// Piecewise inverse function
density_inds[density_inds < K_density] = K_density;
inds.fitnessScaling = inds.fitnessScaling * (1/(density_inds - (K_density - 1))); //K_density - 1 to get maximum fitness at K_density value
}
modifyChild() {
// draw a child position near the first parent, within bounds
do child.x = parent1.x + rnorm(1, 0, 0.02);
while ((child.x < 0.0) | (child.x > 1.0));
do child.y = parent1.y + rnorm(1, 0, 0.02);
while ((child.y < 0.0) | (child.y > 1.0));
return T;
}
late() {
// Creating log file
i3.evaluate(p1);
indall = p1.individualCount;
indtotal = i3.neighborCountOfPoint(c(outx, outy), p1);
if (indtotal == 0){
indnear = "NA";
mean_pheno = "NA";
mean_env = "NA";
} else {
// Get all individuals within a distance from point
indnear = i3.nearestNeighborsOfPoint(c(outx, outy), p1, indtotal);
// Get mean phenotype
mean_pheno = mean(indnear.getValue("phenotype0"));
// Get mean environmental variable
mean_env = mean(OPTIMUM.mapValue(indnear.spatialPosition));
}
line = paste(community.tick, indall, indtotal, mean_pheno, mean_env);
writeFile(logfilename, line, append=T);
}
1000 late() { sim.outputFixedMutations(); }