0001 function [modelRed,hasFlux,maxes,mins] = reduceModel(model,tol,irrevFlag,verbFlag,negFluxAllowedFlag,checkConsistencyFlag,changeBoundsFlag)
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
0031
0032
0033
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
0048 if nargin < 3
0049 irrevFlag = false;
0050 end
0051
0052
0053 if nargin < 4
0054 verbFlag = false;
0055 end
0056
0057
0058
0059 if (nargin < 5)
0060 negFluxAllowedFlag = false;
0061 end
0062
0063
0064 if (nargin < 6)
0065 checkConsistencyFlag = true;
0066 end
0067
0068
0069 if (nargin < 7)
0070 changeBoundsFlag = true;
0071 end
0072
0073
0074 maxes = [];
0075 mins = [];
0076
0077 [nMets,nRxns]= size(model.S);
0078
0079
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
0092 tempModel = changeObjective(model,rxnName);
0093
0094 if (irrevFlag && model.rev(rxnID))
0095
0096 tempModel.lb(rxnID) = -tempModel.ub(rxnID+1);
0097
0098 tempModel.ub(rxnID+1) = 0;
0099 end
0100
0101
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
0116 if abs(maxBound) < tol
0117 maxBound = 0;
0118 end
0119
0120
0121 if abs(minBound) < tol || (minBound < 0 && ~model.rev(rxnID))
0122 minBound = 0;
0123 end
0124
0125
0126 if (irrevFlag && model.rev(rxnID))
0127 if minBound < 0 && maxBound < 0
0128 mins(rxnID) = 0;
0129 mins(rxnID+1) = -maxBound;
0130 maxes(rxnID) = 0;
0131 maxes(rxnID+1) = -minBound;
0132 elseif minBound < 0 && maxBound >= 0
0133 mins(rxnID:rxnID+1) = 0;
0134 maxes(rxnID) = maxBound;
0135 maxes(rxnID+1) = -minBound;
0136 elseif minBound >= 0 && maxBound >= 0
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
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
0169 hasFluxSel = (abs(maxes) > tol | abs(mins) > tol);
0170 hasFlux = find(hasFluxSel);
0171 hasFlux = columnVector(hasFlux);
0172
0173
0174 modelRed = removeRxns(model,model.rxns(~hasFluxSel),irrevFlag,true);
0175
0176
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
0183 nRxnsNew = size(modelRed.S,2);
0184 for rxnID = 1:nRxnsNew
0185 if (~irrevFlag)
0186 if (modelRed.lb(rxnID) >= 0)
0187
0188 modelRed.rev(rxnID) = false;
0189 end
0190 if (modelRed.ub(rxnID) <= 0)
0191
0192
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
0224
0225
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
0247
0248
0249 if (sum(model.c ~= 0) > 0)
0250
0251
0252 solOrigMax = optimizeCbModel(model,'max',[]);
0253 solOrigMin = optimizeCbModel(model,'min',[]);
0254
0255
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