Team:Virginia/Model

TRANSFOAM

Team Attributions
Description Design Experiments Model Safety
Results Future Directions Demonstrate
Parts Collection
Integrated Human Practices Public Engagement Collaborations Mid-Atlantic Meetup
Awards Medals

Introduction

Over the past decade, mathematical and computational modeling have become a key component in the field of synthetic biology. As an engineering discipline, synthetic biology necessitates a significant amount of time and resources through multiple iterations of design and testing. Repetitive simulations and models quicken the engineering design process by predicting the outcome in silico rather than in vivo.

The goal of Virginia iGEM 2019’s model is to identify, analyze, and improve upon the maximum theoretical yield of PHB production by exploring various reaction phenotypes of our bacterial device. The model was additionally explored for potential gene knockout strategies for the purpose of coupling the growth of the device to the production of PHB. To determine optimal PHB producing reaction phenotypes, the device’s metabolic network was analyzed using constraint based Flux Balance Analysis (FBA) and robustness analysis. The theoretical maximal PHB production was found to be 1.50 mol PHB per mol styrene. Trials of robustness analysis showed the device to be limited by oxygen, and also showed different metabolite sensitivities of the device with respect to styrene, oxygen, sulfate, phosphate, and magnesium. Finally, gene knockout strategies were explored using the function OptKnock. The model was unable to produce a knockout set to induce growth coupling suggesting the need for a fed-batch, inducible bioprocess to produce optimal yield of PHBs.

Overall, our constructed genome-wide metabolic model was able to inform and further develop our project by:

Simulating the behavior of our device to determine its theoretical maximum PHB yield.

Simulating PHB production’s sensitivity to several important metabolites in growth media.

Identifying the limitations of device modification which informs our bioprocess' design and implementation.




Genome-Scale Metabolic Model

With over 14 major reactions that make up the device’s biochemical pathway, there are many metabolites which enter and exit the system while interacting with other non-essential reactions in the bacteria.1 Therefore, we needed to consider the entire network of reactions and metabolites in E. coli when determining the maximum PHB yield.

A genome-scale metabolic network (GEM) was obtained from the Biochemical, Genetic, and Genomic (BiGG) database.1 The iML1515 GEM is the latest and most extensive genome-scale reconstruction of the entire metabolic network in E. coli K-12 MG1655, a well-researched and well-documented strain that is assumed to be metabolically equivalent to E. coli K-12 TG1. The network consists of 2,719 known metabolic reactions and 1,192 different metabolites including the reactions and metabolites that make up the endogenous paa pathway.2

While iML1515 contains almost all the essential metabolic reactions that occur in E. coli, it does not contain reactions for cellular mechanisms such as genetic/protein synthesis or regulation.1-3 Instead, the energy required for these reactions is accounted for by the ATP maintenance and Biomass pseudo-reactions, which reflect the bacterium’s ability to survive and its specific growth rate, respectively. With these pseudo-reactions, the network is able to accurately represent all the biochemical reactions in the cell.2,3

Editing the Metabolic Network

To represent the transformed biological device, we edited the iML1515 metabolic network to include new metabolites and reactions in the styrene degradation and PHB production pathways (Table 1). These edits were made using the Constraint-Based Reconstruction and Analysis (COBRA) Toolbox4, an open-source tool that allows users to construct, edit, and analyze metabolic networks through MATLAB, and the final genome-scale metabolic network was then visualized using Escher (Figure 1).5

Because a GEM separates all metabolites based on its location in the bacterium, the transport reactions for styrene across the cell wall and plasma membrane were also added.3 The transport of styrene from the extracellular compartment to the periplasmic space is assumed to be through porins by passive diffusion, and the transport of styrene from the periplasmic space to the cytosol is conducted by the primary active transporter coded by the styE gene.6

Exchange reactions for styrene and PHB were added to allow these metabolites to enter and leave the modelled system. Exchange reactions for other metabolites in M9 minimal medium were already present in the iML1515 GEM.3,7

Additionally, it was assumed that the energy required for the transcription and translation of the two plasmids was negligible, and therefore is accounted for by the given ATP maintenance and Biomass reactions because it is insignificant compared to the energy requirements already present in the cell.

Table 1. added Metabolic, Transport, and Exchange Reactions

Figure 1: visualization of the edited iM1515 genome-scale metabolic model

Flux Balance Analysis

The edited GEM was verified and analyzed using constraint-based Flux Balance Analysis (FBA), which is a computational approach to calculating the flux of metabolites through the network (Figure 2). The FBA model makes it possible to determine the growth rate of the device or the production of PHB under specific environmental and biological constraints provided by the M9 minimal medium with styrene.3,7

Figure 2. the mathematical basis for flux balance analysis. The model represents the reactions and metabolites in the genome-scale metabolic network as a system of linear equations that is defined by the conservation of mass and atoms at steady state. The system is then constrained by only allowing fluxes through exchange reactions with metabolites present in the bacterial growth medium, and by requiring a positive flow of metabolites through the ATP maintenance reaction to ensure cell survival. Finally, the solution space of the linear equations is optimized towards maximizing the flux of metabolites through the target reactions in the objective function.3

To solve the system of equations constructed by the FBA model, the bacterium is assumed to be at steady state. This means that the amount of metabolites flowing into the system has to equal the amount of metabolites flowing out. The FBA model does not show how the metabolite concentrations or rates change over time. Consequently, the fluxes have units of mmol/(gDW*hr), which is defined to be the amount of the metabolite in millimoles in a gram of bacterial cells (dry weight) over the time span of an hour. Negative flux values indicate metabolites flowing into the system while positive flux values indicate metabolites flowing out of the system. These fluxes were standardized to the specific growth rate of the bacteria to determine the theoretical yield as a molar ratio between styrene and PHB.3

GEM Verification

To confirm that the edited iML1515 GEM follows the laws of conservation of mass and energy, the model was constrained such that there was no influx of metabolites into the cell before running an FBA that is optimized towards the Biomass reaction. Then, an FBA was run after only adding styrene to the system, followed by similar FBAs with other important substrates such as oxygen, ammonium, and phosphate. Regardless of the metabolite that was added, the bacterial cell did not survive indicating that the model was correctly constructed with no mass or energy generating reactions.7

Theoretical Maximum Yield

With the FBA model constructed and verified, we analyzed our biological device to determine the maximum PHB yield. Obtaining a theoretical maximum yield will provide a reference to calculate and analyze the percent yield of PHB production from experimentation. It will also be a reliable marker of success for different iterations of our device’s design and experimentation.

To determine the molar ratio of styrene to PHB, the model was constrained to allow influx of all metabolites in M9 minimal medium and to force styrene uptake of 1 mmol/(gDW*hr). The maximum allowable flux for oxygen uptake was constrained to 20 mmol/(gDW*hr) to represent the actual oxygen uptake rate of E. coli.8 An FBA was then ran on the constrained model that was optimized towards the PHB exchange reaction. The maximum flux of 1.17 mmol/(gDW*hr) through the PHB reaction indicated that the molar ratio of styrene consumption to PHB production was 1:1.17.

However, when the model was constrained to a styrene uptake of 2 mmol/(gDW*hr), the molar ratio increased to 1:1.37, which suggests that the theoretical maximum yield of PHB is also dependent on other factors in the system. The change in ratio also showed that our simple method of using basic optimization was not effective for determining theoretical maximum yield. Therefore we improved our max-production method by conducting a robustness analysis between styrene and PHB.

A clear piece-wise function was observed from the robustness analysis graph (Figure 3). The function up to about 6 mmol/(gDW*hr) is linear with a non zero y-intercept as a result of the ATPM reaction always requiring a minimum amount of carbon. While the relationship is linear, the non zero y-intercept prevents us from assuming the ratio between PHB and styrene is always constant. Therefore, the ratio at the peak of the graph, 1.50 molar PHB per styrene, is the actual theoretical maximum production of the device.

Figure 3. robustness analysis of styrene consumption. Robustness analysis forces the consumption of styrene from its minimum to maximum possible values. Multiple FBAs were ran while optimizing toward PHB production at every value of styrene uptake.3 Under minimal medium conditions, PHB is produced most optimally at an influx of about 6 mmol/(gDW*hr). PHB is not produced at styrene uptake levels below 1 mmol/(gDW*hr) as the cells cannot survive. Forced styrene uptake beyond maximum production causes PHB production to decrease due to the energy required to manage the excess styrene in the system.


The piece-wise nature of the robustness graph indicates there is a metabolite other than styrene which is the primary limiting reactant. Oxygen was suspected to be the limiting metabolite in the system because the device’s biochemical pathway requires two moles of O2 to convert one mole of styrene to two moles of PHB. To test this hypothesis, double robustness analysis was run using styrene and O2 uptake as the two control reactions (Figure 4). The result of this analysis shows that maximum PHB production occurs at maximum oxygen uptake, supporting the hypothesis that oxygen is a limiting reactant of the biological device.

Figure 4: double robustness analysis for styrene and oxygen. Conducted similarly to single-variable robustness analysis, the uptake rates of two metabolites are forced and multiple FBAs are run while optimizing for PHB production.3 a) Under minimal medium conditions, PHB is produced most optimally at an influx of about 6 mmol/(gDW*hr) styrene and 20 mmol/(gDW*hr) oxygen, which is the maximum influx for oxygen in the constrained system. b) As shown in the vector plot, PHB production is primarily sensitive to variation in styrene and only a little sensitive to oxygen, yet oxygen is still the limiting metabolite due to the clear separation of phase planes.

Identifying oxygen as the limiting metabolite provides insight that increasing PHB production demands concomitant increases in oxygen uptake. The theoretical maximum yield was therefore revised to a molar sensitivity ratio of three variables: 1.50 : 3.30 : 1.00 PHB/oxygen/styrene. Using this new ratio, we will measure the observed consumption of styrene and oxygen as well as the production of PHB during experimentation. The addition of oxygen to this ratio will allow us to experimentally adjust the uptake of oxygen and/or styrene toward the maximum production ratio to improve the percent yield of our device more efficiently.

Sensitivity Analysis

Next, we determined the sensitivity of PHB production to certain metabolites found in the growth media. Previous studies showed that PHBs are synthesized maximally when the cell is growth limited by metabolites other than its carbon source. These studies found that the rate of PHB production by Pseudomonas putida is increased by deficiencies in nitrogen, phosphorus, magnesium, or sulfur.9 Using double robustness analysis, we tested the sensitivity of each metabolite on the production of PHB (Figure 5).

Figure 5: double robustness analysis to show sensitivity to metabolites found in the growth media. i, ii, iii, iv a)The y-component of the gradient at maximum styrene uptake is the criterion of sensitivity for each metabolite. i b) Nitrogen: Not very sensitive with a y-component of gradient at molar -0.2110 PHB/NH4. ii b) Phosphorus: Fairly sensitive with a y-component of gradient at molar -1.4328 PHB/Pi. iii b) Magnesium: Extremely sensitive with a y-component of gradient at molar -1,782 PHB/Mg2+. iv b) Sulfur: Not very sensitive with a y-component of gradient at molar -0.7112 PHB/SO4.

The results show production of PHB to be maximized when the influx of each metabolite was at 0, which supports the hypothesis that limiting these metabolites will cause a greater accumulation of PHB in the bacterium. Furthermore, the sensitivity of PHB production to each of these metabolites is useful when testing the actual device through experimentation. While nitrogen and sulfur uptake don’t necessarily have much effect on production, phosphorus and magnesium uptake rates should be precisely measured. Even a slight increase in magnesium uptake will have drastic effects on the percent yield of the device. Overall, knowing the sensitivities to these different metabolites gives us a method for error identification when determining limitations on percent yield.

Growth Coupling Analysis

Our final goal is to determine if PHB production can be coupled to growth of the device. Growth coupling is a synthetic biology process where product synthesis becomes metabolically paired with the device’s growth through genetic modifications that shift the substrate carbon demand toward synthesis of product.10 The main assumption for growth coupling is that all biological devices are evolutionarily primed to optimize their metabolic networks towards biomass production. Growth coupling therefore utilizes well established laboratory evolution methods for producing optimal microbial strains using growth as a selection criterion.11

The growth coupling method works by introducing gene deletions to shut off reactions which feed natural metabolic processes. As a result, the increased demand for byproducts of the target metabolite production reaction causes the bacterium to overproduce the target metabolite. If growth coupling is successful, the reaction phenotype will result in one of three different types as shown in Figure 6.11

Figure 6: production envelopes of three different growth coupling types. Production envelopes encompass all possible growth and production conditions, marked as the blue area. Growth coupling type is reliant on the lower bound, known as the minimum assured production, due to the assumption that cells will always maximize to biomass production. A) Weak growth coupling: the desired metabolite is only produced when growth rate is high. B) Holistic growth coupling: production rate is directly related to the growth rate. C) Strong growth coupling: production of the desired metabolite is a requirement for growth.11

Figure 7: bilevel optimization framework of OptKnock. The inner problem allocates flux distribution based on maximizing a cellular objective (biomass production). The outer problem then maximizes the production of the target metabolite by shutting off essential reactions available to the inner optimization problem.12

The COBRA function OptKnock was used to search for knockout sets. The function of OptKnock is to identify non-intuitive knockout sets which will growth couple a desired metabolite. After constraints are placed, the model is run through a bilevel optimization problem (Figure 7) to acquire a solution. The parameters for OptKnock include the reaction to be maximized (PHB exchange), number of knockouts, and a set of input reactions which a deletion set will be attained from. Please refer to Burgard et. al. for more background on the mathematics of the OptKnock function.12

Several trials of Optknock were run setting the number of deletions to 5, and varying the input reaction sets in order to search for a non-intuitive deletion set. Testing included each of the model’s individual subsystem, every combination of two subsystems, combinations of subsystems involved in the central metabolism of the model, and finally the entire reaction list of the model. Each OptKnock solution was analyzed using the analyzeOptKnock function.12 No trial produced a reaction set which resulted in the growth coupling of PHB.

Next we refined our method to a more intuitive approach by considering the reactions required for converting acetyl-CoA to PHB. Referring to the main premise that growth coupling occurs when the demand for the byproducts of the metabolite producing reaction is increased, we hypothesized that inhibiting the synthesis of CoA and/or NADP+ (the two byproducts of PHB production) could potentially growth couple PHB. To test this hypothesis, we first identified the reactions necessary for synthesizing CoA and NADP+ by running FBAs to maximize each byproduct’s production. The large lists of reactions involved in the synthesis of CoA and NADP+ were then used individually and combined as input reaction sets for OptKnock. The three resulting OptKnock solutions showed PHB production was not able to be growth coupled.

In addition to our results, Table 1 shows PHB production relies mainly on the overproduction of acetyl-CoA. Other studies on growth coupling metabolic models have shown acetyl-CoA as one of the only metabolites for which growth coupling is infeasible. This can be explained by the fact that acetyl-CoA is the main metabolite needed for running the cell’s TCA cycle, the primary metabolic subsystem for the production of growth associated energy.11 Furthermore, looking at the byproducts of PHB production, FBA’s optimization process will force both CoA and NADP+ to be recycled just to produce more acetyl-CoA for consumption. Therefore, trials of OptKnock most likely failed to produce a growth coupled knockout set because acetyl-CoA is a metabolic branch point 一 a place where the cell will either use the metabolite for the production of PHB or for the growth of the cell. While it is possible there exists a yet unidentified knockout set which would cause growth coupling, our results, those from other OptKnock studies, and the analysis of the system’s stoichiometry support the conclusion that PHB production cannot be growth coupled by our device.

Knowing PHB production is incapable of growth coupling provides significant insight into the design and future implementation of our device. First, the time scale for bacterial growth and PHB production has to be small enough such that evolution does not occur. This indicates that our device has to be grown in an unsteady state system provided by a batch or semi-batch reactor rather than the steady state system of a continuous bioreactor.13 To further reduce evolutionary and metabolic strain on the device, the promoter of the PHB production genes should be induced at a specific time point on the bacterial growth curve while the promoter of the styrene degradation genes should be constitutively expressed. The theoretical time point for inducing PHB production will be the time when the bacterial culture reaches carrying capacity due to production at maximum cell density while there are limited metabolites for biomass growth present in the medium.

References

  1. Monk JM, Lloyd CJ, Brunk E, Mih N, Sastry A, King Z, et al. iML1515, a knowledgebase that computes Escherichia coli traits. Nature Biotechnology. 2017 Oct 11;35:904.
  2. Orth JD, Thiele I, Palsson BØ. What is flux balance analysis? Nat Biotechnol. 2010 Mar;28(3):245–8.
  3. Heirendt L, Arreckx S, Pfau T, Mendoza SN, Richelle A, Heinken A, et al. Creation and analysis of biochemical constraint-based models using the COBRA Toolbox v.3.0. Nature Protocols. 2019 Mar 1;14(3):639–702.
  4. King ZA, Dräger A, Ebrahim A, Sonnenschein N, Lewis NE, Palsson BO. Escher: A Web Application for Building, Sharing, and Embedding Data-Rich Visualizations of Biological Pathways. Gardner PP, editor. PLoS Comput Biol>. 2015 Aug 27;11(8):e1004321.
  5. Mooney A, O’Leary ND, Dobson ADW. Cloning and Functional Characterization of the styE Gene, Involved in Styrene Transport in Pseudomonas putida CA-3. Applied and Environmental Microbiology. 2006 Feb 1;72(2):1302–9.
  6. Monk JM, Charusanti P, Aziz RK, Lerman JA, Premyodhin N, Orth JD, et al. Genome-scale metabolic reconstructions of multiple Escherichia coli strains highlight strain-specific adaptations to nutritional environments. Proc Natl Acad Sci USA. 2013 Dec 10;110(50):20338.
  7. Garcia-Ochoa F, Gomez E, Santos VE, Merchuk JC. Oxygen uptake rate in microbial processes: An overview. Biochemical Engineering Journal. 2010 May 15;49(3):289–307.
  8. Braunegg G, Lefebvre G, Genser KF. Polyhydroxyalkanoates, biopolyesters from renewable resources: Physiological and engineering aspects. Journal of Biotechnology. 1998 Oct 27;65(2):127–61.
  9. Van Dien S. From the first drop to the first truckload: commercialization of microbial processes for renewable chemicals. Current Opinion in Biotechnology. 2013 Dec 1;24(6):1061–8.
  10. Alter TB, Ebert BE. Determination of growth-coupling strategies and their underlying principles. BMC Bioinformatics. 2019 Aug 28;20(1):447.
  11. Burgard AP, Pharkya P, Maranas CD. Optknock: A bilevel programming framework for identifying gene knockout strategies for microbial strain optimization. Biotechnol Bioeng. 2003 Dec 20;84(6):647–57.
  12. Blunt W, Levin D, Cicek N. Bioreactor Operating Strategies for Improved Polyhydroxyalkanoate (PHA) Productivity. Polymers. 2018 Oct 26;10(11):1197.