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)

VBSA中文簡介 part I

Sensitivity analysis的用意在為確認數學模型中，輸出（此例為為GFP的濃度）比較容易受到哪個參數（模型中的常數，如GFP production rate等）的影響。而"Variance based sensitivity analysis"指的是利用變異數來估算參數的影響力。例如：GFP的Variance有多少百分比可以被某參數解釋。

VBSA的假設其實非常像統計的ANOVA:
令$X = \{x_i\}|_{i=1}^{n}$為模型中的$n$個參數，另設參數皆位於$[0,1]$區間內（非位於此區間內可以用變數變換的方式達成）。那麼輸出 $y = f(x_1, x_2, \ldots, x_n)$可以表示成多個函數的和：
$$  f(X) = f_0 + \sum_{i=1}^{n} f_i(x_i) + \sum_{i<j}^n f_{i,j}(x_i,x_j) + \ldots  $$
並且假設以上函數對它自己的變數積分起來為零：
$$ \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 $$
(其中$\{i_a\}|_{a=1}^s$是個數列，例如$s=2, i_1=2, i_2=5$那麼上式可以表示為：
$\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 $）

如此一來這些函數就會正交(orthogonal，意味著兩個具有相同變數的函數乘積的積分為零)，所以要計算期望值等的就簡單了，因為除了要看的變數之外，積分統統會零掉：
$$ \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}$$

接下來就可以計算變異數：
$$ Var(y) = E(y^2) - (E(y))^2 = \int_{K^n} f^2 dX - f_0^2 = D  $$
為了方便，將以下積分表示簡寫：
$$ 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} $$
那麼輸出的變異數$D$就會是：(因為兩兩相乘的積分為零，只剩這些項)
$$ D = \sum_{i=1}^n D_i + \sum_{i=1}^n\sum_{j=i+1}^n D_{i,j} + \ldots + D_{1,\ldots , n} $$
所以Sensitivity的重點就在於「$D_i$能解釋多少的$D$」：
$$ S_i = \dfrac{D_i}{D} $$
這就是VBSA的中心概念，而estimators以及我上傳檔案的計算方式會放在part II

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)

部落格成立 (8/25)

成立Dry lab blog
這個blog的用意在於聯繫、紀錄Dry lab目前的進度、未來的規劃以及wiki的草稿，並避免被facebook留言洗到下面去而遍尋不著。

現在大家都可以進來觀看留言，但是如果要「發表文章」請在facebook通知我(Edward Cheng)或是在下面留言，也懇請隊友們一起分工合作！

現在mathjax已經加到這個部落格的<head>裡面了，所以直接打latex會直接排版！