Epidemiology: The SIR model
Simulation of differential equations with turtle graphics using JSXGraph.
Contents
SIR model without vital dynamics
The SIR model measures the number of susceptible, infected, and recovered individuals in a host population. Given a fixed population, let [math]S(t)[/math] be the fraction that is susceptible to an infectious, but not deadly, disease at time t; let [math]I(t)[/math] be the fraction that is infected at time [math]t[/math]; and let [math]R(t)[/math] be the fraction that has recovered. Let [math]\beta[/math] be the rate at which an infected person infects a susceptible person. Let [math]\gamma[/math] be the rate at which infected people recover from the disease.
A single epidemic outbreak is usually far more rapid than the vital dynamics of a population, thus, if the aim is to study the immediate consequences of a single epidemic, one may neglect birth-death processes. In this case the SIR system can be expressed by the following set of differential equations:
- [math] \frac{dS}{dt} = - \beta I S [/math]
- [math] \frac{dR}{dt} = \gamma I [/math]
- [math] \frac{dI}{dt} = -(\frac{dS}{dt}+\frac{dR}{dt}) [/math]
Example Hong Kong flu
- initially 7.9 million people,
- 10 infected,
- 0 recovered.
- estimated average period of infection: 3 days, so [math]\gamma = 1/3[/math]
- infection rate: one new person every other day, so [math]\beta = 1/2[/math]
Thus S(0) = 1, I(0) = 1.27E-6, R(0) = 0, see [1].
The lines in the JSXGraph-simulation below have the following meaning:
* Blue: Rate of susceptible population * Red: Rate of infected population * Green: Rate of recovered population (which means: immune, isolated or dead)
The underlying JavaScript code
var brd = JXG.JSXGraph.initBoard('box', {axis: true, boundingbox: [-5, 1.2, 100, -1.2]});
var S = brd.create('turtle',[],{strokeColor:'blue',strokeWidth:3});
var I = brd.create('turtle',[],{strokeColor:'red',strokeWidth:3});
var R = brd.create('turtle',[],{strokeColor:'green',strokeWidth:3});
var s = brd.create('slider', [[0,-0.3], [30,-0.3],[0,1.27E-6,1]], {name:'s'});
brd.create('text', [40,-0.3, "initially infected population rate (on load: I(0)=1.27E-6)"]);
var beta = brd.create('slider', [[0,-0.4], [30,-0.4],[0,0.5,1]], {name:'β'});
brd.create('text', [40,-0.4, "β: infection rate"]);
var gamma = brd.create('slider', [[0,-0.5], [30,-0.5],[0,0.3,1]], {name:'γ'});
brd.create('text', [40,-0.5, "γ: recovery rate = 1/(days of infection)"]);
var t = 0; // global
brd.create('text', [40,-0.2,
function() {return "Day "+t+": infected="+(7900000*I.Y()).toFixed(1)+" recovered="+(7900000*R.Y()).toFixed(1);}]);
S.hideTurtle();
I.hideTurtle();
R.hideTurtle();
function clearturtle() {
S.cs();
I.cs();
R.cs();
S.hideTurtle();
I.hideTurtle();
R.hideTurtle();
}
function run() {
S.setPos(0,1.0-s.Value());
R.setPos(0,0);
I.setPos(0,s.Value());
delta = 1; // global
t = 0; // global
loop();
}
function turtleMove(turtle,dx,dy) {
turtle.moveTo([dx+turtle.X(),dy+turtle.Y()]);
}
function loop() {
var dS = -beta.Value()*S.Y()*I.Y();
var dR = gamma.Value()*I.Y();
var dI = -(dS+dR);
turtleMove(S,delta,dS);
turtleMove(R,delta,dR);
turtleMove(I,delta,dI);
t += delta;
if (t<100.0) {
active = setTimeout(loop,10);
}
}
function stop() {
if (active) clearTimeout(active);
active = null;
}
function goOn() {
if (t>0) {
if (active==null) {
active = setTimeout(loop,10);
}
} else {
run();
}
}
