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 : 'β' }); 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 (); } }

See also