Package to Facilitate ODE-Based Modeling

Description

This package provides methods to

• import a conceptual ODE-based model stored in tabular form (i.e. as text files or spreadsheets).

• generate source code (either R or Fortran) to be passed to an ODE-solver.

• visualize and export basic information about a model, e.g. for documentation purposes.

Details

Consult the package vignette for details. The concept of writing an ODE system in tabular/matrix form is nicely introduced, e. g., in the book of Reichert, P., Borchardt, D., Henze, M., Rauch, W., Shanahan, P., Somlyody, L., and Vanrolleghem, P. A. (2001): River water quality model No. 1, IWA publishing, ISBN 9781900222822.

The current source code repository is https://github.com/dkneis/rodeo.

Class and class methods

See rodeo-class for the rodeo class and the corresponding class methods.

Non-class methods

Type help(package="rodeo") or see the links below to access the documentation of non-class methods contained in the package.

• buildFromWorkbook Builds and compiles a model fully specified in a workbook (supports .xlsx and .ods format).

• forcingFunctions Generation of forcing functions in Fortran.

• exportDF Export of data frames as TEX or HTML code.

• stoiCreate Creates a stoichiometry matrix from a set of chemical reactions.

• stoiCheck Validates a stoichiometry matrix by checking for conservation of mass.

Author(s)

david.kneis@tu-dresden.de

See Also

Useful links:

• https://github.com/dkneis/rodeo

Build a model from the contents of delimited textfiles

Description

The function builds a rodeo-based model by importing all declarations and equations from tables stored as delimited text.

Usage

buildFromText(
  declarations,
  equations,
  sep = "\t",
  dim = 1,
  set_defaults = TRUE,
  fortran = FALSE,
  sources = NULL,
  ...
)

Arguments

Value

An object of class rodeo.

Note

The delimited text files provided as input are parsed by read.table with header=TRUE and the delimiter specified by sep. The files must contain the following:

• 'declarations' Declares the identifiers of state variables, parameters, and functions used in the model equations. Mandatory columns are 'type', 'name', 'unit', 'description', and 'default'. Entries in the type column must be one of 'variable', 'parameter', or 'function'. If source code is generated for R (fortran=FALSE), any declared functions must be accessible in the environment where the model is run. If fortran=TRUE, the functions must be implemented in the file(s) listed in sources to be included in compilation. Entries in the 'name' column must be unique, valid identifier names (character followed by characters, digits, underscore). Entries in the 'default' column shall be numeric.

• 'equations' Specifies the model equations. Mandatory columns are 'name', 'unit', 'description', 'rate' plus one column for every state variable of the model. The 'rate' columns holds math expressions for the process rates and columns named after state variables contain the corresponding expressions representing stoichiometric factors. All columns are of type character.

The best way to understand the contents of the input files is to study the examples in the folder 'models' shipped with the package. Type system.file("models", package="rodeo") at the R prompt to see where this folder is installed on your system. A full example is given below.

Author(s)

David Kneis david.kneis@tu-dresden.de

See Also

buildFromWorkbook provides similar functionality

Examples

# Build and run a SEIR type epidemic model
decl <- system.file("models/SEIR_declarations.txt", package="rodeo")
eqns <- system.file("models/SEIR_equations.txt", package="rodeo")
m <- buildFromText(decl, eqns)
x <- m$dynamics(times=0:30, fortran=FALSE)
print(head(x))

Build a model from the contents of a workbook

Description

The function builds a rodeo-based model by importing all declarations and equations from a workbook established with common spreadsheet software.

Usage

buildFromWorkbook(
  workbook,
  dim = 1,
  set_defaults = TRUE,
  fortran = FALSE,
  sources = NULL,
  ...
)

Arguments

Value

An object of class rodeo.

Note

The file provided as workbook must contain two sheets:

• 'declarations' Declares the identifiers of state variables, parameters, and functions used in the model equations. Mandatory columns are 'type', 'name', 'unit', 'description', and 'default'. Entries in the type column must be one of 'variable', 'parameter', or 'function'. If source code is generated for R (fortran=FALSE), any declared functions must be accessible in the environment where the model is run. If fortran=TRUE, the functions must be implemented in the file(s) listed in sources to be included in compilation. Entries in the 'name' column must be unique, valid identifier names (character followed by characters, digits, underscore). Entries in the 'default' column shall be numeric.

• 'equations' Specifies the model equations. Mandatory columns are 'name', 'unit', 'description', 'rate' plus one column for every state variable of the model. The 'rate' columns holds math expressions for the process rates and columns named after state variables contain the corresponding expressions representing stoichiometric factors. All columns are of type character.

The best way to understand the contents of a suitable workbook is to study the examples in the folder 'models' shipped with the package. Type system.file("models", package="rodeo") at the R prompt to see where this folder is installed on your system. A full example is given below.

Author(s)

David Kneis david.kneis@tu-dresden.de

See Also

buildFromText provides similar functionality

Examples

# Build and run a SEIR type epidemic model
m <- buildFromWorkbook(
  system.file("models/SEIR.xlsx", package="rodeo")
)
x <- m$dynamics(times=0:30, fortran=FALSE)
print(head(x))

Generate Executable Code

Description

Creates and 'compiles' a function for use with numerical methods from package deSolve or rootSolve.

Arguments

Value

invisible(NULL)

Note

The expected language of the external code passed in sources depends on the value of fortran.

If fortran is FALSE, R code is generated and made executable by eval and parse. Auxiliary code passed via sources is made available via source. The created R function is stored in the object.

If fortran is TRUE, the external code passed in sources must implement a module with the fixed name 'functions'. This module must contain all user-defined functions referenced in process rates or stoichiometric factors. The file names passed to the sources argument must carry an extension which is recognized by the compiler as a source file. Something like ".f95" should work.

If fortran is TRUE, a shared library is created. The library is immediately loaded with dyn.load and it is automatically unloaded with dyn.unload when the object's finalize method is called.

The name of the library (base name without extension) as well as the name of the function to compute the derivatives are stored in the object. These names can be queried with the libName and libFunc methods, respectively. Unless a file path is specified via the lib argument, the library is created in the folder returned by tempdir under a unique random name.

Author(s)

david.kneis@tu-dresden.de

See Also

This method internally calls generate.

Examples

data(vars, pars, funs, pros, stoi)
model <- rodeo$new(vars, pars, funs, pros, stoi, dim=c(1))
# This would trigger compilation assuming that 'functionsCode.f95' contains
# a Fortran implementation of all functions; see vignette for full example
## Not run: 
model$compile(sources="functionsCode.f95")

## End(Not run)

Numerical Integration

Description

Compute a dynamic solution with the numerical algorithms from package deSolve.

Arguments

Value

The matrix returned by the integrator (see ode).

Note

This method can only be used after compile has been called.

The ... argument should not be used to assign values to any of y, parms, times, func. If fortran is TRUE it should also not assign values to dllname, nout, or outnames. All these arguments of ode get their appropriate values automatically.

Author(s)

david.kneis@tu-dresden.de

See Also

Use step for integration over a single time step with a built-in, Fortran-based solver.

Export a Data Frame as HTML/TEX Code

Description

Generates code to include tabular data in a tex document or web site.

Usage

exportDF(
  x,
  tex = FALSE,
  colnames = NULL,
  width = NULL,
  align = NULL,
  funHead = NULL,
  funCell = NULL,
  lines = TRUE,
  indent = 2
)

Arguments

Value

A character string (usually needs to be exported to a file).

Note

The functions funHead and funCell are useful to apply formatting or character replacement. For example, one could use

function(x) {paste0("\\bold{",toupper(x),"}")}

to generate bold, uppercase column names in a TEX table.

Author(s)

David Kneis david.kneis@tu-dresden.de

See Also

The xtable packages provides similar functionality with more sophisticated options. Consider the 'pandoc' software do convert documents from one markup language to another one. Finally, consider the latex package 'datatools' for direct inclusion of delimited text files (e.g. produced by write.table) in tex documents.

Examples

# Create example table
df <- data.frame(stringsAsFactors=FALSE, name= c("growth", "dead"),
  unit= c("1/d","1/d"), expression= c("r * N * (1 - N/K)"," d * N"))

# Export as TEX: header in bold, 1st colum in italics, last column as math
tex <- exportDF(df, tex=TRUE,
  colnames=c(expression="process rate expression"),
  width=c(expression=0.5),
  align=c(expression="p"),
  funHead=setNames(replicate(ncol(df),
    function(x){paste0("\\textbf{",x,"}")}),names(df)),
  funCell=c(name=function(x){paste0("\\textit{",x,"}")},
    expression=function(x){paste0("$",x,"$")})
)
cat(tex,"\n")

# Export as HTML: non-standard colors are used for all columns
tf <- tempfile(fileext=".html")
write(x= exportDF(df, tex=FALSE,
  funHead=setNames(replicate(ncol(df),
    function(x){paste0("<font color='red'>",x,"</font>")}),names(df)),
  funCell=setNames(replicate(ncol(df),
    function(x){paste0("<font color='blue'>",x,"</font>")}),names(df))
), file=tf)
## Not run: 
  browseURL(tf)
  file.remove(tf)

## End(Not run)

Clean-up a rodeo Object

Description

Clean-up method for objects of the rodeo-class.

Value

The method is called implicitly for its side effects when a rodeo object is destroyed.

Note

At present, the method just unloads the object-specific shared libraries created with the compile or initStepper methods.

Author(s)

david.kneis@tu-dresden.de

Generation of Forcing Functions in Fortran

Description

Generates Fortran code to return the current values of forcing functions based on interpolation in tabulated time series data.

Usage

forcingFunctions(x)

Arguments

Value

A character string holding generated Fortran code. Must be written to disk, e.g. using write, prior to compilation.

Note

The fields of the input data frame are interpreted as follows:

• name Name of the forcing function as declared in the table of functions.

• file Name of the text file containing the time series data either as an absolute or relative path. Time information is expected as numeric values in the first column (e.g. as number of seconds after some reference date). The period is used as the decimal character in floating point numbers, numeric values can also be given in scientific format (e.g. as 0.314e+1). Allowed column delimiters are blank, tab, or comma. A sequence of white spaces collapses to a single delimiter but this is not the case for commas. It is strictly recommended to use a consistent delimiter character within a particular file. Blank lines are allowed everywhere in the file, comment lines must start with a '#'. The first non-blank, non-comment line is interpreted as column headers and the name of the first column (holding time info) is essentially ignored).

• column Name of the column in file from which data are to be read.

• mode Integer code to control how the interpolation is performed. Use 0 for constant interpolation with full weight given to the value at the end of a time interval. Use 1 for constant interpolation with full weight given to the value at the begin of a time interval. Any other values (< 0 or > 1) result in linear interpolation with weights being set automatically.

• default Logical. If FALSE, the generated function has the interface 'f(time)'. If TRUE, the generated function has a two-argument interface 'f(time, z)'. If the actual argument 'z' is NaN, the function behaves just like the single-argument version, i.e. interpolation in tabulated data is performed. If 'z' is not NaN, the function returns the value of 'z'.

The generated code provides a single module named 'forcings' which defines as many forcing functions as there are rows in x. The module 'forcings' needs to be made available to the compiler (either at the command line or via inclusion in another file with Fortran's include mechanism). In addition, it must be referenced in the module 'functions' with an appropriate 'use' statement (see example below).

The generated function return scalar values of type double precision. If an error condition is encountered, the return value of a functions equals the largest possible double precision value (generated by Fortran's 'huge' function). In addition, errors trigger calls of the subroutines 'rexit' (at top level) or 'rwarn' (at lower levels). These two functions are available automatically if the Fortran code is compiled using 'R CMD SHLIB'. Otherwise, the two functions need to be defined (see examples below).

In the two-argument version, the second argument is tested against NaN using 'ISNAN'. This function is not part of the Fortran standard but it is supported by most compilers, including gfortan. The Fortran 2003 standard conformal function would be 'IS_IEEE_NAN' which is not yet supported by compiler versions normally installed with R (March 2016).

Author(s)

David Kneis david.kneis@tu-dresden.de

Examples

## Not run: 
  ! Example of a Fortran file to define functions
  include 'forcings.f95'   ! include generated forcings file in compilation
  module functions
  use forcings             ! make forcings available as functions
  implicit none
  contains
  ! ... any non-forcing functions go here ...
  end module

## End(Not run)

## Not run: 
  ! Definition of 'rexit' and 'rwarn' for testing of the generated code
  ! outside of R
  subroutine rexit (x)
    character(len=*), intent(in):: x
    write(*,*) "ERROR: ",trim(adjustl(x))
    stop 1
  end subroutine
  
  subroutine rwarn (x)
    character(len=*), intent(in):: x
    write(*,*) "WARNING: ",trim(adjustl(x))
  end subroutine

## End(Not run)

Clear forcing data of a Fortran-based model

Description

Clears the forcings data of a model written in Fortran. This may be useful when data should be updated through another call to to forcings_init or for the purpose of memory cleanup.

Value

Returns NULL invisibly.

Note

The function takes no arguments and should thus called as forcings_clear().

Author(s)

david.kneis@tu-dresden.de

See Also

forcings_init to initialize the forcing data

Initialize forcings of a Fortran-based model

Description

Initializes the forcings of a model written in Fortran based on the input supplied to forcingFunctions. If the data cannot be initialized for any reason, an error is thrown and a hopefully useful message is provided.

Value

Returns NULL invisibly.

Note

The function takes no arguments and should thus called as forcings_init().

Author(s)

david.kneis@tu-dresden.de

See Also

forcings_clear to clear the data set by forcings_init

Declaration of Functions

Description

Declaration of functions referenced at the ODE's right hand sides of the bacteria growth example model.

Format

A data frame with the following fields:

name :

Name of the function

unit :

Unit of the return value

description :

Short description (text)

Code Generator

Description

This is a low-level method to translate the ODE-model specification into a function that computes process rates and the state variables derivatives (either in R or Fortran). You probably want to use the high-level method compile instead, which uses generate internally.

Arguments

Value

The generated source code as a string. Must be written to disk, e.g. using write, prior to compilation.

Note

Details of this low-level method may change in the future.

Author(s)

david.kneis@tu-dresden.de

See Also

See other methods of the rodeo-class, especially compile which internally uses this method.

Examples

data(vars, pars, funs, pros, stoi)
model <- rodeo$new(vars, pars, funs, pros, stoi, dim=c(1))
fortranCode <- model$generate(lang="f95")
## Not run: 
write(fortranCode, file="")

## End(Not run)

Query Values of Parameters

Description

Query values of parameters of a rodeo-based model.

Arguments

Value

A numeric vector or array.

Author(s)

david.kneis@tu-dresden.de

See Also

The corresponding 'set' method is setPars and examples can be found there. Use getVars to query the values of variables rather than parameters.

Query Values of State Variables

Description

Query values of state variables of a rodeo-based model.

Arguments

Value

A numeric vector or array.

Author(s)

david.kneis@tu-dresden.de

See Also

The corresponding 'set' method is setVars and examples can be found there. Use getPars to query the values of parameters rather than variables.

Initialize Internal ODE Solver

Description

Initializes rodeo's built-in ODE solver. This method must be called prior to using step. Only works with Fortran, not plain R.

Arguments

Value

invisible(NULL)

Note

After this method was called, step can be used to perform the integration.

Author(s)

david.kneis@tu-dresden.de

See Also

To perform integration with the solvers from package deSolve use compile instead of this method.

Initialize a rodeo Object

Description

Initializes an object of the rodeo-class with data frames holding the specification of an ODE system.

Arguments

Value

The method is called implicitly for its side effects when a rodeo object is instantiated with new. It has no accessible return value.

Note

The mandatory fields of the input data frames should be of type character. Additional fields may be present in these data frames and the contents becomes part of the rodeo object. The 'expression' fields of pros and stoi (or the contents of the stoichiometry matrix) should be valid mathematical expressions in R and Fortran. These can involve the names of declared state variables, parameters, and functions as well as numeric constants or basic math operators. Branching or loop constructs are not allowed (but these can appear inside referenced functions). There are currently few reserved words that cannot be used as variable, parameter, function, or process names. The reserved words are 'time', 'left', and 'right'.

Initialization does not assign numeric values to state variables or parameters. Use the decicated methods setVars and setPars for that purpose.

Author(s)

david.kneis@tu-dresden.de

See Also

See the package vignette for examples.

Examples

data(vars, pars, funs, pros, stoi)
model <- rodeo$new(vars, pars, funs, pros, stoi, dim=c(1))
print(model)

Return name of library function

Description

Return the name of the library function for use with deSolve or rootSolve methods.

Value

The name of the function to compute derivatives which is contained in the library built with compile. This name has to be supplied as the func argument of the solver methods in deSolve or rootSolve.

Author(s)

david.kneis@tu-dresden.de

Return library name

Description

Return the pure name of the shared library for use with deSolve or rootSolve methods.

Value

The base name of the shared library file created with compile after stripping of the the platform specific extension. This name has to be supplied as the dllname argument of the solver methods in deSolve or rootSolve.

Author(s)

david.kneis@tu-dresden.de

Declaration of Parameters

Description

Declaration of parameters of the bacteria growth example model.

Format

A data frame with the following fields:

name :

Name of the parameter

unit :

Unit

description :

Short description (text)

Plot Qualitative Stoichiometry Matrix

Description

Visualizes the stoichiometry matrix using standard plot methods. The sign of stoichiometric factors is displayed as upward and downward pointing triangles. Also visualized are dependencies of process rates on variables.

Arguments

Note

The values of state variables and parameters must have been set using the setVars and setPars methods. If the stoichiometric factors are mathematical expressions involving function references, these functions must be defined in R (even if the numerical computations are based on generated Fortran code).

Author(s)

david.kneis@tu-dresden.de

See Also

See other methods of the rodeo-class or stoichiometry for computing the stoichiometric factors only. Alternative options for displaying stoichiometry information are described in the package vignette.

Examples

data(vars, pars, funs, pros, stoi)
model <- rodeo$new(vars, pars, funs, pros, stoi, dim=c(1))
model$setVars(c(bac=0.1, sub=0.5))
model$setPars(c(mu=0.8, half=0.1, yield= 0.1, vol=1000, flow=50, sub_in=1))
monod <- function(c,h) {c/(c+h)}
model$plotStoichiometry(box=c(1))

Declaration of Processes

Description

Definition of processes of the bacteria growth example model.

Format

A data frame with the following fields:

name :

Name of the process

unit :

Unit of the rate expression

description :

Short description (text)

expression :

Mathematical expression (as a string)

rodeo Class

Description

This documents the rodeo class to represent an ODE-based model. See the rodeo-package main page or type help(package="rodeo") for an introduction to the package of the same name.

Fields

prosTbl

A data frame with fields 'name', 'unit', 'description', and 'expression' defining the process rates.

stoiTbl

A data frame with fields 'variable', 'process', and 'expression' reprenting the stoichiometry matrix in data base format.

varsTbl

A data frame with fields 'name', 'unit', 'description' declaring the state variables of the model. The declared names become valid identifiers to be used in the expression fields of prosTbl or stoiTbl.

parsTbl

A data frame of the same structure as vars declaring the parameters of the model. The declared names become valid identifiers to be used in the expression fields of prosTbl or stoiTbl.

funsTbl

A data frame of the same structure as vars declaring any functions referenced in the expression fields of prosTbl or stoiTbl.

dim

Integer vector specifying the spatial dimensions.

vars

Numeric vector, holding values of state variables.

pars

Numeric vector, holding values of parameters.

Class methods

For most of the methods below, a separate help page is available describing its arguments and usage.

• initialize Initialization method for new objects.

• namesVars, namesPars, namesFuns, namesPros Functions returning the names of declared state variables, parameters, functions, and processes, respectively (character vectors). No arguments.

• lenVars, lenPars, lenFuns, lenPros Functions returning the number of declared state variables, parameters, functions and processes, respectively (integer). No arguments.

• getVarsTable, getParsTable, getFunsTable, getProsTable, getStoiTable Functions returning the data frames used to initialize the object. No arguments

• getDim Returns the spatial dimensions as an integer vector. No arguments.

• compile Compiles a Fortran library for use with numerical methods from packages deSolve or rootSolve.

• generate Translate the ODE-model specification into a function that computes process rates and the state variables' derivatives (either in R or Fortran). Consider to use the high-level method compile.

• setVars Assign values to state variables.

• setPars Assign values to parameters.

• getVars Returns the values of state variables.

• getPars Returns the values of parameters.

• stoichiometry Returns the stoichiometry matrix, either evaluated (numeric) or as text.

• plotStoichiometry Plots qualitative stoichiometry information.

See Also

See the rodeo-package main page or type help(package="rodeo") to find the documentation of any non-class methods contained in the rodeo package.

Examples

## Example using high-level functions: Epidemiological SEIR model

# Import model from workbook shipped with this package 
m <- buildFromWorkbook(system.file("models/SEIR.xlsx",
  package="rodeo"))

# Set parameters and initial state (defaults stored in workbook)
p <- m$getParsTable()
m$setPars(setNames(p$default, p$name))

v <- m$getVarsTable()
m$setVars(setNames(v$default, v$name))

# Run a dynamic simulations and print parts of the result
sim <- m$dynamics(time=1:30, fortran=FALSE)
print(head(sim))
print(tail(sim))


### Example using low-level functions: Bacterial growth in bioreactor

# Creation of model object (data frames shipped as package data)
data(vars, pars, funs, pros, stoi)
model <- rodeo$new(vars, pars, funs, pros, stoi, dim=c(1))

# Assignment of parameters and initial values
model$setPars(c(mu=0.8, half=0.1, yield= 0.1, vol=1000,
  flow=50, sub_in=1))
model$setVars(c(bac=0.01, sub=0))

# Implementation of functions declared in 'funs'
monod <- function(c,h) {c/(c+h)}

# Creation of derivatives function in a low-level way; calling
# the 'compile' method is a more convenient alternative
code <- model$generate(name="derivs", lang="r")
derivs <- eval(parse(text=code))

# Explicit call of an integrator from the deSolve package
times <- 0:96
out <- deSolve::ode(y=model$getVars(), times=times, func=derivs,
  parms=model$getPars())
colnames(out) <- c("time", model$namesVars(), model$namesPros())
plot(out)

Run a zero-dimensional model, possibly for multiple scenarios

Description

Triggers dynamic simulation(s) of a 0D model for one to many scenarios. The individual scenarios can differ by the initial state or the values of parameters. Optionally produces basic plots.

Arguments

Value

A data frame with at least three columns and one row for each time requested via the times argument. The first column indicates the scenario, the second column holds the time, and further columns hold the values of state variables and process rates.

Note

If the scenarios argument is used to update initial values and / or parameters, the following applies: For each parameter not being included (by name) in the vector for a particular scenario, the default value will be used (as stored in the workbook from which the model was imported). The same is true for the initial values of variables. See the examples below.

Author(s)

David Kneis david.kneis@tu-dresden.de

See Also

Look at buildFromWorkbook for how to create a suitable model object.

Examples

# build the model
m <- buildFromWorkbook(system.file("models/SEIR.xlsx", package="rodeo"))

# run with defaults
x <- m$scenarios(times=0:30, scenarios=NULL)

# run with updated values 
x <- m$scenarios(times=0:30,
  scenarios=list(default=c(t=1, i=0.2, r=0.4), fast=c(t=1, i=0.5, r=0.1))
)

Assign Values to Parameters

Description

Assign values to parameters of a rodeo-based model.

Arguments

Value

NULL (invisible). The assigned numeric data are stored in the object and can be accessed by the getPars method.

Note

Look at the notes and examples for the setVars method.

Author(s)

david.kneis@tu-dresden.de

See Also

The corresponding 'get' method is getPars. Use setVars to assign values to variables rather than parameters. Consult the help page of the latter function for examples.

Assign Values to State Variables

Description

Assign values to state variables of a rodeo-based model.

Arguments

Value

NULL (invisible). The assigned numeric data are stored in the object and can be accessed by the getVars method.

Note

For a 0-dimensional model (i.e. a model without spatial resolution), x must be a numeric vector whose length equals the number of state variables. The element names of x must match those returned by the object's namesVars method. See the examples for how to bring the vector elements into required order.

For models with a spatial resolution, x must be a numeric array of proper dimensions. The last dimension (cycling slowest) corresponds to the variables and the first dimension (cycling fastest) corresponds to the models' highest spatial dimension. Thus, dim(x) must be equal to c(rev(obj$getDim()), obj$namesVars()) where obj is the object whose variables are to be assigned. The names of the array's last dimension must match the return value of obj$namesVars().

In the common 1-dimensional case, this just means that x must be a matrix with column names matching the return value of obj$namesVars() and as many rows as given by obj$getDim().

Author(s)

david.kneis@tu-dresden.de

See Also

The corresponding 'get' method is getVars. Use setPars to assign values to parameters rather than variables.

Examples

data(vars, pars, funs, pros, stoi)
x0 <- c(bac=0.1, sub=0.5)

# 0-dimensional model
model <- rodeo$new(vars, pars, funs, pros, stoi, dim=c(1))
model$setVars(x0)
print(model$getVars())

# How to sort vector elements
x0 <- c(sub=0.5, bac=0.1)              # doesn't match order of variables
model$setVars(x0[model$namesVars()])

# 1-dimensional model with 3 boxes
nBox <- 3
model <- rodeo$new(vars, pars, funs, pros, stoi, dim=c(nBox))
x1 <- matrix(rep(x0, each=nBox), nrow=nBox, ncol=model$lenVars(),
  dimnames=list(NULL, model$namesVars()))
model$setVars(x1)
print(model$getVars())
print(model$getVars(asArray=TRUE))

# 2-dimensional model with 3 x 4 boxes
model <- rodeo$new(vars, pars, funs, pros, stoi, dim=c(3,4))
x2 <- array(rep(x0, each=3*4), dim=c(4,3,model$lenVars()),
  dimnames=list(dim2=NULL, dim1=NULL, variables=model$namesVars()))
model$setVars(x2)
print(model$getVars())
print(model$getVars(asArray=TRUE))

Numerical Integration Over a Single Time Step

Description

Performs integration over a single time step using a built-in ODE solver. At present, a single solver is implement with limited options. The interface of this method may change when support for other solvers is added.

Arguments

Value

A named numeric vector holding the values of state variables and process rates in all boxes.

Note

This method can only be used after a call to initStepper has been made.

Author(s)

david.kneis@tu-dresden.de

See Also

Use deSolve for advanced solvers with more options and capabilities to handle stiff problems.

Specification of Stoichiometry

Description

Definition of the links between simulated processes and state variables in the bacteria growth example model.

Format

A data frame with the following fields:

variable :

Name of the state variable

process :

Name of the process

expression :

Mathematical expression (as a string)

Validation of a Stoichiometry Matrix

Description

Validates the stoichiometry matrix by checking for conservation of mass (more typically conservation of moles).

Usage

stoiCheck(stoi, comp, env = globalenv(), zero = .Machine$double.eps * 2)

Arguments

Value

A numeric matrix with the following properties:

• There is one row for each process, thus the number and names of rows are the same as in stoi.

• There is one column per checked element, hence column names are equal to the row names of comp.

• A matrix element at position [i,k] represent the mass balance for the process in row i with respect to the element in column k. A value of zero indicates a close balance. Positive (negative) values indicate that mass is gained (lost) in the respective process.

Author(s)

David Kneis david.kneis@tu-dresden.de

See Also

Use stoiCreate to create a stoichiometry matrix from a set of reactions in common notation.

Examples

# Eq. 1 and 2 are from Soetaert et al. (1996), Geochimica et Cosmochimica
# Acta, 60 (6), 1019-1040. 'OM' is organic matter. Constants 'nc' and 'pc'
# represent the nitrogen/carbon and phosphorus/carbon ratio, respectively.
reactions <- c(
  oxicDegrad= "OM + O2 -> CO2 + nc * NH3 + pc * H3PO4 + H2O",
  denitrific= "OM + 0.8*HNO3 -> CO2 + nc*NH3 + 0.4*N2 + pc*H3PO4 + 1.4*H2O",
  dissPhosp1= "H3PO4 <-> H + H2PO4",
  dissPhosp2= "H2PO4 <-> H + HPO4"
)
# Non-evaluated stoichiometry matrix
stoi <- stoiCreate(reactions, toRight="_f", toLeft="_b")
# Parameters ('nc' and 'pc' according to Redfield ratio)
pars <- list(nc=16/106, pc=1/106)
# Elemental composition
comp <- rbind(
  OM=    c(C=1, N="nc", P="pc", H="2 + 3*nc + 3*pc"),
  O2=    c(C=0, N=0,    P=0,    H=0),
  CO2=   c(C=1, N=0,    P=0,    H=0),
  NH3=   c(C=0, N=1,    P=0,    H=3),
  H3PO4= c(C=0, N=0,    P=1,    H=3),
  H2PO4= c(C=0, N=0,    P=1,    H=2),
  HPO4=  c(C=0, N=0,    P=1,    H=1),
  H2O=   c(C=0, N=0,    P=0,    H=2),
  HNO3=  c(C=0, N=1,    P=0,    H=1),
  N2=    c(C=0, N=2,    P=0,    H=0),
  H=     c(C=0, N=0,    P=0,    H=1)
)
# We need the transposed form
comp <- t(comp)
# Mass balance check
bal <- stoiCheck(stoi, comp=comp, env=pars)
print(bal)
print(colSums(bal) == 0)

Stoichiometry Matrix from Reaction Equations

Description

Creates a stoichiometry matrix from a set of reaction equations.

Usage

stoiCreate(
  reactions,
  eval = FALSE,
  env = globalenv(),
  toRight = "_forward",
  toLeft = "_backward"
)

Arguments

Value

A matrix with the following properties:

• The number of columns equals the total number of components present in reactions. The components' names are used as column names.

• The number of rows equals the length of reactions plus the number of reversible reactions. Thus, a single row is created for each non-reversible reaction but two rows are created for reversible ones. The latter represent the forward and backward reaction (in that order). The row names are constructed from the names of reactions, making use of the suffixes toRight and toLeft in the case of reversible reactions.

• The matrix is filled with the stoichiometric factors extracted from reactions. Empty elements are set to zero.

• The type of the matrix (character or numeric) depends on the value of eval.

Note

The syntax rules for reaction equations are as follows (see examples):

• There must be a left hand side and a right hand side. Sides must be separated by one of the arrows '->', '<-', or '<->' with the latter indicating a reversible reaction.

• Names of component(s) must appear at each side of the reaction. These must be legal row/column names in R. If multiple components are consumed or produced, they must be separated by '+'.

• Any stoichiometric factors need to appear before the respective component name using '*' as the separating character. Stoichiometric factors being equal to unity can be omitted.

• A stoichiometric factor is treated as a mathematical expression. In common cases, it is just a numeric constant. However, the expression can also involve references to variables or functions. If such an expression contains math operators '*' or '+' it needs to be enclosed in parenthesis.

Author(s)

David Kneis david.kneis@tu-dresden.de

See Also

Use stoiCheck to validate the mass balance of the generated matrix.

Examples

# EXAMPLE 1: From https://en.wikipedia.org/wiki/Petersen_matrix (July 2016)
#
reactions <- c(
  formS=  "A + 2 * B -> S",
  equiES= "E + S <-> ES",
  decoES= "ES -> E + P"
)
stoi <- stoiCreate(reactions, eval=TRUE, toRight="_f", toLeft="_b")
print(stoi)

# EXAMPLE 2: Decomposition of organic matter (selected equations only)
#
# Eq. 1 and 2 are from Soetaert et al. (1996), Geochimica et Cosmochimica
# Acta, 60 (6), 1019-1040. 'OM' is organic matter. Constants 'nc' and 'pc'
# represent the nitrogen/carbon and phosphorus/carbon ratio, respectively.
reactions <- c(
  oxicDegrad= "OM + O2 -> CO2 + nc * NH3 + pc * H3PO4 + H2O",
  denitrific= "OM + 0.8*HNO3 -> CO2 + nc*NH3 + 0.4*N2 + pc*H3PO4 + 1.4*H2O",
  dissPhosp1= "H3PO4 <-> H + H2PO4",
  dissPhosp2= "H2PO4 <-> H + HPO4"
)
# Non-evaluated matrix
stoi <- stoiCreate(reactions, toRight="_f", toLeft="_b")
print(stoi)
# Evaluated matrix ('nc' and 'pc' according to Redfield ratio)
pars <- list(nc=16/106, pc=1/106)
stoi <- stoiCreate(reactions, eval=TRUE, env=pars, toRight="_f", toLeft="_b")
print(stoi)

Return the Stoichiometry Matrix

Description

Return and optionally evaluate the mathematical expression appearing in the stoichiometry matrix.

Arguments

Value

A matrix of numeric or character type, depending on the value of box.

Note

If the stoichiometric factors are mathematical expressions involving function references, these functions must be defined in R (even if the numerical computations are based on generated Fortran code).

Author(s)

david.kneis@tu-dresden.de

See Also

See other methods of the rodeo-class or plotStoichiometry for a graphical representation of the stoichiometric factors only.

Examples

data(vars, pars, funs, pros, stoi)
model <- rodeo$new(vars, pars, funs, pros, stoi, dim=c(1))
model$setPars(c(mu=0.8, half=0.1, yield= 0.1, vol=1000, flow=50, sub_in=1))
model$setVars(c(bac=0.1, sub=0.5))
monod <- function(c,h) {c/(c+h)}
print(model$stoichiometry(box=NULL))
print(model$stoichiometry(box=c(1)))

Declaration of Variables

Description

Declaration of variables of the bacteria growth example model.

Format

A data frame with the following fields:

name :

Name of the variable

unit :

Unit

description :

Short description (text)