Various implementations of the classical SIR model in Julia
Try the notebooks out in Binder:
GitHub Markdown doesn't parse equations, so here's a description of the underlying SIR model.
- The ordinary differential equation model considers:
- Susceptible, S, with initial condition S(0)=990
- Infected, I, with initial condition, I(0)=10
- Recovered, R, with initial condition R(0)=10
- Total population, N=S+I+R=1000
- Susceptible individuals make contacts with others at rate c (=10.0), with the probability of a contact with an infectious person being I/N. With probability β (=0.05), an infected person will infect a susceptible given a contact.
- Infected individuals recover at a per-capita rate γ (=0.25).
There are two types of parameterization commonly used in this project; the 'standard' version, that considers the number of individuals in the S, I, and R groups, and an alternative version, in which the dynamics of transmission (βSI/N) and recovery (γI) are modelled directly, with S, I, and R being calculated based on these dynamics and the initial conditions for S, I and R.
Simulation with different types of model
The above process can be represented in different kinds of ways:
- Ordinary differential equation using the Euler method
- Ordinary differential equation using DifferentialEquations.jl
- Ordinary differential equation using ModelingToolkit.jl
- Ordinary differential equation using Modia.jl
- Ordinary differential equation using ApproxFun.jl
- Ordinary differential equation with composition using AlgebraicDynamics.jl
- Volterra integral equation using the Adomian decomposition method
- Stochastic differential equation using DifferentialEquations.jl
- Stochastic differential equation using StochasticDiffEq.jl
- Stochastic differential equation using Bridge.jl
- Linear noise approximation (LNA) to the stochastic differential equation
- Multivariate birth process reparameterisation of the stochastic differential equation
- ODEs of means, variances, etc. through moment closure
- Function map
- Function map using DynamicalSystems.jl
- Function map using ModelingToolkit.jl
- Stochastic Markov model
- Stochastic Markov model using Soss.jl
- Jump process (Gillespie) using DifferentialEquations.jl
- Jump process (Gillespie) using reaction networks from Catalyst.jl
- Reaction network conversion to ODEs, SDEs and jump process using ModelingToolkit
- Petri net model to ODEs, SDEs, and jump process using Petri.jl
- Petri net model to ODEs, SDEs, and jump process using AlgebraicPetri.jl
- Jump process (Gillespie) using Gillespie.jl
- Jump process using the Sellke construction
- Jump process using ModelingToolkit.jl
- Discrete event simulation using SimJulia
- Agent-based model using base Julia as well as using DifferentialEquations
- Agent-based model using Agents.jl
Generating simulated data
We usually do not observe the trajectory of susceptible, infected, and recovered individuals. Rather, we often obtain data in terms of new cases aggregated over a particular timescale (e.g. a day or a week).
Use of callbacks
- Changing parameter values at fixed times e.g. lockdown in an SIR model
- Preventing negative populations in stochastic differential equations
- Scheduling recovery times to model a fixed infectious period
In addition to the above examples of simulation, there are also examples of inference of the parameters of the model using counts of new cases. Although these are toy examples, they provide the building blocks for more complex situations.
- Point estimates of parameters of the ODE system using Optim.jl and DiffEqParamEstim.jl
- Bayesian estimates of parameters of the ODE system using Approximate Bayesian Computation
- Bayesian estimates of parameters of the ODE system using Turing.jl
- Bayesian estimates of parameters of the ODE system using NestedSamplers.jl
In conducting inference, it is important to know the extent to which parameters are identifiable from the available data.
Incorporating uncertainty in its many forms is important for using models to make decisions.
When solving continuous-time models like ODEs, the discretization can lead to numerical errors. Probabilistic integration treats this error as a statistical problem to capture the uncertainty in the model outputs generated using the solver.
- Probabilistic integration of an ODE model using Bayesian filtering and
- Probabilistic integration of an ODE model by converting to a SDE using
DiffEqUncertainty.jl(can also handle SDEs and DDEs)
- Global sensitivity analysis of an ODE model using Latin hypercube sampling
- Global sensitivity analysis of an ODE model using multiple algorithms from GlobalSensitivity.jl
- Uncertainty propagation of an ODE model using MonteCarloMeasurements.jl
- Uncertainty propagation of an ODE model using the Koopman expectation and DiffEqUncertainty.jl
- Uncertainty analysis of an ODE model with an uncertain input and an uncertain output using Bayesian melding
- Surrogate models (single input/single output) of an ODE model using Surrogates.jl
- Gaussian process surrogate model (two inputs/one output) of an ODE model using the Python package
mogp-emulator, with an example of history matching
- Bayes linear surrogate model (two inputs/one output) of an ODE model using the R package
hmer, with an example of history matching
- Data-driven models of an ODE using DataDrivenDiffEq.jl
- A neural ODE, which replaces all the derivatives of the model with a neural network
- A universal ODE, modeling force of infection using a neural network
- A partially specified ODE, modeling force of infection using a basus
- Fixed (rather than exponential) distribution of infectious period:
- An Erlang distribution for the infectious period using the method of stages is illustrated in this notebook using AlgebraicDynamics.jl
Comments on implementations
Note that the implementations and choice of parameters may be suboptimal, and are intended to illustrate more-or-less the same underlying biological process with different mathematical representations. Additional optimizations may be obtained e.g. by using
I've also tried to transform parameterisations in discrete time as closely as possible to their continuous counterparts. Please see the great work by Linda Allen for how these different representations compare.
Types of output
Weave.jl, Julia Markdown files (in
tutorials/) are converted into multiple formats.
git clone https://github.com/epirecipes/sir-julia cd sir-julia
julia and run the following.
cd(@__DIR__) import IJulia IJulia.notebook(;dir="notebook")
Adding new examples
Plans for new examples are typically posted on the Issues page.
To add an example, make a new subdirectory in the
tutorials directory, and add a Julia Markdown (
.jmd) document to it. Set the beginning to something like the following:
# Agent-based model using Agents.jl Simon Frost (@sdwfrost), 2020-04-27
- Utility functions
- Time domain
- Initial conditions
- Parameter values
- Random number seed
- Running the model
Change to the root directory of the repository and run the following from within Julia; you will need Weave.jl and any dependencies from the tutorial.
include("build.jl") weave_all() # or e.g. weave_folder("abm") for an individual tutorial
If additional packages are added, then these need to be added to
build_project_toml.jl, which when run, will regenerate
Examples use the following libraries (see the
Project.toml file for a full list of dependencies):
DifferentialEquations.jlecosystem for many of the examples
SimJuliafor discrete event simulations
Agents.jlfor agent-based models
Gillespie.jlfor the Doob-Gillespie process
Petri.jlfor the Petri net models
AlgebraicPetri.jlfor a category theory based modeling framework for creating Petri net models
Turing.jlfor inference using probabilistic programs
NestedSamplers.jlfor nested sampling
GpABCfor inference using Approximate Bayesian Computation
Soss.jlfor Markov models
MomentClosure.jlfor moment closure
Bridge.jlfor stochastic differential equations
Parts of the code were taken from @ChrisRackauckas
DiffEqTutorials, which comes highly recommended.