传染病传播及预防的数学模型

摘要：

随着社会和经济的发展，医学水平能力渐渐得到提高，现今社会的医学水平已经能够有效地预防和控制许多传染病，但是仍然有一些传染病暴发或流行，危害人们的健康和生命。人们也认识到定量地研究传染病的传播规律、为预测和控制传染病蔓延创造条件的重要性。通过建立传染病的传播模型，可以了解传染病的扩散传播规律，为预测和控制传染病提供可靠、足够的信息。

传染病病毒是随时间演变的过程。本文以微分方程的SIR模型为基础，分析传染病的扩散传播规律，建立动态模型。应用传染病动力学模型来描述疾病发展变化的过程和传播规律，预测疾病发生的状态，评估各种控制措施的效果，为预防控制疾病提供最优决策依据, 维护人类健康与社会经济发展。通过人数的规划，建立了传染病的微分方程模型，并用matlab软件拟合出患者人数随着时间的变化的关系曲线，利用控制变量的方法，控制某些变量不变，改变其中某个变量，通过比较找出导致传染病的传染的主要因素，以便做出相应的措施。

本模型的关键在于把确诊患者、疑似患者、治愈者、死亡和正常人划分成可传染者和不可传染者两类人，辅加一些特殊的参数，如：传染率，治愈率等等，构成微分方程组，找出单位时间内正常人人数的变化，确诊患者人数的变化，疑似患者人数的变化，死亡者或治愈者（即退出系统者）的人数的变化，从而建立了微分方程模型。

在模型建立的基础上，通过matlab软件拟合出患者人数随时间变化的曲线关系图，分析图形，得出结果，从而找到解决问题的响应措施。

关键词：动力学模型 微分方程模型 控制变量 matlab软件

一、问题重述

已知某种不完全确知的具有传染性病毒的潜伏期为d1~d2到，病患者的治愈时间为d3天。该病毒可通过直接接触、口腔飞沫进行传播、扩散，该人群的人均每天接触人数为r。为了控制病毒的扩散与传播将该人群分为五类：确诊患者、疑似患者、治愈者、死亡和正常人，可控制参数是隔离措施强度p（潜伏期内的患者被隔离的百分数）。通过合理的假设建立传染病传播的数学模型。

二、问题分析

据题目意思，这是一个传染性病毒随着时间演变的过程，我们要分析、预测、研究它就得建立动态模型，在此我们选用微分方程。因题目中把人群分为五类：确诊患者、疑似患者、治愈者、死亡和正常人，所以我们采用SIR模型。模型中我们找出单位时间内这五类人群人数的变化来建立微分方程，得出模型。再利用matlab画出图形，加以分析，达到得出应对措施的目的。

把考察范围内的人群分为以下种类：

1、健康人群，即易感染（Susceptibles）人群。记其数量为S(t),表示t时刻未感染病但有可能感染该疾病的人数；

2、潜伏期人群，即被感染（Infection）该疾病的人群，记其数量为I(t) 表示t时刻可能感染该疾病的但又不是疑似病患的人数； 3、疑似病患，记其数量为E(t) 表示示t时刻感染该疾病的并是疑似病患的人数； 4、确诊病患,记其数量为Q(t) 表示示t感染该疾病并确诊为患者的人数；

5、恢复人群（Recovered）,记其数量为R(t)，表示t时刻已从感染病者中移出的人数（这部分人数既不是已感染者，也不是非感染者，不具有传染性，也不会再次被感染，他们已经推出了传染系统）。

基于以上的假设，健康人群从潜伏期到移出传染系统的过程图如下：

正常人(易 感染人群) 社会交往 感染 该疾病 疑似患者 确诊患者 医务人员感染 恢复人群 治愈 医院 死亡

三、 模型假设

1. 假设易感人数的变化率与当时的易感人数和感染人数的乘积成正比; 2. 假设从感染数中移除个体的速率与当时的感染人数成正比;

3. 假设考察地区内疾病传播期间忽略人口的出生，死亡，流动等种群动力因素对总人数的影响。即：总人口数不变,记为N。

4. 假设潜伏期人群不会传染健康人，不具有传染性。

5. 假设被隔离的患者无法跟别人接触，不会传染健康人。

6. 假设治愈者已对该病毒有免疫力，不会再被该传染病传染，可以退出系统 7. 假设初始时刻健康人群的总人数为S0=1.1千万,潜伏期的总人数为I0=1，疑似病患的总人数为E0=0，确诊病患的总人数为Q0=0，恢复人群的总人数为

R0=0。

四、符号说明

病毒潜伏期（天） d1~d2 病患者治愈时间（天） d3 病患人均每天接触人数 r 隔离措施强度 p 时刻t内健康人群 S(t) 时刻t内潜伏期人群 I(t) 时刻t内病症疑似人群 E(t) 时刻t内已患病人群 Q(t) 时刻t内治愈或死亡人群 R(t) 传染病传染率 ?

五、建立模型

由模型的假设得到如下关系：S(t)+I(t)+E(t)+Q(t)+R(t)=N

1）根据假设在时刻?t内健康人群变化有：S(t??t)?S(t)???Q(t)(1?p)S(t)?t 2）在时刻?t内治愈或死亡人群的变化有：R(t??t)?R(t)?位时间内患者的恢复率） 3）在时刻?t内病

11为单I(t)?t（

d3d3症疑似人群的变化有：

E(t??t)?E(t)???Q(t)(1?p)[E(t)(1?p)?E(t)p1]?t d34）在时刻?t内已患病人群的变化有（已患病人群等于潜伏期病人转为感染者减去移除人数）：Q(t??t)?Q(t)?5）在

时

刻

?t21I(t)?Q(t)

d1?d2d3内潜伏群期人群的变化有：

12I(t??t)-I(t)?{?Q(t)(1-p)[S(t)?E(t)(1-p)?E(t)p]-I(t)}?td3d1?d2

（

2为单位时间内潜伏期病人转为感染者的比例常数）

d1?d2根据以上变化有

dS????Q(t)(1?p)S(t)?dt?dR1??I(t)dtd3??dE1????Q(t)(1?p)[E(t)(1?p)?E(t)p] ?dtd3?dQ21??I(t)?Q(t)?dtd1?d2d3?dI12???Q(t)(1?p)[S(t)?E(t)(1?p)?E(t)p]?I(t)?d3d1?d2?dt六、模型的求解与验证

模型一分析：

当d1?1,d2?14,d3?30,r?20,p?60%，患者2天后入院治疗，疑似患者2天后被隔离。有初始状态的患者人数为：Q?Q(0)*(者人数随时间变化如图一：

N?Q(0)?E(0)*r)*2，则患

Nt=250,y=115.5535

由上图可以得到：在当d1?1,d2?14,d3?30,r?20,p?60%，患者2天后入院治疗，疑似患者2天后被隔离的条件下。当t??0,13.0971?时，患者的人数是急剧上升的，在t=13.0971达到最大值，此时患者人数为y=228669.8722在采取医疗措施，比如患者入院治疗，隔离疑似患者等后患者人数随着时间的增长呈现下降的趋势，在250天后患者人数为115.5535。 模型二分析：

当d1?1,d2?14,d3?30,r?20,p?60%，患者1.5天后入院治疗，疑似患者1.5天后被隔离。有初始状态的患者人数为：Q?Q(0)*(则患者人数随时间变化如下图二：

N?Q(0)?E(0)*r)*1.5，

Nt=250,y=26.5958

由上图可得：在当d1?1,d2?14,d3?30,r?20,p?60%，患者1.5天后入院治疗，疑似患者1.5天后被隔离的条件下。当t??0,13.3715?时，患者的人数是急剧上升的，在t=13.3715达到最大值，此时患者人数为y=52851.3637在采取医疗措施，比如患者入院治疗，隔离疑似患者等后患者人数随着时间的增长呈现下降的趋势，在250天后患者人数为26.5958。 模型三分析：

当d1?1,d2?14,d3?30,r?20,p?40%，患者2天后入院治疗，疑似患者2天后被隔离。有初始状态的患者人数为：Q?Q(0)*(者人数随时间变化如下图三：

N?Q(0)?E(0)*r)*2，则患

Nt=250,y=126.5086

由上图可得：在当d1?1,d2?14,d3?30,r?20,p?40%，患者2天后入院治疗，疑似患者2天后被隔离的条件下。当t??0,14.324?时，患者的人数是急剧上升的，在t=14.324达到最大值，此时患者人数为y=228870.1396在采取医疗措施，比如患者入院治疗，隔离疑似患者等后患者人数随着时间的增长呈现下降的趋势，在250天后患者人数126.5086。

表1 各个参数对应的数值 问题 最大值 患病人数 隔离措患者入院人均每天接250天后患时间 最大值 施 前 触 病 强度P 天数n 人数r 人数 第二13.0971 228669.8722 0.6 2 20 115.5535 问 第三13.3715 52851.3637 0.6 1.5 20 26.5958 问 第四14.324 228870.1396 0.4 2 20 126.5086 问 从上表可以看出，1.当隔离强度一样的时，患者入院的开始时间将在一定时间内影响到患病人数。明显可以看出，患者2天后入院与1.5天后入院相比，患者的治疗时间延长了，而且患病的人数也增多。因此相关部门应及时将病人隔离并治疗。2.当患者入院开始时间一样时，隔离措施强度将影响到患病人数达到最大时的时间长短。可以看出，隔离措施强度降低后，患者人数相对偏高。因此相关部门应该加强隔离措施强度，提高警惕。从上述两个参数取值变化分析可知，“得病后入院时间”与“隔离措施强度”对于传染病疫情态势发展，具有很大的

敏感性与相关性，其中得病后的患者几时去医院治疗，对于疫情的控制具有更重要的意义。所以，“早发现、早隔离、早治疗”，能够帮助我们有效地、较快地控制传染病的扩散与传播。

当患者入院开始时间一样，隔离强度不一样时，患者人数随时间的变化如下图四：

当隔离强度一样的时，治疗时间不同的患者人数随时间的变化。如图五：

六．模型分析与评估

本模型中采用微分方程的模型，对传染病传播做出合理假设，并对其得过拟合，得出传染病的发展趋势，可以有效预报传染病高潮到来的时刻，对群众接受传染病的预防知识起到很好的警示作用。但模型中的参数都具有随机性，所以得出的结果还是会有误差的存在，不能准确预报每次传染病高潮的到来的准确时间，只能限定在一定时间内。

本模型建立的传染病模型因能有效地预报传染病高潮到来的时刻，所以在现实生活中可以得到进一步的应用，特别是对现今社会中的传染病的爆发跟流行有很好的控制作用。

七、模型应用

根据以上建立的模型可以得到： 病毒传播控制的建议和措施：

根据以上建模的结果分析，我们明确了要制止传染病的蔓延主要手段有：提高医院的医疗水平和卫生水平，加强医疗工作人员的效率，这样可以提高隔离强度，减少入院时间。对此，政府需要积极采取措施来控制传染病的传播，及早发现被传染的人员，将其隔离，切断传染病传播的途径。

可以采取的措施有： 1.控制传染源

防止传染源到处活动排出病原体传播他人，应该将他们隔离看护，直到完全康复。

2.切断传播途径

加强公共场所的管理（如公共厕所，商场，娱乐场所等）； 建筑物通风条件的改善; 3.保护易感人群

加强易感人群的个人卫生意识，号召他们接种疫苗，防止传染病的感染，提高自身的免疫力。

八、参考文献

【1】数学建模简明教程 戴朝寿 孙世良 编著

【2】数学建模方法及其应用 解放军信息工程大学 韩中庚

九、附录

模型一： 图一：

MATLAB的.m文件

function x=illness(t,x)

%S=x(1)I=x(2)Q=x(3)R=x(4)E=x(5); a1=0.1429; p=0.6; d3=30; d2=14;

d1=1;

x=[-a1*x(3)*(1-p)*x(1),a1*x(3)*(1-p)*(x(1)+x(5)*(1-p)+x(5)*p*1/d3)-2/(d1+d2)*x(2),2/(d1+d2)*x(2)-1/d3*x(3),1/d3*x(3),-a1*x(3)*(1-p)*(x(5)*(1-p)+x(5)*p*1/d3)]';

MATLAB源代码：

s0=[900*(0.9997*20)^2,500,900,0,2000] [t,x]=ode23s(@illness,[0,250],s0); [y_max,i_max]=max(x(:,3))

t_text=['t=',num2str(t(i_max))]; y_text=['y=',num2str(y_max)];

max_text=char('maximum',t_text,y_text);%生成标志最大值点的字符串 y=num2str(x(end,3)) plot(t,x(:,3)); hold on

plot(t(i_max),y_max,'r','MarkerSize',20); text(t(i_max)+0.3,y_max+0.05,max_text); xlabel('t'),ylabel('y'),hold off 模型二： 图二：

function x=ill3(t,x)

%S=x(1)I=x(2)Q=x(3)R=x(4)E=x(5); a1=0.1429; p=0.6; d3=30; d2=14; d1=1;

x=[-a1*x(3)*(1-p)*x(1),a1*x(3)*(1-p)*(x(1)+x(5)*(1-p)+x(5)*p*1/d3)-2/(d1+d2)*x(2),2/(d1+d2)*x(2)-1/d3*x(3),1/d3*x(3),-a1*x(3)*(1-p)*(x(5)*(1-p)+x(5)*p*1/d3)]';

MATLAB源程序：

s0=[900*(0.9997*20)^1.5,500,900,0,2000] [t3,x3]=ode23s(@ill3,[0,250],s0); [y_max3,i_max3]=max(x3(:,3))

t_text3=['t=',num2str(t3(i_max3))]; y_text3=['y=',num2str(y_max3)]; t_end3=['t=',num2str(t3(end))]; y_end3=['y=',num2str(x3(end,3))];

max_text3=char('maximum',t_text3,y_text3);%生成标志最大值点的字符串 y3=num2str(x3(end,3)) plot(t3,x3(:,3)); hold on

plot(t3(i_max3),y_max3,'r','MarkerSize',20); text(t3(i_max3)+0.3,y_max3+0.05,max_text3); xlabel('t'),ylabel('y'),hold off

模型三： 图三：

MATLAB中.m文件

function x=ill4(t,x)

%S=x(1)I=x(2)Q=x(3)R=x(4)E=x(5); a1=0.1429; p=0.4; d3=30; d2=14; d1=1;

x=[-a1*x(3)*(1-p)*x(1),a1*x(3)*(1-p)*(x(1)+x(5)*(1-p)+x(5)*p*1/d3)-2/(d1+d2)*x(2),2/(d1+d2)*x(2)-1/d3*x(3),1/d3*x(3),-a1*x(3)*(1-p)*(x(5)*(1-p)+x(5)*p*1/d3)]';

MATLAB源程序：

s0=[900*(0.9997*20)^2,500,900,0,2000] [t4,x4]=ode23s(@ill4,[0,250],s0); [y_max4,i_max4]=max(x4(:,3))

t_text4=['t=',num2str(t4(i_max4))]; y_text4=['y=',num2str(y_max4)];

max_text4=char('maximum',t_text4,y_text4);%生成标志最大值点的字符串 y4=num2str(x4(end,3)) plot(t4,x4(:,3)); hold on

plot(t4(i_max4),y_max4,'r','MarkerSize',20); text(t4(i_max4)+0.3,y_max4+0.05,max_text4); xlabel('t'),ylabel('y'),hold off 图四：

s0=[900*(0.9997*20)^2,500,900,0,2000] [t,x]=ode23s(@illness,[0,250],s0); [t4,x4]=ode23s(@ill4,[0,250],s0); plot(t,x(:,3),'k+',t4,x4(:,3),'g+'); 图五：

s0=[900*(0.9997*20)^2,500,900,0,2000] [t,x]=ode23s(@illness,[0,250],s0);

s3=[900*(0.9997*20)^1.5,500,900,0,2000] [t3,x3]=ode23s(@illness,[0,250],s3); plot(t,x(:,3),'k+',t3,x3(:,3),'g+');