Having a kinematic model that mimics the behavior of the ceiling tile manufacturing plant's process water is a valuable design and analysis tool to investigate potential solution strategies such as the proposed synthetic biology solution of engineered yeast. As a first step towards this goal of a simulated design and testing environment, basic mathematical models are developed here to characterize both the kinetics of yeast cell growth and the enzymatic reaction of starch degradation by alpha-amylase. Such a kinetic model for microbial growth and product formation must account for the different phases of growth, most importantly an exponential growth phase during cell seeding,a stationary phase relating to steady-state operation, and a death or cell decline phase associated with conditions of nutrient deficiency or toxic products.
Model Formulation
Several models can be found in the literature [1-4] to describe cell growth, substrate utilization, and product formation. The most widely used models are simpler unstructured models that incorporate the Monod equation to represent the logistic growth characteristics typical of microbials. This model is mathematically derived for cell growth with the following differential equations.
Cell Growth
The rate of cell growth is dependent on cell concentration \(X\). Expressed as
\[ \dfrac{dX}{dt} = \mu_g X \]
the specific growth rate \(\mu_g\), having units of hr -1, is defined as
\[\mu_g = \dfrac{1}{X} \dfrac{dX}{dt}\]
If the rate of cell loss from cell death and endogenous metabolism is \(K_d\), the cell death rate \(R_d\) and the net cell growth can be written as
\[R_d = K_d X\]
\[\dfrac{dX}{dt} = \mu_g X - K_d X\]
Rearranging terms, the net specific growth \(\mu\) is
\[\mu = \mu_g - K_d = \dfrac{1}{X} \dfrac{dX}{dt}\]
Taking integrals on both sides
\[ \int_0^t \mu \,dt = \int\dfrac{1}{X}\,dX \]
\[\
\ln X = \mu t + C \qquad \text{where}\qquad C = \ln {X_0}
\]
\[\ \ln \frac{X}{X_0} = \mu t \]
and \[\ X = X_0 e^{\mu t}\]
where \(X_0\) is the initial cell concentration. It follows from the above equation that a graph of \( \ln \frac{X}{X_0}\)plotted against time will display \(\mu\) in its instantaneous slope. When growth is exponential, \(\frac{dX}{dt} = \mu_g X \) and the graph is linear with a slope of \( \mu_g\). At steady state conditions, \(\frac{dx}{dt} = 0 \), and cell growth approaches cell loss, so that \(\mu_g X = K_d X\). During this stationary phase, the graph plateaus with horizontal slope. In the cell death or decline phase where \(K_dX >> \mu_g X\), the graph decreases with negative slope. With first-order kinetics for cell decline,
\[\ \dfrac{1}{X} \dfrac{dX}{dt} = -K_d\qquad \text{ }\qquad \ln \frac{X}{X_0} = -K_d t \]
and \[\ X = X_0 e^{-k_d t}\]
Substrate Utilization
Substrate is utilized for producing new cells and for maintaining cells. If the cell maintenance coefficient is \(m\) and the dimensionless yield factor \(Y_{X / S}\) is the ratio of new cells mass to substrate consumed, that is
\[\ Y_{X / S} = \dfrac {\Delta X}{\Delta S}\]
then similar to the cell growth rate expression of equation (?), the rate of substrate utilization can be expressed as
\[\dfrac {dS}{dt} = -\mu_g \frac {1}{Y_{X / S}} X - mX \]
Product Formation
Similar to cell growth and substrate utilization, products formed can be expressed by a first-order differential equation. The specific rate of product formation, \(\frac{1}{X}\frac {dP}{dt}\) is proportional to the specific rate of cell growth \( \mu_g \) during the exponential growth phase whereas in the non-growth steady-state phase it is a constant, \(\beta\). If the yield ratio \(Y_{P / X}\) is defined as
\[\ Y_{P / X} = \dfrac {\Delta P}{\Delta X}\]
then the rate of product formation is expressed as
\[\dfrac {dP}{dt} = -\mu_g Y_{P / X} X + \beta X \]
Monod Equation for Limited-Substrate Growth
More often than not, cell growth saturates when a substrate becomes limited. If this growth-limiting substrate is S, the kinetics of limited substrate growth are commonly described by the Monod equation [1-4]:
\[ \mu_g = \frac {\mu_{m} S}{K_s + S}\]
where \(K_s\) is called the saturation coefficient. The Monod equation has the same form as another equation that describes the kinetics of enzymatic reactions and called the Michaelis-Menten equation, which is derived in the next section. The Monod equation is suited for microbial growth such as yeast where the growth is not rapid and when cell concentrations are low.
Michaelis-Menten Kinetics for Enzymatic Reactions
For enzyme-catalyzed reactions such as starch degradation by alpha-amylase, reaction kinetics are modeled with the Michaelis-Menten rate expression. The Michaelis-Menten kinetics applies to a single substrate \(S\) and single enzyme \(E\) reaction. If \(ES\) is the enzyme-substrate complex and \(P\) is the product, the reaction is shown by
\[S + E\mathrel{\mathop{\rightleftharpoons}^{k_{1}}_{k_{-1}}}ES\mathop{\rightarrow}^{k_{2}}P + E\]
By writing three coupled differential equations (not shown here) for the rate of product formation, \(\dfrac{d}{dt}C_p\), the rate of enzyme-substrate formation, \(\dfrac{d}{dt}C_{ES}\), and the rate of substrate consumption,
\(\dfrac{d}{dt}C_s\), and solving these equation
the velocity of the forward reaction defined as \(\nu\) according to
\[ \nu = rate \: of \: S \mathop{\rightarrow} P = \dfrac{d}{dt}\left [ C_p \right ] \]
one arrives at the well known Michaelis-Menten rate equation
\[ \nu = \frac{\nu_{max} C_s}{K_m + C_s} \]
where \(K_m\) is an equilibrium concentration given by
\[\ K_m=\frac{k_{-1}+k_2}{k_1} \]
\(C_s\) is the substrate concentration and \( \nu_{max}\) is the maximum rate given by
\[ \nu_{max} = K_2 C_E(0)\]
It should be noted that of the two model parameters, while \(K_m\) is an intrinsic constant, \(\nu_{max}\)depends on the initial enzyme concentration \(C_E(0)\).Since the rate of product formation \(\nu\), equals the rate of substrate consumption.
\[ \nu = \dfrac{d}{dt}C_p\ = -\dfrac{d}{dt}C_s \]
It can be seen from the equation that when
\[ C_s = K_m \qquad \text{ }\qquad \nu = \frac {1}{2} \nu_{max} \]
Model Summary
The two parts of the model describing the kinetics for cell growth and the enzymatic reaction is summarized below. The first three differential equations are related to cell growth and are obtained by substituting for \(\mu_g\) in equations 13 and 4 with equation 16, the Monod equation. The fourth differential equation is the Micaelis-Menten equation for the enzyme reaction.
\[\dfrac{dX}{dt} = \frac {\mu_{m} S}{K_s + S} X - K_d X\]
\[\dfrac{dS}{dt} = - \frac {\mu_{m} S}{(K_s + S)Y_{X / S}} X - M X \]
\[\dfrac {dP}{dt} = \frac {\mu_{m} Y_{P / X} S}{K_s + S} X + \beta X \]
\[\dfrac {dC_s}{dt} = - \frac{\nu_{max} C_s}{K_m + C_s} \]
Estimating Model Parameters from Experimental Measurements
Two parameters relating to the Monod equation, the maximum specific growth rate μ, and the saturation constant Ks are estimated graphically. These graphical methods are based on our readings of the literature [1, 3]. First, the cell growth described in equation 16, the Monod equation, is linearized similar to a linearization of the Michaelis-Menten equation.
\[\dfrac {1}{\mu} = \frac {K_s}{\mu_m}.\frac{1}{S}+\frac{1}{\mu_m} \]
Where
\[ \frac{1}{\mu} = \frac{X}{\frac{dX}{dt}} \]
This takes the form of a double-reciprocal Line-Weaver Burk plot of 1/μ vs 1/S. From the slope and intercept of this plot, model parameters μ and Ks are estimated. The graph below shows the estimated parameters derived from the experimental data for the genetically modified yeast and wild type yeast in a glucose substrate (YPD yeast media).
Figure 1 Double-reciprocal Line-Weaver Burk plot for Cell Growth To check our work, we also used a simpler method, a plot of the log ratios of cell concentrations against time. This method applies to the exponential phase of growth and is useful for getting an estimate of μ from the slope of the graph, according to Equation 8. The log plots are shown below. The figure also captures a summary Table of model parameters.
Figure 2 Log plot of exponential cell growth to estimate the maximum specific cell growth rates for the genetically modified and wild type yeast in a glucose substrate, and the genetically modified yeast in starch substrate. The Table summarizes the model parameters from the two graphical methods. The growth yield Yx/s is obtained from the experimental measurements; calculated as the ratio of the differences between cell densities and glucose concentration for two time points, at time zero and the time point corresponding to the end of exponential growth
MATLAB for Parameter Estimation and Model Simulation
Since the graphical analysis shown above allowed only some of the model parameters to be estimated, we went to MATLAB to accomplish estimating all the parameters. The problem of parameter estimation is a nonlinear optimization problem. It is an iterative process where one starts with an initial estimate of the parameters, runs the forward problem of solving the set of ordinary differential equations (ODE) with these parameters, and then compare the output of the ODE solver with the corresponding experimental measures in a least squares optimizing cost function, after which a new set of model parameters are generated and the process repeats itself until an error criteria is met for the differences between the experimental and model generated outputs. Implementing this in MATLAB caused difficulties for us as we were still new to MATLAB. We reached out to John Marken of the William and Mary iGEM team. With his great help he offered, we made some progress, however, the delay waiting for our final experimental cell growth data to run the optimization, set us back in our modeling schedule. Our attempts with MATLAB to solve the iterative optimization problem is below. We will continue to work on implementing this nonlinear optimization problem in MATLAB.
Future Modeling
To model the ceiling tile plant closer to reality, the model must incorporate both the yeast and the bacteria that are in the process water and are competing for the same nutrients and oxygen. For this, an extended model will be developed. However, the cell growth model developed here also works for bacteria with the only difference being in the model parameters. Bacteria model parameters will be estimated from experimental data to represent the bacteria prevalent in the ceiling tile process water and wastewater. Other factors not included in this model are variations in pH, temperature, and DO. Moving the model to include these variable will be a next step. The ultimate goal of the model is to support the design and development of an implementable solution.
References
1. Shuler, M. L., & Kargi, F. (2002). Bioprocess engineering: Basic concepts, Chapter 4 and Chapter 6. Upper Saddle River, NJ: Prentice Hall.
2. Zangirolami, T.C., Carlsen, M., Nielsen, J., & Jørgensen, S.B.. (2002). Growth and enzyme production during continuous cultures of a high amylase-producing variant of Aspergillus Oryzae. Brazilian Journal of Chemical Engineering, 19(1), 55-68
3. Wang, L., D. Ridgway, T. Gu and M.Y. Murray, 2009. Kinetic modeling of cell growth and
product formation in submerged culture of recombinant Aspergillus niger. Chem. Eng. Com.,
196:481-490.
4. Akpa J (2012) Modeling of a Bioreactor for the Fermentation of Palm wine
by Saccaharomyce cerevisiae (yeast) and lactobacillus (bacteria). Bioresource
Technology 3: 231-240.
