Team:Slovenia/PROJECT/oscillator/details

From 2010.igem.org

(Difference between revisions)
(New page: {{Team:Slovenia/templatetop}} {{Team:Slovenia/templateproject}} <html> <style> #subgumb1{display:none}#subgumb2{display:none}#subgumb2a{display:none}#subgumb2b{display:none}#subgumb3{disp...)
 
(3 intermediate revisions not shown)
Line 14: Line 14:
</style>
</style>
<div id="overhead">
<div id="overhead">
-
<div id="naslov">oscillator</div>
+
<div id="naslov">oscillator - details</div>
</div>
</div>
<div id="besedilo">
<div id="besedilo">
Line 23: Line 23:
__TOC__
__TOC__
-
<div id="okvircek" style="background-color:#CECDCD;border:1px solid #AEABAE;margin-top:9px;padding:11px">
+
<h2>MATLAB script and how it works</h2>
-
<html>
+
 
-
<h3>State of the art:</h3>
+
In order to perform the stochastic simulations stochastic simulation algorithm (SSA) is commonly used . Our MATLAB script file is an implementation of the SSA algorithm for MATLAB environment. The basic algorithm we used is as follows:<br>
<ul>
<ul>
-
<li>Genetic oscillators have been experimentally demonstrated and extensively modeled but their number remains very small</li>
+
<li>1. Initialization: Initialize the number and initial concentrations of species in the system, reactions used, reactions constants and kinetics.</li>
-
<li>Problem with wider application of genetic oscillators is that there  is a very small number of available parts with matched properties</li>
+
<li>2. Execution step: Generate random numbers to determine the next reaction to occur. The probability of a given reaction to be chosen is proportional to the concentration of reagents (some of the algorithms also determine the time step, time step we used was constant but small, to make the execution step efficient - but at the cost of simulation taking longer to complete).</li>
-
<li>Selection of oscillator properties is often fortuitous and time constants are adjustable only in a narrow range</li>
+
<li>3. Update: Update the species concentrations based on the reaction chosen in step 2.</li>
-
<li>Construction of each new genetic oscillator is laborious due to the difficulties of matching the properties of its constituents</li>
+
<li>4. Iterate: Go back to Step 2 unless the simulation time has been exceeded.</li>
</ul>
</ul>
-
 
<br>
<br>
 +
Note that using constant volume concentration can also mean number of molecules.
 +
 +
<h3>Main variables used in our MATLAB script are as follows:</h3>
-
<h3>Aims:</h3>
 
<ul>
<ul>
-
<li>Define the universal type of parts for the construction of  genetic oscillators that would be available in large number of variants  with matched properties</li>
+
<li>S is the number of species</li>
-
<li>Identify the orthogonal system that could function in prokaryotic as well as in eukaryotic cells</li>
+
<li>R is the number of reactions</li>
-
<li>Explore the new potential features of such oscillators</li>
+
<li>v(S,R) is a matrix, called stoichiometric matrix; each row represents one reaction, while the columns are the species; the numbers in the row tell us what species is used and which one produced (and how much of it) in one reaction</li>
-
<li>Demonstrate the functionality of basic building blocks</li>
+
<li>k is a matrix of kinetic rates; rows are different genes (ie gene for NICTAL/RNA polymerase/GFP), while the columns are kinetic rates used in different reactions</li>
 +
<li>x is a vector of concentration (or molecule count) of each species</li>
 +
<li>t vector of change - calculated in every step for every reaction</li>
</ul>
</ul>
<br>
<br>
-
<h3>Conclusions:</h3>
+
Each step of the simulation (for i =1:N loop) first evaluates the probabilities of all reactions (propensity function ajx), again proportional to the concentration (or molecule count) of the species and kinetic rates.
-
<ul>
+
<br>
-
<li>Zinc fingers or other designed DNA-binding proteins (e.g. TAL  effectors) provide the large repertoire&nbsp; (hundreds) of readily available  artificial repressors with very similar properties</li>
+
Change of the system state is calculated. This is where the randomness occurs. For every reaction (for j = 1:R loop) we calculate t as t = v(:,j)'*poissrnd(a(j)*tau), where v(:,j) presents the stoichiometric vector for reaction j, a(j) presents the result of the propensity function for this reaction and poissrnd(x) a random number somewhere close to x according to Poisson distribution.<br>
-
<li>T7 RNA polymerase provides the orthogonal platform for oscillators for both prokaryotic and eukaryotic systems</li>
+
-
<li>Increasing the number of repressors in the repressilator topology  increases the stability and provides oscillators with desired time  constant</li>
+
-
<li>Addition of a low constitutive level of transcription of T7 RNA  polymerase in a Smolen oscillator initiates oscillations even in the deterministic system without of having to rely on stochastic  fluctuations</li>
+
-
<li>We demonstrated that only repressilators with odd-number of repressors result in sustained oscillations</li>
+
-
<li>Even monomeric repressor can lead to oscillations, i.e. cooperativity is not essential</li>
+
-
</ul>
+
-
</html>
+
-
</div>
+
-
<h2>Background and overview</h2>
+
The validation of calculated change needs to be examined. If it causes negative values of concentrations (x), we do not take it as valid. In other case we simply add vector t to the vector change from previous reaction. Once we finish the calculation for all reactions, we add the change vector to the species vector x, updating the values to reflect reaction execution.<br>
-
An oscillator is a key device in any sequential logic circuit. Its repetitive signal serves for synchronization of events and therefore insuring the correct results of any logical or control operation inside a digital device. While this digitally-accurate concept cannot be directly transferred to biological systems, the nature of this device can still be used to our benefit. As we strive to control processes inside a single or multiple cells this tool could offer help in triggering the selected events in any given order or at any given moment. For ensuring that kind of control, we would need a biological oscillator that would be self-sustaining and would have adjustable frequency.<br><br>
+
If we want to know how the concentration (or molecule count) was changing throughout the simulation, we need to save the x value in every step. The code under %% Plotting results simply turns the saved values of the simulations into a graph comprehensible by the person sitting behind the computer.
-
<h2>Modeling</h2>
+
<h2>MATLAB Simbiology tool and how it works</h2>
-
While experimental work is what counts at the end of the day, modeling of any biological system can help with roughly determining what experimental parameters should be used to achieve the best possible results.<br>
+
To import a simple example of our model in sbml format we could follow these simple steps:
-
We created two different models of a biological oscillator for that purpose, i.e. Smolen-type oscillator (Smolen et al., 1998) and represillator (Elowitz et al., 2000). Both models have been experimentally tested and extensively modeled based on respective biological processes. Our aim was to investigate in which manner the use of ZNF-based panel of repressors could extend the current experimentally realizable set of oscillators and what additional properties could be achieved. We investigated the validity of some premises of genetic oscillators, such as the requirement for the oligomeric repressors to support nonlinear binding profiles. Additionally, we performed stochastic and deterministic simulations.<br><br>
+
-
<h2>How is deterministic approximation different from stochastic?</h2>
+
<ul>
 +
<li>open the simbiology tool by writting "simbiology" in the Matlab console</li>
 +
<li>create a new project from menu File ( File -&gt; New Project )</li>
 +
<li>in the new project wizard window in the left panel menu select  "Add model"</li>
 +
<li>in the right panel "Add model", using scrolled bar, select item "Load model from file"</li>
 +
<li>select the source folder where the sbml file is located</li>
 +
<li>click "Finish" to build selected model</li>
 +
</ul>
 +
<br><br>
 +
Several variable changes of the existing default configuration parameters of simbiology are required to set-up configuration for simulations. The following settings should be used:<br><br>
-
Deterministic approximation is a mathematical model that follows certain rules and equations without any randomness. It is determined by one or usually more differential equations which produce the same output for a given input every time we run the simulation, meaning that the results for of the same starting parameters will always be the same. Because deterministic system is determined by ordinary differential equations the calculation of the results is actually just a problem of solving this system, which was solved in our case with MATLAB, using the SimBiology toolkit.<br><br>
+
<ul>
 +
<li>in the "Project explorer", expand the "Model session" tree</li>
 +
<li>click on "Configuration Settings"</li>
 +
<li>set the simulation time to 30.000 (recommended) or more</li>
 +
<li>set the solver type to "ode15s"</li>
 +
<li>click on the "Data logging" tab</li>
 +
<li>deselect all species logs you do not want in your result figure</li>
 +
<li>select the species you want in your result figure</li>
 +
</ul><br><br>
-
On the other hand, the key feature of a stochastic process is certain level of randomness. The evolution of stochastic system is usually determined by a probability distribution. Regardless of the distribution used it is highly unlikely that two results from the same input (initial parameters) would be exactly the same.<br><br>
+
<h2>SGN Sim and how to create the input file</h2>
-
All our stochastic models used Poisson distribution for executing reactions at each step of the simulation, which produces slightly different results every time the simulation is run even if the input is the same, making them more natural and representative. In a complex system such as biological oscillator it would be too complex (or even impossible) to take into account every variable present, however using a probability distribution for some of the processes is a good enough approximation of the natural situation. While many different probability distributions are known, the Poisson distribution is the most commonly used distribution when it comes to molecular interaction.<br>
+
SGN Sim is easy to use and versatile program that offers a great deal of tools for stochastic simulation. Its base is SSA algorithm described in the 'MATLAB script' section, but with certain additions. For example, we can delay the production of translation with a simple fix in the input file. Two great guides come along with the program explaining how to create an input file, what it should contain and how to use certain tricks. You can find the guides <html><a href="http://www.cs.tut.fi/~sanchesr/SGN/SGNSim.html">here</a></html>. The program then parses your input file and runs the simulation. Here is a quick recap of what we used in our input file:
-
For the calculation of discrete stochastic simulation the Gillespie algorithm (stochastic simulation algorithm - SSA) was used (Gillespie D. T. 2007 Stochastic simulation of chemical kinetics). Two different tools were used for stochastic model, a script written for MATLAB and a freeware program SGNSim (runs on Microsoft windows), that serves for simulation of gene networks ("SGN Sim, a Stochastic Genetic Networks Simulator", A.S.Ribeiro, J.Lloyd-Price, 2007; <html><a href="http://www.cs.tut.fi/~sanchesr/SGN/SGNSim.html">http://www.cs.tut.fi/~sanchesr/SGN/SGNSim.html</a></html> ). MATLAB script used a version of SSA algorithm, while SGNSim uses a delayed SSA (SSA that includes delays for transcription and translation).<br><br>
+
<ul>
-
 
+
<li>"time" and "stop _time" parameters determine when the simulation starts and ends, the "read_out interval" tells the program the saving frequency of calculation results and "output_file" tells the program what output file to use (it can save the results directly as an Excel file, a txt file, etc.).</li>
-
A distinction was made in simulations of prokaryotic (bacterial) and eukaryotic (mammalian) systems. The main difference is the delay between transcription and translation. While this process is coupled in bacteria there is a large delay in eukaryotic cells due to the RNA maturation and translocation&nbsp; from the nucleus into the cytoplasm.<br><br>
+
<li>Segment:
 +
<br>
 +
lua !{
 +
...
 +
}
 +
<br>
 +
tells the program what are the constants used in our simulation.</li>
 +
<li>Segment:<br>
 +
population {
 +
...
 +
}
 +
<br>determines starting concentrations (or molecule counts) for all species.</li>
 +
<li>Segment:<br>
 +
reaction {
 +
...
 +
}
 +
<br>
 +
is the part where the reactions are determined.</li>
 +
</ul><br>
 +
To describe a reaction we simply use the notation from population { } segment and follow the syntax provided by the two guides:
 +
A + B --[k]--&gt; C;
-
Deterministic models seem to be more robust. Wide variety of starting parameter values were simulated and most of the results were an oscillating system. On the other hand, stochastic system quickly reaches an equilibrium state if we do not use realistic parameter values. That is why deterministic model was used to compare the results with stochastic simulations, however extensive testing, including difference in prokaryotic and eukaryotic cells, were only done with stochastic model.<br><br>
+
where A, B and C are the reagents and k is the kinetic constant used for this reaction. If we use a letter for the rate, we must declare it in the lua !{ } section, however the parser allows using a number as well (A + B --[0.02]--&gt; C). Degradation is simply noted as:
-
Both script file for MATLAB and input file for SGNSim were written by Slovenian iGEM team members.<br><br><br><br>
+
A --[k_deg]--&gt; ;
-
[1] Paul Smolen, Douglas A. Baxter and John H. Byrne, Department of Neurobiology and Anatomy, University of Texas, Medical School, Houston, Texas, 1998. Frequency selectivity, multistability, and oscillations emerge from models of genetic regulatory system.<br>
+
The most important part is the time delay we can add to products of the reactions. If I use the example from our input file, where we had to delay the production of translation:
-
Am J Physiol Cell Physiol 274: 531-542.<br><br>
+
-
[2] Michael B. Elowitz, Stanislas Leibler, Departments of Molecular Biology and Physics, Princeton University, USA, 2000. A synthetic oscillatory network of transcriptional regulators.<br>
+
A --[transcription_rate]--&gt; B(translation_delay);
-
Nature Vol 403: 335-338<br>
+
 +
We can again replace translation_delay by a number, the parser even allows mathematical expressions (such as 700/200 - 700 divided by 200). A great feature of this program is also that you can use a probability distribution inside the delay. We used exponential distribution for the translation delay:
 +
A--[transcription_rate]--&gt;B(exp:(1/translation_delay));
 +
Simple, yet very effective and useful.<br><br>
 +
As we can see setting up a biological model and running a stochastic simulation is easy and very accurate with this program. We were also surprised how fast it calculates the result of simulation, compared to our MATLAB script.
 +
<br><br>
 +
<html><a href="https://2010.igem.org/Team:Slovenia/PROJECT/oscillator">Back to Oscillator.</a></html>
<!--STOP BESEDILO-->
<!--STOP BESEDILO-->

Latest revision as of 20:21, 27 October 2010

Chuck Norris facts:

oscillator - details


Contents


MATLAB script and how it works

In order to perform the stochastic simulations stochastic simulation algorithm (SSA) is commonly used . Our MATLAB script file is an implementation of the SSA algorithm for MATLAB environment. The basic algorithm we used is as follows:

  • 1. Initialization: Initialize the number and initial concentrations of species in the system, reactions used, reactions constants and kinetics.
  • 2. Execution step: Generate random numbers to determine the next reaction to occur. The probability of a given reaction to be chosen is proportional to the concentration of reagents (some of the algorithms also determine the time step, time step we used was constant but small, to make the execution step efficient - but at the cost of simulation taking longer to complete).
  • 3. Update: Update the species concentrations based on the reaction chosen in step 2.
  • 4. Iterate: Go back to Step 2 unless the simulation time has been exceeded.


Note that using constant volume concentration can also mean number of molecules.

Main variables used in our MATLAB script are as follows:

  • S is the number of species
  • R is the number of reactions
  • v(S,R) is a matrix, called stoichiometric matrix; each row represents one reaction, while the columns are the species; the numbers in the row tell us what species is used and which one produced (and how much of it) in one reaction
  • k is a matrix of kinetic rates; rows are different genes (ie gene for NICTAL/RNA polymerase/GFP), while the columns are kinetic rates used in different reactions
  • x is a vector of concentration (or molecule count) of each species
  • t vector of change - calculated in every step for every reaction


Each step of the simulation (for i =1:N loop) first evaluates the probabilities of all reactions (propensity function ajx), again proportional to the concentration (or molecule count) of the species and kinetic rates.
Change of the system state is calculated. This is where the randomness occurs. For every reaction (for j = 1:R loop) we calculate t as t = v(:,j)'*poissrnd(a(j)*tau), where v(:,j) presents the stoichiometric vector for reaction j, a(j) presents the result of the propensity function for this reaction and poissrnd(x) a random number somewhere close to x according to Poisson distribution.

The validation of calculated change needs to be examined. If it causes negative values of concentrations (x), we do not take it as valid. In other case we simply add vector t to the vector change from previous reaction. Once we finish the calculation for all reactions, we add the change vector to the species vector x, updating the values to reflect reaction execution.

If we want to know how the concentration (or molecule count) was changing throughout the simulation, we need to save the x value in every step. The code under %% Plotting results simply turns the saved values of the simulations into a graph comprehensible by the person sitting behind the computer.

MATLAB Simbiology tool and how it works

To import a simple example of our model in sbml format we could follow these simple steps:

  • open the simbiology tool by writting "simbiology" in the Matlab console
  • create a new project from menu File ( File -> New Project )
  • in the new project wizard window in the left panel menu select "Add model"
  • in the right panel "Add model", using scrolled bar, select item "Load model from file"
  • select the source folder where the sbml file is located
  • click "Finish" to build selected model



Several variable changes of the existing default configuration parameters of simbiology are required to set-up configuration for simulations. The following settings should be used:

  • in the "Project explorer", expand the "Model session" tree
  • click on "Configuration Settings"
  • set the simulation time to 30.000 (recommended) or more
  • set the solver type to "ode15s"
  • click on the "Data logging" tab
  • deselect all species logs you do not want in your result figure
  • select the species you want in your result figure


SGN Sim and how to create the input file

SGN Sim is easy to use and versatile program that offers a great deal of tools for stochastic simulation. Its base is SSA algorithm described in the 'MATLAB script' section, but with certain additions. For example, we can delay the production of translation with a simple fix in the input file. Two great guides come along with the program explaining how to create an input file, what it should contain and how to use certain tricks. You can find the guides here. The program then parses your input file and runs the simulation. Here is a quick recap of what we used in our input file:

  • "time" and "stop _time" parameters determine when the simulation starts and ends, the "read_out interval" tells the program the saving frequency of calculation results and "output_file" tells the program what output file to use (it can save the results directly as an Excel file, a txt file, etc.).
  • Segment:
    lua !{ ... }
    tells the program what are the constants used in our simulation.
  • Segment:
    population { ... }
    determines starting concentrations (or molecule counts) for all species.
  • Segment:
    reaction { ... }
    is the part where the reactions are determined.

To describe a reaction we simply use the notation from population { } segment and follow the syntax provided by the two guides: A + B --[k]--> C;

where A, B and C are the reagents and k is the kinetic constant used for this reaction. If we use a letter for the rate, we must declare it in the lua !{ } section, however the parser allows using a number as well (A + B --[0.02]--> C). Degradation is simply noted as:

A --[k_deg]--> ;

The most important part is the time delay we can add to products of the reactions. If I use the example from our input file, where we had to delay the production of translation:

A --[transcription_rate]--> B(translation_delay);

We can again replace translation_delay by a number, the parser even allows mathematical expressions (such as 700/200 - 700 divided by 200). A great feature of this program is also that you can use a probability distribution inside the delay. We used exponential distribution for the translation delay:

A--[transcription_rate]-->B(exp:(1/translation_delay));

Simple, yet very effective and useful.

As we can see setting up a biological model and running a stochastic simulation is easy and very accurate with this program. We were also surprised how fast it calculates the result of simulation, compared to our MATLAB script.

Back to Oscillator.


Retrieved from "http://2010.igem.org/Team:Slovenia/PROJECT/oscillator/details"