Modified simple cell model (補完

We referenced and modified the model made by 2007 imperial iGEM team to our needs. This model simulate basic behavior of single cell, and assumptions are as follows:

1. All molecules, including proteins, protein complexes and small molecules are uniformly distributed in cell bodies.
2.  The diffusion rate constant $k16$ of AHL is determined only by the AHL gradient between cytoplasm (denoted as $[AHL_i]$) and extracellular matrix (denoted as $[AHL_e]$).
3. Initial $[AHL_i] = 0,$ and $[AHL_e] = \mbox{constant.}$
4. LuxR protein ($[LuxR]$) is produced either by housekeeping gene which is assumed to have a constant transcription rate $k1,$ or by positive feedback system discussed in 7. LuxR degrades with constant $k12.$
5. AHL binds to LuxR protein, thus forming complex $[C]$ at rate $k5.$ The complex degrades into AHL and LuxR at rate $k13.$ 
6. The aforementioned complex dimerize into dimer $[D]$ at rate $k6,$ and the dimer dissociates and forms two complex at rate $k14.$
7. The complex binds to inducible promoter of LuxR and GFP ($[GFP]$), which has three characteristics:
1. The promoter complies to Hill's equation with cooperativity 1.
2. Maximum transcription rate of this promoter is $k2.$
3. Hill's dissociation constant is $k3.$
8. The degradation rate of GFP is $k18.$ 
9. Translation and degradation rate of mRNA is $k4$ and $k11,$ respectively.
10. Concentration of every species of molecules is adjusted for cell growth (dilution):
\[\dfrac{dx}{dt} = - \dfrac{k8*x5}{k9+x5}*x\]
11. Parameters concerning bacterial growth:
1. $k7$  is growth yield of bacteria
2. $k8$ is maximum growth rate 
3. $k9$  is half-saturation constant
4. $k10$ is bacterial death fraction per time
5. [S] is nutrient in medium
6. [N] is cell density

\[\dfrac{d[S]}{dt} = -\dfrac{1}{k7}\cdot [N]\cdot\dfrac{k8*[S]}{k9+[S]}\]

\[\dfrac{d[N]}{dt} = [N]\cdot\dfrac{k8*[S]}{k9+[S]} - k10*[N]\]

Implementation:

We implement the model in MATLAB. Users can simulate this deterministic model many times with different rate constant (by specifying the standard derivation of constants). 

Components of unknown variables:

x1 : [RHLR]

x2 : [Complex (AHL-RHLR)]

x3 : [dimer_of_complex]

x4 : [AHL_intracellular]

x5 : [S] (Nutrient in medium)

x6 : [N] (cell density)

x7 : [GFP]

System of ODEs:

\[\dfrac{d[LuxR]}{dt} =  (k4/k11)*(k1+(k2*x3/(k3+x3)))- k5*x1*x4 + k13*x2 - k12*x1 - (k8*x5/(k9+x5))*x1\]

\[\dfrac{d[Complex]}{dt} =   k5*x1*x4           - k13*x2           - 2*k6*(x2^2)          + 2*k14*x3          - (k8*x5/(k9+x5))*x2\]

\[\dfrac{d[dimer]}{dt} =   k6*(x2^2)- k14*x3- (k8*x5/(k9+x5))*x3\]

\[\dfrac{d[AHL_i]}{dt} =   k16*(AHLe - x4)  + k13*x2- k5*x1*x4 - k15*x4 - (k8*x5/(k9+x5))*x4\]

\[\dfrac{d[Nutrient]}{dt} =   (-1/k7)*x6*(k8*x5/(k9+x5))\]

\[\dfrac{d[Cell density]}{dt} =   x6*((k8*x5/(k9+x5)) - k10)\]

\[\dfrac{d[GFP]}{dt} =   (k17/k11)*((k2*x3/(k3+x3)))- k18*x7 - (k8*x5/(k9+x5))*x7\]

2D diffusion with cell division

We have extended simple cell model into 2D diffusion model as mentioned before, however both systems reach saturated concentration quickly, and it is far from reality (figure 1.)

If cell division occurred, more cells can receive the AHL and trigger the reaction, yields more GFP. The time reach GFP saturation will be much more later than the time which original static model did.
We build the 2D cell division model following these assumptions and constraints:
1. Collision between cells never occurred. Each cell will detect the neighbor cell and go away from the neighborhood. We can define the moving speed of each cell, the default speed is 1grid/15sec.
2. The cell grow horizontally, the cell division will occur after the cell reach two times of initial length. In fact we use 30 state (not including initial state) to describe the behavior of cell growing and division.

States of growing cell. Red represents intracellular region, green represents cell membrane region and blue represents extracellular region

    We can define the duration of each states of each cell. We set 60sec / states so that the cell will finish the division in 1800sec (30min) if no neighbor or obstacle exists.
3. Protein cannot diffuse to extra-cellular region during cell growing, moving and division.

Here is the simulation result ([AHL_Initial] = 10^-5 M)

Finally we can compare this model with real data and other two more simple model.  The simple model reach saturation within 500 seconds and the 2D diffusion model without cell division reach saturation within 1500 seconds.(Fig.1) However the real data and 2D model with cell division reach saturation much more later. Re scale the figure we can find the obvious relationship between simulation data and model data without numerical transformation and curve fitting (Fig.2):

Figure 1.  Time series from the simulation of three different models and the real data(solid line with marker)

Figure 2. Time series from simulation of 2D cell division model (dashed line) and real data (solid line with marker)

Introduction to stochastic modelling (by 彰引學長)

(By 彰引學長)
As we know, deterministic model for a biology system uses sets of ODEs (ordinary differential equation) to describe the system dynamic behaviour. These differential equation are treated under continuous manner, which is a approximation of large quantity of molecule number in the environment. Moreover, following the theory that each reaction is the collision with proper direction and energy of molecules, we can assume the molecules are spread uniformly about the environment and each reaction rate can be evaluated by classical chemical kinetic law. However, such assumption failed if molecule number is few.

When the molecule number is few, the effect of probability emerges. We should consider the probability of every effective collision. In this model, we define that a \emph{state} is one or more species with a particular number of molecule for each species. Each reaction involves at most two molecules as reactant. The reactions lead to \emph{state transitions}. Take our model into consideration. The species involving in this model are:

• DNA
• DNA with TF bound
• mRNA of LuxR
• mRNA of LuxR and GFP
• LuxR
• GFP
• AHL (intra-cell or extra-cell)
• AHL-LuxR complex

Assume that DNA is stable and AHL remains constant during the time we concern. All the reaction that will change the number of each species are:

Transcription of LuxR mRNA from DNA
\[DNA_{bound}+DNA_{unbound} \longrightarrow DNA_{bound}+DNA_{unbound}+mRNA_{LuxR}\]
Transcription of LuxR and GFP mRNA from bound DNA
\[DNA_{bound}\longrightarrow DNA_{bound}+mRNA_{LuxR-GFP}\]
Degradation of mRNA
\[mRNA_{LuxR}\longrightarrow \emptyset\]
\[mRNA_{LuxR-GFP}\longrightarrow \emptyset\]
Translation of LuxR
\[mRNA_{LuxR}\longrightarrow mRNA_{LuxR}+LuxR\]
Translation of GFP
\[mRNA_{LuxR-GFP}\longrightarrow LuxR+GFP\]
Degradation of LuxR
\[LuxR\longrightarrow \emptyset\]
Degradation of GFP
\[GFP\longrightarrow \emptyset\]
LuxR-AHL complex association
\[LuxR+AHL\longrightarrow complex\]
LuxR-AHL complex dissociation
\[complex\longrightarrow LuxR+AHL\]
Binding of complex and unbound DNA
\[complex+DNA_{unbound} \longrightarrow DNA_{bound} \]
Unbinding of complex and unbound DNA
\[DNA_{bound} \longrightarrow complex+DNA_{unbound}\]
Degradation of AHL
\[AHL\longrightarrow \emptyset\]

Some of you might wonder "where is the Hill equation?" In fact, Hill equation is embedded in this model. Assume the binding rate of $complex$ and $DNA_{unbound}$ is $k_1$ and unbinding rate is $k_2$. According to the reaction listed above, if $complex$ is constant, the equilibrium conditions are:
\[k_1\times  complex \times DNA_{unbound} =k_2 DNA_{bound}\]
\[\frac{DNA_{bound}}{DNA_{unbound}}=\frac{k_1\times  complex }{k_2}\]
Consider the total amount of DNA is a constant,
\[DNA_{bound}=DNA_{total}\frac{k_1\times complex}{k_2+k_1\times complex}\]
The generation of $mRNA$ is proportion to $DNA_{bound}$. Hence, the equation becomes a Hill equation with order $n=1$.

Introduction to Variance-based sensitivity analysis (VBSA)

The aim of sensitivity analysis is to examine which parameter in a model affects output most. In our model, there are 18 parameters (shown below) that have impact on GFP concentrations. The mathematical basis of variance-based sensitivity analysis is to estimate the effect of a parameter using variance. That is, how much portion of variance of GFP could the variance of a particular parameter explain?

Let $X = \{x_i\}|_{i=1}^{n}$ be $n$ parameters in out model, and that these parameters are in interval $[0,1]$ by transformation. It can be assumed that output $y = f(x_1, x_2, \ldots, x_n)$ can be decomposed as summation of several functions:

$$  f(X) = f_0 + \sum_{i=1}^{n} f_i(x_i) + \sum_{i<j}^n f_{i,j}(x_i,x_j) + \ldots  $$

By assumption of independence, one can derive:
$$ \int_0^1 f_{i_1,i_2,\ldots , i_s} (x_{i_1}, x_{i_2}, \ldots , x_{i_s})  dx_{i_k} = 0 \mbox{ for } k = 1 \mbox{ to } s $$
(where $\{i_a\}|_{a=1}^s$ is a sequence. For example, if $s=2, i_1=2, i_2=5$ then the former formula is equivalent to
$\int_0^1 f_{2,5} (x_{2}, x_{5})  dx_{2} = 0 \mbox{ and} \int_0^1 f_{2,5} (x_{2}, x_{5})  dx_{5} = 0 $）

This assumption implies orthogonality, which makes calculations of expected value much easier due to orthogonality:
$$ \begin{cases}
E(y) = f_0
\\
E(y|x_i) = f_0 + f_i
\\
E(y|x_i,x_j) = f_0 + f_i + f_j + f_{i,j} \mbox{ and so on}
\end{cases}$$

Therefore, the variance can be derived as follows:
$$ Var(y) = E(y^2) - (E(y))^2 = \int_{K^n} f^2 dX - f_0^2 = D  $$
Let D be abbreviation of integral defined below
$$ D_{i_1, i_2 \ldots , i_s} = \int_0^1\int_0^1\ldots\int_0^1 f_{i_1, i_2 \ldots , i_s}^2 dx_{i_1} \ldots dx_{i_s}. $$
Then the variance of output, denoted as $D,$ will be:
$$ D = \sum_{i=1}^n D_i + \sum_{i=1}^n\sum_{j=i+1}^n D_{i,j} + \ldots + D_{1,\ldots , n} $$
The index of how much portion of variance of GFP could the variance of a particular parameter explain is:
$$ S_i = \dfrac{D_i}{D} $$

Bibliography

1. Andrea Saltelli, Paola Annoni, Ivano Azzini, Francesca Campolongo, Marco Ratto, and Stefano Tarantola. Variance based sensitivity analysis of model output. Design and estimator for the total sensitivity index. Computer Physics Communications, 181(2):259{270, 2010
2. Karen Chan, Andrea Saltelli, Stefano Tarantola. SENSITIVITY ANALYSIS OF MODEL OUTPUT: VARIANCE-BASED METHODS MAKE THE DIFFERENCE. Proceedings of the 1997 Winter Simulation Conference (http://www.informs-sim.org/wsc97papers/0261.PDF)

iGEM衣服團購

iGEM衣服團購九月五日截止喔！
如要參賽或紀念請早洽facebook軮哲

2D deterministic model

If we construct the grids and simulate the simple model in each grid, and take the diffusion between grids into considerations, we can get 2D or 3D space models.

 We can set different parameters in each grid to simulate different environments in 2D space. For example, set all rate constant (except diffusion constant) to zero to simulate the extracellular environment, or set the low diffusion constant to  simulate the obstacle of diffusion in cell membrane. Hence we can define the extracellular/intracellular/cell membrane area by setting different rate or diffusion constant.
The following movie is the simulation results of 100*100 grids which contains 100 cells (the size of each cell is 6*6 grids), we can observe the GFP is generated only in cytoplasm:

Also we can observe the diffusion process of AHL:

Finally to get the different behavior of each grid, we can set the random rate constant on them, then the molecules will not be uniformly distributed in the cytoplasm:

Simple Deterministic Model

We modify the model which created in 2007 by the imperial iGEM team. The basic model simulate the behavior of single cell, and follow these assumptions:

1. All molecules, including proteins, protein complexes and small molecules are uniformly distributed in cell bodies.
2. The diffusion rate of AHL is determined only by the AHL gradient between cytoplasm and extracellular matrix..
3. The mechanism of GFP formation is shown below:

    a. At the beginning, the AHL receptor, LuxR protein, is produced by housekeeping gene and has a constant producing rate:

   b. AHL molecules bind LuxR and then the complex "LuxR-AHL" are formed:

  Where k2 is the rate constant of binding of LuxR and AHL, and k3 is the rate constant of dissociation of complexes "LuxR-AHL"

   c. The complexes "LuxR-AHL" bind the promoter and trigger the transcription and translation reaction to produce the GFP protein. The rate of GFP protein formation is determined by Hill function.

      To magnify the signal, we add the luxR gene at the downstream of  the gfp gene. The more AHL are detected, the more GFP and LuxR protein are generated and form a positive feedback loop.

      Where k4,Ke,n are the constant of Hill Function, and [LuxR_AHL] is the concentration of the complexes "LuxR-AHL"

4. All molecules, including LuxR protein, AHL, LuxR-AHL complex and GFP are degraded on corresponded rate constant.

Implementation:

We implement the model in MATLAB. Users can simulate this deterministic model many times with different rate constant (by specifying the standard derivation of constants). 

Figure. We simulate the model  100 times with different rate constant (by specifying the std derivation)