甲流H1N1传播数学模型

1、摘要甲型H1N1流感的迅速蔓延引起了世界各国的广泛关注。本文结合当前出现的新形势，对初始模型不断改进，根据患病后机体存在免疫性这一实际情况，建立模型来分析HINI的传播特点。然后，选取三个代表性国家以定性与定量分析相结合探讨H1N1的传播特点。最终，为相关部门提出有效的防控对策。首先，我们建立了针对一般的传染病传播的微分方程模型一，得出了患病人数无限增长的结论。在此模型的基础上，通过对易感染者和未感染者加以区分，将其改进为Logistic 模型，模型二可以预测出传染病高潮到来的时刻。在假设该传染病无免疫性的前提下，结合病人可以治愈的实际情况，建立模型三，在此模型中得到的接触数s =1 是一个阈

2、值，得出控制传染病的关键在于感染期内有效接触感染的健康者人数不超过病人数的结论。进而，根据H1N1患者在患病后就存在免疫性即可移出感染系统的特点，建立模型四。为了避免不能得到解析解的情况，我们运用MATLAB软件，根据龙格库塔方法求解，并根据相轨图分析传播特点。考虑到各国都采取了隔离的措施来延缓H1N1的传播，我们通过模型五来描述这种影响。由于模型五较为复杂，本文未做详细讨论。然后，结合WTO公布的数据，针对这次甲型H1N1流感的传播的特点采用了模型四，通过定量地分析发现H1N1流感疫情在全球范围内得到了明显控制但总趋势上仍在蔓延。随后选取墨西哥、美国、日本这三个疫情较为严重的国家对H1N1流

3、行进行了分析。发现美国的疫情得到了部分控制，墨西哥和日本的疫情不容乐观。最近，H1N1病毒变异后出现耐药性，这是新情况，应引起重视。最后针对目前的情况，我们认为通过增强集体免疫的意识，提高卫生水平，完善医疗结构可以有效防控疫情。关键字： 微分方程模型 龙格库塔方法 相轨图 MATLAB 目 录1.问题的重述与分析32.模型的假设33.符号说明44.模型建立5 4.1 模型演化的基本思路.5 4.2 模型一的建立.54.3 模型二的建立.64.4 模型三的建立.74.5 模型四的建立.84.6 模型五的建立.95.模型的求解及结果的分析95.1 模型四的求解95.2 相轨迹的分析.115.3 H

4、1N1在全球的传播特点分析.145.4 H1N1在墨西哥的传播特点分析.175.5 H1N1在日本的传播特点分析.205.6 H1N1在美国的传播特点分析215.7 H1N1的流行特点.226我们的建议.237.模型优缺点的分析.248.参考文献.249.附录 .241问题的重述与分析 在2009年4月下旬，世界卫生组织宣布出现一种新的甲型流感病毒，即甲型H1N1流感。甲型H1N1流感（influenza A (H1N1)）又为A（H1N1）型流感，人感染猪流感。是一种急性、传染性呼吸器官疾病。其特征为突发，咳嗽，呼吸困难，发热及迅速转归。它携带有H1N1亚型猪流感病毒毒株，包含有禽流感、猪流

5、感和人流感三种流感病毒的核糖核酸基因片断，同时拥有亚洲猪流感和非洲猪流感病毒特征。这种病毒是全新的，以前未曾有过人间传播，病毒具有传染性，很容易在人与人和国与国之间传播。流感蔓延已给世界的经济发展和人民生活带来了很大影响。由于各国进行了仔细的监测、彻底的调查和坦诚的报告，我们已对病毒的传播以及可造成的一系列病症有了一定的初步了解，为预测和控制传染病蔓延创造了重要条件。我们首先对一般的传染病的传播建立模型。分析传染病蔓延的条件和控制传染病蔓延的措施。然后结合WTO公布的数据，针对这次甲型H1N1流感的传播的特点建立数学模型，定量地分析在世界范围的传播情况。再选取墨西哥、美国、日本这三个有特点的国

6、家进行H1N1流行的分析。最终根据最新的情况分析本次H1N1流行的新特点和新情况。本次H1N1流行的爆发时刻进行预测，并针对各个不同的地区的情况提出防控建议。2.模型假设1.WTO提供的全国疫情统计真实可信。2.将H1N1所有可能的传播途径都视为与病源的直接接触。3.在疾病传播期内所考察的地区的总人数N视为常数,即认为本地区流入的人口与流入的人口数相等,时间以天为计量单位。4.设每个病人单位时间有效接触的人数(所谓“有效接触”是指病人与健康者接触时,足以使健康者受到感染而成为病人)可视为常数。5.根据资料可知未出现症状的感染H1N1的人处于潜伏期并且不具有传染性。6.假设潜伏期为一常数。7根据

7、目前的医学调查资料，对于H1N1一个康复者，他势必会更注重自己的个人卫生习惯并主动远离H1N1传染源；从社会心理学的角度来看，其身边的人会主动远离他。因此，我们可以假设一个H1N1康复者二度感染SARS的概率为0，这些人归为“退出者”。8.流入和流出的人群中的带菌者处于潜伏期。9. 被隔离的人群完全断绝与外界的接触,不再具有传染性3.符号说明i：已感染H1N1病毒并已经出现明显症状的人所占的比例s: 未感染H1N1病毒的健康人r：移出感染系统的人称为移出者l：每个病人每天有效接触的人数，有效接触是指足以使人致病。N：总人数s：表示的意义为在一个感染期内每个病人的有效接触人数，称为接触数。m：感

8、染H1N1病人每天治愈的比例Q：退出率，为H1N1患者的日死亡率和日治愈率之和。G：隔离者在人群中的比例W：未隔离者在人群中的比例：接触病源的人的发病率。：每天由可控人群和不可控人群转化为病人的日转化率。：不可控人群(在后面的分析中可得到)在发病后到被隔离前平均每天接触的人的数目。4.模型的建立 4.1模型演化的基本思路 我们从最一般的情况开始讨论，由浅入深逐步改进，共建立出五个模型逐步逼近现实的情况。首先只考虑每个病人每天有效接触人数l，建立起基本的微分方程，得出于实际不相符的情况，即病人人数无限增加。原因是若有效接触的是病人，则不能使病人数增加。针对模型一的缺点，我们区分了易感染者和未感染

9、者，我们建立了模型二，发现所得的方程为Logistic 模型，仍然是一个无限增长的模型并可以预测出传染病高潮到来时刻。模型预测出最终所有的人都会感染与实际不相符。考虑病人可以治愈，另外假设这种传染病无免疫性，我们可以改进得到模型三，在模型中得到的接触数s =1 是一个阈值，控制传染病的关键在于感染期内有效接触感染的健康者人数不超过病人数。针对H1N1这种得病后就有免疫性的疾病。在传染病有免疫即病人治愈后即移出感染系统的条件下，建立模型四。移出感染系统的人称为移出者。这个方程不能得到解析解，我们根据相轨图分析传播特点。在发现H1N1流行的过程中各国都采取了有效地措施，考虑到隔离者对模型的影响，我

10、们建立了模型五。4.2模型一的建立首先假设已感染人数 (病人) i(t)，每个病人每天有效接触(足以使人致病)人数为l。则在一个小时间段内新增的感染人数应全部是这段时间内已感染的病人传染给健康人的，由此可得方程： （1）在无限小的时候可以将方程（1）化简为微分方程： （2）设这个微分方程的初值为: ，可以得到这个微分方程的解解这个方程可得： 方程（1）得出的题存在明显的问题：当时，说明感染的病人数无限增多，传染病不可控制。这明显是不合理的结论。考虑到若有效接触的是病人，则不能使病人数增加，模型存在改进的空间。为了解决这个不合理的因素必须区分已感染者(病人)和未感染者(健康人)。4.3模型二的建

11、立针对模型一的缺点，模型二改进部分在于区分了已感染者和未感染者即健康人和病人。首先假设总人数N不变，病人和健康人的 比例分别为。每个病人每天有效接触人数为l, 且使接触的健康人致病。将l成为日治愈率。我们仍然研究时间内的情况。这段时间内新增病例仍然是已感染的病人传染的，建立模型方程如下： （3）同模型一同样，在无限小的时候可以将方程（1）化简为微分方程： （4）已感染人数的比例和未感染人数的比例之和应该为1，即 （5）联立（4）（5）这两方程可以得到方程组： （6） 对方程（6）求解可得: （7）我们发现所得的方程为Logistic 模型，做出方程的图像。1/2tmii010t 由图形知： （

12、8）当t取tm时di/dt取最大值，说明这个时候是传染病高潮到来时刻,应该引起我们的重视。从图像中还可以看出当日接触率l下降的时候，图像中的tm 的位置将会右移增大。说明控制日接触率可以延缓传染病爆发的到来。另外当 t趋近于时，i趋近于1 。说明经过无限长的时间，所以人都会得病，这显然与实际不符。产生错误的原因在于没有考虑病人可以治愈这个因素。说明这个模型仍然有改进的空间。 4.4 模型三的建立针对模型二的缺点，考虑病人可以治愈这个因素我们又建立了模型三，另外还要假设这种传染病无免疫性，也就是说，病人治愈成为健康人，健康人可再次被感染。同时总人数N不变，病人和健康人的 比例分别为，每个病人每天

13、有效接触人数为l, 且使接触的健康人致病，将l称为日接触率。病人每天治愈的比例为m称为日治愈率。 仍然在一段时间新增的病例数为在这段时间内被传染的人数减去已近治愈的人，我们建立数学方程： （9）可以将这个表示为微分方程的形式： （10）为了化简方程，我们引入了变量，其中l 为日接触率，1/m为感染期。s表示的意义为在一个感染期内每个病人的有效接触人数，称为接触数。将s代入微分方程可以将方程表示为： （11）根据方程（11）画出di/dt于i的图像： di/dti01s >11-1/s 接着我们做出it的图像，发现当s取不同值的时候it 的图像不同。如下图所示：i0ts £1i0

14、0tis >11-1/si0i0从图像中可以看出： （12） 由此可见接触数s =1 是一个阈值，只有s<1病人才可以全部治愈，传染病不再流行。所以控制传染病的关键在于感染期内有效接触感染的健康者人数不超过病人数.4.5 模型四的建立这个模型近似于经典的S-I-R模型。S(susceptible)是易感人群，I(Infected)是感染人群，R(Recovered)是已经康复的人群。感染人群有一定几率传染易感人群，使其转变成为感染人群。而感染人群也可能得到治愈成为已康复的人群。感染的强度和恢复天数的长短可以由参数控制。在模型三中我们假设了传染病无免疫性，但考虑到H1N1这种治愈后就

15、有免疫性的疾病，应该假设传染病有免疫性，免疫即病人治愈后即移出感染系统的条件下。在传染病有免疫性的条件下，我们建立模型四。移出感染系统的人称为移出者。 总人数N不变，病人、健康人和移出者的比例分别为。病人的日接触率l , 日治愈率m, 接触数 s = l / m。同样在一段时间内建立方程：病人、健康人和移出者比例之和为1： （13）在一段时间新增的病例数为在这段时间内被传染的人数减去已近治愈的人： （14）健康人减少的人数等于被已感染的人传染的人数： （15） 上述方程组用微分方程组表示如下： （16）对于这个方程组无法求出 的解析解4.6 模型五的建立因为在H1N1流行的过程中，各个国家都采

16、取了一定的措施，大部分国家都采取了隔离的制度，所以在模型四的基础到，考虑到隔离人数比例G和未隔离人数比例W，还有接触过发病后没有及时隔离治疗的人的人数 ，我们又建立了改进的模型，方程组如下： （17） （18） （19） （20） （21）以上方程组的初值分别代表各个字母代表的初值。因为这个方程组的分析过于复杂，所以这里就不对这个模型的建立多做讨论了。5. 模型的求解及结果分析 5.1 模型四的求解 在模型建立的部分中我们一共建立四个模型，其中模型五是最接近实际情况的，但由于这个模型设计的参数过多，数据难于确定，并不利于对H1N1的传播流行进行分析。所以我们采用模型四进行分析求解，试图找出H1

17、N1传播的规律的特点。模型四建立的方程组如下: (16)s(t) ， i(t)的求解极度困难，在此我们先做数值计算来预估计s(t) ， i(t)的一般变化规律。在方程（16）中设=1，=0.3，i（0）= 0.02，s（0）=0.98，我们求助于matlab中的龙格库塔方法来求出它们的数值解。程序见附件一。输出的简明计算结果列入表1。i(t) , s(t)的图形以下两个图形，is图形称为相轨线,初值i(0)=0.02,s(0)=0.98相当于图2中的P0点,随着t的增,(s,i)沿轨线自右向左运动.由表1、图1、图2可以看出,i(t)由初值增长至约t=7时达到最大值,然后减少,t,i0,s(t

18、)则单调减少,t,s0.0398. 并分析i(t),s(t)的一般变化规律.t 0 1 2 3 4 5 6 7 8i(t)0.02000.03900.07320.12850.20330.27950.33120.34440.3247s(t)0.98000.95250.90190.81690.69270.54380.39950.28390.2027 t 9 10 15 20 25 30 35 40 45i(t)0.28630.24180.07870.02230.00610.00170.00050.00010s(t)0.14930.11450.05430.04340.04080.04010.0399

19、0.03990.0398 5.2相轨线的分析我们在数值计算和图形观察的基础上，利用相轨线讨论解i（t）,s（t）的性质。 i s平面称为相平面，相轨线在相平面上的定义域（s，i）D D = （s，i）| s0，i0 ， s + i 1 （22） 在方程16消去并注意到的定义，可得， （23） 利用积分特性容易求出方程(23)的解为： （24）在定义域D内,24式表示的曲线即为相轨线,如图所示.其中箭头表示了随着时间t的增加s(t)和i(t)的变化趋向.从相轨线可以得到以下几个结论：1.不论初始条件s0,i0如何,病人消失将消失,即： 其证明如下： 首先,由(3) 而 故 存在; 由(2) 而

20、故 存在;再由(1)知存在。其次,若,对于充分大的t 有 , 这将导致，与 存在相矛盾.从图形上看,不论相轨线从P1或从P2点出发,它终将与s轴相交(t充分大).2. 被感染的健康者的比例是,在图形上是相轨线与s轴在(0,1/)内交点的横坐标。在式(24)中令i=0得到, 是方程 （25）在(0,1/)内的根.3. 1/是一个阈值，当>1/时传染病会蔓延，1/时传染病就不会蔓延若>1/,则开始有，i(t)先增加, 令=0，可得当s=1/时,i(t)达到最大值： （26） 然后s<1/时，有 ，所以i(t)减小且趋于零,s(t)则单调减小至,如图中由P1(，)出发的轨线.若 1

21、/,则恒有，i(t)单调减小至零,s(t)单调减小至,如图中由P2(s0,i0)出发的轨线 可以看出,如果仅当病人比例i(t)有一段增长的时期才认为传染病在蔓延,那么1/是一个阈值,当>1/(即>1/s0)时传染病就会蔓延.而减小传染期接触数,即提高阈值1/使得1/(即 1/),传染病就不会蔓延(健康者比例的初始值是一定的,通常可认为接近1)。 并且,即使>1/, 减小时, 增加(通过作图分析), 降低,也控制了蔓延的程度.我们注意到在=中,人们的卫生水平越高,日接触率越小;医疗水平越高,日治愈率越大,于是越小,所以提高卫生水平和医疗水平有助于控制传染病的蔓延. 从另一方面看

22、, 是传染期内一个病人传染的健康者的平均数,称为交换数,其含义是一病人被个健康者交换.所以当 即时必有 .既然交换数不超过1,病人比例i(t)绝不会增加,传染病不会蔓延。5.3 H1N1在全球的传播特点分析这个初夏，甲型H1N1流感袭来。我们知道，传染病流行是有一定规律可循的，比如，一些疾病有着固定的易传染时间段，从感染到发作的时间比较固定，传染能力、致死率等因素也可以被我们获知。这意味着，我们可以建立模型四来描述传染病传播的特点，预测传播的规模、速度等。在这个模型中，最重要的因素之一是流行病的传播能力，也就是一个患者平均可以传染几个人，这个数值叫做接触数。如果接触数小于1，那么流行病就能被控

23、制住。如果大于1，就有流行的风险。根据推算，1918年的西班牙流感期间，再生数大于2.5，所以造成了大面积的杀伤。这一次甲型H1N1流感，几个研究组通过开始阶段的观察，初步估计是这个数值大约在1.42.5之间。 根据我们搜集到数据，首先用MATLAB画出5月18日到11月15日这个时间段内H1N1在全球的发病数,如下图所示：根据得出的图形，可以看出H1N1全球蔓延趋势如下：在H1N1流行的初期，就是在开始统计的45天的时间段病例总数增长并不快，增长率近似于1。甚至在一定阶段显现出缓慢增加的形式，由于初期病例数量小造成的慢速的传染现象。一至两周后，总病例数开始急剧上升，增长趋势类似于指数型爆发增

24、长，这有可能是由于患病人数增加增大了传染几率，短时间内有大量个体被传染并且发病，同时由于该病在传染给人后会有持续近一周的潜伏期，大量前期被传染的个体发病，造成了总病例数急剧攀升。大约在一个半月后，总病例增长趋势放缓，开始平稳增长，伴随小幅振荡，这时该病毒的传染进入了平稳期，染病人数增长率稳定。但是，从已知数据和对病毒的研究分析看来，至少从统计开始的200天的时间段内，感染H1N1总人数没有平缓下来的迹象，可见该病毒传染性强，应该增强防控手段。根据我们搜集到数据，首先用MATLAB画出5月18日到11月15日这个时间段内H1N1在全球的发病数,如下图所示：全球因为感染H1N1而死亡的病例趋势图如

25、下：在对H1N1全球死亡病例趋势图的分析中，我们可以看到死亡趋势与病毒的感染发病总人数趋势有一致的变化，但整体上变化向后推迟了2周。这两周时间大致等于重症患者从发病到死亡的平均时间周期，与医学数据相统一。全球新增病例趋势图如下： 从图中我们看到在H1N1传播的初期新增病例总数很小，在10天之后出现了较大的增幅，这可看做是H1N1的快速增长期，在传播了70天左右从图像上可以明显的看出，新增病例达到了一个高峰。随后，起伏的波峰趋于稳定，这是采取了有效的措施而产生的结果，可以看做疫情的传播得到了初步的控制。5.4 H1N1在墨西哥的传播特点分析墨西哥是全球H1N1疫情比较严重的国家之一，发现疫情早的

26、病人数中，所以选择墨西哥作为代表国家之一今年四月，墨西哥首次爆发疫情时，政府如临大敌，全国停班停课，公共场所停止营业，许多地方陷入空城状态，民众上街都以口罩掩面。但这种场面不会再出现，墨 西哥民众已经知道，在新流感在症状出现后立即服药就可以有效缓和病情。墨西哥当局还在研究制定未来停课的相关疫情标准。有关卫生单位预测，墨西哥在即将到来的冬天流感高峰中，新流感病例可能达到500万。到目前为止，墨西哥累计感染人数已经接近3万，其中有226例死亡病例。而全球的新流感累计病例已经超过30万，更有3900多人因为新流感丧命。根据我们搜集到数据，首先用MATLAB画出5月18日到7月24日这个时间段内H1N

27、1在墨西哥的发病数,如下图所示：根据得出的图形，可以看出H1N1在墨西哥蔓延趋势如下：在H1N1传播的初期就墨西哥就已经拥有了很多的病例，并且一直增长速度保持在一个较高的水平，不但没有减缓的趋势，在开始统计的第30天和第47天之后甚至出现了更大的增长率，这是一个小的爆发，和全球的趋势无法对应，说明这是局部的疫情爆发。从数据中分析墨西哥的H1N1的疫情防控工作做的不够到位，没能有效地控制住疫情的扩散。墨西哥感染H1N1而死亡的人数趋势图如下所示： 从图中我们可以看出感染发病总人数趋势有基本一致的变化，但整体上变化向后推迟了2周，这两周时间大致等于重症患者从发病到死亡的平均时间周期，与医学数据相统

28、一。疫情流行的初期死亡人数较多，随后人数变少，说明后来疫情得到了及时的控制，重症患者减少。H1N1在墨西哥的流行的新增病例如下图所示： 在对H1N1墨西哥死亡病例趋势图的分析中，我们可以发现图中的起伏非常大，这是由消息的交流不及时造成的，数据中很多天没有新增病例的及时通报，我们把视为0.另外后面新增病例越来越多，再一次验证了墨西哥的疫情没有得到有效控制。5.5 H1N1在日本的传播特点分析 根据我们搜集到数据，首先用MATLAB画出5月18日到7月24日这个时间段内H1N1在日本的发病数,如下图所示 . 从日本确诊病例趋势图中可以看出H1N1型流感在日本的流行起初较慢，甚至在出现一周后出现传播

29、停滞的情况，这主要归功于日本较为发达的医疗卫生保健系统以及日本政府对此次流感的重视，因此，在人口密度如此大的国家，病毒的传播并没有呈现爆发性的增长。随着时间的继续，日本确诊病例趋势图显示病例数在半个月后出现较大增长，并且呈现出指数型上升趋势，流感爆发。分析原因，主要归咎于日本的人口密度，人与人之间紧密接触造成流感急剧爆发。5.6 H1N1在美国的传播特点分析根据我们搜集到数据，首先用MATLAB画出5月18日到7月24日这个时间段内H1N1在美国的发病病例,如下图所示 从美国确诊病例趋势图中可以看出H1N1型流感在美国的流行今本呈指数函数变化。在病毒传播的初期病例的增加较为缓慢，随后增长率越来

30、越大，在开始统计的40天和45天，有较大的斜率值，这是爆发期，受H1N1感染的人数迅速增加，但在开始统计的约50天处，增长率明显放缓，说明美国政府的防控措施产生了效果。5.7 H1N1流行的新特点11月23日世界卫生组织最新公布的统计数据显示，全球已经出现至少6770例甲型H1N1流感患者死亡病例。报道前的一周，全球已知就有520人死于甲型H1N1流感。挪威公共健康研究所于11月19日宣布，在两名甲型H1N1流感死者和一名重症患者身上，发现变种的病毒株。11月21日又有媒体披露：英、美两国出现抗药性甲型H1N1流感病毒，其中，美国北卡罗来纳州报告的4例对“达菲”呈抗药性的甲型H1N1流感确诊病例已经有3例死亡，这是美国迄今人数最多的群发性抗药性甲流病例报告。可以看出抗药性甲流病例的出现对我们的防控体系提出了挑战，模型中本来治愈的病人会有可能再次被传染，这就为新一次的爆发创造了条件，所以我们应该加紧新型疫苗的研制。6.我们的建议根据对SIR模型的分析,