SPARSE-MOD SEIR Model

JR Mihaljevic

December 2020

The SEIR Example Model

Here we present a walk-through of using the SPARSE-MOD SEIR Model. SEIR stands for the numbers (not fractions) of Susceptible-Exposed-Infectious-Removed hosts, which define the state variables of the model. The differential equations of this model are: \[\begin{align} \frac{dS}{dt} &= \overbrace{\mu N}^{\text{Reproduction}} - \overbrace{\hat{\beta} S I}^{\text{Transmission}} - \overbrace{\mu S}^{\text{Natural Mortality}} \\ \frac{dE}{dt} &= \hat{\beta} S I - \overbrace{(\delta + \mu) E}^{\text{Latency and Mortality}} \\ \frac{dI}{dt} &= \delta E - \overbrace{(\gamma + \mu) I}^{\text{Recovery and Mortality}} \\ \frac{dR}{dt} &= \gamma I - \mu R \end{align}\]

In this classical model, Susceptible individuals become exposed to the pathogen, moving to the Exposed class. There is a period of latency or incubation (\(1/\delta\)) before the hosts become Infectious. Infectious hosts can then recover from the pathogen at an average rate of \(1/\gamma\), in which case we assume the Recovered individuals are immune to the pathogen for life. The model includes host demography, in which we assume that birth rate is equal to death rate which is equal to \(\mu\). This ensures that the local population size remains constant through time (not accounting for temporary migration events in the SPARSEMODr framework).1 The model also assumes that all host classes contribute to reproduction and that offspring are fully susceptible. Finally, the model assumes that mortaltity is natural and not pathogen-induced. In other words, the pathogen is not meaningfully virulent in terms of host life-expectancy. This model has been typically used for childhood diseases that confer life-long immunity. Note thta we also use the notation \(\hat{\beta}\) to emphasize that the transmission rate can be modeled as frequency- or density-dependent.2

In this example, we will demonstrate how we can use the time-windows feature to implement seasonality in transmission rates that can generate sustained oscillations in pathogen prevalence in this SEIR model with host demography.

Generating a synthetic meta-population

First, we will generate some data that describes the meta-population3 of interest. Note that this is essentially the same set-up as compared to the SPARSEMODr COVID-19 model vignette.

# Set seed for reproducibility
set.seed(2)

# Number of focal populations:
n_pop = 15

# Population sizes + areas
## Draw from neg binom:
pop_N = rnbinom(n_pop, mu = 50000, size = 1)
census_area = rnbinom(n_pop, mu = 50, size = 3)

# Identification variable for later
pop_ID = c(1:n_pop)

# Assign coordinates, plot for reference
lat_temp = runif(n_pop, 32, 37)
long_temp = runif(n_pop, -114, -109)

pop_local_df =
  data.frame(pop_ID = pop_ID,
             pop_N = pop_N,
             census_area,
             lat = lat_temp,
             long = long_temp) %>%
  # Assign regions by quadrant
  ## Used later for aggregation
  mutate(region = case_when(
    lat >= 34.5 & long <= -111.5 ~ "1",
    lat >= 34.5 & long >  -111.5 ~ "2",
    lat <  34.5 & long >  -111.5 ~ "3",
    lat <  34.5 & long <= -111.5 ~ "4"
  ))

# Plot the map:
ggplot(pop_local_df) +
  geom_point(aes(x = long, y = lat, color = region),
             shape = 19) +
  scale_color_viridis_d(direction = -1) +
  # Map coord
  coord_quickmap() +
  theme_classic() +
  theme(
    axis.line = element_blank(),
    axis.title = element_blank(),
    plot.margin = unit(c(0, 0.1, 0, 0), "cm")
  )

# Calculate pairwise dist
## in meters so divide by 1000 for km
dist_mat = geosphere::distm(cbind(pop_local_df$long, pop_local_df$lat))/1000
hist(dist_mat, xlab = "Distance (km)", main = "")

# We need to determine how many Exposed individuals
# are present at the start in each population
E_pops = vector("numeric", length = n_pop)
# We'll assume a total number of exposed across the
# full meta-community, and then randomly distribute these hosts
n_initial_E = 10
# (more exposed in larger populations)
these_E <- sample.int(n_pop,
                      size = n_initial_E,
                      replace = TRUE,
                      prob = pop_N)
for(i in 1:n_initial_E){
  E_pops[these_E[i]] <- E_pops[these_E[i]] + 1
}

pop_local_df$E_pops = E_pops

Setting up the \(\texttt{time_windows}\) object

One of the benefits of the SPARSEMODr design is that the user can specify how certain parameters of the model change over time (see our vignette on key features of SPARSEMODr for more details). We demonstrate this below. In this particular example, we allow the time-varying R0 to change daily in out model, which forces \(\hat{\beta}\) to fluctuate daily. In this particular example, we apply sinusoidal-forcing, allowing R0 to peak in the fall-winter months and wane in the spring-summer months. We’ll use this particular forcing equation; \[R_{0}(t) = \bar{R}_{0} (1 + cos(\frac{2 \pi t}{t_{\text{mode}}}))\] In this case we have an average R0 and \(t_{\text{mode}}\) controls the number of R0 cycles per year.

We’ll use the \(\texttt{time_windows()}\) function to generate a pattern of time-varying R0 that looks like the following:

We’ll change the R0, but leave the migration parameters constant.

Running the spatial SEIR model in parallel

Now we have all of the input elements needed to run SPARSEMODr’s SEIR model. Below we demonstrate a workflow to generate stochastic realizations of the model in parallel. Notice that in this host species, breeding (and mortality) occurs on average once every two years.

# How many realizations of the model?
n_realz = 30

# Need to assign a distinct seed for each realization
## Allows for reproducibility
input_realz_seeds = c(1:n_realz)

# Run the model in parallel

model_output =
  model_parallel(
      # Necessary inputs
      input_dist_mat = dist_mat,
      input_census_area = pop_local_df$census_area,
      input_tw = tw,
      input_realz_seeds = input_realz_seeds,
      # OTHER MODEL PARAMS
      trans_type = 1, # freq-dependent trans
      stoch_sd = 0.75,  # stoch transmission sd,
      control = seir_control    # data structure with seir-specific params
  )

glimpse(model_output)
#> Rows: 1,642,500
#> Columns: 12
#> $ pops.seed         <int> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, …
#> $ pops.pop          <int> 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, …
#> $ pops.time         <int> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 2, 2, …
#> $ pops.S_pop        <int> 3289, 1011, 6182, 5501, 71650, 13455, 80433, 9139, …
#> $ pops.E_pop        <int> 0, 0, 0, 0, 2, 0, 1, 0, 0, 2, 0, 4, 0, 0, 1, 0, 0, …
#> $ pops.I_pop        <int> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
#> $ pops.R_pop        <int> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
#> $ events.birth      <int> 4, 2, 14, 9, 81, 13, 108, 16, 0, 17, 144, 242, 94, …
#> $ events.exposed    <int> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
#> $ events.infectious <int> 0, 0, 0, 0, 0, 0, 0, 0, 0, 2, 0, 2, 0, 0, 1, 0, 0, …
#> $ events.recov      <int> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
#> $ events.death      <int> 2, 2, 8, 2, 94, 19, 118, 11, 0, 17, 137, 189, 125, …

Plotting the output

First we need to manipulate and aggregate the output data. Here we show an example to plot the number of infectious individuals in the populations over time.

Now we’ll create a simple figure to visualize the stochastic realizations over time.

  1. See our vignette on key features of SPARSEMODr, including migration dynamics for more general details of the SPARSEMODr package.↩

  2. See our vignette on key features of SPARSEMODr, including migration dynamics for more general details of the SPARSEMODr package.↩

  3. A set of distinct, focal populations that are connected by migration↩