In cmhoove14/DDNTD: Density Dependence in Neglected Tropical Diseases

knitr::opts_chunk$set(echo = TRUE, warnings = FALSE)

require(deSolve)
require(rootSolve)
require(adaptivetau)
require(ggplot2)
require(tidyverse)
require(gganimate)
require(magick)

devtools::load_all("../../DDNTD")

Effects of MDA on $R_{eff}$

Basic schistosomiasis model simulation with MDA administered to treated segment of the population and untreated segment of the population left alone

#Base time and starting values for state variables
base_start <- c(S=5000, E=0, I=0, Wt=10, Wu=10)
years <- 20
base_time <- c(1:(365*years))

#Run to equibrium with base parameter set
base_eqbm <- runsteady(y = base_start, func = DDNTD::schisto_base_mod,
                       parms = DDNTD::base_pars)[["y"]]

#simulate annual MDA 
eff = 0.93 #93% efficacy
base_pars["cvrg"] <- 0.5 # 50 percent coverage

mda.events = data.frame(var = rep('Wt', years/2),
                        time = c(1:(years/2))*365,
                        value = rep((1 - eff), years/2),
                        method = rep('mult', years/2))

schisto_base_sim <- sim_schisto_mod(nstart = base_eqbm, 
                                    time = base_time, 
                                    model = schisto_base_mod,
                                    pars = base_pars,
                                    events_df = mda.events) %>% 
  mutate(W = Wt*base_pars["cvrg"] + Wu*(1-base_pars["cvrg"]),
         kap = sapply(W, k_from_log_W),
         kap_Wt = sapply(Wt, k_from_log_W),
         kap_Wu = sapply(Wu, k_from_log_W),
         Reff = sapply(W, Reff_W,
                       pars = base_pars)[1,],
         Reff_Wt = sapply(Wt, Reff_W,
                          pars = base_pars)[1,],
         Reff_Wu = sapply(Wu, Reff_W,
                          pars = base_pars)[1,],
         Reff2 = Reff_Wt * base_pars["cvrg"] + Reff_Wu * (1-base_pars["cvrg"]),
         Coverage = base_pars["cvrg"])

schisto_base_plot <- schisto_base_sim %>% 
  gather("Treatment_Group", "Worm_Burden", Wt:W) %>% 
    ggplot(aes(x = time, y = Worm_Burden, lty = Treatment_Group)) +
    annotate("rect", xmin = 365, xmax = max(mda.events$time), ymin = -Inf, ymax = Inf,
             alpha = .2) +
      geom_line(size = 1.2, col = "purple") +
      scale_linetype_manual(values = c("W" = 1,
                                       "Wt" = 2,
                                       "Wu" = 3),
                            labels = c("Mean", "Treated", "Untreated")) +
      scale_x_continuous(breaks = c(0:years)*365,
                         labels = c(-1:(years-1))) +
      theme_classic() +
      ggtitle("Human infection dynamics", 
              subtitle = paste0("Anual MDA for ", years/2, " years at ", base_pars["cvrg"]*100, " % coverage"))

schisto_base_plot

Look at influence of increasing MDA coverage on dynamics

base_pars -> base_pars_cvrg99 ; base_pars_cvrg99["cvrg"] = .99
base_pars -> base_pars_cvrg90 ; base_pars_cvrg90["cvrg"] = .90
base_pars -> base_pars_cvrg80 ; base_pars_cvrg80["cvrg"] = .80
base_pars -> base_pars_cvrg70 ; base_pars_cvrg70["cvrg"] = .70
base_pars -> base_pars_cvrg60 ; base_pars_cvrg60["cvrg"] = .60

schisto_cvrg_fx <- function(pars){
  schisto_sim <- sim_schisto_mod(nstart = base_eqbm, 
                                      time = base_time, 
                                      model = schisto_base_mod,
                                      pars = pars,
                                      events_df = mda.events) %>% 
  mutate(W = Wt*pars["cvrg"] + Wu*(1-pars["cvrg"]),
         kap = sapply(W, k_from_log_W),
         kap_Wt = sapply(Wt, k_from_log_W),
         kap_Wu = sapply(Wu, k_from_log_W),
         Reff = sapply(W, Reff_W,
                       pars = base_pars),
         Reff_Wt = sapply(Wt, Reff_W,
                          pars = base_pars),
         Reff_Wu = sapply(Wu, Reff_W,
                          pars = base_pars),
         Reff2 = Reff_Wt * pars["cvrg"] + Reff_Wu * (1-pars["cvrg"]),
         Coverage = pars["cvrg"])

  return(schisto_sim)
}

schisto_cvrg60_sim <- schisto_cvrg_fx(base_pars_cvrg60)
schisto_cvrg70_sim <- schisto_cvrg_fx(base_pars_cvrg70)
schisto_cvrg80_sim <- schisto_cvrg_fx(base_pars_cvrg80)
schisto_cvrg90_sim <- schisto_cvrg_fx(base_pars_cvrg90)
schisto_cvrg99_sim <- schisto_cvrg_fx(base_pars_cvrg99)

schisto_cvrgs_df <- rbind(schisto_base_sim,
                          schisto_cvrg60_sim, 
                          schisto_cvrg70_sim,
                          schisto_cvrg80_sim,
                          schisto_cvrg90_sim,
                          schisto_cvrg99_sim)

schisto_cvrgs_plot <- schisto_cvrgs_df %>% 
  ggplot(aes(x = time, y = W, col = as.factor(Coverage))) +
    annotate("rect", xmin = 365, xmax = max(mda.events$time), ymin = -Inf, ymax = Inf,
             alpha = .2) +
      geom_line(size = 1.2) +
      scale_x_continuous(breaks = c(0:years)*365,
                         labels = c(-1:(years-1))) +
      theme_classic() +
      labs(title = "Human infection dynamics", 
           subtitle = paste0("Annual MDA for ", years/2, " years at variable coverage"),
           x = "Time (years)",
           color = "Coverage")

schisto_cvrgs_plot

Animate worm burden simulations through time at different coverages

w_anim <- schisto_cvrgs_plot +
  transition_reveal(time)

w_anim

Baseline $R_{eff}$ curve with base parameter set

test_ws <- exp_seq(1e-4,200,200)

#Get values of mean worm burden to plot R effective curve over
Reffs <- data.frame(test_ws = test_ws) %>% 
  mutate(w_det_ks = sapply(test_ws, k_from_log_W),
         Reff = sapply(test_ws, Reff_W, pars = base_pars)[1,])

base_reff <- Reffs %>% 
  ggplot(aes(x = test_ws, y = Reff)) +
  geom_line(size = 1.2) +
  theme_classic() +
  scale_x_continuous(trans = "log",
                     breaks = c(0.001, 0.01, 0.1, 1, 10, 100),
                     labels = c("0.001", "0.01", "0.1", "1", "10", "100"),
                     limits = c(0.001, 200)) +
  scale_y_continuous(breaks = c(0:3),
                     limits = c(0,3)) +
  geom_hline(yintercept = 1, lty = 2) +
  labs(x = expression(Mean~Worm~Burden~(italic(W))),
       y = expression(italic(R[eff]))) 

base_reff

Map $R_{eff}$ estimates from model simulations through time onto $R_{eff}$ curve to observe changes through time as MDA pushes population towards the breakpoint

reff_anim <- base_reff +
  geom_point(data = schisto_cvrgs_df, aes(x = W, y = Reff, col = as.factor(Coverage)),
             size = 4) +
  theme(legend.position = "none") +
  transition_time(time)

reff_anim

Combine two animations in single image

w_gif <- image_read(animate(w_anim, width = 720, height = 360))
reff_gif <- image_read(animate(reff_anim, width = 720, height = 360))

comb_gif <- image_append(c(w_gif[1], reff_gif[1]), stack = TRUE)

for(i in 2:100){
  combined <- image_append(c(w_gif[i], reff_gif[i]), stack = TRUE)
  comb_gif <- c(comb_gif, combined)
}

comb_gif

Stochastic schisto model

Annual MDA at 80% coverage

set.seed(10)
sim_schisto_stoch <- function(sim){
  sim_schisto_stoch_mod(nstart = round(base_eqbm),
                        params = base_pars_cvrg80,
                        tf = max(base_time),
                        events_df = mda.events) %>% 
    mutate(W = Wt*base_pars_cvrg80["cvrg"] + Wu*(1-base_pars_cvrg80["cvrg"]),
           kap = map_dbl(W, k_from_log_W),
           Reff = pmap_dbl(list(parameters = list(base_pars_cvrg80),
                                W = W,
                                kap = kap), getReff),
           sim = sim)
}

schisto_stoch_sims <- bind_rows(map_df(c(1:10), sim_schisto_stoch)) %>% 
  mutate(W = if_else(is.nan(Reff), 0.005, W),
         Reff = if_else(is.nan(Reff), 
                        getReff(base_pars_cvrg80, 0.005, k_from_log_W(0.005)), 
                        Reff))

schisto_stoch_sims_plot <- schisto_stoch_sims %>% 
  ggplot(aes(x = time, y = W, col = as.factor(sim))) +
    annotate("rect", xmin = 365, xmax = max(mda.events$time), ymin = -Inf, ymax = Inf,
             alpha = .2) +
    geom_line(size = 1.1) +
    scale_x_continuous(breaks = c(0:years)*365,
                       labels = c(-1:(years-1))) +
    theme_classic() +
    theme(legend.position = "none") +
    labs(x = "time (years)",
         y = expression(mean~worm~burden~(italic(W))),
         title = "Human infection dynamics from stochastic model", 
         subtitle = paste0("Anual MDA for ", years/2, " years at ", base_pars_cvrg80["cvrg"]*100, " % coverage"))

schisto_stoch_sims_plot +
  transition_reveal(time)

base_reff +
  geom_point(data = schisto_stoch_sims, aes(x = W, y = Reff, col = as.factor(sim)),
             size = 4) +
  theme(legend.position = "none") +
  transition_time(time) +
  labs(title = "Effective reproduction number year {round(frame_time/365,2)}")

Effects of snail habitat reduction on $R_{eff}$

Now look into interventions other than MDA that "push the curve down" rather than acting on the state variable $W$

Example: Amount of snail habitat (modeled as snail environmental carrying capacity)

base_pars -> base_pars_C.90 ; base_pars_C.90["C"] = base_pars["C"]*.90
base_pars -> base_pars_C.80 ; base_pars_C.80["C"] = base_pars["C"]*.80
base_pars -> base_pars_C.70 ; base_pars_C.70["C"] = base_pars["C"]*.70
base_pars -> base_pars_C.60 ; base_pars_C.60["C"] = base_pars["C"]*.60
base_pars -> base_pars_C.50 ; base_pars_C.50["C"] = base_pars["C"]*.50

Reffs_C_red <- Reffs %>% 
  mutate(Reff_C.90 = pmap_dbl(list(parameters = list(base_pars_C.90),
                                   W = test_ws,
                                   kap = w_det_ks), getReff),
         Reff_C.80 = pmap_dbl(list(parameters = list(base_pars_C.80),
                                   W = test_ws,
                                   kap = w_det_ks), getReff),
         Reff_C.70 = pmap_dbl(list(parameters = list(base_pars_C.70),
                                   W = test_ws,
                                   kap = w_det_ks), getReff),
         Reff_C.60 = pmap_dbl(list(parameters = list(base_pars_C.60),
                                   W = test_ws,
                                   kap = w_det_ks), getReff),
         Reff_C.50 = pmap_dbl(list(parameters = list(base_pars_C.50),
                                   W = test_ws,
                                   kap = w_det_ks), getReff))

Reffs_C_red_plot <- Reffs_C_red %>% 
  gather("Snail Habitat", "Reff", Reff:Reff_C.50) %>% 
  ggplot(aes(x = test_ws, y = Reff, col = `Snail Habitat`)) +
    geom_line(size = 1.2) +
    theme_classic() +
    scale_x_continuous(trans = "log",
                       breaks = c(0.001, 0.01, 0.1, 1, 10, 100),
                       labels = c("0.001", "0.01", "0.1", "1", "10", "100"),
                       limits = c(0.001, 200)) +
    scale_y_continuous(breaks = c(0:3),
                       limits = c(0,3)) +
    geom_hline(yintercept = 1, lty = 2) +
    scale_color_manual(breaks = c("Reff", "Reff_C.90", "Reff_C.80", "Reff_C.70", "Reff_C.60", "Reff_C.50"),
                       labels = c("Base",
                                  "10% red",
                                  "20% red",
                                  "30% red",
                                  "40% red",
                                  "50% red"),
                       values = c("black", "grey80", "grey60", "grey40", "grey20", "grey5")) +
    labs(x = expression(Mean~Worm~Burden~(italic(W))),
         y = expression(italic(R[eff]))) 

Reffs_C_red_plot

What if we have an annual intervention that affects the snail habitat like vegetation removal?

#Generate forcing function that changes carrying capacity parameter over time
C_force_fx <- gen_force_fx(burn_in = 365, 
                           t_int = max(mda.events$time),
                           t_max = max(base_time), 
                           freq = 365, 
                           base_pars, "C", 0.75)

#Generate copies of base parameter set in dataframe and replace with C that changes through time
C_force_pars <- as.data.frame(as.list(base_pars)) %>% slice(rep(1:n(), each = max(base_time))) %>% 
  mutate(C = map_dbl(c(1:max(base_time)), C_force_fx))


#New schisto model object with C forcing function
schisto_base_mod_C_force_fx <- function(t, n, parameters) {

  f_N <- parameters["f_N"]
  C <- C_force_fx(t)
  mu_N <- parameters["mu_N"]
  sigma <- parameters["sigma"]
  mu_I <- parameters["mu_I"]
  mu_W <- parameters["mu_W"]
  H <- parameters["H"]
  mu_H <- parameters["mu_H"]
  lambda <- parameters["lambda"]
  beta <- parameters["beta"]
  cvrg <- parameters["cvrg"]
  gamma <- parameters["gamma"]
  xi <- parameters["xi"]

    S=n[1]
    E=n[2]
    I=n[3]
    Wt=n[4]
    Wu=n[5]

    N=S+E+I

    W=(cvrg*Wt) + ((1-cvrg)*Wu) #weighting treated and untreated populations

  #Clumping parameters based on worm burden
    k_Wt = k_from_log_W(Wt)
    k_Wu = k_from_log_W(Wu)

  #Miracidial estimate from treated and untreated populations assuming 1:1 sex ratio, mating probability, density dependence
    M = (0.5*Wt*H*cvrg)*phi_Wk(Wt, k_Wt)*f_Wgk(Wt, gamma, k_Wt) +
        (0.5*Wu*H*(1-cvrg))*phi_Wk(Wu, k_Wu)*f_Wgk(Wu, gamma, k_Wu)


    dSdt= f_N*(1-(N/C))*(S+E) - mu_N*S - beta*M*S #Susceptible snails

    dEdt= beta*M*S - (mu_N+sigma)*E #Exposed snails

    dIdt= sigma*E - (mu_N+mu_I)*I #Infected snails

    #worm burden in human
    dWtdt= (lambda*I*R_Wv(Wt, xi)) - ((mu_W+mu_H)*Wt)
    dWudt= (lambda*I*R_Wv(Wu, xi)) - ((mu_W+mu_H)*Wu)



    return(list(c(dSdt,dEdt,dIdt,dWtdt,dWudt)))
}

schisto_C_force_sim <- sim_schisto_mod(nstart = base_eqbm, 
                                       time = base_time, 
                                       model = schisto_base_mod_C_force_fx,
                                       parameters = base_pars,
                                       events_df = NA) %>% 
  mutate(W = Wt*base_pars["cvrg"] + Wu*(1-base_pars["cvrg"]),
         kap = map_dbl(W, k_from_log_W))

schisto_C_force_sim_Reffs <- numeric()

for(i in 1:max(base_time)){
  schisto_C_force_sim_Reffs[i] <- getReff(as.list(C_force_pars[i,]), schisto_C_force_sim$W[i], schisto_C_force_sim$kap[i])
}

schisto_C_force_sim <- schisto_C_force_sim %>% 
  mutate(Reff = schisto_C_force_sim_Reffs)

schisto_C_force_sim_plot <- schisto_C_force_sim %>% 
    ggplot(aes(x = time, y = W)) +
    annotate("rect", xmin = 365, xmax = max(mda.events$time), ymin = -Inf, ymax = Inf,
             alpha = .2) +
      geom_line(size = 1.2, col = "purple") +
      scale_linetype_manual(values = c("W" = 1,
                                       "Wt" = 2,
                                       "Wu" = 3),
                            labels = c("Mean", "Treated", "Untreated")) +
      scale_x_continuous(breaks = c(0:years)*365,
                         labels = c(-1:(years-1))) +
      scale_y_continuous(limits = c(0, base_eqbm["Wt"]+1)) +
      theme_classic() +
      ggtitle("Human infection dynamics", 
              subtitle = paste0("Annual reduction in snail habitat e.g. due to vegetation removal"))

schisto_C_force_sim_plot

Animate worm burden simulations through time with snail habitat reduction intervention

schisto_C_force_anim <- schisto_C_force_sim_plot +
  transition_reveal(time)

schisto_C_force_anim

Animate $R_{eff}$ over the course of the snail habitat intervention

reff_anim_C_force <- Reffs_C_red_plot +
  geom_point(data = schisto_C_force_sim, aes(x = W, y = Reff),
             size = 4, col = "red") +
  transition_time(time)

reff_anim_C_force

Combine two animations in single image

schisto_C_force_gif <- image_read(animate(schisto_C_force_anim, width = 720, height = 360))
reff_C_force_gif <- image_read(animate(reff_anim_C_force, width = 720, height = 360))

comb_C_force_gif <- image_append(c(schisto_C_force_gif[1], reff_C_force_gif[1]), stack = TRUE)

for(i in 2:100){
  combined <- image_append(c(schisto_C_force_gif[i], reff_C_force_gif[i]), stack = TRUE)
  comb_C_force_gif <- c(comb_C_force_gif, combined)
}

comb_C_force_gif

Effects of combined interventions on $R_{eff}$

MDA + Snail habitat control

schisto_cvrg_C_red_fx <- function(pars){
  schisto_sim <- sim_schisto_mod(nstart = base_eqbm, 
                                      time = base_time, 
                                      model = schisto_base_mod_C_force_fx,
                                      parameters = pars,
                                      events_df = mda.events) %>% 
  mutate(W = Wt*pars["cvrg"] + Wu*(1-pars["cvrg"]),
         kap = map_dbl(W, k_from_log_W))

  schisto_sim_Reffs <- numeric()

for(i in 1:max(base_time)){
  schisto_sim_Reffs[i] <- getReff(as.list(C_force_pars[i,]), schisto_sim$W[i], schisto_sim$kap[i])
}

schisto_sim <- schisto_sim %>% 
  mutate(Reff = schisto_sim_Reffs,
         Coverage = pars["cvrg"])

  return(schisto_sim)
}

schisto_cvrg60_C_red_sim <- schisto_cvrg_C_red_fx(base_pars_cvrg60)
schisto_cvrg70_C_red_sim <- schisto_cvrg_C_red_fx(base_pars_cvrg70)
schisto_cvrg80_C_red_sim <- schisto_cvrg_C_red_fx(base_pars_cvrg80)
schisto_cvrg90_C_red_sim <- schisto_cvrg_C_red_fx(base_pars_cvrg90)
schisto_cvrg99_C_red_sim <- schisto_cvrg_C_red_fx(base_pars_cvrg99)

schisto_cvrgs_C_red_df <- rbind(schisto_C_force_sim %>% mutate(Coverage = 0),
                                schisto_cvrg60_C_red_sim, 
                                schisto_cvrg70_C_red_sim,
                                schisto_cvrg80_C_red_sim,
                                schisto_cvrg90_C_red_sim,
                                schisto_cvrg99_C_red_sim)

schisto_cvrgs_C_red_plot <- schisto_cvrgs_C_red_df %>% 
  ggplot(aes(x = time, y = W, col = as.factor(Coverage))) +
    annotate("rect", xmin = 365, xmax = max(mda.events$time), ymin = -Inf, ymax = Inf,
             alpha = .2) +
      geom_line(size = 1.2) +
      scale_x_continuous(breaks = c(0:years)*365,
                         labels = c(-1:(years-1))) +
      theme_classic() +
      labs(title = "Human infection dynamics", 
           subtitle = paste0("Anual MDA for ", years/2, " years at variable coverage AND annual snail habitat reduction"),
           x = "Time (years)",
           color = "Coverage")

schisto_cvrgs_C_red_plot

reff_MDA_C_red_anim <- base_reff +
  geom_point(data = schisto_cvrgs_C_red_df, aes(x = W, y = Reff, col = as.factor(Coverage)),
             size = 4) +
  theme(legend.position = "none") +
  transition_time(time)

reff_MDA_C_red_anim

cmhoove14/DDNTD documentation built on Nov. 23, 2019, 7:04 p.m.