Team:SJTU-BioX-Shanghai/Modelling

From 2009.igem.org

Revision as of 17:05, 18 October 2009 by Pnstontale (Talk | contribs)

Team Logo ButtonTeam Logo ButtonProject ButtonTeam Logo Button

Project introduction. Inspired by the natural regulator of circadian bioclock exhibited in most eukaryotic organisms, our team has designed an E.coli-based genetic network with the toxin-antitoxin system so that the bacterium oscillates between two states of dormancy and activity (more...)

Mathematical Modeling

Contents

Motivation

From the background of our system, one can see that our model is designed to regulate the bio-clock of E.coli. However, this system is more complicated than one can think in that it consists of a variety of reactions, which is quite hard for us to analysis without mathematics. So it is indispensible to justify that the concentration of relE will oscillate with a specific periodicity. Moreover, understanding the important parameters of the system allows one to choose appropriate promoters and rbses to make our system work. Therefore, we have created the mathematical models to prove our project is feasible, to guide our experiment, and to predict the length of E.coli’s sleeping time.


Introduction

First of all, let's pick out the three key genes, which names relE, relB and Lon, from the original network. Then the interaction between them can be expressed as a graph as follows.

File:SJTU09 Modeling html 1.jpg
Interaction between relE, relB and Lon

Actually, there are four types of reactions in these systems, which are highlighted by different colors.

  • The arrow marked by green is a single gene expression
  • The arrow marked by yellow is the relE cleavage process
  • The arrow marked by blue is the relB-Lon hydrolysis
  • the arrow marked by purple is the relE and relB polymerization

From the diagram mentioned above, obviously almost each two genes or its mRNAs will have a bio-chemical reaction.


Assumptions and hypothesis

  1. The number of ribosome in E.coli is sufficiently large at all times
  2. The cleavage reaction caused by relE is a enzyme process
  3. Any protein and its related mRNA generated at a constant rate
  4. Either protein or mRNA degrades in direct ratio to its concentration

In order to make the quantitative analysis, we should simplify our model by making the above-mentioned assumptions. The first one is the vital assumption in our system. Since a cell may have many thousands of ribosomes, the number being positively correlated with growth rate[1], we consider this assumption reasonable. In fact we can overlook the opposite situation, which will lead us to consider the piecewise ODEs (ordinary differential equations).Unfortunately, such kind of equations cannot be solved by Matlab. The second assumption is supported by The Bacterial Toxin RelE Displays Codon-Specific Cleavage of mRNAs in the Ribosomal A Site[2] while the last two assumptions are illustrated by Stochastic simulations Application to molecular networks[3].


Model development

Single gene expression

According to assumption 3 and 4, we have made, we will have the following six equations:

SJTU09 Modeling html1.gif

Where α1, α3and α5 represent the strength of promoters. α2, α4 and α6represent the strength of rbs. βi(i=1,2...6) represent the degradation rate of the associated proteins or mRNAs.


relE cleavage process

Based on assumption 1 and 2, we choose the classical Michaelis-Menten equations to describe the relE cleavage process:

SJTU09 Modeling html2.gif

Where km is the Michaelis constant, which is only related with the affinity of the enzyme (relE) for the substrate. kcati(i=1,2,3) is related with the cleavage rate of relE.


relB-Lon hydrolysis

For relB-Lon hydrolysis, we have the following equations:

SJTU09 Modeling html3.gif

Where km is the Michaelis constant, which is related with the affinity of the Lon. kcat0 is related with both relB and Lon.


relE and relB polymerization

For relE and relB polymerization, we choose the second-order reaction:

SJTU09 Modeling html4.gif

where k represents the reaction coefficient.


Modeling the full system

Based on the former analysis, we can generalize our model into the following six ordinary differential equations:

SJTU09 Modeling html5.gif

The following table shows the description of the parameters.

SJTU09 Modeling html6.gif

However, since we proposed these equations, we began to realize some problems for us before we can get satisfactory results from them:

  1. From the table above, we can see that there are heaps of parameters for us to find out. And some of them are even impossible for us to get.
  2. We are only interested in the amounts of some proteins, mainly the relB and relE. However, there are many other things, such as the various mRNA, for us to cope within the equations.
  3. The complexity of the relations among different elements adds great difficulties for us to calculate and analyze.

So, we decided to only focus on the relations between two of the three proteins, which mostly interests us.

Having the above notion in our minds, we proposed our simplified new model. We keep the majority of the terms in our new model, except the relE cleavage term. Since we are not pretty sure about the principle of relE cleavage process, we choose the classical hill equations. Now, we can write down our new model in the following three ordinary differential equations:

SJTU09 Modeling html7.gif

Now, we can see the new model is much simpler for Matlab to cope with. However, we still should find the values of some parameters from a variety of databases. They are listed as below:

SJTU09 Modeling html8.gif

Since most parameters can be found in several papers, what we should do is to choose the appropriate parameters, such as the promoter strength and the degradation rate, to make the system oscillate. Fortunately, Matlab can handle all these works, which will be presented in the next section.


Simulation

At the beginning of the simulation part, let's replace the parameters with simple characters.

  • a: relE promoter strength
  • b: relE degradation rate
  • c: Reaction coefficient(relE and relB)
  • d: relB promoter strength
  • e: relB degradation rate
  • f: reaction coefficient(Lon and relB)
  • g: Lon promoter strength
  • h: Lon degradation rate

Now, we are ready to calculate the right value of these eight parameters by Matlab. Remember our target is to select the proper parameters to make the system oscillate. If so, the bacteria will fall asleep when relE rises to a high level and the bacteria will be active when relE falls down to a low level.

We set the following value for the parameters, to see if they satisfy our requirement.

a = {11,16,21,26,31};b = {1,6};c = {1,6,11};d = {21,26,31};

e = {1,6};f = 46;g = {21,26,31};h = {1,6}


At first we find that the majority of the systems with specific parameters are asymptotic stable. In order to figure out the reason why they cannot show us oscillations, we consider the oscillation criterion.

How to judge a system that cannot oscillate?

With the help of the criterion, we can remove many situations which will lead to asymptotic stable. This will reduce the burden of matlab to some degree. Matlab tells us the following values parameter value can help the system to reach an oscillation state.

SJTU09 Modeling html9.gif

And these parameter value listed below will lead to a damping oscillation state.

SJTU09 Modeling html10.gif

After 5 hours work by Matlab, one can see three situations with specific parameter values.

oscillate with periodicity

oscillate with periodicity


damped oscillations

damped oscillations


equilibrium states

equilibrium states


Conclusion

Since we have three types of promoters (j23110,j23118,j23116) with their relative strength 1:0.6:0.07, we are now able to choose the proper one to make our system oscillate. All the well-performing situation are listed as below:

SJTU09 Modeling html11.gif

From these pictures, it is obvious that a larger value of “a” means a higher probability of oscillation. In addition, a less value of “d” means a higher probability of oscillation. Moreover, the value of “g” plays an important role in oscillation system. We find that when we set the value of “g” lower than 13, then this system will fall into chaos. When we set the value of “g” from 18 to 22, we observed the oscillations. When we set the value of “g” higher than 25, the system will reduce to damping oscillations and reach the equilibrium state immediately. In conclusion, we can set g=22 to have the best effect. According to our former discussion, we are now able to choose the appropriate promoters to make our system oscillate.


References & Attachment

[1] Michael T. Madigan, John M. Martinko, Brock biology of microorganism 11th edition, Pearson Education International 2006, ISBN:0-13-196893-9

[2] Kim Pedersen, Andrey V. Zavialov,Michael Yu. Pavlov, Johan Elf,Kenn Gerdes, and Mans Ehrenberg. The Bacterial Toxin RelE Displays Codon-Specific Cleavage of mRNAs in the Ribosomal A Site. Cell, Vol. 112, 131–140, January 10, 2003

[3] Didier Gonze, Stochastic simulations Application to molecular networks March 27, 2007




Since this site doesn't support Tex format for formulas, we provide our
source document SJTU09_Modeling.zip for you to download.