|
|
由于编者水平以及时间的原因，本文必定存在不合理的地方，还望知友批评指正。
2019-nCoV新型冠状病毒传染分析(数学建模)
电子显微镜下的2019-nCoV
摘要
本文首先采用抽样检测法对2019-nCoV早期的模型的合理性及实用性进行了评价，然后我们通过对传染病的共性及2019-nCoV的特性的分析。得出三个基本假设并且把人群理想化为三类(S类, I类, R类)，建立起基本的SIR模型，再对SIR模型中的三类人群间的相互转化关系的分析，结合马氏链得出三种人群间变化率的矩阵T，由于2019-nCoV的特性，可知SIR模型中的两个参数a(t), b(t)是以时间为变量的函数。我们根据武汉疫情的数据，通过多现实的数据拟合法分别得到a(t), b(t)及T结合，从而建立出模型。由于医疗条件的逐步改善，一定会制定出一套有效的治疗方案，甚至到后期的2019-nCoV疫苗的研发。于是我们在不改变人群分类的情况下，增加了一个系数c(c表示有效治疗方案的成功率,考虑到人群年龄问题以及由于肺炎引起的并发症导致的死亡, 故c为常数)，来进一步完善此模型。
本文利用数学软件(MATLAB)很好的实现了模型运算，并结合实际数据得出了人群与实践的关系图。从图中可以很好的反映出各类人群的变化规律，他们的变化规律与实际变化相吻合，从而证明了我们的模型基本符合要求。
一、 问题的提出
2019新型冠状病毒，即“2019-nCoV”，因2019年武汉病毒性肺炎病例而被发现，2020年1月12日被世界卫生组织命名。目前疫情仍未得到控制，并且已经对美好的春节造成了巨大的破坏。由于在家闲来无事，因此，有必要根据2019-nCoV流行的特点，建立数学模型预测其传染，从而采取措施预防和控制其发展。而建立该模型我们要综合各方面的因素才能使模型合理化。
二、 问题的分析
主要通过分析湖北(武汉)地区的受感染人数的变化规律，我们对该地区预测流行病的变化趋势提出以下模型假设：
1.将人群分为三类：
易感染人数：用S表示；
病人数(已受感染者，即确诊者)：用I表示；
移出者人数(包括“被治愈者”和“死亡者”)：用R表示。
2.该地区人口不流动(考虑到湖北已经封省，该假设是合理的)，设最初易感染人数为N，此时I，R均为0.
3.被隔离人群完全断绝与外界接触，不再具有传染性(考虑到现有的医疗条件，该假设也是合理的)。
三、模型的分析与建立
(1)初期数据的模拟
传染病早期可以采用指数模型进行模拟: N(t)=N_{0}(1+k)^{t} ，当然，此情况是在社会来不及防备以及群众不重视的基础上导致的。由于前期武汉市政府的不重视、人民对本次疫情的忽视以及春节带来的附加作用，我们可以认为，新型冠状病毒前期是可以满足该指数模型的。
湖北省人数统计
湖北省肺炎疫情人数预测图
由MATLAB拟合出的人数曲线如上所示。湖北未在23日提供确诊人数，故将该点去除后拟合。此外，MATLAB给出的拟合结果如下：
General model: f(x) = a*(1+b)^x
Coefficients (with 95% confidence bounds):
a = 0.7421 (-0.05313, 1.537)
b = 0.3169 (0.2601, 0.3738)
Goodness of fit:
SSE: 6.771e+04
R-square: 0.949
Adjusted R-square: 0.9468
RMSE: 54.26
显然，此模型只能适用于早期的预测。(此模型对于SARS早期的预测具有极高的吻合度，但当疫情进入中期之后，差别会越来越大)
(2)模型的建立
建立 SIR 模型。易感染者，感染者，移出者之和是个恒量即 N=S+I+R 。假设病人康复后具有免疫力，人与人之间有相同的接触率。最终由如下两种假设决定状态之间的转变率：
① 感染者的增长率是和感染者 I 与易感染者 S 的乘积成正比的
② 感染者 I 到移出者 R 的变化率是与感染者成正比。
基于以上两条得出模型的微分方程(具体细节容易理解)：
\begin{equation} \left\{ \begin{array}{**lr**} \frac{dS}{dt}=-aSI & \\ \frac{dI}{dt}=aSI-bI &\\ \frac{dR}{dt}=bI & \end{array} \right. \end{equation}
其中a，b都是以时间为变量的参数， a(t) 为日感染率， b(t) 为日移出率。
但是因为此疫情目前仍处于上升期，且相关数据较少，故按照以上微分方程组无法求出
a(t),b(t) 的解析解，因此我们先作数值计算。
参考多方资料后，我们设 a=0.0000003,b=0.0077266 , I(0)=1，S(0)=1000000 (假设其中十分之一的人口受到此次疫情的影响)
用MATLAB进行编程
MATLAB代码
上图是由MATLAB求解微分方程后得出的结果。可以看到，21天到25天的数据，也就是截止到1月26日24时，预测的数据都是符合实际情况的。
1/22-1/26湖北省感染人数统计
但是，预测数据给出的结果显然是不符合实际情况的，随着疫情的扩张，感染率势必降低，移出率势必提高。因此，感染率a和移出率b不会是一个常数，该模型仍然有需要改进的地方。
(3)模型的改进与分析
1.相轨线分析
首先定义 \sigma=\frac{a}{b} ，将原微分方程组化简：
\frac{dI}{dS}=\frac{1}{\sigma*S}-1,I|_{S=S_{0}}=I_{0}
容易求出此微分方程的解为：
I=(S_{0}+I_{0})-S+\frac{1}{\sigma}ln\frac{S}{S_{0}}
显然当 S=\frac{1}{\sigma} 时 I 最大。对于基本模型来说，其 \frac{1}{\sigma}=\frac{0.0077266}{0.0000003}\approx25755
由相图可知， 1/\sigma 所代表的是一个阈值，当 S>1/\sigma 时，传染病在蔓延; 当 S<1/\sigma 时，传染病则不会蔓延。
2.模型的改进
由(二)分析可知，此模型对于前期疫情的预测比较准确，但是对于后期疫情的发展则显然不符合实际结果。其原因是因为感染率a和移出率b不会是一个常数。
参考与本次2019-nCoV疫情非常相似的2003年非典疫情的数据可知，后期传染率应呈指数型下降。并将移出率用sigmoid函数进行优化。
① 首先保持移出率不变，优化感染率函数
令 a=0.0000003((stepfun(t,0)-stepfun(t,25))+stepfun(t,25)*exp(-0.02(t-25)))
感染率曲线(第25天时武汉市采取响应措施，故从该日起感染率下降)
由改进后的模型建立的预测图，前期预测仍然符合预期，且大大缓解了高峰时期患病的人数。
下面改变感染率函数下降的斜率，由MATLAB拟合出结果
增大感染率下降的梯度，可以有效的降低感染人数。
②在①的基础上改进移出率函数
a=0.0000003((stepfun(t,0)-stepfun(t,25))+stepfun(t,25)*exp(-0.05(t-25)))
b=0.0077266(stepfun(t,0)+stepfun(t,25)*(5/(1+exp(-t+28))))
移出率b是参考SARS后期治愈率即死亡率的数据之后得出的结论。
在改进模型的情况下，我们发现，相比于①的结果，肺炎疫情在60天左右即得到了有效的控制，且感染人数呈数十倍的下降。最后总的预测感染人数会在5万左右。当然，预测的数据仍然具有一定的问题，稍后将在模型问题中具体阐述。改进后的模型考虑了疫情控前和控后的情况，具体在参数a与参数b中体现，因此更加符合真实情况。
③灵敏度分析
对a做灵敏度分析
如图可见，感染人数对a的变化还是比较灵敏的。这对于现实中的疫情防控具有很好的指导作用，尤其是控后，如何快速的降低感染率以便快速的控制疫情是防止疫情蔓延的重点。
对b做灵敏度分析
同样，如何提高移出率(主要是指治愈率)也是疫情防控后的重点。
四、模型存在的问题
(1)问题分析
虽然改进后的模型能够对疫情做出比较理想的趋势分析，但是对于疫情后期的处理仍然与真实情况有所偏差。主要是由三个原因导致的，一是SIR模型过于精简，将真实情况过度理想化。本次疫情的感染不光是由感染者传染的，对于一部分易感人群，只要携带病原体，均可以作为传染源，但是该部分人群在模型的建立当中并没有体现出来。二是本次疫情仍处于上升期，且真实数据不足，加上湖北省政府前期的怠慢导致的数据异常，我们没有办法对控前做一个很好的拟合，同样，对控后亦无从而知。三是武汉市人口基数众多，且正逢春节，按照官方的说法，在武汉封城之前，有约500万人离开武汉。这对于模型的建立是一个极大的挑战，因为一般的模型都要求人口固定且无人口交流。
下面主要从模型建立的角度阐述本案例存在的问题。参考由SARS疫情建立的比较合理的数学模型。
(2)比较合理的模型建立过程
①基本假设
1) 假设所考察人群的总数恒定，且无病源的输入和输出
2) 将所考察人群分为正常人(易感人群)、确诊患者(传染者)、退出者(治愈和死亡)、疑似患者(被隔离但是还没有确诊的人群)
3) 假设已治愈的患者不会二次感染
4) 假设所有患者均为“他人输入型患者”，即不考虑个体自身发病
5) 假设各类人群在人群中总体分布均匀
6) 假设已被隔离的人群之间不会发生交叉感染
7) 不考虑隐性的病毒携带者，即只要携带病毒即可发病
②符号约定
I ：确诊病人
E ：潜伏期病人(感染了但处于潜伏期且具有传染性的人)
R ：退出者
S ：普通易感者
a_{1} ：病人的传染系数
a_{2} ：潜伏期病人的传染系数
d_{1}-d_{2} ：病毒的潜伏期
d_{3} ：病患者的治愈时间
r ：该人群的人均每天接触人数
p ：可控制参数是隔离措施强度(潜伏期内的患者被隔离的百分数)
\omega ：确诊患者中仍有机会能够传染的人数比例
③模型建立
\begin{equation} \left\{ \begin{array}{**lr**} \frac{dS}{dt}=-a_{1}{\omega}Ir-a_{2}(1-p)Er & \\ \frac{dE}{dt}=a_{1}{\omega}Ir+a_{2}(1-p)Er-\frac{2}{d_{1}+d_{2}}E & \\ \frac{dI}{dt}=\frac{2}{d_{1}+d_{2}}E-\frac{1}{d_{3}}I &\\ \frac{dR}{dt}=\frac{1}{d_{3}}I & \end{array} \right. \end{equation}
④模型分析
相比于之前的模型，此模型使用了更多的参数，符合真实情况中的疫情传播。
四种人群转换示意图
目前官方统计的四项数据
可见，此模型对真实情况具有良好的匹配性。但由于此前数据不足，不能够做出完整的预测，该模型适用于疫情结束之后的分析，有助于政府、医疗机构、人群充分认识疫情，并能够总结经验教训，在下次疫情出现的时候可以迅速做出反应。
五、总结
截止2020/1/27 19:37，武汉肺炎疫情全国总计确诊2840例，疑似5794例，死亡81例，治愈55例。且1月26日的疑似病例仍处于快速增长状态，形式可谓相当严峻，本文主要从数学模型的角度分析了疫情的大致走势，帮助大家对新型冠状病毒疫情有更好的认识。相信我们能够在党中央的带领下，在一线医生的努力下战胜2019-nCoV!
六、引用
[1] 姜启源 谢金星 叶俊，《数学模型》(第三版)，北京：高等教育出版社，2003
[2] Mark M. Meerschaert，《数学建模方法与分析》(原书第四版)，机械工业出版社，2015
[3] 百度文库：https://wenku.baidu.com/view/c565c516fc4ffe473368ab1d.html
[4] 百度文库：https://wenku.baidu.com/view/76544379e009581b6bd9eba2.html
<hr/>注：对本文进行了部分修改，于2020/2/6日完成文稿编辑，不再以图片形式展示。
<hr/>修改
- 感谢知友@堂庭多棪树的指正。原文表述中，将疑似病例当做易感人群 S 不妥，故删除“疑似病例”的表述。
说明
- 相关代码在文中已经给出，需要的请细看文章后自取。
- word和pdf文档在网盘中被识别为违法信息，无法再提供相关文档。
度盘和谐 |
|
发表于 2023-3-31 21:36:17