In dswdejonge/fwstability: Infer Food-Web Stability
knitr::opts_chunk$set(
collapse = TRUE,
comment = "#>",
eval = FALSE
)
The fwstability package is developed to infer food-web stability from energy-flux models. Basically, an energy-flux food-web model shows the transfer of energy or materials from one species to another when they eat each other. These models are steady-state and mass-balanced. Mathematical methods have been developed to transform these steady-state models to the dynamic realm in order to infer stability. This package implements these mathematical methods, so that you can easily infer stability from your energy-flux food-web model.
We are currently preparing a paper introducing this package and an in-depth explanation of the science behind energy-flux food-web stability. Access to the manuscript is available upon reasonable request until publication:
Daniëlle S.W. de Jonge, Peter C. de Ruiter, Johan van de Koppel, Dick van Oevelen. Inferring stability of energy-flux food-web models. In preparation.
This vignette provides you a quick start, and a more detailed explanation about the methods. Reading the vignette is highly recommended so that you can interpret your results correctly.
1. Quick start
1.1. Install and load
To install the development version of the package:
# if you don't have the devtools package yet:
install.packages("devtools")
# Install the fwstability package, development version:
devtools::install_github("dswdejonge/fwstability", build_vignettes = TRUE)
There is also the fwmodels package available from github, that contains several published energy-flux food-web models that are used as examples in this vignette.
# To use the example food-webs, install the fwmodels package as well
devtools::install_github("dswdejonge/fwmodels")
If the installation was successful, we load the fwstability package (and the fwmodels package for examples).
library(fwstability)
library(fwmodels)
1.2. Prepare data
First, you must prepare your energy-flux food-web model to be used by the package. See section 2 for the requirements in the physiology of the food-web model.
The flowmatrix always has source compartments in rows, and sink compartments in columns. In other words, the value in each cell is the flow moving from the row compartment toward to the column compartment.
A case that perhaps requires some extra explanation is the 'frac_matrix'. This matrix has the same size as the flow matrix and indicates the fraction of each flow that comprises defecation. Take the example where there is a material flow from Species A to Detritus of 2 units due to defecation, and another of 8 units due to mortality of the population. The total flow between Species A and Detritus is 10 units (this is the value present in the flow matrix). Of those total 10 units, 2 units are caused by defecation. Therefore, the frac_matrix will have the value 2/10 = 0.2. If there is only a flux between Species A and Detritus due to defecation, the fraction will be 1. For all flows which are not comprised of defecation at all, for example when Species A eats Species B, the fraction in the frac_matrix will be 0.
Dead compartment are all compartments that are not alive, e.g. detritus, nutrients, etc. Information about these compartments should be included in a model element called 'dead'.
External compartments are compartments from which material is imported into and exported from the model, and which can be disregarded for inferring interaction strengths and stability. If present, their names can be included in a model element called 'externals'.
A hypothetical example of what your data can look like:
# 1. Biomasses for each food web compartment
# (e.g. abundance per unit area).
compartments <- c("detritus", "speciesA", "speciesB")
biomasses <- c(30, 20, 10) # e.g. grams per square meter
names(biomasses) <- compartments
# 2. Quantified fluxes between compartments
# (e.g. abundance per unit area per unit time).
detritus_to_speciesA <- 10
speciesA_to_detritus <- 3
speciesA_to_speciesB <- 7
speciesB_to_detritus <- 7
flowmatrix <- matrix(data = 0,
nrow = length(compartments),
ncol = length(compartments))
row.names(flowmatrix) <- compartments
colnames(flowmatrix) <- compartments
flowmatrix["detritus", "speciesA"] <- detritus_to_speciesA
flowmatrix["speciesA", "detritus"] <- speciesA_to_detritus
flowmatrix["speciesA", "speciesB"] <- speciesA_to_speciesB
flowmatrix["speciesB", "detritus"] <- speciesB_to_detritus
# 3. Energy conversion efficiencies for all living compartments
# (assimilation and production efficiencies, fraction 0-1).
assimilation_eff <- c(NA, 0.8, 0.7)
growth_eff <- c(NA, 0.6, 0.4)
names(assimilation_eff) <- compartments
names(growth_eff) <- compartments
# 4. Information on dead compartments, if any
# 4.a. The names of all dead compartments
dead_names <- "detritus"
# 4.b. The fraction of each flow which comprises defecation
# frac_matrix needs to be the same size and the flow matrix
frac_matrix <- matrix(data = 0,
nrow = length(compartments),
ncol = length(compartments))
row.names(frac_matrix) <- compartments
colnames(frac_matrix) <- compartments
# Defecation is the part of ingested food (food source -> species)
# which is NOT assimilated (1 - assimilation efficiency)
# For more explanation about physiology, see Section 2.
defecation_by_speciesA <- flowmatrix["detritus", "speciesA"] *
(1-assimilation_eff["speciesA"])
defecation_by_speciesB <- flowmatrix["speciesA", "speciesB"] *
(1-assimilation_eff["speciesB"])
# You can simply divide the defecation flow by the total flow
# between compartments to obtain the value for the frac_matrix
frac_matrix["speciesA", "detritus"] <- defecation_by_speciesA /
flowmatrix["speciesA", "detritus"]
frac_matrix["speciesB", "detritus"] <- defecation_by_speciesB /
flowmatrix["speciesB", "detritus"]
Combine all information into a named list as follows:
mymodel = list(
type = "EF", # EF stands for energy flux
FM = flowmatrix,
BM = biomasses,
AE = assimilation_eff,
GE = growth_eff,
dead = list(
names = dead_names,
frac = frac_matrix))
1.3. Jacobian matrix
Now we can estimate a Jacobian matrix from the energy-flux model (see section 3 and our article for more detail).
The Jacobian matrix shows the per unit biomass effect of species i (in rows) on the population growth rate of species j (in columns) per unit time. In other words, the Jacobian matrix shows you the interaction strengths between species. Thus, the unit of the JM is 'per unit time'.
The diagonal values are intraspecific effects, thus the density-dependent effects a compartment has on itself (works for both species and detritus).
We use an example food web from the fwmodels package. We recommend using the default diagonal, i.e. calculated from the model (example JM1 below):
library(fwstability)
library(fwmodels)
model <- LovinkhoeveCP
# Documentation for the getJacobian function
?getJacobian
# Create a Jacobian matrix, with a diagonal calculated from the model (default).
# The diagonal values are the upper-limit of self-dampening effects in the system.
JM1 <- getJacobian(model, diagonal = "model")
# Create a Jacobian matrix, with an all-zero diagonal.
JM2 <- getJacobian(model, diagonal = 0)
# Create a Jacobian matrix, with -1 on the diagonal.
JM3 <- getJacobian(model, diagonal = -1)
# Create a Jacobian matrix, with a numeric vector as diagonal
diagonal <- runif(dim(model$FM)[1], min = -1, max = 0)
JM4 <- getJacobian(model, diagonal = diagonal)
# The preview of the first 5 rows and columns of the Jacobian matrix:
print(JM1[1:5, 1:5])
1.4. Infer stability
Now we can infer stability from the Jacobian matrix. All methods are based on the concept that a Jacobian matrix is only stable if all real parts of the eigenvalues are negative.
If the system is stable the system will asymptotically return to equilibrium after a small pulse perturbation. If the system is unstable the system will not return to equilibrium after a small pulse perturbation, but move away from the equilibrium state.
Method s1 (below) finds the dominant eigenvalue of the Jacobian matrix. If positive, the system is NOT stable. If negative the system is stable. In a stable system the dominant eigenvalue represent the return rate to equilibrium. Method s1 is NOT suitable to compare relative stability between different systems if these systems operate on different time-scales.
Method s2 (below) also finds the dominant eigenvalue, but from a normalized Jacobian matrix with an all-zero diagonal (except for detritus). A normalized JM is unitless, and can therefore be used to compare relative stability between systems. The value of s2 is expected to lie between 0 and 1, and represents the minimum fraction of self-dampening necessary to stabilize the system. A lower value of s2 means a more stable system. If s2 > 1, more than available self-dampening is necessary for a stable system. Hence, the system is inherently unstable. We recommend this method.
Method s3 (below) is the scalar method, used in early publication (e.g. de Ruiter et al. 1997). The value of s3 can be interpreted the same as s2 (Neutel & Thorne 2014), however, the resolution of s2 is probably better than the resolution of s3, and s2 is quicker to calculate.
Both method s2 and s3 assume the JM without self-dampening effects on the diagonal is unstable.
# --- Method s1 ---
# Calculate stability with the eigenvalue method (default),
# from a Jacobian matrix with unit per time.
s1 <- getStability(JM1)
# If s1 is positive, the system is NOT stable.
# If s1 is negative, the system is stable.
# If s1 is stable, the asymptotic return rate is:
return_rate <- abs(s1) # per unit time (e.g. per day)
# If s1 is stable, the time for the system to return
# to about 37% of its original state is:
return_time <- 1/return_rate # unit of time (e.g. days or years)
# BEWARE: s1 and the derived return rates and return times
# are NOT suitable to compare relative stability among different systems.
# --- Method s2 ---
# Calculate stability with the eigenvalue method (default),
# from a normalized Jacobian matrix with an all-zero diagonal.
# Normalize JM (JM must have quantified diagonal):
JMnorm <- normalizeJacobian(JM1, dead_names = model$dead$names)
# Find dominant eigenvalue:
s2 <- getStability(JMnorm)
# s2 is expected to lie between 0 and 1
# s2 represents the minimum fraction of available self-dampening
# necessary to stabilize the system.
# --- Method 3 ---
# Calculate stability with the scalar method.
# Natural mortality rates (MR in per unit time) and the names of all
# dead (i.e. detritus) compartments are necessary.
s3 <- getStability(JM2, method = "scalar",
MR = model$MR, dead_names = model$dead$names)
# If the diagonal is calculated from the model,
# the diagonal can be used as MR
s3 <- getStability(JM1, method = "scalar",
MR = abs(diag(JM1)), dead_names = model$dead$names)
2. Preparing a food web model for stability analysis
Your energy-flux model should always stick to the physiological concepts explained in this diagram:
2.1. A linear inverse model
In a linear inverse model (LIM) all processes in a food web are captured by a set of linear equations. These linear equations can be either equalities (the outcome of an equation should equal another value) or inequalities (the outcome of an equation should be equal or greater than another value). The equations can parameterized with known data, but often not all processes are studied well enough to quantify. The unknowns in the equations can then be found with linear programming, thereby fully quantifying your model. A simple example of a linear optimization problem can be found on \link{Wikipedia}[https://en.wikipedia.org/wiki/Linear_programming#Example].
A linear inverse food web model (LIM) can be prepared with the R-package LIM (https://CRAN.R-project.org/package=LIM). More information on using LIM for food web quantification can be found in:
van Oevelen, D., van den Meersche, K., Meysman, F.J.R., Soetaert, K., Middelburg, J.J., Vézina, A.F., 2010. Quantifying food web flows using linear inverse models. Ecosystems 13, 32–45. https://doi.org/10.1007/s10021-009-9297-6.
A LIM model produced by the LIM package can be used directly by the fwstability package if it is formatted in a certain way, as explained below. If your LIM is not formatted according to the requirements of the fwstability package, you will need to extract certain data from the model yourself before using the package. For the latter, please refer to 2.2. Energy-flux model.
In order to obtain a Jacobian matrix, certain physiological processes must be calculated from your LIM. The fwstability can automatically extract this data if you labelled them with a unique tag. In the model you specify:
- which compartments are dead compartments (e.g. detritus and carrion, default tag is "dead")
- which flows specifically represent defecation (default tag is "def")
- which flows specifically represent mortality (default tag is "mort")
- which variables represent total assimilation (default tag is "ass", beware to not tag assimilation efficiency)
- which variables represent total growth (secondary production, default tag is "growth", beware to not tag growth efficiency)
You can use the default tags, but it is also possible to use your own tags. In the latter case you will need to provide the package with the name of your chosen tags in order to extract the correct information. The package will warn you if it cannot find the tag in your model.
Rules and tips for using tags:
- The official component name is the name given in the 'stocks' section. Sometimes, stock values are already given in the 'parameter' section and then assigned to the defined stocks in the 'stocks' section. Be aware not to tag the compartments in the 'parameter' section as these are not the official stock names used throughout the rest of the model analysis.
- Flows that represent defecation and mortality must be named in the 'flows' section and include the defecation or mortality tag. If defecation and mortality are combined in the same flow, they should still be declared separately (so you get the same flow twice, but with a different name).
- The variables representing total defecation, mortality, assimilation, and growth should be named by combining the component name and tag only. The code subtracts the tag from the name of the variable and uses it to place the data with the right compartment. So, a variable containing total rabbit mortality should be named "rabbitMort". If you would name it "thisisRabbitMort" it is not recognized due to the superfluous text. It is possible to specify multiple mortality flows, e.g. if the carcass of the rabbit consists of both high-energy and low-energy material which you want to separate into two separate carcass compartments. In this case it is necessary to define rabbitMort as the sum of these flows in the variable section of the model.
- Tags are not case sensitive. So, RABBITmort and RaBbItMoRt are both valid.
- The tag may occur anywhere in the string (so, it does not have to be added specifically at the beginning or ending of the name of the compartment/flow/variable). Ensure that the rest of names do not accidentally contain a tag! For example, if you include the Indian wild ass (Equus hemionus khur) in your foodweb, a variable representing its respiration called respirationAss will be interpreted as the assimilation variable of an unknown compartment 'respiration' because it contains the default assimilation tag "ass".
- Dead compartments, defecation, and mortality can be omitted from your model if you don't want to include them (for example make them all export to an external compartment). However, including assimilation and growth is required for all non-dead compartments. If you have no growth or assimilation data (for example for bacteria) you can assume the full ingestion is used for assimilation and/or growth and define your variable accordingly.
General tips for setting up a LIM:
- It is recommended to be consistent in your naming of flows and variables so later you can easily extract data. For naming flows for example: sourceProcessSink.
- In order to prevent certain processes to be quantified as zero, include minimum constraints. This is especially important if you later want to normalize the Jacobian matrix, because this process does not allow zeroes on the diagonal. For example, if you don't set a minimum constraint for non-predatory mortality, the parsimonious solution might set this flow to zero. This is not realistic and it the occurrence of these zeroes may hamper your future calculations.
- The inclusion of biomasses in the LIM is not required to be able to solve the model. However, biomasses are required for all internal compartments in order to calculate the Jacobian matrix. If there are compartments that do not have biomasses in the LIM, they must either be quantified anyway or excluded from the getJacobian() function by including them as 'externals'.
Specifying your compartments in a LIM input file can look like this:
## EXTERNALS
CO2
EXPORT
## END EXTERNALS
## COMPARTMENTS
Plant = 700 {g C m-2}
Animal = 50 {g C m-2}
deadDetritus = 1000 {g C m-2} ! contains the default "dead" tag
xxxDetritus = 1000 {g C m-2} ! chosen tag "xxx" is equally valid if you specify this for the fwstability package.
## END COMPARTMENTS
Specifying flows can look like this:
## FLOWS
NetPrimProd : CO2 -> Plant
animalGrazPlant : Plant -> Animal
animalGrazDetritus: deadDetritus -> Animal
animalResp : Animal -> CO2
animalMort : Animal -> EXPORT
animalYYY : Animal -> EXPORT ! equally valid mortality tag, if specified for the package later.
plantMort : Plant -> deadDetritus
animalDef : Animal -> deadDetritus
detritusMineralization: deadDetritus -> CO2
## END FLOWS
Make sure you use the compartment names as specified in "COMPARTMENT", so including the tag for dead compartments. The name of defecation and mortality flows should only contain the name of the compartment and the tag. If it contains more words, like "thisAnimalMort" the package ignores this flow. If, for some reason, you want to include multiple parallel defecation or mortality flows, you need to define a variable (in the VARIABLE section) that includes the sum of all defecation and mortality flows. Like this:
## FLOWS
animalDefOne : Animal -> deadDetritusOne
animalDefTwo : Animal -> deadDetritusTwo
## END FLOWS
## VARIABLES
animalIngestion = animalGrazPlant + animalGrazDetritus
animalIngestion = flowto(Animal) ! Useful feature of LIM package
animalDef = animalDefOne + animalDefTwo ! Sum all defecation flows
animalAss = animalIngestion - animalDef ! Default assimilation tag "ass"
animalGrowth = animalAss * animalGrowthEfficiency ! Default growht tag "growth"
## END VARIABLES
When your LIM input file is done, you can read it into R, and optionally solve it with the LIM package.
# Use the LIM package for importing your model
# install.packages("LIM")
library(LIM)
lim <- Read(system.file("extdata", "foodweb.lim", package = "fwstability"))
# This example is only one of many ways to solve a LIM.
lim_solved <- Ldei(Setup(lim))
The model should be included in a list that can later be parsed into the getJacobian function. The getJacobian matrix will need a list with at least the elements "type" and "model". The element "type" must be set to "LIM", and the element "model" is the read-in LIM (beware, not yet set up!). If there is no solution specified, the parsimonious solution is used (minimizing the sum of squares).
model <- list(
type = "LIM",
model = lim
)
# Ready to use in the fwstability package!
If you have solved the model, you can include it in your list under the element "web" so that it can be used by the package. This element must be a named vector with flow values and their names.
model <- list(
type = "LIM",
model = lim,
web = lim_solved$X
)
# Ready to use in the fwstability package!
If you use your own tags in the model, for example "RIP" in stead of "MORT", you need to specify this in the list under the elements "aTag", "gTag", "mTag", "deadTag", and/or "defTag". If one or more of these elements are omitted, the default tags are used.
# Example when I want to use my own tags for growth, mortality and defecation.
model <- list(
type = "LIM",
model = lim,
gTag = "production",
mTag = "RIP",
defTag = "poop"
)
If setting up your model takes a long time, you can include your set-up model in the element "setup" so that it is not setup within the function.
setupLIM <- Setup(lim)
mymodel <- list(
type = "LIM",
model = lim,
setup = setupLIM
)
For a working example, please refer to 5.3 A complex application.
2.2. Energy flux model
If you don't use a linear inverse model (LIM) to mass-balance your system, or your LIM is not formatted correctly to be used by the fwstability package, you will need to extract the necessary data first. In order to obtain a Jacobian matrix you need:
A named flow matrix containing the total flow of material/energy between compartments.
A named vector with biomasses for all compartments (including the dead ones).
A named vector with assimilation efficiencies for all non-dead compartments. An organism assimilates (A) the part of all ingested food (I) that is not defecated (D): $A = I-D$. So, assimilation efficiency (AE) is the assimilated amount of food (A) divided by the total of ingested food (I): $AE = A/I$. (see diagram under quick start)
A named vector with growth efficiencies for all non-dead compartments. The part of assimilated food (A) that is not respired (R) is used for growth or secondary production (G): $G = A-R$. Therefore, the growth efficiency (GE) is the total growth (G) divided by the assimilated part of food (A): $GE = G/A$. (see diagram under quick start)
If you have dead compartments in your system, you should specify a list containing two elements with the names of dead compartments (element "names"), and a matrix with the fraction of each flow that constitutes defecation (element "frac"). This matrix has the same size as the flow matrix and indicates the fraction of each flow that comprises defecation. Take the example where there is a material flow from Species A to Detritus of 2 units due to defecation, and another of 8 units due to mortality of the population. The total flow between Species A and Detritus is 10 units (this is the value present in the flow matrix). Of those total 10 units, 2 units are caused by defecation. Therefore, the frac_matrix will have the value 2/10 = 0.2. If there is only a flux between Species A and Detritus due to defecation, the fraction will be 1. For all flows which are not comprised of defecation at all, for example when Species A eats Species B, the fraction in the frac_matrix will be 0.
A character vector containing the names of external compartments (these are removed from the flowmatrix before the calculations).
The above information should be combined in a list that also contains the element "type" which is set to "EF" (short for energy flux) so that it can be parsed into the getJacobian function. It will look like this:
model = list(
# EF stands for energy flux
type = "EF",
# a named square flowmatrix with fluxes between all compartments
FM = FM,
# a named numeric vector with biomasses for all compartments
BM = BM,
# a named numeric vector with assimilation efficiencies for all compartments
AE = AE,
# a named numeric vector with growth efficiencies for all compartments
GE = GE,
# info on the dead (detritus) compartments
dead = list(
# character vector with names of all dead (detritus) compartments
names = dead_names,
# a matrix with the fraction defecation of each flow in the flowmatrix
frac = frac_matrix
),
# optional: a character vector with the names of external compartments
# that should not be included in the calculatinons,
externals = externals,
# optional: a named numeric vector with mortality rates (per unit time);
MR = MR
)
Lets take a look at the data in an example food web from the fwmodel package: a food web of soil treated with conventional farming practices.
# Information on the Lovinkhoeve food web model
?Lovinkhoeve
# This food web model is an energy flux ("EF") model
LovinkhoeveCP$type
# This food web model consists of 17 compartments.
dim(LovinkhoeveCP$FM)
rownames(LovinkhoeveCP$FM)
# The biomasses (kg C / ha) of all compartments are
LovinkhoeveCP$BM
# The assimilation and production (growth) efficiency of all compartments are
LovinkhoeveCP$AE
LovinkhoeveCP$GE
# The natural death (mortality) rates ( /yr ) of all compartments are
LovinkhoeveCP$MR
# The flow matrix contains calculated flow sizes between compartments (kg C / ha / yr).
# The first five rows and columns of the flow matrix can be viewed as:
LovinkhoeveCP$FM[1:5, 1:5]
# Mortality and defecation fluxes can be distinguished using a matrix that
# specifies the fraction of each flow comprised of defecation.
# The first five rows and columns of this matrix can be viewed as:
LovinkhoeveCP$dead$frac[1:5, 1:5]
If you have a LIM with the wrong data format, the functions getFlowMatrix and getVariables may still be of use to easily calculate the needed data. The function getFlowMatrix is similar to the function Flowmatrix from the LIM package, but enhanced so that it also works if parallel flows exist between the same two compartments. The function getVariables returns a named list with the values of all defined variables calculated from the food web flows. So, for example, if you have defined variables related to assimilation but for some reason didn't use the tag system, you can easily extract the variables values and find your assimilation efficiency.
# Use the LIM package for importing your model
require(LIM)
lim <- Read(system.file("extdata", "foodweb.lim", package = "fwstability"))
FM <- getFlowMatrix(lim)
print(FM)
vars <- getVariables(lim)
print(vars)
3. Obtaining a Jacobian matrix
Mathematically speaking, the Jacobian matrix contains the first order partial derivatives of a set of differential equations. If the differential equations describe the change in biomass over time for a set of food web compartments, the Jacobian elements can be interpreted as interaction strengths: the effect of population i on population j.
In the fwstability package population i that exerts an effect is in rows, and population j that receives the effect is in columns. If you want to, you can change this by transposing the matrix. If the differential equations describe the total change in biomass per time unit $\frac{dBiomass}{dt}$, then the unit of the Jacobian matrix is per unit time and is also called the Community matrix. If the differential equations describe something else (like the per capita change in biomass per unit time), the unit of the Jacobian matrix will also be different (for example per unit biomass per unit time) and is therefore more difficult to ecologically interpret.
An example of working with units: If you have a food web model with flows in the unit g Carbon $m^{-2} d^{-1}$ and biomass in g Carbon $m^{-2}$, the interaction strengths will have the unit $d^{-1}$.
In order to translate the energy-flux food-web model from the static domain to the dynamic domain, the system is described in Lotka-Volterra type equations and fluxes are equated to elements of these differential equations. For full explanation, please refer to the paper. In the getJacobian function, the effect of a predator $j$ on its prey $i$ ($\alpha_{ij}$) is calculated as:
$$\alpha_{ij}=-\frac{F_{ij}}{x{j}}$$
and the effect of a prey $i$ on its predator $j$ ($\alpha_{ji}$, which is also valid for the effect of detritus on its consumer) is calculated as:
$$\alpha_{ji}=\frac{a_jp_jF_{ij}}{x{i}}$$
where $F_{ij}$ is the flux between resource population $i$ and consumer population $j$, $a_j$ is the assimilation efficiency of the consumer, $p_j$ is the production efficiency of the consumer, and $x_{i,j}$ is the biomass of the respective population (de Ruiter et al., 1995; Neutel et al., 2002).
The effect of a population $j$ on detritus $d$ ($\alpha_{dj}$) is calculated as:
$$\alpha_{dj}=\frac{1}{x_j}(\sum_{i=1}^n (1-a_j)F_{ij}\epsilon_{jd} + \sum_{i=1}^n (1-a_i)F_{ji}\epsilon_{id} + (\sum_{i=1}^na_jp_jF_{ij}-\sum_{i=1}^nF_{ji})\delta_{jd} - F_{dj})$$
where $n$ is the total number of compartments in the food web, $\epsilon_{jd}$ is the fraction of total egestion by $j$ that flows into detritus $d$, and $\delta_{jd}$ is the fraction of total mortality of $j$ that flows into detritus $d$.
The effect of another detritus compartment $o$ on detritus $d$ ($\alpha_{do}$) is calculated as:
$$\alpha_{do}=\frac{1}{x_o}(\sum_{i=1}^n(1-a_i)F_{oi}\epsilon_{id}+F_{od}$$
The detritus functions used to construct Jacobian matrices in previous literature (de Ruiter 1995, Neutel et al. 2002, Neutel et al. 2007, Neutel & Thorne 2014) were only taking into consideration one detritus compartment. The presented equations used in this package allow the inclusion of multiple detritus compartments in the food-web, including linear ontogeny of detritus. For details, review our paper.
The diagonal of a Jacobian matrix represents intraspecific interactions. The default diagonal of the function getJacobian() is calculated from the model. Diagonal values can also be set to a single chosen value or numeric vector. Diagonal values for fauna are found as their non-predatory mortality (Neutel & Thorne, 2014):
$$\alpha_{ii}=-\frac{\sum_{j=1}^n a_i p_i F_{ji}-\sum_{j=1}^n F_{ij}}{x_i}$$
and for detritus as:
$$\alpha_{dd}=-\frac{1}{x_d} \sum_{j=1}^n (a_j F_{dj}+\sum_{o=1}^nF_{do})$$
The function for detritus diagonal values is adjusted from Neutel & Thorne (2014) to allow density-dependent detritus ontogeny. The diagonal values found with these functions represent the upper-limit of self-dampening effect possible in the system. Therefore, Jacobian matrixes with model-based diagonals are often stable. If you want to compare stability between two food web models that are both stable, you can redefine stability as the minimum fraction of natural mortality needed to stabilize the matrix (see 3. Stability analysis).
The getJacobian function should be used like this:
# If you have an energy-flux model:
JM <- getJacobian(
model = list(
type = "EF",
# a named square flowmatrix with fluxes between all compartments
FM = FM,
# a named numeric vector with biomasses for all compartments
BM = BM,
# a named numeric vector with assimilation efficiencies for all compartments
AE = AE,
# a named numeric vector with growth efficiencies for all compartments
GE = GE,
# info on the dead (detritus) compartments
dead = list(
# character vector with names of all dead (detritus) compartments
names = dead_names,
# a matrix with the fraction defecation of each flow in the flowmatrix
frac = frac_matrix
),
# a character vector with the names of external compartments
# that should not be included in the calculatinons,
externals = externals,
# a named numeric vector with mortality rates (per unit time);
# only required if you want the diagonal to be "model"
MR = MR
),
# how to treat the diagonal. Default is "model", can also be a number or numeric vector
diagonal = "model"
)
# If you have a Linear Inverse Model formatted correctly with tags
getJacobian(
model = list(
type = "LIM",
# A read-in LIM model
LIM = Read(path_to_lim_file),
# Optionally a setup LIM model
setup = Setup(Read(path_to_lim_file)),
# A named numeric vector with all flow solutions.
# If not given the least sum of squares solution is found.
web = flow_solutions,
# tag in the LIM assigned to the total assimilation of each compartment
aTag = "ass",
# tag in the LIM assigned to the total growth (production) of each compartment
gTag = "growth",
# tag in the LIM assigned to the total mortality of each compartment
mTag = "mort",
# tag in the LIM included in the compartment name of dead (detritus) compartments
deadTag = "dead"
),
# how to treat the diagonal. Default is "model", can also be a number or numeric vector
diagonal = "model",
)
4. Stability analysis
A Jacobian matrix is only stable if all real parts of the eigenvalues are negative. Therefore, the stability of a Jacobian matrix can be measured as the maximum real value of the eigenvalues, also called the dominant eigenvalue. A negative value indicates a stable matrix, a positive value indicates an unstable matrix. If the dominant eigenvalue is negative (thus stable), the absolute value of the dominant eigenvalue is the rate of return of the slowest trajectory back to equilibrium. If the Jacobian matrix has the unit per unit time, the rate of return is also per unit time. The return rate can be used to calculate return time.
# Find the dominant eigenvalue
s <- getStability(JM)
return_rate <- abs(s) # per unit time (e.g. per day or per year)
return_time <- 1/return_rate # unit of time (e.g. days or years)
The diagonal of the Jacobian matrix is strongly correlated to stability (the sum of the diagonal values equals the sum of eigenvalues). Therefore, it is important to know what values are on the diagonal before you interpret your results. If you calculated the diagonal values from the model, they represent upper-limits of self-dampening effects. In other words, it is the maximum amount of buffering effect the system can provide to stabilize. If the system is unstable even with upper-limit diagonal values, the natural state of the system is unstable and will collapse.
There are two important pieces of information to consider if you are interpreting return times:
1) The return time does NOT indicate the return time back to equilibrium, but the time to return to about 37% ($1/e$) of the equilibrium state. The reason for this is mathematical and has to do with natural logarithms, and is explained more in depth in our paper.
2) Although return rates and return times can be useful information to review, they are NOT suitable to compared relative stability between ecosystems. Ecosystems can operate on different time-scales, which makes a comparison between return rates and return times not 'fair'. In other words, system A can be much slower to return to equilibrium than system B, but only because system A simply includes many slow-living organisms. If time-dependency is removed from system A and B and stability is evaluated again, system A might turn out to actually be more stable than B despite its longer absolute return times.
In ecology, it is useful and powerful to compare stability between two or more systems. In order to allow such comparison of relative stability, time dependency must be removed. In this case stability is not defined as the slowest return rate to equilibrium (per unit time), but redefined as the dimensionless minimum fraction of self-dampening effect necessary to stabilize the matrix. There are two ways to achieve this:
1) Option 1: Stability can be found as the minimum scalar of natural mortality rates needed to stabilize the system. This method is included as the 'scalar' method in the getStability function, and is based on iteratively increasing the fraction of natural mortality rates used as the Jacobian diagonal with a certain step-size until the Jacobian matrix is found to be stable. This method requires a vector with natural mortality rates and the names of all dead (detritus) compartments.
2) Option 2: The Jacobian matrix is normalized by dividing each interaction strength by absolute value of the the corresponding diagonal value (per unit time divided by per unit time is dimensionless). Then, all diagonal values - except for detritus compartments - are set to zero (self-dampening effects are removed) and the dominant eigenvalue is found. The matrix can be normalized with the normalizeJacobian function. This matrix requires a fully quantified diagonal with negative, non-zero, values.
In both cases it is assumed detritus feedbacks are constant and cannot be scaled.
Both these values can be interpreted as the fraction of maximum possible self-dampening effects needed to stabilize the matrix (Neutel & Thorne, 2014). In other words, we expect this value to lie between 0 and 1. If the value is larger than 1 it means more than available self-dampening is needed to stabilize the matrix, thus the system as it is will collapse. These time-independent stability measure can be interpreted as a biological tipping point, and can be used to compare relative stability between systems. If the system has less self-dampening than the fraction self-dampening as indicated by the stability measure, it will collapse.
The scalar method (1) is used in the first papers reviewing stability of different systems. However, the iterative algorithmic implementation is more time-consuming and the resolution of the stability value depends on the step-size of the iterations. Therefore, we recommend finding the dominant eigenvalue of the normalized zero-diagonal Jacobian matrix (2) instead of the scalar method.
Finding time-independent measures of stability:
# Using the scalar method
# You need to provide mortality rates (MR) and
# a character vector with names of dead compartments.
s <- getStability(JM, method = "scalar", MR, dead_names)
# Finding the dominant eigenvalue of a normalized zero-diagonal Jacobian
# This method requires a quantified non-zero diagonal and
# a character vector with names of dead compartments.
JMnorm <- normalizeJacobian(JM, dead_names)
s <- getStability(JMnorm)
5. Analyze food web for (de)stabilizing characteristics
It can be useful for your research to review other network characteristics, to try to understand why stability between systems differ. This package does NOT provide a comprehensive set of functions to do such analysis, but it does provide a starting point.
Effect individual compartments and links
Functions are included to review the effect on stability when individual compartments are removed from the system, or when certain interaction strengths are altered.
The function assessComps removes individual compartments and reassesses the stability. It returns a dataframe with a row for every compartment, and the absolute and relative change in stability if that compartment is removed. If the change in stability (delta) is negative, the system becomes more stable if the respective food web compartment is removed from the Jacobian matrix. If delta is positive the system becomes less stable. You need to specify what stability method you want to use (default is the eigenvalue method, but the scalar method is also possible).
df <- assessComps(JM)
It is also possible to review how altering the interaction strengths impacts the stability measure. There two functions included for this:
1) assessLinksFixed. This function alters each individual interaction strength (one cell of the Jacobian) once by a fixed function (default is times 2) and review the effect of stability. Thus, it returns a dataframe with as many rows as nrows(JM) times ncols(JM), with relative and absolute changes in stability. Again, if the change in stability (delta) is negative, the system becomes more stable if the respective food web compartment is removed from the Jacobian matrix. If delta is positive the system becomes less stable.
2) assessLinksPerm. This function alters a pair of interaction strengths (i on j, AND j on i) in a number of permutations to be between 0 and 2 times the initial interaction strength. It returns a dataframe that is half of the dataframe created by assessLinksFixed, and included the probability that the matrix becomes unstable if the respective interaction pair is altered.
In both functions you can specify which stability measure you would like to use, the default is the eigenvalue method.
df <- assessLinksFixed(JM, func = function(x){return(x*2)})
df <- assessLinksPerm(JM, perms = 100, threshold = 0.01)
Feedback loops
Different scientific publications have shown that looking at individual interaction strengths might not be very insightful to understand stability of systems. It is often better to look at the patterning of interaction strengths, specifically in feedback loops. Therefore, the function assessFeedback is added that wraps multiple functions:
- First it finds all loops (of a certain length k if specified) in your network. The indices of each loop (eg 1 > 2 > 3 > 1) are stored in a textfile. Finding these loops is done with a recursive depth-first-search function. Beware that it can take a really long time to find all loops in a large network! In that case it might be better to only review loops of length 2 and 3, that are known to most heavily influence stability (see works by Neutel).
- Then for each loop it finds the feedback i.e. the product of all interaction strengths, and the loop weight i.e. the geometric mean of all absolute interaction strengths in the loop.
The function returns a dataframe with all loops noted as 'compName1->compName2->compNamek' (column "loop"), with feedbacks (column "fdb") and loop weights (column "lw") of those loops.
How to use this function:
# The following statements finds ALL loops in the system,
# stores these in the textfile "allLoops.txt" (default output name),
# and finds each loop feedback and scaled loop weight.
loops.df <- assessFeedback(JM, findLoops = TRUE)
# If you only want to review loops of a certain length
loops.df <- assessFeedback(JM, k = 3, findLoops = TRUE)
# If you don't want to scale your loop weight
# (e.g when you also don't scale your stability values)
loops.df <- assessFeedback(JM, findLoops = TRUE, scale = FALSE)
If you have a really large network, it might not be useful to find all loops in the system, and to only focus on loops of length 2 and 3. You can do this by using the depth-first search algorithm (dfs) separately from the assessFeeback function, and then run the assessFeedback function by reading in the files created by the dfs function calls:
# Create an adjacency matrix that indicates which interactions exist,
# i.e. only 0 and 1 in stead of quantified interaction strengths.
AM <- abs(JM)
AM[which(AM > 0)] <- 1
# Find loops of length 2 and 3 and store in text files
# "allLoops_k=2.txt" and "allLoops_k=3.txt".
file_name_loop2 <- dfs(AM, k = 2, verbose)
file_name_loop3 <- dfs(AM, k = 3, verbose)
# Evaluate loop feedback and weight for these loops of length 2 and 3.
loops2 <- assessFeedback(JM, findLoops = FALSE, file = file_name_loop2)
loops3 <- assessFeedback(JM, findLoops = FALSE, file = file_name_loop3)
It is also possible to assess the skewness of fluxes, which is known to be related to stability, through calculating weighted connectance (van Altena et al. 2016). A lower value indicates a more skewed flux distribution, which is found to be correlated to less stable systems.
Cw <- getCw(FM)
References
Allesina, S., & Bondavalli, C. (2003). Steady state of ecosystem flow networks: A comparison between balancing procedures. Ecological Modelling, 165(2–3), 221–229. https://doi.org/10.1016/S0304-3800(03)00075-9
Allesina, S., Grilli, J., Barabás, G., Tang, S., Aljadeff, J., & Maritan, A. (2015). Predicting the stability of large structured food webs. Nature Communications, 6, 1–6. https://doi.org/10.1038/ncomms8842
Allesina, S., & Tang, S. (2012). Stability criteria for complex ecosystems. Nature, 483(7388), 205–208. https://doi.org/10.1038/nature10832
Allesina, S., & Tang, S. (2015). The stability–complexity relationship at age 40: a random matrix perspective. Population Ecology, 57(1), 63–75. https://doi.org/10.1007/s10144-014-0471-0
Arnoldi, J. F., Loreau, M., & Haegeman, B. (2016). Resilience, reactivity and variability: A mathematical comparison of ecological stability measures. Journal of Theoretical Biology, 389, 47–59. https://doi.org/10.1016/j.jtbi.2015.10.012
Berlow, E. L., Neutel, A. M., Cohen, J. E., De Ruiter, P. C., Ebenman, B., Emmerson, M., … Petchey, O. (2004). Interaction Strengths in Food Webs: Issues and Opportunities. Journal of Animal Ecology, 73(3), 585–598.
Boit, A., & Gaedke, U. (2014). Benchmarking successional progress in a quantitative food web. PLoS ONE, 9(2). https://doi.org/10.1371/journal.pone.0090404
Botton, S., Van Heusden, M., Parsons, J. R., Smidt, H., & Van Straalen, N. (2006). Resilience of microbial systems towards disturbances. Critical Reviews in Microbiology, 32(2), 101–112. https://doi.org/10.1080/10408410600709933
Christensen, V., & Pauly, D. (1992). Ecopath II – a software for balancing steady-state ecosystem models and calculating network characteristics. Ecological Modelling, 61, 169–185.
Davis, R. C. (1981). Structure and Function of Two Antarctic Terrestrial Moss Communities. Ecological Monographs, 51(2), 125–143. https://doi.org/10.2307/2937260
de Ruiter, P. C., Neutel, A. M., & Moore, J. C. (1995). Energetics, Patterns of Interaction Strengths, and Stability in Real Ecosystems. Science, 269(5228), 1257–1260. https://doi.org/10.1126/science.269.5228.1257
de Ruiter, P. C., Van Veen, J. A., Moore, J. C., Brussaard, L., & Hunt, H. W. (1993). Calculation of nitrogen mineralization in soil food webs. Plant and Soil, 157(2), 263–273. https://doi.org/10.1007/BF00011055
Drazen, J. C., & Sutton, T. T. (2017). Dining in the Deep: The Feeding Ecology of Deep-Sea Fishes. Annual Review of Marine Science, 9(1), 337–366. https://doi.org/10.1146/annurev-marine-010816-060543
Fath, B. D., Scharler, U. M., Ulanowicz, R. E., & Hannon, B. (2007). Ecological network analysis: network construction. Ecological Modelling, 208(1), 49–55. https://doi.org/10.1016/j.ecolmodel.2007.04.029
Heal, O. W., & MacLean, S. F. (1975). Comparative productivity in ecosystems - secondary productivity. In W. H. Van Dobben & R. H. Lowe-McConnell (Eds.), Unifying Concepts in Ecology (pp. 89–108). The Hague: Dr. W. Junk.
Heymans, J. J., Coll, M., Libralato, S., Morissette, L., & Christensen, V. (2014). Global patterns in ecological indicators of marine food webs: A modelling approach. PLoS ONE, 9(4). https://doi.org/10.1371/journal.pone.0095845
Holling, C. S. (1959). Some characteristics of simple types of predation and parasitism. Canadian Entomology, 91, 385–398.
Hoving, H.-J. T., Perez, J. A. A., Bolstad, K. S. R., Braid, H. E., Evans, A. B., Fuchs, D., … Xavier, J. C. C. (2014). The study of deep-sea cephalopods. In Advances in Marine Biology (Vol. 67, pp. 235–359). Elsevier Ltd. https://doi.org/10.1016/B978-0-12-800287-2.00003-2
Hunt, H. W., Coleman, D. C., Ingham, E. R., Ingham, R. E., Elliott, E. T., Moore, J. C., … Morley, C. R. I. (1987). The detrital food web in a shortgrass prairie. Biol, 3, 57–68.
Jacquet, C., Moritz, C., Morissette, L., Legagneux, P., Massol, F., Archambault, P., & Gravel, D. (2016). No complexity–stability relationship in empirical ecosystems. Nature Communications, 7, 12573. https://doi.org/10.1038/ncomms12573
Kuiper, J. J., Van Altena, C., de Ruiter, P. C., Van Gerven, L. P. A. A., Janse, J. H., & Mooij, W. M. (2015). Food-web stability signals critical transitions in temperate shallow lakes. Nature Communications, 6, 1–7. https://doi.org/10.1038/ncomms8727
Landi, P., Minoarivelo, H. O., Brännström, Å., Hui, C., & Dieckmann, U. (2018). Complexity and stability of ecological networks: a review of the theory. Population Ecology, 60(4), 319–345. https://doi.org/10.1007/s10144-018-0628-3
Latham, L. G. (2006). Network flow analysis algorithms. Ecological Modelling, 192(3–4), 586–600. https://doi.org/10.1016/j.ecolmodel.2005.07.029
Logofet, D. O. (1993). Matrices and Graphs: Stability Problems in Mathematical Ecology. Boca Raton, FL: CRCPress.
May, R. M. (1972). Will a large network be stable? Nature, 238, 37–38. https://doi.org/10.1038/238413a0
Montoya, J. M., Woodward, G., Emmerson, M. C., & Solé, R. V. (2009). Press Perturbations and Indirect Effects in Real Food Webs. Ecology, 90(9), 2426–2433.
Moore, J. C., Berlow, E. L., Coleman, D. C., De Ruiter, P. C., Dong, Q., Hastings, A., … Wall, D. H. (2004). Detritus, trophic dynamics and biodiversity. Ecology Letters, 7(7), 584–600. https://doi.org/10.1111/j.1461-0248.2004.00606.x
Moore, J. C., & de Ruiter, P. C. (2012). Energetic Food Webs (First). Oxford, UK: Oxford University Press.
Moore, J. C., Ruiter, P. C. De, & Hunt, H. W. (1993). Influence of Productivity on the Stability of Real and Model Ecosystems. Science, 261, 906–908.
Neubert, M. G., & Caswell, H. (1997). Alternatives to resilience for measuring the responses of ecological systems to perturbations. Ecology, 78(3), 653.
Neutel, A. M., Heesterbeek, J. A. P., & Ruiter, P. C. De. (2002). Stability in Real Food Webs: Weak Links in Long Loops. Science, 296(5570), 1120–1123. https://doi.org/10.1126/science.1068326
Neutel, A. M., Heesterbeek, J. A. P., Van De Koppel, J., Hoenderboom, G., Vos, A., Kaldeway, C., … De Ruiter, P. C. (2007). Reconciling complexity with stability in naturally assembling food webs. Nature, 449(7162), 599–602. https://doi.org/10.1038/nature06154
Neutel, A. M., & Thorne, M. A. S. (2014). Interaction strengths in balanced carbon cycles and the absence of a relation between ecosystem complexity and stability. Ecology Letters, 17(6), 651–661. https://doi.org/10.1111/ele.12266
Nilsson, K. A., & McCann, K. S. (2016). Interaction strength revisited—clarifying the role of energy flux for food web stability. Theoretical Ecology, 9(1), 59–71. https://doi.org/10.1007/s12080-015-0282-8
Novak, M., Yeakel, J. D., Noble, A. E., Doak, D. F., Emmerson, M., Estes, J. A., … Wootton, J. T. (2016). Characterizing Species Interactions to Understand Press Perturbations: What Is the Community Matrix? Annual Review of Ecology, Evolution, and Systematics, 47(1), 409–432. https://doi.org/10.1146/annurev-ecolsys-032416-010215
Pimm, S. L., Lawton, J. H., & Cohen, J. E. (1991). Food web patterns and their consequences. Nature, 350, 669–674.
Ramirez-Llodra, E., Brandt, A., Danovaro, R., De Mol, B., Escobar, E., German, C. R., … Vecchione, M. (2010). Deep, diverse and definitely different: Unique attributes of the world’s largest ecosystem. Biogeosciences, 7(9), 2851–2899. https://doi.org/10.5194/bg-7-2851-2010
Soetaert, K. (2009a). Package NetIndices , network indices and food web descriptors in R. Ecology, 1, 1–10.
Soetaert, K. (2009b). rootSolve: Nonlinear root finding, equilibrium and steady-state analysis of ordinary differential equations.
Soetaert, K., & Herman, P. M. J. (2009). A Practical Guide to Ecological Modelling. Using R as a Simulation Platform. Springer.
Soetaert, K., & Van Oevelen, D. (2010). Package LIM, implementing linear inverse models in R, 37.
van Altena, C., Hemerik, L., & de Ruiter, P. C. (2016). Food web stability and weighted connectance: the complexity-stability debate revisited. Theoretical Ecology, 9(1), 49–58. https://doi.org/10.1007/s12080-015-0291-7
van Oevelen, D., Soetaert, K., Garcia, R., de Stigter, H. C., Cunha, M. R., Pusceddu, A., & Danovaro, R. (2011). Canyon conditions impact carbon flows in food webs of three sections of the Nazaré canyon. Deep-Sea Research Part II: Topical Studies in Oceanography, 58(23–24), 2461–2476. https://doi.org/10.1016/j.dsr2.2011.04.009
van Oevelen, D., Soetaert, K., & Heip, C. (2012). Carbon flows in the benthic food web of the Porcupine Abyssal Plain: The (un)importance of labile detritus in supporting microbial and faunal carbon demands. Limnology and Oceanography, 57(2), 645–664. https://doi.org/10.4319/lo.2012.57.2.0645
van Oevelen, D., van den Meersche, K., Meysman, F. J. R., Soetaert, K., Middelburg, J. J., & Vézina, A. F. (2010). Quantifying food web flows using linear inverse models. Ecosystems, 13(1), 32–45. https://doi.org/10.1007/s10021-009-9297-6
Vézina, A. F., & Platt, T. (1988). Food web dynamics in the ocean. I. Best-estimates of flow networks using inverse methods. Marine Ecology Progress Series, 42, 269–287.
Wilkinson, J. H. (1965). The Algebraic Eigenvalue Problem. Oxford, UK: Clarendon Press.
Wilson, E. E., & Wolkovich, E. M. (2011). Scavenging: How carnivores and carrion structure communities. Trends in Ecology and Evolution, 26(3), 129–135. https://doi.org/10.1016/j.tree.2010.12.011
Zwart, K. B., Burgers, S. L. G. E., Bloem, J., Bouwman, L. A., Brussaard, L., Lebbink, G., … de Ruiter, P. C. (1994). Population dynamics in the belowground food webs in two different agricultural systems. Agriculture, Ecosystems & Environment, 51(1), 187–198. https://doi.org/10.1016/0167-8809(94)90043-4
dswdejonge/fwstability documentation built on Dec. 7, 2022, 7:24 p.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
knitr::opts_chunk$set( collapse = TRUE, comment = "#>", eval = FALSE )
The fwstability package is developed to infer food-web stability from energy-flux models. Basically, an energy-flux food-web model shows the transfer of energy or materials from one species to another when they eat each other. These models are steady-state and mass-balanced. Mathematical methods have been developed to transform these steady-state models to the dynamic realm in order to infer stability. This package implements these mathematical methods, so that you can easily infer stability from your energy-flux food-web model.
We are currently preparing a paper introducing this package and an in-depth explanation of the science behind energy-flux food-web stability. Access to the manuscript is available upon reasonable request until publication:
Daniëlle S.W. de Jonge, Peter C. de Ruiter, Johan van de Koppel, Dick van Oevelen. Inferring stability of energy-flux food-web models. In preparation.
This vignette provides you a quick start, and a more detailed explanation about the methods. Reading the vignette is highly recommended so that you can interpret your results correctly.
1. Quick start
1.1. Install and load
To install the development version of the package:
# if you don't have the devtools package yet: install.packages("devtools") # Install the fwstability package, development version: devtools::install_github("dswdejonge/fwstability", build_vignettes = TRUE)
There is also the fwmodels package available from github, that contains several published energy-flux food-web models that are used as examples in this vignette.
# To use the example food-webs, install the fwmodels package as well devtools::install_github("dswdejonge/fwmodels")
If the installation was successful, we load the fwstability package (and the fwmodels package for examples).
library(fwstability) library(fwmodels)
1.2. Prepare data
First, you must prepare your energy-flux food-web model to be used by the package. See section 2 for the requirements in the physiology of the food-web model.
The flowmatrix always has source compartments in rows, and sink compartments in columns. In other words, the value in each cell is the flow moving from the row compartment toward to the column compartment.
A case that perhaps requires some extra explanation is the 'frac_matrix'. This matrix has the same size as the flow matrix and indicates the fraction of each flow that comprises defecation. Take the example where there is a material flow from Species A to Detritus of 2 units due to defecation, and another of 8 units due to mortality of the population. The total flow between Species A and Detritus is 10 units (this is the value present in the flow matrix). Of those total 10 units, 2 units are caused by defecation. Therefore, the frac_matrix will have the value 2/10 = 0.2. If there is only a flux between Species A and Detritus due to defecation, the fraction will be 1. For all flows which are not comprised of defecation at all, for example when Species A eats Species B, the fraction in the frac_matrix will be 0.
Dead compartment are all compartments that are not alive, e.g. detritus, nutrients, etc. Information about these compartments should be included in a model element called 'dead'.
External compartments are compartments from which material is imported into and exported from the model, and which can be disregarded for inferring interaction strengths and stability. If present, their names can be included in a model element called 'externals'.
A hypothetical example of what your data can look like:
# 1. Biomasses for each food web compartment # (e.g. abundance per unit area). compartments <- c("detritus", "speciesA", "speciesB") biomasses <- c(30, 20, 10) # e.g. grams per square meter names(biomasses) <- compartments # 2. Quantified fluxes between compartments # (e.g. abundance per unit area per unit time). detritus_to_speciesA <- 10 speciesA_to_detritus <- 3 speciesA_to_speciesB <- 7 speciesB_to_detritus <- 7 flowmatrix <- matrix(data = 0, nrow = length(compartments), ncol = length(compartments)) row.names(flowmatrix) <- compartments colnames(flowmatrix) <- compartments flowmatrix["detritus", "speciesA"] <- detritus_to_speciesA flowmatrix["speciesA", "detritus"] <- speciesA_to_detritus flowmatrix["speciesA", "speciesB"] <- speciesA_to_speciesB flowmatrix["speciesB", "detritus"] <- speciesB_to_detritus # 3. Energy conversion efficiencies for all living compartments # (assimilation and production efficiencies, fraction 0-1). assimilation_eff <- c(NA, 0.8, 0.7) growth_eff <- c(NA, 0.6, 0.4) names(assimilation_eff) <- compartments names(growth_eff) <- compartments # 4. Information on dead compartments, if any # 4.a. The names of all dead compartments dead_names <- "detritus" # 4.b. The fraction of each flow which comprises defecation # frac_matrix needs to be the same size and the flow matrix frac_matrix <- matrix(data = 0, nrow = length(compartments), ncol = length(compartments)) row.names(frac_matrix) <- compartments colnames(frac_matrix) <- compartments # Defecation is the part of ingested food (food source -> species) # which is NOT assimilated (1 - assimilation efficiency) # For more explanation about physiology, see Section 2. defecation_by_speciesA <- flowmatrix["detritus", "speciesA"] * (1-assimilation_eff["speciesA"]) defecation_by_speciesB <- flowmatrix["speciesA", "speciesB"] * (1-assimilation_eff["speciesB"]) # You can simply divide the defecation flow by the total flow # between compartments to obtain the value for the frac_matrix frac_matrix["speciesA", "detritus"] <- defecation_by_speciesA / flowmatrix["speciesA", "detritus"] frac_matrix["speciesB", "detritus"] <- defecation_by_speciesB / flowmatrix["speciesB", "detritus"]
Combine all information into a named list as follows:
mymodel = list( type = "EF", # EF stands for energy flux FM = flowmatrix, BM = biomasses, AE = assimilation_eff, GE = growth_eff, dead = list( names = dead_names, frac = frac_matrix))
1.3. Jacobian matrix
Now we can estimate a Jacobian matrix from the energy-flux model (see section 3 and our article for more detail).
The Jacobian matrix shows the per unit biomass effect of species i (in rows) on the population growth rate of species j (in columns) per unit time. In other words, the Jacobian matrix shows you the interaction strengths between species. Thus, the unit of the JM is 'per unit time'.
The diagonal values are intraspecific effects, thus the density-dependent effects a compartment has on itself (works for both species and detritus).
We use an example food web from the fwmodels package. We recommend using the default diagonal, i.e. calculated from the model (example JM1 below):
library(fwstability) library(fwmodels) model <- LovinkhoeveCP # Documentation for the getJacobian function ?getJacobian # Create a Jacobian matrix, with a diagonal calculated from the model (default). # The diagonal values are the upper-limit of self-dampening effects in the system. JM1 <- getJacobian(model, diagonal = "model") # Create a Jacobian matrix, with an all-zero diagonal. JM2 <- getJacobian(model, diagonal = 0) # Create a Jacobian matrix, with -1 on the diagonal. JM3 <- getJacobian(model, diagonal = -1) # Create a Jacobian matrix, with a numeric vector as diagonal diagonal <- runif(dim(model$FM)[1], min = -1, max = 0) JM4 <- getJacobian(model, diagonal = diagonal) # The preview of the first 5 rows and columns of the Jacobian matrix: print(JM1[1:5, 1:5])
1.4. Infer stability
Now we can infer stability from the Jacobian matrix. All methods are based on the concept that a Jacobian matrix is only stable if all real parts of the eigenvalues are negative.
If the system is stable the system will asymptotically return to equilibrium after a small pulse perturbation. If the system is unstable the system will not return to equilibrium after a small pulse perturbation, but move away from the equilibrium state.
Method s1 (below) finds the dominant eigenvalue of the Jacobian matrix. If positive, the system is NOT stable. If negative the system is stable. In a stable system the dominant eigenvalue represent the return rate to equilibrium. Method s1 is NOT suitable to compare relative stability between different systems if these systems operate on different time-scales.
Method s2 (below) also finds the dominant eigenvalue, but from a normalized Jacobian matrix with an all-zero diagonal (except for detritus). A normalized JM is unitless, and can therefore be used to compare relative stability between systems. The value of s2 is expected to lie between 0 and 1, and represents the minimum fraction of self-dampening necessary to stabilize the system. A lower value of s2 means a more stable system. If s2 > 1, more than available self-dampening is necessary for a stable system. Hence, the system is inherently unstable. We recommend this method.
Method s3 (below) is the scalar method, used in early publication (e.g. de Ruiter et al. 1997). The value of s3 can be interpreted the same as s2 (Neutel & Thorne 2014), however, the resolution of s2 is probably better than the resolution of s3, and s2 is quicker to calculate.
Both method s2 and s3 assume the JM without self-dampening effects on the diagonal is unstable.
# --- Method s1 --- # Calculate stability with the eigenvalue method (default), # from a Jacobian matrix with unit per time. s1 <- getStability(JM1) # If s1 is positive, the system is NOT stable. # If s1 is negative, the system is stable. # If s1 is stable, the asymptotic return rate is: return_rate <- abs(s1) # per unit time (e.g. per day) # If s1 is stable, the time for the system to return # to about 37% of its original state is: return_time <- 1/return_rate # unit of time (e.g. days or years) # BEWARE: s1 and the derived return rates and return times # are NOT suitable to compare relative stability among different systems. # --- Method s2 --- # Calculate stability with the eigenvalue method (default), # from a normalized Jacobian matrix with an all-zero diagonal. # Normalize JM (JM must have quantified diagonal): JMnorm <- normalizeJacobian(JM1, dead_names = model$dead$names) # Find dominant eigenvalue: s2 <- getStability(JMnorm) # s2 is expected to lie between 0 and 1 # s2 represents the minimum fraction of available self-dampening # necessary to stabilize the system. # --- Method 3 --- # Calculate stability with the scalar method. # Natural mortality rates (MR in per unit time) and the names of all # dead (i.e. detritus) compartments are necessary. s3 <- getStability(JM2, method = "scalar", MR = model$MR, dead_names = model$dead$names) # If the diagonal is calculated from the model, # the diagonal can be used as MR s3 <- getStability(JM1, method = "scalar", MR = abs(diag(JM1)), dead_names = model$dead$names)
2. Preparing a food web model for stability analysis
Your energy-flux model should always stick to the physiological concepts explained in this diagram:
2.1. A linear inverse model
In a linear inverse model (LIM) all processes in a food web are captured by a set of linear equations. These linear equations can be either equalities (the outcome of an equation should equal another value) or inequalities (the outcome of an equation should be equal or greater than another value). The equations can parameterized with known data, but often not all processes are studied well enough to quantify. The unknowns in the equations can then be found with linear programming, thereby fully quantifying your model. A simple example of a linear optimization problem can be found on \link{Wikipedia}[https://en.wikipedia.org/wiki/Linear_programming#Example].
A linear inverse food web model (LIM) can be prepared with the R-package LIM (https://CRAN.R-project.org/package=LIM). More information on using LIM for food web quantification can be found in:
van Oevelen, D., van den Meersche, K., Meysman, F.J.R., Soetaert, K., Middelburg, J.J., Vézina, A.F., 2010. Quantifying food web flows using linear inverse models. Ecosystems 13, 32–45. https://doi.org/10.1007/s10021-009-9297-6.
A LIM model produced by the LIM package can be used directly by the fwstability package if it is formatted in a certain way, as explained below. If your LIM is not formatted according to the requirements of the fwstability package, you will need to extract certain data from the model yourself before using the package. For the latter, please refer to 2.2. Energy-flux model.
In order to obtain a Jacobian matrix, certain physiological processes must be calculated from your LIM. The fwstability can automatically extract this data if you labelled them with a unique tag. In the model you specify:
- which compartments are dead compartments (e.g. detritus and carrion, default tag is "dead")
- which flows specifically represent defecation (default tag is "def")
- which flows specifically represent mortality (default tag is "mort")
- which variables represent total assimilation (default tag is "ass", beware to not tag assimilation efficiency)
- which variables represent total growth (secondary production, default tag is "growth", beware to not tag growth efficiency)
You can use the default tags, but it is also possible to use your own tags. In the latter case you will need to provide the package with the name of your chosen tags in order to extract the correct information. The package will warn you if it cannot find the tag in your model.
Rules and tips for using tags:
- The official component name is the name given in the 'stocks' section. Sometimes, stock values are already given in the 'parameter' section and then assigned to the defined stocks in the 'stocks' section. Be aware not to tag the compartments in the 'parameter' section as these are not the official stock names used throughout the rest of the model analysis.
- Flows that represent defecation and mortality must be named in the 'flows' section and include the defecation or mortality tag. If defecation and mortality are combined in the same flow, they should still be declared separately (so you get the same flow twice, but with a different name).
- The variables representing total defecation, mortality, assimilation, and growth should be named by combining the component name and tag only. The code subtracts the tag from the name of the variable and uses it to place the data with the right compartment. So, a variable containing total rabbit mortality should be named "rabbitMort". If you would name it "thisisRabbitMort" it is not recognized due to the superfluous text. It is possible to specify multiple mortality flows, e.g. if the carcass of the rabbit consists of both high-energy and low-energy material which you want to separate into two separate carcass compartments. In this case it is necessary to define rabbitMort as the sum of these flows in the variable section of the model.
- Tags are not case sensitive. So, RABBITmort and RaBbItMoRt are both valid.
- The tag may occur anywhere in the string (so, it does not have to be added specifically at the beginning or ending of the name of the compartment/flow/variable). Ensure that the rest of names do not accidentally contain a tag! For example, if you include the Indian wild ass (Equus hemionus khur) in your foodweb, a variable representing its respiration called respirationAss will be interpreted as the assimilation variable of an unknown compartment 'respiration' because it contains the default assimilation tag "ass".
- Dead compartments, defecation, and mortality can be omitted from your model if you don't want to include them (for example make them all export to an external compartment). However, including assimilation and growth is required for all non-dead compartments. If you have no growth or assimilation data (for example for bacteria) you can assume the full ingestion is used for assimilation and/or growth and define your variable accordingly.
General tips for setting up a LIM:
- It is recommended to be consistent in your naming of flows and variables so later you can easily extract data. For naming flows for example: sourceProcessSink.
- In order to prevent certain processes to be quantified as zero, include minimum constraints. This is especially important if you later want to normalize the Jacobian matrix, because this process does not allow zeroes on the diagonal. For example, if you don't set a minimum constraint for non-predatory mortality, the parsimonious solution might set this flow to zero. This is not realistic and it the occurrence of these zeroes may hamper your future calculations.
- The inclusion of biomasses in the LIM is not required to be able to solve the model. However, biomasses are required for all internal compartments in order to calculate the Jacobian matrix. If there are compartments that do not have biomasses in the LIM, they must either be quantified anyway or excluded from the getJacobian() function by including them as 'externals'.
Specifying your compartments in a LIM input file can look like this:
## EXTERNALS
CO2
EXPORT
## END EXTERNALS
## COMPARTMENTS
Plant = 700 {g C m-2}
Animal = 50 {g C m-2}
deadDetritus = 1000 {g C m-2} ! contains the default "dead" tag
xxxDetritus = 1000 {g C m-2} ! chosen tag "xxx" is equally valid if you specify this for the fwstability package.
## END COMPARTMENTS
Specifying flows can look like this:
## FLOWS
NetPrimProd : CO2 -> Plant
animalGrazPlant : Plant -> Animal
animalGrazDetritus: deadDetritus -> Animal
animalResp : Animal -> CO2
animalMort : Animal -> EXPORT
animalYYY : Animal -> EXPORT ! equally valid mortality tag, if specified for the package later.
plantMort : Plant -> deadDetritus
animalDef : Animal -> deadDetritus
detritusMineralization: deadDetritus -> CO2
## END FLOWS
Make sure you use the compartment names as specified in "COMPARTMENT", so including the tag for dead compartments. The name of defecation and mortality flows should only contain the name of the compartment and the tag. If it contains more words, like "thisAnimalMort" the package ignores this flow. If, for some reason, you want to include multiple parallel defecation or mortality flows, you need to define a variable (in the VARIABLE section) that includes the sum of all defecation and mortality flows. Like this:
## FLOWS
animalDefOne : Animal -> deadDetritusOne
animalDefTwo : Animal -> deadDetritusTwo
## END FLOWS
## VARIABLES
animalIngestion = animalGrazPlant + animalGrazDetritus
animalIngestion = flowto(Animal) ! Useful feature of LIM package
animalDef = animalDefOne + animalDefTwo ! Sum all defecation flows
animalAss = animalIngestion - animalDef ! Default assimilation tag "ass"
animalGrowth = animalAss * animalGrowthEfficiency ! Default growht tag "growth"
## END VARIABLES
When your LIM input file is done, you can read it into R, and optionally solve it with the LIM package.
# Use the LIM package for importing your model # install.packages("LIM") library(LIM) lim <- Read(system.file("extdata", "foodweb.lim", package = "fwstability")) # This example is only one of many ways to solve a LIM. lim_solved <- Ldei(Setup(lim))
The model should be included in a list that can later be parsed into the getJacobian function. The getJacobian matrix will need a list with at least the elements "type" and "model". The element "type" must be set to "LIM", and the element "model" is the read-in LIM (beware, not yet set up!). If there is no solution specified, the parsimonious solution is used (minimizing the sum of squares).
model <- list( type = "LIM", model = lim ) # Ready to use in the fwstability package!
If you have solved the model, you can include it in your list under the element "web" so that it can be used by the package. This element must be a named vector with flow values and their names.
model <- list( type = "LIM", model = lim, web = lim_solved$X ) # Ready to use in the fwstability package!
If you use your own tags in the model, for example "RIP" in stead of "MORT", you need to specify this in the list under the elements "aTag", "gTag", "mTag", "deadTag", and/or "defTag". If one or more of these elements are omitted, the default tags are used.
# Example when I want to use my own tags for growth, mortality and defecation. model <- list( type = "LIM", model = lim, gTag = "production", mTag = "RIP", defTag = "poop" )
If setting up your model takes a long time, you can include your set-up model in the element "setup" so that it is not setup within the function.
setupLIM <- Setup(lim) mymodel <- list( type = "LIM", model = lim, setup = setupLIM )
For a working example, please refer to 5.3 A complex application.
2.2. Energy flux model
If you don't use a linear inverse model (LIM) to mass-balance your system, or your LIM is not formatted correctly to be used by the fwstability package, you will need to extract the necessary data first. In order to obtain a Jacobian matrix you need:
A named flow matrix containing the total flow of material/energy between compartments.
A named vector with biomasses for all compartments (including the dead ones).
A named vector with assimilation efficiencies for all non-dead compartments. An organism assimilates (A) the part of all ingested food (I) that is not defecated (D): $A = I-D$. So, assimilation efficiency (AE) is the assimilated amount of food (A) divided by the total of ingested food (I): $AE = A/I$. (see diagram under quick start)
A named vector with growth efficiencies for all non-dead compartments. The part of assimilated food (A) that is not respired (R) is used for growth or secondary production (G): $G = A-R$. Therefore, the growth efficiency (GE) is the total growth (G) divided by the assimilated part of food (A): $GE = G/A$. (see diagram under quick start)
If you have dead compartments in your system, you should specify a list containing two elements with the names of dead compartments (element "names"), and a matrix with the fraction of each flow that constitutes defecation (element "frac"). This matrix has the same size as the flow matrix and indicates the fraction of each flow that comprises defecation. Take the example where there is a material flow from Species A to Detritus of 2 units due to defecation, and another of 8 units due to mortality of the population. The total flow between Species A and Detritus is 10 units (this is the value present in the flow matrix). Of those total 10 units, 2 units are caused by defecation. Therefore, the frac_matrix will have the value 2/10 = 0.2. If there is only a flux between Species A and Detritus due to defecation, the fraction will be 1. For all flows which are not comprised of defecation at all, for example when Species A eats Species B, the fraction in the frac_matrix will be 0.
A character vector containing the names of external compartments (these are removed from the flowmatrix before the calculations).
The above information should be combined in a list that also contains the element "type" which is set to "EF" (short for energy flux) so that it can be parsed into the getJacobian function. It will look like this:
model = list( # EF stands for energy flux type = "EF", # a named square flowmatrix with fluxes between all compartments FM = FM, # a named numeric vector with biomasses for all compartments BM = BM, # a named numeric vector with assimilation efficiencies for all compartments AE = AE, # a named numeric vector with growth efficiencies for all compartments GE = GE, # info on the dead (detritus) compartments dead = list( # character vector with names of all dead (detritus) compartments names = dead_names, # a matrix with the fraction defecation of each flow in the flowmatrix frac = frac_matrix ), # optional: a character vector with the names of external compartments # that should not be included in the calculatinons, externals = externals, # optional: a named numeric vector with mortality rates (per unit time); MR = MR )
Lets take a look at the data in an example food web from the fwmodel package: a food web of soil treated with conventional farming practices.
# Information on the Lovinkhoeve food web model ?Lovinkhoeve # This food web model is an energy flux ("EF") model LovinkhoeveCP$type # This food web model consists of 17 compartments. dim(LovinkhoeveCP$FM) rownames(LovinkhoeveCP$FM) # The biomasses (kg C / ha) of all compartments are LovinkhoeveCP$BM # The assimilation and production (growth) efficiency of all compartments are LovinkhoeveCP$AE LovinkhoeveCP$GE # The natural death (mortality) rates ( /yr ) of all compartments are LovinkhoeveCP$MR # The flow matrix contains calculated flow sizes between compartments (kg C / ha / yr). # The first five rows and columns of the flow matrix can be viewed as: LovinkhoeveCP$FM[1:5, 1:5] # Mortality and defecation fluxes can be distinguished using a matrix that # specifies the fraction of each flow comprised of defecation. # The first five rows and columns of this matrix can be viewed as: LovinkhoeveCP$dead$frac[1:5, 1:5]
If you have a LIM with the wrong data format, the functions getFlowMatrix and getVariables may still be of use to easily calculate the needed data. The function getFlowMatrix is similar to the function Flowmatrix from the LIM package, but enhanced so that it also works if parallel flows exist between the same two compartments. The function getVariables returns a named list with the values of all defined variables calculated from the food web flows. So, for example, if you have defined variables related to assimilation but for some reason didn't use the tag system, you can easily extract the variables values and find your assimilation efficiency.
# Use the LIM package for importing your model require(LIM) lim <- Read(system.file("extdata", "foodweb.lim", package = "fwstability")) FM <- getFlowMatrix(lim) print(FM) vars <- getVariables(lim) print(vars)
3. Obtaining a Jacobian matrix
Mathematically speaking, the Jacobian matrix contains the first order partial derivatives of a set of differential equations. If the differential equations describe the change in biomass over time for a set of food web compartments, the Jacobian elements can be interpreted as interaction strengths: the effect of population i on population j.
In the fwstability package population i that exerts an effect is in rows, and population j that receives the effect is in columns. If you want to, you can change this by transposing the matrix. If the differential equations describe the total change in biomass per time unit $\frac{dBiomass}{dt}$, then the unit of the Jacobian matrix is per unit time and is also called the Community matrix. If the differential equations describe something else (like the per capita change in biomass per unit time), the unit of the Jacobian matrix will also be different (for example per unit biomass per unit time) and is therefore more difficult to ecologically interpret.
An example of working with units: If you have a food web model with flows in the unit g Carbon $m^{-2} d^{-1}$ and biomass in g Carbon $m^{-2}$, the interaction strengths will have the unit $d^{-1}$.
In order to translate the energy-flux food-web model from the static domain to the dynamic domain, the system is described in Lotka-Volterra type equations and fluxes are equated to elements of these differential equations. For full explanation, please refer to the paper. In the getJacobian function, the effect of a predator $j$ on its prey $i$ ($\alpha_{ij}$) is calculated as: $$\alpha_{ij}=-\frac{F_{ij}}{x{j}}$$ and the effect of a prey $i$ on its predator $j$ ($\alpha_{ji}$, which is also valid for the effect of detritus on its consumer) is calculated as: $$\alpha_{ji}=\frac{a_jp_jF_{ij}}{x{i}}$$ where $F_{ij}$ is the flux between resource population $i$ and consumer population $j$, $a_j$ is the assimilation efficiency of the consumer, $p_j$ is the production efficiency of the consumer, and $x_{i,j}$ is the biomass of the respective population (de Ruiter et al., 1995; Neutel et al., 2002).
The effect of a population $j$ on detritus $d$ ($\alpha_{dj}$) is calculated as:
$$\alpha_{dj}=\frac{1}{x_j}(\sum_{i=1}^n (1-a_j)F_{ij}\epsilon_{jd} + \sum_{i=1}^n (1-a_i)F_{ji}\epsilon_{id} + (\sum_{i=1}^na_jp_jF_{ij}-\sum_{i=1}^nF_{ji})\delta_{jd} - F_{dj})$$
where $n$ is the total number of compartments in the food web, $\epsilon_{jd}$ is the fraction of total egestion by $j$ that flows into detritus $d$, and $\delta_{jd}$ is the fraction of total mortality of $j$ that flows into detritus $d$.
The effect of another detritus compartment $o$ on detritus $d$ ($\alpha_{do}$) is calculated as:
$$\alpha_{do}=\frac{1}{x_o}(\sum_{i=1}^n(1-a_i)F_{oi}\epsilon_{id}+F_{od}$$
The detritus functions used to construct Jacobian matrices in previous literature (de Ruiter 1995, Neutel et al. 2002, Neutel et al. 2007, Neutel & Thorne 2014) were only taking into consideration one detritus compartment. The presented equations used in this package allow the inclusion of multiple detritus compartments in the food-web, including linear ontogeny of detritus. For details, review our paper.
The diagonal of a Jacobian matrix represents intraspecific interactions. The default diagonal of the function getJacobian() is calculated from the model. Diagonal values can also be set to a single chosen value or numeric vector. Diagonal values for fauna are found as their non-predatory mortality (Neutel & Thorne, 2014): $$\alpha_{ii}=-\frac{\sum_{j=1}^n a_i p_i F_{ji}-\sum_{j=1}^n F_{ij}}{x_i}$$ and for detritus as: $$\alpha_{dd}=-\frac{1}{x_d} \sum_{j=1}^n (a_j F_{dj}+\sum_{o=1}^nF_{do})$$ The function for detritus diagonal values is adjusted from Neutel & Thorne (2014) to allow density-dependent detritus ontogeny. The diagonal values found with these functions represent the upper-limit of self-dampening effect possible in the system. Therefore, Jacobian matrixes with model-based diagonals are often stable. If you want to compare stability between two food web models that are both stable, you can redefine stability as the minimum fraction of natural mortality needed to stabilize the matrix (see 3. Stability analysis).
The getJacobian function should be used like this:
# If you have an energy-flux model: JM <- getJacobian( model = list( type = "EF", # a named square flowmatrix with fluxes between all compartments FM = FM, # a named numeric vector with biomasses for all compartments BM = BM, # a named numeric vector with assimilation efficiencies for all compartments AE = AE, # a named numeric vector with growth efficiencies for all compartments GE = GE, # info on the dead (detritus) compartments dead = list( # character vector with names of all dead (detritus) compartments names = dead_names, # a matrix with the fraction defecation of each flow in the flowmatrix frac = frac_matrix ), # a character vector with the names of external compartments # that should not be included in the calculatinons, externals = externals, # a named numeric vector with mortality rates (per unit time); # only required if you want the diagonal to be "model" MR = MR ), # how to treat the diagonal. Default is "model", can also be a number or numeric vector diagonal = "model" ) # If you have a Linear Inverse Model formatted correctly with tags getJacobian( model = list( type = "LIM", # A read-in LIM model LIM = Read(path_to_lim_file), # Optionally a setup LIM model setup = Setup(Read(path_to_lim_file)), # A named numeric vector with all flow solutions. # If not given the least sum of squares solution is found. web = flow_solutions, # tag in the LIM assigned to the total assimilation of each compartment aTag = "ass", # tag in the LIM assigned to the total growth (production) of each compartment gTag = "growth", # tag in the LIM assigned to the total mortality of each compartment mTag = "mort", # tag in the LIM included in the compartment name of dead (detritus) compartments deadTag = "dead" ), # how to treat the diagonal. Default is "model", can also be a number or numeric vector diagonal = "model", )
4. Stability analysis
A Jacobian matrix is only stable if all real parts of the eigenvalues are negative. Therefore, the stability of a Jacobian matrix can be measured as the maximum real value of the eigenvalues, also called the dominant eigenvalue. A negative value indicates a stable matrix, a positive value indicates an unstable matrix. If the dominant eigenvalue is negative (thus stable), the absolute value of the dominant eigenvalue is the rate of return of the slowest trajectory back to equilibrium. If the Jacobian matrix has the unit per unit time, the rate of return is also per unit time. The return rate can be used to calculate return time.
# Find the dominant eigenvalue s <- getStability(JM) return_rate <- abs(s) # per unit time (e.g. per day or per year) return_time <- 1/return_rate # unit of time (e.g. days or years)
The diagonal of the Jacobian matrix is strongly correlated to stability (the sum of the diagonal values equals the sum of eigenvalues). Therefore, it is important to know what values are on the diagonal before you interpret your results. If you calculated the diagonal values from the model, they represent upper-limits of self-dampening effects. In other words, it is the maximum amount of buffering effect the system can provide to stabilize. If the system is unstable even with upper-limit diagonal values, the natural state of the system is unstable and will collapse.
There are two important pieces of information to consider if you are interpreting return times:
1) The return time does NOT indicate the return time back to equilibrium, but the time to return to about 37% ($1/e$) of the equilibrium state. The reason for this is mathematical and has to do with natural logarithms, and is explained more in depth in our paper.
2) Although return rates and return times can be useful information to review, they are NOT suitable to compared relative stability between ecosystems. Ecosystems can operate on different time-scales, which makes a comparison between return rates and return times not 'fair'. In other words, system A can be much slower to return to equilibrium than system B, but only because system A simply includes many slow-living organisms. If time-dependency is removed from system A and B and stability is evaluated again, system A might turn out to actually be more stable than B despite its longer absolute return times.
In ecology, it is useful and powerful to compare stability between two or more systems. In order to allow such comparison of relative stability, time dependency must be removed. In this case stability is not defined as the slowest return rate to equilibrium (per unit time), but redefined as the dimensionless minimum fraction of self-dampening effect necessary to stabilize the matrix. There are two ways to achieve this:
1) Option 1: Stability can be found as the minimum scalar of natural mortality rates needed to stabilize the system. This method is included as the 'scalar' method in the getStability function, and is based on iteratively increasing the fraction of natural mortality rates used as the Jacobian diagonal with a certain step-size until the Jacobian matrix is found to be stable. This method requires a vector with natural mortality rates and the names of all dead (detritus) compartments.
2) Option 2: The Jacobian matrix is normalized by dividing each interaction strength by absolute value of the the corresponding diagonal value (per unit time divided by per unit time is dimensionless). Then, all diagonal values - except for detritus compartments - are set to zero (self-dampening effects are removed) and the dominant eigenvalue is found. The matrix can be normalized with the normalizeJacobian function. This matrix requires a fully quantified diagonal with negative, non-zero, values.
In both cases it is assumed detritus feedbacks are constant and cannot be scaled.
Both these values can be interpreted as the fraction of maximum possible self-dampening effects needed to stabilize the matrix (Neutel & Thorne, 2014). In other words, we expect this value to lie between 0 and 1. If the value is larger than 1 it means more than available self-dampening is needed to stabilize the matrix, thus the system as it is will collapse. These time-independent stability measure can be interpreted as a biological tipping point, and can be used to compare relative stability between systems. If the system has less self-dampening than the fraction self-dampening as indicated by the stability measure, it will collapse.
The scalar method (1) is used in the first papers reviewing stability of different systems. However, the iterative algorithmic implementation is more time-consuming and the resolution of the stability value depends on the step-size of the iterations. Therefore, we recommend finding the dominant eigenvalue of the normalized zero-diagonal Jacobian matrix (2) instead of the scalar method.
Finding time-independent measures of stability:
# Using the scalar method # You need to provide mortality rates (MR) and # a character vector with names of dead compartments. s <- getStability(JM, method = "scalar", MR, dead_names) # Finding the dominant eigenvalue of a normalized zero-diagonal Jacobian # This method requires a quantified non-zero diagonal and # a character vector with names of dead compartments. JMnorm <- normalizeJacobian(JM, dead_names) s <- getStability(JMnorm)
5. Analyze food web for (de)stabilizing characteristics
It can be useful for your research to review other network characteristics, to try to understand why stability between systems differ. This package does NOT provide a comprehensive set of functions to do such analysis, but it does provide a starting point.
Effect individual compartments and links
Functions are included to review the effect on stability when individual compartments are removed from the system, or when certain interaction strengths are altered.
The function assessComps removes individual compartments and reassesses the stability. It returns a dataframe with a row for every compartment, and the absolute and relative change in stability if that compartment is removed. If the change in stability (delta) is negative, the system becomes more stable if the respective food web compartment is removed from the Jacobian matrix. If delta is positive the system becomes less stable. You need to specify what stability method you want to use (default is the eigenvalue method, but the scalar method is also possible).
df <- assessComps(JM)
It is also possible to review how altering the interaction strengths impacts the stability measure. There two functions included for this:
1) assessLinksFixed. This function alters each individual interaction strength (one cell of the Jacobian) once by a fixed function (default is times 2) and review the effect of stability. Thus, it returns a dataframe with as many rows as nrows(JM) times ncols(JM), with relative and absolute changes in stability. Again, if the change in stability (delta) is negative, the system becomes more stable if the respective food web compartment is removed from the Jacobian matrix. If delta is positive the system becomes less stable.
2) assessLinksPerm. This function alters a pair of interaction strengths (i on j, AND j on i) in a number of permutations to be between 0 and 2 times the initial interaction strength. It returns a dataframe that is half of the dataframe created by assessLinksFixed, and included the probability that the matrix becomes unstable if the respective interaction pair is altered.
In both functions you can specify which stability measure you would like to use, the default is the eigenvalue method.
df <- assessLinksFixed(JM, func = function(x){return(x*2)}) df <- assessLinksPerm(JM, perms = 100, threshold = 0.01)
Feedback loops
Different scientific publications have shown that looking at individual interaction strengths might not be very insightful to understand stability of systems. It is often better to look at the patterning of interaction strengths, specifically in feedback loops. Therefore, the function assessFeedback is added that wraps multiple functions:
- First it finds all loops (of a certain length k if specified) in your network. The indices of each loop (eg 1 > 2 > 3 > 1) are stored in a textfile. Finding these loops is done with a recursive depth-first-search function. Beware that it can take a really long time to find all loops in a large network! In that case it might be better to only review loops of length 2 and 3, that are known to most heavily influence stability (see works by Neutel).
- Then for each loop it finds the feedback i.e. the product of all interaction strengths, and the loop weight i.e. the geometric mean of all absolute interaction strengths in the loop.
The function returns a dataframe with all loops noted as 'compName1->compName2->compNamek' (column "loop"), with feedbacks (column "fdb") and loop weights (column "lw") of those loops.
How to use this function:
# The following statements finds ALL loops in the system, # stores these in the textfile "allLoops.txt" (default output name), # and finds each loop feedback and scaled loop weight. loops.df <- assessFeedback(JM, findLoops = TRUE) # If you only want to review loops of a certain length loops.df <- assessFeedback(JM, k = 3, findLoops = TRUE) # If you don't want to scale your loop weight # (e.g when you also don't scale your stability values) loops.df <- assessFeedback(JM, findLoops = TRUE, scale = FALSE)
If you have a really large network, it might not be useful to find all loops in the system, and to only focus on loops of length 2 and 3. You can do this by using the depth-first search algorithm (dfs) separately from the assessFeeback function, and then run the assessFeedback function by reading in the files created by the dfs function calls:
# Create an adjacency matrix that indicates which interactions exist, # i.e. only 0 and 1 in stead of quantified interaction strengths. AM <- abs(JM) AM[which(AM > 0)] <- 1 # Find loops of length 2 and 3 and store in text files # "allLoops_k=2.txt" and "allLoops_k=3.txt". file_name_loop2 <- dfs(AM, k = 2, verbose) file_name_loop3 <- dfs(AM, k = 3, verbose) # Evaluate loop feedback and weight for these loops of length 2 and 3. loops2 <- assessFeedback(JM, findLoops = FALSE, file = file_name_loop2) loops3 <- assessFeedback(JM, findLoops = FALSE, file = file_name_loop3)
It is also possible to assess the skewness of fluxes, which is known to be related to stability, through calculating weighted connectance (van Altena et al. 2016). A lower value indicates a more skewed flux distribution, which is found to be correlated to less stable systems.
Cw <- getCw(FM)
References
Allesina, S., & Bondavalli, C. (2003). Steady state of ecosystem flow networks: A comparison between balancing procedures. Ecological Modelling, 165(2–3), 221–229. https://doi.org/10.1016/S0304-3800(03)00075-9
Allesina, S., Grilli, J., Barabás, G., Tang, S., Aljadeff, J., & Maritan, A. (2015). Predicting the stability of large structured food webs. Nature Communications, 6, 1–6. https://doi.org/10.1038/ncomms8842
Allesina, S., & Tang, S. (2012). Stability criteria for complex ecosystems. Nature, 483(7388), 205–208. https://doi.org/10.1038/nature10832
Allesina, S., & Tang, S. (2015). The stability–complexity relationship at age 40: a random matrix perspective. Population Ecology, 57(1), 63–75. https://doi.org/10.1007/s10144-014-0471-0
Arnoldi, J. F., Loreau, M., & Haegeman, B. (2016). Resilience, reactivity and variability: A mathematical comparison of ecological stability measures. Journal of Theoretical Biology, 389, 47–59. https://doi.org/10.1016/j.jtbi.2015.10.012
Berlow, E. L., Neutel, A. M., Cohen, J. E., De Ruiter, P. C., Ebenman, B., Emmerson, M., … Petchey, O. (2004). Interaction Strengths in Food Webs: Issues and Opportunities. Journal of Animal Ecology, 73(3), 585–598.
Boit, A., & Gaedke, U. (2014). Benchmarking successional progress in a quantitative food web. PLoS ONE, 9(2). https://doi.org/10.1371/journal.pone.0090404
Botton, S., Van Heusden, M., Parsons, J. R., Smidt, H., & Van Straalen, N. (2006). Resilience of microbial systems towards disturbances. Critical Reviews in Microbiology, 32(2), 101–112. https://doi.org/10.1080/10408410600709933
Christensen, V., & Pauly, D. (1992). Ecopath II – a software for balancing steady-state ecosystem models and calculating network characteristics. Ecological Modelling, 61, 169–185.
Davis, R. C. (1981). Structure and Function of Two Antarctic Terrestrial Moss Communities. Ecological Monographs, 51(2), 125–143. https://doi.org/10.2307/2937260
de Ruiter, P. C., Neutel, A. M., & Moore, J. C. (1995). Energetics, Patterns of Interaction Strengths, and Stability in Real Ecosystems. Science, 269(5228), 1257–1260. https://doi.org/10.1126/science.269.5228.1257
de Ruiter, P. C., Van Veen, J. A., Moore, J. C., Brussaard, L., & Hunt, H. W. (1993). Calculation of nitrogen mineralization in soil food webs. Plant and Soil, 157(2), 263–273. https://doi.org/10.1007/BF00011055
Drazen, J. C., & Sutton, T. T. (2017). Dining in the Deep: The Feeding Ecology of Deep-Sea Fishes. Annual Review of Marine Science, 9(1), 337–366. https://doi.org/10.1146/annurev-marine-010816-060543
Fath, B. D., Scharler, U. M., Ulanowicz, R. E., & Hannon, B. (2007). Ecological network analysis: network construction. Ecological Modelling, 208(1), 49–55. https://doi.org/10.1016/j.ecolmodel.2007.04.029
Heal, O. W., & MacLean, S. F. (1975). Comparative productivity in ecosystems - secondary productivity. In W. H. Van Dobben & R. H. Lowe-McConnell (Eds.), Unifying Concepts in Ecology (pp. 89–108). The Hague: Dr. W. Junk.
Heymans, J. J., Coll, M., Libralato, S., Morissette, L., & Christensen, V. (2014). Global patterns in ecological indicators of marine food webs: A modelling approach. PLoS ONE, 9(4). https://doi.org/10.1371/journal.pone.0095845
Holling, C. S. (1959). Some characteristics of simple types of predation and parasitism. Canadian Entomology, 91, 385–398.
Hoving, H.-J. T., Perez, J. A. A., Bolstad, K. S. R., Braid, H. E., Evans, A. B., Fuchs, D., … Xavier, J. C. C. (2014). The study of deep-sea cephalopods. In Advances in Marine Biology (Vol. 67, pp. 235–359). Elsevier Ltd. https://doi.org/10.1016/B978-0-12-800287-2.00003-2
Hunt, H. W., Coleman, D. C., Ingham, E. R., Ingham, R. E., Elliott, E. T., Moore, J. C., … Morley, C. R. I. (1987). The detrital food web in a shortgrass prairie. Biol, 3, 57–68.
Jacquet, C., Moritz, C., Morissette, L., Legagneux, P., Massol, F., Archambault, P., & Gravel, D. (2016). No complexity–stability relationship in empirical ecosystems. Nature Communications, 7, 12573. https://doi.org/10.1038/ncomms12573
Kuiper, J. J., Van Altena, C., de Ruiter, P. C., Van Gerven, L. P. A. A., Janse, J. H., & Mooij, W. M. (2015). Food-web stability signals critical transitions in temperate shallow lakes. Nature Communications, 6, 1–7. https://doi.org/10.1038/ncomms8727
Landi, P., Minoarivelo, H. O., Brännström, Å., Hui, C., & Dieckmann, U. (2018). Complexity and stability of ecological networks: a review of the theory. Population Ecology, 60(4), 319–345. https://doi.org/10.1007/s10144-018-0628-3
Latham, L. G. (2006). Network flow analysis algorithms. Ecological Modelling, 192(3–4), 586–600. https://doi.org/10.1016/j.ecolmodel.2005.07.029
Logofet, D. O. (1993). Matrices and Graphs: Stability Problems in Mathematical Ecology. Boca Raton, FL: CRCPress.
May, R. M. (1972). Will a large network be stable? Nature, 238, 37–38. https://doi.org/10.1038/238413a0
Montoya, J. M., Woodward, G., Emmerson, M. C., & Solé, R. V. (2009). Press Perturbations and Indirect Effects in Real Food Webs. Ecology, 90(9), 2426–2433.
Moore, J. C., Berlow, E. L., Coleman, D. C., De Ruiter, P. C., Dong, Q., Hastings, A., … Wall, D. H. (2004). Detritus, trophic dynamics and biodiversity. Ecology Letters, 7(7), 584–600. https://doi.org/10.1111/j.1461-0248.2004.00606.x
Moore, J. C., & de Ruiter, P. C. (2012). Energetic Food Webs (First). Oxford, UK: Oxford University Press.
Moore, J. C., Ruiter, P. C. De, & Hunt, H. W. (1993). Influence of Productivity on the Stability of Real and Model Ecosystems. Science, 261, 906–908.
Neubert, M. G., & Caswell, H. (1997). Alternatives to resilience for measuring the responses of ecological systems to perturbations. Ecology, 78(3), 653.
Neutel, A. M., Heesterbeek, J. A. P., & Ruiter, P. C. De. (2002). Stability in Real Food Webs: Weak Links in Long Loops. Science, 296(5570), 1120–1123. https://doi.org/10.1126/science.1068326
Neutel, A. M., Heesterbeek, J. A. P., Van De Koppel, J., Hoenderboom, G., Vos, A., Kaldeway, C., … De Ruiter, P. C. (2007). Reconciling complexity with stability in naturally assembling food webs. Nature, 449(7162), 599–602. https://doi.org/10.1038/nature06154
Neutel, A. M., & Thorne, M. A. S. (2014). Interaction strengths in balanced carbon cycles and the absence of a relation between ecosystem complexity and stability. Ecology Letters, 17(6), 651–661. https://doi.org/10.1111/ele.12266
Nilsson, K. A., & McCann, K. S. (2016). Interaction strength revisited—clarifying the role of energy flux for food web stability. Theoretical Ecology, 9(1), 59–71. https://doi.org/10.1007/s12080-015-0282-8
Novak, M., Yeakel, J. D., Noble, A. E., Doak, D. F., Emmerson, M., Estes, J. A., … Wootton, J. T. (2016). Characterizing Species Interactions to Understand Press Perturbations: What Is the Community Matrix? Annual Review of Ecology, Evolution, and Systematics, 47(1), 409–432. https://doi.org/10.1146/annurev-ecolsys-032416-010215
Pimm, S. L., Lawton, J. H., & Cohen, J. E. (1991). Food web patterns and their consequences. Nature, 350, 669–674.
Ramirez-Llodra, E., Brandt, A., Danovaro, R., De Mol, B., Escobar, E., German, C. R., … Vecchione, M. (2010). Deep, diverse and definitely different: Unique attributes of the world’s largest ecosystem. Biogeosciences, 7(9), 2851–2899. https://doi.org/10.5194/bg-7-2851-2010
Soetaert, K. (2009a). Package NetIndices , network indices and food web descriptors in R. Ecology, 1, 1–10.
Soetaert, K. (2009b). rootSolve: Nonlinear root finding, equilibrium and steady-state analysis of ordinary differential equations.
Soetaert, K., & Herman, P. M. J. (2009). A Practical Guide to Ecological Modelling. Using R as a Simulation Platform. Springer.
Soetaert, K., & Van Oevelen, D. (2010). Package LIM, implementing linear inverse models in R, 37.
van Altena, C., Hemerik, L., & de Ruiter, P. C. (2016). Food web stability and weighted connectance: the complexity-stability debate revisited. Theoretical Ecology, 9(1), 49–58. https://doi.org/10.1007/s12080-015-0291-7
van Oevelen, D., Soetaert, K., Garcia, R., de Stigter, H. C., Cunha, M. R., Pusceddu, A., & Danovaro, R. (2011). Canyon conditions impact carbon flows in food webs of three sections of the Nazaré canyon. Deep-Sea Research Part II: Topical Studies in Oceanography, 58(23–24), 2461–2476. https://doi.org/10.1016/j.dsr2.2011.04.009
van Oevelen, D., Soetaert, K., & Heip, C. (2012). Carbon flows in the benthic food web of the Porcupine Abyssal Plain: The (un)importance of labile detritus in supporting microbial and faunal carbon demands. Limnology and Oceanography, 57(2), 645–664. https://doi.org/10.4319/lo.2012.57.2.0645
van Oevelen, D., van den Meersche, K., Meysman, F. J. R., Soetaert, K., Middelburg, J. J., & Vézina, A. F. (2010). Quantifying food web flows using linear inverse models. Ecosystems, 13(1), 32–45. https://doi.org/10.1007/s10021-009-9297-6
Vézina, A. F., & Platt, T. (1988). Food web dynamics in the ocean. I. Best-estimates of flow networks using inverse methods. Marine Ecology Progress Series, 42, 269–287.
Wilkinson, J. H. (1965). The Algebraic Eigenvalue Problem. Oxford, UK: Clarendon Press.
Wilson, E. E., & Wolkovich, E. M. (2011). Scavenging: How carnivores and carrion structure communities. Trends in Ecology and Evolution, 26(3), 129–135. https://doi.org/10.1016/j.tree.2010.12.011
Zwart, K. B., Burgers, S. L. G. E., Bloem, J., Bouwman, L. A., Brussaard, L., Lebbink, G., … de Ruiter, P. C. (1994). Population dynamics in the belowground food webs in two different agricultural systems. Agriculture, Ecosystems & Environment, 51(1), 187–198. https://doi.org/10.1016/0167-8809(94)90043-4
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.