疫情 | 用python仿真SIR模型 正确认识传染病

这两个月来世界人民都在关注肺炎疫情，新冠病毒的传播时刻牵动着大家的神经，一时之间传染病让人们谈虎色变，其实早在上个世纪，就已经开始了对传染病的研究和分析，SIR模型作为早期提出来的经典模型一直作为人们分析传染病早期的工具。

这个公式就是病毒传播的SIR模型，对比疫情，其中参数的意思是

It ：t 时刻感染的人数

Pd  ： 感染率

Pr   ： 治愈率

N : 可以被感染的人群

St ： 可能被感染但还未感染的人数

而其中的 Pd / Pr 就是一直被提起的R0，基础繁殖数。理论上R0 > 1，这个病毒就可以传播到整个人群中；R0 < 1，病毒就会慢慢自行消亡。S曲线就是扩散的传播曲线。病毒主要靠人传人，在早期感染人数很少时传播速度很慢，但随着感染人数达到一定规模，就会启动增速飞轮，这时候疫情就很难控制了，但随着时间增加，治愈人数的增多会使得扩散到一定值后放缓，整个过程是开始和结束慢，中间快。

SIR模型是在SI模型上的改进，但引入了R，R不仅指治愈的人，也指死亡的人，这两种情况造成的影响都是相同的，都是去除（remove）的效果，整个SIR模型是建立在几个理想化的条件下的：

1、总人数不变，忽略出生于正常死亡

2、假设 t 时刻单位时间内，一个病人能传染的易感者数目与此环境内易感者总数s(t)成正比，设定比例系数为β，从而在t时刻单位时间内被所有病人传染的人数为βs(t)i(t)

3、t 时刻，单位时间内从染病者中移出的人数与病人数量成正比，比例系数为γ，单位时间内移出者的数量为γi(t)

4、移除者不再被感染

一个简单的SIR模型如下，通过scipy中的odeint完成微分方程的求解，diff_eqs为对应的SIR模型的微分方程，它所对应的微分方程为

import scipy.integrate as spi
import numpy as np
import pylab as pl
alpha=1
beta=0.1
#观察时间
TS=1.0 
#观察结束时间
ND=50.0 
#初始易感人数
S0=1-1e-6 
#初始感染人数
I0=1e-6 
INPUT = (S0, I0, 0.0)
def diff_eqs(INP,t):
    Y=np.zeros((3))
    V = INP
    Y[0] = - alpha * V[0] * V[1]
    Y[1] = alpha * V[0] * V[1] - beta * V[1]
    Y[2] = beta * V[1]
    return Y

t_range = np.arange(0.0, ND+TS, TS) 
#odeint 数值求解微分方程
RES = spi.odeint(diff_eqs,INPUT,t_range)

pl.subplot(111)
pl.plot(RES[:,0], '-og', label='Susceptible')
pl.plot(RES[:,1], '-r', label='Infective')
pl.plot(RES[:,2], '-k', label='Removal')
pl.legend(loc=0)
pl.title('SIR')
pl.xlabel('Time')
pl.ylabel('Numbers')
pl.xlabel('Time')
pl.show()

程序运行结果为

可以看到即使在初始感染人数如此小的情况下，只要经过几天就会快速积累起庞大的患者，这对于医疗系统和医疗资源都是巨大的挑战，这也是为什么钟南山院士一直提出的早发现、早隔离、早诊断、早治疗。

除了SIR模型之外，还有几个其他经典的模型，如SIRS（康复者的免疫力是暂时的，康复者会转化为易感者）、SEIR（加入传染病潜伏期）、SEIRS（同时有潜伏期且免疫暂时的条件）等模型，可以去阅读Albert-Laszlo Barabasi的Network Science中传染病模型的一章，除了这种模型之外，还有利用网格研究病毒在城市中的蔓延情况的模型，可以看一下http://www.45fan.com/article.php?aid=20021636528147093950613304，介绍的非常详实且通俗易懂，介绍了在多种干预措施下的影响，如减少公共交通、隔离场所等。新冠病毒在我们国家的有利措施下 必然会被消灭，就像钟南山院士说的武汉是一座英雄的城市，中国也是一个英雄的国家。