SARS-CoV-2 introduction and spread in connected systems

Introduction

In the previous vignettes, we have focused on simulating outbreaks for a particular scenario focused on wild or captive deer populations. In this vignette, we can test how connected systems (wild and captive separated by a fence line) differ in outbreak characteristics, compared to isolated scenarios. When these scenarios are connected, there may be differences in prevalence, incidence proportion, and persistence due to fence line interactions. This analysis corresponds with Objective 5 of Rosenblatt et al. In Prep.

The systems considered here are all combinations of wild and captive scenarios detailed in the vignette simulating outbreaks, including:

Captive deer in outdoor ranch facilitates connected to wild deer in a rural setting;
Captive deer in outdoor ranch facilitates connected to wild deer in a suburban setting;
Captive deer in intensive, predominantly indoor facilitates connected to wild deer in a rural setting;
Captive deer in intensive, predominantly indoor facilitates connected to wild deer in a suburban setting;

To do this, we will use the same framework used in previous vignettes:

Set the number of simulations
Draw a random sample of parameters used to estimate infection probabilities, contact rates, and other important parameters.
Derive parameters used the SIRS ODE equation set
Define initial SIR compartment sizes
Solve the SIRS ODE equations and calculate daily compartment sizes
Calculate average daily prevalence, probability of persistence, and incidence proportion.

Note that we do not calculate FOI or \(R_0\) in this vignette, as these metrics are specific to the scenarios, not the overall system.

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

set.seed(23)
nsamples <- 1000

2. Sample parameters for simulations

We load the dataset and remind R how many iterations were run in each context (nsamples). We set a seed for consistent results, identify human prevalence and duration of projection (days). We then set parameters that are used consistently across all simulations in the vignette.

We bring in expert-elicited parameters used in previous vignettes:

#Bring in estimates from expert elicitation
elicitation_data <- draw_elicitation_samples(nsamples = nsamples)

We define common parameter values for all systems:

#Parameters commonly used across projections
nWild <- rpois(nsamples,1000) #Abundance
A_w <- 100 #Area
habitat <- "med" #Habitat classification

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

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

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

t_contact_deer_deer_null <- get_EE_param_vals(data = elicitation_data, my_param = "Deer Proximity Duration (minutes)") #Estimate duration of deer-deer proximity event...

We also fix parameters key to the SIRS ODE, including recovery rate to infection, rate of immunity loss, and human prevalence.

#Fixed parameters
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

System parameters

Below, we combine a wild scenario and a captive scenario, with a fence line separating these populations. For each combination, we define new parameters that have not been defined previously. As we work through these combinations, fewer parameters are added as they exist in earlier combinations.

Captive deer in outdoor ranch facilties connected to wild, rural deer

#Ranch side of the fence...

#Deer-to-deer proximity rate in ranch context, mimicking wild proximity rates with the influence of baiting (events per day).
omega_cc_ranch <- calc_contact_rate(nsamples = nsamples, type_contact = habitat, N_w = nWild)*(get_EE_param_vals(data = elicitation_data, my_param = "Proximity rate with baiting (17 events without baiting)")/17) 

#Probability of infection given that duration of proximity event.
sigma_aero_deer_deer_null <- calc_sigma_aero(C_nu = C_nu_deer,
                                        t_contact = t_contact_deer_deer_null / 60,
                                        r = r_deer, nsamples = nsamples) 

#Estimate duration of human-deer proximity event in ranch facility context.
t_contact_deer_human_ranch <- get_EE_param_vals(data = elicitation_data, my_param = "Deer-Human Proximity Duration, Captive (minutes)") 

#Human-to-deer proximity rate in ranch context, mimicking suburban proximitty rates (events per day).
omega_hc_ranch <- get_EE_param_vals(data = elicitation_data, my_param = "Deer-Human Proximity Rate, Suburban (per 120 days)") /120 

#Infection probability given the duration of a human-deer proximity event.
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)

#Deer-to-deer proximity rate along fenceline (events per day).
omega_cw <- 0.00072 / get_EE_param_vals(data = elicitation_data, my_param = "Direct Contact Probability") 

#On the wild side of the fence

#Deer-to-deer proximity rate
omega_ww <- calc_contact_rate(nsamples = nsamples, type_contact = habitat, N_w = nWild)

#Duration of human-deer proximity event
t_contact_deer_human_rural <- get_EE_param_vals(data = elicitation_data, my_param = "Deer-Human Proximity Duration, Rural (minutes)")

#Human-to-deer proximity rate
omega_hw_rural <- get_EE_param_vals(data = elicitation_data, my_param = "Deer-Human Proximity Rate, Rural (per 120 days)") /120

#Infection probability given the duration of a human-deer proximity event.
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))

Captive deer in outdoor ranch facilties connected to wild, suburban deer

Below we modify what’s happening on the wild side of the fence to a suburban context, and can keep the captive side parameterized to describe the outdoor ranch facility scenario.

#Duration of human-deer proximity event
t_contact_deer_human_suburban <- get_EE_param_vals(data = elicitation_data, my_param = "Deer-Human Proximity Duration, Suburban (minutes)")

#Human-to-deer proximity rate
omega_hw_suburban <- get_EE_param_vals(data = elicitation_data, my_param = "Deer-Human Proximity Rate, Suburban (per 120 days)") /120

#Infection probability given the duration of a human-deer proximity event.
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))

Captive deer in intensive facilties connected to wild, rural deer

Since we have defined parameter sets for both wild scenarios, below we modify what’s happening on the captive side of the fence.

#Deer-to-deer proximity rate in ranch context, mimicking wild proximity rates with the influence of baiting (events per day).
omega_cc_intensive <- get_EE_param_vals(data = elicitation_data, my_param = "Deer-Deer Proximity Rate, Captive (per day)")

#Probability of infection given that duration of proximity event.
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))

#Human-to-deer proximity rate in intensive facilities.
omega_hc_intensive <- get_EE_param_vals(data = elicitation_data, my_param = "Deer-Human Proximity Rate, Suburban (per 120 days)") /120

#Duration of human-deer proximity event
t_contact_deer_human_intensive <- get_EE_param_vals(data = elicitation_data, my_param = "Deer-Human Proximity Duration, Captive (minutes)") 

#Infection probability given the duration of a human-deer proximity event.
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))

Captive deer in intensive facilties connected to wild, suburban deer

The final combination of intensive captive facilities and wild, suburban scenarios does not require any additional parameters, so we can go ahead and derive parameters for each system.

3. Derive parameters used the SIRS ODE equation set

ranch_rural.params <- alternative(alpha_immunity = alpha_immunity_null,
                            omega_ww = omega_ww, omega_cw = omega_cw, omega_cc = omega_cc_ranch,
                            omega_hw = omega_hw_rural, omega_hc = omega_hc_ranch,
                            sigma_aero_deer_deer_wild = sigma_aero_deer_deer_null, sigma_aero_deer_deer_captive = sigma_aero_deer_deer_null, sigma_aero_deer_human_wild = sigma_aero_deer_human_rural, 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))

ranch_suburban.params <- alternative(alpha_immunity = alpha_immunity_null,
                            omega_ww = omega_ww, omega_cw = omega_cw, omega_cc = omega_cc_ranch,
                            omega_hw = omega_hw_suburban, omega_hc = omega_hc_ranch,
                            sigma_aero_deer_deer_wild = sigma_aero_deer_deer_null, sigma_aero_deer_deer_captive = sigma_aero_deer_deer_null, sigma_aero_deer_human_wild = sigma_aero_deer_human_suburban, 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_rural.params <- alternative(alpha_immunity = alpha_immunity_null,
                            omega_ww = omega_ww, omega_cw = omega_cw, omega_cc = omega_cc_intensive,
                            omega_hw = omega_hw_rural, omega_hc = omega_hc_intensive,
                            sigma_aero_deer_deer_wild = sigma_aero_deer_deer_null, sigma_aero_deer_deer_captive = sigma_aero_deer_deer_intensive, sigma_aero_deer_human_wild = sigma_aero_deer_human_rural, sigma_aero_deer_human_capt = sigma_aero_deer_human_intensive, 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_suburban.params <- alternative(alpha_immunity = alpha_immunity_null,
                            omega_ww = omega_ww, omega_cw = omega_cw, omega_cc = omega_cc_intensive,
                            omega_hw = omega_hw_suburban, omega_hc = omega_hc_intensive,
                            sigma_aero_deer_deer_wild = sigma_aero_deer_deer_null, sigma_aero_deer_deer_captive = sigma_aero_deer_deer_intensive, sigma_aero_deer_human_wild = sigma_aero_deer_human_suburban, sigma_aero_deer_human_capt = sigma_aero_deer_human_intensive, 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))

4. Define initial SIR compartment sizes

Next, we define initial compartment sizes that will be used in all systems.

inits.fall <- initial_compartments(draws = nsamples) #Calculate initial distributions of animals across SIR compartments
inits.steady <- initial_compartments(draws = nsamples, steady = TRUE) #Calculate initial distributions of animals across SIR compartments, excluding solving cumulative cases.

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

We now define the duration of the simulation, and solve the SIRS ODE equations based on our parameter estimates.

times <- seq(0, 120, by = 1)

proj_ranch_rural <- run(iter = nsamples, initial_compartments = inits.fall, initial_compartments_steady = inits.steady, params = ranch_rural.params, times = times, name = "Outdoor ranch and rural system")

proj_ranch_suburban <- run(iter = nsamples, initial_compartments = inits.fall, initial_compartments_steady = inits.steady, params = ranch_suburban.params, times = times, name = "Outdoor ranch and suburban system")

proj_intensive_rural <- run(iter = nsamples, initial_compartments = inits.fall, initial_compartments_steady = inits.steady, params = intensive_rural.params, times = times, name = "Intensive facility and rural system")

proj_intensive_suburban <- run(iter = nsamples, initial_compartments = inits.fall, initial_compartments_steady = inits.steady, params = intensive_suburban.params, times = times, name = "Intensive facility and suburban system")

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

As we did in whitetailedSIRS::SIRS_analysis_by_context, we compute the average prevalence, cumulative infections, and persistence for each system (captive:wild scenario combination).

persist.threshold <- 0.001 #Set prevalence threshold

#Ranch and rural system
proj_ranch_rural %>%
  mutate(ode_df = map(ode_proj, as.data.frame)) %>%
  pull(ode_df) %>%
  list_rbind(names_to = "run_id") %>% 
  mutate(., Context = c(rep(unique(proj_ranch_rural$Context)[1],length(times)*nsamples))) %>% 
  group_by(Context, run_id) %>% 
  summarize(Wild = mean(I_wild), Captive = mean(I_captive), Cumulative_wild_infections = last(I_wild_cumulative), Cumulative_captive_infections = last(I_captive_cumulative)) %>% 
   mutate(., Context = factor(Context, levels = c("Outdoor ranch and rural system"))) %>% 
  arrange(., run_id) -> Prev_cumulative_ranch_rural_system_df
#> `summarise()` has grouped output by 'Context'. You can override using the
#> `.groups` argument.

proj_ranch_rural %>%
  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(proj_ranch_rural$Context)[1],nsamples))) %>%
  mutate(., Persist_wild = I_wild > persist.threshold, Persist_captive = I_captive > persist.threshold) %>%
   mutate(., Context = factor(Context, levels = c("Outdoor ranch and rural system"))) %>% 
  select(., run_id, Context, Persist_wild, Persist_captive) %>% 
   merge(Prev_cumulative_ranch_rural_system_df,.) %>% 
      arrange(., run_id)-> Prev_cumulative_persist_ranch_rural_system_df

#Ranch and suburban system
proj_ranch_suburban %>%
  mutate(ode_df = map(ode_proj, as.data.frame)) %>%
  pull(ode_df) %>%
  list_rbind(names_to = "run_id") %>% 
  mutate(., Context = c(rep(unique(proj_ranch_suburban$Context)[1],length(times)*nsamples))) %>% 
  group_by(Context, run_id) %>% 
  summarize(Wild = mean(I_wild), Captive = mean(I_captive), Cumulative_wild_infections = last(I_wild_cumulative), Cumulative_captive_infections = last(I_captive_cumulative)) %>% 
   mutate(., Context = factor(Context, levels = c("Outdoor ranch and suburban system"))) %>% 
  arrange(., run_id) -> Prev_cumulative_ranch_suburban_system_df
#> `summarise()` has grouped output by 'Context'. You can override using the
#> `.groups` argument.

proj_ranch_suburban %>%
  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(proj_ranch_suburban$Context)[1],nsamples))) %>%
  mutate(., Persist_wild = I_wild > persist.threshold, Persist_captive = I_captive > persist.threshold) %>%
   mutate(., Context = factor(Context, levels = c("Outdoor ranch and suburban system"))) %>% 
  select(., run_id, Context, Persist_wild, Persist_captive) %>% 
   merge(Prev_cumulative_ranch_suburban_system_df,.) %>% 
      arrange(., run_id)-> Prev_cumulative_persist_ranch_suburban_system_df

#Intensive facility and rural system
proj_intensive_rural %>%
  mutate(ode_df = map(ode_proj, as.data.frame)) %>%
  pull(ode_df) %>%
  list_rbind(names_to = "run_id") %>% 
  mutate(., Context = c(rep(unique(proj_intensive_rural$Context)[1],length(times)*nsamples))) %>% 
  group_by(Context, run_id) %>% 
  summarize(Wild = mean(I_wild), Captive = mean(I_captive), Cumulative_wild_infections = last(I_wild_cumulative), Cumulative_captive_infections = last(I_captive_cumulative)) %>% 
   mutate(., Context = factor(Context, levels = c("Intensive facility and rural system"))) %>% 
  arrange(., run_id) -> Prev_cumulative_intensive_rural_system_df
#> `summarise()` has grouped output by 'Context'. You can override using the
#> `.groups` argument.

proj_intensive_rural %>%
  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(proj_intensive_rural$Context)[1],nsamples))) %>%
  mutate(., Persist_wild = I_wild > persist.threshold, Persist_captive = I_captive > persist.threshold) %>%
   mutate(., Context = factor(Context, levels = c("Intensive facility and rural system"))) %>% 
  select(., run_id, Context, Persist_wild, Persist_captive) %>% 
   merge(Prev_cumulative_intensive_rural_system_df,.) %>% 
      arrange(., run_id)-> Prev_cumulative_persist_intensive_rural_system_df

#Ranch and suburban system
proj_intensive_suburban %>%
  mutate(ode_df = map(ode_proj, as.data.frame)) %>%
  pull(ode_df) %>%
  list_rbind(names_to = "run_id") %>% 
  mutate(., Context = c(rep(unique(proj_intensive_suburban$Context)[1],length(times)*nsamples))) %>% 
  group_by(Context, run_id) %>% 
  summarize(Wild = mean(I_wild), Captive = mean(I_captive), Cumulative_wild_infections = last(I_wild_cumulative), Cumulative_captive_infections = last(I_captive_cumulative)) %>% 
   mutate(., Context = factor(Context, levels = c("Intensive facility and suburban system"))) %>% 
  arrange(., run_id) -> Prev_cumulative_intensive_suburban_system_df
#> `summarise()` has grouped output by 'Context'. You can override using the
#> `.groups` argument.

proj_intensive_suburban %>%
  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(proj_intensive_suburban$Context)[1],nsamples))) %>%
  mutate(., Persist_wild = I_wild > persist.threshold, Persist_captive = I_captive > persist.threshold) %>%
   mutate(., Context = factor(Context, levels = c("Intensive facility and suburban system"))) %>% 
  select(., run_id, Context, Persist_wild, Persist_captive) %>% 
   merge(Prev_cumulative_intensive_suburban_system_df,.) %>% 
      arrange(., run_id)-> Prev_cumulative_persist_intensive_suburban_system_df

7. Compare outbreak metrics to simulations with isolated scenarios

Below, we take the results from whitetailedSIRS::SIRS_analysis_by_contexts and calculate differences in prevalence, cumulative infections, and persistence when fence line interactions can occur. First, we read in the scenario results and format to match the combined systems.

#Read in scenario results for comparison with system results
df <- whitetailedSIRS::scenario_results

df %>% 
   mutate(., run_id = case_when(Context == "Outdoor ranch" ~ run_id-0,
      Context == "Intensive facility" ~ run_id-(nsamples*1),
      Context == "Wild, rural" ~ run_id-(nsamples*2),
      Context == "Wild, suburban" ~ run_id-(nsamples*3),)) -> df_to_match

Next, we calculate differences in prevalence, persistence, and incidence proportions between scenarios and broader systems. Below, we go system by system:

#Ranch, rural
Prev_cumulative_persist_ranch_rural_system_df %>% 
   merge(., df_to_match[which(df_to_match$Context== "Outdoor ranch"),c("run_id","Prevalence", "Persist", "Cumulative_infections")], by = "run_id") %>% 
   rename(., "System_wild_prevalence" = "Wild", "System_captive_prevalence" = "Captive", "System_wild_cumulative" = "Cumulative_wild_infections", "System_captive_cumulative" = "Cumulative_captive_infections","System_wild_persistence" = "Persist_wild", "System_captive_persistence" = "Persist_captive", "Context_captive_prevalence" = "Prevalence", "Context_captive_persistence" = "Persist", "Context_captive_cumulative" = "Cumulative_infections") %>% 
   merge(., df_to_match[which(df_to_match$Context== "Wild, rural"),c("run_id","Prevalence", "Persist", "Cumulative_infections")], by = "run_id") %>% 
   rename(., "Context_wild_prevalence" = "Prevalence", "Context_wild_persistence" = "Persist", "Context_wild_cumulative" = "Cumulative_infections") %>% 
   mutate(., Change_Prev_Wild = System_wild_prevalence - Context_wild_prevalence, 
          Change_Prev_Captive = System_captive_prevalence - Context_captive_prevalence, 
          Proportional_Change_Prev_Wild = Change_Prev_Wild/Context_wild_prevalence, 
          Proportional_Change_Prev_Captive = Change_Prev_Captive/Context_captive_prevalence, 
          Change_Cumulative_Wild = System_wild_cumulative - Context_wild_cumulative, 
          Change_Cumulative_Captive = System_captive_cumulative - Context_captive_cumulative, 
          Change_Persist_Wild = System_wild_persistence - Context_wild_persistence, 
          Change_Persist_Captive = System_captive_persistence - Context_captive_persistence) %>% 
   select(., Context, Change_Prev_Wild, Change_Prev_Captive, Proportional_Change_Prev_Wild, Proportional_Change_Prev_Captive, Change_Cumulative_Wild, Change_Cumulative_Captive, Change_Persist_Wild, Change_Persist_Captive) %>% 
   rename(., System = Context) -> ranch_rural_results

#Ranch, suburban
Prev_cumulative_persist_ranch_suburban_system_df %>% 
   merge(., df_to_match[which(df_to_match$Context== "Outdoor ranch"),c("run_id","Prevalence", "Persist", "Cumulative_infections")], by = "run_id") %>% 
   rename(., "System_wild_prevalence" = "Wild", "System_captive_prevalence" = "Captive", "System_wild_cumulative" = "Cumulative_wild_infections", "System_captive_cumulative" = "Cumulative_captive_infections","System_wild_persistence" = "Persist_wild", "System_captive_persistence" = "Persist_captive", "Context_captive_prevalence" = "Prevalence", "Context_captive_persistence" = "Persist", "Context_captive_cumulative" = "Cumulative_infections") %>% 
   merge(., df_to_match[which(df_to_match$Context== "Wild, suburban"),c("run_id","Prevalence", "Persist", "Cumulative_infections")], by = "run_id") %>% 
   rename(., "Context_wild_prevalence" = "Prevalence", "Context_wild_persistence" = "Persist", "Context_wild_cumulative" = "Cumulative_infections") %>% 
   mutate(., Change_Prev_Wild = System_wild_prevalence - Context_wild_prevalence, 
          Change_Prev_Captive = System_captive_prevalence - Context_captive_prevalence, 
          Proportional_Change_Prev_Wild = Change_Prev_Wild/Context_wild_prevalence, 
          Proportional_Change_Prev_Captive = Change_Prev_Captive/Context_captive_prevalence, 
          Change_Cumulative_Wild = System_wild_cumulative - Context_wild_cumulative, 
          Change_Cumulative_Captive = System_captive_cumulative - Context_captive_cumulative, 
          Change_Persist_Wild = System_wild_persistence - Context_wild_persistence, 
          Change_Persist_Captive = System_captive_persistence - Context_captive_persistence) %>% 
   select(., Context, Change_Prev_Wild, Change_Prev_Captive, Proportional_Change_Prev_Wild, Proportional_Change_Prev_Captive, Change_Cumulative_Wild, Change_Cumulative_Captive, Change_Persist_Wild, Change_Persist_Captive) %>% 
   rename(., System = Context) -> ranch_suburban_results

#Intensive, rural
Prev_cumulative_persist_intensive_rural_system_df %>% 
   merge(., df_to_match[which(df_to_match$Context== "Intensive facility"),c("run_id","Prevalence", "Persist", "Cumulative_infections")], by = "run_id") %>% 
   rename(., "System_wild_prevalence" = "Wild", "System_captive_prevalence" = "Captive", "System_wild_cumulative" = "Cumulative_wild_infections", "System_captive_cumulative" = "Cumulative_captive_infections","System_wild_persistence" = "Persist_wild", "System_captive_persistence" = "Persist_captive", "Context_captive_prevalence" = "Prevalence", "Context_captive_persistence" = "Persist", "Context_captive_cumulative" = "Cumulative_infections") %>% 
   merge(., df_to_match[which(df_to_match$Context== "Wild, rural"),c("run_id","Prevalence", "Persist", "Cumulative_infections")], by = "run_id") %>% 
   rename(., "Context_wild_prevalence" = "Prevalence", "Context_wild_persistence" = "Persist", "Context_wild_cumulative" = "Cumulative_infections") %>% 
   mutate(., Change_Prev_Wild = System_wild_prevalence - Context_wild_prevalence, 
          Change_Prev_Captive = System_captive_prevalence - Context_captive_prevalence, 
          Proportional_Change_Prev_Wild = Change_Prev_Wild/Context_wild_prevalence, 
          Proportional_Change_Prev_Captive = Change_Prev_Captive/Context_captive_prevalence, 
          Change_Cumulative_Wild = System_wild_cumulative - Context_wild_cumulative, 
          Change_Cumulative_Captive = System_captive_cumulative - Context_captive_cumulative, 
          Change_Persist_Wild = System_wild_persistence - Context_wild_persistence, 
          Change_Persist_Captive = System_captive_persistence - Context_captive_persistence) %>% 
   select(., Context, Change_Prev_Wild, Change_Prev_Captive, Proportional_Change_Prev_Wild, Proportional_Change_Prev_Captive, Change_Cumulative_Wild, Change_Cumulative_Captive, Change_Persist_Wild, Change_Persist_Captive) %>% 
   rename(., System = Context) -> intensive_rural_results

#Intensive, suburban
Prev_cumulative_persist_intensive_suburban_system_df %>% 
   merge(., df_to_match[which(df_to_match$Context== "Intensive facility"),c("run_id","Prevalence", "Persist", "Cumulative_infections")], by = "run_id") %>% 
   rename(., "System_wild_prevalence" = "Wild", "System_captive_prevalence" = "Captive", "System_wild_cumulative" = "Cumulative_wild_infections", "System_captive_cumulative" = "Cumulative_captive_infections","System_wild_persistence" = "Persist_wild", "System_captive_persistence" = "Persist_captive", "Context_captive_prevalence" = "Prevalence", "Context_captive_persistence" = "Persist", "Context_captive_cumulative" = "Cumulative_infections") %>% 
   merge(., df_to_match[which(df_to_match$Context== "Wild, suburban"),c("run_id","Prevalence", "Persist", "Cumulative_infections")], by = "run_id") %>% 
   rename(., "Context_wild_prevalence" = "Prevalence", "Context_wild_persistence" = "Persist", "Context_wild_cumulative" = "Cumulative_infections") %>% 
   mutate(., Change_Prev_Wild = System_wild_prevalence - Context_wild_prevalence, 
          Change_Prev_Captive = System_captive_prevalence - Context_captive_prevalence, 
          Proportional_Change_Prev_Wild = Change_Prev_Wild/Context_wild_prevalence, 
          Proportional_Change_Prev_Captive = Change_Prev_Captive/Context_captive_prevalence, 
          Change_Cumulative_Wild = System_wild_cumulative - Context_wild_cumulative, 
          Change_Cumulative_Captive = System_captive_cumulative - Context_captive_cumulative, 
          Change_Persist_Wild = System_wild_persistence - Context_wild_persistence, 
          Change_Persist_Captive = System_captive_persistence - Context_captive_persistence) %>% 
   select(., Context, Change_Prev_Wild, Change_Prev_Captive, Proportional_Change_Prev_Wild, Proportional_Change_Prev_Captive, Change_Cumulative_Wild, Change_Cumulative_Captive, Change_Persist_Wild, Change_Persist_Captive) %>% 
   rename(., System = Context) -> intensive_suburban_results

We can then combine and summarize the differences in outbreak metrics. The resulting table follows the format of table 2 in Rosenblatt et al. In Prep.

sirs_results_systems <- rbind(ranch_suburban_results,intensive_suburban_results,ranch_rural_results,intensive_rural_results)

sirs_results_systems %>% 
   group_by(., System) %>% 
   summarise(., Wild_Prevalence_Change = paste0(round(median(Change_Prev_Wild)*100,4), " (", round(quantile(Change_Prev_Wild,probs = 0.1)*100,4), "-",round(quantile(Change_Prev_Wild,probs = 0.9)*100,4),")"), Captive_Prevalence_Change = paste0(round(median(Change_Prev_Captive)*100,4), " (", round(quantile(Change_Prev_Captive,probs = 0.1)*100,4), "-",round(quantile(Change_Prev_Captive,probs = 0.9)*100,4),")"), 
             Wild_Prevalence_Proportional_Change = paste0(round(median(Proportional_Change_Prev_Wild)*100,4), " (", round(quantile(Proportional_Change_Prev_Wild,probs = 0.1)*100,4), "-",round(quantile(Proportional_Change_Prev_Wild,probs = 0.9)*100,4),")"), Proportional_Captive_Prevalence_Change = paste0(round(median(Proportional_Change_Prev_Captive)*100,4), " (", round(quantile(Proportional_Change_Prev_Captive,probs = 0.1)*100,4), "-",round(quantile(Proportional_Change_Prev_Captive,probs = 0.9)*100,4),")"),
             Wild_Incidence_Change = paste0(round(median(Change_Cumulative_Wild)*100,4), " (", round(quantile(Change_Cumulative_Wild,probs = 0.1)*100,4), "-",round(quantile(Change_Cumulative_Wild,probs = 0.9)*100,4),")"), Captive_Incidence_Change = paste0(round(median(Change_Cumulative_Captive)*100,4), " (", round(quantile(Change_Cumulative_Captive,probs = 0.1)*100,4), "-",round(quantile(Change_Cumulative_Captive,probs = 0.9)*100,4),")"), 
             Wild_Persist_Change = paste0(round(binom.confint(sum(Change_Persist_Wild),nsamples,methods = "exact")$mean,2), " (", round(binom.confint(sum(Change_Persist_Wild),nsamples,methods = "exact", conf.level = 0.80)$lower,4), "-", round(binom.confint(sum(Change_Persist_Wild),nsamples,methods = "exact", conf.level = 0.80)$upper,4), ")"),
             Captive_Persist_Change = sum(Change_Persist_Captive)) ->fence_effect

fence_effect %>%
  kbl() %>%
  kable_paper("hover", full_width = T)

System	Wild_Prevalence_Change	Captive_Prevalence_Change	Wild_Prevalence_Proportional_Change	Proportional_Captive_Prevalence_Change	Wild_Incidence_Change	Captive_Incidence_Change	Wild_Persist_Change
Outdoor ranch and suburban system	0.0022 (0-0.1499)	0 (0-8e-04)	0.4618 (0.0143-103.8241)	0.0028 (3e-04-0.0294)	0.0441 (0-3.2064)	5e-04 (0-0.0163)	0 (1e-04-0.0039)
Intensive facility and suburban system	0.0107 (0-0.5219)	0 (0-1e-04)	11.2484 (0.1707-539.1921)	1e-04 (0-0.0025)	0.2065 (2e-04-10.6701)	0 (0-0.0017)	0.01 (0.0032-0.0105)
Outdoor ranch and rural system	0.0037 (0-0.5565)	0 (0-8e-04)	4.4704 (0.0849-3094.426)	9e-04 (0-0.019)	0.0813 (0-10.954)	1e-04 (0-0.0163)	0.01 (0.0103-0.0212)
Intensive facility and rural system	0.0141 (0-1.0539)	0 (0-1e-04)	122.2533 (0.7449-19925.6897)	0 (0-0.0013)	0.2784 (2e-04-20.4194)	0 (0-0.0014)	0.02 (0.0137-0.0258)

Next Steps

This vignette detailed the simulations of connected wild and captive scenarios, refered to as systems. We compared outbreaks in these systems to outbreaks in discrete scenarios. By doing so, we better understand the role of fence line interactions between wild and captive deer populations in the spread and establishment of SARS-CoV-2 in white-tailed deer. The next vignette runs a similar process, but runs through multiple parameter sets that mimic a range of management alternatives that could potentially mitigate the introduction, spread, and establishment of SARS-CoV-2 in deer.

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.

Analysis referenced in Rosenblatt et al. In Prep

Elias Rosenblatt