0001 function [minFlux,maxFlux,Vmin,Vmax] = fluxVariability(model,optPercentage,osenseStr,rxnNameList,verbFlag, allowLoops)
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
0035
0036
0037 if (nargin < 2)
0038 optPercentage = 100;
0039 end
0040 if (nargin < 3)
0041 osenseStr = 'max';
0042 end
0043 if (nargin < 4)
0044 rxnNameList = model.rxns;
0045 end
0046 if (nargin < 5)
0047 verbFlag = false;
0048 end
0049 if (nargin < 6)
0050 allowLoops = true;
0051 end
0052 if (isempty(optPercentage))
0053 optPercentage = 100;
0054 end
0055 if (isempty(osenseStr))
0056 osenseStr = 'max';
0057 end
0058 if (isempty(rxnNameList))
0059 rxnNameList = model.rxns;
0060 end
0061
0062
0063 global CBT_LP_PARAMS
0064 if (exist('CBT_LP_PARAMS', 'var'))
0065 if isfield(CBT_LP_PARAMS, 'objTol')
0066 tol = CBT_LP_PARAMS.objTol;
0067 else
0068 tol = 1e-6;
0069 end
0070 if isfield(CBT_LP_PARAMS, 'minNorm')
0071 minNorm = CBT_LP_PARAMS.minNorm;
0072 else
0073 minNorm = 0;
0074 end
0075 else
0076 tol = 1e-6;
0077 minNorm = 0;
0078 end
0079
0080
0081 if (sum(model.c ~= 0) > 0)
0082 hasObjective = true;
0083 optSol = optimizeCbModel(model,osenseStr, 0, allowLoops);
0084 if (optSol.stat > 0)
0085 objRxn = model.rxns(model.c~=0);
0086 if (strcmp(osenseStr,'max'))
0087 objValue = floor(optSol.f/tol)*tol*optPercentage/100;
0088 else
0089 objValue = ceil(optSol.f/tol)*tol*optPercentage/100;
0090 end
0091 else
0092 error('Infeasible problem - no optimal solution!');
0093 end
0094 else
0095 hasObjective = false;
0096 end
0097
0098 if (verbFlag == 1)
0099 h = waitbar(0,'Flux variability analysis in progress ...');
0100 end
0101 if (verbFlag > 1)
0102 fprintf('%4s\t%4s\t%10s\t%9s\t%9s\n','No','Perc','Name','Min','Max');
0103 end
0104
0105 if (~isfield(model,'b'))
0106 model.b = zeros(size(model.S,1),1);
0107 end
0108
0109 [nMets,nRxns] = size(model.S);
0110 rxnListFull = model.rxns;
0111 LPproblem.c = model.c;
0112 LPproblem.lb = model.lb;
0113 LPproblem.ub = model.ub;
0114 LPproblem.csense(1:nMets) = 'E';
0115 LPproblem.csense = LPproblem.csense';
0116 if hasObjective
0117 LPproblem.A = [model.S;columnVector(model.c)'];
0118 LPproblem.b = [model.b;objValue];
0119 if (strcmp(osenseStr,'max'))
0120 LPproblem.csense(end+1) = 'G';
0121 else
0122 LPproblem.csense(end+1) = 'L';
0123 end
0124 else
0125 LPproblem.A = model.S;
0126 LPproblem.b = model.b;
0127 end
0128
0129
0130
0131 LPproblem.osense = -1;
0132 tempSolution = solveCobraLP(LPproblem);
0133 LPproblem.basis = tempSolution.basis;
0134
0135
0136 maxFlux = zeros(length(rxnNameList), 1);
0137 minFlux = zeros(length(rxnNameList), 1);
0138
0139 if length(minNorm)> 1 || minNorm > 0
0140 Vmin=zeros(nRxns,nRxns);
0141 Vmax=zeros(nRxns,nRxns);
0142
0143
0144 allowLoops=1;
0145 else
0146 Vmin=[];
0147 Vmax=[];
0148 end
0149
0150 solutionPool = zeros(length(model.lb), 0);
0151 if ~exist('matlabpool') || (matlabpool('size') == 0)
0152 m = 0;
0153 for i = 1:length(rxnNameList)
0154 if mod(i,10) == 0, clear mex, end
0155 if (verbFlag == 1),fprintf('iteration %d. skipped %d\n', i, round(m));end
0156 LPproblem.c = zeros(nRxns,1);
0157 rxnBool=strcmp(rxnListFull,rxnNameList{i});
0158 LPproblem.c(rxnBool) = 1;
0159
0160 LPproblem.osense = -1;
0161 LPsolution = solveCobraLP(LPproblem);
0162
0163 maxFlux(i) = LPsolution.full(LPproblem.c~=0);
0164
0165
0166
0167 if length(minNorm)> 1 || minNorm > 0
0168 QPproblem=LPproblem;
0169 QPproblem.lb(LPproblem.c~=0)=maxFlux(i)-1e-12;
0170 QPproblem.ub(LPproblem.c~=0)=maxFlux(i)+1e12;
0171 QPproblem.c(:)=0;
0172
0173 if length(minNorm)==1
0174 minNorm=ones(nRxns,1)*minNorm;
0175 end
0176 QPproblem.F = spdiags(minNorm,0,nRxns,nRxns);
0177
0178 solution = solveCobraQP(QPproblem);
0179 if isempty(solution.full)
0180 pause(eps)
0181 end
0182 Vmax(:,rxnBool)=solution.full(1:nRxns,1);
0183 end
0184
0185 LPproblem.osense = 1;
0186 LPsolution = solveCobraLP(LPproblem);
0187
0188 minFlux(i) = LPsolution.full(LPproblem.c~=0);
0189
0190
0191
0192
0193
0194 if length(minNorm)> 1 || minNorm > 0
0195 QPproblem=LPproblem;
0196 QPproblem.lb(LPproblem.c~=0)=maxFlux(i)-1e-12;
0197 QPproblem.ub(LPproblem.c~=0)=maxFlux(i)+1e12;
0198 QPproblem.c(:)=0;
0199 QPproblem.F = spdiags(minNorm,0,nRxns,nRxns);
0200
0201 if length(minNorm)==1
0202 minNorm=ones(nRxns,1)*minNorm;
0203 end
0204 QPproblem.F = spdiags(minNorm,0,nRxns,nRxns);
0205
0206 solution = solveCobraQP(QPproblem);
0207 Vmin(:,rxnBool)=solution.full(1:nRxns,1);
0208 end
0209
0210
0211 if ~allowLoops
0212 if any( abs(LPproblem.c'*solutionPool - maxFlux(i)) < tol)
0213
0214 m = m+.5;
0215 else
0216 LPproblem.osense = -1;
0217 LPsolution = solveCobraMILP(addLoopLawConstraints(LPproblem, model));
0218 maxFlux(i) = LPsolution.obj/1000;
0219 end
0220 if any( abs(LPproblem.c'*solutionPool - minFlux(i)) < tol)
0221 m = m+.5;
0222
0223 else
0224 LPproblem.osense = 1;
0225 LPsolution = solveCobraMILP(addLoopLawConstraints(LPproblem, model));
0226 minFlux(i) = LPsolution.obj/1000;
0227 end
0228 end
0229 if (verbFlag == 1)
0230 waitbar(i/length(rxnNameList),h);
0231 end
0232 if (verbFlag > 1)
0233 fprintf('%4d\t%4.0f\t%10s\t%9.3f\t%9.3f\n',i,100*i/length(rxnNameList),rxnNameList{i},minFlux(i),maxFlux(i));
0234 end
0235 end
0236 else
0237 parfor i = 1:length(rxnNameList)
0238
0239
0240 c = zeros(nRxns,1);
0241 c(strcmp(rxnListFull,rxnNameList{i})) = 1000;
0242 if allowLoops
0243 LPsolution = solveCobraLP(struct(...
0244 'A', LPproblem.A,...
0245 'b', LPproblem.b,...
0246 'lb', LPproblem.lb,...
0247 'ub', LPproblem.ub,...
0248 'csense', LPproblem.csense,...
0249 'c',c,...
0250 'osense',-1, ...
0251 'basis', LPproblem.basis ...
0252 ));
0253
0254 maxFlux(i) = LPsolution.full(c~=0);
0255
0256 LPsolution = solveCobraLP(struct(...
0257 'A', LPproblem.A,...
0258 'b', LPproblem.b,...
0259 'lb', LPproblem.lb,...
0260 'ub', LPproblem.ub,...
0261 'csense', LPproblem.csense,...
0262 'c',c,...
0263 'osense',1, ...
0264 'basis', LPproblem.basis ...
0265 ));
0266 minFlux(i) = LPsolution.full(c~=0);
0267 else
0268 LPsolution = solveCobraMILP(addLoopLawConstraints(struct(...
0269 'A', LPproblem.A,...
0270 'b', LPproblem.b,...
0271 'lb', LPproblem.lb,...
0272 'ub', LPproblem.ub,...
0273 'csense', LPproblem.csense,...
0274 'c',c,...
0275 'osense',-1 ...
0276 ), model));
0277 maxFlux(i) = LPsolution.obj/1000;
0278
0279 LPsolution = solveCobraMILP(addLoopLawConstraints(struct(...
0280 'A', LPproblem.A,...
0281 'b', LPproblem.b,...
0282 'lb', LPproblem.lb,...
0283 'ub', LPproblem.ub,...
0284 'csense', LPproblem.csense,...
0285 'c',c,...
0286 'osense',1 ...
0287 ), model));
0288 minFlux(i) = LPsolution.obj/1000;
0289 end
0290 end
0291 end
0292
0293
0294 if (verbFlag == 1)
0295 if ( regexp( version, 'R20') )
0296 close(h);
0297 end
0298 end
0299
0300 maxFlux = columnVector(maxFlux);
0301 minFlux = columnVector(minFlux);