ftINITInternalAlg
This function runs the MILP for a step in ftINIT.
model a reference model structure
rxnScores a vector of scores for the reactions in the model.
Positive scores are reactions to keep and negative
scores are reactions to exclude. Rxns set to 0 are excluded
from the problem.
metData boolean matrix with mets as rows and rxns as columns
saying which reaction produces each detected met (opt, default [])
essentialRxns cell array of reactions that are essential and that
have to be in the resulting model. This is normally
used when fitting a model to task (see fitTasks)
prodWeight a score that determines the value of having
net-production of metabolites. This is a way of having
a more functional network as it provides a reason for
including bad reactions for connectivity reasons. This
score is for each metabolite, and the sum of these weights
and the scores for the reactions is what is optimized
allowExcretion true if excretion of all metabolites should be allowed.
This results in fewer reactions being considered
dead-ends, but all reactions in the resulting model may
not be able to carry flux. If this is "false" then the
equality constraints are taken from model.b. If the
input model lacks exchange reactions then this should
probably be "true", or a large proportion of the model
would be excluded for connectivity reasons
remPosRev If true, the positive reversible reactions are removed from the problem.
This is used in step 1 of ftINIT
params parameters for the MILP, for example MIPGap and TimeLimit
startVals Start values for the MILP, typically used when rerunning
with a higher MIPGap, to use the results from the previous
run
fluxes Fluxes from the last run.
verbose If true, the MILP progression will be shown.
deletedRxns reactions which were deleted by the algorithm (only
rxns included in the problem)
metProduction array that indicates which of the
metabolites in presentMets that could be
produced
0: metabolite could not be produced
1: metabolite could be produced
res The result from the MILP
turnedOnRxns The reactions determined to be present (only
rxns included in the problem)
fluxes The fluxes from the MILP
This function is the actual implementation of the algorithm. See
ftINIT for a higher-level function for model reconstruction.
Usage: [deletedRxns,metProduction,res,turnedOnRxns,fluxes]=runINIT9(model,...
rxnScores,presentMets,essentialRxns,prodWeight,allowExcretion,...
remPosRev,params)0001 function [deletedRxns,metProduction,res,turnedOnRxns,fluxes]=ftINITInternalAlg(model,rxnScores,metData,essentialRxns,prodWeight,allowExcretion,remPosRev,params,startVals,fluxes,verbose) 0002 % ftINITInternalAlg 0003 % This function runs the MILP for a step in ftINIT. 0004 % 0005 % model a reference model structure 0006 % rxnScores a vector of scores for the reactions in the model. 0007 % Positive scores are reactions to keep and negative 0008 % scores are reactions to exclude. Rxns set to 0 are excluded 0009 % from the problem. 0010 % metData boolean matrix with mets as rows and rxns as columns 0011 % saying which reaction produces each detected met (opt, default []) 0012 % essentialRxns cell array of reactions that are essential and that 0013 % have to be in the resulting model. This is normally 0014 % used when fitting a model to task (see fitTasks) 0015 % prodWeight a score that determines the value of having 0016 % net-production of metabolites. This is a way of having 0017 % a more functional network as it provides a reason for 0018 % including bad reactions for connectivity reasons. This 0019 % score is for each metabolite, and the sum of these weights 0020 % and the scores for the reactions is what is optimized 0021 % allowExcretion true if excretion of all metabolites should be allowed. 0022 % This results in fewer reactions being considered 0023 % dead-ends, but all reactions in the resulting model may 0024 % not be able to carry flux. If this is "false" then the 0025 % equality constraints are taken from model.b. If the 0026 % input model lacks exchange reactions then this should 0027 % probably be "true", or a large proportion of the model 0028 % would be excluded for connectivity reasons 0029 % remPosRev If true, the positive reversible reactions are removed from the problem. 0030 % This is used in step 1 of ftINIT 0031 % params parameters for the MILP, for example MIPGap and TimeLimit 0032 % startVals Start values for the MILP, typically used when rerunning 0033 % with a higher MIPGap, to use the results from the previous 0034 % run 0035 % fluxes Fluxes from the last run. 0036 % verbose If true, the MILP progression will be shown. 0037 % 0038 % deletedRxns reactions which were deleted by the algorithm (only 0039 % rxns included in the problem) 0040 % metProduction array that indicates which of the 0041 % metabolites in presentMets that could be 0042 % produced 0043 % 0: metabolite could not be produced 0044 % 1: metabolite could be produced 0045 % res The result from the MILP 0046 % turnedOnRxns The reactions determined to be present (only 0047 % rxns included in the problem) 0048 % fluxes The fluxes from the MILP 0049 % 0050 % This function is the actual implementation of the algorithm. See 0051 % ftINIT for a higher-level function for model reconstruction. 0052 % 0053 % Usage: [deletedRxns,metProduction,res,turnedOnRxns,fluxes]=runINIT9(model,... 0054 % rxnScores,presentMets,essentialRxns,prodWeight,allowExcretion,... 0055 % remPosRev,params) 0056 0057 if isempty(essentialRxns) 0058 essentialRxns={}; 0059 end 0060 essentialRxns=essentialRxns(:); 0061 if isempty(prodWeight) 0062 prodWeight=0.5; 0063 end 0064 0065 0066 %The model should be in the reversible format and all relevant exchange 0067 %reactions should be open 0068 if isfield(model,'unconstrained') 0069 EM='Exchange metabolites are still present in the model. Use simplifyModel if this is not intended'; 0070 dispEM(EM,false); 0071 end 0072 0073 0074 essential = ismember(model.rxns,essentialRxns); 0075 0076 %Some nice to have numbers 0077 nMets=numel(model.mets); 0078 nRxns=numel(model.rxns); 0079 nRxnsWithOnOff = nRxns; 0080 0081 %Reactions with score 0 will just be left in the model, and will not be part of the problem (but can carry flux). 0082 %It is possible to set score = 0 for e.g. spontaneous reactions, exchange rxns, etc., which may not be that interesting 0083 %to remove. 0084 0085 %if makeIrrev is on, we can just as well skip all positive reversible rxns - they will be turned on without carrying any flux 0086 %since they can form a loop within themself (fwd-rev) 0087 if remPosRev 0088 rxnScores(rxnScores > 0 & model.rev~=0) = 0; 0089 end 0090 0091 %Handle metabolomics 0092 %A problem with the metabolomics is that some of the producer reactions 0093 %for a metabolite could be excluded from the problem. We solve this by 0094 %adding them with score 0. They can still be seen as positive, since they 0095 %either don't matter (there is another producer) or they can contribute to an 0096 %increased score. 0097 revRxns=model.rev~=0; 0098 essRevRxns = find(revRxns & essential); 0099 essIrrevRxns = find(~revRxns & essential); 0100 0101 if ~isempty(metData) 0102 %remove any metData rows that is connected to an essential rxn - 0103 %these will be on regardless and will cause problems below. 0104 containsEssential = any(metData(:,essential),2); 0105 metData = metData(~containsEssential,:); 0106 end 0107 0108 if ~isempty(metData) 0109 metRxns = any(metData,1).'; 0110 posRxns = rxnScores > 0 | ((rxnScores == 0) & metRxns); 0111 else 0112 posRxns = rxnScores > 0; 0113 end 0114 negRxns = rxnScores < 0; 0115 0116 posRevRxns = find(posRxns & revRxns & ~essential); 0117 negRevRxns = find(negRxns & revRxns & ~essential); 0118 posIrrevRxns = find(posRxns & ~revRxns & ~essential); 0119 negIrrevRxns = find(negRxns & ~revRxns & ~essential); 0120 nPosRev = numel(posRevRxns); 0121 nNegRev = numel(negRevRxns); 0122 nPosIrrev = numel(posIrrevRxns); 0123 nNegIrrev = numel(negIrrevRxns); 0124 nEssRev = numel(essRevRxns); 0125 nEssIrrev = numel(essIrrevRxns); %not used, but left for symmetry 0126 nMetabolMets = size(metData,1); 0127 0128 if ~isempty(metData) 0129 metNegRev = negRxns & revRxns & ~essential & metRxns; 0130 metNegIrrev = negRxns & ~revRxns & ~essential & metRxns; 0131 nMetNegRev = sum(metNegRev); 0132 nMetNegIrrev = sum(metNegIrrev); 0133 end 0134 0135 0136 milpModel = model; 0137 0138 %These six categories need to be handled separately: 0139 % 0140 %PosIrrev (positive reaction score, irreversible): 0141 %flux >= 0.1*Yi, flux - 0.1*Yi - VPI == 0, 0 <= Yi(onoff) <= 1, 0 <= VPI <= Inf 0142 %The nice thing with rxns with positive score is that they do not require a boolean. The optimizer will 0143 %strive for maximizing the Yi here, so it can be a continuous variable, where we just force a flux if 0144 %the variable is on. 0145 %PosRev: 0146 % 1: Split up the flux into positive and negative: flux - vnrp + vnrn == 0, 0 <= vprp,vprn <= Inf. 0147 % 2. Force one of the fluxes on: vprp + vprn >= 0.1*Yi, vprp + vprn - 0.1*Yi - vprv1 == 0. 0 <= Yi(onoff) <= 1, vprv1 >= 0 0148 % 3. We need a bool to make sure that one of vprp and vprn are always zero: 0149 % vprb E {0,1}: vprp <= 100*vprb (if bool is 0, vprp is zero): vprp - 100*vprb + vprv2 == 0, vprv2 >= 0 0150 % vprn <= (1-vprb)*100 (if bool is one, vprn is zero): vprn + 100*vprb + vprv3 == 0, -100 <= vprv3 <= inf 0151 %Neg Irrev: 0152 %For negative scores, we need to use an integer, I see no way around that. 0153 % flux < 100*Yi, Yi E {0,1}. Flux - 100*Yi + vni == 0 0154 %Neg Rev: 0155 %Two cases: 0156 %A: Without irrev model 0157 % 1: Split up the flux into positive and negative: flux - vnrp + vnrn == 0, 0 <= vprp,vprn <= Inf. 0158 % 2: Force the Yi (on/off) var on if the flux is on: vnrp + vnrn <= 0.1 * Yi, Yi E {0,1}: vnrp + vnrn - 0.1 * Yi + vnrv1 == 0, vnrv1 >= 0. 0159 %B: Irrev model (Not used for now) 0160 % 1. We now have two reactions, but still just boolean variable. We know that 0161 % the flux can only be positive, so we just say that flux1 + flux2 <= 0.1 * Yi, Yi E {0,1}: flux1 + flux2 - 0.1 * Yi + vnrv1 == 0, vnrv1 >= 0. 0162 % We don't care if there is any loop - there is just no benefit for the objective to generate one, and it doesn't matter. 0163 %Ess Irrev (essential rxn, i.e. forced to carry flux, irreversible): 0164 %flux >= 0.1 - solved with lb = 0.1. 0165 %Ess Rev (these are not really used): 0166 % 1: Split up the flux into positive and negative: flux - verp + vern == 0, 0 <= verp,vern <= Inf. 0167 % 2. Force one of the fluxes on: verp + vern >= 0.1, verp + vern - verv1 == 0. verv1 >= 0.1 0168 % 3. We need a bool to make sure that one of verp and vern are always zero: 0169 % verb E {0,1}: verp <= 100*verb (if bool is 0, verp is zero): verp - 100*verb + verv2 == 0, verv2 >= 0 0170 % vern <= (1-verb)*100 (if bool is one, vern is zero): vern + 100*verb + verv3 == 0, -100 <= verv3 <= inf 0171 0172 % Now metabolomics 0173 % The basic idea is that the bonus is gotten if any of the producer reactions are on. We can 0174 % therefore use the "on" variables directly for the reactions that are included. We therefore add one variable per met, mon, which can be 0175 % continuous and between 0 and 1. We then say that mon <= v1 + v2 + ... + vn, i.e. that 0176 % mon + mv1 - v1 - v2 ... - vn == 0, 0 <= mon <= 1, 0 <= mv1 <= inf 0177 % mon gets -prodWeight in the c vector (i.e. objective function) 0178 % A tricky part is that some of the reactions producing a metabolite are left outside the problem 0179 % This is solved by moving them into the problem, but with the score 0. They can still be treated as positive 0180 % reactions. In the end, we are not interested if they are on though, but they may 0181 % allow for production of a metabolite, giving no benefit of turning on another producer reaction. 0182 % Another tricky thing is that some metabolite producers have negative score. There is 0183 % no automatic mechanism for forcing flux on if the variable is on - this is simply 0184 % not needed for negative variables - unless they are a metabolite producer. We therefore need 0185 % to add an extra constraint there, similar to the case of positive. It gets a bit complicated. 0186 % For reversible, we need a bool to make sure that one of vnrp and vnrn stays zero: 0187 % vnrbm E {0,1}: vnrp <= 100*vnrbm (if bool is 0, vnrp is zero): vnrp - 100*vnrbm + vnrvm1 == 0, vnrvm1 >= 0 0188 % vnrn <= (1-vnrbm)*100 (if bool is one, vnrn is zero): vnrn + 100*vnrbm + vnrvm2 == 0, -100 <= vnrvm2 <= inf 0189 % We then also say that vnrp + vnrn >= 0.1*Yi, vnrp + vnrn - 0.1*Yi - vnrvm3 == 0, vnrvm3 >= 0 0190 % For irreversible, it is fairly straight-forward. We just say that flux >= 0.1*Yi, i.e. flux - 0.1*Yi - vnim == 0, 0 <= vnim <= Inf. 0191 0192 0193 %The total matrix then looks like this. Note that we have ordered the categories, two reactions each, 0194 %in the S matrix to make it easier to follow the figure. In practice, they come in random order, which 0195 %is why the SEye variable is used further down. 0196 %all categories below have two columns - the two column starts where label starts 0197 %since the label cannot fit in two chars, they are on multiple lines 0198 % 0199 %The met variables and constraints are not in this figure for practical reasons. 0200 % 0201 % vprb vprv3 vnrn vern verv2 mv1 0202 % pi ni ei Ypi Yni vpi vprn vprv2 vnrp verp verb mon 0203 % pr nr er Ypr Ynr vprp vprv1 vni vnrv1 verv1 verv3 0204 % SSSSSSSSSSSS0000000000000000000000000000000000000000000000 0205 % SSSSSSSSSSSS0000000000000000000000000000000000000000000000 S block 0206 % SSSSSSSSSSSS0000000000000000000000000000000000000000000000 0207 % 1 N - Pos irrev block 0208 % 1 N - 0209 % 1 - 1 Pos rev block 1 0210 % 1 - 1 0211 % N 1 1 - Pos rev block 2 0212 % N 1 1 - 0213 % 1 M 1 Pos rev block 3 0214 % 1 M 1 0215 % 1 C 1 Pos rev block 4 0216 % 1 C 1 0217 % 1 M 1 Neg irrev block 0218 % 1 M 1 0219 % 1 - 1 Neg rev block 1 0220 % 1 - 1 0221 % M 1 1 1 Neg rev block 2 0222 % M 1 1 1 0223 % 1 - 1 Ess rev block 1 0224 % 1 - 1 0225 % 1 1 - Ess rev block 2 0226 % 1 1 - 0227 % 1 M 1 Ess rev block 3 0228 % 1 M 1 0229 % 1 C 1 Ess rev block 4 0230 % 1 C 1 0231 % -------- 1 1 Met block - Here, we assume that all variables support each met. 0232 % -------- 1 1 In practice, fewer of the "on" variables with -1 should be included 0233 % 0234 % 0235 %M = -100 0236 %N = -0.1 0237 %D = 0.1 0238 %C = 100 0239 %- = -1 0240 % 0241 %It looks a bit different for the irrev model, but still quite similar 0242 0243 %build the A matrix 0244 %S row 0245 %When forcing on essential rxns, use the flux value of the previous run (set to 0.1 the first time) 0246 %Don't set it above 0.1, may starve something else out. Leave a margin of 1% from the last run. 0247 forceOnLim = 0.1; 0248 forceOnLimEss=min(abs(fluxes)*0.99,0.1); 0249 0250 varsPerNegRev = 3; 0251 0252 0253 %Figure out the number of variables needed for metabolomics 0254 if ~isempty(metData) 0255 nMetVars = 2*size(metData,1) + 4*nMetNegRev + nMetNegIrrev; %mon,mv1; vnrbm,vnrvm1,vnrvm2,vnrvm3; vnim 0256 else 0257 nMetVars = 0; 0258 end 0259 0260 0261 sRow = [milpModel.S sparse(nMets, nPosIrrev*2+nPosRev*7+nNegIrrev*2+nNegRev*(1+varsPerNegRev)+nEssRev*6 + nMetVars)]; 0262 sEye = speye(nRxns); 0263 nYBlock = nPosIrrev+nPosRev+nNegIrrev+nNegRev; 0264 piRows = [sEye(posIrrevRxns,:) speye(nPosIrrev)*-forceOnLim sparse(nPosIrrev,nPosRev+nNegIrrev+nNegRev) speye(nPosIrrev)*-1 sparse(nPosIrrev,nPosRev*6+nNegIrrev+nNegRev*varsPerNegRev+nEssRev*6+nMetVars)]; 0265 0266 prRows1 = [sEye(posRevRxns,:) sparse(nPosRev,nYBlock + nPosIrrev) speye(nPosRev)*-1 speye(nPosRev) sparse(nPosRev,nPosRev*4+nNegIrrev+nNegRev*varsPerNegRev+nEssRev*6+nMetVars)]; 0267 prRows2 = [sparse(nPosRev,nRxns + nPosIrrev) speye(nPosRev)*-forceOnLim sparse(nPosRev,nNegIrrev + nNegRev + nPosIrrev) speye(nPosRev) speye(nPosRev) sparse(nPosRev,nPosRev) speye(nPosRev)*-1 sparse(nPosRev,nPosRev*2+nNegIrrev+nNegRev*varsPerNegRev+nEssRev*6+nMetVars)]; 0268 prRows3 = [sparse(nPosRev,nRxns + nYBlock + nPosIrrev) speye(nPosRev) sparse(nPosRev,nPosRev) speye(nPosRev)*-100 sparse(nPosRev,nPosRev) speye(nPosRev) sparse(nPosRev,nPosRev+nNegIrrev+nNegRev*varsPerNegRev+nEssRev*6+nMetVars)]; 0269 prRows4 = [sparse(nPosRev,nRxns + nYBlock + nPosIrrev + nPosRev) speye(nPosRev) speye(nPosRev)*100 sparse(nPosRev,nPosRev*2) speye(nPosRev) sparse(nPosRev, nNegIrrev+nNegRev*varsPerNegRev+nEssRev*6+nMetVars)]; 0270 0271 niRows = [sEye(negIrrevRxns,:) sparse(nNegIrrev,nPosIrrev + nPosRev) speye(nNegIrrev)*-100 sparse(nNegIrrev,nNegRev + nPosIrrev + nPosRev*6) speye(nNegIrrev) sparse(nNegIrrev,nNegRev*varsPerNegRev + nEssRev*6+nMetVars)]; 0272 0273 nrRows1 = [sEye(negRevRxns,:) sparse(nNegRev,nYBlock + nPosIrrev + nPosRev*6 + nNegIrrev) speye(nNegRev)*-1 speye(nNegRev) sparse(nNegRev,nNegRev + nEssRev*6+nMetVars)]; 0274 nrRows2 = [sparse(nNegRev,nRxns + nPosIrrev + nPosRev + nNegIrrev) speye(nNegRev)*-100 sparse(nNegRev,nPosIrrev + nPosRev*6 + nNegIrrev) speye(nNegRev) speye(nNegRev) speye(nNegRev) sparse(nNegRev,nEssRev*6+nMetVars)]; 0275 nrRows = [nrRows1;nrRows2]; 0276 0277 erRows1 = [sEye(essRevRxns,:) sparse(nEssRev,nYBlock + nPosIrrev + nPosRev*6 + nNegIrrev + nNegRev*varsPerNegRev) speye(nEssRev)*-1 speye(nEssRev) sparse(nEssRev, nEssRev*4+nMetVars)]; 0278 erRows2 = [sparse(nEssRev,nRxns + nYBlock + nPosIrrev + nPosRev*6 + nNegIrrev + nNegRev*varsPerNegRev) speye(nEssRev) speye(nEssRev) speye(nEssRev)*-1 sparse(nEssRev, nEssRev*3+nMetVars)]; 0279 erRows3 = [sparse(nEssRev,nRxns + nYBlock + nPosIrrev + nPosRev*6 + nNegIrrev + nNegRev*varsPerNegRev) speye(nEssRev) sparse(nEssRev, nEssRev*2) speye(nEssRev)*-100 speye(nEssRev) sparse(nEssRev, nEssRev+nMetVars)]; 0280 erRows4 = [sparse(nEssRev,nRxns + nYBlock + nPosIrrev + nPosRev*6 + nNegIrrev + nNegRev*varsPerNegRev + nEssRev) speye(nEssRev) sparse(nEssRev, nEssRev) speye(nEssRev)*100 sparse(nEssRev, nEssRev) speye(nEssRev) sparse(nEssRev, nMetVars)]; 0281 0282 %now the mets 0283 if ~isempty(metData) 0284 %Order of rxn "on" vars: Ypi Ypr Yni Ynr. This is a bit messy, since they are not in the 0285 %same order as the reactions in the S matrix or metData. We therefore resort the vars in 0286 %metData to match the order of the vars. 0287 srtMetData = sparse([metData(:,posIrrevRxns) metData(:,posRevRxns) metData(:,negIrrevRxns) metData(:,negRevRxns)]); 0288 0289 %First the setup for giving bonus if the met is included. 0290 %The mon vars come first followed by the mv1 vars 0291 metRows1 = [sparse(nMetabolMets,nRxns) -srtMetData sparse(nMetabolMets, nPosIrrev + nPosRev*6 + nNegIrrev + nNegRev*varsPerNegRev + nEssRev*6) speye(nMetabolMets) speye(nMetabolMets) sparse(nMetabolMets, 4*nMetNegRev + nMetNegIrrev)]; 0292 %Then negative rev: 0293 % Variable order: eye(vnrbm),eye(vnrvm1),eye(vnrvm2),eye(vnrvm3) 0294 % vnrp,vnrn are the two first variables among the neg rev vars. 0295 % For reversible, we need a bool to make sure that one of vnrp and vnrn stays zero: 0296 % vnrbm E {0,1}: vnrp <= 100*vnrbm (if bool is 0, vnrp is zero): vnrp - 100*vnrbm + vnrvm1 == 0, vnrvm1 >= 0 (metRows2) 0297 % vnrn <= (1-vnrbm)*100 (if bool is one, vnrn is zero): vnrn + 100*vnrbm + vnrvm2 == 0, -100 <= vnrvm2 <= inf (metRows3) 0298 % We then also say that vnrp + vnrn >= 0.1*Yi, -0.1*Yi + vnrp + vnrn - vnrvm3 == 0, vnrvm3 >= 0 (metRows4) 0299 nrEye = speye(nNegRev); 0300 %vnrp - 100*vnrbm + vnrvm1 == 0 0301 metRows2 = [sparse(nMetNegRev,nRxns + nYBlock + nPosIrrev + nPosRev*6 + nNegIrrev) ... %zeros up to vnrp 0302 nrEye(metNegRev(negRevRxns),:) ... %vnrp 0303 sparse(nMetNegRev, nNegRev*(varsPerNegRev-1) + nEssRev*6 + nMetabolMets*2) ... %zeros up to vnrbm 0304 speye(nMetNegRev)*-100 ... %vnrbm 0305 speye(nMetNegRev) ... % vnrvm1 0306 sparse(nMetNegRev, nMetNegRev*2 + nMetNegIrrev)];%fill up the rest with zeros 0307 %vnrn + 100*vnrbm + vnrvm2 == 0 0308 metRows3 = [sparse(nMetNegRev,nRxns + nYBlock + nPosIrrev + nPosRev*6 + nNegIrrev + nNegRev) ... %zeros up to vnrn 0309 nrEye(metNegRev(negRevRxns),:) ... %vnrp 0310 sparse(nMetNegRev, nNegRev*(varsPerNegRev-2) + nEssRev*6 + nMetabolMets*2) ... %zeros up to vnrbm 0311 speye(nMetNegRev)*100 ... %vnrbm 0312 sparse(nMetNegRev, nMetNegRev) ... %zeros up to vnrvm2 0313 speye(nMetNegRev) ... % vnrvm2 0314 sparse(nMetNegRev, nMetNegRev + nMetNegIrrev)];%fill up the rest with zeros 0315 %-0.1*Yi + vnrp + vnrn - vnrvm3 == 0 0316 metRows4 = [sparse(nMetNegRev,nRxns + nPosIrrev + nPosRev + nNegIrrev) ... %zeros up to Yi for neg rev 0317 nrEye(metNegRev(negRevRxns),:)*-0.1 ... %Yi for met neg rev 0318 sparse(nMetNegRev, nPosIrrev + nPosRev*6 + nNegIrrev) ... %zeros up to vnrp 0319 nrEye(metNegRev(negRevRxns),:) ... %vnrp 0320 nrEye(metNegRev(negRevRxns),:) ... %vnrn 0321 sparse(nMetNegRev, nNegRev*(varsPerNegRev-2) + nEssRev*6 + nMetabolMets*2 + nMetNegRev*3) ...%zeros up to vnrvm3 0322 speye(nMetNegRev)*-1 ... % vnrvm3 0323 sparse(nMetNegRev, nMetNegIrrev)];%fill up the rest with zeros 0324 0325 %Then negative irrev, i.e. flux - 0.1*Yi - vnim == 0 0326 niEye = speye(nNegIrrev); 0327 metRows5 = [sEye(metNegIrrev,:) ... %flux 0328 sparse(nMetNegIrrev, nPosIrrev + nPosRev) ... % zeros up to Yi for neg irrev 0329 niEye(metNegIrrev(negIrrevRxns),:)*-0.1 ... %Yi for met neg irrev 0330 sparse(nMetNegIrrev, nNegRev + nPosIrrev + nPosRev*6 + nNegIrrev + nNegRev*varsPerNegRev + nEssRev*6 + nMetabolMets*2 + 4*nMetNegRev) ... %zeros up to vnim 0331 -speye(nMetNegIrrev) ]; %vnim 0332 metRows = [metRows1;metRows2;metRows3;metRows4;metRows5]; 0333 metVarC = [ones(nMetabolMets,1)*-prodWeight;zeros(nMetVars - nMetabolMets,1)]; 0334 %vnrbm is boolean, so between 0 and 1 0335 %vnrvm1 >= 0 0336 %-100 <= vnrvm2 <= inf 0337 % vnrvm3 >= 0 0338 %0 <= vnim <= Inf 0339 metLb = [zeros(nMetabolMets*2 + 2*nMetNegRev,1);ones(nMetNegRev,1)*-100;zeros(nMetNegRev + nMetNegIrrev,1)]; 0340 metUb = [ones(nMetabolMets,1);inf(nMetabolMets,1);ones(nMetNegRev,1);inf(3*nMetNegRev + nMetNegIrrev,1)]; 0341 metVartype = [repmat('C', 1, nMetabolMets*2), ... 0342 repmat('B', 1, nMetNegRev), ... 0343 repmat('C', 1, 3*nMetNegRev + nMetNegIrrev)]; 0344 0345 else 0346 metRows = []; 0347 metVarC = []; 0348 metLb = []; 0349 metUb = []; 0350 metVartype = []; 0351 end 0352 0353 prob.a = [sRow;piRows;prRows1;prRows2;prRows3;prRows4;niRows;nrRows;erRows1;erRows2;erRows3;erRows4;metRows]; 0354 prob.A = prob.a; 0355 prob.b = zeros(size(prob.A,1),1); 0356 prob.c = [zeros(nRxns,1);rxnScores(posIrrevRxns)*-1;rxnScores(posRevRxns)*-1;rxnScores(negIrrevRxns)*-1;rxnScores(negRevRxns)*-1;zeros(nPosIrrev + nPosRev*6 + nNegIrrev + nNegRev*varsPerNegRev + nEssRev*6,1);metVarC]; 0357 prob.lb = [milpModel.lb;zeros(nYBlock + nPosIrrev + 5*nPosRev,1);ones(nPosRev,1)*-100;zeros(nNegIrrev+nNegRev*varsPerNegRev+nEssRev*2,1);ones(nEssRev,1)*forceOnLim;zeros(nEssRev*2,1);ones(nEssRev,1)*-100;metLb]; 0358 prob.lb(essIrrevRxns) = forceOnLimEss(essIrrevRxns);%force flux for the ess irrev rxns - this is the only way these are handled 0359 prob.ub = [milpModel.ub;ones(nYBlock,1);inf(nPosIrrev+nPosRev*2,1);ones(nPosRev,1);inf(3*nPosRev+nNegIrrev+varsPerNegRev*nNegRev+3*nEssRev,1);ones(nEssRev,1);inf(2*nEssRev,1);metUb]; 0360 0361 prob.vartype = [repmat('C', 1, nRxns + nPosIrrev + nPosRev), ... 0362 repmat('B', 1, nNegIrrev + nNegRev), ... 0363 repmat('C', 1, nPosIrrev + 2*nPosRev), ... 0364 repmat('B', 1, nPosRev), ... 0365 repmat('C', 1, 3*nPosRev + nNegIrrev + varsPerNegRev*nNegRev + 3*nEssRev), ... 0366 repmat('B', 1, nEssRev), ... 0367 repmat('C', 1, 2*nEssRev), ... 0368 metVartype]; 0369 0370 onoffVarInd = (1:(nPosIrrev + nPosRev + nNegIrrev + nNegRev)) + nRxns; 0371 onoffPosIrrev = onoffVarInd(1:nPosIrrev); 0372 onoffPosRev = onoffVarInd((1:nPosRev)+nPosIrrev); 0373 onoffNegIrrev = onoffVarInd((1:nNegIrrev)+nPosIrrev+nPosRev); 0374 onoffNegRev = onoffVarInd((1:nNegRev)+nPosIrrev+nPosRev+nNegIrrev); 0375 0376 metVarInd = (1:nMetabolMets) + (length(prob.vartype) - nMetVars); 0377 0378 if allowExcretion 0379 prob.csense = [repmat('L', 1, length(milpModel.mets)), ... 0380 repmat('E', 1, length(prob.b) - length(milpModel.mets))]; 0381 else 0382 prob.csense = '='; 0383 end 0384 0385 params.intTol = 10^-7; %This value is very important. If set too low 0386 %there is a risk that gurobi fails due to numerical 0387 %issues - this happened for Gurobi v. 10.0 with TestModelL. 0388 %On the other hand, it shouldn't be too large 0389 %either. With this value, fluxes of 10-7 can slip 0390 %through, which should be fine. Another option if 0391 %this becomes a problem is to set NumericFocus=2, 0392 %which makes the solver slower but fixes the issue 0393 %with TestModelL. 0394 0395 prob.osense = 1; %minimize 0396 0397 if ~isempty(startVals) 0398 prob.start = startVals; %This doesn't work... 0399 end 0400 0401 res=optimizeProb(prob,params,verbose); 0402 0403 0404 if ~checkSolution(res) 0405 if strcmp(res.origStat, 'TIME_LIMIT') 0406 EM='Time limit reached without finding a solution. Try increasing the TimeLimit parameter.'; 0407 else 0408 EM='The problem is infeasible'; 0409 end 0410 dispEM(EM); 0411 end 0412 0413 %get the on/off vals 0414 onoff = ones(nRxnsWithOnOff,1);%standard is on, i.e. all reactions not included in the problem + ess is on 0415 onoff(posIrrevRxns) = res.full(onoffPosIrrev); 0416 onoff(posRevRxns) = res.full(onoffPosRev); 0417 onoff(negIrrevRxns) = res.full(onoffNegIrrev); 0418 onoff(negRevRxns) = res.full(onoffNegRev); 0419 0420 onoff2 = zeros(nRxnsWithOnOff,1);%standard is off, i.e. all reactions not included in the problem + ess is off 0421 onoff2(posIrrevRxns) = res.full(onoffPosIrrev); 0422 onoff2(posRevRxns) = res.full(onoffPosRev); 0423 onoff2(negIrrevRxns) = res.full(onoffNegIrrev); 0424 onoff2(negRevRxns) = res.full(onoffNegRev); 0425 0426 0427 %investigate a bit - this code could be used for fault finding. We have 0428 %this code commented out, but it is definitely worth testing this if for example 0429 %the solver is changed. 0430 %problematic = onoff < 0.99 & onoff > 0.01 & (rxnScores ~= 0).'; 0431 %if sum(problematic) > 0 0432 % disp('There are reactions that are in between on and off in the MILP. This may be due to the MILP parameter settings, the tolerances may be too large.') 0433 % disp(['No problematic reactions: ', num2str(sum(problematic))]) 0434 %end 0435 0436 %so, it is only positive that are problematic. Likely because the 0.1 is too large. Test to change to 0.01 0437 %unique(onoff(onoff < 0.99 & onoff > 0.01)) 0438 0439 %Get all reactions used in the irreversible model 0440 %The reaction score check is for metabolomics - we add those reactions as 0441 %positive even though the score is 0. And, we don't want them included in the 0442 %results. 0443 deletedRxns=(onoff < 0.5).' & (rxnScores ~= 0).'; 0444 turnedOnRxns=(onoff2 >= 0.5).' & (rxnScores ~= 0).'; 0445 0446 fluxes = res.full(1:nRxns); 0447 0448 %extract the met data 0449 metProduction = logical(round(res.full(metVarInd))); 0450 0451 end