# Computer Simulation of Epidemics with the SIR Model

In this notebook we will explore a simple numerical simulation of how an epidemic spreads through a population. The model is called "SIR", because every individual in the population is assigned to one of three *compartments*:
\begin{align}
S &= \textrm{Susceptible} \\
I &= \textrm{Infected} \\
R &= \textrm{Removed (or Recovered)}
\end{align}



The code snippets in this notebook are written in the [Julia programming language](https://julialang.org/). You can run each cell by placing your cursor inside the cell and either clicking the "Run" button at the top of the page, or typing `Shift+Enter`. Make sure you run the code cells in order, as some of the later cells depend on the earlier ones.


## Numerical Models

Before we begin modelling the epidemic, let's look at a simpler numerical model for the motion of an object that is accelerating. The system has a number of properties of interest:

- $t$, the time (in seconds) since the object started moving
- $x$, the distance (in metres) of the object from its starting point
- $v$, the object's velocity
- $a$, the object's acceleration

These properties are related by the following system of equations:

\begin{align}
v &= \frac{dx}{dt} & \text{velocity is equal to the rate of change of distance with respect to time}\\
a &= \frac{dv}{dt} & \text{acceleration is equal to the rate of change of velocity with respect to time}\\
\text{or equivalently, }
a &= \frac{d^2x}{dt^2}
\end{align}

Given the four properties $[t,x,v,a]$, we can work out what the system will look like at some later time $(t+\Delta t)$ by updating each property according to these equations. We read "$\Delta t$" as "delta t" and it means "a small change in time".

The following function updates distance and velocity according to the above equations, given an acceleration and the timestep $\Delta t$. Note that the update uses a midpoint formula based on half of the "old" velocity (the velocity at the beginning of the timestep) and half of the "new" velocity (the velocity at the end of the timestep).


```julia
function update(x, velocity, acceleration, delta_t)
    new_v = velocity + acceleration * delta_t
    new_x = x + 0.5(velocity * delta_t) + 0.5(new_v * delta_t)
    return new_x, new_v
end
```




    update (generic function with 1 method)





We can use this update function to simulate the motion in the $x$ direction of an object that accelerates from a standing start at 1 metre per second per second*, for example, a shopping trolley rolling downhill along a street where the distance along the street is $x$.

\* That wasn't a stutter. An acceleration of one metre per second per second ($1 \text{ms}^{-2}$) means that every second the trolley increases its speed by one metre per second.


```julia
delta_t = 1.0 # time step is one second
nsteps = 10 # total time ten seconds

### TRY CHANGING THESE VALUES! ###
x = 0.0 # starting position (a distance of 0)
v = 0.0 # starting velocity (0 means object is stationary at time t=0)
a = 1.0 # accelerating at 1m/s^2

distance = Array{Float64}(undef, nsteps) # how far has the object travelled at each timestep?
for t = 1:nsteps
    x, v = update(x, v, a, delta_t)
    distance[t] = x
end

using Plots
gr()
plot(distance[1:nsteps], xlabel="time", ylabel="distance", legend=false)
```




    

    



## Computing the SIR Model

Back to the SIR model: the proportion of individuals in each compartment of the population changes as the epidemic spreads, according to a system of differential equations:

\begin{align}
\frac{dS}{dt} &= -\beta S(t) I(t) & \text{rate of change in susceptible } S \\
\frac{dI}{dt} &= \beta S(t) I(t) - \gamma I(t)  & \text{rate of change in infected } I \\
\frac{dR}{dt} &= \gamma I(t) & \text{rate of change in removed } R
\end{align}

The model parameters $\beta$ (beta) and $\gamma$ (gamma) represent the infection rate and the recovery rate, respectively. Note that $\frac{dS}{dt} + \frac{dI}{dt} + \frac{dR}{dt} = 0$.



There's one other number that is not a parameter to this model, but is often mentioned in discussing epidemics: the *basic reproduction number* $R_0$, which represents the expected number of people who would be directly infected by a single infected person.

We can simulate the change in each compartment over time by repeatedly updating each compartment based on the corresponding differential equation. In the SIR model, this can be calculated as:

\begin{align}
R_0 &= \frac{\beta N}{\gamma},
\end{align}

where $N = S+I+R$ is the total size of the population.

Here is update function to calculate the new proportion of individuals in each of the compartments $[S,I,R]$ after a timestep $\Delta t$:


```julia
function update_SIR(s, i, r, beta, gam, delta_t)
    new_s = s - (beta * s * i * delta_t)
    new_i = i + (beta * s * i * delta_t) - (gam * i * delta_t)
    new_r = r + (gam * i * delta_t)
    return [new_s, new_i, new_r]
end
```




    update_SIR (generic function with 1 method)



We can use this function to calculate the change in each group over time for an imaginary outbreak in a group of 100 people over a period of six months. Initially only one person is infected. Try changing the initial numbers in each compartment, as well as the model parameters beta and gamma, and see how it changes the spread over time.


```julia
num_days = 180
delta_t = 1.0 # timestep of one day
n_steps = round(Int64, num_days / delta_t)
result = Array{Float64}(undef, n_steps, 3)

### TRY CHANGING THESE VALUES! ###
s = 99 # susceptible
i =  1 # infectious
r =  0 # removed
beta = 0.0021 # probability of a infection
gam = 0.1 # 1/recovery time e.g. average 10 days to recover

time = 1:n_steps
for step = time
    result[step,:] = [s, i, r]
    s, i, r = update_SIR(s, i, r, beta, gam, delta_t)
end

N = s+i+r
R0 = beta * N / gam

using Printf
plot(time*delta_t, result[1:n_steps,:], title="Example SIR results", xlabel="Day", ylabel="Number of cases", labels = ["Susceptible" "Infected" "Removed"])
annotate!(150, 50, @sprintf("R0 = %.1f", R0), :blue)
```




    

    



## Validating Predictions Against Real-World Data

We can also look at real world data and see how well the actual progression of COVID-19 is simulated by our SIR model. The file [au_covid.csv](au_covid.csv) contains data for the first 230 days of the COVID-19 outbreak in Australia, taken from the first 230 days of the [2019 Novel Coronavirus COVID-19 (2019-nCoV) Data Repository by Johns Hopkins CSSE](https://github.com/CSSEGISandData/COVID-19/).


```julia
using CSV
x = CSV.read("au_covid.csv", transpose=true)
states = convert(Matrix, x)
state_names = reshape(String.(names(x)), 1, :)
plot(states, title="COVID-19 Cases per State or Territory", xlabel="Day", ylabel="Number of cases", labels=state_names, legend=:topleft)
```




    

    



We'll just look at the national total, as the numbers for most states are too small to accurately model.


```julia
national=sum(states,dims=2)
plot(national,title="COVID-19 Cases in Australia", xlabel="Day", ylabel="Number of cases", legend=false)
```




    

    



For a population of the size of Australia, we can choose arbitrary values of beta and gamma to see how closely the SIR model matches the observed data. Try changing the values in the code below to see how it affects the fit.


```julia
popsize=24600000 # population of Australia
num_days = size(states,1)
delta_t = 1.0 # timestep of one day
n_steps = round(Int64, num_days / delta_t)
result = Array{Float64}(undef, n_steps, 3)

### TRY CHANGING THESE VALUES! ###
infected = 4
beta = 3.5e-9
gam = 0.05 # 20-day infectious period

s = popsize-infected
i = infected
r = 0
time = 1:n_steps
for step = time
    result[step,:] = [s,i,r]
    s, i, r = update_SIR(s, i, r, beta, gam, delta_t)
end
predicted=result[time,2]+result[time,3]

R0 = beta * popsize / gam

plot(time*delta_t, [national[time],predicted], title="Actual vs. Observed Cases", xlabel="Day", ylabel="Number of Cases", labels = ["Actual" "Predicted"])
annotate!(50, 10000, @sprintf("R0 = %.1f", R0), :blue)
```




    

    



Rather than guessing the values of $\beta$ and $\gamma$, it would be better to fit them to the data to minimize the mean squared error between the predicted and actual values for each day. The [DifferentialEquations.jl](https://docs.sciml.ai/v1.6/analysis/parameter_estimation.html) package allows us to perform this fitting in Julia - but that's a subject for a future notebook.

14839b38d182f5a5b694c0d795cfb62356cec67c
ipynb
Jupyter Notebook
sir_model.ipynb
milthorpe/sir
60c5d1c4d9313dca9c2ff455e1f361371db591fd
CC-BY-4.0
sir_model.ipynb
milthorpe/sir
60c5d1c4d9313dca9c2ff455e1f361371db591fd
CC-BY-4.0
sir_model.ipynb
milthorpe/sir
60c5d1c4d9313dca9c2ff455e1f361371db591fd
CC-BY-4.0
Qwen/Qwen-72B
1. YES
2. YES