%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%% Figures 4 and S3 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % Chisholm et al., Implications of asymptomatic carriers % % for infectious disease transmission and control %%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% set(0,'DefaultTextFontname', 'Times') set(0,'DefaultTextFontSize', 20) set(0,'DefaultAxesFontName', 'Times') set(0,'DefaultAxesFontSize', 20) close all clear all T = 2000; N = 1000; gamma = 0.1; %%%%%%%%%%%%%%%%%%%%%%%%%%% %%% 4(a) and S3(a)--(c) %%% %%%%%%%%%%%%%%%%%%%%%%%%%%% RR = [0.75 1.2]; omega0 = [0.1 0.3]; alpha = 0.1; beta = 0.12; tau = 0.1; xi = 0.6; eta0 = RR * xi; for j = 1 : 2 eta = eta0(j); R = RR(j); for i = 1 : 2 omega = omega0(i); f = @(t,x) [beta * (x(2) + eta * x(1)) * (1 - alpha) * (N - x(1) - x(2)) / N + ... omega * x(2) - x(1) * (xi * gamma + tau); ... beta * (x(2) + eta * x(1)) * alpha * (N - x(1) - x(2)) / N + ... tau * x(1) - x(2) * (gamma + omega)]; if i == 1 [t,xa] = ode45(@(t,x) f(t,x),[0 T],[10 0]); xa(xa<0)=0; Chatpre(j) = xa(end,1)/N; Ihatpre(j) = xa(end,2)/N; Phatpre(j) = Chatpre(j) + Ihatpre(j); xaa = xa; ta = t; else [t,xa] = ode45(@(t,x) f(t,x),[T 2*T],[xa(end,1) xa(end,2)]); xa(xa<0)=0; Chatpost(j) = xa(end,1)/N; Ihatpost(j) = xa(end,2)/N; Phatpost(j) = Chatpost(j) + Ihatpost(j); xaa = [xaa; xa]; ta = [ta; t]; end end if j == 1 figure(1) line(ta,(xaa(:,2)+xaa(:,1))/N,'LineWidth',2); hold on else figure(1) line(ta,(xaa(:,2)+xaa(:,1))/N,'LineWidth',2,'LineStyle', ':') hold on end end xlabel('Time') title('Change \omega') axis([1 4000 0 0.35]) ylabel('P') legend(['R = ',num2str(RR(1))],[' R = ',num2str(RR(2))]) Y = [Chatpre, Ihatpre, Phatpre; Chatpost, Ihatpost, Phatpost]; figure(2) subplot(1,3,1) bar([Chatpre(1) Chatpost(1) 0 Chatpre(2) Chatpost(2)]) Labels = {'B', 'A','','B','A'}; set(gca, 'XTick', 1:5, 'XTickLabel', Labels); ylabel('P_C') axis([0 6 0 0.35]) subplot(1,3,2) bar([Ihatpre(1) Ihatpost(1) 0 Ihatpre(2) Ihatpost(2)]) Labels = {'B', 'A','','B','A'}; set(gca, 'XTick', 1:5, 'XTickLabel', Labels); ylabel('P_I') title('Change \omega') axis([0 6 0 0.35]) subplot(1,3,3) bar([Phatpre(1) Phatpost(1) 0 Phatpre(2) Phatpost(2)]) Labels = {'B', 'A','','B','A'}; set(gca, 'XTick', 1:5, 'XTickLabel', Labels); ylabel('P') axis([0 6 0 0.35]) %%%%%%%%%%%%%%%%%%%%%%%%%%% %%% 4(b) and S3(d)--(f) %%% %%%%%%%%%%%%%%%%%%%%%%%%%%% RR = [0.7 1.2]; tau0 = [0.06 0.02]; alpha = 0.1; beta = 0.13; omega = 0.1; xi = 0.6; eta0 = RR * xi; for j = 1 : 2 eta = eta0(j); R = RR(j); for i = 1 : 2 tau = tau0(i); f = @(t,x) [beta * (x(2) + eta * x(1)) * (1 - alpha) * (N - x(1) - x(2)) / N + ... omega * x(2) - x(1) * (xi * gamma + tau); ... beta * (x(2) + eta * x(1)) * alpha * (N - x(1) - x(2)) / N + ... tau * x(1) - x(2) * (gamma + omega)]; if i == 1 [t,xa] = ode45(@(t,x) f(t,x),[0 T],[10 0]); xa(xa<0)=0; Chatpre(j) = xa(end,1)/N; Ihatpre(j) = xa(end,2)/N; Phatpre(j) = Chatpre(j) + Ihatpre(j); xaa = xa; ta = t; else [t,xa] = ode45(@(t,x) f(t,x),[T 2*T],[xa(end,1) xa(end,2)]); xa(xa<0)=0; Chatpost(j) = xa(end,1)/N; Ihatpost(j) = xa(end,2)/N; Phatpost(j) = Chatpost(j) + Ihatpost(j); xaa = [xaa; xa]; ta = [ta; t]; end end if j == 1 figure(3) line(ta,(xaa(:,2)+xaa(:,1))/N,'LineWidth',2); hold on else figure(3) line(ta,(xaa(:,2)+xaa(:,1))/N,'LineWidth',2,'LineStyle', ':') hold on end end xlabel('Time') title('Change \tau') axis([1 4000 0 0.35]) ylabel('P') legend(['R = ',num2str(RR(1))],[' R = ',num2str(RR(2))]) Y = [Chatpre, Ihatpre, Phatpre; Chatpost, Ihatpost, Phatpost]; figure(4) subplot(1,3,1) bar([Chatpre(1) Chatpost(1) 0 Chatpre(2) Chatpost(2)]) Labels = {'B', 'A','','B','A'}; set(gca, 'XTick', 1:5, 'XTickLabel', Labels); ylabel('P_C') axis([0 6 0 0.35]) subplot(1,3,2) bar([Ihatpre(1) Ihatpost(1) 0 Ihatpre(2) Ihatpost(2)]) Labels = {'B', 'A','','B','A'}; set(gca, 'XTick', 1:5, 'XTickLabel', Labels); ylabel('P_I') title('Change \tau') axis([0 6 0 0.35]) subplot(1,3,3) bar([Phatpre(1) Phatpost(1) 0 Phatpre(2) Phatpost(2)]) Labels = {'B', 'A','','B','A'}; set(gca, 'XTick', 1:5, 'XTickLabel', Labels); ylabel('P') axis([0 6 0 0.35]) %%%%%%%%%%%%%%%%%%%%%%%%%%% %%% 4(c) and S3(g)--(i) %%% %%%%%%%%%%%%%%%%%%%%%%%%%%% RR = [0.70 1.4]; alpha0 = [0.7 0.2]; omega = 0.1; beta = 0.12; tau = 0.06; xi = 0.5; eta0 = RR * xi; for j = 1 : 2 eta = eta0(j); R = RR(j); for i = 1 : 2 alpha = alpha0(i); f = @(t,x) [beta * (x(2) + eta * x(1)) * (1 - alpha) * (N - x(1) - x(2)) / N + ... omega * x(2) - x(1) * (xi * gamma + tau); ... beta * (x(2) + eta * x(1)) * alpha * (N - x(1) - x(2)) / N + ... tau * x(1) - x(2) * (gamma + omega)]; if i == 1 [t,xa] = ode45(@(t,x) f(t,x),[0 T],[10 0]); xa(xa<0)=0; Chatpre(j) = xa(end,1)/N; Ihatpre(j) = xa(end,2)/N; Phatpre(j) = Chatpre(j) + Ihatpre(j); xaa = xa; ta = t; else [t,xa] = ode45(@(t,x) f(t,x),[T 2*T],[xa(end,1) xa(end,2)]); xa(xa<0)=0; Chatpost(j) = xa(end,1)/N; Ihatpost(j) = xa(end,2)/N; Phatpost(j) = Chatpost(j) + Ihatpost(j); xaa = [xaa; xa]; ta = [ta; t]; end end if j == 1 figure(5) line(ta,(xaa(:,2)+xaa(:,1))/N,'LineWidth',2); hold on else figure(5) line(ta,(xaa(:,2)+xaa(:,1))/N,'LineWidth',2,'LineStyle', ':') hold on end end title('Change \alpha') xlabel('Time') axis([1 4000 0 0.35]) ylabel('P') legend(['R = ',num2str(RR(1))],[' R = ',num2str(RR(2))]) Y = [Chatpre, Ihatpre, Phatpre; Chatpost, Ihatpost, Phatpost]; figure(6) subplot(1,3,1) bar([Chatpre(1) Chatpost(1) 0 Chatpre(2) Chatpost(2)]) Labels = {'B', 'A','','B','A'}; set(gca, 'XTick', 1:5, 'XTickLabel', Labels); ylabel('P_C') axis([0 6 0 0.35]) subplot(1,3,2) bar([Ihatpre(1) Ihatpost(1) 0 Ihatpre(2) Ihatpost(2)]) Labels = {'B', 'A','','B','A'}; set(gca, 'XTick', 1:5, 'XTickLabel', Labels); ylabel('P_I') title('Change \alpha') axis([0 6 0 0.35]) subplot(1,3,3) bar([Phatpre(1) Phatpost(1) 0 Phatpre(2) Phatpost(2)]) Labels = {'B', 'A','','B','A'}; set(gca, 'XTick', 1:5, 'XTickLabel', Labels); ylabel('P') axis([0 6 0 0.35])