8CM00 Systems Biology
Assignment 1 Part 2: Pathways and graphs, September 19, 2016
Peter Hilbers
The assignment of the part of Hilbers consists of 2 parts. Below the description of part 2 is given. Previous week part 1 about sequence alignments has been given.
The intake of food of different composition yields different amounts of energy. For instance, the energy yield from a gram of fatty acids is approximately 37 kJ, compared to 17 kJ for carbohydrates. In this assignment we are studying what the maximum amount of energy is that can be produced because of a certain food intake. Next to it we also compare the major differences in metabolic pathways between two different food intakes.
On March 3, 2013, in a Nature Biotechnology paper ([1]) Recon 2 has been launched. Recon 2, the second version of the Reconstruction of The Human Genome is one of the largest constraints-based metabolic models to date with 7440 re- actions, 2169 genes and 5764 compartmented metabolites. Recently, July 2016, an improved version of Recon 2 has been published. This model, Recon2.2 ([2]) is available in SBML format and stored on Canvas. To process this file both Python as Matlab can be used. For the Matlab usage a file ’Instructions to in- stall and initialize the COBRA Toolbox in Matlab’ has been added. (The SBML files can be also translated into a Matlab structure by using the SBML toolbox http://sbml.org/software/sbmltoolbox. Detailed description of the toolbox can be found at http://sbml.org/Software/SBMLToolbox/docs/ManualV3.pdf).
For the exercise description we use Python and a Python program, ’recon2 2.py’, is available in which the basic elements are given. To assist you in your task a detailed description of the relevant parts of ’recon2 2.py’ are given. Although they are given in Python the description should be easily transferable to other programming tools and languages (including Matlab).
1. Reading the file and constructing the model
1-1
1-2
Assignment 1 Part 2: Pathways and graphs, September 19, 2016
By executing run_model with ’recon 2.2.xml’ in the d.models-folder the fol- lowing piece of code is run:
sbml = read_sbml(model_filename)
model_stats(sbml, display_errors=True)
In the first line the sbml-file is read and an sbml object constructed. Next in model_stats(sbml, display_errors=True) the model is analyzed in terms of the number and type of species and reactions.
Exercise 2.1: Write a program to read the sbml-file and to construct the model.
The output should be:
recon_2.2
6047 species
5324 variable
723 fixed
0 non-reactants
So in total the model consists of 6047 metabolites of which 723 have a boundary condition. A list of all these species can be obtained by calling model.getListOfSpecies(). Each list entry is a species-object. Such a species- object has a large number of attributes. Most of them can be accessed by calling the appropriate get-method. For example
sID=’M_glc_D_b’
myspecies=model.getSpecies(sID)
myspecies.getId() # M_glc_D_b is the id for “D-glucopyranose”
myspecies.getName() # ’D-glucopyranose’
myspecies.getBoundaryCondition() # True
Also in model_stats the reactions of the model are analysed. Each reaction converts a number of reactants into a number of prod- ucts, where both the reactants and products are references to species. Similarly as for the species all reactions can be obtained by model.getListOfReactions(). The reactants of a reaction in this list can be accessed by for reactant in reaction.getListOfReactants() and to ob- tain the species’ id of a reactant one can use sID = reactant.getSpecies().
Assignment 1 Part 2: Pathways and graphs, September 19, 2016 1-3
Having the species’ id, all methods described above for species can be applied. For each of the reactants and products, the species reference also contains a field called stoichiometry that can be accessed by reactant.getStoichiometry(). It states the net number of reactants and products that are consumed and produced, respectively, by the reaction. For instance, for the biomass reaction a reactant is the M_biomass_protein_c thas a stoichiometry of 0.706.
For each reaction the flux is limited by a lower and upper bound. The smallest lower bound is −∞, the largest upper bound being +∞. When the lower bound is negative and the upper bound positive, the reaction is reversible. Otherwise the reaction is irreversible. Notice that in flux calculations the sign of the flux value determines the direction.
Reactions are divided into 2 types. A reaction for which either all reactants or all products have a boundary condition is considered a source and sink reaction, respectively. For recon2.2 we have:
7785 reactions
7038 non-source/sink
747 source/sink
Exercise 2.2: Write a program that determines and prints to output the total number of reactions and the number of source and sink reactions.
2. Model balancing
The next step in running the model is the model_balancing. For all non- source/sink reactions it is checked whether the law of conservation of mass holds meaning that the total mass of the reactants equals the total mass of the products. In this calculation the stoichiometric number of each of the species is of course used. In our case recon 2.2 turns out to have 5 unbalanced all related to the artifical biomass reaction.
Exercise 2.3: Write a program that determines and prints the unbalanced reac- tions.
3. Flow calculations
The last step in running the model is the execution of max_fluxes(sbml). This method forms the core to construct the graphs needed for the pathway analysis and by this method the maximal amount of ATP (in flux units i.e. moles per second) is calculated that is produced by given food intake. This is done by setting the objective function for the max-flux calculation to ATP production:
1-4
Assignment 1 Part 2: Pathways and graphs, September 19, 2016
objective = ’DM_atp_c_’
print ’’
In the current version a diversity of carbon sources are given and a maximal flux calculation is performed for each carbon source in isolation, for instance to calculate the maximum amount of ATP production for a source consisting of glucose only ’EX glc(e)’ is set to 1 and all other carbon sources are set to 0.
Note: In recon 2.2.xml reactions have as prefix R and metabolites M . When a reaction is denoted as ’EX glc(e)’ it is written in recon 2.2.xml as ’R EX glc(e)’. In the old notation that is used because of legacy problems a ( is denoted as LPAREN and ) as RPAREN , so reaction ’EX glc(e)’s id is “R EX glc LPAREN e RPAREN ”
Apart from carbon sources there are also other kinds of sources, called media. such as Ca2+ and H2O. In the calculations it is assumed that there is an unlimited supply of these media sources, i.e. their flux upper bound is set to +∞.
The set of media is:
media = [ ’EX_ca2(e)’, ’EX_cl(e)’, ’EX_fe2(e)’, ’EX_fe3(e)’,
’EX_h(e)’, ’EX_h2o(e)’, ’EX_k(e)’, ’EX_na1(e)’,
’EX_nh4(e)’, ’EX_so4(e)’, ’EX_pi(e)’ ]
while the carbon sources are:
carbon_sources=[ # sugars
’EX_glc(e)’, ’EX_fru(e)’,
# fatty acids
’EX_ppa(e)’,
’EX_but(e)’,
’EX_hx(e)’,
’EX_octa(e)’,
’EX_HC02175(e)’,
’EX_HC02176(e)’,
’EX_ttdca(e)’,
’EX_hdca(e)’,
’EX_ocdca(e)’,
’EX_arach(e)’,
’EX_docosac_’,
# C3:0
# C4:0
# C6:0
# C8:0
# C10:0
# C12:0
# C14:0
# C16:0
# C18:0
# C20:0
# C22:0
Assignment 1 Part 2: Pathways and graphs,
September 19, 2016 1-5 # C24:0
# C26:0
’EX_lgnc(e)’,
’EX_hexc(e)’,
# amino acids
’EX_ala_L(e)’, ’EX_arg_L(e)’, ’EX_asn_L(e)’,
’EX_asp_L(e)’, ’EX_cys_L(e)’, ’EX_gln_L(e)’,
’EX_glu_L(e)’, ’EX_gly(e)’, ’EX_his_L(e)’,
’EX_ile_L(e)’, ’EX_leu_L(e)’, ’EX_lys_L(e)’,
’EX_met_L(e)’, ’EX_phe_L(e)’, ’EX_pro_L(e)’,
’EX_ser_L(e)’, ’EX_thr_L(e)’, ’EX_trp_L(e)’,
’EX_tyr_L(e)’, ’EX_val_L(e)’
]
In the max-flux calculation not only the maximum amount of ATP, f opt, is calculated but also for each reaction the flux. These fluxes are stored in the list v:
v, f_opt = max_flux(sbml, carbon_source, objective, normoxic, media)
Exercise 2.4: Change this code such that flux calculation can be performed with each carbon source having an arbitrary input value (so that the calculation can be done not for a single carbon source as food intake, but for a mix of carbon sources as input). Run the calculation for the ’standard’ food intake of each carbon source having the input value 1.
Exercise 2.5: Construct based upon the output fluxes of the ’standard’ food intake max-flux calculation (see previous exercise) the ’default food’ directed graph DF . The nodes of DF are the metabolites of ’recon2.2’ and two arbitrary metabolites R and P have an edge from R to P if there is reaction with a positive flux in which R occurs as reactant and P as product (for a negative flux the direction has to be reversed). How many metabolites have a non-zero degree?
Exercise 2.6: Add to graph DF an additional node, called ’food intake’ with an edge from ’food intake’ to each carbon source. The flux on this edge should be the input value of the carbon source. The next step is to omit from this new graph all nodes with zero degree (so a node having no incoming and no outgoing edges). This resulting graph is called DFfi0.
Calculate for this graph DFfi0 a ’maximum spanning tree’ MST0 with as root ’food intake’. When a node is reachable from ’food intake’ its distance is the largest sum of the fluxes along a (directed) path from ’food intake’ to this node.
1-6
Assignment 1 Part 2: Pathways and graphs, September 19, 2016
Exercise 2.7: By email you will also receive another food intake pattern. Repeat the previous 2 exercises for this new food intake. The new graph should be called DFfi1 and you are being asked to discuss the differences between DFfi0 and DFfi1. These differences should include:
• the degree distribution
• the average degree
• the variance of degree
• rank by degree: hubs
• the distance density function of the ’food intake’ node, i.e, the maximum spanning tree for the ’food intake’ node as root.
References
[1] Thiele et al, A community-driven global reconstruction of human metabolism, Nature Biotechnology, 3 March 2013; doi:10.1038/nbt.2488
[2] Swainston et al, Recon 2.2: from reconstruction to model of human metabolism, Metabolomics (2016) 12:109 DOI 10.1007/s11306-016-1051-4
Submission by email: p.a.j.hilbers@tue.nl
Deadline for submission part 1 and part 2: October 17, 2016, 23:59:59