doublerobustnessanalysis.m

分析微生物化学计量模型中代谢通量双鲁棒性的matlab函数

function [controlFlux1, controlFlux2, objFlux] = doubleRobustnessAnalysis(model, controlRxn1, controlRxn2, nPoints, plotResFlag, objRxn,objType)
%doubleRobustnessAnalysis Performs robustness analysis for a pair of reactions of
% interest and an objective of interest
%
% [controlFlux1, controlFlux2, objFlux] = doubleRobustnessAnalysis(model, controlRxn1, controlRxn2, nPoints, plotResFlag, objRxn, objType)
%
% model         Model structure
% controlRxn    Reaction of interest whose value is to be controlled
% nPoints       Number of flux values per dimension (opt, default = 20)
% plotResFlag   Plot results (opt, default true)
% objRxn        Objective reaction to be maximized (opt, default whatever
% is defined in model)
% objType       Maximize ('max') or minimize ('min') objective (opt,
% default 'max')
%
% controlFlux   Flux values within the range of the maximum and minimum for
%               reaction of interest  
% objFlux       Optimal values of objective reaction at each control
% reaction flux value
%
% Monica Mo and Markus Herrgard 8/20/07

if (nargin < 4)
    nPoints = 20;
end
if (nargin < 5)
    plotResFlag = true;
end
if (nargin > 6)
    baseModel = changeObjective(model,objRxn);
else
    baseModel = model;
end
if (nargin <7)
    objType = 'max';
end

if (findRxnIDs(model,controlRxn1))
    tmpModel = changeObjective(model,controlRxn1);
    solMin1 = optimizeCbModel(tmpModel,'min');
    solMax1 = optimizeCbModel(tmpModel,'max');
else
    error('Control reaction 1 does not exist!');
end
if (findRxnIDs(model,controlRxn2))
    tmpModel = changeObjective(model,controlRxn2);
    solMin2 = optimizeCbModel(tmpModel,'min');
    solMax2 = optimizeCbModel(tmpModel,'max');
else
    error('Control reaction 2 does not exist!');
end

objFlux = [];
controlFlux1 = linspace(solMin1.f,solMax1.f,nPoints)';
controlFlux2 = linspace(solMin2.f,solMax2.f,nPoints)';

h = waitbar(0,'Double robustness analysis in progress ...');
for i=1:nPoints
    for j = 1:nPoints
        waitbar(((i-1)*nPoints+j)/nPoints^2,h);
        modelControlled = changeRxnBounds(baseModel,controlRxn1,controlFlux1(i),'b');
        modelControlled = changeRxnBounds(modelControlled,controlRxn2,controlFlux2(j),'b');
        solControlled = optimizeCbModel(modelControlled,objType);
        objFlux(i,j) = solControlled.f;
    end
end
close(h);

if (plotResFlag)
    clf
    surf(controlFlux1,controlFlux2,objFlux);
    %shading interp
    xlabel(strrep(controlRxn1,'_','-'));
    ylabel(strrep(controlRxn2,'_','-'));
    axis tight
end