Model

Introduction

In this model, we simplify the actual biology process into a model that only remains input molecule, promotor, transcription gene, mRNA, goal protein and output molecule.

Aim

1. This model is used to explain the phenomenon in our experiments.
2. This model is used to predict the working efficiency of our parts.
3. This model can be used to simulate our experiments and do parameter fitting.
4. This model gives some new perspectives and developed models on QS system.

Symbol

Symbol	Meaning
$v_{generate}$	The generation efficency of mRNA
$[X]$	The concentration of substance $X$
$g_{X}$	The generation rate of substance $X$
$\phi_{X}$	The decay rate of substance $X$
$V_{\max}$	The maximum rate of generation
$C_{saturated}$	The saturated concentration
$N_{\max}$	The maximum population
$r$	Growth rate of E.coli
$[S]_t$	Function of signal molecule decay
$R(t)$	Function of mRNA generate

Assumption

1. mRNA and proteins will decay following Poisson distribution (equivalent to birth-and-death process)
2. All combinations of two proteins are considered as quick reactions (Only control by thermodynamics)
3. The constitutive promoter has a constant rate to transcript proteins.
4. All raw materials inside cells can be considered as constants.

Basic Model

$$ \begin{aligned} \frac{d([mRNA])}{dt}&=v_{generate}-\phi_{mRNA}[mRNA]\\ \frac{d([protein])}{dt}&=g_{protein}[mRNA]-\phi_{protein}[protein] \end{aligned} $$

In these equations, $v_{generate}$ refers to the efficiency of mRNA transcription. $\phi$ refers to the degradation rate of mRNA and protein.

The property of $v_{generate}$ depends on the promoter and the concentration of inducer molecule. If the promoter is pcons, $v_{generate}$ is a constant. Otherwise, it will have a sensitive response to different concentration of inducer molecule. This reponse can be expressed as following form:

$$ v_{generate}([x])=V_{max}·(\frac{(1-\epsilon)·x^n}{k^n+x^n}+\epsilon) $$

$k$ refers to the dissociation constant and $x$ refers to the concentration of inducer concentration. $\epsilon$ refers to the leakage of genetic expression.

In comparision, for NOR GATE, the repression of inducer molecule can be expressed as similar form:

$$ v_{generate}([x])=V_{max}·(\frac{1-\epsilon}{1^n+(\frac{x}{k})^n}+\epsilon) $$

For specific concerntration, $v_{generate}$ is a constant, otherwise it is a function of $[x]$

The generated protein is used to produce new signal molecule, which play a role as enzyme. Different from Michaelis-Menten equation, our protein (in other words, enzyme) will degradate while producing new siginal molecule, So this fact should be considered into our fundmental model.

Mathematical expression for producing new signal molecule:

$$ \begin{aligned} \frac{d[EAB]}{dt}&=k_1[E][A][B]-(k_1+k_{-1})[EAB]\\ \frac{d[M_{signal}]}{dt}&=k_2[EAB] \end{aligned} $$

Developed Model

Growth of E.coli

In the developed model, we first take the growth of E.coli into consideration. The growth of E.coli can not only fluctuate the concentration of both reactants and products, but also an important variable in calculate final concentration of products. This model is based on this two fundamental relation:

$$ \begin{aligned} Total&=Concentration·Volumel\\ Volume&=N_{E.coli}·V_{E.coli}\\ \frac{d([protein]·Volume)}{dt}&=g_{protein}[mRNA]·Volume-\phi_{protein}[protein]·Volume \end{aligned} $$

Correspondingly, it is same to equation for mRNA expression:

$$ \frac{d([mRNA])·Volume}{dt}=v_{generate}·Volume-\phi_{mRNA}[mRNA]·Volume\\ $$

$N_{E.coli}$ is a function used to show the population of E.coli, $V_{E.coli}$ refers to the volume of every E.coli, as a constant. So we can divide out the constant $V_{E.coli}$ on both sides of every equations, and take derivative formula:

$$ \frac{d[protein]}{dt}·N_{E.coli}+\frac{dN_{E.coli}}{dt}·[protein]=g_{protein}[mRNA]·N_{E.coli}-\phi_{protein}[protein]·N_{E.coli} $$

Simplify this equation into following form:

$$ \frac{d[protein]}{dt}=g_{protein}[mRNA]-(\phi_{protein}+\frac{N_{E.coli}'}{N_{E.coli}})[protein]\\ N_{E.coli}'=\frac{dN_{E.coli}}{dt} $$

$N_{E.coli}$ is satisfied to following equation:

$$ \begin{aligned} \frac{dN_{E.coli}}{dt}&=rN_{E.coli}(1-\frac{N_{E.coli}}{N_{\max}})\\ N_{E.coli}&=\frac{N_{\max}}{1+(\frac{N_{\max}}{N_{t=0}}-1)·e^{-rt}} \end{aligned} $$

$r$ refers to growth rate of E.coli and $N_{\max}$ refers to the limits of E.coli population. Since $N_{\max}$ and $N_{t=0}$are constants, so we define following parameter:

$$ \frac{N_{\max}}{N_{t=0}}-1=N_{c} $$

And $\frac{N'}{N}$ equals:

$$ \frac{N'}{N}=\frac{N_cre^{-rt}}{1+N_ce^{-rt}} $$

From our experiments, we find there are another two possible factors affecting the production of our system. First one is diffusion of signal molecule at initial time, the other one is the decay of signal molecule with the time flying.

Diffusion of signal molecule at initial time

The concentration of signal is always considered to diffuse into E.coli very rapidly. But from our data, we find that the initial part of our dynamic curve is not fitting to our basic model. Our basic model indicates that the rate of generating will decrease with the time flying, but the experiment shows that the velocity will have a short rise at initial time and then decrease as the way predicted by basic model. Therefore, we take process of diffusion into consideration. Because at very beginning, the concentration of signal in E.coli is very low, and then it will rise by diffusion, so the efficiency of production will rise according to time in a short time period.

We suppose the initial concentration difference between inside of E.coli and outside is $\Delta c(0)$, also we know the time for E.coli to balence this difference:

$$ c(t)= C_{saturated} -\Delta c(0)·e^{-\frac{t}{\tau}} $$

So the generating efficency comes to:

$$ v_{generate} = \frac{V_{\max}}{1+(\frac{k}{ C_{saturated} -\Delta c(0)·e^{-\frac{t}{\tau}}})^n} $$

And we will use this formula to simulate initial state.

The demo is shown above which is a Log linear plot. X-axis refers to the time, Y-axis refers to the generating efficiency. We can easily figure out the concentration will rapidly get to steady state and remains to a constant. Therefore, it will only affect the inital transcription efficiency.

Decay of signal molecule

In basic model, we consider the decay of signal can be neglected because we found there's no significant difference between concentration in vitro. But actually when we meature the rough concentration in the LB with E.coli, we found that the concentration has a linear deacrease through time, which we should take consideration into our model.

The decay can be shown as following equation:

$$ [S]_t=[S]_{initial}-k_{decay}t $$

And the $v_{generate}$ becomes to:

$$ v_{generate}= V_{\max}\frac{([S]_t)^n}{k^n+([S]_t)^n} $$

To illustrate the change taken by the decompose of signal molecule, we can see following simulation curves:

X-axis refers to time. We find the efficiency will not be disturbed greatly at initial time, and will have a rapid decrease when the concentration equals to the half of origin. This property shows that we should control the reaction time otherwise the production will decay without production with the time going by. So the main purpose of this model is to predict when we dilute the input signal solution to obtain the maximum of protein to convert out signal.

Extra model

This part will discuss an interesting model on how the signal molecule affect the growth and population. The reason why we care about this question is that we measured the OD600 under different circumstance and found some special relation between the concentration and the population. In breif, with the rise of concentration, the population will decrease. We wonder the mechanism and propse two hypothesis:

1. The signal molecule is toxic to E.coli, so the population will decrease related to the increase of concentration linearly.

2. The signal molecule induce the synthesis of GFP which occupy the substance that is originally used for growth. It indicates that if the GFP is produced, then the population will be at low level, otherwise the population will be at normal level.

   In our model, we indicates the second hypothesis is more realistic.

   To show the difference between these two hypothesis, we give following equation and diagram:

   For the first hypothesis, we relation between

Model of parameter fitting and simulation

Hill equation

To get the parameter of Hill equation through our data, we tranfer Hill equation to following form:

$$ \begin{aligned} Hill\ equation&:y=V_{max}\times\frac{x^n}{k^n+x^n}\\ New\ form&:\log{\frac{\frac{y}{V_{max}}}{1-\frac{y}{V_{max}}}}=n\log{x}-n\log{k} \end{aligned} $$

In this form, we can get easily get a linear relation between our input concerntration and output GFP. The question is how to find out $V_{max}$ in this equation because this value determine the reprocessed data of output. Another question is, due to the large scale of our data, to ease the workload of proceesing such data. To meet the needs of these two question, first we let each output data substract the minimum among all output data, and define the ratio between each processed output data and the maximum of all output data as the standard output. (NOTICE: The minimum data of this output data set can be the control.)As following shows:

$$ {output}={y_1,y_2,···,y_n} $$ $$ SY_{output}=\{y_1',y_2',···,y_n'\}\quad which\quad y_i=\frac{y_i-\min{Y_{output}}}{\max{Y_{output}}-\min{Y_{output}}} $$

The elements in $SY_{output}$ fit following equation:

$$ \log{\frac{{y_i'}\frac{\max{Youtput}-\min{Y_{output}}}{V_{max}}}{1-{y_i'}\frac{\max{Youtput}-\min{Y_{output}}}{V_{max}}}}=n\log{x_i}-n\log{k} $$

We define the value of $\frac{V_{max}}{\max{Youtput}-\min{Y_{output}}}$ as a parameter $PV_{max}$. So the equation we actually simulate is following one:

$$ \log{\frac{y_i'}{PV_{max}-y_i'}}=n\log{x_i}-n\log{k} $$

​ We use Mathematica as fitting tools, the following code is shown:

outputdata = {output1, output2, output3, output4, output5,output6};
Processeddata = (outputdata - Min[outputdata])/(Max[outputdata] - 
      Min[outputdata]) // N;
data' = {{Log10[10^(-9)], Processeddata[[1]]}, {Log10[10^(-8)], 
    Processeddata[[2]]}, {Log10[10^(-7)], 
    Processeddata[[3]]}, {Log10[10^(-6)], 
    Processeddata[[4]]}, {Log10[10^(-5)], Processeddata[[5]]},{Log10[10^(-4)], Processeddata[[6]]}};
data = {{data'[[1, 1]], data'[[1, 2]]}, {data'[[2, 1]], 
    data'[[2, 2]]}, {data'[[3, 1]], data'[[3, 2]]}, {data'[[4, 1]], 
    data'[[4, 2]]}, {data'[[5, 1]], data'[[5, 2]]}, {data'[[6, 1]], data'[[6, 2]]}};
solu = Flatten[
   Solve[Log10[(y*PVmax)/(1 - (y*PVmax))] == n*x - n*logk, y]];
fitparameter = (FindFit[data, y /. solu, {PVmax, logk, n}, x])
fit = y /. solu /. fitparameter;
Show[ListPlot[data, PlotStyle -> Red], Plot[fit, {x, -10, 0}]]

Example and its output is shown (NOTICE: This example is the fitting curve of the Tra with its limited five data. Actually most of our data, except for tra, has six inputs and outputs, so the original code, which is shown above, has six outputs. When we use this code, we can just import outputs into "outputdata" list and run this programm. ):

outputdata = {16141, 6812, 32977, 362525, 959405};
Processeddata = (outputdata - Min[outputdata])/(Max[outputdata] - 
      Min[outputdata]) // N;
data' = {{Log10[10^(-9)], Processeddata[[1]]}, {Log10[10^(-8)], 
    Processeddata[[2]]}, {Log10[10^(-7)], 
    Processeddata[[3]]}, {Log10[10^(-6)], 
    Processeddata[[4]]}, {Log10[10^(-5)], Processeddata[[5]]}};
data = {{data'[[1, 1]], data'[[1, 2]]}, {data'[[2, 1]], 
    data'[[2, 2]]}, {data'[[3, 1]], data'[[3, 2]]}, {data'[[4, 1]], 
    data'[[4, 2]]}, {data'[[5, 1]], data'[[5, 2]]}};
solu = Flatten[
   Solve[Log10[(y*PVmax)/(1 - (y*PVmax))] == n*x - n*logk, y]];
fitparameter = (FindFit[data, y /. solu, {PVmax, logk, n}, x])
fit = y /. solu /. fitparameter;
Show[ListPlot[data, PlotStyle -> Red], Plot[fit, {x, -10, 0}]]

Then we can get the meaningful parameter from these data quickly and easily.

Simulation of Signal Producing

The Efficiency of Signal Converter

How we can measure the working efficiency of our signal converter is an important question for us. As we all know, the reason why we use GFP to reflect the efficiency of promoter is that we can measure fluoresence easily and establish the quantity relationship between GFP expression and input signal concentration. But when it comes to some other products such as small molecule, they are hard to measure exactly. We use LC-MS to indicate the production of our signal converter roughly, but this data is too rough to instruct our following work. So we will use our model to obtain the parameter of converter indirectly by following experiments and deduction from model.

We symbol $S_1,S_2$ as the concentrations of two signal molecules, signal one and signal two, $GFP$ as the result of fluroesence intensity.

We propose two experiments. First one is using signal two to induce the expression of GFP. We take its results as standard curve. The other experiment is using signal one to obtain signal two, and we use signal two to induce the expression of gene. Also we will have following data:

$$ \begin{aligned}S_1&=\{c_1,c_2,···,c_n\}\\ GFP&=\{F_1,F_2,···,F_n\} \end{aligned} $$

From our model we know the relationship among $S_1,S_2$ and $GFP$ at steady state as following:

$$ \begin{aligned}GFP&=V_{max}·(\frac{(1-\epsilon_1)·{S_2}^n}{k_1^n+{S_2}^n}+\epsilon_1)\\ S_2&=V_{max}·(\frac{(1-\epsilon_2)·{S_1}^m}{k_2^m+{S_1}^m}+\epsilon_2)\end{aligned} $$

From the parameter fitting model, we can determine all parameters in $GFP-S_2$ curve. Therefore, we can use this curve and data of GFP from second experiment to obtain the input signal two concentration.

$$ \begin{aligned} F_i&=V_{max}·(\frac{(1-\epsilon_1)·{[S_2]_i}^n}{k_1^n+{[S_2]_i}^n}+\epsilon_1)\\ F'_i&=\frac{F_i-\epsilon_1}{1-\epsilon_1} \\\Longleftrightarrow\log{[S_2]_i}&=\frac{\log{\frac{F_i'}{V_1-F_i'}}}{n}+\log{k_1} \end{aligned} $$

So we have the data $[S_2]_i$ related to input concentration of signal one, so we can get the relation through using parameter-fitting model would get the parameter of $S_1-S_2$ curve finally.

Rough schematic diagram

This is the concentration curve of protein related to time

This is the concentration curve of protein complex related to time

This is the concentration curve of producing signal molecule related to time

Simulation of NOR GATE

Rough schematic diagram

This is the concentration curve of produced signal molecule related to time

Theoretical Calculation of Modeling

Solution to ODE

The core of our model is to solve following equation and find parameters from experiments:

$$ \frac{dy}{dt}+P(t)y=Q(t) $$

The solution can be decomposed to two parts:

$$ \begin{aligned} \frac{dy}{dt}+P(t)y&=0\\ \frac{dy_{s}}{dt}+P(t)y_{s}&=Q(t) \end{aligned} $$

From fisrt equation we will get:

$$ y=Ce^{-\int P(t)\,dt} $$

How can we use the solution to first equation to solve second equation? The answer is to transfer constant $C$ into a function related to $t$. And the derivative will become to following formula:

$$ \begin{aligned} \frac{dy}{dt}&=C(t)·(-P(t))·e^{-\int P(t)\,dt}+C'(t)·e^{-\int P(t)\,dt}\\ \Longrightarrow \frac{dy}{dt}+P(t)y&=C'(t)·e^{-\int P(t)\,dt}\\ \Longrightarrow Q(t)&=C'(t)·e^{-\int P(t)\,dt}\\ \Longrightarrow C(t)&=\int Q(t)·e^{\int P(t)\,dt}\,dt+C \end{aligned} $$

Therefore, the solution to second equation is:

$$ e^{-\int P(t)\,dt}(\int Q(t)·e^{\int P(t)\,dt}\,dt+C) $$

The difficulty is how we can use such a complex function in next differential equation? Actually we probably cannot get the analytic result of the integral, so it seems impossible to get an exact function for protein concentration. Fortunately, there are still some special properties in our function which wll help us to get a relative solution.

We start from the function of mRNA. Since $P(t)$ is a constant in our first equation, we can directly give the result:

$$ [mRNA]= e^{-\phi_{mRNA}t}(\int Q(t)·e^{\phi_{mRNA}t}dt+C) $$

Now we solve following differetial equation:

$$ \frac{d([protein])}{dt}+\phi_{protein}[protein]=g_{protein}[mRNA] $$

Or for simplicity, we use:

$$ \frac{dy}{dt}+\phi_2·y=g·R(t) $$

According to the differential operator method, we get:

$$ \begin{aligned} (D+\phi_2)y^\ &=g·R(t)\\ \Longleftrightarrow y^*&=\frac{1}{D+\phi_2}·g·R(t)\\ \Longleftrightarrow y^*&=\frac{1}{\phi_2}·(1-(\frac{D}{\phi_2})+(\frac{D}{\phi_2})^2-···)·g·R(t)\\ \Longleftrightarrow y^*&=\frac{1}{\phi_2}·(1-\frac{1}{\phi_2}\frac{d}{dt}+\frac{1}{\phi_2^2}\frac{d^2}{dt^2}-···)·g·R(t) \end{aligned} $$

For $R(t)$, we write the general form:

$$ R(t)=e^{-\phi t}(\int Q(t)e^{\phi t}dt+C) $$

When we take derivation:

$$ R'(t)=(-\phi)·e^{-\phi t}(\int Q(t)e^{\phi t}dt)+e^{-\phi t} Q(t)e^{\phi t}+C·(-\phi)·e^{-\phi t}\\ \Longleftrightarrow R'(t)=(-\phi)R(t)+Q(t) $$

Furthermore:

$$ \begin{aligned} R^{(n)}(t)&=(-\phi)R^{(n-1)}(t)+Q^{(n-1)}(t)\\ R^{(n)}(t)&=(-\phi)^n·R(t)+\sum_{k=1}^{n}k^{n-k}Q^{(k)}(t) \end{aligned} $$

REMARK:

$$ f^{(n)}(t)=\frac{d^nf}{dt^n} $$

Therefore we get:

$$ y^*=\frac{g}{\phi_2}·(\sum_{i=0}^{+\infty}(\frac{\phi_1}{\phi_2})^iR(t)+\sum_{k=1}^{i}(\frac{\phi_1}{\phi_2})^{i}·(\phi_1)^{-k}·Q^{(k)}(t))\\ \Longrightarrow y=\frac{g}{\phi_2}·(\sum_{i=0}^{+\infty}(\frac{\phi_1}{\phi_2})^iR(t)+\sum_{k=1}^{i}(\frac{\phi_1}{\phi_2})^{i}·(\phi_1)^{-k}·Q^{(k)}(t))+A^*·e^{-\phi_2 t} $$

The first summation is simple:

$$ \sum_{i=0}^{+\infty}(\frac{\phi_1}{\phi_2})^iR(t)=\frac{\phi_2}{\phi_2-\phi_1}R(t) $$

Second summation is really complex, so we must do some approximation:

$$ \begin{aligned} &\sum_{i=0}^{+\infty}(\frac{\phi_1}{\phi_2})^{i}\sum_{k=1}^{i}(\phi_1)^{-k}·Q^{(k)}(t)\\ =&\sum_{i=0}^{+\infty}(\frac{\phi_1}{\phi_2})^{i}(\phi_1)^{-1}·\frac{d}{dt}Q(t)\\ =&\frac{\phi_2}{\phi_1(\phi_2-\phi_1)}·\frac{d}{dt}Q(t) \end{aligned} $$

Therefore we get a approximation of protein's concentration:

$$ [protein]=\frac{\phi_2}{\phi_2-\phi_1}e^{-\phi_{mRNA}t}(\int Q(t)·e^{\phi_{mRNA}t}dt+C)+\frac{\phi_2}{\phi_1(\phi_2-\phi_1)}·\frac{d}{dt}Q(t)+A^*·e^{-\phi_2 t} $$

Solution to Our Model

Details of Developed Model

Growth of E.coli

We combine this solution with our equation, and then we get:

$$ \begin{aligned} \left[mRNA\right]&=e^{-\int(\phi_{mRNA}+\frac{r·N_c}{N_c+e^{rt}})\,dt}·(\int v_{generate}·e^{\int(\phi_{mRNA}+\frac{r·N_c}{N_c+e^{rt}})\,dt}\,dt+C_0)\\ \left[protein\right]&=e^{-\int(\phi_{protein}+\frac{r·N_c}{N_c+e^{rt}})\,dt}·(\int g_{protein}\left[mRNA\right]·e^{\int(\phi_{protein}+\frac{r·N_c}{N_c+e^{rt}})\,dt}\,dt+C_0') \end{aligned} $$

We suppose that:

$$ N_c(t)=1+N_c·e^{-rt} $$

Therefore we get:

$$ [mRNA]=C_1·v_{generate}·N_c(t)·e^{-\phi_{mRNA}t}\int\frac{e^{\phi_{mRNA}t}}{N_c(t)}\,dt+C_1·C_0N_c(t)·e^{-\phi_{mRNA}t} $$

As a special case, this is used to decribe if the growth of E.coli is at a steady state:

$$ \lim_{n\to\infty}N_c(t)=1 $$

Then we get a simple formula：

$$ [mRNA]=A·v_{generate}+C·e^{-\phi_{mRNA} t} $$

Further more, we define:

$$ \begin{aligned} A(t)&=C_1·N_c(t)·e^{-\phi_{mRNA}t}\int\frac{e^{\phi_{mRNA}t}}{N_c(t)}\,dt\\ C(t)&=C_1·C_0N_c(t)\\ [mRNA]&=A(t)·v_{generate}+C(t)·e^{-\phi_{mRNA} t} \end{aligned} $$

Consider the inital value of mRNA, we get following relation:

$$ A(0)v_{generate}+C(0)=0 $$

Now let's have a look on this special function and related integration:

$$ A(t)=C_1·N_c(t)·e^{-\phi_{mRNA}t}\int\frac{e^{\phi_{mRNA}t}}{N_c(t)}\,dt\\ N_c(t)=1+N_c·e^{-rt}\\ N_c=\frac{N_{\max}}{N_{t=0}}-1 $$

We can hardly get an analytic solution to this integration theoritically, but we can do some transformation on $N_c(t)$, which helps us solve this problem partly according to this fact:

$$ If\quad|x|<1\\Then\quad\frac{1}{1+x}=\sum_{k=0}^{+\infty}(-x)^k $$

So we suppose:

$$ N_c<1\Longleftrightarrow N_{t=0}>\frac{N_{max}}{2} $$

From the biological perspective, this indicates the initial population of E.coli has been more than the half of maximum population, this assumption roughly fits our experiments. This condition promises following equation:

$$ \because t>0,e^{-rt}<1 \\\therefore N_c·e^{-rt}<1\\ \therefore \frac{1}{N_c(t)}=\sum_{k=0}^{+\infty}(-N_ce^{-rt})^k $$

So we will have:

$$ \int\frac{e^{\phi_{mRNA}t}}{N_c(t)}\,dt=\int\sum_{k=0}^{+\infty}(-N_c)^ke^{(\phi_{mRNA}-kr)t}\,dt\\=\sum_{k=0}^{+\infty}\int(-N_c)^ke^{(\phi_{mRNA}-kr)t}\,dt\\ =\sum_{k=0}^{+\infty}(-N_c)^k(\phi_{mRNA}-kr)^{-1}e^{(\phi_{mRNA}-kr)t} $$

And:

$$ \begin{aligned} A(t)&=C_1·N_c(t)·e^{-\phi_{mRNA}t}\int\frac{e^{\phi_{mRNA}t}}{N_c(t)}\,dt\\ &=C_1N_c(t)·\sum_{k=0}^{+\infty}(-N_c)^k(\phi_{mRNA}-kr)^{-1}e^{-krt}\\ &=C_1N_c(t)·\sum_{k=0}^{+\infty}\frac{(-N_c·e^{-rt})^k}{\phi_{mRNA}-kr} \end{aligned} $$

Therefore:

$$ [mRNA]=C_1N_c(t)·v_{generate}·\sum_{k=0}^{+\infty}\frac{(-N_c·e^{-rt})^k}{\phi_{mRNA}-kr}+C_1C_0N_c(t)·e^{-\phi_{mRNA}t} $$

Before we use this formula to obtain the expression of protein's concentration, we should analyze and simplify it.

Property i :

$$ \exists k_0,\forall k>k_0,|\phi_{mRNA}-kr|>1\\ \Longleftrightarrow k_0>\frac{1+\phi_{mRNA}}{r} \\\therefore \sum_{k=0}^{+\infty}\frac{(-N_c·e^{-rt})^k}{\phi_{mRNA}-kr}=\sum_{k=0}^{k_0}\frac{(-N_c·e^{-rt})^k}{\phi_{mRNA}-kr}+\sum_{k=k_0+1}^{+\infty}\frac{(-N_c·e^{-rt})^k}{\phi_{mRNA}-kr} $$

For the first part:

$$ \begin{aligned} &S_1=\sum_{k=0}^{k_0}\frac{(-N_c·e^{-rt})^k}{\phi_{mRNA}-kr}\\ &=\frac{1}{\phi_{mRNA}}-\frac{N_ce^{-rt}}{\phi_{mRNA}-kr}+\frac{N_c^2e^{-2rt}}{(\phi_{mRNA}-kr)^2}-·····+(\frac{-N_ce^{-rt}}{\phi_{mRNA}-kr})^{k_0}\\ \end{aligned}\\ \lim_{t\to\infty}S_1=\frac{1}{\phi_{mRNA}} $$

For the second part:

$$ |S_2|=\sum_{k=k_0+1}^{+\infty}|\frac{(-N_c·e^{-rt})^k}{\phi_{mRNA}-kr}|<\sum_{k=k_0}^{+\infty}(N_c·e^{-rt})^k=\frac{(N_ce^{-rt})^{k_0+1}}{1-N_ce^{-rt}} \\\ 0\le\lim_{t\to\infty}S_2\le\lim_{t\to\infty}|S_2|\le\lim_{t\to\infty}\frac{(N_ce^{-rt})^{k_0+1}}{1-N_ce^{-rt}}=0\\ \therefore \lim_{t\to\infty}S_2=0 $$

Therefore:

$$ \lim_{t\to\infty}A(t)=\lim_{t\to\infty}C_1·v_{generate}·(1+N_ce^{-rt})(S_1+S_2)\\ =\lim_{t\to\infty}C_1·v_{genrate}·(S_1+S_2+N_ce^{-rt}S_1+N_ce^{-rt}S_2)\\ =\frac{C_1}{\phi_{mRNA}}·v_{generate}+0+0+0\\ =\frac{C_1}{\phi_{mRNA}}·v_{generate} $$

Property ii :

$$ A(t)\approx \frac{C_1·v_{generate}}{\phi_{mRNA}}+{C_1·v_{generate}}(\frac{N_c}{\phi_{mRNA}}-\frac{N_c}{\phi_{mRNA}-r})·e^{-rt}\\+{C_1·v_{generate}}(\frac{N_c^2}{(\phi_{mRNA}-2r)^2}-\frac{N_c^2}{\phi_{mRNA}-r})·e^{-2rt}+o((rt)^3) $$

So we finally get:

$$ [mRNA]= \frac{C_1·v_{generate}}{\phi_{mRNA}}+G·e^{-rt}+H·e^{-2rt}+C(t)·e^{-\phi_{mRNA} t} $$

Now we use this formula to solve following ODE:

$$ [protein]=e^{-\int(\phi_{protein}+\frac{r·N_c}{N_c+e^{rt}})\,dt}·(\int g_{protein}[mRNA]·e^{\int(\phi_{protein}+\frac{r·N_c}{N_c+e^{rt}})\,dt}\,dt+C_0') $$

From previous calculate we could guess the approximate solution to protein's concentration will be following form:

$$ [protein]=A'(t)+B'(t)e^{-rt}+C'(t)e^{-\phi t}+D'(t)e^{-2rt}+E'(t)e^{-(r+\phi)t}+F'(t)e^{-\phi't}\\ \phi = \phi_{mRNA}\\ \phi'= \phi_{protein} $$

Or we can appromixately consider this formula as:

$$ [protein]=S(t)+T(t)·e^{-\kappa t}\\ $$

$\kappa$ is a parameter used to reflect the fact comprehensively.

Finally we get:

$$ \lim_{t\to\infty}[protein]=S=\frac{C_1C_2 g_{protein}}{\phi_{protein}\phi_{mRNA}}v_{generate}\\ \Longrightarrow S\varpropto v_{generate} $$

This result indicates the generated protein concentration has a direct relation with input signal molecule concentration. More importantly, we use Hill equation to describe the final product concentration induced by different concentration of input signal molecule is approperiate.

In our case, after renewing with fresh LB solution, the protein will degradate and never generate new. So another dofferential equation is needed to describe this situation:

$$ \frac{d[protein]}{dt}=-\phi_{protein}[protein] $$

The initial value of this equation is:

$$ [protein]|_{T=t_0}=S $$

Then the function will be:

$$ [protein]=S·e^{-\phi_{protein} t} $$

Diffusion of signal molecule at initial time

Review

We suppose the initial concentration difference between inside of E.coli and outside is $\Delta c(0)$, also we know the time for E.coli to balence this difference:

$$ c(t)= C_{saturated} -\Delta c(0)·e^{-\frac{t}{\tau}} $$

So the generating efficency comes to:

$$ v_{generate} = \frac{V_{\max}}{1+(\frac{k}{ C_{saturated} -\Delta c(0)·e^{-\frac{t}{\tau}}})^n} $$

And we will use this formula to give the initial state.

How to solve?

First we have following relations in mathematics:

$$ \lim_{x\to0}\frac{(1+x)^n}{1+nx}=1\\ \lim_{x \to 0}\frac{1+x^n}{1-x^{2n}} = \lim_{x \to 0}\frac{1}{1-x^{n}}=1 $$

These two equation indicate a group of equivalent infinitesimal, which we can use to do approximation in our problem. The approximation can be done as following way by using two properties:

$$ v_{generate} = \frac{V_{\max}}{1+(\frac{k}{ c(t)})^n}\\ =V_{\max}-\frac{V_{\max}}{1+(\frac{c(t)}{k})^n}\\ =V_{\max}-V_{\max}[1-\frac{c(t)}{k})^n]\\ =V_{\max}(\frac{c(t)}{k})^n\\ =V_{\max}(\frac{C_{saturated}}{k})^n(1-\frac{\Delta c}{k}e^{-\frac{t}{\tau}})^n\\ =V_{\max}(\frac{C_{saturated}}{k})^n(1-n\frac{\Delta c}{k}e^{-\frac{t}{\tau}}) $$

For simplicity, we can rewrite into a simple equation:

$$ v'_{generate} = V'_{\max}-\delta e^{-\frac{t}{\tau}} $$

And the ODE for mRNA can be written into:

$$ \frac{d([mRNA])}{dt}=V'_{\max}-\delta e^{-\frac{t}{\tau}}-\phi_{mRNA}[mRNA] $$

Solution:

$$ [mRNA]=\frac{V_{\max}}{\phi_{mRNA}}-\frac{\tau \delta}{\tau\phi_{mRNA}-1}e^{-\frac{t}{\tau}}+(-\frac{V_{\max}}{\phi_{mRNA}}+\frac{\tau \delta}{\tau\phi_{mRNA}-1})e^{-\phi t} $$

Correspondingly, the function of protein is:

$$ [protein]=\frac{V_{\max}}{\phi_{mRNA}\phi_{protein}}-\frac{\tau^2 \delta}{(\tau\phi_{protein}-1)(\tau\phi_{mRNA}-1)}e^{-\frac{t}{\tau}}\\+\frac{1}{\phi_{protein}-\phi_{mRNA}}(-\frac{V_{\max}}{\phi_{mRNA}}+\frac{\tau \delta}{\tau\phi_{mRNA}-1})e^{-\phi t}+C'e^{-\phi_{protein}t}\\ C'=-\frac{V_{\max}}{\phi_{mRNA}\phi_{protein}}+\frac{\tau^2 \delta}{(\tau\phi_{protein}-1)(\tau\phi_{mRNA}-1)}\\-\frac{1}{\phi_{protein}-\phi_{mRNA}}(-\frac{V_{\max}}{\phi_{mRNA}}+\frac{\tau \delta}{\tau\phi_{mRNA}-1}) $$

The simulation curve for an arbitraty number:

We can see the initial slope of the curve is rasing to a point and then decrease gradually which is highly fixed to the experiment result we get.

Decay of signal molecule

Review

In basic model, we consider the decay of signal can be neglected because we found there's no significant difference between concentration in vitro. But actually when we meature the rough concentration in the LB with E.coli, we found that the concentration has a linear deacrease through time, which we should take consideration into our model.

The decay can be shown as following equation:

$$ [S]_t=[S]_{initial}-k_{decay}t $$

And the $v_{generate}$ becomes to:

$$ v_{generate}= V_{\max}\frac{([S]_t)^n}{k^n+([S]_t)^n}\\ \frac{d}{dt}v_{generate}=-\frac{ V_{\max}·k_{decay}}{k}\frac{n([S]_t)^{n-1}}{(k^n+([S]_t)^n)^2} $$

According to the solution we deduced before, we have:

$$ [protein]=\frac{\phi_2}{\phi_2-\phi_1}e^{-\phi_{mRNA}t}(\int v_{generate}·e^{\phi_{mRNA}t}dt+C)+\frac{\phi_2}{\phi_1(\phi_2-\phi_1)}·\frac{d}{dt}v_{generate}+A^*·e^{-\phi_2 t} $$

Surely the first step is to confirm this equation gives a reasonable result. Through using following mathematics conclution, we can approximately consider the integral as a summation:

$$ \int f(t) dt=F(t)+C=\int_a^tf(\mu)d\mu+F(a)+C $$

We assume $a=0$ which has no effect to the formula but has its biological meaning, which is the starting timepoint. So we have:

$$ \int v_{generate}·e^{\phi_{mRNA}t}dt=\int_0^t v_{generate}(\mu)·e^{\phi_{mRNA}\mu}d\mu\\ \approx\sum_{i=0}^{n}V_{\max}\frac{([S]_{initial}-k_{decay}i\Delta t)^n}{k^n+([S]_{initial}-k_{decay}i\Delta t)^n}·e^{\phi_{mRNA}i\Delta t}·\Delta t $$

Which

$$ n\Delta t=t $$

We use matlab to obtain a rough curve. X-axis refers to time.

This is a important result because it indicates that the production will not always increase with the time going. Actually, there exists a so-called "best time" to process next step in our system. For example, this peak can determine when we dilute input signal to get output signal as much as possible.

Red stars refers to "best time" according to different input concentration from upstream block.

*matlab code:

n = [];
fn = [];
for i=1:T/dt
n = [n i];
t = exp(-a*i*dt);
sum=0;
for j=0:i
sum = sum + (Vm-(j*dt)^n)*exp(a*dt*j)*dt/(k^n+(Vm-(dt*j)^n));
end
y = t*sum+\phi* (Vm - dt*i)^(n-1)/(k^n + (Vm - dt*i)^n)^2;
fn = [fn y];
end
plot(n,fn);
max(fn);

This matlab code shows how we draw the curves and how to find maximum.

Signal Producing

Review

In last part, we gives a approximate value of the protein we will get from our system:

$$ \lim_{t\to\infty}[protein]=S=\frac{C_1C_2 g_{protein}}{\phi_{protein}\phi_{mRNA}}v_{generate} $$

If we consider the decay of molecule, then we rewrite equation above as:

$$ [protein]_{\max}=S|_{t=t_{\max}}(1+\eta(\frac{v'_{generate}}{v_{generate}}+\frac{{v''_{generate}}}{v_{generate}})|_{t=t_{\max}})\\=\frac{C_1C_2 g_{protein}}{\phi_{protein}\phi_{mRNA}}v_{generate}|_{t=t_{\max}}(1+\eta(\frac{dv_{generate}}{dt}+\frac{{d^2v_{generate}}}{dt^2})|_{t=t_{\max}}) $$

Which:

$$ (\frac{dv_{generate}}{dt}+\frac{{d^2v_{generate}}}{dt^2})|_{t=t_{\max}}<0 $$

For simpilicity, we define the final production of goal protein as

$$ [protein]=S\propto v_{generate}|_{t=t_{max}} $$

Analyse

Now we focus on the differential equation related to the signal production:

$$ \frac{d[EAB]}{dt}=k_1[E][A][B]-(k_1+k_{-1})[EAB]\\ \frac{d[M_{signal}]}{dt}=k_2[EAB] $$

In this equation set, $[E]$ equals to the concentration of protein.

$$ [E]=[protein] $$

Finally we have:

$$ [EAB]=\frac{k_1[A][B]S}{-\phi_{protein}+k_{-1}+k_{1}}·e^{-\phi_{protein}t}+C_2·e^{-(k_{-1}+k_{1})t} $$

And the initial state:

$$ [EAB]|_{t=0}=0\\ \Longrightarrow C_2=\frac{k_1[A][B]S}{\phi_{protein}-k_{-1}-k_1} $$

Therefore:

$$ [EAB]=\frac{k_1[A][B]S}{-\phi_{protein}+k_{-1}+k_{1}}·(e^{-\phi_{protein}t}-e^{-(k_{-1}+k_{1})t})\\ \lambda_1=\phi_{protein}\\ \lambda_2=k_1+k_{-1}\\ \Longrightarrow [EAB]=C_2(e^{-\lambda_1t}-e^{-\lambda_2t}) $$

Finally we get:

$$ [M_{signal}]=k_2C_2(-\frac{e^{-\lambda_1t}}{\lambda_{1}}+\frac{e^{-\lambda_2t}}{\lambda_{2}})+C_3\\ \lim_{t\to\infty}[M_{signal}]=\frac{k_1k_2[A][B]}{\lambda_1\lambda_2}S $$

NOR Gate

This result can easily transit to NOR gate, because from mathematical perspective, the different places are initial values and another form of Hill equation. To describe the mechanism of NOR gate, we supposed that the whole system remains at steady state. (In other word, all concentrations remain as constants.)

$$ v_{inhibition}=V_{\max}-v_{generate}\\ \frac{d([mRNA])}{dt}=v_{inhibition}-\phi_{mRNA}[mRNA]=V_{\max}-v_{generate}-\phi_{mRNA}[mRNA]\\ \frac{d([protein])}{dt}=g_{protein}[mRNA]-\phi_{protein}[protein]\\ $$

We get:

$$ [mRNA]=\frac{v_{inhibition}}{\phi_{mRNA}}+C·e^{-\phi_{mRNA} t}\\ \Longrightarrow[mRNA]=A+B·e^{-\phi_{mRNA}t}\\ \lim_{t\to\infty}[mRNA]=\frac{v_{inhibition}}{\phi_{mRNA}} $$ $$ [protein]=g_{protein}\phi_{protein}^{-1}A+g_{protein}B(\phi_{protein}-\phi_{mRNA})^{-1}e^{-\phi_{mRNA} t}+C'e^{-\phi_{protein} t}\\ \Longrightarrow [protein]=A'+B'e^{-\phi_{mRNA} t}+C'e^{-\phi_{protein} t} $$

Furthermore we get:

$$ [M_{signal}]=A'·(k_1+k_2)^{-1}+B'·(-\phi_{mRNA}+k_1+k_2)^{-1}·e^{-\phi_{mRNA}t}\\+C'·(-\phi_{protein}+k_1+k_2)^{-1}·e^{-\phi_{protein}t}+D·e^{-(k_1+k_2)t} $$

With the relation:

$$ D=[M_{signal}]|_{t=0}-\{A'·(k_1+k_2)^{-1}+B'·(-\phi_{mRNA}+k_1+k_2)^{-1}\\+C'·(-\phi_{protein}+k_1+k_2)^{-1}\} $$

Also we have:

$$ \lim_{t\to\infty}[M_{signal}]=\frac{A'}{(k_1+k_2)}\\=\frac{g_{protein}}{(k_1+k_2)\phi_{protein}\phi_{mRNA}}·v_{inhibition} $$

Extra model

Review

Model of E.coli population:

$N_{E.coli}$ is satisfied to following equation:

$$ \frac{dN_{E.coli}}{dt}=rN_{E.coli}(1-\frac{N_{E.coli}}{N_{\max}})\\ N_{E.coli}=\frac{N_{\max}}{1+(\frac{N_{\max}}{N_{t=0}}-1)·e^{-rt}} $$

$r$ refers to growth rate of E.coli and $N_{\max}$ refers to the limits of E.coli population. Since $N_{\max}$ and $N_{t=0}$are constants, so we define following parameter:

$$ \frac{N_{\max}}{N_{t=0}}-1=N_{c} $$

Two hypothesis:

1. The signal molecule is toxic to E.coli, so the population will decrease related to the increase of concentration linearly.
2. The signal molecule induce the synthesis of GFP which occupy the substance that is originally used for growth. It indicates that if the GFP is produced, then the population will be at low level, otherwise the population will be at normal level.

Analyse

To show the difference between these two hypothesis, we give following equation:

Hypothesis 1:

$$ \frac{dN_{E.coli}}{dt}=rN_{E.coli}(1-\frac{N_{E.coli}}{N_{\max}})-\gamma N_{E.coli} $$

$\gamma$ refers to the death rate caused by toxic substance,$[S]$ refers to the concentration of signal molecule and $[S]_{critical}$ refers to the critical point which means all E.coli are dead:

$$ \gamma =r·(1-\frac{[S]}{[S]_{critical}}) $$

Therefore:

$$ \lim_{t\to+\infty}N_{E.coli}=N_{\max}·(1-\frac{[S]}{[S]_{critical}}) $$

X-axis refers to the time and Y-axis refers to the growth curves. Different curves refer to different concentrations.

Hypothesis 2:

$$ \frac{dN_{E.coli}}{dt}=r·N_{E.coli}(1-\frac{N_{E.coli}}{N_{\max}·\beta}) $$

$\beta$ refers to the ratio of limiting the growth of E.coli, which fits to following equation:

$$ \beta =1-\beta_{\lim}·\frac{[S]^n}{k^n+[S]^n} $$

The reason we use the efficiency of mRNA generation is because this ratio determines how many GFP will be finally produced. For example, if the ratio is high, the production of GFP will be at high level, which also means the most of substance are used to produce GFP instead of growth of E.coli.

Therefore:

$$ N_{E.coli}=\frac{N_{\max}·\beta}{1+(\frac{N_{\max}}{N_{t=0}}-1)·e^{-r t}} $$

X-axis refers to the time and Y-axis refers to the growth curves. Different curves refer to different concentrations. Low concentration refers to high population and high concentration refers to low population.

From our data, we found the result showed that hypothesis two was more realistic.

The experiment shows an obvious difference between low concentration and high concentration which fitts to the hypothesis two.

But we also cannot eliminate the hypthesis one, because the curves of low concentration go to steady state but the high concentration go slightly down. If we use hypothesis two to explain this phenomemon, that is: The production of GFP highly occupy the resource and leave little resource for the growth of E.coli even cannot mantain the population at the steady state. If we use hypothesis one, then the result is obvious that signal molecule is toxic to E.coli which causes unavoidable death of E.coli. So further study is required.

Model for Our Project

$$ \begin{align} \frac{{\mathrm{d}\left( {QS1R} \right)}}{{\mathrm{d}t}}&= {C_{QS1R}} + H\left( {{{\left[ M \right]}_e}} \right)\\ \frac{{\mathrm{d}\left( {\left[ {mRN{A_1}} \right]} \right)}}{{\mathrm{d}t}}&= \frac{{\mathrm{d}\left( {QS1R} \right)}}{{\mathrm{d}t}} - {\emptyset _1}\left[ {mRNA} \right]\\ \frac{{\mathrm{d}\left( {\left[ {protei{n_1}} \right]} \right)}}{{\mathrm{d}t}}&= {g_1}\left[ {mRN{A_1}} \right] - {\emptyset _2}\left[ {protei{n_1}} \right] - \frac{{\mathrm{d}\left( {\left[ M \right]} \right)}}{{\mathrm{d}t}}\\ {K_1}&= \frac{{\left[ {protei{n_1}} \right]\left[ {QS1 - AHL} \right]}}{{\left[ M \right]}}\\ H\left( {\left[ x \right]} \right)&= \frac{{{V_{\max }}{{\left[ x \right]}^m}}}{{{{\left[ x \right]}^m} + {K_a}^m}}\\ {\left[ M \right]_e}&= {P_e}\left[ M \right]\\ {P_{activated}}&= {P_e} + {P_e}'= P\left[ {QS1{R_{translated}}|{M_{combined}}} \right]\\ 1&= {P_{activated}} \cdot {P_{combined}} + {P_{cons}} \cdot \overline {{P_{combined}}} + {P_{inactivated}}\\ \frac{{\mathrm{d}\left( {CI} \right)}}{{\mathrm{d}t}}&= {C_{CI}} + H\left( {{{\left[ M \right]}_{e'}}} \right)\\ \frac{{\mathrm{d}\left( {\left[ {mRN{A_2}} \right]} \right)}}{{\mathrm{d}t}}&= \frac{{\mathrm{d}\left( {CI} \right)}}{{\mathrm{d}t}} - {\emptyset _3}\left[ {mRN{A_2}} \right]\\ \frac{{\mathrm{d}\left( {\left[ {protei{n_2}} \right]} \right)}}{{\mathrm{d}t}}&= {g_2}\left[ {mRN{A_2}} \right] - {\emptyset _4}\left[ {protei{n_2}} \right] - H\left( {{{\left[ {protei{n_2}} \right]}_e}} \right)\\ {\left[ {protei{n_2}} \right]_e}&= {P_{activated}}\left[ {protei{n_2}} \right]\\ 1&= {P_{activated}} \cdot {P_{combined}} + {P_{cons}} \cdot \overline {{P_{combined}}} + {P_{inactivated}}\\ \frac{{\mathrm{d}\left( {QS2I} \right)}}{{\mathrm{d}t}}&= {C_{QS2I}} + H'\left( {{{\left[ {protei{n_2}} \right]}_e}} \right)\\ \frac{{\mathrm{d}\left( {\left[ {mRN{A_3}} \right]} \right)}}{{\mathrm{d}t}}&= \frac{{\mathrm{d}\left( {QS2I} \right)}}{{\mathrm{d}t}} - {\emptyset _5}\left[ {mRN{A_3}} \right]\\ \frac{{\mathrm{d}\left( {\left[ {protei{n_3}} \right]} \right)}}{{\mathrm{d}t}}&= {g_3}\left[ {mRN{A_3}} \right] - {\emptyset _6}\left[ {protei{n_3}} \right] - \frac{{\mathrm{d}\left( {\left[ N \right]} \right)}}{{\mathrm{d}t}}\\ \frac{{\mathrm{d}\left( {CI} \right)}}{{\mathrm{d}t}}&= {C_{CI}} + H\left( {{{\left[ M \right]}_{e'}}} \right)\\ \frac{{\mathrm{d}\left( {\left[ {mRN{A_2}} \right]} \right)}}{{\mathrm{d}t}}&= \frac{{\mathrm{d}\left( {CI} \right)}}{{\mathrm{d}t}} - {\emptyset _3}\left[ {mRN{A_2}} \right]\\ \frac{{\mathrm{d}\left( {\left[ {protei{n_2}} \right]} \right)}}{{\mathrm{d}t}}&= {g_2}\left[ {mRN{A_2}} \right] - {\emptyset _4}\left[ {protei{n_2}} \right] - H\left( {{{\left[ {protei{n_2}} \right]}_e}} \right)\\ {\left[ {protei{n_2}} \right]_e}&= {P_{activated}}\left[ {protei{n_2}} \right]\\ 1&= {P_{activated}} \cdot {P_{combined}} + {P_{cons}} \cdot \overline {{P_{combined}}} + {P_{inactivated}}\\ \frac{{\mathrm{d}\left( {QS2I} \right)}}{{\mathrm{d}t}}&= {C_{QS2I}} + H'\left( {{{\left[ {protei{n_2}} \right]}_e}} \right)\\ \frac{{\mathrm{d}\left( {\left[ {mRN{A_3}} \right]} \right)}}{{\mathrm{d}t}}&= \frac{{\mathrm{d}\left( {QS2I} \right)}}{{\mathrm{d}t}} - {\emptyset _5}\left[ {mRN{A_3}} \right]\\ \frac{{\mathrm{d}\left( {\left[ {protei{n_3}} \right]} \right)}}{{\mathrm{d}t}}&= {g_3}\left[ {mRN{A_3}} \right] - {\emptyset _6}\left[ {protei{n_3}} \right] - \frac{{\mathrm{d}\left( {\left[ N \right]} \right)}}{{\mathrm{d}t}} \end{align} $$