Home > RAVEN > qMOMA.m

qMOMA

PURPOSE ^

qMOMA

SYNOPSIS ^

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

DESCRIPTION ^

 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)

   Rasmus Agren, 2010-12-16

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SOURCE CODE ^

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 %   Rasmus Agren, 2010-12-16
0022 %
0023 
0024 if nargin<3
0025     fluxMinWeight=1;
0026 end
0027 
0028 %Match the reactions and metabolites in the small model (modelA) to the
0029 %larger model
0030 [rxnExists,mapRxns]=ismember(modelA.rxns,modelB.rxns);
0031 
0032 %Check that the smaller model is a subset of the larger one
0033 if any(~rxnExists)
0034     throw(MException('','All reactions in the test model must exist in the reference model'));
0035 end
0036 
0037 %In order to make the calculations a little easier to formulate I reshape
0038 %modelA.S so that it is the same dimension as modelB.S
0039 S=modelA.S;
0040 lb=modelA.lb;
0041 ub=modelA.ub;
0042 c=modelA.c;
0043 modelA.S=sparse(numel(modelA.mets),numel(modelB.rxns));
0044 modelA.lb=zeros(numel(modelB.rxns),1);
0045 modelA.ub=modelA.lb;
0046 modelA.c=modelA.lb;
0047 modelA.S(:,mapRxns)=S;
0048 modelA.lb(mapRxns)=lb;
0049 modelA.ub(mapRxns)=ub;
0050 modelA.c(mapRxns)=c;
0051 
0052 %Create the H matrix for the quadratic problem. Note that there are no
0053 %linear terms in the equation
0054 
0055 %Equation to solve is:(xA1 - xB1)^2 + (xA2 - xB2)^2 ....
0056 %Is equal to xA1^2 - xA1*xB1 - xB1*xA1 + xB1^2 + xA2^2 - xA2*xB2 - xB2*xA2 + xB2^2...
0057 
0058 %For three fluxes this would give the following H matrix (first three rows
0059 %for A and last three rows for B)
0060 
0061 %   A1 A2 A3 B1 B2 B3
0062 %A1 1  0  0  -1  0  0
0063 %A2 0  1  0  0  -1  0
0064 %A3 0  0  1  0  0  -1
0065 %B1 -1  0  0  1  0  0
0066 %B2 0  -1  0  0  1  0
0067 %B3 0  0  -1  0  0  1
0068 
0069 %Create stoichiometric matrix and bounds that contain both models
0070 fullS=[modelA.S sparse(size(modelA.S,1),size(modelA.S,2));sparse(size(modelB.S,1),size(modelB.S,2)) modelB.S];
0071 fullLB=[modelA.lb;modelB.lb];
0072 fullUB=[modelA.ub;modelB.ub];
0073 fullB=zeros(size(modelA.S,1)+size(modelB.S,1),1);
0074 
0075 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];
0076 
0077 [x,fval,flag,output]=quadprog(H,zeros(size(H,1),1),[],[],fullS,fullB,fullLB,fullUB);
0078 
0079 if any(x)
0080     fluxA=x(1:numel(modelB.rxns));
0081     fluxA=fluxA(mapRxns);
0082     fluxB=x(numel(modelB.rxns)+1:end);
0083     flag=1; %Since it never converges good enough
0084 else
0085     fluxA=zeros(numel(modelA.rxns),1); %Still the old number
0086     fluxB=zeros(numel(modelB.rxns),1);
0087     flag=-1;
0088 end
0089 end

Generated on Tue 16-Jul-2013 21:50:02 by m2html © 2005