rescaleModelForINIT

 rescaleModelForINIT

 The idea with this function is to rescale the MILP problem in ftINIT to avoid large differences
 in flux magnitudes between reactions. Such differences cause among other things
 difficulties regarding tolerances for integer variables.
 For now it just scales down all reactions with high stoichiometric coefficients
 There is room for improvement here - the best would be to convert mets such as albumin
 to instead represent 1/100 albumin - that would create much less extreme coefficients.
 This type of improvement is known as scaling in the literature around LPs and MILPs.

 model         the model to be modified (input and output)
 maxStoichVal  all reactions with stoichiometric coefficent higher than this 
               will be scaled down. (opt, default 250)

0001 function model=rescaleModelForINIT(model, maxStoichDiff)
0002 % rescaleModelForINIT
0003 %
0004 % The idea with this function is to rescale the MILP problem in ftINIT to avoid large differences
0005 % in flux magnitudes between reactions. Such differences cause among other things
0006 % difficulties regarding tolerances for integer variables.
0007 % For now it just scales down all reactions with high stoichiometric coefficients
0008 % There is room for improvement here - the best would be to convert mets such as albumin
0009 % to instead represent 1/100 albumin - that would create much less extreme coefficients.
0010 % This type of improvement is known as scaling in the literature around LPs and MILPs.
0011 %
0012 % model         the model to be modified (input and output)
0013 % maxStoichVal  all reactions with stoichiometric coefficent higher than this
0014 %               will be scaled down. (opt, default 250)
0015 
0016 if (nargin < 2)
0017     maxStoichDiff = 25;
0018 end
0019 
0020 %Define maxMinRatio to be the ratio between the largest and smallest
0021 %stoichiometric coefficients in each reaction.
0022 %Scale all rxns with maxMinRatio > maxStoichDiff - just set all coeffs that are
0023 %larger than maxStoichDiff*minVal to that value, and center the mean coeff to 1 for all rxns.
0024 SAbs = abs(model.S);
0025 for i = 1:numel(model.rxns)
0026     tmp = SAbs(:,i);
0027     tmp = tmp(tmp ~= 0);
0028     mn = min(tmp);
0029     %modify matrix
0030     sign = ones(numel(model.mets),1);
0031     sign(model.S(:,i) < 0) = -1;
0032     sel = SAbs(:,i) > maxStoichDiff*mn;
0033     model.S(sel,i) = sign(sel) .* maxStoichDiff*mn;
0034 end
0035 %center around 1
0036 boolMat = model.S ~= 0;
0037 absMat = abs(model.S);
0038 rxnScales = sum(boolMat,1)./sum(absMat,1);
0039 
0040 model.S = model.S .* rxnScales;
0041 
0042 %{
0043 if (nargin < 2)
0044     maxStoichVal = 250;
0045 end
0046 
0047 
0048 %find all reactions containing high stoichiometry coefficients
0049 maxCoeffs = max(model.S, [], 1);
0050 rxnsToCheck = maxCoeffs > maxStoichVal;
0051 
0052 %for debugging:
0053 %constructEquations(model, model.rxns(rxnsToCheck))
0054 
0055 rxnsToCheckInd = find(rxnsToCheck);
0056 
0057 %We just scale all reactions so they don't have a coefficient > maxStoichVal for now
0058 
0059 %do this with a loop, doesn't really matter
0060 for rxnInd = rxnsToCheckInd
0061     %find the lowest index which is not a pointless metabolite, such as H+, H2O etc.
0062     rxnS = model.S(:,rxnInd);
0063     %metInd = find(rxnS ~= 0);
0064     
0065     largestCoeff = max(abs(rxnS));
0066     scaleFactor = maxStoichVal/largestCoeff;
0067     model.S(:,rxnInd) = model.S(:,rxnInd) .* scaleFactor;
0068 end
0069 %}
0070 end