Team:UCAS/Model
Overview
Modeling is the art of description, estimation, prediction and optimization with limited knowledge and experimental data, which is exactly what we have to face. Our goal is to build engineering bacteria with antibiotics degradation and procedural suicide genes obtained via screening, and many parameters are not available or not easy to search. The industrial appliance of our idea would be seriously constrained if too many parameters have to be measured. Therefore, we employ ingenious tactics to build a practicable model of antibiotics degradation protein expression, bacteria growth under toxin-antitoxin coupling, and the overall circuit by in-depth analysis with least experiments.
We divide our model into 3 sections. In part 1, we employ numerical algorithm to describe bacteria growth with accuracy, and perform a deep-going analysis, using statistical methods, to show the influence of toxin and toxin-antitoxin coupling. Part 2 deals with the expression of tetracycline resistance protein, we use piecewise linear method for simplification, eigen decomposition method to extract the dominant sector. In part 3 we integrate the two parts into a complete model, with which we do simulations to reach our conclusions.
Bacteria growth and protein expression are the most fundamental issues in synthetic biology. Our model minimizes the necessary experiments to get all necessary parameters, which is conducive to the industrial appliance and also propagable to modeling of other projects.
Growth With TA
Dynamics of Free Growth
To describe growth of bacteria, many growth curve have been proposed. Under substrate-sufficient conditions, the Logistic curve, Gompertz equation, polynomial fitting, Stannard model, and their modifications[1] are widely applied. Meanwhile, Monod equation[2], Baranyi & Roberts model and many computational models[3][4] are more suitable in substrate-limited conditions.
Fig. 1
OD600 (red) and growth rate (blue) curve with error bar. We add 1.5μL bacteria solution into 150μL LB culture medium at 310K and measure the optical density at 600nm every ten minutes, with medium being shaking all along. Growth rate are calculated by five-point method. Error bars show the standard deviation of fifteen parallel samples.
Fig. 1 shows the curves of OD and growth rate. Since the quantity of substrate is limited in our experiments and the parameters should have clear physical meaning for further analysis, we adopt Baranyi & Roberts model[5] as follows：
\frac{\text{dN}}{\text{dt}}=\frac{v\left( t \right)}{1+v\left( t \right)}{{\mu }_{m}}\times \left( 1-\frac{N}{{{N}_{m}}} \right)\text{N}\quad(1)
Parameter v(t) reflects the physiological state of the cells, which is determined by initial conditions and is independent of the history of bacteria and previous environments [6]. It can be assumed that the derivative of v is a constant under invariant conditions. Therefore v(t) can be solved:
v\left( t \right)={{v}_{0}}\exp qt\quad(2)
Optical density(OD) we measure in experiments is proportional to the cell density. Then the growth rate equation can be expressed as followings.
\text{ }\!\!\mu\!\!\text{ }=\frac{1}{OD}\frac{\text{dOD}}{\text{dt}}=\left( 1-\frac{\text{OD}}{\text{O}{{\text{D}}_{m}}} \right)\left( 1-\frac{1}{1+{{v}_{0}}\exp qt} \right){{\mu }_{m}}\quad(3)
The parameters of equation (3) cannot be fitted by simple regression, so we develop another approach. For any given parameters, we use Runge-Kutta method [7] to calculate how OD evolve with time after the initial value, and use Particle Swarm Optimization(PSO) algorithm [8] to find the suitable parameters that minimize the standard deviation between evolution and measurement.
The fitting results are shown in Fig. 2. Dots are experimental values, and red curve is the fitted results. In the first three to four hours, the physiological factor $\left( 1-\frac{1}{1+{{v}_{0}}\exp qt} \right)$ limits the growth rate, as is shown in green curve. The later stage is governed by the environment factor $\left( 1-\frac{\text{OD}}{\text{O}{{\text{D}}_{m}}} \right)$ shown in brown.
Fig. 2
Experimental values and regression curve. Dots are experimental values, and regression results are shown in red curve. Brown line shows the relative value of environment factor $\left( 1-\frac{\text{OD}}{\text{O}{{\text{D}}_{m}}} \right)$ and green line show the physiological factor $\left( 1-\frac{1}{1+{{v}_{0}}\exp qt} \right)$ respectively.
The fitted values for parameters are shown in the table 1. Since the bacteria in our experiments are of same species and are prepared by a same procedure, we regard the initial value of v as a constant in all these experiments. Other parameters ${\rm{O}}{{\rm{D}}_{\rm{m}}}$, ${{\rm{\mu }}_{\rm{m}}}$ and q are influenced by experimental conditions.
| Symbol | Description | Unit | Value |
|---|---|---|---|
| {\rm{O}}{{\rm{D}}_{\rm{m}}} | Maximum OD | [--] | 1.166 |
| {{\rm{\mu }}_{\rm{m}}} | Maximum growth rate | [1/min] | 0.0155 |
| q | Time derivative of v at free conditions | [1/min] | 0.0213 |
| {{\rm{v}}_0} | Initial value for v | [--] | 0.0282 |
| N | Overall number of bacteria in our solution | [cell] | Unnecessary |
| {{\rm{N}}_{\rm{m}}} | Environmental capacity of bacteria | [cell] | Unnecessary |
| OD | Optical density at 600 nm | [--] | Measured |
| v | Parameter of physiological state of the cells | [--] | Calculated |
Symbols, descriptions and values in sub model of free growth.
Effect of Toxin
Toxin and antitoxin (TA) modules frequently coexist in bacterial genomes. Toxins may cause growth inhibition or death by targeting at DNA replication, mRNA stability, protein synthesis, cell-wall biosynthesis, and ATP synthesis [9].
We study the effects of toxin alone on cell growth in our experiments. We prepare five kinds of plasmid: without toxin, with toxin 134, toxin 136, toxin 1204 and toxin 6249, all of which have RGP_Ptec promoter and are therefore regulated by tetracycline antibiotics. These plasmids are introduced into bacteria cells. We respectively measured the bacteria growth under different concentration of anhydrotetracycline(aTc), which induce the expression of toxin.
Fig. 3
Bacteria growth with toxins. Toxins are under regulation of aTc. (a)OD measured in experiments. (b) Growth rate
Fig. 3 shows the growth of bacteria with toxin 1204 plasmid. The big gap between $0.8{\rm{\;ug}}/{\rm{ml}}$ and $1.0{\rm{\;ug}}/{\rm{ml}}$ can be noticed at the first glance, and this switch phenomenon is analyzed later. We fit the data with methods that mentioned previously, toquantify how parameters in growth function change with toxin. The results are shown in Fig. 4
Fig. 4
Effect of toxin on the parameters in growth equation. Parameter ${{\rm{N}}_{\rm{m}}},{{\rm{\mu }}_{\rm{m}}}\;{\rm{and}}\;q$ of bacteria with different toxin plasmids in different concentration of aTc are estimated by Runge-Kutta and PSO methods. GFP stands for the control group.
The switch phenomenon is merely a speculation so far; we need to carry out scientific examinations. The problem is equivalent to test the hypothesis that there is significant difference in the values of parameters ${{\rm{N}}_{\rm{m}}},{\mu _m},q$ between experimental groups and control groups. Since the measurements are not accurate, we think that a maximum difference of 5% is appropriate. The hypothesis we are going to test is as follows:
{{\rm{H}}_0}:\left| {{{\rm{N}}_{\rm{m}}}\left( {{\rm{toxin}}} \right) - {{\rm{N}}_{\rm{m}}}\left( {{\rm{RGP}}} \right)} \right| < 5{\rm{\% }}{{\rm{N}}_{\rm{m}}}\left( {{\rm{RGP}}} \right)
{{\rm{H}}_0}':\left| {{{\rm{\mu }}_{\rm{m}}}\left( {{\rm{toxin}}} \right) - {{\rm{\mu }}_{\rm{m}}}\left( {{\rm{RGP}}} \right)} \right| < 5{\rm{\% }}{{\rm{\mu }}_{\rm{m}}}\left( {{\rm{RGP}}} \right)
{{\rm{H}}_0}'':\left| {{\rm{q}}\left( {{\rm{toxin}}} \right) - {\rm{q}}\left( {{\rm{RGP}}} \right)} \right| < 5{\rm{\% q}}\left( {{\rm{RGP}}} \right)
We use two-sample heteroscedasticity Student’s t-test method [10] to test these hypothesis, the results are shown in Fig. 5
Fig. 5
Test result of the hypothesis for each type of bacteria in different concentration of aTc. White blocks represent that the hypothesis cannot be rejected at the confidence of 90%, and colored blocks represent rejection at the same confidence.
From Fig. 5, it can be convinced that the experimental group is the same as the control group under low aTc concentrations at the confidence of 90%, and distinct from control group under high aTc concentrations. At a concentration of around $0.8\text{ }\!\!~\!\!\text{ }\!\!\mu\!\!\text{ g}/\text{ml}$ to $1..5\text{ }\!\!~\!\!\text{ }\!\!\mu\!\!\text{ g}/\text{ml}$, a rough boundary can be seen, which confirmed that a switch effect occurs. The in-depth analysis of Tc’s switch effect will be discussed in section 3.
As for the quantitative descriptions, there is significant error in the calculated values of each parameter, so we can only approximately think that, the $\text{O}{{\text{D}}_{\text{m}}}$ as well as$~q$ are reduced by half and ${{\text{ }\!\!\mu\!\!\text{ }}_{\text{m}}}$ decay in exponential form when toxins are expressed:
{\rm{O}}{{\rm{D}}_{\rm{m}}} = \left\{ {\begin{array}{*{20}{c}} {O{D_{m0}},}&{c \le {c_0}}\\ {1/2O{D_{m0}},}&{c > {c_0}} \end{array}} \right.
{\rm{q}} = \left\{ {\begin{array}{*{20}{c}} {{q_0},}&{c \le {c_0}}\\ {1/2{q_0},}&{c > {c_0}} \end{array}} \right.\quad\quad\quad\quad\quad\quad\quad(4)
{{\rm{\mu }}_{\rm{m}}} = \left\{ {\begin{array}{*{20}{c}} {{\mu _{m0}},}&{c \le {c_0}}\\ {{\mu _{m0}}\exp {\alpha _i}\left( {c - {c_0}} \right),}&{c > {c_0}} \end{array}} \right.
The estimated values of parameters are shown in table 2. As can be seen, bacteria is most sensitive with toxin 1204. Further researches are needed to test these conjectures.
| Symbol | Description | Unit | Value | |
|---|---|---|---|---|
| {\rm{O}}{{\rm{D}}_{{\rm{m}}0}} | Maximum OD under toxin-free conditions | [--] | 1.121 | |
| {{\rm{q}}_0} | Time derivative of v under toxin-free conditions | [cell/min] | 0.0198 | |
| {{\rm{\mu }}_{{\rm{m}}0}} | Maximum growth rate under toxin-free conditions | [1/min] | 0.0146 | |
| {{\rm{\alpha }}_{\rm{i}}} | Influence factor of toxin with aTc inducer | Toxin 134 | [ml/ug] | 1.32 |
| Toxin 136 | [ml/ug] | 0.67 | ||
| Toxin 1204 | [ml/ug] | 1.83 | ||
| Toxin 6249 | [ml/ug] | 1.45 | ||
| c | Concentration of inducer aTc | [ml/ug] | experimental | |
Symbols, descriptions and values in sub model of growth with toxins alone.
Toxin-antitoxin Coupling
There are three types of toxin-antitoxin systems known by far, among which only type II is considered in our experiments and discussion. In type II systems, toxins and antitoxins are proteins encoded by small neighboring genes on a same operon[11]. The general reaction is:
\ \\ mRNA_{toxin}\xrightarrow{k_{pt}}Toxin\\ Toxin\xrightarrow{k_{dt}}\phi\\ mRNA_{anti}\xrightarrow{k_{pa}}Antitoxin\\ Antitoxin\xrightarrow{k_{da}}\phi\\ mToxin+nAntitoxin\xrightarrow{k_1}TAcomplex\\ mToxin+nAntitoxin\xleftarrow{k_{-1}}TAcomplex\\ Toxin+Target\xrightarrow{k_{2}}Disfunction\\ Toxin+Target\xleftarrow{k_{-2}}Disfunction
These reactions can be modelled in a set of ordinary differential equations as followings:
\frac{{{\rm{d}}\left[ {{\rm{Toxin}}} \right]}}{{{\rm{dt}}}} = {k_{pt}}\left[ {mRN{A_{toxin}}} \right] - {k_{dt}}\left[ {Toxin} \right] - {k_1}\left[ {Anti} \right]\left[ {Toxin} \right] + {k_{ - 1}}\left[ {TAComplex} \right] - {k_2}\left[ {{\rm{Disfunction}}} \right]\left[ {Toxin} \right] + {k_{ - 2}}\left[ {{\rm{Disfunction}}} \right]
\frac{{{\rm{d}}\left[ {{\rm{Anti}}} \right]}}{{{\rm{dt}}}} = {k_{pa}}\left[ {mRN{A_{anti}}} \right] - {k_{dt}}\left[ {Anti} \right] - {k_1}\left[ {Anti} \right]\left[ {Toxin} \right] + {k_{ - 1}}\left[ {TAComplex} \right] - {k_2}\left[ {{\rm{Disfunction}}} \right]\left[ {Toxin} \right] + {k_{ - 2}}\left[ {{\rm{Disfunction}}} \right]
\frac{{{\rm{d}}\left[ {{\rm{TAComplex}}} \right]}}{{{\rm{dt}}}} = {{\rm{k}}_1}\left[ {Anti} \right]\left[ {Toxin} \right] - {{\rm{k}}_{ - 1}}\left[ {TAComplex} \right]
\frac{{{\rm{d}}\left[ {{\rm{Disfunction}}} \right]}}{{{\rm{dt}}}} = {{\rm{k}}_2}\left[ {Target} \right]\left[ {Toxin} \right] - {{\rm{k}}_{ - 2}}\left[ {Disfuntion} \right]
\quad(5)
From reference [12], toxins are relatively stae bases and antitoxins are instable acids; they can produce a stable TA complex by neutralization to block the function of the toxin. In mathematical form, the numerical relation is:
{{\rm{k}}_{{\rm{pt}}}} \ll {{\rm{k}}_{{\rm{pa}}}}
{{\rm{k}}_{dt}} \ll {k_{da}}
{{\rm{k}}_1} \gg {k_{ - 1}}
{{\rm{k}}_2} \gg {k_{pt}}
So, we can assume that the existing toxin molecules are either in binding with target molecules, or in the toxin-antitoxin complex. The antitoxin is produced to neutralize toxin and the excessive part is degraded. The accurate solution of equation set (5) is not necessary because the expression of toxin only affect the terminal concentration of tetracycline degradation.
We do a series of experiments to understand the coupling of toxin and antitoxins. Bacteria with genes of aTc induced toxin 134, 136, 1204, 6249 and their corresponding IPTG induced antitoxin are prepared.
Firstly, we add aTc in gradient concentration and excessive IPTG(8mmol/L) at the same initial time and measure their growth, as is shown in Fig. 6.
Fig. 6
Growth rate of bacteria with TA system in aTc of gradient concentration and excessive IPTG.
It can be clearly seen from Fig. 6 that, the growth rate curve of toxin 134, 1204, 6249 are of the same shape and same height as the control group, but with a delay. In addition, the growth rate curve of toxin 136 shows that the performance of antitoxin is not as good as others, thus we only focus on toxin 134,1204 and 6429. The maximums of growth rate is slightly affected by the toxin of aTc, but it will do not change the discussion above.
To further understand the process of neutralization, we analyze the variation of physiological state v. Parameter v can be solved from equation (3):
{\rm{v}} = \frac{{\mu /{\mu _m}}}{{1 - OD/O{D_m} - \mu /{\mu _m}}}\quad(6)
${{\text{ }\!\!\mu\!\!\text{ }}_{\text{m}}}~$and$~O{{D}_{m}}$ is fitted by the control group, as is shown in table 3.
| Symbol | Description | Unit | Value | |
|---|---|---|---|---|
| {\rm{O}}{{\rm{D}}_{\rm{m}}} | Maximum OD | {\rm{aTc\;}}5.0{\rm{\mu g}}/{\rm{mL}} | [--] | 1.110 |
| {\rm{aTc\;}}2.5{\rm{\mu g}}/{\rm{mL}} | 1.111 | |||
| {\rm{aTc\;}}1.7{\rm{\mu g}}/{\rm{mL}} | 1.116 | |||
| {{\rm{\mu }}_{\rm{m}}} | Maximum growth rate | {\rm{aTc\;}}5.0{\rm{\mu g}}/{\rm{mL}} | [1/min] | 0.0131 |
| {\rm{aTc\;}}2.5{\rm{\mu g}}/{\rm{mL}} | 0.0171 | |||
| {\rm{aTc\;}}1.7{\rm{\mu g}}/{\rm{mL}} | 0.0194 | |||
Fitted values of $\text{O}{{\text{D}}_{\text{m}}}$ and ${{\text{ }\!\!\mu\!\!\text{ }}_{\text{m}}}$.
The variation of v is shown in Fig. 7.
Fig. 7
Variation of parameter q in Fig. 6
As expected, the shape and maximum are the same within error limits, and a delay can be clearly seen. Delay increase with concentration. The reason for the delay is that, the way across the cytomembrane is different for aTc and IPTG. ATc enter the cell by passive diffusion, and IPTG by active transportation, Therefore the expression of toxin is faster. During the delay, antitoxin are synthesized and form binding with toxin. As long as the toxin is neutralized by its antitoxin, bacteria return to the normal state and grow as if in free conditions. This can be described in following forms:
{\rm{v}} = \left\{ {\begin{array}{*{20}{c}} {{v_0},}&{t \le {t_0}}\\ {{v_0}\exp qt,}&{t > {t_0}} \end{array}} \right.\quad(7)
To calculate ${{\text{t}}_{0}}$, we use correlation coefficient to measure the similarity of growth rate of control group ${{\text{v}}_{0}}\left( \text{t} \right)$ and experiment group $\text{v}\left( \text{t}+{{\text{t}}_{0}} \right)$, then ${{\text{t}}_{0}}$ should be the interval that maximum the correlation coefficient:
\mathop {\max }\limits_{{t_0}} \left| {\frac{{E\left( {\left( {v\left( {t + {t_0}} \right) - Ev} \right)\left( {{v_0}\left( t \right) - E{v_0}} \right)} \right)}}{{\sqrt {D\left( {v\left( t \right)} \right)D\left( {{v_0}\left( t \right)} \right)} }}} \right|\quad(8)
The calculated ${{\text{t}}_{0}}$ is shown in table 4:
| {\rm{aTc\;\mu g}}/{\rm{mL}} | Toxin 134 /min | Toxin 1204 /min | Toxin 6249 /min |
|---|---|---|---|
| 5.0 | 170 | 130 | 100 |
| 2.5 | 50 | 50 | 50 |
| 1.7 | 30 | 50 | 30 |
Delay ${{\rm{t}}_0}$ for bacteria with different TA system at different concentration of inducer.
At a same inducer concentration, a shorter delay ${{\rm{t}}_0}$ indicates a higher efficiency of antitoxin. From table 3, we can conclude that the corresponding antitoxin of toxin 6249 is most efficient.
Secondly, we add aTc at a concentration of $1.7\text{ }\!\!\mu\!\!\text{ g}/\text{mL}$ and excessive IPTG(8mmol/L) at different time and measure their growth, as is shown in Fig. 8.
Fig. 8
Growth rate of bacteria with TA system induced by aTc and excessive IPTG at different time.
The delay is computed with the same method mentioned above, as is shown in table 5.
| Conditions | Toxin 134 /min | Toxin 1204 /min | Toxin 6249 /min |
|---|---|---|---|
| IPTG 4h after | 140 | 360 | 510 |
| Simultaneous | 40 | 50 | 30 |
| IPTG 4h before | 10 | 20 | 10 |
Delay ${{\rm{t}}_0}$ when inducer is added at different time.
The result shows that a much longer time is needed for bacteria to overcome the effect of toxin, if antitoxin is induced later than toxin. On the contrary, bacteria are not influenced by the toxin inducer when antitoxins are expressed earlier.
From another perspective, we can see that toxin 134 is inactivated before the antitoxin is induced by IPTG, which is because the leakage of antitoxin promoter or the instability of toxin. It takes the longest time for bacteria cells to overcome the effect of toxin 6249 after antitoxin is induced, therefore toxin 6249 is most toxic and also stable.
From the analysis mentioned above, we can conclude that toxin 6249 is most toxic and stable, its antitoxin is most efficient in neutralization as well. Thus, toxin 6249 is our final choice. The demonstration above shows how to select an appropriate TA block with the mathematical model.
Data processing are performed with MATLAB 2014a. Most figures are plotted with MATLAB 2014a and R 3.3.1. Download a complete collection of all files here.
Protein Expression
Dominant Sectors
Proteins are expressed through transcription, translation process. Transcription involves RNA polymerase binding, RNA formation; translation involves the binding of ribosome complex to RBS site, elongation and termination. The bindings of first two processes are the limiting and rate-defining step.All of these substances can be degraded, therefore each has a degradation rate. It can be seen how many parameters there are!
Models usually relay on searching the parameters from previous works, but the values are measured under different conditions by different methods so the errors cannot be simply ignored. To make things worse, the ordinary differential equations are sensitive to parameters, and the errors will be amplified with time.
To exemplify how error can be amplified, we would like to use model and parameter values from previous iGEM competitions. We appreciate the work done by 2015 TU Delft iGEM Team [13]. They measured the fluorescence which is similar to our methods, therefore we choose their results for illustration.
We add a normal-distributed random relative error with a maximum of 5% to each of the 9 parameters in their model, which is far within the significant digit of these parameters. The procedure is repeated for 100 times, with the results of which we estimate how error in parameters will influence the predictions. Figure 9 is the calculation result, and the red area is the 70% confidence interval. The mean relative error is around 10%.
Fig. 9
Prediction of GFP and errors of 2015 TU Delft iGEM Team’s model by adding random relative error with a maximum of 5%. The green curve is the mean value, and red area is the range mean value added and subtracted by the standard error.
The example above clearly demonstrates the fact that, ODE model is sensitive to the values of parameters. Therefore, a model with as many as 9 parameters collected from different sources are not only unnecessary but also unreliable.
Can there be a model that fulfill the same functions as the traditional models, but with less parameters and with the least experiments? Actually, it is possible.
Generally, the processes of protein can be precisely modeled as equation (9).
\ \\ \frac{dM}{dt} = PoPs(t)·n(t)-\gamma_M(t)M(t)\\ \frac{dI}{dt}=RiPS(t)·M(t)-(\gamma_P(t)+k_m(t))I(t)\quad\quad\quad (9)\\ \frac{dP}{dt}=k_m(t)I(t)-\gamma_P(t)P(t)
| Symbol | Description | Unit |
|---|---|---|
| M | Total number of mRNA | [--] |
| I | Total number of immature peptide chain | [--] |
| P | Total number of protein molecules | [--] |
| n(t) | Total number of plasmid in all cells | [--] |
| PoPS | Number of transcription polymerases per second | \left[ {\frac{{Poly.}}{{plasmid \times s}}} \right] |
| RiPS | Number of translation Ribosomes per second | \left[ {\frac{{Ribosome}}{{mRNA \times s}}} \right] |
| {k_m} | Mature rate of peptide chain | \left[ {\frac{{{\rm{Protein}}}}{{{\rm{chain}} \times {\rm{s}}}}} \right] |
| {\gamma _P} | Degradation rate of peptide chain and protein | \left[ {{{\min }^{ - 1}}} \right] |
Symbols, descriptions and units in equation set (9).
Further analysis is unpractical when all these parameters are time-verifying, so we obtain the solution of equation set(9) at first, and reach the strict solution thereafter.
Lee and Bailey established a detailed model of plasmid replication [15]. The form of the model is too complex, but the result shows that the replication period is converged to be synchronized with cell cycle. Accordingly, we can assume that the growth rate of total number of plasmids is equal to that of cells:
\frac{{dn}}{{dt}} = \mu \cdot n\quad(10)
Firstly, we discuss the solution of equation set (9) when μ is constant. The equations can thus be rewritten in a matrix form:
\frac{d}{{dt}}\left( {\begin{array}{*{20}{c}} n\\ M\\ I\\ P \end{array}} \right) = \left( {\begin{array}{*{20}{c}} \mu &0&0&0\\ {PoPS}&{ - {\gamma _M}}&0&0\\ 0&{RiPS}&{ - ({\gamma _P} + {k_m})}&0\\ 0&0&{{k_m}}&{ - {\gamma _P}} \end{array}} \right)\left( {\begin{array}{*{20}{c}} n\\ M\\ I\\ P \end{array}} \right)\quad(11)
Let ${\rm{x}}$ denote the vector ${\rm{x}}$. For any given initial value ${x_0}$, there exists a unique eigen decomposition of ${x_0}$:
{{\bf{x}}_0} = \sum\limits_i^{} {{a_i}{{\bf{v}}_i}} \quad(12)
The ${\alpha _{\rm{i}}}$,i=1,2,3,4 are coefficients that are determined by the initial conditions. In equation (12), the eigen vectors ${v_i}$ and their eigen values ${\lambda _i}$ are listed below:
| Serial number i | Eigen Vector {{\rm{v}}_i} | Eigen Value {\lambda _i} |
|---|---|---|
| 1 | \left( {\begin{array}{*{20}{c}} {({\gamma _M} + \mu )({\gamma _P} + \mu )({\gamma _P} + {k_m} + \mu )}\\ {PoPS({\gamma _P} + \mu )({\gamma _P} + {k_m} + \mu )}\\ {PoPS \cdot RiPS \cdot ({\gamma _P} + \mu )}\\ {PoPS \cdot RiPS \cdot {k_m}} \end{array}} \right) | \mu |
| 2 | \left( {\begin{array}{*{20}{c}} 0\\ {({\gamma _P} - {\gamma _M})({\gamma _P} - {\gamma _M} + {k_m})}\\ {({\gamma _P} - {\gamma _M})RiPS}\\ {RiPS \cdot {k_m}} \end{array}} \right) | - {\gamma _M} |
| 3 | {\left( {\begin{array}{*{20}{c}} 0&0&{ - ({\gamma _P} - {\gamma _M} + {k_m})}&{{k_m}} \end{array}} \right)^T} | - ({\gamma _P} + {k_m}) |
| 4 | {\left( {\begin{array}{*{20}{c}} 0&0&0&1 \end{array}} \right)^T} | - {\gamma _P} |
Eigen decomposition of equation (11).
Then, the solution of equation (6) is:
{\bf{x}}(t) = \sum\limits_i^{} {{a_i}{{\bf{v}}_i}\exp {\lambda _i}t} \quad(13)
By a brief estimation, we can have a concept of the contribution proportion of each term in equation (13). The typical values of parameters are selected as shown in table 8 only for an estimation. The contribution proportion is independent of the initial value of n[13].
| Parameter | Unit | Value | Reference |
|---|---|---|---|
| μ | \left[ {{{\min }^{ - 1}}} \right] | 6 \times {10^{ - 3}} | Experiment |
| PoPS | \left[ {\frac{{Poly.}}{{plasmid \times s}}} \right] | 5 \times {10^{ - 4}} | [16] |
| {\gamma _M} | \left[ {{{\min }^{ - 1}}} \right] | 0.29 | [16] |
| RiPS | \left[ {\frac{{Ribosome}}{{mRNA \times s}}} \right] | 3 \times {10^{ - 4}} | [17] |
| {\gamma _P} | \left[ {{{\min }^{ - 1}}} \right] | 4 \times {10^{ - 3}} | [18] |
| {k_M} | \left[ {{{\min }^{ - 1}}} \right] | 0.04 | [13] |
| n | [--] | Unneccessary | - |
Typical values of parameters in equation (9).
Fig. 10
Each colored area shows the contribution proportion of the sum of absolute value of each term in equation (13). The initial value of x is chosen to be ${{\rm{x}}_0} = {(\begin{array}{*{20}{c}}1&0&0&0\end{array})^T}$
The dominant items can be clearly seen from figure 10. The term 1 and term 4 of equation (12) take the vast majority of the sum of absolute values. Term 2 is very little and term 3 die out as time goes on. Thus, we only retain term 1 and term 4 in our equation.
x(t) = \alpha {{\bf{v}}_1}\exp \mu t + \beta {{\bf{v}}_4}\exp - {\gamma _P}t\quad(14)
Term 1 is the equilibrium sector of protein expression, at which state the production rate and degradation rate of mRNA, peptide chain and protein are equal. Term 4 is the degradation sector. It cannot be ignored because the lifetime of protein is rather stable that mRNA and peptide chain. Now, the two dominant sectors are extracted.
General Solution
Note that equations (11) to (14) hold only under constant growth rate condition. It can be expected that equilibrium sector and the degradation sector will still be the dominant ones when parameters are time-dependent.
\left( {\begin{array}{*{20}{c}} {{\rm{n}}(t)}\\ {P(t)} \end{array}} \right) = \alpha (t)\left( {\begin{array}{*{20}{c}} {({\gamma _M} + \mu (t))({\gamma _P} + \mu (t))({\gamma _P} + {k_m} + \mu (t))}\\ {PoPS(t) \cdot RiPS(t) \cdot {k_m}} \end{array}} \right) + \beta (t)\left( {\begin{array}{*{20}{c}} 0\\ 1 \end{array}} \right)\quad(15)
$\alpha (t)$ and $\beta (t)$ are unique for any given time t, because the two vectors are linear independent. From equation (13) and (15), we can express the situation at time t+dt in two ways:
\left( {\begin{array}{*{20}{c}} {{\rm{n}}(t + dt)}\\ {P(t + dt)} \end{array}} \right) = \alpha (t + dt)\left( {\begin{array}{*{20}{c}} {({\gamma _M} + \mu (t + dt))({\gamma _P} + \mu (t + dt))({\gamma _P} + {k_m} + \mu (t + dt))}\\ {PoPS(t + dt) \cdot RiPS(t + dt) \cdot {k_m}} \end{array}} \right) + \beta (t + dt)\left( {\begin{array}{*{20}{c}} 0\\ 1 \end{array}} \right)
= \alpha (t)\left( {\begin{array}{*{20}{c}} {({\gamma _M} + \mu (t))({\gamma _P} + \mu (t))({\gamma _P} + {k_m} + \mu (t))}\\ {PoPS(t) \cdot RiPS(t) \cdot {k_m}} \end{array}} \right)\exp \mu (t)dt + \beta (t)\left( {\begin{array}{*{20}{c}} 0\\ 1 \end{array}} \right)\exp ( - {\gamma _P}dt)\quad(16)
So, the following equation must hold:
\begin{array}{l} \alpha (t + dt) = \frac{{({\gamma _M} + \mu (t))({\gamma _P} + \mu (t))({\gamma _P} + {k_m} + \mu (t))}}{{({\gamma _M} + \mu (t + dt))({\gamma _P} + \mu (t + dt))({\gamma _P} + {k_m} + \mu (t + dt))}}\alpha (t)\exp \mu (t)dt\\ = \frac{{\alpha (t)\exp \mu (t)dt}}{{\left( {1 + \frac{{d\mu }}{{{\gamma _M} + \mu (t)}}} \right)\left( {1 + \frac{{d\mu }}{{{\gamma _P} + \mu (t)}}} \right)\left( {1 + \frac{{d\mu }}{{{\gamma _P} + {k_m} + \mu (t)}}} \right)}} \end{array}\quad(17)
\begin{array}{l} d(\ln \alpha (t)) = \mu (t)dt - \ln \left( {1 + \frac{{d\mu }}{{{\gamma _M} + \mu (t)}}} \right) - \ln \left( {1 + \frac{{d\mu }}{{{\gamma _P} + \mu (t)}}} \right) - \ln \left( {1 + \frac{{d\mu }}{{{\gamma _P} + {k_m} + \mu (t)}}} \right)\\ = \left( {\mu (t) - \left( {\frac{1}{{{\gamma _M} + \mu (t)}} + \frac{1}{{{\gamma _P} + \mu (t)}} + \frac{1}{{{\gamma _P} + {k_m} + \mu (t)}}} \right)\frac{{d\mu }}{{dt}}} \right)dt\\ = \left( {\mu (t) - \left( {\frac{1}{{{\gamma _M} + \mu (t)}} + \frac{1}{{{\gamma _P} + \mu (t)}} + \frac{1}{{{\gamma _P} + {k_m} + \mu (t)}}} \right)\frac{{d\mu }}{{dt}}} \right)dt \end{array}\quad(18)
From table 8, the relative magnitude among parameters is:
\frac{{d\mu }}{{d{\rm{s}}}} < \mu \approx {\gamma _P} < {k_{\rm{m}}} < {\gamma _M}\quad(19)
Therefore, only the third fraction is reserved. The solution for parameter is complicated, but it is very close to an exponential decay. Then the quantity of protein can be expressed as:
P(t) = E\exp \int_0^t {\left( {\mu (s) - \frac{1}{{{\gamma _P} + \mu (s)}}\frac{{d\mu }}{{ds}}} \right)} ds + D\exp - {\gamma _P}t\quad(20)
If growth rate is really low compared to degradation rate, and no protein molecule exists at time 0, then 0=P(0)=E+D. Equation (20) can be further simplified as:
P(t) = {P_0}\left( {\exp \int_0^t {\mu (s)ds} - \exp - {\gamma _P}t} \right)\quad(21)
Equation (20) and (21) are called Protein Expression Formula (PEF).
Induced Protein Expression
In this section, we compare Protein Expression Formula with experimental data, and discuss how parameter ${P_0}$ verifies with the concentration of tetracycline.
From table 7 and equation (21), ${P_0}$ can be solve as:
{P_0} = \frac{{PoPS \cdot RiPS \cdot {k_m} \cdot {n_0}}}{{({\gamma _M} + {\mu _0})({\gamma _P} + {\mu _0})({\gamma _P} + {k_m} + {\mu _0})}}\quad(22)
The fluorescence of GFP is in linear relation with P(t). Parameter ${n_0}$, the total number of plasmids, is in proportion with the number of cells and hence in proportion with OD600. The remaining terms in equation (22) is a function of concentration. Due to the limited accuracy of instruments, derivative of growth rate $\frac{{d\mu }}{{ds}}$ cannot be calculated. Consequently, the experimental equation of formula (20) can be expressed as:
\frac{{F(t)}}{{O{D_0}}} = B(c)\exp \int\limits_0^t {\mu (s)ds} - a\exp ( - {\gamma _P}t) + b\quad(23)
B(c) in equation (23) is as follows:
B(c) = \frac{{PoPS(c) \cdot RiPS(c) \cdot {k_m}}}{{({\gamma _M} + {\mu _0})({\gamma _P} + {\mu _0})({\gamma _P} + {k_m} + {\mu _0})}}\quad(24)
Coefficient a and b are decided by the initial state and is unimportant in our discussion. For any given ${\gamma _P}$ , equation (23) can be fitted by multiple linear regression. Then we assign ${\gamma _P}$ point by point with a step of 0.000001 and find the best ${\gamma _P}$ that maximize the linear coefficient.
In experiments, plasmid with tetR-ptet-tetX-GFP sequence is imported to E.coli. The regression results of equation (18) turn out to be excellent. The fitted parameters are shown in table 9. B(c) curve is shown in Figure 11.
| {\rm{Tc}}{\rm{. Conc}}{\rm{. /}}\mu g \cdot m{L^{ - 1}} | B(c) \cdot {10^{ - 4}} | {\gamma _P}/{\min ^{ - 1}} | {R^2} |
|---|---|---|---|
| 40 | 2.8752 | 0.0042 | 0.9858 |
| 20 | 3.0742 | 0.0039 | 0.9863 |
| 10 | 3.3639 | 0.0050 | 0.9942 |
| 5 | 3.7352 | 0.0043 | 0.9855 |
| 2.5 | 5.0076 | 0.0048 | 0.9877 |
| 1.25 | 4.5055 | 0.0050 | 0.9147 |
| 1.0 | 1.5259 | 0.0058 | 0.9780 |
| 0.5 | 0.1379 | 0.0060 | 0.9593 |
| 0 | 0.1236 | 0.0063 | 0.9671 |
Regression results of tetR-ptet-tetX-GFP E.coli with equation (23).
Fig. 11
Parameter B(c) under different concentrations.
As is shown in figure 11, B(c) is low but not zeros when the concentration of tetracycline is zeros, this phenomenon indicates that leakage of promoter happens. The protein expression remains depressed when concentration is below 0.8 mg/L, which is another demonstration of switch effect discussed in part 1. Later, B(c) increase rapidly and reach its peak at around 2.5 mg/L, at which concentration the promoter is activated to the most. After that, excessive tetracycline molecules bind with mRNA and reduce the RiPS, thus B(c) decreases in an exponential pattern.
B(c) = \left\{ {\begin{array}{*{20}{c}} {{B_0},}&{c \le {c_0}}\\ {{{\rm{B}}_m}\frac{{c - {c_0}}}{{{c_m} - {c_0}}},}&{{c_0} < c \le {c_m}}\\ {{{ \rm{B}}_m}\left( {b\exp \lambda (c - {c_m}) + 1 - b} \right),}&{c > {c_m}} \end{array}} \right.\quad(25)
The degradation rate ${\gamma _P}$ is not much influenced by the concentration. It is worth mentioning that ${\gamma _P}$ from regression is close to that measured in other methods. Symbols, descriptions and values are listed in table 10.
| Symbol | Symbol | Unit | Value |
|---|---|---|---|
| {B_0} | B when Tc concentration equals to zeros | \left[ {\frac{{{\rm{protein }}}}{{{\rm{plasmid}}}}} \right] | 1.236 \times {10^3} |
| {B_m} | Maximum value of B | 5.008 \times {10^4} | |
| {c_0} | Switch threshold concentration of tetracycline | [mg/L] | 0.8 |
| {c_m} | Tc concentration at the peak of B | [mg/L] | 2.5 |
| b | Parameter in equation (25) | [--] | 0.54 |
| \lambda | Parameter in equation (25) | [L/mg] | 0.069 |
| {\gamma _P} | Protein degradation rate fitted from data | [{\min ^{ - 1}}] | 0.005 |
Symbols, descriptions and values in sub model of free growth.
Overall Circuit
Diffusion and Degradation Rate
There are only two details needed to be discussed : diffusion and degradation rate.
In the tetracycline uptake process of E.coli, two systems are involved: the rapid concentration-driven system which occurs initially and is insensitive to metabolic inhibitors, and the slow energy-driven system whose rate is only one-tenth of the rapid system[19]. The latter is not as important as the former, hence we only consider the concentration-driven factors. Though it is still unknown whether the concentration-driven system can be categorized as passive diffusion or facilitated transport, the transportation rate is in proportion to the gradient of tetracycline concentration. So, we just take it as diffusion process.
Diffusion through cell membranes can be described by permeability coefficient, which is an important is defined by the following formula:
Papp = \frac{Q}{{A \times {\rm{ }}\Delta {\rm{C}} \times T}}\quad(26)
The parameters in equation (26) are listed as below.
| Symbol | Description | Unit | Value | Reference |
|---|---|---|---|---|
| Papp | Permeability | [cm/s] | 0.8\times10^{-6} | [20] |
| A_{0} | Mean area of E.coli cytomembrane | [m^{2}] | 3.9\times10^{-12} | [20] |
| V_{0} | Mean volume of an E.coli cell | [m^{3}] | 5.6\times10^{-19} | [21] |
| \Delta C | Concentration difference across membrane | [mol/L] | -- | -- |
Values of parameters in equation (26).
Then the differential of tetracycline concentration within the cell caused by diffusion can be calculated directly. $G=Papp\times \frac{A_{0}}{V_{0}}=0.67 s^{-1}$.
\frac{{d[Tet]}}{{dt}} = Papp \times \frac{{{A_0}}}{{{V_0}}} \times \Delta C = G\Delta C\quad(27)
Protein degradation can be described by Hill equation:
\theta = {k_{cat}}P\frac{{{{[Tet]}^n}}}{{{K_d} + {{[Tet]}^n}}}\quad(28)
We measured the parameters as follows:
Where P is the number of degradation proteins, $k_{cat}$ is the turn over number of protein.
Every modules in our engineering bacteria have been discussed right now, the last thing to do is to put them together.
| Symbol | Description | Unit | Value |
|---|---|---|---|
| P | Number of degradation proteins | [--] | -- |
| k_{cat} | Turnover number of protein | [s^{-1}] | 0.21 |
| n | Hill coefficient | [--] | 1 |
| k_{d} | Apparent dissociation constant | [\mu M] | 48.9 |
| \Delta C | Concentration difference across membrane | -- | caculated |
Complete Model
The complete model is as below:
{\mu = \frac{1}{N}\frac{{dN}}{{dt}} = \left\{ {\begin{array}{*{20}{c}} {\left( {1 - \frac{N}{{{N_m}}}} \right)\left( {1 - \frac{1}{{1 + {v_0}\exp qt}}} \right){\mu _{m0}}\exp \alpha ([Tet] - {c_0}),}&{\frac{P}{N} \ge thres}\\ {0,}&{\frac{P}{N} < thres} \end{array}} \right.}
{P(t) = B({\rm{C}}){n_0}\exp \int\limits_0^t {\mu (s)ds} - D\exp - {\gamma _P}t}
{\frac{{d[Tet]}}{{dt}} = G\left( {C - [Tet]} \right) - \frac {{{ k_{cat}}[Tet]}}{{{K_d} + [Tet]}}\frac{P}{N}}
| Symbol | Description | Unit | Value |
|---|---|---|---|
| V | Volume of the whole solution | {m^3} | Given |
| C | Tetracycline concentration outside the cells | mg/L | Given |
| thres | Toxin functions when protein in a cell is below this threshold. | [--] | Unknown |
References
[7] Timothy Sauer. Numerical Analysis (Second Edition). New York: Pearson Education, 2012. Chapter 5.
[8] Kennedy J. Particle Swarm Optimization[M]. Springer US, 2011.
