Numerical Simulation of Dynamic Systems with R: An Introduction

background-size: cover

# Numerical Simulation of Dynamic Systems With R: An Introduction

#### Thomas Petzoldt and Karline Soetaert

Contribution to the "Open Data and Open Science in the Aquatic Sciences" workshop<br>
https://aquaticdatasciopensci.github.io/
<br>
July 16th 2021.

.small[Version 2021-07-13, source code freely available from [https://github.com/dynamic-R](https://github.com/dynamic-R/hacking-limnology)]

---

Use cursor keys for navigation, press .red["O"] for a slide .red[O]verview

---

## Summary

Dynamical systems are found everywhere, in mathematics, physics, chemistry, engineering and business - ecology and aquatic sciences are no exception.

One common approach to describe such systems is by means of differential equations. Following the ideas of Forrester, differential equations appear quite naturally if we describe changes in a system in terms of growth and decay, which together make up a mass balance. A somewhat bigger challenge is to solve these differential equations - some of us may still remember the challenges of differential calculus in high school courses.

Fortunately, computer algorithms allow to solve even complex differential equation systems numerically. This has opened a world of practical applications to be accessible for anyone that has a basic knowledge in scientific computing, or is willing to acquire this knowledge. R is such a scientific computing language that offers powerful methods to solve differential equations. Moreover, some R-packages are especially designed to also solve spatially variable problems that are often found of importance in aquatic sciences.

The introduction will demonstrate selected examples: growth of organisms, predator-prey interaction, spread of diseases, and transport-reaction problems. The formulation of such models in R can be surprisingly compact and close to the mathematical equations.

---
class: center, middle, title-slide
background-image: url("img/water-light.jpg")

# Introduction

---
J.W. Forrester (1918-2016)

## Traditional linear thinking

<br>

.center[
<div id="htmlwidget-391f312db952e03994f6" style="width:1008px;height:108px;" class="grViz html-widget"></div>
<script type="application/json" data-for="htmlwidget-391f312db952e03994f6">{"x":{"diagram":"digraph causal {\n         graph [rankdir = \"LR\", bgcolor=\"none\"]\n           node [shape = box, penwidth=2, fontname = \"Helvetica\"]\n             Information Action Result\n           edge [penwidth=1.5]\n             Information -> Action -> Result\n}","config":{"engine":"dot","options":null}},"evals":[],"jsHooks":[]}</script>
]

---
J.W. Forrester (1918-2016)

## Dynamic thinking: Feedback loops

<br>

.center[
<div id="htmlwidget-65fde5b6bf6b9feea383" style="width:1008px;height:288px;" class="grViz html-widget"></div>
<script type="application/json" data-for="htmlwidget-65fde5b6bf6b9feea383">{"x":{"diagram":"digraph feedback1 {\n         graph [rankdir = \"LR\", bgcolor=\"none\"]\n           node [shape = box, penwidth=2, fontname = \"Helvetica\", minlen=3]\n             Information Action Result\n           edge [penwidth=1.5]\n             Information -> Action -> Result\n           edge [penwidth=0.7, tailport = \"s\", headport = \"s\", constraint = false, color=tomato, label=\"feedback\"]\n             Result -> Information\n           \n}","config":{"engine":"dot","options":null}},"evals":[],"jsHooks":[]}</script>
]

---

## Jay W. Forrester: System Dynamics

* flow of information, flow of matter
* thinking in feedback loops
* easy to understand, if we think in pools and changes

From: [Forrester, J.W. (2009) Some Basic Concepts in System Dynamics](https://web.mit.edu/sysdyn/sd-intro/D-4165-1.pdf)

---

## Dynamic Systems are Everywhere

.center[
![:scale 40%](img/batch.jpg)
![:scale 40%](img/small-stream.jpg)
![:scale 40%](img/marmots.jpg)
![:scale 40%](img/lake-dynamic.jpg)
]

---

## Example 1: Predator-Prey model

![Lotka-Volterra time series](img/predprey_ts.gif)![Lotka-Volterra state diagram](img/predprey_state.gif)

---

## Example 2: Nitrification in a River (1D)
.center[
![:scale 90%](img/river.gif)
]

---

## Example 3: Covid-19: Flatten the Curve!
.center[
<iframe src="https://weblab.hydro.tu-dresden.de/apps/sir-app/" width=800px, height=500px></iframe>
]
---

## Example 4: An Epidemiological Model in 2D

---
class: center, middle, title-slide
background-image: url("img/water-light.jpg")

# .darkblue[Description of Dynamic Systems]

---

## Common Properties of Dynamic Systems

### State at a given time

- vector of real numbers
- point in state space
]

### Evolution rule

- function that describes <br> .red[what future states follow from the current state]
]

<br> <br>
.gray[*) Wikipedia, [Dynamical_system](https://en.wikipedia.org/wiki/Dynamical_system)]

---
class: page-font-large

## Different Ways of Description

* static, mass-balanced snapshot of the system, or
* time series of state variables (pools)
]
#### .blue[Dynamic approach]

* description of .green[changes]
* .green[deterministic] .gray[or stochastic]
* .gray[individual-based models]
* .gray[discrete time: difference equations]
* equations in .green[continuous time] `$\rightarrow$` .green[differential equations] (ODEs)

---

## Differential Calculus

![:scale 35.1%](img/newton.jpg)

Isaac Newton, 1643 - 1727
]]

![:scale 38%](img/leibniz.jpg)

Gottfried Wilhelm Leibniz, 1646 - 1716
]]

--- 
.center[
##### `$\Rightarrow$` describes .red[changes] of pools
]

--- 
.small[.gray[
Newton: https://commons.wikimedia.org/wiki/File:IsaacNewton_1642-1727.jpg (CC-BY-SA 4.0) <br>
Leibniz: https://commons.wikimedia.org/wiki/File:Gottfried_Wilhelm_Leibniz,_Bernhard_Christoph_Francke.jpg (public domain)
]]

---
class: page-font-huge

## Dynamic systems: how and why?

"Many believe that system dynamics has helped them become skilled at
inventing the future, either by sketching out causal loops on the 
.green[back of an envelope], or by assembling equations of cause and effect in a
.green[computer model]. Both approaches work."

## Dynamic systems: how and why?

- .gray[discrete time: <br> another interesting story ...]
]

.pull-right[
#### Why numerical integration?
* Analytical solution<br> does not exist for all systems.

* Flexibility to simulate <br>time-dependent inputs.

.small[.green["... students can deal with high-order dynamic systems
without ever discovering that their elders consider such to be very difficult."]
.gray[J. W. Forrester 2009]]
]

---
class: center, middle, title-slide

# Differential Equation Systems in R

---
class: page-font-large

## Import and Export

`\begin{align}
\frac{dN}{dt} = \mbox{import} - \mbox{export}\\
\end{align}`

<div id="htmlwidget-675062afa5c58d169f2a" style="width:1008px;height:288px;" class="grViz html-widget"></div>
<script type="application/json" data-for="htmlwidget-675062afa5c58d169f2a">{"x":{"diagram":"digraph feedback {\n         graph [rankdir = \"LR\", bgcolor=\"none\"]\n           node [shape = box, penwidth=2, fontname = \"Helvetica\"]\n             N\n           node [shape = octagon, penwidth=0.5, style=\"rounded\", fixedsize=25, fontsize=8]\n             Source Sink\n           node [shape = none, fontsize=10]\n             import export\n           edge [penwidth=1.5]\n             Source -> import -> N -> export -> Sink\n         \n}","config":{"engine":"dot","options":null}},"evals":[],"jsHooks":[]}</script>

`$N$` = .gray[any quantity: matter, energy, information], e.g. **N**umber of individuals

---

## Exponential Growth

`\begin{align}
\frac{dN}{dt} &= \mbox{birth} - \mbox{death}\\
              & = \mbox{birthrate} \cdot N - \mbox{deathrate}\cdot N\\
              &= b \cdot N - d \cdot N = (b-d) \cdot N \\
\frac{dN}{dt} &= r \cdot N\\
\end{align}`

<div id="htmlwidget-68ee93127015cd9c35cc" style="width:1008px;height:288px;" class="grViz html-widget"></div>
<script type="application/json" data-for="htmlwidget-68ee93127015cd9c35cc">{"x":{"diagram":"digraph feedback {\n         graph [rankdir = \"LR\", bgcolor=\"none\"]\n           node [shape = box, penwidth=2, fontname = \"Helvetica\"]\n             \"Population N\"\n           node [shape = octagon, penwidth=0.5, style=\"rounded\", fixedsize=25, fontsize=8]\n             Source Sink\n           node [shape = none, fontsize=10]\n             birth death\n           edge [penwidth=1.5]\n             Source -> birth -> \"Population N\" -> death -> Sink\n           edge [penwidth=0.7, tailport = \"n\", headport = \"n\", constraint = false, color=tomato]\n             \"Population N\" -> birth \n             \"Population N\" -> death\n}","config":{"engine":"dot","options":null}},"evals":[],"jsHooks":[]}</script>

---

## Analytical Integration

.pull-left[
![](index_files/figure-html/expgrowth-1.png)
```r
r    <- 0.5
N0   <- 10
dt   <- 0.1
time <- seq(0, 10, dt)

N <- N0 * exp(r * time)
plot(time, N, type="l")
```
]

.pull-right[
`\begin{align}
\frac{dN}{dt} &= r\cdot N\\
\int_0^t \frac{dN}{dt} &= \int_0^t r\cdot N\\
\int_0^t \frac{1}{N} dN &= \int_0^t r dt\\
\ln(N) &= rt + c\\
N_t &= N_0 e^{rt}\\
\end{align}`

* Integral calculus = "Gold Standard" ...
* ... but not always possible
]

---

## Numerical Integration with the Euler Method

![](index_files/figure-html/expgrowth-euler-1.png)

```r
N <- numeric(length(time))

N[1] <- N0
for (i in 2:length(time)) {
  N[i] <- N[i-1] + r * N[i-1] * dt
}
plot(time, N, type = "l", col = "red")
```
]

.pull-right[
`\begin{align}
\frac{dN}{dt} &= r\cdot N\\
\frac{\Delta N}{\Delta t} &= r\cdot N\\
N_t &= N_{t-1} + r\cdot N_{t-1} \cdot \Delta t
\end{align}`

.page-font-large[
* .red[discrete time steps]
* approximate differential equation<br> by `$\rightarrow$` difference equation
* easy, but integration error
]
]

---

## Numerical Integration Made Easy with deSolve

#### Code: close to mathematical notation

```r
library("deSolve")
model <- function (time, y, parms) {
  with(as.list(c(y, parms)), {
*   dN_dt <- r * N * (1 - N / K)
    list(c(dN_dt))
  })
}
y0 <- c(N = 0.1)
parms <- c(r = 0.1, K = 10)
times <- seq(from = 0, to = 100, by = 1)
```
]

```r
*out <- ode(y0, times, model, parms)
plot(out)
```

![](index_files/figure-html/ode-naive2-1.png)

* better performance, reduced integration error
]

---
## R Package deSolve

* **initial value problems**, .gray[ODE, PDE, DAE, DDE]
* popular .green[industry class solvers] from ODEPACK<sup>1</sup><br>
  .gray[and classical Euler and Runge-Kutta solvers]
* .green[convenience tools] (e.g. `plot`, forcings, events) 
]

* packages: rootSolve, bvpSolve, FME, ReacTran, ...
* .green[compatible interfaces]
* See [CRAN Task View](https://cran.r-project.org/web/views/DifferentialEquations.html) for related packages
* .green[extensive documentation:] papers, books, websites,
[mailing list](https://stat.ethz.ch/mailman/listinfo/r-sig-dynamic-models), 
[StackOverflow](https://stackoverflow.com/search?tab=newest&q=%5br%5d%20desolve?ode)
]

.gray[
<br><br><br>
`$^1$` Thanks to the great ODEPACK authors (ctb=contributors) who made their code publicly available:

Peter N. Brown [ctb] (files ddaspk.f, dvode.f, zvode.f), George D. Byrne [ctb] (files dvode.f, zvode.f), Ernst Hairer [ctb] (files radau5.f, radau5a), Alan C. Hindmarsh [ctb] (files ddaspk.f, dlsode.f, dvode.f, zvode.f, opdkmain.f, opdka1.f), Cleve Moler [ctb] (file dlinpck.f), Linda R. Petzold [ctb] (files ddaspk.f, dlsoda.f), Youcef Saad [ctb] (file dsparsk.f), Clement W. Ulrich [ctb] (file ddaspk.f)
]

---
class: page-font-large

## What is `ode`?

* **`ode`** is used to call one of the solvers
* then the solver calls the model function at adequate time steps
* the model function contains the **derivatives** of the system

## Solvers
* default: .red[**`lsoda`**]
* Livermore solver for Ordinary Differential Equations, automatic
* switches between stiff and non-stiff solver

... other solvers, fine-tuning and pitfalls, see [docs](https://cran.r-project.org/package=deSolve/)

---
## Automatic Step Size

Most solvers adapt integration steps automatically,<br> to maintain pre-defined precision.

User can specify tolerances **`atol`**, **`rtol`**, default: `$10^{-6}$`

Outputs collected at pre-specified time points.

#### Fixed-step solvers: `euler`, `rk2`, `rk4`

still quite popular, because easy to implement

`$\rightarrow$` problems with precision, CPU effort, stability<br>
`$\rightarrow$` useful in special cases

]

```r
out <- ode(y0, 
* times = 1:10,    # output time steps
  model, 
  parms, 
  atol = 1e-6,     # absolute tolerance
  rtol = 1e-6      # relative tolerance
)
```

<br><br>
* use .red[**`lsoda`**]
* no need to use very small output steps
* keep tolerances in mind,<br> consider to rescale huge or tiny values
]

---
class: page-font-huge
background-image: url("img/water-light.jpg")

# Some Practical Examples

1. Resource limited growth

1. A basic Covid 19 model

1. Basic idea of a lake: the chemostat

1. Lotka-Volterra's predator-prey interactions

1. Reaction and transport in a river

---
background-image: url("img/batch.jpg")

## .bg-contrast[Example 1: Resource Limited Growth]

---
## Resource Limited Growth

<br>

<div id="htmlwidget-b22d9321ec32079dfd79" style="width:1008px;height:216px;" class="grViz html-widget"></div>
<script type="application/json" data-for="htmlwidget-b22d9321ec32079dfd79">{"x":{"diagram":"digraph feedback {\n         graph [rankdir = \"LR\", bgcolor=\"none\"]\n           node [shape = box, penwidth=2, fontname = \"Helvetica\"]\n             Nutrient, Algae\n           node [shape = none, fontsize=10]\n             growth\n           edge [penwidth=1.5]\n             Nutrient -> growth -> Algae \n           edge [penwidth=0.7, tailport = \"n\", headport = \"n\", constraint = false, color=tomato]\n             Algae -> growth \n}","config":{"engine":"dot","options":null}},"evals":[],"jsHooks":[]}</script>

* two state variables
* functional response: growth depends on nutrient, e.g. Holling, Monod, ...
* feedback: more algae, more growth

---

### System of equations

<div id="htmlwidget-13b39a11f84cdfb39b12" style="width:1008px;height:144px;" class="grViz html-widget"></div>
<script type="application/json" data-for="htmlwidget-13b39a11f84cdfb39b12">{"x":{"diagram":"digraph feedback {\n         graph [rankdir = \"LR\", bgcolor=\"none\"]\n           node [shape = box, penwidth=2, fontname = \"Helvetica\"]\n             Nutrient, Algae\n           node [shape = none, fontsize=10]\n             growth\n           edge [penwidth=1.5]\n             Nutrient -> growth -> Algae \n           edge [penwidth=0.7, tailport = \"n\", headport = \"n\", constraint = false, color=tomato]\n             Algae -> growth \n}","config":{"engine":"dot","options":null}},"evals":[],"jsHooks":[]}</script>

----

`\begin{align}
\frac{d\mbox{Algae}}{dt} &=       & &r \cdot f(\mbox{Nutrient}) \cdot \mbox{Algae}\\
\\
\frac{d\mbox{Nutrient}}{dt} &=   &- & r \cdot f(\mbox{Nutrient}) \cdot  \mbox{Algae} \cdot \frac{1}{Y} 
\end{align}`

with:

$$
f(\mbox{Nutrient}) =  \text{functional response, e.g. Holling I, II, III}
$$

---

## Resource Limited Growth

#### The Model

```r
model <- function (time, y, parms) {
  with(as.list(c(y, parms)), {
    f <- P/(kP + P)
    dAlg_dt <- r * f * Alg
    dP_dt   <- - r * 1/Y * f * Alg
    list(c(dAlg_dt, dP_dt))
  })
}
y <- c(
  Alg = 10,  # algae,      mol/m3 carbon
  P   =  5   # phosphorus, mol/m3 P
) # in mg/L
parms <- c(r = 0.1, kP = 0.5, Y = 106) # Y = C:P - Redfield ratio
```

#### Where do parameters come from?

* `$r$`, `$k_P$` derived from lab experiments
* `$Y$` first guess from stoichiometry

```r
library("marelac")
redfield(1, species="P")
```

```
##     C   H   O  N P
## 1 106 263 110 16 1
```

#### Simulate the Model

```r
out <- ode(y, times, model, parms)
plot(out, las = 1)
```

![](index_files/figure-html/unnamed-chunk-5-1.png)

#### Compare Scenarios

```r
out1 <- ode(y, times, model, parms =  c(r = 0.2, kP = 0.5, Y = 106))
out2 <- ode(y, times, model, parms =  c(r = 0.1, kP = 1.0, Y = 120))

plot(out1, out2, las = 1)
legend("topright", legend=c("Scenario 1", "Scenario 2"), lty=1:2, col=1:2)
```

![](index_files/figure-html/unnamed-chunk-6-1.png)

---
background-image: url("img/1280px-SARS-CoV-2-CDC-23312.jpg")

## .bg-contrast[Example 2: A Covid-19 Model]

---

## A standard SIR model

<div id="htmlwidget-04af221b426bc0303f25" style="width:1008px;height:288px;" class="grViz html-widget"></div>
<script type="application/json" data-for="htmlwidget-04af221b426bc0303f25">{"x":{"diagram":"digraph feedback {\n         graph [rankdir = \"LR\", bgcolor=\"none\"]\n           node [shape = box, penwidth=2, fontname = \"Helvetica\", style=\"filled\", color=\"dodgerblue\"]\n             S, I, R\n           node [shape = none, fontsize=10, style=\"\"]\n             infection, recovery, mortality\n            node [shape = octagon, penwidth=0.5, style=\"rounded\", fixedsize=25, fontsize=8]\n             X\n           edge [penwidth=1.5]\n             S -> infection -> I\n             I -> recovery -> R\n             I -> mortality\n             mortality -> X\n}","config":{"engine":"dot","options":null}},"evals":[],"jsHooks":[]}</script>

#### Total population is subdivided into 3 subpopulations

* `$S$`: susceptible, `$I$`: infected, `$R$`: recovered

This is the general scheme from which different modifications can be derived,
e.g. an additional state variable `$E$` (SEIR model) of exposed, or of deceased individuals `$X$` (SIRX). It is also possible to consider multiple groups of `$S$`, `$I$`, `$R$`, their spatial
distribution or influence of stochasticity.

---

### Equations describing spread of the disease:

`\begin{align}
\frac{dS}{dt} &= -infection\\
\frac{dI}{dt} &= infection - recovery - mortality\\
\frac{dR}{dt} &= recovery
\end{align}`

where

`\begin{align}
infection &= b \cdot I \cdot S\\
recovery  &= g \cdot I\\
mortality &= m \cdot I
\end{align}`

| Name      | Value      | Description         | Unit                        |
| --------- | ---------- | ------------------- | --------------------------- |
| b         | 0.00000002 | infection parameter | `$\mathrm{ind^{-1}~d^{-1}}$` |
| g         | 0.07       | recovery parameter  | `$\mathrm{d^{-1}}$`           |
| m         | 0.007      | mortality parameter | `$\mathrm{d^{-1}}$`           |

---

### Code

```r
SIR <- function(t, state, parameters) {
  with (as.list(c(state, parameters)), {
    infection <- b * S * I
    recovery  <- g * I
    mortality <- m * I

dS_dt     <- -infection 
    dI_dt     <-  infection - recovery - mortality
    dR_dt     <-  recovery

list(c(dS_dt, dI_dt, dR_dt),  # the time derivatives
         Population = S+I+R)      # extra output variable
  })
}

parms <- c(
  b = 0.00000002,    # [1/ind/d], infection parameter 
  g = 0.07,          # [1/d],     recovery rate of infected individuals
  m = 0.007          # [1/d],     mortality rate of infected individuals
)
y0 <- c(S = 11500000 - 1000, I = 1000, R = 0)
times <- seq(from = 0, to = 365, by = 1)   # time sequence, in days
out   <- ode(y = y0, times = times, func = SIR, parms = parms)
plot(out, las = 1, mfrow = c(1, 4))
```

![](index_files/figure-html/SIR-solution-1.png)

---

## External Interventions

* Let's assume society implements social distancing **after** realizing exponential growth.
* Implementation makes use of .green[forcings].
* This allows to consider any time-dependent inputs.

```r
SIR_distancing <- function(t, state, parameters) {
  with (as.list(c(state, parameters)), {
*   contact <- f_contact(t)
    Infection <- contact * b * S * I
    Recovery  <- g * I
    Mortality <- m * I
    
    dS_dt     <- -Infection 
    dI_dt     <-  Infection - Recovery - Mortality
    dR_dt     <-  Recovery

list(c(dS_dt, dI_dt, dR_dt),  # the time derivatives
         Population = S+I+R)      # extra output variable
  }) 
}

*f_contact <- approxfun(x = c(0, 50,  51, 365), y = c(1,  1, 0.5, 0.5), rule = 2)

y0 <- c(S = 11500000 - 1000, I = 1000, R = 0)

parms <- c(
  b = 0.00000002,    # [1/ind/d], infection parameter 
  g = 0.07,          # [1/d],     recovery rate of infected individuals
  m = 0.007          # [1/d],     mortality rate of infected individuals
)

times <- seq(from = 0, to = 365, by = 1)   # time sequence, in days
out2   <- ode(y = y0, times = times, func = SIR_distancing, parms = parms)
plot(out, out2, las = 1, which = c("S", "I", "R"), mfrow = c(1, 4))
plot(times, rep(1, length(times)), 
     ylim = c(0,1), type="l", xlab = "time", ylab = "", 
     main = "social contact", las = 1, lty = 2)
lines(times, f_contact(times), lty = 2, col="red")
```

![](index_files/figure-html/SIR-distancing-1.png)

---
background-image: url("img/neunzehnhain.jpg")

## .bg-contrast[Example 3: A Lake with Inflow and Outflow]

---

## The Lake as a Bioreactor

![Chemostat](img/chemostat.svg)

* Consider a chemostat bioreactor as a basic model of a lake:
* Phosphorus (P) is delivered from a source with a constant *dilution* rate **D**.
* Algae (Alg) take up phosphorus and grow.
* Part of algae and remaing phosphorus are flushed away.

---

## A Chemostat with Algae and Phosphorus

`\begin{align}
\frac{dAlg}{dt}  & = growth - export \\
\frac{dP}{dt}    & = (import - export) - 1/Y \cdot growth  
\\
\text{with:}\\
growth & = f(P)
\end{align}`

---

`\begin{align}
\frac{dAlg}{dt}  & = r \cdot Alg - D \cdot Alg\\
\frac{dP}{dt}    & =  D \cdot (P_0 - P) - 1/Y \cdot r \cdot Alg\\
\\
\text{with:}\\
r   & = r_{max} \cdot \frac{P}{kp + P}  \qquad\text{(Monod equation)}\\
\end{align}`

---

## The Chemostat Code

```r
library("deSolve")
library("rootSolve")

chemostat <- function(time, init, parms) {
  with(as.list(c(init, parms)), {
    r       <- r_max * P/(kp + P)  # Monod equation
    dAlg_dt <- r * Alg - D * Alg
    dP_dt   <- D * (P0 - P) - 1/Y * r * Alg
    list(c(dAlg_dt, dP_dt), r = r)
   })
}
parms <- c(
  r_max = 0.5,    # 1/d
  kp    = 0.5,    # half saturation constant, P (mol/m3)
  Y     = 106,    # yield coefficient (stoichiometric C:P ratio)
  D     = 0.1,    # 1/d
  P0    = 5       # P in inflow (mol/m3)
)
times <- seq(0, 40, 0.1)        # (d)
init  <- c(Alg = 10, P = 5)    # Phytoplankton C and Phosphorus P (mol/m3)

## =============================================================================
## Dynamic simulation
## =============================================================================
out <- ode(init, times, chemostat, parms)
plot(out, mfrow=c(1, 3), las = 1)
```

![](index_files/figure-html/Chemostat-dynamic-1.png)

---

### Equilibrium

#### Steady state solution

```r
state <- data.frame(
  D = seq(0, 0.6, length.out = 100),
  X = 0,
  S = 0
)

for (i in 1:nrow(state)) {
  parms["D"] <- state$D[i]
  times <- c(0, Inf)
  out <- runsteady(init, times, chemostat, parms)
  state[i, 2:3] <- out$y
}

par(mfrow = c(3, 1))
plot(S ~ D, data = state, type = "l")
plot(X ~ D, data = state, type = "l")
plot(S * X ~ D, data = state, type = "l")
```

![](index_files/figure-html/chemostat-equilibrium-1.png)

**Note:** this can also be solved analytically.

---
background-image: url("img/marmots.jpg")

## .bg-contrast[Example 4: Predator-Prey Model]

---

## Interaction Between Populations

### Lotka-Volterra's Predator and Prey

<div id="htmlwidget-112d70e144d5b14540ca" style="width:1008px;height:288px;" class="grViz html-widget"></div>
<script type="application/json" data-for="htmlwidget-112d70e144d5b14540ca">{"x":{"diagram":"digraph lotka {\n         graph [rankdir = \"LR\", bgcolor=\"none\"]\n           node [shape = box, penwidth=2, fontname = \"Helvetica\"]\n             Predator, Prey\n           node [shape = octagon, penwidth=0.5, style=\"rounded\", fixedsize=25, fontsize=8]\n             Source Sink\n           node [shape = none, fontsize=10]\n             growth death grazing\n           edge [penwidth=1.5]\n             Source -> growth -> Prey -> grazing -> Predator -> death -> Sink\n           edge [penwidth=0.7, tailport = \"n\", headport = \"n\", constraint = false, color=tomato]\n             Prey -> growth\n             Predator -> grazing\n           edge [penwidth=0.7, tailport = \"s\", headport = \"s\", constraint = false, color=tomato]\n             Prey -> grazing\n             Predator -> death\n\n}","config":{"engine":"dot","options":null}},"evals":[],"jsHooks":[]}</script>

* state diagram, cycles
* damped cycle if we introduce a Monod term

---

### Lotka-Volterra Model

<div id="htmlwidget-5e3ee158686e707e54b2" style="width:1008px;height:144px;" class="grViz html-widget"></div>
<script type="application/json" data-for="htmlwidget-5e3ee158686e707e54b2">{"x":{"diagram":"digraph lotka {\n         graph [rankdir = \"LR\", bgcolor=\"none\"]\n           node [shape = box, penwidth=2, fontname = \"Helvetica\"]\n             Predator, Prey\n           node [shape = octagon, penwidth=0.5, style=\"rounded\", fixedsize=25, fontsize=8]\n             Source Sink\n           node [shape = none, fontsize=10]\n             growth death grazing\n           edge [penwidth=1.5]\n             Source -> growth -> Prey -> grazing -> Predator -> death -> Sink\n           edge [penwidth=0.7, tailport = \"n\", headport = \"n\", constraint = false, color=tomato]\n             Prey -> growth\n             Predator -> grazing\n           edge [penwidth=0.7, tailport = \"s\", headport = \"s\", constraint = false, color=tomato]\n             Prey -> grazing\n             Predator -> death\n\n}","config":{"engine":"dot","options":null}},"evals":[],"jsHooks":[]}</script>

`\begin{align}
\frac{dPrey}{dt} &= \mbox{growth} - \mbox{grazing}\\
\\
\frac{dPred}{dt} &= g \cdot \mbox{grazing} - \mbox{mortality}
\end{align}`

With:

`\begin{align}
\mbox{growth}    & = a \cdot Prey & \qquad \text{(exponential growth)}\\
\mbox{grazing}   & = b \cdot Prey \cdot Pred & \qquad \text{(interaction)}\\
\mbox{mortality} & = -e \cdot Pred & \qquad \text{(exponential decay)}\\
                 &  g              & \text{(grazing efficiency)}
\end{align}`

---

### Lotka-Volterra Model: Implementation

```r
require(deSolve)   
y0    <- c(Prey=300, Pred=10)                       # state variable initial conditions
parms <- c(a=0.05, K=500, b=0.0002, g=0.8, e=0.03)  # parameter values

LV <- function(t, state, parameters) { 
  with (as.list(c(state, parameters)), {
    growth     <- a * Prey
    grazing    <- b * Prey * Pred
    mortality  <- e * Pred
    dPrey_dt <- growth - grazing
    dPred_dt <- g * grazing - mortality
    return (list(c(dPrey_dt, dPred_dt),  # vector of derivatives
                 sum = Prey + Pred))     # output variable
  })
}
times <- 1:1000
out   <- ode(y = y0, func = LV, times = times, parms = parms)
plot(out, mfrow = c(1, 3), las = 1)
```

![](index_files/figure-html/LV-basic-1.png)

---

## Efficient Model Formulation in Matrix Style

#### Two preys and two preadors `$\rightarrow$` 4 equations

```r
model <- function(t, n, parms) {
  with(as.list(c(n, parms)), {
    dN1_dt <-  r1 * N1      - a13 * N1 * N3  # Prey1
    dN2_dt <-  r2 * N2      - a24 * N2 * N4  # Prey2
    dN3_dt <- a13 * N1 * N3 -  r3 * N3       # Predator1
    dN4_dt <- a24 * N2 * N4 -  r4 * N4       # Predator2
    return(list(c(dN1_dt, dN2_dt, dN3_dt, dN4_dt)))
  })
}

times <- seq(from=0, to=500, by = 0.1)
n0    <- c(N1=1, N2=1, N3=2, N4=2)            # Number of individuals
parms <- c(r1 = 0.1, r2=0.1, r3=0.1, r4=0.1,  # net growth = birth - death
          a13=0.2, a24 = 0.1                  # interaction
)

out <- ode(n0, times, model, parms)
```

#### Matrix formulation: appears as one single equation

```r
model <- function(t, N, parms) {
  with(parms, {
*   dN_dt <- r * N  + N * (A %*% N)
    return(list(dN_dt))
  })
}

parms <- list(
  r = c(r1 = 0.1, r2 = 0.1, r3 = -0.1, r4 = -0.1),
  ## pairwise interactions:
  A = matrix(c(0.0, 0.0, -0.2,  0.0,    # prey 1
               0.0, 0.0,  0.0, -0.1,    # prey 2
               0.2, 0.0,  0.0,  0.0,    # predator 1; eats prey 1
               0.0, 0.1,  0.0,  0.0),   # predator 2; eats prey 2
               nrow = 4, ncol = 4, byrow = TRUE)
)

out <- ode(n0, times, model, parms)
plot(out, las = 1)
```

![](index_files/figure-html/unnamed-chunk-8-1.png)

---
background-color: #00305d
background-image: url("img/polluted-river3.jpg")

## .orange[Example 5:<br> Reaction<br> and Transport<br> in a<br> Polluted River]

---

## PDE's in 1D: Transport in a River

* A 1D river model has in fact two dimensions, time and space.
* The model has derivatives in both dimensions, so we get partial differential equations (PDE).
* Different methods exist to solve such systems -- **deSolve** supports the **method of lines**:
    * space is discretisized into fixed grid cells
    * time is handled by standard ODE solvers with variable time step
    * functions `ode.1D`, `ode.2D` and `ode.3D` provide pre-defined Jacobi matrices
* Transport with own matrix functions or with **R** package **ReacTran**
* System of equations for all grid cells forms a state vector with, e.g. several 1000 equations,
that are then solved quasiparallel by the solver.

Fast solvers + efficient estimation the Jacobians + R's vectorization<br> make multidimensional models surprinsingly fast.

---

## Example: Nitrification in a River

<div id="htmlwidget-52ad96ed96f954092277" style="width:1008px;height:216px;" class="grViz html-widget"></div>
<script type="application/json" data-for="htmlwidget-52ad96ed96f954092277">{"x":{"diagram":"digraph river_1d {\n         graph [rankdir = \"LR\", bgcolor=\"none\"]\n           node [shape = box, penwidth=2, fontname = \"Helvetica\"]\n             O2, NH4, NO3\n           node [shape = octagon, penwidth=0.5, style=\"rounded\", fixedsize=25, fontsize=8]\n             Source\n           node [shape = none, fontsize=10]\n             nitrification aeration\n           edge [penwidth=1.5]\n             Source -> aeration -> O2\n             NH4 -> nitrification\n             O2 -> nitrification\n\n             nitrification -> NO3\n           edge [penwidth=0.7, tailport = \"s\", headport = \"s\", constraint = false, color=tomato]\n             O2 -> aeration\n}","config":{"engine":"dot","options":null}},"evals":[],"jsHooks":[]}</script>

---

## Model in 0D: Batch

```r
parms <- c(
  r_nitri = 2.0,   # 1/day   nitrification rate
  kO2     = 0.001, # mol/m3  half sat. constant
  k2      = 2.8,   # 1/day  re-aeration rate
  O2sat   = 0.3    # mol/m3 saturation concentration
)

Nitrification <- function(t, y, parms) {
  with(as.list(c(y, parms)), {

aeration      <- k2 * (O2sat - O2)
    nitrification <- r_nitri * O2/(O2 + kO2) * NH4

dNH4 <- -nitrification
    dNO3 <- +nitrification
    dO2  <-  aeration - 2 * nitrification
    list(c(dNH4, dNO3, dO2))
  })
}

y0 <- c(NH4 = 0.5, NO3 = 0, O2 = 0.3) # in mol/m3

times <- seq(from = 0, to = 5, length.out=101)

out <- ode(y=y0, times = times, func=Nitrification, parms = parms)
plot(out, mfrow=c(1, 3), lwd = 2, las = 1, ylab="mol/m3")
```

![](index_files/figure-html/river0d-1.png)

---

## River Geometry

<div id="htmlwidget-650ee0c33979d9506094" style="width:1008px;height:288px;" class="grViz html-widget"></div>
<script type="application/json" data-for="htmlwidget-650ee0c33979d9506094">{"x":{"diagram":"digraph river_nh4 {\n         graph [rankdir = \"LR\", bgcolor=\"none\"]\n           node [shape = box, penwidth=1.5, width=0.25, fontname = \"Helvetica\", fontsize=16]\n           1 2 3 4 5 6 7  N\n           node [shape = none]\n           up \"...\" down\n           edge [penwidth=1.2, len=1]\n           up -> 1 -> 2 -> 3 -> 4 -> 5 -> 6 -> 7  -> \"...\" -> N -> down\n     \n}","config":{"engine":"dot","options":null}},"evals":[],"jsHooks":[]}</script>

* 100km river stretch, 100 grid cells, each 1km
* upper and lower boundary conditions for all substances
* for simplicity: constant river geometry, flow velocity and rate parameters

---

## Implementation of River Geometry and Transport

```r
library("ReacTran")
Nitrification <- function(t, y, parameters) {
  with(as.list(c(parameters)),{
    NH4 <- y[1:N]
    NO3 <- y[(N+1):(2*N)]
    O2  <- y[(2*N+1):(3*N)]

aeration      <- k2 * (O2sat - O2)
    nitrification <- r_nitri * O2/(O2 + kO2) * NH4

dNH4 <- -nitrification +
      tran.1D(C = NH4, D = D, v = v, C.up = NH4up, C.down = NH4dwn, A = A, dx = Grid)$dC

dNO3 <- +nitrification +
      tran.1D(C = NO3, D = D, v = v, C.up = NO3up, C.down = NO3dwn, A = A, dx = Grid)$dC

dO2  <-  aeration - 2 * nitrification +
      tran.1D(C = O2 , D = D, v = v, C.up = O2up,  C.down = O2dwn,  A = A, dx = Grid)$dC

list(c(dNH4, dNO3, dO2))
  })
}
```

* `N` grid cells
* state vector: `$\mathrm{NH4_{1..100}, NO3_{101..200}, O2_{201..300}}$`
* upper and lower boundary conditions
* river transect area `A` 
* longitudinal Diffusion `D`
* grid cell length `dx`

---

## Full Source Code .gray[(with source code for animations)]

```r
library("ReacTran")
library("deSolve")
library("rootSolve")
library("gifski")

## Measurement units used
# Time   = day
# Mass   = mol
# Space  = m
# Volume = m3

N    <- 100  # number of grid cells
cell <- 1000 # 1000m = 1km / grid cell
Grid <- setup.grid.1D(N = N, L = N * cell)

#=============================================================================
# Model parameters
#=============================================================================

## river characteristics
v  <- 40000       # m/day
w  <- 15          # m
h0 <- 0.35        # m
A  <- w * h0      # m2

parameters <- c(
  r_nitri       = 2.0,   # 1/day   nitrification rate
  kO2           = 0.001, # mol/m3  half sat. constant
  NH4up         = 0.5,   # mol/m3  upper boundary
  NO3up         = 0.0,   # mol/m3  -- " --
  O2up          = 0.3,   # mol/m3  -- " --
  NH4dwn        = 0.0,   # mol/m3  lower boundary
  NO3dwn        = 0.0,   # mol/m3  -- " --
  O2dwn         = 0.3,   # mol/m3  -- " --
  D             = 100,   # m2/day dispersion coefficient
  k2            = 2.8,   # 1/day  re-aeration rate
  O2sat         = 0.3    # mol/m3 saturation concentration
)

#------------------------------------------------------------------------------
# Model definition:
#------------------------------------------------------------------------------

Nitrification <- function(t, y, parameters) {
  with(as.list(c(parameters)),{
    NH4 <- y[1:N]
    NO3 <- y[(N+1):(2*N)]
    O2  <- y[(2*N+1):(3*N)]

aeration      <- k2 * (O2sat - O2)
    nitrification <- r_nitri * O2/(O2 + kO2) * NH4

dNH4 <- -nitrification +
      tran.1D(C = NH4, D = D, v = v, C.up = NH4up, C.down = NH4dwn, A = A, dx = Grid)$dC

dNO3 <- +nitrification +
      tran.1D(C = NO3, D = D, v = v, C.up = NO3up, C.down = NO3dwn, A = A, dx = Grid)$dC

dO2  <-  aeration - 2 * nitrification +
      tran.1D(C = O2 , D = D, v = v, C.up = O2up,  C.down = O2dwn,  A = A, dx = Grid)$dC

list(c(dNH4, dNO3, dO2))
  })
}

#=============================================================================
# Steady-state solution
#=============================================================================

std <- steady.1D(y = runif(3*Grid$N), names=c("NH4", "NO3", "O2"),
   func = Nitrification, parms = parameters, nspec = 3, pos = TRUE)

plot(std, mfrow = c(3, 1), las=1, xlab="km", ylab="mol/m3")

#=============================================================================
# Dynamic solution
#=============================================================================

# State variables and initial conditions
NH4_0  <- 0.0 # mol/m3
NO3_0  <- 0.0 # mol/m3
O2_0   <- 0.3 # mol/m3

state <- c(rep(NH4_0, N), rep(NO3_0, N), rep(O2_0, N))

# Time sequence
time_seq <- seq(from = 0, to = 5, length.out=101)

out1D  <- ode.1D(y=state, times=time_seq, func=Nitrification, parms=parameters, nspec=3)
image(out1D, xlab = "time, days", ylab = "Distance, m", grid= Grid$x.mid,
      main =c("NH4","NO3","O2"), add.contour=TRUE, mfrow = c(1, 3))

#=============================================================================
# Animated figure
#=============================================================================
plot_poly <- function(data, time) {
  poly <- function(x, y, pcol, ylab, ylim) {
    plot(x, y, type="n", las=1, xlab="", ylab="", ylim=ylim)
    polygon(x, y, col=pcol, lty="blank")
    mtext(side=1, line=1.5, text="km")
    mtext(side=2, line=3.5, text=ylab, las=3)

}
  y <- data[time,-1]
  NH4 <- y[1:N]
  NO3 <- y[(N+1):(2*N)]
  O2  <- y[(2*N+1):(3*N)]
  x   <- c(1, 1:N, N)

poly(x, c(0, NH4, 0), ylim=c(0, 0.5), pcol="#a6cee3", ylab="NH4 (mol/m3)")
  poly(x, c(0, NO3, 0), ylim=c(0, 0.5), pcol="#b2df8a", ylab="NO3 (mol/m3)")
  poly(x, c(0,  O2, 0), ylim=c(0, 0.3), pcol="#1f78b4", ylab="O2 (mol/m3)")
}

for (tt in 1:length(time_seq))  {
  png(paste0("river", 1000 + tt, ".png"), width=1600, height=800, pointsize = 18)
  par(mfrow=c(3,1))
  par(mar=c(3,5,.5,0), las=1, cex.axis=1.4, cex.lab=1.4, cex.main=2)
  plot_poly(out1D, tt)
  dev.off()
}
gifski(dir(pattern="^river.*png$"), gif_file = "river.gif", width=1600, height=800, delay=1/25)
```

---

### Nitrification in a River (1D): Result
.center[
![:scale 90%](img/river.gif)
]

---
class: middle, title-slide, page-font-large
background-image: url("img/water_tud_nioz_r_oer_by.jpg")
background-size: cover

# Conclusions

* User-contributed packages for **R** to solve differential equations numerically

* Easy to step in ...

* ... supports also complex applications.

* Part of a comprehensive "package ecosystem".

---
background-image: url("img/outlook.jpg")

# .bg-contrast[Outlook]

---

## Further reading ...

The details are described in the package vignettes of [deSolve](https://CRAN.R-project.org/package=deSolve) and [ReacTran](https://CRAN.R-project.org/package=ReacTran) 
and in the following papers: 
<a name=cite-Soetaert2010b></a>[Soetaert, Petzoldt, and Setzer (2010b)](https://doi.org/10.18637/jss.v033.i09), 
<a name=cite-Soetaert2010_rjournal></a>[Soetaert, Petzoldt, and Setzer (2010c)](https://doi.org/10.32614/RJ-2010-013), 
<a name=cite-ReacTran></a>[Soetaert and Meysman (2010)](#bib-ReacTran) and books: 
<a name=cite-Soetaert2008></a>[Soetaert and Herman (2009)](#bib-Soetaert2008), 
<a name=cite-Soetaert2012></a>[Soetaert, Cash, and Mazzia (2012)](#bib-Soetaert2012).

The package relies on the source code and publications of the original algorithms, especially:
<a name=cite-Hindmarsh1983></a>[Hindmarsh (1983)](#bib-Hindmarsh1983), 
<a name=cite-Petzold1983></a>[Petzold (1983)](#bib-Petzold1983), 
<a name=cite-Cash1990></a>[Cash and Karp (1990)](#bib-Cash1990), 
<a name=cite-Hairer2009></a>[Hairer, Norsett, and Wanner (2009)](#bib-Hairer2009) and 
<a name=cite-Hairer2010></a>[Hairer and Wanner (2010)](#bib-Hairer2010).

### Online Resources

* The CRAN Task View [Differential Equations](https://cran.r-project.org/web/views/DifferentialEquations.html) gives an overview over other packages with differential equation solvers and related tools.

* The [deSolve development page](https://desolve.r-forge.r-project.org) contains links to more detailed tutorials, books and papers.

A bunch of free docs and examples can be found in the internet, just search for `deSolve` and `R`.

---
class: center, middle, title-slide
background-image: url("img/water-light.jpg")

## Check out our github repos at:

# https://github.com/dynamic-R

## and thanks for watching.

---

## Copyright

This resource was created by [tpetzoldt](github.com/tpetzoldt) and 
[karlines](https://github.com/karlines). It is provided 
as is without warranty.

If not otherwise stated, images are own photos and drawings.

The original drawings of <a name=cite-Forrester1990></a>[Forrester (1971b)](https://web.mit.edu/sysdyn/sd-intro/D-4165-1.pdf) are cited, the [Newton](https://commons.wikimedia.org/wiki/File:IsaacNewton_1642-1727.jpg) (CC-BY-SA 4.0), 
[Leibniz](https://commons.wikimedia.org/wiki/File:Gottfried_Wilhelm_Leibniz,_Bernhard_Christoph_Francke.jpg) (public domain)
and 
[SARS-CoV-2](https://commons.wikimedia.org/wiki/File:SARS-CoV-2-CDC-23312.png) (public domain)
images were provided by Wikimedia Commons under a creative commons license. 
Please follow the link for details.

---

## Bibliography

<a name=bib-Cash1990></a>[Cash, J. R. and A. H. Karp](#cite-Cash1990)
(1990). "A Variable Order Runge-Kutta Method for Initial Value Problems
With Rapidly Varying Right-Hand Sides". In: _ACM Transactions on
Mathematical Software_ 16, pp. 201-222.

<a name=bib-Forrester1990></a>[Forrester, J.](#cite-Forrester1990)
(1971b). _Some Basic Concepts in System Dynamics_. Sloan School of
Management, Massachusetts Institute of Technology. URL:
[https://web.mit.edu/sysdyn/sd-intro/D-4165-1.pdf](https://web.mit.edu/sysdyn/sd-intro/D-4165-1.pdf).

<a name=bib-Hairer2009></a>[Hairer, E., S. P. Norsett, and G.
Wanner](#cite-Hairer2009) (2009). _Solving Ordinary Differential
Equations I: Nonstiff Problems. Second Revised Edition_. Heidelberg:
Springer-Verlag.

<a name=bib-Hairer2010></a>[Hairer, E. and G. Wanner](#cite-Hairer2010)
(2010). _Solving Ordinary Differential Equations II: Stiff and
Differential-Algebraic Problems. Second Revised Edition_. Heidelberg:
Springer-Verlag.

<a name=bib-Hindmarsh1983></a>[Hindmarsh, A. C.](#cite-Hindmarsh1983)
(1983). "'ODEPACK', a Systematized Collection of ODE Solvers". In:
_Scientific Computing, Vol. 1 of IMACS Transactions on Scientific
Computation_. Ed. by R. Stepleman. Amsterdam: IMACS / North-Holland,
pp. 55-64.

<a name=bib-Petzold1983></a>[Petzold, L. R.](#cite-Petzold1983) (1983).
"Automatic Selection of Methods for Solving Stiff and Nonstiff Systems
of Ordinary Differential Equations". In: _SIAM Journal on Scientific
and Statistical Computing_ 4, pp. 136-148.

<a name=bib-Soetaert2012></a>[Soetaert, K., J. Cash, and F.
Mazzia](#cite-Soetaert2012) (2012). _Solving Differential Equations in
R_. Springer Verlag.

<a name=bib-Soetaert2008></a>[Soetaert, K. and P. M. J.
Herman](#cite-Soetaert2008) (2009). _A Practical Guide to Ecological
Modelling. Using R as a Simulation Platform_. ISBN 978-1-4020-8623-6.
Springer, p. 372.

<a name=bib-ReacTran></a>[Soetaert, K. and F. Meysman](#cite-ReacTran)
(2010). _ReacTran: Reactive transport modelling in 1D, 2D and 3D_. R
package version 1.3.

<a name=bib-Soetaert2010_rjournal></a>[Soetaert, K., T. Petzoldt, and
R. W. Setzer](#cite-Soetaert2010_rjournal) (2010c). "Solving
Differential Equations in R". In: _The R Journal_ 2.2, pp. 5-15. DOI:
[10.32614/RJ-2010-013](https://doi.org/10.32614%2FRJ-2010-013).

<a name=bib-Soetaert2010b></a>[Soetaert, K., T. Petzoldt, and R. W.
Setzer](#cite-Soetaert2010b) (2010b). "Solving Differential Equations
in R: Package deSolve". In: _Journal of Statistical Software_ 33.1, pp.
1-25. DOI:
[10.18637/jss.v033.i09](https://doi.org/10.18637%2Fjss.v033.i09).

---

## Appendix

Things that should be mentioned somewhere, but not yet in the introductory
lecture

* Solver functions in deSolve
* C and Fortran interface

---

## Performance

Performance varies, depending on the type of the model. In general, simulation
can be fast if models can be written as matrices.

In case the performance appears as too slow, models can be written in **C** or 
**Fortran**. In such cases, control is transferred from R to compiled code and both
the solver and the model communicate directly. The interface is described in the [Compiled Code Package Vignette](https://cran.r-project.org/web/packages/deSolve/vignettes/compiledCode.pdf).

In addition, several third parpy packages exist that allow on-the-fly code 
generation from within R, especially 
package [cOde](https://CRAN.R-project.org/package=cOde), 
package [rodeo](https://github.com/dkneis/rodeo) and 
package [odin](https://mrc-ide.github.io/odin/).

---

## Solver Overview

![:scale 75%](img/solver-overview.png)