1226-SIR epidemic model (Kermack-McKendrick).sagews

SIR epidemic model (Kermack-McKendrick)

In this model, there are three classes of individuals: Susceptible, Infected, and Recovered (or Removed). Susceptible individuals are infected at a rate of β times the number of infected individuals, I. Infected individuals recover at a rate of ν. Vital dynamics have been added using the parameter μ, which is the rate at which new susceptibles are introduced and susceptibles, infecteds, and recovereds are depleted (birth and death rate).

This demo plots the number of infected (red) and susceptible (blue) individuals over time, as well as a phase plot (fraction of susceptibles vs fraction of infecteds) in green. You can set β, ν, and the initial number of susceptible and infected individuals.

# SIR epidemic model (Kermack-McKendrick model)
import numpy
import scipy
from scipy.integrate import odeint
from sage.plot.line import Line

@interact
def euler_method(beta = input_box(2, label = 'β: '), nu = input_box(0.8, label = 'ν: '), mu = input_box(0, label = 'μ: '), tmax = input_box(20, label = 'max t: '), Ifrac = slider(0,1,0.05,label='Initial fraction infected: ', default=0.05)):
  html(r"$\displaystyle\frac{dS}{dt} = \mu - \mu S - \beta S I$")
  html(r"$\displaystyle\frac{dI}{dt} = \beta S I - \mu I - \nu I$")
  
  def dXdt(X,t):
    dX1dt = (mu - mu*X[0] - beta*X[0]*X[1])
    dX2dt = (beta*X[0]*X[1] - mu*X[1] - nu*X[1])
    return numpy.array([dX1dt, dX2dt], dtype=float)
    
  t = numpy.arange(0,tmax,0.1)
  pts = odeint(dXdt,[1.0-Ifrac, Ifrac],t)
  pts_S = [[t[i],pts[i][0]] for i in range(len(t))]
  pts_I = [[t[i],pts[i][1]] for i in range(len(t))]
  show(line(pts_S, rgbcolor=(0,0,1)) + line(pts_I, rgbcolor=(1,0,0)) + text("S", (0.1,1-Ifrac), rgbcolor=(0,0,1)) + text("I", (0.1,Ifrac), rgbcolor=(1,0,0)), axes_labels=['time', ''], ymax=1)

  html("<h2>Phase plane</h2>")
  P1 = line([(0,1),(1,0)], rgbcolor=(0,0,0))
  t = numpy.arange(0,10,0.1)
  for i in range(1,10):
      I0 = i/10.0
      S0 = 1.0-I0
      pts = odeint(dXdt,[S0, I0], t) 
      pts_xy = [[x[0],x[1]] for x in pts]
      P1 = P1+line(pts_xy, rgbcolor=(0,1,0))
      
  show(P1,xmin=0,ymin=0,xmax=1,ymax=1, axes_labels=['Susceptible fraction', 'Infected fraction'])

Interact: please open in CoCalc