埃博拉病毒的传播

西安工业大学数学建模竞赛承 诺 书我们仔细阅读了中国大学生数学建模竞赛的竞赛规则.我们完全明白，在竞赛开始后参赛队员不能以任何方式（包括电话、电子邮件、网上咨询等）与队外的任何人（包括指导教师）研究、讨论与赛题有关的问题。我们知道，抄袭别人的成果是违反竞赛规则的, 如果引用别人的成果或其他公开的资料（包括网上查到的资料），必须按照规定的参考文献的表述方式在正文引用处和参考文献中明确列出。我们郑重承诺，严格遵守竞赛规则，以保证竞赛的公正、公平性。如有违反竞赛规则的行为，我们将受到严肃处理。我们参赛选择的题号是（从A/B/C中选择一项填写）： B 所属学校（请填写完整的全名）： 西安工业大学 参赛队员 (打印并签名) ：1. 陈文兴 2. 闫丽萍 3. 魏栩 指导教师或指导教师组负责人 (打印并签名)： 日期： 2015 年 8 月 1 日埃博拉病毒传播及控制分析摘要埃博拉病毒是能引起人类和灵长类动物产生埃博拉出血热的烈性传染病病毒，有很高的死亡率。本文根据研究人员统计所给出的前四十周人类和猩猩的发病数量和死亡数量等信息，对该病毒的传播、预测与控制进行研究并建立模型，并分析了隔离措施的严格执行和药物治疗效果的提高等措施对控制疫情的作用。 针对问题一，在了解埃博拉病毒的传播情况后，根据猩猩的发病情况建立了马尔萨斯模型：。在此模型中，较好地描述病毒在“虚拟猩猩种群”中的传播情况；根据“虚拟猩猩种群”中的数据，用matlab拟合出不同状态下猩猩数量的变化曲线，并以发病状态为例建立灰色预测模型，从而较准确的预测出接下来第80、120、200周的猩猩发病状态的数据。针对问题二，为描述埃博拉病毒在“虚拟种群“中的相互传播规律及人和猩猩的疫情发展状况，建立SEIR模型模型求解时，通过对模型的推导，我们发现不能给出每个函数的解析解，因此考虑利用matlab中的ode45函数进行求解。得出了患者数量随时间的变化规律。同样利用灰色预测模型预测出“虚拟人类种群”在第80、120、200周的相关数据。针对问题三，在问题已建立的模型之上画图分析两个因素：通过某种特效药改变治愈率到80%，控制患者和健康人群的接触即控制隔离强度。对埃博拉病毒传播的影响，并通过图中控制后的患者人数，利用模型二中的关系表达式，计算出，45，50，55周的潜伏期人数，治愈人数，死亡人数。针对问题四，对问题三进一步讨论改变隔离强度和治愈率对病毒传播的影响，分别用matlab作出患者人数随时间的变化曲线，对比分析，可得出：当降低患者与健康者的接触率和使用特殊药物提高治愈率时，随着时间的延迟，患者人数急剧下降。所以实际生活中，改变通过这两个指标可以有效的控制病情的传播。关键词：马尔萨斯模型 SEIR模型 灰色预测 隔离强度一、问题重述埃博拉病毒是能引起人类和灵长类动物产生埃博拉出血热的烈性传染病病毒，埃博拉病毒感染者有很高的死亡率（在50%至90%之间），主要通过接触而非空气传播。迄今为止，已有多次疫情爆发的记录，最近的一次在2014年，截至2014年9月25日，此次在西非爆发的埃博拉疫情已经导致逾3000人死亡，另有6500被确诊为感染。本文假设某地区有20万居民和3000只猩猩。研究人员统计了前40周人类和猩猩的发病数量和死亡数量等信息（见附件一、附件二），根据相关信息研究回答以下问题：1、根据猩猩的发病数量和死亡数量，建立病毒传播模型，动态描述病毒在“虚拟猩猩种群”中的传播，并预测猩猩接下来的疫情变化，并给出“虚拟猩猩种群”在第80周、第120周、第200周的相关数据；2、建立“虚拟种群”相互感染的疾病传播模型，综合描述人和猩猩疫情的发展，并预测接下来疫情在这两个群体中的发展情况，并给出 “虚拟人类种群”在第80周、第120周、第200周的相关数据；3、由于41周外界专家的介入，及严格控制了人类与猩猩的接触，且通过某种特效药物将隔离治疗人群的治愈率提高到了80%。预测接下来疫情在“虚拟人类种群”的发展情况，并对比第2问的预测结果说明其作用和影响，给出“虚拟人类种群”在第45周、第50周、第55周的相关数据；4、依据前述数学模型，分析各种疫情控制措施的严格执行和药物（包括防疫药物、检疫药物和治疗药物等）效果的提高等措施对控制疫情的作用。2、 问题分析 根据题意，这是一个传染性病毒随时间蔓延的过程，需要研究传染病在传播过程中各类种群的数量变化，特别是通过研究患者和疑似患者、患者和治愈者的数量变化，预测传染病的传染的高峰期和持续时间长度，从而我们可以认识到相应的疫情控制措施对控制传染病传播所能达到的效果。针对问题一，根据附件一“虚拟猩猩种群”的数据，初步观察到发病状态、累计自愈及累计死亡的猩猩数量，并对数据做定量分析得到截至每周累计发病的数量，利用matlab编程得出病毒传播速度的散点图，针对病毒的传播过程，首先，我们用表示t时刻的猩猩发病个数，用表示每天每个猩猩有效接触的个数，考虑到时刻发病个数的增加，建立微分方程 ，通过马尔萨斯模型求解得：。利用灰色预测，预测出后期猩猩每个状态的数量，并以发病状态为例建立灰色预测模型，接着运用matlab编程假设 ，再对其作一次性累加生成运算得到新的生成数列，紧接着对作紧邻均值生成得出数据阵和数据向量，再对参数列进行最小二乘估计最后建立出了灰色模型（GM(1,1)模型）。我们又经过对GM(1,1)模型的残差检验，最终得出了预测结果。针对问题二，我们把人群分为五类：健康者（易感染者）、确诊患者、疑似患者、潜伏期感染者、死亡者和治愈者，采用SIER模型，并将死亡者和治愈者都归于系统移出者统称为恢复人群，关系如下图：S(t)：健康者E(t)：疑似患者Q(t)：潜伏期感染者I(t)：确诊患者S(t)：健康者R(t)：恢复人群12在此基础上，找出单位周内这五类群体数量的变化来建立微分方程，得出模型。再利用matlab编程画出图形，改变其隔离强度药物治愈率后重新作图进行比较，对结果进行分析，并利用SEIR模型对埃博拉病毒的传播规律进行定性分析描述，对未来种群数量变化用灰色预测模型进行定量预测，并分析各种疫情控制措施对控制疫情的作用。针对问题三，在问题已建立的模型之上画图分析两个因素：通过某种特效药改变治愈率，控制患者和健康人群的接触即控制隔离强度。对埃博拉病毒传播的影响。针对问题四，对问题三进一步讨论改变隔离强度和治愈率对病毒传播的影响，分别用matlab作出患者人数随时间的变化曲线，对比分析，给出有效控制疫情传播的建议。三、符号说明符号解释说明S(t)t时刻正常者（易受感染）数量E(t)t时刻疑似患者的数量Q(t)t时刻处于潜伏期的数量I(t)t时刻确诊患者的数量T(t)t时刻退出传染系统的数量（包括治愈者和死亡者）1 潜伏期的人数中转化为确诊患病的数量占潜伏期数量的比例2 每日退出传染系统的数量比例a3 确诊患者的治愈时间患者的人均周接触个体数量因接触被感染的概率潜伏期内的患者被隔离的强度a患者被治愈的概率四、模型假设1、 假设单位时间内感染病毒的数量与现有的感染者成比例；2、 假设单位时间内治愈数量与现有感染者成比例；3、 假设单位时间内死亡数量与现有的感染者成比例；4、 假设患者治愈恢复后不会再被感染同种病毒，有很强的免疫能力，即被移除出此传染系统；5、假设正常者被传染后，进入一段时间的潜伏期，处于潜伏期的群体不会表现症状，不可传染健康者，不具有传染性；6、假设患者入院即表示患者被隔离治疗，被视为无法跟别人接触，故不会传染健康者；7、假设实际治愈周期过后，如果患者没有治愈，则认为患者死亡，即实际治愈周期过后，患者都被移出此感染系统；8、假设考察地区内疾病传播期间忽略个体的出生，死亡，流动等种群动力因素对总种群数量的影响。即：总种群数量不变，记为N=20万；9、假设人能以一定的概率接触到猩猩，当接触有传播能力的猩猩后有一定的概率感染病毒，而人发病后与猩猩的接触可以忽略。10、假设可以及时发现疑似患者并隔离治疗，并且剩下一部分未被隔离的感染者变成患者后感染正常人。11、假设人类发病后与猩猩的接触可以忽略，即人类不作为猩猩的感染源。五、模型的建立与求解5.1问题一模型建立与求解据问题一，由附件“虚拟猩猩种群”中的数据,利用Matlab作累计发病个数-时间的散点图如下：（程序参见附件1）从图可看出：埃博拉病毒传播的速度在前40周始终呈上升的趋势，但上升的斜率有减小的趋势。5.1.1马尔萨斯模型的建立与求解“虚拟猩猩种群”的埃博拉病毒传播预测模型类似于人口增长的预测模型，故首先采用马尔萨斯模型（Malthusian 模型）进行建模。设时刻的确诊患者数量是连续、可微函数，并且每天每个患者有效接触（足以使人致病的接触）的人数为常数，考察到+病人人数的增加，则有 再设时有个病人，即得微分方程 解之可得： 其中， 为常数。根据“猩猩种群”疫情数据中的确诊患者的数据散点图(图1)，考虑利用马尔萨斯模型来预测埃博拉病毒的传播情况。用matlab求得，。即得马尔萨斯模型如下（程序参见附件2）： 图：结果表明，随着的增加，猩猩感染的个数持续增长。马尔萨斯拟合及预测图线与猩猩在前40周发病情况图线拟合程度较为符合。由图分析知：马尔萨斯模型是关于人口或种群增长的模型，它发现人口或种群成指数增长。即在该模型中可引意为，猩猩感染病毒的个数随着时间的增长呈指数增长变化。但现实生活中，由于猩猩数量的有限，所以该模型在预测短期的数据上是可以的，但对未来长期预测情况跟实际显然是不太相符合的，因此暂不考虑用该模型进行数据预测。因此对问题一中数据的预测采用灰色预测模型。5.1.2 利用灰色系统预测结果 我们建立的马尔萨斯较为合理的模型仅仅只是对埃博拉病毒的传播规律进了动态描述，而不能定量的进行数据预测。因此，我们先将附件中“虚拟猩猩种群”的数据用excel进行了处理：分别计算出附件中前40周猩猩的潜伏群体、处于发病状态、累计自愈、累计因病死亡的数量，然后在matlab中做出上述四个状态下猩猩数量关于周数的曲线拟合。（程序参见附件3）1、 潜伏群体 2、 处于病状态 3、 累计自愈 4、 累计因病死亡 为了预测接下来80周，120周，200周的情况，我们分别建立灰色预测模型来预测四种状态的情况，由于模型一致，本文以猩猩发病状态为例建立模型，首先取均值将前40周分为八组（具体情况见附件excel表格）。得出如下表（单位：只）：时间1-56-1011-1516-2021-2526-3031-3536-40猩猩发病数量50.254.451.248.445.442.439.236假设.5.1.2.1 GM(1,1)模型的建立为了使其成为有规律的时间序列数据，对其作一次累加生成运算，即令从而得到新的生成数列.对做紧邻均值生成. 则数据阵和数据向量为 对参数列进行最小二乘估计，可得 （其中，为发展系数，反映的发展趋势；为灰色作用量，反映数据间的变化关系. ）从而可得出GM(1,1)模型： 其中，为时间响应函数形式。5.1.2.2 GM(1,1)模型的残差检验残差大小检验，即对模型值和实际值的残差进行逐点检验. （1）根据预测公式，计算，得 （2）累减生成序列， 原始序列： （3）计算绝对残差和相对残差序列 绝对残差序列： 相对残差：GM(1,1)模型的残差检验结果：相对残差不超过0.05%，精确度高。5.1.2.3GM(1,1)模型求解利用Matlab进行预测，分别取预测的个数为8,16,32（程序参见附件4），得到实际值与预测值如下表（单位：只）：时间1-56-1011-1516-2021-2526-3031-3536-40 实际值50.254.451.248.445.442.439.236预测值50.254.8451.2947.9844.8741.9739.2336.71 时间75-80115-120195-200实际值预测值5.593.060.92本文用75-80,115-120,195-200均值近似代替第80,120,200周的情况，问题一的结果预测见下表：潜伏群体处于发病状态累计自愈累计因病死亡第80周622286600第120周313314660第200周143146605.2 问题二模型的建立与求解5.2.1 SEIR模型建立与求解针对问题二，我们充分利用附件一、二，并注意到埃博拉病毒传播过程中每一个群体都处于动态的变化中。对S来说，一部分未被隔离的潜伏期感染者能感染正常者，使其成为潜伏期感染者流出S；对于E来说，流入者包括一部分潜伏期的感染者和一部分正常者，流出者包括一部分没有被感染的正常者和隔离后被确诊患者；对于I来说，它既有从包括隔离和未被隔离的H中确诊的流入者，也有已经治愈的流出者；对于R来说，它只有从I中治愈转化而来的流入者。以上过程在传染的每一时刻都是相同的。在某一时刻对S、E、I、R取其对时间的微分，假设患者可以随意接触和感染正常人。分析阶段时间内，疫情的发展与变化，这样既可建立传染病传播模型的微分方程组如下： （1）正常人-疑似患者：病人尚未被隔离，所以疫情发展比较迅速，此时病人人均每天接触个正常人，假设t时刻病人人数为，则新增疑似患者人数为，。 （2）疑似患者-潜伏期： 疑似患者中包括病毒携带者和非病毒携带者，病毒携带者会进入潜伏期，而非病毒携带者最终还是正常人。 设疑似患者中病毒携带者占疑似患者的比例为，假设时刻疑似患者人数为，潜伏期患者人数为，则，故新增潜伏期人数为。 （3）潜伏期-确诊患者：因为每日潜伏期病人变为确诊患者的数量呈指数增长，用表示这一特性。那么新增确诊患者人数为，现在要确定，如果潜伏期天数为到，假设其变化到了一个稳定阶段，那么随着天数的增加潜伏期的病人越来越多，其概率分布呈指数稳步增长，则每天有概率的人变为埃博拉病毒患者，即。所以新增患者人数：。 （4）确诊患者-治愈、死亡：设T为退出系统人数（治愈者和死亡者），如果治愈天数设为，那么天后病人要么死亡要么被治愈，而被治愈的人产生抗体，不再会被传染，所以被治愈的人和死亡的人都算作退出系统的人。设系统退出率为，则有退出人数。的求解方法与相同，即随着天数的增加退出传染系统的人数也越来越多，则。故新退出传染系统的人数。 根据上述的式子可进一步得出：所以得出以下：5.2.2 SEIR模型求解通过对模型的推导，我们发现不能给出每个函数的解析解，因此考虑利用matlab中的ode45函数进行求解。首先，对传染病模型进行标准化，再带入参数，并由此建立微分方程组函数M文件，随后用ode函数对该文件进行调用，即可得到微分方程组的解向量，然后利用plot函数画出此解向量即可得到虚拟种群患者随周数变化的曲线图如下（程序详见附件5） 由上图分析可知:人类和猩猩患者的数量持续增多，由于人与猩猩数量差别较大，对20万人类来说，随着人类和猩猩发病者增多，感染源群体增多，易感人群发病机率增大，则人类发病数量增加较快；考虑到人发病之后与猩猩的接触可以忽略，而对于3000只猩猩来说，死亡和自愈之后立即退出传染系统，则在总量一定的情况下随着猩猩的传染源减少，猩猩的发病者数量相对于人类发病者缓慢增加。SEIR模型综合描述了40周人类与猩猩相互传染的疫情发展状况，前20周“虚拟种群系统”患者数量增加缓慢，后期感染源增多，在缺乏特效药物及隔离强度较低的情况下，患者总量增加较快，符合实际情况，可见模型的合理性。5.2.3 利用灰色系统进行数据预测（见附件Excel中的灰色预测数据及结果表格） 我们先将附件中“虚拟患者种群”的数据用excel进行了处理：分别计算出附件中前40周人类的潜伏群体、处于发病状态、隔离治疗、累计自愈、累计因病死亡的数量，然后在matlab中做出上述五个状态下人类数量关于周数的曲线拟合。（程序参见附件6）1、潜伏群体 2、处于病状态3、隔离治疗 4、累计自愈 5、累计因病死亡 利用灰色预测模型得到问题二的预测结果如下：（单位：人）潜伏人群处于发病状态隔离治疗累计治愈累计因病死亡第80周54362413223670第120周45191818455531第200周311911226595375.3 问题三模型的建立与求解5.3.1 模型建立与求解 考虑到人类与猩猩的接触率和药物治愈率两个人为控制因素，在问题二中所建立的SEIR模型的基础上进行了完善。通过降低疑似患者与易感人群的接触，即虚拟种群中患者人数I(t)与易感人群S(t)的接触率r变小，增加隔离强度P，同时依题将治愈率a提高到80%（）作图如下： 由图分析知：在未使用药物之前，40周的患者数量已达到2623人，在41周时，专家介入，有效控制易感人群（健康者）与疑似患者的接触率，即对疑似患者进行部分隔离，使得新进入潜伏期的人数在减少。因此，由于时间的延迟，患者人数短期内迅速增长，并达到峰值5436人，其后由于使用某种特效药将隔离治疗人群的治愈率提高了80%，大量的患者被治愈，且受感染的人数越来越少，导致患者人数急剧下降，最后平缓的趋于0，在80周将基本没有患者。和问题二中分别预测的第80、120、200周患者数量36、19、19人对比知，这种特效药物的使用和对患者有效隔离对埃博拉病毒的传播可以进行有效地控制。所以，在实际对埃博拉病毒的控制中，对于特效药物的研发与使用，以及加强隔离防护措施是非常有必要的。5.3.2 利用改进后的SIER模型进行数据预测1、处于发病状态人数根据上图，可以得到在45、50、55周处于发病状态的人数，即对应周数的坐标I（45）=5418, I（50）=4295, I（55）=2003。2、退出系统人数（治愈者和死亡者）在问题二所建立的（4）模型进行改进：退出系统人数，其中a3=2，P=0.05,利用1中所预测的I(t)数据，即可得到T（45）、T（50）、T（55）。在此基础上，由于人类患者的治愈率a提高到a=%80，得到：累计治愈者人数=aT（t）；累计死亡者人数=(1-a)T（t）。3、潜伏期人数因在疾病的早起阶段埃博拉病毒可能不具有：“高度的”传染性，所以对问题而中的（2）模型进行改进：取定=0.25，=43，p=0.05,潜伏期人数Q（t）=I(t),利用1 中I(t)的数据，即可得出Q（t）。4、隔离治疗人数在此问题中我们取定隔离强度P=0.05，则得到:隔离治疗人数=PQ（t）,利用3中Q（t）的预测数据，即可得到45、50、55周的隔离治疗人数。综上，所得预测结果如下表：潜伏人群处于发病状态隔离治疗累计治愈累计因病死亡第45周6115541830540041414第50周3790 429518932321063第55周1376 20036815284755.4 问题四的模型建立与求解 本问只讨论两个因素对控制埃博拉病毒传播的作用 1.易感人群与疑似患者的隔离控制即改变隔离强度 2.通过防疫.检疫.治疗功能的药物使用改变治愈率 研究结果如下图显示（程序见附件7）：由上图分析可知：四条曲线比较可知，当隔离强度不同时，对患者人数最高峰出现的时间和病毒传播的持续时间（即患者全部痊愈没有再出现患者）有极大的影响。在隔离强度较小时，患者人数的最高峰出现时间靠后，传染病持续的传播时间较长；在隔离强度较大时，患者人数能较快的出现最高峰再