0001 function [optKnockSol,bilevelMILPproblem] = OptKnock(model,selectedRxnList,options,constrOpt,prevSolutions,verbFlag,solutionFileNameTmp)
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
0038
0039
0040
0041
0042
0043
0044
0045
0046
0047
0048
0049
0050
0051
0052
0053
0054
0055
0056
0057
0058
0059 global MILPproblemType;
0060 global selectedRxnIndIrrev;
0061 global rxnList;
0062 global irrev2rev;
0063 global solutionFileName;
0064 global biomassRxnID;
0065 global OptKnockKOrxnList;
0066 global OptKnockObjective;
0067 global OptKnockGrowth;
0068 global solID;
0069
0070
0071 if (~exist('options','var') || isempty(options) )
0072 error('OptKnock: No target reaction specified')
0073 else
0074 if ~isfield(options,'vMax'), options.vMax = 1000; end
0075 if ~isfield(options,'numDel'), options.numDel = 5; end
0076 if ~isfield(options,'numDelSense'), options.numDelSense = 'L'; end
0077 if ~isfield(options,'solveOptKnock'), options.solveOptKnock = true; end
0078 end
0079
0080 if ~exist('constrOpt','var')
0081 constrOpt.rxnInd = [];
0082 constrOpt.values = [];
0083 constrOpt.sense = [];
0084 end
0085
0086 if (nargin < 5)
0087 prevSolutions = [];
0088 end
0089 if (nargin < 6)
0090 verbFlag = false;
0091 end
0092 if (nargin < 7)
0093 solutionFileName = 'optKnockSolutions.mat';
0094 else
0095 solutionFileName = solutionFileNameTmp;
0096 end
0097
0098
0099 [modelIrrev,matchRev,rev2irrev,irrev2rev] = convertToIrreversible(model);
0100
0101
0102
0103 selPrevSolIrrev = [];
0104 for i = 1:length(prevSolutions)
0105 prevSolRxnList = prevSolutions{i};
0106 selPrevSol = ismember(model.rxns,prevSolRxnList);
0107 selPrevSolIrrev(:,i) = selPrevSol(irrev2rev);
0108 end
0109
0110 [nMets,nRxns] = size(modelIrrev.S);
0111
0112
0113
0114 selSelectedRxn = ismember(model.rxns,selectedRxnList);
0115 selSelectedRxnIrrev = selSelectedRxn(irrev2rev);
0116 selectedRxnIndIrrev = find(selSelectedRxnIrrev);
0117 cnt = 0;
0118 prevRxnID = -10;
0119 nSelected = length(selectedRxnIndIrrev);
0120 selRxnCnt = 1;
0121 while selRxnCnt <= nSelected
0122 rxnID = selectedRxnIndIrrev(selRxnCnt);
0123 if (matchRev(rxnID)>0)
0124 cnt = cnt + 1;
0125 selectedRxnMatch(cnt,1) = selRxnCnt;
0126 selectedRxnMatch(cnt,2) = selRxnCnt+1;
0127 selRxnCnt = selRxnCnt + 1;
0128 end
0129 selRxnCnt = selRxnCnt + 1;
0130 end
0131
0132
0133 constrOptIrrev = setConstraintsIrrevModel(constrOpt,model,modelIrrev,rev2irrev);
0134
0135
0136 cLinear = zeros(nRxns,1);
0137 cInteger = zeros(sum(selSelectedRxnIrrev),1);
0138
0139
0140 targetRxnID = find(ismember(model.rxns,options.targetRxn));
0141 targetRxnIDirrev = rev2irrev{targetRxnID}(1);
0142 cLinear(targetRxnIDirrev) = 1;
0143
0144
0145 bilevelMILPproblem = createBilevelMILPproblem(modelIrrev,cLinear,cInteger,selSelectedRxnIrrev,...
0146 selectedRxnMatch,constrOptIrrev,[],options,selPrevSolIrrev);
0147
0148
0149
0150 if isfield(options,'initSolution')
0151 if (length(options.initSolution) > options.numDel | ~all(ismember(options.initSolution,selectedRxnList)))
0152 warning('Initial solution not valid - starting from a random initial solution')
0153 bilevelMILPproblem.x0 = [];
0154 else
0155
0156 selInitRxn = ismember(model.rxns,options.initSolution);
0157 selInitRxnIrrev = selInitRxn(irrev2rev);
0158 initRxnIndIrrev = find(selInitRxnIrrev);
0159 initIntegerSol = ~ismember(selectedRxnIndIrrev,initRxnIndIrrev);
0160 selInteger = bilevelMILPproblem.vartype == 'B';
0161 [nConstr,nVar] = size(bilevelMILPproblem.A);
0162 bilevelMILPproblem.x0 = nan(nVar,1);
0163 bilevelMILPproblem.x0(selInteger) = initIntegerSol;
0164
0165
0166
0167
0168
0169
0170
0171
0172
0173
0174
0175 end
0176 else
0177 bilevelMILPproblem.x0 = [];
0178 end
0179
0180
0181 bilevelMILPproblem.osense = -1;
0182
0183 if (verbFlag)
0184 [nConstr,nVar] = size(bilevelMILPproblem.A);
0185 nInt = length(bilevelMILPproblem.intSolInd);
0186 fprintf('MILP problem with %d constraints %d integer variables and %d continuous variables\n',...
0187 nConstr,nInt,nVar);
0188 end
0189
0190 bilevelMILPproblem.model = modelIrrev;
0191
0192
0193 MILPproblemType = 'OptKnock';
0194 rxnList = model.rxns;
0195 biomassRxnID = find(modelIrrev.c==1);
0196 solID = 0;
0197 OptKnockObjective = [];
0198 OptKnockGrowth = [];
0199 OptKnockKOrxnList = {};
0200
0201
0202 if (options.solveOptKnock)
0203 optKnockSol = solveCobraMILP(bilevelMILPproblem,'printLevel',verbFlag);
0204 if (~isempty(optKnockSol.cont))
0205 optKnockSol.fluxes = convertIrrevFluxDistribution(optKnockSol.cont(1:length(matchRev)),matchRev);
0206 end
0207 if (~isempty(optKnockSol.int))
0208
0209 optKnockRxnInd = selectedRxnIndIrrev(optKnockSol.int < 1e-4);
0210 optKnockSol.rxnList = model.rxns(unique(irrev2rev(optKnockRxnInd)));
0211 end
0212 else
0213 optKnockSol.rxnList = {};
0214 optKnockSol.fluxes = [];
0215 end
0216
0217
0218