optimizeTwoCbModels

PURPOSE ^

optimizeTwoCbModels Simultaneously solve two flux balance problems and

SYNOPSIS ^

function [solution1,solution2,totalFluxDiff] = optimizeTwoCbModels(model1,model2,osenseStr,minFluxFlag,verbFlag)

DESCRIPTION ^

optimizeTwoCbModels Simultaneously solve two flux balance problems and
minimize the difference between the two solutions

 [solution1,solution2,totalFluxDiff] =
   optimizeTwoCbModels(model1,model2,osenseStr,minFluxFlag,verbFlag)

INPUTS
 model1        The first COBRA model
 model2        The second COBRA model
 model (the following fields are requires - others can be supplied)
   S             Stoichiometric matrix
   b             Right hand side = 0
   c             Objective coefficients
   lb            Lower bounds
   ub            Upper bounds

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
 solution1     Solution for the 1st model
 solution2     Solution for the 2nd model
 totalFluxDiff 1-norm of the difference between the flux vectors sum|v1-v2| 

 solution
   f         Objective value
   x         Primal (flux vector)


 First solves two separate FBA problems:
                                 f1 = max/min c1'v1
                                 subject to S1*v1 = b1
                                            lb1 <= v1 <= ub1
                                 f2 = max/min c2'v2
                                 subject to S2*v2 = b2
                                            lb2 <= v2 <= ub2

 Then solves the following LP to obtain the two flux vectors with the
 smallest possible 1-norm difference between them

                                 min |v1-v2|
                                   s.t. S1*v1 = b1
                                        c1'v1 = f1
                                        lb1 <= v1 <= ub1
                                        S2*v2 = b2
                                        c2'v2 = f2
                                        lb2 <= v2 <= ub2
 
 Finally optionally minimizes the 1-norm of the flux vectors

 Markus Herrgard 1/4/07

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SOURCE CODE ^

0001 function [solution1,solution2,totalFluxDiff] = optimizeTwoCbModels(model1,model2,osenseStr,minFluxFlag,verbFlag)
0002 %optimizeTwoCbModels Simultaneously solve two flux balance problems and
0003 %minimize the difference between the two solutions
0004 %
0005 % [solution1,solution2,totalFluxDiff] =
0006 %   optimizeTwoCbModels(model1,model2,osenseStr,minFluxFlag,verbFlag)
0007 %
0008 %INPUTS
0009 % model1        The first COBRA model
0010 % model2        The second COBRA model
0011 % model (the following fields are requires - others can be supplied)
0012 %   S             Stoichiometric matrix
0013 %   b             Right hand side = 0
0014 %   c             Objective coefficients
0015 %   lb            Lower bounds
0016 %   ub            Upper bounds
0017 %
0018 %OPTIONAL INPUTS
0019 % osenseStr     Maximize ('max')/minimize ('min') (Default = 'max')
0020 % minFluxFlag   Minimize the absolute value of fluxes in the optimal MOMA
0021 %               solution (Default = false)
0022 % verbFlag      Verbose output (Default = false)
0023 %
0024 %OUTPUTS
0025 % solution1     Solution for the 1st model
0026 % solution2     Solution for the 2nd model
0027 % totalFluxDiff 1-norm of the difference between the flux vectors sum|v1-v2|
0028 %
0029 % solution
0030 %   f         Objective value
0031 %   x         Primal (flux vector)
0032 %
0033 %
0034 % First solves two separate FBA problems:
0035 %                                 f1 = max/min c1'v1
0036 %                                 subject to S1*v1 = b1
0037 %                                            lb1 <= v1 <= ub1
0038 %                                 f2 = max/min c2'v2
0039 %                                 subject to S2*v2 = b2
0040 %                                            lb2 <= v2 <= ub2
0041 %
0042 % Then solves the following LP to obtain the two flux vectors with the
0043 % smallest possible 1-norm difference between them
0044 %
0045 %                                 min |v1-v2|
0046 %                                   s.t. S1*v1 = b1
0047 %                                        c1'v1 = f1
0048 %                                        lb1 <= v1 <= ub1
0049 %                                        S2*v2 = b2
0050 %                                        c2'v2 = f2
0051 %                                        lb2 <= v2 <= ub2
0052 %
0053 % Finally optionally minimizes the 1-norm of the flux vectors
0054 %
0055 % Markus Herrgard 1/4/07
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     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 % Figure out objective sense
0078 if (strcmp(osenseStr,'max'))
0079     osense = -1;
0080 else
0081     osense = +1;
0082 end
0083 
0084 % Find model dimensionalities
0085 [nMets1,nRxns1] = size(model1.S);
0086 [nMets2,nRxns2] = size(model2.S);
0087 
0088 % Match model reaction sets
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 % Fill in the RHS vector if not provided
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 % Solve original FBA problems
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 % Minimize the difference between flux solutions
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     % Set up the optimization problem
0129     % min sum(delta+ + delta-)
0130     % 1: S1*v1 = 0
0131     % 2: S2*v2 = 0
0132     % 3: delta+ >= v1-v2
0133     % 4: delta- >= v2-v1
0134     % 5: c1'v1 >= f1 (optimal value of objective)
0135     % 6: c2'v2 >= f2
0136     %
0137     % delta+,delta- >= 0
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     % Re-solve the problem
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         % Construct the RHS vector
0197         b = [zeros(nMets1+nMets2+2*nCommon+2*nRxns1+2*nRxns2,1);f1;f2;ceil(totalFluxDiff/tol)*tol];
0198 
0199         % Construct the objective (sum of all delta+ and delta-)
0200         c = [zeros(nRxns1+nRxns2+2*nCommon,1);ones(2*nRxns1+2*nRxns2,1)];
0201 
0202         % Construct the ub/lb
0203         % delta+ and delta- are in [0 10000]
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

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