elasticities.m

％The Metabolic Networks Toolbox contains functions to create, ％modify, display, and simulate bioche

%[epsilon_1,pi_1,parameter_names,epsilon_2,rho_2,pi_2] = elasticities(network,s,split);
%
% reaction elasticities for mass action kinetics and metabolite concentrations s
%
% The parameters comprise all kinetic parameters, followed by all external concentrations
%
% For mass--action kinetics:
% analytically determined. the parameter vector is supposed to
% contain all forward rate constants, followed by all backward rate
% constants (including dummy values for the irreversible reactions)
%
% For other kinetics:
%
% see also 'numerical_elasticities'

function [epsilon_1,pi_1,parameter_names,epsilon_2,rho_2,pi_2] = elasticities(network,s,split);

if ~exist('split','var'), split = 0; end

if isfield(network,'used'),   used=network.used;
else,                         used=ones(length(network.actions),1); 
end

p  = network.kinetics;
N  = network.N;
Nf = abs(N.*(N<0));
Nb =     N.*(N>0);
[n_S,n_A]     = size(N);
if split, n_A = 2*n_A; end

switch p.type
  
  case 'mass-action',

    p.k_bwd = p.k_bwd .* network.reversible;
    
    s = s + (10^-10).*(abs(s)<10.^-10);
    epsilon_1 =  diag( p.k_fwd)*(diag(1./s)*Nf*diag(prod( repmat(s,1,n_A) .^ Nf)))' ...
	-diag( p.k_bwd)*(diag(1./s)*Nb*diag(prod( repmat(s,1,n_A) .^ Nb)))' ;
    
    pi_1 = [   diag(prod( repmat(s,1,n_A) .^ Nf)), ...
	     - diag(prod( repmat(s,1,n_A) .^ Nb))      ] ;
    
    pi_1 = [pi_1, epsilon_1(:,find(network.external))];
    
    parameter_names = cellstr([repmat('k',2*n_A,1), ...
		  repmat(num2str((1:n_A)'),2,1),...
		  [ repmat('+',n_A,1);  repmat('-', n_A,1)] ]);
    parameter_names = [parameter_names; network.metabolites(find(network.external))];
  
    % second order elasticities  
    if nargout>3,
     epsilon_2 = zeros(n_A,n_S,n_S);

     for i=1:n_A, 
      epsilon_2(i,:,:) = p.k_fwd(i) * prod(s.^Nf(:,i)) * ( diag(1./s) * Nf(:,i)...
                          *  Nf(:,i)' * diag(1./s) - diag(1./s.^2.*Nf(:,i)) )  ...
	               - p.k_bwd(i) * prod(s.^Nb(:,i)) * ( diag(1./s) * Nb(:,i) *...
                           Nb(:,i)' * diag(1./s) - diag(1./s.^2.*Nb(:,i)) );
    end
 
    n_parameters = length(parameter_names);   
    rho_2    = zeros(n_A,n_S,n_parameters);
    dum = [    (    prod( repmat(s,1,n_A) .^ Nf)' * (1./s'))'.*Nf , ...
           - (    prod( repmat(s,1,n_A) .^ Nb)' * (1./s'))'.*Nb ];
    for i = 1:n_A, 
      rho_2(i,:,i)     = dum(:,i);
      rho_2(i,:,i+n_A) = dum(:,i+n_A);
    end
    
    pi_2     = zeros(n_A,n_parameters,n_parameters);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%   second elasticities w.r.t. external metabolites must be fixed !!!!
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

    pi_1(:,n_A+find(network.reversible==0))    = 0;
    rho_2(:,:,n_A+find(network.reversible==0)) = 0;
    end
    
  otherwise,
    
    % compute elasticities numerically
    
    delta = 0.0001;

    v   = network_velocities(s,network,network.kinetics,split);

    S_ext       = s(find(network.external));
    S_ext_names = network.metabolites(find(network.external));
    
    
    epsilon_1 = nan*zeros(n_A,n_S);
    for i=1:n_S,
      s_pert = s; s_pert(i) = s(i)*(1+delta);
      if s_pert(i)==0, s_pert(i)=10^-6; end
      v_pert = network_velocities(s_pert,network,network.kinetics,split);
      epsilon_1(:,i) = (v_pert-v)/(s_pert(i)-s(i));
    end

   if nargout>1,    
 
   [par,parameter_names] = parameters2vector(network.kinetics,S_ext,S_ext_names);
    n_par = length(par);
 
   pi_1 = nan*zeros(n_A,n_par);
    for i=1:n_par,
      par_pert  = par; par_pert(i) = par(i)*(1+delta);
      if par_pert(i)==0, par_pert(i)=10^-6; end
      [nkinetics,Sext_pert] = vector2parameters(network.kinetics,par_pert,sum(network.external));
      s_pert    = s; s_pert(find(network.external)) = Sext_pert;
      v_pert    = network_velocities(s_pert,network,nkinetics,split);
      pi_1(:,i) = (v_pert-v)/(par_pert(i)-par(i));
    end
    end

    if nargout>3,
      
      epsilon_2 = nan*zeros(n_A,n_S,n_S);
    for i=1:n_S,
      for k=1:i,
	s_pert = s; s_pert(i) = s(i)*(1+delta); s_pert(k) = s(k)*(1+delta); 
        if s_pert(i)==0, s_pert(i)=10^-6; end 
	if s_pert(k)==0, s_pert(k)=10^-6; end
        v_pert = network_velocities(s_pert,network,network.kinetics,split);
	epsilon_2(:,i,k) = 0.5 * ( v_pert - v - epsilon_1*(s_pert-s)) / ((s_pert(i)-s(i))* (s_pert(k)-s(k)) );
	epsilon_2(:,k,i) = epsilon_2(:,i,k);
      end
    end

    rho_2     = nan*zeros(n_A,n_S,n_par);
    for i=1:n_S,
      for k=1:n_par,
	s_pert   = s;       s_pert(i) =   s(i)*(1+delta);  
	par_pert = par;   par_pert(k) = par(k)*(1+delta);  if par_pert(k)==0, par_pert(k)=10^-6; end
	[nkinetics,Sext_pert]    = vector2parameters(network.kinetics,par_pert,sum(network.external));
	s_pert(find(network.external)) = Sext_pert;
		if   s_pert(i)==s(i),   s_pert(i)=s(i)+10^-6; end 
	v_pert       = network_velocities(s_pert,network,nkinetics,split);
        rho_2(:,i,k) = 0.5 * ( v_pert - v - epsilon_1*(s_pert-s)-pi_1*(par_pert-par)) ...
	                   / ( (s_pert(i)-s(i))* (par_pert(k)-par(k)) );
      end
    end
    
    pi_2      = nan*zeros(n_A,n_par,n_par);
    for i=1:n_par,
      for k=1:i,
	par_pert    = par; par_pert(i) = par(i)*(1+delta); par_pert(k) = par(k)*(1+delta);
        if par_pert(i)==0, par_pert(i)=10^-6; end 
	if par_pert(k)==0, par_pert(k)=10^-6; end
	[nkinetics,Sext_pert]   = vector2parameters(network.kinetics,par_pert,sum(network.external));
	s_pert    = s; s_pert(find(network.external)) = Sext_pert;
	v_pert      = network_velocities(s_pert,network,nkinetics,split);
	pi_2(:,i,k) = 0.5 * ( v_pert - v - pi_1*(par_pert-par)) / ( (par_pert(i)-par(i))* (par_pert(k)-par(k)) );
  	pi_2(:,k,i) =   	pi_2(:,i,k);
      end
    end
    end
end

epsilon_1(find(1-used),:)   = 0;
pi_1(find(1-used),:)   = 0;

epsilon_1 = sparse(epsilon_1);
pi_1 = sparse(pi_1);

if nargout>3,
  epsilon_2(find(1-used),:,:) = 0;
  pi_2(find(1-used),:,:) = 0;
  rho_2(find(1-used),:,:) = 0;
end