SARS-CoV-2 outbreak characteristics across wild and captive settings, with an initial spillover
Analysis referenced in Rosenblatt et al. In Prep
Elias Rosenblatt
Source:vignettes/SIRS_analysis_by_contexts_initialspill.Rmd
SIRS_analysis_by_contexts_initialspill.RmdIntroduction
This vignette follows the same blueprint as the SIRS_analysis_by_context
vignette. The primary difference here is that initial infected
compartment sizes are set non-zero values and human prevalence is set to
zero. This mimics an initial spillover event of a given magnitude, to
test how outbreak dynamics differ from continuous spillover detailed in
whitetailedSIRS::SIRS_analysis_by_context. These
differences are visualized in the vignette Visualize_by_context. Most of the
simulation code is suppresed in the rendered vignette, but can be viewed
in the .Rmd file. This analysis corresponds with Objective 4 of
Rosenblatt et al. In Prep.
Again, general approach with simulating SARS-CoV-2 outbreaks here is to:
Set the number of simulations
Draw a random sample of parameters used to estimate infection probabilities, contact rates, and other important parameters.
Derive parameters used the SIRS ODE equation set
Define initial SIR compartment sizes
Solve the SIRS ODE equations and calculate daily compartment sizes
Calculate risk of introduction and magnitude of spread given parameters
Calculate average daily prevalence, probability of persistence, and incidence proportion.
1. Set number of simulations to run and context of simulation
The same scenarios are considered in this vignette as detailed in the vignette simulating continuous human introduction. The major difference here is that there is an initial infection in the population and that human prevalence is set to 0%. By comparing results between an initial spillover and continual spillover, we can get a sense of how clustered exposures to a single event differ in outbreak dynamics compared to continuous spillover. Below are the scenario descriptions again.
Wild deer in rural conditions: In this scenario, wild deer are free-ranging in an area with a rural human density (3.1 humans/km2). We assumed that deer in this context were subject to regulated hunting either using still-hunting, or ground blind or treestand approaches. We also assumed that baiting and backyard feeding were illegal but may still occur. We again used Habib et al.’s (2011) proximity rate model. We estimated the duration of deer-deer proximity events, and rates and duration of human-deer proximity events using expert elicitation.
Wild deer in suburban conditions: Wild deer are free-ranging in an area of suburban human density (100 humans/km2). Deer-to-deer proximity rates were derived using the same parameters as used in the rural setting. The primary difference between this context and the wild, rural context are higher human-to-deer proximity rates and longer duration of proximity events as estimated by expert elicitation.
Outdoor captive ranch: We modeled a population of captive deer housed in an outdoor ranching facility. We assumed proximity rates between deer in this context was the same as estimated in wild contexts, with the increase of these proximity rates due to supplemental feeding. We assume that proximity rates between humans and deer in this context are the same as those estimated for the wild, suburban context, but the duration of these proximity events are longer, reflecting those typical of a captive facility.
Intensive captive facility: We considered captive deer held in a facility focused on captive breeding or exposition. Deer in this facility were predominantly indoors with high stocking densities and low indoor air exchange rates (AER; 1-hr). We estimated all proximity rates and duration of proximity using expert elicitation.
We set a seed so the results here are repeatable, and set the number of samples to the desired amount (nsamples). For this example, we ran 1000 simulations.
set.seed(23)
nsamples <- 1002. Sample parameters for simulations
Next, we draw an assortment of parameters to solve the SIRS ODE equations for a SARS-CoV-2 outbreak. If you have not read them, see the SIR model description and example vignettes to understand the workflow here.
elicitation_data <- whitetailedSIRS::draw_elicitation_samples(nsamples = nsamples)Following Rosenblatt et al. (In Prep.), four contexts were defined, including outdoor captive ranch, intensive captive facility, wild deer in rural conditions, and wild deer in suburban conditions. In this RMD file, these scenarios will be simulated independently of each other.
Outdoor captive ranch
Proximity Rates
Proximity rates were derived from a rate model developed by Habib et al. (2011), as a function of habitat conditions (habitat) and density (nwild/A_W). For all uses of the Habib et al. (2011) model in this analysis, we fixed density at 10 deer per sq km and habitat at 26% wooded habitat (classified in the calc_contact_rate function as type_contact = “med”). We also used expert elicitation estimates of the probability of direct contact for fluid transmission.
nWild <- rpois(nsamples,1000) #Abundance
A_w <- 100 #Area
habitat <- "med" #Habitat classification
sigma_season <- 1 #Season adjustment for proximity rate
epsilon_dc <- get_EE_param_vals(data = elicitation_data, my_param = "Direct Contact Probability") #Probability of direct contact between deer, given proximity.Next, we define various proximity rates for the ranch context. Object names reflect other parameters that must be included in the SIRS ODE solver (e.g. omega_ww_ranch, set to 0).
omega_ww_ranch <- rep(0, nsamples) #Deer-to-deer proximity rate in wild (set to 0; events per day).
omega_cw_null <- rep(0, nsamples) #Deer-to-deer proximity rate along fenceline (set to 0; events per day).
omega_cc_ranch <- calc_contact_rate(nsamples = nsamples, type_contact = habitat, N_w = nWild, rho_attractant = get_EE_param_vals(data = elicitation_data, my_param = "Proximity rate with baiting (17 events without baiting)")/17) #Deer-to-deer proximity rate in ranch context, mimicing wild proximity rates with the influence of baiting (events per day).
omega_hw_ranch <- rep(0, nsamples) #Human-to-deer proximity rate in wild (set to 0; events per day).
omega_hc_ranch <- get_EE_param_vals(data = elicitation_data, my_param = "Deer-Human Proximity Rate, Suburban (per 120 days)") /120 #Human-to-deer proximity rate in ranch context, mimicing suburban proximitty rates (events per day).Infection probabilities
We calculate infection probabilities for three transmission pathways - aerosol transmission from deer-to-deer, fluid transmission from deer-to-deer, and aerosol transmission from human-to-deer. Below, we will use various package functions to define parameters to estimate the infection probability given proximity or contact.
C_nu_human <- rnorm(n = nsamples, mean = 10^5.6, sd = 10^1.2)#viral load in humans (genomic copies per ml)
C_nu_deer <- 10^5.6 * get_EE_param_vals(data = elicitation_data, my_param = "Viral Load") #viral load in deer saliva, relative to humans (genomic copies per ml)
r_deer <- get_EE_param_vals(data = elicitation_data, my_param = "Dose-Response")# Dose response coefficient for deer and SARS-CoV-2
#Infection probability calculation of aerosol transmission from deer-to-deer
t_contact_deer_deer_null <- get_EE_param_vals(data = elicitation_data, my_param = "Deer Proximity Duration (minutes)") #Estimate duration of proximity event...
sigma_aero_deer_deer_ranch <- calc_sigma_aero(C_nu = C_nu_deer,
t_contact = t_contact_deer_deer_null / 60,
r = r_deer, nsamples = nsamples) #...and estimate probability of infection given that duration of proximity event.
sigma_aero_deer_deer_wild_null <- rep(0, nsamples) #Estimate infection probability in out in the wild as 0 (needs to be included for SIRS ODE equations)
#Infection probability of 0.1 ml of saliva being transferred between deer on contact
sigma_dc_deer_deer_null <- calc_sigma_dc(C_nu = C_nu_deer, nsamples = nsamples) #Calculate infection probability
#Infection probability calculation of aerosol transmission from humans-to-deer
t_contact_deer_human_ranch <- get_EE_param_vals(data = elicitation_data, my_param = "Deer-Human Proximity Duration, Captive (minutes)") #Estimate duration of human-deer proximity event in ranch facility context...
sigma_aero_deer_human_ranch <- calc_sigma_aero(ER = 0.53, C_nu = C_nu_human,
t_contact = t_contact_deer_human_ranch / 60,
r = r_deer, nsamples = nsamples)#... and calculate infection probability given the duration of a human-deer proximity event.
sigma_aero_deer_human_wild_null <- rep(0, nsamples)#Estimate human-to-deer infection probability out in the wild as 0 (needs to be included for SIRS ODE equations)Intensive captive facility
Below, we list additional parameters used for the “Intensive captive facility” scenario. Note that this scenario will use objects defined earlier in this RMarkdown.
omega_ww_intensive <- rep(0, nsamples)
omega_cc_intensive <- get_EE_param_vals(data = elicitation_data, my_param = "Deer-Deer Proximity Rate, Captive (per day)")
omega_hw_intensive <- rep(0, nsamples)
omega_hc_intensive <- get_EE_param_vals(data = elicitation_data, my_param = "Deer-Human Proximity Rate, Suburban (per 120 days)") /120
sigma_aero_deer_deer_intensive <- calc_sigma_aero(C_nu = C_nu_deer,
t_contact = t_contact_deer_deer_null / 60,
r = r_deer, nsamples = nsamples, AER = rep(1, nsamples))
t_contact_deer_human_intensive <- get_EE_param_vals(data = elicitation_data, my_param = "Deer-Human Proximity Duration, Captive (minutes)")
sigma_aero_deer_human_intensive <- calc_sigma_aero(ER = 0.53, C_nu = C_nu_human,
t_contact = t_contact_deer_human_intensive / 60,
r = r_deer, nsamples = nsamples, AER = rep(1, nsamples))Wild deer in rural conditions
Again, we list additional parameters used for the “Wild deer in rural conditions” scenario. Note that this scenario will use objects defined earlier in this RMarkdown.
omega_ww_rural <- calc_contact_rate(nsamples = nsamples, type_contact = habitat, N_w = nWild)
omega_cc_rural <- rep(0, nsamples)
omega_hw_rural <- get_EE_param_vals(data = elicitation_data, my_param = "Deer-Human Proximity Rate, Rural (per 120 days)") /120
omega_hc_rural <- rep(0, nsamples)
t_contact_deer_human_rural <- get_EE_param_vals(data = elicitation_data, my_param = "Deer-Human Proximity Duration, Rural (minutes)")
sigma_aero_deer_deer_rural <- calc_sigma_aero(C_nu = C_nu_deer,
t_contact = t_contact_deer_deer_null / 60,
r = r_deer, nsamples = nsamples, AER = rep(4, nsamples))
sigma_aero_deer_deer_captive_null <- rep(0, nsamples)
sigma_aero_deer_human_rural <- calc_sigma_aero(ER = 0.53, C_nu = C_nu_human,
t_contact = t_contact_deer_human_rural / 60,
r = r_deer, nsamples = nsamples, AER = rep(4, nsamples))
sigma_aero_deer_human_captive_null <- rep(0, nsamples)Wild deer in suburban conditions
Again, we list additional parameters used for the “Wild deer in suburban conditions” scenario. Note that this scenario will use objects defined earlier in this RMarkdown.
omega_ww_suburban <- calc_contact_rate(nsamples = nsamples, type_contact = habitat, N_w = nWild)
omega_cc_suburban <- rep(0, nsamples)
omega_hw_suburban <- get_EE_param_vals(data = elicitation_data, my_param = "Deer-Human Proximity Rate, Suburban (per 120 days)") /120
omega_hc_suburban <- rep(0, nsamples)
t_contact_deer_human_suburban <- get_EE_param_vals(data = elicitation_data, my_param = "Deer-Human Proximity Duration, Suburban (minutes)")
sigma_aero_deer_deer_suburban <- calc_sigma_aero(C_nu = C_nu_deer,
t_contact = t_contact_deer_deer_null / 60,
r = r_deer, nsamples = nsamples, AER = rep(4, nsamples))
sigma_aero_deer_human_suburban <- calc_sigma_aero(ER = 0.53, C_nu = C_nu_human,
t_contact = t_contact_deer_human_suburban / 60,
r = r_deer, nsamples = nsamples, AER = rep(4, nsamples))Recovery and Immunity Loss Rates
For all simulations, animals recover from their infection at a given rate (gamma_recov), and lose their immunity and become susceptible at a given rate (alpha_immunity). We use a recovery rate based on a 6 day infectious period (Palmer et al. 2021), and immunity loss rate based on expert elicited estimates. The recovery rate is repeated across all iterations of a simulation (nsamples). The immunity loss rate is derived from [nsamples] draws from an error distribution estimated by the expert elicitation process. Note that human prevalence is fixed at 0.
gamma_recov <- rep(1/6, nsamples)
alpha_immunity_null <- 1 / get_EE_param_vals(data = elicitation_data, my_param = 'Temporary Immunity')
I_human_null <- 03. Derive parameters used in the SIRS ODE euqation set.
Next, we create derived parameters for each scenario:
ranch.params <- alternative(alpha_immunity = alpha_immunity_null,
omega_ww = omega_ww_ranch, omega_cw = omega_cw_null,omega_cc = omega_cc_ranch,
omega_hw = omega_hw_ranch, omega_hc = omega_hc_ranch,
sigma_aero_deer_deer_wild = sigma_aero_deer_deer_wild_null, sigma_aero_deer_deer_captive = sigma_aero_deer_deer_ranch, sigma_aero_deer_human_wild = sigma_aero_deer_human_wild_null, sigma_aero_deer_human_capt = sigma_aero_deer_human_ranch, epsilon_dc = epsilon_dc, sigma_dc_deer_deer = sigma_dc_deer_deer_null, gamma_recov = gamma_recov, I_human = rep(I_human_null, nsamples), boost = rep(0, nsamples))
intensive.params <- alternative(
omega_ww = omega_ww_intensive, omega_cw = omega_cw_null, omega_cc = omega_cc_intensive,
omega_hw = omega_hw_intensive, omega_hc = omega_hc_intensive,
sigma_aero_deer_deer_wild = sigma_aero_deer_human_wild_null, sigma_aero_deer_deer_captive = sigma_aero_deer_deer_intensive, sigma_aero_deer_human_wild = sigma_aero_deer_human_wild_null, sigma_aero_deer_human_capt = sigma_aero_deer_human_intensive, sigma_dc_deer_deer = sigma_dc_deer_deer_null,
alpha_immunity = alpha_immunity_null, epsilon_dc = epsilon_dc, gamma_recov = gamma_recov, I_human = rep(I_human_null, nsamples), boost = rep(0, nsamples))
rural.params <- alternative(omega_ww = omega_ww_rural, omega_cw = omega_cw_null, omega_cc = omega_cc_rural, omega_hw = omega_hw_rural, omega_hc = omega_hc_rural, sigma_aero_deer_deer_wild = sigma_aero_deer_deer_rural, sigma_aero_deer_deer_captive = sigma_aero_deer_deer_captive_null, sigma_aero_deer_human_wild = sigma_aero_deer_human_rural, sigma_aero_deer_human_capt = sigma_aero_deer_human_captive_null, sigma_dc_deer_deer = sigma_dc_deer_deer_null,alpha_immunity = alpha_immunity_null, epsilon_dc = epsilon_dc, gamma_recov = gamma_recov, I_human = rep(I_human_null, nsamples), boost = rep(0, nsamples))
suburban.params <- alternative(omega_ww = omega_ww_suburban, omega_cw = omega_cw_null, omega_cc = omega_cc_suburban, omega_hw = omega_hw_suburban, omega_hc = omega_hc_suburban, sigma_aero_deer_deer_wild = sigma_aero_deer_deer_suburban, sigma_aero_deer_deer_captive = sigma_aero_deer_deer_captive_null, sigma_aero_deer_human_wild = sigma_aero_deer_human_suburban, sigma_aero_deer_human_capt = sigma_aero_deer_human_captive_null, sigma_dc_deer_deer = sigma_dc_deer_deer_null,alpha_immunity = alpha_immunity_null, epsilon_dc = epsilon_dc, gamma_recov = gamma_recov, I_human = rep(I_human_null, nsamples), boost = rep(0, nsamples))4. Define initial SIR compartment sizes
Next, we define three sets of compartment sizes for both captive and wild scenarios that differ in the size of the initial exposure. These range from a widespread exposure (1 in 1000 deer), a minimal exposure (1 in 1 million), and a extremely rare exposure (1 in 1 billion).
captive.inits.fall.1in1000 <- initial_compartments(S_wild_prop = 0, S_captive_prop = 0.999, I_captive_prop = 0.001, draws = nsamples)
captive.inits.steady.1in1000 <- initial_compartments(S_wild_prop = 0, S_captive_prop = 0.999, I_captive_prop = 0.001, draws = nsamples, steady = TRUE)
wild.inits.fall.1in1000 <- initial_compartments(S_captive_prop = 0, S_wild_prop = 0.999, I_wild_prop = 0.001, draws = nsamples)
wild.inits.steady.1in1000 <- initial_compartments(S_captive_prop = 0, S_wild_prop = 0.999, I_wild_prop = 0.001, draws = nsamples, steady = TRUE)
captive.inits.fall.1in1mil <- initial_compartments(S_wild_prop = 0, S_captive_prop = 0.999999, I_captive_prop = 0.000001, draws = nsamples)
captive.inits.steady.1in1mil <- initial_compartments(S_wild_prop = 0, S_captive_prop = 0.999999, I_captive_prop = 0.000001, draws = nsamples, steady = TRUE)
wild.inits.fall.1in1mil<- initial_compartments(S_captive_prop = 0, S_wild_prop = 0.999999, I_wild_prop = 0.000001, draws = nsamples)
wild.inits.steady.1in1mil <- initial_compartments(S_captive_prop = 0, S_wild_prop = 0.999999, I_wild_prop = 0.000001, draws = nsamples, steady = TRUE)
captive.inits.fall.1in1bil <- initial_compartments(S_wild_prop = 0, S_captive_prop = 0.999999999, I_captive_prop = 0.000000001, draws = nsamples)
captive.inits.steady.1in1bil <- initial_compartments(S_wild_prop = 0, S_captive_prop = 0.999999999, I_captive_prop = 0.000000001, draws = nsamples, steady = TRUE)
wild.inits.fall.1in1bil <- initial_compartments(S_captive_prop = 0, S_wild_prop = 0.999999999, I_wild_prop = 0.000000001, draws = nsamples)
wild.inits.steady.1in1bil <- initial_compartments(S_captive_prop = 0, S_wild_prop = 0.999999999, I_wild_prop = 0.000000001, draws = nsamples, steady = TRUE)5. Solve the SIRS ODE equations and calculate daily compartment sizes
We solve the SIRS ODE equations based on various initial exposures.
times <- seq(0, 120, by = 1)
#1 in 1000 deer
proj.ranch.1in1000 <- run(iter = nsamples, initial_compartments = captive.inits.fall.1in1000, initial_compartments_steady = captive.inits.steady.1in1000, params = ranch.params, times = times, name = "Outdoor ranch")
proj.intensive.1in1000 <- run(iter = nsamples, initial_compartments = captive.inits.fall.1in1000, initial_compartments_steady = captive.inits.steady.1in1000, params = intensive.params, times = times, name = "Intensive facility")
proj.rural.1in1000 <- run(iter = nsamples, initial_compartments = wild.inits.fall.1in1000, initial_compartments_steady = wild.inits.steady.1in1000 ,params = rural.params, times = times, name = "Wild, rural")
proj.suburban.1in1000 <- run(iter = nsamples, initial_compartments = wild.inits.fall.1in1000, initial_compartments_steady = wild.inits.steady.1in1000, params = suburban.params, times = times, name = "Wild, suburban")
sirs_results_contexts.1in1000 <- rbind(proj.ranch.1in1000,proj.intensive.1in1000, proj.rural.1in1000, proj.suburban.1in1000)
#1 in 1 million deer
proj.ranch.1in1mil <- run(iter = nsamples, initial_compartments = captive.inits.fall.1in1mil, initial_compartments_steady = captive.inits.steady.1in1mil, params = ranch.params, times = times, name = "Outdoor ranch")
proj.intensive.1in1mil <- run(iter = nsamples, initial_compartments = captive.inits.fall.1in1mil, initial_compartments_steady = captive.inits.steady.1in1mil, params = intensive.params, times = times, name = "Intensive facility")
proj.rural.1in1mil <- run(iter = nsamples, initial_compartments = wild.inits.fall.1in1mil, initial_compartments_steady = wild.inits.steady.1in1mil ,params = rural.params, times = times, name = "Wild, rural")
proj.suburban.1in1mil <- run(iter = nsamples, initial_compartments = wild.inits.fall.1in1mil, initial_compartments_steady = wild.inits.steady.1in1mil, params = suburban.params, times = times, name = "Wild, suburban")
#1 in 1 billion deer
proj.ranch.1in1bil <- run(iter = nsamples, initial_compartments = captive.inits.fall.1in1bil, initial_compartments_steady = captive.inits.steady.1in1bil, params = ranch.params, times = times, name = "Outdoor ranch")
proj.intensive.1in1bil <- run(iter = nsamples, initial_compartments = captive.inits.fall.1in1bil, initial_compartments_steady = captive.inits.steady.1in1bil, params = intensive.params, times = times, name = "Intensive facility")
proj.rural.1in1bil <- run(iter = nsamples, initial_compartments = wild.inits.fall.1in1bil, initial_compartments_steady = wild.inits.steady.1in1bil ,params = rural.params, times = times, name = "Wild, rural")
proj.suburban.1in1bil <- run(iter = nsamples, initial_compartments = wild.inits.fall.1in1bil, initial_compartments_steady = wild.inits.steady.1in1bil, params = suburban.params, times = times, name = "Wild, suburban")6 and 7: Calculate risk of introduction and magnitude of spread given parameters, average daily prevalence, probability of persistence, and incidence proportion.
We are going to calculate several metrics for each iteration to summarize across contexts. We combine the projections for each context into one list object, and summarize the average prevalence (Infected) and incidence proportion (Cumulative_infections). We
We save these resulting summary data frame for future use.
1 in 1000 deer initially infected
First, we record prevelence and incidence proportion:
#Next, we record SARS-CoV-2 persistence:
We then calculate \(R_0\) and FOI for each scenario, and combine with outbreak dynamics above to create one object for export.
1 in 1 million deer initially infected
First, we record prevelence and incidence proportion:
#Next, we record SARS-CoV-2 persistence:
We then calculate \(R_0\) and FOI for each scenario, and combine with outbreak dynamics above to create one object for export.
1 in 1 billion deer initially infected
First, we record prevelence and incidence proportion:
#Next, we record SARS-CoV-2 persistence:
We then calculate \(R_0\) and FOI for each scenario, and combine with outbreak dynamics above to create one object for export.
Results from these simulations are stored as
whitetailedSIRS::initial_infection_results_1_in_1000,
whitetailedSIRS::initial_infection_results_1_in_1mil, and
whitetailedSIRS::initial_infection_results_1_in_1bil.
initial_infection_results_1_in_1000 <- df_initial_infections_1_in_1000
initial_infection_results_1_in_1mil <- df_initial_infections_1_in_1mil
initial_infection_results_1_in_1bil <- df_initial_infections_1_in_1bilNext Steps:
This vignette detailed the simulations run in Rosenblatt et al. (In Prep) to simulate outbreaks from an initial exposure event. The next vignette produce the figures published in Rosenblatt et al. In Prep..
Click here to return to the vignette simulating outbreaks from continous exposure to infectious humans.
Literature Cited
- Rosenblatt, E., J.D. Cook, G.V. Direnzo, E.H.C. Grant, F. Arce, K. Pepin, F.J. Rudolph, M.C. Runge, S. Shriner, D. Walsh, B.A. Mosher. In Prep. Epidemiological modeling of SARS-CoV-2 in white-tailed deer (Odocoileus virginianus) reveals conditions for introduction and widespread transmission.