Introduction

The purpose of the soilfoodwebs package is to help analyze and simulate soil food webs. The following five functions are the core of the package:

Calculate the fluxes through a food web given biomass, parameters, and food web structure.
Calculate the direct and indirect contribution of each trophic species (i.e., node) to carbon and nitrogen mineralization.
Calculate food web stability and smin.
Simulate the food web away from equilibrium.
Conduct a detritus decomposition simulation and calculate the decomposition constant.

The package can complete the following tasks using functions built to work with the communities that are input:

Modify the fluxes to balance carbon and nitrogen demands.
Modify the structure of the food web.

# Load the package
library(soilfoodwebs)

How to build a community

The basis of the soilfoodwebs package is a community, which contains two elements: a feeding matrix and a properties database.

This vignette will contain the code to build a simple 5 trophic species (i.e., node) community from scratch and then run it through all the analyses included in the package. A different community is included as an example in the package, called intro_comm.

Building the properties database:

First, we will build the properties database. This database must contain all the columns listed below in this example.

properties = data.frame(ID = c("Predator", "Prey1", "Prey2", "Microbe", "Detritus"), # Name
                        d = c(1,3,0.5,1.2,0), # Death rate
                    a = c(0.61,0.65,0.45,1,1), # Assimilation efficiency
                    p = c(0.5,0.4,0.3,0.3,1), # Production efficiency for carbon (nitrogen production efficiency is assumed to be 1)
                    B = c(0.1,8,5,9000,1380000), # Biomass
                    CN = c(4.5,4.8,5,8,30), # Carbon to nitrogen ratio
                    DetritusRecycling = c(0,0,0,0,1), # proportion of detritus recycling
                    isDetritus = c(0,0,0,0,1), # Boolean: Is this pool detritus?
                    isPlant = c(0,0,0,0,0), # Boolean: Is this pool a plant?
                    canIMM = c(0,0,0,1,0)) # Boolean: Can the pool immobilize inorganic nitrogen?

The properties a, p, d, B, and CN are used throughout the soil food web literature (Moore and de Ruiter, 2012).

The units of death rate d are \(time^{-1}\), which sets the time step for any simulations. Typically, death rates are included as the inverse of turnover time in years.

The units for biomass B set the unit for the model analyses. These must be in some units of carbon if both carbon and nitrogen analyses are being used. If B is not in units of carbon, then the C:N ratio calculations are nonsensical.

The units should also agree with the units of CN. If C:N ratio is a molar ratio (i.e., \(mol_{Carbon}/mol_{Nitrogen}\)), then the biomass should also be \(mol_{Carbon}\).

The area unit is more flexible and should be defined by the data itself. An example input for biomass data would be \(grams_{Carbon}\) \(hectare^{-1}\).

The conversion efficiency is set by assimilation efficiency a and production efficiency p. Currently, the model assumes that production efficiency for nitrogen is perfect, so production efficiency sets the relative efficiency of carbon versus nitrogen use. These are ratio values that are unitless and should be between 0 and 1.

The DetritusRecycling vector should sum to 1 and is the proportion of dead organic matter either from death or egestion that is returned to the detritus pool. In communities with only 1 detritus pool, this is just a vector with one 1 as shown here. Otherwise, you can split detritus recycling into fractions.

The final properties are Boolean supplied to the model as 1 = TRUE or 0 = FALSE. isDetritus identifies detritus pools. Detrital pools must have a==1, p==1, d==0, and is the only pool type where DetritusRecycling can be non-zero. isPlant identifies plant nodes and allows them to gain carbon at a per capita rate. This feature is poorly developed and may not be reliable. Finally, canIMM identifies trophic species that can immobilize inorganic nitrogen and can therefore have net negative nitrogen mineralization.

Building the interaction matrix:

The interaction matrix can be build in one of two ways. First, it can be built with a data frame of predator-prey interactions and the predator preference for each prey item. This data frame is added to the function build_foodweb, which compiles the interaction matrix and combines it into a list with the properties data frame. Second, you can build the interaction matrix by hand and create the community manually. The second option is described below.

To build a community using the build_foodweb function, you need two data frames: the feeding and the properties. The properties data frame was build in the last section. The feeding data frame always has three columns: Predator, Prey, and Preferences.

Predator and Prey columns contain the names–as they appear in the properties data frame ID column–of all predator-prey pairs.

The Preferences column contains a numerical value weighting the relative feeding preferences of each node after differences in prey abundance are taken into account. If you want predators to feed in proportion to prey abundance, set all of these to 1 (Moore and de Ruiter, 2012). This assumption could be modified by setting the diet ratios using values other than 1. For example, to make the Predator eat Prey1 twice as as Prey2 relative to their biomass you would put a 2 in the row with Predator and Prey1 and a 1 in the row with Predator and Prey2. Note that this does not mean that the predator will eat twice as much of Prey1 as Prey2 unless they have the same biomass.

In the example community, the Predator eats both prey. Prey1 eats both Microbes and Detritus, while Prey2 only eats Detritus. Microbes eat Detritus.

# Create a list of feeding relationships:
feedinglist <- data.frame(Predator = c("Predator", "Predator","Prey1","Prey1", "Prey2", "Microbe"),
           Prey = c("Prey1", "Prey2", "Microbe", "Detritus", "Detritus", "Detritus"),
           Preference = c(1,1,1,1,1,1))

Building the community:

Now you can assemble the community using the build_foodweb function:

# Make the community

our_comm <- build_foodweb(feeding = feedinglist,
                          properties = properties)

# Clean up our working files
rm(feedinglist)

Manual community building:

We can also build the interaction matrix and the community manually.

First, the interaction matrix contains feeding preferences for trophic species on each other. It is formatted with the species names in the row and column names. The rows are the predators and the columns are their prey,so the first row of the matrix below notes that the Predator eats Prey1 and Prey2 in proportion to their biomass.

sp_names = c("Predator", "Prey1", "Prey2", "Microbe", "Detritus")
feedingmatrix = matrix(c(0,1,1,0,0,
                         0,0,0,1,1,
                         0,0,0,0,1,
                         0,0,0,0,1,
                         0,0,0,0,0), 
                       dimnames = list(sp_names, sp_names),
                       ncol = 5, nrow = 5, byrow = TRUE)

To build the community manually, you combine the two features together into a list with the names imat and prop.

# Make the community
our_comm <- list(imat = feedingmatrix,
                 prop = properties)
# Clean up our working files
rm(feedingmatrix, sp_names, properties)

Community analysis (comana)

The comana function is the central function of the package. It analyses the community fluxes and can produce plots of the results. Conceptually, the fluxes are calculated from the top of the food web down, so predation rates at higher trophic levels determine the amount of consumption that is necessary for lower trophic levels to maintain equilibrium (i.e., net zero growth; Moore and de Ruiter, 2012). Technically, the fluxes are calculated by solving the system of linear equations defined by setting each species` change in biomass to zero and calculating their consumption rates. The matrix equation solved is \(\textbf{Ax}=\textbf{b}\), where \(\textbf{A}\) defines feeding interactions, \(\textbf{x}\) is the target vector of total consumption rates for each species, and \(\textbf{b}\) is a vector of loss to natural death. Solving the system using matrix methods at once, instead of sequentially from the top down, allows cannibalism and mutual feeding to be calculated simultaneously.

Inside the comana function, the community is checked for quality using the checkcomm(our_comm). This function adds an extra column to the properties data frame called MutualPred to identify trophic species that are cannibals and eat each other. This is will explored more below.

Analysis of our community produces the following results:

# Analysis
ana1 <- comana(our_comm)

#Returns the following results:

# 1. Consumption rates for each species or input rates into the system for any species/pool that has no prey
ana1$consumption
#>     Predator        Prey1        Prey2      Microbe     Detritus 
#> 3.278689e-01 9.308371e+01 1.945262e+01 3.600201e+04 2.524394e+04

# 2. Carbon mineralized by each species
ana1$Cmin
#>     Predator        Prey1        Prey2      Microbe     Detritus 
#>     0.100000    36.302648     6.127575 25201.407313     0.000000

# 3. Nitrogen mineralization by each species
ana1$Nmin
#>      Predator         Prey1         Prey2       Microbe      Detritus 
#>    0.01880342   -2.98928393   -0.23343141 -150.00837686    0.00000000

# 3.1. Contribution of each prey species to nitrogen mineralized when feeding on each of their prey
ana1$Nminmat
#>          Predator      Prey1       Prey2    Microbe     Detritus
#> Predator        0 0.01196581 0.006837607 0.00000000    0.0000000
#> Prey1           0 0.00000000 0.000000000 0.01633488   -3.0056188
#> Prey2           0 0.00000000 0.000000000 0.00000000   -0.2334314
#> Microbe         0 0.00000000 0.000000000 0.00000000 -150.0083769
#> Detritus        0 0.00000000 0.000000000 0.00000000    0.0000000

# 4. Consumption rates for each species on each of their prey items
ana1$fmat
#>          Predator     Prey1     Prey2   Microbe    Detritus
#> Predator        0 0.2017654 0.1261034 0.0000000     0.00000
#> Prey1           0 0.0000000 0.0000000 0.6031342    92.48058
#> Prey2           0 0.0000000 0.0000000 0.0000000    19.45262
#> Microbe         0 0.0000000 0.0000000 0.0000000 36002.01045
#> Detritus        0 0.0000000 0.0000000 0.0000000     0.00000

# 5. The community returned back
ana1$usin
#> $imat
#>          Predator Prey1 Prey2 Microbe Detritus
#> Predator        0     1     1       0        0
#> Prey1           0     0     0       1        1
#> Prey2           0     0     0       0        1
#> Microbe         0     0     0       0        1
#> Detritus        0     0     0       0        0
#> 
#> $prop
#>         ID   d    a   p        B   CN DetritusRecycling isDetritus isPlant
#> 1 Predator 1.0 0.61 0.5 1.00e-01  4.5                 0          0       0
#> 2    Prey1 3.0 0.65 0.4 8.00e+00  4.8                 0          0       0
#> 3    Prey2 0.5 0.45 0.3 5.00e+00  5.0                 0          0       0
#> 4  Microbe 1.2 1.00 0.3 9.00e+03  8.0                 0          0       0
#> 5 Detritus 0.0 1.00 1.0 1.38e+06 30.0                 1          1       0
#>   canIMM MutualPred
#> 1      0       <NA>
#> 2      0       <NA>
#> 3      0       <NA>
#> 4      1       <NA>
#> 5      0       <NA>

We can plot the results of the community analysis using the plotting built into the function. Options include the web, nitrogen mineralization, and carbon mineralization. The workhorse for the food web plot is the function diagram::plotmat, so check there to understand the options.

# Produce a plot if you want
old.par <- par(no.readonly = TRUE)
par(mfrow = c(2,2), mar = c(4,4,2,2))
ana1 <- comana(our_comm, mkplot = T, BOX.SIZE = 0.15, BOX.CEX = 0.7, edgepos = c(0.15, 0.85))
par(old.par)

Parameter uncertainty

Drawing uncertain parameters from a distribution can be easily accomplished using the function parameter_uncertainty. This function currently allows the user to draw any parameter from either a uniform, normal, or gamma distribution. It requires information on the distribution, some measure of uncertainty (e.g., standard deviation, min/max).

For example, we can draw biomass values from a gamma distribution given the mean values listed in our community and some measure of uncertainty, perhaps from replicate samples (Koltz et al., 2018). For this example we will assume a coefficient of variation of 0.2 for all trophic species. Our example focuses on the direct and indirect effects of each organism on mineralization, whomineralizes as the output of interest. Other options are comana and CNsim.

# Draw parameter uncertainty:

# Build function inputs to indicate distribution, error measure, and error type using an empty matrix.
guidemat <- matrix(NA, nrow = length(our_comm$prop$ID), ncol = 2, dimnames = list(our_comm$prop$ID, c("B","a")))

# Replicate this matrix across the inputs.
distributionin = guidemat
errormeasurein = guidemat
errortypein = guidemat

# Test two most useful distributions (normal often produces negative values).
distributionin[,1] = "gamma"
distributionin[,2] = "uniform"

# Gamma error uses coefficient of variation. Uniform distribution MUST be the minimum (maximum is the parameter value in the input community).
errortypein[,1] = "CV"
errortypein[,2] = "Min"

# Set these values.
errormeasurein[,1] = 0.2
errormeasurein[,2] = 0

# Run the calculation.
oot <- parameter_uncertainty(usin = our_comm, distribution = distributionin, errormeasure = errormeasurein, errortype = errortypein, fcntorun = "whomineralizes", parameters = c("B", "a"), replicates = 10)

# Bind together the result
oot <- do.call("rbind", oot)

# Summarize the direct and indirect effects of Predators across 10 parameter draws:
summary(subset(oot, ID == "Predator")[,-1])
#>     DirectC             DirectN             IndirectC        
#>  Min.   :2.576e-06   Min.   :-2.288e-04   Min.   :7.995e-05  
#>  1st Qu.:3.117e-06   1st Qu.:-1.357e-04   1st Qu.:2.189e-04  
#>  Median :3.880e-06   Median :-1.231e-04   Median :3.745e-04  
#>  Mean   :4.028e-06   Mean   :-1.279e-04   Mean   :4.906e-04  
#>  3rd Qu.:4.258e-06   3rd Qu.:-9.918e-05   3rd Qu.:5.210e-04  
#>  Max.   :7.197e-06   Max.   :-8.110e-05   Max.   :1.729e-03  
#>    IndirectN        
#>  Min.   :4.462e-05  
#>  1st Qu.:1.412e-04  
#>  Median :2.445e-04  
#>  Mean   :3.077e-04  
#>  3rd Qu.:4.060e-04  
#>  Max.   :9.477e-04
# NOTE: You can subset to either include or exclude any trophic species from the result by filtering on the ID column.

Correcting stoichiometry

The community we analyzed above has an important error in stoichiometry. The net nitrogen mineralization for the two prey species is negative even though they cannot immobilize nitrogen. The function corrstoich fixes this by first trying to modify their diet preferences and then reducing their production efficiency to make sure they are getting enough nitrogen (Buchkowski and Lindo, 2021).

our_comm2 <- corrstoich(our_comm)
our_comm2
#> $imat
#>          Predator Prey1 Prey2 Microbe Detritus
#> Predator        0     1     1       0        0
#> Prey1           0     0     0     184        1
#> Prey2           0     0     0       0        1
#> Microbe         0     0     0       0        1
#> Detritus        0     0     0       0        0
#> 
#> $prop
#>         ID   d    a         p        B   CN DetritusRecycling isDetritus
#> 1 Predator 1.0 0.61 0.5000000 1.00e-01  4.5                 0          0
#> 2    Prey1 3.0 0.65 0.4000000 8.00e+00  4.8                 0          0
#> 3    Prey2 0.5 0.45 0.1666667 5.00e+00  5.0                 0          0
#> 4  Microbe 1.2 1.00 0.3000000 9.00e+03  8.0                 0          0
#> 5 Detritus 0.0 1.00 1.0000000 1.38e+06 30.0                 1          1
#>   isPlant canIMM MutualPred
#> 1       0      0       <NA>
#> 2       0      0       <NA>
#> 3       0      0       <NA>
#> 4       0      1       <NA>
#> 5       0      0       <NA>

The feeding preference of Prey1 is modified because they can switch to Microbes over Detritus in this example. Prey2 only has one food source, so it reduces its production efficiency from 0.3 to 0.16. The resulting community has acceptable nitrogen mineralization rates. Because the problem is solved using quadratic programming, small negative numbers are possible from rounding error.

# Print the new nitrogen mineralization rates.
comana(our_comm2)$Nmin
#>      Predator         Prey1         Prey2       Microbe      Detritus 
#>  1.880342e-02 -2.220446e-16  0.000000e+00 -1.507052e+02  0.000000e+00

Sometimes the extreme diet corrections are not reasonable. In this case, a user-defined matrix mirroring imat can be supplied with the maximum diet contributions of each trophic species in the diet. Notice how the leftover correction necessary to balance Prey1’s diet is now accomplished by reducing production efficiency.

# Set diet limits
DL <- our_comm$imat 
# Note: This only works because all feeding is currently set as 1--the same as 100% of the diet can be from that resource if necessary.
DL["Prey1", "Microbe"] = 0.2 # Microbes only 20% of the diet.
corrstoich(our_comm, dietlimits = DL)
#> $imat
#>          Predator Prey1 Prey2  Microbe Detritus
#> Predator        0     1     1  0.00000        0
#> Prey1           0     0     0 38.33333        1
#> Prey2           0     0     0  0.00000        1
#> Microbe         0     0     0  0.00000        1
#> Detritus        0     0     0  0.00000        0
#> 
#> $prop
#>         ID   d    a         p        B   CN DetritusRecycling isDetritus
#> 1 Predator 1.0 0.61 0.5000000 1.00e-01  4.5                 0          0
#> 2    Prey1 3.0 0.65 0.2480000 8.00e+00  4.8                 0          0
#> 3    Prey2 0.5 0.45 0.1666667 5.00e+00  5.0                 0          0
#> 4  Microbe 1.2 1.00 0.3000000 9.00e+03  8.0                 0          0
#> 5 Detritus 0.0 1.00 1.0000000 1.38e+06 30.0                 1          1
#>   isPlant canIMM MutualPred
#> 1       0      0       <NA>
#> 2       0      0       <NA>
#> 3       0      0       <NA>
#> 4       0      1       <NA>
#> 5       0      0       <NA>

Modifying the community structure

The community can be modified to add, combine, or remove trophic species. The function for combining trophic species has the default functionality to combine the two trophic species that have the most similar feeding relationships (i.e., predators and prey), but user-defined inputs are also acceptable (Buchkowski and Lindo, 2021).

Combining species outputs the biomass weighted mean of their parameters and sums their biomass. The combined trophic species are named first/second.

# Combine the most similar trophic species
comtrosp(our_comm)
#> $imat
#>             Predator Prey1/Prey2      Microbe     Detritus
#> Predator           0  0.07692308 0.000000e+00 0.000000e+00
#> Prey1/Prey2        0  0.00000000 3.599712e-07 7.222900e-07
#> Microbe            0  0.00000000 0.000000e+00 7.246377e-07
#> Detritus           0  0.00000000 0.000000e+00 0.000000e+00
#> 
#> $prop
#>            ID        d         a         p        B        CN DetritusRecycling
#> 1    Predator 1.000000 0.6100000 0.5000000 1.00e-01  4.500000                 0
#> 2 Prey1/Prey2 2.038462 0.5730769 0.3615385 1.30e+01  4.876923                 0
#> 4     Microbe 1.200000 1.0000000 0.3000000 9.00e+03  8.000000                 0
#> 5    Detritus 0.000000 1.0000000 1.0000000 1.38e+06 30.000000                 1
#>   isDetritus isPlant canIMM MutualPred
#> 1          0       0      0       <NA>
#> 2          0       0      0       <NA>
#> 4          0       0      1       <NA>
#> 5          1       0      0       <NA>

# Combine the Predator and Prey1
comtrosp(our_comm, selected = c("Predator", "Prey1"))
#> $imat
#>                Predator/Prey1      Prey2      Microbe     Detritus
#> Predator/Prey1      0.0379867 0.03846154 3.599712e-07 3.599712e-07
#> Prey2               0.0000000 0.00000000 0.000000e+00 7.246377e-07
#> Microbe             0.0000000 0.00000000 0.000000e+00 7.246377e-07
#> Detritus            0.0000000 0.00000000 0.000000e+00 0.000000e+00
#> 
#> $prop
#>               ID        d         a         p        B        CN
#> 1 Predator/Prey1 2.975309 0.6495062 0.4012346 8.10e+00  4.796296
#> 3          Prey2 0.500000 0.4500000 0.3000000 5.00e+00  5.000000
#> 4        Microbe 1.200000 1.0000000 0.3000000 9.00e+03  8.000000
#> 5       Detritus 0.000000 1.0000000 1.0000000 1.38e+06 30.000000
#>   DetritusRecycling isDetritus isPlant canIMM     MutualPred
#> 1                 0          0       0      0 Predator/Prey1
#> 3                 0          0       0      0           <NA>
#> 4                 0          0       0      1           <NA>
#> 5                 1          1       0      0           <NA>

The combination of Predator and Prey1 creates two new features that need to be considered. First, it creates a cannibalistic interaction noted in the properties data frame. It also has non-unitary feeding preferences, because these are biomass weighted preferences from the consumption rates of the original prey items. We can leave these default features or remove them using the following options.

# Delete cannibalism and reset feeding preferences
comtrosp(our_comm, selected = c("Predator", "Prey1"), deleteCOMBOcannibal = T, allFEEDING1 = T)
#> $imat
#>                Predator/Prey1 Prey2 Microbe Detritus
#> Predator/Prey1              0     1       1        1
#> Prey2                       0     0       0        1
#> Microbe                     0     0       0        1
#> Detritus                    0     0       0        0
#> 
#> $prop
#>               ID        d         a         p        B        CN
#> 1 Predator/Prey1 2.975309 0.6495062 0.4012346 8.10e+00  4.796296
#> 3          Prey2 0.500000 0.4500000 0.3000000 5.00e+00  5.000000
#> 4        Microbe 1.200000 1.0000000 0.3000000 9.00e+03  8.000000
#> 5       Detritus 0.000000 1.0000000 1.0000000 1.38e+06 30.000000
#>   DetritusRecycling isDetritus isPlant canIMM MutualPred
#> 1                 0          0       0      0       <NA>
#> 3                 0          0       0      0       <NA>
#> 4                 0          0       0      1       <NA>
#> 5                 1          1       0      0       <NA>

You can also remove trophic species from the community by name.

# Delete trophic species
removenodes(our_comm, "Predator")
#> $imat
#>          Prey1 Prey2 Microbe Detritus
#> Prey1        0     0       1        1
#> Prey2        0     0       0        1
#> Microbe      0     0       0        1
#> Detritus     0     0       0        0
#> 
#> $prop
#>         ID   d    a   p       B   CN DetritusRecycling isDetritus isPlant
#> 2    Prey1 3.0 0.65 0.4       8  4.8                 0          0       0
#> 3    Prey2 0.5 0.45 0.3       5  5.0                 0          0       0
#> 4  Microbe 1.2 1.00 0.3    9000  8.0                 0          0       0
#> 5 Detritus 0.0 1.00 1.0 1380000 30.0                 1          1       0
#>   canIMM
#> 2      0
#> 3      0
#> 4      1
#> 5      0

You can also add a new node.

# Add a new trophic species
newnode(our_comm, "NewNode", prey = c(Detritus = 1), predator = c(Predator = 2, Prey1 = 0.1), newprops = c(d = 1, a = 0.1, p = 0.1, B = 10, CN = 10, DetritusRecycling = 0, isDetritus = 0, isPlant = 0, canIMM = 0))
#> $imat
#>          Predator Prey1 Prey2 Microbe Detritus NewNode
#> Predator        0     1     1       0        0     2.0
#> Prey1           0     0     0       1        1     0.1
#> Prey2           0     0     0       0        1     0.0
#> Microbe         0     0     0       0        1     0.0
#> Detritus        0     0     0       0        0     0.0
#> NewNode         0     0     0       0        1     0.0
#> 
#> $prop
#>         ID   d    a   p        B   CN DetritusRecycling isDetritus isPlant
#> 1 Predator 1.0 0.61 0.5 1.00e-01  4.5                 0          0       0
#> 2    Prey1 3.0 0.65 0.4 8.00e+00  4.8                 0          0       0
#> 3    Prey2 0.5 0.45 0.3 5.00e+00  5.0                 0          0       0
#> 4  Microbe 1.2 1.00 0.3 9.00e+03  8.0                 0          0       0
#> 5 Detritus 0.0 1.00 1.0 1.38e+06 30.0                 1          1       0
#> d  NewNode 1.0 0.10 0.1 1.00e+01 10.0                 0          0       0
#>   canIMM MutualPred
#> 1      0       <NA>
#> 2      0       <NA>
#> 3      0       <NA>
#> 4      1       <NA>
#> 5      0       <NA>
#> d      0       <NA>

Calibrating the fluxes to a specific ecosystem

soilfoodwebs calculates the fluxes of carbon and nitrogen based on the community demands, which means that it is very likely that the initial predictions of total mineralization do not match empirical data. This because generic parameter for conversion efficiencies (i.e., a and p) are unlikely to match a given study system.

The user adjust the parameters to deal with this mismatch using two sources: (1) empirical measurements of mineralization or (2) known relationships with abitoic variables.

Using empirical measurements of mineralization

It is not possible to adjust all the food web parameters with overall ecosystem flux data (e.g., carbon mineralization), because there are too many parameters in the model to constrain. Instead, Holtkamp et al. (2011) suggests modifying the microbial production efficiency to match carbon and nitrogen mineralization rates. Modifying only microbial efficiencies is justified because they represent the vast majority of soil heterotrophic respiration and mineralization. We suggest the same strategy.

We will demonstrate this with the intro_comm because it has two microbial pools. We assume the baseline parameters are the correct ones, so we can calculate true carbon and ntirogen mineralization. In reality, these would be measured empirically or estimated for the target ecosystem using available data bases (Jian et al. 2021).

# Create function to modify fungal production efficiency
ccxf <- function(p, ccx){
  
  ccx$prop$p[4] = p[1]
  ccx$prop$p[5] = p[2]
  
  return(c(Cmin = sum(comana(corrstoich(ccx))$Cmin), Nmin = sum(comana(corrstoich(ccx))$Nmin)))
}

# Return carbon and nitrogen mineralization across these gradients
res1 = expand.grid(1:10/10,1:10/10)
res2 = res1

for(i in 1:100){
  res2[i,] = ccxf(as.numeric(res1[i,]), ccx = intro_comm)
}

res = cbind(res1, res2)
colnames(res) = c("p1", "p2", "Cmin", "Nmin")

# Create color gradient to work with
res$p1col = palette.colors(10, "Polychrome 36")[res$p1*10]
res$p2col = palette.colors(10, "Polychrome 36")[res$p2*10]

# Plot the results
old.par <- par(no.readonly = TRUE)
par(mfrow = c(2,2))
plot(Cmin~p1, data = res, col = p2col, xlab = "Prod. eff. Fungi 1", pch = 19)
abline(h = sum(comana(corrstoich(intro_comm))$Cmin), lty = 2)

plot(Nmin~p1, data = res, col = p2col, xlab = "Prod. eff. Fungi 1", pch = 19)
abline(h = sum(comana(corrstoich(intro_comm))$Nmin), lty = 2)

plot(Nmin~Cmin, data = res, col = p2col, bg = p1col, pch = 21)
abline(h = sum(comana(corrstoich(intro_comm))$Nmin), lty = 2)
abline(v = sum(comana(corrstoich(intro_comm))$Cmin), lty = 2)
plot.new()
legend("topright", legend = seq(0.1, 1, by = 0.1), 
       col = palette.colors(10, "Polychrome 36"),
       pch = 19, cex = 0.5, title = "Production efficiency scale", ncol = 2, xpd = T)

par(old.par)

The results above demonstrate that a single measure of carbon mineralization (Cmin) or nitrogen mineralization(Nmin) is not sufficient to fit two food web parameters like microbial production efficiency (x-axis and color in top two figures) because multiple combinations can produce the true value (dashed lines).

When the user has data on both carbon and nitrogen mineralization, two microbial production efficiencies can be resolved accurately (bottom figure; color and fill show the two production efficiencies).

In a applied situation, optimization, for example stats::optim, could be used to match the parameters and the data. Gauzens et al. 2019 provide excellent examples for how to use metabolic scaling relationships to fit loss rates in their R package fluxweb. Their package places metabolic losses in a different location than soilfoodwebs, but the user could predict metabolic losses using body size and then adjust production efficiencies accordingly to get the desired carbon mineralization rates.

More detailed relationships between abiotic variables and microbial assimilation and production efficiency are available. These can be helpful because the relationship between body size and metabolism are less useful for microorganisms.

Using defined relationships

Another option is to use defined relationships between respiration (Manzoni et al. 2012). For example, Xu et al. (2014) provide some relationships for microbial death rate, respiration rate, and production efficiency across temperature, moisture, and substrate quality gradients. We implement these for the community here to demonstrate how to adjust microbial death rate and conversion efficiency.


temp_moist_microbe <- function(temp, moist, comm){
  
  # Set functional range:
  if(temp < 0 | moist < 0.05) stop("Function not coded for these conditions.")
  
  # Functions and parameters from Xu et al. 2014
  fmm = ifelse(moist < 0.6,log10(0.05/moist)/log10(0.05/moist),1)
  
  fmt = 2.5^((temp - 25)/10)

  # New microbial death rate
  comm$prop[comm$prop$ID == "Microbe", "d"] = 
    comm$prop[comm$prop$ID == "Microbe", "d"]*fmm*fmt
  
  # New microbial assimilation efficiency
  comm$prop[comm$prop$ID == "Microbe", "a"] = 
    (0.43 - (-0.012*(temp-15)))*
    (comm$prop[comm$prop$ID == "Microbe", "CN"]/comm$prop[comm$prop$ID == "Detritus", "CN"])^0.6
  
  # Ignoring changes in maintenance respiration here. See the cited paper for the details.

  return(comm)
}

# Explore the carbon mineralization rate across temperature as an example

res <- data.frame(temp = seq(1,30, length = 10))

res$Cmin = NA

for(i in 1:10){
  res$Cmin[i] = 
    sum( # Sum all C min rates
      comana( # Calculate Cmin
        corrstoich( # Correct stoichiometry
      temp_moist_microbe(res[i, "temp"], 0.7, our_comm)
      )
      )$Cmin
      )
}

plot(Cmin~temp, data = res, xlab = "Temperature", ylab = "C. mineralization", type = "l", lwd = 2)

Contributions to carbon and nitrogen mineralization

A core feature of soilfoodwebs is to calculate the contribution of each trophic species to carbon and nitrogen mineralization. This is achieved using the function whomineralizes. The direct contribution sums up each trophic species’ carbon and nitrogen mineralization from comana, while their indirect contribution calculates the change in fluxes without that trophic species and removes their direct contributions. Similar calculations are outlined by Holtkamp et al. (2011), but our version of them produces a different result. Specifically, our calculation uses the following formula for indirect effects \((IE_X)\) on mineralization of element \(X\):

\[IE_X = (X_{with} - X_{without}- X_{mineralization})/X_{with} \]

where \(X_{with}\) is the total mineralization with the trophic species,\(X_{without}\) is the mineralization without the trophic species, and \(X_{mineralization}\) is the mineralization coming directly from the trophic species (e.g., respiration rate).

# Calculate the mineralization contributions
whomineralizes(our_comm)
#>         ID      DirectC       DirectN     IndirectC     IndirectN
#> 1 Predator 3.961347e-06 -0.0001227279  2.410963e-05  2.362742e-04
#> 2    Prey1 1.438074e-03  0.0195107322  3.709909e-05 -6.907725e-05
#> 3    Prey2 2.427345e-04  0.0015235815 -7.783569e-06 -9.776173e-05
#> 4  Microbe 9.983152e-01  0.9790884142 -8.584855e-17 -2.345553e-04
#> 5 Detritus 0.000000e+00  0.0000000000 -8.548112e-03  1.016577e+00

Remember that the units of these fluxes are set by the units of biomass B and death rate d.

Calculate inputs

You can also calculate the inputs and outputs of the food web using the function calculate_inputs. This function saves a data frame, but also prints the results on the console.

inout <- calculate_inputs(our_comm)
#> [1] "------------------------------------------------"
#> [1] "Equilibrium is maintained by these inputs:"
#> [1] "------------------------------------------------"
#> [1] "The system gains organic compounds from the detritus pool:"
#> [1] "Gains of C from detritus:25368"
#> [1] "Gains of N from detritus:-151"
#> [1] "Gains detritus at a C:N ratio of:-168"
#> [1] "------------------------------------------------"
#> [1] "The web mineralizes 25368 units of carbon"
#> [1] "The web must immobilize 151 units of nitrogen"
inout
#> $IN
#> Predator    Prey1    Prey2  Microbe Detritus 
#>        0        0        0        0        0 
#> 
#> $r_i
#> Predator    Prey1    Prey2  Microbe Detritus 
#>        0        0        0        0        0 
#> 
#> $Cmin
#>    Predator       Prey1       Prey2     Microbe    Detritus 
#>     0.10000    36.30265    13.13052 25318.47018     0.00000 
#> 
#> $Nmin
#>      Predator         Prey1         Prey2       Microbe      Detritus 
#>    0.01880342    0.00000000    0.00000000 -150.70517965    0.00000000

Notice that in this example the equilibrium can only be maintained by a gain of carbon and loss of nitrogen. This is because the dead microbes and animals are high nitrogen items that must be lost from the community.

We can add a necromass pool with the right C:N ratio to fix this problem if we don’t want to mix our high C:N detrital inputs with our low C:N ratio detrital recycling.

# Add the new node
our_comm_necro = newnode(our_comm, "Necromass", prey = NA, predator = c(Prey1 = 1, Prey2 = 1, Microbe = 1), newprops = c(d = 0, a = 1, p = 1, B = 100, CN = 5, DetritusRecycling = 1, isDetritus = 1, isPlant = 0, canIMM = 0))
#> Warning in checkcomm(COMM): Rescaling Detritus Recycling to sum to 1.
# Modify the DetritusRecycling of the old detritus pool
our_comm_necro$prop$DetritusRecycling = c(0,0,0,0,0,1)

# Check out the community:
our_comm_necro
#> $imat
#>           Predator Prey1 Prey2 Microbe Detritus Necromass
#> Predator         0     1     1       0        0         0
#> Prey1            0     0     0       1        1         1
#> Prey2            0     0     0       0        1         1
#> Microbe          0     0     0       0        1         1
#> Detritus         0     0     0       0        0         0
#> Necromass        0     0     0       0        0         0
#> 
#> $prop
#>          ID   d    a   p        B   CN DetritusRecycling isDetritus isPlant
#> 1  Predator 1.0 0.61 0.5 1.00e-01  4.5                 0          0       0
#> 2     Prey1 3.0 0.65 0.4 8.00e+00  4.8                 0          0       0
#> 3     Prey2 0.5 0.45 0.3 5.00e+00  5.0                 0          0       0
#> 4   Microbe 1.2 1.00 0.3 9.00e+03  8.0                 0          0       0
#> 5  Detritus 0.0 1.00 1.0 1.38e+06 30.0                 0          1       0
#> d Necromass 0.0 1.00 1.0 1.00e+02  5.0                 1          1       0
#>   canIMM MutualPred
#> 1      0       <NA>
#> 2      0       <NA>
#> 3      0       <NA>
#> 4      1       <NA>
#> 5      0       <NA>
#> d      0       <NA>

# New inputs
inout2 = calculate_inputs(our_comm_necro)
#> [1] "------------------------------------------------"
#> [1] "Equilibrium is maintained by these inputs:"
#> [1] "36086 units of C to Detritus"
#> [1] "1203 units of N to Detritus"
#> [1] "------------------------------------------------"
#> [1] "The system looses organic compounds from the detritus pool:"
#> [1] "Looses of C from detritus:10838"
#> [1] "Looses of N from detritus:1352"
#> [1] "Looses detritus at a C:N ratio of:8"
#> [1] "------------------------------------------------"
#> [1] "The web mineralizes 25248 units of carbon"
#> [1] "The web must immobilize 150 units of nitrogen"

Stability

The stability of the community at equilibrium can be calculated using an eigen decomposition of the Jacobian matrix. Several features are available:

Basic analysis of the Jacobian.
Analysis of the Jacobian adding density-dependent population regulation to selected trophic species.
Soil food web model also use a modified stability calculation, where smin is calculated. smin is the value multiplied by the diagonal terms of the Jacobian necessary to make the food web stable (de Ruiter et al. 1995).

Two stability functions return these calculations. stability calculates the Jacobian or smin algebraically using the structure of soil food web models. stability2 estimates the Jacobian using the function rootSolve::jacobian.full.

stab1 <- stability(our_comm)

# The Jacobian:
stab1$J
#>               Predator        Prey1         Prey2       Microbe      Detritus
#> Predator  5.768275e-12  0.007692308  7.692308e-03  0.000000e+00  0.000000e+00
#> Prey1    -2.017654e+00  0.000000000  0.000000e+00  1.466774e-03  7.971596e-06
#> Prey2    -1.261034e+00  0.000000000 -1.276756e-14  0.000000e+00  1.902973e-06
#> Microbe   0.000000e+00 -6.346616813  0.000000e+00  2.220446e-16  7.862879e-03
#> Detritus  2.278689e+00  1.793401176 -2.641488e+00 -2.816830e+00 -2.624094e-02

# The eigen decomposition:
stab1$eigen
#> eigen() decomposition
#> $values
#> [1] -0.007044524+0.2010952i -0.007044524-0.2010952i -0.005289296+0.1263491i
#> [4] -0.005289296-0.1263491i -0.001573304+0.0000000i
#> 
#> $vectors
#>                            [,1]                       [,2]
#> [1,]  9.857894e-06+5.10678e-05i  9.857894e-06-5.10678e-05i
#> [2,] -1.026049e-03+1.47423e-04i -1.026049e-03-1.47423e-04i
#> [3,] -3.180124e-04+6.35186e-05i -3.180124e-04-6.35186e-05i
#> [4,] -7.144527e-03-7.11323e-02i -7.144527e-03+7.11323e-02i
#> [5,]  9.974407e-01+0.00000e+00i  9.974407e-01+0.00000e+00i
#>                             [,3]                        [,4]             [,5]
#> [1,]  6.505374e-06-5.351043e-05i  6.505374e-06+5.351043e-05i -3.452614e-08+0i
#> [2,]  3.392426e-04+1.161714e-04i  3.392426e-04-1.161714e-04i  1.237190e-03+0i
#> [3,]  5.352138e-04+2.747616e-05i  5.352138e-04-2.747616e-05i -1.237183e-03+0i
#> [4,] -7.710979e-03-4.480378e-02i -7.710979e-03+4.480378e-02i -6.809188e-03+0i
#> [5,]  9.989658e-01+0.000000e+00i  9.989658e-01+0.000000e+00i  9.999753e-01+0i

# The largest eignvalue:
stab1$rmax # The community is stable if this is negative.
#> [1] -0.001573304

# Use the function calc_smin to recover smin:
calc_smin(our_comm)
#> [1] 0.004

# Add density-dependence using stability2:
stability2(our_comm, densitydependence = c(1,1,1,0,0))
#> $J
#>      PredatorC       Prey1C       Prey2C      MicrobeC     DetritusC
#> [1,] -1.000000  0.007692308  0.007692308  0.000000e+00  0.000000e+00
#> [2,] -2.017654 -2.999999937  0.000000000  1.466774e-03  7.971596e-06
#> [3,] -1.261034  0.000000000 -0.500000006  0.000000e+00  1.902973e-06
#> [4,]  0.000000 -6.346613191  0.000000000  2.021099e-08  7.846155e-03
#> [5,]  3.278728  4.793355401 -2.141496225 -2.816830e+00 -2.622422e-02
#> 
#> $eigen
#> eigen() decomposition
#> $values
#> [1] -2.9890965+0.0000000i -0.9878663+0.0000000i -0.5199380+0.0000000i
#> [4] -0.0146618+0.1480693i -0.0146618-0.1480693i
#> 
#> $vectors
#>                 [,1]           [,2]            [,3]                        [,4]
#> [1,] -0.001621502+0i  0.07500995+0i  0.004734199+0i  6.243183e-09+2.040386e-07i
#> [2,]  0.420113768+0i -0.07556810+0i -0.003884991+0i  6.170475e-07+2.564393e-05i
#> [3,] -0.000821620+0i  0.19388862+0i  0.299337225+0i -3.744874e-06+6.123591e-07i
#> [4,]  0.891564517+0i -0.47874413+0i -0.061790166+0i  4.102973e-03+5.253630e-02i
#> [5,]  0.169156230+0i -0.84963182+0i  0.952124832+0i -9.986106e-01+0.000000e+00i
#>                             [,5]
#> [1,]  6.243183e-09-2.040386e-07i
#> [2,]  6.170475e-07-2.564393e-05i
#> [3,] -3.744874e-06-6.123591e-07i
#> [4,]  4.102973e-03-5.253630e-02i
#> [5,] -9.986106e-01+0.000000e+00i
#> 
#> 
#> $rmax
#> [1] -0.01466177

# stability and stability2 should produce similar results when used in default:
stability(our_comm)$rmax
#> [1] -0.001573304
stability2(our_comm)$rmax
#> [1] -0.001575481

Simulate the food web away from equilibrium

The function CNsim wraps a number of underlying functions together to retrieve the model parameters and simulate the model (de Ruiter et al. 1994). Feeding occurs using Type-I (i.e., linear) functional responses, while density-dependence is user defined. Because the model is at equilibrium to start, this function is interesting when the starting conditions are modified or it is used for a detritus simulation experiment.

# Simulate equilibirum over 100 time steps
sim1 <- CNsim(our_comm)
# Modify predator biomass to 80% of equilibrium value
sim2 <- CNsim(our_comm, start_mod = c(0.8,1,1,1,1))

# Modify predator biomass to 80% of equilibrium value with density-dependence.
sim3 <- CNsim(our_comm, start_mod = c(0.8,1,1,1,1), densitydependence = c(1,0,0,0,0))

# Plot the results:
old.par <- par(no.readonly = TRUE)
par(mfrow= c(2,2), mar = c(5,4,2,2))
plot(Predator_Carbon~Day, data = sim2, type = "l", col = "blue")
points(Predator_Carbon~Day, data = sim3, type = "l", col = "orange")
points(Predator_Carbon~Day, data = sim1, type = "l")

plot(Prey1_Carbon~Day, data = sim2, type = "l", col = "blue")
points(Prey1_Carbon~Day, data = sim3, type = "l", col = "orange")
points(Prey1_Carbon~Day, data = sim1, type = "l")

plot(Prey2_Carbon~Day, data = sim2, type = "l", col = "blue")
points(Prey2_Carbon~Day, data = sim3, type = "l", col = "orange")
points(Prey2_Carbon~Day, data = sim1, type = "l")

plot(Microbe_Carbon~Day, data = sim2, type = "l", col = "blue")
points(Microbe_Carbon~Day, data = sim3, type = "l", col = "orange")
points(Microbe_Carbon~Day, data = sim1, type = "l")
legend("bottomright", legend = c("Base", "80% predator", "w/ density-dependence"), col = c("black", "blue", "orange"), lty = 1)

par(old.par)

Decomposition constant

The model can calculate the decomposition constant and the effect of each soil organism on the decomposition rate. The function `decompexpt can calculate decomposition constants (k) with the option overtime enabled to return a vector of decomposition.

Besides returning the decomposition rate, the function returns a table of the effect of each trophic species on the decomposition constant. The direct effects set each species consumption to zero, but do not recalculate the food web consumption rates. The indirect effects recalculate the consumption rates and subtract the direct effects.

Overtime produces a list of simulations showing how each trophic species affects overall decomposition rate over time. These results can be easily plotted. The user must select the detritus pool, because the list is multilayered when there is more than one detritus pool.

# The function decompexpt 
decompres = decompexpt(our_comm, overtime = 10)

# The decomposition constants for each detritus pool identified in the community
decompres$basedecomp
#>   Detritus 
#> 0.02616952

# A table of the direct and indirect effects 
decompres$decompeffects
#>         ID DetritusID       Direct      Indirect
#> 1 Predator   Detritus 0.0000000000  4.767834e-05
#> 2    Prey1   Detritus 0.0025607998  1.428502e-05
#> 3    Prey2   Detritus 0.0005386456 -1.363312e-05
#> 4  Microbe   Detritus 0.9969005546 -1.670087e-05

# Plot the decomposition with and without Microbe
decompdata = decompres$overtime$Detritus
plot(Original~Day, data = decompdata, type = "l", ylim = c(0.5,1))
points(Microbe~Day, data = decompdata, col = "red", type = "l")
legend("bottomleft", legend = c("Original", "No Microbe"), col = c("black", "red"), lty = 1)

Decomposition experiments in silico

The model can be used to simulate the decomposition of detritus in an in silico experiment wherein the surrounding food web is not affected. Furthermore, the surrounding food web can be held at equilibrium or modified during the experiment. The inspiration for this functionality comes from the work of Zheng et al. (1997).

# Simulate detritus experiment at equilibirum over 100 time steps.
# NOTE: This is equivalent to the overtime option presented for decompexpt for the original community although numerical errors in the solver makes it diverge slightly.
sim1 <- CNsim(our_comm, DETEXPT = "Detritus", TIMES = 1:50)
#> Warning in CNsim(our_comm, DETEXPT = "Detritus", TIMES = 1:50): Making the
#> starting detritus experiment 100 units of C.

# Modify microbial biomass to 80% of equilibrium value
sim2 <- CNsim(our_comm, start_mod = c(1,1,1,0.5,1), DETEXPT = "Detritus", TIMES = 1:50)
#> Warning in CNsim(our_comm, start_mod = c(1, 1, 1, 0.5, 1), DETEXPT =
#> "Detritus", : Making the starting detritus experiment 100 units of C.

# Modify microbial biomass to 80% of equilibrium value with density-dependence.
sim3 <- CNsim(our_comm, DETEXPT = "Detritus", start_mod = c(1,1,1,0.5,1), densitydependence = c(0,0,0,1,0), TIMES = 1:50)
#> Warning in CNsim(our_comm, DETEXPT = "Detritus", start_mod = c(1, 1, 1, :
#> Making the starting detritus experiment 100 units of C.

# Plot the results:
plot(DetExpt~Day, data = sim2, type = "l", col = "blue")
points(DetExpt~Day, data = sim3, type = "l", col = "orange")
points(DetExpt~Day, data = sim1, type = "l")
legend("topright", legend = c("Base", "80% microbial biomass", "w/ density-dependence"), col = c("black", "blue", "orange"), lty = 1)

References

Buchkowski, R. W., and Z. Lindo. 2021. Stoichiometric and structural uncertainty in soil food web models. Functional Ecology 35:288–300.

Gauzens, B., Barnes, A., Giling, D. P., Hines, J., Jochum, M., Lefcheck, J. S., Rosenbaum, B., Wang, S., & Brose, U. (2019). fluxweb: An R package to easily estimate energy fluxes in food webs. Methods in Ecology and Evolution, 10(2), 270–279.

Holtkamp, R., A. van der Wal, P. Kardol, W. H. van der Putten, P. C. de Ruiter, and S. C. Dekker. 2011. Modelling C and N mineralisation in soil food webs during secondary succession on ex-arable land. Soil Biology & Biochemistry 43:251–260.

Jian, J., Vargas, R., Anderson-Teixeira, K., Stell, E., Herrmann, V., Horn, M., Kholod, N., Manzon, J., Marchesi, R., Paredes, D., & Bond-Lamberty, B. (2021). A restructured and updated global soil respiration database (SRDB-V5). Earth System Science Data, 13(2), 255–267.

Koltz, A. M., A. Asmus, L. Gough, Y. Pressler, and J. C. Moore. 2018. The detritus-based microbial-invertebrate food web contributes disproportionately to carbon and nitrogen cycling in the Arctic. Polar Biology 41:1531–1545.

Manzoni, S., Taylor, P., Richter, A., Porporato, A., & Ågren, G. I. (2012). Environmental and stoichiometric controls on microbial carbon-use efficiency in soils: Research review. New Phytologist, 196(1), 79–91.

Moore, J.C. and de Ruiter, P. C. 2012. Energetic Food Webs: An Analysis of Real and Model Ecosystems. Oxford University Press, Oxford.

de Ruiter, P. C., A.-M. Neutel, and J. C. Moore. 1994. Modelling food webs and nutrient cycling in agro-ecosystems. Trends in Ecology & Evolution 9:378–383.

de Ruiter, P. C., A.-M. Neutel, and J. C. Moore. 1995. Energetics, Patterns of Interaction Strengths, and Stability in Real Ecosystems. Science 269:1257–1260.

Xu, X., Schimel, J. P., Thornton, P. E., Song, X., Yuan, F., & Goswami, S. (2014). Substrate and environmental controls on microbial assimilation of soil organic carbon: A framework for Earth system models. Ecology Letters, 17(5), 547–555.

Zheng, D. W., J. Bengtsson, and G. I. Ågren. 1997. Soil food webs and ecosystem processes: decomposition in donor-control and Lotka-Volterra systems. American Naturalist 149:125–148.

soilfoodweb_vignette

Robert W. Buchkowski

2023-05-08