network_integrate_kin.m

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

function [t,x,v] = network_integrate_kin(x0,p,T1,T2,N,velocity_fct,delta_par)

%[t,x,v] = network_integrate_kin(x0,p,T1,T2,N,velocity_fct,delta_par)
%
%time integration for metabolic systems with oscillatory parameter
%given kinetics function as m-file and structure array of parameters
%
%ARGUMENTS
% x0           initial concentrations (vector)
% p            parameters (structure array)
% T1, T2       integration times for tuning and simulation
% N            stoichiometric matrix
% velocity_fct function handle to kinetics function
% delta_par    structure array describing the parameter perturbation
%              (optional)
%
%FIELDS OF delta_par: 
% type (optional)  'complex', 'cos'
% omega            circular frequency
% name             parameter name
% value            amplitude of oscillation
%
%For extracting a single oscillation period, see 'choose_period'

odeoptions=[];

if exist('delta_par','var'),
  if delta_par.omega~=0,
    odeoptions1=odeset('MaxStep',2*pi/delta_par.omega/10);
    odeoptions2=odeset('MaxStep',2*pi/delta_par.omega/100);
  end

if ~isfield(delta_par,'type'), delta_par.type='complex'; end

end

if exist('delta_par','var'),
  switch delta_par.type,
    case 'cos',
      [t,x] = ode15s(@derivatives_pert_cos,[ -T1, 0],x0,odeoptions1,N,p,velocity_fct,delta_par);
      [t,x] = ode15s(@derivatives_pert_cos,[0:T2/100:T2],x(end,:).',odeoptions2,N,p,velocity_fct,delta_par);
    case 'complex',
      [t,x] = ode15s(@derivatives_pert_complex,[ -T1, 0],x0,odeoptions1,N,p,velocity_fct,delta_par);
      [t,x] = ode15s(@derivatives_pert_complex,[0:T2/100:T2],x(end,:).',odeoptions2,N,p,velocity_fct,delta_par);
  end  
else
  if T1>0,
    [t,x] = ode15s(@derivatives,[-T1, 0],x0,odeoptions,N,p,velocity_fct);
    [t,x] = ode15s(@derivatives,[ 0:T2/100:T2],x(end,:).',odeoptions,N,p,velocity_fct);
  else,
    [t,x] = ode15s(@derivatives,[ 0:T2/1000:T2],x0,odeoptions,N,p,velocity_fct);
  end
end

if exist('delta_par','var'),
  switch delta_par.type,
    case 'cos',
      for j=1:size(t),
	delta_pt = real(delta_par.value*exp(i*delta_par.omega*t(j)));
	this_p   = setfield(p,delta_par.name,getfield(p,delta_par.name) + delta_pt);
	v(j,:)   = feval(velocity_fct,x(j,:).',this_p).';
      end
    case 'complex',
      for j=1:size(t),
	delta_pt = delta_par.value*exp(i*delta_par.omega*t(j));
	this_p   = setfield(p,delta_par.name,getfield(p,delta_par.name) + delta_pt);
	v(j,:)   = feval(velocity_fct,x(j,:).',this_p).';
      end
  end
  
else
    for j=1:size(t),
      v(j,:)= feval(velocity_fct,x(j,:).',p).';
    end
end

function x_dot = derivatives(t,x,N,p,velocity_fct)

x_dot = N * feval(velocity_fct,x,p);

function x_dot = derivatives_pert_cos(t,x,N,p,velocity_fct,delta_par)

delta_pt = real(delta_par.value*exp(i*delta_par.omega*t));
this_p   = setfield(p,delta_par.name,getfield(p,delta_par.name) + delta_pt);

x_dot = N * feval(velocity_fct,x,this_p);


function x_dot = derivatives_pert_complex(t,x,N,p,velocity_fct,delta_par)

delta_pt = delta_par.value*exp(i*delta_par.omega*t);
this_p   = setfield(p,delta_par.name,getfield(p,delta_par.name) + delta_pt);

x_dot = N * feval(velocity_fct,x,this_p);