--- title: "A model of wombat-mite dynamics among burrows" author: "Shane A. Richards" date: "12/11/2018" output: html_document: toc: true toc_float: true editor_options: chunk_output_type: console --- # Summary This document provides the host-disease modelling results presented in *Population-scale treatment informs solutions for control of environmentally transmitted wildlife disease* (Martin et al.). # The model See Supplememtary Materials for full details of the host-disease model. The code below provides the parameters for which we have some reasonable estimates. ```{r message=FALSE, warning=FALSE} rm(list = ls()) # clear memory # load useful packages library(tidyverse) library(DT) library(gridExtra) library(ggpubr) DaysSimulate <- as.integer(2.0*365) # days of simulation # Model parameters d <- 30 # mean days infected wombat to move to high infection [d] f <- 10 # mean days for infected burrow to lose infection [d] p <- 0.13 # probability wombat changes burrow each day [d-1] tau <- 547.5 # mean days between independence = births[d] L <- 10*365 # mean life-expectancy of wombats [d] LL <- L # mean life-expectancy with low infection [d] LH <- 90 # mean life-expectancy with high infection [d] = 90 vt <- 78 # total duration when treatment is applied [d] v <- 0.90 # proportion of burrows treated on a treatment day dv <- 1*7 # days between treatment days w <- 5 # days that the treatment remains effective pV <- 0.50 # daily probability treatment decays if not triggered qL <- 0.00 # probability L-wombat infects burrow # poorly known model parameters # qB <- 0.2 # probability burrow infects wombat # qH <- 0.25 # probability H-wombat infects burrow zS <- 0.333 # probability treatment is successfully applied to susceptible zL <- zS # probability treatment is successfully applied to low zH <- zS # probability treatment is successfully applied to high # derived model parameters pH <- 1.0/d # daily probability of low to high [d-1] pD <- 1.0/f # daily probability of disease loss per burrow [d-1] m0 <- 1.0/L # base-line mortality probability [d-1] mL <- 1.0/LL # low infection mortality probability [d-1] mH <- 1.0/LH # high infection mortality probability [d-1] pT <- 1.0/w # daily probability wombat loses effect of treatment [d-1] # set up vector of treatment fractions for each day v_day <- rep(0, DaysSimulate) v_day[1:vt] <- v*ifelse((0:(vt-1) %% dv) == 0, 1, 0) last_treatment <- max(which(v_day == v)) ``` Simulations are provided for three scenarios that are defined by lesser known parameter values. ```{r} # create a data frame with three base-line scenarios df_baseline <- tibble( scenario = c("A", "B", "C"), qB = c(0.500, 0.333, 0.200), # prob burrow transmits wombat qH = c(0.600, 0.500, 0.250), # prob high infects burrow zS = c(0.333, 0.250, 0.333), # prob treatment is applied f = c(10.00, 10.00, 20.00) # mean days disease persists in burrow ) datatable(df_baseline) # display paramter values that distinguish scenarios ``` In order to check that our model can provide reasonable predictions, we now compare them with observed wombats counts. ```{r message = FALSE} # read in the raw counts on day one of each survey df_plot_obs <- read_csv("WombatCounts.csv") df_plot_obs$Date <- as.Date(df_plot_obs$Date, format = "%d/%m/%y") ``` ```{r} # set the default parameter values for each scenario SetScenario <- function(scenario) { # Fixed model parameters d <<- 30 # mean days infected wombat to move to high infection [d] p <<- 0.13 # probability wombat changes burrow each day [d-1] tau <<- 547.5 # mean days between independence = births[d] L <<- 10*365 # mean life-expectancy of wombats [d] LL <<- L # mean life-expectancy with low infection [d] LH <<- 90 # mean life-expectancy with high infection [d] = 90 vt <<- 78 # total duration when treatment is applied [d] v <<- 0.90 # proportion of burrows treated on a treatment day dv <<- 1*7 # days between treatment days w <<- 5 # days that the treatment remains effective pV <<- 0.50 # daily probability treatment decays if not triggered if (scenario == "A") { i <- 1 } else if (scenario == "B") { i <- 2 } else { i <- 3 } # set some global parameters qB <<- df_baseline$qB[i] qH <<- df_baseline$qH[i] zS <<- df_baseline$zS[i] zL <<- df_baseline$zS[i] zH <<- df_baseline$zS[i] f <<- df_baseline$f[i] pD <<- 1.0/f # daily probability of disease loss per burrow [d-1] # set derived global derived model parameters pH <<- 1.0/d # daily probability of low to high [d-1] pD <<- 1.0/f # daily probability of disease loss per burrow [d-1] m0 <<- 1.0/L # base-line mortality probability [d-1] mL <<- 1.0/LL # low infection mortality probability [d-1] mH <<- 1.0/LH # high infection mortality probability [d-1] pT <<- 1.0/w # daily probability wombat loses effect of treatment [d-1] # set up global vector of treatment fractions for each day v_day <<- rep(0, DaysSimulate) v_day[1:vt] <<- v*ifelse((0:(vt-1) %% dv) == 0, 1, 0) last_treatment <<- max(which(v_day == v)) } # set the treatment transfer probabilities (equal for all wombat states) SetzS <- function(z_val) { zS <<- z_val zL <<- z_val zH <<- z_val } # set the treatment transfer probabilities (equal for all wombat states) Setw <- function(w_val) { w <<- w_val pT <<- 1.0/w_val # daily probability wombat loses effect of treatment [d-1] } # set the treatment duaration SetV <- function(wks) { vt <<- (7*wks-1) + 1 # 12 -> 78 # set up vector of treatment fractions for each day v_day <<- rep(0, DaysSimulate) v_day[1:vt] <<- v*ifelse((0:(vt-1) %% dv) == 0, 1, 0) last_treatment <<- max(which(v_day == v)) } # set the treatment duaration SetCoverage <- function(v) { vt <<- 78 # set up vector of treatment fractions for each day v_day <<- rep(0, DaysSimulate) v_day[1:vt] <<- v*ifelse((0:(vt-1) %% dv) == 0, 1, 0) last_treatment <<- max(which(v_day == v)) } ``` ```{r} # ================================================================ # The main function that takes the global parameters and simulates # wombat and disease dynamics # ================================================================ SimulateDynamics <- function() { # create a data frame to store the simulation predictions df_P <- tibble( # day-day time = 1:(DaysSimulate+1), Udt = 0.0, UdT = 0.0, UDt = 0.0, UDT = 0.0, Sdt = 0.0, SdT = 0.0, SDt = 0.0, SDT = 0.0, Idt = 0.0, IdT = 0.0, IDt = 0.0, IDT = 0.0, LDt = 0.0, LdT = 0.0, Ldt = 0.0, LDT = 0.0, Hdt = 0.0, HdT = 0.0, HDt = 0.0, HDT = 0.0) # within-day dynamics (8 steps) # 1: Death, 2: Burrows lose disease, 3: Wombat disease progression # 4: Wombats lose immunity, 5: Burrows lose tretament viability, # 6: Burrows gain treatment, 7: Wombats forage (may switch burrows) # 8: Wombat offspring gain independence df_p <- tibble( # within-day step = 1:8, Udt = 0.0, UdT = 0.0, UDt = 0.0, UDT = 0.0, Sdt = 0.0, SdT = 0.0, SDt = 0.0, SDT = 0.0, Idt = 0.0, IdT = 0.0, IDt = 0.0, IDT = 0.0, LDt = 0.0, LdT = 0.0, Ldt = 0.0, LDT = 0.0, Hdt = 0.0, HdT = 0.0, HDt = 0.0, HDT = 0.0) # Initial State df_P$Udt[1] <- 28/40 df_P$UdT[1] <- 0.00 df_P$UDt[1] <- 0.00 df_P$UDT[1] <- 0.00 df_P$Sdt[1] <- 3/40 df_P$SdT[1] <- 0.00 df_P$SDt[1] <- 0.00 df_P$SDT[1] <- 0.00 df_P$Idt[1] <- 0.00 df_P$IdT[1] <- 0.00 df_P$IDt[1] <- 0.00 df_P$IDT[1] <- 0.00 df_P$Ldt[1] <- 0.00 df_P$LdT[1] <- 0.00 df_P$LDt[1] <- 3/40 df_P$LDT[1] <- 0.00 df_P$Hdt[1] <- 0.00 df_P$HdT[1] <- 0.00 df_P$HDt[1] <- 6/40 df_P$HDT[1] <- 1.00 - sum(df_P[1,2:20]) # ensure probs sum to 1 # perform daily time-steps for (t in 1:DaysSimulate) { # (1) wombat dies df_p$Udt[1] <- df_P$Udt[t] + m0*df_P$Sdt[t] + m0*df_P$Idt[t] + mL*df_P$Ldt[t] + mH*df_P$Hdt[t] df_p$UdT[1] <- df_P$UdT[t] + m0*df_P$SdT[t] + m0*df_P$IdT[t] + mL*df_P$LdT[t] + mH*df_P$HdT[t] df_p$UDt[1] <- df_P$UDt[t] + m0*df_P$SDt[t] + m0*df_P$IDt[t] + mL*df_P$LDt[t] + mH*df_P$HDt[t] df_p$UDT[1] <- df_P$UDT[t] + m0*df_P$SDT[t] + m0*df_P$IDT[t] + mL*df_P$LDT[t] + mH*df_P$HDT[t] df_p$Sdt[1] <- (1.0 - m0)*df_P$Sdt[t] df_p$SdT[1] <- (1.0 - m0)*df_P$SdT[t] df_p$SDt[1] <- (1.0 - m0)*df_P$SDt[t] df_p$SDT[1] <- (1.0 - m0)*df_P$SDT[t] df_p$Idt[1] <- (1.0 - m0)*df_P$Idt[t] df_p$IdT[1] <- (1.0 - m0)*df_P$IdT[t] df_p$IDt[1] <- (1.0 - m0)*df_P$IDt[t] df_p$IDT[1] <- (1.0 - m0)*df_P$IDT[t] df_p$Ldt[1] <- (1.0 - mL)*df_P$Ldt[t] df_p$LdT[1] <- (1.0 - mL)*df_P$LdT[t] df_p$LDt[1] <- (1.0 - mL)*df_P$LDt[t] df_p$LDT[1] <- (1.0 - mL)*df_P$LDT[t] df_p$Hdt[1] <- (1.0 - mH)*df_P$Hdt[t] df_p$HdT[1] <- (1.0 - mH)*df_P$HdT[t] df_p$HDt[1] <- (1.0 - mH)*df_P$HDt[t] df_p$HDT[1] <- (1.0 - mH)*df_P$HDT[t] # cat(paste("Step 1: ", sum(df_p[1, 2:21]), "\n", sep = "")) # (2) disease persistence df_p$Udt[2] <- df_p$Udt[1] + pD*df_p$UDt[1] df_p$UDt[2] <- (1.0 - pD)*df_p$UDt[1] df_p$UdT[2] <- df_p$UdT[1] + pD*df_p$UDT[1] df_p$UDT[2] <- (1.0 - pD)*df_p$UDT[1] df_p$Idt[2] <- df_p$Idt[1] + pD*df_p$IDt[1] df_p$IDt[2] <- (1.0 - pD)*df_p$IDt[1] df_p$IdT[2] <- df_p$IdT[1] + pD*df_p$IDT[1] df_p$IDT[2] <- (1.0 - pD)*df_p$IDT[1] df_p$Sdt[2] <- df_p$Sdt[1] + pD*df_p$SDt[1] df_p$SDt[2] <- (1.0 - pD)*df_p$SDt[1] df_p$SdT[2] <- df_p$SdT[1] + pD*df_p$SDT[1] df_p$SDT[2] <- (1.0 - pD)*df_p$SDT[1] df_p$Ldt[2] <- df_p$Ldt[1] df_p$LdT[2] <- df_p$LdT[1] df_p$LDt[2] <- df_p$LDt[1] df_p$LDT[2] <- df_p$LDT[1] df_p$Hdt[2] <- df_p$Hdt[1] df_p$HdT[2] <- df_p$HdT[1] df_p$HDt[2] <- df_p$HDt[1] df_p$HDT[2] <- df_p$HDT[1] # cat(paste("Step 2: ", sum(df_p[2, 2:21]), "\n", sep = "")) # (3) within-host disease progression df_p$Udt[3] <- df_p$Udt[2] df_p$UDt[3] <- df_p$UDt[2] df_p$UdT[3] <- df_p$UdT[2] df_p$UDT[3] <- df_p$UDT[2] df_p$Idt[3] <- df_p$Idt[2] df_p$IDt[3] <- df_p$IDt[2] df_p$IdT[3] <- df_p$IdT[2] df_p$IDT[3] <- df_p$IDT[2] df_p$Sdt[3] <- df_p$Sdt[2] df_p$SdT[3] <- df_p$SdT[2] df_p$SDt[3] <- df_p$SDt[2] df_p$SDT[3] <- df_p$SDT[2] df_p$Ldt[3] <- (1.0 - pH)*df_p$Ldt[2] df_p$Hdt[3] <- df_p$Hdt[2] + pH*df_p$Ldt[2] df_p$LdT[3] <- (1.0 - pH)*df_p$LdT[2] df_p$HdT[3] <- df_p$HdT[2] + pH*df_p$LdT[2] df_p$LDt[3] <- (1.0 - pH)*df_p$LDt[2] df_p$HDt[3] <- df_p$HDt[2] + pH*df_p$LDt[2] df_p$LDT[3] <- (1.0 - pH)*df_p$LDT[2] df_p$HDT[3] <- df_p$HDT[2] + pH*df_p$LDT[2] # cat(paste("Step 3: ", sum(df_p[3, 2:21]), "\n", sep = "")) # (4) treatment loss for wombats df_p$Udt[4] <- df_p$Udt[3] df_p$UDt[4] <- df_p$UDt[3] df_p$UdT[4] <- df_p$UdT[3] df_p$UDT[4] <- df_p$UDT[3] df_p$Idt[4] <- (1.0 - pT)*df_p$Idt[3] df_p$Sdt[4] <- df_p$Sdt[3] + pT*df_p$Idt[3] df_p$IdT[4] <- (1.0 - pT)*df_p$IdT[3] df_p$SdT[4] <- df_p$SdT[3] + pT*df_p$IdT[3] df_p$IDt[4] <- (1.0 - pT)*df_p$IDt[3] df_p$SDt[4] <- df_p$SDt[3] + pT*df_p$IDt[3] df_p$IDT[4] <- (1.0 - pT)*df_p$IDT[3] df_p$SDT[4] <- df_p$SDT[3] + pT*df_p$IDT[3] df_p$Ldt[4] <- df_p$Ldt[3] df_p$LdT[4] <- df_p$LdT[3] df_p$LDt[4] <- df_p$LDt[3] df_p$LDT[4] <- df_p$LDT[3] df_p$Hdt[4] <- df_p$Hdt[3] df_p$HdT[4] <- df_p$HdT[3] df_p$HDt[4] <- df_p$HDt[3] df_p$HDT[4] <- df_p$HDT[3] # cat(paste("Step 4: ", sum(df_p[4, 2:21]), "\n", sep = "")) # (5) treatment loss for burrows df_p$Udt[5] <- df_p$Udt[4] + pV*df_p$UdT[4] df_p$UdT[5] <- (1.0 - pV)*df_p$UdT[4] df_p$UDt[5] <- df_p$UDt[4] + pV*df_p$UDT[4] df_p$UDT[5] <- (1.0 - pV)*df_p$UDT[4] df_p$Sdt[5] <- df_p$Sdt[4] + pV*df_p$SdT[4] df_p$SdT[5] <- (1.0 - pV)*df_p$SdT[4] df_p$SDt[5] <- df_p$SDt[4] + pV*df_p$SDT[4] df_p$SDT[5] <- (1.0 - pV)*df_p$SDT[4] df_p$Ldt[5] <- df_p$Ldt[4] + pV*df_p$LdT[4] df_p$LdT[5] <- (1.0 - pV)*df_p$LdT[4] df_p$LDt[5] <- df_p$LDt[4] + pV*df_p$LDT[4] df_p$LDT[5] <- (1.0 - pV)*df_p$LDT[4] df_p$Hdt[5] <- df_p$Hdt[4] + pV*df_p$HdT[4] df_p$HdT[5] <- (1.0 - pV)*df_p$HdT[4] df_p$HDt[5] <- df_p$HDt[4] + pV*df_p$HDT[4] df_p$HDT[5] <- (1.0 - pV)*df_p$HDT[4] df_p$Idt[5] <- df_p$Idt[4] + pV*df_p$IdT[4] df_p$IdT[5] <- (1.0 - pV)*df_p$IdT[4] df_p$IDt[5] <- df_p$IDt[4] + pV*df_p$IDT[4] df_p$IDT[5] <- (1.0 - pV)*df_p$IDT[4] # cat(paste("Step 5: ", sum(df_p[5, 2:21]), "\n", sep = "")) # (6) treatment gain for burrows: v_day[t] vv <- v_day[t] # fraction of burrows treated today df_p$UdT[6] <- df_p$UdT[5] + vv*df_p$Udt[5] df_p$Udt[6] <- (1.0 - vv)*df_p$Udt[5] df_p$UDT[6] <- df_p$UDT[5] + vv*df_p$UDt[5] df_p$UDt[6] <- (1.0 - vv)*df_p$UDt[5] df_p$SdT[6] <- df_p$SdT[5] + vv*df_p$Sdt[5] df_p$Sdt[6] <- (1.0 - vv)*df_p$Sdt[5] df_p$SDT[6] <- df_p$SDT[5] + vv*df_p$SDt[5] df_p$SDt[6] <- (1.0 - vv)*df_p$SDt[5] df_p$LdT[6] <- df_p$LdT[5] + vv*df_p$Ldt[5] df_p$Ldt[6] <- (1.0 - vv)*df_p$Ldt[5] df_p$LDT[6] <- df_p$LDT[5] + vv*df_p$LDt[5] df_p$LDt[6] <- (1.0 - vv)*df_p$LDt[5] df_p$HdT[6] <- df_p$HdT[5] + vv*df_p$Hdt[5] df_p$Hdt[6] <- (1.0 - vv)*df_p$Hdt[5] df_p$HDT[6] <- df_p$HDT[5] + vv*df_p$HDt[5] df_p$HDt[6] <- (1.0 - vv)*df_p$HDt[5] df_p$IdT[6] <- df_p$IdT[5] + vv*df_p$Idt[5] df_p$Idt[6] <- (1.0 - vv)*df_p$Idt[5] df_p$IDT[6] <- df_p$IDT[5] + vv*df_p$IDt[5] df_p$IDt[6] <- (1.0 - vv)*df_p$IDt[5] # cat(paste("Step 6: ", sum(df_p[6, 2:21]), "\n", sep = "")) # (7) Wombat burrow switching # calc burrow states after non-switching wombats return to burrow p_SdT <- 0.0 p_SDT <- 0.0 p_IdT <- 0.0 p_IDT <- 0.0 p_LdT <- 0.0 p_LDT <- 0.0 p_HdT <- 0.0 p_HDT <- 0.0 p_UdT <- df_p$UdT[6] p_UDT <- df_p$UDT[6] p_Udt <- df_p$Udt[6] + p*( df_p$Sdt[6] + df_p$Idt[6] + df_p$Ldt[6] + df_p$Hdt[6] + df_p$SdT[6] + df_p$IdT[6] + df_p$LdT[6] + df_p$HdT[6]) p_UDt <- df_p$UDt[6] + p*( df_p$SDt[6] + df_p$IDt[6] + df_p$LDt[6] + df_p$HDt[6] + df_p$SDT[6] + df_p$IDT[6] + df_p$LDT[6] + df_p$HDT[6]) p_Sdt <- (1.0 - p)*((1.0 - zS)*df_p$SdT[6] + df_p$Sdt[6]) p_SDt <- (1.0 - p)*(1.0 - qB)*((1.0 - zS)*df_p$SDT[6] + df_p$SDt[6]) p_Ldt <- (1.0 - p)*(1.0 - qL)*((1.0 - zL)*df_p$LdT[6] + df_p$Ldt[6]) p_LDt <- (1.0 - p)*(qL*df_p$Ldt[6] + df_p$LDt[6]) + (1.0 - p)*(1.0 - zL)*(qL*df_p$LdT[6] + df_p$LDT[6]) + (1.0 - p)*qB*((1.0 - zS)*df_p$SDT[6] + df_p$SDt[6]) p_Hdt <- (1.0 - p)*(1.0 - qH)*((1.0 - zH)*df_p$HdT[6] + df_p$Hdt[6]) p_HDt <- (1.0 - p)*(qH*df_p$Hdt[6] + df_p$HDt[6]) + (1.0 - p)*(1.0 - zH)*(qH*df_p$HdT[6] + df_p$HDT[6]) p_Idt <- (1.0 - p)*(df_p$Idt[6] + df_p$IdT[6]) + (1.0 - p)*(zS*df_p$SdT[6] + zL*df_p$LdT[6] + zH*df_p$HdT[6]) p_IDt <- (1.0 - p)*(df_p$IDt[6] + df_p$IDT[6]) + (1.0 - p)*(zS*df_p$SDT[6] + zL*df_p$LDT[6] + zH*df_p$HDT[6]) sum_all <- p_Udt + p_UdT + p_UDt + p_UDT + p_Sdt + p_SdT + p_SDt + p_SDT + p_Idt + p_IdT + p_IDt + p_IDT + p_Ldt + p_LdT + p_LDt + p_LDT + p_Hdt + p_HdT + p_HDt + p_HDT # cat(paste(" P prime sum: ", sum_all, "\n", sep = "")) # number of switchers looking for new burrow NS <- p*(df_p$Sdt[6] + df_p$SDt[6] + (1.0 - zS)*(df_p$SdT[6] + df_p$SDT[6])) NL <- p*(df_p$Ldt[6] + df_p$LDt[6] + (1.0 - zL)*(df_p$LdT[6] + df_p$LDT[6])) NH <- p*(df_p$Hdt[6] + df_p$HDt[6] + (1.0 - zH)*(df_p$HdT[6] + df_p$HDT[6])) NI <- p*(df_p$Idt[6] + df_p$IDt[6] + df_p$IdT[6] + df_p$IDT[6]) + p*(zS*df_p$SdT[6] + zL*df_p$LdT[6] + zH*df_p$HdT[6] + zS*df_p$SDT[6] + zL*df_p$LDT[6] + zH*df_p$HDT[6]) N <- NS + NL + NH + NI # total density of switchers # cat(paste(" Total switchers: ", N, "\n", sep = "")) # which burrow states do switchers encounter phi_dt <- p_Udt / (p_Udt + p_UdT + p_UDt + p_UDT) phi_dT <- p_UdT / (p_Udt + p_UdT + p_UDt + p_UDT) phi_Dt <- p_UDt / (p_Udt + p_UdT + p_UDt + p_UDT) phi_DT <- p_UDT / (p_Udt + p_UdT + p_UDt + p_UDT) # reduce unoccupied burrows df_p$Udt[7] <- p_Udt - phi_dt*N df_p$UdT[7] <- p_UdT - phi_dT*N df_p$UDt[7] <- p_UDt - phi_Dt*N df_p$UDT[7] <- p_UDT - phi_DT*N # increase occupied burrows by same amount df_p$Sdt[7] <- p_Sdt + NS*(phi_dt + phi_dT*(1.0 - zS)) df_p$SDt[7] <- p_SDt + NS*(1.0 - qB)*(phi_Dt + phi_DT*(1.0 - zS)) df_p$Ldt[7] <- p_Ldt + NL*(1.0 - qL)*(phi_dt + phi_dT*(1.0 - zL)) df_p$LDt[7] <- p_LDt + NS*qB*(phi_Dt + phi_DT*(1.0 - zS)) + NL*(phi_dt*qL + phi_dT*qL*(1.0 - zL) + phi_Dt + phi_DT*(1.0 - zL)) df_p$Hdt[7] <- p_Hdt + NH*(1.0 - qH)*(phi_dt + phi_dT*(1.0 - zH)) df_p$HDt[7] <- p_HDt + NH*(phi_dt*qH + phi_dT*qH*(1.0 - zH) + phi_Dt + phi_DT*(1.0 - zH)) df_p$Idt[7] <- p_Idt + NI*(phi_dt + phi_dT) + phi_dT*(zS*NS + zL*NL + zH*NH) df_p$IDt[7] <- p_IDt + NI*(phi_Dt + phi_DT) + phi_DT*(zS*NS + zL*NL + zH*NH) # not triggered and occupied should be zero df_p$SdT[7] <- p_SdT df_p$SDT[7] <- p_SDT df_p$IdT[7] <- p_IdT df_p$IDT[7] <- p_IDT df_p$LdT[7] <- p_LdT df_p$LDT[7] <- p_LDT df_p$HdT[7] <- p_HdT df_p$HDT[7] <- p_HDT # cat(paste("Step 7: ", sum(df_p[7, 2:21]), "\n", sep = "")) # (8) Wombat independence # densities of healthy and diseased independent wombats bU <- (0.5/tau) * ( df_p$Sdt[7] + df_p$SdT[7] + df_p$SDt[7] + df_p$SDT[7] + df_p$Idt[7] + df_p$IdT[7] + df_p$IDt[7] + df_p$IDT[7]) bD <- (0.5/tau) * ( df_p$Ldt[7] + df_p$LdT[7] + df_p$LDt[7] + df_p$LDT[7] + df_p$Hdt[7] + df_p$HdT[7] + df_p$HDt[7] + df_p$HDT[7]) # survival probability U <- df_p$Udt[7] + df_p$UdT[7] + df_p$UDt[7] + df_p$UDT[7] st <- min(bU + bD, U) / (bU + bD) # which burrow states do new independent wombats encounter? phi_dt <- df_p$Udt[7] / U phi_dT <- df_p$UdT[7] / U phi_Dt <- df_p$UDt[7] / U phi_DT <- df_p$UDT[7] / U df_p$Sdt[8] <- df_p$Sdt[7] + st*bU*(phi_dt + phi_dT*(1.0 - zS)) df_p$SDt[8] <- df_p$SDt[7] + st*bU*(1.0 - qB)*(phi_Dt + phi_DT*(1.0 - zS)) df_p$Ldt[8] <- df_p$Ldt[7] + st*bD*(1.0 - qL)*(phi_dt + phi_dT*(1.0 - zL)) df_p$LDt[8] <- df_p$LDt[7] + st*bU*qB*(phi_Dt + phi_DT*(1.0 - zS)) + st*bD*(phi_Dt + phi_DT*(1.0 - zL) + qL*phi_dt + qL*phi_dT*(1.0 - zL)) df_p$Idt[8] <- df_p$Idt[7] + st*phi_dT*(bU*zS + bD*zL) df_p$IDt[8] <- df_p$IDt[7] + st*phi_DT*(bU*zS + bD*zL) df_p$Udt[8] <- df_p$Udt[7] - st*phi_dt*(bU + bD) df_p$UdT[8] <- df_p$UdT[7] - st*phi_dT*(bU + bD) df_p$UDt[8] <- df_p$UDt[7] - st*phi_Dt*(bU + bD) df_p$UDT[8] <- df_p$UDT[7] - st*phi_DT*(bU + bD) df_p$SdT[8] <- df_p$SdT[7] df_p$SDT[8] <- df_p$SDT[7] df_p$IdT[8] <- df_p$IdT[7] df_p$IDT[8] <- df_p$IDT[7] df_p$LdT[8] <- df_p$LdT[7] df_p$LDT[8] <- df_p$LDT[7] df_p$Hdt[8] <- df_p$Hdt[7] df_p$HdT[8] <- df_p$HdT[7] df_p$HDt[8] <- df_p$HDt[7] df_p$HDT[8] <- df_p$HDT[7] # cat(paste("Step 8: ", sum(df_p[8, 2:21]), "\n", sep = "")) # transfer day values to results df_P[t+1, 2:21] <- df_p[8, 2:21] } # generate a data frame for plotting state dynamics df_plot <- gather(df_P, key = State, value = Probability, 2:21) df_plot$Occupant <- str_sub(df_plot$State, start = 1, end = 1) df_plot$Diseased <- str_sub(df_plot$State, start = 2, end = 2) df_plot$Treatment <- str_sub(df_plot$State, start = 3, end = 3) df_plot$State <- factor(df_plot$State, levels = c( "Udt", "UdT", "UDt", "UDT", "Sdt", "SdT", "SDt", "SDT", "Idt", "IdT", "IDt", "IDT", "Ldt", "LdT", "LDt", "LDT", "Hdt", "HdT", "HDt", "HDT")) df_plot$Occupant <- factor(df_plot$Occupant, levels = c("U", "S", "I", "L", "H")) df_plot_state <- df_plot pl_state <- ggplot(df_plot_state, aes(x = (time - 1)/365, y = Probability, color = Treatment)) + geom_line() + labs(x = "Years") + facet_grid(Occupant ~ Diseased) + theme_bw() # calculate some other states df_macro <- df_P %>% mutate( U = Udt + UdT + UDt + UDT, S = Sdt + SdT + SDt + SDT, I = Idt + IdT + IDt + IDT, L = Ldt + LdT + LDt + LDT, H = Hdt + HdT + HDt + HDT) %>% select(time, U, S, I, L, H) df_plot <- gather(df_macro, key = State, value = Probability, 2:6) df_plot$State <- factor(df_plot$State, levels = c("U", "S", "I", "L", "H")) df_plot_wombat <- df_plot pl_wombat <- ggplot(df_plot_wombat, aes(x = (time - 1)/365, y = Probability, color = State)) + geom_line() + scale_colour_manual(values=c("black", "forestgreen", "blue", "orange", "red")) + labs(x = "Years", color = "Burrow\nstate") + theme_bw() df_figure <- df_macro %>% mutate(Occupied = 1 - U) %>% select(-U) df_figure <- gather(df_figure, key = State, value = Probability, 2:6) df_figure$State <- factor(df_figure$State, levels = c("Occupied", "S", "I", "L", "H")) start_date <- "2015-09-16" # "2015-09-16" TreatmentEnd <- "2015-12-02" # last day of treatment last_recovered <- "2016-03-17" # last day observed recovered animals start_day <- 0 # "2015-09-16" TreatmentEnd_day <- 77 # last day of treatment last_recovered_day <- 146 # last day observed recovered animals df_plot <- df_macro %>% mutate(N = 1 - U) %>% dplyr::select(time, N, H, I) df_plot <- gather(df_plot, key = Parameter, value = Probability, 2:4) df_plot$Parameter <- factor(df_plot$Parameter, levels = c("N", "H", "I")) df_plot$Date <- as.Date(df_plot$time - 1, origin = start_date) scalar <- 1/50 df_plot_obs$Total_scaled <- scalar*df_plot_obs$Total df_plot_obs$Infected_scaled <- scalar*df_plot_obs$Infected pl_compare <- ggplot(filter(df_plot, Parameter %in% c("N", "H")), aes(x = Date, y = Probability, color = Parameter)) + annotate("rect", fill = "blue", alpha = 0.25, xmin = as.Date(start_date), xmax = as.Date(TreatmentEnd), ymin = 0, ymax = Inf) + geom_line() + scale_colour_manual(values=c("black", "red")) + geom_point(data = df_plot_obs, aes(x = Date, y = Total_scaled), color = "black", inherit.aes = FALSE) + geom_point(data = df_plot_obs, aes(x = Date, y = Infected_scaled), color = "red", inherit.aes = FALSE) + labs(x = "Years", color = "Abundance", shape = "State") + guides(color=FALSE) + theme_bw() df_figure$Date <- as.Date(df_figure$time - 1, origin = start_date) df_figure$Day <- df_figure$time - 1 pl_burrow <- ggplot(df_figure, aes(x = Day, y = Probability, color = State)) + annotate("rect", fill = "blue", alpha = 0.25, xmin = start_day, xmax = TreatmentEnd_day, ymin = 0, ymax = Inf) + annotate("rect", fill = "grey87", xmin = TreatmentEnd_day, xmax = last_recovered_day, ymin = 0, ymax = Inf) + geom_line() + scale_colour_manual(values=c("black", "forestgreen", "blue", "orange", "red")) + geom_point(data = df_plot_obs, aes(x = Day, y = Total_scaled), color = "black", inherit.aes = FALSE) + labs(x = "Date", color = "Burrow\nstate") + scale_x_continuous( breaks=c(0, 175, 350, 525, 700), labels=c("16-Sept-16\nWeek 0", "09-Mar-16\nWeek 25", "31-Aug-16\nWeek 50", "22-Feb-17\nWeek 75", "16-Aug-17\nWeek 100")) + theme_bw() + theme( axis.text = element_text(size=6), # axis.text.x = element_text(angle=35, hjust=1, vjust=1), axis.title = element_text(size=9), plot.subtitle = element_text(size=12), panel.grid = element_blank()) # calculate some other states df_disease <- df_P %>% mutate( d = Udt + UdT + Sdt + SdT + Ldt + LdT + Hdt + HdT + Idt + IdT, D = UDt + UDT + SDt + SDT + LDt + LDT + HDt + HDT + IDt + IDT ) %>% select(time, d, D) df_probs <- df_macro %>% mutate( N = 1.0 - U, # wombat population density fracWDiseased = (L+H)/(1.0-U), # fraction of wombats diseased fracWTreated = I / N # fraction of wombats immune ) df_probs$fracBDiseased = df_disease$D df_probs <- df_probs %>% select(time, N, fracWDiseased, fracWTreated, fracBDiseased) # extract the disease prevalence aming burrows at end of treatment disease_prevalence <- df_probs$fracBDiseased[last_treatment] df_plot <- gather(df_probs, key = State, value = Probability, 2:5) df_plot$State <- factor(df_plot$State, levels = c("N", "fracWDiseased", "fracWTreated", "fracBDiseased")) df_plot$Date <- as.Date(df_plot$time - 1, origin = start_date) pl_disease <- ggplot(df_plot, aes(x = Date, y = Probability, color = State)) + annotate("rect", fill = "blue", alpha = 0.25, xmin = as.Date(start_date), xmax = as.Date(TreatmentEnd), ymin = 0, ymax = Inf) + geom_line() + scale_colour_manual( labels = c("Wombat density", "Disease prevalence\namong wombats", "Treatment prevalence\namong wombats", "Disease prevalence\namong burrows"), values=c("black", "red", "blue", "orange") ) + labs(x = "Years", color = "Variable") + theme_bw() + theme(panel.grid = element_blank()) # store all the relevant output into a list return (list( df_P = df_P, pl_state = pl_state, pl_wombat = pl_wombat, pl_compare = pl_compare, pl_burrow = pl_burrow, pl_disease = pl_disease, prevalence = disease_prevalence) ) } ``` # Results First, we check that our model can make host-disease predictions that are consistent with the observed data on wombat abundance. Next, we ask what would happen to disease prevalence among burrows at the end of the treatment program if four of the model's parameters are changed; changes to these parameters reflect a change in the management strategy. ## Observations and predictions ```{r fig.width=9, fig.height=3, eval = TRUE} SetScenario("A") l_A1 <- SimulateDynamics() SetScenario("B") l_B1 <- SimulateDynamics() SetScenario("C") l_C1 <- SimulateDynamics() ggarrange( l_A1$pl_burrow + theme(legend.position=c(0.99,0.99), legend.justification=c(1,1), legend.text=element_text(size=8), legend.title=element_text(size=10), legend.spacing.y=unit(0.25,"line")), l_B1$pl_burrow + guides(color=FALSE), l_C1$pl_burrow + guides(color=FALSE), ncol = 3, nrow = 1, labels = "AUTO", label.y = 1) ``` These three scenarios result in simulated host numbers that are consistent with the observations on wombat abundance in the field. ## Predicting management outcomes Below we simulate disease dynamics when key management parameters are manipulated. ### Treatment delivery success Change the probability that treatment dose is successfully applied during flap triggering. ```{r fig.width=6, fig.height=3, eval = TRUE} vec_zS <- c(0.1, 0.15, 0.2, 0.25, 0.3, 0.333, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1.0) prev_A <- rep(0, length(vec_zS)) prev_B <- rep(0, length(vec_zS)) prev_C <- rep(0, length(vec_zS)) SetScenario("A") for (i in 1:length(vec_zS)) { SetzS(vec_zS[i]) prev_A[i] <- SimulateDynamics()$prevalence } SetScenario("B") for (i in 1:length(vec_zS)) { SetzS(vec_zS[i]) prev_B[i] <- SimulateDynamics()$prevalence } SetScenario("C") for (i in 1:length(vec_zS)) { SetzS(vec_zS[i]) prev_C[i] <- SimulateDynamics()$prevalence } df_results <- tibble( Scenario = sort(rep(c("A", "B", "C"), length(vec_zS))), zS = rep(vec_zS, 3), prev = c(prev_A, prev_B, prev_C) ) df_reference <- df_results %>% filter(zS %in% c(0.25, 0.333)) df_reference <- df_reference[c(2,3,6), ] plot_zS <- ggplot(df_results, aes(x = zS, y = prev, color = Scenario)) + geom_point() + geom_line() + geom_point(data = df_reference, aes(x = zS, y = prev), color = "black", size = 2, inherit.aes = FALSE) + labs( x = "Probability treatment is successfully applied, zS", y = "Disease prevalence") + theme_bw() + theme( panel.grid = element_blank(), axis.title = element_text(size = 11), axis.text = element_text(size = 9), legend.position=c(0.95,0.95), legend.justification=c(1,1) ) plot_zS df_reference ``` ### Treatment longevity Change the mean number of days that a wombat is immune to mite infection due to the treatment. ```{r fig.width=6, fig.height=3, eval = TRUE} # set range of treatment durations vec_w <-c(c(2.5, 5.0, 7.5, 10, 15), seq(from = 20, to = 60, by = 10)) # where to store the disease prevalances at end of treatment prev_A <- rep(0, length(vec_w)) prev_B <- rep(0, length(vec_w)) prev_C <- rep(0, length(vec_w)) SetScenario("A") # simulate scenario A for (i in 1:length(vec_w)) { Setw(vec_w[i]) prev_A[i] <- SimulateDynamics()$prevalence } SetScenario("B") for (i in 1:length(vec_w)) { Setw(vec_w[i]) prev_B[i] <- SimulateDynamics()$prevalence } SetScenario("C") for (i in 1:length(vec_w)) { Setw(vec_w[i]) prev_C[i] <- SimulateDynamics()$prevalence } # create a data frame with the resuls in long format df_results <- tibble( Scenario = sort(rep(c("A", "B", "C"), length(vec_w))), w = rep(vec_w, 3), prev = c(prev_A, prev_B, prev_C) ) df_reference <- df_results %>% filter(w == 5) plot_w <- ggplot(df_results, aes(x = w, y = prev, color = Scenario)) + geom_point() + geom_line() + geom_point(data = df_reference, aes(x = w, y = prev), color = "black", size = 2, inherit.aes = FALSE) + labs( x = "Treatment efficacy, w (days)", y = "Disease prevalence") + guides(color = FALSE) + theme_bw() + theme( panel.grid = element_blank(), axis.title = element_text(size = 11), axis.text = element_text(size = 9) ) plot_w df_reference ``` ### Treatment duration Change the number of weeks that the treatment program is applied. ```{r eval = TRUE} # set range of treatment durations vec_weeks <- seq(from = 6, to = 20, by = 2) # where to store the disease prevalances at end of treatment prev_A <- rep(0, length(vec_weeks)) prev_B <- rep(0, length(vec_weeks)) prev_C <- rep(0, length(vec_weeks)) SetScenario("A") # simulate scenario A for (i in 1:length(vec_weeks)) { SetV(vec_weeks[i]) prev_A[i] <- SimulateDynamics()$prevalence } SetScenario("B") for (i in 1:length(vec_weeks)) { SetV(vec_weeks[i]) prev_B[i] <- SimulateDynamics()$prevalence } SetScenario("C") for (i in 1:length(vec_weeks)) { SetV(vec_weeks[i]) prev_C[i] <- SimulateDynamics()$prevalence } # create a data frame with the resuls in long format df_results <- tibble( Scenario = sort(rep(c("A", "B", "C"), length(vec_weeks))), wks = rep(vec_weeks, 3), prev = c(prev_A, prev_B, prev_C) ) df_reference <- df_results %>% filter(wks == 12) plot_wks <- ggplot(df_results, aes(x = wks, y = prev, color = Scenario)) + geom_point() + geom_line() + geom_point(data = df_reference, aes(x = wks, y = prev), color = "black", size = 2, inherit.aes = FALSE) + labs( x = "Weeks treatment is applied", y = "Disease prevalence") + # xlim(0,20) + guides(color = FALSE) + theme_bw() + theme( panel.grid = element_blank(), axis.title = element_text(size = 11), axis.text = element_text(size = 9) ) plot_wks df_reference ``` ### Treatment coverage Change the proportion of active burrows that have treatment delivering flaps applied. ```{r eval = TRUE} # set range of treatment durations vec_v <- seq(from = 0.4, to = 1.0, by = 0.1) # where to store the disease prevalances at end of treatment prev_A <- rep(0, length(vec_v)) prev_B <- rep(0, length(vec_v)) prev_C <- rep(0, length(vec_v)) SetScenario("A") # simulate scenario A for (i in 1:length(vec_v)) { SetCoverage(vec_v[i]) prev_A[i] <- SimulateDynamics()$prevalence } SetScenario("B") for (i in 1:length(vec_v)) { SetCoverage(vec_v[i]) prev_B[i] <- SimulateDynamics()$prevalence } SetScenario("C") for (i in 1:length(vec_v)) { SetCoverage(vec_v[i]) prev_C[i] <- SimulateDynamics()$prevalence } # create a data frame with the resuls in long format df_results <- tibble( Scenario = sort(rep(c("A", "B", "C"), length(vec_v))), v = rep(vec_v, 3), prev = c(prev_A, prev_B, prev_C) ) df_reference <- df_results %>% filter(v > 0.85 & v < 0.95) plot_v <- ggplot(df_results, aes(x = v, y = prev, color = Scenario)) + geom_point() + geom_line() + geom_point(data = df_reference, aes(x = v, y = prev), color = "black", size = 2, inherit.aes = FALSE) + labs( x = "Proportion of burrows treated", y = "Disease prevalence") + guides(color = FALSE) + theme_bw() + theme( panel.grid = element_blank(), axis.title = element_text(size = 11), axis.text = element_text(size = 9) ) plot_v df_reference ``` ## All management options Combine results from all four management modifications into a single figure. ```{r eval = TRUE} p_all <- ggarrange( plot_zS + ylim(0,0.16) + xlim(0.15,1), plot_w + ylim(0,0.16) + xlim(0,60), plot_wks + ylim(0,0.16) + xlim(5,20), plot_v + ylim(0,0.16) + xlim(0.4,1), ncol = 2, nrow = 2, # common.legend = FALSE, legend = "none", labels = "AUTO", label.y = 1) p_all ```