network_velocities.m

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

% v = network_velocities(s,network,kinetics,split)
%
% calculate flux velocities for a network, given a data structure
% describing the kinetics type and containing the kinetic parameters
%
% the field 'network.kinetics' is overridden by 'kinetics', if
% given as an argument
%
% There are alternative ways to specify the kinetics, indicated by kinetics.type:
%
% (i)   'numeric': the kinetics of all reactions are defined by a single m-file
%       indicated by field 'velocity_function'
%       with parameters from field 'parameters'
%
% (ii)  'mass-action': (all reactions are of mass action type): 
%       the parameters are given as vectors in additional fields
%       k_fwd (forward rate constants), k_bwd rate constants, exponents
%       (for reactions of type 2A -> B, optional)
%
% (iii) 'saturated': kinetics contains field 'reactions', a cell structure of
%       structure arrays, each describing one reaction: the kind of kinetics is given in the field
%
%       'type'                    parameters (further fields)
%        'mass-action'           k, km
%        'michaelis_menten_rev'  Vm_S, Km_S, Vm_P, Km_P
%        'multiple_michaelis_menten_rev'  Vm_S, Km_S, Vm_P, Km_P
%        'influx'                v_in         
%        'efflux'                k_out
%        'rev_uni_uni'           Pk, Pkm, k, km   (k, km: vectors of length 4 )
%        'rev_bi_uni',           Pk, Pkm, k, km   (k, km: vectors of length 4 )
%        'rev_uni_bi'            Pk, Pkm, k, km   (k, km: vectors of length 4 )
%        'rev_bi_bi'             Pk, Pkm, k, km   (k, km: vectors of length 4 )
%        'rev_multiple'          Pk, Pkm, k, km   (k, km: vectors of length 1+n_substrates+n_products )
%        'irrev_uni'             Pk, k
%        'irrev_bi'              Pk, k
%
% (iv)   'formula_strings': kinetics contains field 'reactions', a cell structure of
%        structure arrays, each describing one reaction. parameter
%        names and values are given in a column cell array of strings
%        and a column vector, respectively.
%
% s: metabolite concentrations
%
% Note: this function is slow. For integration steps, use 'integrate_network_der' instead.

function v = network_velocities(s,network,kinetics,split)

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

if size(s,1)==1, s=s'; end

KINETICS   = kinetics;
N   = network.N;

n_A = size(N,2);

switch KINETICS.type
      
  case 'formula_strings',
    for it = 1:length(network.metabolites),
%      fprintf([ network.metabolites{it} ' = ' num2str(s(it)) ';\n']);
      eval([ network.metabolites{it} ' = ' num2str(s(it)) ';']);
    end
    for i=1:length(KINETICS.reactions),
      r = KINETICS.reactions{i};
      for itt = 1:length(r.parameters),
	%	fprintf([r.parameters{itt}.name ' = '
              %	num2str(KINETICS.parameter_values(r.parameters{itt}.index)) ';']);
	eval([r.parameters{itt}.name ' = ' num2str(KINETICS.parameter_values(r.parameters{itt}.index)) ';']);
      end
      %fprintf(['v(i) = ' r.string ';\n']);
      eval(['v(i) = ' r.string ';']);
    end
    v=v';
%    pause;
    
  case 'numeric',
    if isfield(KINETICS,'use_only_met'),
      dummy = ones(KINETICS.n_met_tot,1);
      dummy(KINETICS.use_only_met)=s;
      s=dummy;
    end

    v = feval(KINETICS.velocity_function,s,KINETICS.parameters);
    if isfield(KINETICS,'use_only_act'),
       v=v(KINETICS.use_only_act);
    end

  case 'mass-action',
    if isfield(KINETICS,'exponents'),
      Nf = KINETICS.exponents.*(N<0);
      Nb = KINETICS.exponents.*(N>0);
    else,
      Nf  = abs(N.*(N<0));
      Nb  =     N.*(N>0);      
    end
    v  = (  (                 KINETICS.k_fwd'  .* prod( repmat(s,1,n_A) .^ Nf) ) ...
	      - ( ( network.reversible' .* KINETICS.k_bwd') .* prod( repmat(s,1,n_A) .^ Nb) ) )';
  
  case 'saturated',
    
    if ~split,
      v = nan*ones(n_A,1);
    
      for i=1:length(network.actions),
	
	R   = KINETICS.reactions{i};
	S   = s(find(network.N(:,i)<0));
	P   = s(find(network.N(:,i)>0));
	
	switch R.type,
	  
	  case 'mass-action'	  
	    v(i) = prod(s.^Nf)* R.k - prod(s.^Nb) * R.km * network.reversible(i);	    
	  case 'michaelis_menten_rev',
	    v(i)= ( R.Vm_S/R.Km_S * S -R.Vm_P/R.Km_P * P)/(1+S/R.Km_S+P/R.Km_P);	    
	  case 'multiple_michaelis_menten_rev',
	    v(i)= ( R.Vm_S/R.Km_S * prod(S) -R.Vm_P/R.Km_P * prod(P))/(1+S/R.Km_S+P/R.Km_P);
	  case 'multiple_michaelis_menten_irrev',
	    v(i)= ( R.Vm_S/R.Km_S * prod(S) )/(1+S/R.Km_S);	    
	  case 'influx',
	    v(i) = R.v_in;	    
	  case 'efflux',
	    v(i)= prod(S) * R.k_out;	    
	  case 'rev_uni_uni',	    
	    v(i) = ( R.Pk * S - R.Pkm * P ) / ( S * R.k(1) + R.km(1) +  R.k(2) + P * R.km(2));	    
	  case 'rev_bi_uni',
	    v(i) = ( R.Pk * S(1) * S(2) - R.Pkm * P ) /...
		 ( ( R.km(2) + S(2)*R.k(2))*(S(1)*R.k(1)+ P*R.km(3)) + R.k(3)*(S(1)*R.k(1) + S(2)*R.k(2)) ...
		   + R.km(1)*(R.km(2) + R.k(3) + P*R.km(3)) );	    
	  case 'rev_uni_bi',
	    v(i) = ( R.Pk * S - R.Pkm * P(1) * P(2) )/...
		   ( ( R.k(2) + P(1)*R.km(2))*(P(2)*R.km(3) + S*R.k(1)) + R.km(1)*(P(2)*R.km(3) + P(1)*R.km(2)) ...
		     + R.k(3)*(R.k(2) + R.km(1) + S*R.k(1)) );	    
	  case 'rev_bi_bi',
	    v(i) = ( R.Pk * S(1) * S(2)- R.Pkm * P(1) * P(2))/...
		   ( (R.km(2)+R.k(3))*(R.k(4)*R.km(1) + R.k(1)*R.k(4)*S(1) + R.km(4)*R.km(1)*P(2)) ...
		     + S(2)*R.k(2) *(R.k(3) *(R.k(4)+R.km(4)*P(2)) + S(1)*R.k(1)*(R.k(3)+R.k(4)))...
		     + P(1)*R.km(3)*(R.km(2)*(R.km(1)+R.k(1)*S(1)) + P(2)*R.km(4)*(R.km(1)+R.km(2)))...
		     + (R.k(1)*S(1)+R.km(4)*P(2))*R.k(2)*R.km(3)*S(2)*P(1) );	    
	  case 'rev_multiple',
	    v(i) = ( R.Pk * prod(S) - R.Pkm * prod(P) ) / ( prod(S) * R.k(1) + R.km(1) +  R.k(2) + prod(P) * R.km(2));	    
	  case 'irrev_uni',
	    v(i) =  R.Pk * S / (S * R.k(1) + R.k(2) );	    
	  case 'irrev_bi',
	    v(i) = ( R.Pk * S(1) * S(2) ) /...
		   (  S(2)*R.k(2) * S(1)*R.k(1)  + R.k(3)*(S(1)*R.k(1) + S(2)*R.k(2)) );
	end
      end	
    else,
      v = nan*ones(2*n_A,1);
      
      for i=1:length(network.actions),
	
	R   = KINETICS.reactions{i};
	S   = s(find(network.N(:,i)<0));
	P   = s(find(network.N(:,i)>0));
	
	switch R.type,
	  
	  case 'mass-action'	  
	    v(2*i-1:2*i) = [prod(s.^Nf)* R.k;  prod(s.^Nb) * R.km * network.reversible(i)];	    
	  case 'michaelis_menten_rev',
	    v(2*i-1:2*i)= [ R.Vm_S/R.Km_S * S ; R.Vm_P/R.Km_P * P]/(1+S/R.Km_S+P/R.Km_P);	    
	  case 'multiple_michaelis_menten_rev',
	    v(2*i-1:2*i)= [ R.Vm_S/R.Km_S * prod(S);  R.Vm_P/R.Km_P * prod(P)]/(1+S/R.Km_S+P/R.Km_P);
	  case 'multiple_michaelis_menten_irrev',
	    v(2*i-1:2*i)= [ R.Vm_S/R.Km_S * prod(S) /(1+S/R.Km_S); 0];	    
	  case 'influx',
	    v(2*i-1:2*i) = [R.v_in;0];	    
	  case 'efflux',
	    v(2*i-1:2*i)= [prod(S) * R.k_out; 0];
	  case 'rev_uni_uni',	    
	    v(2*i-1:2*i) = [ R.Pk * S; R.Pkm * P ] / ( S * R.k(1) + R.km(1) +  R.k(2) + P * R.km(2));	    
	  case 'rev_bi_uni',
	    v(2*i-1:2*i) = [ R.Pk * S(1) * S(2) ;  R.Pkm * P ] /...
		 ( ( R.km(2) + S(2)*R.k(2))*(S(1)*R.k(1)+ P*R.km(3)) + R.k(3)*(S(1)*R.k(1) + S(2)*R.k(2)) ...
		   + R.km(1)*(R.km(2) + R.k(3) + P*R.km(3)) );	    
	  case 'rev_uni_bi',
	    v(2*i-1:2*i) = [ R.Pk * S ; R.Pkm * P(1) * P(2) ]/...
		   ( ( R.k(2) + P(1)*R.km(2))*(P(2)*R.km(3) + S*R.k(1)) + R.km(1)*(P(2)*R.km(3) + P(1)*R.km(2)) ...
		     + R.k(3)*(R.k(2) + R.km(1) + S*R.k(1)) );	    
	  case 'rev_bi_bi',
	    v(2*i-1:2*i) = [ R.Pk * S(1) * S(2);  R.Pkm * P(1) * P(2)]/...
		   ( (R.km(2)+R.k(3))*(R.k(4)*R.km(1) + R.k(1)*R.k(4)*S(1) + R.km(4)*R.km(1)*P(2)) ...
		     + S(2)*R.k(2) *(R.k(3) *(R.k(4)+R.km(4)*P(2)) + S(1)*R.k(1)*(R.k(3)+R.k(4)))...
		     + P(1)*R.km(3)*(R.km(2)*(R.km(1)+R.k(1)*S(1)) + P(2)*R.km(4)*(R.km(1)+R.km(2)))...
		     + (R.k(1)*S(1)+R.km(4)*P(2))*R.k(2)*R.km(3)*S(2)*P(1) );	    
	  case 'rev_multiple',
	    v(2*i-1:2*i) = [ R.Pk * prod(S); R.Pkm * prod(P)] / ( prod(S) * R.k(1) + R.km(1) +  R.k(2) + prod(P) * R.km(2));	    
	  case 'irrev_uni',
	    v(2*i-1:2*i) =  [R.Pk * S / (S * R.k(1) + R.k(2) );0];	    
	  case 'irrev_bi',
	    v(2*i-1:2*i) = [ ( R.Pk * S(1) * S(2) ) /...
		   (  S(2)*R.k(2) * S(1)*R.k(1)  + R.k(3)*(S(1)*R.k(1) + S(2)*R.k(2)) ); 0];
	end
	
      end
    end
end