H7N9禽流感传播模型

-.z.H7N9禽流感传播模型摘要本文以2013年3月底再我国中东部地区爆发的H7N9新型禽流感为研究对象，通过对其传染和发展过程的调查和研究，来建立人类感染禽流感的微分方程模型。由于并未证实此类禽流感病毒是否具有人传人的特性，可以把人类患禽流感的通道归结为两类：一是带有禽流感病毒的禽传染人；二是感染禽流感病毒的人传人。我们在建立禽类传染病SEIR模型的基础上再建立人患有H7N9禽流感的SIR模型关键词H7N9禽流感微分方程禽类SEIR模型SIR模型一问题重述H7N9型禽流感是一种新型禽流感，于2013年3月底在和两地率先发现。H7N9型禽流感是全球首次发现的新亚型流感病毒，尚未纳入我国法定报告传染病监测报告系统，并且至2013年4月初尚未有疫苗推出。被该病毒感染均在早期出现发热等症状，至2013年4月尚未证实此类病毒是否具有人传染人的特性。2013年4月经调查，H7N9禽流感病毒基因来自于东亚地区野鸟和中国、、鸡群的基因重配。截至2013年5月6日16时，全国已确诊130人，31人死亡，42人痊愈。病例分布于、、、、、、、、等地。严重急性呼吸综合征（SevereAcuteRespiratorySyndromes），又称传染性非典型肺炎，简称SARS，是一种因感染SARS冠状病毒引起的新的呼吸系统传染性疾病。主要通过近距离空气飞沫传播，以发热，头痛，肌肉酸痛，乏力，干咳少痰等为主要临床表现，严重者可出现呼吸窘迫。本病具有较强的传染性，在家庭和医院有显著的聚集现象。首发病例，也是全球首例。于2002年11月出现在，并迅速形成流行态势。2002年11月－2003年8月5日，29个国家报告临床诊断病例病例8422例，死亡916例。报告病例的平均死亡率为9.3%。问题1：试根据H7N9型禽流感的传播规律，建立数学模型研究H7N9型禽流感的感染人数和死亡人数，并预测我国今后三个月的感染人数和死亡人数。问题2：根据你们建立的数学模型，假设在2002年底发生传染的病毒不是SARS，而是H7N9型禽流感，则情况会怎么样？问题3：H7N9型禽流感的出现对我国家禽业的影响巨大，试建立数学模型对我国家禽业的影响进行分析。问题4：根据上述模型和分析结果，从社会的稳定和提高人民的生活质量方面考虑，给我国政府写一篇建议文稿（1000字以），说明应该采取的手段。二问题分析问题1：由于尚未证实H7N9是否具有人传人的特性，我们可以假设人类感染禽流感来自两个方面：禽类感染人和人传人。设两者的传染系数分别为和。根据传染病的SEIR模型建立禽类禽流感传染模型，在此基础上来建H7N9禽流感传染模型。问题2：吧三模型建立与求解3.1问题1：建立禽类禽流感SEIR模型和人感染H7N9禽流感模型3.1.1禽类禽流感SEIR模型：3.1.1.1模型基本假设（1）2013年4月调查，H7N9禽流感病毒基因来自于东亚地区野鸟和中国，，鸡群的基因重配。本模型主要考虑鸡和野鸟，不考虑对其他动物的传染且假设鸡群和鸟类是在疫区均匀分布，由H7N9的特性，可以假设感染病毒后禽的死亡率为100%；（2）在禽流感传播期间，所考虑地区的总禽数不变，不考虑生死和鸟类的迁移；全部禽类可以分为对禽流感没有免疫力的易感染禽只（Susceptible），已感染病毒而未被发现的病毒潜伏禽只（E*posed），已经感染病毒并且被发现的受感染禽只（Infections）和已死亡禽只（Died）。任意t时刻4类禽只占总数的比列为，，，；（3）记时刻4种类别的数量分别为，，，；（4）患有禽流感的鸡和鸟类其传染病毒的能力一致；（5）每只已感染禽只每天有效接触易感染禽只的平均数量为，病毒潜伏禽只有效接触易感染禽只的平均数为，已死亡的禽只有效接触易感染禽只的平均数为；（6）易感染的禽只接触已感染禽只，死亡禽只，病毒潜伏禽只后有部分转为了已感染禽只有（1-）部分转化为了病毒潜伏禽只；（7）病毒潜伏禽只经过潜伏期后会转为已感染禽只，其占病毒潜伏禽只总数的比例为；（8）已感染禽只在一段时间后死亡，死亡数占已感染禽只数的比例为。如图所示四类禽只的转化关系：3.1.1.2SEIR模型的建立由模型中禽只总数不变可得（1）易感染禽只一部分转为病毒潜伏禽只一部分转为已感染禽只，由此我们可以得到（2）易感染禽只一部分转为病毒潜伏禽只，而病毒潜伏禽只一部分转为了已感染禽只，得到以下方程（3）病毒潜伏禽只一部分转为已感染禽只，而一部分已感染禽只死亡，可以得到（4）对于死亡禽只（5）联立以上方程3.1.1.3H7N9禽流感SIR模型仿真我们同样假设一下的参数值：禽对人的有效接触系数满足，病人对健康人的有效解除系数为满足。日治愈率，日死亡率。，。同样运用simulink对微分方程求解。3.1.2人感染H7N9禽流感模型的建立3.1.2.1模型假设与分析（1）在禽流感传播期间不考虑人口自然是出生率和自然死亡率，并假设疫区的人口总数为，人群可以分为健康者，病人，治愈者，死亡者。t时刻四类人在总人数所占比例分别为，，，。初始时刻四类人数分别为，，，；（2）禽流感患者治愈后会产生免疫作用；（3）同时考虑病禽传人和人传人是，禽流感的传播途径如图（4）假设每个病人每天有效接触健康者的平均数量为，每天被治愈的人数占病人总数的倍，每天死亡的病人数占病人总数的倍；（5）为了更好的表示感人禽只对人患禽流感的影响，这里引入疫情总指标函数，表示对其总体情况对人类的影响。其中其中，，。取影响系数，，。3.1.2.2模型建立有总人数的不变可得（1）健康者由禽只感染和病人的感染而变为病人，由此得到方程（2）一部分健康者转为了病人，而一部分病人被治愈或者死亡，有（3）被治愈者人数和死亡人数分别为，（4）联立上面的方程得到3.2问题2：3.2.1根据全国非典大事件时间表的不完全统计：2002年11月16日有5人丧命，300人被感染。2003年1月两名肺炎病人，该院接触过上述两名病人的医务人员中有八人发生同样疾病。2003年2月9日，305例非典型肺炎病例，死亡5例。2003年2月26日非典型肺炎首度发生在越南首都河。2003年2月21日，7名旅客感染2003年3月4日染病的大学第二附属医院退休教授剑伦不治去世。2003年3月25日省中医院护士长叶欣殉职2003年2月26日，20人感染。2003年3月14日，死亡1人2003年3月5日，在加拿大多伦多，一名姓关（译音）的妇女从回来后，因为患有非典型肺炎不治身亡。2003年3月6日，接报第一例输入性非典病例。2003年3月10日，22人感染。2003年3月29日，死亡1人。2003年3月11日，感染1人。2003年3月13日，感染的人数是47人。死亡1人。2003年3月15日，死亡1人，感染41人。2003年3月23日，感染1人。2003年3月25日，省中医院护士长叶欣殉职，她是抗“非典”战斗中第一位被患者传染而牺牲的医护人员。2003年3月26日，新华社首次报道疫情。2003年4月2日，中国省3月份有361起新病例，9人死亡。2003年4月3日，感染40名病人，死亡6例。2003年4月4日，感染100名病人，死亡7例。2003年4月6日，国际劳工组织芬兰籍男子派克"阿罗因患非典在京去世，成为首名在华得非典去世的外国人。2003年4月09日，可能感染100人。2003年4月13日，首批8位非典患者康复出院；2003年4月20日，卫生部承认疫情统计存在较大疏漏。2003年4月21日至4月底，非典疫情处于较高水平，连续十几天每天都有90至100个病例，最高一天达到150多人。2003年4月22日，最新疫情为非典病例588例，死亡28例，疑似666例。2003年4月23日，发现首例非典疑似病例。2003年4月24日，台北和平医院爆发封院事件，是第一间因SARS感染而封院的事件；而和平医院的封院方法犯下了“常识性的错误”。2003年4月26日，市委、市政府发布《关于保障医务人员全力以赴做好防治非典型肺炎工作的若干决定》。2003年4月27日，市委、市政府发布《关于加强防治非典型肺炎工作的决定》。全面建立非典患者的社会救助机制。暂停娱乐场所经营，并开始公布各区县疫情和隔离区信息。2003年4月28日，市政府发布《关于做好非典患者密切接触人员医学观察的通告》和《关于加强防治非典型肺炎工作的决定》。截至28日10时，累计收治非典病例2474人，其中确诊病例1199人，疑似病例1275人，治愈出院78人，死亡59人。2003年5月初开始，疫情出现小幅回落，但仍处于高发平台期。2003年5月4日，SARS疫情恶化，可能及疑似病例增加人数，在一天之均创下历史新高！世界卫生组织已派出专家组去台2003年5月6日《工伤保险条例》颁布。据此，医务人员在救治患者的过程中感染非典应算做工伤。2003年5月8日中国地报告有疫情的省份达25个。、、、、、**6个省区未发现疫情。三级甲等医院中日友好医院作为非典定点医院投入使用。2003年5月9日新增病例数首次跌破50例。市防治非典联合工作小组根据总书记的批示，决定打好五月攻坚战。同日宣布，医务人员的感染比例已经呈明显下降趋势。从4月21日到5月1日，每天平均是15.81人。从5月2日开始到5月8日，平均每天6.3人。2003年5月10日，6名外籍非典患者全部痊愈出院。2003年5月11日，疫情大幅下降已有两天。有人担忧，、等地疫情有向农村蔓延的趋势2003年5月13日，新发病例首次为零。2003年5月15日，小汤山非典定点医院第一批7名病人痊愈出院。2003年5月19日，通报的新增非典确诊病例首次降至个位数，总数为7例。省首次没有新增临床诊断病例。2003年5月21日，地坛医院送走第100位非典康复者。2003年5月22日起，8万名高三年级学生开始返校进行考前复习，其它年级的中小学生也将陆续分期、分批、分区域复课。体温超标者严禁入校。2003年5月23日和的科研人员联合宣布，从果子狸等野生动物体找到非典病毒前体。2003年5月24日，向、、、、、等11省市捐赠一批急需的防治非典物资。2003年5月29日，新收治非典确诊病例首次为零。当日确诊与疑似病例之和也首次降至个位数。2003年5月30日，淇宣布：防治非典“五月攻坚战”取得阶段性成果。2003年6月1日，首都高校出现首批返校的应届毕业生。2003年6月2日，非典疫情出现三个“零的突破”：当日新收治确诊病例、疑似病例转确诊病例、非典病例死亡人数均为零。2003年6月8日首次迎来新增非典病例零纪录。3.2.2模型的假设地总人数N可视为常数，即流入人口等于流出人口。（此模型取N=13000000，初始感染人数取以上统计的300，死亡数5）据人口所处的健康状态，将人群分为：健康者，SARS病人，退出者（被治愈者、免疫者和死亡者）。3．在政府的强制措施下，人口基本不流动，故无病源的流入和流出，避免了交叉感染，降低了感染基数。隔离的人断绝了与外界的联系，不具有传染性。SARS康复者二度感染的概率为0。国家完善了监控手段，加强了对SARS病毒监控的力度，故可假设所有感染SARS病毒的人群都进入了SARS病人类和疑似类。由于对SARS病原体的研究不够深入，无有效药物可以使人体免疫，同时SARS病毒感染后，大量繁殖，破坏免疫系统，故不可免疫。3.2.3参数的设定和符号说明s(t)：t时刻健康者在总体人群中的比例i(t)：t时刻SARS病人在总体人群中的比例l(t)：t时刻疑似病人在总体人群中的比例r(t)：t时刻被治愈者、死亡者和免疫者在总体人群中的比例之和。：SARS病人日接触率。为每个病人每天有效接触（足以使健康者受感染变为病人）的平均人数。：日治愈率。为每天被治愈的病人占病人总数的比例。：日转化率。为每天危险群体中的疑似病人被确诊为SARS患者的比例。：日死亡率。为每天SARS病人死亡的数量和当天病人总数量的比值。：疑似感染率。为每天感染为疑似病人的比例。3.2.4模型建立模型一感染为SARS患者情况由假设，每个病人每天可使个健康者变为病人，因为病人人数为，所以每天共有个健康者被感染，于是就是病人数的增加率，又因为每天被治愈率为，死亡率为，所以每天有个病人被治愈，有个病人死亡。则病人的感染为由于对于退出者()由假设可知：故SARS患者率模型一的方程建立如下：(3)模型二疑似患者的变化情况与前面同样的分析，得到疑似患者率模型二：(5)3.2.5模型求解3.2.5.1参数的确定：1.的确定=，=，=从统计结果给数据得：=0.055076，=0.038183，=0.002443。已确诊病例累计300现有疑似病例300死亡累计5治愈出院累计0当天退出数当天病人数0当天病例305退出率0治愈率03394021833174310.0394430.076566482610254365200.0115380.0826925886662846166190.0258480.0743136937823555136840.0190060.0804097748633964127740.0155040.082687877954427398730.0103090.0836298810934876109900.0101010.076768111412555678310650.0028170.0732391199127559781212100.0099170.0644631347135866831612910.0123930.0642911440140875901713880.0122480.06484115531415821001814540.012380.06877616361468911091115410.0071380.0707331741149396115715920.0043970.07223618031537100118616790.0035740.07028189715101031211717360.0097930.0697196015231071341018080.0055310.074115204915141101411318850.0068970.074801213614861121521819130.0094090.07945621771425114168919450.0046270.086375222713971161751519740.0075990.088652226514111201863119980.0155160.093093230413781292084120100.0203980.103483234713381342441319920.0065260.1224923701308139252619970.0030050.126189238813171402571720080.0084660.127988240512651412733820060.0189430.136092242012501453072719820.0136230.154894243412501473322019580.0102150.169561243712491503495019450.0257070.179434244412251543955418950.0284960.208443244412211564478318530.0447920.24123245612051585285617790.0314780.296796246511791605828817480.0503430.332952249011341636674116690.0245660.399641249911051677044416330.0269440.431108250410691687478515970.0532250.467752251210051728284115140.0270810.54689625149411758666314760.0426830.58672125178031769287914160.0557910.655367252076017710068513380.0635280.7518682521747181108797512530.7781320.86751825217391902053672780.2410077.38489225217341902120352110.16587710.0473925217241912154171760.09659112.2386425217181912171181590.11320813.6540925217161912189421410.29787215.524822521713191223126990.26262622.535352521550191225720730.27397330.9178125214511912277-116354-21.53742.16667252235118111243312170.0271160.92358325227118111573211840.0270270.9771962522418111897411520.0642361.0321182522318112635810780.0538031.171614252266818113218410200.0823531.295098252225718314031419360.1506411.498932252215518415431107950.1383651.940881252231841653966850.1401462.4131392522518617471985890.3361632.966044252241871944523910.1329924.971867252231891994213390.0619475.882006252231892015-575319-1.802516.3166142523218314463788940.4228191.61745252321861821565160.1085273.52907252321871876-0.289883.1184972确定很明显从我们建立的模型是无法得到s、i、、的解析解。为了解决这个问题我们用MATLAB软件中龙格—库塔方法求出他们的数值解。先通过实际统计数据算出每一天的s、i、、初步确定的围为1到2，通过调试。我们发现当1.4时，理论图形与实际图形有最佳的吻合。如图：matlab程序:functiony=ill(t,*)w=1.4;z=0.0575;v=0.00000003;u=0.01;y=[w.**(1).**(2)-z.**(1)-v**(1),-w.**(1).**(2)]'ts=0:0.01:70;*0=[300/13000000,1-300/13000000];[t,*]=ode45('ill2',ts,*0);[t,*];plot(t,*(:,1)),grid,pause（*(1):i(t);*(2):s(t)）3.2.5.1参数的确定：对比h7n9模型若2002至2003年全国出现禽流感，可以假设sars模型中的日接触率，即=1.4为禽流感禽与人之间的日接触率（因为到目前为止主要是H7N9只要是禽传染），即。代入到我们已经建好的H7N9禽流感传染模型除外，其他参数按照H7N9模型计算的假设给值。基于simulink得到仿真结果：对比两种仿真结果我们可以得到：（1）左图中病人数增长达最大值所需的时间少，健康者人数减少的更快，即假如在2002年底发生的是H7N9禽流感，则疾病爆发的时间更快。（2）左图中病人数达到最大值是其在人口总数总所占的比列更大说明假如在2002年底发生H7N9禽流感，在疾病爆发时刻患病人数更多。（3）作图中死亡者人数增长率比右图，且左图中最后的死亡人数比例比右图高，说明假如在2002年底发生H7N9禽流感死亡率比2013年高，且总死亡人也会多。3.33.3.1模型假设假设一：我国家禽业的发展是独立的，不受其他国家疫情、进出口量的影响。假设二：如果不考虑H7N9型禽流感的影响，我国家禽业的发展规律是自然延续的，不会发生特变化，其他因素的影响远小于H7N9型禽流感的影响。3.3.2符号说明：符号说明符号说明2013年全国禽肉产量发生禽流感时2013年家禽总产量当季无禽流感时家禽产量当季政府捕杀家禽数量无禽流感时一季度家禽产量发生禽流感时一季度家禽产量无禽流感时二季度家禽产量发生禽流感时一季度家禽产量无禽流感时三季度家禽产量发生禽流感时一季度家禽产量无禽流感时四季度家禽产量发生禽流感时一季度家禽产量一季度禽流感病毒致家禽死亡率一季度政府捕杀家禽数量二季度禽流感病毒致家禽死亡率二季度政府捕杀家禽数量三季度禽流感病毒致家禽死亡率三季度政府捕杀家禽数量四季度禽流感病毒致家禽死亡率四季度政府捕杀家禽数量一季度减少养殖规模比例三季度减少养殖规模比例二季度减少养殖规模比例四季度减少养殖规模比例3.3.3问题分析从网络资料来看，H7N9型禽流感的出现对我国家禽业造成巨大影响，损失巨大，这半个月疫情的到了有效控制，没有增加病毒感染者，家禽业有所恢复，在近期有望达到疫前水平。从历年的我国禽蛋产量的数据来看，今年来家禽业保持着良好的势头，中国家禽业的发展是有目共睹的。如果没有H7N9型禽流感，2013年的禽蛋产量仍然会保持持续增长。H7N9型禽流感只是一个突发事件，不会对旅游业造成长期的影响，不过在近几个月还有一个恢复过程。因此认为今后一段时期家禽业的影响因素有：（1）家禽业的在发展动力。如随着经济的发展，人民生活水平的提高，对禽蛋的需求会不断增加；随着我国养殖业的不断规壮大，产量会扩大。（2）H7N9型禽流感发生期间对家禽业的影响。两种因素的影响特点有：第一种是家禽业进步的在本质因素，具有根本性、长期性。第二种是突发事件，是偶然因素，具有暂时性，作用效果看有突出重要性。3.3.4模型的建立与求解3.3.4.1不考虑H7N9型禽流感对2013年禽蛋产量的预测。从2003-2012年的禽蛋产量数据，可以看到各年的分布有明显的规律性，见表1规律一：除2006年外，各年都在不断增长，从2003到2012年增加528万，平均增幅52.8万，见图1。规律二：各年增长率有一些波动，其中2006出现了负增长。其实有点波动是很正常的，因为事物都是在曲折中前进，螺旋式上升的，见图2。通过以上数据分析可以得到：2007年、2008年增长速度达到5%左右，但后几年发展比较稳健，市场机制逐步健全，年增长稳定在1.5%左右。故预测，2013年的禽蛋产量为2900万吨。根据历年禽蛋产量，现在需要预测2013年禽蛋产量：模型：线性回归模型用一条直线来近似地表示它们之间的关系，这条直线就叫做禽蛋产量m对运量*的回归直线，方程就叫做禽蛋产量m对时间的回归方程，记作

．

其中，称为回归系数．为了和禽蛋产量m的实际观测值区别，用m表示回归方程中的禽蛋产量，由的每一个取值所得到的的值叫做总成本的回归值．

一般地，如果随机变量与变量之间呈现*种线性关系，则与之间一元线性回归模型为

也称为变量对变量的一元线性回归方程，，称为回归系数．

建立模型与求解：

1、最小二乘法

根据样本值确定方程中系数，．

设在一次试验中，取得对数据，这对数据就是一组样本值，根据这一组样本值可以寻求一对系数，．但由于是一个随机变量，所以如果通过另一组试验又可得到一对，的值．也就是说，我们通过一组数据所得到的是系数，的估计值，记作，，通过一组试验数据所求出的回归方程为

，（2）

称为经验回归方程．又称为经验公式，，叫做经验回归系数．

为了求系数，的估计值，，我们常用的方法是最小二乘法．利用matlab求解求解程序为>>*=[12345678910];>>y=[2333237124382424252927022743276328112861];>>result=polyfit(*,y,1)result=1.0e+003*0.06422.2446故可以预测2013年的禽蛋产量为模型分析与检验：依据以上模型，估计近十年的禽蛋产量如表2所下。年份禽蛋产量/万吨增长率200323333.0200423711.6200524382.820062424-0.6200725294.3200827026.8200927431.5201027630.7201128111.8201228611.8表12003-2012年中国禽蛋产量及增长率图22003-2012年中国禽蛋产量图32003-2012年中国禽蛋产量增长率年份实际禽蛋产量预测禽蛋产量200323332308.8200423712373200524382437.2200624242501.4200725292565.6200827022629.8200927432694201027632758.2201128112822.4201228612886.62013--2950.8表22003-2012年中国禽蛋产量及增长率3.3.4.2考虑H7N9型禽流感对2013年禽肉、禽蛋产量的预测。假设2013年全国禽肉产量为，按分季产肉量计算，一季度约占35%，二季度和三季度各占20%，四季度占25%。正常情况，一季度可生产禽肉，二季度可生产禽肉，三季度生产禽肉，四季度可生产禽肉。但从各地区调查情况看，受禽流感的影响，禽流感病毒会导致家禽死亡率为a。在疫情严重的地区政府会有计划的捕杀家禽，设捕杀量为q。由于目前消费者对家禽安全性的担忧，禽产品消费需求大幅度下降，造成短期出现禽产品严重供大于求，禽产品价格大幅度下降。据各地反映，目前，肉禽价格下降为20-30%，高者达50%以上；禽蛋价格下降10-20%。这种情况将影响家禽生产后劲，降低养殖的积极性，设农户减少家禽苗数占养殖规模的比例为b。则按照以上分析得到，家禽数目为由于H7N9型禽流感发生在4月份，故对第一季度没有影响，第一季度家禽产量为在第二季度，尤其是4月底和五月初，疫情比较严重，影响家禽产量很大。病毒致死率a2=0.4，减少养殖规模比例为b2=0.2,捕杀量为q2设在5月底疫情得到控制时，第二季的家禽产量为在第三季，由于疫情基本消失，市场开始回暖。病毒致死率，减少养殖规模比例为,捕杀量为，第三季的家禽产量为在第四季，由于疫情消失，人们对家禽安全性的信任提高，对家禽的消费恢复到H7N9型禽流感发生前，养殖企业开始扩大养殖规模。病毒致死率，减少养殖规模比例为,捕杀量为，第三季的家禽产量为综上对2013年家禽产量建立如下数学模型根据（1）中模型，设，。则,,,;故2013年家禽产量为2568万吨用matlab分析如图3所示，各季度产量对比如图4所示。图3各季度产量图4各季度产量对比表3.3.5模型分析该模型的优点是充分考虑了H7N9型禽流感发生的初始期、严重期、消灭期等不同时期对我国家禽业的影响。该模型中各时期的影响系数是不同的，更贴近现实。该模型的不足之处是由于H7N9型禽流感还没完全消失，统计数据还不完善，对于各系数的取值还需进一步提高精度。3.4建议文稿根据我们所建立的模型可以知道，禽流感对人民身体健康造成了较大的威胁，造成了一定的恐慌，还严重影响了我国家禽业的正常发展，而且这种影响在一定时间还会持续存在。疫情发生后党和政府高度重视，采取了一系列的措施防控疫情，并且取得了较好的效果。但目前抗击H7N9禽流感病毒过程中仍然存在一些薄弱环节和不足之处：1.部分人对H7N9禽流感病毒缺乏科学的认识和应对措施。2.禽类产销受阻，损失严重。受H7N9禽流感疫情影响，家禽饲养和流通存在销售不畅，压栏严重，成本增加，损失严重等状况。3.抵御风险能力不足，自救乏力。家禽业主要以出售活禽和鲜蛋为主，自身加工、储