| type | name | unit | description | default |
|---|---|---|---|---|
| variable | LB | mg / mL | medium concentration | 1.5e+01 |
| variable | Ecoli | cells / mL | bacterial density | 1.0e+05 |
| parameter | mu_max | 1 / h | max. growth rate | 2.0e+00 |
| parameter | h | mg / mL | half saturation const. | 1.8e-01 |
| parameter | y | cells / mg | yield of Ecoli per LB | 1.0e+06 |
| parameter | q | mL / h | flow rate | 1.0e+01 |
| parameter | v | mL | reactor volume | 1.0e+02 |
| parameter | LB_in | mg / mL | input medium conc. | 1.5e+01 |
System analysis: Bioreactor model (using rodeo)
1 Study system
In this example, we consider a continuous, fully mixed bioreactor (CFSTR) where heterotrophic E. coli bacteria grow on a rich organic medium known as LB. The system is inoculated with bacteria once and thereafter fed with sterile medium only (Figure 1).
This particular reactor is not overly useful, but the principle has many applications, for example
in the industry, where microbial metabolites are exploited (e.g. in drug or food production),
in waste water treatment targeted at microbial degradation rather than production (with the inflow usually being non-sterile) .
2 Model specification
Using the matrix-based notation introduced earlier, the model is fully specified by the Table 1 and Table 2.
Table 1 does not need much explanation. It declares the two state variables and provides a set of default initial values. Very similar for the parameters. All other information is primarily meant as a documentation. A model without documentation is essentially useless!
Table 2 implements the ODE right hand sides. Essentially, the column with name rates holds the vector of process rates. The columns named like the two state variables make up the stoichiometry matrix. A minimum documentation is provided through columns description and unit. The name column is primarily needed to extract the value of the process rates from simulation output for, e.g., the purpose of plotting.
| name | unit | description | rate | Ecoli | LB |
|---|---|---|---|---|---|
| flow | 1 / h | im-/export | q / v | -Ecoli | LB_in - LB |
| growth | cells / mL / h | growth | mu_max * LB / (LB + h) * Ecoli | 1 | -1 / y |
3 Organizing the model with spreadsheets
Create a new workbook using spreadsheet software of your choice.
Add a sheet
declarationsto the workbook and fill in the contents of Table 1 (using copy-and-paste).Add a sheet
equationsto the workbook and fill in the contents of Table 2 (using copy-and-paste).Safe the workbook as
bioreactor.xlsxin a folder of your choice.
Alternative workflow: Save the contents of Table 1 and Table 2 as two separate text files named, e.g., declarations.txt and equations.txt, respectively. In the following code sections, replace any calls to buildFromWorkbook by buildFromText and provide the two text files as inputs.
4 Simulation of system dynamics
4.1 Model run with default inputs and outputs
Now, let’s run a dynamic simulation with that model. Unless we take special action, the defaults from sheet declarations (Table 1) will automatically be used as initial values and parameters.
rm(list=ls())
library("rodeo") # load package
workbook <- "models/bioreactor.xlsx" # adjust as necessary
m <- buildFromWorkbook(workbook) # assemble model
x <- m$dynamics(times=seq(0, 12, 0.1)) # run simulation
plot(x, log="y", las=2) # plot calls plot.deSolve here
Comparing the few lines of code above against previous code blocks hopefully convinces you of the benefit of the matrix-based notation and subsequent automatic code generation.
4.2 Customizing output graphics
The object returned by the dynamics method is a plain matrix holding columns for time, state variables, and process rates. The matrix is actually an object of class deSolve which controls the behavior of plot (Figure 2; see the help of plot.deSolve). For custom plotting, you may want to cast the matrix as a data.frame for flexibility. This is illustrated by the example below (Figure 3).
x <- as.data.frame(x) # matrix to data.frame
plot(Ecoli ~ time, data=x, type="l", log="y", yaxt="n", ylab="", col="red")
axis(2, col="red")
par(new=TRUE)
plot(LB ~ time, data=x, type="l", bty="n", yaxt="n", ylab="", col="blue")
axis(4, col="blue")
legend("right", bty="n", lty=1, col=c("red", "blue"),
legend=c("E. coli", "LB"))
In the context of plotting, it may be convenient to use the model’s namesVars and namesPros methods to distinguish state variables and process rates in the output matrix returned by dynamics.
4.3 Altering input values
To run a simulation with altered parameters or initial values, one could, of course, simple change the defaults stored in the workbook from which the model is loaded, save the workbook, and re-run the R code.
Alternatively (and faster), you can overwrite the defaults directly in the R code. For that purpose, you want to make use of the respective get?ars and set?ars methods as illustrated below (replace ? by P for parameters or V for variables, respectively).
p <- m$getPars() # query current parameters using *g*etPars
p["LB_in"] <- p["LB_in"] / 10 # overwrite a parameter (access by name)
m$setPars(p) # make change effective by calling *s*etPars
x <- m$dynamics(times=seq(0, 12, 0.1)) # re-run simulation
plot(x, log="y", las=2, which=c("LB", "Ecoli"))
4.4 Quick comparison of alternative dynamics
Very often, we want to compare predictions of a model for different initial values or parameters. For this purpose, the rodeo package provides as class method scenarios. This method can trigger a series of dynamic simulations where, in each run, the values of particular parameters and/or initial values are modified. All values not being set explicitly will remain unchanged. The scenarios method is demonstrated below. The method attempts to produce graphics for a quick visual comparison of the scenarios (Figure 5). These graphics are not necessarily ideal with regard to scaling, legend position, etc., and you typically want to use the returned data frame for making your own custom plots.
# scenarios must be defined as a named list
sc <- list(
default = c(m$getPars(), m$getVars()), # unchanged
clean_start = c(LB=0, Ecoli=1e3), # altered initial values
high_flow = c(q=10 * m$getPars()[["q"]]) # altered parameter
)
x <- m$scenarios(times=seq(0, 12, 0.1), scenarios=sc)
5 Finding steady-state solutions
5.1 Relevance
The control of biomass production in a continuous flow bio-reactor is best understood by looking at steady-state solutions. Then, by definition, biomass and resource concentration remain constant and so do the reactor’s output loads.
Imagine you are responsible for the operation of the bioreactor. Specifically, you are supposed to maximize steady-state production (i.e. the output load of bacteria) by choosing optimum values for the flow rate and the resource concentration in the input. How can you do that? Instead of going by trial-and-error, you certainly want support from modeling.
5.2 The quick way: Numerical long-term integration
Since we have the technology ready to predict the dynamics of bacteria in the bioreactor, we can exploit it for the calculation of steady-state conditions. All we need to do is to run the model with constant external forcings (flow, input concentration) over a very long time interval. While this strategy is not very elegant from a mathematical point of view, it is essentially fail-proof.
# perform long-term integration to find the steady state
times <- c(0, 24*365*20, 24*365*21) # 20 and 21 years !
x <- m$dynamics(times=times) # simulation
x <- x[,m$namesVars()] # only keep state variables
if (max(abs(x[2,]-x[3,])) > 1e-9) # check convergence
stop("not in steady state yet")
print(x[nrow(x),]) # final state = steady state LB Ecoli
1.800e-01 1.482e+07 5.3 The elegant way: Finding the analytic solution
Reverting the information from Table 1 into plain math notation yields the following set of ODE (Equation 1) where \(q/v\) was replaced by the inverse of the residence time, \(1/\tau\), for convenience.
\[ \begin{aligned} \tfrac{d}{dt}Ecoli &= & \mu_{max} \cdot \frac{LB}{LB + h} \cdot Ecoli & - \frac{1}{\tau} \cdot Ecoli \\ \tfrac{d}{dt}LB &= -\frac{1}{y} \cdot & \mu_{max} \cdot \frac{LB}{LB + h} \cdot Ecoli & + \frac{1}{\tau} \cdot (LB_{in} - LB) \end{aligned} \tag{1}\]
Fortunately, the steady-state solution for this system can be found using simple algebra. After setting the left-hand sides to zero (recall the definition of steady-state), the first ODE can be divided by \(Ecoli\) and then solved for the concentration of the medium (\(LB\)). The resulting expression is (Equation 2):
\[ LB = \frac{h}{\tau \cdot \mu_{max} - 1} \tag{2}\]
Similarly, the second ODE can be rearranged to find the steady-state concentration of \(Ecoli\) (Equation 3).
\[ Ecoli = \dfrac{y \cdot \left( LB_{in} - LB \right)}{\tau \cdot \mu_{max} \cdot \dfrac{LB}{LB + h}} \tag{3}\]
To solve Equation 3, you need to substitute all instances of \(LB\) on the right hand side by the previously found steady-state solution (Equation 2) and you will end up with the somewhat lengthy expression of Equation 4.
\[ Ecoli = \dfrac{y \cdot \left( LB_{in} - \dfrac{h}{\tau \cdot \mu_{max} - 1} \right)}{\tau \cdot \mu_{max} \cdot \left( \dfrac{\dfrac{1}{\tau \cdot \mu_{max} - 1}}{\dfrac{1}{\tau \cdot \mu_{max} - 1} + 1} \right)} \tag{4}\]
Fortunately, the entire puzzling denominator of Equation 4 turns out to vanish (i.e. equals unity) upon further simplification. Consequently, the steady-state solution for \(Ecoli\) is given by Equation 5.
\[ Ecoli = y \cdot \left( LB_{in} - \dfrac{h}{\tau \cdot \mu_{max} - 1} \right) \tag{5}\]
It is important to notice that the steady-state solution is only dependent on parameters, but it is independent of the system’s initial state. For example, the amount of bacteria used to inocculate the system does not matter.
Let’s validate Equation 5 against the solution found by long-term integration.
Ecoli_steady <- with(as.list(m$getPars()),
y * (LB_in - h / (v / q * mu_max - 1)) # analytic steady-state solution
)
print(Ecoli_steady)[1] 148200006 Sensitivity of the steady-state solution
6.1 Relevance
By means of sensitivity analysis, we can check how the bioreactor’s state and output responds to altered operating conditions. This is certainly valuable knowledge if we were responsible for operation under variable boundary conditions (e.g. temperatures, flow rate). Ultimately, sensitivity analysis is the first step to optimizing the reactor’s performance.
6.2 Sensitivity of biomass and output load with regard to flow rates
Here, we will analyze the response of the reactor’s steady-state to altered flow rates. Specifically, we focus on the bacterial concentration in the reactor and the corresponding output load.
For the purpose of sensitivity analysis, we could employ the numerical approach to finding the steady state. However, use of the analytical solution found earlier (Equation 5) is way more efficient from a computational point of view.
The following code section inspects both the concentration and output load of bacterial for a continuous, dense array of flow rates (Figure 6).
m <- buildFromWorkbook(workbook) # reset to defaults
q <- 1:250 # flow rates to investigate
Ecoli_steady <- function(q, pars) { # q: flow, pars: other parameters
with(as.list(pars[names(pars) != "q"]), # drop "q" from "pars", if present
y * (LB_in - h / (v / q * mu_max - 1)) # analytic steady-state solution
)
}
EC_conc <- sapply(q, Ecoli_steady, pars=m$getPars()) # concentrations
EC_load <- EC_conc * q # output loads
par(mfrow=c(1,2))
plot(q, EC_conc, type="l", xlab="q (mL/h)", ylab="E. coli (cells/mL)")
plot(q, EC_load, type="l", xlab="q (mL/h)", ylab="Output load (cells/h)")
par(mfrow=c(1,1))
The sensitivity analysis suggest that the concentration of bacteria remains almost constant over a wide range of flow rates (left panel of Figure 6). At the same time, the reactor’s output increases in an almost linear manner (right panel of Figure 6). However, when the flow rate exceeds a critical point, biomass production apparently breaks down. What is the cause of this crash?
Simple answer: As the flow rate increases, the input load of resources also increases which is reflected by the rise in bacterial biomass. However, a higher flow rate also means higher losses as bacteria are flushed from the reactor in shorter time. Beyond a certain flow rate, losses can no longer be compensated by growth as bacteria are already growing at their intrinsic maximum rate. At this point, the system switches from resource-limited to loss-controlled.
The critical flow rate where biomass production breaks down even if resources were provided in unlimited amounts (very high \(LB_{in}\)) can be identified from Equation 5 directly. Because of the relation \(\tau = v / q\), the denominator of the fraction is a function of the flow rate and, when \(\tau = 1 / \mu_{max}\), we get a zero division. This explains the suspicious discontinuity visible in Figure 6.
To summarize: Even if resources were available in huge excess, the actual bacterial growth rate cannot exceed \(\mu_{max}\). But the loss rate given by \(1/\tau\) (or \(q/v\)) does not have a fixed upper limit. So, whenever \(q/v > \mu_{max}\), a sustainable bacterial population cannot exist.
7 Finding optimum operating conditions
7.1 Objective
In real-world applications, we typically want bioreactors to reach the best possible performance (e.g., in terms production or degradation). Here, after seeing the outcome of the sensitivity analysis (Figure 6), we may be interested in finding the particular value of the flow rate which maximizes the reactor’s steady-state output load. This is the flow rate corresponding to the peak in the right hand side panel of Figure 6. Again, we can chose between a numeric or analytic solution strategy.
7.2 The elegant way: Finding the optimum analytically
With Equation 5 we have a closed-form expression for the steady-state concentration of bacteria but also for the reactor’s output load (concentration times flow rate). Now, what we need to find is a root of the derivative of the reactor’s output load with regard to the flow rate. Remember: Where the derivative is zero, the original function has a maximum (or minimum) according to this well-known theorem.
So, lets multiply Equation 5 by the flow rate \(q\) to obtain the expression for the output load of bacteria, \(L\), in steady state (Equation 6)
\[ L = q \cdot y \cdot \left( LB_{in} - \dfrac{h}{\tau \cdot \mu_{max} - 1} \right) \tag{6}\]
which can be simplified to Equation 7 considering the definition of the residence time (\(\tau = v/q\)).
\[ L = q \cdot y \cdot LB_{in} - \dfrac{q^2 \cdot y \cdot h}{\mu_{max} \cdot v - q} \tag{7}\]
Next, we need to derive Equation 7 with respect to the flow rate (\(dL/dq\)). Making use of the quotient rule of calculus, we get Equation 8
\[ \frac{dL}{dq} = y \cdot LB_{in} - \dfrac{2 \cdot q \cdot y \cdot h \cdot (\mu_{max} \cdot v - q) - q^2 \cdot y \cdot h \cdot (-1)}{(\mu_{max} \cdot v - q)^2} \tag{8}\]
and after a series of rearrangements (involving binomial formulas, particularly the equality of \(2ab - b^2\) and \(a^2 - (a-b)^2\)) we end up with the expression from Equation 9.
\[ \frac{dL}{dq} = y \cdot LB_{in} + y \cdot h - \dfrac{y \cdot h \cdot (\mu_{max} \cdot v)^2}{(\mu_{max} \cdot v - q)^2} \tag{9}\]
Since we want to find the root of the derivative, we set the left hand side of Equation 9 to zero and after another series of rearrangements, we obtain the reduced quadratic expression of Equation 10.
\[ 0 = q^2 - 2 \cdot \mu_{max} \cdot v \cdot q + \dfrac{(\mu_{max} \cdot v)^2}{1 + \frac{h}{LB_{in}}} \tag{10}\]
Equation 10 can be solved with the well-known formula for reduced quadratic expressions where the leading coefficient (the one in front of the quadratic term) equals one. Here, the valid solution of the two possible alternatives is the one with the negative sign before the square root (Equation 11).
\[ q = \mu_{max} \cdot v - \sqrt{(\mu_{max} \cdot v)^2 \cdot \left(1 - \dfrac{1}{1 + \frac{h}{LB_{in}}} \right)} \tag{11}\]
For the purpose of verification, we redraw Figure 6 with the critical flow rate predicted by Equation 11 marked by a dashed line (Figure 7).
plot(q, EC_load, type="l", ylim=c(0, 5e9), xlab="q (mL/h)", ylab="Output load (cells/h)")
abline(v=with(as.list(m$getPars()),
mu_max * v - sqrt((mu_max * v)^2 * (1 - 1 / (1 + h/LB_in) ))
), lty=2)
7.3 The quick way: Numerical solution with optimize
For more complex systems, deriving a closed-form expression for optimum operation conditions may not be possible. Moreover, calculus and fiddling about lengthy mathematical expressions is not everybody’s joy. Being able to use numerical optimization is thus a quite useful skill.
In the section on sensitivity analysis, we already developed a function to compute the steady state concentration of the bacteria for a given flow rate. We can instantly re-use that function to compute the corresponding output load.
Load_steady <- function(q, pars) {
Ecoli_steady(q, pars) * q # load = concentration * flow
}We can then ask the one-dimensional optimize function to carry out the search for the optimum flow rate within a prescribed interval. By default, optimize searches for a minimum of the provided function, so we want to set the maximum argument (or multiply our function by -1).
opt <- optimize(Load_steady, interval=c(0,250), pars=m$getPars(),
maximum=TRUE)
print(opt$maximum)[1] 178.2214Note that we could use this approach even in the absence of the closed-form steady-state solution (Equation 5). Then, we would simply make use of numerical computation twice: (1) for finding the steady state and (2) for optimization.
Load_steady_numeric <- function(q, model) {
parms <- model$getPars()
parms["q"] <- q # overwrite flow parameter
model$setPars(parms)
times <- c(0, 24*365*20, 24*365*21) # 20 and 21 years !
x <- model$dynamics(times=times) # simulation
x <- x[,model$namesVars()] # only keep state variables
if (max(abs(x[2,]-x[3,])) > 1e-9) # check convergence
stop("not in steady state yet")
x[nrow(x), "Ecoli"] * q # output load
}
opt <- optimize(Load_steady_numeric, interval=c(150,200), model=m,
maximum=TRUE)
print(opt$maximum)[1] 178.22148 The bioreactor from a state-space perspective
8.1 Basics
The state of a system at a particular point in time can be regarded as a point in the so-called state space (or phase space). The number of dimensions of the state space equals the number of state variables. In the bioreactor example, there would be just two dimensions: bacteria and substrate concentration.
If we would add, for example, another bacterial species and let both the reactor’s volume and temperature vary over time, we would end up with a 5-dimensional state space. A large number of dimensions is not a problem in general but it makes visualization more challenging.
8.2 Attractors
As long as external forcings (\(q\), \(LB_{in}\)) were kept constant, the bioreactor’s state variables always converged towards a steady state. Such points in state space toward which a system tends to evolve are called attractors.
Attractors can be imagined like magnets. If the the initial state of a system is close enough to the attractor, future states are pulled toward that point (or shape) in state space. The region around the attractor where the apparent magnetic force is noticeable is called a basin of attraction.
8.3 Attractor plot for the bioreactor
In a two dimensional state-space plot, the system’s state at a particular time point appears as a dot. When we connect the dots for consecutive time points, we obtain so-called trajectories. Plotting trajectories is a useful technique to visualize system dynamics. In a system with one or more steady states, trajectory plots can also be used to visualize attractors.
In the following code, we assume a set of contrasting initial states, each defined by a pair of values for the resource concentration (LB) and the bacterial inoculum (Ecoli). Starting from each of the hypothetical initial states, we predict the future system dynamics and instantly add the corresponding trajectory to a plot (Figure 8).
set.seed(2) # for reproducible results
m <- buildFromWorkbook(workbook) # reset model to defaults
plot(x=c(1e5, 2e9), y=c(1e-4, 50), type="n", # empty plot
log="xy", xlab="E. coli", ylab="LB")
initial <- as.matrix(expand.grid(LB=10^(-3:1), # assumed initial states
Ecoli=10^(5:9)))
for (i in 1:nrow(initial)) { # start from each ...
m$setVars(initial[i,]) # set initial values
x <- m$dynamics(times=c(0, 1.2^(-3:32))) # predict into future
lines(LB ~ Ecoli, data=x) # add trajectory
points(x[nrow(x),"Ecoli"], x[nrow(x),"LB"], # highlight final state
pch=20, cex=2)
}
What does this puzzling plot tell us?
Whenever we start with a very large inoculum, we see a sudden marked drop in resource concentration. Thereafter, bacterial density only declines slowly while the resource concentration rises again.
If both resource and bacterial concentration are initially low, the trajectories follow a curve. Only after resources have reached a higher level (through the incoming load), we experience massive proliferation of bacteria.
Finally, all trajectories lead to a unique point, the attractor, representing the steady state.
8.4 Classification of attractors (selection)
Fixed point attractors: A particular point in state space, as in Figure 8.
Limit cycles: A closed trajectory in state space. Other trajectories starting at close-by positions spiral into that trajectory as time passes. Limit cycles are usually connected to systems with stable oscillations.
Strange attractors: If the attractor cannot be described by a simple geometric shape (e.g. a point, line, surface, or sphere), it is called “strange”. Strange attractors are usually related to systems whose dynamics are hard to predict, because the solution is extremely sensitive to perturbations. Strange attractors are linked to deterministic chaos. A well-known example is the so-called Lorenz attractor.
8.5 Relevance of studying trajectories and attractors
Knowledge about attractors and its counterpart (so-called repellers) is essential for understanding the vulnerability of systems to perturbations. The simple bioreactor is an example of a very stable system with a single attractor (steady state) to which the system returns under all circumstances (Figure 8).
However, more complex systems often have multiple steady states. It is then highly relevant to know under which conditions we must expect a switch from one steady state to the other. Such critical conditions are often referred to as tipping points.