linearMOMA

PURPOSE ^

linearMOMA Performs a linear version of the MOMA (minimization of

SYNOPSIS ^

function [solutionDel,solutionWT,totalFluxDiff,solStatus] =linearMOMA(modelWT,modelDel,osenseStr,minFluxFlag,verbFlag)

DESCRIPTION ^

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

 [solutionDel,solutionWT,totalFluxDiff,solStatus] = 
       linearMOMA(modelWT,modelDel,osenseStr,minFluxFlag,verbFlab)

INPUTS
 modelWT           Wild type model
 modelDel          Deletion strain model

OPTIONAL INPUTS
 osenseStr         Maximize ('max')/minimize ('min') (Default = 'max')
 minFluxFlag       Minimize the absolute value of fluxes in the optimal MOMA
                   solution (Default = false)
 verbFlag          Verbose output (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|
 solStatus         Solution status

 Solves the problem

 min sum|v_wt - v_del|
     S_wt*v_wt = 0
     lb_wt <= v_wt <= ub_wt
     c_wt'*v_wt = f_wt
     S_del*v_del = 0
     lb_del <= v_del <= ub_del

 Here f_wt is the optimal wild type objective value found by FBA

 Notes:

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

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