%Produces game theory plots. clear clc pop = 100000; t_MAX = 1000; beta1 = 0.37; beta2 = 0.37; mu1 = 1/(365.25*76.4); mu2 = 1/20; alpha1 = pop*mu1; phi = 0.97; R0 = 4.072727272727272727272727272727; lambda1 = 0.5; lambda2 = 0.5; ggamma = 1/7; N = alpha1/mu1; alpha2 = (R0^2)*(mu2^2)*alpha1*(ggamma+mu1)*(lambda1+mu1)*(lambda2+mu2)/(mu1*beta1*beta2*lambda1*lambda2); y_MAX_h = alpha1/mu1; y_MAX_m = alpha2/mu2; rHI = 1-((alpha1*mu2^2*(ggamma+mu1)*(lambda1+mu1)*(lambda2+mu2))/(alpha2*beta1*beta2*lambda1*lambda2*mu1)); rNE = 0:rHI/299:rHI; C = zeros(length(rNE)); R0var2 = zeros(length(rNE)); for i=1:length(rNE) beta1r = (1-rNE(i))*beta1; E = (alpha1*alpha2*beta1r*beta2*lambda1*lambda2-mu1*mu2*mu2*N*N*(lambda1+mu1)*(lambda2+mu2)*(ggamma+mu1))/(alpha2*beta1r*beta2*lambda1*lambda2*(lambda1+mu1)+beta2*lambda1*mu1*mu2*N*(lambda1+mu1)*(lambda2+mu2)); Z = (alpha2*beta2*lambda1*lambda2*E)/(mu2*(lambda2+mu2)*(beta2*lambda1*E+mu2*N*(ggamma+mu1))); C(i) = beta1*Z/(mu1*N+beta1*Z); R0var2(i) = (1/mu2)*sqrt(alpha2*beta1r*beta2*lambda1*lambda2*mu1/(alpha1*(ggamma+mu1)*(lambda1+mu1)*(lambda2+mu2))); end %What cost provides herd immunity? beta1r = (1-rHI)*beta1; E = (alpha1*alpha2*beta1r*beta2*lambda1*lambda2-mu1*mu2*mu2*N*N*(lambda1+mu1)*(lambda2+mu2)*(ggamma+mu1))/(alpha2*beta1r*beta2*lambda1*lambda2*(lambda1+mu1)+beta2*lambda1*mu1*mu2*N*(lambda1+mu1)*(lambda2+mu2)); Z = (alpha2*beta2*lambda1*lambda2*E)/(mu2*(lambda2+mu2)*(beta2*lambda1*E+mu2*N*(ggamma+mu1))); CHI = beta1*Z/(mu1*N+beta1*Z); %plot C against rNE figure box on hold on set(gcf,'Units','Inches','Position',[0, 0, 3, 2.5]); plot(C(2:length(C)),rNE(2:length(rNE)),'k','LineWidth',1.5,'MarkerSize',3) plot(linspace(0,1),rHI*ones(length(linspace(0,1))),'k--','LineWidth',1,'MarkerSize',1) xlabel('relative cost of protection $C$','Interpreter','latex') xticks([0 0.2 0.4 0.6 0.8 max(max(C)) 1]) xticklabels({'$0$','$0.2$','$0.4$','$0.6$','$0.8$','$C_\mathrm{max}$',''}) set(gca, 'FontSize', 10) ylabel('protection level $r_\mathrm{NE}$','Interpreter','latex') yticks([0 0.2 0.4 0.6 0.8 rHI 1]) yticklabels({'$0$','$0.2$','$0.4$','$0.6$','$0.8$','$r_\mathrm{HI}$','$1$'}) set(gca,'TickLabelInterpreter', 'latex'); plot([max(max(C)) 1],[0 0],'k','LineWidth',1.5,'MarkerSize',3) %plot R0 against rpop with rHI rpop = 0:(1/299):1; R0var = zeros(length(rpop)); for i = 1:length(rpop) R0var(i) = (1/mu2)*sqrt(alpha2*beta1*(1-rpop(i))*beta2*lambda1*lambda2*mu1/(alpha1*(ggamma+mu1)*(lambda1+mu1)*(lambda2+mu2))); end figure set(gcf,'Units','Inches','Position',[0, 0, 3, 2.5]); box on hold on plot(rpop,R0var(:,1),'b','LineWidth',1.5,'MarkerSize',3) plot(rHI*ones(length(linspace(0,1))),linspace(0,1),'r--','LineWidth',1,'MarkerSize',1) plot(linspace(0,1),ones(length(linspace(0,1))),'r--','LineWidth',1,'MarkerSize',1) xlabel('protection level $r_\mathrm{pop}$','Interpreter','latex') xticks([0 0.2 0.4 0.6 0.8 rHI 1]) xticklabels({'$0$', '$0.2$', '$0.4$', '$0.6$', '$0.8$', '$r_\mathrm{HI}$', ''}) set(gca, 'FontSize', 10) ylabel('basic reproduction number $R_0$','Interpreter','latex') yticks([0 1 2 3 4 R0]) yticklabels({'$0$','$1$','$2$','$3$','$4$'}) set(gca,'TickLabelInterpreter', 'latex'); axis([0 1 0 R0]) %plot R0 against C figure set(gcf,'Units','Inches','Position',[0, 0, 3, 2.5]); box on hold on plot(C,R0var2,'b','LineWidth',1.5,'MarkerSize',3) plot(linspace(0,1),ones(length(linspace(0,1)),1),'r--','LineWidth',1,'MarkerSize',1) xlabel('relative cost of protection $C$','Interpreter','latex') xticks([0 0.2 0.4 0.6 0.8 max(max(C)) 1]) xticklabels({'$0$','$0.2$','$0.4$','$0.6$','$0.8$','$C_\mathrm{max}$',''}) set(gca, 'FontSize', 10) ylabel('basic reproduction number $R_0$','Interpreter','latex') yticks([0 1 2 3 4 R0]) yticklabels({'$0$','$1$','$2$','$3$','$4$',''}) set(gca,'TickLabelInterpreter', 'latex'); axis([0 1 0 R0]) plot([max(max(C)) 1],[R0 R0],'b','LineWidth',1.5,'MarkerSize',3)