Difference between revisions of "Epidemiology: The SIR model"

From JSXGraph Wiki
Jump to navigationJump to search
 
(24 intermediate revisions by 2 users not shown)
Line 13: Line 13:
  
 
:<math> \frac{dI}{dt} = -(\frac{dS}{dt}+\frac{dR}{dt}) </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 [http://www.cs.princeton.edu/introcs/94diffeq/].
  
 
The lines in the JSXGraph-simulation below have the following meaning:
 
The lines in the JSXGraph-simulation below have the following meaning:
Line 19: Line 28:
 
  * <span style="color:green">Green: Rate of recovered population (which means: immune, isolated or dead)
 
  * <span style="color:green">Green: Rate of recovered population (which means: immune, isolated or dead)
 
<html>
 
<html>
<link rel="stylesheet" type="text/css" href="http://jsxgraph.uni-bayreuth.de/distrib/jsxgraph.css" />
+
<form><input type="button" value="clear and run a simulation of 100 days" onClick="clearturtle();run()">
<script type="text/javascript" src="http://jsxgraph.uni-bayreuth.de/distrib/prototype.js"></script>
+
<input type="button" value="stop" onClick="stop()">
<script type="text/javascript" src="http://jsxgraph.uni-bayreuth.de/distrib/jsxgraphcore.js"></script>
+
<input type="button" value="continue" onClick="goOn()"></form>
<form><input type="button" value="clear and run" onClick="clearturtle();run()"></form>
+
</html>
<div id="box" class="jxgbox" style="width:600px; height:450px;"></div>
+
 
<script language="JavaScript">
+
<jsxgraph box="box" width="600" height="450">
var brd = JXG.JSXGraph.initBoard('box', {originX: 20, originY: 300, unitX: 6, unitY: 250});
+
var brd = JXG.JSXGraph.initBoard('box', {axis: true, boundingbox: [-5, 1.2, 100, -1.2]});
  
var S = brd.createElement('turtle',[],{strokeColor:'blue',strokeWidth:3});
+
var S = brd.create('turtle',[],{strokeColor:'blue',strokeWidth:3});
var I = brd.createElement('turtle',[],{strokeColor:'red',strokeWidth:3});
+
var I = brd.create('turtle',[],{strokeColor:'red',strokeWidth:3});
var R = brd.createElement('turtle',[],{strokeColor:'green',strokeWidth:3});
+
var R = brd.create('turtle',[],{strokeColor:'green',strokeWidth:3});
 
              
 
              
var xaxis = brd.createElement('axis', [[0,0], [1,0]], {});
+
var s = brd.create('slider', [[0,-0.3], [30,-0.3],[0,1.27E-6,1]], {name:'s'});
var yaxis = brd.createElement('axis', [[0,0], [0,1]], {});
+
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:'&beta;'});
var s = brd.createElement('slider', [[0,-0.3], [10,-0.3],[0,0.03,1]], {name:'s'});
+
brd.create('text', [40,-0.4, "&beta;: infection rate"]);
brd.createElement('text', [12,-0.3, "initially infected population rate"]);
+
var gamma = brd.create('slider', [[0,-0.5], [30,-0.5],[0,0.3,1]], {name:'&gamma;'});
var beta = brd.createElement('slider', [[0,-0.4], [10,-0.4],[0,0.5,1]], {name:'&beta;'});
+
brd.create('text', [40,-0.5, "&gamma;: recovery rate = 1/(days of infection)"]);
brd.createElement('text', [12,-0.4, "&beta;: infection rate"]);
+
 
var gamma = brd.createElement('slider', [[0,-0.5], [10,-0.5],[0,0.3,1]], {name:'&gamma;'});
+
var t = 0; // global
brd.createElement('text', [12,-0.5, "&gamma;: recovery rate"]);
 
  
brd.createElement('text', [12,-0.2,  
+
brd.create('text', [40,-0.2,  
         function() {return "S(t)="+brd.round(S.pos[1],3) +", I(t)="+brd.round(I.pos[1],3) +", R(t)="+brd.round(R.pos[1],3);}]);
+
         function() {return "Day "+t+": infected="+(7900000*I.Y()).toFixed(1)+" recovered="+(7900000*R.Y()).toFixed(1);}]);
 
              
 
              
 
S.hideTurtle();
 
S.hideTurtle();
Line 61: Line 69:
 
   S.setPos(0,1.0-s.Value());
 
   S.setPos(0,1.0-s.Value());
 
   R.setPos(0,0);
 
   R.setPos(0,0);
   I.setPos(0,s.X());
+
   I.setPos(0,s.Value());
 
                  
 
                  
 
   delta = 1; // global
 
   delta = 1; // global
   t = 0.0;  // global
+
   t = 0;  // global
 
   loop();
 
   loop();
 
}
 
}
 
              
 
              
 
function turtleMove(turtle,dx,dy) {
 
function turtleMove(turtle,dx,dy) {
   turtle.lookTo([1.0+turtle.pos[0],dy+turtle.pos[1]]);
+
   turtle.moveTo([dx+turtle.X(),dy+turtle.Y()]);
  turtle.fd(dx*Math.sqrt(1+dy*dy));
 
 
}
 
}
 
              
 
              
 
function loop() {
 
function loop() {
   var dS = -beta.Value()*S.pos[1]*I.pos[1];
+
   var dS = -beta.Value()*S.Y()*I.Y();
   var dR = gamma.Value()*I.pos[1];
+
   var dR = gamma.Value()*I.Y();
 
   var dI = -(dS+dR);
 
   var dI = -(dS+dR);
 
   turtleMove(S,delta,dS);
 
   turtleMove(S,delta,dS);
Line 82: Line 89:
 
                  
 
                  
 
   t += delta;
 
   t += delta;
   if (t<30.0) {
+
   if (t<100.0) {
     setTimeout(loop,10);
+
     active = setTimeout(loop,10);
 
   }
 
   }
 
}
 
}
           
 
</script>
 
</html>
 
  
====Example Hong Kong flu====
+
function stop() {
* initially 7.9 million people,  
+
  if (active) clearTimeout(active);
* 10 infected,
+
  active = null;
* 0 recovered.
+
}
* estimated average period of infection: 3 days, so <math>\gamma = 1/3</math>
+
function goOn() {
* infection rate: one new person every day, so <math>\beta = 1/2</math>
+
  if (t>0) {
 +
    if (active==null) {
 +
      active = setTimeout(loop,10);
 +
    }
 +
  } else {
 +
    run();
 +
  }
  
Thus S(0) = 1, I(0) = 1.27E-6, R(0) = 0, see See [http://www.cs.princeton.edu/introcs/94diffeq/].
+
}
 +
</jsxgraph>
  
 
===The underlying JavaScript code===
 
===The underlying JavaScript code===
<source lang="html4strict">
 
<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" onClick="clearturtle();run()"></form>
 
<div id="box" class="jxgbox" style="width:600px; height:450px;"></div>
 
</source>
 
 
 
<source lang="javascript">
 
<source lang="javascript">
var brd = JXG.JSXGraph.initBoard('box', {originX: 20, originY: 300, unitX: 20, unitY: 250});
+
var brd = JXG.JSXGraph.initBoard('box', {axis: true, boundingbox: [-5, 1.2, 100, -1.2]});
  
var S = brd.createElement('turtle',[],{strokeColor:'blue',strokeWidth:3});
+
var S = brd.create('turtle',[],{strokeColor:'blue',strokeWidth:3});
var I = brd.createElement('turtle',[],{strokeColor:'red',strokeWidth:3});
+
var I = brd.create('turtle',[],{strokeColor:'red',strokeWidth:3});
var R = brd.createElement('turtle',[],{strokeColor:'green',strokeWidth:3});
+
var R = brd.create('turtle',[],{strokeColor:'green',strokeWidth:3});
 
              
 
              
var xaxis = brd.createElement('axis', [[0,0], [1,0]], {});
+
var s = brd.create('slider', [[0,-0.3], [30,-0.3],[0,1.27E-6,1]], {name:'s'});
var yaxis = brd.createElement('axis', [[0,0], [0,1]], {});
+
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:'&beta;'});
var s = brd.createElement('slider', [[0,-0.3], [10,-0.3],[0,0.03,1]], {name:'s'});
+
brd.create('text', [40,-0.4, "&beta;: infection rate"]);
brd.createElement('text', [12,-0.3, "initially infected population rate"]);
+
var gamma = brd.create('slider', [[0,-0.5], [30,-0.5],[0,0.3,1]], {name:'&gamma;'});
var beta = brd.createElement('slider', [[0,-0.4], [10,-0.4],[0,0.5,1]], {name:'&beta;'});
+
brd.create('text', [40,-0.5, "&gamma;: recovery rate = 1/(days of infection)"]);
brd.createElement('text', [12,-0.4, "&beta;: infection rate"]);
+
 
var gamma = brd.createElement('slider', [[0,-0.5], [10,-0.5],[0,0.3,1]], {name:'&gamma;'});
+
var t = 0; // global
brd.createElement('text', [12,-0.5, "&gamma;: recovery rate"]);
 
  
brd.createElement('text', [12,-0.2,  
+
brd.create('text', [40,-0.2,  
         function() {return "S(t)="+brd.round(S.pos[1],3) +", I(t)="+brd.round(I.pos[1],3) +", R(t)="+brd.round(R.pos[1],3);}]);
+
         function() {return "Day "+t+": infected="+(7900000*I.Y()).toFixed(1)+" recovered="+(7900000*R.Y()).toFixed(1);}]);
 
              
 
              
 
S.hideTurtle();
 
S.hideTurtle();
Line 145: Line 147:
 
   S.setPos(0,1.0-s.Value());
 
   S.setPos(0,1.0-s.Value());
 
   R.setPos(0,0);
 
   R.setPos(0,0);
   I.setPos(0,s.X());
+
   I.setPos(0,s.Value());
 
                  
 
                  
   delta = 0.3; // global
+
   delta = 1; // global
   t = 0.0;  // global
+
   t = 0;  // global
 
   loop();
 
   loop();
 
}
 
}
 
              
 
              
 
function turtleMove(turtle,dx,dy) {
 
function turtleMove(turtle,dx,dy) {
   turtle.lookTo([1.0+turtle.pos[0],dy+turtle.pos[1]]);
+
   turtle.moveTo([dx+turtle.X(),dy+turtle.Y()]);
  turtle.fd(dx*Math.sqrt(1+dy*dy));
 
 
}
 
}
 
              
 
              
 
function loop() {
 
function loop() {
   var dS = -beta.Value()*S.pos[1]*I.pos[1];
+
   var dS = -beta.Value()*S.Y()*I.Y();
   var dR = gamma.Value()*I.pos[1];
+
   var dR = gamma.Value()*I.Y();
 
   var dI = -(dS+dR);
 
   var dI = -(dS+dR);
 
   turtleMove(S,delta,dS);
 
   turtleMove(S,delta,dS);
Line 166: Line 167:
 
                  
 
                  
 
   t += delta;
 
   t += delta;
   if (t<30.0) {
+
   if (t<100.0) {
     setTimeout(loop,10);
+
     active = setTimeout(loop,10);
 
   }
 
   }
}</source>
+
}
 +
 
 +
function stop() {
 +
  if (active) clearTimeout(active);
 +
  active = null;
 +
}
 +
function goOn() {
 +
  if (t>0) {
 +
    if (active==null) {
 +
      active = setTimeout(loop,10);
 +
    }
 +
  } else {
 +
    run();
 +
  }
 +
 
 +
}
 +
</source>
 +
 
 +
===See also===
 +
* [[Epidemiology: The SEIR model]]
 +
* [[Population growth models]]
 +
* [[Lotka-Volterra equations]]
  
 
===References===
 
===References===

Latest revision as of 15:59, 20 February 2013

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]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:'&beta;'});
brd.create('text', [40,-0.4, "&beta;: infection rate"]);
var gamma = brd.create('slider', [[0,-0.5], [30,-0.5],[0,0.3,1]], {name:'&gamma;'});
brd.create('text', [40,-0.5, "&gamma;: 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();
   }

}

See also

References