reduceModel

PURPOSE ^

reduceModel Removes from the model all of the reactions that are never used (max and

SYNOPSIS ^

function [modelRed,hasFlux,maxes,mins] = reduceModel(model,tol,irrevFlag,verbFlag,negFluxAllowedFlag,checkConsistencyFlag,changeBoundsFlag)

DESCRIPTION ^

reduceModel Removes from the model all of the reactions that are never used (max and
 min are < tol). Finds the minimal bounds for the flux through each reaction.
 Also returns the results for flux variability analysis (maxes, mins).

 [modelRed,hasFlux,maxes,mins] = reduceModel(model,tol,irrevFlag,verbFlag,negFluxAllowedFlag,checkConsistencyFlag,changeBoundsFlag)

INPUT
 model                 COBRA model structure

OPTIONAL INPUTS
 tol                   Tolerance for non-zero bounds - bounds smaller in absolute
                       value than this value will be set to zero (Default = 1e-6)
 irrevFlag             Determines if the models should be treated using 
                       the irreversible form. (Default = false)
 verbFlag              Verbose output (Default = false)
 negFluxAllowedFlag    Allow negative fluxes through irrev reactions
                       (Default = false)
 checkConsistencyFlag  Do a consistency check of the optimal solution
                       (Default = true)
 changeBoundsFlag      Change upper/lower bounds to the minimal bounds
                       (Default = true)

OUTPUTS
 modelRed              Reduced model
 hasFlux               The indexes of the reactions that are not blocked 
                       in the model
 maxes                 Maximum fluxes
 mins                  Minimum fluxes

 Gregory Hannum and Markus Herrgard 7/20/05

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SUBFUNCTIONS ^

SOURCE CODE ^

0001 function [modelRed,hasFlux,maxes,mins] = reduceModel(model,tol,irrevFlag,verbFlag,negFluxAllowedFlag,checkConsistencyFlag,changeBoundsFlag)
0002 %reduceModel Removes from the model all of the reactions that are never used (max and
0003 % min are < tol). Finds the minimal bounds for the flux through each reaction.
0004 % Also returns the results for flux variability analysis (maxes, mins).
0005 %
0006 % [modelRed,hasFlux,maxes,mins] = reduceModel(model,tol,irrevFlag,verbFlag,negFluxAllowedFlag,checkConsistencyFlag,changeBoundsFlag)
0007 %
0008 %INPUT
0009 % model                 COBRA model structure
0010 %
0011 %OPTIONAL INPUTS
0012 % tol                   Tolerance for non-zero bounds - bounds smaller in absolute
0013 %                       value than this value will be set to zero (Default = 1e-6)
0014 % irrevFlag             Determines if the models should be treated using
0015 %                       the irreversible form. (Default = false)
0016 % verbFlag              Verbose output (Default = false)
0017 % negFluxAllowedFlag    Allow negative fluxes through irrev reactions
0018 %                       (Default = false)
0019 % checkConsistencyFlag  Do a consistency check of the optimal solution
0020 %                       (Default = true)
0021 % changeBoundsFlag      Change upper/lower bounds to the minimal bounds
0022 %                       (Default = true)
0023 %
0024 %OUTPUTS
0025 % modelRed              Reduced model
0026 % hasFlux               The indexes of the reactions that are not blocked
0027 %                       in the model
0028 % maxes                 Maximum fluxes
0029 % mins                  Minimum fluxes
0030 %
0031 % Gregory Hannum and Markus Herrgard 7/20/05
0032 
0033 % Sets the tolerance for zero flux determination
0034 if nargin < 2
0035     global CBT_LP_PARAMS
0036     if (exist('CBT_LP_PARAMS', 'var'))
0037         if isfield(CBT_LP_PARAMS, 'objTol')
0038             tol = CBT_LP_PARAMS.objTol;
0039         else
0040             tol = 1e-6
0041         end
0042     else
0043         tol = 1e-6;
0044     end
0045 end
0046 
0047 % Sets the irrevFlag to default
0048 if nargin < 3
0049     irrevFlag = false;
0050 end
0051 
0052 % Print out more stuff
0053 if nargin < 4
0054     verbFlag = false;
0055 end
0056 
0057 % Allow negative irreversible fluxes (default: reverse the reaction
0058 % direction)
0059 if (nargin < 5)
0060     negFluxAllowedFlag = false;
0061 end
0062 
0063 % Check if the reduced model produces consistent results
0064 if (nargin < 6)
0065     checkConsistencyFlag = true;
0066 end
0067 
0068 % Change to minimal bounds
0069 if (nargin < 7)
0070     changeBoundsFlag = true;
0071 end
0072 
0073 %declare some variables
0074 maxes = [];
0075 mins = [];
0076 %modelRed = model;
0077 [nMets,nRxns]= size(model.S);
0078 
0079 %obtain maxes and mins for the fluxes
0080 rxnID = 1;
0081 h = waitbar(0,'Model reduction in progress ...');
0082 while rxnID <= nRxns
0083     if mod(rxnID,10) == 0
0084         waitbar(rxnID/nRxns,h);
0085     end
0086     rxnName = model.rxns{rxnID};
0087     if (verbFlag)
0088         fprintf('%s\t',rxnName);
0089     end
0090 
0091     % Set the objective function to the current reactiom
0092     tempModel = changeObjective(model,rxnName);
0093 
0094     if (irrevFlag && model.rev(rxnID))
0095         % Make the forward reaction reversible temporarily
0096         tempModel.lb(rxnID) = -tempModel.ub(rxnID+1);
0097         % Disable the reverse reaction
0098         tempModel.ub(rxnID+1) = 0;
0099     end
0100 
0101     %solve for the minimum and maximum for the current reaction
0102     sol = optimizeCbModel(tempModel,'max',[]);
0103     if (sol.stat > 0)
0104         maxBound = sol.f;
0105     else
0106         maxBound = model.ub(rxnID);
0107     end
0108     sol = optimizeCbModel(tempModel,'min',[]);
0109     if (sol.stat > 0)
0110         minBound = sol.f;
0111     else
0112         minBound = model.lb(rxnID);
0113     end
0114 
0115     %eliminate very small boundaries and set predetermined reversible boundaries
0116     if abs(maxBound) < tol
0117         maxBound = 0;
0118     end
0119 
0120     % Ignore negative lower bounds for irrev reactions
0121     if abs(minBound) < tol || (minBound < 0 && ~model.rev(rxnID))
0122         minBound = 0;
0123     end
0124 
0125     %set the new appropriate bounds
0126     if (irrevFlag && model.rev(rxnID))
0127         if minBound < 0 && maxBound < 0 % Negative flux
0128             mins(rxnID) = 0;
0129             mins(rxnID+1) = -maxBound;
0130             maxes(rxnID) = 0;
0131             maxes(rxnID+1) = -minBound;
0132         elseif minBound < 0 && maxBound >= 0 % Reversible flux
0133             mins(rxnID:rxnID+1) = 0;
0134             maxes(rxnID) = maxBound;
0135             maxes(rxnID+1) = -minBound;
0136         elseif minBound >= 0 && maxBound >= 0 % Positive flux
0137             mins(rxnID) = minBound;
0138             mins(rxnID+1) = 0;
0139             maxes(rxnID) = maxBound;
0140             maxes(rxnID+1) = 0;
0141         end
0142 
0143         if (verbFlag)
0144             fprintf('%g\t%g\n',mins(rxnID),maxes(rxnID));
0145             fprintf('%s\t',model.rxns{rxnID+1});
0146             fprintf('%g\t%g\n',mins(rxnID+1),maxes(rxnID+1));
0147         end
0148         % Jump over the reverse direction
0149         rxnID = rxnID + 1;
0150     else
0151         maxes(rxnID)=maxBound;
0152         mins(rxnID)=minBound;
0153         if (verbFlag)
0154             fprintf('%g\t%g\n',minBound,maxBound);
0155         end
0156     end
0157 
0158     rxnID = rxnID + 1;
0159 end
0160 if ( regexp( version, 'R20') )
0161         close(h);
0162 end
0163 
0164 if (verbFlag)
0165     fprintf('\n');
0166 end
0167 
0168 % Create a list of flux indexes that have non-zero flux (hasFlux)
0169 hasFluxSel = (abs(maxes) > tol | abs(mins) > tol);
0170 hasFlux = find(hasFluxSel);
0171 hasFlux = columnVector(hasFlux);
0172 
0173 % Remove reactions that are blocked
0174 modelRed = removeRxns(model,model.rxns(~hasFluxSel),irrevFlag,true);
0175 
0176 % Update bounds
0177 if (changeBoundsFlag)
0178     modelRed.lb = columnVector(mins(hasFlux));
0179     modelRed.ub = columnVector(maxes(hasFlux));
0180     selInconsistentBounds = (modelRed.ub < modelRed.lb);
0181     modelRed.ub(selInconsistentBounds) = modelRed.lb(selInconsistentBounds);
0182     %update the reversible list with new bounds
0183     nRxnsNew = size(modelRed.S,2);
0184     for rxnID = 1:nRxnsNew
0185         if (~irrevFlag)
0186             if (modelRed.lb(rxnID) >= 0)
0187                 % Only runs in positive direction
0188                 modelRed.rev(rxnID) = false;
0189             end
0190             if (modelRed.ub(rxnID) <= 0)
0191 
0192                 % Only runs in negative direction -> reverse the reaction
0193                 modelRed.rev(rxnID) = false;
0194                 if (~negFluxAllowedFlag)
0195                     ubTmp = modelRed.ub(rxnID);
0196                     lbTmp = modelRed.lb(rxnID);
0197                     modelRed.S(:,rxnID) = -modelRed.S(:,rxnID);
0198                     modelRed.ub(rxnID) = -lbTmp;
0199                     modelRed.lb(rxnID) = -ubTmp;
0200                     modelRed.c(rxnID) = -modelRed.c(rxnID);
0201                     modelRed.rxns{rxnID} = [modelRed.rxns{rxnID} '_r'];
0202                 end
0203             end
0204         end
0205     end
0206 
0207     if (checkConsistencyFlag)
0208         fprintf('Perform model consistency check\n');
0209         modelOK = checkConsistency(model,modelRed,tol);
0210         if (~modelOK)
0211             modelRed = expandBounds(model,modelRed,tol);
0212         end
0213     end
0214 else
0215     if (checkConsistencyFlag)
0216         fprintf('Perform model consistency check\n');
0217         modelOK = checkConsistency(model,modelRed,tol);
0218     end
0219 end
0220 
0221 %%
0222 function modelRed = expandBounds(model,modelRed,tol)
0223 % Expand bounds to achieve the desired objective value
0224 %
0225 % modelRed = expandBounds(model,modelRed,tol)
0226 %
0227 
0228 modelOK = false;
0229 cushion = tol;
0230 tempModel = modelRed;
0231 while (~modelOK)
0232     narrowInd = find(modelRed.ub-modelRed.lb < cushion & modelRed.ub ~= modelRed.lb);
0233     tempModel.lb(narrowInd) = tempModel.lb(narrowInd) - cushion;
0234     narrowIrrevInd =intersect(narrowInd,find(~tempModel.rev));
0235     tempModel.lb(narrowIrrevInd) = max(tempModel.lb(narrowIrrevInd),0);
0236     tempModel.ub(narrowInd) = tempModel.ub(narrowInd) + cushion;
0237     modelRed.lb(narrowInd) = tempModel.lb(narrowInd);
0238     modelRed.ub(narrowInd) = tempModel.ub(narrowInd);
0239     cushion = cushion*2;
0240     modelOK = checkConsistency(model,tempModel,tol);
0241 end
0242 
0243 %%
0244 function modelOK = checkConsistency(model,modelRed,tol)
0245 %
0246 % modelOK = checkConsistency(model,modelRed,tol)
0247 %
0248 
0249 if (sum(model.c ~= 0) > 0)
0250 
0251     % Original model
0252     solOrigMax = optimizeCbModel(model,'max',[]);
0253     solOrigMin = optimizeCbModel(model,'min',[]);
0254 
0255     % Reduced model
0256     solRedMax = optimizeCbModel(modelRed,'max',[]);
0257     solRedMin = optimizeCbModel(modelRed,'min',[]);
0258 
0259     diffMax = abs(solRedMax.f - solOrigMax.f);
0260     diffMin = abs(solRedMin.f - solOrigMin.f);
0261 
0262     if (diffMax > tol || diffMin > tol)
0263         fprintf('Inconsistent objective values %g %g %g %g\n',solOrigMax.f,solRedMax.f,solOrigMin.f,solRedMin.f);
0264         modelOK = false;
0265     else
0266         fprintf('Model is consistent\n');
0267         modelOK = true;
0268     end
0269 
0270 else
0271     modelOK = true;
0272 end

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