qMOMA

 qMOMA
   Uses quadratic programming to minimize the sum((fluxAi - fluxBi)^2)

   modelA        a model structure for the test case. This model must be a
                 subset of modelB (no reactions that are not in modelB)
   modelB        a model structure for the reference case
   fluxMinWeight a double >=1 that determines whether minimization of the
                 sum of fluxes should also be taken into account in the
                 optimization. A value of 2.0 means that sum(fluxAi)^2 +
                 sum(fluxBi)^2 has equal weight as sum((fluxAi - fluxBi)^2).
                 Values of around 1.01 should be enough to get rid of loops
                 (opt, default 1)

   fluxA         the resulting flux distribution in the test model
   fluxB         the resulting flux distribution in the reference model
   flag          1 if the optimization terminated successfully

   Usage: [fluxA,fluxB, flag]=qMOMA(modelA,modelB,fluxMinWeight)

0001 function [fluxA,fluxB, flag]=qMOMA(modelA,modelB,fluxMinWeight)
0002 % qMOMA
0003 %   Uses quadratic programming to minimize the sum((fluxAi - fluxBi)^2)
0004 %
0005 %   modelA        a model structure for the test case. This model must be a
0006 %                 subset of modelB (no reactions that are not in modelB)
0007 %   modelB        a model structure for the reference case
0008 %   fluxMinWeight a double >=1 that determines whether minimization of the
0009 %                 sum of fluxes should also be taken into account in the
0010 %                 optimization. A value of 2.0 means that sum(fluxAi)^2 +
0011 %                 sum(fluxBi)^2 has equal weight as sum((fluxAi - fluxBi)^2).
0012 %                 Values of around 1.01 should be enough to get rid of loops
0013 %                 (opt, default 1)
0014 %
0015 %   fluxA         the resulting flux distribution in the test model
0016 %   fluxB         the resulting flux distribution in the reference model
0017 %   flag          1 if the optimization terminated successfully
0018 %
0019 %   Usage: [fluxA,fluxB, flag]=qMOMA(modelA,modelB,fluxMinWeight)
0020 
0021 if nargin<3
0022     fluxMinWeight=1;
0023 end
0024 
0025 %Match the reactions and metabolites in the small model (modelA) to the
0026 %larger model
0027 [rxnExists,mapRxns]=ismember(modelA.rxns,modelB.rxns);
0028 
0029 %Check that the smaller model is a subset of the larger one
0030 if any(~rxnExists)
0031     dispEM('All reactions in the test model must exist in the reference model');
0032 end
0033 
0034 %In order to make the calculations a little easier to formulate I reshape
0035 %modelA.S so that it is the same dimension as modelB.S
0036 S=modelA.S;
0037 lb=modelA.lb;
0038 ub=modelA.ub;
0039 c=modelA.c;
0040 modelA.S=sparse(numel(modelA.mets),numel(modelB.rxns));
0041 modelA.lb=zeros(numel(modelB.rxns),1);
0042 modelA.ub=modelA.lb;
0043 modelA.c=modelA.lb;
0044 modelA.S(:,mapRxns)=S;
0045 modelA.lb(mapRxns)=lb;
0046 modelA.ub(mapRxns)=ub;
0047 modelA.c(mapRxns)=c;
0048 
0049 %Create the H matrix for the quadratic problem. Note that there are no
0050 %linear terms in the equation
0051 
0052 %Equation to solve is:(xA1 - xB1)^2 + (xA2 - xB2)^2 ....
0053 %Is equal to xA1^2 - xA1*xB1 - xB1*xA1 + xB1^2 + xA2^2 - xA2*xB2 - xB2*xA2 + xB2^2...
0054 
0055 %For three fluxes this would give the following H matrix (first three rows
0056 %for A and last three rows for B)
0057 
0058 %   A1 A2 A3 B1 B2 B3
0059 %A1 1  0  0  -1  0  0
0060 %A2 0  1  0  0  -1  0
0061 %A3 0  0  1  0  0  -1
0062 %B1 -1  0  0  1  0  0
0063 %B2 0  -1  0  0  1  0
0064 %B3 0  0  -1  0  0  1
0065 
0066 %Create stoichiometric matrix and bounds that contain both models
0067 fullS=[modelA.S sparse(size(modelA.S,1),size(modelA.S,2));sparse(size(modelB.S,1),size(modelB.S,2)) modelB.S];
0068 fullLB=[modelA.lb;modelB.lb];
0069 fullUB=[modelA.ub;modelB.ub];
0070 fullB=zeros(size(modelA.S,1)+size(modelB.S,1),1);
0071 
0072 H=[eye(size(fullS,2)/2)*fluxMinWeight eye(size(fullS,2)/2)*-1;eye(size(fullS,2)/2)*-1 eye(size(fullS,2)/2)*fluxMinWeight];
0073 
0074 x=quadprog(H,zeros(size(H,1),1),[],[],fullS,fullB,fullLB,fullUB);
0075 
0076 if any(x)
0077     fluxA=x(1:numel(modelB.rxns));
0078     fluxA=fluxA(mapRxns);
0079     fluxB=x(numel(modelB.rxns)+1:end);
0080     flag=1; %Since it never converges good enough
0081 else
0082     fluxA=zeros(numel(modelA.rxns),1); %Still the old number
0083     fluxB=zeros(numel(modelB.rxns),1);
0084     flag=-1;
0085 end
0086 end