0001 function [growthRates,shadowPrices1,shadowPrices2] = phenotypePhasePlane(model,controlRxn1,controlRxn2,nPts,range1,range2)
0002
0003
0004
0005
0006
0007
0008
0009
0010
0011
0012
0013
0014
0015
0016
0017
0018
0019
0020
0021
0022
0023
0024
0025
0026
0027
0028
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
0041 rxnID1 = findRxnIDs(model,controlRxn1);
0042 metID1 = find(model.S(:,rxnID1));
0043 rxnID2 = findRxnIDs(model,controlRxn2);
0044 metID2 = find(model.S(:,rxnID2));
0045
0046
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
0054 h = waitbar(0,'generating PhPP');
0055 global CBT_LP_PARAMS
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
0061 for j = 1:nPts
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
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
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