In SWotherspoon/QPress: Qualitative Network Analysis

knitr::opts_chunk$set(collapse = TRUE, warning = FALSE, message = FALSE, comment = "#>",
                      fig.height = 6, fig.width = 6)
figrefs <- c()
figref <- function(z) paste0("Figure ", ifelse(z %in% figrefs, which(figrefs == z), { warning ("unknown figure reference: ", z); "??" }))
tablerefs <- c()
tableref <- function(z) paste0("Table ", ifelse(z %in% tablerefs, which(tablerefs == z), { warning ("unknown table reference: ", z); "??" }))

The QPress package provides facilities for qualitatively modelling the impact of a press perturbation upon a network model. This document illustrates the comparison of alternative models in terms of their responses to perturbation. Two examples are provided for experimental mesocosms in (I) a lake environment [@Hulot2000], and (II) Antarctica (Davidson et al. unpublished).

Introduction

An experimental mesocosm is an ecological tool that enables the experimenter to assess assemblage-level responses to environmental change. Mesocosm experiments are performed in the "natural environment" rather than in a laboratory, but in an enclosure that is small enough that key variables can be directly controlled. Because the responses of component species to perturbation can be observed directly in a mesocosm setting, the results from such studies lend themselves to analysis using qualitative network modelling. In particular, these analyses can be used to make inferences regarding the trophic and competitive interactions that drive assemblage-level responses to perturbation in a mesocosm setting and potentially more broadly.

In this document we use the QPress package to analyse two sets of models corresponding to two separate mesocosm experiments. Example I is for experimental mesocosms in a lake in France that were subjected to high and low nutrient treatments [@Hulot2000]. Example II examines the effects of elevated concentrations of CO$_2$ on the Southern Ocean microbial loop, studied in environmentally controlled tanks (referred to as "minicosms" by the authors) in Antarctica (Davidson et al. unpublished).

Example I: Lake mesocosm

r figrefs <- c(figrefs, "Lake")

Hulot et al. [-@Hulot2000] describe an eight-variable lake mesocosm model comprising phosphorous, three algal groups, two size-classes of herbivores, invertebrate carnivores and carnivorous fish (r figref("Lake")). Unambiguous, directional responses of model components to increased phosphorous observed in the mesocosm experiment were increases in each of: phosphorous (Phos), large herbivores (Herb2), fish (Carn2) and periphyton (AlgP). Melbourne-Thomas et al. [-@JMT2012] compare the ability of Hulot et al.'s original model (r figref("Lake")a) and a variant of the model with one edge added and one edge deleted (r figref("Lake")b) to reproduce these changes in response to a positive perturbation to phosphorous. Specifically, the variant model (r figref("Lake")b) includes a self-limitation effect for invertebrate carnivores, and the predator-prey edge from fish to invertebrate carnivores has been deleted.

**`r figref("Lake")`:** Two versions of the eight-variable lake mesocosm model analyzed by Melbourne-Thomas et al. [-@JMT2012] and Hosack et al. [-@Hosack2008], and originally presented by Hulot et al. [-@Hulot2000]. Abbreviations are: Alg1, edible algae; Alg2, protected algae; AlgP, periphyton; Carn1, invertebrate carnivores; Carn2, fish; Herb1, small herbivores; Herb2, large herbivores; Phos, phosphorus.

The models shown in r figref("Lake") were created using Dia and so are read using model.dia. Here we assign the two model objects to lake.a and lake.b.

library(QPress)
lake.a <- model.dia("Lake-a.dia")
lake.b <- model.dia("Lake-b.dia")

r figrefs <- c(figrefs, "Selector1", "Perturb1")

We then simulate for model (a) and use impact.barplot to view simulation outcomes corresponding with observations from the mesocosm experiment (r figref("Selector1")).

simlake.a <- system.simulate(1000, lake.a)
impact.barplot(simlake.a)

Doing the same for model (b) gives r figref("Perturb1").

**`r figref("Selector1")`:** Node selector for the lake mesocosm example.

simlake.a <- system.simulate(1000, lake.a)
simlake.b <- system.simulate(1000, lake.b)
opar <- par(mfrow = c(1, 2))
impact.barplot0(simlake.a, c(0, 0, 0, 0, 0, 0, 0, 1), c(NA, NA, 1, NA, 1, NA, 1, 1), main = "Model (a)")
impact.barplot0(simlake.b, c(0, 0, 0, 0, 0, 0, 0, 1), c(NA, NA, 1, NA, 1, NA, 1, 1), main = "Model (b)")
par(opar)

Approximately 50% of the simulations match with observations for model (a), whereas for model (b) this increases to over 80%. This suggests that model (b) provides a more parsimonious representation of the lake mesocosm system. Melbourne-Thomas et al. [-@JMT2012] provide a direct comparison of the two models using Bayes factors.

Example II: Antarctic mesocosm

r figrefs <- c(figrefs, "Antarctic-a")

This example is based on experiments conducted at Australia's Davis Station in Antarctica, which assessed the responses of a natural community of Antarctic marine microbes from near-shore waters to elevated concentrations of CO$_2$ using environmentally controlled tanks (mesocosms). Davidson et al. (in prep) describe the main components of the Southern Ocean microbial loop represented in these mesocosms and the ecological interactions between them. Specifically:

Small and large phytoplankton cells consume macronutrients and iron (Fe), although iron is potentially less important for growth of small phytoplankton;
Small phytoplankton are consumed by heterotrophic nanoflagellates and microzooplankton;
Large phytoplankton (diatoms) are consumed by microzooplankton;
Bacteria break down DOC (dissolved organic carbon) and are consumed by heterotrophic nanoflagellates;
Phytoplankton, microzooplankton and heterotrophic nanoflagellates supply the DOC pool.

These interactions are represented in a model that was built using Dia (black-coloured edges in r figref("Antarctic-a")).

When this system was subjected to elevated CO$_2$ concentrations in experimental mesocosms, the following responses were observed:

Small phytoplankton increased;
Large phytoplankton decreased;
Heterotrophic nanoflagellates decreased; and
Bacteria increased.

In the model shown in r figref("Antarctic-a") we assume that the direct effects of CO$_2$ correspond with these observed changes.

**`r figref("Antarctic-a")`:** Model (a) representing the effects of CO$_2$ on the Southern Ocean microbial loop, determined directly from observed responses in experimental mesocosms. These direct effects are shown in red.

We then analyse the response this model to a positive press perturbation to CO$_2$. In this case we enforce a specific ordering of the nodes, and also enforce self-limitation for each model variable.

r figrefs <- c(figrefs, "Selector2", "Perturb2")

labels <- c("CO2", "Fe", "Macronutrients", "Small phytoplankton", "Large phytoplankton",
            "Heterotrophic nanoflagellates", "Microzooplankton", "Bacteria", "DOC")
antarctic.a <- model.dia("Antarctic-a.dia", labels = labels)
antarctic.a <- enforce.limitation(antarctic.a)

As for Example I, we then simulate and use impact.barplot (with the selections shown in r figref("Selector2")) to produce r figref("Perturb2").

simantarctic.a <- system.simulate(1000, antarctic.a)
impact.barplot(simantarctic.a)

**`r figref("Selector2")`:** Node selector for the Antarctic mesocosm example..

simantarctic.a <- system.simulate(1000, antarctic.a)
impact.barplot0(simantarctic.a, c(1, 0, 0, 0, 0, 0, 0, 0, 0), c(NA, NA, NA, 1, -1, -1, NA, 1, NA), main = "Model (a)")

r figrefs <- c(figrefs, "Antarctic-b")

For this version of the model, less than 20% of the simulations match with our "Monitor" criteria, suggesting that it does not provide a particularly good representation of the Antarctic microbial loop system or its response to increased CO$_2$ concentrations. Given that increases in dissolved CO$_2$ are generally expected to have negative impacts of the growth of marine organisms, we could assume a direct negative effect of CO$_2$ on each of large phytoplankton, small phytoplankton, heterotrophic nanoflagellates, and also the pool of available iron (r figref("Antarctic-b")).

**`r figref("Antarctic-b")`:** Alternative model (b) representing the effects of CO$_2$ on the Southern Ocean microbial loop.

Reading this alternative representation from a Dia file

antarctic.b <- model.dia("Antarctic-b.dia", labels = labels)
antarctic.b <- enforce.limitation(antarctic.b)

r figrefs <- c(figrefs, "Perturb3")

and exploring simulation outcomes with impact.barplot (with the same "Monitor" criteria as before) produces r figref("Perturb3").

simantarctic.b <- system.simulate(1000, antarctic.b)
impact.barplot(simantarctic.b)

r figrefs <- c(figrefs, "Antarctic-c")

In this case, we still only see approximately 20% of simulations that meet our criteria, suggesting that we need to re-think the structure of our model. Davidson et al. (in prep) indicate that bacteria "transform dissolved organic carbon (DOC)... thereby repackaging carbon and making it available to higher trophic levels". Looking at our original model, we haven't quite captured this in our representation of interactions between bacteria, DOC and macronutrients. In a third version of our microbial loop model (model (c) --- r figref("Antarctic-c")) we therefore assume a negative effect of bacteria on DOC and a positive effect of bacteria on macronutrients (which includes carbon).

r figrefs <- c(figrefs, "Perturb4")

Reading a third Dia version of the model

antarctic.c <- model.dia("Antarctic-c.dia", labels = labels)
antarctic.c <- enforce.limitation(antarctic.c)

and repeating the simulation and impact.barplot steps from above gives r figref("Perturb4").

simantarctic.b <- system.simulate(1000, antarctic.b)
impact.barplot0(simantarctic.b, c(1, 0, 0, 0, 0, 0, 0, 0, 0), c(NA, NA, NA, 1, -1, -1, NA, 1, NA), main = "Model (b)")

**`r figref("Antarctic-c")`:** Alternative model (c) representing the effects of CO$_2$ on the Southern Ocean microbial loop.

simantarctic.c <- system.simulate(1000, antarctic.c)
impact.barplot0(simantarctic.c, c(1, 0, 0, 0, 0, 0, 0, 0, 0), c(NA, NA, NA, 1, -1, -1, NA, 1, NA), main = "Model (c)")