Flux balance analysis on yeast metabolic models for data reuse

1. Goal
To verify the state-of-the-art systems biology model for Saccharomyces cerevisiae metabolic networks using flux balance analysis (FBA), with the experimental data produced from a formalized functional genomics study [King09]. 

2. Introduction to FBA
FBA is a constraint-based approach using linear programming to identify a flux distribution that optimize the given objective function [Edwards01]. FBA applies mass and energy balance constraints to model steady state behavior of reconstructed networks, assuming the production and utilization rate of a metabolite must be a balance in a steady state. Further utilized constraints involve irreversibility of reactions and the maximum flux through any reaction or transporter. The network of reactions is uniquely defined by a stoichiometric matrix S, whose element represents the stoichiometry of metabolite in the reactions. Exchange fluxes for boundary metabolites that are allowed to be in and out of the cell boundary are also incorporated in the S matrix. The steady state constraints for FBA is thus expressed as: Sv = 0 and vmin_i ≤ v_i ≤ vmax_i, i = 1,..., N, where v is the vector of N reaction fluxes. vmin_i and vmax_i correspond to the lower and upper bound for individual flux v_i. 

3. Method
3.1. Reference model
We used the published genome scale metabolic network model for S. cerevisiae, namely iND750 model, by Duarte et al. [Duarte04] as our reference model, which is downloadable from website http://systemsbiology.ucsd.edu/In_Silico_Organisms/Yeast. The biomass formation reaction was also taken from Duarte et al.. 

3.2. Model constraints for optimization
To simulate the model, proper model constraints have to be defined. This involves in setting up the constraints for maximum uptake rate for nutrients according to the growth media and lower and upper bound for reactions according to different genotypes. 

The maximum uptake rates for the exchange metabolites for rich media (YPD, link: ) and minimal media (MM, link: ), were both adapted from Snitkin et al.. Methionione was added to the minimal media in order to match the our real experimental condition. 

For each of the 20 genes in the mutants investigated, flux balance analysis was done to simulate the change of the flux distribution with deletion. To account for the hypothesised effect of deletion of a gene, the associated reaction(s) were disallowed by setting lower and upper bound for those fluxes to 0. In a standard gene deletion simulation using FBA, 1) in case of enzyme complex encoded by a set of genes, reaction will be all disallow if one of the genes is deleted; 2) in case of iso-enzymes encoded by different individual genes, the reaction will be intact unless all the iso-enzymes have been deleted. In this simulation study, all enzymes annotated with multiple genes have been treated as enzyme complexes here. 

To model the genome of the wild-type the reactions catalysed by the genes: YEL021W (URA3), YLR303W (MET15), YCL018W (LEU2), and YOR202W (HIS3) were disallowed for both the wild-type and deletion simulations. Among the 13 enzymes (and associated reactions) with hypothesised annotation in the study, 11 are present in iND750 model and 10 of them have different or no annotations at all, two are absent in iND750 and not being able to verify with the model. Thus these annotations iND750 model has been modified according to the hypotheses. 

3.3. Flux range analysis
Flux range analysis was performed to check the possible flux variability for the studied reactions satisfying two slightly different set of constraints: (1) structural constraints only, (2) with one additional constraint of a minimum growth rate of 0.05. The flux range was then obtained by finding the minimum and maximum fluxes for the reaction of interest individually using the wild-type model under the two sets of constraints. The analysis was conducted for both the minimal defined and YPD media. 

3.4. Growth rate prediction for gene deletion mutants
We then performed standard FBA simulation to predict the maximum growth rate for the wild-type and deletion mutants by maximising the flux for the biomass formation reaction while satisfying all the constraints. The difference of the maximal growth rate between the deletants and wild-type was calculated and compared with our experimental data. According to the FBA model, the flux in general have unit of mmol/gDW/hr; the flux value for biomass formation reaction can be interpreted as the growth rate with unit of 1/hr, whose absolute value is not directly comparable to the growth rate estimated from the growth curves in our experiments. However, as we focus on comparing the change of direction (increase or decrease) in the growth rate between the wild-type and mutants only, a normalisation from the predicted growth rate is not necessary here. The experimental difference was estimated by the median growth rate changes, i.e. the growth rate for deletion mutants minus the growth rate for the wild-type from different experiments. 

3.5. Computational tools for FBA simulation
The FBA simulation was conducted in Python, utilizing the open source constraint-based optimization packages such as lpSolver and optimization framework OpenOpt. The python tool for reading the SBML model and FBA simulation can be found following this link xxx. The python script for FBA simulation of iND750 model for data reuse study can be found here xxx. To run the script, installation of the following packages are required: 
	•	Numpy (http://numpy.scipy.org/)
	•	libSBML and its python binding (http://sbml.org/Software/libSBML)
	•	lpSolver and its python binding (http://lpsolve.sourceforge.net/5.5/)
	•	OpenOpt (http://openopt.org/)

4. Results
Table 1 shows the FBA model's range analysis. The analysis shows that the expected effect of perturbing the flux of the reactions examined on growth will be small. Four of these reactions are blocked in the model and are unable to carry any flux at all. For five of the reactions the flux range for minimal growth is large and encompasses zero. Therefore the model predicts that elimination these two types of reactions will not prevent the cell from growing. The model also predicts that the removal of two of the reactions (involving 4 gene annotations) will stop the cell growing as the minimum growth requires their flux to be positive.

Table 2 shows the comparison of the FBA predicted growth rate changes compared with the experimental measured data. The kinetic details are not yet available for those reactions of interest that involve in multiple enzymes (encoded by different genes). 

Assuming that the enzymes all play the same role for the particular reaction that the multiple genes. 
The model thus predicts no change for the single gene deletion mutant that actually. Thus here we report the model predictions assuming the enzymes involve in multiple genes are all enzyme complexes, that is, all relevant genes need to be present in the cell in order to have the associated reaction catalysed. 

As we can see that a standard FBA simulation is not possible to predict the subtle change of these deletion mutants over wild-type. 
There are several related reasons for the inability of the FBA model to model the observed growth:
	•	With FBA models it is impossible to obtain a higher maximal growth rate through the removal of reactions from the model. 
	•	The reactions investigated are mostly non-essential based on the flux range analysis.
	•	The FBA model is known to be both incomplete and incorrect: with expected mistakes in the network structure, the definition of biomass formation, and the limits in nutrient uptake rates.
	•	Standard FBA models cannot model the perturbation of removing isoenzymes. 

We have also tried using the minimisation of metabolic adjustment (MOMA) approach to modelling deletion growth rate (see further information). However, similar results were obtained as with standard FBA, and the model still failed to explain the experimental observation that most deletant strains grow faster than wild-type. 

         
References

[King09] King, R. D., Rowland, J., Oliver, S. G., Young, M., Aubrey, W., Byrne, E., Liakata, M., Markham, M., Pir, P., Soldatova, L. N., Sparkes, A., Whelan, K. E., Clare, A. The Automation of Science. Science, 2009, 324(5923):85-89.

[Edwards01] Edwards JS, Ibarra RU, Palsson BO: In silico predictions of Escherichia coli metabolic capabilities are consistent with experimental data. Nat Biotechnol 2001, 19:125-130.

[Duarte04] Duarte, N. C., Herrgard, M. J. & Palsson, B. Ø., Reconstruction and validation of Saccharomyces cerevisiae iND750, a fully compartmentalized genome-scale metabolic model. Genome Res. 14, 1298–1309 (2004).

[Snitkin08] Evan S. Snitkin, Aimee M. Dudley, Daniel M. Janse, Kaisheen Wong, George M. Church, and Daniel Segrè: Model-driven analysis of experimentally determined growth phenotypes for 465 yeast gene deletion mutants under 16 different conditions, Genome Biology (2008) 9:R140.

Table 1: Flux range analysis for the reactions involved in the hypothesised gene annotations. 
		Media		Minimal	Minimal	YPD	YPD	
studied enzyme	Reaction ID in iND750	Hypothesised associated genes	Reaction reversible 	Range for structural validity	Range for minimal growth	Range for structural validity	Range for minimal growth	
2-aminoadipate transaminase (2.6.1.39)	R_AATA	YER152C, 
YGL202W, 
YJL060W	Yes	[0, 7.4]	[0.01, 7.0]	[0, 8.1]	[0, 7.7]	
1-acylglycerol-3-phosphate O-acyltransferase (2.3.1.51)	R_AGAT_SC	YDL052C	No	[0, 0.02]	[0.001, 0.02]	[0, 0.02]	[0.001, 0.02]	
5,10-methenyltetrahydrofolate synthetase (6.3.3.2)	R_FTHFCLm	YER183C	No	[0, 36.1]	[0, 34.7]	[0, 36.9]	[0, 35.6]	
glucosamine-6-phosphate deaminase (3.5.99.6)	R_G6PDA	YHR163W, 
YNR034W, 
YGR248	No	[0, 61.5]	[0, 57.3]	[0, 62.7]	[0, 59.0]	
glutaminase (3.5.1.2)	R_GLUN	YIL033C	No	[0, 61.5]	[0, 57.3]	[0, 62.7]	[0, 59.0]	
mannitol-1-phosphate 5-dehydrogenase (1.1.1.17)	R_M1PD	YNR073C	Yes	0	0	0	0	
maleylacetoacetate isomerase (5.2.1.2)	R_MACACI	YLL060C	No	0	0	0	0	
polyamine oxidase (1.5.3.11)	R_POLYAO2	YMR020W	No	0	0	0	0	
purine-nucleoside phosphorylase (2.4.2.1)	R_PUNP1	YLR209C, 
YLR017W	Yes	[-61.5, 30.7]	[-57.3, 28.7]	[-62.7, 31.4]	[-59.0, 29.5]	
pyridoxal kinase (2.7.1.35)	R_PYDXK	YNR027W, YPR121W	No	0	0	0	0	
serine O-acetyltransferase (2.3.1.30)	R_SERATi	YJL218W	No	[0, 11.6]	[0, 10.9]	[0, 11.9]	[0, 11.3]	
L-threonine 3-dehydrogenase (1.1.1.103)	N/A*	YDL168W, 
YLR070C	-	-	-	-	-	
L-aspartate oxidase (1.4.3.16)	N/A*	YJL045W	-	-	-	-	-	
*N/A means the reactions not being present in iND750 model and thus the consequent flux range analysis results are not available. 

Table 2: Comparison of the predicted maximum growth rate change (mutant - wild-type) with the measured experimental growth rate change for the 20 manually studied gene deletion mutants. 
Reaction ID in iND750	Deletion genes	Experiment_MinimalMedia	Simulation_MinimalMedia	Experiment_YPDMedia	Simulation_YPDMedia	
R_AATA	YER152C	0.009	-0.733	0.019	-0.222	
R_AATA	YGL202W	-0.024	-0.733	0.024	-0.222	
R_AATA	YJL060W	0.013	-0.733	0.024	-0.222	
R_AGAT_SC	YDL052C	0.009	-0.733	0.034	-0.805	
R_FTHFCLm	YER183C	0.022	0.000	0.014	0.000	
R_G6PDA	YGR248W	0.017	0.000	0.007	0.000	
R_G6PDA	YHR163W	-0.222	0.000	0.005	0.000	
R_G6PDA	YNR034W	0.023	0.000	0.028	0.000	
R_GLUN	YIL033C	-0.079	0.000	-0.205	0.000	
R_M1PD	YNR073C	0.016	0.000	0.024	0.000	
R_MACACI	YLL060C	0.011	0.000	0.014	0.000	
R_POLYAO2	YMR020W	0.016	0.000	0.023	0.000	
R_PUNP1	YLR017W	0.003	0.000	0.008	0.000	
R_PUNP1	YLR209C	0.017	0.000	0.004	0.000	
R_PYDXK	YNR027W	0.013	0.000	0.023	0.000	
R_PYDXK	YPR121W	0.036	0.000	0.025	0.000	
R_SERATi	YJL218W	0.015	-0.733	0.038	0.000	
N/A*	YDL168W	0.018	0.000	0.024	0.000	
N/A*	YJL045W	0.016	0.000	0.000	0.000	
N/A*	YLR070C	0.012	0.000	0.019	0.000	
*N/A means the reactions not being present in iND750 model and thus simulated flux distribution of knockout will show no difference with the wildtype.