System analysis: Simultaneous ODE

Author

David Kneis (firstname.lastname @ tu-dresden.de)

Published

November 29, 2024

\(\leftarrow\) BACK TO TEACHING SECTION OF MY HOMEPAGE

1 Introduction

1.1 The meaning of “simultaneous”

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 ambient P, but it also alters the 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 solve just one of the ODE without simultaneously solving the other ones too, carefully keeping track of the dependencies.

1.2 Example: 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 1).

Figure 1: A bacterial culture (represented by a single cell) taking up organic matter and oxygen from the medium. The culture vessel allows for the replenishment of consumed oxygen by atmospheric exchange. Organic matter is represented by a glucose molecule.

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 1

\[ \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{1}\]

where

\(k_d\) is a rate constant to describe the efficiency of degradation (unit: 1/time). The value of that constant must account for the (constant) density of the bacterial population and the respiration 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 depends on temperature and salinity.
\(s\) is a stoichiometric coefficient specifying the amount of oxygen consumed per amount of degraded \(OM\) (units discussed below in Table 1).

Like the classical model of Streeter & Phelps published in 1925, Equation 1 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. In fact, \(DO\) can even go negative.
Further, the model does not acount for the settling of organic particles, the possible role of surface-attached bacterial biofilms, a possible accumulation of recalcitrant \(OM\), and many other details …

Let’s spent some thoughts on the value of the stoichiometric parameter \(s\) appearing in Equation 1. As it specifies how much oxygen is consumed per amount of degraded substrate, we should be able to infer a reasonable value from the corresponding chemical reaction. Let us adopt the reaction formula

\[CH_2O + O_2 \rightarrow CO_2 + H_2O \]

where the organic matter has the same stoichiometry as glucose (\(C_6H_{12}O_6\)) but, for simplicity, all quantities have been divided by six. Apparently, one mole of \(O_2\) is consumed per mole of organic carbon (\(C\)) which results in \(s=1\).

However, unlike marine scientists, freshwater biologists tend to use mass concentrations instead of molar concentrations which necessarily affects the value of \(s\). In addition, the concentration of organic matter is not necessarily specificed in units of carbon but it could also be, for example, glucose equivalents. A few typical cases are illustrated in Table 1.

Table 1: Appropriate values of parameter \(s\) depending on chosen units.
Unit of OM	Unit of DO	Value of parameter s
mol C mL^-1	mol O₂ mL^-1	1
g C mL^-1	g O₂ mL^-1	32/12 = 2.67
mol C₆H₁₂O₆ mL^-1	mol O₂ mL^-1	6
g C₆H₁₂O₆ mL^-1	g O₂ mL^-1	192/180 = 1.067

To find the values of \(s\) listed in Table 1, you need to know

the stoichiometry of the reaction (given by the reaction formula)
the molar masses of the involved elements as specified in the periodic table
the rule for converting between molar and mass concentrations typically written as \(n = m / M\) with \(n\): amount of moles, \(m\): mass in grams, \(M\): molar mass of the substance in grams per mole.

2 Approaches to solving simultaneous ODE

2.1 General features

Besides a few special cases, systems of simultaneous ODE cannot be solved by analytical methods. Consequently, we need to adopt numerical algorithms to find approximate solutions in an iterative manner. Obviously, at this point, computers become really essential.

Like in the case of individual ODE, there are two possible goals when solving simultaneous ODE: We can be interested in either the system’s dynamics or a possible steady-state.

2.2 Dynamic solutions

2.2.1 Theory

Dynamic solutions of the initial value problem require integration. To get a basic idea of how numerical integration works, consider Figure 2. The latter illustrates the idea behind Euler’s method applied to a single ODE. The method is equivalent to linear extrapolation into the future using the known derivative(s). Actual numerical algorithms are still rooted in this basic idea but they employ sophisticated iteration schemes to achieve convergence towards the true solution. The Runge-Kutta methods are good examples of such algorithms with moderate complexity.

Figure 2: The simplest approach to numerical integration known as Euler’s method. The idea is to start from a known initial value and estimate the function’s value at a future point in time using the known derivative. Essentially, the tangent is used for extrapolation. In the case of non-linear functions, the estimate is necessarily very inaccurate.

2.2.2 Computing dynamic solutions in R

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 and the corresponding package vignettes. We will use this package heavily in the future.

My favorite solver from the deSolve package is called lsoda as it is very efficient and tackles even difficult problems where simpler solvers may struggle. You may want to inspect the respective help page in detail later.

install.packages("deSolve")        # install package (once)

library("deSolve")                 # load package

?lsoda                             # help for recommended solver

We will apply this solver to Equation 1 shortly in Section 3.

2.3 Steady-state solutions

2.3.1 Theory

A steady-state solution of a system of simultaneous ODE can be obtained by setting all left-hand sides to zero, indicating no change in any of the variables (Equation 2).

\[ \begin{aligned} \frac{d}{dt}y_1 &= 0 \\ \frac{d}{dt}y_2 &= 0 \\ \dots \\ \frac{d}{dt}y_n &= 0 \end{aligned} \tag{2}\]

The resulting algebraic system of equations is often non-linear and can only be tackled with numerical root-finding algorithms. An initial guess of the solution is generally needed to start such algorithms.

Remember that a non-trivial steady-state solution does not necessarily exist. For example, Equation 1 apparently has only a trivial steady state at \(OM = 0\) and \(DO = DO_{sat}\).

2.3.2 Computing steady-state solutions using R

To directly obtain steady-steady solutions, one can use the rootSolve package available through this CRAN page. It is in fact a companion-package to deSolve. Although the direct approach is elegant and often fast, the algorithm may be “fragile”. Specifically, the root-solving algorithm often requires a good initial guess to converge towards a solution. A good guess may not be easy to obtain.

As an alternative, one can simply use the same approach as for initial value problems (i.e., deSolve) and run the integration over a very long period of time while keeping all boundary conditions at constant values. For this, one only needs a guessed initial state but not a guess of the solution. In fact, a very crude approximation of the initial state is sufficient (unless the system has multiple steady states). The only drawback of this very robust approach may be the computational overhead of a very-long-time integration.

3 Dynamic solution of the degradation example

Let’s integrate the simultaneous ODE from Equation 1 using the deSolve package. To do so, we have to provide four elements of code:

A vector with the initial values of all state variables with elements being named. A suitable line of code for the example might read like:

initial <- c(OM=1, DO=8/32)

A vector holding the values of any parameters appearing at the right hand side of the ODE system. A suitable line of code might be:

parameters <- c(kd=0.07, ka=0.1, DO_sat=8/32, s=1)

A vector of times in the future at which the values of the state variables are to be predicted. Of course, time units should be consistent with the units of the parameters.

times <- 0:72

Finally, we must provide a function to compute the right hand sides of the simultaneous ODE. As can be seen in Equation 1, 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 as time, the state variables’ values as vars, and the parameters as pars.

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 elements dydt must be in the same order as the elements in the input vector vars. A safe way to achieve this is a statement like dydt <- dydt[names(vars)].
The vector dydt must 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.

Consequently, a function to return the state variables’ derivatives would look like this:

derivatives <- function(time, vars, pars) {

  # compute the vector of derivatives using 'vars', 'pars' ('time' not used)
  dydt <- c(
    OM = unname(-pars["kd"] * vars["OM"]),
    DO = unname(-pars["kd"] * vars["OM"] * pars["s"] +
      pars["ka"] * (pars["DO_sat"] - vars["DO"]))
  )
  
  dydt <- dydt[names(vars)]    # sort derivatives like the values in 'vars'
  list(dydt)                   # return derivatives wrapped into a list
}

Note the use of unname in the statements above. This function avoids an undesired carry-over of element names from the right hand side (one of the few annoying “features” of R).

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 “OM”, the following code would make it accessible just as OM (instead of vars["OM"]). Although the use of with is very convenient, it is much easier to make bad mistakes.

# function as before, coded differently
derivatives <- function(time, vars, pars) {

  # compute the vector of derivatives using 'vars', 'pars' ('time' not used)
  dydt <- with( as.list(c(vars, pars)), {  # "dissolves" vectors into scalars
    c(
      OM = kd * OM,                               # much better to read ...
      DO = kd * OM * s + ka * (DO_sat - DO)       # ... but more error-prone 
    )
  })
  
  dydt <- dydt[names(vars)]    # sort derivatives like the values in 'vars'
  list(dydt)                   # return derivatives wrapped into a list
}

Now that all required elements of code are ready, let’s call the solver and inspect the predicted dynamics visually (Figure 3).

# call the solver providing it with the initial state, the target times,
# our function to compute the derivatives, and the parameters
x <- lsoda(y=initial, times=times, func=derivatives, parms=parameters)

# outputs from the solver can be plotted instantly
plot(x)

Figure 3: Predicted dynamics of the state variables from Equation 1 for given initial values and parameters.

The predicted dynamics of \(OM\) in Figure 3 is in good agreement with expectation. A first order decay processes should result in convergence towards zero as we have seen earlier.

The predicted dynamics of \(DO\) is a bit more tricky. Apparently, larger amounts of oxygen are consumed initially as the degradation process is very active. Later, as most of \(OM\) has already been degraded, the rate of oxygen consumption declines and concentrations recover due to replenishment from the atmosphere. However, the fact that \(DO\) temporarily goes negative is certainly problematic and reflects a conceptual issue with Equation 1. In fact, the problem is a missing limitation term that prevents the degradation process from continuing once \(DO\) is fully depleted. We will look at approaches to modeling resource limitation in a subsequent session.