medium_modelling_part_one.ipynb

Kernel: Python 3 (ipykernel)

SIR Model

From: https://github.com/hf2000510/infectious_disease_modelling/blob/master/part_one.ipynb

LICENSE: https://github.com/hf2000510/infectious_disease_modelling/blob/master/LICENSE

---

Bibliography

An introduction to mathematical epidemiology, M Martcheva - 2015 - Springer [Google Scholar] [PDF]

from scipy.integrate import odeint
import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline 
#!pip install mpld3
import mpld3
mpld3.enable_notebook()

compartmental model.

A compartmental model separates the population into several compartments, for example:

• Susceptible (can still be infected, “healthy”)

• Infected

• Recovered (were already infected, cannot get infected again)

We require the following parameters

• $N$: total population

• $S(t)$: number of people susceptible on day t

• $I(t)$: number of people infected on day t

• $R(t)$: number of people recovered on day t

• $\beta$: expected amount of people an infected person infects per day

• $D$: number of days an infected person has and can spread the disease

• $\gamma=1/D$: the proportion of infected recovering per day ( = 1/D)

the basic reproduction number $R_0$, is the total number of people an infected person infects. We just used an intuitive formula: $$R_0 = \beta \cdot D\,.$$

Check $R_0$ for Colombia: https://github.com/restrepo/COVID-19/blob/master/covid.ipynb

See Understanding the models that are used to model Coronavirus: Explaining the background and deriving the formulas of the SIR model from scratch. Coding and visualizing the model in Python.

$$\begin{align} \frac{d S}{d t}=&-\beta \cdot I \cdot \frac{S}{N} \\ \frac{d I}{d t}=&\beta \cdot I \cdot \frac{S}{N}-\gamma \cdot I \\ \frac{d R}{d t}=&\gamma \cdot I \end{align}$$

So that ($R$ include recover or died people) $$N=S(t)+I(t)+R(t)$$

If $S\gg I$ and $S\gg R$, then $S/N\approx 1$ and only the second equation is relevant and reduce to the exponential growth $dI/dt=\beta I$.

def deriv(U, t, N=1000, β=1, γ=1/4):
    S, I, R = U
    dSdt = -β * S * I / N
    dIdt = β * S * I / N - γ * I
    dRdt = γ * I
    return [dSdt, dIdt, dRdt]

We want start with an infected, $I_0=I(0)=1$ and a $R_0=R(0)=1.2$. Then, by using $D=4$ $$\beta=R_0/D=0.3\,.$$

N = 1000
β = 0.3  # infected person infects 0.275 other person per day, such that R_0=1.2
D = 4.0 # infections lasts four days
γ = 1.0 / D

S0, I0, R0 = 999, 1, 0  # initial conditions: one infected, rest susceptible

t = np.linspace(0, 199, 200) # Grid of time points (in days)
U0 = [S0, I0, R0] # Initial conditions vector

# Integrate the SIR equations over the time grid, t.
ret = odeint(deriv, U0, t, args=(N, β, γ))
S, I, R = ret.T

def plotsir(t, S, I, R,yscale='linear'):
  f, ax = plt.subplots(1,1,figsize=(10,4))
  ax.plot(t, S, 'b', alpha=0.7, linewidth=2, label='Susceptible')
  ax.plot(t, I, 'y', alpha=0.7, linewidth=2, label='Infected')
  ax.plot(t, R, 'g', alpha=0.7, linewidth=2, label='Recovered')
  plt.yscale(yscale)
  #plt.grid()
  ax.set_xlabel('Time (days)')

  ax.yaxis.set_tick_params(length=0)
  ax.xaxis.set_tick_params(length=0)
  #ax.grid(b=True, which='major', c='w', lw=2, ls='-')
  legend = ax.legend()
  legend.get_frame().set_alpha(0.5)
  for spine in ('top', 'right', 'bottom', 'left'):
      ax.spines[spine].set_visible(False)
  plt.show();

plotsir(t, S, I, R,yscale='log')

Scaling:

Nreal=50E6
sc=Nreal/N
sc

50000.0

print(f'N0={int(I0*sc)}')
print(f'Nmax={int(I.max()*sc)}')

N0=50000
Nmax=778250

I.max()

26.67215538853332

5E4,5*5E4

(50000.0, 250000.0)