cba_flux2model_mrl_constraint_xi

 function [g, geq, Ec_un, beta_M, beta_A, beta_I, c, mu, CJ] = cba_flux2model_mrl_constraint_xi(xi,Mplus,Mminus,Wplus,Wminus,v,N,ind_ext,cba_constraints,kinetic_law)

 choose state parameters (mu0, c, alphas) such that
  o all flux control coefficient on benefit are positive
  o all reaction affinities and fluxes have the same signs

0001  function [g, geq, Ec_un, beta_M, beta_A, beta_I, c, mu, CJ] = cba_flux2model_mrl_constraint_xi(xi,Mplus,Mminus,Wplus,Wminus,v,N,ind_ext,cba_constraints,kinetic_law,h)
0002 
0003 % function [g, geq, Ec_un, beta_M, beta_A, beta_I, c, mu, CJ] = cba_flux2model_mrl_constraint_xi(xi,Mplus,Mminus,Wplus,Wminus,v,N,ind_ext,cba_constraints,kinetic_law)
0004 %
0005 % choose state parameters (mu0, c, alphas) such that
0006 %  o all flux control coefficient on benefit are positive
0007 %  o all reaction affinities and fluxes have the same signs
0008 
0009 eval(default('h','ones(size(v))'));
0010 
0011 [nr, nm] = size(Mplus);
0012 
0013 M       = Mplus + Mminus;
0014 alpha_M = zeros(size(M)); 
0015 alpha_A = zeros(size(M)); 
0016 alpha_I = zeros(size(M)); 
0017 beta_M  = zeros(size(M)); 
0018 beta_A  = zeros(size(M)); 
0019 beta_I  = zeros(size(M)); 
0020 
0021 ind_KM = find(M);
0022 ind_KA = find(Wplus);
0023 ind_KI = find(Wminus);
0024 n_KM   = length(ind_KM);
0025 n_KA   = length(ind_KA);
0026 n_KI   = length(ind_KI);
0027 
0028 mu0  = xi(1:nm);
0029 ln_c = xi(nm+[1:nm]);
0030 c    = exp(ln_c);
0031 
0032 mu   = mu0 + RT*ln_c;
0033 A    = -N'*mu;
0034 A(A==0) = 10^-10; % simple fix to avoid divergence
0035 zeta = exp(h.*A/RT);
0036 
0037 E_sc_T = diag(zeta./[zeta-1]) * Mplus - diag(1./[zeta-1]) * Mminus;
0038 
0039 alpha_M(ind_KM) = 1-xi(2*nm+[1:n_KM]);
0040 beta_M( ind_KM) =   xi(2*nm+[1:n_KM]);
0041 
0042 alpha_A(ind_KA) = 1-xi(2*nm+n_KM + [1:n_KA]);
0043 beta_A( ind_KA) =   xi(2*nm+n_KM + [1:n_KA]);
0044 
0045 alpha_I(ind_KI) = 1-xi(2*nm+n_KM + n_KA + [1:n_KI]);
0046 beta_I( ind_KI) =   xi(2*nm+n_KM + n_KA + [1:n_KI]);
0047 
0048 % -----------------------------------------------------------
0049 % copy from compute_modular_elasticities
0050 
0051 switch kinetic_law,  
0052 
0053   case {'cs','ms'}
0054     alpha_M( find((Mplus==0).*(Mminus==0)) ) = 1;  % irrelevant elements
0055     psi_plus    = prod((1./alpha_M) .^ Mplus, 2);
0056     psi_minus   = prod((1./alpha_M) .^ Mminus,2);
0057 
0058   case {'ds','fd'},
0059     alpha_M( find((Mplus==0).*(Mminus==0)) ) = 1/2; % irrelevant elements
0060     xi_plus  = prod((1./alpha_M-1) .^ Mplus, 2);
0061     xi_minus = prod((1./alpha_M-1) .^ Mminus,2);
0062 
0063   otherwise error('Unknown kinetic law');
0064 
0065 end
0066 
0067 switch kinetic_law,  
0068   case 'cs',
0069     D               = psi_plus + psi_minus - 1;
0070     E_sc_Dc_kinetic = diag(1./D) * beta_M .* [diag(psi_plus) * Mplus + diag(psi_minus) * Mminus];
0071   case 'ms',
0072     D               = psi_plus .* psi_minus;
0073     E_sc_Dc_kinetic = beta_M   .* [ Mplus + Mminus ];
0074   case 'ds',
0075     D               = xi_plus + xi_minus + 1;
0076     E_sc_Dc_kinetic = diag(1./D) * [ diag(xi_plus) * Mplus + diag(xi_minus) * Mminus ];
0077   case 'rp',
0078     D               = ones(nm,1);
0079     E_sc_Dc_kinetic = zeros(nr,nm);
0080   case 'fd',
0081     D               = sqrt(xi_plus .* xi_minus);
0082     E_sc_Dc_kinetic = 1/2 * [ Mplus + Mminus ];
0083 end
0084 
0085 E_sc_c_regulation = alpha_A .* Wplus - beta_I .* Wminus;
0086 
0087 % -----------------------------------------------------------
0088 
0089 Ec_un = diag(v) * [E_sc_T - E_sc_Dc_kinetic + E_sc_c_regulation] * diag(1./c);
0090 
0091 ind_int = setdiff(1:nr,ind_ext);
0092 N_int   = N(ind_int,:);
0093 E_int   = Ec_un(:,ind_int);
0094 
0095 CJ      = eye(nr) - E_int * pinv(N_int*E_int) * N_int;
0096 
0097 g       = [- diag(sign(v)) * CJ' * cba_constraints.zv + cba_constraints.zx_scaled_min; ...
0098            ];%-sign(v .* A) .* double(v~=0) - 0.1 ];
0099 
0100 geq     = [];