Team:Nanjing-China/modeling

Our project is aimed at designing a system to detect and handling pollutants and toxins in the environment especially microcystin. Our system is consisted of several parts, such as the positive feedback circuit, the time sensor, the toxin detector, the suicide circuit and so on. We can find the toxin and measure its concentration by the detector. The time sensor is designed to tell the time to make sure all the measurement is at the same time. With the positive feedback circuit, we can make much more enzyme or binding protein to degrade or catch the toxin. The suicide circuit is used to satisfy the requirement of life safety. We use riboswitches to detect the toxin, which means we need to design several riboswitches and evaluate its efficiency. For the evaluating riboswitch, we use RNAstructure to predict the secondary structure of the sequences and its free energy. With the free energy and some other parameters, we can obtain the kinetics model of the RNA.
The solution of detector model needs the help of MATLAB for solving partial differential equations. We also use Simbiology of MATLAB to simulate all the other circuits.
We compare the modeling results with experiment results. Some prediction can also be found by modeling

The POSITIVE FEEDBACK SYSTEM is used to enlarge the production of enzymes and promote the degradation of toxins.
There are four designs. The first is a blank control. The second one uses luxI to cause positive feedback. The third uses luxR while the fourth uses both luxI and luxR.

There are four designs. The first is a blank control. The second one uses luxI to cause positive feedback. The third uses luxR while the fourth uses both luxI and luxR.
For specific information of modeling, please look at the supplement.

Simulating by MATLAB Simbiology, We get several conclusions which can be compared with the results of the experiment.

From the picture above, we can obtain the conclusion that our positive feedback system works in theory. The binary positive feedback system[ luxI+ and luxR+] maintains much more RFP than the control system in only 300 seconds. The luxI+ system is as good as the binary system. The luxR+ system is only a little better than the control.

As time goes on, the binary system and luxI+ system is still better than luxR+ system and control system, and the differences between the binary system and luxI+ system is two little to be found.
Someone may ask why luxI+ system and binary system is better than luxR+ system, while it is clear that they should be better than the control system. We assume that these E.coli. already remain in a stable state called initial state after cultivating for a long time. In this state, the concentration of luxR is almost stable because the work of the Pcons. As a result, the feedback by luxR is not as useful as that of luxI, which increase more magically[ Results about this will be listed in the supplement.]. In the supplement, we also model about the stability of the system.
In a word, by using positive feedback system, we can produce more enzyme to do the degradation job, and the system is little affected by the strength of Pcons[ For more, please refer to the supplement.].

The time sensor shows the time for reading the position of the red fluorescence in the detection system. It is situated in the AHL producing cell in the detection system following the luxI. This part was tested by using Plac and lacI.

With the help of MATLAB, we have gotten some results.

The diagram above shows that RFP increase as time goes on as our expects. However, we would rather use a stable amount of RFP to make sure that we always read the results of sensor at the same time. A stable concentration of RFP is also a reflect of the environment, which may also affect the result of the sensor..

From this diagram, we can get that RFP reach its stable concentration at the time of 40000s, which will be our time for reading the result of sensor system.
With the help of time senor system, we can always read the results of detector at proper time.

The suicide system is used to wipe off our bacteria to reduce the probability that these E.coli. may escape to the environment, posing a potential threat to the health of human. To realize our aim, we design a circuit that makes the cells dead where there are no lactose or the aimed toxins. To make it simple, we just test whether the cells will die without lactose for this section.

We are convinced by the diagram that the lysozyme will increase with time as long as we remove lactose. The increase of lysozyme will lead the cells to death after work as we hoped.

A model was chosen from a paper[Beisel CL, Smolke CD. Design principles for riboswitch function. 2009. PLoS Comput Biol.
For more information, please consult reference[28]] to test the function of our riboswitches. With the riboswitches, we can detect our aimed toxins, report its concentration and degrade it.
The riboswitch can response to specific ligand, changing the stability of mRNA and in result the concentration of AHL[ It is used in detector.] and enzymes.

Our riboswitch worked in the level of RNA stability. This part is composed of a aptamer
and a hammerhead ribozyme. The ribozyme can degrade the mRNA when the aptamer is not combined with toxin.

Using RNAstructure, we get the free energy of different folding state of the RNA, with which we can work out the important parameter k₁⁹, with which we could estimate the value of k₁ and k₁’. Here we will only list the two fold structures of MC10 [ The folding structure of other riboswitch will be listed wen the attachment file as CT files.] [ The unit of energy results in RNAstructure is kcal/mol. Energy results from different software may be a little different.]. The first can be combined with ligand and is stable, while the second without ligand is not stable with the activation of hammerhead ribozyme.

Using the formula △G^θ = nRTlnK , we can get K1. What’s more, with help of reference [24], value of Kd can also be gotten. With some rational assumption we can obtain k₁, k₁’ , k₂, k₂’. Using all the parameters[ All the parameters will be listed in the attachments. The parameter of riboswitch model will also be listed in the section of Parameter.], we can set up the model and get some important outputs.

Using the equations[ Important equations of this model will be listed. To know more about riboswitch, you may refer to the attachments.] paper mentioned, we get some result that can evaluate the efficiency of our riboswitch.

       Some important equations:

The η is the difference between high and low RFP level. The EC50 is the effective ligand concentration at which a half maximal response is achieved. The P(L=0) is basal RFP level. The P(L∞) is ligand-saturating RFP level. They are all measurements of the efficiency of the riboswitch. We often expect a riboswitch with higher η, lower EC50, lower P(L=0), higher P(L→∞).

From the table, we know MC10 has a very obvious change with the combination of the ligand, while MC31, MC7, MC3 all has a low EC50, a good symbol for low concentration toxin detection.
We also get some similar results by modeling the circuit.[ Details of circuit will be listed in the attachments.]

The result of circuit model is in accordance with results from the formula in the paper.
Using model and formula, we can test the efficiency of designed riboswitches and offer data for the model of detection system.

Using pattern formation mentioned in reference [27], we can realize quantitative measurement of toxins.

From modeling[ The concentration of ligand here is 0.01mol/L.] in the riboswitch section, we know luxI will become stable in about 2×10⁴>s and the stable concentration is about 1.05×10^-4mol/L. AHL will turn stable in about 8 hours, and the stable concentration is 0.0165mol/L. As we didn’t consider the influence of diffusion when handling sender cells, these values may be a just little different form the actual situation. However, the difference will be really little as the diffusion process is really slow. It is a reasonable estimation that speed of the AHL losing in diffusion is about 1.0×10^-10mol/L.

Then we can create a model made of a large circle with a radius of 5cm and a small circle with a radius of 0.1cm. The two circles are concentric with their center on (0 , 0). The small circle stands for the cells that produce AHL, while the large circle stands for the plate we use. AHL will continuously come out from the small circle. The space between two circles are what we will do research in.
On the boundary of small circle:

On the boundary of large circle:

In the space between the two circles:

After simulating for 4000s, we can find the concentration of AHL on the site of radius 0.3cm is 5.3×10-2μM.

Then with results from reference [27], we can conclude that there is a red ring among the place of radius 0.3cm.

So, when the concentrate of ligand is 0.01mol/L, after 4000s, there will be a red belt around the center with radius of about 0.3cm.
As the concentration of ligand increase, the speed of AHL production and degradation will increase. It is reasonable to assume that the amount of AHL going out to diffusion is always the same. In result, there will be more AHL going out as ligand increase, and the red belt will go further out after same amount of time. Besides, if the ligand concentration decreases, the red belt will go in after same amount of time.

The boundary conditions of the inner circle will change, as the ligand concentration decrease.
We may consider another situation.On the boundary of small circle:

All the other conditions are the same as the pervious situation. After 40000s, the result will be as follows. There is a red belt of radius about 0.3cm after 40000s. It is easy to conclude that if we control the time, the radius of belt will move out as the concentration of ligand increase.

Click the picture to see details.

References：
[1]. The page of Anderson promoter:
  http://parts.igem.org/Promoters/Catalog/Anderson
[2]. KENNELL D, RIEZMAN H. TRANSCRIPTION AND TRANSLATION INITIATION FREQUENCIES OF ESCHERICHIA-COLWELAC OPERON. JOURNAL
 OF MOLECULAR BIOLOGY. 1977. Volume 114. Page 1-21.
[3].iGEM Team USTC-China 2011:
  https://2011.igem.org/Team:USTC-China/Drylab/modeling
[4]. Andrzej M. Kierzek, Jolanta Zaim, and Piotr Zielenkiewicz. The Effect of Transcription and Translation Initiation 
Frequencies on the Stochastic Fluctuations in Prokaryotic Gene Expression. THE JOURNAL OF BIOLOGICAL CHEMISTRY. 2001. 
276. 8165-8172.
[5]. iGEM Team Uppsala University 2012:
   https://2012.igem.org/Team:Uppsala_University
[6].iGEM Team Tsinghua-A 2011:
  https://2011.igem.org/Team:Tsinghua-A/Modeling#References
[7].iGEM Team Nanjing-China:
  https://2013.igem.org/Team:Nanjing-China/model
[8]. iGEM Team ETH_Zurich:
  https://2013.igem.org/Team:ETH_Zurich/Parameter
[9].iGEM Team UNIPV-Pavia:
  https://2011.igem.org/Team:UNIPV-Pavia/Project/Modelling
[10]. Makoto Kobayashi, Kyosuke Nagata and Akira Ishihama. Promoter selectivity of Escherichia colWeRNA 
polymerase:effect of base substitutions in the promoter -35 region onpromoter strength. 1990. Nucleic Acids Research. 
18.
[11]. Yue Pan,Tim Durfee, Joseph Bockhorst, Mark Craven. Connecting quantitative regulatory-network models to the
 genome. 2007. Bioinformatics. 23.
[12]. Tomohiro Shimada,Yukiko Yamazaki, Kan Tanaka, Akira Ishihama. The Whole Set of Constitutive Promoters Recognized 
by RNA Polymerase RpoD Holoenzyme of Escherichia coli. 2013. PLOS.
[13]. Jeffrey H. Miller. The lacIGene: Its Role in lac Operon Control and Its Use as a Genetic System. Cold Spring 
Harbor Monograph Archive.
[14]. Marc Taraban,HonglWeZhan,Andrew E. Whitten,David B. Langley,Kathleen S. Matthews,Liskin Swint-Kruse,Jill 
Trewhella. Ligand-induced Conformational Changes and Conformational Dynamics in the Solution Structure of the Lactose
 Repressor Protein. 2008. Journal of Molecular Biology. 376. 466–481.
[15]. Liskin Swint-Kruse,HonglWeZhan, Kathleen Shive Matthews. Integrated Insights from Simulation, Experiment, and 
Mutational Analysis Yield New Details of LacIFunction. 2005. Biochemistry.
[16]. Guillermo Rodrigo, Boris Kirov, ShensWeShen, Alfonso Jaramillo. Theoretical and experimental analysis of the
 forced LacI-AraC oscillator with a minimal gene regulatory model. 2013. AIP. 23.
[17]. Denis Michel. Kinetic approaches to lactose operon induction and bimodality. 2013. Journal of Theoretical Biology. 
325. 62-75.
[18]. Sonia Covaceuszach, Giuliano Degrassi, Vittorio Venturi, Doriano Lamba. Structural Insights into a Novel 
Interkingdom Signaling Circuit by Cartography of the Ligand-Binding Sites of the Homologous Quorum Sensing LuxR-Family.
 2013. International Journal of Molecular Sciences.
[19]. iGEM Team Calgary: 
  https://2012.igem.org/Team:Calgary
[20]. Robert Sidney Cox III, Michael G Surette, Michael B Elowitz. Programming gene expression with combinatorial 
promoters. 2007. Molecular Systems Biology. 3.
[21]. Christopher M. Waters1 and Bonnie L. Bassler. The Vibrio harveyWequorum-sensing system uses shared regulatory 
components to discriminate between multiple autoinducers. 2014. Cold Spring Harbor Laboratory Press.
[22]. The experience page of PlacI:
  http://parts.igem.org/Part:BBa_R0010:Experience
[23]. Maung Nyan Win and Christina D. Smolke. A modular and extensible RNA-based gene-regulatory platform for engineering 
cellular function. 2007. PNAS. Volume 104.
[24]. Gu Kang-ding, Michael Famulok. In vitro selection of specific aptamers against microcystin-LR. 2004. Chin J Prev 
Med. Volume 38.
[25]. WOLFGANG A, PIEKEN DAVID B, OLSEN FRITZ BENSELER, HELLE AURUP, FRrrz ECKSTEIN. Kinetic Characterization of 
Ribonuclease-Resistant 2'-Modified Hammerhead Ribozymes. Science. 1991. Volume253. Pages 314-317.
[26]. Klara R BIRIKH, Paul A HEATON and Fritz ECKSTEIN. The structure, function and application of the hammerhead 
ribozyme. 1996. Eur. J. Biochem.
[27]. Subhayu Basu, Yoram Gerchman1, Cynthia H. Collins, Frances H. Arnold & Ron Weiss. A synthetic multicellular system
 for programmed pattern formation. Nature. 2005. Volume 434.
[28]. Beisel CL, Smolke CD. Design principles for riboswitch function. 2009. PLoS Comput Biol.