ETH 2006 Matlab Ode

From 2006.igem.org

casc3.m

% Run model for engineered 3-step cascade
%
% version 1.0 - 06/05/19
%
% requires: ode_casc3.m

%% assign model parameters

a1  =   0.02;
a2  =   0.02;
a3  =   0.02;
k1  = 100;
k2  =  10;
k3  =   1;
d1  =   0.1;
d2  =   0.1;
d3  =   0.1;
K1  = 100;
K2  = 100;
K3  =  10;
n   =   2;

I   =  100; % inducer concentration

% compile parameter vector

para = [a1;a2;a3;k1;k2;k3;d1;d2;d3;K1;K2;K3;n;I];

NX = 3;            % number of states 
NP = length(para); % number of parameters

% integration parameters

y0      = zeros(NX,1); % initial values
tfinal  = 120;         % final simulation time [min]
tstep   = .1;           % time step 
MaxAbsT = 1.0;         % maximal step size;

opt = odeset('RelTol',1e-3,'AbsTol',1e-6);

%% simulation

[t,y] = ode15s('ode_casc3',[0:tstep:tfinal],y0,opt,para);
  
% plot

figure(1);
clf;
plot(t,y(:,1),'r-',t,y(:,2),'g-',t,y(:,3),'k-');
xlabel('Time [min]');
ylabel('Concentration [-]');
legend('Repressor R_2','Repressor R_3','Output Z');

ode_casc3.m

% Function: calculation of r.h.s. for model 3-step cascade

function varargout = ode_casc3(t,y,flag,para)

switch flag
  case ''
  
    % parameter values    

    a1  = para(1);
    a2  = para(2);
    a3  = para(3);
    k1  = para(4);
    k2  = para(5);
    k3  = para(6);
    d1  = para(7);
    d2  = para(8);
    d3  = para(9);
    K1  = para(10);
    K2  = para(11);
    K3  = para(12);
    n   = para(13);
    I   = para(14);
    
    % state variables
    
    R2 = y(1);
    R3 = y(2);
    Z  = y(3);
    
    % calculate RHS states
    
    dy(1) = a1*k1 + (k1*(I/K1)^n)/(1+(I/K1)^n) - d1*R2;    
    dy(2) = a2*k2 + 1/(1+(R2/K2)^n) - d2*R3;
    dy(3) = a3*k3 + 1/(1+(R3/K3)^n) - d3*Z;
      
   varargout{1} = dy(:);
      
  case 'init'

    y0 = y;
    
  otherwise
    error(['Unknown flag ''' flag '''.']);

end

return