% To generate curves presented in Figure 3 (B) of the main text, choose one of
% the values of delta (the fitness of drug-resistance)

clear all
clc
%% System parameters with units in days
delta=0.38;          % fitness of drug-resistance for first curve
%delta=0.42;         % fitness of drug-resistance for second curve
%delta=0.46;         % fitness of drug-resistance for third curve
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
%% initial conditions 
I0 = 0.001;
Ir0 = 0.0001;
S0 = 1-I0-Ir0;

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;
     if tau(ii)<tau0-0.06
        qqtau(ii)=0.6;
        qtau=qqtau(ii);
        Irstar=(mu/(mu+gamma))*(1-1/R3);
        R1=R0*(1-exp(-(mu+gamma)*tau(ii)));
        R2=qtau*R0*exp(-(mu+gamma)*tau(ii));
        total(ii)=Irstar;
    else
        qqtau(ii)=0.6*exp(-2*(tau(ii)-tau0+0.06));
        qtau=qqtau(ii);
        if tau(ii)<=tau0
        Irstar=(mu/(mu+gamma))*(1-1/R3);
        R1=R0*(1-exp(-(mu+gamma)*tau(ii)));
        R2=qtau*R0*exp(-(mu+gamma)*tau(ii));
        total(ii)=Irstar;
        else
        Irstar=(mu/(mu+gamma))*(R1-1)/((exp((mu+gamma)*tau(ii))*(R1-R3)/qtau+R3));
        R1=R0*(1-exp(-(mu+gamma)*tau(ii)));
        R2=qtau*R0*exp(-(mu+gamma)*tau(ii));
        Istar=((R1-R3)/R2)*Irstar;
        total(ii)=Irstar+Istar;
        end
     end
end
plot(tau,total,'linewidth',3)
hold on

  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  %%
  % Ttau=-gamma*log(1-1/R0)/(mu+gamma)
% EE=1/(1-exp(-(1+gamma/mu)))
% cond=1-1/R3
% qtau=0.3;
% R1=R0*(1-exp(-(mu+gamma)*tau))


  %% 
% 
% for ii=1:ceil(1/gamma/hh)
%         xx(ii)=1-1/R3;
%         tau(ii)=ii*hh;
% %         qqtau(ii)=1*exp(-0.2*tau(ii))/(1+exp(-2*(tau(ii)-taustar)));
% %         qtau=qqtau(ii);
%     if tau(ii)<taustar-0.01
%         qqtau(ii)=0.99;
%         qtau=qqtau(ii);
%         total(ii)=0;
%     else
%         qqtau(ii)=0.99*exp(-10*(tau(ii)-taustar+0.01));
%         qtau=qqtau(ii);
%         if tau(ii)<taustar
%             total(ii)=0;
%         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        
% end
%%
% Ni=3;
% hi=0.1;
% Nj=10;
% hj=1;
% rvec=1:hi:Ni;
% gvec=1:hj:Nj;
% nhi=length(rvec);
% nhj=length(gvec);
% curve=zeros(nhi,nhj);
% bif=zeros(nhi,nhj);
% R0=zeros(nhi,1);
% gg=zeros(nhj,1);
% for ii=1:nhi
%     R0(ii)=1+ii*hi;
%     for jj=1:nhj
%         gg(jj)=1/(jj*hj);
%         curve(ii,jj)=-log(1-1/R0(ii))-(mu+gg(jj))/gg(jj);
%         bif(ii,jj)=-log(1-1/R0(ii))/(mu+gamma);
%     end
% end