6  SIR

Recall the deterministic SIR is this:

\[\begin{aligned} \frac{dS(t)}{dt} &= -\frac{\beta S(t)I(t)}{N} \\ \frac{dI(t)}{dt} &= \frac{\beta S(t)I(t)}{N} - \alpha I(t) \\ \frac{dR(t)}{dt} &= \alpha I(t), \end{aligned}\]

A discrete approximation of this model will be:

\[\begin{aligned} S(t+h) &= S(t) - \frac{\beta S(t)I(t)}{N}h \\ I(t+h) &= I(t) + \frac{\beta S(t)I(t)}{N}h - \alpha I(t)h \\ R(t+h) &= R(t) + \alpha I(t)h, \end{aligned}\]

where \(h > 0\) represents the discrete time step.

\[\begin{aligned} \text{P}(\Delta I(t) = 1|I(t)) &= \frac{\beta S(t)I(t)}{N}h + o(h) \\ \text{P}(\Delta R(t) = 1|R(t)) &= \alpha I(t)h + o(h), \end{aligned}\]

new infections and new recoveries occur at the points of two nonhomogeneous Poisson processes with rates \(\frac{\beta S(t)I(t)}{N}h\) and \(\alpha I(t)\), respectively

6.1 Kolmogorov Forward Equations

6.2 The Gillespie Stochastic Simulation Algorithm (Gillespie SSA)

The algorithm proceeds in two steps:

  1. Simulate time at which the next event will occur; in this case an event refers to the movement of one individual either from \(S\) to \(I\) or from \(I\) to \(R\). The time \(Z\) until the next event occurs follows an exponential distribution with a rate equal to the sum of the rates over all possible events (in this case two events are possible, i.e., movement from \(S\) to \(I\) and \(I\) to \(R\)). Denoting \(c_1=\beta S(t)I(t)/N\) and \(c_2=\alpha I(t)\) (note there are as many \(c\)s as there are reactions), \(Z\) is distributed as:

\[g_Z(z)=\left(\sum_{i=1}^2 c_i\right)\exp\left(-z\sum_{i=1}^2 c_i\right).\]

  1. When the event time has been simulated, next simulate which event occurs. Event rates \(c_i\) are converted into probabilities; one of the events is then selected at random according to:

\[\text{P}(\text{Event}=v)=\frac{c_v}{\sum_{i=1}^2 c_i}.\]

6.2.1 Interactive illustration

The widget below makes both Gillespie steps visible. Move the sliders to set the population state \((S, I, R)\) and the rates \(\beta, \gamma\) — the two distributions on the left react instantly. Hit Step to draw one event: panel ① shows where the random waiting time \(\Delta t\) landed under the exponential density \(\text{Exp}(\lambda)\), panel ② shows where a uniform draw \(u\) landed on the \([0, 1]\) probability bar (left of the divider ⇒ infection, right ⇒ recovery), and the trajectory below records the step.

How to read the panels. Panel ① is the density of the inter-event waiting time \(Z\): a tall thin curve when \(\lambda = c_1 + c_2\) is large (events fire quickly), a wide flat curve when \(\lambda\) is small (events are rare). The dot drops onto the curve at the sampled \(\Delta t\), and the new \(t\) on the trajectory advances by exactly that amount. Panel ② splits the \([0,1]\) interval into two coloured segments whose widths are the event probabilities \(c_1/\lambda\) and \(c_2/\lambda\); the marker is the uniform draw, and which side of the divider it lands on selects the event. Drag \(\beta\) up and the red segment grows; drag \(\gamma\) up and the green segment grows.