Skip to contents

COVID-19 control with a larger NPI space

multiple_NPIs.Rmd

Overview

This vignette demonstrates model predictive optimal control of ICU cases for COVID-19

Load Libraries and Setup

Load the required libraries and initialize the simulation environment.

library(EpiControl)
library(VGAM)
library(parallel)
library(pbapply)
library(zoo)  # For rolling sum operations
library(dplyr)
library(ggplot2)

Initialize the cluster for parallel computing:

cores <- detectCores() - 1
cl <- makeCluster(cores)

clusterSetRNGStream(cl, iseed = 20250501)

Epidemiological and Noise Parameters

Define parameters for pathogens and noise: Can include different pathogens, in this case we consider COVID-19 and Ebola virus.

Epi_pars <- data.frame(
  Pathogen = c("COVID-19", "Ebola"),
  R0 = c(3.5, 2.5),
  gen_time = c(6.5, 15.0),
  gen_time_var = c(2.1, 2.1),
  CFR = c(0.0132, 0.5),
  mortality_mean = c(10.0, 10.0),
  mortality_var = c(1.1, 1.1)
)

Noise_pars <- data.frame(
  repd_mean = 10.5,  # Reporting delay mean
  del_disp = 5.0,    # Reporting delay variance
  ur_mean = 0.3,     # Under-reporting mean
  ur_beta_a = 50.0   # Beta distribution alpha for under-reporting
)

Action Space for Policy Interventions

Define non-pharmaceutical interventions: We have 2 actions, introducing a “lockdown” and no interventions.

# we can use 3, 6 or 10 levels
#Action_space <- data.frame (
#  NPI = c("No restrictions", "Social distancing", "Lockdown"),
#  R_coeff = c(1.0, 0.5, 0.2), #R0_act = R0 * ctrl_states
#  R_beta_a = c(0.0, 5.0, 5.0), #R0_act uncertainty
#  cost_of_NPI = c(0.0, 0.01, 0.15)
#)

#Action_space <- data.frame (
#  NPI = c("No restrictions", "SD1", "SD2", "SD3", "SD4", "Lockdown"),
#  R_coeff = c(1.0, 0.8, 0.6, 0.4, 0.2, 0.1), #R0_act = R0 * ctrl_states
#  R_beta_a = c(0.0, 5.0, 5.0, 5.0, 5.0, 5.0), #R0_act uncertainty
#  cost_of_NPI = c(0.0, 0.02, 0.05, 0.08, 0.15, 0.35)
#)

Action_space <- data.frame (
  NPI = c("No restrictions", "SD1", "SD2", "SD3", "SD4", "SD5", "SD6", "SD7", "SD8", "Lockdown"),
  R_coeff = c(1.0, 0.9, 0.8, 0.7, 0.6, 0.5, 0.4, 0.3, 0.2, 0.1), #R0_act = R0 * ctrl_states
  R_beta_a = c(0.0, 5.0, 5.0, 5.0, 5.0, 5.0, 5.0, 5.0, 5.0, 5.0), #R0_act uncertainty
  cost_of_NPI = c(0.0, 0.01, 0.02, 0.03, 0.04, 0.05, 0.065, 0.08, 0.15, 0.35)
)

Simulation Setup

Set up key parameters and initialize simulation data:

C_target <- 5000
D_target <- 1
r_trans_len <- 7

sim_settings <- list(
  ndays = 121 * 7, #simulation length
  start_day = 1,
  N = 3.5e6, # population size
  I0 = 10, # initial infections
  C_target = C_target, #target cases
  C_target_pen = C_target*1.5, #overshoot penalty threshold
  R_target = 1.0,
  D_target = D_target, #one way to get peaks at 400 is to increase this to 15
  D_target_pen = 50, #max death
  alpha = 1.3/C_target, #~proportional gain (regulates error in cases) covid
  #alpha = 3.25/C_target #~proportional gain (regulates error in cases) ebola
  alpha_d = 0*1.3/D_target,
  ovp = 5.0, #overshoot penalty
  dovp = 0*10.0, #death overshoot penalty
  gamma = 0.95, #discounting factor
  n_ens = 20L, #MC assembly size for 4
  sim_ens = 20L, #assembly size for full simulation
  rf = 7L, #days 14
  R_est_wind = 5L, #rf-2 #window for R estimation
  pred_days = 12L,
  r_trans_steep = 1.5,  # Growth rate
  r_trans_len = r_trans_len,  # Number of days for the transition
  t0 = r_trans_len  / 2, # Midpoint of the transition
  pathogen = 1,
  susceptibles = 1,
  delay = 0,
  ur = 0,
  r_dir = 2
)

column_names <- c("days", "sim_id", "I", "Lambda", "C", "Lambda_C", "S", "Deaths", "Re", "Rew", "Rest", "R0est", "policy", "R_coeff")
episim_data <- data.frame(matrix(0, nrow = (sim_settings$ndays), ncol = length(column_names)))
colnames(episim_data) <- column_names

episim_data$policy <- rep(1, sim_settings$ndays)
episim_data$days <- 1:(sim_settings$ndays)
episim_data[1, ] <- c(1, 1, sim_settings$I0, sim_settings$I0, Noise_pars$ur_mean * sim_settings$I0, Noise_pars$ur_mean * sim_settings$I0, sim_settings$N - sim_settings$I0, 0, Epi_pars[1, "R0"], Epi_pars[1, "R0"], 1, 1, 1, 1)

episim_data_ens <- replicate(sim_settings$sim_ens, episim_data, simplify = FALSE)

for (ii in 1:sim_settings$sim_ens) {
  episim_data_ens[[ii]]$sim_id <- rep(ii, sim_settings$ndays)
}

Running the Simulation

Run simulations in parallel using the pblapply function:


episettings <- list(
  sim_function = Epi_MPC_run_wd,
  reward_function = reward_fun_wd,
  R_estimator = R_estim,
  noise_par = Noise_pars,
  epi_par = Epi_pars,
  actions = Action_space,
  sim_settings = sim_settings,
  parallel = TRUE
)

episettings$cl <- cl

results <- epicontrol(episim_data_ens, episettings)
#> [1] 1
#>   |                                                          |                                                  |   0%

episim_data_ens <- results
stopCluster(cl)

#for (jj in 1:sim_ens) {
#  episim_data_ens[[jj]] <- Epi_MPC_run_wd(episim_data_ens[[jj]], Epi_pars, Noise_pars, Action_space, pred_days = pred_days, n_ens = n_ens, ndays = nrow(episim_data), R_est_wind = R_est_wind, pathogen = 1, susceptibles = 0, delay = 0, ur = 0, r_dir = 2, N = N)
#  setTxtProgressBar(pb,jj)
#}
#close(pb)

for (jj in 1:sim_settings$sim_ens) {
  episim_data_ens[[jj]] <- head(episim_data_ens[[jj]], -sim_settings$pred_days)
  episim_data_ens[[jj]]["D_roll"] <- rollsum(episim_data_ens[[jj]]["Deaths"], 7, fill = NA)
  episim_data_ens[[jj]]["I_roll"] <- rollsum(episim_data_ens[[jj]]["I"], 7, fill = NA)
  episim_data_ens[[jj]]["D_cum"] <- cumsum(episim_data_ens[[jj]]["Deaths"])
  episim_data_ens[[jj]]["I_cum"] <- cumsum(episim_data_ens[[jj]]["I"])
}

# Combine Simulation Results
combined_data <- do.call(rbind, episim_data_ens)

Plotting results 


combined_data <- combined_data %>%
  group_by(sim_id, policy) %>%
  arrange(days) %>%
  mutate(group = cumsum(c(1, diff(days) != 1))) %>%
  ungroup()

combined_data <- combined_data %>%
  arrange(sim_id, days, policy)

policy_labels <- c("1" = "No restrictions", "2" = "SD1", "3" = "SD2", "4" = "SD3", "5" = "SD4", "6" = "SD5", "7" = "SD6", "8" = "SD7")

# find the day when herd immunity reached

# Sample constant value
HI <- 1/Epi_pars[1,"R0"]

# Assuming episim_data_ens[[1]]["S"] is a dataframe or list
# Extract the vector
S_vector <- episim_data_ens[[1]]["S"][[1]]

# Compute the fraction S/N
fraction <- S_vector / sim_settings$N

# Find the index where the fraction is first smaller than the constant
index <- which(fraction < HI)[1]


ggplot(combined_data %>% filter(sim_id == 1)) +
  geom_line(data = subset(combined_data, sim_id != 1), aes(x = days, y = C, color = as.factor(sim_id)), alpha = 0.1) +
  geom_line(aes(x = days, y = C, color = factor(policy, labels = policy_labels), group = 1), alpha = 1.0) +
  geom_line(aes(x = days, y = I, color = factor(policy, labels = policy_labels), group = 1), alpha = 1.0, size=0.25) +
  geom_hline(yintercept = C_target, linetype = "dashed", color = "blue", size=0.25) +
  geom_vline(xintercept = index, linetype = "dashed", color = "black", size=0.25) +
  labs(x = "Days", y = "Reported cases and true infections", color = "Policy") +
  scale_color_manual(values = c("SD1" = "aquamarine", "SD2" = "darkgoldenrod1", "SD3" = "deeppink", "SD4" = "darkorchid1", "SD5" = "purple","SD6" = "red","SD7" = "darkred",  "No restrictions" = "chartreuse3")) +
  #scale_color_manual(values = c("SD1" = "aquamarine", "SD2" = "darkgoldenrod1", "SD3" = "deeppink", "SD4" = "darkorchid1", "SD5" = "purple", "No restrictions" = "chartreuse3")) +
  guides(color = guide_legend(title = "Policy"))+
  ylim(0,10000)