Team:Manchester-Graz/Modeling
Modeling
Key Achievements
- Full genetic setup of DopaDoser control circuit successfully rebuilt, including 40 parameters and over 20 equations, using Matlab and Simbiology
- L-DOPA pathway introduced into the model and expression in the gut predicted
- 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
Figure 1 Schematic representation of our quorum sensing-based regulatory system.
The iGEM Manchester-Graz expression system is designed to regulate single and multi-gene pathways for intestinal expression. The expression system is designed in a flexible and modular manner to control a wide range of pathways. In addition, it can be used not only to regulate pathways in probiotics, but also for industrial protein expression.
The system consists of two quorum sensing (QS) systems: EsaR/I and CepR/I. The EsaR/I system is derived from the plant pathogen Pantoea stewartii [1]. In contrast 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), produced by the EsaI-synthase, is reached, it leads to an allosteric conformational change in EsaR’s structure and 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), produced by CepI, is reached [4].
Our system is designed in such a way that EsaR, EsaI and CepR are constitutively expressed under control of 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 using the expression of cyan fluorescent protein (CFP). In addition to the reporter gene, CepI also 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 a monomer of 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 explore the range of conditions in which our system would be functional. Modeling also provided a way to predict if the system actually would work as we planned it.
The quorum sensing system model was done using Matlab and Simbiology. The genetic set-up was modeled graphically and parameters were taken from the 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 at the DNA level. 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 see 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: the cell (genetic set-up) and the surrounding 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 [5].
Cell growth and division
For simplification some parameters had to be estimated. The cell population was described as a single-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 other 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 cells (No), the volume of a single cell (Vcell), and the total volume of the system (Vtot) 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 see continuous accumulation of molecules in the model. To avoid this, cell division is simulated by decreasing all species by the growth rate (dilution). 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 as well, as some species can diffuse through the cell membrane into the environment. To describe the rate of 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 estimated from the literature to be 10 min-1 [6]. The diffusion back into the cells is more complicated to model. 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].
Transcription and Translation
Transcription speed was estimated to be 45 nt/s [7]. Since the mRNAs have all a different lengths, 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 [8]).
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 converted back to the values mentioned before, because in Simbiology all parameters should be kept in the same unit (molecules). For details, see the list of parameters below (Tab 1).
Parameters
Most of the parameters were estimated based on 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.
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.
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 is driven by the same transcription rate, as they are transcribed into 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 detectable in its free form 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 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, not bound to protein, 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 level of 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 our planned system is capable to handle the expression of various pathways and is flexible if the length of mRNAs or proteins differ.
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 in the model 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 constant in the modeled 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.
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/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 set-up and not including the second quorum sensing system, we can regulate one enzyme and not a whole pathway. If we change the system such that L-DOPA is the only product, we get the following circuit (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 determined 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.
[2] Shong, J. and Collins, C. (2013) Engineering the esaR promoter for tunable quorum sensing-dependent gene expression. American Chemical Society. 2 (10), pp. 568–575.
[3] Shong, J., Huang, Y-M., Bystroff, C. and Collins, C. (2013) Directed evolution of the quorum-sensing regulator EsaR for increased signal sensitivity. American Chemical Society. 8 (4), pp 789–795.
[4] Weingart, C., White, C., Liu, S., Chai, Y., Cho, H., Tsai, C., Wei, Y., Delay, N., Gronquist, M., Eberhard, A. and Winans, S. (2005) Direct binding of the quorum sensing regulator CepR of Burkholderia cenocepacia to two target promoters in vitro. Molecular Microbiology. 7 (2), pp 452-467.
[5] Mosaicoli (2014) Parameters and tools. [online] Available at: https://2014.igem.org/Team:ETH_Zurich/modeling/parameters [Accessed 28 Aug. 2015]
[6] Weber, M. and Buceta, J. (2013) Dynamics of the quorum sensing switch: stochastic and non-stationary effects. BMC Systems Biology, 7(1), pp. 6.
[7] Vogel, U. and Jensen, K. (1994) The RNA chain elongation rate in Escherichia coli depends on the growth rate. Journal of Bacteriology, 176(10), pp.2807-13.
[8] Bremer, H., Dennis, P. P. (1996) Modulation of chemical composition and other parameters of the cell by growth rate
[9] Goryachev, A., Toh, D. and Lee, T. (2006) Systems analysis of a quorum sensing network: Design constraints imposed by the functional requirements, network topology and kinetic constants, Biosystems, 83(2-3), pp. 178–187.
[10] Zhu, J. and Winans, S. (2001) The quorum-sensing transcriptional regulator TraR requires its cognate signaling ligand for protein folding, protease resistance, and dimerization. Proceedings of the National Academy of Sciences, 98(4), pp. 1507–1512.
[11] [11] Prescott, L., Harley, J. and Klein, D. (1996) Microbiology, 3rd edition, WM. C. Brown Company Publishers, Iowa, US.
[12] Roberts, C., Anderson, K., Murphy, E., Projan, S., Mounts, W., Hurlburt, B., Smeltzer, M., Overbeek, R., Disz, T. and Dunman, P. (2006) Characterizing the Effect of the Staphylococcus aureus Virulence Factor Regulator, SarA, on Log-Phase mRNA Half-Lives. Journal of Bacteriology, 188(7), pp. 2593–2603.
[13] Bennett, B., Kimball, E., Gao, M., Osterhout, R., Van Dien, S. and Rabinowitz, J. (2009) Absolute metabolite concentrations and implied enzyme active site occupancy in Escherichia coli. Nature Chemical Biology, 5(8), pp. 593–599.
[14] Capriola, A. (2009) Kinetics Analysis of Tyrosinase. [online] Available at: http://adamcap.com/schoolwork/kinetics-analysis-of-tyrosinase/ [Accessed 19 Sept. 2015]
[15] Wenzel, K. et al. (2014) P-aktuell: Aktuelles zur Pumpentherapie mit intrajejunalem Levodopa-Gel
Flux Balance Analysis
We aimed to use modeling to shape our product design, for example as a tool to optimise compound production under given conditions and to determine an ideal energy source that would optimise our compound of interest. To achieve this, we used the COBRA Toolbox that is 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, which was found on BiGG database [1]. This would serve as our primary model 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 constraints 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. Additional non-native reactions represent the reactions that occur in the bacteria as a result of our genetic engineering. The objective of modelling is to predict the levels of L-DOPA and dopamine that can be produced in various growth conditions.
All metabolites involved in the large metabolic are produced and consumed, such that a steady state condition (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 solutions while optimizing a specified objective function [3]. The media components were all added to the model and were constrained to values that resemble the uptake rate of our 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 use of a dual objective function, also known as Flux Variability Analysis, was essential in order to give a more realistic value for the estimated production rate.
The non-native reactions added were:
For L-DOPA Synthesis:
For Dopamine Synthesis (two step reaction):
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
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 or Ascorbate as substrate, the bacteria have 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 robustly grow on various dietary energy/carbon sources in the human gut and thus the diet could be designed accordingly for adequate growth of bacteria to constantly produce the compound of interest (L-DOPA and Dopamine).
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 on 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 when glucose is used as substrate, there was no positive flux value obtained and hence no L-DOPA was produced on other substrates. These results indicate that apart from glucose no other energy sources are suitable for production of 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, 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) are predicted to be produced even though the enzymes that catalyses both the reaction are different and thus the outcome could be expected to differ. The main constraint added in the non–native reactions is the 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, the same levels of compounds are predicted to be produced by the bacterial model.
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).
Figure 18 L-DOPA (Blue) production as a function of increasing growth rate of E.coli BL21 using glucose as the energy source. A similar trend is also seen when the L-DOPA synthesis rate is 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 can 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 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 (flux balance analysis), and there is only one major constrain added to the non-native compound catalysis reaction, namely the uptake of available oxygen. This model could be further improved by combining constraint-based analysis with enzyme kinetics. This would then give levels and rates of L-DOPA and Dopamine that would have more realistic relevance.
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 could not be captured.
The stoichiometric model used is largely static and more complex multi-species models (combination of human and bacteria model) would be required, using dynamic flux balance analysis [8] in order to achieve results that are predictive for the in vivo behaviour of our system. Furthermore, this model could be further improvised by knocking out specific reactions to further optimise compound synthesis. This could be achieved by running a pre-existing OptKnock script that can determine optimal mutant probiotics strain [9].
[2] Gao, Y.-D., Zhao, Y. and Huang, J. (2014) Metabolic Modeling of Common Escherichia coli Strains in Human Gut Microbiome. BioMedical Research International, pp. 1–11.
[3] Orth, J., Thiele, I. and Palsson, B. (2010) What is flux balance analysis? Nature Biotechnology, 28(3), 245-248.
[4] Flux variability analysis - Constraint-based analysis (2013) [online] Available at: http://cobramethods.wikidot.com/flux-variability-analysis [Accessed: 19 September 2015]
[5] KEGG REACTION: R01815 (2015) [online] Available at: http://www.kegg.jp/dbget-bin/www_bget?rn:R01815 [Accessed: 19 September 2015].
[6] KEGG REACTION: R02383 (2015) [online] Available at: http://www.genome.jp/dbget-bin/www_bget?R02383 [Accessed: 19 September 2015].
[7] Harrison, R., Papp, B., Pal, C., Oliver, S. G. and Delneri, D. (2007) Plasticity of genetic interactions in metabolic networks of yeast. Proceedings of the National Academy of Sciences, 104(7), pp. 2307–2312.
[8] Gonçalo dos Santos Correia (2012) Coupling Metabolic Footprinting and Flux balance Analysis to predict how single gene knockouts perturb microbial metabolism.
[9] Burgard, A., Pharkya, P. and Maranas, C. (2003) Optknock: a bilevel programming framework for identifying gene knockout strategies for microbial strain optimization. Biotechnology and bioengineering, 6(84), pp. 647–657.
