SIR epidemic model

A simple mathematical description of the spread of a disease in a population is the so-called SIR model, which divides the (fixed) population of N individuals into three "compartments" which may vary as a function of time, t:

S(t) are those susceptible but not yet infected with the disease;
I(t) is the number of infectious individuals;
R(t) are those individuals who have recovered from the disease and now have immunity to it.

The SIR model describes the change in the population of each of these compartments in terms of two parameters, β and γ. β describes the effective contact rate of the disease: an infected individual comes into contact with β·N other individuals per unit time (of which the fraction that are susceptible to contracting the disease is S/N). γ is the mean recovery rate: that is, 1/γ is the mean period of time during which an infected individual can pass it on.

The differential equations describing this model were first derived by Kermack and McKendrick [Proc. R. Soc. A, 115, 772 (1927)]:

dS / dt = - β·S·I / N
dI / dt = β·S·I / N - γ·I
dR / dt = γ·I

The following Python code integrates these equations for a disease characterised by parameters β = 0.2, 1/γ = 10days in a population of N=1000 (perhaps flu in a school). The model is started with a single infected individual on day 0: I(0)=1.



In [11]:

    
%matplotlib inline
import numpy as np
from scipy.integrate import odeint
import matplotlib.pyplot as plt

# Total population, N.
N = 1000
# Initial number of infected and recovered individuals, I0 and R0.
I0, R0 = 1, 0
# Everyone else, S0, is susceptible to infection initially.
S0 = N - I0 - R0
# Contact rate, beta, and mean recovery rate, gamma, (in 1/days).
beta, gamma = 0.2, 1./10 
# A grid of time points (in days)
t = np.linspace(0, 200, 200)

# The SIR model differential equations.
def deriv(y, t, N, beta, gamma):
    S, I, R = y
    dSdt = -beta * S * I / N
    dIdt = beta * S * I / N - gamma * I
    dRdt = gamma * I
    return dSdt, dIdt, dRdt

# Initial conditions vector
y0 = S0, I0, R0
# Integrate the SIR equations over the time grid, t.
ret = odeint(deriv, y0, t, args=(N, beta, gamma))
S, I, R = ret.T

# Plot the data on three separate curves for S(t), I(t) and R(t)
plt.figure(figsize=(16, 4), dpi=72)
plt.plot(t, S/1000, 'b', label='Susceptible')
plt.plot(t, I/1000, 'r', label='Infected')
plt.plot(t, R/1000, 'g', label='Recovered with immunity')
plt.xlabel('Time (days)')
plt.ylabel('Percentage of individuals')
plt.xlim(0,200)
plt.ylim(0,1)
plt.legend(loc='center right', shadow=True)
plt.show()

References

Hill, Christian. "The SIR epidemic model." in Hill, Christian. Learning scientific programming with Python. Cambridge University Press, 2016.
Kermack, William O., and Anderson G. McKendrick. "A contribution to the mathematical theory of epidemics." Proceedings of the Royal Society of London A: mathematical, physical and engineering sciences. Vol. 115. No. 772. The Royal Society, 1927.