"

From 2013hs.igem.org

Mathematical modeling

Deterministic model

Our deterministic model represents the change in time of the concentrations of mRNAs and their corresponding proteins. It assumes that the variables behave continuously and obey kinetic rules that can be represented by constants. We are aware that, in practice, the components and variables in the model may not fall under the assumptions of a deterministic model and that there is always the chance that noise effects are being grossly underestimated; however, we propose this model as a general framework for thermoregulator modeling, upon which further work can be done.

Our genetic system is composed by two parts each, one placed in different vectors: the first in the plasmid 2-4B C/C, consists of the transcription factor cI (c0051) under the regulation of a promoter with a riboswitch (J23100 and k115017 ) that can be activated or repressed depending on the temperature; we call this our thermosensor. We found this piece in the iGEM TUDelf-2008 Team and saw we could use it because of its specific translation temperature at 32ºC. The second part located in the vector 2-6B K/K consist of two genes -the insecticide protein VIP3Ca3 and the green fluorescent protein, GFP - that are under the regulation of cI.

The thermosensor is repressed when the temperature is below 32 °C, it does not translate because the mRNA takes a secondary structure just as the image below and cannot enter the ribosome for translation; when the temperatures reaches a range between 32 and 37 °C, the thermosensor is activated, the mRNA unfolds into a straight structure that now have access to the ribosome, start translation, be passed to aminoacids and form a protein .

cI (c0051) is a gene that constitutively represses the promoter r0051. Once our promoter isnt affected by the cI(c0051) the promoter is activated starting the production of Vip3Ac3 and GFP. This is means that a range below 32°C to 37°C ables the production of the Vip and GFP.

Equations

\begin{equation} \large \frac{d[mC]}{dt} = \alpha_{1} - \mu_{1}[mC] \end{equation}

Describe transcription of cI (C0051) mRNA over the change in time; that it is equal to the transcription rate of cI less of the degradation rate of cI(mRNA). Said in other words alpha(a) is how fast mRNA is produced by the mRNA polymerase minus the degradation speed of the mRNA produced (µ[mC]).

\begin{equation} \large \frac{d[C]}{dt} = \alpha_{2} \cdot f_{RBS} \cdot [mC] - \mu_{2}[C] \end{equation}

Describe the translation of cI (C0051) over the change in time, that it is equal to the translation rate of cI(RNA polymerase velocity to make mRNA) by the function RBS of E.coli ribosomes (expressed in equation 5) multiplied by the transcription of cI(equation 1). This part is the amount of cI(c0051)mRNA present in our bacteria, if this is zero (no mRNA produced) we would get no translation results; no material, no product. We need to reduce to the protein translated the degradation rate of cI(protein), because in all organic material they tend to decompose and we need to show this in the equation.

\begin{equation} \large [mC]_{max} = \frac{\alpha _{1}}{\mu _{1}} \end{equation}

In order to get a maximum cI mRNA capacity we described the transcription of cI (C0051) when the transcription of cI (C0051) mRNA over the change in time is 0 because we supposed at that specific time the transcription will be in a maximum production. Then µ[mC] is isolated; getting that µ[mC] is the transcription rate of cI divided by the degradation rate of cI.

\begin{equation} \large [C]_{max} = \frac{\alpha _{1} \cdot \alpha _{2}} {\mu _{1} \cdot \mu _{2}} \end{equation}

It describes the translation of cI (C0051) when the translation of cI (C0051) over the change in time is 0 because we supposed at that specific time the translation will be in the maximum production, such as the third equation. It is solved and the [C]max is equal to the transcription of cI by the translation of cI both divided by the degradation rate of cI(mRNA) multiplied by the degradation rate(protein)

\begin{equation} \large f_{RBS} = \left\{ \begin{array}{rcl} 0 & \mbox{if} & t \geq ON \\ 1 & \mbox{if} & t < OFF \end{array} \right. \end{equation}

The function of RBS of ribosome says that at the range of temperature between 32 ºC and 37 ºC, the production of Vip3Ca3 will be null such as GFP and that´s why we relation with the zero and the term of greater or equal than is ON (cI), and the term 1 express that at below temperature (32ºC) it will start the production of Vip3Ca3 and GFP (ON) and the cI off, the riboswitch. This equation express the primary and secondary structures the riboswitch k115017 can take depending on the temperature. When temperature is below 32ºC the riboswitch have a special structure named hairpin structure that avoids the sequence to entering the ribosome to translate, but at higher temperatures the mRNA is at a normal shape that can enter the ribosome to create a protein.

\begin{equation} \large \frac{d[mV]}{dt} = \alpha_{3} \cdot \frac{K_{D}^h}{K_{D}^h + [C]^h} - \mu_{1}[mV] \end{equation}

It describes the transcription of Vip3Ca3 mRNA over the change in time; that it is equal to the transcription rate of Vip, regulated depending on the production of cI, less the degradation rate of Vip (mRNA)

\begin{equation} \large \frac{d[V]}{dt} = \alpha_{4} \cdot [mV] - \mu_{4}[V] \end{equation}

It describes the translation of Vip3Ca3 over the change in time, that it is equal to the translation rate of Vip multiplied by the transcription of Vip (equation 7) less of degradation rate of Vip (protein).

\begin{equation} \large [mV]_{max} = \frac{\alpha _{3}}{\mu _{3}} \end{equation}

It describes the transcription of Vip3Ca3 when the time is 0 because we supposed at that specific time the transcription will be in the maximum production. Then it is solved the octave equation and is the transcription rate of Vip3Ca3 divided by the degradation rate of VIP.

\begin{equation} \large [V]_{max} = \frac{\alpha _{3} \cdot \alpha _{4}} {\mu _{3} \cdot \mu _{4}} \end{equation}

It describes the translation of Vip3Ca3 when the time is 0 because we supposed at that specific time the translation will be in the maximum production, such as the ninth equation. It is solved and the equation is equal to the transcription of Vip by the translation of Vip is divided by the degradation rate of Vip mRNA by the degradation rate of Vip protein

\begin{equation} \large \frac{d[mG]}{dt} = \alpha_{5} \cdot \frac{K_{D}^h}{K_{D}^h + [C]^h} - \mu_{5}[mG] \end{equation}

It describes the transcription of GFP mRNA over the change in time; that it is equal to the transcription rate of GFP, regulated depending on the production of cI, less the degradation rate of GFP (mRNA)

\begin{equation} \large \frac{d[G]}{dt} = \alpha_{6} \cdot [mG] - \mu_{6}[G] \end{equation}

It describes the translation of GFP over the change in time, that it is equal to the translation rate of GFP multiplied by the transcription of GFP (equation 11) less of degradation rate of GFP (protein).

\begin{equation} \large [mG]_{max} = \frac{\alpha _{5}}{\mu _{5}} \end{equation}

It describes the transcription of GFP when the time is 0 because we supposed at that specific time the transcription will be in the maximum production. Then it is solved the octave equation and is the transcription rate of GFP divided by the degradation rate of GFP.

\begin{equation} \large [G]_{max} = \frac{\alpha _{5} \cdot \alpha _{6}} {\mu _{5} \cdot \mu _{6}} \end{equation}

It describes the translation of GFP when the time is 0 because we supposed at that specific time the translation will be in the maximum production, such as the 13 equation. It is solved and the equation is equal to the transcription of GFP by the translation of Vip is divided by the degradation rate of GFP mRNA by the degradation rate of Vip protein

Parameters and variables

Our gene circuit is made of six different variables: the concentrations of three proteins (cI, VIP and GFP) and their respective mRNA inside a cell. In table 1, the symbols for each variable are shown. Proteins are represented by a single letter and their mRNAs are represented by that same letter with a lowercase "m" before it.

Table 1.- Variables

Symbol	Definition	Gene size in bp	Source
mC, C	Transcription factor cI (mRNA and protein)	775	http://partsregistry.org/wiki/index.php?title=Part:BBa_C0051
mV, V	Insecticide protein VIP (mRNA and protein)	2412	http://www.ncbi.nlm.nih.gov/nuccore/HQ876489
mG, G	Reporter protein GFP (mRNA and protein)	876	http://parts.igem.org/wiki/index.php?title=Part:BBa_E0240

To parameterize our model, we chose to follow the approach of team Beijing 2009; they propose a relationship between the gene length in base pairs and the maximum transcription rate and, similarly, between the protein length in amino acid numbers and the maximum translation rate. Assuming that the number of polymerases and ribosomes are the average values determined for E. coli, the following equations are used:

Parameters

Symbol	Definition	Values	Formula	Source
α1	Transcription rate of cI	5.6	4200/Gene Length (nM/min)	https://2009.igem.org/Team:PKU_Beijing/Modeling/Parameters
α2	Translation rate of cI	9.6	2400RBS/Protein Length	https://2009.igem.org/Team:PKU_Beijing/Modeling/Parameters
α3	Transcription rate of VIP	1.74129353	4200/Gene Length (nM/min)	https://2009.igem.org/Team:PKU_Beijing/Modeling/Parameters
α4	Translation rate of VIP	2.985075	2400RBS/Protein Length	https://2009.igem.org/Team:PKU_Beijing/Modeling/Parameters
α5	Transcription rate of GFP	5.53359684	4200/Gene Length (nM/min)	https://2009.igem.org/Team:PKU_Beijing/Modeling/Parameters
α6	Translation rate of GFP	9.486166	2400RBS/Protein Length	https://2009.igem.org/Team:PKU_Beijing/Modeling/Parameters

For all the variables the degradation rate is expressed by the formula (ln(2)/half life)+(ln(2)/division time), with the same division time of e. coli (30 min), because all the process occurs within it. The only thing that change is the half time; for cI, VIP and GFP (mRNA) is 6.8 minutes and for cI (Selinger, GW, et al., 2003), VIP and GFP protein is more than 10 hours (Varshavsky, (1997) and Tobias et al., 1991).

Degradation

Symbol	Definition	Values	Formula	Source
μ1μ3,μ5,	Degradation rate of cI (mRNA	0.18063836	Half life = 6.8 min, Division time = 30 min	(Selinger, GW, et al., 2003)
μ2,μ4,μ6	Degradation rate of cI (protein)	0.03885825	Half life > 10 h; division time = 30 min	(Varshavsky, (1997) and Tobias et al., 1991)

Simulations

The saturation values for mRNAs and proteins were calculated analytically; but since there are several variables, it becomes complicated to integrate by analytical methods, so we use methods of numerical integration in a computer program by called Simulink. The values of the parameters (rate of transcription, translation, degradation and dissociation) are the ones we have found so far, but we continue researching in order to improve and expand our model.

cI inactive-active simulation

The graphs below shows the change in concentration with respect to time of all the proteins of our circuit in a cI inactive-active state. This simulations are showing our model in a temperature below 32ºC the first 200 minutes making the variables of cI protein production in their minimum value for this initial time, creating the possibility of Vip3Ca3 and GFP proteins to be processed in the bacteria at maximum capacity (maximum capacities shown in the parameter table). Past the half hour the temperature is higher of our ideal parameters of cI production (above 32ºC) allowing the transcription and translation of cI protein in a factor of 1 (represented in the equation number 5 showed above) this is the rise in the blue graph at the time 200 minutes, inhibiting the formation of the other two system parts: Vip3Ca3 and GFP proteins, simulated as the decrease of the purple and green graphics. Whose percents in the E.Coli bacterium drops to the minimum until the cI production stops again and the process restarts.

cI active-inactive simulation

The graphics that are shown represent the change of the concentration of each part in relation to time of all proteins in our cicuit, including all the parameters that we presented and the equation 5 (fRBS), where we represent the repression and activation of the riborswitch making the value of its function 0 or 1. The first graph simulation is showing our model in a temperature between 32 to 37ºC in the first 200 minutes. As you can see, in the graph of cI (c0051), it is active and produced constantly, repressing the Vip3Ca3 and GFP proteins, but when the temperature changes below than the 32 to 37ºC the cI it is now turning off allowing the transcription and translation of Vip3Ca3 and GFP so changes the graphic, but We can appreciate the GFP is producing more than the Vip3Ca3, this because of the bigger base pair size of Vip3Ca3 than GFP that make it longer to be synthesis by the ribosome. This bp length of those proteins are 2412 and 876 respectively.

Click here to download the Simlulations in MathLab

Probability Equation

We tried to modify our model with a probability equation from an article of the National Institute for Mathematical and Biological Synthesis, Knoxville, Tennessee, United States of America published in the internet page PLOSE ONE with the name Is Thermosensing Property of RNA Thermometers Unique? This work looks for the probability of finding a RNA at an openness way.

This could described the chance of finding our riboswitch k115017 open or closed, because it represented based on their use of parameters a more smooth chance in the fold of the thermosensor.

\begin{equation} \large p_{i}\left ( T \right )= \frac{e^{a_{i}+b_{i}T}}{1+e^{a_{i}+b_{i}T}} \end{equation}

"Pi(T) is the probability of finding the window at position “i” in a gene, open at temperature (C), ai and bi are the intercept and slope parameters of how the log-odds of finding an open window at position “i”, log(pi(T)/1- pi(T)), changes with temperature. The ratio –ai/bi indicates the temperature at which the probability of openness of a window is 0.5."http://www.plosone.org/article/info%3Adoi%2F10.1371%2Fjournal.pone.0011308

Our team used the last data of –ai/bi = T when openness is 0.5, to get “ai” expressed in b ( ai= -T*bi) and have our own goods ; to find b the authors of the articles, compared across the 100 genes, found that values were not significantly greater than zero for 24 genes at p value is 0.05, so b value is 0.157. This implies that a small fraction of genes did not show a significant change in openness of its RBS with temperature over the range of temperatures considered. Then we place our control temperature in the isolated equation and get a= -32*0.157 to get the “a” value that was negative 5.024.

But there were some problems using the pi(T) formula in the matlab simulations like malfunction of the equation at the moment of seen graphs or the incorrect insertion of parameters in the simulation, so finally the team decided to leave it apart of the graphics.

We mention this special formula because it was the solution of the fast growth and decrease of the concentrations of the RNAs and proteins in E.coli. We run out of time and it was almost imposible to experiment with it and get better results, but heres the small advance we made with all this.