geometricFBA

PURPOSE ^

geometricFBA finds a unique optimal FBA solution that is (in some sense)

SYNOPSIS ^

function flux = geometricFBA(model,varargin)

DESCRIPTION ^

geometricFBA finds a unique optimal FBA solution that is (in some sense)
central to the range of possible fluxes; as described in
   K Smallbone, E Simeonidis (2009). Flux balance analysis: 
   A geometric perspective. J Theor Biol 258: 311-315
   http://dx.doi.org/10.1016/j.jtbi.2009.01.027

 flux = geometricFBA(model)

INPUT
 model         COBRA model structure

OPTIONAL INPUTS
 Optional parameters can be entered as parameter name followed by
 parameter value: i.e. ...,'epsilon',1e-9)
 printLevel    [default: 1]  printing level
               = 0     silent
               = 1     show algorithm progress and times
 epsilon       [default: 1e-6]    convergence tolerance of algorithm, 
               defined in more detail in paper above
 flexRel       [default: 0] flexibility to flux bounds
               try e.g. 1e-3 if the algorithm has convergence problems

OUTPUT
 flux          unique centered flux

kieran smallbone, 5 May 2010

 This script is made available under the Creative Commons
 Attribution-Share Alike 3.0 Unported Licence (see
 www.creativecommons.org).

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SUBFUNCTIONS ^

SOURCE CODE ^

0001 function flux = geometricFBA(model,varargin)
0002 %geometricFBA finds a unique optimal FBA solution that is (in some sense)
0003 %central to the range of possible fluxes; as described in
0004 %   K Smallbone, E Simeonidis (2009). Flux balance analysis:
0005 %   A geometric perspective. J Theor Biol 258: 311-315
0006 %   http://dx.doi.org/10.1016/j.jtbi.2009.01.027
0007 %
0008 % flux = geometricFBA(model)
0009 %
0010 %INPUT
0011 % model         COBRA model structure
0012 %
0013 %OPTIONAL INPUTS
0014 % Optional parameters can be entered as parameter name followed by
0015 % parameter value: i.e. ...,'epsilon',1e-9)
0016 % printLevel    [default: 1]  printing level
0017 %               = 0     silent
0018 %               = 1     show algorithm progress and times
0019 % epsilon       [default: 1e-6]    convergence tolerance of algorithm,
0020 %               defined in more detail in paper above
0021 % flexRel       [default: 0] flexibility to flux bounds
0022 %               try e.g. 1e-3 if the algorithm has convergence problems
0023 %
0024 %OUTPUT
0025 % flux          unique centered flux
0026 %
0027 %kieran smallbone, 5 May 2010
0028 %
0029 % This script is made available under the Creative Commons
0030 % Attribution-Share Alike 3.0 Unported Licence (see
0031 % www.creativecommons.org).
0032 
0033 param = struct('epsilon',1e-6,'flexRel',0,'printLevel',1);
0034 field = fieldnames(param);
0035 if mod(nargin,2) ~= 1 % require odd number of inputs
0036     error('incorrect number of input parameters')
0037 else
0038     for k = 1:2:(nargin-1)
0039         param.(field{strcmp(varargin{k},field)}) = varargin{k+1};
0040     end
0041 end
0042 param.flexTol = param.flexRel * param.epsilon; % absolute flexibility
0043 
0044 % determine optimum
0045 FBAsolution = optimizeCbModel(model);
0046 ind = find(model.c);
0047 if length(ind) == 1
0048     model.lb(ind) = FBAsolution.f;
0049 end
0050 
0051 A = model.S;
0052 b = model.b;
0053 L = model.lb;
0054 U = model.ub;
0055 
0056 % ensure column vectors
0057 b = b(:); L = L(:); U = U(:);
0058 
0059 % Remove negligible elements
0060 J = any(A,2); A = A(J,:); b = b(J);
0061 
0062 % presolve
0063 v = nan(size(L));
0064 J = (U-L < param.epsilon);
0065 v(J) = (L(J)+U(J))/2;
0066 J = find(isnan(v));
0067 
0068 if param.printLevel
0069     fprintf('%s\t%g\n\n%s\t@%s\n','# reactions:',length(v),'iteration #0',datestr(now,16)); 
0070 end
0071 
0072 L0 = L; U0 = U;
0073 for k = J(:)'
0074     f = zeros(length(v),1); f(k) = -1;
0075     [dummy,opt,conv] = easyLP(f,A,b,L0,U0);
0076     if conv
0077         vL = max(-opt,L(k));
0078     else
0079         vL = L(k);
0080     end
0081     [dummy,opt,conv] = easyLP(-f,A,b,L0,U0);
0082     if conv
0083         vU = min(opt,U(k));
0084     else vU = U(k);
0085     end
0086     if abs(vL) < param.epsilon
0087         vL = 0;
0088     end
0089     if abs(vU) < param.epsilon
0090         vU = 0;
0091     end
0092     vM = (vL + vU)/2;
0093     if abs(vM) < param.epsilon
0094         vM = 0;
0095     end
0096     if abs(vU - vL) < param.epsilon
0097         vL = (1-sign(vM)* param.flexTol)*vM; 
0098         vU = (1+sign(vM)* param.flexTol)*vM; 
0099     end
0100     L(k) = vL; U(k) = vU;
0101 end
0102 
0103 v = nan(size(L)); 
0104 J = (U-L < param.epsilon);
0105 v(J) = (L(J)+U(J))/2; v = v.*(abs(v) > param.epsilon);
0106 
0107 if param.printLevel
0108     fprintf('%s\t\t%g\n%s\t\t%g\n\n','fixed:',sum(J),'@ zero:',sum(v==0));
0109 end
0110 
0111 % iterate
0112 J  = find(U-L >= param.epsilon); 
0113 n   = 1;
0114 mu  = [];
0115 Z   = [];
0116 
0117 while ~isempty(J)
0118     
0119     if param.printLevel
0120         fprintf('%s #%g\t@%s\n','iteration',n,datestr(now,16)); 
0121     end
0122     
0123     if n == 1
0124         M = zeros(size(L));
0125     else
0126         M = (L+U)/2;
0127     end
0128     
0129     mu(:,n) = M;                                                %#ok<AGROW>
0130     allL = L; allU = U; allA = A; allB = b;    
0131     [a1,a2] = size(A);
0132     
0133     % build new matrices
0134     for k = 1:(n-1)
0135         [b1,b2] = size(allA);
0136         f = sparse(b2+2*a2,1); f((b2+1):end) = -1;
0137         opt = -Z(k);
0138         allA = [allA,sparse(b1,2*a2);
0139             speye(a2,a2),sparse(a2,b2-a2),-speye(a2),speye(a2);
0140             f(:)'];                                             %#ok<AGROW>
0141         allB = [allB;mu(:,k);opt];                              %#ok<AGROW>
0142         allL = [allL;zeros(2*a2,1)];                            %#ok<AGROW>
0143         allU = [allU;inf*ones(2*a2,1)];                         %#ok<AGROW>
0144     end
0145     
0146     [b1,b2] = size(allA);
0147     f = zeros(b2+2*a2,1); f((b2+1):end) = -1;
0148     allA = [allA,sparse(b1,2*a2);
0149         speye(a2,a2),sparse(a2,b2-a2),-speye(a2),speye(a2)];    %#ok<AGROW>
0150     allB = [allB;M];                                            %#ok<AGROW>
0151     allL = [allL;zeros(2*a2,1)];                                %#ok<AGROW>
0152     allU = [allU;inf*ones(2*a2,1)];                             %#ok<AGROW>
0153     
0154     [v,opt,conv] = easyLP(f,allA,allB,allL,allU);
0155     if ~conv, disp('error: no convergence'); flux = (L+U)/2; return; end    
0156     
0157     opt = ceil(-opt/eps)*eps;
0158     Z(n) = opt;                                                 %#ok<AGROW>
0159     allA = [allA; sparse(f(:)')];                               %#ok<AGROW>
0160     allB = [allB; -opt];                                        %#ok<AGROW>
0161     
0162     for k = J(:)'        
0163         f = zeros(length(allL),1); f(k) = -1;
0164         [dummy,opt,conv] = easyLP(f,allA,allB,allL,allU);
0165         if conv
0166             vL = max(-opt,L(k));
0167         else
0168             vL = L(k); 
0169         end
0170         [dummy,opt,conv] = easyLP(-f,allA,allB,allL,allU);
0171         if conv
0172             vU = min(opt,U(k));
0173         else
0174             vU = U(k);
0175         end
0176         if abs(vL) < param.epsilon
0177             vL = 0;
0178         end
0179         if abs(vU) < param.epsilon
0180             vU = 0;
0181         end
0182         vM = (vL + vU)/2;
0183         if abs(vM) < param.epsilon
0184             vM = 0;
0185         end
0186         if abs(vU - vL) < param.epsilon
0187             vL = (1-sign(vM)* param.flexTol)*vM; 
0188             vU = (1+sign(vM)* param.flexTol)*vM; 
0189         end
0190         L(k) = vL; 
0191         U(k) = vU;
0192     end
0193     
0194     v = nan(size(L)); 
0195     J = (U-L < param.epsilon);
0196     v(J) = (L(J)+U(J))/2; v = v.*(abs(v) > param.epsilon);
0197 
0198     if param.printLevel
0199         fprintf('%s\t\t%g\n%s\t\t%g\n\n','fixed:',sum(J),'@ zero:',sum(v==0));
0200     end
0201     
0202     n = n+1;
0203     J = find(U-L >= param.epsilon);   
0204     
0205     flux = v;
0206 end
0207 
0208 function [v,fOpt,conv] = easyLP(c,A,b,lb,ub)
0209 %easyLP
0210 %
0211 % solves the linear programming problem:
0212 %   max c'x subject to
0213 %   A x = b
0214 %   lb <= x <= ub.
0215 %
0216 % Usage: [v,fOpt,conv] = easyLP(c,A,b,lb,ub)
0217 %
0218 %   c           objective coefficient vector
0219 %   A           LHS matrix
0220 %   b           RHS vector
0221 %   lb         lower bound
0222 %   ub         upper bound
0223 %
0224 %   v           solution vector
0225 %   fOpt        objective value
0226 %   conv        convergence of algorithm [0/1]
0227 %
0228 % the function is a wrapper for the "solveCobraLP" script.
0229 %
0230 %kieran smallbone, 5 may 2010
0231 
0232 csense(1:length(b)) = 'E';
0233 model = struct('A',A,'b',b,'c',full(c),'lb',lb,'ub',ub,'osense',-1,'csense',csense);
0234 solution = solveCobraLP(model);
0235 v = solution.full;
0236 fOpt = solution.obj;
0237 conv = solution.stat == 1;

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