Difference between revisions of "Team:Manchester-Graz/Modeling"

Line 318: Line 318:
 
<p>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.</p>
 
<p>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.</p>
  
<p>In collaboration with UCL iGEM and Norwich iGEM we also introduced other single and multi gene pathways into the model and observed the production of the corresponding products</p>
+
<p>In collaboration with <a href="https://2015.igem.org/Team:Manchester-Graz/Collaborations/UCL" target="_blank">UCL iGEM</a> and <a href="https://2015.igem.org/Team:Manchester-Graz/Collaborations/NRP-UEA-Norwich" target="_blank">Norwich iGEM</a> we also introduced other single and multi gene pathways into the model and observed the production of the corresponding products.</p>
  
 
<div style="height:50px; width:900px;"></div>
 
<div style="height:50px; width:900px;"></div>

Revision as of 18:46, 18 September 2015

iGEM Manchester Header

iGEM Manchester - Modeling

Modeling

Key Achievements

  • Full genetic setup successfully rebuild, including 40 parameters and over 20 equations, using Matlab and Simbiology
  • L-DOPA pathway introduced into the model and expression in the gut observed
  • Simulated gut energy/oxygen conditions in existing E. coli lab strain model (BL21).
  • L-DOPA biosynthesis achievable from different diets, provided growth rate of probiotics is regulated or constrained.

Kinetic Model - The System

System_new
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 PesaRC 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 PesaS-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, PaidA, 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 PesaRC initiates, while the PesaS-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 PaidA 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.

Manchester-Graz_Modeling_Fig2
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].

Manchester-Graz_Modeling_Reaction1

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

Manchester-Graz_Modeling_Reaction2

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 (No), the volume of a single cell (Vcell) and the total volume of the system (Vtot) are determined in the beginning of the simulation.

Manchester-Graz_Modeling_Formulas1
Manchester-Graz_Modeling_Fig3
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 (Vctot) and environment volume (Vext) is set to a constant volume over the simulation [6].

Manchester-Graz_Modeling_Formulas2

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*1023 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
Manchester-Graz_Modeling_Tab1

All Parameters in detail:

Manchester-Graz_Modeling_Parameters

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
Manchester-Graz_Modeling_Tab2

Reactions

Overview :

Manchester-Graz_Modeling_Formulas3

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

Manchester-Graz_Modeling_Fig4
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

Manchester-Graz_Modeling_Fig5
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 PesaS 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 PesaRC 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).

Manchester-Graz_Modeling_Fig6
Figure 6 mRNA transcription over time.

Expression

Manchester-Graz_Modeling_Fig7
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.

Manchester-Graz_Modeling_Fig8
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).

Manchester-Graz_Modeling_Fig9
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.

Manchester-Graz_Modeling_Fig10
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
Manchester-Graz_Modeling_Tab3
Table 4 New Species for the pathway
Manchester-Graz_Modeling_Tab4

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.

Manchester-Graz_Modeling_Tab5

Results

Manchester-Graz_Modeling_Fig11
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).

Manchester-Graz_Modeling_Fig12
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).

Manchester-Graz_Modeling_Fig13
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).

Manchester-Graz_Modeling_Fig14
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.

In collaboration with UCL iGEM and Norwich iGEM we also introduced other single and multi gene pathways into the model and observed the production of the corresponding products.

[1] http://www.ipm.iastate.edu/ipm/info/plant-diseases/stewarts-wilt
[2] Shong et al (2013) Engineering the esaR Promoter for Tunable Quorum Sensing- Dependent Gene Expression
[3] Shong et al (2013) Directed Evolution of the Quorum-Sensing Regulator EsaR for Increased Signal Sensitivity
[4] Weingart et al (2005) Direct binding of the quorum sensing regulator CepR of Burkholderia cenocepacia to two target promoters in vitro
[5] https://2014.igem.org/Team:ETH_Zurich/modeling/parameters
[6] Weber M., Buceta J. (2013) Dynamics of the quorum sensing switch: stochastic and non-stationary effects, BMC Systems Biology
[7] Bremer, H., Dennis, P. P. (1996) Modulation of chemical composition and other parameters of the cell by growth rate
[8] Vogel U, Jensen KF (1994) The RNA chain elongation rate in Escherichia coli depends on the growth rate
[9] Goryachev, A.B., D.J. Toh and T. Lee (2006) System analysis of a quorum sensing network: Design constraints imposed by the functional requirements, network topology and kinetic constant
[10] Zhu J, Winans SC (2001) The quorum-sensing transcriptional regulator TraR requires its cognate signaling ligand for protein folding, protease resistance, and dimerization
[11] Prescott LM, Harley JP, Klein DA (1996) Microbiology
[12] Roberts C et al (2006) Characterizing the effect of the Staphylococcus aureus virulence factor regulator, SarA, on log-phasemRNA half-lives
[13] Bennett BD et al (2009) Absolute metabolite concentrations and implied enzyme active site occupancy in Escherichia coli
[14] http://adamcap.com/schoolwork/kinetics-analysis-of-tyrosinase/
[15] Wenzel, K. et al. (2014) P-aktuell: Aktuelles zur Pumpentherapie mit intrajejunalem Levodopa-Gel

Flux Balance Analysis

We aimed to use modelling to determine numerous things to help shape our product design, such as optimise compound productions under given conditions and obtained an ideal energy source that would optimise our compound of interest. To achieve this, we used COBRA Toolbox that implemented in MATLAB to conduct Flux Balance analysis on various gut strains of E.coli. This involved sourcing existing models of our particular strain BL21, was found on BiGG database [1]. This would serve as our primary strain as the same particular strain was used for transformation and protein expression in the experimental outcome.

The concept was to optimize the non-native reaction by adding constrains and prioritising a particular reaction, such that a theoretical value is presented that defines the possibility of the reaction to take place [2]. In the models, there are numerous reactions that make up a large metabolic network. Addition of non-native reaction represents the reaction that occurs in the bacteria as a result of manipulating the genomic structure of the bacteria. This expression would then lead to synthesis of the compound of interest in the bacteria. Thus, the objective of modelling is to provide additional support to the lab experiments and we aim to quantify levels of L-DOPA and dopamine via modelling.

All metabolites involved in the large metabolic are produced and consumed, such that a steady state model (the net flux of each metabolite is zero) is maintained [3]. The flux values are obtained by a linear program solver established in MATLAB that provides possible solution such that a steady state model is maintained [3]. The media components were all added to the model and were constrained such that it resembles the uptake rate of lab strain.

To calculate the growth rates of E.coli BL21 the reaction that was made the objective function was 'Ec_biomass_iJO1366_core_53p95M’. This reaction in the model defines the survival of the bacteria and the flux through this reaction shows the growth rate of the bacteria on different energy source which are typically found in diet. This reaction was then further constrained such that flux values are obtained, that represents the synthesis of L-DOPA and Dopamine by the bacterial strain while the survival reaction also occurs in the model [4]. This dual objective function also known as the Flux Variability Analysis and it was essential in order to give a more realistic synthesis value or production rate. All the data obtained were tabulated and various Figures were produced for further analysis.

The non-native reactions added were:

For L-DOPA Synthesis:

Manchester-Graz_Model_Man_Reaction1

For Dopamine Synthesis (two step reaction):

Manchester-Graz_Model_Man_Reaction2

The media components added were:

Acetate, Adenine, L-alanine, L-arginine, L-asparagine, L-aspartate, Biotin Calcium, Cobalamin, Chloride, Cobalt, Copper, L-cysteine, Iron (II), Iron (III), Folate, L-glutamine, L-glutamate, Glycine, Sulphide, L-histidine, L-isoleucine, Potassium, L-cysteine, L-leucine, L-lysine, L-methionine, Magnesium, Sodium, Niacin, Ammonia, L-phenylalanine, Phosphate, Panthothenate, L-proline, Propionate, Pyridoxal, Riboflavin, L-serine, Sulphate, Thiamine, L-threonine, L-tryptophan, L-tyrosine, Uracil, L-valine, Xanthine [7]

All the components were constrained to an uptake rate of 1 mmol gram of dry weight per hour (mmolgDW -1h-1). Carbon/Energy source and Oxygen levels were constrained to an uptake rate of 10 mmolgDW -1h-1 [7].

Results and Discussion

Manchester-Graz_Model_Man_Fig%201
Figure 15 Maximum E.coli BL21 growth rates on different Energy Sources.

Figure 15 shows the growth rates of E.coli BL21 on different energy/carbon Sources. From the Figure it can be seen that using Glucose and Ascorbate substrate, the bacteria has high flux values that represent the growth rate (hr-1). These varying growth rates are obtained because each substrate is metabolized specifically, giving varying reactant starting values for the survival reaction to use. These results indicate that the bacteria could grow on various diet sources and thus the diet could be designed accordingly for adequate growth of bacteria to constantly produce the compound of interest (L-DOPA and Dopamine).

Manchester-Graz_Model_Man_Fig%202
Figure 16 L-DOPA (Blue) and Dopamine (Red) production at maximum growth rate of E.coli BL21 at different Energy Sources. The values shown are the synthesis rate of compounds in mmol gram of dry cell weight per hour (mmolgDW -1h-1).

Figure 16 shows the results of L-DOPA synthesized by the bacteria over different carbon sources. These results are dependent on Figure 15 as the survival reaction was then constrained to the maximum growth rate that yields out of particular carbon source. From this Figure, it could be seen that, apart from glucose used as substrate, there was no flux value obtained and hence no L-DOPA was produced on other substrates. From these results, it seems that apart from glucose; no other energy sources produced any significant levels of L-DOPA and dopamine. This is because the necessary metabolites (L-Tyrosine which is also important substrate to produce L-DOPA and dopamine) were used up for the survival reaction. As this reaction was constrained to the flux that resembles the maximum growth rate, so more weight is given to the survival reaction and thus, no L-DOPA or Dopamine was produced in the BL21 model. Since glucose as energy source produced the highest growth rate, there were enough metabolites in the model required to L-produce L-DOPA and Dopamine. This was not the case for other sources and hence, there was no _DOPA or dopamine produced.

From the same Figure, it could be seen that same levels (same flux values of L-DOPA and Dopamine) was produced even though the enzymes that catalyses both the reaction are different and thus resulting outcome should differ. The main constrain added in the non–native reaction is uptake rate of oxygen and energy source and due to the similarity of the added reactions (in both cases the levels of oxygen is essential) the model fails to differentiate between the two reaction and thus, same levels of compounds are predicted to be produced by the bacterial model.

Manchester-Graz_Model_Man_Fig%203
Figure 17 L-DOPA (Red) and Dopamine (Blue) production at half the maximum growth rate of E.coli BL21 using different Energy Sources. The values shown are the synthesis rate of compounds in mmol gram of dry cell weight per hour (mmolgDW -1h-1).
Manchester-Graz_Model_Man_Fig%204
Figure 18 L-DOPA (Blue) production as a function of increasing growth rate of E.coli BL21 using glucose as the energy source. Similar trend is also seen on when the L-DOPA synthesis rate were obtained as a function of increasing bacterial growth rate using different Energy Sources.

Figure 17 shows the flux values where the survival reactions are constrained to half their growth rate. Unlike Figure 16 it could be seen that the compounds are produced using various substrates as energy source. This is because there are essential metabolites available for the non-native reactions to take place since not all of the metabolites would be used for the survival reaction. This is further shown by the Figure 18 where the bacterial growth potential was constrained to increasing level of percentage growth. The Figure shows a negative trend, as the percentage of growth rate increases, the flux value for L-DOPA production decreases.

The growth of the bacteria could be constrained by manipulating the genome of the bacteria, not only to induce the necessary enzymes and cofactors within them, but also to constrain their growth to maintain a certain density/population. The necessary compounds (L-DOPA) could be produced as a function of bacterial growth rate/ constant cell density and thus could be altered as per the dietary and dosage requirements of the Parkinson’s patients. Thus, this holds great potential in terms of being a novel therapeutic delivery system for the treatment of Parkinson’s.

Limitations and Future Outlook

There was no difference in the levels (flux values) of L-DOPA and Dopamine produced despite the fact that both are synthesized independent of each other and the enzymes responsible for their catalysis also differ. This is because the enzyme kinetics and catalysis are not part of this particular analysis and there is only one major constrain added to the non-native compound catalysis reaction is the uptake of available oxygen. This model could be further improved by incorporating constrain-based analysis with enzyme kinetics. This would then give levels and rate of L-DOPA and Dopamine that would have more realistic relavence.

Furthermore, the model is largely limited to any activities occurring in the extracellular and intracellular regions of the bacteria and fails to capture any external factors outside the model, which could have an adverse effect on the synthesis of the compound. Since the existing model used for analysis just takes into account the effect of what is already present in the model, it fails to capture any components that might have a potential effect on our pathway. For example, there are other energy sources that are absent to the model and therefore any effects of those components couldn’t be captured.

The physical values obtained from this model lack any physiological relevance and thus more substantial modelling is required in order to attain data with significant in vivo relevance. The model used is largely static and thus more complex co-model (combination of human and bacteria model) is required and dynamic flux balance analysis [8] is required in order to achieve results with in vivo cues. Furthermore, this model could be further improvised by knocking out specific reaction to further optimise compound synthesis. Thus could be achieved by running a pre-existing opt knock script that can determine mutant probiotics strain. This would knockout any unnecessary reaction that could have potential side effects when this therapeutics system is applied in clinical practice [9].

[1] http://bigg.ucsd.edu/models/iB21_1397
[2] Y. Z. a. J. H. Yue-Dong Gao, "Metabolic Modeling of Common Escherichia coli Strains in Human Gut Microbiome," Hindawi Publishing Corporation, pp. 1-11, 2014.
[3] http://www.ncbi.nlm.nih.gov/pmc/articles/PMC3108565/
[4] http://cobramethods.wikidot.com/flux-variability-analysis
[5] http://www.kegg.jp/dbget-bin/www_bget?rn:R01815
[6] http://www.genome.jp/dbget-bin/www_bget?R02383
[7] http://www.pnas.org/content/104/7/2307/F1.expansion.html - rates
[8] Gonçalo dos Santos Correia (2012) COUPLING METABOLIC FOOTPRINTING AND FLUX BALANCE ANALYSIS TO PREDICT HOW SINGLE GENE KNOCKOUTS PERTURB MICROBIAL METABOLISM.
[9] http://www.ncbi.nlm.nih.gov/pubmed/14595777