% exercise 2
% First In - first out temporal order
% x=[X Z1 Z2 Z3]';

function [T,y]=exercise2_KattrinArning(tstart,tend)

tspan=[tstart,tend];% time interval for solution 
ICList=[0;0;0;0];% initial conditions

[T,y]=ode23(@ODEsys,tspan,ICList);%solve ODE system
figure;subplot(2,1,1);plot(T,y(:,1));title('X(t)');xlabel('t')
subplot(2,1,2);plot(T,y(:,2),'b',T,y(:,3),'r',T,y(:,4),'g');title('Z1(t), Z2(t), Z3(t)');xlabel('t')


function dxdt=ODEsys(t,x)

%constants
beta_X=1;
beta_Z1=1;
beta_Z2=1;
beta_Z3=1;
alpha_X=1;
alpha_Z1=1;
alpha_Z2=1;
alpha_Z3=1;
K_XZ1=0.5;
K_Z1Z2=0.5;
K_Z2Z3=0.5;
%K_XZ2=0.4;
%K_XZ3=0.7;
%K_Z1Z3=0.5;

% ODE system
% if regulator X reaches a certain threshold (K_XZ1), gene Z1 is turned on,
% if X falls below K_XZ, gene Z1 is turned off again
% if Z1 reaches a certain threshold (K_Z1Z2), gene Z2 is turned on, if Z1
% falls below K_Z1Z2, gene Z2 is turned off again
% if Z2 reaches a certain threshold (K_Z2Z3), gene Z3 is turned on, if Z2
% falls below K_Z2Z3, gene Z3 is turned off again

dxdt=[beta_X*(1-stepfunc(t-5))*(stepfunc(t-1))-alpha_X*x(1);...
      beta_Z1*stepfunc(x(1)-K_XZ1)-alpha_Z1*x(2);...
      beta_Z2*stepfunc(x(2)-K_Z1Z2)-alpha_Z2*x(3);...
      beta_Z3*stepfunc(x(3)-K_Z2Z3)-alpha_Z3*x(4)];


function teta=stepfunc(t)
teta=(t >= 0);