MOMA

PURPOSE ^

MOMA Performs a quadratic version of the MOMA (minimization of

SYNOPSIS ^

function [solutionDel,solutionWT,totalFluxDiff,solStatus] =MOMA(modelWT,modelDel,osenseStr,verbFlag,minNormFlag)

DESCRIPTION ^

MOMA Performs a quadratic version of the MOMA (minimization of
metabolic adjustment) approach 

 [solutionDel,solutionWT,totalFluxDiff,solStatus] = MOMA(modelWT,modelDel,osenseStr,verbFlag,minNormFlag)

INPUTS
 modelWT           Wild type model
 modelDel          Deletion strain model

OPTIONAL INPUTS
 osenseStr         Maximize ('max')/minimize ('min') (Default = 'max')
 verbFlag          Verbose output (Default = false)
 minNormFlag       Work with minimum 1-norm flux distribution for the FBA
                   problem (Default = false)
 
OUTPUTS
 solutionDel       Deletion solution structure
 solutionWT        Wild-type solution structure
 totalFluxDiff     Value of the linear MOMA objective, i.e.
                   sum(v_wt-v_del)^2
 solStatus         Solution status

 Solves two different types of MOMA problems:

 1) MOMA that avoids problems with alternative optima (this is the
 default)

    First solve:
    
    max c_wt'*v_wt0
     lb_wt <= v_wt0 <= ub_wt
     S_wt*v_wt0 = 0
    
    Then solve: 

    min sum(v_wt - v_del)^2
     S_wt*v_wt = 0
     S_del*v_del = 0
     lb_wt <= v_wt <= ub_wt
     lb_del <= v_del <= ub_del
     c_wt'*v_wt = f_wt

   Here f_wt is the optimal wild type objective value found by FBA in the
   first problem. Note that the FBA solution v_wt0 is not used in the second
   problem. This formulation avoids any problems with alternative optima

 2) MOMA that uses a minimum 1-norm wild type FBA solution (this approach
 is used if minNormFlag = true)

    First solve

    max c_wt'*v_wt0
     lb_wt <= v_wt0 <= ub_wt
     S_wt*v_wt0 = 0

    Then solve

    min |v_wt|
     S_wt*v_wt = b_wt
     c_wt'*v_wt = f_wt
     lb_wt <= v_wt <= ub_wt
    
    Here f_wt is the objective value obtained in the 1st optimization.

    Finally solve:

    min sum(v_wt - v_del)^2
     S_del*v_del = 0
     lb_del <= v_del <= ub_del

 Notes:

 1) These formulation allows for selecting for more appropriate
 optimal wild type FBA solutions as the starting point as opposed to
 picking an arbitrary starting point (original MOMA implementation).

 2) The reaction sets in the two models do not have to be equal as long as
 there is at least one reaction in common

 Markus Herrgard 11/7/06

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SOURCE CODE ^

0001 function [solutionDel,solutionWT,totalFluxDiff,solStatus] = ...
0002     MOMA(modelWT,modelDel,osenseStr,verbFlag,minNormFlag)
0003 %MOMA Performs a quadratic version of the MOMA (minimization of
0004 %metabolic adjustment) approach
0005 %
0006 % [solutionDel,solutionWT,totalFluxDiff,solStatus] = MOMA(modelWT,modelDel,osenseStr,verbFlag,minNormFlag)
0007 %
0008 %INPUTS
0009 % modelWT           Wild type model
0010 % modelDel          Deletion strain model
0011 %
0012 %OPTIONAL INPUTS
0013 % osenseStr         Maximize ('max')/minimize ('min') (Default = 'max')
0014 % verbFlag          Verbose output (Default = false)
0015 % minNormFlag       Work with minimum 1-norm flux distribution for the FBA
0016 %                   problem (Default = false)
0017 %
0018 %OUTPUTS
0019 % solutionDel       Deletion solution structure
0020 % solutionWT        Wild-type solution structure
0021 % totalFluxDiff     Value of the linear MOMA objective, i.e.
0022 %                   sum(v_wt-v_del)^2
0023 % solStatus         Solution status
0024 %
0025 % Solves two different types of MOMA problems:
0026 %
0027 % 1) MOMA that avoids problems with alternative optima (this is the
0028 % default)
0029 %
0030 %    First solve:
0031 %
0032 %    max c_wt'*v_wt0
0033 %     lb_wt <= v_wt0 <= ub_wt
0034 %     S_wt*v_wt0 = 0
0035 %
0036 %    Then solve:
0037 %
0038 %    min sum(v_wt - v_del)^2
0039 %     S_wt*v_wt = 0
0040 %     S_del*v_del = 0
0041 %     lb_wt <= v_wt <= ub_wt
0042 %     lb_del <= v_del <= ub_del
0043 %     c_wt'*v_wt = f_wt
0044 %
0045 %   Here f_wt is the optimal wild type objective value found by FBA in the
0046 %   first problem. Note that the FBA solution v_wt0 is not used in the second
0047 %   problem. This formulation avoids any problems with alternative optima
0048 %
0049 % 2) MOMA that uses a minimum 1-norm wild type FBA solution (this approach
0050 % is used if minNormFlag = true)
0051 %
0052 %    First solve
0053 %
0054 %    max c_wt'*v_wt0
0055 %     lb_wt <= v_wt0 <= ub_wt
0056 %     S_wt*v_wt0 = 0
0057 %
0058 %    Then solve
0059 %
0060 %    min |v_wt|
0061 %     S_wt*v_wt = b_wt
0062 %     c_wt'*v_wt = f_wt
0063 %     lb_wt <= v_wt <= ub_wt
0064 %
0065 %    Here f_wt is the objective value obtained in the 1st optimization.
0066 %
0067 %    Finally solve:
0068 %
0069 %    min sum(v_wt - v_del)^2
0070 %     S_del*v_del = 0
0071 %     lb_del <= v_del <= ub_del
0072 %
0073 % Notes:
0074 %
0075 % 1) These formulation allows for selecting for more appropriate
0076 % optimal wild type FBA solutions as the starting point as opposed to
0077 % picking an arbitrary starting point (original MOMA implementation).
0078 %
0079 % 2) The reaction sets in the two models do not have to be equal as long as
0080 % there is at least one reaction in common
0081 %
0082 % Markus Herrgard 11/7/06
0083 
0084 if (nargin <3 || isempty(osenseStr))
0085     osenseStr = 'max';
0086 end
0087 if (nargin < 4)
0088     verbFlag = false;
0089 end
0090 if (nargin < 5)
0091     minNormFlag = false;
0092 end
0093 
0094 % LP solution tolerance
0095 global CBT_LP_PARAMS
0096 if (exist('CBT_LP_PARAMS', 'var'))
0097     if isfield(CBT_LP_PARAMS, 'objTol')
0098         tol = CBT_LP_PARAMS.objTol;
0099     else
0100         tol = 1e-6;
0101     end
0102 else
0103     tol = 1e-6;
0104 end
0105 
0106 [nMets1,nRxns1] = size(modelWT.S);
0107 [nMets2,nRxns2] = size(modelDel.S);
0108 
0109 % Match model reaction sets
0110 selCommon1 = ismember(modelWT.rxns,modelDel.rxns);
0111 nCommon = sum(selCommon1);
0112 if (nCommon == 0)
0113     error('No common rxns in the models');
0114 end
0115 
0116 solutionWT.f = [];
0117 solutionWT.x = [];
0118 solutionWT.stat = -1;
0119 solutionDel.f = [];
0120 solutionDel.x = [];
0121 solutionDel.stat = -1;
0122 
0123 if (verbFlag)
0124     fprintf('Solving wild type FBA: %d constraints %d variables ',nMets1,nRxns1);
0125 end
0126 % Solve wt problem
0127 if minNormFlag
0128     solutionWT = optimizeCbModel(modelWT,osenseStr,true);
0129 else
0130     solutionWT = optimizeCbModel(modelWT,osenseStr);
0131 end
0132 
0133 if (verbFlag)
0134     fprintf('%f seconds\n',solutionWT.time);
0135 end
0136 % Round off solution to avoid numerical problems
0137 
0138 if (strcmp(osenseStr,'max'))
0139     objValWT = floor(solutionWT.f/tol)*tol;
0140 else
0141     objValWT = ceil(solutionWT.f/tol)*tol;
0142 end
0143 
0144 % Variables in the following problem are
0145 % x = [v1;v2;delta]
0146 % where v1 = wild type flux vector
0147 %       v2 = deletion strain flux vector
0148 %       delta = v1 - v2
0149 
0150 if (solutionWT.stat > 0)
0151     
0152     if minNormFlag
0153 
0154         b = zeros(nMets2,1);
0155         A = modelDel.S;
0156         c = -2*solutionWT.x;
0157         F = 2*eye(nRxns2);
0158         lb = modelDel.lb;
0159         ub = modelDel.ub;
0160         csense(1:nMets2) = 'E';
0161 
0162     else
0163 
0164         % Construct the LHS matrix
0165         % Rows:
0166         % 1: Swt*v1 = 0 for the wild type
0167         % 2: Sdel*v2 = 0 for the deletion strain
0168         % 5: c'v1 = f1 (wild type)
0169         deltaMat = createDeltaMatchMatrix(modelWT.rxns,modelDel.rxns);
0170         deltaMat = deltaMat(1:nCommon,1:(nRxns1+nRxns2+nCommon));
0171         A = [modelWT.S sparse(nMets1,nRxns2+nCommon);
0172             sparse(nMets2,nRxns1) modelDel.S sparse(nMets2,nCommon);
0173             deltaMat;
0174             modelWT.c' sparse(1,nRxns2+nCommon)];
0175 
0176         % Construct the RHS vector
0177         b = [zeros(nMets1+nMets2+nCommon,1);objValWT];
0178 
0179         % Linear objective = 0
0180         c = zeros(nRxns1+nRxns2+nCommon,1);
0181 
0182         % Construct the ub/lb
0183         % delta [-10000 10000]
0184         lb = [modelWT.lb;modelDel.lb;-10000*ones(nCommon,1)];
0185         ub = [modelWT.ub;modelDel.ub;10000*ones(nCommon,1)];
0186 
0187         % Construct the constraint direction vector (G for delta's, E for
0188         % everything else)
0189         csense(1:(nMets1+nMets2+nCommon)) = 'E';
0190         if (strcmp(osenseStr,'max'))
0191             csense(end+1) = 'G';
0192         else
0193             csense(end+1) = 'L';
0194         end
0195 
0196         % F matrix
0197         F = [sparse(nRxns1+nRxns2,nRxns1+nRxns2+nCommon);
0198             sparse(nCommon,nRxns1+nRxns2) 2*eye(nCommon)];
0199 
0200     end
0201     
0202     if (verbFlag)
0203         fprintf('Solving MOMA: %d constraints %d variables ',size(A,1),size(A,2));
0204     end
0205     
0206     % Solve the linearMOMA problem
0207     [QPproblem.A,QPproblem.b,QPproblem.F,QPproblem.c,QPproblem.lb,QPproblem.ub,QPproblem.csense,QPproblem.osense] = deal(A,b,F,c,lb,ub,csense,1);
0208     %QPsolution = solveCobraQP(QPproblem,[],verbFlag-1);
0209     QPsolution = solveCobraQP(QPproblem, 'printLevel', verbFlag-1);
0210 
0211     if (verbFlag)
0212         fprintf('%f seconds\n',QPsolution.time);
0213     end
0214 
0215     % Get the solution(s)
0216     if (QPsolution.stat > 0)
0217         if minNormFlag
0218             solutionDel.x = QPsolution.full;
0219         else
0220             solutionDel.x = QPsolution.full((nRxns1+1):(nRxns1+nRxns2));
0221             solutionWT.x = QPsolution.full(1:nRxns1);
0222         end
0223         solutionDel.f = sum(modelDel.c.*solutionDel.x);
0224         totalFluxDiff = sum((solutionWT.x-solutionDel.x).^2);
0225     end
0226     solutionDel.stat = QPsolution.stat;
0227     solStatus = QPsolution.stat;
0228     solutionDel.solver = QPsolution.solver;
0229     solutionDel.time = QPsolution.time;
0230     
0231 else
0232     warning('Wild type FBA problem is infeasible or unconstrained');
0233 end
0234 
0235

Generated on Thu 21-Jun-2012 15:39:23 by m2html © 2003