传染病数学模型
病毒扩散与传播的控制模型
摘要
本文基于传统的传染病模型，以微分方程的方法作为理论基础，结合采取的措施不同的情况，用 MATLAB 软件拟合出患者人数与时间的曲线关系，从中得出应采取的相应的应对措施。
在考虑地区总人数不变，人群被分为五类：确诊患者、疑似患者、治愈者、死亡和正常人，再将这几类分为可传染性和不可传染性两种。我们找出 单位时间内正常人数的变化、单位时间内潜伏期病人数的变化、单位时间内确诊患者人数的变化、单位时间内退出的人数的变化、单位时间内疑似患者人数的变化等关系 建立微分方程模型，得到病毒扩散与传播的控制模型。
在此基础上，我们将所要求的问题带入模型得到 患者人数随时间变化的曲线图，根据这图形得出模型结果的变化。 这样一来就可根据这结果的变化得出相应的应对措施。
此外对该传染病的潜伏期及治愈期进行了灵敏度分析，发现潜伏期的变化会对整个模型的结果产生较大影响，而治愈期的变化只会使传染病的持续时间缩短，但对累积的患病人数影响不大。
应尽量避免患者与正常人接触，减少正常人患病的可能性；加大隔离措施强度；减少拖延患者去住院的时间，让患者及时住院治疗。养成良好的卫生习惯，保证科学睡眠，适当锻炼，减少压力，保证营养，增强个人抵抗力，降低被病毒感染的危险。
关键词： 病毒扩散与传播 微分方程模型 曲线拟合
一、问题重述
已知某种不完全确知的具有传染性病毒的潜伏期为天，病患者的治愈时间为天。该病毒可通过直接接触、口腔飞沫进行传播、扩散，该人群的人均每天接触人数为。为了控制病毒的扩散与传播将该人群分为五类：确诊患者、疑似患者、治愈者、死亡和正常人，可控制参数是隔离措施强度（潜伏期内的患者被隔离的百分数）。
要求：
1 、在合理的假设下试建立该病毒扩散与传播的控制模型；
2 、利用你所建立的模型针对如下数据进行模拟
条件1：；
条件2：已经知道的初始发病人数为890、疑似患者为2000；
条件3：隔离措施强度；
条件4：患者2天后入院治疗，疑似患者2天后被隔离,试给出患者人数随时间变化的曲线图，并明确标识图中的一些特殊点的具体数据，分析结果的合理性。
3 、若将2中的条件4改为条件：患者1.5天后入院治疗，疑似患者1.5天后被隔离,模拟结果有何变化？
4 、若仅将2中的条件3改为条件：隔离措施强度，模拟结果有何变化？
5 、若仅将2中的条件1改为条件：，模拟结果有何变化？
6 、分析问题中的参数对计算结果的敏感性。
7 、针对如上数据给政府部门写一个不超过400字的建议报告。
二、问题分析
在考虑地区总人数不变，人群被分为五类：确诊患者、疑似患者、治愈者、死亡和正常人，我们可知，治愈者、死亡和正常人不可能传染病毒，我们把问题转化为如何找出正确的关系表达式来表达出每天病人增加的总数的问题，找出 单位时间内正常人数的变化、单位时间内潜伏期病人数的变化、单位时间内确诊患者人数的变化、单位时间内退出的人数的变化、单位时间内疑似患者人数的变化等关系 建立微分方程模型，得到病毒扩散与传播的控制模型。
在问题一的已求出得到病毒扩散与传播的控制模型的基础上，我们将后几问所要求的问题带入模型就可得到 患者人数随时间变化的曲线图，我们可以根据这些图得出模型结果的变化。 这样一来就可根据这些模型结果的变化得出相应的应对措施。
三、模型假设
1 、将病毒所有传播途径都视为与病原的接触；
2 、在疾病传播期间内所考察地区的总人数 N 视为常数，即认为本地区流入的人数与流出的人数相等，时间以天为计时单位；
3 、该病毒处于潜伏期的病毒不具有传染性；
4 、治愈者二度感染的概率为 0 ，他们以退出传染体系，因此将他们归为“退出者”；
5 、不考虑这段时间内人口出生率和自然死亡率，而对于由病毒引起的死亡人数，也将其归为“退出者”；
6 、被隔离的人群完全断绝与外界的接触，不再具有传染性；
7 、不考虑被隔离而实际又未被感染者，因为这部分人没有自由活动，对疾病的传播（感染和被感染）基本不造成任何影响；
8 、将人群分为以下四类：
正常人：易感染者
确诊患者：传染者
退出者：“治愈者”和“死亡者”统称；
疑似患者：被隔离但还没有确诊或者排除的人员；
四、符号约定
：确诊病人；
：潜伏期病人（感染了但处于潜伏期没有传染性的人)；
：类似病人（症状类似感染但其实没有感染的人)；
：退出者（痊愈和死亡的确诊病人)；
：普通易感者；
：病人的传染系数；
：潜伏期病人的传染系数（假设潜伏期病人也有传染性，但小于）；
：传染性病毒的潜伏期；
：病患者的治愈时间；
：该人群的人均每天接触人数；
：可控制参数是隔离措施强度（潜伏期内的患者被隔离的百分数）；
五、模型的建立与求解
5 ．1 模型一
5.1.1 模型准备
根据人口守恒原理,可建立如下模型：
模型将疫区的总人口数看成不变(不考虑流动人口) ,将疫区所有的人(假设人口的自然出生率和死亡率在疫期相等)分为:
：确诊病人
：潜伏期病人（感染了但处于潜伏期没有传染性的人)
：类似病人（症状类似感染但其实没有感染的人)
：退出者（痊愈和死亡的确诊病人)
：普通易感者
5.1.1 .1 单位时间内正常人数的变化：
……………⑴
5.1.1 .2 单位时间内潜伏期病人数的变化：
…………⑵
5.1.1 .3 单位时间内确诊患者人数的变化：
………………………⑶
5.1.1 .4 单位时间内退出的人数的变化：
……………………………………⑷
5.1.1 .5 单位时间内疑似患者人数的变化：
……⑸
其中，，， ， 为初值
(1) 、传染病毒的平均潜伏期为，即单位时间内潜伏期病人以比例常数，转为感染者；
(2) 、确诊病人平均死亡或痊愈的疗程为 ,即单位时间感染者的恢复率为；
(3) 、疑似患者平均疗程为，即单位时间内疑似患者的恢复率为；
(4) 、单位时间内每个易感者与病人的接触率参数为；
(5) 、易感者与疑似患者的接触率参数；
(6) 、考虑疑似患者感染病菌转为潜伏期病人,但潜伏期病人不会转为疑似患者；
(7) 、为隔离措施强度；
(8) 、为痊愈的被解除隔离的疑似患者；
(9) 、、和 均为被隔离对象；
(10) 、为疑似患者；
初值的设定: （这些数据是一个根据人口总数和医学常识的估计值。）
，
千万,不考虑流动人口;
；
；
；
参数的设定:
，设传染病平均潜伏期为5;
， 设确诊病人平均死亡或痊愈的疗程为20;
，设疑似病人平均疗程为20;
，疑似病人与易感者的接触率参数也假设固定。
5.1.2 模型的建立
， 千万,不考虑流动人口，，，。
5.2 模型二
5.2.1 模型建立
当，患者两天后入院治疗、疑似患者两天后被隔离时。
这样可以得到患者人数随时间变化的曲线图：（如下）
5.2.2 结果分析
从图中我们可以看出患者人数随时间变化先是急剧升高，这说明这是病毒传播初期的发展趋势，然后可以看到最高点(第12.37天)时患者人数达到最大值6803000人，通过采取患者入院治疗和疑似患者的隔离措施，我们从图中可以明显看出患者的人数呈下降趋势，并且在100天后患者人数降低到540800人。
5.3 模型三
5.3.1 模型建立
当，患者1.5天后入院治疗、疑似患者1.5天后被隔离时。
这样可以得到患者人数随时间变化的曲线图：（如下）
5.3.2 结果分析
从图中我们可以看出在最高点(第12.39天)时患者人数达到最大值6776000人，通过采取患者入院治疗和疑似患者的隔离措施，我们从图中可以明显看出患者的人数呈下降趋势，并且在100天后患者人数降低到516500人。
与问题二相比较我们可以知道提前入院治疗可以减少患者的人数，同时也可以更好的控制病毒的传播，更好的预防正常人患病。
5.4 模型四
5.4.1 模型建立
当，患者两天后入院治疗、疑似患者两天后被隔离时。
这样可以得到患者人数随时间变化的曲线图：（如下）
5.4.2 结果分析
从图中我们可以看到最高点(第12.12天)时患者人数达到最大值6802000人，通过采取患者入院治疗和疑似患者的隔离措施，我们从图中可以明显看出患者的人数呈下降趋势，并且在100天后患者人数降低到540800人。
与问题二相比较我们可以知道降低隔离措施会增加患者的人数，同时会对控制病毒的传播带来负面影响，会导致更多的正常人患病。所以我们建议医院加大病毒的控制强度。
5.5 模型五
5.5.1 模型建立
当，患者两天后入院治疗、疑似患者两天后被隔离时。
这样可以得到患者人数随时间变化的曲线图：（如下）
5.5.2 结果分析
从图中我们可以看到最高点(第12.48天)时患者人数达到最大值6803000人，通过采取患者入院治疗和疑似患者的隔离措施，我们从图中可以明显看出患者的人数呈下降趋势，并且在100天后患者人数降低到540800人。
与问题二相比较我们可以知道确诊患者的人均接触人数增大时，患病高峰期会延迟，同时会对控制病毒的传播带来负面影响，会导致更多的正常人患病。所以我们建议这些确诊患者减少与外界正常人的接触，这样减少正常人患病的可能性。
5.6 模型五
5.6.1 模型建立
通过建立的模型，我们对问题二、三、四、五进行了定量的计算，问题二、三、四、五是改变模型中的一些参数的值，得到不同参数下的结果，并分析了参数的改变对患者数量最大值，达到这个最大值的时间以及疫情得到控制的时间的变化，通过对这些数据结果的分析，可以得到某一参数的改变对病毒传播过程的影响，通过数据前后的对比，可以分析参数对计算结果的敏感性，针对这个问题，通过以上几问的求解，可以得到以下表格（表 1 ）：
表 1 各个参数对应的数值
| 问题 | 最大值时间 | 患病人数最大值 | 隔离措施强度 P | 患者入院前天数 n | 人均每天接触人数 R |
|---|---|---|---|---|---|
| 第二问 | 12.374 | 6803000 | 0.6 | 2 | 10 |
| 第三问 | 12.39 | 6776000 | 0.6 | 1.5 | 10 |
| 第四问 | 12.12 | 6802000 | 0.4 | 2 | 10 |
| 第五问 | 12.48 | 6803000 | 0.6 | 2 | 250 |
（ 1 ）、对患者 m 天后入院治疗的灵敏度分析
通过问题二与问题三的解答可知，我们通过改变患者入院治疗的时间，即将 m 分别取值为 2 和 1.5 天，然后得到患者人数随时间变化曲线（见图 1 和 2 ），通过分析对比可知：当 m=2 时，患者人数的最高峰达到了 6803000 人；当 m=1.5 时，患者人数的最高峰达到了 6776000 人。因此通过控制患者发病后入院治疗的时间就可适当的减少病情持续时间。
（ 2 ）、 对隔离强度 p 的灵敏度分析
通过对问题二与问题四的图象的观察，我们通过改变隔离措施强度 分析数据的变化，当 P 减小后患者人数最大值相比第二问来说增加了，同时达到这个最大值的时间也增大了，我们还可以看出病情消退的时间也稍稍的增长了，这说明隔离措施 P 减小时患者人数的最高峰增大了，同时达到这个数值的时间也相应的增大了，所以我们应尽量地增大隔离措施强度 P ，这样来椒使病毒消退时间减小，所以政府应尽量地增大隔离措施 P
（ 3 ）、对人均每天接触人数 R 的 灵敏度分析
通过对问题二与问题五的图像观察，我们通过改变人均接触人数 R 分析数据的变化，当 R 增大时患病人数最大值在增加，同时达到这个最大值的时间也在增加，病毒的消退时间也相应地增加，所以我们应尽量控制患者人均接触人数，应尽量地让大家少接触人，这样控制了人均接触人数，病毒的消退时间也会大大减小，疫情容易得到控制
5.7 病毒扩散与传播控制的建议报告
随着社会的进步，科学技术的发展，传统的传染病得到了有效的控制，但同时新发的传染病却不断出现，对于新发的某种不完全确知的具有传染性的病毒的突袭，我们首先要了解该病毒的传播方式，做好相应的防范措施。
通过本模型结果可知，被感染的人数有很大的差距，r越大被感染的人就越多，所以我们应该尽量避免与病人接触，因而要尽量少去人多的地方。
隔离措施强度相比较可知，p越小病情很越难控制，所以政府要加强隔离措施强度。且要减少拖延患者去住院的时间，让患者及时住院治疗。而且减少患者与外界人的接触，减少正常人患病的可能性。
最根本的实质还是平时要养成良好的卫生习惯，保证科学睡眠，适当锻炼，减少压力，保证营养，增强个人抵抗力，才可以降低被病毒感染的危险。
六、模型的评价与推广
6.1 模型的评价
6.1.1 模型的优点
⑴ 、将医学领域的问题转化到数学领域上进行分析和讨论，可以定量地得出传染病的发展趋势以及对未来的预测结果，具有很强的理论性和可靠性。
⑵ 、模型中涉及到的参变量都有相应的数据来源，结合一定的数据可以很方便计算出，而且各变量间关系明确，易于模型的求解。
⑶ 、由于本文的数学模型基于连续的微分方程，不会得出准确的解析解，我们在合理参数确定的前提下，将参数进行拟合，准确模拟出传染病发展趋势走向的曲线，从宏观的角度上给社会一个明晰的概念，易于被社会接受，对政府和人们采取措施起到了指导作用，具有一定的实用价值和直观性。
⑷ 、针对传染病对我们各方面的影响，我们给政府相应的政策，在现实生活中很具有实用性。
6.1.2 模型的不足
⑴ 、采用微分方程方法建立数学模型，易受外界因素变化的影响，其稳定性具有相对性，这就提出了外界干扰对该模型的影响程度的研究，从而建立传染病模型的稳定性理论，这点需在模型推广中进一步讨论。
⑵ 、模型中的参数变量有其自身的随机性，虽然我们采取对已知数据进行统计平均的处理方法，但在计算结果上仍存在一定的误差。
⑶ 、模型中涉及到的参数较多，在实际生活当中很难确定各参数。
6.2 模型的推广
我们建立该传染病模型的方法和思想对其他类似的问题也很适用，可广泛应用于人口、交通、肿瘤、战争、几何、物理、化学、体育、社会、经济等方面。对于微分模型受外界因素的干扰情况，我们可以借助一个称为李雅普诺夫函数的辅助函数和扰动微分方程所计算出来的全导数的符号性质来直接推断方程组的稳定性问题，亦称为李雅普诺夫直接法。
七、参考文献
[1] 姜启源谢金星 叶俊，《数学模型》（第三版），北京：高等教育出版社， 2003。
[2] 甘筱青，《数学建模教育及竞赛》，江西：江西高校出版社，2004。
[3] 赵静 但琦，《数学建模与数学实验》（第二版），北京：高等教育出版社，2003。
[4] 吴建国，《数学模型案例精编》，北京：中国水利水电出版社，2005。
八、附件
问题二程序：
function x=ill(t,x)
%s1=x(1) e=x(2) i=x(3) r=x(4) a=x(5);
m1=10;
m2=1*10.^(-11);
w=0.6;
d3=30;
d2=11;
d1=1;
x=[-m1x(3)(1-w)x(1)-m2(x(5)(1-w)+x(5)w1/d3)x(1),m1x(3)(1-w)(x(1)+x(5)(1-w)+x(5)w1/d3)-2/(d1+d2)x(2),2/(d1+d2)x(2)-1/d3x(3),1/d3x(3),m2*(x(5)(1-w)+x(5)w1/d3)x(1)-m1x(3)(1-w)(x(5)(1-w)+x(5)w1/d3)]';
s0=[10000000,500,890,0,2000];
[t,x]=ode23s(@ill,[0,100],s0)
plot(t,x( : ,3));
hold on
text(0,890,' (0,890)','color','r')
text(12.37,6.803E+006,' (12.37,6.803E+006)','color','r')
text(74,5.408E+005,' (100,5.408E+005)','color','r')
plot(0,890,'g+',12.37,6.803E+006,'g+',100,5.408E+005,'g+')
问题三程序：
function x=ill(t,x)
%s1=x(1) e=x(2) i=x(3) r=x(4) a=x(5);
m1=10;
m2=1*10.^(-11);
w=0.6;
d3=30;
d2=11;
d1=1;
x=[-m1x(3)(1-w)x(1)-m2(x(5)(1-w)+x(5)w1/d3)x(1),m1x(3)(1-w)(x(1)+x(5)(1-w)+x(5)w1/d3)-2/(d1+d2)x(2),2/(d1+d2)x(2)-1/d3x(3),1/d3x(3),m2*(x(5)(1-w)+x(5)w1/d3)x(1)-m1x(3)(1-w)(x(5)(1-w)+x(5)w1/d3)]';
s0=[10000000,500,890,0,2000];
[t,x]=ode23s(@ill,[0,100],s0)
plot(t,x( : ,3));
hold on
text(0,890,' (0,890)','color','r')
text(12.39,6.776E+006,' (12.39,6.776E+006)','color','r')
text(80,5.165E+005,' (100,5.165E+005)','color','r')
plot(0,890,'g+',12.39,6.776E+006,'g+',100,5.165E+005,'g+')
问题四程序：
function x=ill(t,x)
%s1=x(1) e=x(2) i=x(3) r=x(4) a=x(5);
m1=10;
m2=1*10.^(-11);
w=0.4;
d3=32;
d2=11;
d1=1;
x=[-m1x(3)(1-w)x(1)-m2(x(5)(1-w)+x(5)w1/d3)x(1),m1x(3)(1-w)(x(1)+x(5)(1-w)+x(5)w1/d3)-2/(d1+d2)x(2),2/(d1+d2)x(2)-1/d3x(3),1/d3x(3),m2*(x(5)(1-w)+x(5)w1/d3)x(1)-m1x(3)(1-w)(x(5)(1-w)+x(5)w1/d3)]';
s0=[10000000,500,890,0,2000];
[t,x]=ode23s(@ill,[0,100],s0)
plot(t,x( : ,3));
hold on
text(0,890,' (0,890)','color','r')
text(12.12,6.802E+006,' (12.12,6.802E+006)','color','r')
text(74,5.408E+005,' (100,5.408E+005)','color','r')
plot(0,890,'g+',12.12,6.802E+006,'g+',100,5.408E+005,'g+')
问题五程序：
function x=ill(t,x)
%s1=x(1) e=x(2) i=x(3) r=x(4) a=x(5);
m1=250;
m2=1*10.^(-11);
w=0.6;
d3=32;
d2=11;
d1=1;
x=[-m1x(3)(1-w)x(1)-m2(x(5)(1-w)+x(5)w1/d3)x(1),m1x(3)(1-w)(x(1)+x(5)(1-w)+x(5)w1/d3)-2/(d1+d2)x(2),2/(d1+d2)x(2)-1/d3x(3),1/d3x(3),m2*(x(5)(1-w)+x(5)w1/d3)x(1)-m1x(3)(1-w)(x(5)(1-w)+x(5)w1/d3)]';
s0=[10000000,500,890,0,2000];
[t,x]=ode23s(@ill,[0,100],s0)
plot(t,x( : ,3));
hold on
text(0,890,' (0,890)','color','r')
text(12.48,6.803E+006,' (12.48,6.803E+006)','color','r')
text(74,5.408E+005,' (100,5.408E+005)','color','r')
plot(0,890,'g+',12.48,6.803E+006,'g+',100,5.408E+005,'g+')