% To generate curves presented in Figure 3 (A) of the main text, choose one of
% the values of the parameter 'a'
clear all
clc
%% System parameters with units in days
delta=0.4;          % fitness of drug-resistance
mu=1/(70*365);
gamma=1/(5);
R0=3;
beta=R0*(mu+gamma);
betar=delta*beta;
R3=delta*R0;
hh=0.000005;        % steps of increments for delay
%% parameter 'a'
%a=1.5;             % first curve in Fig 3 (A)
a=0.5;              % second curve in Fig 3 (A)


taustar=-log(1-1/R0)/(mu+gamma);         % delay when R1=1
tau0=-log(1-delta)/(mu+gamma);           % delay when R1=R3

%% curves of endemic state 
for ii=1:ceil(1/gamma/hh)
   xx(ii)=1-1/R3;
    tau(ii)=ii*hh;
    qqtau(ii)=2*exp(-a*tau(ii))/(1+exp(-2*(tau(ii)-1)));
    qtau=qqtau(ii);
     if tau(ii)<tau0
        Irstar=(mu/(mu+gamma))*(1-1/R3);
        total(ii)=Irstar;
     else
        R1=R0*(1-exp(-(mu+gamma)*tau(ii)));
        R2=qtau*R0*exp(-(mu+gamma)*tau(ii));
        Irstar=(mu/(mu+gamma))*(R1-1)/((exp((mu+gamma)*tau(ii))*(R1-R3)/qtau+R3));
        Istar=((R1-R3)/R2)*Irstar;
        total(ii)=Irstar+Istar;
     end
end
plot(tau,total,'linewidth',3)
hold on