禽流感H7N9传播模型的动力学分析

周飘飘1, 祝光湖1,2

(1. 桂林电子科技大学 数学与计算科学学院, 广西 桂林 541004; 2. 广东省疾病预防控制中心 广东省公共卫生研究院, 广州 511430)

摘要 H7N9型禽流感严重威胁人类健康和生命安全为研究H7N9病毒的传播规律,提出了一个结合人群、家禽和环境中病毒之间相互作用的SI-V-SEIR禽流感传染病模型通过动力学分析,给出基本再生数R0的表达式,并证明无病平衡点和地方病平衡点的稳定性接着应用模型分析广东省2016年—2017年的H7N9疫情,获得疫情初期R0=18.8,此时禽类的接种率需达到94.7%才能控制病毒在禽类和环境中的传播,而采取措施后R0=0.14结果表明,降低环境中的病毒载量、和禽类之间以及禽到人的传染率能有效地减少染病人数

H7N9; 动力学模型; 稳定性; 基本再生数

引 言

自2013年2月华东地区发现首例人感染新型H7N9流感病毒病例以来[1],H7N9流感已经在中国出现5波疫情截至2017年6月,全国已报告1 533人确诊感染H7N9流感病毒,其中592人死亡疫情迅速蔓延,感染和死亡人数不断增加,引起了全球公共卫生界的警惕广东省是H7N9流感传播的高危地区[2],截止2017年6月,广东省共有258例人感染H7N9流感病例,主要发生在每年的10月份到次年4月份由于目前尚无有效治愈H7N9的方法,死亡率较高,对人群健康和社会经济造成了严重影响,因此研究其传播机制具有十分重要的现实意义

目前已有部分数学模型从不同的侧重点对H7N9进行研究例如,2014年,Zhang等[3]建立了一种具有候鸟、居民禽、家禽和人的动力学模型,发现候鸟很可能为长江三角洲地区2013年—2014年疫情的最初感染源,估计基本再生数为4.11~6.022016年,Lin等[4]建立了几种传播模型,发现有环境类感染的不含周期效应的模型最适合模拟H7N9的传播模式同年Lu等[5]采用Bayes(贝叶斯)模型分析人群发病率的时空变化,评估了关闭和重开活禽市场对疫情的影响2017年,Xing等[6]建立了一种具有候鸟、居民禽、家禽和人的动力学模型,探索温度和周期性关闭活禽市场对H7N9病毒的影响很多研究表明,接触家禽或者环境中的H7N9病毒是人类感染禽流感的最主要途径[7-8]为此,我们将进一步考虑禽类、环境中的病毒和人的相互作用,提出一个新的SI-V-SEIR模型

本文将从理论和数值上分析此模型的动力学行为,并把模型应用于广东省2016年—2017年的禽流感疫情,分析该地区在关闭活禽市场后的基本再生数变化,并讨论不同的模型参数对传播动力学的影响本文第1节详细叙述模型,给出了基本再生数的表达式,并证明模型的无病平衡点和地方病平衡点的稳定性;第2节给出模型的应用,通过选取广东地区的疫情数据进行数值分析;最后给出结论

1 数 学 模 型

建立一个简单的H7N9传染病模型其中,Np(t)和Nh(t)分别表示禽类和人类的总数量,Np(t)分为两个子类:易感禽类Sp(t)和感染禽类Ip(t);Nh(t)分为4个子类:易感者Sh(t)、潜伏期患者Eh(t)、感染者Ih(t)和康复者Rh(t)人类和家禽可通过与受感染家禽接触或环境中的病毒(V)感染H7N9病毒变量V的单位是病毒粒子[9],用于模拟被感染的家禽向环境中传播H7N9病毒的演变V源于受感染家禽,存活于环境中[3,9]疾病传播的流程图如图1示

图1 H7N9传播模型的流程图
Fig. 1 The flowchart of the H7N9 transmission model

对应的动力学模型可由下列高维微分方程组表示:

(1)

模型参数定义见表1(具体值对应广东省情况)感染力λ11λ12λ21λ22定义如下:

因目前没有直接证据证明H7N9可以人传人[6,10-11],故假设人类不参与病原体的传播模型动力学行为主要由下列子系统决定:

(2)

注意到禽类总数Np满足故禽类总数的平衡点是全局渐近稳定,由此获得式(2)的极限系统:

(3)

系统(3)的稳定性决定了系统(1)和(2)的稳定性,因此下面的理论分析主要针对系统(3)根据模型对应的生物意义及简单计算,方程(3)的正向不变集为

表1 模型参数及取值

Table 1 Parameters and their values

parameterdefinitionvaluenoteArecruitment rate of poultry14 580 000[a]bhbirth rate of human1/(75×52) week-1[b]dhnatural death rate of human(1/(75×52))×0.98 week-1[b]λpbaseline transmission rate among poultryvariableλhbaseline transmission rate from poultry to humanvariableβpbaseline transmission rate from virus in VE to poultryvariableβhbaseline transmission rate from VE to humansvariableα1rescaled constant of environmental infectiousness to poultry0.034 6 week-1estimationα2rescaled constant of environmental infectiousness to humans0.024 4 week-1estimationdprate of screening and culling in poultry1/6 week-1[c]dvelimination rate of VE in natural environment0.961 0 week-1estimationrproduction rate of VE from poultry0.143 3 week-1estimationγreciprocal of the incubation period7/3.3 week-1ref. [3]ηdisease-related death of human0.96 week-1ref. [3]δrecovery rate of infectious human1 week-1ref. [9]

[a]为由家禽出栏率估计所得

[b]为假设人均寿命为75岁

[c]为由现实中实测估计所得

Note [a] estimation from the marketing rate of domestic poultry.

[b] with an average life time of human taken as 75 years.

[c] estimation from fild survey.

2 基本再生数与稳定性

基本再生数R0作为衡量疾病爆发或消亡的阈值,有重要的生物意义[12]根据文献[13]提出的下一代生成矩阵的方法,得

从而

其中基本再生数是矩阵FV-1的谱半径,因此

2.1 无病平衡点的全局渐近稳定性

令系统(2)的右边为零,求得模型的无病平衡点(Sp,Ip,V)=(A/dp,0,0),注意到系统(3)的等号右边可看作一个向量值函数:

(4)

定理1R0<1时,模型(1)的无病平衡点局部渐近稳定

证明 系统(3)在无病平衡点的Jacobi矩阵为

特征方程为

(5)

下面只需证明式(5)的解有负实部即可

则特征方程(5)可简写为

λ2++v=0

(6)

方程(6)的解为从而J0的所有特征值均为实数若要使得特征值都是负数, 应有u2<Δ,而u2-Δ=4(mdv-α1βpNp)另外R0<1等价于因此当R0<1时,Jacobi矩阵的特征值全为负实数,无病平衡点是局部渐近稳定的

进一步地有下面的结论

定理2R0<1时,模型(1)的无病平衡点全局渐近稳定,此时疾病灭绝

证明 由于R0<1,所以

(7)

由于

因此系统(7)只有平衡点(0,0)由Dulac定理可知,系统(7)是全局渐进稳定的,则当t→∞时,有(Ip1(t),V1(t))→(0,0)由于Ip(t)≤Ip1(t),V1(t)≤V(t),所以当t→∞时,有(Ip(t),V(t))→(0,0)那么系统(2)中Sp(t)→A/dp,即无病平衡点是全局吸引的,同时无病平衡点又是局部稳定的,所以系统(2)中无病平衡点是全局渐近稳定的证毕

此时,环境中的病毒和感染禽类的数量均趋向于零,代入系统(1)可知,模型中的感染病人也趋向于零,疾病将灭绝

2.2 地方病平衡点的全局渐近稳定性

为了分析地方病平衡点的稳定性,需要一个引理(文献[14]中的推论3.2)

引理1 设函数连续可微,满足

1) f(x)是协作的,即它的Jacobi矩阵是Metzler矩阵,且∂f(x)不可约;

2) f(0)=0,且对所有的i=1,2,…,n,当xi=0时fi(x)≥0;

3) f(x)严格次线性

s(∂f(0))是矩阵∂f(0)的特征值的最大实部,如果s(∂f(0))>0,那么对系统而言,要么系统的任意解趋向无穷;要么存在唯一的正平衡点,全局渐近稳定

定理3R0>1,模型(2)存在唯一的地方病平衡点且全局渐近稳定,此时病毒将长期存在于禽类、环境和人群中

证明 1)系统(3)的Jacobi矩阵为

这里的f(x)由式(4)定义,其中

容易看出它的非对角元均非负,因此是Metzler矩阵A=df(x)为邻接矩阵,由于A所对应的可达性矩阵Q中全部元素均为1,所以∂f(x)不可约[15]

2)f(0)=0是显然成立的;当xi=0时,即当Ip,V分别为0时,则

易验证fi(x)≥0,i=1,2

3) 即对∀α∈(0,1),有

f2(αx)-αf2(x)=rIp(1-α)>0

f(αIp,αV)>αf(Ip,V),因此是严格次线性的

注意到∂f(0)=F-V,而文献[13]指出R0>1等价于s(F-V)>0因为系统(3)有紧正向不变集Ω意味着它的解不会趋向于无穷,故由引理1,系统(3)存在唯一的正平衡点全局渐近稳定,对应地,系统(2)存在唯一的地方病平衡点它全局渐近稳定设地方病平衡点稳定时人类的稳态人口为于是由平衡点关系获得

其中

由此可见的值被唯一确定下来

以上结论表明,当R0>1时,无病平衡点不稳定,系统(1)存在唯一正平衡点,H7N9病毒感染者将一直存在,并且感染人数逐渐稳定,形成地方病

3 模 型 应 用

3.1 参数选取

选取广东省2016年11月—2017年6月H7N9的周病例数(从香港特别行政区政府卫生防护中心获取),在此期间总发病数为63例,发病率为0.057 3/10 000初值由广东区疫情数据估计而得,模型以周为单位根据广东省统计信息网,估计禽类的扑杀率dp=1/6和禽类载量A=14 580 000,再利用假设禽类初始总数Np(0)=17 948 000从2017年1月开始,政府采取了多种控制措施,特别是加强了对活禽市场管理(如清理和关闭),从而降低了疾病的传染率,因此将传染率(λpλhβpβh)分成两段代入模型,并利用Markov(马尔科夫)链Monte-Carlo(蒙特卡罗)估计未知参数,得到图2模型的拟合图可看出, H7N9疫情在2017年1月左右最为严重,随后感染人数逐渐下降,疫情逐渐消失

3.2 数值分析

运用MATLAB软件进行数值模拟,进一步分析动力学性质分别对初值Ip(0)和V0(0)以及参数α1,α2,rdv取不同值,研究初始时刻染病禽类数量和环境病毒载量、以及环境病毒的传染率、染病家禽对病毒生产率和环境病毒的消除率对人群感染的影响,模型其他参数的取值基于广东实际情况进行选取和估计

分别将初始时刻的病毒数量初值V(0)和感染病毒数Ip(0)视为变量,其他参数不变,感染病人的时间演化图如图3(a)所示结果发现模型峰值随着Ip(0)和V(0)的不断增大而升高,即感染人数逐渐增多,在2017年1月初达到峰值,之后病例数迅速下降且逐渐趋于零通过计算估计在2016年11月—2017年1月初期间,模型的基本再生数R0=18.8,这期间疾病流行并迅速蔓延此时需要对家禽的接种率为V/N=1-1/R0=94.7%,从而可控制H7N9在禽类和环境中的传播疫情的蔓延使得广东省各地迅速采取了各类措施,如关闭市场和清理环境等,对应的降低了传染力,而后的基本再生数下降为R0=0.140 4<1,此时H7N9疫情随着时间的推移逐渐消亡

图2 2016年11月—2017年6月广东省H7N9的周病例数和模型的拟合图
Fig. 2 The fitting results of the H7N9 human cases from Nov. 2016 to Jun. 2017 in Guangdong with the model

(a) 不同的感染禽类初值的影响 (b) 不同的环境中病毒初值的影响
(a) Different initial values of infected poultry (b) Different initial values of virus in the environment
图3 在不同的初值条件下感染人数的时间演化图
Fig. 3 The time evolution of the infected human number with different initial values

比较图3发现,感染禽类和环境中的病毒初值对感染H7N9人数都有很大的影响,感染禽类比环境中病毒对发病率影响更加显著,这与Breban等[9]分析2009年疫情得到的环境中病毒传播水平较低的结果相符合

(a) 环境病毒对禽类的传播参数α1对病例数的影响 (b) 环境病毒对人类的传播参数α2对病例数的影响
(a) The effects of parameter α1of rescaled environmental (b) The effects of parameter α2 of rescaled environmental infectiousness to poultry on the human incidence infectiousness to human on the human incidence

(c) 家禽对环境中病毒的生产率r对病例数的影响 (d) 病毒死亡率dv对病例数的影响
(c) The effects of parameter r of VE production rate (d) The effects of parameter dv of VE elimination rate in from poultry on the human incidence natural environment on the human incidence
图4 模型参数对感染人群的影响
Fig. 4 Effects of model parameters on the human incidence

保持初值以及其他参数不变,分别改变α1α2rdv,得到的人群感染数量演化图如图4所示,发现模型峰值随着α1α2r的增大及dv的减少而升高,此时感染病例数逐渐增多

4 结 论

本文提出并分析了一个含有环境中病毒的H7N9传染病模型,运用了微分方程稳定性理论和数值模拟对模型进行讨论证明了当基本再生数小于1时系统(1)存在无病平衡点且全局渐近稳定,此时疾病随时间的推进将逐渐消亡若基本再生数大于1,系统(1)存在全局稳定的正平衡点,H7N9传染病将发展为地方病

通过数值分析发现,降低初始时刻环境中的病毒和感染禽类数量能有效减少人类的感染,但仅仅降低初值还不足够,后期还要通过降低家禽间的传染率(例如对染病禽类进行扑杀),降低禽类到人类的传染率(例如实施关闭活禽市场),以及减少环境中病毒对禽类的传染率和对人类的传染率(例如减少去比较密集的地方,离养殖场较近的地方),减少感染的家禽对环境中病毒的生产率(例如对染病家禽进行焚烧和掩埋,防止病毒进入环境),均可以有效地减弱H7N9的传播广东省对控制禽流感实施了许多防控措施,例如,“1110政策”“集中屠宰、冷链配送、生鲜上市”和“定期关闭活禽交易市场”等,使得禽流感疫情得到了有效的控制,但H7N9仍有卷土重来的风险,仍需及时进行风险评估和多加防范

文中的模型还有许多因素没有考虑,例如温度和湿度等环境因素对H7N9病例数的影响,我们将在以后的工作中做进一步研究与讨论

参考文献

[1] GAO R, CAO B, HU Y, et al. Human infection with a novel avian-origin influenza a (H7N9) virus[J]. The New England Journal of Medicine, 2013, 368(20): 1888-1897.

[2] 徐继承, 黄水平, 肖伟伟, 等. 中国大陆地区438例人感染H7N9禽流感空间聚集性分析[J]. 中华流行病学杂志, 2014, 35(11): 1270-1274. (XU Jicheng, HUANG Shuiping, XIAO Weiwei, et al. Spatial aggregation of 438 human infections with avian influenza A (H7N9) in the mainland of China[J]. Chinese Journal of Epidemiology, 2014, 35(11): 1270-1274.(in Chinese))

[3] ZHANG J, JIN Z, SUN G Q, et al. Determination of original infection source of H7N9 avian influenza by dynamical model[J]. Scientific Reports, 2014, 4: 4846. DOI: 10.1038/srep04846.

[4] LIN Q, LIN Z, CHIU A P Y, et al. Seasonality of influenza A (H7N9) virus in China: fitting simple epidemic models to human cases[J]. PloS One, 2016, 11(3): e0151333. DOI: 10.1371/journal.pone.0151333.

[5] LU J, LIU W, XIA R, et al. Effects of closing and reopening live poultry markets on the epidemic of human infection with avian influenza A virus[J]. Journal of Biomedical Research, 2016, 30(2): 112-119.

[6] XING Y, SONG L, SUN G Q, et al. Assessing reappearance factors of H7N9 avian influenza in China[J]. Applied Mathematics and Computation, 2017, 309: 192-204.

[7] ZHOU X, LI Y, WANG Y, et al. The role of live poultry movement and live bird market biosecurity in the epidemiology of influenza A (H7N9): a cross-sectional observational study in four eastern China provinces[J]. Journal of Infection, 2015, 71(4): 470-479.

[8] CHEN L J, LIN X D, GUO W P, et al. Diversity and evolution of avian influenza viruses in live poultry markets, free-range poultry and wild wetland birds in China[J]. Journal of General Virology, 2016, 97(4): 844-854.

[9] BREBAN R, DRAKE J M, STALLKNECHT D E, et al. The role of environmental transmission in recurrent avian influenza epidemics[J]. PLoS Computational Biology, 2009, 5(4): e1000346. DOI: 10.1371/journal.pcbi.1000346.

[10] WEN Y M, KLENK H D. Erratum: H7N9 avian influenza virus: search and re-search[J]. Emerging Microbes & Infections, 2013, 2(5): e26. DOI: 10.1038/emi.2013.26.

[11] ZHOU P, MA J, LAI A, et al. Avian influenza A (H7N9) virus and mixed live poultry-animal markets in Guangdong province: a perfect storm in the making? [J]. Emerging Microbes & Infections, 2015, 4(10): e63. DOI: 10.1038/emi.2015.63.

[12] 马知恩, 周义仓, 王稳地, 等. 传染病动力学的数学建模与研究[M]. 北京: 科学出版社, 2004.(MA Zhien, ZHOU Yicang, WANG Wendi, et al. The Mathematical Modeling and Study of Infectious Disease Dynamics[M]. Beijing: Science Press, 2004.(in Chinese))

[13] VAN DEN DRIESSCHE P, WATMOUGH J. Reproduction numbers and sub-threshold endemic equilibria for compartmental models of disease transmission[J]. Mathematical Biosciences, 2002, 180(1): 29-48.

[14] ZHAO X Q, JING Z J. Global asymptotic behavior in some cooperative systems of functional differential equations[J]. Canadian Applied Mathematics Quarterly, 1996, 4(4): 421-444.

[15] 王伟贤. 矩阵不可约的充要条件[J]. 数学的实践与认识, 1988(4): 42-46.(WANG Weixian. The necessary condition for the matrix irreducibility[J]. Mathematics in Practice and Theory, 1988(4): 42-46.(in Chinese))

Dynamic Analysis of the Avian Influenza A (H7N9) Transmission Model

ZHOU Piaopiao1, ZHU Guanghu1,2

(1. School of Mathematics and Computing Science, Guilin University of Electronic Technology, Guilin, Guangxi 541004, P.R.China; 2. Guangdong Provincial Institute of Public Health, Guangdong Provincial Center for Disease Control and Prevention, Guangzhou 511430, P.R.China)

Abstract: Avian influenza A(H7N9) is always a big threat to human health and safety. Aimed at the transmission patterns of A(H7N9), a new SI-V-SEIR epidemic model was put forward, which incorporated the viral interactions among humans, poultry and environment. Through dynamic analysis, the expression of basic reproduction number R0 was given, and the stability of disease-free and endemic equilibrium points was proved. The proposed model was further applied to study the 2016—2017 A(H7N9) outbreaks in Guangdong province. It is found that R0=18.8 in the early outbreak, which indicates 94.7% of poultry to be vaccinated for the control of the virus transmission in poultry and environment. After control, R0 will fall down to 0.14. The results show that, reduction of the viral load in environment and the infection ratios among poultry and from poultry to humans could effectively lower human infections.

Key words: H7N9; dynamic model; stability; basic reproduction number

中图分类号 O29; O175.1

文献标志码:A

DOI: 10.21656/1000-0887.390175

收稿日期 2018-06-25;

修订日期:2018-09-19

基金项目 国家自然科学基金(11661026);广州市科技计划(201804010383);广西自然科学基金(2017GXNSFAA198235)

作者简介

周飘飘(1993—),女,硕士生(E-mail: 1150732811@qq.com);

祝光湖(1979-),男,副教授,博士(通讯作者. E-mail: ghzhu@guet.edu.cn).

文章编号1000-0887(2019)03-0311-10

ⓒ 应用数学和力学编委会,ISSN 1000-0887

Foundation item: The National Natural Science Foundation of China(11661026)

引用本文/Cite this paper:周飘飘, 祝光湖. 禽流感H7N9传播模型的动力学分析[J]. 应用数学和力学, 2019,40(3): 311-320.ZHOU Piaopiao, ZHU Guanghu. Dynamic analysis of the avian influenza A (H7N9) transmission model[J]. Applied Mathematics and Mechanics, 2019, 40(3): 311-320.