Evaluation of management alternatives across systems

Introduction:

In this vignette, we run outbreak simulations to estimate the effects of various management alternatives on the dynamics of a SARS-CoV-2 outbreak in wild and captive white-tailed deer. We considered captive:wild systems that are separated by a fence. These management alternatives may be implemented in both wild and captive populations, or in one scenario. We focus on 11 alternatives, ranging across agriculture, public health, and wildlife sectors. These alternatives are detailed in Cook et al. In Prep.. These alternatives (and parameter changes) include:

Local reductions in wild deer populations through agency culling or hunter harvest (Reduce wild deer densities by 10% and 25%)
Prohibit deer feeding and baiting to reduce congregations of wild deer (Remove inflation factor of proximity rates due to bait and supplemental feeding)
Pause/reduce research permits for work on wild deer or other susceptible wildlife/animal species (Reduce human-wild deer interaction)
Encourage limiting human-deer interactions including participating in conduct that attracts wildlife (Proportionally reduce human-deer proximity rates reflecting rural:suburban rates)
Encourage (and train on the proper) use of enhanced personal protective equipment during interactions with captive and wild deer (Integrate effects of masks in aerosol transmission from humans to deer)
Animal vaccination requirement (All captive deer start simulations as recovered; All captive deer start simulations as recovered and are boosted when immunity wanes)
Upgrade fences to double fences that prevent transmission between captive and wild deer (Transmission between captive and wild deer is set to zero)
Increase ventilation in indoor, intensive facilities (Increase air exchange rate)

We detail the alternatives and the systems where they can be implemented:

ActionTable <- data.frame(Alternative = c("Reduce wild density by 10%", "Reduce wild density by 25%", "Eliminate baiting", "Pause Research", "Restrict human interactions with deer","PPE, Wild" , "PPE, Captive", "PPE, Both", "Vaccinate captive deer", "Vaccinate and boost captive deer", "Double-fence captive facilities", "Increase air exchange rate in intensive facilities"), Sector = c("Wildlife", "Wildlife", "Wildlife", "Wildlife", "Public Health","Public Health","Public Health","Public Health", "Agriculture", "Agriculture","Agriculture", "Agriculture"), Ranch_Rural_System = c(TRUE,TRUE,TRUE,TRUE,TRUE,TRUE,TRUE,TRUE,TRUE,TRUE,TRUE,FALSE), Ranch_Suburban_System = c(TRUE,TRUE,TRUE,TRUE,TRUE,TRUE,TRUE,TRUE,TRUE,TRUE,TRUE,FALSE), Intensive_Facility_Rural_System = c(TRUE,TRUE,TRUE,TRUE,TRUE,TRUE,TRUE,TRUE,TRUE,TRUE,TRUE,TRUE), Intensive_Facility_Suburban_System = c(TRUE,TRUE,TRUE,TRUE,TRUE,TRUE,TRUE,TRUE,TRUE,TRUE,TRUE,TRUE))

ActionTable %>%
  kbl() %>%
  kable_paper("hover", full_width = F)

Alternative	Sector	Ranch_Rural_System	Ranch_Suburban_System	Intensive_Facility_Rural_System	Intensive_Facility_Suburban_System
Reduce wild density by 10%	Wildlife	TRUE	TRUE	TRUE	TRUE
Reduce wild density by 25%	Wildlife	TRUE	TRUE	TRUE	TRUE
Eliminate baiting	Wildlife	TRUE	TRUE	TRUE	TRUE
Pause Research	Wildlife	TRUE	TRUE	TRUE	TRUE
Restrict human interactions with deer	Public Health	TRUE	TRUE	TRUE	TRUE
PPE, Wild	Public Health	TRUE	TRUE	TRUE	TRUE
PPE, Captive	Public Health	TRUE	TRUE	TRUE	TRUE
PPE, Both	Public Health	TRUE	TRUE	TRUE	TRUE
Vaccinate captive deer	Agriculture	TRUE	TRUE	TRUE	TRUE
Vaccinate and boost captive deer	Agriculture	TRUE	TRUE	TRUE	TRUE
Double-fence captive facilities	Agriculture	TRUE	TRUE	TRUE	TRUE
Increase air exchange rate in intensive facilities	Agriculture	FALSE	FALSE	TRUE	TRUE

To evaluate the effects of these interventions, we follow the same analytical framework as in other vignettes. Because of we are looking at multiple systems, our work flow will differ a little in this vignette. Steps include:

Set number of simulations to run
Sample parameters for simulations
Establish systems and test alternatives
Combine and summarize results across alternatives and systems
Impact of cross-sector collaborations (multiple alternatives implemented in concert).

1. Set number of simulations to run

Like in other vignettes, we set a seed and sample size. We use a small sample size for ease of vignette loading. We also define the length of the simulation here.

2. Sample parameters for simulations

Common Parameters

We bring in elicitation data for various parameters included in these simulations.

We can then define a number of parameters that will be used in every system in which we want to test alternatives.

Scenario parameters

Here, we define parameter values used for the baseline projection for each system:

Captive deer in outdoor ranch facilities

#Ranch Parameters
#Deer-to-deer proximity rate in ranch context, mimicing 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) 

#Human-to-deer proximity rate in ranch context, mimicing 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 

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

#Estimate probability of infection from deer-deer proximity given that 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) 

#Calculate 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)

Captive deer in intensive facilities

#Deer-to-deer proximity rate in intensive context (events per day).
omega_cc_intensive <- get_EE_param_vals(data = elicitation_data, my_param = "Deer-Deer Proximity Rate, Captive (per day)") 

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

#Estimate 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)")

#Estimate probability of infection from deer-deer proximity 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))

#Calculate 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))

Wild deer in a suburban setting

#Deer-to-deer proximity rate (events per day).
omega_ww_suburban <- 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)

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

#Estimate 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)")

#Estimate probability of infection from deer-deer proximity given that duration of proximity event.
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))

#Calculate 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))

Wild deer in a rural setting

#Deer-to-deer proximity rate (events per day).
omega_ww_rural <- 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)

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

#Estimate 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)")

#Estimate probability of infection from deer-deer proximity given that duration of proximity event.
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))

#Calculate 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))

3.Establish Systems and Test Alternatives

Below we define a captive-wild system, simulate SARS-CoV-2 outbreaks with no interventions, and then apply the 12 candidate interventions by altering a relevant parameter. All 13 simulations (1 no action + 12 interventions) are combined into one data frame at the end of each code chunk.

Combined system of outdoor ranch facilities and wild deer in a suburban setting

ranch_suburban.inits.fall <- initial_compartments(draws = nsamples)
ranch_suburban.inits.steady <- initial_compartments(draws = nsamples, steady = TRUE)


ranch_suburban.params <- alternative(alpha_immunity = alpha_immunity_null,
                            omega_ww = omega_ww_suburban, 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_suburban, sigma_aero_deer_deer_captive = sigma_aero_deer_deer_ranch, 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))

proj_ranch_suburban_null <- run(iter = nsamples, initial_compartments = ranch_suburban.inits.fall, initial_compartments_steady = ranch_suburban.inits.steady, params = ranch_suburban.params, times = times, name = "Baseline")

#Alternatives
#1. 10% density reduction
omega_ww_alt1_density.red.10 <- whitetailedSIRS::calc_contact_rate(nsamples = nsamples, type_contact = habitat, N_w = nWild*0.9)*(get_EE_param_vals(data = elicitation_data, my_param = "Proximity rate with baiting (17 events without baiting)")/17) 

ranch_suburban.params_alt1 <- alternative(alpha_immunity = alpha_immunity_null,
                            omega_ww = omega_ww_alt1_density.red.10, 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_suburban, sigma_aero_deer_deer_captive = sigma_aero_deer_deer_ranch, 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))

proj_ranch_suburban_alt1 <- run(iter = nsamples, initial_compartments = ranch_suburban.inits.fall, initial_compartments_steady = ranch_suburban.inits.steady, params = ranch_suburban.params_alt1, times = times, name = "Reduce wild density by 10%")

#2. 25% density reduction
omega_ww_alt2_density.red.25 <- whitetailedSIRS::calc_contact_rate(nsamples = nsamples, type_contact = habitat, N_w = nWild*0.75)*(get_EE_param_vals(data = elicitation_data, my_param = "Proximity rate with baiting (17 events without baiting)")/17)

ranch_suburban.params_alt2 <- alternative(alpha_immunity = alpha_immunity_null,
                            omega_ww = omega_ww_alt2_density.red.25, 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_suburban, sigma_aero_deer_deer_captive = sigma_aero_deer_deer_ranch, 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))

proj_ranch_suburban_alt2 <- run(iter = nsamples, initial_compartments = ranch_suburban.inits.fall, initial_compartments_steady = ranch_suburban.inits.steady, params = ranch_suburban.params_alt2, times = times, name = "Reduce wild density by 25%")

#3. Eliminate baiting
omega_ww_alt3_no.baiting <- whitetailedSIRS::calc_contact_rate(nsamples = nsamples, type_contact = habitat, N_w = nWild*0.75)

ranch_suburban.params_alt3 <- alternative(alpha_immunity = alpha_immunity_null,
                            omega_ww = omega_ww_alt3_no.baiting, 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_suburban, sigma_aero_deer_deer_captive = sigma_aero_deer_deer_ranch, 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))

proj_ranch_suburban_alt3 <- run(iter = nsamples, initial_compartments = ranch_suburban.inits.fall, initial_compartments_steady = ranch_suburban.inits.steady, params = ranch_suburban.params_alt3, times = times, name = "Eliminate baiting")

#4. Pause research (assuming random 1-50% of human-wild deer interaction)
omega_hw_alt4_pause.research <- omega_hw_suburban*runif(nsamples, min = 0.5, max = 0.99) #reduce human deer proximity rate in the wild up to 50%

ranch_suburban.params_alt4 <- alternative(alpha_immunity = alpha_immunity_null,
                            omega_ww = omega_ww_suburban, omega_cw = omega_cw, omega_cc = omega_cc_ranch,
                            omega_hw = omega_hw_alt4_pause.research, omega_hc = omega_hc_ranch,
                            sigma_aero_deer_deer_wild = sigma_aero_deer_deer_suburban, sigma_aero_deer_deer_captive = sigma_aero_deer_deer_ranch, 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))

proj_ranch_suburban_alt4 <- run(iter = nsamples, initial_compartments = ranch_suburban.inits.fall, initial_compartments_steady = ranch_suburban.inits.steady, params = ranch_suburban.params_alt4, times = times, name = "Pause Research")

#5. Require PPE for humans in wild settings
sigma_aero_deer_human_wild_alt5_PPE <- calc_sigma_aero(ER = 0.53, C_nu = C_nu_human, t_contact = t_contact_deer_human_suburban/60, r = r_deer, V_d = 0.009*runif(nsamples, min = (1-0.994), max = (1-0.509)), nsamples = nsamples)

ranch_suburban.params_alt5 <- alternative(alpha_immunity = alpha_immunity_null,
                            omega_ww = omega_ww_suburban, 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_suburban, sigma_aero_deer_deer_captive = sigma_aero_deer_deer_ranch, sigma_aero_deer_human_wild = sigma_aero_deer_human_wild_alt5_PPE, 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))

proj_ranch_suburban_alt5 <- run(iter = nsamples, initial_compartments = ranch_suburban.inits.fall, initial_compartments_steady = ranch_suburban.inits.steady, params = ranch_suburban.params_alt5, times = times, name = "PPE, Wild")

#6. Require PPE for humans in captive settings
sigma_aero_deer_human_capt_alt6_PPE <- calc_sigma_aero(ER = 0.53, C_nu = C_nu_human, t_contact = t_contact_deer_human_ranch/60, r = r_deer, V_d = 0.009*runif(nsamples, min = (1-0.994), max = (1-0.509)), nsamples = nsamples)

ranch_suburban.params_alt6 <- alternative(alpha_immunity = alpha_immunity_null,
                            omega_ww = omega_ww_suburban, 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_suburban, sigma_aero_deer_deer_captive = sigma_aero_deer_deer_ranch, sigma_aero_deer_human_wild = sigma_aero_deer_human_suburban, sigma_aero_deer_human_capt = sigma_aero_deer_human_capt_alt6_PPE, 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))

proj_ranch_suburban_alt6 <- run(iter = nsamples, initial_compartments = ranch_suburban.inits.fall, initial_compartments_steady = ranch_suburban.inits.steady, params = ranch_suburban.params_alt6, times = times, name = "PPE, Captive")

#7. Require PPE for humans in wild and captive settings (use derived parameters from 5 and 6)

ranch_suburban.params_alt7 <- alternative(alpha_immunity = alpha_immunity_null,
                            omega_ww = omega_ww_suburban, 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_suburban, sigma_aero_deer_deer_captive = sigma_aero_deer_deer_ranch, sigma_aero_deer_human_wild = sigma_aero_deer_human_wild_alt5_PPE, sigma_aero_deer_human_capt = sigma_aero_deer_human_capt_alt6_PPE, 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))

proj_ranch_suburban_alt7 <- run(iter = nsamples, initial_compartments = ranch_suburban.inits.fall, initial_compartments_steady = ranch_suburban.inits.steady, params = ranch_suburban.params_alt7, times = times, name = "PPE, Both")

#8. Vaccinate captive deer (start with r_captive = 1)
ranch_suburban.inits_alt8 <- initial_compartments(draws = nsamples, S_captive_prop = 0, R_captive_prop = 1)
ranch_suburban.inits_alt8.steady <- initial_compartments(draws = nsamples, S_captive_prop = 0, R_captive_prop = 1, steady = TRUE)

proj_ranch_suburban_alt8 <- run(iter = nsamples, initial_compartments = ranch_suburban.inits_alt8, initial_compartments_steady = ranch_suburban.inits_alt8.steady,params = ranch_suburban.params, times = times, name = "Vaccinate captive deer")

#9. Vaccinate and boost captive deer (start with r_captive = 1)
ranch_suburban.inits_alt9 <- initial_compartments(draws = nsamples, S_captive_prop = 0, R_captive_prop = 1)
ranch_suburban.inits_alt9.steady <- initial_compartments(draws = nsamples, S_captive_prop = 0, R_captive_prop = 1, steady = TRUE)

boost_alt9 <- rep(0.9,nsamples)

ranch_suburban.params_alt9 <- alternative(alpha_immunity = alpha_immunity_null,
                            omega_ww = omega_ww_suburban, 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_suburban, sigma_aero_deer_deer_captive = sigma_aero_deer_deer_ranch, 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 = boost_alt9)


proj_ranch_suburban_alt9 <- run(iter = nsamples, initial_compartments = ranch_suburban.inits_alt9, initial_compartments_steady = ranch_suburban.inits_alt9.steady, params = ranch_suburban.params_alt9, times = times, name = "Vaccinate and boost captive deer")

#10. Double fence captive facilities
omega_cw_alt10 <- rep(0, nsamples)

ranch_suburban.params_alt10 <- alternative(alpha_immunity = alpha_immunity_null,
                            omega_ww = omega_ww_suburban, omega_cw = omega_cw_alt10, omega_cc = omega_cc_ranch,
                            omega_hw = omega_hw_suburban, omega_hc = omega_hc_ranch,
                            sigma_aero_deer_deer_wild = sigma_aero_deer_deer_suburban, sigma_aero_deer_deer_captive = sigma_aero_deer_deer_ranch, 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))


proj_ranch_suburban_alt10 <- run(iter = nsamples, initial_compartments = ranch_suburban.inits.fall, initial_compartments_steady = ranch_suburban.inits.steady, params = ranch_suburban.params_alt10, times = times, name = "Double-fence captive facilities")

#11. Restrict human interaction (using proportional reduction from suburban to rural)
human_restrictions <- (get_EE_param_vals(data = elicitation_data, my_param = "Deer-Human Proximity Rate, Rural (per 120 days)") / 120)/(get_EE_param_vals(data = elicitation_data, my_param = "Deer-Human Proximity Rate, Suburban (per 120 days)") / 120)

ranch_suburban.params_alt11 <- alternative(alpha_immunity = alpha_immunity_null,
                            omega_ww = omega_ww_suburban, omega_cw = omega_cw, omega_cc = omega_cc_ranch,
                            omega_hw = omega_hw_suburban*human_restrictions, omega_hc = omega_hc_ranch,
                            sigma_aero_deer_deer_wild = sigma_aero_deer_deer_suburban, sigma_aero_deer_deer_captive = sigma_aero_deer_deer_ranch, 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))


proj_ranch_suburban_alt11 <- run(iter = nsamples, initial_compartments = ranch_suburban.inits.fall, initial_compartments_steady = ranch_suburban.inits.steady, params = ranch_suburban.params_alt11, times = times, name = "Restrict human interactions")

#12. Improve indoor facility air exchange (does not apply to outdoor ranch facilities)

proj_ranch_suburban <- rbind(proj_ranch_suburban_null, proj_ranch_suburban_alt1, proj_ranch_suburban_alt2, proj_ranch_suburban_alt3, proj_ranch_suburban_alt4, proj_ranch_suburban_alt5, proj_ranch_suburban_alt6, proj_ranch_suburban_alt7, proj_ranch_suburban_alt8, proj_ranch_suburban_alt9, proj_ranch_suburban_alt10, proj_ranch_suburban_alt11)

Calcualte prevalence, cumulative infections, and persistence for interventions in the outdoor ranch : wild, suburban system

persist.threshold <- 0.001

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(., persistWild = I_wild > persist.threshold, persistCaptive = I_captive > persist.threshold) %>%
  select(., run_id, persistWild, persistCaptive) %>%
  arrange(., run_id) -> Persist_ranch_suburban_df

proj_ranch_suburban %>%
  mutate(ode_df = map(ode_proj, as.data.frame)) %>%
  pull(ode_df) %>%
  list_rbind(names_to = "run_id") %>% 
  mutate(., Complex = rep("Outdoor ranch and suburban deer",nrow(.)), Alternative = c(rep(unique(proj_ranch_suburban$Context)[1],length(times)*nsamples),rep(unique(proj_ranch_suburban$Context)[2],length(times)*nsamples), rep(unique(proj_ranch_suburban$Context)[3],length(times)*nsamples),rep(unique(proj_ranch_suburban$Context)[4],length(times)*nsamples), rep(unique(proj_ranch_suburban$Context)[5],length(times)*nsamples), rep(unique(proj_ranch_suburban$Context)[6],length(times)*nsamples), rep(unique(proj_ranch_suburban$Context)[7],length(times)*nsamples),rep(unique(proj_ranch_suburban$Context)[8],length(times)*nsamples),rep(unique(proj_ranch_suburban$Context)[9],length(times)*nsamples),rep(unique(proj_ranch_suburban$Context)[10],length(times)*nsamples),rep(unique(proj_ranch_suburban$Context)[11],length(times)*nsamples),rep(unique(proj_ranch_suburban$Context)[12],length(times)*nsamples))) %>% 
  group_by(Complex,Alternative, run_id) %>% 
  summarize(meanWild = mean(I_wild), meanCaptive = mean(I_captive), 
            cumulativeWild = last(I_wild_cumulative), cumulativeCaptive = last(I_captive_cumulative)) %>% 
   merge(.,Persist_ranch_suburban_df, by = "run_id")-> alternatives_for_ranch_suburban_complex

We can see the result of this processing, is a dataframe listing outbreak metrics for all iterations of each management alternative (including the default baseline). We will repeat this for the other systems, and then use summary functions to compare the effect of each alternative to the baseline, for each system.

head(alternatives_for_ranch_suburban_complex)
#>   run_id                         Complex Alternative     meanWild  meanCaptive
#> 1      1 Outdoor ranch and suburban deer    Baseline 6.676796e-02 0.0688039960
#> 2      2 Outdoor ranch and suburban deer    Baseline 9.735172e-05 0.0000827675
#> 3      3 Outdoor ranch and suburban deer    Baseline 5.937537e-02 0.0598487381
#> 4      4 Outdoor ranch and suburban deer    Baseline 3.225770e-05 0.0002175466
#> 5      5 Outdoor ranch and suburban deer    Baseline 4.790617e-02 0.0479719483
#> 6      6 Outdoor ranch and suburban deer    Baseline 2.432780e-02 0.0266367681
#>   cumulativeWild cumulativeCaptive persistWild persistCaptive
#> 1    1.397087341       1.435637721        TRUE           TRUE
#> 2    0.002060018       0.001751416       FALSE          FALSE
#> 3    1.236692239       1.246979083        TRUE           TRUE
#> 4    0.000686242       0.004626696       FALSE          FALSE
#> 5    0.967220064       0.968339379        TRUE           TRUE
#> 6    0.537436553       0.573908352        TRUE           TRUE

Combined system of outdoor ranch facilities and wild deer in a rural setting

ranch_rural.inits.fall <- initial_compartments(draws = nsamples)
ranch_rural.inits.steady <- initial_compartments(draws = nsamples, steady = TRUE)

ranch_rural.params <- alternative(alpha_immunity = alpha_immunity_null,
                            omega_ww = omega_ww_rural, 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_rural, sigma_aero_deer_deer_captive = sigma_aero_deer_deer_ranch, 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))

proj_ranch_rural_null <- run(iter = nsamples, initial_compartments = ranch_rural.inits.fall, initial_compartments_steady = ranch_rural.inits.steady, params = ranch_rural.params, times = times, name = "Baseline")

#Alternatives
#1. 10% density reduction
omega_ww_alt1_density.red.10 <- whitetailedSIRS::calc_contact_rate(nsamples = nsamples, type_contact = habitat, N_w = nWild*0.9)*(get_EE_param_vals(data = elicitation_data, my_param = "Proximity rate with baiting (17 events without baiting)")/17) 

ranch_rural.params_alt1 <- alternative(alpha_immunity = alpha_immunity_null,
                            omega_ww = omega_ww_alt1_density.red.10, 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_rural, sigma_aero_deer_deer_captive = sigma_aero_deer_deer_ranch, 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))

proj_ranch_rural_alt1 <- run(iter = nsamples, initial_compartments = ranch_rural.inits.fall, initial_compartments_steady = ranch_rural.inits.steady, params = ranch_rural.params_alt1, times = times, name = "Reduce wild density by 10%")

#2. 25% density reduction
omega_ww_alt2_density.red.25 <- whitetailedSIRS::calc_contact_rate(nsamples = nsamples, type_contact = habitat, N_w = nWild*0.75)*(get_EE_param_vals(data = elicitation_data, my_param = "Proximity rate with baiting (17 events without baiting)")/17)

ranch_rural.params_alt2 <- alternative(alpha_immunity = alpha_immunity_null,
                            omega_ww = omega_ww_alt2_density.red.25, 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_rural, sigma_aero_deer_deer_captive = sigma_aero_deer_deer_ranch, 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))

proj_ranch_rural_alt2 <- run(iter = nsamples, initial_compartments = ranch_rural.inits.fall, initial_compartments_steady = ranch_rural.inits.steady, params = ranch_rural.params_alt2, times = times, name = "Reduce wild density by 25%")

#3. Eliminate baiting
omega_ww_alt3_no.baiting <- whitetailedSIRS::calc_contact_rate(nsamples = nsamples, type_contact = habitat, N_w = nWild*0.75)

ranch_rural.params_alt3 <- alternative(alpha_immunity = alpha_immunity_null,
                            omega_ww = omega_ww_alt3_no.baiting, 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_rural, sigma_aero_deer_deer_captive = sigma_aero_deer_deer_ranch, 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))

proj_ranch_rural_alt3 <- run(iter = nsamples, initial_compartments = ranch_rural.inits.fall, initial_compartments_steady = ranch_rural.inits.steady, params = ranch_rural.params_alt3, times = times, name = "Eliminate baiting")

#4. Pause research (assuming random 1-50% of human-wild deer interaction)
omega_hw_alt4_pause.research <- omega_hw_rural*runif(nsamples, min = 0.5, max = 0.99) #reduce human deer proximity rate in the wild up to 50%

ranch_rural.params_alt4 <- alternative(alpha_immunity = alpha_immunity_null,
                            omega_ww = omega_ww_rural, omega_cw = omega_cw, omega_cc = omega_cc_ranch,
                            omega_hw = omega_hw_alt4_pause.research, omega_hc = omega_hc_ranch,
                            sigma_aero_deer_deer_wild = sigma_aero_deer_deer_rural, sigma_aero_deer_deer_captive = sigma_aero_deer_deer_ranch, 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))

proj_ranch_rural_alt4 <- run(iter = nsamples, initial_compartments = ranch_rural.inits.fall, initial_compartments_steady = ranch_rural.inits.steady, params = ranch_rural.params_alt4, times = times, name = "Pause Research")

#5. Require PPE for humans in wild settings
sigma_aero_deer_human_wild_alt5_PPE <- calc_sigma_aero(ER = 0.53, C_nu = C_nu_human, t_contact = t_contact_deer_human_rural/60, r = r_deer, V_d = 0.009*runif(nsamples, min = (1-0.994), max = (1-0.509)), nsamples = nsamples)

ranch_rural.params_alt5 <- alternative(alpha_immunity = alpha_immunity_null,
                            omega_ww = omega_ww_rural, 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_rural, sigma_aero_deer_deer_captive = sigma_aero_deer_deer_ranch, sigma_aero_deer_human_wild = sigma_aero_deer_human_wild_alt5_PPE, 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))

proj_ranch_rural_alt5 <- run(iter = nsamples, initial_compartments = ranch_rural.inits.fall, initial_compartments_steady = ranch_rural.inits.steady, params = ranch_rural.params_alt5, times = times, name = "PPE, Wild")

#6. Require PPE for humans in captive settings
sigma_aero_deer_human_capt_alt6_PPE <- calc_sigma_aero(ER = 0.53, C_nu = C_nu_human, t_contact = t_contact_deer_human_ranch/60, r = r_deer, V_d = 0.009*runif(nsamples, min = (1-0.994), max = (1-0.509)), nsamples = nsamples)

ranch_rural.params_alt6 <- alternative(alpha_immunity = alpha_immunity_null,
                            omega_ww = omega_ww_rural, 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_rural, sigma_aero_deer_deer_captive = sigma_aero_deer_deer_ranch, sigma_aero_deer_human_wild = sigma_aero_deer_human_rural, sigma_aero_deer_human_capt = sigma_aero_deer_human_capt_alt6_PPE, 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))

proj_ranch_rural_alt6 <- run(iter = nsamples, initial_compartments = ranch_rural.inits.fall, initial_compartments_steady = ranch_rural.inits.steady, params = ranch_rural.params_alt6, times = times, name = "PPE, Captive")

#7. Require PPE for humans in wild and captive settings (use derived parameters from 5 and 6)

ranch_rural.params_alt7 <- alternative(alpha_immunity = alpha_immunity_null,
                            omega_ww = omega_ww_rural, 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_rural, sigma_aero_deer_deer_captive = sigma_aero_deer_deer_ranch, sigma_aero_deer_human_wild = sigma_aero_deer_human_wild_alt5_PPE, sigma_aero_deer_human_capt = sigma_aero_deer_human_capt_alt6_PPE, 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))

proj_ranch_rural_alt7 <- run(iter = nsamples, initial_compartments = ranch_rural.inits.fall, initial_compartments_steady = ranch_rural.inits.steady, params = ranch_rural.params_alt7, times = times, name = "PPE, Both")

#8. Vaccinate captive deer (start with r_captive = 1)
ranch_rural.inits_alt8 <- initial_compartments(draws = nsamples, S_captive_prop = 0, R_captive_prop = 1)
ranch_rural.inits_alt8_steady <- initial_compartments(draws = nsamples, S_captive_prop = 0, R_captive_prop = 1, steady = T)

proj_ranch_rural_alt8 <- run(iter = nsamples, initial_compartments = ranch_rural.inits_alt8, initial_compartments_steady = ranch_rural.inits_alt8_steady, params = ranch_rural.params, times = times, name = "Vaccinate captive deer")

#9. Vaccinate and boost captive deer (start with r_captive = 1)
ranch_rural.inits_alt9 <- initial_compartments(draws = nsamples, S_captive_prop = 0, R_captive_prop = 1)
ranch_rural.inits_alt9_steady <- initial_compartments(draws = nsamples, S_captive_prop = 0, R_captive_prop = 1, steady = T)


boost_alt9 <- rep(0.9,nsamples)

ranch_rural.params_alt9 <- alternative(alpha_immunity = alpha_immunity_null,
                            omega_ww = omega_ww_rural, 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_rural, sigma_aero_deer_deer_captive = sigma_aero_deer_deer_ranch, 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 = boost_alt9)


proj_ranch_rural_alt9 <- run(iter = nsamples, initial_compartments = ranch_rural.inits_alt9, initial_compartments_steady = ranch_rural.inits_alt9_steady, params = ranch_rural.params_alt9, times = times, name = "Vaccinate and boost captive deer")

#10. Double fence captive facilities
omega_cw_alt10 <- rep(0, nsamples)

ranch_rural.params_alt10 <- alternative(alpha_immunity = alpha_immunity_null,
                            omega_ww = omega_ww_rural, omega_cw = omega_cw_alt10, omega_cc = omega_cc_ranch,
                            omega_hw = omega_hw_rural, omega_hc = omega_hc_ranch,
                            sigma_aero_deer_deer_wild = sigma_aero_deer_deer_rural, sigma_aero_deer_deer_captive = sigma_aero_deer_deer_ranch, 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))


proj_ranch_rural_alt10 <- run(iter = nsamples, initial_compartments = ranch_rural.inits.fall, initial_compartments_steady = ranch_rural.inits.steady, params = ranch_rural.params_alt10, times = times, name = "Double-fence captive facilities")

#11. Restrict human interaction (using proportional reduction from suburban to rural)
human_restrictions <- (get_EE_param_vals(data = elicitation_data, my_param = "Deer-Human Proximity Rate, Rural (per 120 days)") / 120)/(get_EE_param_vals(data = elicitation_data, my_param = "Deer-Human Proximity Rate, Suburban (per 120 days)") / 120)

ranch_rural.params_alt11 <- alternative(alpha_immunity = alpha_immunity_null,
                            omega_ww = omega_ww_rural, omega_cw = omega_cw, omega_cc = omega_cc_ranch,
                            omega_hw = omega_hw_rural*human_restrictions, omega_hc = omega_hc_ranch,
                            sigma_aero_deer_deer_wild = sigma_aero_deer_deer_rural, sigma_aero_deer_deer_captive = sigma_aero_deer_deer_ranch, 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))


proj_ranch_rural_alt11 <- run(iter = nsamples, initial_compartments = ranch_rural.inits.fall, initial_compartments_steady = ranch_rural.inits.steady, params = ranch_rural.params_alt11, times = times, name = "Restrict human interactions")

#12. Improve indoor facility air exchange (does not apply to outdoor ranch facilities)

proj_ranch_rural <- rbind(proj_ranch_rural_null, proj_ranch_rural_alt1, proj_ranch_rural_alt2, proj_ranch_rural_alt3, proj_ranch_rural_alt4, proj_ranch_rural_alt5, proj_ranch_rural_alt6, proj_ranch_rural_alt7, proj_ranch_rural_alt8, proj_ranch_rural_alt9, proj_ranch_rural_alt10, proj_ranch_rural_alt11)

Calculate prevalence, cumulative infections, and persistence for interventions in the outdoor ranch : wild, rural system

persist.threshold <- 0.001

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(., persistWild = I_wild > persist.threshold, persistCaptive = I_captive > persist.threshold) %>%
  select(., run_id, persistWild, persistCaptive) %>%
  arrange(., run_id) -> Persist_ranch_rural_df

proj_ranch_rural %>%
  mutate(ode_df = map(ode_proj, as.data.frame)) %>%
  pull(ode_df) %>%
  list_rbind(names_to = "run_id") %>% 
  mutate(., Complex = rep("Outdoor ranch and rural deer",nrow(.)), Alternative = c(rep(unique(proj_ranch_rural$Context)[1],length(times)*nsamples),rep(unique(proj_ranch_rural$Context)[2],length(times)*nsamples), rep(unique(proj_ranch_rural$Context)[3],length(times)*nsamples),rep(unique(proj_ranch_rural$Context)[4],length(times)*nsamples), rep(unique(proj_ranch_rural$Context)[5],length(times)*nsamples), rep(unique(proj_ranch_rural$Context)[6],length(times)*nsamples), rep(unique(proj_ranch_rural$Context)[7],length(times)*nsamples),rep(unique(proj_ranch_rural$Context)[8],length(times)*nsamples),rep(unique(proj_ranch_rural$Context)[9],length(times)*nsamples),rep(unique(proj_ranch_rural$Context)[10],length(times)*nsamples),rep(unique(proj_ranch_rural$Context)[11],length(times)*nsamples),rep(unique(proj_ranch_rural$Context)[12],length(times)*nsamples))) %>% 
  group_by(Complex,Alternative, run_id) %>% 
  summarize(meanWild = mean(I_wild), meanCaptive = mean(I_captive), 
            cumulativeWild = last(I_wild_cumulative), cumulativeCaptive = last(I_captive_cumulative)) %>% 
   merge(.,Persist_ranch_rural_df, by = "run_id")-> alternatives_for_ranch_rural_complex

Combined system of intensive captive facilities and wild deer in a suburban setting

intensive_suburban.inits.fall <- initial_compartments(draws = nsamples)
intensive_suburban.inits.steady <- initial_compartments(draws = nsamples, steady = TRUE)


intensive_suburban.params <- alternative(alpha_immunity = alpha_immunity_null,
                            omega_ww = omega_ww_suburban, 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_suburban, 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))

proj_intensive_suburban_null <- run(iter = nsamples, initial_compartments = intensive_suburban.inits.fall, initial_compartments_steady = intensive_suburban.inits.steady, params = intensive_suburban.params, times = times, name = "Baseline")

#Alternatives
#1. 10% density reduction
omega_ww_alt1_density.red.10 <- whitetailedSIRS::calc_contact_rate(nsamples = nsamples, type_contact = habitat, N_w = nWild*0.9)*(get_EE_param_vals(data = elicitation_data, my_param = "Proximity rate with baiting (17 events without baiting)")/17) 

intensive_suburban.params_alt1 <- alternative(alpha_immunity = alpha_immunity_null,
                            omega_ww = omega_ww_alt1_density.red.10, 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_suburban, 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))

proj_intensive_suburban_alt1 <- run(iter = nsamples, initial_compartments = intensive_suburban.inits.fall, initial_compartments_steady = intensive_suburban.inits.steady, params = intensive_suburban.params_alt1, times = times, name = "Reduce wild density by 10%")

#2. 25% density reduction
omega_ww_alt2_density.red.25 <- whitetailedSIRS::calc_contact_rate(nsamples = nsamples, type_contact = habitat, N_w = nWild*0.75)*(get_EE_param_vals(data = elicitation_data, my_param = "Proximity rate with baiting (17 events without baiting)")/17)

intensive_suburban.params_alt2 <- alternative(alpha_immunity = alpha_immunity_null,
                            omega_ww = omega_ww_alt2_density.red.25, 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_suburban, 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))

proj_intensive_suburban_alt2 <- run(iter = nsamples, initial_compartments = intensive_suburban.inits.fall, initial_compartments_steady = intensive_suburban.inits.steady, params = intensive_suburban.params_alt2, times = times, name = "Reduce wild density by 25%")

#3. Eliminate baiting
omega_ww_alt3_no.baiting <- whitetailedSIRS::calc_contact_rate(nsamples = nsamples, type_contact = habitat, N_w = nWild*0.75)

intensive_suburban.params_alt3 <- alternative(alpha_immunity = alpha_immunity_null,
                            omega_ww = omega_ww_alt3_no.baiting, 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_suburban, 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))

proj_intensive_suburban_alt3 <- run(iter = nsamples, initial_compartments = intensive_suburban.inits.fall, initial_compartments_steady = intensive_suburban.inits.steady, params = intensive_suburban.params_alt3, times = times, name = "Eliminate baiting")

#4. Pause research (assuming random 1-50% of human-wild deer interaction)
omega_hw_alt4_pause.research <- omega_hw_suburban*runif(nsamples, min = 0.5, max = 0.99) #reduce human deer proximity rate in the wild up to 50%

intensive_suburban.params_alt4 <- alternative(alpha_immunity = alpha_immunity_null,
                            omega_ww = omega_ww_suburban, omega_cw = omega_cw, omega_cc = omega_cc_intensive,
                            omega_hw = omega_hw_alt4_pause.research, omega_hc = omega_hc_intensive,
                            sigma_aero_deer_deer_wild = sigma_aero_deer_deer_suburban, 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))

proj_intensive_suburban_alt4 <- run(iter = nsamples, initial_compartments = intensive_suburban.inits.fall, initial_compartments_steady = intensive_suburban.inits.steady, params = intensive_suburban.params_alt4, times = times, name = "Pause Research")

#5. Require PPE for humans in wild settings
sigma_aero_deer_human_wild_alt5_PPE <- calc_sigma_aero(ER = 0.53, C_nu = C_nu_human, t_contact = t_contact_deer_human_suburban/60, r = r_deer, V_d = 0.009*runif(nsamples, min = (1-0.994), max = (1-0.509)), nsamples = nsamples)

intensive_suburban.params_alt5 <- alternative(alpha_immunity = alpha_immunity_null,
                            omega_ww = omega_ww_suburban, 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_suburban, sigma_aero_deer_deer_captive = sigma_aero_deer_deer_intensive, sigma_aero_deer_human_wild = sigma_aero_deer_human_wild_alt5_PPE, 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))

proj_intensive_suburban_alt5 <- run(iter = nsamples, initial_compartments = intensive_suburban.inits.fall, initial_compartments_steady = intensive_suburban.inits.steady, params = intensive_suburban.params_alt5, times = times, name = "PPE, Wild")

#6. Require PPE for humans in captive settings
sigma_aero_deer_human_capt_alt6_PPE <- calc_sigma_aero(ER = 0.53, C_nu = C_nu_human, t_contact = t_contact_deer_human_intensive/60, r = r_deer, V_d = 0.009*runif(nsamples, min = (1-0.994), max = (1-0.509)), nsamples = nsamples)

intensive_suburban.params_alt6 <- alternative(alpha_immunity = alpha_immunity_null,
                            omega_ww = omega_ww_suburban, 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_suburban, 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_capt_alt6_PPE, 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))

proj_intensive_suburban_alt6 <- run(iter = nsamples, initial_compartments = intensive_suburban.inits.fall, initial_compartments_steady = intensive_suburban.inits.steady, params = intensive_suburban.params_alt6, times = times, name = "PPE, Captive")

#7. Require PPE for humans in wild and captive settings (use derived parameters from 5 and 6)

intensive_suburban.params_alt7 <- alternative(alpha_immunity = alpha_immunity_null,
                            omega_ww = omega_ww_suburban, 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_suburban, sigma_aero_deer_deer_captive = sigma_aero_deer_deer_intensive, sigma_aero_deer_human_wild = sigma_aero_deer_human_wild_alt5_PPE, sigma_aero_deer_human_capt = sigma_aero_deer_human_capt_alt6_PPE, 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))

proj_intensive_suburban_alt7 <- run(iter = nsamples, initial_compartments = intensive_suburban.inits.fall, initial_compartments_steady = intensive_suburban.inits.steady, params =  intensive_suburban.params_alt7, times = times, name = "PPE, Both")

#8. Vaccinate captive deer (start with r_captive = 1)
intensive_suburban.inits_alt8.fall <- initial_compartments(draws = nsamples, S_captive_prop = 0, R_captive_prop = 1)
intensive_suburban.inits_alt8.steady <- initial_compartments(draws = nsamples, S_captive_prop = 0, R_captive_prop = 1, steady = TRUE)

proj_intensive_suburban_alt8 <- run(iter = nsamples, initial_compartments = intensive_suburban.inits_alt8.fall, initial_compartments_steady = intensive_suburban.inits_alt8.steady, params = intensive_suburban.params, times = times, name = "Vaccinate captive deer")

#9. Vaccinate and boost captive deer (start with r_captive = 1)
intensive_suburban.inits_alt9.fall <- initial_compartments(draws = nsamples, S_captive_prop = 0, R_captive_prop = 1)
intensive_suburban.inits_alt9.steady <- initial_compartments(draws = nsamples, S_captive_prop = 0, R_captive_prop = 1, steady = TRUE)

boost_alt9 <- rep(0.9,nsamples)

intensive_suburban.params_alt9 <- alternative(alpha_immunity = alpha_immunity_null,
                            omega_ww = omega_ww_suburban, 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_suburban, 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 = boost_alt9)


proj_intensive_suburban_alt9 <- run(iter = nsamples, initial_compartments = intensive_suburban.inits_alt9.fall, initial_compartments_steady = intensive_suburban.inits_alt9.steady, params = intensive_suburban.params_alt9, times = times, name = "Vaccinate and boost captive deer")

#10. Double fence captive facilities
omega_cw_alt10 <- rep(0, nsamples)

intensive_suburban.params_alt10 <- alternative(alpha_immunity = alpha_immunity_null,
                            omega_ww = omega_ww_suburban, omega_cw = omega_cw_alt10, omega_cc = omega_cc_intensive,
                            omega_hw = omega_hw_suburban, omega_hc = omega_hc_intensive,
                            sigma_aero_deer_deer_wild = sigma_aero_deer_deer_suburban, 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))


proj_intensive_suburban_alt10 <- run(iter = nsamples, initial_compartments = intensive_suburban.inits.fall, initial_compartments_steady = intensive_suburban.inits.steady, params = intensive_suburban.params_alt10, times = times, name = "Double-fence captive facilities")

#11. Restrict human interaction (using proportional reduction from suburban to rural)
human_restrictions <- (get_EE_param_vals(data = elicitation_data, my_param = "Deer-Human Proximity Rate, Rural (per 120 days)") / 120)/(get_EE_param_vals(data = elicitation_data, my_param = "Deer-Human Proximity Rate, Suburban (per 120 days)") / 120)

intensive_suburban.params_alt11 <- alternative(alpha_immunity = alpha_immunity_null,
                            omega_ww = omega_ww_suburban, omega_cw = omega_cw, omega_cc = omega_cc_intensive,
                            omega_hw = omega_hw_suburban*human_restrictions, omega_hc = omega_hc_intensive,
                            sigma_aero_deer_deer_wild = sigma_aero_deer_deer_suburban, 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))


proj_intensive_suburban_alt11 <- run(iter = nsamples, initial_compartments = intensive_suburban.inits.fall, initial_compartments_steady = intensive_suburban.inits.steady, params = intensive_suburban.params_alt11, times = times, name = "Restrict human interactions")

#12. Improve indoor facility air exchange
sigma_aero_deer_human_intensive_alt12 <- calc_sigma_aero(ER = 0.53, C_nu = C_nu_human, t_contact = t_contact_deer_human_intensive/60, r = r_deer, nsamples = nsamples)

sigma_aero_deer_deer_intensive_alt12 <- calc_sigma_aero(C_nu = C_nu_deer, t_contact = t_contact_deer_deer_null/60, r = r_deer, nsamples = nsamples)

intensive_suburban.params_alt12 <- alternative(alpha_immunity = alpha_immunity_null,
                            omega_ww = omega_ww_suburban, 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_suburban, sigma_aero_deer_deer_captive = sigma_aero_deer_deer_intensive_alt12, sigma_aero_deer_human_wild = sigma_aero_deer_human_suburban, sigma_aero_deer_human_capt = sigma_aero_deer_human_intensive_alt12, 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))


proj_intensive_suburban_alt12 <- run(iter = nsamples, initial_compartments = intensive_suburban.inits.fall, initial_compartments_steady = intensive_suburban.inits.steady, params = intensive_suburban.params_alt12, times = times, name = "Improve indoor facility air quality")

proj_intensive_suburban <- rbind(proj_intensive_suburban_null, proj_intensive_suburban_alt1, proj_intensive_suburban_alt2, proj_intensive_suburban_alt3, proj_intensive_suburban_alt4, proj_intensive_suburban_alt5, proj_intensive_suburban_alt6, proj_intensive_suburban_alt7, proj_intensive_suburban_alt8, proj_intensive_suburban_alt9, proj_intensive_suburban_alt10, proj_intensive_suburban_alt11, proj_intensive_suburban_alt12)

Calculate prevalence, cumulative infections, and persistence for interventions in the intensive facility : wild, suburban system

persist.threshold <- 0.001

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(., persistWild = I_wild > persist.threshold, persistCaptive = I_captive > persist.threshold) %>%
  select(., run_id, persistWild, persistCaptive) %>%
  arrange(., run_id) -> Persist_intensive_suburban_df


proj_intensive_suburban %>%
  mutate(ode_df = map(ode_proj, as.data.frame)) %>%
  pull(ode_df) %>%
  list_rbind(names_to = "run_id") %>% 
  mutate(., Complex = rep("Intensive facility and suburban deer",nrow(.)), Alternative = c(rep(unique(proj_intensive_suburban$Context)[1],length(times)*nsamples),rep(unique(proj_intensive_suburban$Context)[2],length(times)*nsamples), rep(unique(proj_intensive_suburban$Context)[3],length(times)*nsamples),rep(unique(proj_intensive_suburban$Context)[4],length(times)*nsamples), rep(unique(proj_intensive_suburban$Context)[5],length(times)*nsamples), rep(unique(proj_intensive_suburban$Context)[6],length(times)*nsamples), rep(unique(proj_intensive_suburban$Context)[7],length(times)*nsamples),rep(unique(proj_intensive_suburban$Context)[8],length(times)*nsamples),rep(unique(proj_intensive_suburban$Context)[9],length(times)*nsamples),rep(unique(proj_intensive_suburban$Context)[10],length(times)*nsamples),rep(unique(proj_intensive_suburban$Context)[11],length(times)*nsamples),rep(unique(proj_intensive_suburban$Context)[12],length(times)*nsamples), rep(unique(proj_intensive_suburban$Context)[13],length(times)*nsamples))) %>% 
  group_by(Complex,Alternative, run_id) %>% 
  summarize(meanWild = mean(I_wild), meanCaptive = mean(I_captive), 
            cumulativeWild = last(I_wild_cumulative), cumulativeCaptive = last(I_captive_cumulative)) %>% 
  merge(.,Persist_intensive_suburban_df, by = "run_id") -> alternatives_for_intensive_suburban_complex

Combined system of intensive captive facilities and wild deer in a rural setting

intensive_rural.inits.fall <- initial_compartments(draws = nsamples)
intensive_rural.inits.steady <- initial_compartments(draws = nsamples, steady = TRUE)


intensive_rural.params <- alternative(alpha_immunity = alpha_immunity_null,
                            omega_ww = omega_ww_rural, 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_rural, 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))

proj_intensive_rural_null <- run(iter = nsamples, initial_compartments = intensive_rural.inits.fall, initial_compartments_steady = intensive_rural.inits.steady, params = intensive_rural.params, times = times, name = "Baseline")

#Alternatives
#1. 10% density reduction
omega_ww_alt1_density.red.10 <- whitetailedSIRS::calc_contact_rate(nsamples = nsamples, type_contact = habitat, N_w = nWild*0.9)*(get_EE_param_vals(data = elicitation_data, my_param = "Proximity rate with baiting (17 events without baiting)")/17) 

intensive_rural.params_alt1 <- alternative(alpha_immunity = alpha_immunity_null,
                            omega_ww = omega_ww_alt1_density.red.10, 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_rural, 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))

proj_intensive_rural_alt1 <- run(iter = nsamples, initial_compartments = intensive_rural.inits.fall, initial_compartments_steady = intensive_rural.inits.steady, params = intensive_rural.params_alt1, times = times, name = "Reduce wild density by 10%")

#2. 25% density reduction
omega_ww_alt2_density.red.25 <- whitetailedSIRS::calc_contact_rate(nsamples = nsamples, type_contact = habitat, N_w = nWild*0.75)*(get_EE_param_vals(data = elicitation_data, my_param = "Proximity rate with baiting (17 events without baiting)")/17)

intensive_rural.params_alt2 <- alternative(alpha_immunity = alpha_immunity_null,
                            omega_ww = omega_ww_alt2_density.red.25, 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_rural, 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))

proj_intensive_rural_alt2 <- run(iter = nsamples, initial_compartments = intensive_rural.inits.fall, initial_compartments_steady = intensive_rural.inits.steady, params = intensive_rural.params_alt2, times = times, name = "Reduce wild density by 25%")

#3. Eliminate baiting
omega_ww_alt3_no.baiting <- whitetailedSIRS::calc_contact_rate(nsamples = nsamples, type_contact = habitat, N_w = nWild*0.75)

intensive_rural.params_alt3 <- alternative(alpha_immunity = alpha_immunity_null,
                            omega_ww = omega_ww_alt3_no.baiting, 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_rural, 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))

proj_intensive_rural_alt3 <- run(iter = nsamples, initial_compartments = intensive_rural.inits.fall, initial_compartments_steady = intensive_rural.inits.steady, params = intensive_rural.params_alt3, times = times, name = "Eliminate baiting")

#4. Pause research (assuming random 1-50% of human-wild deer interaction)
omega_hw_alt4_pause.research <- omega_hw_rural*runif(nsamples, min = 0.5, max = 0.99) #reduce human deer proximity rate in the wild up to 50%

intensive_rural.params_alt4 <- alternative(alpha_immunity = alpha_immunity_null,
                            omega_ww = omega_ww_rural, omega_cw = omega_cw, omega_cc = omega_cc_intensive,
                            omega_hw = omega_hw_alt4_pause.research, omega_hc = omega_hc_intensive,
                            sigma_aero_deer_deer_wild = sigma_aero_deer_deer_rural, 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))

proj_intensive_rural_alt4 <- run(iter = nsamples, initial_compartments = intensive_rural.inits.fall, initial_compartments_steady = intensive_rural.inits.steady, params = intensive_rural.params_alt4, times = times, name = "Pause Research")

#5. Require PPE for humans in wild settings
sigma_aero_deer_human_wild_alt5_PPE <- calc_sigma_aero(ER = 0.53, C_nu = C_nu_human, t_contact = t_contact_deer_human_rural/60, r = r_deer, V_d = 0.009*runif(nsamples, min = (1-0.994), max = (1-0.509)), nsamples = nsamples)

intensive_rural.params_alt5 <- alternative(alpha_immunity = alpha_immunity_null,
                            omega_ww = omega_ww_rural, 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_rural, sigma_aero_deer_deer_captive = sigma_aero_deer_deer_intensive, sigma_aero_deer_human_wild = sigma_aero_deer_human_wild_alt5_PPE, 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))

proj_intensive_rural_alt5 <- run(iter = nsamples, initial_compartments = intensive_rural.inits.fall, initial_compartments_steady = intensive_rural.inits.steady, params = intensive_rural.params_alt5, times = times, name = "PPE, Wild")

#6. Require PPE for humans in captive settings
sigma_aero_deer_human_capt_alt6_PPE <- calc_sigma_aero(ER = 0.53, C_nu = C_nu_human, t_contact = t_contact_deer_human_intensive/60, r = r_deer, V_d = 0.009*runif(nsamples, min = (1-0.994), max = (1-0.509)), nsamples = nsamples)

intensive_rural.params_alt6 <- alternative(alpha_immunity = alpha_immunity_null,
                            omega_ww = omega_ww_rural, 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_rural, 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_capt_alt6_PPE, 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))

proj_intensive_rural_alt6 <- run(iter = nsamples, initial_compartments = intensive_rural.inits.fall, initial_compartments_steady = intensive_rural.inits.steady, params = intensive_rural.params_alt6, times = times, name = "PPE, Captive")

#7. Require PPE for humans in wild and captive settings (use derived parameters from 5 and 6)

intensive_rural.params_alt7 <- alternative(alpha_immunity = alpha_immunity_null,
                            omega_ww = omega_ww_rural, 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_rural, sigma_aero_deer_deer_captive = sigma_aero_deer_deer_intensive, sigma_aero_deer_human_wild = sigma_aero_deer_human_wild_alt5_PPE, sigma_aero_deer_human_capt = sigma_aero_deer_human_capt_alt6_PPE, 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))

proj_intensive_rural_alt7 <- run(iter = nsamples, initial_compartments = intensive_rural.inits.fall, initial_compartments_steady = intensive_rural.inits.steady, params = intensive_rural.params_alt7, times = times, name = "PPE, Both")

#8. Vaccinate captive deer (start with r_captive = 1)
intensive_rural.inits.fall_alt8 <- initial_compartments(draws = nsamples, S_captive_prop = 0, R_captive_prop = 1)
intensive_rural.inits.steady_alt8 <- initial_compartments(draws = nsamples, S_captive_prop = 0, R_captive_prop = 1, steady = TRUE)

proj_intensive_rural_alt8 <- run(iter = nsamples, initial_compartments = intensive_rural.inits.fall_alt8, initial_compartments_steady = intensive_rural.inits.steady_alt8, params = intensive_rural.params, times = times, name = "Vaccinate captive deer")

#9. Vaccinate and boost captive deer (start with r_captive = 1)
intensive_rural.inits.fall_alt9 <- initial_compartments(draws = nsamples, S_captive_prop = 0, R_captive_prop = 1)
intensive_rural.inits.steady_alt9 <- initial_compartments(draws = nsamples, S_captive_prop = 0, R_captive_prop = 1, steady = TRUE)

boost_alt9 <- rep(0.9,nsamples)

intensive_rural.params_alt9 <- alternative(alpha_immunity = alpha_immunity_null,
                            omega_ww = omega_ww_rural, 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_rural, 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 = boost_alt9)


proj_intensive_rural_alt9 <- run(iter = nsamples, initial_compartments = intensive_rural.inits.fall_alt9, initial_compartments_steady = intensive_rural.inits.steady_alt9, params = intensive_rural.params_alt9, times = times, name = "Vaccinate and boost captive deer")

#10. Double fence captive facilities
omega_cw_alt10 <- rep(0, nsamples)

intensive_rural.params_alt10 <- alternative(alpha_immunity = alpha_immunity_null,
                            omega_ww = omega_ww_rural, omega_cw = omega_cw_alt10, omega_cc = omega_cc_intensive,
                            omega_hw = omega_hw_rural, omega_hc = omega_hc_intensive,
                            sigma_aero_deer_deer_wild = sigma_aero_deer_deer_rural, 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))


proj_intensive_rural_alt10 <- run(iter = nsamples, initial_compartments = intensive_rural.inits.fall, initial_compartments_steady = intensive_rural.inits.steady, params = intensive_rural.params_alt10, times = times, name = "Double-fence captive facilities")

#11. Restrict human interaction (using proportional reduction from suburban to rural)
human_restrictions <- (get_EE_param_vals(data = elicitation_data, my_param = "Deer-Human Proximity Rate, Rural (per 120 days)") / 120)/(get_EE_param_vals(data = elicitation_data, my_param = "Deer-Human Proximity Rate, Suburban (per 120 days)") / 120)

intensive_rural.params_alt11 <- alternative(alpha_immunity = alpha_immunity_null,
                            omega_ww = omega_ww_rural, omega_cw = omega_cw, omega_cc = omega_cc_intensive,
                            omega_hw = omega_hw_rural*human_restrictions, omega_hc = omega_hc_intensive,
                            sigma_aero_deer_deer_wild = sigma_aero_deer_deer_rural, 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))


proj_intensive_rural_alt11 <- run(iter = nsamples, initial_compartments = intensive_rural.inits.fall, initial_compartments_steady = intensive_rural.inits.steady, params = intensive_rural.params_alt11, times = times, name = "Restrict human interactions")

#12. Improve indoor facility air exchange
sigma_aero_deer_human_intensive_alt12 <- calc_sigma_aero(ER = 0.53, C_nu = C_nu_human, t_contact = t_contact_deer_human_intensive/60, r = r_deer, nsamples = nsamples)

sigma_aero_deer_deer_intensive_alt12 <- calc_sigma_aero(C_nu = C_nu_deer, t_contact = t_contact_deer_deer_null/60, r = r_deer, nsamples = nsamples)

intensive_rural.params_alt12 <- alternative(alpha_immunity = alpha_immunity_null,
                            omega_ww = omega_ww_rural, 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_rural, sigma_aero_deer_deer_captive = sigma_aero_deer_deer_intensive_alt12, sigma_aero_deer_human_wild = sigma_aero_deer_human_rural, sigma_aero_deer_human_capt = sigma_aero_deer_human_intensive_alt12, 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))


proj_intensive_rural_alt12 <- run(iter = nsamples, initial_compartments = intensive_rural.inits.fall, initial_compartments_steady = intensive_rural.inits.steady, params = intensive_rural.params_alt12, times = times, name = "Improve indoor facility air quality")

proj_intensive_rural <- rbind(proj_intensive_rural_null, proj_intensive_rural_alt1, proj_intensive_rural_alt2, proj_intensive_rural_alt3, proj_intensive_rural_alt4, proj_intensive_rural_alt5, proj_intensive_rural_alt6, proj_intensive_rural_alt7, proj_intensive_rural_alt8, proj_intensive_rural_alt9, proj_intensive_rural_alt10, proj_intensive_rural_alt11, proj_intensive_rural_alt12)

Calculate prevalence, cumulative infections, and persistence for interventions in the intensive facility : wild, rural system

persist.threshold <- 0.001

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(., persistWild = I_wild > persist.threshold, persistCaptive = I_captive > persist.threshold) %>%
  select(., run_id, persistWild, persistCaptive) %>%
  arrange(., run_id) -> Persist_intensive_rural_df

proj_intensive_rural %>%
  mutate(ode_df = map(ode_proj, as.data.frame)) %>%
  pull(ode_df) %>%
  list_rbind(names_to = "run_id") %>% 
  mutate(., Complex = rep("Intensive facility and rural deer",nrow(.)), Alternative = c(rep(unique(proj_intensive_rural$Context)[1],length(times)*nsamples),rep(unique(proj_intensive_rural$Context)[2],length(times)*nsamples), rep(unique(proj_intensive_rural$Context)[3],length(times)*nsamples),rep(unique(proj_intensive_rural$Context)[4],length(times)*nsamples), rep(unique(proj_intensive_rural$Context)[5],length(times)*nsamples), rep(unique(proj_intensive_rural$Context)[6],length(times)*nsamples), rep(unique(proj_intensive_rural$Context)[7],length(times)*nsamples),rep(unique(proj_intensive_rural$Context)[8],length(times)*nsamples),rep(unique(proj_intensive_rural$Context)[9],length(times)*nsamples),rep(unique(proj_intensive_rural$Context)[10],length(times)*nsamples),rep(unique(proj_intensive_rural$Context)[11],length(times)*nsamples),rep(unique(proj_intensive_rural$Context)[12],length(times)*nsamples), rep(unique(proj_intensive_rural$Context)[13],length(times)*nsamples))) %>% 
  group_by(Complex,Alternative, run_id) %>% 
  summarize(meanWild = mean(I_wild), meanCaptive = mean(I_captive), 
            cumulativeWild = last(I_wild_cumulative), cumulativeCaptive = last(I_captive_cumulative)) %>% 
  merge(.,Persist_intensive_rural_df, by = "run_id") -> alternatives_for_intensive_rural_complex

4.Combine and summarize results

Below, we combine the results for all alternatives, for all combinations of captive and wild scenarios (systems). These results are saved in the package’s data folder, and can be accessed through whitetailedSIRS::alternatives_across_systems.

alternatives_across_systems <- rbind(alternatives_for_ranch_suburban_complex, alternatives_for_ranch_rural_complex, alternatives_for_intensive_suburban_complex, alternatives_for_intensive_rural_complex)

alternatives_across_systems %>% 
   group_by(., Complex, Alternative) %>% 
   summarize(., Wild_Cumulative = median(cumulativeWild), Captive_Cumulative = median(cumulativeCaptive))
#> `summarise()` has grouped output by 'Complex'. You can override using the
#> `.groups` argument.
#> # A tibble: 50 × 4
#> # Groups:   Complex [4]
#>    Complex                        Alternative Wild_Cumulative Captive_Cumulative
#>    <chr>                          <chr>                 <dbl>              <dbl>
#>  1 Intensive facility and rural … Baseline            0.861                 1.47
#>  2 Intensive facility and rural … Double-fen…         0.494                 1.47
#>  3 Intensive facility and rural … Eliminate …         0.00270               1.47
#>  4 Intensive facility and rural … Improve in…         0.832                 1.14
#>  5 Intensive facility and rural … PPE, Both           0.856                 1.46
#>  6 Intensive facility and rural … PPE, Capti…         0.857                 1.46
#>  7 Intensive facility and rural … PPE, Wild           0.861                 1.47
#>  8 Intensive facility and rural … Pause Rese…         0.861                 1.47
#>  9 Intensive facility and rural … Reduce wil…         0.783                 1.47
#> 10 Intensive facility and rural … Reduce wil…         0.521                 1.47
#> # ℹ 40 more rows

5.Impact of Cross-Sector Collaboration

Finally, we simulated SARS-CoV-2 outbreaks in each system with multiple alternatives applied. Here, we implement the most effective alternative for each sector (agriculture, public health, and wildlife) in concert to understand the cumulative effects of coordinated actions across sectors. The results of this final simulation are stored in whitetailedSIRS::cross_sector_results.

#Vaccinate..
cross_sector.inits.fall <- initial_compartments(draws = nsamples, S_captive_prop = 0, R_captive_prop = 1)
cross_sector.inits.steady <- initial_compartments(draws = nsamples, S_captive_prop = 0, R_captive_prop = 1, steady = TRUE)

#...and boost
boost_alt9 <- rep(0.9,nsamples)

#PPE for all interactions
sigma_aero_deer_human_capt_alt6_PPE <- calc_sigma_aero(ER = 0.53, C_nu = C_nu_human, t_contact = t_contact_deer_human_ranch/60, r = r_deer, V_d = 0.009*runif(nsamples, min = (1-0.994), max = (1-0.509)), nsamples = nsamples)

sigma_aero_deer_human_wild_alt5_PPE <- calc_sigma_aero(ER = 0.53, C_nu = C_nu_human, t_contact = t_contact_deer_human_rural/60, r = r_deer, V_d = 0.009*runif(nsamples, min = (1-0.994), max = (1-0.509)), nsamples = nsamples)

#Eliminate baiting
omega_ww_alt3_no.baiting <- whitetailedSIRS::calc_contact_rate(nsamples = nsamples, type_contact = habitat, N_w = nWild*0.75)

#Ranch and suburban complex
ranch_suburban_cross_sector.params <- alternative(alpha_immunity = alpha_immunity_null,
                            omega_ww = omega_ww_alt3_no.baiting, 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_suburban, sigma_aero_deer_deer_captive = sigma_aero_deer_deer_ranch, sigma_aero_deer_human_wild = sigma_aero_deer_human_wild_alt5_PPE, sigma_aero_deer_human_capt = sigma_aero_deer_human_capt_alt6_PPE, 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 = boost_alt9)

proj_ranch_suburban_cross_sector <- run(iter = nsamples, initial_compartments = cross_sector.inits.fall, initial_compartments_steady = cross_sector.inits.steady, params =  ranch_suburban_cross_sector.params, times = times, name = "Cross-Sector Collaboration: Vaccinate and Boost, PPE for all interactions, and eliminating baiting")

proj_ranch_suburban_cross_sector %>%
  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(., persistWild = I_wild > persist.threshold, persistCaptive = I_captive > persist.threshold) %>%
  select(., run_id, persistWild, persistCaptive) %>%
  arrange(., run_id) -> Persist_ranch_suburban_cross_sector_df

proj_ranch_suburban_cross_sector %>%
  mutate(ode_df = map(ode_proj, as.data.frame)) %>%
  pull(ode_df) %>%
  list_rbind(names_to = "run_id") %>% 
  mutate(., Complex = rep("Ranch and suburban deer",nrow(.)), Alternative = c(rep(unique(proj_ranch_suburban_cross_sector$Context)[1],length(times)*nsamples))) %>% 
  group_by(Complex,Alternative, run_id) %>% 
  summarize(meanWild = mean(I_wild), meanCaptive = mean(I_captive), 
            cumulativeWild = last(I_wild_cumulative), cumulativeCaptive = last(I_captive_cumulative)) %>% 
   merge(.,Persist_ranch_suburban_cross_sector_df, by = "run_id") -> cross_sector_for_ranch_suburban_complex
#> `summarise()` has grouped output by 'Complex', 'Alternative'. You can override
#> using the `.groups` argument.

#Ranch and rural complex
ranch_rural_cross_sector.params <- alternative(alpha_immunity = alpha_immunity_null,
                            omega_ww = omega_ww_alt3_no.baiting, 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_rural, sigma_aero_deer_deer_captive = sigma_aero_deer_deer_ranch, sigma_aero_deer_human_wild = sigma_aero_deer_human_wild_alt5_PPE, sigma_aero_deer_human_capt = sigma_aero_deer_human_capt_alt6_PPE, 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 = boost_alt9)

proj_ranch_rural_cross_sector <- run(iter = nsamples, initial_compartments = cross_sector.inits.fall, initial_compartments_steady = cross_sector.inits.steady, params = ranch_rural_cross_sector.params, times = times, name = "Cross-Sector Collaboration: Vaccinate and Boost, PPE for all interactions, and eliminating baiting")

proj_ranch_rural_cross_sector %>%
  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(., persistWild = I_wild > persist.threshold, persistCaptive = I_captive > persist.threshold) %>%
  select(., run_id, persistWild, persistCaptive) %>%
  arrange(., run_id) -> Persist_ranch_rural_cross_sector_df

proj_ranch_rural_cross_sector %>%
  mutate(ode_df = map(ode_proj, as.data.frame)) %>%
  pull(ode_df) %>%
  list_rbind(names_to = "run_id") %>% 
  mutate(., Complex = rep("Ranch and rural deer",nrow(.)), Alternative = c(rep(unique(proj_ranch_rural_cross_sector$Context)[1],length(times)*nsamples))) %>% 
  group_by(Complex,Alternative, run_id) %>% 
  summarize(meanWild = mean(I_wild), meanCaptive = mean(I_captive), 
            cumulativeWild = last(I_wild_cumulative), cumulativeCaptive = last(I_captive_cumulative)) %>% 
  merge(.,Persist_ranch_rural_cross_sector_df, by = "run_id")-> cross_sector_for_ranch_rural_complex
#> `summarise()` has grouped output by 'Complex', 'Alternative'. You can override
#> using the `.groups` argument.

#Intensive facility and suburban complex
intensive_suburban_cross_sector.params <- alternative(alpha_immunity = alpha_immunity_null,
                            omega_ww = omega_ww_alt3_no.baiting, 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_suburban, sigma_aero_deer_deer_captive = sigma_aero_deer_deer_intensive, sigma_aero_deer_human_wild = sigma_aero_deer_human_wild_alt5_PPE, sigma_aero_deer_human_capt = sigma_aero_deer_human_capt_alt6_PPE, 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 = boost_alt9)

proj_intensive_suburban_cross_sector <- run(iter = nsamples, initial_compartments = cross_sector.inits.fall, initial_compartments_steady = cross_sector.inits.steady, params = intensive_suburban_cross_sector.params, times = times, name = "Cross-Sector Collaboration: Vaccinate and Boost, PPE for all interactions, and eliminating baiting")

proj_intensive_suburban_cross_sector %>%
  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(., persistWild = I_wild > persist.threshold, persistCaptive = I_captive > persist.threshold) %>%
  select(., run_id, persistWild, persistCaptive) %>%
  arrange(., run_id) -> Persist_intensive_suburban_cross_sector_df

proj_intensive_suburban_cross_sector %>%
  mutate(ode_df = map(ode_proj, as.data.frame)) %>%
  pull(ode_df) %>%
  list_rbind(names_to = "run_id") %>% 
  mutate(., Complex = rep("Intensive facility and suburban deer",nrow(.)), Alternative = c(rep(unique(proj_intensive_suburban_cross_sector$Context)[1],length(times)*nsamples))) %>% 
  group_by(Complex,Alternative, run_id) %>% 
  summarize(meanWild = mean(I_wild), meanCaptive = mean(I_captive), 
            cumulativeWild = last(I_wild_cumulative), cumulativeCaptive = last(I_captive_cumulative)) %>% 
   merge(.,Persist_intensive_suburban_cross_sector_df, by = "run_id")-> cross_sector_for_intensive_suburban_complex
#> `summarise()` has grouped output by 'Complex', 'Alternative'. You can override
#> using the `.groups` argument.

#Intensive facility and rural complex
intensive_rural_cross_sector.params <- alternative(alpha_immunity = alpha_immunity_null,
                            omega_ww = omega_ww_alt3_no.baiting, 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_rural, sigma_aero_deer_deer_captive = sigma_aero_deer_deer_intensive, sigma_aero_deer_human_wild = sigma_aero_deer_human_wild_alt5_PPE, sigma_aero_deer_human_capt = sigma_aero_deer_human_capt_alt6_PPE, 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 = boost_alt9)

proj_intensive_rural_cross_sector <- run(iter = nsamples, initial_compartments = cross_sector.inits.fall, initial_compartments_steady = cross_sector.inits.steady, params = intensive_rural_cross_sector.params, times = times, name = "Cross-Sector Collaboration: Vaccinate and Boost, PPE for all interactions, and eliminating baiting")

proj_intensive_rural_cross_sector %>%
  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(., persistWild = I_wild > persist.threshold, persistCaptive = I_captive > persist.threshold) %>%
  select(., run_id, persistWild, persistCaptive) %>%
  arrange(., run_id) -> Persist_intensive_rural_cross_sector_df

proj_intensive_rural_cross_sector %>%
  mutate(ode_df = map(ode_proj, as.data.frame)) %>%
  pull(ode_df) %>%
  list_rbind(names_to = "run_id") %>% 
  mutate(., Complex = rep("Intensive facility and rural deer",nrow(.)), Alternative = c(rep(unique(proj_intensive_rural_cross_sector$Context)[1],length(times)*nsamples))) %>% 
  group_by(Complex,Alternative, run_id) %>% 
  summarize(meanWild = mean(I_wild), meanCaptive = mean(I_captive), 
            cumulativeWild = last(I_wild_cumulative), cumulativeCaptive = last(I_captive_cumulative)) %>% 
   merge(.,Persist_intensive_rural_cross_sector_df, by = "run_id")-> cross_sector_for_intensive_rural_complex
#> `summarise()` has grouped output by 'Complex', 'Alternative'. You can override
#> using the `.groups` argument.


cross_sector_results <- rbind(cross_sector_for_ranch_suburban_complex, cross_sector_for_ranch_rural_complex, cross_sector_for_intensive_suburban_complex, cross_sector_for_intensive_rural_complex)

cross_sector_results %>% 
   group_by(Complex) %>% 
   summarize(., prevalenceWild = median(meanWild), prevalenceCaptive = median(meanCaptive), cumulativeWild = median(cumulativeWild), cumulativeCaptive = median(cumulativeCaptive), meanPersistWild = (sum(persistWild)/nsamples), meanPersistCaptive= (sum(persistCaptive)/nsamples))
#> # A tibble: 4 × 7
#>   Complex      prevalenceWild prevalenceCaptive cumulativeWild cumulativeCaptive
#>   <chr>                 <dbl>             <dbl>          <dbl>             <dbl>
#> 1 Intensive f…     0.00000147       0.000000661      0.0000313        0.0000140 
#> 2 Intensive f…     0.00000523       0.000000682      0.000113         0.0000143 
#> 3 Ranch and r…     0.00000147       0.000000101      0.0000313        0.00000226
#> 4 Ranch and s…     0.00000522       0.000000107      0.000113         0.00000226
#> # ℹ 2 more variables: meanPersistWild <dbl>, meanPersistCaptive <dbl>

We can take these results and compare to baseline results, as well as single management alternatives.

Literature Cited

Cook, J.D., E. Rosenblatt, G.V. Direnzo, E.H.C. Grant, B.A. Mosher, F. Arce, S. Christensen, R. Ghai, M.C. Runge. In Prep. Using decision science to evaluate the risk and management of SARS-CoV-2 zoonotic transmission between humans and white-tailed deer.

Analysis referenced in Cook et al. In Prep