library(slimr)
#> Welcome to the slimr package for forward population genetics simulation in SLiM. For more information on SLiM please visit https://messerlab.org/slim/ .
#>
#> Attaching package: 'slimr'
#> The following object is masked from 'package:methods':
#>
#> initialize
#> The following object is masked from 'package:base':
#>
#> interaction
This example is based on recipe 19.5 (spatial host-parasite model)
and 19.6 (evolutionary trait-matching host-parasite mode) from the SLiM
Manual, a spatial trait-matching host parasite model. We will use all of
the main slimr
verbs in this model. We will then show how
slimr
can be used in conjunction with SLiMGUI to aid in
model development.
For the purposes of modifying existing SLiM models in
slimr
, we can use some neat slimr
functionality. slimr
includes all of the recipes for SLiM
scripts found in the SLiM Manual. to look at one of the recipes you can
simply access the slim_recipes
. Here we simply print out
the recipe from section 5.3.4 of the SLiM Manual:
cat(slim_recipes$`5.3.4`)
#> // Keywords: migration, dispersal, SLiMgui
#>
#> initialize() {
#> initializeMutationRate(1e-7);
#> initializeMutationType("m1", 0.5, "f", 0.0);
#> initializeMutationType("m2", 0.5, "f", 0.3);
#> initializeGenomicElementType("g1", m1, 1.0);
#> initializeGenomicElement(g1, 0, 99999);
#> initializeRecombinationRate(1e-8);
#> }
#> 1 late() {
#> mSide = 10; // number of subpops along one side of the grid
#> for (i in 1:(mSide * mSide))
#> sim.addSubpop(i, 500);
#>
#> subpops = sim.subpopulations;
#> for (x in 1:mSide)
#> for (y in 1:mSide)
#> {
#> destID = (x - 1) + (y - 1) * mSide + 1;
#> ds = subpops[destID - 1];
#> if (x > 1) // left to right
#> ds.setMigrationRates(destID - 1, runif(1, 0.0, 0.05));
#> if (x < mSide) // right to left
#> ds.setMigrationRates(destID + 1, runif(1, 0.0, 0.05));
#> if (y > 1) // top to bottom
#> ds.setMigrationRates(destID - mSide, runif(1, 0.0, 0.05));
#> if (y < mSide) // bottom to top
#> ds.setMigrationRates(destID + mSide, runif(1, 0.0, 0.05));
#>
#> // set up SLiMgui's population visualization nicely
#> xd = ((x - 1) / (mSide - 1)) * 0.9 + 0.05;
#> yd = ((y - 1) / (mSide - 1)) * 0.9 + 0.05;
#> ds.configureDisplay(c(xd, yd), 0.4);
#> }
#>
#> // remove 25% of the subpopulations
#> subpops[sample(0:99, 25)].setSubpopulationSize(0);
#>
#> // introduce a beneficial mutation
#> target_subpop = sample(sim.subpopulations, 1);
#> sample(target_subpop.genomes, 10).addNewDrawnMutation(m2, 20000);
#> }
#> 10000 late() { sim.outputFixedMutations(); }
SLiM recipes can be converted into slimr_script
object
using the function as_slimr_script()
.
test_script <- as_slimr_script(slim_recipes$`5.3.4`)
test_script
#> <slimr_script[3]>
#> block_init_1:initialize() {
#> initializeMutationRate(1e-07);
#> initializeMutationType("m1", 0.5, "f", asFloat(0));
#> initializeMutationType("m2", 0.5, "f", 0.3);
#> initializeGenomicElementType("g1", m1, asFloat(1));
#> initializeGenomicElement(g1, 0, 99999);
#> initializeRecombinationRate(1e-08);
#> }
#>
#> block_2:1 late() {
#> mSide = 10;
#> for (i in 1:(mSide * mSide)) sim.addSubpop(i, 500);
#> subpops = sim.subpopulations;
#> for (x in 1:mSide) for (y in 1:mSide) {
#> destID = (x - 1) + (y - 1) * mSide + 1;
#> ds = subpops[destID - 1];
#> if (x > 1) ds.setMigrationRates(destID - 1, runif(1, asFloat(0), 0.05));
#> if (x < mSide) ds.setMigrationRates(destID + 1, runif(1, asFloat(0), 0.05));
#> if (y > 1) ds.setMigrationRates(destID - mSide, runif(1, asFloat(0), 0.05));
#> if (y < mSide) ds.setMigrationRates(destID + mSide, runif(1, asFloat(0), 0.05));
#> xd = ((x - 1)/(mSide - 1)) * 0.9 + 0.05;
#> yd = ((y - 1)/(mSide - 1)) * 0.9 + 0.05;
#> ds.configureDisplay(c(xd, yd), 0.4);
#> }
#> subpops[sample(0:99, 25)] %.% setSubpopulationSize(0);
#> target_subpop = sample(sim.subpopulations, 1);
#> sample(target_subpop.genomes, 10) %.% addNewDrawnMutation(m2, 20000);
#> }
#>
#> block_3:10000 late() {
#> sim.outputFixedMutations();
#> }
However, we cannot edit or modify a slimr_script
once it
has been created. To edit and mix and match between different
slimr_script
objects we need to be able to reconstruct the
slimr_script()
call that could recreate the above.
slimr
provides the reconstruct
function for
this purpose:
cat(reconstruct(test_script))
#> slim_script(
#>
#> slim_block(initialize(), {
#> initializeMutationRate(1e-07)
#> initializeMutationType("m1", 0.5, "f", asFloat(0))
#> initializeMutationType("m2", 0.5, "f", 0.3)
#> initializeGenomicElementType("g1", m1, asFloat(1))
#> initializeGenomicElement(g1, 0, 99999)
#> initializeRecombinationRate(1e-08)
#> }),
#>
#> slim_block(1, late(), {
#> mSide = 10
#> for (i in 1:(mSide * mSide)) sim.addSubpop(i, 500)
#> subpops = sim.subpopulations
#> for (x in 1:mSide) for (y in 1:mSide) {
#> destID = (x - 1) + (y - 1) * mSide + 1
#> ds = subpops[destID - 1]
#> if (x > 1) ds.setMigrationRates(destID - 1, runif(1, asFloat(0), 0.05))
#> if (x < mSide) ds.setMigrationRates(destID + 1, runif(1, asFloat(0), 0.05))
#> if (y > 1) ds.setMigrationRates(destID - mSide, runif(1, asFloat(0), 0.05))
#> if (y < mSide) ds.setMigrationRates(destID + mSide, runif(1, asFloat(0), 0.05))
#> xd = ((x - 1)/(mSide - 1)) * 0.9 + 0.05
#> yd = ((y - 1)/(mSide - 1)) * 0.9 + 0.05
#> ds.configureDisplay(c(xd, yd), 0.4)
#> }
#> subpops[sample(0:99, 25)] %.% setSubpopulationSize(0)
#> target_subpop = sample(sim.subpopulations, 1)
#> sample(target_subpop.genomes, 10) %.% addNewDrawnMutation(m2, 20000)
#> }),
#>
#> slim_block(10000, late(), {
#> sim.outputFixedMutations()
#> })
#> )
Now we can simply copy and paste that output into an R script, which
we can edit as we like. We can do the two steps above in a single step
using as_slimr_code()
, which takes a SLiM script as a
character vector and returns the slimr
code equivalent, and
printing it to the console.
as_slimr_code(slim_recipes$`5.3.4`)
#> slim_script(
#>
#> slim_block(initialize(), {
#> initializeMutationRate(1e-07)
#> initializeMutationType("m1", 0.5, "f", asFloat(0))
#> initializeMutationType("m2", 0.5, "f", 0.3)
#> initializeGenomicElementType("g1", m1, asFloat(1))
#> initializeGenomicElement(g1, 0, 99999)
#> initializeRecombinationRate(1e-08)
#> }),
#>
#> slim_block(1, late(), {
#> mSide = 10
#> for (i in 1:(mSide * mSide)) sim.addSubpop(i, 500)
#> subpops = sim.subpopulations
#> for (x in 1:mSide) for (y in 1:mSide) {
#> destID = (x - 1) + (y - 1) * mSide + 1
#> ds = subpops[destID - 1]
#> if (x > 1) ds.setMigrationRates(destID - 1, runif(1, asFloat(0), 0.05))
#> if (x < mSide) ds.setMigrationRates(destID + 1, runif(1, asFloat(0), 0.05))
#> if (y > 1) ds.setMigrationRates(destID - mSide, runif(1, asFloat(0), 0.05))
#> if (y < mSide) ds.setMigrationRates(destID + mSide, runif(1, asFloat(0), 0.05))
#> xd = ((x - 1)/(mSide - 1)) * 0.9 + 0.05
#> yd = ((y - 1)/(mSide - 1)) * 0.9 + 0.05
#> ds.configureDisplay(c(xd, yd), 0.4)
#> }
#> subpops[sample(0:99, 25)] %.% setSubpopulationSize(0)
#> target_subpop = sample(sim.subpopulations, 1)
#> sample(target_subpop.genomes, 10) %.% addNewDrawnMutation(m2, 20000)
#> }),
#>
#> slim_block(10000, late(), {
#> sim.outputFixedMutations()
#> })
#> )
Here is the slimr
code to create recipes 19.5 and
19.6
as_slimr_code(slim_recipes$`19.5`)
#> slim_script(
#>
#> slim_block(all = initialize(), {
#> defineConstant("K", 100)
#> defineConstant("R", log(20))
#> defineConstant("A", 0.015)
#> defineConstant("SIDE", 10)
#> defineConstant("S", SIDE * SIDE)
#> defineConstant("N0_host", asInteger((135.6217 + 0.01) * S))
#> defineConstant("N0_parasitoid", asInteger((109.301 + 0.01) * S))
#> defineConstant("S_P", 0.5)
#> defineConstant("S_H", 0.2)
#> defineConstant("D_H", 0.2)
#> defineConstant("CROSS_SCRIPT", "subpop.addCrossed(individual, mate);")
#> initializeSLiMModelType("nonWF")
#> initializeInteractionType(1, "xy", maxDistance = S_P)
#> i1.setInteractionFunction("l", asFloat(1))
#> initializeInteractionType(2, "xy", maxDistance = S_H)
#> i2.setInteractionFunction("l", asFloat(1))
#> }),
#>
#> slim_block(host = initialize(), {
#> initializeSpecies(avatar = "🐛", color = "cornflowerblue")
#> initializeSLiMOptions(dimensionality = "xy")
#> }),
#>
#> slim_block(parasitoid = initialize(), {
#> initializeSpecies(avatar = "🦟", color = "red")
#> initializeSLiMOptions(dimensionality = "xy")
#> }),
#>
#> slim_block(all = 2, 250, first(), {
#> host_pop = host.subpopulations
#> hosts = host_pop.individuals
#> parasitoid_pop = parasitoid.subpopulations
#> parasitoids = parasitoid_pop.individuals
#> i1.evaluate(c(host_pop, parasitoid_pop))
#> i2.evaluate(host_pop)
#> parasitoid_density_byhost = i1.localPopulationDensity(hosts, parasitoid_pop)
#> parasitoids.tag = 0
#> parasitoids.setValue("PREY_POS", NULL)
#> P_parasitized_byhost = 1 - exp(-A * parasitoid_density_byhost)
#> killed = (runif(hosts.size()) < P_parasitized_byhost)
#> hosts.tag = asInteger(killed)
#> preys = hosts[killed]
#> for (prey in preys) {
#> hunter = i1.drawByStrength(prey, 1, parasitoid_pop)
#> preyPos = prey.spatialPosition
#> preyPos = c(hunter.getValue("PREY_POS"), preyPos)
#> hunter.tag = hunter.tag + 1
#> hunter.setValue("PREY_POS", preyPos)
#> }
#> unhunted = hosts[!killed]
#> host_density_by_unhunted = i2.localPopulationDensity(unhunted, host_pop)
#> P_survives_by_unhunted = exp(-host_density_by_unhunted/K)
#> survived = (runif(unhunted.size()) < P_survives_by_unhunted)
#> dead = unhunted[!survived]
#> dead.tag = 1
#> }),
#>
#> slim_block(host = reproduction(), {
#> if (individual.tag != 0) return
#> mate = i2.drawByStrength(individual)
#> if (mate.size()) {
#> litterSize = rpois(1, exp(R))
#> if (litterSize > 0) {
#> offspring = sapply(seqLen(litterSize), CROSS_SCRIPT)
#> positions = rep(individual.spatialPosition, litterSize)
#> positions = positions + rnorm(litterSize * 2, 0, D_H)
#> positions = p1.pointReflected(positions)
#> offspring.setSpatialPosition(positions)
#> }
#> }
#> }),
#>
#> slim_block(parasitoid = reproduction(), {
#> litterSize = individual.tag
#> if (litterSize > 0) {
#> mate = i1.drawByStrength(individual)
#> if (mate.size()) {
#> offspring = sapply(seqLen(litterSize), CROSS_SCRIPT)
#> offspring.setSpatialPosition(individual.getValue("PREY_POS"))
#> }
#> }
#> }),
#>
#> slim_block(all = 1, early(), {
#> host.addSubpop("p1", N0_host)
#> p1.setSpatialBounds(c(0, 0, SIDE, SIDE))
#> p1.individuals.setSpatialPosition(p1.pointUniform(N0_host))
#> parasitoid.addSubpop("p2", N0_parasitoid)
#> p2.setSpatialBounds(c(0, 0, SIDE, SIDE))
#> p2.individuals.setSpatialPosition(p2.pointUniform(N0_parasitoid))
#> }),
#>
#> slim_block(host = survival(), {
#> return((individual.age == 0))
#> }),
#>
#> slim_block(parasitoid = survival(), {
#> return((individual.age == 0))
#> }),
#>
#> slim_block(all = 250, late(), {
#>
#> })
#> )
as_slimr_code(slim_recipes$`19.6`)
#> slim_script(
#>
#> slim_block(all = initialize(), {
#> defineConstant("K", 100)
#> defineConstant("R", log(20))
#> defineConstant("A", 0.015)
#> defineConstant("S", 10^2)
#> defineConstant("N0_host", asInteger((135.6217 + 0.01) * S))
#> defineConstant("N0_parasitoid", asInteger((109.301 + 0.01) * S))
#> defineConstant("S_M", asFloat(1))
#> defineConstant("S_S", asFloat(2))
#> initializeSLiMModelType("nonWF")
#> }),
#>
#> slim_block(host = initialize(), {
#> initializeSpecies(avatar = "🐛", color = "cornflowerblue")
#> initializeMutationType("m1", 0.5, "n", asFloat(0), 0.1)
#> initializeGenomicElementType("g1", m1, asFloat(1))
#> initializeGenomicElement(g1, 0, 9999)
#> initializeMutationRate(1e-07)
#> initializeRecombinationRate(1e-08)
#> }),
#>
#> slim_block(parasitoid = initialize(), {
#> initializeSpecies(avatar = "🦟", color = "red")
#> initializeMutationType("m2", 0.5, "n", asFloat(0), 0.1)
#> initializeGenomicElementType("g2", m2, asFloat(1))
#> initializeGenomicElement(g2, 0, 9999)
#> initializeMutationRate(1e-07)
#> initializeRecombinationRate(1e-08)
#> }),
#>
#> slim_block(host = mutationEffect(m1), {
#> return(asFloat(1))
#> }),
#>
#> slim_block(parasitoid = mutationEffect(m2), {
#> return(asFloat(1))
#> }),
#>
#> slim_block(all = 2, 10000, first(), {
#> hosts = host.subpopulations.individuals
#> parasitoids = parasitoid.subpopulations.individuals
#> x1 = hosts.size()/S
#> x2 = parasitoids.size()/S
#> host_values = hosts.tagF
#> parasitoid_values = parasitoids.tagF
#> mean_parasitoid = mean(parasitoid_values)
#> scale = dnorm(asFloat(0), asFloat(0), S_M)
#> host_match = dnorm(host_values, mean_parasitoid, S_M)/scale
#> parasitoids.tag = 0
#> P_parasitized_byhost = 1 - exp(-A * x2 * host_match)
#> killed = rbinom(hosts.size(), 1, P_parasitized_byhost)
#> hosts.tag = killed
#> hunters = sapply(hosts[killed == 1], "sample(parasitoids, 1, " + "weights=dnorm(applyValue.tagF - parasitoid_values, asFloat(0.0), S_M));")
#> for (hunter in hunters) hunter.tag = hunter.tag + 1
#> survivors = hosts[killed == 0]
#> P_survives = exp(-x1/K)
#> survived = rbinom(survivors.size(), 1, P_survives)
#> dead = survivors[survived == 0]
#> dead.tag = 1
#> }),
#>
#> slim_block(host = reproduction(), {
#> if (individual.tag != 0) return
#> litterSize = rpois(1, exp(R))
#> if (litterSize > 0) {
#> mate = subpop.sampleIndividuals(1, exclude = individual)
#> for (i in seqLen(litterSize)) subpop.addCrossed(individual, mate)
#> }
#> }),
#>
#> slim_block(parasitoid = reproduction(), {
#> litterSize = individual.tag
#> if (litterSize > 0) {
#> mate = subpop.sampleIndividuals(1, exclude = individual)
#> for (i in seqLen(litterSize)) subpop.addCrossed(individual, mate)
#> }
#> }),
#>
#> slim_block(all = 1, early(), {
#> host.addSubpop("p1", N0_host)
#> parasitoid.addSubpop("p2", N0_parasitoid)
#> log = community.createLogFile("host-parasite log.txt", logInterval = 1)
#> log.addTick()
#> log.addPopulationSize(host)
#> log.addMeanSDColumns("host", "p1.individuals.tagF;")
#> log.addPopulationSize(parasitoid)
#> log.addMeanSDColumns("parasitoid", "p2.individuals.tagF;")
#> }),
#>
#> slim_block(all = early(), {
#> scale = dnorm(asFloat(0), asFloat(0), S_S)
#> hosts = host.subpopulations.individuals
#> phenotypes = hosts.sumOfMutationsOfType(m1)
#> hosts.fitnessScaling = dnorm(phenotypes, asFloat(0), S_S)/scale
#> hosts.tagF = phenotypes
#> parasitoids = parasitoid.subpopulations.individuals
#> phenotypes = parasitoids.sumOfMutationsOfType(m2)
#> parasitoids.fitnessScaling = dnorm(phenotypes, asFloat(0), S_S)/scale
#> parasitoids.tagF = phenotypes
#> hosts[hosts.age > 0] %.% fitnessScaling = asFloat(0)
#> parasitoids[parasitoids.age > 0] %.% fitnessScaling = asFloat(0)
#> }),
#>
#> slim_block(all = 10000, late(), {
#>
#> })
#> )
Combining aspects of the above two slimr
scripts, and
adding slimr
verb functionality gives a final full
model:
n_gen <- 1000
slim_script(
slim_block(all = initialize(), {
r_template_constant("K", 100)
r_template_constant("R", log(20))
r_template_constant("A", 0.015)
r_template_constant("SIDE", 10)
defineConstant("S", SIDE * SIDE)
defineConstant("N0_host", asInteger((135.6217 + 0.01) * S))
defineConstant("N0_parasitoid", asInteger((109.301 + 0.01) * S))
r_template_constant("S_P", 0.5)
r_template_constant("S_H", 0.2)
r_template_constant("D_H", 0.2)
defineConstant("CROSS_SCRIPT", "subpop.addCrossed(individual, mate);")
initializeSLiMModelType("nonWF")
initializeInteractionType(1, "xy", maxDistance = S_P)
i1.setInteractionFunction("l", asFloat(1))
initializeInteractionType(2, "xy", maxDistance = S_H)
i2.setInteractionFunction("l", asFloat(1))
}),
slim_block(host = initialize(), {
initializeSpecies(avatar = "🐛", color = "cornflowerblue")
initializeSLiMOptions(dimensionality = "xy")
}),
slim_block(parasitoid = initialize(), {
initializeSpecies(avatar = "🦟", color = "red")
initializeSLiMOptions(dimensionality = "xy")
}),
slim_block(all = 2, !!n_gen, first(), {
host_pop = host.subpopulations
hosts = host_pop.individuals
parasitoid_pop = parasitoid.subpopulations
parasitoids = parasitoid_pop.individuals
i1.evaluate(c(host_pop, parasitoid_pop))
i2.evaluate(host_pop)
parasitoid_density_byhost = i1.localPopulationDensity(hosts, parasitoid_pop)
parasitoids.tag = 0
parasitoids.setValue("PREY_POS", NULL)
P_parasitized_byhost = 1 - exp(-A * parasitoid_density_byhost)
killed = (runif(hosts.size()) < P_parasitized_byhost)
hosts.tag = asInteger(killed)
preys = hosts[killed]
for (prey in preys) {
hunter = i1.drawByStrength(prey, 1, parasitoid_pop)
preyPos = prey.spatialPosition
preyPos = c(hunter.getValue("PREY_POS"), preyPos)
hunter.tag = hunter.tag + 1
hunter.setValue("PREY_POS", preyPos)
}
unhunted = hosts[!killed]
host_density_by_unhunted = i2.localPopulationDensity(unhunted, host_pop)
P_survives_by_unhunted = exp(-host_density_by_unhunted/K)
survived = (runif(unhunted.size()) < P_survives_by_unhunted)
dead = unhunted[!survived]
dead.tag = 1
cat("Host Population Size: ", hosts.size())
cat("Parasitoid Population Size: ", parasitoids.size())
}),
slim_block(host = reproduction(), {
if (individual.tag != 0) return
mate = i2.drawByStrength(individual)
if (mate.size()) {
litterSize = rpois(1, exp(R))
if (litterSize > 0) {
offspring = sapply(seqLen(litterSize), CROSS_SCRIPT)
positions = rep(individual.spatialPosition, litterSize)
positions = positions + rnorm(litterSize * 2, 0, D_H)
positions = p1.pointReflected(positions)
offspring.setSpatialPosition(positions)
}
}
}),
slim_block(parasitoid = reproduction(), {
litterSize = individual.tag
if (litterSize > 0) {
mate = i1.drawByStrength(individual)
if (mate.size()) {
offspring = sapply(seqLen(litterSize), CROSS_SCRIPT)
offspring.setSpatialPosition(individual.getValue("PREY_POS"))
}
}
}),
slim_block(all = 1, early(), {
host.addSubpop("p1", N0_host)
p1.setSpatialBounds(c(0, 0, SIDE, SIDE))
p1.individuals.setSpatialPosition(p1.pointUniform(N0_host))
parasitoid.addSubpop("p2", N0_parasitoid)
p2.setSpatialBounds(c(0, 0, SIDE, SIDE))
p2.individuals.setSpatialPosition(p2.pointUniform(N0_parasitoid))
}),
slim_block(host = survival(), {
return((individual.age == 0))
}),
slim_block(parasitoid = survival(), {
return((individual.age == 0))
}),
slim_block(all = !!n_gen, late(), {
community.simulationFinished()
})
) -> hp_script
hp_script
#> <slimr_script[10]>
#> block_init:species all initialize() {
#> defineConstant("K", ..K..);
#> defineConstant("R", ..R..);
#> defineConstant("A", ..A..);
#> defineConstant("SIDE", ..SIDE..);
#> defineConstant("S", SIDE * SIDE);
#> defineConstant("N0_host", asInteger((135.6217 + 0.01) * S));
#> defineConstant("N0_parasitoid", asInteger((109.301 + 0.01) * S));
#> defineConstant("S_P", ..S_P..);
#> defineConstant("S_H", ..S_H..);
#> defineConstant("D_H", ..D_H..);
#> defineConstant("CROSS_SCRIPT", "subpop.addCrossed(individual, mate);");
#> initializeSLiMModelType("nonWF");
#> initializeInteractionType(1, "xy", maxDistance = S_P);
#> i1.setInteractionFunction("l", asFloat(1));
#> initializeInteractionType(2, "xy", maxDistance = S_H);
#> i2.setInteractionFunction("l", asFloat(1));
#> }
#>
#> block_init:species host initialize() {
#> initializeSpecies(avatar = "🐛", color = "cornflowerblue");
#> initializeSLiMOptions(dimensionality = "xy");
#> }
#>
#> block_init:species parasitoid initialize() {
#> initializeSpecies(avatar = "🦟", color = "red");
#> initializeSLiMOptions(dimensionality = "xy");
#> }
#>
#> block_04:ticks all 2:1000 first() {
#> host_pop = host.subpopulations;
#> hosts = host_pop.individuals;
#> parasitoid_pop = parasitoid.subpopulations;
#> parasitoids = parasitoid_pop.individuals;
#> i1.evaluate(c(host_pop, parasitoid_pop));
#> i2.evaluate(host_pop);
#> parasitoid_density_byhost = i1.localPopulationDensity(hosts, parasitoid_pop);
#> parasitoids.tag = 0;
#> parasitoids.setValue("PREY_POS", NULL);
#> P_parasitized_byhost = 1 - exp(-A * parasitoid_density_byhost);
#> killed = (runif(hosts.size()) < P_parasitized_byhost);
#> hosts.tag = asInteger(killed);
#> preys = hosts[killed];
#> for (prey in preys) {
#> hunter = i1.drawByStrength(prey, 1, parasitoid_pop);
#> preyPos = prey.spatialPosition;
#> preyPos = c(hunter.getValue("PREY_POS"), preyPos);
#> hunter.tag = hunter.tag + 1;
#> hunter.setValue("PREY_POS", preyPos);
#> }
#> unhunted = hosts[!killed];
#> host_density_by_unhunted = i2.localPopulationDensity(unhunted, host_pop);
#> P_survives_by_unhunted = exp(-host_density_by_unhunted/K);
#> survived = (runif(unhunted.size()) < P_survives_by_unhunted);
#> dead = unhunted[!survived];
#> dead.tag = 1;
#> cat("Host Population Size: ", hosts.size());
#> cat("Parasitoid Population Size: ", parasitoids.size());
#> }
#>
#> block_05:species host 1:1000 reproduction() {
#> if (individual.tag != 0) return;
#> mate = i2.drawByStrength(individual);
#> if (mate.size()) {
#> litterSize = rpois(1, exp(R));
#> if (litterSize > 0) {
#> offspring = sapply(seqLen(litterSize), CROSS_SCRIPT);
#> positions = rep(individual.spatialPosition, litterSize);
#> positions = positions + rnorm(litterSize * 2, 0, D_H);
#> positions = p1.pointReflected(positions);
#> offspring.setSpatialPosition(positions);
#> }
#> }
#> }
#>
#> block_06:species parasitoid 1:1000 reproduction() {
#> litterSize = individual.tag;
#> if (litterSize > 0) {
#> mate = i1.drawByStrength(individual);
#> if (mate.size()) {
#> offspring = sapply(seqLen(litterSize), CROSS_SCRIPT);
#> offspring.setSpatialPosition(individual.getValue("PREY_POS"));
#> }
#> }
#> }
#>
#> block_07:ticks all 1 early() {
#> host.addSubpop("p1", N0_host);
#> p1.setSpatialBounds(c(0, 0, SIDE, SIDE));
#> p1.individuals.setSpatialPosition(p1.pointUniform(N0_host));
#> parasitoid.addSubpop("p2", N0_parasitoid);
#> p2.setSpatialBounds(c(0, 0, SIDE, SIDE));
#> p2.individuals.setSpatialPosition(p2.pointUniform(N0_parasitoid));
#> }
#>
#> block_08:species host 1:1000 survival() {
#> return((individual.age == 0));
#> }
#>
#> block_09:species parasitoid 1:1000 survival() {
#> return((individual.age == 0));
#> }
#>
#> block_10:ticks all 1000 late() {
#> community.simulationFinished();
#> }
#> This slimr_script has templating in block(s) block_init for
#> variables K and R and A and SIDE and S_P
#> and S_H and D_H.
Writing a multispecies model in slimr
is simple. You
only have to use a single named argument in each
slim_block()
call to specify the species to which the block
applies. You use all
to specify that the block applies to
all species. Note in particular the use of
r_template_constant()
, a variation of
r_template()
that inserts the templated variable into a
defineConstant()
statement. defineConstant()
assigns the variable into a constant in the SLiM initialization block so
that it can be uses throughout any block of code in the simulation. It
also makes it easier to adjust parameter values manually within SLiMGUI
if we choose to test the model in SLiMGUI with its interactive
features.
test <- slim_run(hp_script, progress = TRUE)