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