monarchs.slim

species all initialize()
{
	initializeSLiMModelType("nonWF");
	defaults = Dictionary(
		// General params:
		"SEED", getSeed(),
		"RUNTIME", 520,
		"WIDTH", 50.0,
		"HEIGHT", 50.0,  // The migration route consists of an additional 10x this height

		// Monarch params:
		"BUTTERFLY_MILKWEED_DISTANCE", 2.0,
		"CATERPILLAR_MILKWEED_DISTANCE", 0.05,
		"NUM_EGGS", 25,
		"WEEKLY_FLAT_MORTALITY", 0.02,
		"AGE_DEPENDENT_MORTALITY", c(0.0, 0.0, 0.0, 0.0, 0.0, 0.05, 0.15, 0.25, 0.45, 0.65, 0.85, 1.0),
		"SD", 2.0,

		// Milkweed params:
		"MILKWEED_PATCHES", 20,
		"MILKWEED_PER_PATCH", 80,
		"MILKWEED_PATCH_DISTRIBUTION_SIGMA", 0.5,
		"CATERPILLARS_FED_PER_MILKWEED", 2
		);

	// Set up parameters with a user-defined function
	setupParams(defaults);
	setSeed(SEED);

	// Caterpillar foraging interaction.
	initializeInteractionType(1, "xy", reciprocal=T, maxDistance=CATERPILLAR_MILKWEED_DISTANCE);

	// Reproduction part 1 (interaction by which milkweed indexes nearby males).
	initializeInteractionType(2, "xy", reciprocal=T, maxDistance=BUTTERFLY_MILKWEED_DISTANCE);
	i2.setConstraints("exerter", sex="M");

	// Reproduction part 2 (interaction by which females find milkweed).
	initializeInteractionType(3, "xy", reciprocal=T, maxDistance=BUTTERFLY_MILKWEED_DISTANCE);
	i3.setConstraints("receiver", sex="F");
}

species milkweed initialize()
{
	initializeSpecies(avatar="V", color="cornflowerblue");
	initializeSLiMOptions(dimensionality="xy");
}

species monarch initialize()
{
	initializeSpecies(avatar="O", color="red");
	initializeSLiMOptions(dimensionality="xy");
	initializeSex("A");
}

function (object<DataFrame>$)addMilkweed(void)
{
	/*
        Set up the milkweed positions.
        Patches will be distributed according to a uniform random distribution,
        while milkweed within a patch will be distributed according to a normal distribution.
    */
	centroids_x = runif(MILKWEED_PATCHES, 0, WIDTH);
	centroids_y = runif(MILKWEED_PATCHES, 0, HEIGHT);

	count = MILKWEED_PATCHES * MILKWEED_PER_PATCH;
	xs = rep(centroids_x, MILKWEED_PER_PATCH) + rnorm(count, sd=MILKWEED_PATCH_DISTRIBUTION_SIGMA);
	ys = rep(centroids_y, MILKWEED_PER_PATCH) + rnorm(count, sd=MILKWEED_PATCH_DISTRIBUTION_SIGMA);

	return DataFrame("xs", pmax(0.0, pmin(WIDTH, xs)), "ys", pmax(0.0, pmin(HEIGHT, ys)) + HEIGHT * 10);
}

ticks all 1 first()
{
	// Add the milkweed population.
	milkweedData = addMilkweed();
	milkweed.addSubpop("p1", milkweedData.nrow);
	p1.setSpatialBounds(c(0, HEIGHT * 10, WIDTH, HEIGHT * 11));
	p1.individuals.x = milkweedData.getValue("xs");
	p1.individuals.y = milkweedData.getValue("ys");

	// Initialize the monarchs.
	// The monarchs don't have a defined capacity in this model.
	// The simulation will just start with enough caterpillars to eat the milkweed.
	caterpillars_per_milkweed = CATERPILLARS_FED_PER_MILKWEED * 2;
	monarch.addSubpop("p2", milkweedData.nrow * caterpillars_per_milkweed);
	p2.setSpatialBounds(c(0, HEIGHT * 10, WIDTH, HEIGHT * 11));
	p2.individuals.x = rep(milkweedData.getValue("xs"), caterpillars_per_milkweed);
	p2.individuals.y = rep(milkweedData.getValue("ys"), caterpillars_per_milkweed);
	p2.individuals.tagF = 0.0;

	// No pupae or adult butterflies at start.
	monarch.addSubpop("p3", 0); // Pupae.
	p3.setSpatialBounds(c(0, HEIGHT * 10, WIDTH, HEIGHT * 11));

	monarch.addSubpop("p4", 0); // Adult butterflies.
	p4.setSpatialBounds(c(0, HEIGHT * 10, WIDTH, HEIGHT * 11));

	monarch.addSubpop("p5", 0); // Migratory butterflies.
	p5.setSpatialBounds(c(0, 0, WIDTH, HEIGHT * 11));
}

ticks all first()
{
	// Evaluate the two-part spatial interaction for reproduction.
	i2.evaluate(c(p1, p4));
	i3.evaluate(c(p1, p4));
}

species monarch reproduction(NULL, "F")
{
	// Females and males can reproduce if they are both within range of the same milkweed.
	for (plant in p1.individuals)
	{
		nearby_males = i2.nearestInteractingNeighbors(plant, p4.individualCount, p4);
		plant.setValue("nearby_males", nearby_males.index);
	}

	females = p4.subsetIndividuals(sex="F");
	eggs_laid = rpois(females.size(), NUM_EGGS);
	non_zero = (eggs_laid != 0);
	females = females[non_zero];
	eggs_laid = eggs_laid[non_zero];
	mate_pool = p4.individuals;

	for (mother in females, egg_count in eggs_laid)
	{
		// Females not near any milkweed cannot reproduce.
		nearest_milkweeds = i3.nearestInteractingNeighbors(mother, p1.individualCount, p1);
		if (!size(nearest_milkweeds))
			next;

		// Females with no available mates cannot reproduce.
		nearby_males = nearest_milkweeds.getValue("nearby_males");
		if (!size(nearby_males))
			next;

		// Sample a mate and generate offspring.
		mate = mate_pool[sample(nearby_males, 1)];
		nearest_milkweeds.tag = 0;

		for (i in seqLen(egg_count))
		{
			milkweed_for_this_egg = sample(nearest_milkweeds, 1);
			offspring = p2.addCrossed(mother, mate);
			offspring.setSpatialPosition(milkweed_for_this_egg.spatialPosition);
			offspring.tagF = 0.0;

			// Use tag to prevent this female from laying more than 10 eggs per milkweed in range.
			milkweed_for_this_egg.tag = milkweed_for_this_egg.tag + 1;

			if (milkweed_for_this_egg.tag == 10)
			{
				nearest_milkweeds = nearest_milkweeds[nearest_milkweeds.tag < 10];
				if (!size(nearest_milkweeds))
					break;
			}
		}
	}

	self.active = 0;
}

ticks all early()
{
	week = community.tick % 52;

	// After two weeks, caterpillars become pupae.
	// Pupae survival is based on food collected as a caterpillar.
	// Pupae mortality is only applied a single time.
	new_pupae = p2.subsetIndividuals(minAge=2, maxAge=2);

	mortality_flags = runif(size(new_pupae)) > new_pupae.tagF;
	dead = new_pupae[mortality_flags];
	alive = new_pupae[!mortality_flags];
	monarch.killIndividuals(dead);
	p3.takeMigrants(alive);

	// After two weeks as a pupa, survivors become butterflies.
	new_butterflies = p3.subsetIndividuals(minAge=4, maxAge=4);

	if (week < 24)
	{
		p4.takeMigrants(new_butterflies);
	}
	else
	{
		p5.takeMigrants(new_butterflies);
		p5.takeMigrants(p4.individuals);
	}

	// Non-migratory butterfly dispersal and mortality.
	nonmigrants = p4.individuals;

	p4.deviatePositions(NULL, "reprising", INF, "n", SD);
	age_mortality = AGE_DEPENDENT_MORTALITY[nonmigrants.age];
	nonmigrants.fitnessScaling = pmax(0.0, 1.0 - age_mortality - WEEKLY_FLAT_MORTALITY);

	// Caterpillar foraging.
	// Each milkweed evenly provides resources to all caterpillars within range.
	inds_fed_per_milkweed = (week > 24) ? 0 else CATERPILLARS_FED_PER_MILKWEED;
	i1.evaluate(c(p2, p1));

	for (plant in p1.individuals)
	{
		customers = i1.nearestNeighbors(plant, p2.individualCount, p2);
		customers.tagF = customers.tagF + inds_fed_per_milkweed / size(customers);
	}

	// Migratory butterflies.
	migrants = p5.individuals;
	migrants.fitnessScaling = 1.0 - WEEKLY_FLAT_MORTALITY;

	// Migrate south.
	if (week < 40 & week > 24)
	{
		southbound = migrants[migrants.y > HEIGHT];
		southbound.x = runif(size(southbound), 0.0, WIDTH);
		southbound.y = pmax(0.0, southbound.y - rnorm(size(southbound), HEIGHT, SD * 5));
	}

	// Migrate back north. When individuals make it back north,
	// their ages will be set to a lower value and they will start
	// experiencing age-dependent mortality and start reproducing.
	if (week > 42 | week < 24)
	{
		migrants.x = runif(size(migrants), 0.0, WIDTH);
		migrants.y = pmin(HEIGHT * 11, migrants.y + rnorm(size(migrants), HEIGHT, SD * 5));
		returned_north = migrants[migrants.y > HEIGHT * 10];
		returned_north.age = 7;
		p4.takeMigrants(returned_north);
	}

	// Every year, the position of the milkweed is re-randomized.
	// Assuming the simulation starts in summer, this will occur around spring.
	if (week == 32)
	{
		milkweedData = addMilkweed();
		plants = p1.individuals;
		plants.x = milkweedData.getValue("xs");
		plants.y = milkweedData.getValue("ys");
	}
}

ticks all late()
{
	// Output.
	catn("Time: " + community.tick);
	catn("    Caterpillars: " + p2.individualCount);
	catn("    Pupae:        " + p3.individualCount);
	catn("    Butterflies:  " + p4.individualCount);
	catn("    Migrators:    " + p5.individualCount);
}

ticks all RUNTIME late()
{
	catn("End of simulation (run time reached)");
	community.simulationFinished();
}

function (void)setupParams(object<Dictionary>$ defaults)
{
	if (!exists("PARAMFILE")) defineConstant("PARAMFILE", "./params.json");
	if (!exists("OUTDIR")) defineConstant("OUTDIR", ".");
	defaults.addKeysAndValuesFrom(Dictionary("PARAMFILE", PARAMFILE, "OUTDIR", OUTDIR));

	if (fileExists(PARAMFILE)) {
		defaults.addKeysAndValuesFrom(Dictionary(readFile(PARAMFILE)));
		defaults.setValue("READ_FROM_PARAMFILE", PARAMFILE);
	}

	defaults.setValue("OUTBASE", OUTDIR + "/out_" +	defaults.getValue("SEED"));
	defaults.setValue("OUTPATH", defaults.getValue("OUTBASE") + ".trees");

	for (k in defaults.allKeys) {
		if (!exists(k))
			defineConstant(k, defaults.getValue(k));
		else
			defaults.setValue(k, executeLambda(k + ";"));
	}

	// print out default values
	catn("===========================");
	catn("Model constants: " + defaults.serialize("pretty"));
	catn("===========================");
}