Genome Scale Modeling
One of our main objectives for this project is to build a comprehensive, genome-scale model that details the function of the Lux system and its interaction with the host organism, Saccharomyces Cerevisiae (s288c).
Genome-scale models (M-Models) represent an effort to bridge the gap between isolated systems (e.g. signaling pathways, glycolysis) to modeling whole systems of organisms1. Metabolic network reconstructions input the full range of known biochemical reactions in an organism to give comprehensive insight into how the components of the genome result in a particular phenotype. These elaborate models are usually encoded in the Systems Biology Markup Language (SBML).
There are currently many M-Models available in literature, mostly detailing the networks of different strains of bacteria and the simpler eukaryotes2. The yeast strain of Saccharomyces Cerevisiae (s288c) was one of the first to be converted into a network model. For this and for reasons of experimental simplicity, we chose this strain as our base model.
The UC system-wide biochemically, genetically, and genomically (BiGG) database features two different networks for the genome of Saccharomyces Cerevisiae. The latest model reconstruction, iND750, consists of 1059 metabolites, 1266 reactions, and 750 genes2. In early 2007, a community of researchers came together to standardize the metabolic map of Saccharomyces Cerevisiae. They collaborated to periodically update the model and make it available in the public domain (yeast.sf.net) for the years that followed. The latest version to date, version 7, includes reactions published by Aung et al. in 20133.
Ultimately, we decided to model the Lux system in both available models – BiGG and v7. After evaluating the strengths and weaknesses of both models, we were not convinced that either model provided a complete and accurate representation of the yeast metabolism. The v7 model would seem to be the logical choice, because it is frequently used and updated by the yeast community; however, we were concerned when we discovered many reactions and metabolites missing from the native fatty acid biosynthesis pathway. This pathway is crucial to the Lux system, because luxD interacts with myristic-ACP, also known as tetradecanoyl-ACP4,5. With the v7 model terminating fatty acid biosynthesis at octanoyl-ACP, we were forced to add the necessary reactions with the assistance of the KEGG database4,5.
Metabolic Control Analysis and Flux Balance Analysis
The term metabolic flux is used to describe the rate of turnover of metabolites in these large-scale networks. Flux for an individual reaction is the difference between the rates of the forward reaction and the reverse reaction6.
The main objective of Flux Balance Analysis (FBA) is to find a solution for a system of linear equations, (S*v = b). Using the coefficients -1,0, and 1, we can identify which metabolites actively participate in each reaction and form the stoichiometric matrix, S. b is assumed to be a zero matrix, because at steady-state, the rate of consumption is equal to the rate of production. Hence, flux balance analysis returns v, the column vector of the flux of each reaction solved at steady state7.
However, FBA has many limitations. FBA cannot predict metabolic concentrations, because it does not utilize kinetic parameters. Also, it cannot account for regulatory agents such as those involved in transcription or translation of the genome7,8.
Constraint-Based Modeling
We used the COBRA Toolbox plugin in MATLAB to edit and analyze our SBML model9. The steady-state assumption that is employed in FBA is one example of a constraint. Constraints define the allowable solution space in FBA and accounts for many reaction properties that cannot be encoded by simple stoichiometry. Metabolic constraints, such as flux coefficients, define the physico-chemical boundaries of a reaction. The lower bound of the flux coefficient also tells us if the reaction is reversible or irreversible10.
Space constraints, such as the use of compartments, define the physical boundaries on which metabolites can react with one another. There are many other types of constraints that can be embedded into a genome-scale model such as regulatory constraints, thermodynamic constraints, and temporal constraints. All parameters seek to eliminate infeasible solutions and define the solution space as closely as possible1,10.
The most important constraint is the objective function. The COBRA Toolbox allows us to maximize or minimize the flux through certain reactions of interest1. For instance, setting the objective to biomass production allows us to predict the growth rate of a yeast cell once it reaches steady-state. For our purposes, maximizing the flux through the light reaction was set as the main objective.
Results
In both models, we first defined a minimum media to mimic experimental conditions. The components of Yeast Extract Peptone Dextrose (YPD) were listed and linked to exchange reactions in both the v7 and the BIGG model. Exchange reactions are defined to describe the metabolite composition in the immediate surroundings of the cell. The upper and lower bounds for all other exchange reactions were then set to ‘0’.
Once the Lux reactions were inserted into both models, we analyzed the metabolic burden on the cell using a two-step process. First, we ran the built-in COBRA function robustnessAnalysis. This function plots the flux through a pathway against the flux through the optimized pathway1. In this case, the flux through the light pathway was negatively correlated with the flux through the growth pathway. We concluded that this result could only have been obtained if the light pathway was redirecting a metabolite that would normally have lent itself to growth.
To confirm this hypothesis, we needed to identify which reactions in the native genome were coupled to light. We developed a script that set light as the objective function, and then constrained growth in increments until it reached its max value. At each growth value, we took the vector of all the fluxes through every reaction. In the end, each reaction fell into one of four categories: incr, decr, zero, and neither. In incr, we had a list of reactions strongly coupled to growth. The flux through these reactions increased at every point that growth increased. Likewise, in decr, we held a list of reactions strongly coupled to light. In zero, we stored a list of reactions that were uncoupled to growth and light. In neither, we stored reactions, whose flux varied with no discernible pattern.
The BiGG model returned 34 reactions that were strongly coupled to light and the v7 model returned 56 reactions. The transport reactions in this category were then removed to better ascertain which subsystems (Glycolysis, Citric Acid cycle etc.) these reactions belonged to4,5. The diagram below traces the light pathway to its origin and shows its progression into the fatty acid biosynthesis pathway. This diagram allowed us to understand how the light pathway redirected metabolites from the growth pathway.
Finally, to confirm our results, we ran a robustnessAnalysis on flux through the light pathway versus that through transketolase, an enzyme in the pentose phosphate pathway. We were expecting to see the flux of light monotonically increase, as predicted by the previous sorting algorithm. However, we observed a peak in the graph of the flux through light and a subsequent decline in light with increasing transketolase activity. This effect can only be explained by the lack of constraints on the flux through the growth pathway. From this, we conclude that past a threshold value, the metabolite produced through this pathway is redirected into the growth reaction.
Future Directions
We can use experimental data to further validate and refine predictive capabilities of our model. A measured growth curve, as well as values for nutrient uptake and secretion can be used as constraints on the flux of their respective reactions. We can also compare the flux value for light units given by the FBA function to that of experimental results.
Dynamic FBA is a technique that solves for the mathematical model at certain snapshots in time instead of at steady-state. It can also factor in concentrations of reactants and kinetic rate constants, which can make the predictive results of the model closer to experimental values8.
Scripts
Link to Constraint-Based Modeling Scripts
References
[1] Schellenberger J, Que R, Fleming RMT, Thiele I, Orth JD, Feist AM, Zielinski DC, Bordbar A, Lewis NE, Rahmanian S, Kang J, Hyduke DR, Palsson BØ.
[2] Duarte, Natalie C., Markus J. Herrgård, and Bernhard Ø. Palsson. “Reconstruction and Validation of Saccharomyces Cerevisiae iND750, a Fully Compartmentalized Genome-Scale Metabolic Model.” Genome Research 14.7 (2004): 1298–1309. PMC. Web. 18 Sept. 2015.
[3] Herrgard, M et al. “A consensus yeast metabolic network reconstruction obtained from a community approach to systems biology.” Nat Biotechnol 26 (2008) : 1155-1160.
[4] Kanehisa, M., Goto, S., Sato, Y., Kawashima, M., Furumichi, M., and Tanabe, M.; Data, information, knowledge and principle: back to metabolism in KEGG. Nucleic Acids Res. 42, D199–D205 (2014).
[5] Kanehisa, M. and Goto, S.; KEGG: Kyoto Encyclopedia of Genes and Genomes. Nucleic Acids Res. 28, 27-30 (2000).
[6] Voet, Donald, and Judith G. Voet. Biochemistry. New York: J. Wiley & Sons, 1995. 439. Print.
[7] Orth, Jeffrey D., Ines Thiele, and Bernhard Ø. Palsson. "What is flux balance analysis?." Nature biotechnology 28.3 (2010): 245-248.
[8] Vincent J. Henry; Anita E. Bandrowski; Anne-Sophie Pepin; Bruno J. Gonzalez; Arnaud Desfeux
Database 2014; doi: 10.1093/database/bau069
[9] Becker, S. et al., "Quantitative prediction of cellular metabolism with constraint-based models: The COBRA Toolbox", Nat. Protoc 2, 727-738 (2007).
[10] Lewis NE, Nagarajan H, Palsson BO. Constraining the metabolic genotype-phenotype relationship using a phylogeny of in silico methods. Nature Reviews Microbiology. (2012)