Description of ftINITFillGapsMILP

0001 function [x,I,exitFlag]=ftINITFillGapsMILP(model, toMinimize, params, scores, verbose)
0002 % ftINITFillGapsMILP
0003 %   Returns the minimal set of fluxes that satisfy the model using
0004 %   mixed integer linear programming. This is an optimized variant of the
0005 %   old function "getMinNrFluxes" that is adapted to ftINIT.
0006 %   It does not need to make an irrev model, which takes time. The problem
0007 %   also becomes smaller (fewer integers but larger matrix). Only tested with
0008 %   Gurobi.
0009 %
0010 %    model         a model structure
0011 %   toMinimize    either a cell array of reaction IDs, a logical vector
0012 %                 with the same number of elements as reactions in the model,
0013 %                 of a vector of indexes for the reactions that should be
0014 %                 minimized (optional, default model.rxns)
0015 %   params        parameter structure as used by getMILPParams (optional)
0016 %   scores        vector of weights for the reactions. Negative scores
0017 %                 should not have flux. Positive scores are not possible in this
0018 %                 implementation, and they are changed to max(scores(scores<0)).
0019 %                 Must have the same dimension as toMinimize (find(toMinimize)
0020 %                 if it is a logical vector) (optional, default -1 for all reactions)
0021 %   verbose       if true, the MILP progression will be shown.
0022 %
0023 %   x             the corresponding fluxes for the full model
0024 %   I             the indexes of the reactions in toMinimize that were used
0025 %                 in the solution
0026 %   exitFlag      1: optimal solution found
0027 %                -1: no feasible solution found
0028 %                -2: optimization time out
0029 %
0030 %   NOTE: Uses 1000 mmol/gDW/h as an arbitary large flux. Could possibly
0031 %   cause problems if the fluxes in the model are larger than that.
0032 %
0033 % Usage: [x,I,exitFlag]=getMinNrFluxes(model, toMinimize, params, scores)
0034 
0035 % glpk solver as implemented by COBRA does not work well for MILP.
0036 global CBT_MILP_SOLVER
0037 if strcmp(getpref('RAVEN','solver'),'cobra') && strcmp(CBT_MILP_SOLVER,'glpk')
0038     dispEM('The current solver is set to ''cobra'', while in COBRA the MILP solver has been set to ''glpk''. The COBRA implementation of glpk is not well suitable for solving MILPs. Please install the Gurobi or an alternative MILP solver.',true);
0039 end
0040 
0041 exitFlag=1;
0042 
0043 if nargin<2
0044     toMinimize=model.rxns;
0045 else
0046     if ~iscell(toMinimize)
0047         toMinimize=model.rxns(toMinimize);
0048     end
0049 end
0050 
0051 %For passing parameters to the solver
0052 if nargin<3
0053     params=struct();
0054 end
0055 
0056 if nargin<4
0057     %It says that the default is -1, but that is to fit with other code
0058     scores=ones(numel(toMinimize),1)*1;
0059 else
0060     if numel(scores)~=numel(toMinimize)
0061         EM='The number of scores must be the same as the number of reactions to minimize';
0062         dispEM(EM);
0063     end
0064     
0065     %Change positive scores to have a small negative weight. This is a
0066     %temporary solution.
0067     scores(scores>=0)=max(scores(scores<0));
0068     
0069     %It says that the default is -1, but that is to fit with other code
0070     scores=scores*-1;
0071 end
0072 
0073 %The trick to make this possible using a reversible model is that for
0074 %reversible reactions, we do the following. Set the flux Vi == Vipos - Vineg.
0075 %What happens then is that if the flux is positive, Vipos will have a nonzero value,
0076 %and when the flux is negative, Vineg will have a positive value. In addition
0077 %we can add an arbitrary constant to both Vipos and Vineg. For example, if the
0078 %flux Vi is -1, Vineg can be 4 and Vipos 3. This is however not a problem,
0079 %because we can be sure that Vineg + Vipos >= abs(Vi). Since we are interested
0080 %in forcing the ints to be on when there is a flux, it doesn't matter if we overestimate
0081 %the flux! So, we can simply constrain the boolean Yi to Yi*Maxflux >= Vineg + Vipos.
0082 
0083 %The matrix then becomes as this:
0084 %         S       p n int b     var
0085 %         SSSSSSSS        -
0086 %         SSSSSSSS         -
0087 %         SSSSSSSS          -
0088 %         SSSSSSSS           -
0089 %         SSSSSSSS            -
0090 %         SSSSSSSS             -
0091 %         -       1 -
0092 %               -  1 -
0093 %           -          M         -
0094 %             -         M         -
0095 %                 - - M         -
0096 %                  - -   M         -
0097 % An example with 8 rxns and 6 metabolites. - means -1, M max flux, and S is the S matrix.
0098 % 4 rxns are to be minimized(1,3,5,7) and 1,7 are reversible. The p and n
0099 % are the Vipos and Vineg variables (2 rxns of each). The ints are the Yi for
0100 % the variables that are to be minimized (the rest of the rxns doesn't have any).
0101 % The mets here are the constraints, so right under the S matrix, you have
0102 % Vi == Vipos - Vineg for the reactions 1 and 7 while the two next rows represent
0103 % the non-reversible rxns 3 and 5, where we simply say that yi*M >= Vi. The last
0104 % 2 rows are the reactions yi*M >= Vipos + Vineg. To the right, we first have a
0105 % -I matrix for setting the b constraints, and under that we have rxns that are
0106 % just variables (rxns) between 0 and Inf to complete the constraints mentioned above.
0107 %Ex: yi*M >= Vipos + Vineg is impl. as yi*M - Vipos - Vineg - var == 0, 0 <= var <= Inf.
0108 %
0109 %All rows should be equal to zero, so we don't set the b vector in the problem
0110 %The reactions should be constrained as follows
0111 %S - as given in model.lb and model.ub
0112 %pos and neg - between 0 and inf
0113 %ints - between 0 and 1
0114 %b - as stated in the model.b vector - if this has one column, that is lb and ub (fixed value), if two columns, that is lb and ub
0115 %var - between zero and inf
0116 
0117 [minLog, I]=ismember(model.rxns,toMinimize);
0118 indexes=find(minLog);
0119 revIndexes = find(minLog & (model.rev == 1));
0120 irrevIndexes = find(minLog & (model.rev == 0));
0121 revIndexesInInd = find(model.rev(indexes) == 1);
0122 irrevIndexesInInd = find(model.rev(indexes) == 0);
0123 
0124 %Add binary constraints in the following manner: -  Add one unique
0125 %"metabolite" for each integer reaction as a substrate.
0126 %   These metabolites can have net production
0127 %-  Add reactions for the production of each of those metabolites. The
0128 %   amount produced in one reaction unit must be larger than the largest
0129 %   possible flux in the model (but not too large to avoid bad scaling)
0130 
0131 maxFlux=1000;
0132 %we build the total matrix as blocks: [S pos neg int b var]
0133 
0134 %s block
0135 intArray=speye(numel(model.rxns))*-1;
0136 intArrayRev=intArray(revIndexes,:);
0137 intArrayIrrev = intArray(irrevIndexes,:);
0138 sBlock=[model.S;intArrayRev;intArrayIrrev;sparse(numel(revIndexes), numel(model.rxns))]; %the S matrix and what is below
0139 
0140 %pos/neg blocks
0141 revposorneg1 = sparse(numel(model.mets), numel(revIndexes));
0142 revpos2 = speye(numel(revIndexes));
0143 revneg2 = -revpos2;
0144 revposneg3 = sparse(numel(irrevIndexes), numel(revIndexes));
0145 revposneg4 = revneg2;
0146 posBlock = [revposorneg1;revpos2;revposneg3;revposneg4];
0147 negBlock = [revposorneg1;revneg2;revposneg3;revposneg4];
0148 
0149 %int block
0150 int1 = sparse(numel(model.mets), numel(indexes));
0151 int2 = sparse(numel(revIndexes), numel(indexes));
0152 tmpEye = speye(numel(indexes))*maxFlux;
0153 int3 = tmpEye(irrevIndexesInInd,:);%we here select only the irrev indexes from the indexes
0154 int4 = tmpEye(revIndexesInInd,:);
0155 intBlock = [int1;int2;int3;int4];
0156 
0157 %b block
0158 b1 = speye(numel(model.mets))*-1;
0159 b2 = sparse(numel(indexes) + numel(revIndexes), numel(model.mets));
0160 bBlock = [b1;b2];
0161 
0162 %var block
0163 var1 = sparse(numel(model.mets), numel(indexes));
0164 var2 = sparse(numel(revIndexes), numel(indexes));
0165 tmpEye = speye(numel(indexes))*-1;
0166 var3 = tmpEye(irrevIndexesInInd,:);
0167 var4 = tmpEye(revIndexesInInd,:);
0168 varBlock = [var1;var2;var3;var4];
0169 
0170 prob.A = [sBlock posBlock negBlock intBlock bBlock varBlock];
0171 prob.a = prob.A;%I think this is needed as well
0172 
0173 prob.c=[zeros(numel(model.rxns),1);zeros(numel(revIndexes)*2,1); scores(:);zeros(numel(model.mets) + numel(indexes),1)]; %Minimize the sum of reaction scores for reactions that are on
0174 
0175 %ub and lb, text copied from above
0176 %S - as given in model.lb and model.ub
0177 %pos and neg - between 0 and inf
0178 %ints - between 0 and 1
0179 %b - as stated in the model.b vector - if this has one column, that is lb and ub (fixed value), if two columns, that is lb and ub
0180 %var - between zero and inf
0181 
0182 if size(model.b,2)==2
0183     bub = model.b(:,2);
0184 else
0185     bub = model.b(:,1);
0186 end
0187 
0188 prob.lb = [model.lb;zeros(numel(revIndexes)*2,1);zeros(numel(indexes),1);model.b(:,1);zeros(numel(indexes),1)];
0189 prob.ub = [model.ub;inf(numel(revIndexes)*2,1);ones(numel(indexes),1);bub;inf(numel(indexes),1)];
0190 
0191 prob.b=zeros(size(prob.a,1), 1);
0192 
0193 %Use the output from the linear solution as starting point. Only the values
0194 %for the integer variables will be used, but all are supplied.
0195 intsIndexes = find(prob.c ~= 0);
0196 %The start point is not important (solved quickly anyway), so just skip it.
0197 %prob.sol.int.xx=zeros(numel(prob.c),1);
0198 %prob.sol.int.xx(intsIndexes(sol.x(indexes)>10^-3))=1;%these doesn't work for gurobi anyways...
0199 prob.x0=[];
0200 prob.vartype=repmat('C', 1, size(prob.A,2));
0201 prob.vartype(intsIndexes) = 'B';
0202 prob.csense = repmat('E', 1, size(prob.A,1));
0203 prob.osense=1; %minimize the objective
0204 
0205 %prob=rmfield(prob,{'blx','bux','blc','buc'});
0206 params.intTol = 10^-9; %experment with this value
0207 params.TimeLimit = 300;
0208 params.Seed = 26;%This is weird - although it says "optimal solution found", we can get different results with different
0209                  %values of the objective function, where one is more optimal than the other (pretty big difference...)
0210 %params.CSClientLog = 3;%generates a warning in gurobi, but may be of interest for other solvers
0211 
0212 % Optimize the problem
0213 res = optimizeProb(prob,params,verbose);
0214 isFeasible=checkSolution(res);
0215 
0216 if ~isFeasible
0217     x=[];
0218     I=[];
0219     exitFlag=-1;
0220     return;
0221 end
0222 
0223 x=res.full(1:numel(model.rxns));%the fluxes
0224 I=res.full(intsIndexes) > 10^-3;%The margin for integers in gurobi is 10^-5, not 10^-12 that was previously used! Use 10^-3 to have some margin!
0225 
0226 tmp = res.full(intsIndexes);
0227 sel = (tmp > 10^-12) & (tmp < 0.5);
0228 if sum(sel) > 0
0229     %This may indicate that there is a problem with the tolerances in the solver
0230     disp(['ftINITFillGapsMILP: Some variables meant to be boolean in the MILP have intermediate values. Num vars: ' num2str(sum(sel))])
0231 end
0232 
0233 end
ftINITFillGapsMILP

PURPOSE

SYNOPSIS

DESCRIPTION

CROSS-REFERENCE INFORMATION

SOURCE CODE
