Model changing booster coverage over time using set_pev_epi

[lmhaile]

We would like to implement a vaccine strategy where coverage of the booster dose varies over time. set_pev_epi is currently set up to permit varying coverage of the initial three vaccine doses over time, but only flat coverage of the booster dose.

In this example scenario, I set booster coverage to kick in 12 months after the initial series:

# Set colour palette:
cols <- c("#E69F00", "#56B4E9", "#009E73", "#F0E442", "#0072B2", "#D55E00", "#CC79A7")

plot_doses <- function(){
  output$month <- ceiling(output$timestep / month)
  doses <- output[, c(grep("n_pev" , names(output)), grep("month", names(output)))]
  doses <- aggregate(cbind(doses[1:4]),
                     by = list(doses$month), 
                     FUN = sum)
  doses <- as.matrix(t(doses[, -1]))
  
  barplot(doses, xlab = "Month",
          ylab = "Number of doses",
          col = cols[1:6], space = 0, axes = T,
          beside = FALSE, xaxs = "i", yaxs = "i",
          ylim = c(0, 230))
  grid(lty = 2, col = "grey80", lwd = 0.5);box()
  axis(side = 1, lty = 1, col = "black", pos = 0)
  legend("topleft", box.lty = 0, legend = c("Dose 1","Dose 2","Dose 3","Booster 1"),
         fill = cols[1:6], bg="transparent", cex = 0.8, y.intersp = 1.5)
}


plot_doses_multiple_boosters <- function(){
  output$month <- ceiling(output$timestep / month)
  doses <- output[, c(grep("n_pev" , names(output)), grep("month", names(output)))]
  doses <- aggregate(cbind(doses[1:5]),
                     by = list(doses$month), 
                     FUN = sum)
  doses <- as.matrix(t(doses[, -1]))
  
  barplot(doses, xlab = "Month",
          ylab = "Number of doses",
          col = cols[1:6], space = 0, axes = T,
          beside = FALSE, xaxs = "i", yaxs = "i",
          ylim = c(0, 230))
  grid(lty = 2, col = "grey80", lwd = 0.5);box()
  axis(side = 1, lty = 1, col = "black", pos = 0)
  legend("topleft", box.lty = 0, legend = c("Dose 1","Dose 2","Dose 3","Booster 1", "Booster 2"),
         fill = cols[1:6], bg="transparent", cex = 0.8, y.intersp = 1.5)
}

# initial vaccine scenario  ----------------------------------------------------
year <- 365
month <- 30
sim_length <- 3 * year
human_population <- 10000
starting_EIR <- 20

simparams <- get_parameters(list(
  human_population = human_population,
  clinical_incidence_rendering_min_ages = 0,
  clinical_incidence_rendering_max_ages = 5 * year,
  individual_mosquitoes = FALSE
)
)

simparams <- set_equilibrium(parameters = simparams, init_EIR = starting_EIR)

static_booster_coverage <- set_pev_epi(
  simparams,
  profile = rtss_profile, # We will model implementation of the RTSS vaccine.
  timesteps = 1 * year, # Vaccination will begin at 1 year into the simulation.
  coverages = 1, # Vaccine coverage is 100%.
  min_wait = 0, # There is no minimum wait since the last vaccination.
  age = 5 * month, # Individuals will be vaccinated once they reach 5 months of age.
  booster_timestep = 12 * month, # The booster is administered 12 months following the third dose. 
  booster_coverage = 0.95, # 95% of those vaccinated with the primary series will be boosted.
  booster_profile = list(rtss_booster_profile) # We will model implementation of the RTSS booster.
)

output <- run_simulation(timesteps = sim_length * 2, parameters = static_booster_coverage)
plot_doses()

I attempted to change booster coverage over time using multiple set_pev_epi() calls, but this updates booster coverage to a flat value for the time period:

# try changing booster coverage at a specified timestep ------------------------
changed_coverage <- set_pev_epi(
  static_booster_coverage,
  profile = rtss_profile, 
  timesteps = 1 * year,  
  coverages = 1,         
  min_wait = 0,          
  age = 5 * month,        
  booster_timestep =  48 * month, #change booster coverage to 20%, 4 years after initial series
  booster_coverage = 0.20, 
  booster_profile = list(rtss_booster_profile)
)


output <- run_simulation(timesteps = sim_length * 2, parameters = changed_coverage)
plot_doses()

If I specify two different coverage values in the same set_pev_epi() call, this manifests as multiple boosters:

# multiple boosters  -----------------------------------------------------------
rtssepiparams2 <- set_pev_epi(
  simparams,
  profile = rtss_profile, 
  timesteps = 1 * year, 
  coverages = 1, 
  age = 5 * month, 
  min_wait = 0, 
  booster_timestep = c(12 * month, 24 * month), # Here, we are testing a strategy with 2 boosters, one at 1 year after the 3rd dose and the second 2 years after the 3rd dose.
  booster_coverage = c(1, 1), # For each of the two boosters, coverage is 100%.
  booster_profile = list(rtss_booster_profile, rtss_booster_profile) 
)

output <- run_simulation(timesteps = sim_length * 2, parameters = rtssepiparams2)
plot_doses_multiple_boosters()

Given that vaccine coverage can vary over time for the initial three doses, I imagine it would not be a huge lift to implement similar for the booster dose. I think it would involve some change to these functions among others:

malariasimulation/R/pev.R

Lines 154 to 179 in 397a7b7

	create_pev_booster_listener <- function(
	variables,
	coverage,
	booster_number,
	pev_profile_index,
	next_booster_event,
	next_booster_delay,
	renderer,
	strategy
	) {
	render_name <- paste0("n_pev_", strategy, "_booster_", booster_number)
	renderer$set_default(render_name, 0)
	force(next_booster_event) # because R lazy evaluation is rubbish
	force(next_booster_delay)
	force(coverage)
	function(timestep, target) {
	target <- sample_bitset(target, coverage)
	variables$pev_timestep$queue_update(timestep, target)
	variables$pev_profile$queue_update(pev_profile_index, target)
	renderer$render(render_name, target$size(), timestep)
	if (!is.null(next_booster_event)) {
	next_booster_event$schedule(target, next_booster_delay)
	}
	}
	}

malariasimulation/R/pev.R

Lines 210 to 306 in 397a7b7

	attach_pev_dose_listeners <- function(
	variables,
	parameters,
	dose_events,
	booster_events,
	booster_delays,
	booster_coverages,
	pev_profile_indices,
	strategy,
	renderer
	) {
	# set up dosing
	for (d in seq_along(dose_events)) {
	dose_events[[d]]$add_listener(
	create_dosage_renderer(renderer, strategy, d)
	)
	if (d == length(dose_events)) {
	dose_events[[d]]$add_listener(
	create_pev_efficacy_listener(
	variables,
	pev_profile_indices[[1]]
	)
	)
	if (length(booster_events) > 0) {
	seasonal_boosters <- FALSE
	if (!is.null(parameters$pev_epi_seasonal_boosters)) {
	seasonal_boosters <- parameters$pev_epi_seasonal_boosters
	}
	if (seasonal_boosters) {
	dose_events[[d]]$add_listener(
	create_seasonal_booster_scheduler(
	booster_events[[1]],
	booster_delays[[1]],
	parameters
	)
	)
	} else {
	dose_events[[d]]$add_listener(
	individual::reschedule_listener(
	booster_events[[1]],
	booster_delays[[1]]
	)
	)
	}
	}
	}
	}
	# set up boosters
	for (b in seq_along(booster_events)) {
	if (b == length(booster_events)) {
	next_booster_event <- NULL
	next_booster_delay <- NULL
	} else {
	next_booster_event <- booster_events[[b + 1]]
	next_booster_delay <- diff(
	booster_delays[c(b, b + 1)]
	)
	}
	booster_events[[b]]$add_listener(
	create_pev_booster_listener(
	variables = variables,
	coverage = booster_coverages[[b]],
	booster_number = b,
	pev_profile_index = pev_profile_indices[[b + 1]],
	next_booster_event = next_booster_event,
	next_booster_delay = next_booster_delay,
	renderer = renderer,
	strategy = strategy
	)
	)
	}
	}
	create_seasonal_booster_scheduler <- function(
	booster_event,
	booster_delay,
	parameters
	) {
	function(timestep, target) {
	delay <- booster_delay - timestep %% 365
	if (delay < 0) {
	delay <- delay + 365
	}
	if (delay <= parameters$pev_epi_min_wait) {
	delay <- delay + 365
	}
	booster_event$schedule(target, delay)
	}
	}
	sample_pev_param <- function(profile_index, profile_list, param_name) {
	mu <- vnapply(profile_list, function(p) p[[param_name]][[1]])
	sigma <- vnapply(profile_list, function(p) p[[param_name]][[2]])
	rnorm(length(profile_index), mu[profile_index], sigma[profile_index])
	}

We can chat about this more on Thursday @giovannic !