Modeling

Kinetic Model - The System

Figure 1 Schematic representation of our quorum sensing based regulatory system.

The iGEM Manchester-Graz expression system is designed to regulated single and multi-gene pathways for an intestine expression. For controlling a wide range of pathways it is designed in a flexible and modular manner. In addition it can not only be used to regulate pathways in probiotics but is also applicable for industrial protein expression.

The system consists of two quorum sensing (QS) systems EsaR/I and CepR/I. The EsaR/I system belongs to the plant pathogen Pantoea stewartii [1]. Contrary to common QS-systems EsaR/I uses a repressor rather than an activator. EsaR, the repressor, binds to its corresponding binding sites on the P_esaRC promoter and represses the expression of the genes under the promoter’s control [2,3]. When a certain concentration of 3-oxohexanoyl-homoserinelactone (3OC6-HSL) that is produced by the EsaI-synthase, is reached, it leads to an allosteric conformation change in EsaR’s structure that inhibits its repressor function. When positioned in the -60 region of the P_esaS-promoter EsaR can also work as an activator, by facilitating RNA-polymerase recruitment [2].

The second QS-System we are using, CepR/I, belongs to the opportunistic pathogen Burkholderia cenocepacia. Similar to the LuxR/I system, CepR acts as an activator of its corresponding promoter, P_aidA, when a certain level of octanoyl-homoserinelactone (C8-HSL) is reached [4]. C8-HSL is produced by CepI.

Our system is designed in a way that EsaR, EsaI and CepR are constitutively expressed by the same promoter. As long as the 3OC6-HSL concentration is low enough, EsaR will additionally increase its own transcription, creating a positive feedback loop (Fig. 1). When the 3OC6-HSL threshold is reached, transcription of the P_esaRC initiates, while the P_esaS-feedback loop is turned off. The activation of the promoter is shown and measured on the expression of cyan fluorescent protein (CFP). Additionally to the reporter gene also CepI gets expressed, resulting in the time-shifted activation of our second QS-system. When the C8-HSL threshold is reached, CepR can work as an activator of the P_aidA promoter that transcribes monomer red fluorescent protein (mRFP) as a second reporter gene (Fig. 1).

CFP and mRFP are used to test the system and can afterwards be replaced by genes of a desired pathway. The system is very flexible and can be optimized to get the best results depending on the introduced genes. In this model the system was optimized to regulate the Dopamine/L-DOPA pathway.

The Model

While planning our project we decided to model the system with Matlab to see what we expect to get. Modeling provides also a way to test if the system works as we planned it.

The quorum sensing system model was done using Matlab and Simbiology. The genetic setup was modeled graphically and parameters were taken from literature or estimated. A Simbiology-Model consists of species (e.g. DNA, mRNA, Proteins) that interact in reactions (e.g. DNA-> DNA + mRNA). Each reaction has parameters that affect the velocity of the reaction.

Figure 2 Diagram view of the system in Simbiology

The modeling was started on the DNA-basis. The DNA is transcribed into mRNA which gets translated into Proteins. The DNA can be bound by two proteins (Fig. 2), EsaR or CepRHSL2 (dimer of CepR bound to C8-HSL). Protein binding leads to either a basal or a maximal expression of the three different mRNAs depending if the corresponding promoter is affected by the binding (for more details look into the equations). e.g. the expression of EsaR, EsaI and CepR gets induced further by the binding of EsaR, whereas the expression of CFP and CepI gets repressed.

Our model consists of two compartments: cell (genetic setup), surrounded by the environment. The homoserine lactones are transported through the cell membrane into the environment due to diffusion. The diffusion into the environment and back are influenced by the size of the cell.

Cell growth and division

For simplification some parameters had to be assumed. Cell population was estimated as a one cell organism surrounded by the environment, so the output of the model is based on one cell in the whole system. Since the amount of all species in the model increase during the simulation a division time had to be assumed (τ=45 min). From the division time we get the rate the amount of DNA increases, which is μ = ln2/τ, the growth rate [6].

This is also valid for the DNA currently bound to some species during the simulation.

The number of cells in the model (N) had also to be defined as increasing with the fixed parameter μ, the growth rate. The number of cells defines further the total volume of cells in the model (Fig. 3). The starting number of cell (N_o), the volume of a single cell (V_cell) and the total volume of the system (V_tot) are determined in the beginning of the simulation.

Figure 3 Schematic representation of the cell growth and included parameters. The total volume is constant. The cell population grows over time and affects the volume of the environment.

Without cell division we would get aggregation of species in the model. To get no aggregation of species cell division is simulated by decreasing all species by the growth rate. Additionally all species are constantly decreasing with a specific rate (=degradation rate) to simulate natural degradation of mRNA, Proteins etc.

Diffusion through the cell wall was modeled because some species can diffuse through the cell membrane into the environment. To get the rate of the diffusion a new parameter “D” had to be introduced. D is the diffusion rate by which the species get through the cell membrane. D was found in literature to be 10 min^-1 [6]. The diffusion back into the cells is more complicated. It depends on the total cell density, which changes over time. So a new parameter “r” which depends on the relation of a single cell volume to environment volume has to be defined. The total volume which is the sum of cell volume (V_ctot) and environment volume (V_ext) is set to a constant volume over the simulation [6].

Transcription and Translation

Transcription speed was found to be 45 nt/s [8]. Since the mRNAs have all a different length, the transcription rate depends on how many nucleotides have to be transcribed per mRNA.

e.g.:
mRNAEsaR_EsaI_CepR: 2400 bp -> ~1 mRNA/min
mRNARFP:                         850 bp -> ~3 mRNA/min

Same is for translation which depends on the length of the amino acid sequence of the Proteins translated (20AA/s [7]).

e.g.:
EsaR: 250AA -> ~5 molecules/min
mRFP: 225AA -> ~6 molecules/min

Binding of Proteins and DNA

The Binding of Proteins to other Proteins or DNA in this model was measured in molecules/time related to the cell volume. Most references showed a value of e.g. 0.1 1/(nM*min). Those values were calculated back to the values before mentioned, because of Simbiology where all parameters should be kept in the same unit (molecules). Using the Avogadro constant (6.022*10²³ molecules/mole) and the cell volume the parameters were estimated to values and units to fit into the model. For details, look into the listed Parameters below (Tab 1).

Parameters

Most of the parameters were estimated after literature values of the LuxR/I system due to the lack of literature about the EsaR/I and CepR/I quorum sensing systems. For all parameters standardized prefixes were used to keep the model clear.

Table 1 Values and prefixes of parameters used in the model

All Parameters in detail:

Species

All the species used in the model are listed in table 2 with a short description. The species names were kept as short as possible but clear enough not to get confused.

Table 2 Species used in the model

Reactions

Overview :

All Reactions in detail :

DNA -> DNA + mRNAEsaR_EsaI_CepR (basal)
DNA -> DNA + mRNARFP (basal)
DNA -> DNA + mRNACFP_CepI (max)
DNA_EsaR -> DNA_EsaR + mRNAEsaR_EsaI_CepR (max)
DNA_EsaR -> DNA_EsaR + mRNACFP_CepI (basal)
DNA_EsaR -> DNA_EsaR + mRNARFP (basal)
mRNAEsaR_EsaI_CepR -> mRNAEsaR_EsaI_CepR + EsaR + EsaI + CepR
EsaR + DNA <-> DNA_EsaR
EsaI -> EsaI + 3OC6-HSL
CepI -> CepI + C8-HSL
EsaR + 3OC6-HSL <-> EsaRHSL
2*EsaRHSL <-> EsaRHSL2
mRNACFP_CepI -> mRNACFP_CepI + CFP + CepI
CepR + C8-HSL <-> CepRHSL
2*CepRHSL <-> CepRHSL2
CepRHSL2 + DNA <-> DNA_ CepRHSL2
DNA_CepRHSL2 -> DNA_CepRHSL2 + mRNARFP (max)
DNA_CepRHSL2 -> DNA_CepRHSL2 + mRNAEsaR_EsaI_CepR (basal)
DNA_CepRHSL2 -> DNA_CepRHSL2 + mRNACFP_CepI (max)
mRNARFP -> mRNARFP + mRFP
8C-HSL <-> 8C-HSLext
3OC6-HSL <-> 3OC6-HSLext

EsaI -> Ø
CepI -> Ø
EsaR -> Ø
CepR -> Ø
mRNAEsaR_EsaI_CepR-> Ø
mRNACFP_CepI-> Ø
mRNARFP-> Ø
CFP-> Ø
mRFP-> Ø
8C-HSL -> Ø
3OC6-HSL -> Ø
8C-HSL_ext -> Ø
3OC6-HSL_ext -> Ø
EsaRHSL-> Ø
EsaRHSL2-> Ø
CepRHSL-> Ø
CepRHSL2-> Ø
DNA ->DNA + DNA (DNA duplication)
DNA-> Ø (DNA loss due to cell division)
DNA_CepRHSL2 -> DNA_CepRHSL2 + DNA
DNA_EsaR -> DNA_EsaR + DNA

Equations

d([3OC6_HSL])/dt = f3OC6HSL*EsaI – (dHSL1*[3OC6_HSL]+ μ*[3OC6_HSL]) – (D*[3OC6_HSL]) + r*D*[3OC6_HSLext] – (bEsaRHSL1*EsaR*[3OC6_HSL]-ubEsaRHSL1*EsaRHSL)

d(C8_HSL)/dt = f8CHSL*CepI – (dHSL2*C8_HSL+ μ*C8_HSL) – (D*C8_HSL) + r*D*C8_HSLext – (bCepR*CepR*C8_HSL-ubCepR*CepRHSL)

d(CepI)/dt = - (dCepI*CepI+ μ*CepI) + tlCFP*mRNACFP_CepI

d(CepR)/dt = - (dCepR*CepR+ μ*CepR) + tlRNA1*mRNAEsaI_EsaR_CepR – (bCepR*CepR*C8_HSL-ubCepR*CepRHSL)

d(EsaR)/dt = - (dEsaR*EsaR+ μ*EsaR) + tlRNA1*mRNAEsaI_EsaR_CepR – (bEsaRHSL1*EsaR*[3OC6_HSL]-ubEsaRHSL1*EsaRHSL) – (bDNA_EsaR*EsaR*DNA-ubDNA_EsaR*DNA_EsaR)

d(EsaI)/dt = - (dEsaI*EsaI+ μ*EsaI) + tlRNA1*mRNAEsaI_EsaR_CepR

d(EsaRHSL2)/dt = bEsaRHSL*EsaRHSL*EsaRHSL-ubEsaRHSL*EsaRHSL2 – (dEsaRHSL2*EsaRHSL2+ μ*EsaRHSL2)

d(CepRHSL2)/dt = bCepRHSL*CepRHSL*CepRHSL-ubCepRHSL*CepRHSL2 – (dCepRHSL*CepRHSL2+ μ*CepRHSL2) – (bCepRHSL2*CepRHSL2*DNA-ubCepRHSL2*DNA_CepRHSL2)

d(DNA)/dt = μ*DNA – (bCepRHSL2*CepRHSL2*DNA-ubCepRHSL2*DNA_CepRHSL2) – (bDNA_EsaR*EsaR*DNA-ubDNA_EsaR*DNA_EsaR) – (μ *DNA) + μ *DNA_EsaR + μ*DNA_CepRHSL2

d(CFP)/dt = tlCFP*mRNACFP_CepI – (dCFP*CFP+ μ*CFP)

d(mRNARFP)/dt = -( dmRNAmRFP*mRNARFP+ μ*mRNARFP) + tRFP*DNA_CepRHSL2 + tRFP2*DNA + tRFP2*DNA_EsaR

d(mRFP)/dt = tlRFP*mRNARFP – (dRFP*mRFP+ μ*mRFP)

d(EsaRHSL)/dt = -( -2* bEsaRHSL*EsaRHSL*EsaRHSL-ubEsaRHSL*EsaRHSL2) + bEsaRHSL1*EsaR*[3OC6_HSL]-ubEsaRHSL1*EsaRHSL – (dEsaRHSL*EsaRHSL+ μ*EsaRHSL)

d(CepRHSL)/dt = -( -2* bCepRHSL*CepRHSL*CepRHSL-ubCepRHSL*CepRHSL2) + bCepR*CepR*C8_HSL-ubCepR*CepRHSL – (dCepRHSL*CepRHSL+ μ*CepRHSL)

d(DNA_CepRHSL2)/dt = bCepRHSL2*CepRHSL2*DNA-ubCepRHSL2*DNA_CepRHSL2

d(DNA_EsaR)/dt = bDNA_EsaR*EsaR*DNA-ubDNA_EsaR*DNA_EsaR

d(mRNAEsaI_EsaR_CepR)/dt = -( dmRNAesaI*mRNAEsaI_EsaR_CepR + μ *mRNAEsaI_EsaR_CepR) + tDNA1*DNA + tDNA1_2*DNA_EsaR + tDNA1*DNA_CepRHSL2

d(mRNACFP_CepI)/dt = -( dmRNA2*mRNACFP_CepI+ μ*mRNACFP_CepI) + tDNA2*DNA + tDNA2_2*DNA_EsaR + tDNA2*DNA_CepRHSL2

d([3OC6_HSLext])/dt =D*[3OC6_HSL] – (r*D*[3OC6_HSLext]) – (d3OC6_HSL*[3OC6_HSLext])

d(C8_HSLext)/dt = D*C8_HSL – (r*D*C8_HSLext) – (dC8_HSL*C8_HSLext)

Results

Overview

Figure 4 Quantity of the different free, unbound species in the cell over time. CepR gets bound by C8-HSL immediately and is therefore not found in free form.

The expression of EsaR, EsaI and CepR are driven by the same transcription rate, as they are transcribed onto the same mRNA. Whereas EsaR gets bound to 3OC6-HSL and dimerised to EsaRHSL2, EsaI increases during time and produces 3OC6-HSL (Fig. 4). CepR is not really there at any time, because it binds to C8-HSL fast, which is produced by CepI. This is due to the transcription of mRNACFP_CepI before EsaR gets expressed and can repress transcription of mRNACFP_CepI. The first gene of interest (CFP) gets expressed just as CepI from the beginning. The formation of CepRHSL2 from CepR and C8-HSL induces further the transcription of the second gene of interest (mRFP).

DNA binding

Figure 5 Binding of the DNA to different Proteins over time.

If we simulate DNA binding over time we get a graphical trend that shows DNA gets bound to different proteins at different time points (Fig. 5). While free DNA (blue) decreases the other species increase. After about 2 minutes DNA gets first bound to EsaR (yellow), which is expressed constitutively and its expression is positively influenced by a feedback loop by binding in front of P_esaS which facilitates polymerase recruitment. However, EsaR binds to the corresponding homoserine lactone (HSL) 3OC6-HSL (Fig.4, red), which is produced by the synthase EsaI. The transcription of the mRNAEsaR_EsaI_CepR decreases after binding of EsaR to 3OC6-HSL (Fig. 6), turning off the positive feedback loop. Simultaneously repression of P_esaRC is stopped. Thus free DNA, bound to no Protein complex leads to a fast transcription of CFP and CepI, whereas the binding to EsaR represses its transcription. The second gene of interest (mRFP) gets expressed first with a small basal expression but after the DNA binds to CepRHSL2 the transcription gets induced (Fig. 6).

Figure 6 mRNA transcription over time.

Expression

Figure 7 Expression of the genes of interests over time.

Expression of the gene of interest is also observed over time (Fig. 7, 8). We can see the expression of the second gene of interest (mRFP) starts ~12min later than the expression of the first gene of interest (CFP). At first the expression is shifted ~30min to reach the same amount of molecules, but after ~2.5 hours the expression of mRFP is increasing and reaching the same level as CFP (Fig. 8). After 3 hours mRFP is expressed more than CFP and is increasing strongly. This is due to the stronger transcription rate of mRNARFP, which is way shorter than mRNACFP_CepI.

The model showed that it is capable to handle the expression of various pathways and is flexible if the length of mRNAs or Proteins differ. It is based on molecules because of the easier handling of Simbiology. However, the amount of molecules can also be recalculated to a concentration value. For CFP we get for example a value of 1.25 μmol/(min*L) per cell.

Figure 8Long term observation of the expression of the two genes of interest.

Incorporation of the Dopamine/L-DOPA pathway

The Dopamine and L-DOPA pathway is a good example to test the expression system because it consists of three enzymes: Tyrosine hydroxylase, CYP2D6 (Cytochrome P450 enzyme 2D6) and the AADC (Aromatic L-amino acid decarboxylase). With these enzymes tyrosine is converted intracellular into dopamine (Fig. 9).

Figure 9 L-Dopa and Dopamine pathway. Tyrosinase (=Tyrosine Hydroxylase), AADC and CYP2D6 are the enzymes used to incorporate the pathway to the system.

The incorporation of the pathway was done by changing the mRNA transcription rate and the translation rate of the different proteins. This was necessary, because the length of the DNA sequences and the proteins change the rate of the transcription and translation.

New species had to be added to the model to introduce the pathway. Dopamine and L-Dopa are the desired end products of the pathway and can also diffuse through the cell membrane (Fig. 10). The tyrosine amount in the cell was estimated as 280 molecules, which is 1% of the total tyrosine in the cytoplasm [13]. The available tyrosine in the cell is held by dosing tyrosine constantly to the cell, since the cell produces constantly tyrosine.

Figure 10 Overview of incorporating the pathway to the model. New species and reactions were added to model the pathway
which is controlled by the previous expression system.

To get an overview of the model all new parameters and species are listed in the following tables.

Table 3 New Parameters to introduce the Pathway

Table 4 New Species for the pathway

Enzyme velocity

AADC is known to be the rate limiting step in the pathway. The enzyme velocity of the Tyrosinase was found in literature to be 0.3 conversions 1/min [14]. So the velocities of the other enzymes were estimated according to this value.

Results

Figure 11 Expression of the three genes of the pathway. The graph of the Tyrosinase is the same as CYP2D6 (yellow).

The three enzymes of the pathway were modeled in a way all enzymes are expressed in an equally amount. This is the case if the first gene of interest in the expression system (CFP) is replaced by the Tyrosinase and CYP2D6, whereas the second gene of interest (mRFP) gets replaced by AADC (Fig. 11).

Figure 12 Dopamine and L-Dopa yield in the environment/gut.

The model showed a high yield of dopamine and a lower yield of L-DOPA in the gut (Fig. 12). If we recalculate the molecule values into concentrations we get a yield of 401.75 μmol/h/L or 61.54 mg/h/L. This means, it takes 1.6 hours to get an amount of 100mg Dopamine in a 1L approach. By changing the whole setup and not including the second quorum sensing system, we can regulate one enzyme and not a whole pathway. If we change the system in a way L-DOPA is the only product we get the following setup (Fig. 13).

Figure 13 Simbiology overview of the changed genetic setup

The simulation gives us a high L-DOPA yield of 415.14 µmol/h/L or 81.86 mg/h/L (Fig. 14).

Figure 14 L-DOPA production and yield in the environment.

The 81.86 mg/h/L we calculated back from the molecules we get from the simulation are a promising result. Current Parkinson’s disease treatment involves either relatively high doses of L-DOPA by taking pills at a time or a constant flux of 20-200mg L-DOPA per hour using the Duodopa ® therapy [15]. Since our project is designed to replace currently used treatments it is really promising that one day our DopaDoser is being used in modern medicine as a treatment.

The modeling of the Dopamine and L-DOPA pathway showed a wide range of possible application for the expression system. It is feasible to control single and multi-gene pathways and by its modular design it is easy to get exactly the results needed.

Flux Analysis