The package ATNr defines the differential equations and parametrisation of different versions of the Allometric Trophic Network (ATN) model. It is structured around a model object that contains a function implementing the ordinary differential equations (ODEs) of the model and various attributes defining the different parameters to run the ODEs.
Two different versions of the model are implemented:
Scaled version: Delmas et al. (2017)
Unscaled version incorporating nutrient dynamics: Schneider et al. (2016)
Unscaled version without nutrients: Binzer et al. (2016)
The version without nutrients from Delmas et
al. (2017) is scaled, meaning that the biological rates
controlling the growth rate of the species are normalised by the growth
rate of the smallest basal species. For more details on the three
models, see the specific vignette:
vignette(model_descriptions, package = "ATNr")
.
The definition of ATN is based on a model object (formally a S4 class in R). The model object is initialised specifying a fixed set of parameters: the number of species, the number of basal species, species body masses, a matrix defining the trophic interactions and, for the version including the nutrient dynamics, the number of nutrients. The first thing to do is therefore to create the corresponding R variables. While one can use an empirical food web for its analysis, it is also possible to generate synthetic food webs using the niche model from Williams and Martinez (2000) or using allometric scaling as defined in Schneider et al. (2016).
ATNr has two functions to generate synthetic food webs,
create_niche_model()
for the niche model (Williams and Martinez (2000)) and
create_Lmatrix()
for the allometric scaling model (Schneider et al. (2016)). The niche model
requires information on the number of species and connectance of the
desired food web:
library(ATNr)
set.seed(123)
<- 20 # number of species
n_species <- 0.3 # connectance
conn <- create_niche_model(n_species, conn)
fw # The number of basal species can be calculated:
<- sum(colSums(fw) == 0) n_basal
As the niche model does not rely on allometry, it is possible to
estimate species body masses from their trophic levels, which can be
calculated form the TroLev
function of the package. For
instance:
= TroLev(fw) #trophic levels
TL <- 1e-2 * 10 ^ (TL - 1) masses
The allometric scaling model generate links based on species body masses. Therefore, it requires as an input a vector containing the body mass of species, as well as a parameter informing on the number of basal species desired. It produces a so-called L matrix which formally quantifies the probability for a consumer to successfully attack and consumer an encountered resource:
<- 20
n_species <- 5
n_basal <- sort(10^runif(n_species, 2, 6)) #body mass of species
masses <- create_Lmatrix(masses, n_basal) L
This L matrix can then be transformed into a binary food web:
<- L
fw > 0] <- 1 fw[fw
More details about the generative models and the the usage precaution around them can be found in the section “The food web generative functions”
As soon as a food web is stored in a matrix, it is possible to create a model object that refers to the desired specific model
# initialisation of the model object. It is possible to create a ode corresponding to
# Schneider et al. 2016, Delmas et al. 2016 or Binzer et al. 2016:
# 1) Schneider et al. 2016
<- 3
n_nutrients <- create_model_Unscaled_nuts(n_species, n_basal, n_nutrients, masses, fw)
model_unscaled_nuts # 2) Delmas et al. 2016:
<- create_model_Scaled(n_species, n_basal, masses, fw)
model_scaled # 3) Binzer et al. 2016
<- create_model_Unscaled(n_species, n_basal, masses, fw) model_unscaled
Once created, it is possible to access to the methods and attributes of the object to initialise or update them:
# updating the hill coefficient of consumers in the Unscaled_nuts model:
$q <- rep(1.4, model_unscaled_nuts$nb_s - model_unscaled_nuts$nb_s)
model_unscaled_nuts# Changing the assimilation efficiencies of all species to 0.5 in the Scaled model:
$e = rep(0.5, model_scaled$nb_s)
model_scaled# print the different fields that can be updated and their values:
# str(model_unscaled_nuts)
It is important to keep in mind that some rules apply here:
The order of the species in the different fields must be
consistent: the first species in the $BM
object corresponds
to the first species in the $fw
object and in the
$e
object.
The objects that are specific to a species type (i.e. basal
species or consumers) are dimensioned accordingly: the handling time
($h
) sets the handling time of consumers on resources.
Therefore, the h
matrix has a number of rows equal to the
number of species and a number of columns equal to the number of
consumers (as non consumer species do not have a handling time by
definition). In that case, the first row correspond to the first species
and the first column to the first consumer.
The object describing the interactions between plants and
nutrients ($K
or $V
) are matrices for which
the number of rows equals to the number of nutrients and a number of
columns matches the number of basal species (this point is specific to
the Schneider model which is the only one including an explicit dynamics
of the nutrient pool).
To run the population dynamics, all the parameters must be defined. It is possible to automatically load a by default parametrisation using the dedicated functions:
# for a model created by create_model_Unscaled_nuts():
<- initialise_default_Unscaled_nuts(model_unscaled_nuts, L)
model_unscaled_nuts # for a model created by create_model_Scaled():
<- initialise_default_Scaled(model_scaled)
model_scaled # for a model created by create_model_Unscaled():
<- initialise_default_Unscaled(model_unscaled) model_unscaled
Importantly, for the unscaled with nutrients model of Schneider et al. (2016), the calculation of
consumption rates rely on the L matrix created above or, in case of
empirical networks, by a matrix that defines the probability of a
consumer to successfully attack and consume an encountered prey. The
default initialisation of the Unsacled_nuts
and
Unscaled
models can also include temperature effects (20° C
by default).
Once all the parameters are properly defined, the ODEs can be
integrated by any solver. We present here a solution based on
lsoda
from library deSolve
(Soetaert, Petzoldt, and Setzer (2010)), but
other solutions exist (sundialr
is also a possibility). The
package propose a direct wrapper to lsoda
with the function
lsoda_wrapper
:
<-runif(n_species, 2, 3) # starting biomasses
biomasses <- append(runif(3, 20, 30), biomasses) # nutrient concentration
biomasses # defining the desired integration time
<- seq(0, 1500, 5)
times <- lsoda_wrapper(times, biomasses, model_unscaled_nuts) sol
To have more control of the integration, it is however possible to
not use the wrapper proposed in the package and directly work with the
lsoda
function. Here is an example:
# running simulations for the Schneider model
$initialisations()
model_unscaled_nuts<- deSolve::lsoda(
sol
biomasses,
times,function(t, y, params) {
return(list(params$ODE(y, t)))
},
model_unscaled_nuts )
Not that the call of
model_unscaled_nuts$initialisation()
is important here as
it pre-computes some variables to optimise code execution. This function
is normally internally called by lsoda_wrapper
. In case of
an integration that does not rely on this wrapper function, the call to
$initialisation()
is needed for ALL the model types.
The package also contains a simple function to plot the time series obtained: plot_odeweb. The colours only differentiate the species using their ranks in the food web matrix (from blue to red).
par(mar = c(4, 4, 1, 1))
plot_odeweb(sol, model_unscaled_nuts$nb_s)
It is possible to create a model object using empirical food webs, however synthetic ones can be valuable tools to explore different theoretical questions. To allow this possibility, two different models are available in the package: the niche model (Williams and Martinez (2000)) or the allometric scaling model (Schneider et al. (2016)). Thereafter, we use the following function to visualise the adjacency matrices (where rows correspond to resources and columns to consumers) of the food webs:
# function to plot the fw
<- function(mat, title = NULL) {
show_fw par(mar = c(.5, .5, 2, .5))
<- nrow(mat)
S <- mat[nrow(mat):1, ]
mat <- t(mat)
mat image(mat, col = c("goldenrod", "steelblue"),
frame = FALSE, axes = FALSE)
title(title)
grid(nx = S, ny = S, lty = 1, col = adjustcolor("grey20", alpha.f = .1))
}
The niche model orders species based on their trophic niche, randomly sampled from a uniform distribution. For each species \(i\), a diet range (\(r_i\)) is then drawn from a Beta distribution and a diet center \(c_i\) from a uniform distribution. For each species \(i\), all species that have trophic niche within the interval \([c_i - r_i / 2, c_i + r_i / 2]\) are considered to be prey of species \(i\). In this package, we followed the modification to the niche model of Williams and Martinez (2000) as specified in Allesina, Alonso, and Pascual (2008).
Generating a food web from the niche model is made by a simple call to the corresponding functions:
<- 50 # number of species
S <- 0.2 # connectance
C <- create_niche_model(S, C) fw
The function ensure that the food web returned are not composed of disconnected sub networks (i.i several connected components).
The allometric scaling model assumes an optimal consumer/resource body mass ratio (Ropt, default = 100) for attack rates, i.e. the probability that when a consumer encounter a species it will predate on it. In particular, each attack rate is calculated using a Ricker function:
\[ a_{ij} = \left( \frac{m_i}{m_j \cdot Ropt} \cdot e^{(1 - \frac{m_i}{m_j \cdot Ropt})} \right) ^\gamma \]
where \(m_i\) is the body mass of species \(i\) and \(\gamma\) sets the width of the trophic niche.
Generating a food web with the allometric scaling model necessitate
few more steps. The trophic niche of species is defined by a body mass
interval and is quantified (see fig. 2 and 3 from Schneider et al.,
2016). This quantified version actually return the probabilities of a
successful attack event to occur when a consumer encounter a prey. These
probabilities are estimated with a Ricker function of 4 parameters: the
body masses of the resource and of the consumer, the optimal
predator-prey body mass ratio Ropt
and the width of the
trophic niche gamma
. A threshold (th
) filters
out links with very low probabilities of attack success. The
probabilities are stored in a matrix obtained from:
# number of species and body masses
<- 20
n_species <- 5
n_basal # body mass of species. Here we assume two specific rules for basal and non basal species
<- c(sort(10^runif(n_basal, 1, 3)), sort(10^runif(n_species - n_basal, 2, 6)))
masses <- create_Lmatrix(masses, n_basal, Ropt = 100, gamma = 2, th = 0.01) L
Then, a food web is a binary version of the L matrix that can be stored either using booleans (FALSE/TRUE) or numeric values (0/1):
# boolean version
<- L > 0
fw # 0/1 version:
<- L
fw > 0] <- 1
fw[fw show_fw(fw, title = "L-matrix model food web")
ATNr makes it relatively easy to vary one parameter to assess its effect on the population dynamics. For example, we can study how changes in temperatures from 15 to 30 Celsius degrees affects the number species to go extinct.
set.seed(12)
# 1) define number of species, their body masses, and the structure of the
# community
<- 50
n_species <- 20
n_basal <- 2
n_nut # body mass of species
<- 10 ^ c(sort(runif(n_basal, 1, 3)),
masses sort(runif(n_species - n_basal, 2, 9)))
# 2) create the food web
# create the L matrix
<- create_Lmatrix(masses, n_basal, Ropt = 50, gamma = 2, th = 0.01)
L # create the 0/1 version of the food web
<- L
fw > 0] <- 1
fw[fw # 3) create the model
<- create_model_Unscaled_nuts(n_species, n_basal, n_nut, masses, fw)
model # 4) define the temperature gradient and initial conditions
<- seq(4, 22, by = 2)
temperatures <- rep(NA, length(temperatures))
extinctions # defining biomasses
<- runif(n_species + n_nut, 2, 3)
biomasses # 5) define the desired integration time.
<- seq(0, 100000, 100)
times # 6) and loop over temperature to run the population dynamics
<- 0
i for (t in temperatures){
# initialise the model parameters for the specific temperature
# Here, no key parameters (numbers of species or species' body masses) are modified
# Therefore, it is not needed to create a new model object
# TO reinitialise the different parameters is enough
<- initialise_default_Unscaled_nuts(model, L, temperature = t)
model # updating the value of q, same for all consumers
$q = rep(1.4, n_species - n_basal)
model$S <- rep(10, n_nut)
model# running simulations for the Schneider model:
<- lsoda_wrapper(times, biomasses, model, verbose = FALSE)
sol # retrieve the number of species that went extinct before the end of the
# simulation excluding here the 3 first columns: first is time, 2nd and 3rd
# are nutrients
<- i + 1
i <- sum(sol[nrow(sol), 4:ncol(sol)] < 1e-6)
extinctions[i] }
plot(temperatures, extinctions,
pch = 20, cex = 0.5, ylim = c(0,50), frame = FALSE,
ylab = "Number of Extinctions", xlab = "Temperature (°C)")
lines(temperatures, extinctions, col = 'blue')
Predator-prey body mass ratio and environment temperature have been shown to affect persistence of species in local communities, e.g. Binzer et al. (2016). Here, we use the ATNr model (name here) to replicate the results from Binzer et al. (2016). In particular, we compute the fraction of species species that persist for predator-prey body mass ratio values in \(\left[ 10^{-1}, 10^4 \right]\) and temperature values in \(\{0, 40\}\) °C.
First, we create a food web with 30 species and initialize within a for loop the model with a given value of body mass ratio and temperature. Species persistence is calculate as the fraction of species that are not extinct at the end of the simulations.
# set.seed(142)
# number of species
<- 30
S
# vector containing the predator prey body mass ratios to test
<- 10 ^ seq(-1, 4, by = .5)
scaling
# vectors to store the results
<- c()
persistence0 <- c()
persistence40
# create the studied food web
<- create_niche_model(S = S, C = 0.1)
fw # calculating trophic levels
= TroLev(fw)
TL <- runif(S, 2, 3)
biomasses
# run a loop over the different pred-prey body mass ratios
for (scal in scaling) {
# update species body masses following the specific body mass ratio
<- 0.01 * scal ^ (TL - 1)
masses
# create the models with parameters corresponding to 0 and 40 degrees Celcius
<- create_model_Unscaled(nb_s = S,
mod0 nb_b = sum(colSums(fw) == 0),
BM = masses,
fw = fw)
<- initialise_default_Unscaled(mod0, temperature = 0)
mod0 $c <- rep(0, mod0$nb_s - mod0$nb_b)
mod0$alpha <- diag(mod0$nb_b)
mod0
<- create_model_Unscaled(nb_s = S,
mod40 nb_b = sum(colSums(fw) == 0),
BM = masses,
fw = fw)
<- initialise_default_Unscaled(mod40, temperature = 40)
mod40 $c <- rep(0, mod40$nb_s - mod40$nb_b)
mod40$alpha <- diag(mod40$nb_b)
mod40
<- seq(1, 1e9, by = 1e7)
times
# run the model corresponding to the 0 degree conditions
<- lsoda_wrapper(times, biomasses, mod0, verbose = FALSE)
sol <- append(persistence0, sum(sol[nrow(sol), -1] > mod0$ext) / S)
persistence0 # run the model corresponding to the 40 degrees conditions
<- lsoda_wrapper(times, biomasses, mod40, verbose = FALSE)
sol <- append(persistence40, sum(sol[nrow(sol), -1] > mod40$ext) / S)
persistence40 }
Similarly to Binzer et al. (2016), species persistence increases with increasing values of predator-prey body mass ratios, but temperature effect differs depending n this ratio: when predator prey body mass ratio is low, high temperature lead to more persistence while increasing predator prey body mass ratio tend to reduce the effects of temperature.
plot(log10(scaling), persistence40,
xlab = expression("Body mass ratio between TL"[i + 1]* " and TL"[i]),
ylab = "Persistence",
ylim = c(0, 1),
frame = FALSE, axes = FALSE, type = 'l', col = "red")
lines(log10(scaling), persistence0, col = "blue")
axis(2, at = seq(0, 1, by = .1), labels = seq(0, 1, by = .1))
axis(1, at = seq(-1, 4, by = 1), labels = 10 ^ seq(-1, 4, by = 1))
legend(0.1, 0.9, legend = c("40 \u00B0C", "0 \u00B0C"), fill = c("red", "blue"))
The paradox of enrichment states that by increasing the carrying capacity of basal species may destabilize the population dynamics (Rosenzweig (1971)). Here, we show how this can be studied with the ATNr; we used the model from Delmas et al. (2017), but similar results can be obtained using the other two models in the package.
First, we create a food web with 10 species and initialize the model
set.seed(1234)
<- 10
S <- NULL
fw <- NULL
TL <- create_niche_model(S, C = .15)
fw <- TroLev(fw)
TL
<- 0.01 * 100 ^ (TL - 1) #body mass of species
masses
<- create_model_Scaled(nb_s = S,
mod nb_b = sum(colSums(fw) == 0),
BM = masses,
fw = fw)
<- initialise_default_Scaled(mod)
mod <- seq(0, 300, by = 2)
times <- runif(S, 2, 3) # starting biomasses biomasses
Then, we solve the system specifying the carrying capacity of basal
species equal to one (mod$K <- 1
) and then increased
this to ten (mod$K <- 10
)
$K <- 1
mod<- lsoda_wrapper(times, biomasses, mod, verbose = FALSE)
sol1 $K <- 10
mod<- lsoda_wrapper(times, biomasses, mod, verbose = FALSE) sol10
As shown in the plot below, for K = 1 the system reaches a stable equilibrium, whereas when we increase the carrying capacity (K = 10) the system departs from this stable equilibrium and periodic oscillations appear.
par(mfrow = c(2, 1))
plot_odeweb(sol1, S)
title("Carrying capacity = 1")
plot_odeweb(sol10, S)
title("Carrying capacity = 10")
The building block of this package are the C++ classes to solve ODEs for the ATN model. Model parameters are stored in such classes and can be changed only by addressing them within the respective objects. For instance:
set.seed(1234)
<- 20
nb_s <- 5
nb_b <- 2
nb_n <- sort(10 ^ runif(nb_s, 2, 6)) #body mass of species
masses = runif(nb_s + nb_n, 2, 3)
biomasses <- create_Lmatrix(masses, nb_b, Ropt = 50)
L 1:nb_b] <- 0
L[, <- L
fw > 0] <- 1
fw[fw <- create_model_Unscaled_nuts(nb_s, nb_b, nb_n, masses, fw)
model_unscaled_nuts <- initialise_default_Unscaled_nuts(model_unscaled_nuts, L)
model_unscaled_nuts <- 30 #this does not change the model parameter
nb_s $nb_s #this is the model parameter
model_unscaled_nuts#> [1] 20
Changing parameters that are used in the create_model
functions without creating a new model object is almost always a bad
idea. Those parameters are important structural parameters, and changing
one of them implies changes in most of the other variables contained in
the model object. For instance, the example above, changing the number
of species in the model object will lead to inconsistencies in the
different variables: the dimensions of objects storing attack rates,
body masses and so on won’t match the updated number of species. Some
basic checks are made before starting the integration in the
lsoda_wrapper
function, based on the
run_checks
procedure also available in the package.
<- seq(0, 15000, 150)
times $nb_s = 40
model_unscaled_nuts# this will return an error :
# sol <- lsoda_wrapper(times, biomasses, model_schneider)
However, some modification can remain undetected. For instance, modifying species’ body masses only won’t raise any errors. However, a change in species body mass should be associated to a change in all the associated biological rates. The following code won’t raise any errors, but will produce results relying on a model with an incoherent set of parameters and therefore wrong result:
set.seed(1234)
<- 20
nb_s <- 5
nb_b <- 2
nb_n <- sort(10 ^ runif(nb_s, 2, 6)) #body mass of species
masses = runif(nb_s + nb_n, 2, 3)
biomasses <- create_Lmatrix(masses, nb_b, Ropt = 50)
L 1:nb_b] <- 0
L[, <- L
fw > 0] <- 1
fw[fw <- create_model_Unscaled_nuts(nb_s, nb_b, nb_n, masses, fw)
model_unscaled_nuts <- initialise_default_Unscaled_nuts(model_unscaled_nuts, L)
model_unscaled_nuts $BM <- sqrt(model_unscaled_nuts$BM) # we change body masses within the model
model_unscaled_nuts<- lsoda_wrapper(seq(1, 5000, 50), biomasses, model_unscaled_nuts)
sol par(mar = c(4, 4, 1, 1))
plot_odeweb(sol, model_unscaled_nuts$nb_s)
In general, each time one of these key parameters is modified
(nb_s
, nb_b
, nb_n
for the
Schneider model, BM
, fw
), it is strongly
recommended to create a new model objects with the updated
parameters:
<- 30
nb_s <- 2
nb_n <- sort(10 ^ runif(nb_s, 2, 6)) #body mass of species
masses <- runif(nb_s + nb_n, 2, 3)
biomasses <- create_Lmatrix(masses, nb_b, Ropt = 50)
L 1:nb_b] <- 0
L[, <- L
fw > 0] <- 1
fw[fw # create a new object:
<- create_model_Unscaled_nuts(nb_s, nb_b, nb_n, masses, fw)
model_unscaled_nuts <- initialise_default_Unscaled_nuts(model_unscaled_nuts, L)
model_unscaled_nuts # safely run the integration:
<- lsoda_wrapper(times, biomasses, model_unscaled_nuts) sol
Specifically to the Schneider model, changing the ‘Lmatrix’ requires
to update the feeding rate ($b
).
Changing the dimensions of a vector or matrix object somehow implies
a change in the number of species (see section above). For instance,
decreasing the length of the assimilation efficiencies vector
$e
should imply a consistent update of all the other
parameters (handling times, attack rates, body masses, etc depending on
the model). The function remove_species
is made to properly
remove species from model objects without having to manually regenerate
all parameters.
As built on Rcpp, the different models are only pointers to C++ objects. It means that this script:
<- 30
nb_s <- 2
nb_n <- sort(10 ^ runif(nb_s, 2, 6)) #body mass of species
masses <- runif(nb_s + nb_n, 2, 3)
biomasses <- create_Lmatrix(masses, nb_b, Ropt = 50)
L 1:nb_b] <- 0
L[, <- L
fw > 0] <- 1
fw[fw # create a new object:
<- create_model_Unscaled_nuts(nb_s, nb_b, nb_n, masses, fw)
model_1 <- initialise_default_Unscaled_nuts(model_1, L)
model_1
# trying to create a new model that is similar to model_1
= model_1 model_2
will not create a new model object. Formally, it creates a new
pointer to the same address, which means that model_1
and
model_2
are in reality the same variable (shallow copy).
Therefore, modifying one modifies the other:
$q = 1.8
model_1# this also updated the value in model_2:
$q
model_2#> [,1]
#> [1,] 1.8
Therefore, to create a new model object based on another one, it is
important to formally create one (either with one of the
create_model_
function, or using new
). More
information on copying variables using pointers here: https://stackoverflow.com/questions/184710/what-is-the-difference-between-a-deep-copy-and-a-shallow-copy
Modifying a R variable inside a *apply
function in does
not modify it:
.3 = function(x, useless) {
plus= x+3
y = useless + 1
useless return(y)
}
= 4:10
useless = useless
useless2 = sapply(1:5, plus.3, useless)
x # the useless variable was not modified:
== useless2
useless #> [1] TRUE TRUE TRUE TRUE TRUE TRUE TRUE
However, this is not the case anymore with a model object. If we consider a model object:
<- 20
n_species <- 5
n_basal = n_species - n_basal
n_cons <- 2
n_nut <- 10 ^ c(sort(runif(n_basal, 0, 3)),
masses sort(runif(n_species - n_basal, 2, 5)))
<- create_Lmatrix(masses, n_basal, Ropt = 100, gamma = 2, th = 0.01)
L #> Warning in create_Lmatrix(masses, n_basal, Ropt = 100, gamma = 2, th = 0.01):
#> Presence of consumer without prey.
<- L
fw > 0] <- 1
fw[fw <- create_model_Unscaled_nuts(n_species, n_basal, n_nut, masses, fw)
model <- initialise_default_Unscaled_nuts(model, L, temperature = 20) model
and a function that takes this model object as an argument, setting the b matrix to 0:
# a function that sets all elements of model$b to 0
<- function(x, model){
a.fun $b = model$b*0
modelreturn(x+1)
}
then, we can see that the global model object is indeed modified when the function is called by *apply:
= c(1,2)
x sum(model$b)
= lapply(x, a.fun, model)
y sum(model$b)
This behaviour is still due to the fact that in a *apply function the model is shallow-copied and each iteration points in fact to the same object in memory. However, this behaviour is not present when using a parallel version of an apply function, as in parallel computations the object is automatically deep-copied and passed to each task separately:
library(parallel)
sum(model$b)
<- initialise_default_Unscaled_nuts(model, L, temperature = 20)
model = mclapply(x, a.fun, model = model, mc.cores=5)
y sum(model$b)