%% Model Parameters
L_wt=1.9; % growth rate of WT (hour^-1) 
K_wt=0.8; % carrying capacity of WT/concentration of saturated culture of WT
D_wt=0.07; % killing rate of WT at the end (lysis-growth, hour^-1)

T=1; % difference in lag time (hours)
D=26; % killing time at the end (hours)

L_m = 1.6;% growth rate of mutant (hour^-1) 
K_m = 1.2;% carrying capacity of mutant
D_m = 0.07;% killing rate of mutant at the end (hour^-1)

p = [L_wt K_wt L_m K_m D_wt D_m]; %define vector with model parameters

%% Evolution Simulation
xx=[]; %initialize empty vector to store time
tt=[]; %initialize empty vector to store OD

xx=[0.8 10^(-7)]';
tt=0;
for i=1:7
    if T>0
        [tt1,xx1]=ode113(@(t,x) eq1(t,x,p),(i-1)*48+[0 T],xx(:,end)'/200);
    else
        tt1=[];
        xx1=[];
    end
   
    [tt2,xx2]=ode113(@(t,x) eq2(t,x,p),(i-1)*48+[T 48-D],xx1(end,:));
    
    if D>0
        [tt3,xx3]=ode113(@(t,x) eq3(t,x,p),(i-1)*48+[48-D 48],xx2(end,:));
    else
        tt3=[];
        xx3=[];
    end
    xx=[xx xx1' xx2' xx3'];
    tt=[tt tt1' tt2' tt3'];
end

tt=tt/24; % convert time to days
xxx=smoothdata(xx(2,:)./(xx(1,:)+xx(2,:)),'movmedian',50); % define mutant fraction

%% Plots
figure
plot(tt,xx(1,:),'g','LineWidth',2)
hold on
plot(tt,xx(2,:),'r','LineWidth',2)
xlabel('days','FontSize',15)
ylabel('OD (each strain)','FontSize',15)
axis([0 tt(end) 0 1.2])
legend('WT','\Delta{\itlasR}','Location','northwest')

figure
plot(tt,xxx,'--k','LineWidth',2)
xlabel('days','FontSize',15)
ylabel('fraction of \Delta{\itlasR}','FontSize',15)

%% Model
% lag phase model
% no competition
% WT population constant, exponential growth of mutants
function dx=eq1(t,x,p)
    dx=zeros(2,1);
    dx(1)=0;
    dx(2)=p(3).*x(2);
end

% growth phase model
% competition 
% logistic growth of both WT and mutants
function dx=eq2(t,x,p)
    dx=zeros(2,1);
    dx(1)=max(0,p(1).*x(1).*(1-(x(1)+(p(3)./p(1)).*x(2))./p(2)));
    dx(2)=max(0,p(1).*x(2).*((p(3)./p(1))-(x(1)+(p(3)./p(1)).*x(2))./p(4)));
end

% death phase model
% exponential decay of both WT and mutants
function dx=eq3(t,x,p)
    dx=zeros(2,1);
    dx(1)=-p(5).*x(1);
    dx(2)=-p(6).*x(2);
end

