Seasonality

Formulating seasonality in compartmental model

Disclaimer

The content of this blog is mostly based on Chapter 5 - Temporally Forced Model of the book “Modeling Infectious Diseases in Human and Animals” by Keeling and Rohani (2008).

The intended purpose of this blog is to summarize the content and provide reproducible code/model implementations in R. For a more detailed explanation on the topic, please refer to the source material.

1 Background

Traditional SIR model assumes that the force of infection depends on the proportion of infectives in the population, formulated as \(\beta \frac{I}{N}\) in the following system of ODEs.

\[ \begin{cases} \frac{dS}{dt} = bN-\beta \frac{I}{N} S\\ \frac{dI}{dt} = \beta \frac{I}{N} S - \gamma I\\ \frac{dR}{dt} = \gamma I\\ \end{cases} \]

Where:

• \(b\) is the per capita birth rate

• \(\beta\) is the transmission rate

• \(\gamma\) is the recovery rate

However, several statistical approaches have revealed that transmission of childhood infections also varies seasonally (Anderson and May 1991). This property can be reflected by using a time varying contact rate \(\beta(t)\).

Different functional forms for \(\beta(t)\) can lead to different disease dynamics. In the following sections, we first focus on sinusoidal \(\beta(t)\) then discuss some alternative forms for \(\beta(t)\).

2 Sinusoidal \(\beta(t)\)

Sinusoidal functions (sine, cosine, etc.) are the most straightforward choices to model \(\beta(t)\) due to their periodic nature. Bailey (1975) formulate \(\beta(t)\) as followed

\[ \beta(t) = \beta_0(1+\beta_1 \text{cos}(\omega t)) \]

Where:

• \(\beta_0\) is the baseline transmission rate (average transmission rate).

• \(\beta_1\) controls the amplitude of seasonality (amplitude = \(\beta_1 *\beta_0\))

• \(\omega\) controls the period of forcing \(\frac{2\pi}{\text{period in time unit}}\) (when period of forcing is 1 year, time unit is in week, \(\omega = \frac{2\pi}{52}\))

How seasonal forcing can amplify small perturbations to the unforced equilibrium

Oscillations in number of infectives under small perturbations, assuming a small amplitude of seasonality have been investigated.

While period of oscillations stays the same (as before seasonal forcing), amplitude is now given by the following equation.

\[ M = \beta_1 \omega \gamma \{ (b* \beta_0 - \omega^2 )^2 + (\omega b R_0)^2 \}^{-\frac{1}{2}} \tag{1}\]

To demonstrate what this function implies, we can try plugging in some values. For measles, the parameters used are \(\frac{1}{\gamma} = 2\), \(b R_0 \sim 0.014\), \(\omega = \frac{\pi}{26}\), which would result in \(M \sim 7.76 \beta_1\). The interpretation for this is: a 10% increase in \(\beta_1\) can lead to 78% increase in case notifications.

Method to derive the amplitude equation

Small perturbations to the unforced equilibrium is modeled by the substitutions \(S = S(1+s)\) and \(I = I(1+i)\) (where \(s\) and \(i\) are small perturbations on the susceptible and infectious population).

The assumption of a small amplitude of seasonality is equivalent to having a small \(\beta_1\) value (i.e. \(\beta_1 << 1\)).

After omitting some intermediate steps, we can obtain the second order differential equation in the small infectious perturbation \(i\).

\[ \frac{d^2 i}{dt^2} + b R_0 \frac{di}{dt} + b \beta_0 i = -\beta_1 \omega \gamma \text{ sin}(\omega t) \]

The integral for this equation gives the period and amplitude of oscillations as driven by the seasonal term.

2.1 Harmonic and Sub-harmonic resonance

Let’s first revisit the familiar SIR model, but now we put our focus on \(S\) and \(I\) compartments

\[ \begin{cases} \frac{dS}{dt} = bN - (\beta(t) \frac{I}{N} + b)S \\ \frac{dI}{dt} = \beta(t) \frac{I}{N} S - (b + \gamma)I \end{cases} \]Note that the death and birth rate are assumed to be the same to keep the total population constant.

\(\beta(t)\) is defined as \(\beta(t) = \beta_0(1- \beta_1 \text{cos}(\omega t))\). Note that minus sign is used so seasonal forcing is at 0 at the start.

In the absence of seasonal forcing (i.e. \(\beta_1 = 0\)), the system fluctuates with frequency \(F\), where.

\[ F^2 = b(\gamma + b)(R_0 - 1) - (\frac{b R_0}{2})^2 \]

These oscillations are hereby referred to as the natural oscillations.

When the natural period of oscillations is approximately the same as seasonal forcing period (i.e., \(F \approx \omega\) ) we observe harmonic resonance, where the amplitude of oscillations may be greatly increased.

Looking at Equation 1, when forcing is relatively small, the amplification of sinusoidal forcing ( \(\frac{M}{\beta_1}\) ) is largest when the forcing frequency \(\omega\) is close the natural oscillation frequency \(F\). (i.e. \(F \approx \sqrt{b \beta_0}\) ).

When natural period of oscillations \(\frac{1}{F}\) is close to an integer multiple of period of forcing \(\frac{1}{\omega}\) (i.e. \(\frac{1}{F} \approx n \frac{1}{\omega}\)), we observe a phenomenon called sub-harmonic resonance. Under this phenomenon, the observed period of oscilations is longer than that of forcing.

Changes in \(\beta_1\) and \(R_0\) can also lead to qualitatively different epidemic patterns.

📌 Introducing time-dependent contact rate generate a variety of dynamical patterns - depending on parameter values - from simple annual epidemics, to multi-annual outbreaks and eventually chaos.

Example: base model (model without seasonality forcing) vs model with harmonic resonance vs model with sub-harmonic resonance

library(odin2) 
library(dust2)

sir_seasonality <- odin2::odin({
  N <- S + I + R
  deriv(S) <- b*N - (beta_t*(I/N) + b)*S
  deriv(I) <-beta_t*(I/N)*S - (b + gamma)*I
  deriv(R) <- gamma*I - R*b

  # seasonality forcing
  beta_0 <- parameter(0.4)
  beta_1 <- parameter(0)
  omega <- parameter(2*3.14/52) # use week as time unit by default
  beta_t <- beta_0*(1 - beta_1*cos(omega*time))
  
  # initialize starting population
  init_S <- parameter(9500)
  init_I <- parameter(500)
  gamma <- parameter(0.05)
  b <- parameter(0.05)
  
  initial(S) <- init_S
  initial(I) <- init_I
  initial(R) <- 0
  
})

library(tidyverse)

# helper function to compute F 
get_f_val <- function(params, R_0=NULL){
  f_val <- NULL
  # --- Compute F -----
  with(params, {
    R_0 <- (beta_0)/(gamma + b)
    f_val <<- sqrt( b*(gamma + b)*(R_0 - 1) - (b*R_0/2)**2 )
  })
  
  f_val
}

# helper function to run the model 
run_mod <- function(mod, pars, duration=100, timestep=0.05){
  # --- initialize simulation time ---- 
  times <- seq(0, duration, timestep)
  
  sys <- dust_system_create(mod, pars)
  dust_system_set_state_initial(sys)
  out <- dust_system_simulate(sys, times)
  out <- dust_unpack_state(sys, out)
  
  out %>% 
    as.data.frame() %>% 
    mutate(
      t = times
    )
}

# helper function to plot the proportion of a specified compartment 
plot_comp_prop <- function(data, comp="I", t_lower = 90, use_log = TRUE, osc_freq=NULL){
  plot <- data %>% 
    mutate(
      N = S + I + R,
      prop = .data[[comp]]/N,
      prop = if(use_log) log(prop) else prop
    ) %>% 
    filter(t >= t_lower) %>%
    ggplot(aes(x = t, y = prop)) +
      geom_line(color = "cornflowerblue") 
  
  # --- annotate start and end of each cycle if oscillation frequency is provided --- 
  if(!is.null(osc_freq)){
    # compute oscillation period
    osc_period <- 2*pi/osc_freq
    
    cycle_df <- data.frame(
        cycle_t = seq(0, by = osc_period, 
                      length.out = round(max(data$t)/osc_period) + 1)
      ) %>% 
      filter(cycle_t >= t_lower)
    
    plot <- plot + geom_vline(
      aes(xintercept = cycle_t),
      data=cycle_df,
      color = "red", linetype = "dashed")
  }
  
  plot + 
    labs(
        title = paste0("Oscillations in ",comp, " compartment"),
        y = if(use_log) paste0("Log(proportion of ",comp, ")") else 
          paste0("proportion of ", comp),
        x = "Time (years)"
      )
}

# helper function to plot the SIR model
plot_sir <- function(data, ylab="Count"){
  data %>% 
    ggplot() +
      geom_line(aes(x = t, y = S, color = "S")) +
      geom_line(aes(x = t, y = I, color = "I")) +
      geom_line(aes(x = t, y = R, color = "R")) +
      scale_color_manual(
        values = c("S" = "cornflowerblue", "I" = "red", "R" = "green")
      )+
      labs(
        title = "Model output",
        x = "Time",
        y = ylab
      )
}

For this example, the following set of parameters is used, note that 1 time unit is equivalent to 1 year

• \(\gamma = \frac{365}{13}\) for 13 days infectious period

• \(b= 0.02\)

• \(\beta_0 = 478\)

• For initial values \(S(0) = 0.06\) and \(I(0) = 0.001\)

• Value for \(\omega\) will be changed to demonstrate resonance and harmonic resonance

The natural frequency of oscillations under this set of parameters is \(F \approx 1.7\) (i.e. the period of natural oscillations is \(\frac{2\pi}{F} \approx 2.098 \text{ (years)}\)). To visualize this, we can try running and plotting the model when there is no forcing ( \(\beta_1 = 0\) ).

no_forcing_pars <- list(
  beta_1 = 0,
  beta_0 = 478,
  gamma = 365/13,
  b = 0.02,
  init_S = 6e-2,
  init_I = 1e-3
)

# compute natural oscillation frequency  
f_val <- get_f_val(no_forcing_pars)

# --- run model + plot proportion of infectives ----
run_mod(sir_seasonality, no_forcing_pars, duration=150) %>% 
  plot_comp_prop(t_lower=140)

Figure 1: Oscillations in I when there is no forcing

To demonstrate how sub-harmonic resonance would affect the proportion of infectives, we can set \(\beta_1 > 0\) ( \(\beta_1 = 0.1\) in this example) to apply seasonal forcing and set \(\omega = F\).

harmonic_pars <- no_forcing_pars
harmonic_pars$beta_1 <- 0.1

# compute natural oscillation frequency  
f_val <- get_f_val(harmonic_pars)
# set omega = F
harmonic_pars$omega <- f_val

# --- run model + plot proportion of infectives ----
run_mod(sir_seasonality, harmonic_pars, duration=150) %>%
  plot_comp_prop(t_lower=140)

Figure 2: Oscillations in I when period of forcing is the same as natural oscillation period

To demonstrate how sub-harmonic resonance would affect the proportion of infectives, we can set \(\beta_1 > 0\) ( \(\beta_1 = 0.1\) in this example) to apply seasonal forcing and set \(\omega = nF\) where \(n\) is an integer greater than 1.

subharmonic_pars <- no_forcing_pars
subharmonic_pars$beta_1 <- 0.1

# compute natural oscillation frequency  
f_val <- get_f_val(subharmonic_pars)
subharmonic_pars$omega <- 2*f_val

# --- run model + plot proportion of infectives ----
run_mod(sir_seasonality, subharmonic_pars, duration=150) %>% 
  plot_comp_prop(t_lower = 140)

Figure 3: Oscillations in I when the natural oscillation period is a multiple of period of forcing

2.2 The anatomy of seasonally forced epidemics

We now take a closer look at the oscillations in the \(I\) compartment by re-visting the condition for an epidemic to grow.

Recall that when the disease incidence increases, \(\frac{dI}{dt} > 0\). After omitting a few steps, we derive the following condition.

\[ \frac{S(t)}{N} > \frac{\gamma + b}{\beta(t)} \]

Deriving the condition

\[ \frac{dI}{dt} = \beta(t) \frac{I}{N} S - (b + \gamma)I > 0 \]

Factoring by \(I\), we get

\[ I(\beta(t) \frac{1}{N} S - b - \gamma) > 0 \]

Since \(I(t) \geq 0\), the condition is true when

\[ \beta(t) \frac{S}{N} - b - \gamma > 0 \]

\[ \Rightarrow \frac{S}{N} > \frac{b + \gamma}{\beta(t)} \]

In other words, \(I(t)\) only increases when the proportion of susceptible is above the threshold \(\frac{\gamma + b}{\beta(t)}\).

This plot below visualizes the observed infectives (log-transformed for clarity) and the proportion of susceptible with the threshold \(\frac{\gamma + b}{\beta(t)}\). Note that both the plot for infectives and susceptible are color coded by whether \(\frac{S(t)}{N} > \frac{b + \gamma}{\beta(t)}\). The parameters are taken from the previous example Figure 2.

library(patchwork)

# compute model output df with threshold 
harmonic_out <- run_mod(sir_seasonality, harmonic_pars, duration=150) %>% 
  filter(t>=140) %>% 
  mutate(
    # compute beta_t and threshold for the filtered timeframe
    beta_t = harmonic_pars$beta_0*(
      1 - harmonic_pars$beta_1*cos(harmonic_pars$omega*t)
    ),
    # threshold for disease to spread
    threshold = (harmonic_pars$gamma + harmonic_pars$b)/beta_t, 
    # check whether prop of S is above threshold
    above_threshold = (S/(S + I + R)) > threshold
  )

# plot for I compartment 
i_prop_plot <- harmonic_out %>% 
  plot_comp_prop(comp = "I", t_lower = 140, use_log = TRUE) +
  geom_point(aes(x = t, y = prop, color = above_threshold), size = 0.8) 

# plot for S compartment compartment with threshold
s_prop_plot <- harmonic_out %>% 
  plot_comp_prop(comp="S", t_lower = 140, use_log = FALSE) +
  geom_point(aes(x = t, y = prop, color = above_threshold), size = 0.8) +
  geom_line(aes(x = t, y = threshold), 
            color = "black", linetype = "dashed")

i_prop_plot/s_prop_plot

We can also generate a similar plot using the parameters from Figure 3.

library(patchwork)

# compute model output df with threshold 
subharmonic_out <- run_mod(sir_seasonality, subharmonic_pars, duration=150) %>% 
  filter(t>=140) %>% 
  mutate(
    # compute beta_t and threshold for the filtered timeframe
    beta_t = subharmonic_pars$beta_0*(
      1 - subharmonic_pars$beta_1*cos(subharmonic_pars$omega*t)
    ),
    # threshold for disease to spread
    threshold = subharmonic_pars$gamma/beta_t, 
    # check whether prop of S is above threshold
    above_threshold = (S/(S + I + R)) > threshold
  )

# plot for I compartment 
i_prop_plot <- subharmonic_out %>% 
  plot_comp_prop(comp = "I", t_lower = 140, use_log = TRUE) +
  geom_point(aes(x = t, y = prop, color = above_threshold), size = 0.8) 

# plot for S compartment compartment with threshold
s_prop_plot <- subharmonic_out %>% 
  plot_comp_prop(comp="S", t_lower = 140, use_log = FALSE) +
  geom_point(aes(x = t, y = prop, color = above_threshold), size = 0.8) +
  geom_line(aes(x = t, y = threshold), 
            color = "black", linetype = "dashed")

i_prop_plot/s_prop_plot

2.3 Bifurcation diagram

//TBD:

3 Other functional forms for \(\beta(t)\)

Formulation for \(\beta(t)\) may varies to better reflect empirical data or prior knowledge. For example, transmission rate for childhood diseases is assumed to be high during school terms and low otherwise, this assumption is formulated as followed.

\[ \beta(t) = \frac{\beta_0}{\frac{1}{365} ((1 + \beta_1)D_+ + (1 - \beta_1)D_-)} (1+\beta_1 \text{Term}(t)) \]

Where:

• \(\text{Term}(t)\) returns \(+1\) when it is school term and \(-1\) otherwise.

• \(\beta_1\) is the amplitude of seasonality forcing

• \(\beta_0\) is the mean transmission rate

• a factor of \(\frac{1}{\frac{1}{365} ((1 + \beta_1)D_+ + (1 - \beta_1)D_-)}\) is multiplied to \(\beta_0\) to ensure that \(\text{mean}(\beta(t)) = \beta_0\), where:

  • \(D_+\) is the number of school days (in a year)

  • \(D_-\) is the number of holidays (in a year)

A closer look at the factor for $\beta_0$

For the un-adjusted model, the formulation for \(\beta(t)\) is as followed

\[ \beta(t) = \beta_0 (1+\beta_1 \text{Term}(t)) \]

The problem of this formulation arises when there is a large discrepancy between number of days where \(\text{Term}(t) = 1\) and \(\text{Term}(t) = -1\), in which case the mean of \(\beta(t)\) will deviate from \(\beta_0\).

The mean for \(\beta(t)\) assuming a period of 1 year can be estimated by

\[ \frac{1}{365} \sum_1^{365} \beta_0 (1 + \beta_1 * \text{Term}(t)) \] \[ \frac{1}{365} \beta_0 \sum_{t=1}^{365} 1 + \beta_1 * \text{Term}(t) \]

Denote \(D_+\) as the number of days where \(\text{Term}(t) = 1\) and \(D_-\) as the number of days where \(\text{Term}(t) = -1\), also recall that \(\sum_{t = 1}^{365} 1 = 365\). We can rewrite the formula into

\[ \frac{1}{365} \beta_0 (365 + \beta_1 D_+ - \beta_1 D_-) \]

For example, if a year have 250 days where \(\text{Term}(t) = 1\), with the settings \(\beta_0 = 0.1\) and \(\beta_1 = 2\), mean of \(\beta(t)\) over the year can be computed by.

\[ \frac{1}{365} \beta_0 (365 + 250\beta_1 - 115\beta_1) \approx 0.17 \]

We can also rewrite that \(365 = D_+ + D_-\)

\[ \frac{1}{365} \beta_0 (D_+ + D_- + \beta_1 D_+ - \beta_1 D_-) = \\ \frac{1}{365} \beta_0 ((1 + \beta_1)D_+ + (1- \beta_1) D_-) \]

Thus, for the mean of \(\beta(t)\) to be equal to \(\beta_0\) we need to multiply a factor \(\frac{1}{\frac{1}{365} ((1 + \beta_1)D_+ + (1 - \beta_1)D_-)}\)

The implementation for SIR model with term-based forcing (using R package odin) is provided below.

sir_term <- odin2::odin({
  N <- S + I + R
  deriv(S) <- b*N - (beta_t*(I/N) + b)*S
  deriv(I) <-beta_t*(I/N)*S - (b + gamma)*I
  deriv(R) <- gamma*I - R*b

  # seasonality forcing
  beta_0 <- parameter(0.4)
  beta_1 <- parameter(0)
  
  school_open <- parameter()
  school_time <- parameter()
  school_data_dim <- parameter()
  dim(school_time) <- school_data_dim
  dim(school_open) <- school_data_dim
  
  term <- interpolate(school_time, school_open, "constant")
  beta_t <- beta_0*(1 + beta_1*term)
  
  # initialize starting population
  init_S <- parameter(9500)
  init_I <- parameter(500)
  gamma <- parameter(0.05)
  b <- parameter(0.05)
  
  initial(S) <- init_S
  initial(I) <- init_I
  initial(R) <- 0

})

get_adjust_factor <- function(b1, formatted_schooldays){
  # compute number of school days for each term
  total_school_days <- sapply(formatted_schooldays$school_days, \(pair){
    pair$max - pair$min + 1
  })
  # compute number of holidays for each term
  total_holidays <- sapply(formatted_schooldays$holidays, \(pair){
    pair$max - pair$min + 1
  })
  
  adjust_factor <- (1/365)*((1+b1)*sum(total_school_days) + (1-b1)*sum(total_holidays))
  1/adjust_factor
}

format_schooldays <- function(school_time, school_open){
  prev <- 0
  school_days <- list()
  holidays <- list()
  sapply(1:length(school_time), \(i){
    if(school_open[i] == -1 && school_time[i]>0){
      school_days[[length(school_days) + 1]] <<- list(min = prev, max = school_time[i]-1)
      prev <<- school_time[i]
    }else if (school_open[i] == 1 && school_time[i]>0){
      holidays[[length(holidays) + 1]] <<- list(min = prev, max = school_time[i]-1)
      prev <<- school_time[i]
    }
  })
  
  # handle cases when the end of school time is not the end of the year
  if(prev < 365){
    if(school_open[length(school_open)] == -1){
      holidays[[length(holidays) + 1]] <- list(min = prev, max = 365)
    }else{
      school_days[[length(school_days) + 1]] <- list(min = prev, max = 365)
    }
  }
  
  list(
    school_days = school_days,
    holidays = holidays
  )
}

# get a list of item, each item is a list of 2, for min x and max x for the annotation
get_annotate_layers <- function(annotate_range, color = "red", alpha = 0.5){
  lapply(annotate_range, \(pair){
    annotate("rect", 
             xmin = pair$min, xmax = pair$max,
             ymin = -Inf, ymax = Inf,
             fill = color, alpha = alpha)
  })
}

Example: we can try running the model using the following set of parameters

• \(\beta_1 = 0.29\) and \(\beta_0 = 0.125\)

• \(\gamma = 1/13\) (mean recovery time of 13 days)

• Initial population: \(S(0) = 499900\) and \(I(0) = 100\)

• The model is simulated for 365 time units, where each unit is 1 day

For \(\text{Term}(t)\), school holidays data for England and Wales is used.

Holiday	Model time
Christmas	356-6
Easter	100-115
Summer	200-251
Autumn Half-term	300-307

The plot below shows how the number of infectives changes over one year, using term-based forcing. Holiday periods are highlighted in red and note that \(\beta_0\) is not adjusted.

school_time <- c(0, 7, 100, 116, 200, 252, 300,308,356)
school_open <- c(-1, 1, -1, 1,   -1,   1,   -1, 1,-1)

term_pars <- list(
  school_time = school_time, 
  school_open = school_open,
  school_data_dim = length(school_time),
  beta_0 = 0.125,
  beta_1 = 0.29,
  gamma = 1/13,
  b = 0.02,
  init_S = 499900,
  init_I = 100
  )

# format school data for generating annotation layer
school_data <- format_schooldays(school_time, school_open)

# run model and plot output
run_mod(sir_term, term_pars, duration = 365, timestep = 1) %>% 
  ggplot() +
    geom_line(aes(x = t, y = I), color = "cornflowerblue") +
  get_annotate_layers(school_data$holiday, alpha = 0.3) +
  labs(
    title = "Number of infectives over a year",
    x = "Time (days)",
    y = "Number of infectives"
  )

Figure 4: Changes in I when beta_0 is not adjusted

Visualization for the same model but with \(\beta_0\) now adjusted by the factor \(\frac{1}{\frac{1}{365} ((1 + \beta_1)D_+ + (1 - \beta_1)D_-)}\) is also provided below.

adjust_factor <- get_adjust_factor(term_pars$beta_1, school_data)

# use adjusted value for beta_0 instead
adjusted_pars <- term_pars
adjusted_pars$beta_0 <- adjust_factor*adjusted_pars$beta_0

# run model again and plot
run_mod(sir_term, adjusted_pars, duration = 365, timestep = 0.5) %>% 
  ggplot() +
    geom_line(aes(x = t, y = I), color = "cornflowerblue") +
  get_annotate_layers(school_data$holiday, alpha = 0.3) +
  labs(
    title = "Number of infectives over a year",
    x = "Time (days)",
    y = "Number of infectives"
  )

Figure 5: Changes in I when beta_0 is adjusted

3.1 Bifurcation diagram

//TBD: ## References {.unnumbered}

References

Anderson, Roy M., and Robert M. May. 1991. Infectious Diseases of Humans: Dynamics and Control. OUP Oxford.

Bailey, N. T. J. 1975. The Mathematical Theory of Infectious Diseases and Its Applications. 2nd Edition. Charles Griffin & Company Ltd, 5a Crendon Street, High Wycombe, Bucks HP13 6LE.

Keeling, Matt J., and Pejman Rohani. 2008. Modeling Infectious Diseases in Humans and Animals. Princeton University Press. https://doi.org/10.2307/j.ctvcm4gk0.