In cmhoove14/DDNTD: Density Dependence in Neglected Tropical Diseases

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

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

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

Background and questions

Relationship between prevalence and clumping parameter, $\kappa$

plot(mapply(est_prev_W_k, 
            (exp_seq(0.00001,1000, 200)), 
            sapply(exp_seq(0.00001,1000, 200), 
                   k_w_fx)), 
     sapply(exp_seq(0.00001,1000, 200), 
            k_w_fx), 
     type = "l", lwd = 2, xlab = "Prevalence", ylab = expression(kappa))

Worm aggregation and mating probability

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

phi_k0.1 <- sapply(test_ws, phi_Wk, k = 0.1)

as.data.frame(cbind(test_ws,
                    phi_k0.1)) %>% 
  ggplot(aes(x = test_ws, y = phi_k0.1)) +
    geom_line(size = 1.2) +
    theme_classic() +
    theme(title = element_text(size = 18),
          axis.title = element_text(size = 15),
          axis.text = element_text(size = 12)) +
    scale_x_continuous(trans = "log",
                       breaks = c(0.0001, 0.001, 0.01, 0.1, 1, 10, 100),
                       labels = c("0.0001", "0.001", "0.01", "0.1", "1", "10", "100"),
                       limits = c(0.01, 200)) +
    ylim(c(0,1)) +
    labs(x = expression(Mean~Worm~Burden~(italic(W))),
         y = expression(paste("Mating Probability, ", italic(Phi(W)))),
         title = "Mating probability across Mean Worm Burden") #+ theme(legend.position = "none")

Worm aggregation and fecundity

rho_zeta1 <- sapply(test_ws, rho_Wk, zeta = 5e-3, k = 0.05)
rho_zeta1_k2 <- sapply(test_ws, rho_Wk, zeta = 5e-3, k = 0.2)
rho_zeta2 <- sapply(test_ws, rho_Wk, zeta = 5e-4, k = 0.05)
rho_zeta2_k2 <- sapply(test_ws, rho_Wk, zeta = 5e-4, k = 0.2)

as.data.frame(cbind(test_ws,
                    rho_zeta1,
                    rho_zeta1_k2,
                    rho_zeta2,
                    rho_zeta2_k2)) %>% 
  gather("zeta_kap", "rho", rho_zeta1:rho_zeta2_k2) %>% 
  ggplot(aes(x = test_ws, y = rho, col = zeta_kap)) +
    geom_line(size = 1.2) +
    theme_classic() +
    theme(title = element_text(size = 18),
          axis.title = element_text(size = 15),
          axis.text = element_text(size = 12)) +
    scale_x_continuous(trans = "log",
                       breaks = c(0.0001, 0.001, 0.01, 0.1, 1, 10, 100),
                       labels = c("0.0001", "0.001", "0.01", "0.1", "1", "10", "100"),
                       limits = c(0.01, 200)) +
    ylim(c(0,1)) +
    labs(x = expression(Mean~Worm~Burden~(italic(W))),
         y = expression(paste("Relative Fecundity, ", italic(rho(W)))),
         title = "Relative fecundity across Mean Worm Burden") #+ theme(legend.position = "none")

test_ws <- seq(1e-2, sqrt(100), 
               by = sqrt(100*2)/1000)^2

phi_k0.1 <- sapply(test_ws, phi_Wk, k = 0.1)
phi_k1 <- sapply(test_ws, phi_Wk, k = 1)
phi_k10 <- sapply(test_ws, phi_Wk, k = 10)
phi_dynak <- mapply(phi_Wk, test_ws, sapply(test_ws, k_w_fx))

as.data.frame(cbind(test_ws,
                    phi_k0.1,
                    phi_k1,
                    phi_k10,
                    phi_dynak)) %>% 
  gather("clump_function", "Mate_Prob", phi_k0.1:phi_dynak) %>% 
  ggplot(aes(x = test_ws, y = Mate_Prob, col = clump_function)) +
    geom_line(size = 1.2) +
    theme_classic() +
    scale_x_continuous(trans = "log",
                       breaks = c(0.0001, 0.001, 0.01, 0.1, 1, 10, 100),
                       labels = c("0.0001", "0.001", "0.01", "0.1", "1", "10", "100"),
                       limits = c(0.01, 200)) +
    ylim(c(0,1)) +
    scale_color_manual(breaks = c("phi_k0.1", "phi_k1", "phi_k10", "phi_dynak"),
                       labels = c(expression(kappa~0.1),
                                  expression(kappa~1),
                                  expression(kappa~10),
                                  expression(kappa~(W))),
                       values = c("green", "blue", "red", "purple")) +
    labs(x = expression(Mean~Worm~Burden~(italic(W))),
         y = expression(Mating~Probability~(italic(Phi))),
         title = "Mating probability across clumping parameter",
         subtitle = "Influence of various formulations of clumping parameter") #+ theme(legend.position = "none")

Influence of the clumping parameter on the mating probability of adult female worms, the source of positive density dependence. When the clumping parameter is allowed to vary dynamically as a function of W, it can be seen that the mating probability remains higher into lower worm burdens, which will decrease the breakpoint population size and therefore make elimination with MDA more difficult

set.seed(10)

w_k_dist <- as.data.frame(expand.grid(samp.size = 1000,
                                      mu = c(1,10,100),
                                      kap = c(0.1,1,10))) %>% 
  mutate(merp = pmap(list("n" = samp.size, "size" = kap, "mu" = mu), rnbinom))

w_k_dist <- data.frame(mu = rep(c(rep(1, 1000),
                                  rep(10, 1000),
                                  rep(100, 1000)), 3),
                       k = c(rep(0.1, 3000),
                             rep(1, 3000),
                             rep(10, 3000)),
                       w = c(rnbinom(1000, mu = 1, size = 0.1),
                             rnbinom(1000, mu = 10, size = 0.1),
                             rnbinom(1000, mu = 100, size = 0.1),
                             rnbinom(1000, mu = 1, size = 1),
                             rnbinom(1000, mu = 10, size = 1),
                             rnbinom(1000, mu = 100, size = 1),
                             rnbinom(1000, mu = 1, size = 10),
                             rnbinom(1000, mu = 10, size = 10),
                             rnbinom(1000, mu = 100, size = 10)))

w_k_dist %>% ggplot(aes(x = w)) +
  geom_density(aes(fill = factor(k))) +
  facet_grid(.~mu) +
  scale_x_continuous(trans = "log",
                     breaks = c(1,10,100,1000)) +
  theme_classic() +
  scale_fill_manual(breaks = c("0.1", "1", "10"),
                    values = c("blue", "red", "purple")) +
  labs(fill = expression(kappa),
       x = expression(Individual~worm~burdens),
       title = "Distribution of individual worm burdens",
       subtitle = expression(Across~values~of~mean~worm~burden~(W)~and~clumping~parameter~(kappa)))

Influence of the clumping parameter on the distribution of individual worm burdens across different values of the mean worm burden. At low mean worm burdens (e.g. W=1, left panel), lower values of kappa implying more dispersed individual worm burdens are favorable for transmission as more individuals are likely to have >= 2 adult worms

Worm aggregation and acquired immunity

gamma_k0.1 <- sapply(test_ws, gam_Wxi, xi = 2.8e-2)

as.data.frame(cbind(test_ws,
                    gamma_k0.1)) %>% 
  ggplot(aes(x = test_ws, y = gamma_k0.1)) +
    geom_line(size = 1.2) +
    theme_classic() +
    theme(title = element_text(size = 15),
          axis.title = element_text(size = 15),
          axis.text = element_text(size = 12)) +
    scale_x_continuous(trans = "log",
                       breaks = c(0.0001, 0.001, 0.01, 0.1, 1, 10, 100),
                       labels = c("0.0001", "0.001", "0.01", "0.1", "1", "10", "100"),
                       limits = c(0.01, 200)) +
    ylim(c(0,1)) +
    labs(x = expression(Mean~Worm~Burden~(italic(W))),
         y = expression(paste("Relative infection probability, ", italic(gamma(W)))),
         title = "Relative infection probability across Mean Worm Burden") #+ theme(legend.position = "none")

Snail FOI as function of miracidia to snail pop density (M/N)

M_seq <- exp_seq(1e-4, 100, 300)

plot(M_seq, sapply(M_seq, function(x){1e-4*x}), type = "l", #ylim = c(0,1),
     lwd = 2, xlab = "Miracidia per snail (M/N)", ylab = "Snail FOI")
  lines(M_seq, sapply(M_seq, function(x){1e-2*(1-exp(-x))}), lwd = 2, col = 2)
  legend("bottomright", bty = "n", lwd = 2, col = c(1,2),
         legend = c("Linear", "Saturating"), title = expression(Snail~FOI~(Lambda)))

Initial Simulations and Checks

Test dynamic model 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(0:(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.8 # 80 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)

schisto_base_sim %>% 
  gather("Snail_Infection", "Density", S:I) %>% 
    ggplot(aes(x = time, y = Density, col = Snail_Infection)) +
      geom_line(size = 1.2) +
      scale_color_manual(values = c("S" = "green",
                                    "E" = "orange",
                                    "I" = "red")) +
      theme_classic() +
      ggtitle("Snail infection dynamics", 
              subtitle = paste0("Anual MDA for ", years/2, " years"))

schisto_base_sim %>% 
  mutate(W = Wt*base_pars["cvrg"] + Wu*(1-base_pars["cvrg"])) %>% 
  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)) +
      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"))

Test stochastic model through time

set.seed(10)
schisto_stoch_sim <- sim_schisto_stoch_mod(nstart = round(base_eqbm),
                                           params = as.list(base_pars),
                                           tf = max(base_time),
                                           events_df = mda.events) %>% 
  mutate(W = Wt*base_pars["cvrg"] + Wu*(1-base_pars["cvrg"]))

schisto_stoch_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.1, col = "purple") +
    scale_x_continuous(breaks = c(0:years)*365,
                       labels = c(-1:(years-1))) +
    theme_classic() +
    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"))

s_haem <- schisto_dat %>% 
  filter(Species == "S. haematobium" & 
           Mean_type == "arithmetic" & 
           !is.na(Intensity) &
           Intensity > 0) %>% 
  mutate(prev_sqrd = Prevalence^2,
         worm_est = (Intensity / 3.6)*2,  #Convert egg output to worm burden estimate assuming 3.6 eggs/mL per mated adult female worm
         clump_est = map2_dbl(worm_est, Prevalence, prev_intensity_to_clump))

s_haem %>% 
  ggplot(aes(x = Prevalence, y = Intensity, size = Sample_size)) + 
    geom_point() + 
    theme_classic() +
    labs(title = expression(Prevalence~intensity~relationship~italic(S.~haematobium)),
         subtitle = "Each point represents a population",
         x = "Prevalence (%)",
         y = "Egg Burden (eggs/10mL)") +
    geom_smooth(method = "lm", formula = y ~ x + I(x^2))

s_haem %>% 
  ggplot(aes(x = worm_est, y = clump_est, size = Sample_size)) + 
    geom_point() + 
    theme_classic() +
    labs(x = "Mean worm burden (W)",
         y = expression(Clumping~Parameter~(italic(kappa)))) +
    geom_smooth(method = "lm") +
    geom_smooth(method = "lm", formula = y~x+I(x^2))

clump_burden_lm_mod <- lm(clump_est ~ worm_est, data = s_haem)
clump_burden_lm_mod2 <- lm(clump_est ~ worm_est + I(worm_est^2), data = s_haem)

anova(clump_burden_lm_mod, clump_burden_lm_mod2)

#Linear is a better/simpler model

k_from_W <- function(W){
  as.numeric(coef(clump_burden_lm_mod)[1] + coef(clump_burden_lm_mod)[2] * W)
}

Generate $R_{eff}$ curve

#Get values of mean worm burden to plot R effective curve over
Reffs <- data.frame(test_ws = test_ws) %>% 
  mutate(kap = sapply(test_ws, k_from_W),
         Absent = pmap_dbl(list(parameters = list(base_pars),
                                     W = test_ws,
                                     kap = kap), Reff_W),
         Present = pmap_dbl(list(parameters = list(base_pars),
                                   W = test_ws,
                                   kap = kap), Reff_W))

Reffs %>% 
  gather("PDD", "Reff", Absent:Present) %>% 
  ggplot(aes(x = test_ws, y = Reff, lty = PDD)) +
  geom_line(size = 1.2) +
  theme_classic() +
  theme(title = element_text(size = 18),
        axis.title = element_text(size = 15),
        axis.text = element_text(size = 12),
        legend.position = c(0.2,0.6)) +
  geom_hline(yintercept = max(Reffs$Absent), lty = 2, col = 2) +
  scale_linetype_manual(values = c("dashed", "solid")) +
  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](W))),
       title = expression(italic(R[eff](W))~Curve)) 

Reffs %>% 
  ggplot(aes(x = test_ws, y = Present)) +
  geom_line(size = 1.2) +
  theme_classic() +
  theme(title = element_text(size = 18),
        axis.title = element_text(size = 15),
        axis.text = element_text(size = 12),
        legend.position = c(0.2,0.6)) +
  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:2),
                     limits = c(0,2)) +
  geom_hline(yintercept = 1, lty = 2) +
  labs(x = expression(Mean~Worm~Burden~(italic(W))),
       y = expression(italic(R[eff](W))),
       title = expression(italic(R[eff](W))~Curve)) 

#Get values of mean worm burden to plot R effective curve over
Reffs_dk <- data.frame(test_ws = test_ws) %>% 
  mutate(k0.1 = 0.1,
         k1 = 1,
         k10 = 10,
         w_det_ks = sapply(test_ws, k_from_log_W),
         Reff_k0.1 = pmap_dbl(list(parameters = list(base_pars),
                                   W = test_ws,
                                   kap = k0.1), getReff),
         Reff_k1 = pmap_dbl(list(parameters = list(base_pars),
                                   W = test_ws,
                                   kap = k1), getReff),
         Reff_k10 = pmap_dbl(list(parameters = list(base_pars),
                                   W = test_ws,
                                   kap = k10), getReff),
         Reff_no_pdd = pmap_dbl(list(parameters = list(base_pars),
                                     W = test_ws,
                                     kap = w_det_ks), getReff_noPDD),
         Reff_k_from_W = pmap_dbl(list(parameters = list(base_pars),
                                       W = test_ws,
                                       kap = w_det_ks), getReff))

Reffs_dk %>% 
  gather("Clumping_Function", "Reff", Reff_k0.1:Reff_k_from_W) %>% 
  ggplot(aes(x = test_ws, y = Reff, col = Clumping_Function)) +
  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)) +
  scale_color_manual(breaks = c("Reff_k0.1", "Reff_k1", "Reff_k10", "Reff_k_from_W", "Reff_no_pdd"),
                     labels = c(expression(kappa==0.1),
                                expression(kappa==1),
                                expression(kappa==10),
                                expression(kappa==f(W)),
                                expression(No~PDD)),
                     values = c("green", "blue", "red", "purple", "black")) +
  geom_hline(yintercept = 1, lty = 2) +
  labs(x = expression(Mean~Worm~Burden~(italic(W))),
       y = expression(italic(R[eff])),
       color = "Clumping\nParameter") #+ theme(legend.position = "none")

Influence of the clumping parameter on $R_{eff}$ curves. The dynamic clumping parameter (green line) makes things much worse, causing both a lower (harder to reach) breakpoint and a higher (more infection) endemic equilibrium)

See how $R_{eff}$ changes over course of deterministic simulation

schisto_base_sim_reff <- schisto_base_sim %>% 
  mutate(W = Wt*base_pars["cvrg"] + Wu*(1-base_pars["cvrg"]),
         kap = map_dbl(W, k_from_log_W),
         Reff = pmap_dbl(list(parameters = list(base_pars),
                              W = W,
                              kap = kap), getReff))

schisto_base_sim_reff %>% 
  ggplot(aes(x = time, y = Reff)) + 
  annotate("rect", xmin = 365, xmax = max(mda.events$time), ymin = -Inf, ymax = Inf,
           alpha = .2) +
  geom_line() + 
    theme_classic() +
    labs(x = "time (years)",
         y = expression(R[eff]),
         title = expression(Effective~reproduction~number~(R[eff])~over~MDA~campaign), 
         subtitle = paste0("Anual MDA for ", years/2, " years")) +
    scale_x_continuous(breaks = c(0:years)*365,
                       labels = c(-1:(years-1)))

See how $R_{eff}$ changes over course of stochastic simulation

schisto_stoch_sim_reff <- schisto_stoch_sim %>% 
  mutate(kap = map_dbl(W, k_from_log_W),
         Reff = pmap_dbl(list(parameters = list(base_pars),
                              W = W,
                              kap = kap), getReff))

schisto_stoch_sim_reff %>% 
  ggplot(aes(x = time, y = Reff)) + 
  annotate("rect", xmin = 365, xmax = max(mda.events$time), ymin = -Inf, ymax = Inf,
           alpha = .2) +
  geom_line() + 
    theme_classic() +
    labs(x = "time (years)",
         y = expression(R[eff]),
         title = "Effective reproduction number over MDA campaign", 
         subtitle = paste0("Anual MDA for ", years/2, " years")) +
    scale_x_continuous(breaks = c(0:years)*365,
                       labels = c(-1:(years-1)))

Clear that at 80% coverage and 93% efficacy, we never get close to the breakpoint with 10 years of annual MDA. This indicated by the rise in $R_{eff}$ with no subsequent decrease (implying traversing the curve down towards the breakpoint). However, this implies that the 20% of the untreated population contributes to transmission just like the 80% that's treated. In reality, MDA is often administered to school-aged children (SAC) while adults are untreated, but SAC also contribute much more to transmission. To exlpore the influence of this, we'll next move to the age-distributed model.

Effects of different interventions on $R_{eff}$ curve

Reffs_int <- data.frame(test_ws = test_ws) %>% 
  mutate(k0.1 = 0.1,
         k1 = 1,
         k10 = 10,
         w_det_ks = sapply(test_ws, k_from_log_W),
         Reff_k0.1 = pmap_dbl(list(parameters = list(base_pars),
                                   W = test_ws,
                                   kap = k0.1), getReff),
         Reff_k1 = pmap_dbl(list(parameters = list(base_pars),
                                   W = test_ws,
                                   kap = k1), getReff),
         Reff_k10 = pmap_dbl(list(parameters = list(base_pars),
                                   W = test_ws,
                                   kap = k10), getReff),
         Reff_no_pdd = pmap_dbl(list(parameters = list(base_pars),
                                     W = test_ws,
                                     kap = w_det_ks), getReff_noPDD),
         Reff_k_from_W = pmap_dbl(list(parameters = list(base_pars),
                                       W = test_ws,
                                       kap = w_det_ks), getReff))

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