埃博拉病毒的传播

西安工业大学数学建模竞赛西安工业大学数学建模竞赛 承承 诺诺 书书 我们仔细阅读了中国大学生数学建模竞赛的竞赛规则. 我们完全明白，在竞赛开始后参赛队员不能以任何方式（包括电话、电子 邮件、网上咨询等）与队外的任何人（包括指导教师）研究、讨论与赛题有关 的问题。 我们知道，抄袭别人的成果是违反竞赛规则的, 如果引用别人的成果或其 他公开的资料（包括网上查到的资料） ，必须按照规定的参考文献的表述方式在 正文引用处和参考文献中明确列出。 我们郑重承诺，严格遵守竞赛规则，以保证竞赛的公正、公平性。如有违 反竞赛规则的行为，我们将受到严肃处理。 我们参赛选择的题号是（从 A/B/C 中选择一项填写）： B 所属学校（请填写完整的全名）： 西安工业大学 参赛队员 (打印并签名) ：1. 陈文兴 2. 闫丽萍 3. 魏栩 指导教师或指导教师组负责人 (打印并签名)： 日期： 2015 年 8 月 1 日 埃博拉病毒传播及控制分析 摘要 埃博拉病毒是能引起人类和灵长类动物产生埃博拉出血热的烈性传染病病 毒，有很高的死亡率。本文根据研究人员统计所给出的前四十周人类和猩猩的 发病数量和死亡数量等信息，对该病毒的传播、预测与控制进行研究并建立模 型，并分析了隔离措施的严格执行和药物治疗效果的提高等措施对控制疫情的 作用。 针对问题一，在了解埃博拉病毒的传播情况后，根据猩猩的发病情况建立 了马尔萨斯模型：。在此模型中，较好地描述病毒在 t etx 0270 . 0 97.154 “虚拟猩猩种群”中的传播情况；根据“虚拟猩猩种群”中的数据，用 matlab 拟合出不同状态下猩猩数量的变化曲线，并以发病状态为例建立灰色预测模型 ，从而较准确的预测出接下来第 897947) 1 (1 0539.600.0669 0124 . 0 )0( 1 1 e a b e a b xkx x dt dx a 80、120、200 周的猩猩发病状态的数据。 针对问题二，为描述埃博拉病毒在“虚拟种群“中的相互传播规律及人和 猩猩的疫情发展状况，建立 SEIR 模型 )(111 )(111)(111 )(111)( 3 312 12 tIeadtdT tIeatQeaadtdI tQeaatIrdtdQ t tt t 模型求解时，通过对模型的推导，我们发现不能给出每个函数的解析解，因此 考虑利用 matlab 中的 ode45 函数进行求解。得出了患者数量随时间的变化规律。 同样利用灰色预测模型预测出“虚拟人类种群”在第 80、120、200 周的相关数 据。 针对问题三，在问题已建立的模型之上画图分析两个因素：通过某种特效 药改变治愈率到 80%，控制患者和健康人群的接触即控制隔离强度。对埃博拉 病毒传播的影响，并通过图中控制后的患者人数，利用模型二中的关系表达式， 计算出，45，50，55 周的潜伏期人数，治愈人数，死亡人数。 针对问题四，对问题三进一步讨论改变隔离强度和治愈率对病毒传播的影 响，分别用 matlab 作出患者人数随时间的变化曲线，对比分析，可得出：当降 低患者与健康者的接触率和使用特殊药物提高治愈率时，随着时间的延迟，患 者人数急剧下降。所以实际生活中，改变通过这两个指标可以有效的控制病情 的传播。 1 关键词关键词：马尔萨斯模型 SEIR 模型 灰色预测 隔离强度 一、问题重述 埃博拉病毒是能引起人类和灵长类动物产生埃博拉出血热的烈性传染病病 毒，埃博拉病毒感染者有很高的死亡率（在 50%至 90%之间） ，主要通过接触而 非空气传播。迄今为止，已有多次疫情爆发的记录，最近的一次在 2014 年，截 至 2014 年 9 月 25 日，此次在西非爆发的埃博拉疫情已经导致逾 3000 人死亡， 另有 6500 被确诊为感染。 本文假设某地区有 20 万居民和 3000 只猩猩。研究人员统计了前 40 周人类 和猩猩的发病数量和死亡数量等信息（见附件一、附件二） ，根据相关信息研究 回答以下问题： 1 1、根据猩猩的发病数量和死亡数量，建立病毒传播模型，动态描述病毒在 “虚拟猩猩种群”中的传播，并预测猩猩接下来的疫情变化，并给出“虚拟猩 猩种群”在第 80 周、第 120 周、第 200 周的相关数据； 2 2、建立“虚拟种群”相互感染的疾病传播模型，综合描述人和猩猩疫情的 发展，并预测接下来疫情在这两个群体中的发展情况，并给出 “虚拟人类种群” 在第 80 周、第 120 周、第 200 周的相关数据； 3 3、由于 41 周外界专家的介入，及严格控制了人类与猩猩的接触，且通过 某种特效药物将隔离治疗人群的治愈率提高到了 80%。预测接下来疫情在“虚 拟人类种群”的发展情况，并对比第 2 问的预测结果说明其作用和影响，给出 “虚拟人类种群”在第 45 周、第 50 周、第 55 周的相关数据； 4 4、依据前述数学模型，分析各种疫情控制措施的严格执行和药物（包括防 疫药物、检疫药物和治疗药物等）效果的提高等措施对控制疫情的作用。 2 2、问题分析 根据题意，这是一个传染性病毒随时间蔓延的过程，需要研究传染病在传 播过程中各类种群的数量变化，特别是通过研究患者和疑似患者、患者和治愈 者的数量变化，预测传染病的传染的高峰期和持续时间长度，从而我们可以认 识到相应的疫情控制措施对控制传染病传播所能达到的效果。 针对问题一，根据附件一“虚拟猩猩种群”的数据，初步观察到发病状态、 累计自愈及累计死亡的猩猩数量，并对数据做定量分析得到截至每周累计发病 的数量，利用 matlab 编程得出病毒传播速度的散点图，针对病毒的传播过程， 首先，我们用表示 t 时刻的猩猩发病个数，用表示每天每个猩猩有效接tx 触的个数，考虑 到时刻发病个数的增加，建立微分方程 ，tttx dt dx ，通过马尔萨斯模型求解得：。利用灰色预测，预测出后期 0 0xx t extx 0 猩猩每个状态的数量，并以发病状态为例建立灰色预测模型，接着运用 matlab 3 编程假设 ，再对其作一次性累加 3642.4,39.2,48.4,45.4,54.4,51.2,50.2 0 ，kX 生成运算得到新的生成数列 ，紧接着对作紧邻 67.292,331.2,3.2,249.6,2,155.8,20450.2,104.6 1 kX kX 1 均值生成得出数据阵和数据向量，再对参数列进行最小二乘估计B n Y T ba, 最后建立出了灰色模型（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 确诊患者的治愈时间 r 患者的人均周接触个体数量 因接触被感染的概率 p 潜伏期内的患者被隔离的强度 a患者被治愈的概率 四、模型假设 1、假设单位时间内感染病毒的数量与现有的感染者成比例； 2、假设单位时间内治愈数量与现有感染者成比例； 3、假设单位时间内死亡数量与现有的感染者成比例； 4、假设患者治愈恢复后不会再被感染同种病毒，有很强的免疫能力，即被移除 出此传染系统； 5、假设正常者被传染后，进入一段时间的潜伏期，处于潜伏期的群体不会表现 症状，不可传染健康者，不具有传染性； 6、假设患者入院即表示患者被隔离治疗，被视为无法跟别人接触，故不会传染 健康者； 7、假设实际治愈周期过后，如果患者没有治愈，则认为患者死亡，即实际治愈 周期过后，患者都被移出此感染系统； 8、假设考察地区内疾病传播期间忽略个体的出生，死亡，流动等种群动力因素 对总种群数量的影响。即：总种群数量不变，记为 N=20万； 9、假设人能以一定的概率接触到猩猩，当接触有传播能力的猩猩后有一定的概 率感染病毒，而人发病后与猩猩的接触可以忽略。 10、假设可以及时发现疑似患者并隔离治疗，并且剩下一部分未被隔离的感染 者变成患者后感染正常人。 11、假设人类发病后与猩猩的接触可以忽略，即人类不作为猩猩的感染源。 五、模型的建立与求解 5.1 问题一模型建立与求解 据问题一，由附件“虚拟猩猩种群”中的数据,利用 Matlab 作累计发病个 数-时间的散点图如下：（程序参见附件 1） 5 从图可看出：埃博拉病毒传播的速度在前 40 周始终呈上升的趋势，但上升 的斜率有减小的趋势。 5.1.1 马尔萨斯模型的建立与求解 “虚拟猩猩种群”的埃博拉病毒传播预测模型类似于人口增长的预测模型， 故首先采用马尔萨斯模型（Malthusian 模型）进行建模。设时刻 的确诊患者t 数量是连续、可微函数，并且每天每个患者有效接触（足以使人致病的接tI 触）的人数为常数，考察 到 +病人人数的增加，则有ttt ttItIttI 再设时有个病人，即得微分方程0t 0 I 0 0II I dt dI 解之可得： t eItI 0 其中， 为常数。, 0 I t etI 0270. 0 97.154 根据“猩猩种群”疫情数据中的确诊患者的数据散点图(图 1)，考虑利用 马尔萨斯模型来预测埃博拉病毒的传播情况。用 matlab 求得 t extx 0 ，。即得马尔萨斯模型如下（程序参见附件 2）：97.154 0 I0270 . 0 t etI 0270. 0 97.154 图： 结果表明，随着 的增加，猩猩感染的个数持续增长。马尔萨斯拟合及预测ttI 图线与猩猩在前 40 周发病情况图线拟合程度较为符合。 由图分析知：马尔萨斯模型是关于人口或种群增长的模型，它发现人口或种群 成指数增长。即在该模型中可引意为，猩猩感染病毒的个数随着时间的增长呈 指数增长变化。但现实生活中，由于猩猩数量的有限，所以该模型在预测短期 的数据上是可以的，但对未来长期预测情况跟实际显然是不太相符合的，因此 暂不考虑用该模型进行数据预测。因此对问题一中数据的预测采用灰色预测模 型。 5.1.2 利用灰色系统预测结果 我们建立的马尔萨斯较为合理的模型仅仅只是对埃博拉病毒的传播规律进 了动态描述，而不能定量的进行数据预测。因此，我们先将附件中“虚拟猩猩 种群”的数据用 excel 进行了处理：分别计算出附件中前 40 周猩猩的潜伏群体、 处于发病状态、累计自愈、累计因病死亡的数量，然后在 matlab 中做出上述四 个状态下猩猩数量关于周数的曲线拟合。 （程序参见附件 3） 1、潜伏群体 2、处于病状态 3、累计自愈 4、累计因病死亡 7 为了预测接下来 80 周，120 周，200 周的情况，我们分别建立灰色预测模 型来预测四种状态的情况，由于模型一致，本文以猩猩发病状态为例建立模型， 首先取均值将前 40 周分为八组（具体情况见附件 excel 表格） 。 得出如下表（单位：只）： 假设. 3642.4,39.2,48.4,45.4,54.4,51.2,50.2 0 ，kX 5.1.2.1 GM(1,1)模型的建立 为了使其成为有规律的时间序列数据，对其作一次累加生成运算，即令 )72 , 1,()()( 1 )0()1( ktkXtX t n 从而得到新的生成数列.对 67.292,331.2,3.2,249.6,2,155.8,20450.2,104.6 1 kX 做紧邻均值生成. 则数据阵和数据向量为 kX 1 B n Y 时间 1-56-1011-1516-2021-2526-3031-3536-40 猩猩发病 数量 50.254.451.248.445.442.439.236 1 349.2- 1 311.6- 1 270.8 - 1 226.9- 1 180- 1 130.2- 177.4- 187 2 1 176 2 1 165 2 1 154 2 1 143 2 1 132 2 1 121 2 1 11 11 11 11 11 11 11 xx xx xx xx xx xx xx B 367.2 331.2 292 249.6 204.2 155.8 104.6 8 7 6 5 4 3 2 0 0 0 0 0 0 0 x x x x x x x Yn 对参数列进行最小二乘估计，可得 T ba, 3483. 0-2044 . 0 -0482. 0-1198 . 0 2993 . 0 4899 . 0 6920 . 0 0022. 0-0016 . 0 -0009. 0-0001 . 0 -0007 . 0 0016 . 0 0025 . 0 1 TT BBB （其中，为发展系数，反映的发展 60.0539 0.0669 1 n TT YBBB b a ax 趋势；为灰色作用量，反映数据间的变化关系. ）b 从而可得出 GM(1,1)模型： 897947) 1 (1 0539.600.0669 0124 . 0 )0( 1 1 e a b e a b xkx x dt dx a 其中，为时间响应函数形式。8970669 . 0 947) 1 (1 )0( e a b e a b xkx a 5.1.2.2 GM(1,1)模型的残差检验 残差大小检验，即对模型值和实际值的残差进行逐点检验. （1）根据预测公式，计算，得 kX 1 8, 2 , 1 7.22,331.2,362,249.6,29155.8,204.104.650.2 1 kkX， （2）累减生成序列， kX 0 82 , 1 k 84,36.71939731,39.254.8755,41.,47.9786,434,51.296350.2,54.84 0 kX 9 原始序列： 39.2,3645.4,42.4,51.2,48.4,50.2,54.4, 0 kX （3）计算绝对残差和相对残差序列 绝对残差序列： 015,0.020.0101,0.087,0.0116,.0019,0.000,0.0082,0 0 相对残差： 055%0.0038%,0.0.0238%0.0255%00179%0037 . 0 %015. 00， GM(1,1)模型的残差检验结果：相对残差不超过 0.05%，精确度高。 5.1.2.3GM(1,1)模型求解 利用 Matlab 进行预测，分别取预测的个数为 8,16,32（程序参见附件 4） ，得到 实际值与预测值如下表（单位：只）： 本文用 75-80,115-120,195-200 均值近似代替第 80,120,200 周的情况，问题一 的结果预测见下表： 潜伏群体处于发病状态累计自愈累计因病死亡 第 80 周 622286600 第 120 周 313314660 第 200 周 14314660 5.2 问题二模型的建立与求解 5.2.1 SEIR 模型建立与求解 针对问题二，我们充分利用附件一、二，并注意到埃博拉病毒传播过程中 每一个群体都处于动态的变化中。对 S 来说，一部分未被隔离的潜伏期感染者 能感染正常者，使其成为潜伏期感染者流出 S；对于 E 来说，流入者包括一部 分潜伏期的感染者和一部分正常者，流出者包括一部分没有被感染的正常者和 隔离后被确诊患者；对于 I 来说，它既有从包括隔离和未被隔离的 H 中确诊的 流入者，也有已经治愈的流出者；对于 R 来说，它只有从 I 中治愈转化而来的 流入者。以上过程在传染的每一时刻都是相同的。在某一时刻对 S、E、I、R 取 其对时间的微分，假设患者可以随意接触和感染正常人。分析阶段时间内，t 疫情的发展与变化，这样既可建立传染病传播模型的微分方程组如下： 时间 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 （1）正常人-疑似患者： 病人尚未被隔离，所以疫情发展比较迅速，此时病人人均每天接触 个正r 常人，假设 t 时刻病人人数为，则新增疑似患者人数为，)(tIE 。ttIrtrtI)()(E （2）疑似患者-潜伏期： 疑似患者中包括病毒携带者和非病毒携带者，病毒携带者会进入潜伏期， 而非病毒携带者最终还是正常人。 设疑似患者中病毒携带者占疑似患者的比例为，假设 时刻疑似患者人数t 为，潜伏期患者人数为，则，故新增潜伏期人数为)(E t)(Q t)()(tEtQ 。EQ （3）潜伏期-确诊患者： 因为每日潜伏期病人变为确诊患者的数量呈指数增长，用表示这一特性。 1 那么新增确诊患者人数为，现在要确定，如果潜伏期天数为ttQI)( 1 1 到，假设其变化到了一个稳定阶段，那么随着天数的增加潜伏期的病人越 1 a 2 a 来越多，其概率分布呈指数稳步增长，则每天有概率的人变 t eaa 12 111 为埃博拉病毒患者，即。所以新增患者人数： t eaa 121 111 。tQeaaI t 12 111 （4）确诊患者-治愈、死亡： 设 T 为退出系统人数（治愈者和死亡者） ，如果治愈天数设为，那么 3 a 天后病人要么死亡要么被治愈，而被治愈的人产生抗体，不再会被传染，所 3 a 以被治愈的人和死亡的人都算作退出系统的人。设系统退出率为，则有退出 3 a 人数。的求解方法与相同，即随着天数的增加退出传染系ttIT 2 )( 2 1 统的人数也越来越多，则。故新退出传染系统的人数 t ea 32 111 。tIeaT t 3 111 根据上述的式子可进一步得出：）（）（41 11 ttIeatTttT ttIeattQeaatIttI ttQeaattIrtQttQ EQ t tt t )(111)()( )(111)(111)()( )(111)()()( 3 312 12 所以得出以下： 5.2.2 SEIR 模型求解 通过对模型的推导，我们发现不能给出每个函数的解析解，因此考虑利用 matlab 中的 ode45 函数进行求解。 首先，对传染病模型进行标准化，再带入参数，并由此建立微分方程组函 数 M 文件，随后用 ode 函数对该文件进行调用，即可得到微分方程组的解向量， 然后利用 plot 函数画出此解向量即可得到虚拟种群患者随周数变化的曲线图如 下（程序详见附件 5） 由上图分析可知:人类和猩猩患者的数量持续增多，由于人与猩猩数量差别 较大，对 20 万人类来说，随着人类和猩猩发病者增多，感染源群体增多，易感 人群发病机率增大，则人类发病数量增加较快；考虑到人发病之后与猩猩的接 触可以忽略，而对于 3000 只猩猩来说，死亡和自愈之后立即退出传染系统，则 )(111 )(111)(111 )(111)( 3 312 12 tIeadtdT tIeatQeaadtdI tQeaatIrdtdQ t tt t 在总量一定的情况下随着猩猩的传染源减少，猩猩的发病者数量相对于人类发 病者缓慢增加。 SEIR 模型综合描述了 40 周人类与猩猩相互传染的疫情发展状况，前 20 周 “虚拟种群系统”患者数量增加缓慢，后期感染源增多，在缺乏特效药物及隔 离强度较低的情况下，患者总量增加较快，符合实际情况，可见模型的合理性。 5.2.3 利用灰色系统进行数据预测（见附件 Excel 中的灰色预测数据及结果表 格） 我们先将附件中“虚拟患者种群”的数据用 excel 进行了处理：分别计算 出附件中前 40 周人类的潜伏群体、处于发病状态、隔离治疗、累计自愈、累计 因病死亡的数量，然后在 matlab 中做出上述五个状态下人类数量关于周数的曲 线拟合。 （程序参见附件 6） 1、潜伏群体 2、处于病状态 3、隔离治疗 4、累计自愈 5、累计因病死亡 利用灰色预测模型得到问题二的预测结果如下：（单位：人） 5.3 问题三模型的建立与求解 5.3.1 模型建立与求解 考虑到人类与猩猩的接触率和药物治愈率两个人为控制因素，在问题二中所 建立的 SEIR 模型的基础上进行了完善。通过降低疑似患者与易感人群的接触， 潜伏人群处于发病状态隔离治疗累计治愈累计因病死亡 第 80 周54362413223670 第 120 周45191818455531 第 200 周31191122659537 13 即虚拟种群中患者人数 I(t)与易感人群 S(t)的接触率 r变小，增加隔离强度 P，同时依题将治愈率 a 提高到 80%（）作图如下：%801 aa 由图分析知：在未使用药物之前，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)数据，)(111t p 3 tIeaT t ）（ 即可得到 T（45） 、T（50） 、T（55） 。在此基础上，由于人类患者的治愈率 a 提 高到 a=%80，得到：累计治愈者人数=aT（t） ；累计死亡者人数=(1-a)T（t） 。 3、潜伏期人数 因在疾病的早起阶段埃博拉病毒可能不具有：“高度的”传染性，所以对 问题而中的（2）模型进行改进：取定=0.25，=43，p=0.05,潜伏期人数r Q（t）=I(t),利用 1 中 I(t)的数据，即可得出 Q（t） 。r t e p 4、隔离治疗人数 在此问题中我们取定隔离强度 P=0.05，则得到:隔离治疗人数=PQ（t）,利 用 3 中 Q（t）的预测数据，即可得到 45、50、55 周的隔离治疗人数。 综上，所得预测结果如下表： 5.4 问题四的模型建立与求解 本问只讨论两个因素对控制埃博拉病毒传播的作用 1.易感人群与疑似患者的隔离控制即改变隔离强度 2.通过防疫.检疫.治疗功能的药物使用改变治愈率 研究结果如下图显示（程序见附件 7）： 由上图分析可知：四条曲线比较可知，当隔离强度不同时，对患者人数最高峰