0001 function [OMNISol,bilevelMILPProblem] = OMNI(model, selectedRxnList, options, constrOpt, measOpt, prevSolutions, verbFlag)
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 global MILPproblemType;
0054 global selectedRxnIndIrrev;
0055
0056 global irrev2rev;
0057
0058
0059
0060
0061
0062
0063
0064 if (nargin < 5)
0065 prevSolutions = [];
0066 end
0067 if (nargin < 6)
0068 verbFlag = false;
0069 end
0070
0071
0072
0073
0074
0075
0076
0077 [modelIrrev,matchRev,rev2irrev,irrev2rev] = convertToIrreversible(model);
0078
0079
0080
0081 selPrevSolIrrev = [];
0082 for i = 1:size(prevSolutions,2)
0083 prevSolRxnList = model.rxns(prevSolutions(:,i)==1);
0084 selPrevSol = ismember(model.rxns,prevSolRxnList);
0085 selPrevSolIrrev(:,i) = selPrevSol(irrev2rev);
0086 end
0087
0088 [nMets,nRxns] = size(modelIrrev.S);
0089
0090
0091
0092 selSelectedRxn = ismember(model.rxns,selectedRxnList);
0093 selSelectedRxnIrrev = selSelectedRxn(irrev2rev);
0094 selectedRxnIndIrrev = find(selSelectedRxnIrrev);
0095 cnt = 0;
0096
0097 nSelected = length(selectedRxnIndIrrev);
0098 selRxnCnt = 1;
0099 while selRxnCnt <= nSelected
0100 rxnID = selectedRxnIndIrrev(selRxnCnt);
0101 if (matchRev(rxnID)>0)
0102 cnt = cnt + 1;
0103 selectedRxnMatch(cnt,1) = selRxnCnt;
0104 selectedRxnMatch(cnt,2) = selRxnCnt+1;
0105 selRxnCnt = selRxnCnt + 1;
0106 end
0107 selRxnCnt = selRxnCnt + 1;
0108 end
0109
0110
0111 constrOptIrrev = setConstraintsIrrevModel(constrOpt,model,modelIrrev,rev2irrev);
0112
0113
0114
0115
0116 cLinear = zeros(nRxns,1);
0117 cInteger = zeros(sum(selSelectedRxnIrrev),1);
0118
0119
0120
0121
0122
0123
0124
0125 sel_meas_rxn = measOpt.rxnSel';
0126 b_meas_rxn = measOpt.values';
0127 wt_meas_rxn = measOpt.weights';
0128 n_m = length(sel_meas_rxn);
0129
0130
0131
0132
0133 sel_m = zeros(nRxns,1);
0134 ord_ir = [];
0135 b_meas_tmp = [];
0136 wt_meas_tmp = [];
0137 for i = 1:n_m
0138 rxn_name = sel_meas_rxn{i};
0139 rxn_id = find(strcmp(model.rxns,rxn_name));
0140 if (~isempty(rxn_id))
0141 b_meas_tmp = [b_meas_tmp;b_meas_rxn(i)];
0142 wt_meas_tmp = [wt_meas_tmp;wt_meas_rxn(i)];
0143
0144 if (model.rev(rxn_id))
0145 rxn_id_ir = rev2irrev{rxn_id}(1);
0146 sel_m(rxn_id_ir) = 1;
0147 sel_m(rxn_id_ir+1) = -1;
0148 else
0149
0150 rxn_id_ir = rev2irrev{rxn_id};
0151 sel_m(rxn_id_ir) = 1;
0152 end
0153
0154 ord_ir = [ord_ir rxn_id_ir];
0155 end
0156 end
0157
0158 [tmp,ord_ind] = sort(ord_ir);
0159
0160 if (sum(wt_meas_rxn) == 0)
0161 measOpts.weights = ones(n_m,1);
0162 else
0163 measOpts.weights = wt_meas_tmp(ord_ind);
0164 end
0165
0166 measOpts.values = b_meas_tmp(ord_ind);
0167
0168 measOpts.rxnSel = sel_m;
0169
0170
0171 bilevelMILPProblem = createBilevelMILPproblem(modelIrrev,cLinear,cInteger,selSelectedRxnIrrev,...
0172 selectedRxnMatch,constrOptIrrev,measOpts,options,selPrevSolIrrev);
0173
0174
0175
0176 if isfield(options,'initSolution')
0177 if (length(options.initSolution) > options.numDel | ~all(ismember(options.initSolution,selectedRxnList)))
0178 warning('Initial solution not valid - starting from a random initial solution')
0179 bilevelMILPProblem.x0 = [];
0180 else
0181
0182 selInitRxn = ismember(model.rxns,options.initSolution);
0183 selInitRxnIrrev = selInitRxn(irrev2rev);
0184 initRxnIndIrrev = find(selInitRxnIrrev);
0185 initIntegerSol = ~ismember(selectedRxnIndIrrev,initRxnIndIrrev);
0186 selInteger = bilevelMILPProblem.vartype == 'B';
0187 [nConstr,nVar] = size(bilevelMILPProblem.A);
0188 bilevelMILPProblem.x0 = nan(nVar,1);
0189 bilevelMILPProblem.x0(selInteger) = initIntegerSol;
0190
0191
0192
0193
0194
0195
0196
0197
0198
0199
0200
0201 end
0202 else
0203 bilevelMILPProblem.x0 = [];
0204 end
0205
0206
0207 bilevelMILPProblem.osense = 1;
0208
0209 if (verbFlag)
0210 [nConstr,nVar] = size(bilevelMILPProblem.A);
0211 nInt = length(bilevelMILPProblem.intSolInd);
0212 fprintf('MILP problem with %d constraints %d integer variables and %d continuous variables\n',...
0213 nConstr,nInt,nVar);
0214 end
0215
0216 bilevelMILPProblem.model = modelIrrev;
0217
0218
0219 MILPproblemType = 'OMNI';
0220
0221
0222
0223
0224
0225
0226
0227
0228 if (options.solveOMNI)
0229 OMNISol = solveCobraMILP(bilevelMILPProblem,'printLevel',0);
0230 if OMNISol.stat~=0
0231 if (~isempty(OMNISol.cont))
0232 OMNISol.fluxes = convertIrrevFluxDistribution(OMNISol.cont(1:length(matchRev)),matchRev);
0233 end
0234 if (~isempty(OMNISol.int))
0235
0236 OMNIRxnInd = selectedRxnIndIrrev(OMNISol.int < 1e-4);
0237 OMNISol.kos = model.rxns(unique(irrev2rev(OMNIRxnInd)));
0238
0239
0240
0241
0242
0243
0244
0245
0246 end
0247 else
0248 OMNISol.fluxes=[];
0249 OMNISol.kos={};
0250 end
0251 else
0252 OMNISol.rxnList = {};
0253 OMNISol.fluxes = [];
0254 end
0255
0256
0257