We are going to expand the example from the previous vignette to include multiple species.

Setting up the simulation

As previously, we start by loading the packages and creating the landscape.

library(metaRange)
library(terra)

# find the file
raster_file <- system.file("ex/elev.tif", package = "terra")

# load it
r <- rast(raster_file)

# scale it
r <- scale(r, center = FALSE, scale = TRUE)

r <- rep(r, 10)
landscape <- sds(r)
names(landscape) <- c("habitat_quality")

# plot the first layer of the landscape
plot(landscape[["habitat_quality"]][[1]], main = "Habitat quality")

Figure 1: The habitat quality of the example landscape. Note: higher value = better habitat quality

As before, we create a simulation and add the landscape to it.

sim <- create_simulation(
    source_environment = landscape,
    ID = "example_simulation",
    seed = 1
)
#> number of time steps: 10
#> time step layer mapping: 1, 2, 3, 4, 5, 6, 7, 8, 9, 10
#> added environment
#> class       : SpatRasterDataset 
#> subdatasets : 1 
#> dimensions  : 90, 95 (nrow, ncol)
#> nlyr        : 10 
#> resolution  : 0.008333333, 0.008333333  (x, y)
#> extent      : 5.741667, 6.533333, 49.44167, 50.19167  (xmin, xmax, ymin, ymax)
#> coord. ref. : lon/lat WGS 84 (EPSG:4326) 
#> source(s)   : memory 
#> names       : habitat_quality
#> 
#> created simulation: example_simulation

Adding more species to the simulation

Instead of only adding one species to the simulation, we can just supply more names in the add_species() call to create more species.

sim$add_species(c("species_1", "species_2"))
#> adding species
#> name: species_1
#> name: species_2

If you are at any point wondering which, or how many species are in the simulation, you can use the species_names() method.

sim$species_names()
#> [1] "species_2" "species_1"

Adding traits to (multiple) species

The add_traits() method is able to add (multiple) traits to multiple species at once, which is useful when setting up a large number of species with the same traits. So instead of only specifying one species argument, we can specify a vector of species names, all of which will get the same traits. Here, we set the initial abundance of the species to be proportional to the habitat quality.

sim$add_traits(
    species = c("species_1", "species_2"),
    population_level = TRUE,
    abundance = 200 * sim$environment$current[["habitat_quality"]]
)
#> adding traits:
#> [1] "abundance"
#> 
#> to species:
#> [1] "species_1" "species_2"
#>

As the name suggests, sim$environment$current[["....."]] always refers to the “current” state of the landscape. So what is the difference to the source SDS we used as input to create the simulation? * The sim$environment$sourceSDS[["....."]] is stored as raster data, has multiple sub-datasets and multiple layer and is potentially stored on the disk, which makes it suitable to store large amounts of data (the whole time series), but makes accessing it more complicated and slow. * The current environment contains a 2 dimensional matrix (i.e. only one layer) with the same name for each of the sub-datasets in the source SDS and is stored in memory. This makes it faster to access and easier to use in calculations (As seen in the code example above). This current environment is automatically updated at the beginning of each time step, and right now (before the simulation has started) stores the condition of the first time step (i.e. the first layer of each sub-dataset of the sourceSDS).

# define a nice color palette
plot_cols <- hcl.colors(100, "BluYl", rev = TRUE)
plot(
    sim,
    obj = "species_1",
    name = "abundance",
    main = "Initial abundance",
    col = plot_cols
)

Figure 2: The initial abundance of the species.

If we would want to add a trait to all species in the simulation, without having to type their names, we could use the already mentioned species_names() method to get a vector of all species names and then use that as the species argument.

sim$add_traits(
    species = sim$species_names(),
    population_level = TRUE,
    reproduction_rate = 1.5,
    carrying_capacity = 1000,
    allee_threshold = 150
)
#> adding traits:
#> [1] "reproduction_rate" "carrying_capacity" "allee_threshold"
#> 
#> to species:
#> [1] "species_2" "species_1"
#>

Since we only have the two species in the simulation this would be equivalent to the previous call.

Adding processes

Until now, we have two species in the simulation that are virtually identical. In order to make them behave differently, we can add different processes to them.

Reproduction

In the case of species 1, we will use the same reproduction model from the previous vignette ricker_reproduction_model and we will let the habitat quality influence the carrying capacity.

sim$add_process(
    species = "species_1",
    process_name = "reproduction",
    process_fun = function() {
        ricker_reproduction_model(
            self$traits$abundance,
            self$traits$reproduction_rate,
            self$traits$carrying_capacity * self$sim$environment$current$habitat_quality
        )
        print(
            paste0(self$name, " mean abundance: ", mean(self$traits$abundance))
        )
    },
    execution_priority = 1
)
#> adding process: reproduction
#> to species:
#> [1] "species_1"
#>

In the case of species 2, we will use a Ricker model with additional Allee effects (via the function: ricker_allee_reproduction_model), which is an adapted version from the model described in: Cabral, J.S. and Schurr, F.M. (2010) [Ref. 1]. This means the populations that are smaller than the Allee threshold will have a negative per-capacity reproduction rate and go extinct over time. The Allee effects is also known as “depensation” or “negative density dependence” and can describe multiple different mechanisms that lead to lower reproduction rates at small population sizes as for example difficulties finding a mate, or increased predation pressure (see: Liermann and Hilborn, 2001) [Ref. 2].

sim$add_process(
    species = "species_2",
    process_name = "reproduction",
    process_fun = function() {
        self$traits$abundance <-
            ricker_allee_reproduction_model(
                self$traits$abundance,
                self$traits$reproduction_rate,
                self$traits$carrying_capacity * self$sim$environment$current$habitat_quality,
                self$traits$allee_threshold
            )
        print(
            paste0(self$name, " mean abundance: ", mean(self$traits$abundance))
        )
    },
    execution_priority = 1
)
#> adding process: reproduction
#> to species:
#> [1] "species_2"
#>

Executing the simulation

Now we can execute the simulation again and compare the results.

set_verbosity(0)
sim$begin()
#> [1] "species_1 mean abundance: 348.709489547274"
#> [1] "species_2 mean abundance: 113.066829172314"
#> [1] "species_1 mean abundance: 577.208810253916"
#> [1] "species_2 mean abundance: 125.016153630418"
#> [1] "species_1 mean abundance: 497.450380346053"
#> [1] "species_2 mean abundance: 143.677696198853"
#> [1] "species_1 mean abundance: 538.40186195848"
#> [1] "species_2 mean abundance: 174.359283926975"
#> [1] "species_1 mean abundance: 518.397874403992"
#> [1] "species_2 mean abundance: 226.71976346029"
#> [1] "species_1 mean abundance: 528.487508675932"
#> [1] "species_2 mean abundance: 310.122736448007"
#> [1] "species_1 mean abundance: 523.467949256324"
#> [1] "species_2 mean abundance: 391.190632359679"
#> [1] "species_1 mean abundance: 525.983603165407"
#> [1] "species_2 mean abundance: 422.940404661898"
#> [1] "species_1 mean abundance: 524.727298281258"
#> [1] "species_2 mean abundance: 451.072431872118"
#> [1] "species_1 mean abundance: 525.355824461575"
#> [1] "species_2 mean abundance: 458.686829593375"

plot(
    sim,
    obj = "species_1",
    name = "abundance",
    main = "Species 1: abundance",
    col = plot_cols
)

Figure 3: The resulting abundance distribution of species 1 after 10 simulation time steps.

plot(
    sim,
    obj = "species_2",
    name = "abundance",
    main = "Species 2: abundance",
    col = plot_cols
)

Figure 4: The resulting abundance distribution of species 2 after 10 simulation time steps.

Note how in areas of lower habitat quality, the populations of species 2 are extinct since their abundance was lower than the Allee threshold. If we would combine this with a dispersal process, this could lead to areas that are colonized by species 2, but are permanent population “sinks” for the species, since they would depend on immigration from other areas and are not self-sustaining.

References:

Cabral, J.S. and Schurr, F.M. (2010) Estimating demographic models for the range dynamics of plant species. Global Ecology and Biogeography, 19, 85–97. doi:10.1111/j.1466-8238.2009.00492.x
Liermann, and Hilborn, (2001), Depensation: evidence, models and implications. Fish and Fisheries, 2: 33–58. doi:10.1046/j.1467-2979.2001.00029.x