Skip to contents

SARS-CoV-2 outbreak characteristics across wild and captive settings

Analysis referenced in Rosenblatt et al. 2023

Elias Rosenblatt

Source: vignettes/SIRS_analysis_by_contexts.Rmd
SIRS_analysis_by_contexts.Rmd

Introduction

This vignette details the analysis used to study differences in outbreak dynamics in four scenarios of SARS-CoV-2 transmission to and among white-tailed deer populations (Rosenblatt et al. In Prep.). Briefly, this study evaluates the risk deer face of direct transmission from humans in different scenarios. Once the disease is introduced, this study also describes the outbreak dynamics within the deer population. The use of whitetailedSIRS functions and workflow mimic the previous vignette demonstrating the versitility of the whitetailedSIRS package. The output of this vignette is used to visualize differences in the vignette whitetailedSIRS::Visualize_by_context.

Again, general approach with simulating SARS-CoV-2 outbreaks here is to:

  1. Set the number of simulations

  2. Draw a random sample of parameters used to estimate infection probabilities, contact rates, and other important parameters.

  3. Derive parameters used the SIRS ODE equation set

  4. Define initial SIR compartment sizes

  5. Solve the SIRS ODE equations and calculate daily compartment sizes

  6. Calculate risk of introduction and magnitude of spread given parameters

  7. Calculate average daily prevalence, probability of persistence, and incidence proportion.

1. Set number of simulations to run and context of simulation

For the purpose of the study, Rosenblatt et al. (In Prep) focused on four scenarios to simulated SARS-CoV-2 introduction and outbreak. These simulations occur in isolation, meaning that an outbreak cannot spread from captive to wild or from wild to captive. Animals in a scenario are either all captive or all wild. Combined systems are address in the “Connected systems” vignette.

Wild deer in rural conditions: In this scenario, wild deer are free-ranging in an area with a rural human density (3.1 humans/km2). We assumed that deer in this context were subject to regulated hunting either using still-hunting, or ground blind or treestand approaches. We also assumed that baiting and backyard feeding were illegal but may still occur. We again used Habib et al.’s (2011) proximity rate model. We estimated the duration of deer-deer proximity events, and rates and duration of human-deer proximity events using expert elicitation.

Wild deer in suburban conditions: Wild deer are free-ranging in an area of suburban human density (100 humans/km2). Deer-to-deer proximity rates were derived using the same parameters as used in the rural setting. The primary difference between this context and the wild, rural context are higher human-to-deer proximity rates and longer duration of proximity events as estimated by expert elicitation.

Outdoor captive ranch: We modeled a population of captive deer housed in an outdoor ranching facility. We assumed proximity rates between deer in this context was the same as estimated in wild contexts, with the increase of these proximity rates due to supplemental feeding. We assume that proximity rates between humans and deer in this context are the same as those estimated for the wild, suburban context, but the duration of these proximity events are longer, reflecting those typical of a captive facility.

Intensive captive facility: We considered captive deer held in a facility focused on captive breeding or exposition. Deer in this facility were predominantly indoors with high stocking densities and low indoor air exchange rates (AER; 1-hr). We estimated all proximity rates and duration of proximity using expert elicitation.

We set a seed so the results here are repeatable, and set the number of samples to the desired amount (nsamples). For this example, we run 100 simulations.

set.seed(23)
nsamples <- 100

2. Sample parameters for simulations

Next, we draw an assortment of parameters to solve the SIRS ODE equations for a SARS-CoV-2 outbreak. If you have not read them, see the SIR model description and example vignettes to understand the workflow here.

elicitation_data <- whitetailedSIRS::draw_elicitation_samples(nsamples = nsamples)

Following Rosenblatt et al. (In Prep), four contexts were defined, including outdoor captive ranch, intensive captive facility, wild deer in rural conditions, and wild deer in suburban conditions. In this RMD file, these scenarios will be simulated independently of each other.

Outdoor captive ranch

Proximity Rates

Proximity rates were derived from a rate model developed by Habib et al. (2011), as a function of habitat conditions (habitat) and density (nwild/A_W). For all uses of the Habib et al. (2011) model in this analysis, we fixed density at 10 deer per sq km and habitat at 26% wooded habitat (classified in the calc_contact_rate function as type_contact = “med”). We also used expert elicitation estimates of the probability of direct contact for fluid transmission.

nWild <- rpois(nsamples,1000) #Abundance
A_w <- 100 #Area
habitat <- "med" #Habitat classification
sigma_season <- 1 #Season adjustment for proximity rate

epsilon_dc <- get_EE_param_vals(data = elicitation_data, my_param = "Direct Contact Probability") #Probability of direct contact between deer, given proximity.

Next, we define various proximity rates for the ranch context. Object names reflect other parameters that must be included in the SIRS ODE solver (e.g. omega_ww_ranch, set to 0).

omega_ww_ranch <- rep(0, nsamples) #Deer-to-deer proximity rate in wild (set to 0; events per day).
omega_cw_null <- rep(0, nsamples) #Deer-to-deer proximity rate along fenceline (set to 0; events per day).
omega_cc_ranch <- calc_contact_rate(nsamples = nsamples, type_contact = habitat, N_w = nWild, rho_attractant = get_EE_param_vals(data = elicitation_data, my_param = "Proximity rate with baiting (17 events without baiting)")/17) #Deer-to-deer proximity rate in ranch context, mimicing wild proximity rates with the influence of baiting (events per day).

omega_hw_ranch <- rep(0, nsamples) #Human-to-deer proximity rate in wild (set to 0; events per day).
omega_hc_ranch <- get_EE_param_vals(data = elicitation_data, my_param = "Deer-Human Proximity Rate, Suburban (per 120 days)") /120 #Human-to-deer proximity rate in ranch context, mimicing suburban proximitty rates (events per day).

Infection probabilities

We calculate infection probabilities for three transmission pathways - aerosol transmission from deer-to-deer, fluid transmission from deer-to-deer, and aerosol transmission from human-to-deer. Below, we will use various package functions to define parameters to estimate the infection probability given proximity or contact.

C_nu_human <- rnorm(n = nsamples, mean = 10^5.6, sd = 10^1.2)#viral load in humans (genomic copies per ml)
C_nu_deer <- 10^5.6 * get_EE_param_vals(data = elicitation_data, my_param = "Viral Load") #viral load in deer saliva, relative to humans (genomic copies per ml)

r_deer <- get_EE_param_vals(data = elicitation_data, my_param = "Dose-Response")# Dose response coefficient for deer and SARS-CoV-2
#Infection probability calculation of aerosol transmission from deer-to-deer
t_contact_deer_deer_null <- get_EE_param_vals(data = elicitation_data, my_param = "Deer Proximity Duration (minutes)") #Estimate duration of proximity event...
sigma_aero_deer_deer_ranch <- calc_sigma_aero(C_nu = C_nu_deer,
                                        t_contact = t_contact_deer_deer_null / 60,
                                        r = r_deer, nsamples = nsamples) #...and estimate probability of infection given that duration of proximity event.

sigma_aero_deer_deer_wild_null <- rep(0, nsamples) #Estimate infection probability in out in the wild as 0 (needs to be included for SIRS ODE equations)

#Infection probability of 0.1 ml of saliva being transferred between deer on contact
sigma_dc_deer_deer_null <- calc_sigma_dc(C_nu = C_nu_deer, nsamples = nsamples) #Calculate infection probability

#Infection probability calculation of aerosol transmission from humans-to-deer
t_contact_deer_human_ranch <- get_EE_param_vals(data = elicitation_data, my_param = "Deer-Human Proximity Duration, Captive (minutes)") #Estimate duration of human-deer proximity event in ranch facility context...
sigma_aero_deer_human_ranch <- calc_sigma_aero(ER = 0.53, C_nu = C_nu_human, 
                                        t_contact = t_contact_deer_human_ranch / 60,
                                        r = r_deer, nsamples = nsamples)#... and calculate infection probability given the duration of a human-deer proximity event.
sigma_aero_deer_human_wild_null <- rep(0, nsamples)#Estimate human-to-deer infection probability out in the wild as 0 (needs to be included for SIRS ODE equations)

Intensive captive facility

Below, we list additional parameters used for the “Intensive captive facility” scenario. Note that this scenario will use objects defined earlier in this RMarkdown.

omega_ww_intensive <- rep(0, nsamples)
omega_cc_intensive <- get_EE_param_vals(data = elicitation_data, my_param = "Deer-Deer Proximity Rate, Captive (per day)")
omega_hw_intensive <- rep(0, nsamples)
omega_hc_intensive <- get_EE_param_vals(data = elicitation_data, my_param = "Deer-Human Proximity Rate, Suburban (per 120 days)") /120

sigma_aero_deer_deer_intensive <- calc_sigma_aero(C_nu = C_nu_deer,
                                            t_contact = t_contact_deer_deer_null / 60,
                                            r = r_deer, nsamples = nsamples, AER = rep(1, nsamples))

t_contact_deer_human_intensive <- get_EE_param_vals(data = elicitation_data, my_param = "Deer-Human Proximity Duration, Captive (minutes)") 

sigma_aero_deer_human_intensive <- calc_sigma_aero(ER = 0.53, C_nu = C_nu_human, 
                                            t_contact = t_contact_deer_human_intensive / 60,
                                            r = r_deer, nsamples = nsamples, AER = rep(1, nsamples))

Wild deer in rural conditions

Again, we list additional parameters used for the “Wild deer in rural conditions” scenario. Note that this scenario will use objects defined earlier in this RMarkdown.

omega_ww_rural <- calc_contact_rate(nsamples = nsamples, type_contact = habitat, N_w = nWild)
omega_cc_rural <- rep(0, nsamples)
omega_hw_rural <- get_EE_param_vals(data = elicitation_data, my_param = "Deer-Human Proximity Rate, Rural (per 120 days)") /120
omega_hc_rural <- rep(0, nsamples)

t_contact_deer_human_rural <- get_EE_param_vals(data = elicitation_data, my_param = "Deer-Human Proximity Duration, Rural (minutes)")

sigma_aero_deer_deer_rural <- calc_sigma_aero(C_nu = C_nu_deer,
                                        t_contact = t_contact_deer_deer_null / 60,
                                        r = r_deer, nsamples = nsamples, AER = rep(4, nsamples))

sigma_aero_deer_deer_captive_null <- rep(0, nsamples)

sigma_aero_deer_human_rural <- calc_sigma_aero(ER = 0.53, C_nu = C_nu_human, 
                                        t_contact = t_contact_deer_human_rural / 60,
                                        r = r_deer, nsamples = nsamples, AER = rep(4, nsamples))

sigma_aero_deer_human_captive_null <- rep(0, nsamples)

Wild deer in suburban conditions

Again, we list additional parameters used for the “Wild deer in suburban conditions” scenario. Note that this scenario will use objects defined earlier in this RMarkdown.

omega_ww_suburban <- calc_contact_rate(nsamples = nsamples, type_contact = habitat, N_w = nWild)
omega_cc_suburban <- rep(0, nsamples)
omega_hw_suburban <- get_EE_param_vals(data = elicitation_data, my_param = "Deer-Human Proximity Rate, Suburban (per 120 days)") /120
omega_hc_suburban <- rep(0, nsamples)

t_contact_deer_human_suburban <- get_EE_param_vals(data = elicitation_data, my_param = "Deer-Human Proximity Duration, Suburban (minutes)")

sigma_aero_deer_deer_suburban <- calc_sigma_aero(C_nu = C_nu_deer,
                                        t_contact = t_contact_deer_deer_null / 60,
                                        r = r_deer, nsamples = nsamples, AER = rep(4, nsamples))

sigma_aero_deer_human_suburban <- calc_sigma_aero(ER = 0.53, C_nu = C_nu_human, 
                                        t_contact = t_contact_deer_human_suburban / 60,
                                        r = r_deer, nsamples = nsamples, AER = rep(4, nsamples))

Recovery and Immunity Loss Rates

For all simulations, animals recover from their infection at a given rate (gamma_recov), and lose their immunity and become susceptible at a given rate (alpha_immunity). We use a recovery rate based on a 6 day infectious period (Palmer et al. 2021), and immunity loss rate based on expert elicited estimates. The recovery rate is repeated across all iterations of a simulation (nsamples). The immunity loss rate is derived from [nsamples] draws from an error distribution estimated by the expert elicitation process.

gamma_recov <- rep(1/6, nsamples)
alpha_immunity_null <- 1 / get_EE_param_vals(data = elicitation_data, my_param = 'Temporary Immunity')
I_human_null <- 0.05

3. Derive parameters used in the SIRS ODE euqation set.

Next, we create derived parameters for each scenario:

ranch.params <- alternative(alpha_immunity = alpha_immunity_null,
                            omega_ww = omega_ww_ranch, omega_cw = omega_cw_null,omega_cc = omega_cc_ranch,
                            omega_hw = omega_hw_ranch, omega_hc = omega_hc_ranch,
                            sigma_aero_deer_deer_wild = sigma_aero_deer_deer_wild_null, sigma_aero_deer_deer_captive = sigma_aero_deer_deer_ranch, sigma_aero_deer_human_wild = sigma_aero_deer_human_wild_null, sigma_aero_deer_human_capt = sigma_aero_deer_human_ranch, epsilon_dc = epsilon_dc, sigma_dc_deer_deer = sigma_dc_deer_deer_null, gamma_recov = gamma_recov, I_human = rep(I_human_null, nsamples), boost = rep(0, nsamples))

intensive.params <- alternative(
   omega_ww = omega_ww_intensive, omega_cw = omega_cw_null, omega_cc = omega_cc_intensive, 
   omega_hw = omega_hw_intensive, omega_hc = omega_hc_intensive, 
   sigma_aero_deer_deer_wild = sigma_aero_deer_human_wild_null, sigma_aero_deer_deer_captive =  sigma_aero_deer_deer_intensive, sigma_aero_deer_human_wild = sigma_aero_deer_human_wild_null, sigma_aero_deer_human_capt = sigma_aero_deer_human_intensive, sigma_dc_deer_deer = sigma_dc_deer_deer_null,
   alpha_immunity = alpha_immunity_null, epsilon_dc = epsilon_dc, gamma_recov = gamma_recov, I_human = rep(I_human_null, nsamples), boost = rep(0, nsamples))

rural.params <- alternative(omega_ww = omega_ww_rural, omega_cw = omega_cw_null, omega_cc = omega_cc_rural, omega_hw = omega_hw_rural, omega_hc = omega_hc_rural, sigma_aero_deer_deer_wild = sigma_aero_deer_deer_rural, sigma_aero_deer_deer_captive = sigma_aero_deer_deer_captive_null, sigma_aero_deer_human_wild = sigma_aero_deer_human_rural,  sigma_aero_deer_human_capt = sigma_aero_deer_human_captive_null, sigma_dc_deer_deer = sigma_dc_deer_deer_null,alpha_immunity = alpha_immunity_null, epsilon_dc = epsilon_dc, gamma_recov = gamma_recov, I_human = rep(I_human_null, nsamples), boost = rep(0, nsamples))

suburban.params <- alternative(omega_ww = omega_ww_suburban, omega_cw = omega_cw_null, omega_cc = omega_cc_suburban, omega_hw = omega_hw_suburban, omega_hc = omega_hc_suburban, sigma_aero_deer_deer_wild = sigma_aero_deer_deer_suburban, sigma_aero_deer_deer_captive = sigma_aero_deer_deer_captive_null, sigma_aero_deer_human_wild = sigma_aero_deer_human_suburban,  sigma_aero_deer_human_capt = sigma_aero_deer_human_captive_null, sigma_dc_deer_deer = sigma_dc_deer_deer_null,alpha_immunity = alpha_immunity_null, epsilon_dc = epsilon_dc, gamma_recov = gamma_recov, I_human = rep(I_human_null, nsamples), boost = rep(0, nsamples))

4. Define initial SIR compartment sizes 

captive.inits.fall <- initial_compartments(S_wild_prop = 0, draws = nsamples)
captive.inits.steady <- initial_compartments(S_wild_prop = 0, draws = nsamples, steady = TRUE)

wild.inits.fall <- initial_compartments(S_captive_prop = 0, draws = nsamples)
wild.inits.steady <- initial_compartments(S_captive_prop = 0, draws = nsamples, steady = TRUE)

5. Solve the SIRS ODE equations and calculate daily compartment sizes 

times <- seq(0, 120, by = 1)
proj.ranch <- run(iter = nsamples, initial_compartments = captive.inits.fall, initial_compartments_steady = captive.inits.steady, params = ranch.params, times = times, name = "Outdoor ranch")

proj.intensive <- run(iter = nsamples, initial_compartments = captive.inits.fall, initial_compartments_steady = captive.inits.steady, params = intensive.params, times = times, name = "Intensive facility")

proj.rural <- run(iter = nsamples, initial_compartments = wild.inits.fall, initial_compartments_steady = wild.inits.steady ,params = rural.params, times = times, name = "Wild, rural")

proj.suburban <- run(iter = nsamples, initial_compartments = wild.inits.fall, initial_compartments_steady = wild.inits.steady, params = suburban.params, times = times, name = "Wild, suburban")

sirs_results_contexts <- rbind(proj.ranch,proj.intensive, proj.rural, proj.suburban)

6 and 7: Calculate risk of introduction and magnitude of spread given parameters, average daily prevalence, probability of persistence, and incidence proportion.

We are going to calculate several metrics for each iteration to summarize across contexts. First, we combine the projections for each context into one list object, and summarize the average prevalence (Infected) and incidence proportion (Cumulative_infections). We save this resulting summary dataframe for future use.

sirs_results_contexts %>%
  mutate(ode_df = map(ode_proj, as.data.frame)) %>%
  pull(ode_df) %>%
  list_rbind(names_to = "run_id") %>% 
  mutate(., Context = c(rep(unique(sirs_results_contexts$Context)[1],length(times)*nsamples),rep(unique(sirs_results_contexts$Context)[2],length(times)*nsamples), rep(unique(sirs_results_contexts$Context)[3],length(times)*nsamples),rep(unique(sirs_results_contexts$Context)[4],length(times)*nsamples))) %>% 
  group_by(Context, run_id) %>% 
   summarize(Wild = mean(I_wild), Captive = mean(I_captive), Prevalence = Wild + Captive, Cumulative_wild_infections = last(I_wild_cumulative), Cumulative_captive_infections = last(I_captive_cumulative), Cumulative_infections = sum(Cumulative_wild_infections,Cumulative_captive_infections)) %>% 
  mutate(., Context = factor(Context, levels = c("Outdoor ranch", "Intensive facility", "Wild, rural","Wild, suburban"))) %>% 
   select(., -Wild, -Captive, -Cumulative_wild_infections, -Cumulative_captive_infections) %>% 
  group_by(.,Context) %>% 
  arrange(., Context, run_id) -> Prev_cumulative_df

head(Prev_cumulative_df)
#> # A tibble: 6 × 4
#> # Groups:   Context [1]
#>   Context       run_id Prevalence Cumulative_infections
#>   <fct>          <int>      <dbl>                 <dbl>
#> 1 Outdoor ranch      1 0.00000782              0.000165
#> 2 Outdoor ranch      2 0.0892                  1.85    
#> 3 Outdoor ranch      3 0.0647                  1.34    
#> 4 Outdoor ranch      4 0.0000217               0.000459
#> 5 Outdoor ranch      5 0.000142                0.00300 
#> 6 Outdoor ranch      6 0.0661                  1.35

Next, we will join information from the steady-state equilibrium, using a threshold value of 0.1% prevalence at equilibrium for persistence.


persist.threshold <- 0.001

sirs_results_contexts %>%
  mutate(steady_sir = map(steady_state, "y"),
         steady_sir = map(steady_sir, as_tibble_row)) %>%
  pull(steady_sir) %>%
  list_rbind(names_to = "run_id") %>%
  mutate(., Context = c(rep(unique(sirs_results_contexts$Context)[1],nsamples),rep(unique(sirs_results_contexts$Context)[2],nsamples), rep(unique(sirs_results_contexts$Context)[3],nsamples),rep(unique(sirs_results_contexts$Context)[4],nsamples))) %>%
  mutate(., Context = factor(Context, levels = c("Outdoor ranch", "Intensive facility", "Wild, rural","Wild, suburban"))) %>% 
  mutate(., Persist.I_wild = I_wild > persist.threshold, Persist.I_captive = I_captive > persist.threshold, Persist = as.logical(Persist.I_wild + Persist.I_captive)) %>%
  select(., run_id,Context, Persist) %>%
  arrange(., Context, run_id) %>% 
   merge(Prev_cumulative_df,.) -> Prev_cumulative_persist_df

Next, we step through each context’s parameters to estimate strength of introduction (quantified as force-of-infection) and spread (quantified as basic reproductive number \(R_0\)). We then combine with the prevalence and persistence summary above to create the dataset that we will visualized in the visualization vignette.

ranch.df <- list_cbind(map(ranch.params, as_data_frame))
colnames(ranch.df) <-  names(ranch.params)
ranch.df %>% 
  mutate(., run_id = min(Prev_cumulative_df[which(Prev_cumulative_df$Context=="Outdoor ranch"),"run_id"])+0:(nsamples-1),r0 = unlist((beta_aero_ww+beta_aero_cc+beta_dc_ww+beta_dc_cc)/gamma_recov),
         FOI = unlist((beta_aero_hw+beta_aero_hc)*I_human),
         Context = "Outdoor ranch", Setting = "Captive") %>% 
  select(., run_id, Context,  Setting, r0, FOI) -> ranch.r0.FOI

intensive.df <- list_cbind(map(intensive.params, as_data_frame))
colnames(intensive.df) <-  names(intensive.params)
intensive.df %>% 
  mutate(., run_id = min(Prev_cumulative_df[which(Prev_cumulative_df$Context=="Intensive facility"),"run_id"])+0:(nsamples-1),r0 = unlist((beta_aero_ww+beta_aero_cc+beta_dc_ww+beta_dc_cc)/gamma_recov),
         FOI = unlist((beta_aero_hw+beta_aero_hc)*I_human),
         Context = "Intensive facility", Setting = "Captive") %>% 
  select(., run_id, Context,  Setting, r0, FOI) -> intensive.r0.FOI

rural.df <- list_cbind(map(rural.params, as_data_frame))
colnames(rural.df) <-  names(rural.params)
rural.df %>% 
  mutate(., run_id = min(Prev_cumulative_df[which(Prev_cumulative_df$Context=="Wild, rural"),"run_id"])+0:(nsamples-1),r0 = unlist((beta_aero_ww+beta_aero_cc+beta_dc_ww+beta_dc_cc)/gamma_recov),
         FOI = unlist((beta_aero_hw+beta_aero_hc)*I_human),
         Context = "Wild, rural", Setting = "Wild") %>% 
  select(., run_id, Context, Setting,  r0, FOI) -> rural.r0.FOI

suburban.df <- list_cbind(map(suburban.params, as_data_frame))
colnames(suburban.df) <-  names(suburban.params)
suburban.df %>% 
  mutate(., run_id = min(Prev_cumulative_df[which(Prev_cumulative_df$Context=="Wild, suburban"),"run_id"])+0:(nsamples-1),r0 = unlist((beta_aero_ww+beta_aero_cc+beta_dc_ww+beta_dc_cc)/gamma_recov),
         FOI = unlist((beta_aero_hw+beta_aero_hc)*I_human),
         Context = "Wild, suburban", Setting = "Wild") %>% 
  select(., run_id, Context, Setting, r0, FOI) -> suburban.r0.FOI

r0.FOI <- rbind(ranch.r0.FOI,intensive.r0.FOI, rural.r0.FOI,suburban.r0.FOI)

#Merge with average prevalence (step 1) and persistence (step 2)
merge(r0.FOI,Prev_cumulative_persist_df[,c("run_id","Prevalence", "Persist", "Cumulative_infections")], by = "run_id", all.x = TRUE) %>% 
  mutate(Context = factor(Context,  levels = c("Outdoor ranch", "Intensive facility", "Wild, rural","Wild, suburban"))) -> df

We will develop plots and tables of this results data frame. Here is a simple summary to end this vignette, with median spread (R0), probability introduction from humans (FOI, log10-transformed), average prevalence for the 120-day fall projection (Prevalence), and if there are any infections at the end of the 120-day fall projection (Persistence).

df %>% 
  group_by(., Context) %>% 
  reframe(., R0.Median = round(quantile(r0,probs = c(0.5)),2), 
          R0.80CI = paste0(round(quantile(r0,probs = c(0.1)),2)," - ",round(quantile(r0,probs = c(0.9)),2)), 
          FOI.Median = round((quantile(FOI*100,probs = c(0.5))),2), 
          FOI.80CI = paste0(round((quantile(FOI*100,probs = c(0.1))),2)," - ",round((quantile(FOI*100,probs = c(0.9))),2)), 
          Percent.Prob.1in1000.Infected.Median = round((1-(exp(-last(times)*quantile(FOI, probs = c(0.5)))^1000))*100, 1),
          Percent.Prob.1in1000.80CI = paste0(round((1-(exp(-last(times)*quantile(FOI, probs = c(0.1)))^1000))*100, 1)," - ",round((1-(exp(-last(times)*quantile(FOI, probs = c(0.9)))^1000))*100, 1)), 
          Prevalence.Median = round(quantile(Prevalence, probs = c(0.5)),6), 
          Prevalence.80CI = paste0(round(quantile(Prevalence,probs = c(0.1)),6)," - ",round(quantile(Prevalence,probs = c(0.9)),6)), 
          Persistence.Mean = round(sum(Persist)/nsamples,2), 
          Persistence.80CI = paste0(round(binom.confint(x = sum(Persist),n = nsamples, conf.level = 0.95, methods = "exact")$lower,2)," - ", round(binom.confint(x = sum(Persist),n = nsamples, conf.level = 0.95, methods = "exact")$upper,2)), 
          Incidence.Proportion.Median  = median(Cumulative_infections),
          Incidence.Proportion.80CI  = paste0(round(quantile(Cumulative_infections,probs = c(0.1)),2)," - ",round(quantile(Cumulative_infections,probs = c(0.9)),2)), ) -> results

results %>%
  kbl() %>%
  kable_paper("hover", full_width = F)
Context R0.Median R0.80CI FOI.Median FOI.80CI Percent.Prob.1in1000.Infected.Median Percent.Prob.1in1000.80CI Prevalence.Median Prevalence.80CI Persistence.Mean Persistence.80CI Incidence.Proportion.Median Incidence.Proportion.80CI
Outdoor ranch 2.35 0.36 - 17.88 0 0 - 0 48.9 7.7 - 97.1 0.049872 2.2e-05 - 0.094106 0.68 0.58 - 0.77 1.0108415 0 - 1.97
Intensive facility 8.15 1.04 - 112.35 0 0 - 0.01 82.4 18.9 - 100 0.073934 0.001248 - 0.10853 0.90 0.82 - 0.95 1.5355313 0.03 - 2.23
Wild, rural 1.17 0.17 - 9.31 0 0 - 0 1.2 0 - 11.2 0.000053 0 - 0.07232 0.53 0.43 - 0.63 0.0011257 0 - 1.49
Wild, suburban 1.17 0.17 - 9.31 0 0 - 0 16.2 4 - 62.9 0.000405 5e-06 - 0.074064 0.55 0.45 - 0.65 0.0094693 0 - 1.55

Next Steps

This vignette detailed the simulations of continuous exposure to humans infected with SARS-CoV-2. The next vignette runs a similar process, but simulates outbreaks caused by initial infections in a deer population, rather than continual introduction from humans. Results from both this and the next vignette are collated into another vignette to produce the figures similar to those in Rosenblatt et al. In Prep..

Click here to return to the vignette demonstrating an example use of whitetailedSIRS.

Literature Cited

  • Rosenblatt, E., J.D. Cook, G.V. Direnzo, E.H.C. Grant, F. Arce, K. Pepin, F.J. Rudolph, M.C. Runge, S. Shriner, D. Walsh, B.A. Mosher. In Prep. Epidemiological modeling of SARS-CoV-2 in white-tailed deer (Odocoileus virginianus) reveals conditions for introduction and widespread transmission.