Team:KU Leuven/Modeling/Internal
Internal Model
1. 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
2. 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.
3. 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 chose 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:
- The effects of cell division can be neglected
- The substrate pool can not be depleted and the concentration (or amount of molecules) of substrate in the cell is constant
- The exterior of the cell contains no leucine at t=0 and is perfectly mixed
- Diffusion happens independent of cell movement and has a constant rate
4. System
After this extensive literature search, we can finally set up our complete system of ODEs for every cell.
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:
Figure 1: The internal network of cell A visualized by Simbiology
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} $$
We visualize these ODE's in the Simbiology toolbox. This gives us the following diagrams:
Figure 2: The internal network of cell A visualized by Simbiology
5. 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.
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.
Figure 4: Simulation of all processes in Cell B in OFF state with and without AHL induction
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.
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.
