Epidemiology: The SIR model

From JSXGraph Wiki
Revision as of 12:27, 22 January 2009 by A WASSERMANN (talk | contribs)
The printable version is no longer supported and may have rendering errors. Please update your browser bookmarks and please use the default browser print function instead.

Simulation of differential equations with turtle graphics using JSXGraph.

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]\displaystyle{ S(t) }[/math] be the fraction that is susceptible to an infectious, but not deadly, disease at time t; let [math]\displaystyle{ I(t) }[/math] be the fraction that is infected at time [math]\displaystyle{ t }[/math]; and let [math]\displaystyle{ R(t) }[/math] be the fraction that has recovered. Let [math]\displaystyle{ \beta }[/math] be the rate at which an infected person infects a susceptible person. Let [math]\displaystyle{ \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]\displaystyle{ \frac{dS}{dt} = - \beta I S }[/math]
[math]\displaystyle{ \frac{dR}{dt} = \gamma I }[/math]
[math]\displaystyle{ \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]\displaystyle{ \gamma = 1/3 }[/math]
  • infection rate: one new person every other day, so [math]\displaystyle{ \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

<link rel="stylesheet" type="text/css" href="http://jsxgraph.uni-bayreuth.de/distrib/jsxgraph.css" />
<script type="text/javascript" src="http://jsxgraph.uni-bayreuth.de/distrib/prototype.js"></script>
<script type="text/javascript" src="http://jsxgraph.uni-bayreuth.de/distrib/jsxgraphcore.js"></script>
<form><input type="button" value="clear and run a simulation of 100 days" onClick="clearturtle();run()"></form>
<div id="box" class="jxgbox" style="width:600px; height:450px;"></div>
var brd = JXG.JSXGraph.initBoard('box', {originX: 20, originY: 300, unitX: 20, unitY: 250});

var S = brd.createElement('turtle',[],{strokeColor:'blue',strokeWidth:3});
var I = brd.createElement('turtle',[],{strokeColor:'red',strokeWidth:3});
var R = brd.createElement('turtle',[],{strokeColor:'green',strokeWidth:3});
            
var xaxis = brd.createElement('axis', [[0,0], [1,0]], {});
var yaxis = brd.createElement('axis', [[0,0], [0,1]], {});
            
var s = brd.createElement('slider', [[0,-0.3], [30,-0.3],[0,0.03,1]], {name:'s'});
brd.createElement('text', [40,-0.3, "initially infected population rate"]);
var beta = brd.createElement('slider', [[0,-0.4], [30,-0.4],[0,0.5,1]], {name:'&beta;'});
brd.createElement('text', [40,-0.4, "&beta;: infection rate"]);
var gamma = brd.createElement('slider', [[0,-0.5], [30,-0.5],[0,0.3,1]], {name:'&gamma;'});
brd.createElement('text', [40,-0.5, "&gamma;: recovery rate = 1/(days of infection)"]);

var t = 0; // global

brd.createElement('text', [40,-0.2, 
        function() {return "Day "+t+": infected="+brd.round(7900000*I.pos[1],1)+" recovered="+brd.round(7900000*R.pos[1],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.X());
                
  delta = 1; // global
  t = 0;  // global
  loop();
}
             
function turtleMove(turtle,dx,dy) {
  turtle.lookTo([1.0+turtle.pos[0],dy+turtle.pos[1]]);
  turtle.fd(dx*Math.sqrt(1+dy*dy));
}
             
function loop() {
  var dS = -beta.Value()*S.pos[1]*I.pos[1];
  var dR = gamma.Value()*I.pos[1];
  var dI = -(dS+dR);
  turtleMove(S,delta,dS);
  turtleMove(R,delta,dR);
  turtleMove(I,delta,dI);
                
  t += delta;
  if (t<100.0) {
    setTimeout(loop,10);
  }
}

References