Drawing contour plots of R0 vs two given parameters

Question

I have been working on drawing the contour plots of the basic reproduction number (R0) of a model with respect to two parameters. For practice, I am currently working on replicating the graph attached below from a paper by Kifle et al, https://www.sciencedirect.com/science/article/pii/S2211379722000122.

The MATLAB codes that I used to calculate R0 of the model are given below.

%Parameters 

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

%Basic Reproduction Number

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
K_5 = gamma_3 + delta_3 + 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)

This is what I tried so far and it doesn't seem right.

[beta,eta_1] = meshgrid(0.1:0.001:1,0.1:0.001:1);

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)


%Drawing the plot

surf(beta,eta_1,R_0)
hold on
z2 = 0*beta + 1
surf(beta,eta_1,z2,'MarkerFaceColor','red')

I would be very much grateful, if someone could please help me draw the contour plots using MATLAB or R.

Answer 1

I think if I were writing the function in R, I would allow the other arguments as parameters to allow tweaking:

R0 <- function(beta, eta_1, theta = 141302, mu = 0.001229, tau = 0.45,
               zeta = 1/14, alpha = 0.75214, q = 0.31167, eta_2 = 0.02557,
               eta_3 = 1/7, delta_1 = 0.16673, delta_2 = 0.00147, 
               delta_3 = 0.00038, gamma_1 = 0.00827, gamma_2 = 0.00787, 
               gamma_3 = 0.20186) {
  
    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
    K_5 <- gamma_3 + delta_3 + mu

    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)
}

You can then recreate the above plot like this:

library(dplyr)
library(geomtextpath)

expand.grid(beta = seq(0, 1, 0.01), eta1 = seq(0, 1, 0.01)) %>%
  mutate(R0 = apply(., 1, function(x) R0(x[1], x[2]))) %>%
  ggplot(aes(beta, eta1)) +
  geom_contour_filled(aes(z = R0)) +
  geom_textcontour(aes(z = R0, label = after_stat(level)), size = 6) +
  scale_x_continuous(name = bquote(beta), expand = c(0, 0)) +
  scale_y_continuous(name = bquote(eta[1]), expand = c(0, 0)) +
  coord_equal() +
  theme_classic(base_size = 20) +
  scale_fill_manual(values = c("red", "yellow", "blue", "gray", "#00bfbf",
                               "purple", "red", "yellow"), name = bquote(R[0])) +
  guides(fill = guide_colorsteps(key_height = unit(10, "mm"))) +
  theme(legend.key.height = unit(15, "mm"))