0001 function [solution1,solution2,totalFluxDiff] = optimizeTwoCbModels(model1,model2,osenseStr,minFluxFlag,verbFlag)
0002
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
0048
0049
0050
0051
0052
0053
0054
0055
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 end
0063 else
0064 tol = 1e-6;
0065 end
0066
0067 if (nargin < 3)
0068 osenseStr = 'max';
0069 end
0070 if (nargin < 4)
0071 minFluxFlag = false;
0072 end
0073 if (nargin < 5)
0074 verbFlag = false;
0075 end
0076
0077
0078 if (strcmp(osenseStr,'max'))
0079 osense = -1;
0080 else
0081 osense = +1;
0082 end
0083
0084
0085 [nMets1,nRxns1] = size(model1.S);
0086 [nMets2,nRxns2] = size(model2.S);
0087
0088
0089 commonRxns = ismember(model1.rxns,model2.rxns);
0090 nCommon = sum(commonRxns);
0091 if (nCommon == 0)
0092 error('No common rxns in the models');
0093 end
0094
0095
0096 if (~isfield(model1,'b'))
0097 model1.b = zeros(size(model1.S,1),1);
0098 end
0099 if (~isfield(model2,'b'))
0100 model2.b = zeros(size(model2.S,1),1);
0101 end
0102
0103 csense = [];
0104
0105 if (verbFlag)
0106 fprintf('Solving original FBA problems: %d constraints %d variables ',nMets1+nMets2,nRxns1+nRxns2);
0107 end
0108
0109
0110 FBAsol1 = optimizeCbModel(model1,osenseStr);
0111 FBAsol2 = optimizeCbModel(model2,osenseStr);
0112
0113 if (verbFlag)
0114 fprintf('%f seconds\n',FBAsol1.time+FBAsol2.time);
0115 end
0116
0117
0118 if (FBAsol1.stat > 0 && FBAsol1.stat > 0)
0119 f1 = FBAsol1.f;
0120 f2 = FBAsol2.f;
0121 if (strcmp(osenseStr,'max'))
0122 f1 = floor(f1/tol)*tol;
0123 f2 = floor(f2/tol)*tol;
0124 else
0125 f1 = ceil(f1/tol)*tol;
0126 f2 = ceil(f2/tol)*tol;
0127 end
0128
0129
0130
0131
0132
0133
0134
0135
0136
0137
0138 A = [model1.S sparse(nMets1,nRxns2+2*nCommon);
0139 sparse(nMets2,nRxns1) model2.S sparse(nMets2,2*nCommon);
0140 createDeltaMatchMatrix(model1.rxns,model2.rxns)
0141 model1.c' sparse(1,nRxns2+2*nCommon);
0142 sparse(1,nRxns1) model2.c' sparse(1,2*nCommon);];
0143 c = [zeros(nRxns1+nRxns2,1);ones(2*nCommon,1)];
0144 lb = [model1.lb;model2.lb;zeros(2*nCommon,1)];
0145 ub = [model1.ub;model2.ub,;10000*ones(2*nCommon,1)];
0146 b = [model1.b;model2.b;zeros(2*nCommon,1);f1;f2];
0147 clear csense;
0148 csense(1:(nMets1+nMets2)) = 'E';
0149 csense(end+1:end+2*nCommon) = 'G';
0150 if (strcmp(osenseStr,'max'))
0151 csense(end+1:end+2) = 'G';
0152 else
0153 csense(end+1:end+2) = 'L';
0154 end
0155
0156
0157 if (verbFlag)
0158 fprintf('Minimize difference between solutions: %d constraints %d variables ',size(A,1),size(A,2));
0159 end
0160
0161 [LPproblem.A,LPproblem.b,LPproblem.c,LPproblem.lb,LPproblem.ub,LPproblem.csense,LPproblem.osense] = deal(A,b,c,lb,ub,csense,1);
0162 LPsol = solveCobraLP(LPproblem);
0163
0164 if (verbFlag)
0165 fprintf('%f seconds\n',LPsol.time);
0166 end
0167
0168 if (LPsol.stat > 0)
0169 totalFluxDiff = LPsol.obj;
0170 solution1.f = f1;
0171 solution2.f = f2;
0172 solution1.x = LPsol.full(1:nRxns1);
0173 solution2.x = LPsol.full(nRxns1+1:nRxns1+nRxns2);
0174 else
0175 totalFluxDiff = [];
0176 solution1.f = [];
0177 solution2.f = [];
0178 solution1.x = [];
0179 solution2.x = [];
0180 end
0181
0182 if (LPsol.stat > 0 & minFluxFlag)
0183 A = [model1.S sparse(nMets1,nRxns2+2*nCommon+2*nRxns1+2*nRxns2);
0184 sparse(nMets2,nRxns1) model2.S sparse(nMets2,2*nCommon+2*nRxns1+2*nRxns2);];
0185 A = [A;
0186 createDeltaMatchMatrix(model1.rxns,model2.rxns) sparse(2*nCommon,2*nRxns1+2*nRxns2)];
0187 A = [A;
0188 speye(nRxns1,nRxns1) sparse(nRxns1,nRxns2) sparse(nRxns1,2*nCommon) speye(nRxns1,nRxns1) sparse(nRxns1,nRxns1+2*nRxns2);
0189 -speye(nRxns1,nRxns1) sparse(nRxns1,nRxns2) sparse(nRxns1,2*nCommon) sparse(nRxns1,nRxns1) speye(nRxns1,nRxns1) sparse(nRxns1,2*nRxns2);
0190 sparse(nRxns2,nRxns1) speye(nRxns2,nRxns2) sparse(nRxns2,2*nCommon) sparse(nRxns2,2*nRxns1) speye(nRxns2,nRxns2) sparse(nRxns2,nRxns2);
0191 sparse(nRxns2,nRxns1) -speye(nRxns2,nRxns2) sparse(nRxns2,2*nCommon) sparse(nRxns2,2*nRxns1) sparse(nRxns2,nRxns2) speye(nRxns2,nRxns2);];
0192 A = [A;
0193 model1.c' sparse(1,nRxns2+2*nCommon+2*nRxns1+2*nRxns2);
0194 sparse(1,nRxns1) model2.c' sparse(1,2*nCommon+2*nRxns1+2*nRxns2);
0195 sparse(1,nRxns1+nRxns2) ones(1,2*nCommon) sparse(1,2*nRxns1+2*nRxns2)];
0196
0197 b = [zeros(nMets1+nMets2+2*nCommon+2*nRxns1+2*nRxns2,1);f1;f2;ceil(totalFluxDiff/tol)*tol];
0198
0199
0200 c = [zeros(nRxns1+nRxns2+2*nCommon,1);ones(2*nRxns1+2*nRxns2,1)];
0201
0202
0203
0204 lb = [model1.lb;model2.lb;zeros(2*nCommon+2*nRxns1+2*nRxns2,1)];
0205 ub = [model1.ub;model2.ub;10000*ones(2*nCommon+2*nRxns1+2*nRxns2,1)];
0206 csense(1:(nMets1+nMets2)) = 'E';
0207 csense((nMets1+nMets2)+1:(nMets1+nMets2+2*nCommon+2*nRxns1+2*nRxns2)) = 'G';
0208 if (strcmp(osenseStr,'max'))
0209 csense(end+1:end+2) = 'G';
0210 else
0211 csense(end+1:end+2) = 'L';
0212 end
0213 csense(end+1) = 'L';
0214
0215 if (verbFlag)
0216 fprintf('Minimizing flux distribution norms: %d constraints %d variables ',size(A,1),size(A,2));
0217 end
0218
0219 [LPproblem.A,LPproblem.b,LPproblem.c,LPproblem.lb,LPproblem.ub,LPproblem.csense,LPproblem.osense] = deal(A,b,c,lb,ub,csense,1);
0220 LPsol = solveCobraLP(LPproblem);
0221
0222 if (verbFlag)
0223 fprintf('%f seconds\n',LPsol.time);
0224 end
0225
0226 if (LPsol.stat > 0)
0227 solution1.x = LPsol.full(1:nRxns1);
0228 solution2.x = LPsol.full(nRxns1+1:nRxns1+nRxns2);
0229 end
0230
0231 end
0232 end
0233
0234 solution1.stat = LPsol.stat;
0235 solution2.stat = LPsol.stat;
0236