Team:Grenoble/Projet/Modelling/Stochastic

From 2011.igem.org

(Difference between revisions)
Line 66: Line 66:
and reactions happen randomly in the medium. On a macroscopic scale, the reactions can be seen  
and reactions happen randomly in the medium. On a macroscopic scale, the reactions can be seen  
as deterministic, and the statistical properties are summarized by constants in classical ODEs.  
as deterministic, and the statistical properties are summarized by constants in classical ODEs.  
 +
 +
      </p>
 +
      <p>
 +
<img src="https://static.igem.org/mediawiki/2011/5/5f/Mouvement_brownien.png" class="centerwide" style="height: 300px; width: auto;"/>
 +
      </p>
 +
      <p>
However, at a cell's scale, the influence of the randomness of the reactions is no longer  
However, at a cell's scale, the influence of the randomness of the reactions is no longer  
negligible, especially for biosensors. For an efficient measurement biosensor systems must provide  
negligible, especially for biosensors. For an efficient measurement biosensor systems must provide  
Line 71: Line 77:
      </p>
      </p>
      <p>
      <p>
-
IMAGEEUH
+
In his algorithm Gillespie uses propensity theory to describe the behaviour of such a medium. Each
 +
reaction occuring in the cell, like construction or destruction of a protein, has a certain
 +
propensity (see it as the probability that the event occurs in a coming amount of time).
      </p>
      </p>
 +
<img src="https://static.igem.org/mediawiki/2011/2/2b/Stochastic_illustration.png" class="centerwide" style="height: 300px; width: auto;"/>
      <p>
      <p>
-
Each reaction has a certain probability to occur in an interval of time, and this probability
+
The propensities are calculated at each time step. Then, one random number is generated and its value
-
depends, for example, on the concentration of reactant species in the medium and on the affinity
+
will set the reaction that will occur for the current time step. Of course, a reaction with a relatively high
-
constants of the reactions.
+
propensity will have more chance to occur than another one.
      </p>
      </p>
      <p>
      <p>
-
In his algorithm Gillespie uses propensity theory to describe the behaviour of such a medium. Each
+
After this step the concentrations are changed according to the reaction, and the same process is repeated until
-
reaction occuring in the cell, like construction or destruction of a protein, has a certain
+
the ending time of the simulation is reached.
-
propensity. These propensities depend on the reactants concentrations and on other molecules in the  
+
-
cell (e.g. : The production of a protein from a gene placed after a pLac promoter depends on the  
+
-
concentrations of lacI and/or IPTG molecules in the cell).
+
      </p>
      </p>
 +
     
 +
<p><strong>Reference :</strong> Daniel T. Gillespie (1977). "Exact Stochastic Simulation of Coupled Chemical Reactions".
 +
The Journal of Physical Chemistry 81 (25): 2340–2361</p>
 +
<p>Daniel T. Gillespie (1976). "A General Method for Numerically Simulating the Stochastic Time Evolution
 +
of Coupled Chemical Reactions". Journal of Computational Physics 22 (4): 403–434</p>
 +
 +
      <h2>Our genetical network</h2>
 +
 
      <ol>
      <ol>
 +
      <strong>
<li>Propensity functions and the Gillespie Algorithm</li>
<li>Propensity functions and the Gillespie Algorithm</li>
 +
      </strong>
<p>
<p>
  Our stochastic model only describes the stochastic behaviour of the Toggle switch genetical network.  
  Our stochastic model only describes the stochastic behaviour of the Toggle switch genetical network.  
Line 92: Line 108:
  the precision, the behaviour and the limits of our system.
  the precision, the behaviour and the limits of our system.
</p>
</p>
 +
<p>
<p>
  The propensity functions used in our models are derived from the ODEs we have already written for  
  The propensity functions used in our models are derived from the ODEs we have already written for  
Line 102: Line 119:
</tr>
</tr>
<tr>
<tr>
-
 
-
 
-
 
  <td>
  <td>
-
    <math >
+
    $ \phi \arrow TetR $
-
      <mrow>
+
-
<mi>&Phi;</mi>
+
-
<mo>&rarr;</mo>
+
-
<mi>TetR</mi>
+
-
      </mrow>
+
-
    </math>
+
  </td>
  </td>
-
 
-
 
-
 
  <td>
  <td>
-
  <math >
+
    $ \frac{k_{pLac}P_{Lac total}}{1+(\frac{LacI}{1 + \frac{IPTG}{K_{LacI - IPTG}}})^{n_{n_pLac}}} $
-
  <mrow>
+
  </td>
-
  <mfrac>
+
</tr>
-
  <mrow>
+
<tr>
-
    <msub>
+
  <td>
-
    <mi>k</mi>
+
    $ TetR \arrow \phi$
-
  <mi>plac</mi>
+
  </td>
-
    </msub>
+
  <td>
-
   
+
    $ \delta_{TetR} TetR $
-
      <mrow>
+
  </td>
-
<msub>
+
</tr>
-
<mi>P</mi>
+
<tr>
-
      <mi>lac total</mi>
+
  <td>
-
</msub>
+
    $ \phi \arrow LacI $
-
      </mrow>
+
  <td>
-
   
+
    $ \frac{k_{pTet}P_{Tet total}}{1+(\frac{TetR}{1 + \frac{aTc}{K_{TetR - aTc}}})^{n_{n_pTet}}} $
-
  </mrow>
+
  </td>
-
  <mrow>
+
</tr>
-
    <mi>1 +</mi>
+
<tr>
-
    <msup>
+
  <td>
-
  <mrow>
+
    $ LacI \arrow \phi$
-
    <mfenced open="(" close=")" separators=",">
+
  </td>
-
      <mrow>
+
  <td>
-
<mfrac>
+
    $ \delta_{LacI} LacI $
-
  <mrow>
+
  </td>
-
   
+
</tr>
-
      <mrow>
+
 
-
<mi>lacI</mi>
+
-
      </mrow>
+
-
   
+
-
  </mrow>
+
-
  <mrow>
+
-
    <mi>1 +</mi>
+
-
    <mfrac>
+
-
      <mrow>
+
-
+
-
  <mrow>
+
-
    <mi>IPTG</mi>
+
-
  </mrow>
+
-
+
-
      </mrow>
+
-
      <mrow>
+
-
<msub>
+
-
<mi>K</mi>
+
-
      <mi>lacI - IPTG</mi>
+
-
</msub>
+
-
      </mrow>
+
-
    </mfrac>
+
-
  </mrow>
+
-
</mfrac>
+
-
      </mrow>
+
-
    </mfenced>
+
-
  </mrow>
+
-
  <msub>
+
-
  <mi>n</mi>
+
-
<mi>plac</mi>
+
-
  </msub>
+
-
    </msup>
+
-
  </mrow>
+
-
</mfrac>
+
-
      </mrow>
+
-
    </math>
+
-
    </td>
+
-
 
+
-
 
+
-
  </tr>
+
-
  <tr>
+
-
    <td>
+
-
      <math >
+
-
      <mrow>
+
-
<mi>TetR</mi>
+
-
<mo>&rarr;</mo>
+
-
<mi>&Phi;</mi>
+
-
      </mrow>
+
-
      </math>
+
-
    </td>
+
-
     
+
-
   
+
-
    <td>
+
-
      <math >
+
-
      <mrow>
+
-
<msub>
+
-
<mi>&delta;</mi>
+
-
      <mi>TetR</mi>
+
-
</msub>
+
-
<mi>TetR</mi>
+
-
      </mrow>
+
-
      </math>
+
-
    </td>
+
-
 
+
-
 
+
-
    </tr>
+
-
    <tr>
+
-
      <td>
+
-
<math >
+
-
  <mrow>
+
-
    <mi>&Phi;</mi>
+
-
    <mo>&rarr;</mo>
+
-
    <mi>lacI</mi>
+
-
  </mrow>
+
-
</math>
+
-
      </td>
+
-
 
+
-
 
+
-
 
+
-
      <td>
+
-
<math >
+
-
  <mrow>
+
-
    <mfrac>
+
-
      <mrow>
+
-
<msub>
+
-
<mi>k</mi>
+
-
      <mi>pTet</mi>
+
-
</msub>
+
-
+
-
  <mrow>
+
-
    <msub>
+
-
    <mi>P</mi>
+
-
  <mi>Tet total</mi>
+
-
    </msub>
+
-
  </mrow>
+
-
+
-
      </mrow>
+
-
      <mrow>
+
-
<mi>1 +</mi>
+
-
<msup>
+
-
      <mrow>
+
-
<mfenced open="(" close=")" separators=",">
+
-
  <mrow>
+
-
    <mfrac>
+
-
      <mrow>
+
-
+
-
  <mrow>
+
-
    <mi>TetR</mi>
+
-
  </mrow>
+
-
+
-
      </mrow>
+
-
      <mrow>
+
-
<mi>1 +</mi>
+
-
<mfrac>
+
-
  <mrow>
+
-
   
+
-
      <mrow>
+
-
<mi>aTc</mi>
+
-
      </mrow>
+
-
   
+
-
  </mrow>
+
-
  <mrow>
+
-
    <msub>
+
-
    <mi>K</mi>
+
-
  <mi>TetR - aTc</mi>
+
-
    </msub>
+
-
  </mrow>
+
-
</mfrac>
+
-
      </mrow>
+
-
    </mfrac>
+
-
  </mrow>
+
-
</mfenced>
+
-
      </mrow>
+
-
      <msub>
+
-
      <mi>n</mi>
+
-
    <mi>pTet</mi>
+
-
      </msub>
+
-
</msup>
+
-
      </mrow>
+
-
    </mfrac>
+
-
  </mrow>
+
-
</math>
+
-
      </td>
+
-
     
+
-
 
+
-
 
+
-
    </tr>
+
-
    <td>
+
-
      <math >
+
-
      <mrow>
+
-
<mi>lacI</mi>
+
-
<mo>&rarr;</mo>
+
-
<mi>&Phi;</mi>
+
-
      </mrow>
+
-
      </math>
+
-
    </td>
+
-
     
+
-
   
+
-
    <td>
+
-
      <math >
+
-
      <mrow>
+
-
<msub>
+
-
<mi>&delta;</mi>
+
-
      <mi>lacI</mi>
+
-
</msub>
+
-
<mi>lacI</mi>
+
-
      </mrow>
+
-
      </math>
+
-
    </td>
+
-
+
-
   
+
</table>
</table>
<p>
<p>
Line 321: Line 156:
  In our Matlab code the propensities are computed at each time step in the file Stochastic_model.m.
  In our Matlab code the propensities are computed at each time step in the file Stochastic_model.m.
</p>
</p>
-
<p>
+
-
  Then a pseudo-random number is generated in the interval [0;1] with Matlab function rand(). This
+
<strong>
-
  number will set which one of the reactions will occur during this iteration, or time step. The interval
+
-
  of time that separates each iteration is set by the sum of all propensities. If all reactions have very
+
-
  high propensities, many iterations will happen in a certain amount of time – And many iterations means
+
-
  many reactions as well. By contrary, if all propensities are low, few reactions will occur during the same
+
-
  amount of time. The time interval is then inversely proportionnal to the sum of all propensities.
+
-
</p>
+
-
<p>
+
-
  This is the core of the algorithm. It can be found in the file Gillespie.m
+
-
  The species concentrations are finally changed according to the reaction fixed
+
-
  by the random number and the loop keeps running until the ending time of the
+
-
  simulation is reached (ending time fixed by the user).
+
-
</p>
+
-
<p>Reference : Daniel T. Gillespie (1977). "Exact Stochastic Simulation of Coupled Chemical Reactions".
+
-
The Journal of Physical Chemistry 81 (25): 2340–2361</p>
+
-
<p>Daniel T. Gillespie (1976). "A General Method for Numerically Simulating the Stochastic Time Evolution
+
-
of Coupled Chemical Reactions". Journal of Computational Physics 22 (4): 403–434</p>
+
<li>Runs and statistical properties</li>
<li>Runs and statistical properties</li>
 +
</strong>
<p>
<p>
  In a Gillespie simulation, the information brought by a few instances of the Gillespie simulation is  
  In a Gillespie simulation, the information brought by a few instances of the Gillespie simulation is  

Revision as of 15:21, 25 October 2011

Grenoble 2011, Mercuro-Coli iGEM


Modelling - Stochastic

With precedent modelling, we get the deterministic behavior of our system. All the parameters are precisely known and the solution obtained is always the same whatever the number of simulation is. However, in the switch area, the choice between one of the two state is randomly made by bacterias. That's why, we need a modelling which is taking into account the randomness of the choice of the state.

Sensitivity to noise

In a first time, we just used the deterministic model. The gradient of IPTG and the homogenous concentration of aTc are modeled by a normal distribution with for standard deviation a predetermined percentage of the distribution average.


Figure 1: Logarithmic gradient of IPTG and aTc repartition on the plate with deterministic modelling
Figure 2: Logarithmic gradient of IPTG and aTc repartition on the plate with random deterministic modelling
Figure 3: Concentration of both repressor on the plate with deterministic modelling.
Figure 4: Concentration of both repressor on the plate with random deterministic modelling.

With the random model, fluctuations in the system are more present. This affects the quality of the switch. That's why is important to take into account the stochastic aspect (random and probalistic studies) of the system. In our stochastic model, we used gillespie algorithm. An often used one in stochastic modelling.

Gillespie algorithm

During a chemical reaction, the molecules move at random in the medium obeying brownian motion, and reactions happen randomly in the medium. On a macroscopic scale, the reactions can be seen as deterministic, and the statistical properties are summarized by constants in classical ODEs.

However, at a cell's scale, the influence of the randomness of the reactions is no longer negligible, especially for biosensors. For an efficient measurement biosensor systems must provide the expected precision of the measure.

In his algorithm Gillespie uses propensity theory to describe the behaviour of such a medium. Each reaction occuring in the cell, like construction or destruction of a protein, has a certain propensity (see it as the probability that the event occurs in a coming amount of time).

The propensities are calculated at each time step. Then, one random number is generated and its value will set the reaction that will occur for the current time step. Of course, a reaction with a relatively high propensity will have more chance to occur than another one.

After this step the concentrations are changed according to the reaction, and the same process is repeated until the ending time of the simulation is reached.

Reference : Daniel T. Gillespie (1977). "Exact Stochastic Simulation of Coupled Chemical Reactions". The Journal of Physical Chemistry 81 (25): 2340–2361

Daniel T. Gillespie (1976). "A General Method for Numerically Simulating the Stochastic Time Evolution of Coupled Chemical Reactions". Journal of Computational Physics 22 (4): 403–434

Our genetical network

  1. Propensity functions and the Gillespie Algorithm
  2. Our stochastic model only describes the stochastic behaviour of the Toggle switch genetical network. The toggle switch is the core of our system, it is the most sensitive part of the network and sets the precision, the behaviour and the limits of our system.

    The propensity functions used in our models are derived from the ODEs we have already written for deterministic modelling :

    Chemical reaction Propensity
    $ \phi \arrow TetR $ $ \frac{k_{pLac}P_{Lac total}}{1+(\frac{LacI}{1 + \frac{IPTG}{K_{LacI - IPTG}}})^{n_{n_pLac}}} $
    $ TetR \arrow \phi$ $ \delta_{TetR} TetR $
    $ \phi \arrow LacI $ $ \frac{k_{pTet}P_{Tet total}}{1+(\frac{TetR}{1 + \frac{aTc}{K_{TetR - aTc}}})^{n_{n_pTet}}} $
    $ LacI \arrow \phi$ $ \delta_{LacI} LacI $

    The parameters are the same as those used in ODEs and can be found on the parameters page. In our Matlab code the propensities are computed at each time step in the file Stochastic_model.m.

  3. Runs and statistical properties
  4. In a Gillespie simulation, the information brought by a few instances of the Gillespie simulation is enough to get the general behaviour of the genetical system. However, it is not sufficient when the information that we want is the mean or the variance of the concentrations in each species.

    In this case a great number of runs is necessary to have a correct estimate of the expected values. This is why we had to write a Matlab code to iterate a great number of runs. This part of the code can be found in the file Main_gillespie.m.

    Once the runs are computed through Gillespie algorithm, we have to extract the information from them. In this purpose we wrote the Hist.m, test.m and Dynamicdistros.m files. These files are specific to our system, we hope they can give an idea of how to analyse the results obtained via our Gillespie code, but they are not as easily understandable as other files.

    Once the datasets are obtained we have to extract its statistical properties. Refer to the next section for more information.

IMPORTANT NOTE:

We tried to write a Matlab code that is as easily adaptable to any other system as possible. However, because of the lack of time and the great amount of work it requires, we could not build a completely generic MATLAB function handler for Gillespie simulations. We provide the source codes here of the Matlab stochastic scripts for our simulations and tried to comment them as much as possible. Note that, if you want to adapt our code to a completely different system, only the Stochastic_model.m and parameters.mat files need to be changed, but a good understanding of the whole code is necessary.

Mean, standard deviation and statistical properties

To get the number of bacteria and the minimal step for the IPTG gradient we get, we need to compute the mean and the standard deviation of each of the two species of the toggle switch.

X1 and X2 are here the matrices which will represent the two toggle switch ways in each bacterium on the whole plate. On eahc point of the plate are wells containing a great number of bacteria. Each way in each bacterium is a random variable. The X matrices are then matrices of nbcells*nbcell/well random variables.

Thanks to stochastic modelling we can obtain the mean and variance in each of the nbwells wells on the plate.

To design our final device we need to know the width of the interface between the two ways of the toggle switch. We also need to know the number of bacteria needed in the wells to have a proper measurement.

The interface which will be the colored part of our plate will turn red when populations in way 1 and populations in way 2 are in presence on the same point on the plate. We then need to know the statistical properties of the (X1X2)well random variable.

X1X2(well) is of course not continuous but discrete, we just want to highlight the deviation problem caused by σX1X2(well))

We want to know µX1X2(well) and σX1X2(well) to obtain respectively :

  1. The width of the "gaussian" function of µX1X2(well) to set the minimal definition (the ΔIPTG between wells) of our final device
  2. The minimum number of bacteria we want in the wells.
  3. We have nbcell/well independant random variables with the same probability density function. According to central limit theorem, the mean of X1X2 = (X1X2 cell1 + X1X2 cell2 + ... + X1X2 celln) / n is µX1X2(well) and its standard deviation is σX1X2(well)/n. The width of the gaussian is therefore easily calculable (µX1X2(well) = µX1(well) + µX2(well)), but it's not that easy for the standard deviation σX1X2(well)

    If we consider Y1 and Y2 two random variables correlated, with known mean and variance (µ1, µ2, σ1, σ2)
    Var(Y1Y2) = E[(Y1Y2 - E[(Y1Y2])2]
    = E[(Y1Y2)2 -2E[Y1Y2]Y1Y2 + E[Y1Y2]2]
    = E[(Y1Y2)] -2(µ1µ2)2 + µ1µ2

    and E[(Y1Y2)] is not reductible to function of µ1 and µ2. To obtian this value we then needed to create a composite random variable X3 wich will be calculated during for each run of our MATLAB stochastic algorithm (see previous section).
    x3 = x1 * x2
    We can thus get the variance and mean of X1X2 (X3) through simulation. If we want to get an error inferior to 10% for example around the µX1X2(well) curve, the number of bacteria n needed will be n so that : σ X 3 well n / µ X1X2 well < 10% With such a precision we can then calculate the IPTG definition between the wells.