Background

In ecology, a metacommunity is a network of communities of organisms that are interconnected by dispersal and colonization dynamics (Leibold et al. 2004). Here I describe how to use a metacommunity simulation package (MCSim) for the R statistical language to explore how disturbance can interact with metacommunity characteristics to result in different emergent patterns in biodiversity. For an introudction to MCSim, see the tutorial “An introduction to MetaCommunity Simulations (MCSim) for R”. For more information about the theoretical underpinnings of MCSim, please see Sokol et al. (2017) for a sensitivity analysis of how biodiversity metrics respond to changes in model parameters, and Sokol et al. (2015) and Sokol et al. (2020) for applications of the simulation platform to create predictions that can be tested with empirically collected field data.

This simulation is a zero-sum, individual-oriented, iterative, lottery-based simulation. That means that the number of individuals in the simulated metacommunity remains constant through time and each time step is a complete, synchronous generation for the entire metacommunity. A lottery process moves the simulation forward in time with each time step, allowing for stochastic metacommunity dynamics to occur. However, environmental filtering can also be built into the simulation, allowing for deterministic species sorting to occur.

To simulate disturbance, we can string together multiple MCSim simulation scenarios in series to create changing environmental conditions through time. Examples of how disturbance can be simulated include altering environmental filters to create a press distrubance, or creating a temporary decrease in carrying capacity to simulate increased mortality and a temporary reduction in the population. Disturbance frequency can be adjusted based on the number of timesteps each landscape state is assigned, and the number of scenarios that are strung together in series. Note that there are limitations in how this approach can be applied because of the zero-sum assumption in the model, however, the modularity and flexibility of MCSim provides many options for exploring how metacommunity characteristics can interact with perturbations in the landscape.

Please find the most recent releases of MCSim here, and the current version under development at https://github.com/sokole/MCSim.

Getting started

In order to install the current version of MCSim from GitHub, you will need to have the devtools package installed. install.packages("devtools"). The remotes package provides another option for installing packages from GitHub.

Once devtools is installed, you can install the current version of MCSim

# -- Install the current dev version of MCSim
devtools::install_github('sokole/MCSim')

# -- Install the version used in this tutorial
devtools::install_github('sokole/MCSim@v0.5')
# v0.5 is used in this demo

After installing the package, make sure you can load it in your R environment.

library(MCSim)

Note that we will also use functions from the tidyverse, so be sure to make sure you have installed tidyverse from CRAN and load the library:

## -- Attaching packages --------------------------------------- tidyverse 1.3.0 --
## v ggplot2 3.3.5     v purrr   0.3.4
## v tibble  3.0.3     v dplyr   1.0.1
## v tidyr   1.1.1     v stringr 1.4.0
## v readr   1.3.1     v forcats 0.5.0
## -- Conflicts ------------------------------------------ tidyverse_conflicts() --
## x dplyr::filter() masks stats::filter()
## x dplyr::lag()    masks stats::lag()

Running simulations with dynamic landscapes

There are two steps to creating a metacommunity simulation in MCSim:

  1. Make a “landscape” using make.landscape
  2. Run a simulation using metasim

Note that multiple landscapes can be concatenated in series to create a landscape with properties that shift overtime. In our first example, we will concatenate two landscapes to create a change in the environmental characteristics for a subset of sites in the landscape.

Example 1. Change in environmental state at a subset of sites

Step 1. Create the initial species pool

First, we will set up the characteristics of the taxa in the species pool for the simulation, including their habitat preferences (niche.positions) and flexibility in the diversity of habitats in which they can survive (niche.breadths).

set.seed(1234) #set random seed

# niche positions, niche breadths, and relative abundances for three species
niche.positions <-  c(-.8, 0, .9)
niche.breadths <- c(.3, .3, .3)

Here, we create a site by species abundance data.frame for three species at five sites.

my.J.t0 <- data.frame(
  a = c(10000, 1000, 500, 500, 100),
  b = c(10, 30, 10, 10, 40),
  c = c(100, 500, 500, 1000, 10000))

Step 2. Create the landscapes

Next, we will create the initial and “disturbed” landscapes that we will concatenate to create a simulation in which the environmental characteristics (Ef) becomes homogenized half way through the simulation.

# calculate site carrying capacities as the row totals from the site by species data.frame
my.JL.vals <- rowSums(my.J.t0)

# make the initial landscape, with variable environmental conditions across the five sites
my.landscape.1 <- make.landscape(
  site.coords = data.frame(
    x = c(1, 2, 3, 4, 5),
    y = c(1, 3, 1, 5, 2)),
  Ef = c(-.8, -.6, 0, .25, .9),
  m = 0.01,
  JL = my.JL.vals)
## [1] "success!!"
# make the "disturbed" landscape, in which the environmental filter is homogenized by assigning the same value to Ef across all sites.
my.landscape.2 <- make.landscape(
  site.coords = data.frame(
    x = c(1, 2, 3, 4, 5),
    y = c(1, 3, 1, 5, 2)),
  Ef = rep(.9, 5), # Note homogenization of Ef across sites
  m = 0.01,
  JL = my.JL.vals)
## [1] "success!!"

Step 3. Run the simulation

Here we use the metasim.disturbance wrapper to run MCSim simulations over the component landscapes concatenated in a landscape.list. Note that the time.interval.durations argument determines how many timesteps each landscape in the landscape.list is used.

sim.result <- metasim.disturbance(
  scenario.name = "niche_shift",
  landscape.list = list(my.landscape.1,my.landscape.2),
  time.interval.durations = c(10,10),
  trait.Ef = niche.positions,
  trait.Ef.sd = niche.breadths,
  J.t0.initial = my.J.t0,
  nu = 0, 
  W.r = 10,
  output.dir.path = "my_disturbance_sim_output_directory")

Step 4. Visualize the simulation outcomes

Here are some simple visualizations you can use for a simulation:

Plot the sites in two-dimensional space:

# view map of sites
my.landscape.1$site.coords %>%
  ggplot(aes(x,y,size = my.landscape.1$site.info$JL)) +
  geom_point()

Use the plot.coenoclines function in MCSim to visualize species’ habitat preferences and niche breadths relative to the position of the sites along the environmental gradient:

# Landscape 1 coenoclines
plot.coenoclines(landscape = my.landscape.1,
                 trait.Ef = niche.positions,
                 trait.Ef.sd = niche.breadths)

# Landscape 2 coenoclines -- note the shift of site positions on the x-axis
plot.coenoclines(landscape = my.landscape.2,
                 trait.Ef = niche.positions,
                 trait.Ef.sd = niche.breadths)

Use a dot plot to visualize how metacommunity structure changes over the course of the simulation:

# look at change in composition over time
plot.dot.plots(sim.result, timesteps = c(1:20))
## WARNING: this sim result includes a changing landscape, this plotting function is has not been optimized for changing landscapes and the plot will be based on the initial landscape configuration

Note that the homogenization of the environment (Ef) homogenizes the species pool after timestep 10.

Example 2. Change site carrying capacities (JL)

Step 1. Create the initial species pool

Note that we will use the same species pool as in the previous example.

Step 2. Create the landscapes

In this example, the initial and “disturbed” landscapes will have the same environmental characteristics (Ef), however, landscape 2 will have site-level carrying capacities reduced by three orders of magnitude.

set.seed(1234) #set random seed

# make a landscapes
my.landscape.1 <- make.landscape(
  site.coords = data.frame(
    x = c(1, 2, 3, 4, 5),
    y = c(1, 3, 1, 5, 2)),
  Ef = c(-.8, -.6, 0, .25, .9),
  m = 0.5,
  JL = rep(10000, 5))
## [1] "success!!"
my.landscape.2 <- make.landscape(
  site.coords = data.frame(
    x = c(1, 2, 3, 4, 5),
    y = c(1, 3, 1, 5, 2)),
  Ef = c(-.8, -.6, 0, .25, .9),
  m = 0.5,
  JL = rep(10, 5))
## [1] "success!!"

Step 3. Run the simulation

Again, use the metasim.disturbance wrapper to run MCSim simulations over the component landscapes concatenated in a landscape.list.

sim.result <- metasim.disturbance(
  scenario.name = "JL_reduction",
  landscape.list = list(my.landscape.1,my.landscape.2),
  time.interval.durations = c(10,10),
  trait.Ef = niche.positions,
  trait.Ef.sd = niche.breadths,
  J.t0.initial = my.J.t0,
  nu = 0,
  W.r = 0,
  output.dir.path = "my_disturbance_sim_output_directory")

Step 4. Visualize the simulation outcomes

Use the plot.coenoclines function in MCSim to visualize species’ habitat preferences and niche breadths relative to the position of the sites along the environmental gradient. Note that in this scenario, the ceonoclines are the same before and after the disturbance, because we are not manipulating Ef.

# Landscape 1 coenoclines
plot.coenoclines(landscape = my.landscape.1,
                 trait.Ef = niche.positions,
                 trait.Ef.sd = niche.breadths)

# Landscape 2 coenoclines -- note the shift of site positions on the x-axis
plot.coenoclines(landscape = my.landscape.2,
                 trait.Ef = niche.positions,
                 trait.Ef.sd = niche.breadths)

The dot plot shows that locally rare taxa are lost at each site after the “disturbance”, which was a reduction in site carrying capacity:

# look at change in composition over time
plot.dot.plots(sim.result, timesteps = c(1:20))
## WARNING: this sim result includes a changing landscape, this plotting function is has not been optimized for changing landscapes and the plot will be based on the initial landscape configuration