0001 function [gdlsSolution, bilevelMILPProblem, gdlsSolutionStructs] = GDLS(model, targetRxns, varargin)
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 MAXFLUX = 1000;
0036 MAXDUAL = 1000;
0037 EPS = 1e-4;
0038 gdlsSolutionStructs = [];
0039
0040 if nargin < 2
0041 error('Model and target reaction(s) must be specified')
0042 elseif mod((nargin-2),2)==0
0043 options.targetRxns = targetRxns;
0044 for i=1:2:(nargin-2)
0045 if ismember(varargin{i},{'nbhdsz','M','maxKO','selectedRxns','koType','koCost','minGrowth', ...
0046 'timeLimit','iterationLimit'})
0047 options.(varargin{i}) = varargin{i+1};
0048 else
0049 display(['Unknown option ' varargin{i}]);
0050 end
0051 end
0052 elseif isstruct(targetRxns)
0053 options = varargin{1};
0054 else
0055 error('Invalid number of entries')
0056 end
0057
0058
0059 if ~isfield(options,'koType'), options.koType = 'rxns'; end
0060 if ~isfield(options,'nbhdsz'), options.nbhdsz=1; end
0061 if ~isfield(options,'M'), options.M=1; end
0062 if ~isfield(options,'maxKO'), options.maxKO=50; end
0063 if ~isfield(options,'iterationLimit'), options.iterationLimit=70; end
0064 if ~isfield(options,'timeLimit'), options.timeLimit=252000 ; end
0065 if isfield(options,'targetRxns')
0066 selTargetRxns = logical(ismember(model.rxns,options.targetRxns));
0067 if ~any(selTargetRxns)
0068 error([options.targetRxns ' not found. Double check spelling.']);
0069 end
0070 else
0071 error('No target reaction specified')
0072 end
0073 if isfield(options,'selectedRxns')
0074 selSelectedRxns = logical(ismember(model.rxns,options.selectedRxns));
0075 else
0076 selSelectedRxns = true(length(model.rxns),1);
0077 end
0078
0079
0080 switch lower(options.koType)
0081 case 'rxns'
0082
0083 model.selRxnMatrix = selMatrix(selSelectedRxns)';
0084 possibleKOList = model.rxns(selSelectedRxns);
0085
0086 case 'genesets'
0087
0088
0089 possibleKOList = unique(model.grRules(selSelectedRxns));
0090 if isempty(possibleKOList{1}), possibleKOList = possibleKOList(2:end); end
0091 for i = 1:length(possibleKOList)
0092 model.selRxnMatrix(:,i) = double(strcmp(possibleKOList{i},model.grRules));
0093 end
0094
0095 case 'genes'
0096
0097 model.selRxnMatrix = model.rxnGeneMat;
0098 possibleKOList = model.genes;
0099 otherwise
0100 error('Unrecognized KO type')
0101 end
0102
0103
0104 if ~isfield(options,'koCost')
0105 if isfield(model,'koCost')
0106 if length(model.koCost) == 1
0107 options.koCost = ones(length(possibleKOList,1)) * model.koCost;
0108 else
0109 options.koCost = model.koCost;
0110 end
0111 else
0112 options.koCost = ones(length(possibleKOList),1);
0113 end
0114 elseif length(model.koCost) == 1
0115 options.koCost = ones(length(possibleKOList,1)) * model.koCost;
0116 else
0117 options.koCost = model.koCost;
0118 end
0119
0120
0121 [selExc] = findExcRxns(model,true,true);
0122
0123
0124 model.ub(isinf(model.ub)) = MAXFLUX;
0125 model.ub(model.ub>MAXFLUX) = MAXFLUX;
0126 model.lb(isinf(model.lb)) = -MAXFLUX;
0127 model.lb(model.lb<-MAXFLUX) = -MAXFLUX;
0128 model.rxnsPresent = ones(length(model.rxns),1);
0129
0130
0131 solution0 = fluxBalance(model,selTargetRxns);
0132
0133
0134
0135 [nMets nRxns] = size(model.S);
0136 nInt = size(model.selRxnMatrix,2);
0137
0138 jMu = find(~isinf(model.lb));
0139 nMu = length(jMu);
0140 jNu = find(~isinf(model.ub));
0141 nNu = length(jNu);
0142
0143 y0 = (model.selRxnMatrix' * ~model.rxnsPresent) ~= 0;
0144 y0 = repmat(y0,1,options.M);
0145 y0p = y0;
0146 iiter = 1;
0147 change = true;
0148 t1 = clock;
0149 while change
0150 change = false;
0151
0152 y = false(nInt, 0);
0153 fmax = zeros(1, 0);
0154 for istart = 1:size(y0, 2)
0155 for irun = 1:options.M
0156
0157 c = [selTargetRxns;
0158 zeros(nMets, 1);
0159 zeros(nMu, 1);
0160 zeros(nNu, 1);
0161 zeros(nRxns, 1);
0162 zeros(nInt, 1)];
0163
0164 A = [ model.S sparse(nMets, nMets) sparse(nMets, nMu) sparse(nMets, nNu) sparse(nMets, nRxns) sparse(nMets, nInt);
0165 sparse(nRxns, nRxns) model.S' -sparse(jMu, 1:nMu, ones(nMu, 1), nRxns, nMu) sparse(jNu, 1:nNu, ones(nNu, 1), nRxns, nNu) speye(nRxns) sparse(nRxns, nInt);
0166 speye(nRxns) sparse(nRxns, nMets) sparse(nRxns, nMu) sparse(nRxns, nNu) sparse(nRxns, nRxns) model.selRxnMatrix .* repmat(model.lb, 1, nInt);
0167 speye(nRxns) sparse(nRxns, nMets) sparse(nRxns, nMu) sparse(nRxns, nNu) sparse(nRxns, nRxns) model.selRxnMatrix .* repmat(model.ub, 1, nInt);
0168 sparse(nRxns, nRxns) sparse(nRxns, nMets) sparse(nRxns, nMu) sparse(nRxns, nNu) speye(nRxns) model.selRxnMatrix * MAXDUAL;
0169 sparse(nRxns, nRxns) sparse(nRxns, nMets) sparse(nRxns, nMu) sparse(nRxns, nNu) speye(nRxns) -model.selRxnMatrix * MAXDUAL;
0170 model.c' sparse(1, nMets) model.lb(jMu)' -model.ub(jNu)' sparse(1, nRxns) sparse(1, nInt);
0171 sparse(1, nRxns) sparse(1, nMets) sparse(1, nMu) sparse(1, nNu) sparse(1, nRxns) ((y0(:, istart) == 0) - (y0(:, istart) == 1))';
0172 sparse(size(y, 2), nRxns) sparse(size(y, 2), nMets) sparse(size(y, 2), nMu) sparse(size(y, 2), nNu) sparse(size(y, 2), nRxns) ((y == 0) - (y == 1))';
0173 sparse(1, nRxns) sparse(1, nMets) sparse(1, nMu) sparse(1, nNu) sparse(1, nRxns) options.koCost'; ];
0174 b = [ zeros(nMets, 1);
0175 model.c;
0176 model.lb;
0177 model.ub;
0178 zeros(nRxns, 1);
0179 zeros(nRxns, 1);
0180 0;
0181 options.nbhdsz - nnz(y0(:, istart));
0182 ones(size(y, 2), 1) - sum((y ~= 0), 1)';
0183 options.maxKO; ];
0184 csense = char(['E' * ones(nMets, 1);
0185 'E' * ones(nRxns, 1);
0186 'G' * ones(nRxns, 1);
0187 'L' * ones(nRxns, 1);
0188 'G' * ones(nRxns, 1);
0189 'L' * ones(nRxns, 1);
0190 'E';
0191 'L';
0192 'G' * ones(size(y, 2), 1);
0193 'L'; ]);
0194 lb = [ model.lb;
0195 -Inf * ones(nMets, 1);
0196 zeros(nMu, 1);
0197 zeros(nNu, 1);
0198 -Inf * ones(nRxns, 1);
0199 zeros(nInt, 1) ];
0200 ub = [ model.ub;
0201 Inf * ones(nMets, 1);
0202 Inf * ones(nMu, 1);
0203 Inf * ones(nNu, 1);
0204 Inf * ones(nRxns, 1);
0205 ones(nInt, 1) ];
0206 vartype = char(['C' * ones(nRxns, 1);
0207 'C' * ones(nMets, 1);
0208 'C' * ones(nMu, 1);
0209 'C' * ones(nNu, 1);
0210 'C' * ones(nRxns, 1);
0211 'B' * ones(nInt, 1); ]);
0212 osense = -1;
0213
0214 if isfield(options,'minGrowth')
0215 A = [A; model.c' sparse(1, nMets + nMu + nNu + nRxns + nInt)];
0216 b = [b; options.minGrowth];
0217 csense = [csense; 'G'];
0218 end
0219
0220 [bilevelMILPProblem.c, bilevelMILPProblem.A,...
0221 bilevelMILPProblem.b, bilevelMILPProblem.lb,...
0222 bilevelMILPProblem.ub, bilevelMILPProblem.csense,...
0223 bilevelMILPProblem.vartype, bilevelMILPProblem.osense,...
0224 bilevelMILPProblem.x0] = ...
0225 deal(c, A, b, lb, ub, csense, vartype, osense, []);
0226
0227
0228 solution1 = solveCobraMILP(bilevelMILPProblem);
0229
0230
0231 if solution1.stat~=1
0232 continue;
0233 end
0234 yt = solution1.full((end - nInt + 1):end) > EPS;
0235
0236 model.rxnsPresent = ~(model.selRxnMatrix * yt);
0237 solution2 = fluxBalance(model,selTargetRxns,false);
0238 if abs(solution2.obj - solution1.obj) > EPS
0239 continue;
0240 end
0241
0242 fmax(:, end + 1) = solution1.obj;
0243 y(:, end + 1) = yt;
0244 end
0245 end
0246
0247 if size(y, 2) == 0
0248 continue;
0249 end
0250
0251 [fmaxsort, ifmaxsort] = sort(fmax);
0252 y = y(:, ifmaxsort);
0253 y = y(:, max([1 size(y, 2) - options.M + 1]):end);
0254 y = shrinkKnockouts(y, model, selTargetRxns);
0255 if size(y, 2) ~= size(y0, 2) || any(any(y~=y0))&&any(any(y~=y0p))
0256 y0p=y0;
0257 y0 = y;
0258 change = true;
0259 end
0260
0261 fprintf('Iteration %d\n', iiter);
0262 fprintf('----------%s\n', char('-' * ones(1, floor(log10(iiter)) + 1)));
0263 for iend = 1:size(y0, 2)
0264 model.rxnsPresent = ~(model.selRxnMatrix * y0(:, iend));
0265 [solSynMax solSynMin solBiomass] = fluxBalance(model,selTargetRxns);
0266 fprintf('Knockout cost: %d\n', options.koCost' * y0(:, iend));
0267 if nnz(y0) > 0
0268 fprintf('Knockouts:\n%s', sprintf('\t%s\n', possibleKOList{y0(:, iend)}));
0269 end
0270 printLabeledData(model.rxns(selExc),solSynMax.full(selExc),true);
0271 fprintf('\n');
0272
0273
0274 gdlsSolutionStructs.(sprintf('Iteration_%d',iiter)).(sprintf('solution_%2',i)).solBiomass = solBiomass;
0275 gdlsSolutionStructs.(sprintf('Iteration_%d',iiter)).(sprintf('solution_%2',i)).solSynMin = solSynMin;
0276 gdlsSolutionStructs.(sprintf('Iteration_%d',iiter)).(sprintf('solution_%2',i)).solSynMax = solSynMax;
0277 end
0278 elapsed_time = etime(clock,t1)
0279 iiter = iiter + 1;
0280 if (elapsed_time >= options.timeLimit) || (iiter >= options.iterationLimit)
0281 break;
0282 end
0283 end
0284
0285
0286 fprintf('\nGenerating Output\n');
0287 gdlsSolution.int = y0;
0288 gdlsSolutions.KOs = cell(max(sum(y0,2)),size(y0,2));
0289 for i = 1:size(y0,2)
0290 gdlsSolution.KOs(1:nnz(y0(:,i)),i) = possibleKOList(y0(:,i));
0291 [solSynMax solSynMin solBiomass] = fluxBalance(model,selTargetRxns,false);
0292 gdlsSolution.biomass(1,i) = solBiomass.obj;
0293 gdlsSolution.minTargetProd(1,i) = solSynMin.obj;
0294 gdlsSolution.maxTargetProd(1,i) = solSynMax.obj;
0295 end
0296
0297 end
0298 function y = shrinkKnockouts(y, model, selTargetRxns)
0299
0300 for iycol = 1:size(y, 2)
0301 model.rxnsPresent = ~(model.selRxnMatrix * y(:, iycol));
0302 solution1 = fluxBalance(model, selTargetRxns, false);
0303 for i = find(y(:, iycol))'
0304 yt = y(:, iycol);
0305 yt(i) = 0;
0306
0307 model.rxnsPresent = ~(model.selRxnMatrix * yt);
0308 solution2 = fluxBalance(model, selTargetRxns, false);
0309
0310 if solution2.obj >= solution1.obj
0311 y(:, iycol) = yt;
0312
0313 y(:, iycol) = shrinkKnockouts(y(:, iycol), model, selTargetRxns);
0314 end
0315 end
0316 end
0317
0318 y = unique(y', 'rows')';
0319 end
0320
0321 function [solSynMax solSynMin solBiomass] = fluxBalance(model,selTargetRxns,verbFlag)
0322 if nargin < 3
0323 verbFlag = true;
0324 end
0325 model.x0 = [];
0326 modelb = model;
0327 model_syn = model;
0328
0329 [nMets nRxns] = size(model.S);
0330
0331 yt = model.rxnsPresent;
0332
0333 modelb.A = [ model.S;
0334 sparse(1:nnz(~yt), find(~yt), ones(nnz(~yt), 1), nnz(~yt), nRxns) ];
0335 modelb.b = [ zeros(nMets, 1);
0336 zeros(nnz(~yt), 1) ];
0337 modelb.csense = char('E' * ones(1, nMets + nnz(~yt)));
0338 modelb.vartype = char('C' * ones(1, nRxns));
0339 modelb.osense = -1;
0340 solBiomass = solveCobraMILP(modelb);
0341
0342
0343 model_syn.A = [ model.S;
0344 sparse(1:nnz(~yt), find(~yt), ones(nnz(~yt), 1), nnz(~yt), nRxns);
0345 model.c' ];
0346 model_syn.b = [ zeros(nMets, 1);
0347 zeros(nnz(~yt), 1);
0348 solBiomass.obj ];
0349 model_syn.c = selTargetRxns;
0350 model_syn.csense = char('E' * ones(1, nMets + nnz(~yt) + 1));
0351 model_syn.vartype = char('C' * ones(1, nRxns));
0352 model_syn.osense = -1;
0353 solSynMax = solveCobraMILP(model_syn);
0354 model_syn.osense = 1;
0355 solSynMin = solveCobraMILP(model_syn);
0356
0357 if verbFlag
0358 fprintf('Biomass flux: %f\n', solBiomass.obj);
0359 fprintf('Synthetic flux: [%f, %f]\n', solSynMin.obj, solSynMax.obj);
0360 end
0361 end