In cmhoove14/DDNTD: Density Dependence in Neglected Tropical Diseases

knitr::opts_chunk$set(echo = TRUE)

devtools::load_all()

require(tidyverse)
require(viridis)
require(RColorBrewer)

W_eq <- 60
Ip_eq <- 0.1

fit_pars <- base_pars
fit_alpha_beta_Neq <- fit_pars_from_eq_vals(beta_guess = 2e-4, alpha_guess = 2e-4, Neq_guess = 1000, 
                                            W = W_eq, Ip = Ip_eq, pars = base_pars)

fit_pars["beta"] <- fit_alpha_beta_Neq[1]
fit_pars["alpha"] <- fit_alpha_beta_Neq[2]
N_eq = fit_alpha_beta_Neq[3]
I_eq = N_eq*Ip_eq

Reff_W(W = W_eq, pars = fit_pars)  
R_eff_W_I(W = W_eq, I = I_eq, pars = fit_pars)  

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

Reff_R0 <- data.frame(W = W_seq) %>% 
  mutate(dWdt = sapply(W, dWdt_W, pars = fit_pars),
         Reff = sapply(W, R_eff_W_I, pars = fit_pars, I = I_eq,
                       PDD = phi_Wk, DDF = rho_Wk, DDI = nil_1),
         R0 = get_R0(fit_pars),
         Reff_W = sapply(W, Reff_W, pars = fit_pars)[1,],
         I_P = sapply(W, Reff_W, pars = fit_pars)[3,],
         Lambda = sapply(W, Reff_W, pars = fit_pars)[4,])

reff_crv <- Reff_R0 %>% 
  ggplot(aes(x = W, y = Reff_W)) +
    geom_line(size = 1.2) +
    geom_hline(yintercept = 1, lty = 2) +
    theme_classic() +
    theme(axis.text = element_text(size = 12),
          axis.title = element_text(size = 14))  +
    #ylim(c(0,2)) +
    scale_x_continuous(trans = "log",
                       breaks = c(1e-4, 1e-3, 1e-2, 0.1, 1, 10, 100),
                       labels = c("0.0001", "0.001", "0.01", "0.1", "1", "10", "100"),
                       limits = c(1e-4, 200)) +
    labs(x = "Mean Worm Burden (W)",
         y = expression(R[eff]),
         title = expression(The~Effective~Reproductive~Rate))

reff_crv

Reff_R0 %>% 
  gather("R", "Value", Reff:Reff_W) %>% 
  ggplot(aes(x = W, y = Value, col = R)) +
    geom_line(size = 1.2) +
    geom_hline(yintercept = 1, lty = 2) +
    theme_classic() +
    #ylim(c(0,10)) +
        scale_x_continuous(trans = "log",
                       breaks = c(1e-4, 1e-3, 1e-2, 0.1, 1, 10, 100),
                       labels = c("0.0001", "0.001", "0.01", "0.1", "1", "10", "100"),
                       limits = c(1e-4, 200)) +
    labs(x = "Mean Worm Burden (W)",
         y = "Reproduction Number",
         title = expression(R[eff]~R[0]~Comparison))

If you don't allow the infected snail population to change as W changes, Reff -> big numbers as W is reduced

Simulate treated/untreated worm burden compartment model dynamics through time

#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 = fit_pars)[["y"]]

#simulate annual MDA 
eff = 0.94 #94% efficacy
fit_pars["cvrg"] <- 0.6 # 60 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 = fit_pars,
                                    events_df = mda.events) %>% 
  mutate(W = Wt*fit_pars["cvrg"] + Wu*(1-fit_pars["cvrg"]),
         N = S+E+I,
         Reff = sapply(W, Reff_W,
                       pars = fit_pars)[1,],
         Reff_Wt = sapply(Wt, Reff_W,
                          pars = fit_pars)[1,],
         Reff_Wu = sapply(Wu, Reff_W,
                          pars = fit_pars)[1,],
         Reff2 = Reff_Wt * fit_pars["cvrg"] + Reff_Wu * (1-fit_pars["cvrg"]),
         reff_w_i = map2_dbl(W, I, R_eff_W_I, pars = fit_pars),
         Coverage = fit_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 ", fit_pars["cvrg"]*100, " % coverage"))

schisto_base_plot

schisto_snail_plot <- schisto_base_sim %>%  
  dplyr::select(time, S, E, I, N) %>% 
  gather("Infection_class", "Pop_Size", S:N) %>% 
    ggplot(aes(x = time, y = Pop_Size/fit_pars["K"], col = Infection_class)) +
    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(y = "Pop size scaled to carrying capacity, K",
           x = "time (years)") +
      ggtitle("Snail infection dynamics", 
              subtitle = paste0("Anual MDA for ", years/2, " years at ", fit_pars["cvrg"]*100, " % coverage"))

schisto_snail_plot

schisto_base_sim %>% 
  ggplot(aes(x = time, y = reff_w_i)) +
    geom_line(size = 1.2) 

Simulate single worm burden compartment model through time

With annual MDA

#Base time and starting values for state variables
mod_start <- c(S=5000, E=0, I=0, W=60)

#Run to equibrium with base parameter set
mod_eqbm <- runsteady(y = mod_start, func = DDNTD::schisto_mod,
                       parms = fit_pars)[["y"]]

mda.events2 = data.frame(var = rep('W', years/2),
                        time = c(1:(years/2))*365,
                        value = rep((1 - eff)*fit_pars["cvrg"]+(1-fit_pars["cvrg"]), years/2),
                        method = rep('mult', years/2))

schisto_mod_sim <- sim_schisto_mod(nstart = mod_eqbm, 
                                    time = base_time, 
                                    model = schisto_mod,
                                    pars = fit_pars,
                                    events_df = mda.events2) %>% 
  mutate(N = S+E+I,
         Reff = sapply(W, Reff_W,
                       pars = fit_pars)[1,],
         reff_w_i = map2_dbl(W, I, R_eff_W_I, pars = fit_pars),
         Coverage = fit_pars["cvrg"], 
         Intervention = "Annual MDA")

schisto_mod_plot <- schisto_mod_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_x_continuous(breaks = c(0:years)*365,
                         labels = c(-1:(years-1))) +
      theme_classic() +
      ggtitle("Human infection dynamics", 
              subtitle = paste0("Annual MDA for ", years/2, " years at ", fit_pars["cvrg"]*100, " % coverage"))

schisto_mod_plot

schisto_mod_sim %>% 
  ggplot(aes(x = time, y = reff_w_i)) +
    geom_line(size = 1.2) +
    theme_classic() +
    geom_hline(yintercept = 1, lty = 2) +
    scale_x_continuous(breaks = c(0:years)*365,
                       labels = c(-1:(years-1)))

With annual snail control

snail_eff = 0.85

snail.events = data.frame(var = rep(c('S', "E", "I"), years/2),
                          time = rep(c(1:(years/2))*365, each = 3),
                          value = rep(1-snail_eff, (years/2)*3),
                          method = rep('mult',( years/2)*3))

schisto_snail_sim <- sim_schisto_mod(nstart = mod_eqbm, 
                                     time = base_time, 
                                     model = schisto_mod,
                                     pars = fit_pars,
                                     events_df = snail.events) %>% 
  mutate(N = S+E+I,
         Reff = sapply(W, Reff_W,
                       pars = fit_pars)[1,],
         reff_w_i = map2_dbl(W, I, R_eff_W_I, pars = fit_pars),
         Coverage = fit_pars["cvrg"], 
         Intervention = "Annual Snail Control")

schisto_snail_plot <- schisto_snail_sim %>% 
    ggplot(aes(x = time, y = W)) +
    annotate("rect", xmin = 365, xmax = max(snail.events$time), ymin = -Inf, ymax = Inf,
             alpha = .2) +
      geom_line(size = 1.2, col = "purple") +
      scale_x_continuous(breaks = c(0:years)*365,
                         labels = c(-1:(years-1))) +
      theme_classic() +
      ylim(c(0,62)) +
      ggtitle("Human infection dynamics", 
              subtitle = paste0("Anual Snail Control for ", years/2, " years at ", snail_eff*100, " % efficiency"))

schisto_snail_plot

schisto_snail_sim %>% 
  ggplot(aes(x = time, y = reff_w_i)) +
    geom_line(size = 1.2) +
    theme_classic() +
    geom_hline(yintercept = 1, lty = 2) +
    scale_x_continuous(breaks = c(0:years)*365,
                       labels = c(-1:(years-1)))

Combined snail and MDA intervention

mda.snail.events <- rbind(mda.events2, snail.events) %>% 
  arrange(time)

schisto_comb_sim <- sim_schisto_mod(nstart = mod_eqbm, 
                                    time = base_time, 
                                    model = schisto_mod,
                                    pars = fit_pars,
                                    events_df = mda.snail.events) %>% 
  mutate(N = S+E+I,
         Reff = sapply(W, Reff_W,
                       pars = fit_pars)[1,],
         reff_w_i = map2_dbl(W, I, R_eff_W_I, pars = fit_pars),
         Coverage = fit_pars["cvrg"], 
         Intervention = "Annual MDA & Snail Control")

schisto_comb_plot <- schisto_comb_sim %>% 
    ggplot(aes(x = time, y = W)) +
    annotate("rect", xmin = 365, xmax = max(mda.snail.events$time), ymin = -Inf, ymax = Inf,
             alpha = .2) +
      geom_line(size = 1.2, col = "purple") +
      scale_x_continuous(breaks = c(0:years)*365,
                         labels = c(-1:(years-1))) +
      theme_classic() +
      ylim(c(0,62)) +
      ggtitle("Human infection dynamics", 
              subtitle = paste0("Anual Snail Control for ", years/2, " years at ", snail_eff*100, " % efficiency", 
                                "/nand Annual Snail Control for ", years/2, " years at ", snail_eff*100, " % efficiency"))

schisto_comb_plot

schisto_comb_sim %>% 
  ggplot(aes(x = time, y = reff_w_i)) +
    geom_line(size = 1.2) +
    theme_classic() +
    geom_hline(yintercept = 1, lty = 2) +
    scale_x_continuous(breaks = c(0:years)*365,
                       labels = c(-1:(years-1)))

Check out $R_{eff}$ values over a range of I and W

W_I_grid <- expand.grid(W = exp_seq(0.001, 200, 200),
                        I = exp_seq(0.001, 200, 200))

#Get reff across all values of W and I
w_i_reff <- mapply(R_eff_W_I, W_I_grid$W, W_I_grid$I, MoreArgs = list(pars = fit_pars))

W_I_Reff_df <- as.data.frame(cbind(W_I_grid, w_i_reff)) %>% 
  mutate(W_I_ratio = W/I)

reff1_vals <- W_I_Reff_df %>% 
  filter(round(w_i_reff,2) == 1.00)

#Get reff just from W (assume fast snail infection dynamics)
w_reff <- as.data.frame(t(sapply(W_I_grid$W, Reff_W, pars = fit_pars))) %>% 
  mutate(W = W_I_grid$W,
         I = N*I_P)

Reff_W_I_heat <- W_I_Reff_df %>% 
  ggplot(aes(x = W, y = I)) +
    geom_tile(aes(fill = log(w_i_reff))) +
    theme_classic() +
    theme(axis.title = element_text(size = 14),
          axis.text = element_text(size = 12)) +
    scale_fill_distiller(type = "div", palette = "RdBu") +
    scale_x_continuous(trans = "log",
                       breaks = c(1e-3, 1e-2, 0.1, 1, 10, 100),
                       labels = c("0.001", "0.01", "0.1", "1", "10", "100"),
                       limits = c(0.01, 200))  +
    scale_y_continuous(trans = "log",
                       breaks = c(1e-3, 1e-2, 0.1, 1, 10, 100),
                       labels = c("0.001", "0.01", "0.1", "1", "10", "100"),
                       limits = c(0.01, 200)) +
    geom_line(data = w_reff, aes(x = W, y = I), col = "white", size = 1.2) +
    geom_point(data = data.frame(W_endemic = W_eq,
                                 I_endemic = I_eq),
               aes(x = W_endemic,
                   y = I_endemic),
               col = "red", size = 3) +
    geom_point(data = w_reff[which(round(w_reff$Reff,2) == 1.00),][1,],
               aes(x = W,
                   y = I),
               col = "blue", size = 3) +
    geom_line(data = reff1_vals, aes(x = W, y = I), col = "black", lty = 2) +
    labs(x = expression(Mean~Worm~Burden~(W)),
         y = expression(Infected~snail~population~size~(I)),
         fill = expression(paste("log(", R[eff], "(W,I))")))

Reff_W_I_heat

Add intervention simulations through time to plot

int_df <- rbind(schisto_mod_sim, schisto_comb_sim)

Reff_W_I_heat +
  geom_line(data = int_df,# %>% filter(time <= (years/2)*365), 
            aes(x = W, y = I, col = time, group = Intervention), size = 1.2) +
    scale_color_viridis(option = "viridis", direction = -1)

Different version showing the Reff curve with the MDA intervention through time

reff_crv + ylim(c(0,5)) +
  geom_line(data = schisto_mod_sim,
            aes(x = W, y = reff_w_i, col = time),
            size = 1.2) +
    scale_color_distiller(type = "div", palette = "RdBu")

Worm burden breakpoint population size across key parameters

low_alpha <- fit_pars_from_eq_vals(2e-4, 2e-4, fit_pars["K"], 10, 0.1, fit_pars)[2]

high_alpha <- fit_pars_from_eq_vals(2e-4, 2e-4, fit_pars["K"], 100, 0.1, fit_pars)[2]

pars_grid <- expand.grid(K = seq(500, 1500, length.out = 5),
                         mu_W = seq(1/(2*365), 1/(5*365), length.out = 8),
                         omega = seq(0.01, 0.1, length.out = 5),
                         alpha = seq(low_alpha, high_alpha, length.out = 20))

#Determine W_bp across parameter and W ranges
  pars_grid$W_bp <- apply(pars_grid, 1, 
                           FUN = function(x){use_pars <- fit_pars ; 
                           use_pars["K"] <- x[["K"]] ;
                           use_pars["mu_W"] <- x[["mu_W"]] ;
                           use_pars["omega"] <- x[["omega"]] ;
                           use_pars["alpha"] <- x[["alpha"]] ;
                           get_breakpoint_W_N(0.1, use_pars["K"]*0.7, use_pars)[1]})

  pars_grid %>% 
    ggplot(aes(x = alpha, y = mu_W)) +
      geom_tile(aes(fill = W_bp)) +
      theme_classic() +
      facet_grid(K~omega) +
      scale_fill_viridis(option = "magma") +
      scale_y_continuous(name = expression(Mean~Adult~Worm~Lifespan),
                         breaks = c(1/(2*365),1/(3*365),1/(4*365),1/(5*365)),
                         labels = c(2:5))

W_get_N <- function(W, pars){
    ##standard snail parameters
    r=pars["r"]             # recruitment rate (from sokolow et al)
    K=pars["K"]          # carrying capacity corresponding to 50 snails per square meter
    mu_N=pars["mu_N"]          # Mean mortality rate of snails (assuming lifespan of 60days)
    sigma=pars["sigma"]         # Transition rate from exposed to infected (assuming pre-patency period of ~4 weeks) doi:10.4269/ajtmh.16-0614
    mu_I=pars["mu_I"]          # Increased mortality rate of infected snails
    theta=pars["theta"]          # mean cercarial shedding rate per adult snail doi:10.4269/ajtmh.16-0614

  #Adult Worm, Miracidia and Cercariae Parameters
    mu_W = pars["mu_W"]   # death rate of adult worms
    mu_H = pars["mu_H"] # death rate of adult humans
    m = pars["m"]             # mean eggs shed per female worm per 10mL urine (truscott et al)
    v = pars["v"]           # mean egg viability (miracidia per egg)

  #Density dependence parameters
    zeta = pars["zeta"]       # parameter of fecundity reduction function
    xi = pars["xi"]        # parameter for acquired immunity function http://doi.wiley.com/10.1111/j.1365-3024.1992.tb00029.x

  #Human parameters
    H = pars["H"]
    U = pars["U"]          # mL urine produced per child per day /10mL https://doi.org/10.1186/s13071-016-1681-4
    omega = pars["omega"]          #  infection risk/contamination of SAC  (related to sanitation/education/water contact) 10.1186/s13071-016-1681-4

  #Transmission parameters
    alpha=pars["alpha"]       # Cercarial infection probability
    beta=pars["beta"]       # Man to snail trnamission probability for linear FOI

  # Get miracidial density as function of worm burden
    #Update clumping parameter
      kap = k_w_fx(W)

    #Estimate mating probability
      phi = phi_Wk(W = W, k = kap)

    #Estimate density dependent fecundity
      rho = rho_Wk(W, zeta, k = kap)

    # Estimate total miracidia entering snail habitat
      M_tot = 0.5*H*omega*v*m*W*phi*rho*U

  # Get man-to-snail FOI and snail population sizeas solution given M_tot and other parameters
      #Estimate N
        N <- uniroot.all(function(N) K*(1-(mu_N+(beta*M_tot/N))/(r*(1+(beta*M_tot/(N*(mu_N+sigma))))))-N,
                         interval = c(0,K))

    I_P <- (N*sigma)/((N*mu_I*(mu_N+sigma)/beta*M_tot) + mu_I + sigma)    

    return(c(N, I_P, M_tot) )   
}

N_from_W <- data.frame(W = W_seq) %>% 
  mutate(N_W = sapply(W, W_get_N, pars = fit_pars)[1,],
         I_P = sapply(W, W_get_N, pars = fit_pars)[2,],
         M = sapply(W, W_get_N, pars = fit_pars)[3,])

N_from_W %>% 
  ggplot(aes(x = W, y = N_W)) + geom_line() + theme_classic()

N_from_W %>% 
  ggplot(aes(x = W, y = I_P)) + geom_line() + theme_classic() +
          scale_x_continuous(trans = "log",
                       breaks = c(1e-4, 1e-3, 1e-2, 0.1, 1, 10, 100),
                       labels = c("0.0001", "0.001", "0.01", "0.1", "1", "10", "100"),
                       limits = c(1e-4, 200)) 

N_from_W %>% 
  ggplot(aes(x = W, y = M)) + geom_line() + theme_classic() +
    ylim(c(0,10)) +
          scale_x_continuous(trans = "log",
                       breaks = c(1e-4, 1e-3, 1e-2, 0.1, 1, 10, 100),
                       labels = c("0.0001", "0.001", "0.01", "0.1", "1", "10", "100"),
                       limits = c(1e-4, 200)) 

main_plot <- Reff_R0 %>% 
  ggplot(aes(x = W, y = dWdt*365)) +
    geom_line(size = 1.2) +
    geom_hline(yintercept = 0, lty = 2) +
    theme_classic() +
    theme(axis.text = element_text(size = 12),
          axis.title = element_text(size = 14),
          axis.title.y = element_text(angle = 0, vjust = 0.5)) +
    scale_x_continuous(trans = "log",
                       breaks = c(1e-4, 1e-3, 1e-2, 0.1, 1, 10, 100),
                       labels = c("0.0001", "0.001", "0.01", "0.1", "1", "10", "100"),
                       limits = c(1e-4, 200)) +
    labs(x = "Mean Worm Burden (W)",
         y = expression(italic(frac(dW, dt))))

main_plot

inset <- Reff_R0 %>% 
  ggplot(aes(x = W, y = dWdt*365)) +
    geom_line(size = 1.2) +
    geom_hline(yintercept = 0, lty = 2) +
    theme_classic() +
    theme(axis.text = element_text(size = 12),
          axis.title = element_blank(),
          axis.title.y = element_blank()) +
    ylim(c(-1e-3,1e-3)) +
    scale_x_continuous(trans = "log",
                       breaks = c(1e-4, 1e-3, 1e-2, 0.1, 1),
                       labels = c("0.0001", "0.001", "0.01", "0.1", "1"),
                       limits = c(1e-4, 1))

inset

dd_product <- data.frame(W = W_seq) %>% 
  mutate(phi_rho_gam = sapply(W, 
                              function(W) phi_Wk(W, k_w_fx(W)) * rho_Wk(W, age_strat_pars["zeta"], k_w_fx(W)) * gam_Wxi(W, age_strat_pars["xi"])),
         phi_rho_z10_gam = sapply(W, 
                                  function(W) phi_Wk(W, k_w_fx(W)) * rho_Wk(W, age_strat_pars["zeta"]*10, k_w_fx(W)) * gam_Wxi(W, age_strat_pars["xi"])),
         M = sapply(W, 
                    function(W) phi_Wk(W, k_w_fx(W)) * rho_Wk(W, age_strat_pars["zeta"], k_w_fx(W)) * W*fit_pars["H"]*fit_pars["m"]*fit_pars["v"]*fit_pars["omega"]*fit_pars["U"]*0.5))

  dd_product %>% 
    ggplot(aes(x = W, y = phi_rho_gam)) + 
      geom_line() +
      geom_line(aes(x = W, y = phi_rho_z10_gam), col = 2) +
      theme_classic() +
      labs(y = expression(Phi~rho~gamma))

  dd_product %>% 
    ggplot(aes(x = W, y = M)) + 
      geom_line() +
      theme_classic() +
      labs(y = "M") +
    ylim(c(0,10)) +
      scale_x_continuous(trans = "log",
                         breaks = c(1e-4, 1e-3, 1e-2, 0.1, 1),
                         labels = c("0.0001", "0.001", "0.01", "0.1", "1"),
                         limits = c(1e-4, 1))

Reff_IP <- data.frame(W = W_seq) %>% 
  mutate(Reff_I01 = sapply(W, R_eff_W_I_P, pars = fit_pars, I_P = 0.1,
                       PDD = phi_Wk, DDF = rho_Wk, DDI = nil_1)[1,],
         Reff_I005 = sapply(W, R_eff_W_I_P, pars = fit_pars, I_P = 0.05,
                       PDD = phi_Wk, DDF = rho_Wk, DDI = nil_1)[1,],
         Reff_I001 = sapply(W, R_eff_W_I_P, pars = fit_pars, I_P = 0.01,
                       PDD = phi_Wk, DDF = rho_Wk, DDI = nil_1)[1,])

Reff_IP %>% 
  gather("I Prev", "Reff", Reff_I01:Reff_I001) %>% 
  ggplot(aes(x = W, y = Reff, col = `I Prev`)) +
    geom_line(size = 1.2) +
    geom_hline(yintercept = 1, lty = 2) +
    theme_classic() +
        scale_x_continuous(trans = "log",
                       breaks = c(1e-4, 1e-3, 1e-2, 0.1, 1, 10, 100),
                       labels = c("0.0001", "0.001", "0.01", "0.1", "1", "10", "100"),
                       limits = c(1e-4, 200)) +
    labs(x = "Mean Worm Burden (W)",
         y = expression(R[eff]),
         title = expression(R[eff]~I[P]~Relationship))

W_bp_IP <- data.frame(I_P = seq(0.001,0.1,0.001)) %>% 
  mutate(Wbp = sapply(I_P, W_bp_I_P, pars = fit_pars,
                      PDD = phi_Wk, DDF = rho_Wk, DDI = nil_1))

W_bp_IP %>% 
  mutate(Wbp = if_else(Wbp > 100, NA_real_, Wbp)) %>% 
  ggplot(aes(x = I_P, y = Wbp)) +
    geom_line(size = 1.2) +
    theme_classic() +
    theme(axis.text = element_text(size = 12),
          axis.title = element_text(size = 14))  +
    scale_y_continuous(trans = "log",
                       breaks = c(1e-4, 1e-3, 1e-2, 0.1, 1, 10),
                       labels = c("0.0001", "0.001", "0.01", "0.1", "1", "10"),
                       limits = c(1e-4, 10)) +
    scale_x_continuous(breaks = c(0,0.025,0.05,0.075,0.1),
                       labels = c("0", "2.5", "5.0", "7.5", "10")) +    
    labs(y = expression(Worm~Burden~Breakpoint~(W[bp])),
         x = expression(paste("Infected Snail Prevalence (", I[P], ", %)")),
         title = "Breakpoint as function of snail prevalence")

W_bp_R0 <- data.frame(R0 = seq(0.95,5,0.05)) %>% 
  mutate(Wbp = sapply(R0, W_bp_R_0, 
                      PDD = phi_Wk, DDF = rho_Wk, DDI = nil_1,
                      zeta = age_strat_pars["zeta"], xi = age_strat_pars["xi"]))

W_bp_R0 %>% 
  ggplot(aes(x = R0, y = Wbp)) +
    geom_line(size = 1.2) +
    theme_classic() +
    theme(axis.text = element_text(size = 12),
          axis.title = element_text(size = 14))  +
    scale_x_continuous(breaks = c(1:5),
                       limits = c(1,5)) +    
    scale_y_continuous(breaks = seq(0.2, 1.6, by = 0.2),
                       limits = c(0.2,1.6)) +    
    labs(y = expression(italic(W[bp])),
         x = expression(italic(R[0])))

wt_wu <- expand.grid(W_t = exp_seq(0.01, 100, 100),
                     W_u = exp_seq(0.01, 100, 100))

tau <- mapply(FUN = function(wt, wu){(0.1-wu)/(wt-0.94*wt-wu)},
              wt_wu$W_t, wt_wu$W_u)*100

tau <- case_when(tau > 0 & tau < 100 ~ tau,
                 tau < 0 ~ NA_real_,
                 tau > 100 ~ NA_real_)

cvrg_sum <- cbind(wt_wu, tau)

as.data.frame(cvrg_sum) %>% 
  ggplot(aes(x = W_t, y = W_u, fill = tau, z = tau)) +
    scale_x_continuous(trans = "log",
                       breaks = c(0.1, 1,10,100),
                       labels = c("0.1", "1", "10", "100")) +
    scale_y_continuous(trans = "log",
                       breaks = c(0.1, 1,10,100),
                       labels = c("0.1", "1", "10", "100")) +
    geom_tile()


wt_tau <- sapply(exp_seq(0.01, 100, 100), function(Wt){(0.1-Wt)/(-0.94*Wt)}) * 100
wt_tau <- case_when(wt_tau > 0 & wt_tau < 100 ~ wt_tau,
                 wt_tau < 0 ~ NA_real_,
                 wt_tau > 100 ~ NA_real_)
plot(exp_seq(0.01, 100, 100), wt_tau, type = "l", xlim = c(0,5))

W_bp_R0 <- W_bp_R0 %>% 
  mutate(bp_cvrg = sapply(Wbp, cvrg_from_W_bp, W_t = W, epsilon = 0.94))

W_bp_R0 %>% 
  ggplot(aes(x = R0, y = bp_cvrg)) +
    geom_line(size = 1.2) +
    theme_classic() +
    theme(axis.text = element_text(size = 12),
          axis.title = element_text(size = 14))  +
    scale_x_continuous(breaks = c(1:4),
                       limits = c(1,4)) +    
    labs(y = expression(MDA~Coverage~to~reach~breakpoint),
         x = expression(R[0]),
         title = "Coverage to reach breakpoint as function of basic reproductive rate")

test_omegas <- seq(0.5,0.99,0.01)
result_Wbps <- numeric(length = length(test_omegas))

for(i in test_omegas){
  fit_pars -> use_pars ; use_pars["omega"] = fit_pars["omega"]*i

  result_Wbps[which(test_omegas == i)] <- W_bp(use_pars, PDD = phi_Wk, DDF=rho_Wk, DDI = nil_1)
}

W_bp_omega <- data_frame(omega_red = 100 - test_omegas*100,
                         W_bp = result_Wbps)

W_bp_omega %>% 
  ggplot(aes(x = omega_red, y = W_bp)) +
    geom_line(size = 1.2) +
    theme_classic() +
    theme(axis.text = element_text(size = 12),
          axis.title = element_text(size = 14))  +
    scale_x_continuous(breaks = c(0, 10, 20,30),
                       limits = c(0, 31)) +    
    labs(y = expression(italic(W[bp])),
         x = expression(paste("% reduction in exposure/contamination ", (italic(omega)))))

W_bp_omega <- W_bp_omega %>% 
  mutate(bp_cvrg = sapply(W_bp, cvrg_from_W_bp, W_t = W, epsilon = 0.94))

W_bp_omega %>% 
  ggplot(aes(x = omega_red, y = bp_cvrg*100)) +
    geom_line(size = 1.2) +
    theme_classic() +
    geom_hline(yintercept = 100, lty = 2) +
    theme(axis.text = element_text(size = 12),
          axis.title = element_text(size = 14))  +
    scale_x_continuous(breaks = c(0, 10, 20,30),
                       limits = c(0, 31)) +    
    labs(y = expression(MDA~Coverage~to~reach~W[bp]),
         x = expression(paste("% reduction in exposure/contamination ", (italic(omega))))) +
    annotate("rect", xmin = 0, xmax = 42, ymin = 100.01, ymax = Inf, 
             alpha = 0.2, col = "red", fill = "red")

W_bp_omega_Wt <- expand.grid(W_bp = W_bp_omega$W_bp,
                             W_t = exp_seq(1,100,30)) %>% 
  mutate(omega_red = rep(W_bp_omega$omega_red, times = 30),
         bp_cvrg = mapply(cvrg_from_W_bp, W_t, W_bp, MoreArgs = list(epsilon = 0.94))*100,
         bp_cvrg = case_when(bp_cvrg > 0 & bp_cvrg <100 ~ bp_cvrg,
                             bp_cvrg < 0 ~ NA_real_,
                             bp_cvrg > 100 ~ NA_real_))


W_bp_omega_Wt %>% 
  ggplot(aes(y = W_t, x = omega_red, z = bp_cvrg)) + 
    geom_tile(aes(fill = bp_cvrg)) + 
    scale_x_continuous(breaks = c(0, 10, 20,30),
                       limits = c(0, 31)) +
    scale_y_continuous(trans = "log",
                       breaks = c(1,10,100)) +
    scale_fill_viridis(option = "magma", limits = c(0,100), direction = -1) +
    theme_bw() +
    theme(axis.title = element_text(size = 14),
          axis.text = element_text(size = 12)) +
    labs(y = expression(paste("Pre-MDA mean worm burden, ", W[t])),
         x = expression(paste("% reduction in exposure/contamination ", (italic(omega)))),
         fill = expression(italic(cvrg))) #+ annotate("text", x = 15, y = 18, label = expression(paste(W[bp], " not achievable \nwith MDA")), col = "red", size = 6)

cvrg_df <- expand.grid(W_t = exp_seq(1e-3,100,10),
                       W_u = exp_seq(1e-3,100,10),
                       K = seq(fit_pars["K"], fit_pars["K"]*0.7, length.out = 3),
                       omega = seq(fit_pars["omega"], fit_pars["omega"]*0.7, length.out = 3),
                       epsilon = 0.94) %>% 
  mutate(cvrg = mapply(cvrg_from_pars, W_t, W_u, epsilon, K, omega, MoreArgs = list(pars = fit_pars))*100)

cvrg_df_lt0 <- cvrg_df %>% 
  filter(cvrg < 0 | is.na(cvrg))

cvrg_df_gt100 <- cvrg_df %>% 
  filter(cvrg > 100)

cvrg_df %>% 
  mutate(cvrg = case_when(cvrg > 0 & cvrg < 100 ~ cvrg,
                          cvrg < 0 ~ NA_real_,
                          cvrg > 100 ~ NA_real_)) %>% 
  ggplot(aes(y = W_t, x = W_u, z = cvrg)) + 
    geom_tile(aes(fill = cvrg)) + 
    scale_x_continuous(trans = "log",
                       breaks = c(1e-3, 1,100),
                       labels = c("0.001", "1",  "100")) +
    scale_y_continuous(trans = "log",
                       breaks = c(1e-3, 1e-2, 0.1, 1,10,100),
                       labels = c("0.001", "0.01", "0.1", "1", "10", "100")) +
    scale_fill_viridis(limits = c(0,100), direction = -1) +
#    scale_fill_gradientn(colours = c("black", "blue", "red", "white")) +
    geom_tile(data = cvrg_df_lt0, aes(y = W_t, x = W_u, fill = cvrg), fill = "white") +
    geom_tile(data = cvrg_df_gt100, aes(y = W_t, x = W_u, fill = cvrg), fill = "black") +
    theme_bw() +
    facet_grid(K ~ omega, labeller = label_bquote("K"==.(K/fit_pars["K"]),
                                                  omega==.(omega/fit_pars["omega"]))) +
    theme(axis.title = element_text(size = 14),
          axis.text = element_text(size = 12), 
          strip.background = element_rect(fill = "white"),
          strip.text = element_text(size = 12)) +
    labs(y = expression(paste("Pre-MDA mean worm burden, ", W[T])),
         x = expression(paste("Untreated mean worm burden, ", W[U])),
         fill = expression(italic(Tau))) #+ annotate("text", x = 15, y = 18, label = expression(paste(W[bp], " not achievable \nwith MDA")), col = "red", size = 6)

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