0001 function [solutionDel,solutionWT,totalFluxDiff,solStatus] = ...
0002 linearMOMA(modelWT,modelDel,osenseStr,minFluxFlag,verbFlag)
0003
0004
0005
0006
0007
0008
0009
0010
0011
0012
0013
0014
0015
0016
0017
0018
0019
0020
0021
0022
0023
0024
0025
0026
0027
0028
0029
0030
0031
0032
0033
0034
0035
0036
0037
0038
0039
0040
0041
0042
0043
0044
0045
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
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
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
0090 solutionWT = optimizeCbModel(modelWT,osenseStr);
0091
0092 if (verbFlag)
0093 fprintf('%f seconds\n',solutionWT.time);
0094 end
0095
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
0103
0104
0105
0106
0107
0108
0109 if (solutionWT.stat > 0)
0110
0111
0112
0113
0114
0115
0116
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
0123 b = [zeros(nMets1+nMets2+2*nCommon,1);objValWT];
0124
0125
0126 c = [zeros(nRxns1+nRxns2,1);ones(2*nCommon,1)];
0127
0128
0129
0130 lb = [modelWT.lb;modelDel.lb;zeros(2*nCommon,1)];
0131 ub = [modelWT.ub;modelDel.ub;10000*ones(2*nCommon,1)];
0132
0133
0134
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
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
0173 b = [zeros(nMets1+nMets2+2*nCommon+2*nRxns1+2*nRxns2,1);objValWT;ceil(totalFluxDiff/tol)*tol];
0174
0175
0176 c = [zeros(nRxns1+nRxns2+2*nCommon,1);ones(2*nRxns1+2*nRxns2,1)];
0177
0178
0179
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;