cba_homogeneous_cost

 CBA_HOMOGENEOUS_COST - Determine economic potentials from the principle of even cost

 [w, delta_w, y, zx] = cba_homogeneous_cost(network, v, cba_constraints, y_given, method)

 Assumptions: 
  o Balanced distribution of enzyme cost (min = sum(u^2))
  o Benefit arises only from fluxes and external production (not from concentrations)

 Input
   y_given: (optional) scalar or vector (#reactions x 1) of enzyme costs to be approximated

   method: {'cost', 'force'}  if method == 'force', the vector y_given does not represent
                                    the cost y, but the force y./v; in this case, it has to 
                                    be explicitly given as a vector

0001 function [w, delta_w, y, zx] = cba_homogeneous_cost(network,v,cba_constraints,y_given,method)
0002 
0003 % CBA_HOMOGENEOUS_COST - Determine economic potentials from the principle of even cost
0004 %
0005 % [w, delta_w, y, zx] = cba_homogeneous_cost(network, v, cba_constraints, y_given, method)
0006 %
0007 % Assumptions:
0008 %  o Balanced distribution of enzyme cost (min = sum(u^2))
0009 %  o Benefit arises only from fluxes and external production (not from concentrations)
0010 %
0011 % Input
0012 %   y_given: (optional) scalar or vector (#reactions x 1) of enzyme costs to be approximated
0013 %
0014 %   method: {'cost', 'force'}  if method == 'force', the vector y_given does not represent
0015 %                                    the cost y, but the force y./v; in this case, it has to
0016 %                                    be explicitly given as a vector
0017 
0018 % THE LATTER OPTION HAS NOT BEEN TESTED YET
0019 
0020 
0021 eval(default('y_given','[]','method','''cost'''));
0022 
0023 ind_act = find(v ~=0 );
0024 ind_int = find(network.external ==0);
0025 ind_ext = find(network.external ==1);
0026 
0027 cba_constraints = cba_update_constraints(cba_constraints,network.N(find(network.external),:));
0028 
0029 switch method, 
0030   case 'cost',
0031     %% if cost unknown -> choose according to benefit function
0032     if isempty(y_given), y_given = [cba_constraints.zv' * v ] / length(ind_act); end
0033 
0034     %% cost: scalar ->  vector
0035     if length(y_given) == 1, y_given = y_given * double(v~=0); end
0036 
0037   case 'force',
0038     y_given(ind_act) = y_given(ind_act) ./ v(ind_act);
0039 end
0040 
0041 if sum(y_given(ind_act) <= 0), error('Impossible enzyme costs'); end 
0042 
0043 % rescale benefit function to match the predefined costs
0044 
0045 cba_constraints.z_ext = cba_constraints.z_ext * [sum(y_given)] / [cba_constraints.zv'*v];
0046 cba_constraints.z_int = cba_constraints.z_int * [sum(y_given)] / [cba_constraints.zv'*v];
0047 cba_constraints.zv    = cba_constraints.zv    * [sum(y_given)] / [cba_constraints.zv'*v];
0048 
0049 v_act        = v(ind_act);
0050 y_given_act  = y_given(ind_act);
0051 N_int        = network.N(ind_int,ind_act);
0052 
0053 
0054 % quadratic function min = wc' * A * wc + a' * wc
0055 % s.t.                 B * wc <= b;
0056 %(without constraints the solution would be simply ww_int = - pinv(A) * a;)
0057 
0058 epsilon = 10^-3;
0059 
0060 A  = N_int * diag( [v_act.^2] ./ y_given_act ) * N_int';
0061 a  = N_int * diag( [v_act.^2] ./ y_given_act ) * cba_constraints.zv(ind_act);
0062 B  = - [diag(v_act) * N_int'];
0063 nr = length(v_act);
0064 b  = - [epsilon * ones(nr,1) - diag(v_act) * cba_constraints.zv(ind_act)];
0065 
0066 % regularisation
0067 
0068 if find(eig(A)==0),
0069   alpha = 10^-8; A = A + alpha * eye(size(A));
0070 end
0071 
0072   if exist('cplexqp','file'),
0073     opt = cplexoptimset('Display','off','Algorithm','interior-point-convex');
0074     w_int = cplexqp(A,a,[],[],[],[],[],[],[],opt); 
0075   else,
0076     opt = optimset('Display','off','Algorithm','interior-point-convex');
0077     w_int = quadprog(A,a,[],[],[],[],[],[],[],opt); 
0078   end
0079 
0080 if(sum(double(B * w_int > b))),
0081   %% constraints violated
0082   if exist('cplexqp','file'),
0083     opt = cplexoptimset('Algorithm', 'interior-point-convex','Display','off');
0084     w_int = cplexqp(A,a,B,b,[],[],[],[],[],opt);
0085   else
0086     opt = optimset('Algorithm', 'interior-point-convex','Display','off');
0087     w_int = quadprog(A,a,B,b,[],[],[],[],[],opt);
0088   end
0089 end
0090 
0091 
0092 % ---------------------
0093 
0094 w(ind_ext,1) = cba_constraints.z_ext;
0095 w(ind_int)   = w_int;
0096 
0097 delta_w = network.N'*w;
0098 
0099 y           = v .* [delta_w + cba_constraints.z_int];
0100 zx          = zeros(size(w)); 
0101 zx(ind_ext) = cba_constraints.z_ext;
0102 
0103 % ---------------------
0104 
0105 if find(y(v~=0) == 0), warning('Zero enzyme cost encountered'); end