Stability

Living organisms have genetic and biochemical tools to respond to stress. These tools normally up-regulate some stress-specific genes, and, directly or indirectly, down-regulate the rest of normal genes (most of them housekeeping genes). On the other hand, every E. coli component works worse under conditions different from their optimum (for example, nutrients uptake or ATP production), so that, every E. coli function will be reduced. We model the relative expression of biobricks by variations of temperature, pH and salt concentration.

Protein production

The mRNA and protein production of a cell can be described by the kinetic equations:

$$ \begin{eqnarray} \frac{d N_{mRNA} }{dt}&=& P \cdot r_T - N_{mRNA} \cdot \lambda_{mRNA} \nonumber ,\\ \frac{d N_{Pr} }{dt}&=& N_{mRNA} \cdot r_P - N_{Pr} \cdot \lambda_{Pr}, \nonumber \end{eqnarray} $$

where $N_{\it{mRNA}} (N_{Pr} )$ is the number of biobrick mRNA (protein) molecules per cell, $\lambda_{mRNA} $ ($\lambda_{Pr}$) is the mRNA (protein) decay rate, $r_T$ ($r_P$) is the transcription (translation) rate, P is the mean number of plasmids per cell.

Assuming stationary messenger and protein content during the cell cycle, $\frac{d N_{mRNA} }{dt}=0$ and $\frac{d N_{Pr} }{dt}=0$, we get

This equation describes the amount of total synthesized protein as a function of calculable rates and the number of free plasmids. It can be easily shown that, after many generations, the mean number of proteins per cell is 1.5 times the protein production during a cell cycle.

1) TEMPERATURE

The expression level of our heterologous gene can be calculated by equations 1 and 2. Our main goal is to find which of the terms will change most significantly with the different stress factors. We will be interested in the activity relative to a reference one, which we take as the activity of E. coli bacteria grown on the LB medium at 38.7 ^oC, pH 7, 1% of NaCl, no UV irradiation and 1 atm of pressure.

$$ \begin{equation} \frac{A_{c}}{A^0_{c}} = \prod_i \frac{V_i}{V^0_i} . \end{equation} $$ where $V$ is any of the factors contributing to the activity and the index $0$ refers to the reference case.

Transcription rate

Among several factors, we first discuss the transcription rate. In the range of temperatures (or pH) that we have studied, the polymerase is sufficiently stable to assume that the elongation rate is not affected and the interaction energies do not change significantly (for example, the promoter specific interaction do only change 16% from 33 ^oC to 44 ^oC). We did not find any significant difference between the nucleotide compositions of E. coli under stress [2]. Therefore, main variations in the transcription rate are due to the probability of RNAP-promoter bounds $$ \begin{equation} r_T (T)=\frac{1}{\frac{1}{r_I(\Delta\epsilon)}+\frac{1}{r_E}}Prob_B (\Delta\epsilon, T) \propto Prob_B (T). \nonumber \end{equation} $$ On the other hand, the down-regulation effect affects the expression level of our Biobrick, because it works as a housekeeping gene (with a normal -35 and -10 promoter for the RNA polymerase). For example, E. coli has a common σ factor for housekeeping genes (σ70), and, another one (σ32) for temperature stress. In the case of temperature stress, the bacteria will synthetize σ32, which will bind RNA polymerase (RNAP) cores that won’t be free for σ70, and so they will not be able to synthetize mRNA of the normal housekeeping genes.

Therefore, we just need to find a relation between the amount of σ factors and the expression level, which we derive next with thermodynamic arguments, and also the variation of σ factors with temperature, which can be determined experimentally.

The RNAP is bound to a promoter (unspecific DNA) with $\epsilon_{S}$ ($\epsilon_{NS}$) free energy. There are $N_{NS}$ nonspecific sites on the E. coli genome (sequences of 41bp) and R free E$\sigma$70, with $R \ll N_{NS}$ . The combinatory number of unbound states, with energy R$\epsilon_{NS}$, is $$ \begin{equation} \frac{N_{NS}!}{R!(N_{NS}-R)!} \simeq \frac{N^R_{NS}}{R!} . \nonumber \end{equation} $$ Similarly, the combinatory number with one bound state, with energy $\epsilon_{S}$+(R-1)$\epsilon_{NS}$, is $$ \begin{equation} \frac{N_{NS}!}{(R-1)!(N_{NS}-R-1)!} \simeq \frac{N^{R-1}_{NS}}{(R-1)!} . \nonumber \end{equation} $$ Then, the probability of finding a RNAP bound to a promoter is with the enter $$ \begin{equation} Prob_B = \frac {\frac{R}{N_{NS}} e^{-\frac{\Delta\epsilon}{k_BT }}}{1+\frac{R}{N_{NS}} e^{-\frac{\Delta\epsilon}{k_BT }}}, \end{equation} $$ with $\Delta\epsilon = \epsilon_{S} - \epsilon_{NS}$, calculated with the energy matrix developed in [3] and $k_B$ is the Bolztmann constant. R $\simeq$ 1000 and $N_{NS} \simeq 10^7$ for E. coli.
Using the measured variation of free E$\sigma$70, in arbitrary units, with temperature, in $^o$C, [4] $$ \begin{equation} E\sigma70 = 0.027 T^2 - 2.12 T + 42.2, \nonumber \end{equation} $$ we can determine the relative change of transcription rate $$ \begin{equation} \frac{r_T}{r^0_T} = \frac{Prob_B}{Prob^0_B} = \frac { 1+\frac{R^0}{N_{NS}} e^{-\frac{\Delta\epsilon}{k_BT }} } { \frac {R^0}{R}+\frac{R^0}{N_{NS}} e^{-\frac{\Delta\epsilon}{k_BT }}} , \end{equation} $$ with $$ \begin{equation} \frac{R^0}{R} = \frac {0.027 (T^0)^2 - 2.12 T^0 + 42.2} {0.027 T^2 - 2.12 T + 42.2} . \end{equation} $$

Synthesis rate

Modeling synthesis rate follows the lines discussed in transcription rate, but the working regime is quite different. The energies of the mRNA-ribosome sequences interactions are much larger than in the RNAP-promoter ones. The average difference of free energies (computed with the RBS calculator for random sequences of E. coli genes) is $\Delta\epsilon = -20.79 kT$, which means that, in the equation: $$ \begin{equation} \frac{r_T}{r^0_T} = \frac{Prob_B}{Prob^0_B} = \frac { 1+\frac{R^0}{N_{NS}} e^{-\frac{\Delta\epsilon}{k_BT }} } { \frac {R^0}{R}+\frac{R^0}{N_{NS}} e^{-\frac{\Delta\epsilon}{k_BT }}} , \end{equation} $$ The term $\frac{R^0}{N_{NS}} e^{-\frac{\Delta\epsilon}{k_BT }}>>>\frac {R^0}{R}$. With this, and knowing that there is no effect on the ribosomal protein composition with temperature cahnges [5], we conclude that the most signifficant factor should be the elongation rate.

We discuss next the dependence of elongation rates with temperature. We use the relation of amino acid concentration with the temperature [7]. Using data from [8], we know that small variations on the ternary complex does not affects the elongation rate. The slowests steps are: 1) GTP hydrolysis on ternary complex (0.01 s), 2) dipeptide bond formation (0.03s) and 3) translocation (0.022s). $$ \begin{equation} r_E = \frac{1}{\frac{1}{A_1}e^{\frac{E^1_a}{RT}}+\frac{1}{A_2}e^{\frac{E^2_a}{RT}}+\frac{1}{A_3}e^{\frac{E^3_a}{RT}}}, \end{equation} $$ where $E^i_a$ is (32, 75.2, 29.260) kJ/mol and $A^{-1}_i$ is (2.5 $\cdot 10^7$, 1.6 $\cdot 10^{14}$, 3.9 $\cdot 10^6$).

Summary

Our model includes two main contributions to the changes of the enzymatic activity due to temperature variation: transcription and synthesis rates. In the case of the GFP, the folding free energy term is small and the relative enzymatic activity can be simplified to our master equation for stability $$ \begin{equation} \frac{Ac}{Ac^0} = \frac{r_T}{r^0_T} \frac{r_P}{r^0_P} . \end{equation} $$

pH & SALT CONCENTRATION

In contrast with temperature variations, E. coli internal pH and salt concentration hardly vary. Therefore, we focus on the intrinsic E. coli element that may be affected by these stresses, and influence Biobrick expression. The main factor is the number of free RNAP, due to the synthesis of stress response genes. We estimate the ratio of free RNAP under different conditions by the ratio of their growth rates $$ \begin{equation} \frac{r_T}{r^0_T} = \frac{Prob_B}{Prob^0_B} = \frac { 1+\frac{R^0}{N_{NS}} e^{-\frac{\Delta\epsilon}{k_BT }} } { \frac {R^0}{R}+\frac{R^0}{N_{NS}} e^{-\frac{\Delta\epsilon}{k_BT }}} , \end{equation} $$ with $$ \begin{equation} \frac{R^0}{R} = \frac{growth^0} {growth} . \end{equation} $$

Comparison with wet-lab data

Our model captures the measured shape of the fluorescence per cell as a function of temperature. We have found that transcription rate decreases with temperature as a consequence in the reduction of free RNAP; and synthesis rate grows above the optimal temperature. Unfortunately, this is not the whole story. Our model dramatically fails to fit the amount of variation, which is much larger experimentally than the predicted one. Let us compare theoretical predictions and experimental results. In both studies, strain XL1 Blue transformed with Bb1 or Bb3, we find that the ratio between theoretical prediction and experimental result is the same at lower temperatures than the optimal (in fact, accidentally too close given the error bars). This fact suggests that the temperature dependence not included in our model should come from cell mechanisms; independent from parameters of the Biobricks (promoter strength, ...). Quite on the contrary, the comparison of theoretical predictions and experimental results at higher temperatures show that our model has not yet captured the dependence of the magnitude of the decrease in activity with the Biobrick parameters.

Standardization

We discuss now Biobrick expression in different strains of E. coli. Changes on the number of free RNAP affect the mRNA synthesis rate depending on the energy of the promoter. The relative transcription rate in different strains can be written as $$ \begin{equation} \frac{r_T}{r'_T} = \frac{Prob_B}{Prob'_B} = \frac { 1+\frac{R'}{N_{NS}} e^{-\frac{\Delta\epsilon}{k_BT }} } { \frac {R'}{R}+\frac{R'}{N_{NS}} e^{-\frac{\Delta\epsilon}{k_BT }}} , \end{equation} $$ where we can see that the dominance of the second term over the first one, controlled by the promoter strength, would lead to similar rate independently on the particular strain chosen. Comparing this formula with the relative transcription rate in stability, we find that promoters whose binding energy makes them more standard also implies more stable promoters. Strong promoter sequences experience a decrease in synthesis rates when subjected to an increase of promoter energy. We use these results and the energy matrix derived in [3] to create an application that is able to compute promoter sequences for a given relative (to the maximum) expression rate and associated standard and stable rate. The physical parameters behind are the promoter strength (using the data of Anderson’ promoters) and the binding energy and the tool create a sequence of a constitutive promoter with the appropriated trade-off between strength and energy, therefore, with the appropriated trade-off between stability & standardization (ST$^2$)and rate of expression. We define the ST$^2$ index as $$ F_{st2} = (F - 0.012)/0.51 $$ with $$F = \frac{R'}{N_{NS}} e^{-\frac{\Delta\epsilon}{k_BT }}, $$ where $F_{sty}$ has been scaled to run from 1 to 0. We have used a fitting formula for the rate of expression relative to the maximum rate [3] as a function of the binding energy $\Delta\epsilon$, $$ r_T = 2.29 - 0.28 \frac{\Delta\epsilon}{kT}. $$ With this approach, we can infer the ST$^2$ level and the promoter sequence for a given expression rate. We show here our app in javascript:

Comparison with wet-lab data

Our wet-lab results show that de Biobrick 5 is the more standard one. Our model predicts that systems with less modules or factors will be the more stable and standard. Biobrick 5 has a tetracycline depending promoter, whose limitating step is the binding of tetracycline to TetR. In this scenario, the intracellular concentration of aTc and the binding to the TetR repressor do not depend on the strain, so there is only one factor in which the strain could have an important effect: the TetR intracellular concentration, which has a constitutive expression. Analyzing the other Biobricks, we found that the strains with less expression in the constitutive promoter are the ones with higher expression when we use the Ptet promoter (but JM109 and BB1), which agree with our prediction.

Open License

We model the susceptibility to patent a human creation by $$ P.I. = M \cdot N \cdot A ^{log(\frac{A+I}{2})}, $$ where P.I. is the patentability index, M is the morality, which is a non-legal term accepted in the E.U., and N refers to novelty. Both factors are described by a discrete variable which can be valued 0 or 1 depending on the results of a questionnaire. A is the inventive step and I is the industrial application which can be valued from 1 to 10 depending of the results of a test.

Inventive is the most relevant variable but is regulated by the industrial application. The highest inventive step can be reduced by almost a factor 2 with the lowest industrial application. On the contrary, the lowest inventive step is not affected by the industrial application. After computing all the discussed factors, the patentability index can run from 0 to 10. The reader can find a detailed discussion of the questionnaires, tests, analysis of the index and references in the area Human Practices.

Applications to other projects

With our mathematical model we have developed a tool that allows us to predict how stable a biological system is (in this case, we will study the stability with temperature).

In our project, the system was very simple: E. coli carrying a reporter gene. We study the changes on activity with temperature variations:

$$ \begin{equation} \frac{A_{c}}{A^0_{c}} = \prod_i \frac{V_i}{V^0_i} . \end{equation} $$

Where $V$ is any significantly variable factor, and A_c⁰ is the activity at the reference conditions (in this case, 37 ºC).

In the model of our project, we consider (see Stability):

$$ \begin{equation} \frac{Ac}{Ac^0} = \frac{N_Pr}{N^0_Pr} =\frac{r_T}{r^0_T} \frac{r_P}{r^0_P} . \end{equation} $$

$$ \begin{equation} \frac{r_T}{r^0_T} = \frac{Prob_B}{Prob^0_B} = \frac { 1+\frac{R^0}{N_{NS}} e^{-\frac{\Delta\epsilon}{k_BT }} } { \frac {R^0}{R}+\frac{R^0}{N_{NS}} e^{-\frac{\Delta\epsilon}{k_BT }}} . \end{equation} $$

$$ \begin{equation} \frac{R^0}{R} = \frac {0.027 (T^0)^2 - 2.12 T^0 + 42.2} {0.027 T^2 - 2.12 T + 42.2} . \end{equation} $$

$$ \begin{equation} \frac{r_P}{r^0_P} = \frac{\frac{5}{4}}{\frac{r^0_E}{r_E}+\frac{1}{4}}. \end{equation} $$

$$ \begin{equation} r_E = \frac{1}{\frac{1}{A_1}e^{\frac{E^1_a}{RT}}+\frac{1}{A_2}e^{\frac{E^2_a}{RT}}+\frac{1}{A_3}e^{\frac{E^3_a}{RT}}}, \end{equation} $$

But if we create another, more complex system, more parameters will vary, and more significantly.

For an enzyme, the activity is:

$$ \begin{equation} A_{c} = k e^{-\frac{E_a}{RT } } \frac{N_{Pr}}{1+e^{\frac{\Delta G_F}{RT } } } . \end{equation} $$

Where N_Pr is the number of total proteins, ΔG is the free energy of the folding process, k is the Arrhenius constant of the chemical reaction and E_a is the activation energy of the catalyzed reaction ()

In a first approximation, we can assume that all genes are controlled by the same promoter, and the same folding energy for every protein, and if the proteins are enzymes, we assume the same activation energy for every reaction (for the mean activation energy we used [¹]) is:

$$ \begin{equation} \frac{A_{c}}{A^0_{c}} = {(\frac{r_P}{r^0_P} )^n}{({\frac{r_P}{r^0_P}}{\frac{\frac{1}{1+e^{\frac{\Delta G_F}{RT } } }}{\frac{1}{1+e^{\frac{\Delta G_F}{RT^0 } } }}}{\frac{e^{-\frac{E_a}{RT } }}{e^{-\frac{E_a}{RT^0 } }}})^m} \end{equation} $$

Where:

n: number of different mRNAs,

m: number of proteins,

Ea: activation energy of the catalyzed reaction (we can use the mean value of E. coli reaction, from [1])

ΔG: protein folding energy. Using the equation:

$$ \begin{equation} \Delta G_F (L,T) = \Delta H(L) + \Delta C_p (L)(T-T_h) - T\Delta S(L) - \Delta C_p (L)T\ln\frac{T}{T_S} \nonumber, \end{equation} $$

Where: $$ \begin{eqnarray} \Delta H(L) &=& -5.03 L - 41.6 (kJ/mol)\nonumber, \\ \Delta C_p(L) &=& -0.062 L + 0.53 (kJ/mol)\nonumber, \\ \Delta S(L) &=& -16.8 L -85 (J/mol)\nonumber. \end{eqnarray} $$

T_h = 373.5 K and T_s = 385 K.

Where:

T: temperature (K).

L: Length of the protein (aa). The E. coli mean protein length is 300 aa.

Figure 1: Variation of activity with temperature for a system that have 1, 2, 3, 4 and 5 genes.

Figure 2: Variation of synthetized protein with temperature for a system that have 1, 2, 3, 4 and 5 genes.

Figure 3: Variation of activity with temperature for a system that have 1, 2, 3, 4 and 5 proteins, but only one mRNA.

Conclusions

- Our main conclusion is that stability decreases as the amount modules involved increases. If there are many modules, it will be more stable if just one of these controls the limiting step.

- As the number of different mRNAs decreases, the stability of the system grows at lower temperatures and the activity of these increases at larger.

Other iGEM projects

Team Slovenia 2010. DNA CODING BEYOND TRIPLETS developed a system to use a DNA sequence as a scaffold for enzymes, to optimize a metabolic pathway (better explained on their wiki). We can apply our model as follows:

$$ \begin{equation} \frac{A_{c}}{A^0_{c}} = {(\frac{r_P}{r^0_P} )^n}{({\frac{r_P}{r^0_P}}{\frac{\frac{1}{1+e^{\frac{\Delta G_F}{RT } } }}{\frac{1}{1+e^{\frac{\Delta G_F}{RT^0 } } }}}{\frac{e^{-\frac{E_a}{RT } }}{e^{-\frac{E_a}{RT^0 } }}})^m}{({{\frac{\frac{1}{1+e^{\frac{\Delta G_Z}{RT } } }}{\frac{1}{1+e^{\frac{\Delta G_Z}{RT^0 } } }}}}{{\frac{\frac{1}{1+e^{\frac{\Delta G_I}{RT } } }}{\frac{1}{1+e^{\frac{\Delta G_I}{RT^0 } } }}}})^z}{{\frac{\frac{1}{1+e^{\frac{\Delta G_D}{RT } } }}{\frac{1}{1+e^{\frac{\Delta G_D}{RT^0 } } }}}} \end{equation} $$

Where:

n: number of required transcripts

m: number of required translated proteins

z: number of required zinc-finger—DNA interactions.

E_a: activation energy of every chemical reaction

ΔG_F: Folding energy of every enzyme required for the pathway

ΔG_Z: Folding energy of every zinc-finger domain ([¹⁸])

ΔG_I: Interaction energy of a normal zinc-finger-DNA interaction ([¹⁸]).

They designed it is looking for a more efficient synthesis of violacein. As we do not have any data, we will use the mean activation energy (50 kJ/mol) of every catalyzed reaction. [¹].

Figure 4: Variation of activity with temperature of a metabolic pathway with 1, 2, 3, 4 and 5 enzymes.

As we can see, the result is the same than in the previous analysis, which means that the new modules (zinc-finger and DNA) do not affect significantly to the stability.

Team Heidelberg 2013. PHILOSOPHER’S STONE developed a system to synthesize non ribosomal peptides. One of the simplest mechanisms was the synthesis of indigoidine (see their wiki) :

For this synthesis, we need:

1. Synthesis or acquisition of Glu (assuming that the bacteria acquires from the external medium with a specific transporter).
2. 2. Transform Glu in cGlu
3. 3. Transform cGlu in Ind (assuming that there ir one lower step).

$$ \begin{equation} \frac{A_{c}}{A^0_{c}} = {(\frac{r_P}{r^0_P} )^3}{(\frac{r_P}{r^0_P})^3}{({\frac{\frac{1}{1+e^{\frac{\Delta G_F}{RT } } }}{\frac{1}{1+e^{\frac{\Delta G_F}{RT^0 } } }}}{\frac{e^{-\frac{E_a}{RT } }}{e^{-\frac{E_a}{RT^0 } }}})^2}{\frac{\frac{1}{1+e^{\frac{\Delta G_T}{RT } } }}{\frac{1}{1+e^{\frac{\Delta G_T}{RT^0 } } }}}{\frac{e^{-\frac{E_aT}{RT } }}{e^{-\frac{E_aT}{RT^0 } }}} \end{equation} $$

Where:

E_a: activation energy of every chemical reaction.

E_aT: activation energy of the transport process.

ΔG_F: Folding energy of every enzyme required for the pathway

ΔG_T: Folding energy of the transporter.

As above, we do not have any data, so we will use the mean activation energy (50 kJ/mol) of every catalyzed reaction. From [¹⁹], we can calculate the activation energy for the Lysine transport through membranes: 33.244 kJ/mol.

Figure 5: Variation of the production of Ind with temperature.

In this case, we find variation between this system and a generic described above (the optimum temperature is lower and the maximum value of Ac/Ac ⁰) is also lower than the one with “n,m = 3”. This means that exists a significant effect when we use the module of transport of Gln.

References

1. Ghosh K., Dill K. (2010). Cellular proteomes have broad distributions of protein stability. Biophys J 99: 3996–4002.
2. Soini, J. et al. (2005) Transient increase of ATP as a response to temperature up-shift in Escherichia coli. Microb. Cell Fact. 4, 9.
3. Brewster R.C., Jones D.L., Phillips R.(2012) Tuning promoter strength through RNA polymerase binding site design in Escherichia coli. PLoS Comput Biol 8: e1002811.
4. Skelly S. et al. (1987) Correlation between the 32-kDa s factor levels and in vitro expression of Escherichia coli heat shock genes. Proc. Natl. Acad. Sci. USA 84:8365–8369.
5. Zaritsky, A. (1982). Effects of growth temperature on ribosomes and other physiological properties of Escherichia coli. J. Bacteriol. 151:485-486.
6. Bremer, H., Dennis, P. P. (1996) Modulation of chemical composition and other parameters of the cell by growth rate. Neidhardt, et al. eds. Escherichia coli and Salmonella typhimurium: Cellular and Molecular Biology, 2nd ed. chapter 97.
7. Piperno J.R., Oxender D.L. (1968) Amino acid transport systems in Escherichia coli K-12. J Biol Chem. 243(22):5914–5920.
8. Thompson R.C., Dix D.B., Karim A.M. (1986) The reactions of ribosomes with elongation factor Tu·GTP complexes. Aminoacyl-tRNA-independent reactions in the elongation cycle determine the accuracy of protein synthesis. J Biol Chem, 261:4868–4874.
9. Bilgin N. et al. (1992) Kinetic properties of Escherichia coli ribosomes with altered forms of S12. J Mol Biol. 224 (4):1011-27.
10. Grigorenko B.L. et al. (2008) Mechanism of the chemical step for the guanosine triphosphate (GTP) hydrolysis catalyzed by elongation factor Tu. Biochim Biophys Acta 1784:1908–1917.
11. Gindulyte,A. et al. (2006) The transition state for formation of the peptide bond in the ribosome. Proc. Natl. Acad. Sci. U. S. A., 103:13327–13332.
12. Wilson,K. S., Noller, H.F. (1998) Mapping the position of translational elongation factor EF-G in the ribosome by directed hydroxyl radical probing. Cell, 92: 131–139.
13. Erban R, Chapman SJ and Maini PK (2007) A PRACTICAL GUIDE TO STOCHASTIC SIMULATIONS OF REACTION-DIFFUSION PROCESSES link
14. Klumpp E (2011) Growth-Rate Dependence Reveals Design Principles of Plasmid Copy Number Control PLoS One Vol 6 Issue 5 link
15. Marguet P, Tanouchi Y, Splitz E, Smith C and You L (2010) Oscillations by Minimal Bacterial Suicide Circuits Reveal Hidden Facets of Host-Circuit Phisology PLoS One Vol 5 Issue 7 link
16. Paulsson J and Ehrenberg M (2001) Noise in a minimal regulatory network: plasmid copy number control Cambridge Press link