Skip to contents

Generate time-series with The Self-Organised Instability (SOI) model. Implements a K-leap method for accelerating stochastic simulation.

Usage

simulateSOI(
  n_species,
  x0 = NULL,
  names_species = NULL,
  carrying_capacity = 1000,
  A = NULL,
  k_events = 5,
  t_end = 1000,
  metacommunity_probability = runif(n_species, min = 0.1, max = 0.8),
  death_rates = runif(n_species, min = 0.01, max = 0.08),
  norm = FALSE
)

Arguments

n_species

Integer: number of species

x0

a vector of initial community abundances If (default: x0 = NULL), based on migration rates

names_species

Character: names of species. If NULL, paste0("sp", seq_len(n_species)) is used. (default: names_species = NULL)

carrying_capacity

integer community size, number of available sites (individuals)

A

matrix: interaction matrix defining the positive and negative interactions between n_species. If NULL, powerlawA(n_species) is used. (default: A = NULL)

k_events

integer number of transition events that are allowed to take place during one leap. (default: k_events = 5). Higher values reduce runtime, but also accuracy of the simulation.

t_end

Numeric: the end time of the simulation, defining the modeled time length of the community. (default: t_end = 1000)

metacommunity_probability

Numeric: Normalized probability distribution of the likelihood that species from the metacommunity can enter the community during the simulation. By default, runif(n_species, min = 0.1, max = 0.8) is used. (default: metacommunity_probability = runif(n_species, min = 0.1, max = 0.8))

death_rates

Numeric: death rates of each species. By default, runif(n_species, min = 0.01, max = 0.08) is used. (default: death_rates = runif(n_species, min = 0.01, max = 0.08))

norm

logical scalar indicating whether the time series should be returned with the abundances as proportions (norm = TRUE) or the raw counts (default: norm = FALSE)

Value

simulateSOI returns a TreeSummarizedExperiment class object

Examples

# Generate interaction matrix
A <- miaSim::powerlawA(10, alpha = 1.2)
# Simulate data from the SOI model
tse <- simulateSOI(
    n_species = 10, carrying_capacity = 1000, A = A,
    k_events = 5, x0 = NULL, t_end = 150, norm = TRUE
)