Team:Jilin China/Model/degradation model

Degradition Model

Overview

Our degradation model contained a set of biochemical reactions, as other metabolic models might be. We generated several ordinary differential equations from these reactions with unknown kinetic constants. These equation sets might not have analytic solutions and numeric methods always fail when dealing with unknown parameters in differential equations. One way to fix such problem is to make an approximation. For example, if a reaction X is not highly reversible, the kinetic constant kX for that reaction would be regarded as zero. However, if any of the reaction in a chain process is not that irreversible, such approximations could also fail. This year, we involved a parameter optimizing method and combined the optimization with numeric solvers, which allowed us to get more accurate values for the kinetic constants from an initial guess. As soon as we get the accurate values for the kinetic constants, the degradation model could then predict the degradation rate of our pollutant. We’ve developed this method (we call it a numeric solving – parameter optimizing method) for our future work – since the enzyme could deal with many kinds of pollutant (theoretically), we couldn’t test them all and fit the model with experimental data.

Mathematical modelling

Once the population of the engineered E. coli reached a stable level (and the expression of enzymes would reach a stable level as well), the degradation rate of the pollutant could be measured. The degradation process contains the following progress:


$\mathrm{TfdB}+\Phi\rightleftharpoons\mathrm{TfdB}\Phi$      $k_{1+},k_{1-}$

$\mathrm{TfdB}\Phi\rightleftharpoons\mathrm{TfdB}\Phi'$      $k_{2+},k_{2-}$

$\mathrm{TfdB}\Phi'\rightleftharpoons\mathrm{TfdB}+\Phi'$      $k_{3+},k_{3-}$

$\mathrm{CphA_{-1}}+\Phi'\rightleftharpoons\mathrm{CphA_{-1}}\Phi'$      $k_+',k_-'$

$\mathrm{CphA_{-1}}\Phi'\rightarrow\mathrm{CphA_{-1}}+\Omega$      $k'$

Where Φ represents the pollutant, Φ' represents the intermediate product (the production from TfdB) and Ω represents the final production.

Parameters
VariableExplanation
$[\mathrm{substance}]$concentration of the substance
$t$time
ConstantExplanation
$[\mathrm{TfdB}]$concentration of TfdB
$[\mathrm{CphA_{-1}}]$concentration of CphA-1
$k_{?}$rate constant, see reaction equations
Equations

$\begin{cases} \dfrac{\mathrm{d}[\Phi]}{\mathrm{d}t}=k_{1-}[\mathrm{TfdB}\Phi]-k_{1+}[\mathrm{TfdB}][\Phi]\\\\ \dfrac{\mathrm{d}[\mathrm{TfdB}\Phi]}{\mathrm{d}t}=k_{1+}[\mathrm{TfdB}][\Phi]+k_{2-}[\mathrm{TfdB}\Phi']-(k_{2+}+k_{1-})[\mathrm{TfdB}\Phi]\\\\ \dfrac{\mathrm{d}[\mathrm{TfdB}\Phi']}{\mathrm{d}t}=k_{3-}[\mathrm{TfdB}][\Phi']+k_{2+}[\mathrm{TfdB}\Phi]-(k_{3+}+k_{2-})[\mathrm{TfdB}\Phi']\\\\ \dfrac{\mathrm{d}[\Phi']}{\mathrm{d}t}=k_{3+}[\mathrm{TfdB}\Phi']+k_-'[\mathrm{CphA_{-1}}\Phi']-k_+'[\mathrm{CphA_{-1}}][\Phi']-k_{3-}[\mathrm{TfdB}][\Phi']\\\\ \dfrac{\mathrm{d}[\mathrm{CphA_{-1}}\Phi']}{\mathrm{d}t}=k_+'[\mathrm{CphA_{-1}}][\Phi']-(k'+k_-')[\mathrm{CphA_{-1}}\Phi']\\\\ \dfrac{\mathrm{d}[\Omega]}{\mathrm{d}t}=k'[\mathrm{CphA_{-1}}\Phi'] \end{cases}$
Model fitting

This model could be used on a wide range of substances -- not only various kinds of phenol, but also indole and many other pollutants. For the model fitting part, our aim was not to fit this model for a specific kind of pollutant -- we aimed to involve a fitting method which would be used in our future work. We combined Runge-Kutta method and simulated annealing method. Runge-Kutta method is a widely used numeric method for differential equation solving, this method could only work with given parameter values. Simulated annealing is a optimization method based on a physical process of annealing. We planned to start the fitting with initial guesses for the kinetic constants -- we could get this guess from given homologous protein kinetic rules. The initial values of kinetic constants may not fit our datas well, however, we could randomly change the values and compare them with our experimental datas and calculate the degradition kinetic curve with the RK-4 algorithm. Based on the least squares calculations, we could give up some of the changes (or, we can say 'annealing try'). The annealing tries lead to extremely big least squares values would be given up. The rule to judge whether an annealing try should be abandoned or not is to compare it with the former result. If the former values for the kinetic constants lead to greater least squares, the current try would be accepted to take the place of former values; if the current try could not give smaller least squares, we have to calculate whether it is 'worse enough' -- we might accept it with a 'certain probability':

$p=e^{-\frac{(\Delta L)}{T}}$

where $p$ referes to the possibility, $e$ is the base of natural logs, $\Delta L$ is the difference between the least squares from the current values and the former values, $T$ is the 'current temperature'. Current temperature is involved here to avoid bad annealing tries for the last steps -- in the physical world, though the beginning annealing tries would move particals unexpectedly (to avoid the local best traps -- particals might get through energy barriers to reach their stable positions), particals might not move so far away to an unstable position again in the last several steps. In the last several steps, since $T$ would have a quite small value ($T$ decrease with the progress of annealing), any of positive values of $\Delta L$ will cause a high possibility of rejection.

Algorithm Coding

We designed our whole algorithm with Matlab®, however, to show the whole principle of our algorithm, we didn't use Matlab build-in functions for the key steps in our work. Click here to view our codes.