%%% Tyson (1991) model of cell cycle
%%% Matlab ode solver implementation
%
%


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% Step 1:  Define constants 

global k1norm k2 k3 k4prime k4 k5 k6 k7 k8 k9 ;

k1norm = 0.015 ;       % min-1
k2 = 0 ;           %
k3 = 200 ;         % min-1
k4prime = 0.018 ;  % min-1
k4 = 180 ;         % min-1 (range:  10-1000)
% k5[~P] = 0 ;
k5 = 0 ;           % technically k5*[~P]
k6 = 1 ;           % min-1 (range:  0.1-100)
k7 = 0.6 ;         % min-1
% k8[~P] >> k9
% k9 >> k6
k9 = 100*k6 ;
k8 = 100*k9 ;      % technically k8*[~P] 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% Step 2:  Define simulation time 

tlast = 200 ; % min

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% Step 3:  Initial conditions 

Y = 0.1 ;
YP = 0.1 ;
C2 = 0.5 ;
CP = 0.5 ;
M = 0.1 ;
pM = 0.1 ;

statevar_i = [Y,YP,C2,CP,M,pM] ;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% Step 4:  Run it!

[time,statevars] = ode15s(@dydt_tyson,[0,tlast],statevar_i) ;

whos

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% % % Here you need to set Y equal to the first column of 'statevars', 
% % % YP equal to the second column, etc.
% % % The 'whos' command above should help by showing all the variables
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% Step 5:  Plot/save results

figure
hold on
plot(time,Y,'k') 
plot(time,YP,'b') 
plot(time,C2,'r')
plot(time,CP,'g')
plot(time,M,'m')
plot(time,pM,'c')
legend('cyclin','cyclin-P','cdc2','cdc2-P','MPF','preMPF')

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% % % To reproduce Figure 2 in the paper, you will need to add lines here
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%