Different with decision model, a dynamics model can us a scientific understanding of a particular pathway. In this section, we develop two different dynamics models, which are Population Dynamics Model and System Dynamics Model respectively. The former we conduct mainly focuses on predicting the population change of our E.coli and gives us the optimal induction time. And the latter tries to describe the whole Astaxanthin pathway mathematically and predict the final yield, thereby providing supports to our further experience.

    During our pre-experiment, we face a difficulty that we need to find out the logarithmic growth phase of our E.coli population because only in this period can we induce them with IPTG. In order to tell our experiment team the best time to induce E.coli with IPTG, we, modeling team, use Logistic Regression to fit pre-experiment data and get the optimal inducing time through our model. Therefore, the amount of experiment can be largely decreased and a lot of time can be saved.

2.1 Basic theory

    Logistic Model is the base of Logistic Regression. It is also known as a famous population dynamics model. In this model, the developing tendency of a population is well described. As a result, we try to use MATLAB2018a to fit our pre-experiment data with this model .

    The basic form of the model is

\begin{equation} \frac{dN(t)}{dt}=rN(t)(1-\frac{N(t)}{K}) \end{equation}

    Where, $N(t)$ represents population density, $r$ represents population growth rate, and $K$ means environmental carrying capacity.

    Further, we have

\begin{equation} \frac{dN(t)}{dt}=aN(t)-bN^{2}(t)\quad\quad(a=r,b=\frac{r}{K}) \end{equation} \begin{equation} \int_{N_{0}}^{N(t)}\frac{dx}{ax-bx^{2}}=\int_{t_{0}}^{t}dt=t-t_{0} \end{equation} \begin{equation} \frac{1}{a}\int_{N_{0}}^{N(t)}[\frac{1}{x}+\frac{b}{a-bx}]dx=\int_{t_{0}}^{t}dt=t-t_{0} \end{equation}

    Finally, we can get the solution that

\begin{equation} N(t)=\frac{K}{1+(\frac{K}{x_{0}}-1)e^{-r(t-t_{0})}} \end{equation}

2.2 Model Solving

    We use MATLAB2018a to fit our pre-experiment data. To estimate the parame- ters in the equations above, we develop a First order difference algorithm (MATLAB code:LogisticSimulate.m).

Step 1 Firstly, to calculate the growth rate $r$, we use an one-order difference formula to describe $r$

\begin{equation} r=\frac{N(t)-N(t-1)}{N(t)} \end{equation}

Step 2 Secondly, In order to use ordinary least squares (OLS) to estimate it, we discrete the continuous model by doing:

\begin{equation} \frac{dN(t)}{dt}=\frac{\Delta N(t)}{\Delta t}=\frac{N(t)-N(t-1)}{\Delta t} \end{equation}

To solve it, we use MATLAB to iterate $N(t)$

\begin{equation} N^{'}(k+1)-N^{'}(k)=r^{'}\cdot N^{'}(k)\cdot (1-\frac{N^{'}(k)}{M}) \end{equation}

Where, $N^{'}$ is a empty vector with the same order as $N(t)$, $r^{'}$ is a empty vector with the same order as $r(t)$, and k is number of iterations.

Step 3 Calculate variance and plot.

    The solution of the formulas above are shown in figure 1, in which (a) is Logistic model simulation effect diagram and (b) is Simulate error of Logistic model simulation, and the value of $r$, $K$ is $r=0.8606$ and $K=3.0238$.

    As a result, we obtain the best population concentration to induce is about OD600=$\frac{K}{2}=1.5$, where the population density is in log phase.

3 System Dynamics Modeling

    In this section, we try to simulate our pathway mathematically. We first consider each reaction as an independent part. Further, we assume that each reaction in the six-step pathway is consistent with the Michaelis-Menten equation. Therefore, the basic model of the reaction chain can be obtained according to the Michaelis-Menten equation along with the law of mass action. In addition, the basic model is improved modified due to the competitive effects of different substrates on the same enzyme. By using MATLAB 2018a to solve the model above, the curve of the yield of each product with time can be gotten.

3.1 Basic theory

3.1.1 Michaelis-Menten equation

    In biochemistry, the Michaelis-Menten equation is one of the most famous profound kinetic model. The equation mainly describes the behavior of a typical enzyme in a enzymatic action. For the following enzymatic reactions:

\begin{equation} E+S\stackrel{{k_{f}/k_{r}}}{\longleftrightarrow}ES\stackrel{k_{cat}}{\longrightarrow}E+P \end{equation}

    Where, $k_{f}$ (forward rate), $k_{r}$ (reverse rate) and $k_{cat}$ (catalytic rate) provide rate constants to the equation; the left-right arrow between S (substrate) and ES (enzyme-substrate complex) means that the the binding process of enzyme-substrate complex is reversible; and the right arrow represents the release of P (product) If the postulation that the enzyme concentration is much lower than the substrate concentration is accepted, the reaction speed is:

\begin{equation} v=\frac{d[P]}{dt}=\frac{k_{cat}\cdot[S]\cdot[E]}{k_{m}+[S]} \end{equation}

    Where, [$S$] is substrate concentration, [$E$] is enzyme concentration, and [$P$] is product concentration.

3.1.2 Law of mass reaction

    The law of mass action is: the rate of the elementary reaction is proportional to the power of the concentration of each reactant, where in the power exponent of each reactant concentration is the stoichiometry of the reactant in the reaction equation, and the proportional constant is called the chemical reaction constant. To illustrate the law of mass reaction, we take the following biochemical reactions as an example:

\begin{equation} A+B\stackrel{k}{\longrightarrow}C \end{equation}

    We have:

\begin{equation} \frac{d[C]}{dt}=-\frac{d[A]}{dt}=-\frac{d[B]}{dt}=v \end{equation} \begin{equation} v=k\cdot[A]\cdot[B] \end{equation}

    Where, $v$ is called flux of chemical reaction.

3.2 Basic Model

    According to the basic theory above, along with the postulation that the concentration of enzyme is a constant, our pathway (Figure 2) can be described mathematically as:

\begin{equation} v_{x_{1}}=-\frac{dx_{1}}{dt}=-\frac{k^{crtE}_{cat}\cdot x_{1}\cdot [crtE]}{k_{m}^{crtE}+x_{1}} \end{equation} \begin{equation} v_{x_{2}}=\frac{dx_{2}}{dt}=-v_{x_{1}}-\frac{k^{crtI}_{cat}\cdot x_{2}\cdot [crtI]}{k_{m}^{crtI}+x_{2}} \end{equation} \begin{equation} v_{x_{3}}=\frac{dx_{3}}{dt}=v_{x_{2}}-\frac{k^{crtB}_{cat}\cdot x_{3}\cdot [crtB]}{k_{m}^{crtB}+x_{3}} \end{equation} \begin{equation} v_{x_{4}}=\frac{dx_{4}}{dt}=v_{x_{3}}-\frac{k^{crtY}_{cat}\cdot x_{4}\cdot [crtY]}{k_{m}^{crtY}+x_{4}} \end{equation} \begin{equation} v_{x_{5}}=\frac{dx_{5}}{dt}=v_{x_{4}}-\frac{k^{crtZ}_{cat}\cdot x_{5}\cdot [crtZ]}{k_{m}^{crtZ}+x_{5}} \end{equation} \begin{equation} v_{x_{6}}=\frac{dx_{6}}{dt}=v_{x_{5}}-\frac{k^{BKT}_{cat}\cdot x_{6}\cdot [BKT]}{k_{m}^{BKT}+x_{6}} \end{equation} \begin{equation} v_{x_{7}}=\frac{dx_{7}}{dt}=\frac{k^{BKT}_{cat}\cdot x_{6}\cdot [BKT]}{k_{m}^{BKT}+x_{6}} \end{equation}

    Where, $x_{1}$ to $x_{7}$ represent FPP to Astaxanthin respectively. Canthaxanthin, however, is not included. We will talk about this branch pathway in next section.

3.3 Improved Model

3.3.1 Competitive effects of two substrates on the same enzyme

    Although we used the relationship between the Michaelis-Menten equation and the chain reaction to obtain the expression of Astaxanthin, the model only considered the $(\beta -Carotene)-(Zeaxanthin)-(Astaxanthin)$ pathway. However, there is a branch pathway $(\beta -Carotene)-(Canthaxanthin)-(Astaxanthin)$. In this process, the competitve effect must appear.

\begin{equation} E+S\stackrel{k_{f}/k_{r}}{\longleftrightarrow}ES\stackrel{k_{cat1}}{\longrightarrow}E+P \end{equation} \begin{equation} E+C\stackrel{k_{f}/k_{r}}{\longleftrightarrow}EC\stackrel{k_{cat2}}{\longrightarrow}E+P \end{equation}

    Where, C means competitive substrate and EC represents enzyme-competitive substrate complex. As the equation shown above, two substrates compete the same enzyme, which can be described as the following competitive inhibition kinetic equation:

\begin{equation} [E_{f}]=[E]-[ES]-[EC] \end{equation} \begin{equation} k_{m}=\frac{k_{m}+k_{cat}}{k_{f}} \end{equation} \begin{equation} k_{1}\cdot [E_{f}]\cdot [S]=(k_{r}+k_{cat})\cdot [ES] \end{equation}

    By combining the eqautions above, we can get the improved reaction speed:

\begin{equation} v^{'}_{x_{5}}=v_{x_{4}}-\frac{k^{crtZ}_{cat}\cdot x_{8}\cdot [crtZ]}{k_{m}^{crtZ}(1+\frac{x_{8}}{k_{m}^{crtZ}})+x_{5}} \end{equation} \begin{equation} v^{'}_{x_{6}}=v_{x_{5}}-\frac{k^{crtZ}_{cat}\cdot x_{5}\cdot [crtZ]}{k_{m}^{crtZ}(1+\frac{x_{5}}{k_{m}^{crtZ}})+x_{6}} \end{equation}

    Where, $v^{'}_{x_{5}}$ and $v^{'}_{x_{6}}$ represent improved $v_{x_{5}}$ and $v_{x_{6}}$ respectively.

3.3.2 Product inhibition effect

    Because of the product inhibition effect, the producing speed of Astaxanthin($x_{7}$) can be described as:

\begin{equation} v_{x_{7}}=\frac{dx_{7}}{dt}=\frac{k^{BKT}_{cat}\cdot x_{6}\cdot [BKT]}{k_{m}^{BKT}(1+\frac{x_{7}}{k_{i}})+x_{6}} \end{equation}

    Where, $k_{i}$ is a product inhibition constant.

3.4 Model Solving

    We use MATLAB 2018a to solve the ordinary differential equations above. The following parameters were used:

    The results are shown in Figure 3 where (a) is the curve of the concentration of intermediate products, (b) is the the curve concentration of Astaxanthin.

    From the figure 3, we can expect that the final concentration of Astaxanthin is 5$mg/mL$.

1. Zhengshan Lin. Population dynamics. Beijing:Science Press. 2006
2. BRENDA[DB/OL]. http://www.brenda-enzymes.org/index.php
3. Park SY, Binkley RM, Kim WJ, Lee MH, Lee SY. Metabolic engineering of Escherichia coli for high-level astaxanthin production with high productivity.[J] Metab Eng. 2018.9