ODE-based modeling

Basics and the R package “rodeo”

david.kneis!@!tu-dresden.de

Models …

Fundamental types of models

Fit for extrapolation?

A tree of techniques

Basics of ODE

Quantifying rates of change

  • Differential equations have derivatives on their left-hand side

    \[\frac{dy}{dx} = \ldots\]

  • Expresses the rate of change of a state variable (\(y\)) with respect to another quantity (\(x\))

Derivatives w.r.t. time

Derivatives in the spatial domain

Ordinary vs. partial DE


Derivatives in a single dimension \(\rightarrow\) ODE

\[\frac{dy}{dt} = \ldots\]


Derivatives in multiple dimensions \(\rightarrow\) PDE

\[\frac{\partial y}{\partial t} = \ldots \frac{\partial y}{\partial x} \ldots \]

Elements of a typical ODE

\[\frac{dy}{dt} = f(y, p, t)\]

y State variables (e.g. temperature, O2)

Elements of a typical ODE

\[\frac{dy}{dt} = f(y, p, t)\]

y State variables (e.g. temperature, O2)
p Parameters (e.g. radiation, reflection, heat capacity, …)

Elements of a typical ODE

\[\frac{dy}{dt} = f(y, p, t)\]

y State variables (e.g. temperature, O2)
p Parameters (e.g. radiation, reflection, heat capacity, …)
t Time; appears on right hand side if external forcing varies

Use case: Forward simulation


Prediction and simulation of hypothetical scenarios

Use case: Inverse modeling


Estimation of parameters from observations

Use case: System identification


Analysis of discrepancies between modeled and observed patterns

Simultaneous ODE

Simultaneous ODE


  • Aerobic bacteria (B) grow on organic substrate (S)

  • Well mixed culture

  • Oxygen (X) exchanged with atmosphere

Simultaneous ODE

Exponential growth


\[ \begin{align} \frac{d}{dt} B & = \mu \cdot B \end{align} \]

B Bacterial density (mmol C / L)
µ Growth rate constant (1 / h)

Simultaneous ODE

Accounting for substrate limitation


\[ \begin{align} \frac{d}{dt} B & = \mu \cdot B \cdot \left( \frac{S}{S+h_S} \right) \end{align} \]

S Substrate (mmol C / L)
hS Monod parameter (mmol C / L)

Simultaneous ODE

Accounting for oxygen limitation


\[ \begin{align} \frac{d}{dt} B & = \mu \cdot B \cdot \left( \frac{S}{S+h_S} \right) \cdot \left( \frac{X}{X+h_X} \right) \end{align} \]

X Oxygen (mmol O2 / L)
hX Monod parameter (mmol O2 / L)

Simultaneous ODE

Equivalent for the second state variable


\[ \begin{align} \frac{d}{dt} B & = & \mu \cdot B \cdot \left( \tfrac{S}{S+h_S} \right) \cdot \left( \tfrac{X}{X+h_X} \right) \\ \frac{d}{dt} S & = & \color{orange}{-\frac{1}{f}} \cdot \mu \cdot B \cdot \left( \tfrac{S}{S+h_S} \right) \cdot \left( \tfrac{X}{X+h_X} \right) \end{align} \]

f Yield parameter (mmol B / mmol S)

Simultaneous ODE

\[ \begin{align} \frac{d}{dt} B & = & \mu \cdot B \cdot \left( \tfrac{S}{S+h_S} \right) \cdot \left( \tfrac{X}{X+h_X} \right) \\ \frac{d}{dt} S & = & \color{orange}{-\tfrac{1}{f}} \cdot \mu \cdot B \cdot \left( \tfrac{S}{S+h_S} \right) \cdot \left( \tfrac{X}{X+h_X} \right) \\ \frac{d}{dt} X & = & \color{teal}{\left(1-\tfrac{1}{f} \right)} \cdot \mu \cdot B \cdot \left( \tfrac{S}{S+h_S} \right) \cdot \left( \tfrac{X}{X+h_X} \right) \\ & & + \color{magenta}{k \cdot \left(X_{sat} - X \right)} \end{align} \]

k Exchange rate constant (1 / hour)
Xsat Oxygen saturation level (mmol / L)

Simultaneous ODE

\[ \begin{align} \frac{d}{dt} B & = & \mu \cdot B \cdot \left( \tfrac{S}{S+h_S} \right) \cdot \left( \tfrac{X}{X+h_X} \right) \\ \frac{d}{dt} S & = & \color{orange}{-\tfrac{1}{f}} \cdot \mu \cdot B \cdot \left( \tfrac{S}{S+h_S} \right) \cdot \left( \tfrac{X}{X+h_X} \right) \\ \frac{d}{dt} X & = & \color{teal}{\left(1-\tfrac{1}{f} \right)} \cdot \mu \cdot B \cdot \left( \tfrac{S}{S+h_S} \right) \cdot \left( \tfrac{X}{X+h_X} \right) \\ & & + \color{magenta}{k \cdot \left(X_{sat} - X \right)} \end{align} \]

  • Terms heavily redundant

  • All variables on all RHS \(\rightarrow\) simultaneous solution

Simultaneous ODE in general

\[ \begin{align} \frac{d}{dt} \color{red}{y_1} & = f_1(\left[\color{red}{y_1}, \color{blue}{y_2}, \ldots, \color{green}{y_n}\right], p, t) \\ \frac{d}{dt} \color{blue}{y_2} & = f_2(\left[\color{red}{y_1}, \color{blue}{y_2}, \ldots, \color{green}{y_n}\right], p, t) \\ \ldots & \\ \frac{d}{dt} \color{green}{y_n} & = f_3(\left[\color{red}{y_1}, \color{blue}{y_2}, \ldots, \color{green}{y_n}\right], p, t) \end{align} \]

Solutions in forward modeling

1st case: Simulation of dynamics

  • From a known initial state, we predict future states

  • “Initial Value Problem” (IVP)

  • Requires integration over time (e.g., R package “deSolve”)

2nd case: Steady state

  • Given constant forcings, systems may run into steady state

  • All fluxes are in perfect balance

Option 1: Set LHS to zero and solve for state variables (R package “rootSolve”): Elegant but “fragile”

Option 2: Long-term integration: Fail-proof but slow

Integration algorithms

  • Analytical solutions exist for simple cases only

  • Useful for testing accuracy of numerical algorithms

Integration algorithms

Integration algorithms

Integration algorithms

Widely used

  • Methods of Runge & Kutta

  • Fails on stiff systems

Recommended

  • LSODA from ODEPACK

  • Efficient on both well-behaved and stiff systems

  • Default of deSolve

Concepts of the “rodeo” R package

Simplified lake P budget

Even simpler

\(P\) Total P in water column (g m-3)
\(P_{sed}\) Total sediment P (g m-2)

Set of ODE

\[ \begin{align} \frac{d}{dt} P & = & ... \\ \frac{d}{dt} P_{sed} & = & ... \end{align} \]

  • Expressions for processes that link water and sediment would appear redundantly in both ODE

  • Chemical engineers invented a clever notation …

Matrix-based notation for ODE

\[ \begin{align} \begin{bmatrix} \tfrac{d}{dt}y_1 \\ \tfrac{d}{dt}y_2 \\ \ldots \\ \tfrac{d}{dt}y_n \end{bmatrix} ^T & = & \begin{bmatrix} a \\ b \\ \ldots \\ z \end{bmatrix} ^T & \cdot & \begin{bmatrix} s_{a,1}, & \ldots, & s_{a,n} \\ s_{b,1}, & \ldots, & s_{b,n} \\ \ldots, & \ldots, & \ldots \\ s_{z,1}, & \ldots, & s_{z,n} \\ \end{bmatrix} \end{align} \]

Vector of derivatives

Vector of rate expressions (processes)

Stoichiometry matrix

Matrix-based notation for ODE

\[ \begin{align} \begin{bmatrix} \color{blue}{\tfrac{d}{dt}y_1} \\ \tfrac{d}{dt}y_2 \\ \ldots \\ \tfrac{d}{dt}y_n \end{bmatrix} ^T & = & \begin{bmatrix} \color{teal}{a} \\ \color{orange}{b} \\ \ldots \\ \color{red}{z} \end{bmatrix} ^T & \cdot & \begin{bmatrix} \color{teal}{s_{a,1}}, & \ldots, & s_{a,n} \\ \color{orange}{s_{b,1}}, & \ldots, & s_{b,n} \\ \ldots, & \ldots, & \ldots \\ \color{red}{s_{z,1}}, & \ldots, & s_{z,n} \\ \end{bmatrix} \end{align} \]


\[ \color{blue}{\tfrac{d}{dt}y_1} = \color{teal}{a} \cdot \color{teal}{s_{a,1}} + \color{orange}{b} \cdot \color{orange}{s_{b,1}} + \ldots + \color{red}{z} \cdot \color{red}{s_{z,1}} \]

Application to P budget

\[ \begin{align} \begin{bmatrix} \tfrac{d}{dt}\color{blue}{P} \\ \tfrac{d}{dt}\color{brown}{P_{sed}} \end{bmatrix} ^T = \begin{bmatrix} 1/\tau \cdot P_{in} \\ 1/\tau \cdot P \\ u/d \cdot P \\ k_B \cdot P_{sed} \\ k_R \cdot P_{sed} \end{bmatrix} ^T \cdot \begin{bmatrix} \color{blue}{1} & \color{brown}{0} \\ \color{blue}{-1} & \color{brown}{0} \\ \color{blue}{-1} & \color{brown}{d} \\ \color{blue}{0} & \color{brown}{-1} \\ \color{blue}{1/d} & \color{brown}{-1} \end{bmatrix} \end{align} \]

\(\tau\) Residence time (y)
\(d\) Water depth (m)
\(u\) Settl. velocity (m y-1)
\(P_{in}\) P in inflow (g / m3)
\(k_B\) Burial rate const. (y-1)
\(k_R\) Remob. rate const. (y-1)

Application to P budget

\[ \begin{align} \begin{bmatrix} \color{red}{\tfrac{d}{dt}P} \\ \tfrac{d}{dt}P_{sed} \end{bmatrix} ^T = \begin{bmatrix} \color{blue}{1/\tau \cdot P_{in}} \\ \color{magenta}{1/\tau \cdot P} \\ \color{olive}{u/d \cdot P} \\ k_B \cdot P_{sed} \\ \color{teal}{k_R \cdot P_{sed}} \end{bmatrix} ^T \cdot \begin{bmatrix} \color{blue}{1} & 0 \\ \color{magenta}{-1} & 0 \\ \color{olive}{-1} & d \\ 0 & -1 \\ \color{teal}{1/d} & -1 \end{bmatrix} \end{align} \]

\[\color{red}{\tfrac{d}{dt}P} = \color{blue} {\tfrac{1}{\tau} \cdot P_{in}} \color{magenta}{- \tfrac{1}{\tau} \cdot P} \color{olive}{- \tfrac{u}{d} \cdot P} \color{teal}{+ k_R \cdot P_{sed} \cdot \tfrac{1}{d}}\]

Import, Export, Settling, Burial, Remobilization

Advantages of the notation

\[ \begin{align} \begin{bmatrix} \tfrac{d}{dt}P \\ \tfrac{d}{dt}P_{sed} \end{bmatrix} ^T = \begin{bmatrix} 1/\tau \cdot P_{in} \\ 1/\tau \cdot P \\ u/d \cdot P \\ k_B \cdot P_{sed} \\ k_R \cdot P_{sed} \end{bmatrix} ^T \cdot \begin{bmatrix} 1 & 0 \\ -1 & 0 \\ -1 & d \\ 0 & -1 \\ 1/d & -1 \end{bmatrix} \end{align} \]


Process rate expression do not appear redundantly

  • clarity, consistency
  • faster computations

Advantages of the notation

\[ \begin{align} \begin{bmatrix} \tfrac{d}{dt}P \\ \tfrac{d}{dt}P_{sed} \end{bmatrix} ^T = \begin{bmatrix} 1/\tau \cdot P_{in} \\ 1/\tau \cdot P \\ u/d \cdot P \\ k_B \cdot P_{sed} \\ k_R \cdot P_{sed} \end{bmatrix} ^T \cdot \begin{bmatrix} 1 & 0 \\ -1 & 0 \\ -1 & d \\ 0 & -1 \\ 1/d & -1 \end{bmatrix} \end{align} \]


Straightforward modification

  • to add/drop a process, we add/drop a row

  • similar for state variables

Advantages of the notation

\[ \begin{align} \begin{bmatrix} \tfrac{d}{dt}P \\ \tfrac{d}{dt}P_{sed} \end{bmatrix} ^T = \begin{bmatrix} 1/\tau \cdot P_{in} \\ 1/\tau \cdot P \\ u/d \cdot P \\ k_B \cdot P_{sed} \\ k_R \cdot P_{sed} \end{bmatrix} ^T \cdot \begin{bmatrix} 1 & 0 \\ -1 & 0 \\ -1 & d \\ 0 & -1 \\ 1/d & -1 \end{bmatrix} \end{align} \]


  • The stoichiometry matrix is extremely useful for communicating the model

  • Signs of coefficients can be translated to symbols

Advantages of the notation

\[ \begin{align} \begin{matrix} & P & P_{sed} \\ \text{Import} & \color{red}{\uparrow} & \circ \\ \text{Export} & \color{blue}{\downarrow} & \circ \\ \text{Settling} & \color{blue}{\downarrow} & \color{red}{\uparrow} \\ \text{Burial} & \circ & \color{blue}{\downarrow} \\ \text{Remob.} & \color{red}{\uparrow} & \color{blue}{\downarrow} \end{matrix} \end{align} \]


\(\color{red}{\uparrow}\) increase, \(\color{blue}{\downarrow}\) decline, \(\circ\) no direct effect

More advantages

Vectors and matrices can be stored in tables, alongside with documentation


More advantages

Tables can be processed automatically to …

  • generate source code for simulation

  • perform sanity checks

  • create nicely formatted documentation

R package “rodeo”

Available on CRAN since around 2017

  • Employs the matrix-based notation to separate equations from source code

  • Provides a code generator compliant with integration routines from deSolve

    • native R code

    • compiled Fortran 2000+

  • Supports solving PDE with the “method-of-lines”

Implementation of the P budget model

\[ \begin{align} \begin{bmatrix} \tfrac{d}{dt}P \\ \tfrac{d}{dt}P_{sed} \end{bmatrix} ^T = \begin{bmatrix} 1/\tau \cdot P_{in} \\ 1/\tau \cdot P \\ u/d \cdot P \\ k_B \cdot P_{sed} \\ k_R \cdot P_{sed} \end{bmatrix} ^T \cdot \begin{bmatrix} 1 & 0 \\ -1 & 0 \\ -1 & d \\ 0 & -1 \\ 1/d & -1 \end{bmatrix} \end{align} \]

Tabular representation

\[ \begin{align} \begin{bmatrix} \tfrac{d}{dt}\color{red}{P} \\ \tfrac{d}{dt}\color{red}{P_{sed}} \end{bmatrix} ^T = \begin{bmatrix} 1/\color{teal}{\tau} \cdot \color{teal}{P_{in}} \\ 1/\color{teal}{\tau} \cdot \color{red}{P} \\ \color{teal}{u}/\color{teal}{d} \cdot \color{red}{P} \\ \color{teal}{k_B} \cdot \color{red}{P_{sed}} \\ \color{teal}{k_R} \cdot \color{red}{P_{sed}} \end{bmatrix} ^T \cdot \begin{bmatrix} 1 & 0 \\ -1 & 0 \\ -1 & \color{teal}{d} \\ 0 & -1 \\ 1/\color{teal}{d} & -1 \end{bmatrix} \end{align} \]


1st table Math expressions of RHS
2nd table Declaration & doc. of state variables
3rd table Declaration & doc. of parameters

1st table: Variables

name unit description initial
P g m-3 phosphorus in water column 0.1
Psed g m-3 sediment phosphorus 1.0

2nd table: Parameters

name unit description default
tau years residence time 1.0
d m lake depth 10.0
u m y-1 apparent settling velocity 40.0
kB y-1 rate constant of burial 0.1
kR y-1 remobilization rate const. 0.3
Pin g m-3 inflow concentration 0.1

3rd table: Rates & Stoichiometry

name unit description expression P Psed
imp g m-3 y-1 import 1 / tau * Pin 1 0
exp g m-3 y-1 export 1 / tau * P -1 0
set g m-3 y-1 settling u / d * P -1 d
bur g m-2 y-1 burial kB * Psed 0 -1
rem g m-2 y-1 remobil. kR * Psed 1/d -1

Import sheets from workbook

library("readxl")

read.sheet <- function(sheet) {
  x <- read_excel("models/phosphorus.xlsx", sheet)
  data.frame(x, row.names=x[["name"]])
}

vars <- read.sheet("vars")    # variables
pars <- read.sheet("pars")    # parameters
eqns <- read.sheet("eqns")    # rates & stoich.

The model is an “R6” object

library("rodeo")

PB <- rodeo$new(
  vars= vars,
  pars= pars,
  funs= NULL,
  pros= eqns,
  stoi= as.matrix(eqns[,vars[,"name"]]),
  asMatrix= TRUE,
  dim= 1
)

Code generation

PB$compile(fortran=FALSE)


  • For cheap models, R code is sufficient

  • Compiled Fortran code is way faster

  • Needs “gfortran” (on Windows: “Rtools”)

Supplying initial values

v <- setNames(vars[,"initial"], vars[,"name"])

print(v)
   P Psed 
 0.1  1.0 


PB$setVars(v)

print(PB$getVars())
   P Psed 
 0.1  1.0

Supplying parameter values

p <- setNames(pars[,"default"], pars[,"name"])

PB$setPars(p)

print(PB$getPars())
 tau    d    u   kB   kR  Pin 
 1.0 10.0 40.0  0.1  0.3  0.1

Integration

x <- PB$dynamics(times=0:50, fortran=FALSE)
print(dim(x))
[1] 51  8


print(x[1:3 ,1:5])
     time          P     Psed imp        exp
[1,]    0 0.10000000 1.000000 0.1 0.10000000
[2,]    1 0.03223815 2.066405 0.1 0.03223815
[3,]    2 0.03443036 2.483850 0.1 0.03443036

Visualization

par(mfcol=c(1, PB$lenVars()))
for (v in PB$namesVars())
  plot(x[,c("time",v)], ty="b", xlab="Yr", ylab=v)

Plotting the model structure

PB$plotStoichiometry(box=1, cex=1.5)

Summary

Summary: ODE

  • Nature of simultaneous ODE

  • Numerical solutions: Dynamics or steady-state ?

  • Numerous advantages of a matrix-based notation (Petersen matrix)

Summary: R package rodeo

\[ \begin{align} \begin{bmatrix} \tfrac{d}{dt}P \\ \tfrac{d}{dt}P_{sed} \end{bmatrix} ^T = \begin{bmatrix} 1/\tau \cdot P_{in} \\ 1/\tau \cdot P \\ u/d \cdot P \\ k_B \cdot P_{sed} \\ k_R \cdot P_{sed} \end{bmatrix} ^T \cdot \begin{bmatrix} 1 & 0 \\ -1 & 0 \\ -1 & d \\ 0 & -1 \\ 1/d & -1 \end{bmatrix} \end{align} \]

Adopts the matrix-based notation

ODEs held in tables with basic documentation

deSolve::lsoda( … )

Generated source code passed to numeric solver

Workshop

Contents

  • Examples provided (workbook + R code)

  • Walk through examples

  • Learn by simple modifications

  • Introduction to solving PDE problems (MOL)

Workshop prerequisites

  • Web browser

  • Spreadsheet software, e.g., Libre Office “calc”

  • A version of R with these CRAN packages installed: deSolve, rodeo, readxl

  • The tools to compile source code in R packages (gfortran, “Rtools” on Windows)