clc
clf
clear

%   All equations and values are taken from the article:
%   "Modelling the dynamics of the yeast pheromone pathway." 
%   by Kofahl et. al
%
%   This is the Main file that uses f.m for solving the ODEs.

%intial concentrations for the different species, all are expressed in nM
x0=zeros(41,1);
x0(1) = 100.0; %alpha
x0(2) = 1666.67; %Ste2
x0(3) = 0.0; %Ste2active
x0(4) = 1666.67; %Galphabetagamma
x0(5) = 0.0; %GalphaGTP
x0(6) = 0.0; %Gbetagamma
x0(7) = 0.0; %GaGDP
x0(8) = 235.72; %C
x0(9) = 0.0; %D
x0(10) = 158.33; %Ste5
x0(11) = 158.33; %Ste11
x0(12) = 105.94; %A
x0(13) = 36.4; %Ste7
x0(14) = 686.4; %Fus3
x0(15) = 77.87; %B
x0(16) = 1000; %Ste20
x0(17) = 0.0; %E
x0(18) = 0.0; %F
x0(19) = 0.0; %G
x0(20) = 0.0; %H
x0(21) = 0.0; %I
x0(22) = 0.0; %L
x0(23) = 0.0; %Fus3PP
x0(24) = 0.0; %K
x0(25) = 200; %Ste12
x0(26) = 0.0; %Ste12active
x0(27) = 200; %Bar1
x0(28) = 0.0; %Bar1active
x0(29) = 0.0; %Bar1activeex
x0(30) = 500; %Far1
x0(31) = 0.0; %Far1PP
x0(32) = 0.0; %Far1ubi
x0(33) = 0.0; %M
x0(34) = 0.0; %N
x0(35) = 300; %Cdc28
x0(36) = 0.0; %Sst2
x0(37) = 0.0; %Cre/Fus1
x0(38) = 0.0; %Cas9
x0(39) = 0.0; %Mated diploids
x0(40) = 1000; %Haploids
x0(41) = 0.0; %Cas9-gRNA complex


% %plots time courses for the different species
% voc=[0.011;0.006;0.052];
% x0(1)=voc(1);
% [t,x]=ode23s(@f,[0 100],x0);
% plot(t,x)

time=100;
voc=[0.011;0.006;0.052];

figure('Name','Protein expression for different octanal levels in breath (original yeast parameters, median values)')
subplot(3,2,1)

x0(1)=voc(1);
[t,x]=ode23s(@f,[0 time],x0);
plot(t,x(:,26),'r'),hold on
x0(1)=voc(2);
[t,x]=ode23s(@f,[0 time],x0);
plot(t,x(:,26),'m')
x0(1)=voc(3);
[t,x]=ode23s(@f,[0 time],x0);
plot(t,x(:,26),'b')
xlabel('time (mins)')
ylabel('Concentration of Ste12 (nM)')
legend('Control - 11pM','Smokers - 6pM','Lung cancer patients - 52pM')
legend('Location','northeast')
legend('boxoff')

% title('Cre expression for different octanal levels in breath (original yeast parameters, median values)')
subplot(3,2,2)

x0(1)=voc(1);
[t,x]=ode23s(@f,[0 time],x0);
plot(t,x(:,37),'r'),hold on
xcontrol=x(:,37);
x0(1)=voc(2);
[t,x]=ode23s(@f,[0 time],x0);
plot(t,x(:,37),'m')
x0(1)=voc(3);
[t,x]=ode23s(@f,[0 time],x0);
plot(t,x(:,37),'b')
xcancer=x(:,37);
xlabel('time (mins)')
ylabel('Concentration of Cre/Fus1 (nM)')
legend('Control - 11pM','Smokers - 6pM','Lung cancer patients - 52pM')
legend('Location','southeast')
legend('boxoff')

% ratio=xcancer/xcontrol;
% subplot(2,2,3)
% plot(t,ratio,'g')
% xlabel('time (mins)')
% ylabel('Cre expression ratio between cancer and control')


subplot(3,2,3)

x0(1)=voc(1);
[t,x]=ode23s(@f,[0 time],x0);
plot(t,x(:,38),'r'),hold on

x0(1)=voc(2);
[t,x]=ode23s(@f,[0 time],x0);
plot(t,x(:,38),'m')
x0(1)=voc(3);
[t,x]=ode23s(@f,[0 time],x0);
plot(t,x(:,38),'b')

xlabel('time (mins)')
ylabel('Concentration of Cas9/gRNA (nM)')
legend('Control - 11pM','Smokers - 6pM','Lung cancer patients - 52pM')
legend('Location','southwest')
legend('boxoff')
xlabel('time (mins)')
ylabel('Concentration of Cas9/gRNA (nM)')

% figure('Name','Mating and ADE2 deletion')
time2=1440;

subplot(3,2,4)
x0(1)=voc(1);
[t,x]=ode23s(@f,[0 time2],x0);
plot(t,x(:,39),'r'),hold on
x0(1)=voc(2);
[t,x]=ode23s(@f,[0 time2],x0);
plot(t,x(:,39),'m')
x0(1)=voc(3);
[t,x]=ode23s(@f,[0 time2],x0);
plot(t,x(:,39),'b')
xlabel('time (mins)')
ylabel('Ammount of mated diploids')
legend('Control - 11pM','Smokers - 6pM','Lung cancer patients - 52pM')
legend('Location','northwest')
legend('boxoff')

subplot(3,2,5)
x0(1)=voc(1);
[t,x]=ode23s(@f,[0 time2],x0);
plot(t,x(:,41),'r'),hold on
x0(1)=voc(2);
[t,x]=ode23s(@f,[0 time2],x0);
plot(t,x(:,41),'m')
x0(1)=voc(3);
[t,x]=ode23s(@f,[0 time2],x0);
plot(t,x(:,41),'b')
xlabel('time (mins)')
ylabel('Cas9:gRNA per diploid (nM)')
legend('Control - 11pM','Smokers - 6pM','Lung cancer patients - 52pM')
legend('Location','northwest')
legend('boxoff')