Various implementations of the classical SIR model in Julia
Try the notebooks out in Binder:
The model equations are as follows.
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.
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
- 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
- Jump process (Gillespie) using DifferentialEquations.jl
- Jump process (Gillespie) using reaction networks from Catalyst.jl
- Jump process (Gillespie) using Gillespie.jl
- Jump process using the Sellke construction
- Jump process using ModelingToolkit.jl
- Jump process using Fleck.jl
- Petri net model to ODEs, SDEs, and jump process using Petri.jl
- Petri net model to ODEs, SDEs, and jump process using AlgebraicPetri.jl
- Agent-based model using base Julia as well as using DifferentialEquations
- Agent-based model using Agents.jl
- Transmission network individual-based model using Pathogen.jl
- Optimal control of an SIR epidemic with non-pharmaceutical interventions using ODEs and Optimization.jl
- Optimal control of an SIR epidemic with non-pharmaceutical interventions using ODEs and InfiniteOpt.jl
- Optimal control of an SIR epidemic with non-pharmaceutical interventions using a function map and JuMP.jl
- Flattening the curve of an SIR epidemic with non-pharmaceutical interventions using a function map and JuMP.jl
- Optimal control of an SIR epidemic with vaccination using a function map and JuMP.jl
- Optimal control of an SIR epidemic with non-pharmaceutical interventions using a function map and stochastic dual dynamic programming with SDDP.jl
Building models from smaller, re-usable components make it easier to build complex models quickly, and also makes it easier to document the development of these models.
- Composition of ODE models using ModelingToolkit.jl
- Composition of ODE models using AlgebraicDynamics.jl
- Composition of petri net models using AlgebraicPetri.jl
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).
- Changing parameter values at fixed times e.g. lockdown in an SIR model
- Stopping simulations when infected individuals reach zero in stochastic differential equations
- Scheduling recovery times to model a fixed infectious period
- Preventing out-of-domain errors in a sinusoidally forced ODE model
- A hybrid ODE/jump process model
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 time-varying parameters of an ODE system using Turing.jl
- Bayesian estimates of parameters of the ODE system using NestedSamplers.jl
- Bayesian estimates of parameters of the ODE system using importance sampling, Markov Chain Monte Carlo and Sequential Monte Carlo with Gen.jl
- Bayesian estimates of parameters of the ODE system using message passing using RxInfer.jl
- Bayesian inference of transmission network individual-based model parameters using Pathogen.jl
- Estimating the likelihood of a discrete-time Markov model using a simple particle filter
- Point estimates of parameters of a discrete-time Markov model using Ensemble Kalman Inversion
- Approximate posterior estimates of parameters of a discrete-time Markov model using Ensemble Kalman Sampling
- Bayesian estimates of parameters of a discrete-time Markov model using importance sampling, Markov Chain Monte Carlo and Sequential Monte Carlo with Gen.jl
- Steady state analysis of an SIR model with births and deaths
- Bifurcation/stroboscopic plot of a sinusoidally forced ODE model using brute force simulation
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
ProbNumDiffEq.jl
- 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
- Full likelihood intervals using ProfileLikelihood.jl
- Profile likelihood using ProfileLikelihood.jl
- Profile likelihood using LikelihoodProfiler.jl
- 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 basis
- Fixed (rather than exponential) distribution of infectious period:
- Using a delay differential equation and DelayDiffEq.jl
- Using a fixed delay in a jump system, making use of the integrator
- Using a gamma-distributed delay in a jump system using DelaySSAToolkit.jl
- Using a stochastic delay differential equation in StochasticDelayDiffEq.jl
- Using a fixed delay in a discrete event simulation using SimJulia (see end of file)
- An Erlang distribution for the infectious period using the method of stages:
- Multigroup models
- Fractional differential equations
- A hybrid ODE/jump process model
While this repository is mainly about Julia, it is also possible to use Julia to call code written in other languages. Here are some examples of how to define the vector field of an ODE in C, Python, and R.
Many languages can compile to shared libraries that can be accessed via ccall
. Here are examples of how to define the vector field of an ODE in various languages, and call it using ccall
.
- ODE with derivatives in C
- ODE with derivatives in Fortran 90
- ODE with derivatives in Rust
- ODE with derivatives in FreePascal
- ODE with derivatives in Zig
- ODE with derivatives in Nim
- ODE with derivatives in Python, accessed via PythonCall.jl, also demonstrating Python to Julia code conversion using
modelingtoolkitize
- ODE with derivatives in R, accessed via RCall.jl
- ODE with derivatives defined using the
odin
R package, accessed via RCall.jl
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 StaticArrays
.
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.
Thanks to Weave.jl
, Julia Markdown files (in tutorials/
) are converted into multiple formats.
git clone https://github.com/epirecipes/sir-julia
cd sir-julia
Then launch julia
and run the following.
cd(@__DIR__)
import IJulia
IJulia.notebook(;dir="notebook")
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
Suggested sections:
- Introduction
- Libraries
- Utility functions
- Transitions
- Time domain
- Initial conditions
- Parameter values
- Random number seed
- Running the model
- Post-processing
- Plotting
- Benchmarking
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 Project.toml
.
Examples use the following libraries (see the Project.toml
file for a full list of dependencies):
- The
DifferentialEquations.jl
ecosystem for many of the examples SimJulia
for discrete event simulationsAgents.jl
for agent-based modelsGillespie.jl
for the Doob-Gillespie processPetri.jl
for the Petri net modelsAlgebraicPetri.jl
for a category theory based modeling framework for creating Petri net modelsTuring.jl
for inference using probabilistic programsNestedSamplers.jl
for nested samplingGpABC
for inference using Approximate Bayesian ComputationSoss.jl
for Markov modelsMomentClosure.jl
for moment closureBridge.jl
for stochastic differential equations
Parts of the code were taken from @ChrisRackauckas DiffEqTutorials
, which comes highly recommended.