流行病毒的扩散与传播的控制问题

关于病毒扩散与传播的控制模型 流行病毒的扩散与传播的控制问 一、摘要 本文是针对流行病毒的扩散与传播的控制问题所进行的研究。以猪流感为例，通过建立微分方程来讨论和研究这个问题。最近几天猪流感在墨西哥爆发，引起全世界人的关注。流行病毒的扩散与传播的控制问题得到各国领导人和世界卫生组织的重视。各国都采取各种措施预防猪流感病毒的传播和蔓延。 在问题1中，我们把人群分为五类：健康者，确诊病人、疑似病人、退出者（包括治愈者和死亡者）、自由带菌者，根据流行病毒的发展情况，将会建立两个模型，其中包括流行病毒被控制前和被控制后的模型。 在问题2中，我们利用在问题1中所建立的第二个控制后模型，对满足问题2中的四个条件的模型进行求解，并给出患者人数随时间变化的曲线如图3，并标示图中的一些特殊点的具体数据，通过图像可以看到经过3天后疫情 基本得到控制，12天后无确诊病人，10天后无自由带菌者和疑似病人。 在问题3中，我们在问题2的基础上，对问题2中的第四个条件进行调整，之后就进行模拟求解，并给出患者人数随时间变化的曲线如图4，通过图像可以看到经过3天后疫情 基本得到控制，12天后无确诊病人，10天后无自由带菌者和疑似病人。并且病人达到高峰期时的数量明显低于问题2的数量。 在问题4中，同样是在问题2的基础上，对问题2中的第三个条件进行调整并进行求解，并给出患者人数随时间变化的曲线如图5，通过图像可以看到经过2天后疫情 基本得到控制，11天后无确诊病人，8天后无自由带菌者和疑似病人，并且病人达到高峰期时的数量明显低于问题2的数量。 问题3的结果与问题4相比，问题3所提供的措施在控制病人增长数量上效果比问题4的措施佳，问题4在控制病情所用的时间方面比问题3所提供的措施效果好。 给政府的建议：尽早提出解决治疗强化确诊病人和疑似病人的隔离面对突发性流行病毒的袭击，应当以预防为主，早发现、早治疗、早疏散、早隔离。 关键字：流行病毒扩散、微分方程、控制模型、敏感度 二、问题重述 最近几天猪流感在墨西哥爆发，引起全世界人的关注。流行病毒的扩散与传播的控制问题得到各国领导人和世界卫生组织的重视。各国都采取各种措施预防猪流感病毒的传播和蔓延。假设该病毒的潜伏期为 至 天，得病患者经治疗经过d3天可以治愈，严重的可能引起患者死亡。该病毒可通过直接接触、口腔飞沫进行传播、扩散。设人群中每人每天的接触人数为r。人群中的人可以分为5类：确诊患者、疑似患者、治愈者、死亡人和正常人，可控制参数是隔离措施强度p，即潜伏期内的患者及疑似患者被隔离的百分数。 1.​ 建立流行病病毒扩散与传播的控制模型； 2.​ 利用所建立的模型针对如下数据进行模拟 条件1. =2, =7, =20, =15; 条件2. 已经知道出事发病人数为50人，疑似患者280人； 条件3. 隔离强度p=75%； 条件4. 患者2天后入院治疗；疑似患者2天后被隔离. 试给出患者人数随时间变化的曲线图，并标识图中的一些特殊点的具体数据，结果的合理性. 3.​ 若将2中条件4改为患者1天后入院治疗，疑似患者1天后隔离，模拟结果如何变化？ 4.​ 若将2中条件3改为隔离强度p=90%模拟结果如何变化？ 分析问题中参数对计算结果的敏感性. 针对以上数据给政府部门写一个不超过800字的建议报告. 三、问题分析 针对问题1，为了建立流行病病毒扩散与传播的控制模型，我们把人群分为五类，如下图所示： 结合问题的要求，我们采用微分的思想流行病毒控制前和流行病毒控制后的模型。在流行病毒被控制之前，是按照自然规律扩散与传播的，患者的人数可以视为常数。而流行病毒被控制后，患者人数与隔离措施强度有关。而问题2、3、4是在问题1所建立的第二个模型基础上针对一些数据进行模拟求解，从而找到一个可行的方案来控制疫情。 四、模型的建立与求解 模型的建立 1、控制前阶段（模型Ⅰ） 模型假设 为了模型建立的方便我们提出了2条假设： 假设1 控制前在疾病传播期内所考察地区的总人数N不变，既不考虑生死，也不考虑迁移，不考虑诊治，将人群分为三类：健康者（S’）、疑似病人（R’）、发病者（I’）； 假设2 与发病者接触过的健康者人都会成为疑似病人，由假设人群每人每天平均接触r个正常人，则假设每日新增疑似病人= 假设3随着每天增加发病人也越来越多，其概率分布呈指数稳步增长，则每天有 概率的人变为病患。 从而建立以下微分方程 （程序见附录A）【1】 ： 2、控制后阶段（模型Ⅱ） 模型假设 为了模型建立的方便我们提出了4条假设： 假设1 在疾病传播期内所考察地区的总人数N不变，既不考虑生死，也不考虑迁移。将人群分为五类：健康者，确诊病人、疑似病人、退出者、自由带菌者。 假设2 为了计算方便，假设自由带菌者一经发现就被确定为确诊病人，而不必隔离观察；而被自由带菌者感染的人被发现必须隔离观察，即归为疑似病人；未被发现则归为自由带菌者。 假设3 隔离人群完全断绝与外界接触，不再具有传染性。 假设4 治愈后则产生抗体，不会再被传染，因此，模型中不予考虑。 符号约定 t：时间变量，以天计算 d1：传染性病毒的潜伏期下限 d2：传染性病毒的潜伏期上限 d3：病患者的治愈时间 r：病人人均每天接触r个人 p：每日新增被感染人群中被隔离的概率（即被隔离的强度） S(t): t时刻健康者人数（易感染者） I(t)：t时刻确诊病人人数（己被隔离的病人） R(t)：t时刻疑似病人人数（当日没有发病的隔离者） T(t)：t时刻退出者（包括“死亡者”和“治愈者”） D(t)：t时刻自由带菌者（没有收治的传染病患者） b：疑似病人中每日被排除的人数占疑似病人总数比例，一般每日疑似病人被排除的数量呈指数减少，所以可以假设 d：每日自由带菌者被确诊为病人的概率，一般自由带菌者变为确诊病患的数量呈指数增长，所以可以假设 f:：病人人均每天接触的健康者的人数，f与p有关，设 q：退出率（日死亡率和日治愈率之和），当新增病人越来越多，则每天退出的人也将增多，设 a：疑似病人中每日被确诊人数占疑似病人总数比例,一般每日疑似病人变为确诊病患的数量呈指数增长，所以可以假设 我们可以根据人群的实际流向可以得出以下数学模型：（程序见附录C）【2】 模型求解 （问题二） 根据所建的模型1和模型Ⅱ，针对如下数据进行模拟 条件1. =2, =7, =20, =15; 条件2. 已经知道出事发病人数为50人，疑似患者280人； 条件3. 隔离强度p=75%； 条件4. 患者2天后入院治疗；疑似患者2天后被隔离. 根据条件4，知在出事前两天，没有控制措施，因此我们在对病人人数在前两天发展的情况，用控制前阶段的模型（即模型1）来模拟，得到两天后的疑似患者人数与自由带菌人数分别58621和15963。 （程序见附录C） 图1观察0-2天内人数的变化 根据条件4，知在出事两天后，采取了控制措施，因此我们在对病人人数在两天后发展的情况，用控制后阶段的模型（即模型Ⅱ）来模拟，模型赋的初值是： 运用MATLAB软件求解（程序见附件1），得到患者人数随时间变化的曲线图为：【3】 第3天病人数量达到高峰，12天后，无确诊病人 第2天病人数量达到高峰，10天后，无疑似病人 第2天病人数量达到高峰，10天后，无自由带菌者。 图3 患者人数随时间变化 通过图像可以看到经过3天后疫情 基本得到控制，12天后无确诊病人，10天后无自由带菌者和疑似病人。 此图为在开始作出控制后所作的图，由图可知，首先，可能由于刚开始控制病情，控制措施不佳，所以即前3天左右，病毒的传播速度依然很快。而后有政府迅速找到有效措施，并及时采取，因而该病毒开始得到的控制。 这个图像，符合现实，因此结果合理。 （问题三） 在问题2的基础上，我们将2中条件4改为患者1天后入院治疗，疑似患者1天后隔离， 意思就是政府在病发前1天没有采取任何措施，因此在病情发展的前1天，用控制前阶段的模型（即模型1）来模拟，得到1天后的疑似患者人数与自由带菌人数分别2870和661。 （程序见附录D） 图2 1天内患者人数变化 根据修改后的条件4，知在出事1天后，采取了控制措施，因此我们在对病人人数在两天后发展的情况，用控制后阶段的模型（即模型2）来模拟，模型赋的初值是： 得到患者人数随时间变化的曲线图为： 第3天病人数量达到高峰，12天后，无确诊病人 第2天病人数量达到高峰，10天后，无疑似病人 第2天病人数量达到高峰，10天后，无自由带菌者 图4 患者人数随时间的变化 与问题二进行对比可知，患者越早入院治疗，疑似患者越早被隔离，病情越早得到控制，在第10天的时候病情就已经得到了控制。在病人达高峰时期，病人的数目也有大幅度减少。 （问题四） 在问题2的基础上，我们将将2中条件3改为隔离强度p=90%，得到患者人数随时间变化的曲线图为 第2天病人数量达到高峰，11天后，无确诊病人 1天内病人数量达到高峰，8天后，无疑似病人 1天内病人数量达到高峰，8天后，无自由带菌者 图5 患者人数随时间的变化 由上图可知，当把隔离强度增加到p=90%时，这两天之内，病情就开始得到有效地控制，患者人数逐渐减少，相对问题2中的结果知，在病人达高峰时期，病人的数目也有大幅度减少，并且完全控制用的时间也有所减少。 问题3的结果与问题4相比，问题3所提供的措施在控制病人增长数量上效果比问题4的措施佳，问题4在控制病情所用的时间方面比问题3所提供的措施效果好。 五、模型评价 模型优点 （1）​ 模型1，只考虑了三类人，生病的人得不到治疗，那最终将几乎所以的人都会感染到这种病，这显然不符合实际，因为现实中一旦发现有某种病流行，政府将会采取措施控制病情，并有研制医疗这种病的药物，使病情得到控制，但是在短时期内能用这个模型，因为当一种病刚产生时，政府可能做出的反应没有那么快，来不及采取措施，医院也可能在短时期内没有医疗的方案，所以在短时期内，病人可能显著增加。 （2）​ 模型2是在模型1的基础上建立起来的，模型2是病情开始控制后的模型，所以，人群分的种类就有所增加，把人群分的更细，这种分法也是比较合理的，模型的个股变量关系明确，并且可以比较方便的得到。 （3）​ 模型能够反映政府对传染病的态度和控制力度，如果可能甚至可能进一步考虑控制力度和费用的函数关系，方便政府疫情控制及成本等方面进行更好的决策 模型缺点 （1）​ 模型1和模型2都忽略了流动人口的影响，忽略了交通工具上的传播，实际上从全世界各地区的猪流感疫情来看，交通中的传播有一定的影响。 （2）​ 模型1、2采用常微分方程进行分析，自身就有一定的缺陷：微分方程对流行病的情况而言，其计算结果的准确性、可靠性将受到一定的限制。 六、分析问题中参数对计算结果的敏感性. 根据各个条件分析可知，控制前过程主要受到参数r的影响，而控制后过程则主要由政府采取措施的时t和隔离措施强度p所控制，当控制的时间越早，疫情被控制得越好，当控制时间越晚，疫情被控制得越差，问题3相对问题2只是提早了一天采取措施，高峰期的病人数目大幅度减少。当p越大时，疫情控制得越快，当p的力度不够时，疫情控制得越慢，并且p的改变不但使高峰期病人数目减少，而且控制病情用的时间也减少了。 七、针对以上数据给政府部门写一个不超过800字的建议报告. 随着社会的进步、卫生设施和医疗水平的改善、人类文明的不断进步，虽然诸如霍乱、天花等曾经肆虐全球的传染病已经得到了有效的控制，但是在世界的某些地区，特别是在较贫困的发展中国家，很可能出现传染病流行的情况。因此，找出一个能够预测并为预防和控制提供足够可靠信息的方法是人们所期望的，建立传染病数学模型可以预测、预防以及控制传染病的流行，阻止其蔓延。 传染病模型是根据传染病传播的一般机理，利用以观测到的数据，对传染病的传播规律进行模拟，从而找出他的传播方式和途径，预测在以后较短一段时间内传染病的传播程度，为有关部门提供预防这种传染病的方案，现将有关建议提供如下： 开展我国新发传染病现状调查；建立国家病原微生物和新发传染病实验室；建立以实验室为基础的传染病监测系统；加强新发传染病的诊断技术研究；开展新发传染病的发生机理和预警技术研究；积极引进在我国还没有出现、但在临近国家已经出现的新发传染病的病原学诊断技术和方法；开展有关新发传染病的群众性宣传教育工作，提高公众对新发传染病的认识及如何预防的措施和方法；建立新发传染病的专业疾病控制队伍；开展新发传染病的培训工作，使基层卫生防疫人员能够了解它的临床症状、流行病学特点、病原学特点、治疗和控制措施，控制传染病源，加大隔离力度，尽早安排患者入院治疗。建立专门的医疗卫生体系，尽早提出解决治疗方案强化确诊病人和疑似病人的隔离面对突发性流行病毒的袭击，应当以预防为主，早发现、早治疗、早疏散、早隔离。【4】 参考文献： 【1】赵鸿丽，病毒扩散与传播的控制模型【ER/OL】http://www,docin.com/p-417242.html 2010-6-15 【2】无锡商业技术学院，SARS疫情分析及北京疫情走势预测，http://emuch.net/journal/article.php?id=CJFDTotal-CQZB200501060，2010-6-15 【3】任玉杰，数值分析及其MATLAB实现，北京：高等教育出版社出版社，2007. 【4】本文来自CSDN博客， http://blog.csdn.net/chengrenren/archive/2009/09/04/4519838.aspx，2010-6-15 附录： 附录A function y=ill(t,x) y=[x(2)*(1-0.8*exp(-t)),15*x(1)-x(2)*(1-0.8*exp(-t))]'; ts=0:2; >> x0=[50,280]; >> [t,x]=ode45('ill',ts,x0);[t,x] plot(t,x(:,1),t,x(:,2)) ans = 1.0e+004 * 0 0.0050 0.0280 0.0001 0.0661 0.2870 0.0002 1.5963 5.8621 附录B %说明，在编程的时候为了便于观察，将变量作以下替换： Y（1）=I Y（2）=R Y（3）=T Y（4）=D MATLAB主程序： CT=[50;58621;0;15963];h=0.1; [k,X,Y,wucha,P]=RK4z(@dydx,0,10,CT,h) plot(X,Y(:,1)) xlabel('轴\it X');ylabel('轴\it Y') legend('病人人数') m文件一： function [k,X,Y,wucha,P]=RK4z(dydx,a,b,CT,h) n=fix((b-a)/h); X=zeros(n+1,1); Y=zeros(n+1,length(CT)); X=a:h:b; Y(1,:)= CT'; for k=1:n k1=feval(dydx,X(k),Y(k,:)); x2=X(k)+h/2;y2=Y(k,:)'+k1*h/2; k2=feval(dydx,x2,y2); k3=feval(dydx,x2,Y(k,:)'+k2*h/2); k4=feval(dydx, X(k)+h,Y(k,:)'+k3*h); Y(k+1,:)=Y(k,:)+h*(k1'+2*k2'+2*k3'+k4')/6;k=k+1; end for k=2:n+1 wucha(k)=norm(Y(k)-Y(k-1)); k=k+1; end X=X(1:n+1);Y=Y(1:n+1,:);k=1:n+1;wucha=wucha(1:k,:); P=[k',X',Y,wucha']; m文件二： function dY=dydx(X,Y) dY(1)=Y(2)*(1-0.8*exp(-X))+Y(4)*(1-0.8*exp(-X))-(1-0.95*exp(-X))*Y(1); dY(2)=Y(4)*(15*exp(-0.75*X))*0.75-Y(2)*(1-0.8*exp(-X))-Y(2)*0.8*exp(-X); dY(3)=(1-0.95*exp(-X))*Y(1); dY(4)=Y(4)*(15*exp(-0.75*X))*0.25-Y(4)*(1-0.8*exp(-X)); dY=[dY(1);dY(2);dY(3);dY(4)]; 附录C ts=0:2; >> x0=[50,280]; >> [t,x]=ode45('ill',ts,x0);[t,x] plot(t,x(:,1),t,x(:,2)) 附录D ts=0:1; >> x0=[50,280]; >> [t,x]=ode45('ill',ts,x0);[t,x] plot(t,x(:,1),t,x(:,2))