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).
Figure 1: Positive and negative components controlling the population size within system boundaries.
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).
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 Growth with limited resources
2.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 \(K\), is cells mL-1.
Equation 4 allows to account for resource limitation without introducing the resource itself as a second state variable. This is especially useful if the actual mechanism of limitation is not known 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=K\). An analytic solution of Equation 4 can still be found (Equation 5), but the integration requires knowledge about partial fraction decomposition which is clearly beyond the scope of this course.
Since we already learnt how to integrate ODE numerically, we do not depend on the analytic solution. The R code below illustrates how a dynamic solution of Figure 2 can be obtained using the deSolve package.
rm(list=ls())library("deSolve")dydt <-function(time, vars, pars) {with(as.list(c(vars, pars)),list(mu_max * N * (1- N/K)) # ODE right hand side )}init <-c(N=1e4) # initial valuepars <-c(mu_max=2, K=1e9) # parametersdynam <-lsoda(y=init, parms=pars, # integrationfunc=dydt, times=seq(0, 18, .1))par(mfrow=c(1,2)) # visualizationplot(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))
Figure 2: Dynamic solution of the logistic growth model plotted with a linear y axis (left) and a log-scaled y axis (right). The chosen initial value and parameters are reprentative for bacterial growth in rich medium at close to optimum temperatures.
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}\), reflecting depletion of the resources.
In the microbial context, the population dynamics is typically plotted on a logarithmic scale (Figure 2, right) to cover an increase of cell numbers by 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.
2.2 Monod-type models
2.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 any particular 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 6) 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 organisms once it has been started.
\[
\begin{aligned}
\tfrac{d}{dt} N &= &\mu_{max} \cdot N \cdot \frac{R}{R+h} \\
\tfrac{d}{dt} R &= -\frac{1}{y} &\cdot \mu_{max} \cdot N \cdot \frac{R}{R+h}
\end{aligned}
\tag{6}\]
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 fraction \(R / (R + h)\) appearing in Equation 6 is known as a Monod term with the parameter \(h\) representing a half saturation constant. The meaning becomes obvious when rewriting the term as in Equation 7.
The fraction is essentially a dimensionless multiplier to derive the actual growth rate under resource limitation (\(\mu\)) from the maximum, unlimited rate (\(\mu_{max}\)). Under all circumstances, the term \(R / (R + h)\) takes a value in the interval from zero (\(R=0\)) to one (\(R \rightarrow \infty\)). When \(R=h\), we get \(\mu = \mu_{max} / 2\) and the population grows at half the maximum possible rate (Figure 3).
Figure 3: Value of the term from Equation 7 for different assumptions on the parameter \(h\).
Now that we have a general understanding of the Monod term (Equation 7), let’s implement the model from Equation 6 in R.
Figure 4: Prediction of population growth with resource limitation described by a Monod term (Equation 6). Data are presented with a linear y axis.
Figure 5: Data from Figure 4 replotted using logarithmic y axes.
As illustrated by Figure 4 and Figure 5, the Monod model predicts a similar dynamic as the carrying capacity approach (Figure 2). At the same time, the Monod model is more complex as it needs two state variables and involves two parameters (the half-saturation constant \(h\) and the yield parameter \(y\)). Still, the Monod model is typically the superior approach because
the limitation of growth can be explicitly linked to a specific resource,
Monod terms can be combined to reflect a multitude of potentially limiting resources, including co-limitation (more on that in a subsequent section),
the two parameters have a clear interpretation.
Let’s develop a better understanding for the effects of the two parameters by means of a local sensitivity analysis (i.e. by varying one of the parameters at a time). The following code illustrates the impact of the half saturation constant \(h\) on the predicted dynamics with results displays in Figure 6.
x <-NULL# empty result matrixfor (h inc(0.1, 0.5, 0.9)) { # chose values to be tested s <- deSolve::lsoda(y=vars, # run with updated parametertimes=seq(0, 18, 0.1), func=derivatives,parms=c(h=h, pars[names(pars) !="h"])) s <-cbind(s, h=h) # add parameter value as column x <-rbind(x, s) # append to result matrix}par(mfcol=c(1,2)) # plot resultsfor (v inc("N","R")) {plot(range(x[,"time"]), range(x[,v]), type="n", xlab="Time", ylab=v)for (h inunique(x[,"h"])) { sel <- x[,"h"]==hlines(x[sel,"time"], x[sel,v], col=match(h, unique(x[,"h"]))) }}legend("topright", bty="n", lty=1, col=1:length(unique(x[,"h"])),legend=paste0("h=", unique(x[,"h"])))par(mfcol=c(1,1))
Figure 6: Impact of different values of the half saturation constant on the population dynamics prediced by Equation 6. The model addresses a batch reactor without in- or outflow.
As depicted by Figure 6, the 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.
Note that the final biomass is not affected by the choice of \(h\) at all. The latter depends solely on the yield parameter \(y\). This can be illustrated easily by modifying the code from above such that \(y\) is varied instead of \(h\) (Figure 7).
Figure 7: Impact of different values of the yield parameter on the population dynamics prediced by Equation 6. The model addresses a batch reactor without in- or outflow.
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.
2.2.2 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 8).
Figure 8: The classical Monod model and common alternatives. The occurence of an offset is pretty common because organisms cannot uptake resources at arbitrarily low levels, e.g. for energetic reasons. The application of a power allows for fine-tuning of the shape of the curve and it may also help to circumvent numerical issues during ODE integration.
The Monod model and its close descendants (Figure 8) 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.
2.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, considering Liebig’s law, it may be more appropriate to consider the minimum of all possible limitation terms. For example, Equation 6 could be modified as below to make the growth dependent on the concentrations of two essential resources identified by \(R_a\) and \(R_b\). (Equation 8).
\[
\begin{aligned}
\tfrac{d}{dt} N &= &\mu_{max} \cdot N \cdot min \left( \frac{R_a}{R_a+h_a}, \frac{R_b}{R_b+h_b} \right) \\
\tfrac{d}{dt} R_a &= -\frac{1}{y_a} \cdot &\mu_{max} \cdot N \cdot min \left( \frac{R_a}{R_a+h_a}, \frac{R_b}{R_b+h_b} \right) \\
\tfrac{d}{dt} R_b &= -\frac{1}{y_b} \cdot &\mu_{max} \cdot N \cdot min \left( \frac{R_a}{R_a+h_a}, \frac{R_b}{R_b+h_b} \right)
\end{aligned}
\tag{8}\]
3 Growth in the presence of inhibitors
3.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 1).
Table 1: Major categories of antibiotic drugs with examples. See this paper for a cricical discussion suggesting that the two categories cannot be distinguished under all circumstances.
Category
Effect
Mechanism
Examples
Bacteriostatic
Inhibition of growth
Reversible inhibition of protein systhesis
Tetracyclines
Bactericidal
Killing of cells
Destruction 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_0 \cdot f(c) \cdot N
\tag{9}\]
where the drug concentration is denoted \(c\), the growth rate without exposure is \(\mu_0\), 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_0 - g(c)) \cdot N
\tag{10}\]
If the function \(g(c) = 0\) (lower limit), the drug concentration is too low to have any effect. At \(g(c) = \mu_0\), the effect is the same as for bacteriostatic drugs since \(\tfrac{d}{dt}N = 0\). However, at \(g(c) > \mu_0\), 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.
3.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).
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 9). 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 mu40.195 1.21 1.15 1.09 1.260.39 1.28 1.2 1.19 1.220.78 1.17 1.21 1.2 1.261.56 1.36 1.25 1.24 1.173.125 1.09 0.98 1.18 1.076.25 0.55 0.48 0.62 0.7112.5 0.18 0.34 0.27 0.1125 0.04 0.14 0.15 0.0250 0.02 0.03 0.08 0.08')
Figure 9: Application of the MIC model to P. putida growing with streptomycin. Note that the model only appears non-lineare because of the log-scaled lower axis.
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)
with \(\mu_{max,0}\) being the growth rate for optimum resource supply and no drug exposure.
3.3 Inverted Monod model
The value of the original Monod model used to describe resource limitation (Equation 7) gradually converges from 0 to 1 for increasing resource concentrations (Figure 8). 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)
which can easily be rearranged for convenience (Equation 14).
\[
\mu = \mu_0 \cdot \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 in growth (or another biological process) as compared to a situation where no inhibitor is present (\(X=0\)).
3.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.
