OptKnock

PURPOSE ^

OptKnock Run OptKnock in the most general form

SYNOPSIS ^

function [optKnockSol,bilevelMILPproblem] = OptKnock(model,selectedRxnList,options,constrOpt,prevSolutions,verbFlag,solutionFileNameTmp)

DESCRIPTION ^

OptKnock Run OptKnock in the most general form

 OptKnock(model,selectedRxnList,options,constrOpt,prevSolutions,verbFlag,solutionFileNameTmp)

INPUTS
 model            Structure containing all necessary variables to described a
                  stoichiometric model
   rxns            Rxns in the model
   mets            Metabolites in the model
   S               Stoichiometric matrix (sparse)
   b               RHS of Sv = b (usually zeros)
   c               Objective coefficients
   lb              Lower bounds for fluxes
   ub              Upper bounds for fluxes
   rev             Reversibility of fluxes

 selectedRxnList  List of reactions that can be knocked-out in OptKnock

 options          OptKnock options
   targetRxn       Target flux to be maximized

OPTIONAL INPUTS
 options             OptKnock options
   numDel             # of deletions allowed (Default: 5)
   numDelSense        Direction of # of deletions constraint (G/E/L)
                      (Default: L)
   vMax               Max flux (Default: 1000)
   solveOptKnock      Solve problem within Matlab (Default: true)
   createGams         Create GAMS input file
   gamsFile           GAMS input file name

 constrOpt           Explicitly constrained reaction options
   rxnList            Reaction list
   values             Values for constrained reactions
   sense              Constraint senses for constrained reactions (G/E/L)

 prevSolutions       Previous solutions

 verbFlag            Verbose flag

 solutionFileNameTmp File name for storing temporary solutions

OUTPUTS
 optKnockSol           OptKnock solution structure
   rxnList              Reaction KO list
   fluxes               Flux distribution
 bilevelMILPproblem    optKnock problem structure

NOTES
 OptKnock uses bounds of -vMax to vMax or 0 to vMax for reversible and
 irreversible reactions. If you wish to constrain a reaction, use
 constrOpt.

 Markus Herrgard 3/28/05
 Richard Que (04/27/10) - Added some default parameters.

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SOURCE CODE ^

0001 function [optKnockSol,bilevelMILPproblem] = OptKnock(model,selectedRxnList,options,constrOpt,prevSolutions,verbFlag,solutionFileNameTmp)
0002 %OptKnock Run OptKnock in the most general form
0003 %
0004 % OptKnock(model,selectedRxnList,options,constrOpt,prevSolutions,verbFlag,solutionFileNameTmp)
0005 %
0006 %INPUTS
0007 % model            Structure containing all necessary variables to described a
0008 %                  stoichiometric model
0009 %   rxns            Rxns in the model
0010 %   mets            Metabolites in the model
0011 %   S               Stoichiometric matrix (sparse)
0012 %   b               RHS of Sv = b (usually zeros)
0013 %   c               Objective coefficients
0014 %   lb              Lower bounds for fluxes
0015 %   ub              Upper bounds for fluxes
0016 %   rev             Reversibility of fluxes
0017 %
0018 % selectedRxnList  List of reactions that can be knocked-out in OptKnock
0019 %
0020 % options          OptKnock options
0021 %   targetRxn       Target flux to be maximized
0022 %
0023 %OPTIONAL INPUTS
0024 % options             OptKnock options
0025 %   numDel             # of deletions allowed (Default: 5)
0026 %   numDelSense        Direction of # of deletions constraint (G/E/L)
0027 %                      (Default: L)
0028 %   vMax               Max flux (Default: 1000)
0029 %   solveOptKnock      Solve problem within Matlab (Default: true)
0030 %   createGams         Create GAMS input file
0031 %   gamsFile           GAMS input file name
0032 %
0033 % constrOpt           Explicitly constrained reaction options
0034 %   rxnList            Reaction list
0035 %   values             Values for constrained reactions
0036 %   sense              Constraint senses for constrained reactions (G/E/L)
0037 %
0038 % prevSolutions       Previous solutions
0039 %
0040 % verbFlag            Verbose flag
0041 %
0042 % solutionFileNameTmp File name for storing temporary solutions
0043 %
0044 %OUTPUTS
0045 % optKnockSol           OptKnock solution structure
0046 %   rxnList              Reaction KO list
0047 %   fluxes               Flux distribution
0048 % bilevelMILPproblem    optKnock problem structure
0049 %
0050 %NOTES
0051 % OptKnock uses bounds of -vMax to vMax or 0 to vMax for reversible and
0052 % irreversible reactions. If you wish to constrain a reaction, use
0053 % constrOpt.
0054 %
0055 % Markus Herrgard 3/28/05
0056 % Richard Que (04/27/10) - Added some default parameters.
0057 
0058 % Set these for MILP callbacks
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 %idefault <= 5 deletions; solve OptKnock
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 % Convert to irreversible rxns
0099 [modelIrrev,matchRev,rev2irrev,irrev2rev] = convertToIrreversible(model);
0100 
0101 % Create the index of the previous KO's suggested by OptKnock to avoid obtaining the same
0102 % solution again
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 % Create matchings for reversible reactions in the set selected for KOs
0113 % This is to ensure that both directions of the reaction are knocked out
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 % Set inner constraints for the LP
0133 constrOptIrrev = setConstraintsIrrevModel(constrOpt,model,modelIrrev,rev2irrev);
0134     
0135 % Set objectives for linear and integer parts
0136 cLinear = zeros(nRxns,1);
0137 cInteger = zeros(sum(selSelectedRxnIrrev),1);
0138 
0139 % Set the correct objective coefficient
0140 targetRxnID = find(ismember(model.rxns,options.targetRxn));
0141 targetRxnIDirrev = rev2irrev{targetRxnID}(1);
0142 cLinear(targetRxnIDirrev) = 1;
0143 
0144 % Create the constraint matrices for the bilevel MILP
0145 bilevelMILPproblem = createBilevelMILPproblem(modelIrrev,cLinear,cInteger,selSelectedRxnIrrev,...
0146     selectedRxnMatch,constrOptIrrev,[],options,selPrevSolIrrev);
0147 
0148 % Initial guess (random)
0149 %bilevelMILPproblem.x0 = round(rand(length(bilevelMILPproblem.c),1));
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         % Set initial integer solution
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 %         LPproblem.b = bilevelMILPproblem.b - bilevelMILPproblem.A(:,selInteger)*initIntegerSol;
0166 %         LPproblem.A = bilevelMILPproblem.A(:,bilevelMILPproblem.vartype == 'C');
0167 %         LPproblem.c = bilevelMILPproblem.c(bilevelMILPproblem.vartype == 'C');
0168 %         LPproblem.lb = bilevelMILPproblem.lb(bilevelMILPproblem.vartype == 'C');
0169 %         LPproblem.ub = bilevelMILPproblem.ub(bilevelMILPproblem.vartype == 'C');
0170 %         LPproblem.osense = -1;
0171 %         LPproblem.csense = bilevelMILPproblem.csense;
0172 %         LPsol = solveCobraLP(LPproblem);
0173 %
0174 %         bilevelMILPproblem.x0(~selInteger) = LPsol.full;
0175     end
0176 else
0177     bilevelMILPproblem.x0 = [];
0178 end
0179 
0180 % Maximize
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 % Set these for CPLEX callbacks
0193 MILPproblemType = 'OptKnock';
0194 rxnList = model.rxns;
0195 biomassRxnID = find(modelIrrev.c==1);
0196 solID = 0;
0197 OptKnockObjective = [];
0198 OptKnockGrowth = [];
0199 OptKnockKOrxnList = {};
0200 
0201 % Solve problem
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         % Figure out the KO reactions
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

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