%% %%%%%% MATLAB ODE Simulation of the Budding Yeast Cell Cycle %%%%%% %%
% The ODE model is taken from Chen(2004): 
% Chen, K.C. and Calzone, L. and Csikasz-Nagy, A. and Cross, F.R. and
% Novak, B. and Tyson, J.J. 2004. Integrative analysis of cell cycle 
% control in budding yeast, Molecular Biology of the Cell 15:3841-3862
% The equations and parameters were downloaded from
% http://mpf.biol.vt.edu/research/budding_yeast_model/model_download/bychen04.ode
% and converted into MATLAB code by Kartik Subramanian (VirginiaTech) and
% Kamaludin Dingle (Oxford Univ.). September 2010.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%% main_yeast.m %%
% This is the main file of the ODE simulation. It runs the MATLAB ODE 
% solver and plots the output, which reproduces figure 2 in the paper
% Chen(2004) (NB: only the plot of variables mass,BUD,ORI and SPN).

clc;
clear all;

% This calls the function (variable name) so that the variables are
% recognized by their name rather than y(1), y(2)...- see comments 
% in function for details
variable_name; 

% y0 is a vector containing the initial conditions for the system of ODEs,
% which are defined in the function "initial_conds"
y0 = initial_conds;

% Integration parameters - t0 is the time the simulation runs from, and tf
% the final time. To generate figure 2 of Chen(2004), use t0=0 & tf=210
t0=0;
tf=210;

% "options" is used in the ODE solver below. The event here refers to the 
% reset rules, as defined in Chen(2004) page 3845.  
options=odeset('Events',@event_yeast,'RelTol',1e-4,'AbsTol',1e-6);

% Setting variables for the numerical integration
tout=t0; % tout contains the time values integrated over
yout=y0.'; % yout contains the conc. values for each variable in the ODEs
teout = []; % The "e" means events, te is the time  when the event takes place 
yeout = []; % ye is the y(value of all variables) when the event takes place.
ieout = []; % the "e" refers to the event."i" is the index of event. index ie=1 means one event  
while t0<tf   
    
    % This is the ODE solver
    [t,y,te,ye,ie]=ode15s(@eqns,[t0 tf],y0,options);
    
    nt=length(t);
    
    tout=[tout;t(2:nt)];  
    yout=[yout;y(2:nt, :)];
 
    
    teout = [teout;te]; % The "e" means events. "te" is the time  when the event takes place 
  	yeout = [yeout;ye];
  	ieout = [ieout;ie];
    
    
    
    % If any of the reset rules are applied (detected by the event in "options"
    % above), we refersh the initial conditions for ODEs 
    y0 = y(nt, :);

    	
     % the "e" refers to the event."i" is the index of event. index ie=1
     % means one event
 	if isscalar(ie) == 0 
    		ie = 0;
    end  
    
    % Here we apply the reset rules
 	if ie==1	                
       y0(1,nBUD)=0;
        y0(1,nSPN)=0;
        mdt=90;
        kg=log(2)/mdt;
        MorD=1;
        D=(1.026/kg)-32;
        
        lte1l=0.1;
        f=exp(-kg*D);
       
        y0(1,nMASS) =(abs(MorD+f-1))*y(nt,nMASS);
        y0(1,nlte1)=lte1l;
       
    elseif  ie==2                  
      	y0(1,nORI)=0;
      
    elseif ie==3
        ki20pp=8; 
        Bub2h=1;    
      	y0(1,nVi20)=ki20pp;
        y0(1,nBUB2)=Bub2h;
  
    elseif ie==4  
        ki20p=0.01;
        lte1h=1;
        Bub2l=0.2;
      	y0(1,nVi20)=ki20p;
        y0(1,nlte1)=lte1h;
        y0(1,nBUB2)=Bub2l;
   	
    end
    
    % After every event takes place, reinitialize the start time (t0) to
    % time step at event (t(nt)) and begin integration again
	t0 = t(nt);
    
	if t0 >= tf
		break;
	end

end

%% PLotting
% Here we plot the graphs of the concentrations of the variables over
% time, and reproduce Chen(2004) figure 2 (NB: only the plot of variables
% mass,BUD,ORI and SPN)

figure;

plot(tout,yout(:,nMASS),'b')
hold on
plot(tout,yout(:,nBUD),'r')
hold on
plot(tout,yout(:,nORI),'g')
hold on
plot(tout,yout(:,nSPN),'k')
hold on
legend('MASS','BUD','ORI','SPN')
title('Budding Yeast ODE Model Chen(2004)')
xlabel('Time')
ylabel('Concentration')
ylim([0 4])
xlim([0 tf])