set_formula_strings_kinetics.m

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

function [kinetics,parameter_id] = set_formula_strings_kinetics(network)

% kinetics = set_formula_strings_kinetics(network)
%
% construct a kinetics field (of type 'formula_strings') for a metabolic network

kinetics.type      = 'formula_strings';
kinetics.reactions = cell(length(network.actions),1);

z = 0;
p = {};
p_id = {};
parvalues = [];

 for i=1:length(network.actions)
  type = 'unknown';
  substrate_indices = find(network.N(:,i)<0);
  product_indices   = find(network.N(:,i)>0);
  n_in  = length(substrate_indices) ;
  n_out = length(product_indices) ;
  rev   = network.reversible(i);
  
  switch rev
    case 1,
%      if n_in == 1 & n_out == 1, r.type = 'rev_uni_uni'; end
%      if n_in == 1 & n_out == 2, r.type = 'rev_uni_bi'; end
%      if n_in == 2 & n_out == 1, r.type = 'rev_bi_uni'; end
%      if n_in == 2 & n_out == 2, r.type = 'rev_bi_bi'; end
      type = 'massr';
    case 0,
%      if n_in == 1 , r.type = 'irrev_uni'; end
%      if n_in == 2 , r.type = 'irrev_bi'; end
      type = 'massi';
  end
%  if n_in == 0 & n_out > 0, r.type = 'influx'; end
%  if n_in > 0  & n_out == 0, r.type = 'efflux'; end

   if isfield(network,'inhibition_matrix'), 
     if sum(network.inhibition_matrix(i,:)) + network.inhibited(i),
        type = 'ucii';
     end
   end

   clear r
   switch type,
     case {'massi','massr'},
       kplus_name  = ['kp'];
       kminus_name = ['km'];
       kplus_value = 1;
       kminus_value = 1;
       substrate_names = network.metabolites(substrate_indices);
       substrate_expon = abs(network.N(substrate_indices,i));
       a = [repmat(' * ',length(substrate_names),1)  char(substrate_names) ...
	    repmat(' ^ ',length(substrate_names),1) num2str(substrate_expon) ];
       substrate_string = reshape(a',1,prod(size(a)));
       r.string = [kplus_name substrate_string];
       r.parameters{1} = struct('name','kp','value',kplus_value);
       if strcmp(type,'massr'),
	 product_names = network.metabolites(product_indices);
	 product_expon = abs(network.N(product_indices,i));
	 a=   [repmat(' * ',length(product_names),1) char(product_names) ...
	       repmat(' ^ ',length(product_names),1) num2str(product_expon)  ]; 
	 product_string = reshape(a',1,prod(size(a)));
	 r.string = [r.string ' - ' kminus_name product_string];
	 r.parameters{2} = struct('name',kminus_name,'value',1);
       end
       r.string = strrep(r.string,' ^ 1','');      
       r.parameters = r.parameters';
      
    case 'ucii',
      inhibitor_indices = find(network.inhibition_matrix(i,:));
      if length(inhibitor_indices), 
	I_name = network.metabolites{inhibitor_indices}; 
      else,
	I_name  = ['I']; 
      end
      S_name  = network.metabolites{substrate_indices};
      V_name  = ['V'];
      Km_name = ['Km'];
      Ki_name = ['Ki'];
      r.string = sprintf('(%s * %s / %s )/(1 + %s / %s + %s /%s)',V_name,S_name,Km_name ,S_name,Km_name ,I_name,Ki_name ) ;
      r.parameters = {...
	  struct('name',V_name,'value',1 ),...
	  struct('name',Km_name,'value',1),...
	  struct('name',Ki_name,'value',1)};
      if ~length(inhibitor_indices), r.parameters=[r.parameters { struct('name',I_name,'value',1 )}];
      end
  end
  kinetics.reactions{i}=r;
 end
 
 z=1;
 for it=1:length(network.actions),
   par = kinetics.reactions{it}.parameters;
   for it2 = 1: length(par),
     kinetics.parameters{z}   = ['v' num2str(it) '_' par{it2}.name];
     kinetics.parameter_values(z)   = par{it2}.value;
     kinetics.reactions{it}.parameters{it2}.index = z;
     kinetics.reactions{it}.parameters{it2} = rmfield(kinetics.reactions{it}.parameters{it2},'value');
     if isfield(network,'action_id'),
       parameter_id{z}.object        = [network.action_id{it}.object  '_' par{it2}.name];
       parameter_id{z}.unit         = 'unknown';
       parameter_id{z}.quantity      = 'reaction parameter';
       parameter_id{z}.is_reaction  = 0;
     end
     z=z+1;
   end
 end
 
 kinetics.parameter_values = kinetics.parameter_values';
 kinetics.parameters = kinetics.parameters';
     if isfield(network,'action_id'),
       parameter_id = parameter_id';
 end