tutorial4_solutions

 tutorial4_solutions
   This script contains the solutions for Tutorial 4, see Tutorial 4 in
   "RAVEN tutorials.docx" for more details.

   NOTE: Many of these changes are easier to do in the Excel sheet. They
   are done here in code just to avoid having several model files.

0001 % tutorial4_solutions
0002 %   This script contains the solutions for Tutorial 4, see Tutorial 4 in
0003 %   "RAVEN tutorials.docx" for more details.
0004 %
0005 %   NOTE: Many of these changes are easier to do in the Excel sheet. They
0006 %   are done here in code just to avoid having several model files.
0007 
0008 %Import the Excel model
0009 model=importExcelModel('smallYeastBad.xlsx');
0010 
0011 %Close all uptake and maximize for production
0012 model=setParam(model,'eq',{'glcIN', 'o2IN'},[0 0]);
0013 model=setParam(model,'obj',{'acOUT' 'biomassOUT' 'co2OUT' 'ethOUT' 'glyOUT'},[1 1 1 1 1]);
0014 sol=solveLP(model);
0015 printFluxes(model,sol.x,true); %Nothing produced, good
0016 
0017 %Add some fake reactions
0018 rxns.rxns={'FREE_ATP';'FREE_NADH';'FREE_NADPH'};
0019 rxns.equations={'ATP <=> ADP + phosphate';'NAD(+) <=> NADH';'NADP(+) <=> NADPH'};
0020 model=addRxns(model,rxns,2,'c');
0021 sol=solveLP(model,1);
0022 
0023 %Lots of ethanol produced. Also plot the equations to make the error easier
0024 %to find
0025 printFluxes(model,sol.x,false,[],[],'%rxnID (%rxnName):%flux\n\t%eqn\n');
0026 
0027 %See that ADH1 should only produce one unit of ethanol. Change the reaction
0028 %equation
0029 model=changeRxns(model,'ADH1','acetaldehyde[c] + NADH[c] => ethanol[c] + NAD(+)[c]',3);
0030 sol=solveLP(model,1);
0031 printFluxes(model,sol.x,true); %Nothing produced, good
0032 
0033 %Add excretion of all metabolites
0034 model.b=[model.b inf(numel(model.b),1)];
0035 sol=solveLP(model,1);
0036 printFluxes(model,sol.x,false,10^-5,[],'%rxnID (%rxnName):\n\t%eqn\n\t%flux\n');
0037 
0038 %By looking at the reactions which were unbalanced and that were in the
0039 %flux list one can see that FBP should be changed to result in only one unit
0040 %of F6P
0041 model=changeRxns(model,'FBP','beta-D-fructofuranose 1,6-bisphosphate[c] => beta-D-fructofuranose 6-phosphate[c] + phosphate[c]',3);
0042 sol=solveLP(model,1);
0043 printFluxes(model,sol.x,false,10^-5,[],'%rxnID (%rxnName):\n\t%eqn\n\t%flux\n');
0044 
0045 %The same thing again and one should change PFK to only give one unit of
0046 %F16P
0047 model=changeRxns(model,'PFK','ATP[c] + beta-D-fructofuranose 6-phosphate[c] => ADP[c] + beta-D-fructofuranose 1,6-bisphosphate[c]',3);
0048 sol=solveLP(model,1);
0049 printFluxes(model,sol.x,false,10^-5,[],'%rxnID (%rxnName):\n\t%eqn\n\t%flux\n'); %Now it works
0050 
0051 %Set all uptakes and production to 0
0052 model=setParam(model,'eq',getExchangeRxns(model),0);
0053 
0054 %Since it is checked which metabolites could be consumed without
0055 %production, one can no longer have free production of all metabolites
0056 model.b=model.b(:,1);
0057 I=canConsume(model);
0058 disp(model.mets(I)); %These 12 metabolites can be consumed without any production
0059 
0060 %Allow all uptake
0061 model.b=[ones(numel(model.b),1)*-1000 model.b];
0062 
0063 %Pick CO2 and force uptake of it
0064 model=setParam(model,'eq',{'co2OUT'},-1); %Negative output means input
0065 sol=solveLP(model);
0066 printFluxes(model,sol.x,false,10^-5,[],'%rxnID (%rxnName):\n\t%eqn\n\t%flux\n'); %Now it works
0067 
0068 %See that PDC converts pyruvate (3 carbons) to acetaldehyde (2 carbons)
0069 %without any other products. If one googles, one may realize that CO2 is
0070 %missing. This would be simpler to change in the Excel file (or using
0071 %changeRxns), but one can change it here as an exercise. One therefore
0072 %needs to find the index of the reactions and the index of cytosolic CO2 in
0073 %order to change the reaction
0074 Irxn=ismember(model.rxns,'PDC');
0075 Imet=ismember(model.mets,'CO2_c');
0076 model.S(Imet,Irxn)=1; %The coefficient is 1.0
0077 
0078 %Display the new equation just to be sure
0079 constructEquations(model,Irxn)
0080 
0081 %The solution is now not feasible, meaning that it is no longer possible to
0082 %force uptake of CO2 without any output
0083 sol=solveLP(model);
0084 
0085 %***Second part of tutorial
0086 model=importExcelModel('smallYeastBad2.xlsx',true,false,true); %This has to be loaded with the setting to ignore error or it would find the error
0087 [reducedModel, deletedReactions, deletedMetabolites]=simplifyModel(model,false,false,false,true);
0088 disp(deletedReactions);
0089 disp(deletedMetabolites);
0090 
0091 %It turned out that G15L_c was spelled G15Lc in one reaction. The best
0092 %solution would be just to change it in the reaction list and remove the
0093 %duplicate metabolite, but one can do it here as an exercise. The indexes
0094 %of the two metabolites are needed.
0095 Igood=ismember(model.mets,'G15L_c');
0096 Ibad=ismember(model.mets,'G15Lc');
0097 
0098 %Get all reactions and the coefficients in which the wrong one participates
0099 %move them to be for the right one instead
0100 model.S(Igood,:)=model.S(Igood,:)+model.S(Ibad,:);
0101 
0102 %Delete the bad one
0103 model=removeMets(model,'G15Lc');
0104 [reducedModel, deletedReactions, deletedMetabolites]=simplifyModel(model,false,false,false,true);
0105 disp(deletedReactions);
0106 disp(deletedMetabolites);
0107 
0108 %The only difference was that there were 20 deleted metabolites instead of
0109 %21. Nothing too spectacular. Check production can tell the user what one
0110 %needs to connect.
0111 [notProducedMets, ~, neededForProductionMat, minToConnect]=checkProduction(model,true,model.comps,false);
0112 
0113 %In order to have production of all 54 metabolites one needs to enable
0114 %production of these 12. This small model does not include net synthesis of
0115 %co-factors, so one should concentrate on the other ones. Glycerone
0116 %phosphate allows for connection 18 others, so it seems like a good target.
0117 disp(minToConnect);
0118 
0119 %If one googles around a little bit, and knows the metabolism, one would
0120 %find that DHAP (dihydroxyacetone) and GLYP (glycerone phosphate) are
0121 %actually synonymes. Only use DHAP
0122 Igood=ismember(model.mets,'DHAP_c');
0123 Ibad=ismember(model.mets,'GLYP_c');
0124 
0125 %Get all reactions and the coefficients in which the wrong one participates
0126 %move them to be for the right one instead
0127 model.S(Igood,:)=model.S(Igood,:)+model.S(Ibad,:);
0128 
0129 %Delete the bad one
0130 model=removeMets(model,'GLYP_c');
0131 [reducedModel, deletedReactions, deletedMetabolites]=simplifyModel(model,false,false,false,true);
0132 disp(deletedReactions);
0133 disp(deletedMetabolites);
0134 [notProducedMets, ~, neededForProductionMat, minToConnect]=checkProduction(model,true,model.comps,false);
0135 disp(minToConnect);
0136 
0137 %Still quite a lot of gaps and no immediate way to fix it. One could try
0138 %including reactions from a reference network and see if that helps. Use
0139 %the small yeast model from Tutorial 3
0140 refModel=importExcelModel('smallYeast.xlsx');
0141 [newConnected, cannotConnect, addedRxns, newModel]=fillGaps(model,{refModel},false);
0142 disp(addedRxns);
0143 disp(newConnected);
0144 
0145 %By including the ALD6 reaction from the reference model it was possible to
0146 %connect 21 reactions
0147 [reducedModel, deletedReactions, deletedMetabolites]=simplifyModel(newModel,false,false,false,true);
0148 disp(deletedMetabolites);
0149 disp(deletedReactions);
0150 
0151 %All the model seems to be connected
0152 
0153 %All this stuff may be done in a more automated manner as well
0154 model=importExcelModel('smallYeastBad2.xlsx',true,false,true);
0155 gapReport(model,{refModel});