Plotting contour plots of R0 against two parameters

I have been trying to replicate the following graph from the paper, https://www.sciencedirect.com/science/article/pii/S2211379722000122 for a class project.

The authors draw contour plots of R0 against beta and zeta. The parameter values and the formula for R0 are given below.

%Basic Reproduction Number

theta = 141302; %recruitment rate
mu = 0.001229; %natural death rate
tau = 0.45; %modification factor for A
zeta = 1/14; %influx from Q to S 
beta = 0.88; %transmission coefficient
alpha = 0.75214; %hospitalization rate
q = 0.31167; %influx from Q to I
eta_1 = 0.81692; %influx from E to Q
eta_2 = 0.02557; %influx from E to A
eta_3 = 1/7; %influx from E to I
delta_1 = 0.16673; %disease death rate for A
delta_2 = 0.00147; %disease death rate for I
delta_3 = 0.00038; %disease death rate for J
gamma_1 = 0.00827; %recovery rate for A
gamma_2 = 0.00787; %recovery rate for I
gamma_3 = 0.20186; %recovery rate for J

K_1 = eta_1 + eta_2 + eta_3 + mu
K_2 = zeta + q + mu
K_3 = gamma_1 + delta_1 + mu
K_4 = alpha + gamma_2 + delta_2 + mu

R_0 = beta*(tau*eta_2*K_2*K_4 + K_3*(eta_3*K_2 + eta_1*q))/(K_1*K_2*K_3*K_4)

Any support on drawing the plot would be highly appreciated. Thank you for your time!

Answer 1

[Mathieu NOE]

hello

try this

now beta and zeta must be vectors so that R_0 is a matrix (2D array)

use contourf for the plot

I let you choose which colormap fits best your color taste !

hope it helps

%Basic Reproduction Number
theta = 141302; %recruitment rate
mu = 0.001229; %natural death rate
tau = 0.45; %modification factor for A
% zeta = 1/14; %influx from Q to S 
zeta = linspace(0,1,100);
% beta = 0.88; %transmission coefficient
beta = linspace(-0.001,1,100)'; % NB the min value must be slightly negative to let the contour line "0" show up
alpha = 0.75214; %hospitalization rate
q = 0.31167; %influx from Q to I
eta_1 = 0.81692; %influx from E to Q
eta_2 = 0.02557; %influx from E to A
eta_3 = 1/7; %influx from E to I
delta_1 = 0.16673; %disease death rate for A
delta_2 = 0.00147; %disease death rate for I
delta_3 = 0.00038; %disease death rate for J
gamma_1 = 0.00827; %recovery rate for A
gamma_2 = 0.00787; %recovery rate for I
gamma_3 = 0.20186; %recovery rate for J
K_1 = eta_1 + eta_2 + eta_3 + mu;
K_2 = zeta + q + mu;
K_3 = gamma_1 + delta_1 + mu;
K_4 = alpha + gamma_2 + delta_2 + mu;
R_0 = beta.*(tau.*eta_2*K_2.*K_4 + K_3.*(eta_3.*K_2 + eta_1.*q))./(K_1.*K_2.*K_3.*K_4);

% plot
levels = (0:0.2:1.2);
[c,h] = contourf(beta,zeta,R_0',levels);
clabel(c,h)
h.LineWidth = 2; % contour line width (up to you)
xlabel('\beta');
ylabel('\zeta');
colormap(jet(numel(levels)-1));
colorbar('vert');

[Bas123]

This was perfect, thank you so much!

[Mathieu NOE]

my pleasure !