print(system.file("models", package="rodeoEasy"))System analysis: Limitation and inhibition of growth
\(\leftarrow\) BACK TO TEACHING SECTION OF MY HOMEPAGE
1 Motivation
The budget concept of system theory also applies to populations. While reproduction and immigration add individuals to a population, individuals are lost due mortality and emigration across system boundaries (Figure 1).
In this section, we will take a closer look at the controls of reproduction with a somewhat narrow focus on microorganisms.
What is special about reproduction in microorganisms?
Single-cell organisms like bacteria or eukaryotic phytoplankton proliferate through cell division. Hence, we do not need to distinguish between the growth of individuals and population growth.
Microbial populations generally comprise a very large number of individuals capable of reproducing continuously (e.g. independent of season) under suitable conditions.
Because of these two facts, the growth of microbial populations can be regarded as a continuous process. Consequently, the dynamics of the population size (\(N\)) can be mathematically described by a simple linear ODE (Equation 1) with the growth rate constant \(\mu\) (unit: 1/time).
\[ \tfrac{d}{dt}N = \mu \cdot N \tag{1}\]
As we have seen earlier, integration of Equation 1 for a known initial value leads to an exponential expression (Equation 2). We also learned that Equation 2 implies a simple unique relationship between \(\mu\) and the doubling time \(T\) (Equation 3).
\[ N(t) = N(t_0) \cdot e^{\mu \cdot (t - t_0)} \tag{2}\]
\[ \mu = \frac{ln(2)}{T} \tag{3}\]
Because of the exponential characteristics, the density of microbial populations can increase very fast. The larger the population is already, the more rapid it will grow further. It is obvious that exponential growth cannot continue over long periods of time since population sizes would go out of bounds. Hence, there must be mechanisms that restrict the growth of populations. We will address these mechanisms and the corresponding mathematical formulations in the following sections.
2 Technical note
In the following, we will consider several models implemented in spreadsheet workbooks. They are meant to be processed through the code generation and simulation functions of the R package rodeoEasy. We installed this package earlier.
For simplicity, I distribute the workbooks as part of the rodeoEasy package. To open and check the contents of the workbooks, you need to find them on your system. For that purpose, you should use the following line of R code:
The command will print the location of a folder ‘models’ on your local computer. You can use your file manager to navigate to this folder and check / open the workbooks there.
Warning: If you want to alter these workbooks shipped with the rodeoEasy package, you should make a copy of the original files and save them outside this system folder, e.g. where you normally keep your documents. Reason: The ‘models’ folder including all workbooks will be overwritten (without warning) if you re-install the rodeoEasy package.
3 Growth with limited resources
3.1 Carrying capacity concept
The simplest way to account for resource limitation is to assume a maximum population density, the so-called carrying capacity. As the population size approaches this upper limit, the actual growth rate converges to zero (Equation 4). For suspended microorganisms, a typical unit of the carrying capacity, denoted \(N_{max}\), is cells mL-1.
\[ \tfrac{d}{dt}N = \mu \cdot N \cdot \left( 1 - \frac{N}{N_{max}} \right) \tag{4}\]
Equation 4 allows to account for resource limitation without introducing the resource itself as a second state variable. This is especially useful if the exact mechanism of limitation is not understood in detail.
The formulation of Equation 4 is also known as the logistic growth model (not to be confused with a logistic regression model). The two steady state solutions are \(N=0\) and \(N=N_{max}\). A dynamic solution of Equation 4 can still be found analytically, but the integration requires some knowledge about partial fraction decomposition. I will rather illustrate the numerical solution instead (R code below and Figure 2).
rm(list=ls())
library("deSolve")
dydt <- function(time, vars, pars) {
with(as.list(c(vars, pars)),
list(mu * N * (1 - N/Nmax)) # growth model
)
}
init <- c(N=1e4) # initial value
pars <- c(mu=2, Nmax=1e9) # parameters
dynam <- lsoda(y=init, parms=pars, # integration
func=dydt, times=seq(0, 18, .1))
par(mfrow=c(1,2)) # visualization
plot(dynam[,"time"], dynam[,"N"], type="l",
xlab="Time", ylab="Cells / mL (linear)")
plot(dynam[,"time"], dynam[,"N"], type="l", log="y",
xlab="Time", ylab="Cells / mL (log scale)")
par(mfrow=c(1,1))
The left hand plot of Figure 2 illustrates the typical, S-shaped curve predicted by the logistic growth model. Initially, the population size increases almost negligibly but later, the curve turns sharply upwards indicating the exponential phase. The exponential phase is followed by a transition and the stationary phase where the number of individuals stabilizes as \(N\) approaches \(n_{max}\) (i.e. resources are exhausted).
In the microbial context, the population dynamics is typically plotted on a logarithmic scale (Figure 2, right) to cover several orders or magnitude. In a log-scale plot, the exponential phase is represented by the linear part of the curve. In fact, the exponential phase can only be properly identified in the log-scale plot.
3.2 Monod-type models
3.2.1 Classical Monod model
Although the carrying capacity model from Equation 4 can depict the phenomenon of resource limitation, it does not actually describe the mechanism. A mechanistic understanding is required however, if we like to control population sizes in a bioreactor or an ecosystem. For that purpose, we need to extend the growth equation with an explicit dependency on the resource. At the same time, we need to model the depletion of the resource in response to population growth.
This takes us to a system of two simultaneous ODE (Equation 5) accounting for both the change in population density (\(N\)) and resource concentration (\(R\)). Note that this system describes a batch experiment (e.g. in an Erlenmeyer flask) which does not receive any further resource inputs or organims once it has been started.
\[ \begin{aligned} \tfrac{d}{dt} N &= \mu \cdot N \cdot \frac{R}{R+h} \\ \tfrac{d}{dt} R &= -\frac{1}{y} \cdot \mu \cdot N \cdot \frac{R}{R+h} \end{aligned} \tag{5}\]
The parameter \(y\) essentially controls the relation between alterations in biomass and resource. Specifically, it represents the amount of new biomass being formed per amount of consumed resource. Therefore, \(y\) is often denoted as the yield parameter. Depending on the context, it can either be expressed as a mass ratio or a molar ratio.
The last term of Equation 5 is known as a Monod term (Equation 6) with the parameter \(h\) representing a half saturation constant.
\[ \frac{R}{R + h} \tag{6}\]
Obviously, at \(R=h\), the term takes a value of 0.5 and, when multiplied with a growth term as in Equation 5, the population grows at half of the maximum possible rate.
Remember the matrix-based notation for simultaneous ODE we introduced previously? Even for the simple case of Equation 5, it allows for a clearer, more compact notation (Equation 7; the stoichiometry matrix has a single row only).
\[ \begin{bmatrix} \tfrac{d}{dt} N, & \tfrac{d}{dt} R \end{bmatrix} = \begin{bmatrix} \mu \cdot N \cdot \frac{R}{R+h} \end{bmatrix} \cdot \begin{bmatrix} 1, & -1/y \end{bmatrix} \tag{7}\]
The model presented in Equation 5 / Equation 7 is part of suite of demo models distributed with the rodeoEasy package. To find the respective workbook on your local computer, type the R command
print(system.file("models/growth_batch.xlsx", package="rodeoEasy"))and open that file with the spreadsheet software of your choice. For convenience, the contents of the workbook is also displayed below (Table 1 to Table 4).
| description | name | unit | rate | N | R |
|---|---|---|---|---|---|
| growth | grw | cells / mL / h | mu * N * R / (R + h) | 1 | -1/y |
| name | unit | description | default |
|---|---|---|---|
| N | cells / mL | organism density | 10000 |
| R | µg / mL | resource concentration | 1 |
| name | unit | description | default |
|---|---|---|---|
| mu | 1 / h | growth rate constant | 2e+00 |
| h | µg / mL | half saturation constant | 5e-01 |
| y | cells / µg | yield | 1e+09 |
| name | code |
|---|---|
| dummy | dummy <- function() { NULL } |
R code below triggers a simulation with graphical outputs being displayed in Figure 3.
rm(list=ls())
library("rodeoEasy")
workbook <- system.file("models/growth_batch.xlsx", package="rodeoEasy")
m <- build(workbook)
x <- run.scenarios(m, times=seq(0, 12, .1), plot.vars=FALSE)
par(mfrow=c(2,2))
show.scenarios(x, "N") # display with linear axis
show.scenarios(x, "R")
show.scenarios(x, "N", ylog=TRUE) # display with log-scale axis
show.scenarios(x, "R", ylog=TRUE)
par(mfrow=c(1,1))
In the R code above, you may want to try different values for the half saturation constant \(h\) and the parameter \(y\). What are the individual effects of these two parameters?
As depicted by Figure 4, parameter \(h\) controls the impact of resource availability on the actual growth rate. Low values of \(h\) indicate the organisms’ capability to grow fast even if resources are scarce. Larger values of \(h\) result in slower growth at a given resources concentration and thus result in a delay. By contrast, the yield parameter \(y\) exclusively controls the final biomass reached (Figure 5).
rm(list=ls())
library("rodeoEasy")
workbook <- system.file("models/growth_batch.xlsx", package="rodeoEasy")
m <- build(workbook)
sc <- list(low.h=c(h=0.1), medium.h=c(h=0.5), high.h=c(h=0.9))
x <- run.scenarios(m, times=seq(0, 18, .1), scenarios=sc)
# ... continued from previous code section
sc <- list(low.y=c(y=1e+8), medium.y=c(y=1e+9), high.y=c(y=1e+10))
x <- run.scenarios(m, times=seq(0, 18, .1), scenarios=sc)
To understand the true meaning of yield parameter (\(y\)) in the context of nutrient limitation, we need to distinguish between contrasting metabolisms performed by different types of microorganisms like
photoautotrophic eukaryotes such as, e.g., green algae or diatoms,
chemoautotrophic bacteria,
heterotrophic bacteria.
Photoautotrophs use CO2 as the carbon source for primary production which is not a limiting factor in most environments (important exception: acid lakes). Thus, the limiting factors are typically other constituents like phosphor and nitrogen (e.g. in algae). If the biomass is conventionally expressed in units of organic carbon, the parameter \(y\) mostly reflects the biomass stoichiometry. For example, a standard estimate of phytoplankton stoichiometry is a molar C:N:P ratio of about 106:16:1.
Chemoautotrophs equally use CO2 as the carbon source and they also depend on N, P, and further trace elements for primary production. However, this group of organisms can also be limited by the substrate required for energy production like NH4+ in nitrifiers or HS- in sulfide oxidizing bacteria.
The situation is different in heterotrophs were both biomass and and resource are typically expressed in units of organic carbon. Here, the parameter \(y\) is obviously not a simple stoichiometric ratio but it must also account for the loss of carbon through respiration. In this case, it simply specifies the net amount of new biomass being formed per amount of consumed resource.
3.2.2 Monod model as part of a continuous bioreactor model
In Equation 5 and Equation 5, we employed a Monod term to describe resource limitation in a batch reactor, e.g. a bottle which is inoculated with organisms at time zero and thereafter receives no further inputs (neither resource, nor organisms).
By contrast, we will now use a Monod term to represent resource limitation in a continuous bioreactor model based on the CFSTR concept. We will inoculate the reactor with microorganisms just once at the start. Thereafter, we continuously feed the reactor with sterile resource only to replenish what has been consumed by the organisms.
Again, this model is part of the suite of demo models distributed with the rodeoEasy package. To find the respective workbook on your local computer, type the R command
print(system.file("models/growth_continuous.xlsx", package="rodeoEasy"))and open that file with the spreadsheet software of your choice. For convenience, the contents of the workbook is also displayed below (Table 5 to Table 8).
| description | name | unit | rate | N | R |
|---|---|---|---|---|---|
| growth | growth | cells / mL / h | mu * N * R / (R + h) | 1 | -1/y |
| import | import | depends | 1/tau | N_in | R_in |
| export | export | depends | 1/tau | -N | -R |
| name | unit | description | default |
|---|---|---|---|
| N | cells / mL | organism density | 1e+03 |
| R | µg / mL | resource concentration | 3e-01 |
| name | unit | description | default |
|---|---|---|---|
| mu | 1 / h | growth rate constant | 1e+00 |
| h | µg / mL | half saturation constant | 9e-01 |
| y | cells / µg | yield | 1e+03 |
| tau | h | residence time | 1e+01 |
| R_in | µg / mL | resource in inflow | 1e+00 |
| N_in | cells / mL | organisms in inflow | 0e+00 |
| name | code |
|---|---|
| dummy | dummy <- function() { NULL } |
The R code below triggers a simulation of the continuous bioreactor. Apparently, after an adaptation phase, the chosen settings allow for a constant production of biomass (Figure 6).
rm(list=ls())
library("rodeoEasy")
workbook <- system.file("models/growth_continuous.xlsx", package="rodeoEasy")
m <- build(workbook)
x <- run.scenarios(m, times=0:100)
To fully understand the outcome in Figure 6, I propose to prints some numbers corresponding to the final model time step:
# ... continued from previous code section
last <- nrow(x)
print(x[last,"growth"])[1] 90print(x[last,"export"] * x[last,"N"])[1] 90.00155The numbers illustrate that, in steady-state, the rate of biomass production perfectly balances the losses due to the reactor’s outflow. Since the model parameters were chosen arbitrarily, the system obviously has a tendency for self-regulation.
Can this kind of “self-regulation” be observed under all circumstances? In the context of bioreactor operation (but also lake management), this is clearly an important question. Let’s perform a sensitivity analysis to find out. Figure 7 to Figure 10 illustrate the result of local sensitivity analyses where only one parameter has been varied at a time.
# impact of parameter 'R_in'
sc <- list(XS=c(R_in=0.05), S=c(R_in=0.1), M=c(R_in=0.5),
L=c(R_in=0.75), XL=c(R_in=1))
x <- run.scenarios(m, times=0:100, scenarios=sc)
# impact of parameter 'h'
sc <- list(XS=c(h=0.1), S=c(h=1), M=c(h=2),
L=c(h=10), XL=c(h=50))
names(sc) <- c("XS","S","M","L","XL")
x <- run.scenarios(m, times=0:100, scenarios=sc)
# ... continued from previous code section
# impact of parameter 'tau'
sc <- list(XS=c(tau=1), S=c(tau=2), M=c(tau=5),
L=c(tau=10), XL=c(tau=50))
x <- run.scenarios(m, times=0:100, scenarios=sc)
# impact of parameter 'y'
sc <- list(XS=c(y=1e3), S=c(y=1.5e3), M=c(y=2e3),
L=c(y=2.5e3), XL=c(y=3e3))
x <- run.scenarios(m, times=0:100, scenarios=sc)
What do we learn from the sensitivity analyses?
First of all we can see that a steady biomass production does not happen under all circumstances. If the resource input to the system is too low (Figure 7) or the organisms require a very high ambient concentration to grow optimally (Figure 8), the system converges to a trivial steady state where no organisms are present. The same is true if we increase the flow through the system a lot which is equivalent to drastically shortening the residence time (Figure 9). In those cases, the unconditional losses from the system due to outflow can no longer be compensated by growth.
By contrast, the yield parameter can be varied over a very wide range without ever causing a breakdown of biomass production (Figure 10). Further, altering the yield does not affect the residual resource concentration in any way; it only controls the translation of consumed resources into biomass.
3.2.3 Beyond the classical Monod model
The classical Monod model does not necessarily fit all cases. The little R code below demonstrates common derivatives of the classical model (Figure 11).
rm(list=ls())
r <- seq(0, 10, 0.1) # gradient of resource
half <- seq(1, 5, 1) # half-saturation constants
newPlot <- function(label) {
plot(range(r), 0:1, type="n",
xlab="Resource", ylab="Growth rate multiplier")
legend("bottom", bty="n", legend=label)
}
par(mfrow=c(2,2))
newPlot("Classic model")
for (h in half)
lines(r, r / (h + r), lty=match(h, half))
newPlot("With offset 0.8")
off <- 0.8
for (h in half) {
lines(r, pmax(0, (r-off) / (h + r-off)), lty=match(h, half))
}
newPlot("Using power 2")
for (h in half)
lines(r, r^2 / (h^2 + r^2), lty=match(h, half))
# legend
plot(0, 0, type="n", axes=F, ann=F)
legend("center", bty="n", lty=1:length(half),
legend=paste0("h=",half))
par(mfrow=c(1,1))
The Monod model and its close descendants (Figure 11) propose an immediate impact of the ambient resource concentration on the growth rate. Many organisms, however, are capable of taking up limiting resources in excess when they are available to build a stock for “bad times”. Because of this internal stock, cell division may continue for some time even when ambient resources are apparently depleted. The phenomenon is well known from phytoplankton populations growing in spite of no dissolved inorganic phosphors being measurable. A replacement for the Monod model which explicitly accounts for temporary storage was originally proposed by Droop in 1973 and has been discussed extensively later, for example in this paper.
3.3 Multiple potentially limiting resources
Organisms generally depend on multiple resources multiple of which could potentially be limiting. A simple, intuitive solution would be to concatenate multiple limitation terms in a multiplicative manner. However, depending on the situation, it can also be appropriate to consider the minimum of all possible limitation terms. For example, Equation 7 could be modified as below to make the growth dependent on the concentrations of two essential resources identified by subscripts “a” and “b”. (Equation 8).
\[ \begin{bmatrix} \tfrac{d}{dt} N, & \tfrac{d}{dt} c_a, & \tfrac{d}{dt} c_b \end{bmatrix} = \begin{bmatrix} \mu \cdot N \cdot min \left( \frac{c_a}{c_a+h_a}, \frac{c_b}{c_b+h_b} \right) \end{bmatrix} \cdot \begin{bmatrix} 1, & -1/y_a, & -1/y_b \end{bmatrix} \tag{8}\]
4 Growth in the presence of inhibitors
4.1 Inhibition versus lethal effects
(Micro)organisms are often confronted with unfavorable conditions that prevent population growth irrespective of resource limitation. A very good example are antibiotics used to suppress the proliferation of pathogenic bacteria. In that context, one typically distinguishes two major categories of drugs (Table 9).
| Category | Effect | Mechanism | Examples |
|---|---|---|---|
| Bacteriostatic | Inhibition of growth | Reversible inhibition of protein systhesis | Tetracyclines |
| Bactericidal | Killing of cells | Descruction of cell components | beta-lactams |
Thinking in terms of bacterial growth models, it is obvious that the action of a bacteristatic drug could be represented by an multiplicative expression like Equation 9
\[ \tfrac{d}{dt}N = \mu \cdot f(c) \cdot N \tag{9}\]
where the drug concentration is denoted \(c\) and the function \(f(c)\) takes values from 1 (no effect) to 0 (total inhibition). Hence, at best, application of the drug can suspend microbial growth (\(\tfrac{d}{dt}N = 0\)), hereby increasing the chance for the immune system to eliminate the temporarily suppressed pathogens.
By contrast, a different expression is required to represent bactericidal drugs. It could take the form of Equation 10.
\[ \tfrac{d}{dt}N = (\mu - g(c)) \cdot N \tag{10}\]
Again \(c\) denotes the drug concentration. If the function \(g(c) = 0\) (lower limit), the drug concentration is too low to have any effect. At \(g(c) = \mu\), the effect is the same as for bacteriostatic drugs since \(\tfrac{d}{dt}N = 0\). However, at \(g(c) > \mu\), the drug contributes to the killing of microbes actively.
The concept demonstrated with the example of bactericidal and bacteriostatic drugs is generic and similarly applies to other stressors. Thus, depending on the situation, one has to decide whether a formulation like Equation 9 or Equation 10 is more appropriate.
4.2 MIC model
A the classical model to represent a bacteriostatic effects of drugs is based on the minimum inhibitory concentration (MIC). It is a linear model described by Equation 11 with \(c\) denoting the concentration of the drug and \(f(c)\) being a dimensionless multiplier for the growth rate (like in Equation 9).
\[ f(c) = max \left( 0, 1 - \frac{c}{MIC} \right) \tag{11}\]
Obviously, the parameter \(MIC\) specifies the drug concentration \(c\) satisfying \(f(c) = 0\). Hence, as the name says, the MIC is the lowest concentration to entirely prevent growth.
The concept is illustrated with a small data set of Pseudomonas putida exposed to the antibiotic streptomycin. Growth curves were first recorded using a plate reader (based on optical density). The growth rate constants (mu) were then inferred from the steepest slope of the curves considering four replicates. For simplicity, parameters were “fitted by eye” (Figure 12). For serious fitting, look at the R functions nls or optim.
x <- read.table(header=T, text='
# Growth rate constants (mu; 1/hour) for Pseudomonas putida KT2440
# under exposure to streptomycin (conc; mg/L). The bacteria were
# incubated in Mueller-Hinton broth at 30°C with continuous agitation.
# m1 ... m4 represent measurements on four biological replicates.
conc mu1 mu2 mu3 mu4
0.195 1.21 1.15 1.09 1.26
0.39 1.28 1.2 1.19 1.22
0.78 1.17 1.21 1.2 1.26
1.56 1.36 1.25 1.24 1.17
3.125 1.09 0.98 1.18 1.07
6.25 0.55 0.48 0.62 0.71
12.5 0.18 0.34 0.27 0.11
25 0.04 0.14 0.15 0.02
50 0.02 0.03 0.08 0.08
')library(data.table)
x <- melt(data=as.data.table(x), id.vars="conc", # simplify plotting
variable.name="replicate", value.name="mu")
plot(mu ~ conc, data=x, log="x", # observations
xlab="Streptomycin (mg/L)",
ylab="Growth rate constant (1/hour)")
mu <- function(c, mu0, mic) {mu0 * pmax(0, 1 - c/mic)} # MIC model
conc <- seq(0.1, 100, 0.1)
lines(conc, mu(conc, 1.25, 15)) # fitted by eye
If Equation 11 is inserted into a growth expression (see Equation 9) and a term accounting for resource limitation is added as well (see Equation 4), the full model would read like so (Equation 12):
\[ \tfrac{d}{dt}N = \mu \cdot N \cdot \left( 1 - \frac{N}{K}\right) \cdot max\left( 0, 1 - \frac{c}{MIC} \right) \tag{12}\]
4.3 Inverted Monod model
The value of the original Monod model used to describe resource limitation (Equation 6) gradually converges from 0 to 1 for increasing resource concentrations (Figure 11). To describe inhibition, we want the opposite. That is, we need a term whose value gradually declines from 1 to 0 as the concentration of the inhibitor (X) increases. By keeping the core mathematical structure, this can be achieved simply by subtracting a Monod term from one (Equation 13)
\[ 1 - \frac{X}{X + h} \tag{13}\]
which can easily be rearranged for convenience (Equation 14).
\[ \frac{h}{X + h} \tag{14}\]
Like in the original Monod model, the value of the parameter \(h\) has a simple interpretation in the inverted model as well: \(h\) represents the concentration of the inhibitor \(X\) causing a 50% reduction of a biological process rate as compared to a situation where no inhibitor is present (\(X=0\)).
4.4 Flexible dose-response models
The MIC model or the inverse Monod model are attractive for their simplicity. They are not very flexible, however, namely because the single parameter does not allow for the representation of complex relationships.
If more flexibility is needed, one often wants to fit dose-response curves of the logistic or log-logistic family. These latter are frequently employed in ecotoxicology, for example, and the R package “drc” is specializes on their fitting and analysis.