Difference between revisions of "Team:KU Leuven/Modeling/Internal"

Line 109: Line 109:
 
</ol>
 
</ol>
 
</p>
 
</p>
−
</div>
 
 
</div>
 
</div>
 
</div>
 
</div>

Revision as of 16:41, 18 September 2015

Logo

Internal Model

Introduction

We can think of many relevant questions when implementing a new circuit: how sensitive is the system, how much will it produce and will it affect the growth? As such, it is important to model the effect of the new circuits on the bacteria. This will be done in the Internal Model. We will use two approaches. First we will use a bottom-up approach. This involves building a detailed kinetic model with rate laws. We will use Simbiology and ODEs to study the sensitivity and dynamic processes inside the cell. Afterwards, a top-down model, Flux Balance Analysis (FBA), will be used to study the steady-state values for production flux and growth rate. This part is executed by the iGEM Team of Toulouse as part of a collaboration and can be found here

Simbiology and ODEs

In the next section we will describe our Simbiology model. Simbiology is a toolbox from Matlab designed for the simulation of (bio)chemical reactions. It allows us to calculate systems of ODEs and to visualize the system in a diagram. It also has options to make scans for different parameters, which allows us to study the effect of the specified parameter. We will focus on the production of leucine, Ag43 and AHL in cell A and the changing behavior of cell B due to changing AHL concentration. In this perspective, we will make two models in Simbiology: one for cell A and one for cell B. First we will describe how we made the model and searched for the parameters. Afterwards we check the robustness of the model with a parameter analysis and we do scans to check for the effects of molecular noise.

Quest for parameters

We can divide the different processes that are being executed in the cells in 7 classes: transcription, translation, DNA binding, complexation and dimerization, protein production kinetics, degradation and diffusion. We went on to search the necessary parameters and descriptions for each of these categories. To start making our model we have to pick a unit. We choose to use molecules as unit, because many constants are expressed in this unit and it allows us to drop the dillution terms connected to cell growth. We will also work with a deterministic model instead of a stochastic model. A stochastic model would show us the molecular noise, but we will check this with parameter scans.

The next step is to make some assumptions:

  1. The effects of cell division can be neglected
  2. The substrate pool can not be depleted and the concentration (or amount of molecules) of substrate in the cell is constant
  3. The exterior of the cell contains no leucine at t=0 and is perfectly mixed
  4. Diffusion happens independent of cell movement and has a constant rate

Transcription

Transcription is the first step of gene expression. It involves the binding of RNA polymerase to the promoter region and the formation of mRNA. The transcription rate is dependent on a number of conditions like the promoter strength and the used strain. Our system has both constitutive promoters, which are always active, and inducable promoters, which can be activated or repressed.

Maximum Transcription rate

First we will try to find the maximum transcription rate. The prediction of the transcription rate has been an important hold-back in the past, even though Polymerases per Second was introduced as a unit. This is the amount of Polymerases that passes through a given position in the DNA per time unit and is essentially the transcription rate at a particular location on the DNA (open wetware http://openwetware.org/wiki/PoPS). We can use this value as the transciption rate of the whole gene because it is the value of the slowest and rate-defining step, the binding of polymerase on the promoter. The other step being the movement of polymerase over the gene. It remained difficult to measure in vivo but Kelly et al have introduced a way to measure the activity of promoters using an in vivo standard, promoter J23101. This gave rise to a new unit: Relative activity of promoter (RPU)$=\frac{PoPS_{phi}}{PoPS_{J23101}}$. By using relative units, the variability due to equipment and conditions was drastically decreased. The PoPS of J23101 was found to be 0.03.

As our constitutive promoter, we used J23114. On the iGEM website we found that the whole J constitutive promoter family had been characterized and we found that J23101 had a strength of 1791 au and J23114 had a strength of 256 au with J23112 as standard in this measurement. We conclude that J23114 has a PoPS of 0.00429.

The strength of the pLuxPR promoter has been measured relatively to pLacIQ by iGEM Tokyo_Tech 2010. They found a value for pLuxPR of 550 and for pLacIQ of 390. pLacIQ has also been measured by iGEM Upssala university 2012 and was found to have a RPU of 0.04 giving it a PoPS of 0.0012. From these values we calculate the PoPS of pLuxPR to be 550/390 $\cdot$ 0.0012 =0.00169 PoPS.

Next we sought values for pCI. iGEM NYMU-Taipei 2009 measured the pLux promoter relative to pCI, which allows us to calculate the PoPS for pCI: $\frac{1}{0.43} {\cdot} 0.00169=0.0394$ PoPS. Their results of pPen were not consistent and the team mentioned they lacked experience in the beginning. Results of the earlier used pCI and pLux were consistent though. This is why we only used their last results and took the mean of pPen compared to pTet: $(0.801 {\cdot} 0.537 + 0,191) {\cdot} {\frac{1}{2}}=0.3106$. Using this mean value we were able to calculate the PoPS: $0.3106 {\cdot} 0.0394=0.01224$.

These values all depend on measuring activities and since the strains, media and antibiotic markers won’t be completely the same, the real values will diverge from the values found with these simple calculations. Nevertheless, these values can give us an idea about the relative strength of the different promoters. We could gather the real values of mRNA concentration in our cells by executing a qPCR, but for now we will use the found values for our model. The last important parameter for transcription is the copy number in which the genes are present in the cell. The plasmid that incorporates the used genes, has an ORI with a copy number between 100 and 300. In our model we will use a mean copy number of 200.

Inducable promoters

For inducable promoters the maximum transcription rate has to be multiplied with a factor to integrate the effect of the transcription factors. Our system includes two types of DNA binding proteins: repressors and activators. These proteins bind certain regions in the promoter region of the DNA. This binding can either repress or activate the promoter activity, affecting transcription rate. The repressors used in our system are the cI protein of the lambda phage and the PenI protein. The activator protein is LuxR, which should be activated by AHL. All proteins first need to form homodimers before being able to bind their target DNA.

DNA binding will be simulated using the Hill function : ${\theta} = \frac{{[L]}^n}{{K_d}^n + [L]^n} = \frac{1}{1+(\frac{K_{d}}{[L]})^{n}}$. $\theta$ is the amount of DNA bound by the protein, $L$ is the amount of protein, $K_d$ is the dissociation constant and $n$ is the Hill coefficient. This Hill function has values between 0 and 1 and acts like an ON/OFF switch.
For repressors we are interested in how much of the DNA is still unbound and active: $\frac{1}{1 + ({\frac{[Repressor]}{K_{d}}})^{n}}$
For activators we are interested in how much of the DNA is bound and active: $\frac{1}{1+(\frac{K_{d}}{[Activator]})^{n}}$
The Hill coefficient gives an idea about the strength about the interaction between the DNA and the protein and is sometimes used to estimate the amount of proteins that bind the DNA but caution is needed in using it. Indeed, the Hill function has knwown shortcomings, but it is still very useful because of its simplicity (citate paper).

The constants needed for our model are $K_d$ and $n$. :

Table 1: The constants of the Hill function for cI, LuxR and PenI

Protein Kd (molecules) n Source
cI 20 2 2009 iGEM Aberdeen and Wang et al. (2009)
LuxR/AHL 3.7 1 2007 iGEM ETH Zurich and Basu et al. (2005)
PenI 17 2 Wittman et al. (1993)

${K_d}$ is the value which corresponds with the amount of repressor or activator necessary to repress or activate half of the promoters. Since our copynumber is 200, we multiply the Hill constant with 200 because this higher value seems more logical to us.

Leakiness

A last important value in transcription is the leakiness of the promoter. A shutdown promoter still has a very small transcription which is called the leakiness of the promoter. This value should be as low as possible, because there should be a big difference between ON and OFF promoters. The effect of leakiness of promoters is something that should be checked.

iGEM 2012 Upssala measured the leakiness of the pLuxPR promoter. Without AHL and LuxR, the pLuxPR was active and had a RPU of 10. This means that the leaking promoter has a PoPS of $\frac{10}{390} {\cdot} 0.040 {\cdot} 0.03 = 3.077 {\cdot} 10^{-5} PoPS$. We can use this value to estimate the leakiness of the other promoters.


Translation

The next step in gene expression is the translation of the synthesized mRNA to a protein. This involves the binding of a ribosome complex to the RBS site of the mRNA, elongation and termination. The initiation step is the rate-limiting step of the three. The translation rate is thus tied to the thermodynamics of RBS-ribosome interaction but other sequences (5’-UTR, Shine-Delgarno,repeats, internal startcodons) and possible secondary structures (hairpins) also play a role (Chiam Yu Ng, Salis). (Fong) Salis et al have developed an algorithm to predict the translation initiation rate based upon thermodynamic calculations. This algorithm is used in the Ribosome Binding Site Calculator (https://salislab.net/software/ , Borujeni & Salis). We used this Calculator to predict the translation rates of our mRNAs.

The results are in table 1. The program gives us results in au (arbitrary units). Since we know the translation rate of LuxI and LuxR, we can use these as a base to calculate the other translation rates since the used scale is proportional. (table 1, column 2) LuxI and LuxR have almost the same output from the calculator which corresponds with the paper where they also have the same translation rate. (Ag43 has a very low translation rate, which is not completely illogical, since Ag43 is by far the biggest protein.) Our values are normal values since 1000 is a moderate value and values between 1 and 100 000 are possible. (efficient search, mapping and optimization, salis). We did get warnings about the prediction (NEQ: not at equilibrium) which happens when mRNA may not fold quickly to its equilibrium state.

Table 2: The translation rates for the different proteins as found by the RBS calculator using LuxR and LuxI as standard

Protein Translation rate (a.u.) Translation rate (1/s)
LuxI 261.13 0.016
LuxR 279.61 0.016
cI857 5695.79 0.35
Transaminase B 2195.1 0.13
Ag43-YFP 33.22 0.0020
PenI 6696.73 0.41
RFP 2604.51 0.16
CheZ-GFP 396.95 0.024

For TransaminaseB we also need to add a production term because it is an protein which naturally occurs in E. coli. We choose a value of 0.0033 for this constitutive production rate.

Complexation and oligomerization kinetics

Before the proteins can bind the promoter region, they first have to make complexes. cI and penI form homodimers, while LuxR first forms a heterodimer with AHL and afterwards forms a homodimer with another LuxR/AHL dimer. This next part will describe the kinetics of such complexation. We only need the parameters of LuxR and cI, because the Hill function found for penI already was adapted for the protein in monomer form.

Table 3: Dissociation and association rates for the oligomerization of the transcription factors used in our circuit

Biomolecule Oligomerization rate Source
LuxR/AHL association $1{\cdot}10^{-5}$ 1/(molecule*s) Goryachev et al. (2005)
LuxR/AHL dissociation 0.00333 1/s Goryachev et al. (2005)
LuxRdimer association $1{\cdot}10^{-5}$ 1/(molecule*s) Goryachev et al. (2005)
LuxRdimer dissociation 0.01 1/s Goryachev et al. (2005)
cI dimer association 0.00147 1/(molecule*s) Samiee et al. (2005)
cI dimer dissociation 0.01 1/s Samiee et al. (2005)

Protein production kinetics

Our models comprise two product forming enzymes: LuxI and Transaminase B. LuxI is the enzyme responsible for AHL production and Transaminase B is responsible for Leucine formation.

For LuxI kinetics Schaefer et al. (1996) found a value of 1 molecule/minute for the Vmax. Since we assume that LuxI is fully saturated with substrate, we take this value as the synthesis rate of AHL.

Transaminase B forms leucine from aKIC and glutamate. In this reaction glutamate is converted to aKG. The reaction kinetics is of the ping-pong bi bi form. This model describes a mechanism in which the binding of substrates and release of products is ordered. The enzyme shuttles between a free and a substrate-modified intermediate state.

Ping-Pong Bi-Bi equations

\begin{align} \frac{{\large d}{TB}}{d t}= & \beta_{TB} {\cdot} {m_{ilvE}} - {kf}_{1}{\cdot}{TB}{\cdot}{Glu} + {kf}_{-1}{\cdot}{[TB-GLU]} - {kr}_1{\cdot}{Leucine}_{in}{\cdot}{[TB-GLU]} \\\\ & + {kr}_{-1}{\cdot}{[TB-Leu]} + {kcat2}{\cdot}{[{TBNH}_2-aKIC]} + {kcat4}{\cdot}{[{TBNH}_2-aKG]} - d_{TB}{\cdot}{TB} \end{align} $$\frac{{\large d}{[TB-GLU]}}{d t}= -{kcat1}{\cdot}{[TB-GLU]} + {kf}_{1}{\cdot}{TB}{\cdot}{Glu} - {kf}_{-1}{\cdot}{[TB-GLU]} $$ \begin{align} \frac{{\large d}{[{TBNH}_2]}}{d t}= & {kcat1}{\cdot}{[TB-GLU]} + {kcat3}{\cdot}{[TB-Leu]}- {kf}_{2}{\cdot}{{TBNH}_2}{\cdot}{aKIC} + {kf}_{-2}{\cdot}{[{TBNH}_2-aKIC]} \\\\ & - {kr}_2{\cdot}{{TBNH}_2}{\cdot}{aKG} +{kr}_{2}{\cdot}{[{TBNH}_2-aKG]} \end{align} $$\frac{{\large d}{[{TBNH}_2-aKIC]}}{d t}= -{kcat2}{\cdot}{[{TBNH}_2-aKIC]} + {kf}_{2}{\cdot}{{TBNH}_2}{\cdot}{aKIC} - {kf}_{-2}{\cdot}{[{TBNH}_2-aKIC]} $$ $$\frac{{\large d}{[TB-Leu]}}{d t}= {kr}_1{\cdot}{Leucine}_{in}{\cdot}{[TB-GLU]} - {kr}_{-1}{\cdot}{[TB-Leu]} - {kcat3}{\cdot}{[TB-Leu]} $$ $$\frac{{\large d}{[{TBNH}_2-aKG]}}{d t}= {kr}_2{\cdot}{{TBNH}_2}{\cdot}{aKG} - {kr}_{2}{\cdot}{[{TBNH}_2-aKG]} - {kcat4}{\cdot}{[{TBNH}_2-aKG]}$$


Yang et al. (2005) studied these reactions and they found a method to estimate the constants: the Lambda approximation. Their constants of this reversible Ping-Pong Bi-Bi reactions are put in the next table :

Table 4: Constants of the Ping-Pong Bi-Bi reactions using Transaminase B as found by Yang et al. (2005)

Constant Value
${kf}_{1}$ 0.0098
${kf}_{-1}$ 3300
${kf}_{2}$ 0.0882
${kf}_{-2}$ 5940
${kr}_{1}$ 0.0031184
${kr}_{-1}$ 4620
${kr}_{2}$ 0.0041161
${kr}_{-2}$ 3465
$kcat1$ 33.33
$kcat2$ 60
$kcat3$ 46.67
$kcat4$ 35

Degradation

Proteins, mRNA and other metabolites have a turnover rate. They are degraded over time. The degradation rate will be described with first order kinetics with d as the degradation constant. Not every molecule has the same degradation rate, since some molecules are more stable than others. We can influence the stability of the molecules. For example, in cell B it is important that there is a fast switch between conditions and a fast turnover of CheZ and RFP is necessary. This is why we add a LVA-tag to these proteins. This tag destabilizes the protein and makes them degradade faster. For TransaminaseB we choose a very high degradation rate, because we did not include degradation terms for the TransaminaseB bound to substrate. The degradation rates used in the model are put in the next table:

Table 5: Degradation rates for the different biomolecules used in our circuit

Biomolecule Degradation rate (1/s) Source
LuxI $5 \cdot 10^{-5}$ Goryachev et al. (2005)
luxI-mRNA 0.006 Goryachev et al. (2005)
LuxR $2 \cdot 10^{-4}$ Goryachev et al. (2005)
luxR-mRNA 0.006 Goryachev et al. (2005)
cI857 0.00288 iGEM Aberdeen 2009
cI-mRNA 0.002265 Bernstein et al. (2002)
Transaminase B 4.5 Estimated
ilvE-mRNA 0.00304 Bernstein et al. (2002)
Ag43-YFP $2.84{\cdot}10^{-5}$ Bernstein et al. (2002)
Adhesine mRNA 0.00186 Bernstein et al. (2002)
PenI 0.00333 Estimated
PenI mRNA 0.002265 Bernstein et al. (2002)
RFP-LVA 0.0002814 Andersen et al. (1998) and iGEM KULeuven 2008
RFP mRNA 0.066 Bernstein et al. (2002)
CheZ-GFP-LVA 0.0002814 Andersen et al. (1998)
CheZ-GFP mRNA 0.00333 Estimated
AHL $1.667 \cdot 10^{-4}$ Basu et al. (2005)
AHL external $1.155 \cdot 10^{-5}$ Horswill et al. (2007)
Leucine 0.0167 Estimated
Leucine external 0.000333 Estimated

Diffusion

Our model has 2 types of diffusion. Diffusion from the inside of the cell to the outside over the cell membrane and diffusion in the external medium. The diffusion over the cell membrane is more complicated because some proteins play a role in it and the membrane is not equally permeable for every molecule. For AHL, we found a value given in molecules/second. This unit seems strange because diffusion is usually used with units in concentration. It also leads to strange results since the amount of molecules is being leveled out so the inside amount equals the outside amount even though the concentrations are different. This is why we added a correction for the volume of the cell and the external volume to it.

We assume that the volume of the cell has a shape of a cilinder and is constant. For this simplified volume Goryachev et al. found a value of $5.65 {\cdot} 10^{-16} $ l. We can take this volume as a constant, since cell growth is very small compared to the diffusion. The outside compartment will be modeled as half a sphere with the cell as center and a radius equal to $\sqrt{2{\cdot}D{\cdot}t} + r_0$. We only take half a sphere because there is no upward diffusion. For the initial value of the outside compartment we take a volume (and thus $r_0$ slightly bigger than the cell volume (radius). We choose a value of $5.7 {cdot} 10^{-16}$l which accords to a radius of $5.128{\cdot}10^{-6}$ dm. Because there is also already a volume taken by the cell, the total initial outside compartment has an initial value of $5{\cdot}10^{-18}$ dm. The volume will increase per time step with a $\frac{2{\cdot}\pi{\cdot}D}{\sqrt{2{\cdot}D{\cdot}t}}{\cdot}{(\sqrt{2{\cdot}D{\cdot}t} + r_0)^2}$

With this approximation we get more logical results. We take the diffusion rate for AHL equal to 0.23 1/s. For Leucine we make a difference between inward and outward diffusion. Inward diffusion is facilitated by transporters, so we choose this value to be the largest. The inward diffusion rate is estimated 0.05 and the outward diffusion 0.0005.

Table 6: Diffusion rates of the molecules in our circuit

Biomolecule Diffusion rate (1/s) Source
AHL IN and OUT 0.23 1/s Goryachev et al. (2005)
Leucine IN 0.05 1/s Estimated
Leucine OUT 0.0005 1/s Estimated
AHL extracellular $7.1{\cdot}10^{-9}$ dm²/s Goryachev et al. (2005)
Leucine extracellular $7.3{\cdot}10^{-8}{\cdot}0.09$ dm²/s iGEM Aberdeen 2009

System

After this extensive literature search, we can finally set up our complete system of ODEs for every cell.

Cell A

The designed circuit in Cell A is under control of a temperature sensitive cI repressor. Upon raising the temperature, cI will dissociate from the promoter and the circuit is activated. This leads to the initiation of the production of LuxR and LuxI. LuxI will consecutively produce AHL, which binds with LuxR. The newly formed complex will then activate the production of Leucine and Ag43. Leucine and AHL are also able to diffuse out of the cell into the medium. Ag43 is the adhesine which aids the aggregation of cells A, while Leucine and AHL are necessary to repel cells B.

We can extract the following ODE's from this circuit:

Cell A equations

Symbols:${}$ ${\alpha}$: transcription term, ${\beta}$: translation term, $d$: degradation term,
$D$: diffusion term, ${ K_d}$: dissociation constant, n: Hill coefficient, L: leak term

$$\frac{{\large d} m_{cI}}{d t} = \alpha_1 {\cdot} cI_{gene} - d_{mCI} {\cdot} m_{cI}$$ \begin{align} \frac{{\large d}{cI}}{d t} = \beta_{cI} {\cdot} {m_{cI}} -2 {\cdot} {k_{cI,dim}} {\cdot} {cI}^2 + 2 {\cdot} {k_{-cI,dim}}{\cdot} {[cI]_2} - d_{cI} {\cdot} {cI} \end{align} $$\frac{{\large d}{[cI]_2}}{d t}= k_{cI,dim} {\cdot} {cI}^2 - {k_{-cI,dim}}{\cdot} {[cI]_2} $$ $$\frac{{\large d} m_{LuxI}}{d t} = (L_{lambda} + {\frac{\alpha_{lambda}}{1 + ({\frac{[cI]_2}{K_{d1}}})^{n_{cI}}}}) {\cdot} LuxI_{gene} - d_{mLuxI} {\cdot} m_{LuxI} $$ $$\frac{{\large d} m_{LuxR}}{d t} = (L_{lambda} + {\frac{\alpha_{lambda}}{1 + ({\frac{[cI]_2}{K_{d1}}})^{n_{cI}}}}) {\cdot} LuxR_{gene} - d_{mLuxR} {\cdot} m_{LuxR} $$ $$\frac{{\large d} LuxI}{d t} = \beta_{LuxI} {\cdot} {m_{LuxI}} - d_{LuxI} {\cdot}{LuxI} $$ $$\frac{{\large d} LuxR}{d t} = \beta_{LuxR} {\cdot} {m_{LuxR}} -k_{lux,as} {\cdot}{LuxR}{\cdot}{AHL_{in}} + k_{lux,dis}{\cdot}{[LuxR/AHL]} - d_{LuxR} {\cdot}{LuxR} $$ $$\frac{{\large d} AHL_{in}}{d t} = {k_{luxI}} {\cdot} {luxI} - k_{lux,as} {\cdot}{luxR}{\cdot}{AHL_{in}} + k_{lux,dis}{\cdot}{[luxR/AHL]}+ ( {D_{IN,AHL}} {\cdot} {AHL_{out}} {\cdot}\frac{{V_{cell}}}{V_{external,AHL}} - {D_{OUT,AHL}} {\cdot} {AHL_{in}} ) - d_{AHL,in} {\cdot} {AHL_{in}} $$ $$\frac{{\large d} AHL_{out}}{d t} = ( {D_{OUT,AHL}} {\cdot} {AHL_{in}} {\cdot}\frac{V_{external,AHL}}{{V_{cell}}}- {D_{IN,AHL}}{\cdot}{AHL_{out}} ) -d_{AHL,out}{\cdot}{AHL_{out}} $$ $$\frac{{\large d} [luxR/AHL]}{d t} = k_{lux,as} {\cdot}{luxR}{\cdot}{AHL_{in}} - k_{lux,dis}{\cdot}{[luxR/AHL]} - 2 {\cdot} k_{lux,dim} {\cdot}{[luxR/AHL]^2} + 2 {\cdot}{k_{-lux,dim}}{\cdot}{[luxR/AHL]_{2}} $$ $$\frac{{\large d} [luxR/AHL]_{2}}{d t} = k_{lux,dim} {\cdot}{[luxR/AHL]^2} - k_{- lux,dim} {[luxR/AHL]} $$ $$\frac{{\large d} m_{ilvE}}{d t} = (L_{lux} + \frac{\alpha_{lux}}{1+(\frac{K_{d2}}{[luxR/AHL]_{2}})^{n_{lux}}} ) {\cdot} ilvE_{gene} - d_{milvE} {\cdot} {m_{ilvE}} $$ $$\frac{{\large d} m_{Ag43}}{d t} = (L_{lux} + \frac{\alpha_{lux}}{1+(\frac{K_{d2}}{luxR/AHL]_{2}})^{n_{lux}}} ) {\cdot} Ag43_{gene} - d_{mAg43} {\cdot} {m_{Ag43}} $$ $$\frac{{\large d} Ag43}{d t} = \beta_{Ag43} {\cdot} {m_{Ag43}} - d_{Ag43} {\cdot} {Ag43} $$ $$\frac{{\large d} Transaminase B}{d t} = \beta_{TB} {\cdot} {m_{ilvE}} - kf_1 {\cdot} {Transaminase B} {\cdot}{Glutamate} + kf_{-1}{\cdot}{[TB-GLU]} - kr_1 {\cdot}{Leucine_{in}}{\cdot}{Transaminase B} + kr_{-1}{\cdot}{TB-Leu} + kcat2{\cdot}{[{TBNH}_2-aKIC]} + kcat4{\cdot}{[{TBNH}_2-aKG]} + k_{production} - d_{TB} {\cdot} {Transaminase B}$$ \begin{align} \frac{{\large d}{TB}}{d t}= & \beta_{TB} {\cdot} {m_{ilvE}} - {kf}_{1}{\cdot}{TB}{\cdot}{Glu} + {kf}_{-1}{\cdot}{[TB-GLU]} - {kr}_1{\cdot}{Leucine}_{in}{\cdot}{[TB-GLU]} \\\\ & + {kr}_{-1}{\cdot}{[TB-Leu]} + {kcat2}{\cdot}{[{TBNH}_2-aKIC]} + {kcat4}{\cdot}{[{TBNH}_2-aKG]} - d_{TB}{\cdot}{TB} \end{align} $$\frac{{\large d}{[TB-GLU]}}{d t}= -{kcat1}{\cdot}{[TB-GLU]} + {kf}_{1}{\cdot}{TB}{\cdot}{Glu} - {kf}_{-1}{\cdot}{[TB-GLU]} $$ \begin{align} \frac{{\large d}{[{TBNH}_2]}}{d t}= & {kcat1}{\cdot}{[TB-GLU]} + {kcat3}{\cdot}{[TB-Leu]}- {kf}_{2}{\cdot}{{TBNH}_2}{\cdot}{aKIC} + {kf}_{-2}{\cdot}{[{TBNH}_2-aKIC]} \\\\ & - {kr}_2{\cdot}{{TBNH}_2}{\cdot}{aKG} +{kr}_{2}{\cdot}{[{TBNH}_2-aKG]} \end{align} $$\frac{{\large d}{[{TBNH}_2-aKIC]}}{d t}= -{kcat2}{\cdot}{[{TBNH}_2-aKIC]} + {kf}_{2}{\cdot}{{TBNH}_2}{\cdot}{aKIC} - {kf}_{-2}{\cdot}{[{TBNH}_2-aKIC]} $$ $$\frac{{\large d}{[TB-Leu]}}{d t}= {kr}_1{\cdot}{Leucine}_{in}{\cdot}{[TB-GLU]} - {kr}_{-1}{\cdot}{[TB-Leu]} - {kcat3}{\cdot}{[TB-Leu]} $$ $$\frac{{\large d}{[{TBNH}_2-aKG]}}{d t}= {kr}_2{\cdot}{{TBNH}_2}{\cdot}{aKG} - {kr}_{2}{\cdot}{[{TBNH}_2-aKG]} - {kcat4}{\cdot}{[{TBNH}_2-aKG]}$$ $$\frac{{\large d}{Leucine_{in}}}{d t}= {kcat2}{\cdot}{[{TBNH}_2-aKIC]} - {kr}_1 {\cdot}{Leucine_{in}}{\cdot}{Transaminase B} + {kr}_{-1}{\cdot}{[TB-Leu]} - d_{Leu}{\cdot}{Leucine_{in}} - {D_{OUT,Leu}}{\cdot}{Leucine_{in}} + {D_{IN,Leu}}{\cdot}{Leucine_{out}}{\cdot}\frac{{V_{cell}}}{V_{external,leu}} $$ $$\frac{{\large d} Leucine_{out}}{d t} = (D_{OUT,Leu} {\cdot}{leucine_{in}}{\cdot}\frac{V_{external,leu}}{{V_{cell}}} - D_{IN,Leu} {\cdot}{leucine_{out}}) - d_{Leu,out} {\cdot} {Leucine_{out}} $$ $$\frac{{\large d} V_{external,AHL}}{d t} = \frac{2{\cdot}\pi{\cdot}D_{AHL}}{\sqrt{2{\cdot}D{\cdot}t}}{\cdot}{(\sqrt{2{\cdot}D{\cdot}t} + r_0)^2}$$ $$\frac{{\large d} V_{external,Leu}}{d t} = \frac{2{\cdot}\pi{\cdot}D_{Leu}}{\sqrt{2{\cdot}D{\cdot}t}}{\cdot}{(\sqrt{2{\cdot}D{\cdot}t} + r_0)^2}$$



We visualize these ODE's in the Simbiology Toolbox which results in the following diagram:

Cell A

Figure 1: The internal network of cell A visualized by Simbiology


Cell B

The system of Cell B is also under control of the cI repressor and is activated similar as cell A. The activation by the temperature raise, leads to the production of LuxR. AHL of the medium can diffuse into the cell, binding LuxR and activating the next component of the circuit. This leads to the production of CheZ and PenI. CheZ is the protein responsible for cells to make a directed movement, governed by the repellent Leucine. PenI is a repressor which will shut down the last part of the circuit which was responsible for the production of RFP.

We can extract the following ODEs for Cell B from this sytem:

Cell B equations

Symbols: ${\alpha}$:transcription term, ${\beta}$:translation term, $d$:degradation term,
$D$:diffusion term, ${ K_d}$:dissociation constant, n:Hill coefficient, L:leak term

$$\frac{{\large d} m_{cI}}{d t} = \alpha_1 {\cdot} cI_{gene} - d_1 {\cdot} m_{cI}$$ $$\frac{{\large d}{cI}}{d t} = \beta_1 {\cdot} {cI} -2 {\cdot} {k_{cI,dim}} {\cdot} {cI}^2 + 2 {\cdot}{k_{-cI,dim}}{\cdot} {[cI]_2} - d_{cI} {\cdot} {cI} $$

$$\frac{{\large d}{[cI]_2}}{d t}= k_{cI,dim} {\cdot} {cI}^2 - {k_{-cI,dim}}{\cdot} {[cI]_2} $$ $$\frac{{\large d} m_{luxR}}{d t} = (L_{lambda} + {\frac{\alpha_{lambda}}{1 + ({\frac{cI}{K_{d1}}})^{n_{cI}}}}) {\cdot} luxR_{gene} - d_{mluxR} {\cdot} m_{luxR} $$ $$\frac{{\large d} luxR}{d t} = \beta_{luxR} {\cdot} {m_{luxR}} - d_{luxR} {\cdot}{luxR} $$ $$\frac{{\large d} AHL_{in}}{d t} = ( {D_{IN,AHL}} {\cdot} {AHL_{out}} - {D_{OUT,AHL}} {\cdot} {AHL_{in}} ) - k_{lux,as}{\cdot}{luxR}{\cdot}{AHL_{in}} + k_{lux,dis}{\cdot}{[luxR/AHL]} - d_{AHL,in} {\cdot} {AHL_{in}} $$ $$\frac{{\large d} AHL_{out}}{d t} = ( {D_{OUT,AHL}} {\cdot} {AHL_{in}} - {D_{IN,AHL}}{\cdot}{AHL_{out}} ) -d_{AHL,out}{\cdot}{AHL_{out}} $$ $$\frac{{\large d} [luxR/AHL]}{d t} = k_{lux,as} {\cdot}{luxR}{\cdot}{AHL_{in}} - k_{lux,dis}{\cdot}{[luxR/AHL]} - 2 {\cdot} k_{lux,dim} {\cdot}{[luxR/AHL]^2} + 2 {\cdot}{[luxR/AHL]_{2}} $$ $$\frac{{\large d} [luxR/AHL]_{2}}{d t} = k_{lux,dim} {\cdot}{[luxR/AHL]^2} - k_{- lux,dim} {[luxR/AHL]} $$ $$\frac{{\large d} m_{CheZ}}{d t} = (L_{lux} + \frac{\alpha_{lux}}{1+(\frac{K_{d2}}{[luxR/AHL]_{2}})^{n_{lux}}} ) {\cdot} CheZ_{gene} - d_{mCheZ} {\cdot} {m_{CheZ}} $$ $$\frac{{\large d} m_{PenI}}{d t} = (L_{lux} + \frac{\alpha_{lux}}{1+(\frac{K_{d2}}{luxR/AHL]_{2}})^{n_{lux}}} ) {\cdot} PenI_{gene} - d_{mPenI} {\cdot} {m_{PenI}} $$ $$\frac{{\large d} CheZ}{d t} = \beta_{CheZ} {\cdot} {m_{CheZ}} - d_{CheZ} {\cdot} {CheZ} $$ $$\frac{{\large d} PenI}{d t} = \beta_{PenI} {\cdot} {m_{PenI}} - d_{PenI} {\cdot} {PenI} $$ $$\frac{{\large d} m_{RFP}}{d t} = (L_{Pen} + {\frac{\alpha_{Pen}}{1 + ({\frac{PenI}{K_{d3}}})^{n_{Pen}}}}) {\cdot} RFP_{gene} - d_{mRFP} {\cdot} m_{RFP} $$ $$\frac{{\large d} RFP}{d t} = \beta_{RFP} {\cdot}{m_{RFP}} - d_{RFP} {\cdot}{RFP} $$



We visualize these ODE's in the Simbiology toolbox. This gives us the following diagrams:

Cell B

Figure 2: The internal network of cell A visualized by Simbiology

Results

For cell A we made a simulation with cell A in the ON and OFF mode as visuable in figure 3. When cell A is in the OFF state, the whole designed circuit is in OFF mode. This means that the cI repressor is succesful in repressing the design. If the degradation rate of cI is raised to simulate the temperature rising, we see that all the components of the system show a big increase. For LuxR this increase is only temporary, but this is also explainable. Since LuxI keeps on producing AHL which binds LuxR to form a complex. Indeed all the AHL reacts to form the complex. Some values seem really high (for example the AHL,out en Leucine, extracellular values but they are also in the biggest volume so the concentration is not that high). We also assumed that there is a substrate pool without limiting, which is of course not the real situation. The most important conclusion is that with values backed up by literature, our system qualitatively still shows the desired behavior.


Simulation of all processes in Cell A

Figure 3: Simulation of all processes in Cell A in ON and OFF state

In the OFF state simulations, there is not a big difference between the red and green lines. We do see a very small rise in LuxR but we can ignore this because it is so small. We see a fast equilibration between the external AHL and internal AHL and no drop since there is no LuxR to react with the AHL. The only proteins that are available in high amounts are PenI and RFP. The high amount of PenI are not predicted in our design, but it does not affect the amount of RFP.

In the ON simulations we see a big difference between the red and green lines. When there is AHL available, the production of CheZ and PenI is much bigger and the production of RFP much lower. This is the behavior we wished for. Thus, our system is still qualitatively showing the desired behavior.



Simulation of all processes in Cell B in OFF

Figure 4: Simulation of all processes in Cell B in OFF state with and without AHL induction



Simulation of all processes in Cell B in ON

Figure 5: Simulation of all processes in Cell B in ON state with and without AHL induction



Sensitivities

Now we are going to check which parameters have the biggest effect on the output and are the most important. We can quantify this effect using derivatives: $\frac{\delta {output}}{\delta {parameter}}$. The parameters with the largest sensitivity value, are the parameters that should be best characterized. Furthermore, if they are controllable, the could be varied to our wishes. The sensitivity analysis will be executed in Simbiology. This analysis uses "complex-step approximation" to calculate derivatives of reaction rates. This technique yields accurate results for the vast majority of typical reaction kinetics. We will use full dedimensionalization. This way we can compare the results.

For cell A the output is the Medium Leucine, Adhesine and Medium AHL. We took a time integral of the sensitivity and plotted it in figure 6.

The leucine medium is dependent on Diffusion of Leucine out, leucine consumption, degradation, transcription and translation of ilvE mRNA. We also notice that not all Ping-Pong Bi Bi constants are equally important. kf1, kf-1 and kcat are the most important while kr1, kr-1 and kcat3 are the least important.

AHL medium: is highly sensitive for variations in luxI translation, degradation of LuxI and LuxI mRNA and for the catalytic activity of LuxI.

Adhesine-YFP is sensitive for translation, degradation and transcription.

In Cell A are all outputs dependent on variables closely related to themselves. Indeed, since there is production of AHL and LuxR in the cell, the transcriptional network is always active so the steps concerning the transcriptional network lose their importance. It is thus only a matter of understanding the metabolism, production and degradation terms to correctly model Cell A.

Sensitivity of cell A

Figure 6: Sensitivity analysis of parameters in cell A



We do the same for cell B. In cell B the output is the production of CheZ and RFP and the results are plotted in figure 7. We notice that CheZ-GFP is highly sensitive for the association rates of the LuxR/AHL complex and the LuxR/AHL dimer. RFP production is also sensitive for the association of the LuxR/AHL dimer.

We see that in cell B the dimerisation steps are really important. This is logical, since cell B is dependent on external AHL concentrations to start the transcriptional network. Thus, the steps concerning the binding of AHL to LuxR and making the activated LuxR/AHL dimer, which can start the transcription, are the most important steps. They determine the sensitivity of the network to AHL and as so, form the major component of the network.

Sensitivity of cell B

Figure 7: Sensitivity analysis of parameters in cell B