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
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