phenotypePhasePlane

PURPOSE ^

phenotypePhasePlane Plots three phenotype phase planes for two reactions. The first plot is

SYNOPSIS ^

function [growthRates,shadowPrices1,shadowPrices2] = phenotypePhasePlane(model,controlRxn1,controlRxn2,nPts,range1,range2)

DESCRIPTION ^

phenotypePhasePlane Plots three phenotype phase planes for two reactions.  The first plot is 
 a double robustness analysis, a kind of 3D surface plot.  The second
 two plots show the shadow prices of the metabolites from the two control 
 reactions, which define the phases.  Use the COLORMAP and SHADING 
 functions to change the looks of the plots.

 [growthRates,shadowPrices1,shadowPrices2] = phenotypePhasePlane(model,controlRxn1,controlRxn2,nPts,range1,range2)

INPUTS
 model             COBRA model structure
 controlRxn1       the first reaction to be plotted
 controlRxn2       the second reaction to be plotted

OPTIONAL INPUTS
 nPts              the number of points to plot in each dimension
                   (Default = 50)
 range1            the range of reaction 1 to plot
                   (Default = 20)
 range2            the range of reaction 2 to plot
                   (Default = 20)

OUTPUTS
 growthRates1      a matrix of maximum growth rates
 shadowPrices1     a matrix of rxn 1 shadow prices
 shadowPrices2     a matrix of rxn 2 shadow prices

 Jeff Orth 6/26/08

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SOURCE CODE ^

0001 function [growthRates,shadowPrices1,shadowPrices2] = phenotypePhasePlane(model,controlRxn1,controlRxn2,nPts,range1,range2)
0002 %phenotypePhasePlane Plots three phenotype phase planes for two reactions.  The first plot is
0003 % a double robustness analysis, a kind of 3D surface plot.  The second
0004 % two plots show the shadow prices of the metabolites from the two control
0005 % reactions, which define the phases.  Use the COLORMAP and SHADING
0006 % functions to change the looks of the plots.
0007 %
0008 % [growthRates,shadowPrices1,shadowPrices2] = phenotypePhasePlane(model,controlRxn1,controlRxn2,nPts,range1,range2)
0009 %
0010 %INPUTS
0011 % model             COBRA model structure
0012 % controlRxn1       the first reaction to be plotted
0013 % controlRxn2       the second reaction to be plotted
0014 %
0015 %OPTIONAL INPUTS
0016 % nPts              the number of points to plot in each dimension
0017 %                   (Default = 50)
0018 % range1            the range of reaction 1 to plot
0019 %                   (Default = 20)
0020 % range2            the range of reaction 2 to plot
0021 %                   (Default = 20)
0022 %
0023 %OUTPUTS
0024 % growthRates1      a matrix of maximum growth rates
0025 % shadowPrices1     a matrix of rxn 1 shadow prices
0026 % shadowPrices2     a matrix of rxn 2 shadow prices
0027 %
0028 % Jeff Orth 6/26/08
0029 
0030 if (nargin < 4)
0031     nPts = 50;
0032 end
0033 if (nargin < 5)
0034     range1 = 20;
0035 end
0036 if (nargin < 6)
0037     range2 = 20;
0038 end
0039 
0040 % find rxn and met ID numbers to get shadow prices and reduced costs
0041 rxnID1 = findRxnIDs(model,controlRxn1);
0042 metID1 = find(model.S(:,rxnID1));
0043 rxnID2 = findRxnIDs(model,controlRxn2);
0044 metID2 = find(model.S(:,rxnID2));
0045 
0046 % create empty vectors for the results
0047 ind1 = linspace(0,range1,nPts);
0048 ind2 = linspace(0,range2,nPts);
0049 growthRates = zeros(nPts);
0050 shadowPrices1 = zeros(nPts);
0051 shadowPrices2 = zeros(nPts);
0052 
0053 % calulate points
0054 h = waitbar(0,'generating PhPP');
0055 global CBT_LP_PARAMS  % save the state of the primal only flag.
0056 if isfield( CBT_LP_PARAMS, 'primalOnly')
0057     primalOnlySave = CBT_LP_PARAMS.primalOnly;
0058 end
0059 changeCobraSolverParams('LP', 'primalOnly', false);
0060 for i = 1:nPts %ind1
0061     for j = 1:nPts %ind2
0062         waitbar((nPts*(i-1)+j)/(nPts^2),h);
0063         model1 = changeRxnBounds(model,controlRxn1,-1*ind1(i),'b');
0064         model1 = changeRxnBounds(model1,controlRxn2,-1*ind2(j),'b');
0065                 
0066         fbasol = optimizeCbModel(model1,'max');
0067         growthRates(j,i) = fbasol.f;
0068         try % calculate shadow prices
0069             shadowPrices1(j,i) = fbasol.y(metID1(1));
0070             shadowPrices2(j,i) = fbasol.y(metID2(1));
0071         end
0072     end
0073 end
0074 if ( regexp( version, 'R20') )
0075         close(h);
0076 end
0077 if exist('primalOnlySave', 'var')
0078     changeCobraSolverParams('LP', 'primalOnly', primalOnlySave);
0079 else
0080     changeCobraSolverParams('LP', 'primalOnly', true);
0081 end
0082 
0083 % plot the points
0084 figure(2);
0085 pcolor(ind1,ind2,shadowPrices1);
0086 xlabel(strrep(strcat(controlRxn1,' (mmol/g DW-hr)'),'_','\_')), ylabel(strrep(strcat(controlRxn2,' (mmol/g DW-hr)'),'_','\_')), zlabel('growth rate (1/hr)');
0087 figure(3);
0088 pcolor(ind1,ind2,shadowPrices2);
0089 xlabel(strrep(strcat(controlRxn1,' (mmol/g DW-hr)'),'_','\_')), ylabel(strrep(strcat(controlRxn2,' (mmol/g DW-hr)'),'_','\_')), zlabel('growth rate (1/hr)');
0090 figure(1);
0091 surfl(ind1,ind2,growthRates);
0092 xlabel(strrep(strcat(controlRxn1,' (mmol/g DW-hr)'),'_','\_')), ylabel(strrep(strcat(controlRxn2,' (mmol/g DW-hr)'),'_','\_')), zlabel('growth rate (1/hr)');
0093 
0094 
0095 
0096

Generated on Thu 21-Jun-2012 15:39:23 by m2html © 2003