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)