Epidemiology: The SIR model: Difference between revisions
A WASSERMANN (talk | contribs) No edit summary  | 
				A WASSERMANN (talk | contribs) No edit summary  | 
				||
| (25 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>  | ||
<  | <form><input type="button" value="clear and run a simulation of 100 days" onClick="clearturtle();run()">  | ||
<  | <input type="button" value="stop" onClick="stop()">  | ||
<input type="button" value="continue" onClick="goOn()"></form>  | |||
</html>  | |||
<  | |||
<jsxgraph box="box" width="600" height="450">  | |||
var brd = JXG.JSXGraph.initBoard('box', {  | var brd = JXG.JSXGraph.initBoard('box', {axis: true, boundingbox: [-5, 1.2, 100, -1.2]});  | ||
var S = brd.  | var S = brd.create('turtle',[],{strokeColor:'blue',strokeWidth:3});  | ||
var I = brd.  | var I = brd.create('turtle',[],{strokeColor:'red',strokeWidth:3});  | ||
var R = brd.  | 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:'β'});  | |||
var s = brd.  | brd.create('text', [40,-0.4, "β: infection rate"]);  | ||
brd.  | var gamma = brd.create('slider', [[0,-0.5], [30,-0.5],[0,0.3,1]], {name:'γ'});  | ||
var beta = brd.  | brd.create('text', [40,-0.5, "γ: recovery rate = 1/(days of infection)"]);  | ||
brd.  | |||
var gamma = brd.  | var t = 0; // global  | ||
brd.  | |||
brd.  | brd.create('text', [40,-0.2,    | ||
         function() {return "  |          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.  |    I.setPos(0,s.Value());  | ||
   delta =   |    delta = 1; // global  | ||
   t =   |    t = 0;  // global  | ||
   loop();  |    loop();  | ||
}  | }  | ||
function turtleMove(turtle,dx,dy) {  | function turtleMove(turtle,dx,dy) {  | ||
   turtle.  |    turtle.moveTo([dx+turtle.X(),dy+turtle.Y()]);  | ||
}  | }  | ||
function loop() {  | function loop() {  | ||
   var dS = -beta.Value()*S.  |    var dS = -beta.Value()*S.Y()*I.Y();  | ||
   var dR = gamma.Value()*I.  |    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<  |    if (t<100.0) {  | ||
     setTimeout(loop,10);  |      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();  | |||
   }  | |||
}  | |||
</jsxgraph>  | |||
===The underlying JavaScript code===  | ===The underlying JavaScript code===  | ||
<source lang="javascript">  | <source lang="javascript">  | ||
var brd = JXG.JSXGraph.initBoard('box', {  | var brd = JXG.JSXGraph.initBoard('box', {axis: true, boundingbox: [-5, 1.2, 100, -1.2]});  | ||
var S = brd.  | var S = brd.create('turtle',[],{strokeColor:'blue',strokeWidth:3});  | ||
var I = brd.  | var I = brd.create('turtle',[],{strokeColor:'red',strokeWidth:3});  | ||
var R = brd.  | 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:'β'});  | |||
var s = brd.  | brd.create('text', [40,-0.4, "β: infection rate"]);  | ||
brd.  | var gamma = brd.create('slider', [[0,-0.5], [30,-0.5],[0,0.3,1]], {name:'γ'});  | ||
var beta = brd.  | brd.create('text', [40,-0.5, "γ: recovery rate = 1/(days of infection)"]);  | ||
brd.  | |||
var gamma = brd.  | var t = 0; // global  | ||
brd.  | |||
brd.  | brd.create('text', [40,-0.2,    | ||
         function() {return "  |          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.  |    I.setPos(0,s.Value());  | ||
   delta =   |    delta = 1; // global  | ||
   t =   |    t = 0;  // global  | ||
   loop();  |    loop();  | ||
}  | }  | ||
function turtleMove(turtle,dx,dy) {  | function turtleMove(turtle,dx,dy) {  | ||
   turtle.  |    turtle.moveTo([dx+turtle.X(),dy+turtle.Y()]);  | ||
}  | }  | ||
function loop() {  | function loop() {  | ||
   var dS = -beta.Value()*S.  |    var dS = -beta.Value()*S.Y()*I.Y();  | ||
   var dR = gamma.Value()*I.  |    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<  |    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 14: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]\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
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();
   }
}