Description of tutorial1

0001 % tutorial1
0002 %   This is a short introduction that shows how to load a genome-scale
0003 %   metabolic model (GEM), set reaction constraints, objective function and
0004 %   perform an optimization through flux balance analysis (FBA). The
0005 %   resulting fluxes are visualized and exported to a PDF file.
0006 %   For a more detailed description of the individual functions, see
0007 %   [raven_directory]/doc/index.html.
0008 %   A GEM for the filamentous fungus Penicillium chrysogenum is used in
0009 %   this tutorial. The model can be found in a Microsoft Excel file under
0010 %   the name "iAL1006 v1.00.xlsx" and in SBML file "iAL1006 v1.00.xml".
0011 %   See Tutorial 1 in "RAVEN tutorials.docx" for more details.
0012 
0013 %Import the model from Excel. This function performs a number of checks
0014 %regarding the model structure (such as for incorrectly written equations
0015 %or illegal characters). In this structure there is only one warning; that
0016 %the formula for the metabolite LPE could not be parsed. The "false" flag
0017 %imports a model with exchange reactions in their "closed" form. This makes
0018 %the model unsuited for modelling, but it is useful for some quality
0019 %control steps.
0020 model=importExcelModel('iAL1006 v1.00.xlsx',false);
0021 
0022 %The Excel interface is supposed to work in all the systems (Windows, Unix,
0023 %macOS), but upon any problems, the model can be imported from SBML format
0024 %instead. However, in such case the user would not be able to run Tutorials
0025 %2-4, which involve the editing of Excel files. Run the command below
0026 %(remove "%" sign) instead, if having such problem:
0027 %model=importModel('iAL1006 v1.00.xml',false);
0028 
0029 %The following function prints some properties of the model. The two "true"
0030 %flags say that it should also list potential problems such as dead-end
0031 %reactions or unconnected metabolites.
0032 printModelStats(model,true,true);
0033 
0034 %As can be seen the model contains 1632 reactions, 1395 metabolites and
0035 %1006 genes
0036 
0037 %Most modelling approaches using GEMs are based on the mass balancing
0038 %around the internal metabolites in the system. However, in order for the
0039 %system to uptake or excrete metabolites, some metabolites have been
0040 %defined as "unconstrained". In order to simulate something, those
0041 %metabolites have to be removed from the model. The function simplifyModel
0042 %is a general-purpose function for making models smaller. This includes the
0043 %options such as grouping linear reactions and deleting reactions which
0044 %cannot carry flux. Here it is chosen to delete the exchange metabolites,
0045 %all reactions that are constrained to zero (mainly uptake of non-standard
0046 %carbon sources), and all reactions that cannot carry flux (mainly
0047 %reactions that were dependent on any of those non-standard carbons
0048 %sources).
0049 model=simplifyModel(model,true,false,true,true);
0050 
0051 %As can be seen the model now contains only 1305 reactions, 1037
0052 %metabolites and 1006 genes
0053 
0054 %As a first way of validating the model, calculate the theoretical yield of
0055 %carbon dioxide from glucose. The supplied model already allows for uptake
0056 %of phosphate, sulfate, NH3, O2 and the co-factor precursors thiamin and
0057 %pimelate. The setParam function is used for setting constraints,
0058 %reversibility and objective function coefficients.
0059 %Set the uptake of glucose to be no more than 1 mmol/gDW/h and no uptake of
0060 %ethanol.
0061 model=setParam(model,'ub',{'glcIN' 'etohIN'},[1 0]);
0062 
0063 %Set the objective for the simulation to maximize CO2 production
0064 model=setParam(model,'obj',{'co2OUT'},1);
0065 
0066 %The problem can now be solved using linear programming. The solveLP
0067 %function takes a model and solves the linear programming problem defined
0068 %by the constraints and the objective value coefficients.
0069 sol=solveLP(model);
0070 
0071 %If everything worked fine one should see a structure that contains the
0072 %fields .f which is the negative of the objective value, .stat which is 1
0073 %if the optimization terminated successfully and .x which contains the
0074 %fluxes through each of the reactions in this particular solution of the
0075 %problem. sol.x contains 1305 values which makes it rather difficult to
0076 %interpret. A first approach is to look at only the reactions that
0077 %transports metabolites in to or out from the system. The printFluxes
0078 %function prints (parts of) the flux distribution to the console.
0079 printFluxes(model, sol.x, true, 10^-7); %true means only print exchange fluxes
0080 
0081 %One can see that the system took up one unit of glucose and 6 units of O2
0082 %while producing 6 units of CO2 and 6 units of H2O. To get more detailed
0083 %insight how this happens, one can use the same function to plot all fluxes
0084 %(above 10^-7).
0085 printFluxes(model, sol.x, false, 10^-7);
0086 
0087 %The results show many reactions that have -1000 or 1000 flux. This is
0088 %because there are loops in the solution. In order to clean up the solution
0089 %one can minimize the sum of all the fluxes. This is done by setting the
0090 %third flag to solveLP to true (take a look at solveLP, there are other
0091 %options as well).
0092 sol=solveLP(model,1);
0093 printFluxes(model, sol.x, false, 10^-7);
0094 
0095 %Now there are much fewer reactions that are active. Some of them will
0096 %probably appear to be a little "weird" but that is because the ATP and
0097 %NAD(P)H produced has to go somewhere in order to be mass balanced.
0098 
0099 %If one wants to study the growth instead, change the objective to be the
0100 %production of biomass instead of CO2
0101 model=setParam(model,'obj',{'bmOUT'},1);
0102 sol=solveLP(model,1);
0103 printFluxes(model, sol.x, true, 10^-7);
0104 
0105 %The results show that the growth rate is 0.084/h and that the system now
0106 %also requires sulfate, phosphate, NH3, thiamin and pimelate. Compare this
0107 %to the growth on ethanol instead of glucose. Use three times the molar
0108 %flux of ethanol since it contains 2 carbons rather than 6. 'eq' means
0109 %'equal to' and sets the lower and upper bound to the same value.
0110 modelETH=setParam(model,'eq',{'glcIN' 'etohIN'},[0 3]);
0111 solETH=solveLP(modelETH,1);
0112 printFluxes(modelETH, solETH.x, true, 10^-7);
0113 
0114 %Investigate how metabolism changes between the two carbon sources.
0115 %followChanged takes two flux distributions and lets the user select which
0116 %reactions should be printed. Here the reactions are shown that differ with
0117 %more than 50%, has a flux higher than 0.5 mmol/gDW/h and an absolute
0118 %difference higher than 0.5 mmol/gDW/h.
0119 followChanged(modelETH,sol.x,solETH.x, 50, 0.5, 0.5);
0120 
0121 %There are 65 such reactions. By studying them one can start to get an idea
0122 %about where the key changes occur. Visualization can help a lot in this
0123 %regard. One can investigate how ATP metabolism changes by running the
0124 %following command:
0125 followChanged(modelETH,sol.x,solETH.x, 30, 0.4, 0.4,{'ATP'});
0126 
0127 %See that on glucose ATP is generated in glycolysis but on ethanol it seems
0128 %to have to do with acetate and so on. This allows the user to look further
0129 %and further down until one understands the underlying flux redistributions
0130 %that give rise to different phenotypes.
0131 
0132 %The fluxes can also be visualized on a metabolic map. Green corresponds to
0133 %reactions which are more used for growth on glucose, and red are reactions
0134 %which are more used for growth on ethanol. Open the "GLCvsETH.pdf" file to
0135 %be able to zoom in on individual reactions.
0136 load 'pcPathway.mat' pathway;
0137 drawMap('Glucose vs ethanol',pathway,model,sol.x,solETH.x,modelETH,'GLCvsETH.pdf',10^-5);
tutorial1

PURPOSE

SYNOPSIS

DESCRIPTION

CROSS-REFERENCE INFORMATION

SOURCE CODE
