Simulating stochastic epidemics

This R Shiny app (hosted at https://edbonneville.shinyapps.io/gillespie_shiny/) is an addendum to my master's thesis, which was on the statistical forecasting of infectious disease epidemics using the time between successive infections (infection interevent times).

The app allows the user to simulate and explore susceptible-infectious (SI), susceptible-infectious-susceptible (SIS) and susceptible-infectious-recovered (SIR) type epidemics stochastically with the Gillespie algorithm. All source code is available on this Github repository.

Running it locally

To run the app locally on your personal machine, you need to first install the packages shiny, shinydashboard, plotly, shinyWidgets, tidyverse, and deSolve.

You can then run the following chunk of code:

library(shiny)
shiny::runGitHub("gillespie_shiny", "edbonneville")

Features and usage

The panel on the left allows you to select the:

• Epidemic type: SI, SIS or SIR (without demographics).
• Infection parameter β scaled per year.
• Average infectious period 1/γ in days, in the case of SIS or SIR.
• Shape parameter α of the Gamma distribution of interevent times. By default this is exponential (α = 1), but it can be modified.
• Initial population conditions: number of susceptible S, infectious I and total population size N. Note that N = S + I, and for computational reasons we restrict N ≤ 1000.
• Time frame t_end, i.e. how long we should simulate for. For SI epidemics, the susceptible pool may be depleted before t_end. Analogously for SIR epidemics, the entire population may have recovered before t_end. We restrict t_end to a maximum of 365 days.

On the right panel, in addition to visualising the evolution of the number of susceptible, infectious and recovered individuals over time, you can also:

• Overlay the deterministic (ODE solutions) curves, or view them on their own.
• Visualise the rates of the different events (infections or recoveries) over time.
• View the distribution of infection interevent times.
• Plot infection interevent times against the number of infectious.

The last two features are directly relevant to the thesis. Finally, you can also download a copy of the simulated dataset as a .csv.

Note: in the SIS and SIR case, it is possible after pressing 'Run' that the epidemic will end prematurely. This may be due to a) your chosen parameterisation (e.g. β < γ), or b) due to the initial conditions. If you start with I = 1, the first drawn event (by chance) may be a recovery, leading to extinction. For the latter case, you can either run the simulation again, of increase your initial I to ensure the epidemic takes off.

Downloading data and running a Gamma GLM with inverse link

If you would like to work further along the lines of the thesis, you can model the infection interevent times using a Gamma Generalised Linear Model (GLM) with an inverse link function.

With current app, let us assume you simulate an epidemic under an SIS model and download the data as "SIS_sim". You can then run the following code:

library(tidyverse)

dat <- read.csv("SIS_sim.csv") 

dat_glm <- dat %>% 
  
  # Subset only infections
  filter(transition == "infection") %>% 
  
  # Compute interevent times
  mutate(interev = t - c(0, t[-length(t)])) %>% 
  
  # Remove first row (interev = 0, since first event)
  slice(-1) %>% 
  
  # Add explanatory variables
  mutate(x1 = I, x2 = I^2)

# Run gamma glm with inverse link
glm(interev ~ x1 + x2, family = Gamma(link = "inverse"), data = dat_glm)

Acknowledgements and extra notes

The layout of the app was in part inspired by the "Explore ID Model" Shiny app developed by Sam Abbott at the University of Bristol.

Feel free to use the app and source code for educational and/or personal purposes, and if you have any questions or issues to report feel free to contact me on [email protected].

Edouard Bonneville