2219-SIR Models.sagews

The SIR Model

Randall Pruim

Calvin College

This worksheet is based on a worksheet by Tim McLarnan and Anand Pardhanani (Earlham College).  You can search for Earlham in the published worksheets at sagemath.org/pub.  You will find several worksheets, including one about SIR models.

The Earlham worksheet shows how to use and ODE solver with the SIR model as an example.  I have made several changes to that code and also demonstrate how to fit the same model via stochastic simulation. 

SIR Models

SIR models are a type of compartment model describing the impact of an infectious disease on a population.  This model has three compartments: Susceptible, Infected, and Recovered.  Individuals can move from one compartment to another, and the rates at which they do so are determined by the proportion of the population in each compartment.  For example, when more individuals are infected, the rate of transfer from Susceptible to Infected (and also from Infected to Recovered) increases.  In fact our simple model is entirely driven by this transfer, since it is a model in which Recovered individuals are immune and so cannot become Infected or Susceptible again.

A Differential Equations Model

One of the simplest SIR models is the Kermack-McKendrick model, which is described by the following system of differential equations.

$$

\begin{array}{lcrcr}

S'(t)&=&-\beta \cdot S\cdot I \;  \\

I'(t)&=&\beta \cdot S\cdot I  &-& \gamma I \; \\

R'(t)&=& &+& \gamma I \; .

\end{array}

$$

This model says that Susceptibles become Infected at a rate proportional to the number of potential contacts (roughly the product of $S$ and $I$, and that Infecteds recover at a fixed rate $\gamma$.

A solution using Sage

The function sir_ode takes as arguments the initial numbers of susceptible ($S$), infected ($I$), and recovered ($R$) people, and the length of time over which you want data, and the model parameter ($\beta$ and $\gamma$).  The output is a plot showing $S$ as a solid black line, $I$ as dashed red line, and $R$ as a dash-dot blue line.

def sir_ode(S_init, I_init, R_init,    # initial values
              endtime,                   # time to model
              beta = 2/10^6, gamma=1/14, # model parameters
              colors=['black','red','blue'], thickness=3 # plotting controls
              ):    
    var('t S I R')
    system = [-beta*S*I, beta*S*I- gamma * I, gamma *I]

    # P will be a set of 4-tuples of the form (t,S,I,R)
    P = desolve_system_rk4(system,[S,I,R], ivar=t,
                            ics=[0, S_init, I_init, R_init], 
                            step=1, end_points=endtime)

    plotS = plot(line([(time,s) for time,s,i,r in P],
                 linestyle = '-', rgbcolor=colors[0]), thickness=thickness)
    plotI = plot(line([(time,i) for time,s,i,r in P],
                 linestyle = '--', rgbcolor=colors[1], thickness=thickness))
    plotR = plot(line([(time,r) for time,s,i,r in P], 
                 linestyle = '-.', rgbcolor=colors[2], thickness=thickness))
    return plotS+plotI+plotR

The Earlham Example

Here's an example in which we follow a population that starts with 45400 susceptible people, 2100 infected people, and 2500 recovered people over a period of 300 days using the default values of $\beta$ and $\gamma$.  This example was used in the orignal Earlham SIR worksheet.

sir_ode(45400, 2100, 2500, 300)

Other Examples

We can adjust any of the initial values or parameters of the model and try again.

sir_ode(45400, 2100, 2500, 300, beta=1/50000)

plot1 = sir_ode(45400, 2100, 2500, 150)
plot2 = sir_ode(45400, 2100, 2500, 150, beta=1/10^5, colors=['gray','purple','orange'])
show(plot1 + plot2)

A Stochastic Model

Alternatively, we implement a stochastic model.  In this model we use a binomial distribution to determine how many individualas move from Susceptible to Infected and from Infected to Recovered.

$$

\begin{align*}

\pi_{SI} &= \beta * I \\

\pi_{IR} & = \gamma

\end{align*}

$$

This model is equivalent to the differential equation model above.  The implementation is straightforward.

def sir_stoch(S_init, I_init, R_init,    # initial values
              endtime,                   # time to model
              beta = 2/10^6, gamma=1/14, # model parameters
              colors=['black','red','blue'], thickness=3 # plotting controls
              ):    
             
    from scipy import stats              # use scipy to generate random draws
    P = [(0,S_init, I_init, R_init)]     # a list with one tuple
    for t in range(endtime):
        (time, S, I, R) = P[-1]
        if S > 0 and I > 0: 
            S_to_I = stats.binom.rvs(S, beta * I, size=1)[0]
        else:
            S_to_I = 0
        if I >0:
            I_to_R = stats.binom.rvs(I, gamma, size=1)[0]
        else:
            I_to_R = 0
        P.append( (time+1, S - S_to_I, I + S_to_I - I_to_R, R + I_to_R) )
    
    plotS = plot(line([(time,s) for time,s,i,r in P],
                 linestyle = '-', rgbcolor=colors[0]), thickness=thickness)
    plotI = plot(line([(time,i) for time,s,i,r in P],
                 linestyle = '--', rgbcolor=colors[1], thickness=thickness))
    plotR = plot(line([(time,r) for time,s,i,r in P], 
                 linestyle = '-.', rgbcolor=colors[2], thickness=thickness))
    return plotS+plotI+plotR

sir_ode(15, 2, 0, 300,beta=1/50)

The function below runs both models and is a convenient way to make comparisons.

def sir_both(S_init, I_init, R_init, endtime, beta = 2/10^6, gamma=1/14):   
    p1 = sir_ode(S_init, I_init, R_init, endtime, beta, gamma, colors=['gray','purple', 'green'], thickness=2)
    p2 = sir_stoch(S_init, I_init, R_init, endtime, beta, gamma, thickness=3)
    return p1+p2

The differences are miniscule for a large population.

sir_both(45400, 2100, 2500, 300)

But in small populations, the law of averages doesn't help us out and the stochasitic nature of the second model is easier to see.

sir_both(40, 10, 0, 100,beta=1/500)

sir_both(40, 10, 0, 100,beta=1/500)

sum( [ sir_stoch(40,10,0,100,beta=1/500) for i in range(5) ] )

Reproducible "random" results

We can set the seed of the random number generator (located in the numpy module) if we want reproducable random results.

import numpy as np
np.random.seed(int(12345))
sir_both(40, 10, 0, 100,beta=1/500)

sir_stoch(45400, 2100, 2500, 300)