runPhenotypePhasePlane

 runPhenotypePhasePlane
   Runs phenotype phase plane analysis and plots the results. The first
   plot is a 3D surface plot showing the phenotype phase plane, the other
   two plots show the shadow prices of the metabolites from the two
   control reactions, which define the phases. Modified from the COBRA
   phenotypePhasePlane function.

 Input:
   model           a model structure
   controlRxn1     reaction identifier of the first reaction to be plotted
   controlRxn2     reaction identifier of the second reaction to be plotted
   nPts            the number of points to plot in each dimension (optional,
                   default 50)
   range1          the range [from 0 to range1] of reaction 1 to plot
                   (optional, default 20)
   range2          the range [from 0 to range2] of reaction 2 to plot
                   (optional, default 20)

 Output:
   growthRates1    a matrix of maximum growth rates
   shadowPrices1   a matrix with shadow prices for reaction 1
   shadowPrices2   a matrix with shadow prices for reaction 2

 Modified from COBRA Toolbox phenotypePhasePlane.m

 Usage: [growthRates, shadowPrices1, shadowPrices2] = runPhenotypePhasePlane(model, controlRxn1, controlRxn2, nPts, range1, range2)

0001 function [growthRates, shadowPrices1, shadowPrices2] = runPhenotypePhasePlane(model, controlRxn1, controlRxn2, nPts, range1, range2)
0002 % runPhenotypePhasePlane
0003 %   Runs phenotype phase plane analysis and plots the results. The first
0004 %   plot is a 3D surface plot showing the phenotype phase plane, the other
0005 %   two plots show the shadow prices of the metabolites from the two
0006 %   control reactions, which define the phases. Modified from the COBRA
0007 %   phenotypePhasePlane function.
0008 %
0009 % Input:
0010 %   model           a model structure
0011 %   controlRxn1     reaction identifier of the first reaction to be plotted
0012 %   controlRxn2     reaction identifier of the second reaction to be plotted
0013 %   nPts            the number of points to plot in each dimension (optional,
0014 %                   default 50)
0015 %   range1          the range [from 0 to range1] of reaction 1 to plot
0016 %                   (optional, default 20)
0017 %   range2          the range [from 0 to range2] of reaction 2 to plot
0018 %                   (optional, default 20)
0019 %
0020 % Output:
0021 %   growthRates1    a matrix of maximum growth rates
0022 %   shadowPrices1   a matrix with shadow prices for reaction 1
0023 %   shadowPrices2   a matrix with shadow prices for reaction 2
0024 %
0025 % Modified from COBRA Toolbox phenotypePhasePlane.m
0026 %
0027 % Usage: [growthRates, shadowPrices1, shadowPrices2] = runPhenotypePhasePlane(model, controlRxn1, controlRxn2, nPts, range1, range2)
0028 close all force % Close all existing figure windows (if open)
0029 if nargin < 4
0030     nPts = 50;
0031 end
0032 if nargin < 5
0033     range1 = 20;
0034 end
0035 if nargin < 6
0036     range2 = 20;
0037 end
0038 
0039 rxnID1 = getIndexes(model,controlRxn1,'rxns',true);
0040 metID1 = find(model.S(:,rxnID1));
0041 rxnID2 = getIndexes(model,controlRxn2,'rxns',true);
0042 metID2 = find(model.S(:,rxnID2));
0043 
0044 % Create empty vectors for the results
0045 ind1 = linspace(0,range1,nPts);
0046 ind2 = linspace(0,range2,nPts);
0047 growthRates = zeros(nPts);
0048 shadowPrices1 = zeros(nPts);
0049 shadowPrices2 = zeros(nPts);
0050 [~,hsSol] = solveLP(model);
0051 % Calulate points
0052 PB = ProgressBar2(nPts,'Running PhPP analysis','cli');
0053 for i = 1:nPts %ind1
0054     for j = 1:nPts %ind2
0055         model1 = setParam(model,'eq',controlRxn1,-1*ind1(i));
0056         model1 = setParam(model1,'eq',controlRxn2,-1*ind2(j));
0057         [fbasol,hsSol] = solveLP(model1,0,[],hsSol);
0058         try
0059             growthRates(j,i) = fbasol.x(logical(model1.c));
0060             shadowPrices1(j,i) = fbasol.sPrice(metID1);
0061             shadowPrices2(j,i) = fbasol.sPrice(metID2);   
0062         end
0063     end
0064     count(PB);
0065 end
0066 fprintf('\b\b\b\b\b\b\b\b\b\b\b\b\bCOMPLETE\n');
0067 
0068 % Plot the points
0069 figure(2);
0070 pcolor(ind1,ind2,shadowPrices1);
0071 xlabel(strrep(strcat(controlRxn1,' (mmol/g DW-hr)'),'_','\_')), ylabel(strrep(strcat(controlRxn2,' (mmol/g DW-hr)'),'_','\_')), zlabel('growth rate (1/hr)');
0072 title(['Shadow price ' strrep(model.mets{metID1},'_','-')]);
0073 colorbar();
0074 xticklabels(sprintfc('%d', -xticks))
0075 yticklabels(sprintfc('%d', -yticks))
0076 figure(3);
0077 pcolor(ind1,ind2,shadowPrices2);
0078 xlabel(strrep(strcat(controlRxn1,' (mmol/g DW-hr)'),'_','\_')), ylabel(strrep(strcat(controlRxn2,' (mmol/g DW-hr)'),'_','\_')), zlabel('growth rate (1/hr)');
0079 title(['Shadow price ' strrep(model.mets{metID2},'_','-')]);
0080 colorbar();
0081 xticklabels(sprintfc('%d', -xticks))
0082 yticklabels(sprintfc('%d', -yticks))
0083 figure(1);
0084 surfl(ind1,ind2,growthRates);
0085 xlabel(strrep(strcat(controlRxn1,' (mmol/g DW-hr)'),'_','\_')), ylabel(strrep(strcat(controlRxn2,' (mmol/g DW-hr)'),'_','\_')), zlabel('growth rate (1/hr)');
0086 colormap(hsv(17));
0087 xticklabels(sprintfc('%d', -xticks))
0088 yticklabels(sprintfc('%d', -yticks))
0089 end