0001 function [controlFlux1, controlFlux2, objFlux] = doubleRobustnessAnalysis(model, controlRxn1, controlRxn2, nPoints, plotResFlag, objRxn,objType)
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 if (nargin < 4)
0029 nPoints = 20;
0030 end
0031 if (nargin < 5)
0032 plotResFlag = true;
0033 end
0034 if (nargin > 6)
0035 baseModel = changeObjective(model,objRxn);
0036 else
0037 baseModel = model;
0038 end
0039 if (nargin <7)
0040 objType = 'max';
0041 end
0042
0043 if (findRxnIDs(model,controlRxn1))
0044 tmpModel = changeObjective(model,controlRxn1);
0045 solMin1 = optimizeCbModel(tmpModel,'min');
0046 solMax1 = optimizeCbModel(tmpModel,'max');
0047 else
0048 error('Control reaction 1 does not exist!');
0049 end
0050 if (findRxnIDs(model,controlRxn2))
0051 tmpModel = changeObjective(model,controlRxn2);
0052 solMin2 = optimizeCbModel(tmpModel,'min');
0053 solMax2 = optimizeCbModel(tmpModel,'max');
0054 else
0055 error('Control reaction 2 does not exist!');
0056 end
0057
0058 objFlux = [];
0059 controlFlux1 = linspace(solMin1.f,solMax1.f,nPoints)';
0060 controlFlux2 = linspace(solMin2.f,solMax2.f,nPoints)';
0061
0062 h = waitbar(0,'Double robustness analysis in progress ...');
0063 for i=1:nPoints
0064 for j = 1:nPoints
0065 waitbar(((i-1)*nPoints+j)/nPoints^2,h);
0066 modelControlled = changeRxnBounds(baseModel,controlRxn1,controlFlux1(i),'b');
0067 modelControlled = changeRxnBounds(modelControlled,controlRxn2,controlFlux2(j),'b');
0068 solControlled = optimizeCbModel(modelControlled,objType);
0069 objFlux(i,j) = solControlled.f;
0070 end
0071 end
0072 if ( regexp( version, 'R20') )
0073 close(h);
0074 end
0075
0076 if (plotResFlag)
0077 clf
0078 surf(controlFlux1,controlFlux2,objFlux);
0079
0080 xlabel(strrep(controlRxn1,'_','-'));
0081 ylabel(strrep(controlRxn2,'_','-'));
0082 axis tight
0083 end
0084
0085
0086
0087
0088
0089