In FluvialLandscapeLab/gangsta: The Generalized Algorithm for Nutrient, Growth, Stoichiometric and Thermodynamic Analysis

gangsta basics

The Generalized Algorithm for Nutrient, Growth, Stoichiometric and Thermodynamic Analysis (gangsta) formalizes structured, user-defined conceptual models of linked biogeochemical cycles via instantiation in constraint-based computer code. This vignette describes how to generate a gangsta-derived model of a single time step using the gangsta R package.

Citing this package

Poole, G.C., and Reinhold, A.M. 2018. gangsta v1.0: Generalized algorithm for nutrient, growth, stoichiometric, and thermodynamic analysis. https://github.com/FluvialLandscapeLab/gangstaBuiltPackage

Detailed description of the gangsta and some applications

Reinhold, A.M., G.C. Poole, C. Izurieta, A.M. Helton, and E.S. Bernhardt. 2019. Constraint-based simulation of multiple interactive elemental cycles in biogeochemical systems. Ecological Informatics. 50:102-121. DOI: 10.1016/j.ecoinf.2018.12.008

Formalizing and instantiating a conceptual model

Formalize the conceptual model

In the example in this vignette, we will instantiate the following structured conceptual model:

• Optimization function: maximize biomass of all organisms during the time step
• Tracked elements: carbon (C), nitrogen (N), and oxygen (O)
• Organism set: the two types of organisms (chemoautotrophs [Aut] and methanotrophs [Met]); the processes that they are capable of carrying out (nitrification, methane oxidation, decay, and assimilation of CO$_2$, CH$_4$, NH$_4$^+^, and NO$_3$^-^; and the compounds involved in those processes (CO$_2$, CH$_4$, NH$_4$^+^, NO$_3$^-^, O$_2$, dissolved organic matter [DOM], and Ox [a source/sink for O atoms])
• Biogeochemical constraints: mathematical expressions describing stoichiometry, thermodynamics, molar and energy balances and the interrelationships among each of these

The conceptual model in this vignette is not intended to represent any "real" system. It is simply a demonstration of the R package.

gangsta v1.0 can accept user-defined tracked elements and organism sets. The optimization function and biogeochemical constraints are fixed in this version.

First steps in R

Load the package.

library(gangsta)

For the sake of convenience, we will generate objects to store the values for the model parameters describing:

• the ecological stoichiometry of biomass and dissolved organic matter (DOM). We assume the chemical composition of both biomass and DOM is CN${16/106}$O${110/106}$.

EcoStoichC = 106/106
EcoStoichN = 16/106
EcoStoichO = 110/106

• the respiration rate of the organisms. Here, the units of respiration rate are kJ per $\mu$mol of biomass per hour.

resp = -2.83E-6

• the length of the time step. Here, the units are hours.

timeStep = 24

Importantly, gangsta models use variables where the units are derived using only three base units: number of atoms and molecules (e.g., mol, $\mu$mol), energy (e.g., J or KJ), and time (e.g., s, hr, d, fortnight). The units can be selected by the end-user, so long as all arguments passed to functions in the gangsta package are consistent.

Next, we will begin to create our gangsta objects and store them in an object called myGangstas.

Create gangsta objects describing the tracked elements and organism set

Using the compoundFactory() function, we will first create gangsta objects for the compounds that are neither organisms nor source/sinks. Each compound must have a name (compoundName), a description of the elemental pools contained in the compound (molarRatios), and the initial amount available at the start of the model run (initialMolecules).

The molarRatios describe the chemical composition of the tracked elements contained in compounds. For example, the molarRatios for the elemental pools contained in the CO$_2$ compound are 1 for C and 2 for O because each CO$_2$ molecule contains one atom of C and 2 atoms of O; both C and O are tracked elements. In the cases of NH$_4$^+^ or CH$_4$ compounds, no molar ratio is specified for the hydrogen contained in the NH$_4$^+^ or CH$_4$ compounds because H is not a tracked element in the conceptual model.

The molarRatios serve two important purposes:

• The gangsta reads the tracked elements from the specification of the molarRatios for the elements contained in the complete set of compounds specified by the end user. That is, whichever elements have specified molarRatios are those tracked in the gangsta model.
• The elemental pools associated with the compounds are created from the molarRatios. The names of the pools are created from the names of the objects in the vector, and the ratios are stored in the vector.

The initialMolecules provide the starting values for the compounds at the beginning of the model. The units of initialMolecules are $\mu$mols.

myGangstas =
  c(
    compoundFactory(
      compoundName = "CO2",
      molarRatios = c(C=1, O=2),
      initialMolecules = 0.2
    ),
    compoundFactory(
      compoundName = "CH4" ,
      molarRatios = c(C=1),
      initialMolecules = 0.1
      ),
    compoundFactory(
      compoundName = "NH4",
      molarRatios = c(N=1),
      initialMolecules = 0.3
    ),
    compoundFactory(
      compoundName = "NO3",
      molarRatios = c(N=1, O=3),
      initialMolecules = 0.1
      ),
    compoundFactory(
      compoundName = "O2",
      molarRatios = c(O=2),
      initialMolecules = 0.9
      ),
    compoundFactory(
      compoundName = "DOM",
      molarRatios = c(C=EcoStoichC, N=EcoStoichN, O=EcoStoichO),
      initialMolecules = 0.3
    )
  )

Compounds that are to be modeled as source/sinks should be specified as infiniteCompounds. The initialMolecules for these compounds must always be set to 0.

Next, we will add an infinite compound to myGangstas called "Ox." Ox will serve as a source/sink for O atoms.

myGangstas =
  c(myGangstas,
    compoundFactory(
      compoundName = "Ox",
      molarRatios = c(O=1),
      initialMolecules = 0,
      infiniteCompound = T
      )
  )

The next gangsta objects that we will create are the organisms in our organism set. We will use the compoundFactory() function to do so, because organisms are special types of finite compounds. Organisms have all the properties of finite compounds, and in addition, respire and perform metabolic processes.

Importantly, specifying a respirationRate for a compound results in the gangsta assigning the compound to be of class organism.

The following code creates gangsta objects for the two types of organisms, chemoautotrophs (Aut) and methanotrophs (Met), in our conceptual model.

myGangstas =
  c(myGangstas,
    compoundFactory(
      compoundName = "Aut",
      molarRatios = c(C=EcoStoichC, N=EcoStoichN, O=EcoStoichO),
      initialMolecules = 1,
      respirationRate = resp * timeStep
    ),
    compoundFactory(
      compoundName = "Met",
      molarRatios = c(C=EcoStoichC, N=EcoStoichN, O=EcoStoichO),
      initialMolecules = 1,
      respirationRate = resp * timeStep
    )
  )

Next, we will use processFactory() to create gangsta objects for the processes carried out by organisms. Doing so requires specifying the following:

• the name of the process (i.e., processName);
• the organism carrying out the process (i.e, organismName);
• the chemical affinity associated with the process (i.e., energyTerm; herein, units are KJ);
• the names of the reactant compounds (i.e., fromCompoundNames) and product compounds (toCompoundNames), and the pools of tracked elements associated with those compounds;
• the $\mu$mols of atoms transferred from the reactant compounds to the product compounds (molarTerms);
• the indicies describing the to pools involved in each transfer involved in the process; and
• whether or not to constrain the availability of molecules of compound to those present at the start of the time step.

This subsequent block of code creates gangsta objects for the two dissimilatory processes in our conceptual model and stores these objects in myGangstas.

The energyTerm in our conceptual model is the chemical affinity of the process, which is equal in magnitude but opposite in sign of the Gibbs free energy yield. From the perspective of an organism, a positive energyTerm provides the organism with energy and a negative energyTerm expends energy. Thus, each dissimilatory process has a positive energyTerm.

When calling processFactory(), the names of the objects within the lists inputted to fromCompoundNames, toCompoundNames, and molarTerms are identical and their order must be identical within each of these lists. These lists describe the stoichiometry of the process. For example, in the process of nitrification (AutNitrif) below,

• 1 $\mu$mol of N atoms are transferred from NH$_4$^+^-N to NO$_3$^-^-N,
• 3 $\mu$mol of O atoms are transferred from O$_2$-O to NO$_3$^-^-O, and
• 1 $\mu$mol of O atoms are transferred from O$_2$-O to Ox-O.

myGangstas =
  c(
    myGangstas,
    processFactory(
      myGangstas,
      processName = "AutNitrif",
      energyTerm = 3.485E-4, 
      fromCompoundNames = list(N = "NH4", O = "O2", O = "O2"),
      toCompoundNames = list(N = "NO3", O = "NO3", O = "Ox"),
      molarTerms = list(N = 1, O = 3, O = 1),
      organismName = "Aut"
    ),
    processFactory(
      myGangstas,
      processName = "MetMethaneOxid",
      energyTerm = 8.18E-4, 
      fromCompoundNames = list(C = "CH4", O = "O2", O = "O2"),
      toCompoundNames = list(C = "CO2", O = "CO2", O = "Ox"),
      molarTerms = list(C = 1, O = 2, O = 2),
      organismName = "Met"
    )
  )

This next block of code creates objects describing the assimilatory processes in our conceptual model. Note that the energyTerms are negative for these processes because, from the perspective of an organism, the energy used for assimilation is subtracted from the energy available to an organism to meet its respiratory demands.

The gangsta assumes that all elemental pools containing tracked elements are available for assimilation. However, excess $\mu$mols of atoms are not assimilated. For example, the specification below for autotrophic assimilation of CO$_2$ (AutAssimCO2) indicates that all of the C from the CO$_2$ will be assimilated and as much O from CO$_2$ will be assimilated as needed, but that any excess O will be transferred to Ox. This is handled via molarTerms and transferOptions where the molarTerms = list(C = 1, O = 2, O = 2) and the transferOptions = list(C = 1, O = 2:3). For C, understanding this specification is straightforward; the carbon transfer is from the CO$_2$-C from pool to the Aut-C to pool. However, the O transfer is a bit more complex. The molarTerms for the from pool CO$_2$-O are repeated twice O = 2, O = 2, but there is only one element in the transferOptions list for O where O = 2:3; this tells the gangsta that the two molarTerms for O represent one from pool (CO$_2$-O) and need need to be split between two to pools (Aut-O and Ox-O).

myGangstas =
  c(
    myGangstas,
    processFactory(
      myGangstas,
      processName = "AutAssimCO2",
      energyTerm = -3.5E-03,
      fromCompoundNames = list(C = "CO2", O = "CO2", O = "CO2"),
      toCompoundNames = list(C = "Aut", O = "Aut", O = "Ox"),
      molarTerms = list(C = 1, O = 2, O = 2),
      transferOptions = list(C = 1, O = 2:3),
      organismName = "Aut"
    ),
    processFactory(
      myGangstas,
      processName = "MetAssimCH4",
      energyTerm = -1.09E-03,
      fromCompoundNames = list(C = "CH4"),
      toCompoundNames = list(C = "Met"),
      molarTerms = list(C = 1),
      organismName = "Met"
    ),
    processFactory(
      myGangstas,
      processName = "AutAssimNO3",
      energyTerm = -1.55E-04,
      fromCompoundNames = list(N = "NO3", O = "NO3", O = "NO3"),
      toCompoundNames = list(N = "Aut", O = "Aut", O = "Ox"),
      molarTerms = list(N = 1, O = 3, O = 3),
      transferOptions = list(N = 1, O = 2:3),
      organismName = "Aut"
    ),
    processFactory(
      myGangstas,
      processName = "MetAssimNO3",
      energyTerm = -1.55E-04,
      fromCompoundNames = list(N = "NO3", O = "NO3", O = "NO3"),
      toCompoundNames = list(N = "Met", O = "Met", O= "Ox"),
      molarTerms = list(N = 1, O = 3, O = 3),
      transferOptions = list(N = 1, O = 2:3),
      organismName = "Met"
    ),
    processFactory(
      myGangstas,
      processName = "AutAssimNH4",
      energyTerm = -3.18E-05,
      fromCompoundNames = list(N = "NH4"),
      toCompoundNames = list(N = "Aut"),
      molarTerms = list(N = 1),
      organismName = "Aut"
    ),
    processFactory(
      myGangstas,
      processName = "MetAssimNH4",
      energyTerm = -3.18E-05,
      fromCompoundNames = list(N = "NH4"),
      toCompoundNames = list(N = "Met"),
      molarTerms = list(N = 1),
      organismName = "Met"
    )
  )

Below, we specify the processes of decay for our two types of organisms. In both cases, we assume that decayed organisms become DOM.

myGangstas =
  c(myGangstas,
    processFactory(
      myGangstas,
      processName = "AutDecay",
      energyTerm = 0,
      fromCompoundNames = list(C = "Aut", N = "Aut", O = "Aut"),
      toCompoundNames = list(C = "DOM", N = "DOM", O = "DOM"),
      molarTerms = list(C = 1, N = EcoStoichN, O = EcoStoichO),
      organismName = c("Aut")
    ),
    processFactory(
      myGangstas,
      processName = "MetDecay",
      energyTerm = 0,
      fromCompoundNames = list(C = "Met", N = "Met", O = "Met"),
      toCompoundNames = list(C = "DOM", N = "DOM", O = "DOM"),
      molarTerms = list(C = 1, N = EcoStoichN, O = EcoStoichO),
      organismName = c("Met")
    )
  )

At this point, our conceptual model is completely specified.

Incorporate the goal function and biogeochemical constraints, and generate simulation code

The next line of code generates the gangsta simulation model code. The following are incorporated into the code: the mathematical expressions describing the goal function, tracked, elements, organism set, and biogeochemical constraints associated with the structured conceptual model.

The gangsta simulation model code is written to a file that can be run in lpSolve (http://lpsolve.sourceforge.net/5.5/). Files should be saved with a ".lp" file extension. Executing the next line of code allows the user to either select a .lp file to overwrite or to create a new .lp file. An R interface to lpSolve is contained within the lpSolveAPI package.

writeGangstaModel(gangstaObjects = myGangstas, file = file.choose())

Running lpSolve models from R

The lpSolve model generated by the gangsta can be executed as follows.

Install and load the lpSolveAPI package.

install.packages("lpSolveAPI")
library(lpSolveAPI)

Create the lpSolve model object.

lpModel = read.lp(file.choose(), verbose = "normal")

Solve the lpModel.

solve(lpModel)

The solve function returns a status code and, as a side effect, updates the lpModel object to include the model solution.

Get the results.

results = get.variables(lpModel)
names(results) = dimnames(lpModel)[[2]]
results

For more information on how to interact with lpSolve model objects in R, please see the documentation for the lpSolve and lpSolveAPI packages.

FluvialLandscapeLab/gangsta documentation built on May 6, 2019, 5:05 p.m.