注解

此笔记本可从此处下载： Epidemic_model_SIR.ipynb

流行病模型SIR

import numpy as np
import matplotlib.pyplot as plt
import matplotlib as mpl
import pandas as pd
import numba
from scipy import integrate
%matplotlib notebook

所谓的SIR模型描述了疾病在固定为 \(N\) 随着时间的推移，个人 \(t\) 。

问题描述

中国的人口 \(N\) 个人分为三类(舱室)：

  • 个人 \(S\) 易受感染的；

  • 个人 \(I\) 感染；

  • 从疾病中康复的个人 \(R\) (并且现在已经获得了对它的豁免权)；

哪里 \(S\) ， \(I\) 和 \(R\) 是的函数 \(t\) 。

SIR模型与流行病学中的许多其他分室模型一样，依赖于我们现在介绍的特定参数：

  • \(\beta>0\) 疾病的收缩率(传播参数)

  • \(\gamma>0\) ：平均回收率

个人 \(S\) 在与一种病毒正面接触后被感染 \(I\) 个人的。然而，他对这种疾病产生了免疫力。 : he leaves \(I\) 舱室在a处 \(\gamma\) 治愈率。

模型功能

该模型基于以下假设：

  1. 平均而言，一个人 \(S\) 在人口接触中 \(\beta\) 单位时间内的个人

  2. 被感染的人离开车厢的比率 \(I\) 是 \(\gamma I\) 每单位时间(一旦一个人被感染，他就会对这种疾病产生免疫力)。

  3. 人口规模 \(N = S + I + R\) 是不变的。

这是模型的方程式系统：

\[\begin{split}\Left\{\Begin{array}{ll} \displaystyle\frac{\mathm{d}s}{\mathm{d}t}&=-\displaystyle\frac{\beta S}{N}i，\\ [4mm] \displaystyle\frac{\mathm{d}i}{\mathm{d}t}&=\displaystyle\frac{\beta S}{N}i-\Gamma i，\\ [4mm] \displaystyle\frac{\mathrm{d}R}{\mathrm{d}t}&=\DisplayStyle\Gamma I。 \end{array}\右。\end{split}\]

关于该模型的一些说明：

  • 平均感染期(即感染者可以将其传染的平均时间)等于 \(\displaystyle \frac{{1}}{{\gamma}}\) 。

  • 这是一个确定性模型

  • 平均接触数恒定的假设 \(\beta\) 是一个强有力的约束性假设：它不能适用于所有疾病。

  • 我们可以想象通过考虑以下因素来改进此模型：

    • 新生的婴儿将会对应于 \(S\) 易受影响的个体。我们会引入出生率 \(b\) 。

    • 会离开车厢的死者 \(S\) 或 \(I\) 以同样的速度 \(b\) : this allows to consider a constant population \(N\) 。

有关SIR型号的进一步阅读，请访问此处：https://en.wikipedia.org/wiki/Compartmental_models_in_epidemiology#The_SIR_model

模型图

先生，下图总结了我们研究的模型。

图像0

Source: Wikipedia

问题表述

我们摆姿势

\[\begin{split}X=\BEGIN{pMatrix}S\\ 我\\ r\end{pMatrix}\end{split}\]

我们把我们的微分系统写成：

\[\begin{split}\点{X}= \BEGIN{pMatrix}\DisplayStyle-\frac{\beta S}{N}I\\ \displaystyle\frac{\beta S}{N}i-\Gamma I\\ \γI \end{pMatrix} =f\Left(X\Right)\end{split}\]

设置

N = 350. #Total number of individuals, N
I0, R0 = 1., 0 #Initial number of infected and recovered individuals
S0 = N - I0 - R0 #Susceptible individuals to infection initially is deduced
beta, gamma = 0.4, 0.1 #Contact rate and mean recovery rate
tmax = 160 #A grid of time points (in days)
Nt = 160
t = np.linspace(0, tmax, Nt+1)

def derivative(X, t):
    S, I, R = X
    dotS = -beta * S * I / N
    dotI = beta * S * I / N - gamma * I
    dotR = gamma * I
    return np.array([dotS, dotI, dotR])

X0 = S0, I0, R0 #Initial conditions vector
res = integrate.odeint(derivative, X0, t)
S, I, R = res.T
Seuil = 1 - 1 / (beta/gamma)
Seuil

0.75

plt.figure()
plt.grid()
plt.title("odeint method")
plt.plot(t, S, 'orange', label='Susceptible')
plt.plot(t, I, 'r', label='Infected')
plt.plot(t, R, 'g', label='Recovered with immunity')
plt.xlabel('Time t, [days]')
plt.ylabel('Numbers of individuals')
plt.ylim([0,N])
plt.legend()

plt.show();

<IPython.core.display.Javascript object>

欧拉法

def Euler(func, X0, t):
    dt = t[1] - t[0]
    nt = len(t)
    X  = np.zeros([nt, len(X0)])
    X[0] = X0
    for i in range(nt-1):
        X[i+1] = X[i] + func(X[i], t[i]) * dt
    return X

Nt = 100
Xe = Euler(derivative, X0, t)

plt.figure()

plt.title("Euler method")
plt.plot(t, Xe[:,0], color = 'orange', linestyle = '-.', label='Susceptible')
plt.plot(t, Xe[:,1], 'r-.', label='Infected')
plt.plot(t, Xe[:,2], 'g-.', label='Recovered with immunity')
plt.grid()
plt.xlabel("Time, $t$ [s]")
plt.ylabel("Numbers of individuals")
plt.legend(loc = "best")

plt.show();

<IPython.core.display.Javascript object>

两种方法的比较

plt.figure()

plt.plot(t, Xe[:,0], color = 'orange', linestyle = '-.', label='Euler Susceptible')
plt.plot(t, Xe[:,1], 'r-.', label='Euler Infected')
plt.plot(t, Xe[:,2], 'g-.', label='Euler Recovered with immunity')
plt.plot(t, S, 'orange', label='odeint Susceptible')
plt.plot(t, I, 'r', label='odeint Infected')
plt.plot(t, R, 'g', label='odeint Recovered with immunity')
plt.grid()
plt.xlabel("Time, $t$ [s]")
plt.ylabel("Numbers of individuals")
plt.legend(loc = "best")
plt.title("Comparaison of the two methods")
plt.show();

<IPython.core.display.Javascript object>

是什么 \(R_0\) ？

\(R_0\) 是描述由于患病个体而新感染的平均数的参数。它通常被称为 基本再生数 。这是流行病学中的一个基本概念。

如果 \(R_0 > 1\) 这种流行病会持续下去，否则就会绝迹。

如果一种疾病有一种 \(R_0 = 3\) 例如，平均而言，一个患有这种疾病的人会把它传染给其他三个人。

有关更多详细信息，请参见：https://en.wikipedia.org/wiki/Basic_reproduction_number

import ipywidgets as ipw

def deriv(X, t, beta, gamma):
    S, I, R = X
    dSdt = -beta * S * I / N
    dIdt = beta * S * I / N - gamma * I
    dRdt = gamma * I
    return np.array([dSdt, dIdt, dRdt])

def update(beta = 0.2, gamma = 0.1): #beta = 0.3, gamma = 0.06 -> R0 = 5
    """
    Update function.
    """
    I0, R0 = 1, 0
    S0 = N - I0 - R0
    X0 = S0, I0, R0

    sol = integrate.odeint(deriv, X0, t, args = (beta, gamma))
    line0.set_ydata(sol[:, 0])
    line1.set_ydata(sol[:, 1])
    line2.set_ydata(sol[:, 2])

    txR0 = beta/gamma
    ax.set_title("$N$ = {0} $R_0$ = {1:.2f}".format(str(N).zfill(4), txR0))

Nt = 160
t = np.linspace(0., tmax, Nt)
X = np.zeros((Nt,3))
fig = plt.figure()
ax = fig.add_subplot(1, 1, 1)
plt.grid()
line0, = ax.plot(t, X[:,0], "orange", label = "Susceptible")
line1, = ax.plot(t, X[:,1], "r", label = 'Infected')
line2, = ax.plot(t, X[:,2], "g", label = 'Recovered with immunity')
ax.set_xlim(0., tmax, Nt)
ax.set_ylim(0., N, Nt)
ax.set_xlabel("Time, $t$ [day]")
ax.set_ylabel("Percentage of individuals")
plt.legend()

ipw.interact(update, beta = (0.,1., 0.01),
           gamma = (0.01, 1., 0.01));

<IPython.core.display.Javascript object>

interactive(children=(FloatSlider(value=0.2, description='beta', max=1.0, step=0.01), FloatSlider(value=0.1, d…