enumerateOptimalSolutions

PURPOSE ^

enumerateOptimalSolution returns a set of optimal flux distributions

SYNOPSIS ^

function [solution] = enumerateOptimalSolutions(model)

DESCRIPTION ^

enumerateOptimalSolution returns a set of optimal flux distributions
spanning the optimal set

[solution] enumerateOptimalSolution(model) 

INPUT
 model         COBRA model structure

OUTPUT
 solution      solution strcture
   fluxes      Flux distribution for each iteration
   nonzero     Boolean matrix denoting which fluxes are nonzero for each
               iteration

 Author:  Jan Schellenberger, August 2008
 Based on code by Jennie Reed 
 Reed, J.L. and Palsson, B.O., "Genome-scale in silico models of ''E. coli'' have multiple equivalent phenotypic states: assessment of correlated reaction subsets that comprise network states" , Genome Research, 14:1797-1805(2004).

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SOURCE CODE ^

0001 function [solution] = enumerateOptimalSolutions(model)
0002 %enumerateOptimalSolution returns a set of optimal flux distributions
0003 %spanning the optimal set
0004 %
0005 %[solution] enumerateOptimalSolution(model)
0006 %
0007 %INPUT
0008 % model         COBRA model structure
0009 %
0010 %OUTPUT
0011 % solution      solution strcture
0012 %   fluxes      Flux distribution for each iteration
0013 %   nonzero     Boolean matrix denoting which fluxes are nonzero for each
0014 %               iteration
0015 %
0016 % Author:  Jan Schellenberger, August 2008
0017 % Based on code by Jennie Reed
0018 % Reed, J.L. and Palsson, B.O., "Genome-scale in silico models of ''E. coli'' have multiple equivalent phenotypic states: assessment of correlated reaction subsets that comprise network states" , Genome Research, 14:1797-1805(2004).
0019 
0020 
0021 [m,n] = size(model.S);
0022 
0023 solution.fluxes = zeros(n,0);
0024 solution.nonzero = zeros(n,0);
0025 
0026 sol = optimizeCbModel(model);
0027 maxObjective = sol.f;
0028 
0029 
0030 prevNZ = abs(sol.x) > .0001;
0031 NZ = prevNZ;
0032 solution.fluxes = sol.x;
0033 solution.nonzero = prevNZ;
0034 
0035 while 1
0036     % variables:
0037     %    v's (n), y's (n) w's (n)  3n total variables
0038     
0039     % constriants:
0040     %    m mass balance constraints
0041     A = [model.S, zeros(m,2*n)];
0042     b = zeros(m,1);
0043     csense = '';
0044     for i = 1:m
0045         csense(end+1) = 'E';
0046     end
0047     % constrain UB fluxes w/ integer constraints
0048     A = [A; 
0049         [eye(n,2*n), -diag(model.ub)] ];
0050     b = [b;
0051         zeros(n,1)];
0052     for i = 1:n
0053         csense(end+1) = 'L';
0054     end    
0055     % constrain LB fluxes w/ integer constraints
0056     A = [A; 
0057         eye(n,2*n), -diag(model.lb) ];
0058     b = [b;
0059         zeros(n,1)];
0060     for i = 1:n
0061         csense(end+1) = 'G';
0062     end    
0063     
0064 
0065     % constrain w+y <=1
0066     A = [A; 
0067         zeros(n,n), eye(n,n), eye(n,n) ];
0068     b = [b;
0069         ones(n,1)];
0070     for i = 1:n
0071         csense(end+1) = 'L';
0072     end    
0073 
0074     % constrain with previous zero results
0075     A = [A; 
0076         zeros(1,n), prevNZ', zeros(1,n) ];
0077     b = [b;
0078         1];
0079     csense(end+1) = 'G';
0080 
0081     % constrain with previous results (altbases)
0082     for i = 1:size(NZ,2)
0083         A = [A;
0084             [zeros(1,n), zeros(1,n) NZ(:,i)']];
0085         b(end+1) = sum(NZ(:,i))-1;
0086         csense(end+1) = 'L';
0087     end
0088     
0089     % vartype
0090     vartype = [];
0091     for i = 1:n
0092         vartype(i) = 'C';
0093     end
0094     for i = 1:2*n
0095         vartype(end+1) = 'B';
0096     end
0097 
0098     % lb,ub
0099     lb = [model.lb; zeros(2*n,1)];
0100     ub = [model.ub; ones(2*n,1)];    
0101     % c
0102     c = [model.c; zeros(2*n,1)];
0103 
0104     
0105     % create structure
0106     MILPproblem.A = A;
0107     MILPproblem.b = b;
0108     MILPproblem.c = c;
0109     MILPproblem.csense = csense;
0110     MILPproblem.lb = lb;
0111     MILPproblem.ub = ub;
0112     MILPproblem.osense = -1;
0113     MILPproblem.vartype = vartype;
0114     MILPproblem.x0 = [];%zeros(2*n,1);
0115     %MILPproblem.intSolInd = [];
0116     %MILPproblem.contSolInd = [];
0117     
0118 %    pause;
0119     MILPsol = solveCobraMILP(MILPproblem)
0120 %    MILPsol.full
0121     NZ(:,end+1) = abs(MILPsol.full(1:n))>.000000001;
0122     PrevNZ = NZ(:,end);
0123 
0124 
0125     if (abs(MILPsol.full - maxObjective) > .001)
0126         'done';
0127         return;
0128     end
0129     solution.fluxes = [solution.fluxes,MILPsol.full(1:n)];
0130     solution.nonzero = [solution.nonzero, NZ(:,end)];
0131     solution
0132 end
0133 
0134 return;

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