System analysis: Simultaneous ODE
\(\leftarrow\) BACK TO TEACHING SECTION OF MY HOMEPAGE
1 Motivation
The description of environmental systems generally requires multiple state variables, including, e.g., water levels, concentrations, species abundances. The state of a system is thus expressed as a vector \(\{ y_1, y_2, \dots y_n \}\).
Since all state variables potentially undergo dynamics, their derivatives with respect to time form a corresponding vector \(\{ \frac{d}{dt}y_1, \frac{d}{dt}y_2, \dots \frac{d}{dt}y_n \}\) often shortly written as \(\{ \dot{y_1}, \dot{y_2}, \dots \dot{y_n} \}\).
Very typically, the derivative of one variable is dependent on the value of another variable, resulting in relations like, e.g., \(\frac{d}{dt}y_1 = f(y_2)\). The following examples illustrate such a situation:
The growth of algae results in a change of population density. Since growth requires the uptake of dissolved phosphorus, the change in population density depends on the ambient P concentration.
In a predator-prey system, the growth of the predator population is dependent on prey density. At the same time, the decline in prey population is controlled by the abundance of predators.
If we like to understand the dynamics of those systems, we need to deal with simultaneous ODE. The term simultaneous is used to express coupling between equations. It is impossible to integrate just one of the ODE without simultaneously integrating the other ones too, carefully keeping track of the dependencies.
The coupled nature of the equations must also be recognized in the calculation of steady-states, because a solution must satisfy \(\{ \frac{d}{dt}y_1 = 0, \frac{d}{dt}y_2 = 0, \dots \frac{d}{dt}y_n = 0 \}\). The resulting equation systems are not necessarily tractable by analytical methods because they are often non-linear.
Consequently, systems of simultaneous ODE are generally tackled by numerical methods. Numerical methods generally start with an initial guess of the solution which is refined iteratively until a desired level of accuracy is reached. Numerical integration is not an exception.
To get a basic idea of how numerical integration works, consider Figure 1. The latter illustrates the idea behind Euler’s method applied to a single ODE. The method is equivalent to plain linear extrapolation. Actual numerical algorithms are still rooted in this basic concept but use sophisticated iteration schemes to achieve a very high accuracy.
2 R packages for solving simultaneous ODE
A well-tested and widely used add-on package for numerical integration in R is deSolve. You can find the full package documentation on this CRAN page. The package is especially useful if one wants to simulated system dynamics.
However, the integration methods also help to compute steady-state solutions. For that purpose, the integration is simply performed over a very long period of time with constant boundary conditions. Although this is not very elegant compared to a direct solution, it is a very robust approach that is guaranteed to work. By comparison, the direct solutions available through the rootSolve package (see this CRAN page) are way more fragile.
3 Basics: Example of two connected reservoirs
3.1 Experimental setup and governing equations
To start with, we consider a system of two reservoirs connected by a tube (Figure 2). This system has two state variables: the volume of the bigger tube on the right (\(V_b\)) and the volume of the smaller tube on the left (\(V_s\)). In accordance with experience, the water surfaces in the two reservoirs will level out over time as water is exchanged through the tube (if the latter is filled initially).
If we denote the flow rate in the tube as \(Q\) and assign a positive sign to flow from the big to the small reservoir, we get the following simultaneous ODE (Equation 1).
\[ \begin{aligned} \tfrac{d}{dt}V_s &= Q \\ \tfrac{d}{dt}V_b &= -Q \end{aligned} \tag{1}\]
What we need now is a physics-based expression for \(Q\). Let’s adopt the essence of Darcy’s law which suggests a linear relation between the apparent velocity (\(u\)) and the difference in hydraulic head (\(h\)) over distance \(L\) (Equation 2).
\[ u = k \cdot \frac{\Delta h}{L} \tag{2}\]
Equation 2 originally refers to flow in porous media. To obtain an expression for the flow rate \(Q\), we can combine Equation 2 with the continuity equation (Equation 3) where \(a\) denotes the cross-sectional area through which the liquid is flowing.
\[ Q = u \cdot a \tag{3}\]
For the case depicted by Figure 2, it seems reasonable to collapse all unknown constants into a single parameter denoted \(f\) (Equation 4). We also replace the \(\Delta h\) by the difference in water levels (\(H\)) according to our definition of positive flow direction.
\[ \begin{aligned} Q &= k \cdot \frac{\Delta h}{L} \cdot a \\ &= f \cdot \Delta h \\ &= f \cdot (H_b - H_s) \end{aligned} \tag{4}\]
Using Equation 4, we can now extend the right hand sides in Equation 1. Doing so, we also want to translate water levels (\(H\)) into the actual state variables (\(V\)) based on the individual cylinders’ cross-sectional area \(A\) (Equation 5).
\[ H = \frac{V}{A} \tag{5}\]
Compiling all pieces together, we get the system of simultaneous ODE depicted by Equation 6. It involves the volumes \(V\) and cross-sectional areas \(A\) of the big and small reservoir (subscripts “b” and “s”, respectively) and the resistance parameter \(f\). Just like in Equation 1, the right hand sides solely differ by the sign only.
\[ \begin{aligned} \tfrac{d}{dt}V_s &= f * \left( \frac{V_b}{A_b} - \frac{Vs}{A_s} \right) \\ \tfrac{d}{dt}V_b &= -f * \left( \frac{V_b}{A_b} - \frac{Vs}{A_s} \right) \end{aligned} \tag{6}\]
3.2 Implementing the equations for use with deSolve
We now want to integrate the simultaneous ODE from Equation 6 with numerical methods. In order to use the numerical solvers from deSolve for tackling initial value problems, we have to provide four elements of code:
A vector with the initial values of all state variables with elements being named. Considering Figure 2, an appropriate statement would be
init <- c(Vs=28, Vb=70)instructing the integration algorithm to start from those values.A vector holding the values of any parameters appearing at the right hand side of the ODE system. An appropriate vector could look like
pars <- c(f=0.1, As=3.8, Ab=6.1)based on Equation 6.A vector of times in the future at which the values of the state variables are to be predicted. It could look like, e.g.,
times <- 0:200ortimes <- seq(0, 180, 30)Finally, we must provide a function to compute the right hand sides of the simultaneous ODE. As can be seen in Equation 6, the values depend on the current values of the state variables, the parameters, and (possibly but not here) on the time. Hence, this function’s interface must provide formal arguments for theses three inputs. Consequently, the interface must look like
function (time, vars, pars). When the function is called to actually compute derivatives, the time is passed astime, the state variables’ values asvars, and the parameters aspars.
Next, we need to understand what goes into the body of the function with the characteristic interface function (time, vars, pars). Here is a brief summary:
The function must compute the vector of derivatives. Let’s call it
dydt. Under all circumstances, the elements of elementsdydtmust be in the same order as the elements in the input vectorvars. A safe way to achieve this is a statement likedydt <- dydt[names(vars)].The vector
dydtmust be wrapped into a list before it is being returned from the function. This is technically necessary to separate the computed derivatives from any (optional) auxiliary information to be returned (e.g. for the purpose of debugging).
Consequently, the skeleton of a function to return the state variables’ derivatives would look like this:
derivatives <- function(time, vars, pars) { # time, state v., parameters
# compute the vector of derivatives using 'vars', 'pars', possibly 'time'
dydt <- c(
# statement for 1st variable
# statement for 2nd variable
# ...
)
# ensure that derivatives are ordered like the values in 'vars'
dydt <- dydt[names(vars)]
# return derivatives wrapped into a list
list(dydt)
}When you compose the contents of the vector dydt, a state variable named, e.g., “Vs” should generally be accessed as vars[["Vs"]]. Similarly, a parameter “As” should be accessed as pars[["As"]]. The use of double brackets has two advantages. First, R will complain immediately about attempts to access any undefined item with a clear message. Second, the [[ strips the element name to only yield the bare numerical value. This avoids puzzling name conflicts later on…
If the right hand sides are rather complex, the with statement of R can simplify the access to the values in the vectors vars and pars. If, for example, the vector vars has a named element “Vs”, the following code would make it accessible just as Vs (instead of vars[["Vs"]]). However, while the use of with is convenient, it is much easier to write buggy code compared to using the [[ ]] operator.
dydt <- with( as.list(c(vars, pars)), { # "dissolves" vectors into scalars
c(
# statement for 1st variable
# statement for 2nd variable
# ...
)
}) # end of the with `block3.3 Application to experimental data
The code section below loads a data set obtained with the equipment from Figure 2 and plots the observations (Figure 3).
rm(list=ls())
library("deSolve")
dynam <- read.table(header=TRUE, sep="", text='
# volume in the big cylinder (Vb) and the one small (Vs)
# units: mL
seconds Vb Vs
22 93 4.5
31 88 9.5
35 85 12.5
43 82 16
47 80 19
54 76.5 21
59 75 22.5
65 72.5 24.5
82 68 29
102 65 32.5
129 62.5 35
158 61 36.5
')
# cross-sectional area of the big cylinder (Ab) and the
# small one (As), units: cm squared
const <- c(Ab=6.107, As=3.791)
# let observations start at time zero
dynam[,"seconds"] <- dynam[,"seconds"] - dynam[1,"seconds"]
# plot observations
par(mfcol=c(1,2))
with(dynam, {
plot(seconds, Vs)
plot(seconds, Vb)
})
par(mfcol=c(1,1))
Now, we implement the equations from Equation 6 in a deSolve-compliant manner. After collecting initial values and parameters in the respective vectors, we call the lsoda method to perform the numerical integration. Behind lsoda is a sophisticated algorithm that typically masters even difficult integration problems, namely stiff equations.
Note that the value of the friction parameter \(f\) was found by trial-and-error. If you are keen on learning how to find the optimum value by fitting, look at R’s optimize method for one-dimensional problems like this.
# implementation of the system of ODE
derivs <- function (t, y, p) { # used shorter names
flow <- p[["f"]] * # compute just once
(y[["Vb"]]/p[["Ab"]] - y[["Vs"]]/p[["As"]])
dydt <- c(Vs= flow, Vb= -flow) # set derivatives
dydt <- dydt[names(y)] # ensure element order
list(dydt) # wrap result into list
}
# initial values taken from first observation
init <- unlist(dynam[1,c("Vb", "Vs")])
# parameters; value of 'f' found by trial-and-error
pars <- c(const, f=.055)
# perform integration
simul <- deSolve::lsoda(y=init, times=0:180, func=derivs, parms=pars)
# plot model result together with observation
par(mfcol=c(1,2))
for (v in c("Vs","Vb")) {
plot(dynam[,"seconds"], dynam[,v], xlab="seconds", ylab=v)
lines(simul[,"time"], simul[,v])
if (v == "Vb") {
legend("topright", bty="n", pch=c(1,NA), lty=c(NA,1),
legend=c("Observations","Model"))
}
}
par(mfcol=c(1,1))
4 Generic concept for simultaneous ODE models
4.1 Objectives
In the connected reservoirs example (Section 3), we dealt with just two state variables whose dynamics was due to a single process. Real-world ecological or technical systems, however, consist of many interacting state variables affected by a large number of processes.
Models accounting for both multiple state variables and multiple processes need more care for two reasons:
In spite of more variables and equations, the corresponding program code must remain transparent, easy to maintain, and reusable.
Computational efficiency matters, especially for long-term simulations, parameter estimation, or uncertainty analysis.
Fortunately, transparency and efficiency can achieved at the same time by adopting a standard notation for simultaneous ODE. The latter is used in chemical engineering for decades already. It became more popular in the biological sciences through activated sludge modeling (e.g., Henze et al., 1999), river water quality modeling (e.g., Reichert et al., 2001) or ecological lake modeling (e.g., Omlin et al., 2001).
In the following subsections, I will introduce this notation with an example.
4.2 Example: Aerobic degradation of organic matter
We consider the aerobic decay of organic matter (\(OM\)) according to a first order kinetics and the corresponding impact on dissolved oxygen (\(DO\)). The pool of \(DO\) is simultaneously depleted by microbial consumption and replenished by atmospheric aeration (Figure 5).
For simplicity, we consider a system without external input or outputs, e.g. a bottle. The perspective would also be applicable to a travelling control volume in an idealized river with negligible longitudinal mixing (see the plug-flow reactor concept).
The dynamics of the two state variables \(OM\) and \(DO\) is expressed by the simultaneous ODE from Equation 7
\[ \begin{aligned} \tfrac{d}{dt}OM &= -k_d \cdot OM \\ \tfrac{d}{dt}DO &= -k_d \cdot OM \cdot s + k_a \cdot (DO_{sat} - DO) \end{aligned} \tag{7}\]
where
\(k_d\) is a rate constant to describe the efficiency of degradation (unit: 1/time). The value of that constant must account for both the density of the bacterial population and the activity of the cells.
\(k_a\) is the rate constant of re-aeration (unit: 1/time). The value certainly depends on turbulence.
\(DO_{sat}\) represents the oxygen saturation level in equilibrium with the atmosphere. The value is mostly affected by temperature and salinity.
\(s\) is a stoichiometric coefficient specifying the mass (or moles) of oxygen consumed per mass (or moles) of \(OM\). It can roughly be estimated from the stoichiometry of the underlying chemical reaction.
Like the classical model of Streeter & Phelps published in 1925, Equation 7 is obviously overly simplistic in several aspects:
The rate constant of degradation, \(k_d\), does not vary over time. Thus, in the model, microbes do not reproduce. If reproduction was to be considered, one would need to add bacteria as another state variable!
The equations do not account for oxygen limitation or a possible adaptation of the microbial community to anaerobic conditions. Hence, the model returns non-sense if \(DO\) concentrations fall below about 2 mg/L.
Further, the model does not disentangle oxygen consumption due to the oxidation of carbon (CBOD) or ammonium (NBOD), is does not consider settling and the role of surface-attached biofilms, a possible accumulation of hardly degradable fractions of \(OM\) is neglected, and a constant temperature is assumed.
4.3 Matrix-based notation for simultaneus ODE
Nevertheless, in spite of these “biological shortcomings”, Equation 7 is well suited to demonstrate the standard notation for simultaneous ODE. The fundamental idea is to rewrite Equation 7 as a vector-matrix multiplication like so (Equation 8):
\[ \begin{bmatrix} \tfrac{d}{dt}OM, & \tfrac{d}{dt}DO \end{bmatrix} = \begin{bmatrix} k_d \cdot OM, & k_a \cdot (DO_{sat} - DO) \end{bmatrix} \cdot \begin{bmatrix} -1, & -s \\ 0, & 1 \end{bmatrix} \tag{8}\]
As usual, the state variables’ derivatives appear on the left hand side, but now as a vector. The row vector forming the first factor on the right hand side typically represents the individual processes. The second factor is known as the stoichiometry matrix. The latter connects the processes with the state variables derivatives. Each column of the matrix corresponds to a state variable, each row is associated with a process.
The signs of the stoichiometric coefficients allow for a simple interpretation: If a coefficient at position [\(i,k\)] is positive, the \(i\)-th process triggers an increase in the value of the \(k\)-th state variable. Likewise, negative coefficients indicate a decline of the variable’s value in response to the action of the process. Finally, a coefficient of zero indicates that the process does not have a direct impact on a variable (which does not exclude indirect effects).
So why should we prefer the matrix-based notation (Equation 8) over Equation 7? There are multiple reasons:
the matrix-based notation eliminates redundant terms. In the above example, the term \(-k_d \cdot OM\) appears twice in Equation 7 but only once in Equation 8. Thus, the equations are simpler to write and there is a gain in computation speed.
the stoichiometry matrix helps to validate the model because it summarizes all interactions between processes and variables. It also helps to communicate the model concept (see, e.g., Fig. 2 of this paper).
the standard notation allows the model equations to be stored in tables (e.g. spreadsheets). This simplifies model documentation and allows for the automatic translation into source code. Hence, one can implement models of simultaneous ODE without actual programming.
4.4 Use of the matrix notation with deSolve
The code below illustrates how the matrix-based notation can be used together with the numerical solvers from the deSolve package. Of particular importance is a small detail: the %*% operator. It is the key to matrix multiplication in R. The code’s output is displayed in Figure 6.
rm(list=ls())
# Implementation of the aerobic degradation model
dydt <- function(time, vars, pars) {
with(as.list(c(vars, pars)), { # 'dissolve' input vectors
# Vector of process rates
proc <- c(
degradation = kd * OM,
reaeration = ka * (DOsat - DO)
)
# Stoichiometry matrix
stoi <- rbind(
degradation = c(OM= -1, DO= -s)[names(vars)],
reaeration = c(OM= 0, DO= 1)[names(vars)]
)
stopifnot(identical(names(proc), rownames(stoi)))
# note: The '[names(y)]' in above statements enforces the
# proper order of columns. We also check for a possible
# mismatch in row order.
# Return derivatives, process rates as extra output
# (!) Note the %*% operator for matrix multiplication
# and the t() function performing transposition
list(
derivs= t(proc) %*% stoi, # derivatives
rates=proc # process rates
)
}) # end of 'with' block
}
# Define parameters, initial values, and times of interest
pars <- c(kd=1, ka=0.5, s=2.67, DOsat=10)
vars <- c(OM=1, DO=10)
times <- seq(from=0, to=10, by=1/24)
# Solve and visualize
dyn <- lsoda(y=vars, parms=pars, func=dydt, times=times)
plot(dyn)
In the R code above code, we have used a value of s=2.67 for the amount of oxygen consumed per amount of degraded organic matter. Would you be able to fill the missing values in the table below?
| Unit of OM | Unit of DO | Value of parameter s |
|---|---|---|
| g C mL-1 | g O2 mL-1 | 2.67 |
| g C6H12O6 mL-1 | g O2 mL-1 | ? |
| moles C mL-1 | moles O2 mL-1 | ? |
| moles C6H12O6 mL-1 | moles O2 mL-1 | ? |
4.5 Summary and recommendation
The above example should have made clear that the matrix-based standard notation is simple but helpful. I recommend to always use it because it has proven superiority over loose collections of ODE in almost all applications. The notation is summarized again in generic form in Equation 9 with processes denoted \(a\), \(b\), … \(z\) and state variables dinstinguished by numeric indices. Note the transpose operators (\(T\)) which allow for a more compact notation using column vectors instead of row vectors.
\[ \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}, & s_{a,2}, & \ldots, & s_{a,n} \\ s_{b,1}, & s_{b,2}, & \ldots, & s_{b,n} \\ \ldots, & \ldots, & \ldots, & \ldots \\ s_{z,1}, & s_{z,2}, & \ldots, & s_{z,n} \\ \end{bmatrix} \tag{9}\]