local_adaptation.slim

// initialize genetic architecture of local adaptation
initialize() {
	if (exists("slimgui")) {
		defineConstant("m", 1e-5);   	// migration rate
		defineConstant("mu", 1e-5);  	// mutation rate
		defineConstant("sigsqr", 25);    // width of fitness function
		defineConstant("N", 1000);		// size of populatio
		defineConstant("numpos", 160); 		// num positions on genomic element
		defineConstant("r", 1/numpos);   	// recombination rate
		defineConstant("alpha", 0.1); 	 // mutation effect size

		defineConstant("outputEvery", 5); // how often to output data
		
	}
		
	 // init mutation rate
	initializeMutationRate(mu);
	
	// init mutation type
	initializeMutationType("m1", 0.5, "n", 0.0, alpha);
	
	// init genomic element
	initializeGenomicElementType("g1", m1, 1);
	initializeGenomicElement(g1,  1,  numpos);
	
	initializeRecombinationRate(r);

	m1.convertToSubstitution = F;
	m1.mutationStackPolicy = "l";  //last 
}

// create a subpopulations of size N individuals
// and set migration rate between them
1 {
	sim.addSubpop("p1", N);
	sim.addSubpop("p2", N);
	
	p1.setMigrationRates(p2, m);
	p2.setMigrationRates(p1, m);	
}

1: late() {
	// create vector of all individuals
	inds  = sim.subpopulations.individuals;
	inds.z = inds.sumOfMutationsOfType(m1);
}

fitness(m1) {
	// the QTLs themselves are neutral; their effect is handled below
	return 1.0;
}


//Gaussian fitness function

fitness(NULL, p1) {
	return exp(-1* (1 - individual.z)^2 / (2*sigsqr));
}

fitness(NULL, p2) {
	return exp(-1* (-1 - individual.z)^2 / (2*sigsqr));
}

1:250000 late(){
	if (sim.generation % outputEvery == 0){
		genomes			= sim.subpopulations.genomes;
		muts			= sortBy(unique(c(genomes.mutationsOfType(m1))), "position");
		mut_freqs_p1	= sim.mutationFrequencies(p1, muts);
		mut_freqs_p2	= sim.mutationFrequencies(p2, muts);

		if (sim.generation / outputEvery == 1){
			out_m 			= "m39 " + "position " + "select_coef " + "p1_freq " + "p2_freq " + "origin_gen " +  "migr_rate " + "mut_rate " + "recomb_rate " + "fitness_width " + "n " + "alpha " + "output_gen " + "\n";
		} else {
			out_m = "";
		}
		
		// if mutations are present output
		if (size(muts) > 0 ) {
			for (i in 0:(size(muts)-1)){
				out_m = out_m + "m39 " + muts.position[i] + " " + muts.selectionCoeff[i] + " " + mut_freqs_p1[i] + " " + mut_freqs_p2[i] + " " + muts.originGeneration[i] + " " + m + " " + mu + " " + r + " " + sigsqr + " " + N + " " + alpha + " " + sim.generation + "\n";
			}
		}
   		print(out_m);
	}	
}