乙型病毒性肝炎是由乙肝病毒(hepatitis B virus,HBV)引起的、以肝脏炎性病变为主要表现的一种疾病[1]。据世界卫生组织报道，全球有超过2亿4千万长期慢性乙肝患者，每年超过78万人死于HBV相关性疾病。在中国，HBV是目前3大传染性疾病之一。中国的HBV患者人数接近世界HBV患者人数的1/3[2, 3]，每年有30万人死于HBV相关性疾病，其中18万人死于肝癌[4]。
为了更好地监测和控制传染病疫情，中国国家卫生和计划生育委员会从2003年开始定时发布国家法定传染病疫情分析[5]。针对2002-2009年的中国急性乙肝患者数据，文献[6]构建了数学模型并拟合得到了中国HBV传播的基本再生数R0=2.406，说明中国HBV的流行趋势仍不容乐观。但根据2009年之后的卫计委数据，我们发现2009-2013年中国乙肝患者数的整体趋势是逐渐下降的(图 1)，说明中国乙肝疫情已得到一定缓解，主要原因可能是2002年启动的中国卫生部/全球疫苗免疫联盟(global alliance for vaccines and immunization,GAVI)乙肝疫苗合作项目以及2006年实行的《2006-2010年全国乙型病毒性肝炎防治规划》等防控措施效果的凸显[7, 8]。特别是随着2002年新生儿HBV计划免疫政策的实施，现在婴幼儿感染HBV的比例已经越来越低[9, 10]，并且文献[11]也发现母婴传播对中国HBV传播的影响极小。因此，本研究拟在文献[6]的基础上去除母婴垂直感染的影响，重点关注水平感染并结合新的报道数据对模型进行校对分析，预测中国乙肝疫情的流行趋势，并对现有防控措施进行评估，进而提出下阶段的防控重点，为控制和降低我国HBV流行提供量化支撑。
1 数据与方法 1.1 数据根据国家法定传染病疫情报道，2002-2013年中国大陆(不包含香港、澳门和台湾)的HBV急性患者数据如图 1所示。
 |
| 图 1 2002-2013年中国急性乙肝患者数据分析 |
根据HBV感染的阶段特征，本研究将人群分为6个组：易感人群(S)，潜伏期人群(L)，急性患者(I)，携带者(C)，免疫人群(V)和康复人群(R)，并且各个组别间的传播流程图如图 2所示。由图 2所示，我们发现除了由急性患者和携带者向康复人群流入之外，康复人群与其他组别间并无其他的交互作用，因此HBV传播的数学模型可以用如下5个常微分方程来表示，模型参数的含义及其取值如表 1所示，其中S(2002)、L(2002)、I(2002)、C(2002)和V(2002)分别为HBV易感人群、潜伏期人群、急性患者、携带者和免疫人群数占2002年中国总人口数的比率。
 |
| 图 2 HBV传播流程图 |
| 参数 | 说明 | 数值 | 参考文献 |
| μ | 人口出生率 | 0.012 1 | 16 |
| μ 0 | 自然死亡率 | 0.007 15 | 16 |
| μ 1 | HBV相关疾病死亡率 | 0.2% | 8 |
| β | HBV急性患者传播系数 | 0.95~20 | 未知 |
| ε | HBV携带者传播权重 | 0.16 | 11 |
| σ | HBV潜伏期患者向急性患者转化的系数 | 6 | 12 |
| γ 1 | HBV急性患者的转化系数 | 4 每年 | 11 |
| γ 2 | HBV携带者的康复系数 | 0.025每年 | 11 |
| γ 3 | HBV疫苗接种率 | 0%~100% | 未知 |
| ω | 未接种乙肝疫苗的新生儿比率 | 0%~100% | 未知 |
| q | HBV急性患者未能清除感染的平均概率 | 0.0885 | 13,14,15 |
| ψ | HBV免疫人群的失效率 | 0.1 | 12 |
| S(2002) | HBV易感人群的初始值 | 0~1 | 未知 |
| L(2002) | HBV潜伏期人群的初始值 | 0~1 | 未知 |
| I(2002) | HBV急性患者的初始值 | 5.202 9×10 -4 | 16 |
| C(2002) | HBV携带者的初始值 | 0.090 9 | 16,17 |
| V(2002) | HBV免疫人群的初始值 | 0.394 7 | 16 |
基本再生数表示一个HBV患者在其平均患病周期内所有传染的人数。利用房室模型基本再生数的计算方法[18, 19]，对模型(1)可以计算得到HBV基本再生数为
其中
分别表示1例HBV急性患者在其平均患病周期内所传染的人数，即HBV急性患者的基本再生数，和1例HBV携带者在其平均患病周期内所传染的人数，即HBV携带者的基本再生数。
表示当HBV疫情绝灭时易感人群的数量，z1=μ0+σ,z2=μ0+γ1,z3=μ0+μ1+γ2。
敏感性分析一直被用来检测模型参数对模型输出的影响。本研究首先用拉丁超立方抽样(latin hypercube sampling,LHS)和偏等级相关系数方法(partial rank correlation coefficients,PRCCs)来检测输入参数对基本再生数的影响程度从而确定影响HBV疫情的重要参数。LHS也叫作无重叠的多维分层抽样，抽样的参数值随机分成N个相等的抽样概率间隔。N表示样本量的大小，通常远远大于变化参数的个数以保证准确性[20]。通过用等级转换的LHS矩阵和输出矩阵计算PRCC，PRCC反映了这些参数的不确定性对输出结果的影响，并可据此评估参数变化对输出结果的影响程度[19, 20]。由于对相关参数的实际分布情况知之甚少，对已知参数本研究选用±20%参数值的均匀分布来模拟相应参数的分布情况[20, 21, 22]。根据PRCC的定义，如果某参数的PRCC值为正(负)，那么随着该参数的增加(减小)，对应的基本再生数也将增加(减小)，并且当某参数的PRCC绝对值>0.2表明该参数对输出结果影响显著[20, 21]。
然后，将对基本再生数R0影响显著的参数分别优化20%，并利用如下公式计算每个参数优化时基本再生数的改变百分数
其绝对值的大小顺序即可刻画各参数对基本再生数影响大小的顺序。其中p表示参数，R0(p)表示参数未优化时的基本再生数，R0(
)表示参数优化20%后的基本再生数。
利用遗传算法，以拟合数据与报道数据的误差平方和最小为目标函数，得到模型(1)中待估参数的估计值，再以此估计值为基础，利用自适应Metropolis-Hastings 算法估计出各参数的基线值及其95%的置信区间[23]，结果见表 2。
| 参数 | 基线值 | 95% CI |
| β | 3.192 00 | (3.126 50,3.257 50) |
| γ 3 | 0.454 04 | (0.446 56,0.461 52) |
| ω | 0.552 50 | (0.545 69,0.559 31) |
| S(2002) | 2.75 12×10 -2 | — |
| L(2002) | 7.03 56×10 -4 | — |
根据表 1和表 2中的参数值，我们对模型(1)进行了数值模拟，再对得到的2002-2013各年度模型预测数据利用文献[24]提出的方法计算模型的拟合优度指标，结果得RNL=0.957 5。根据国家卫计委公布的2014年中国乙肝患者数据(I2014=7.284 4×10-4)，同模型预测值(
2014=8.002 5×10-4)比较，我们发现模型的拟合相对误差仅为9.8581%。因此，模型(1)可较好地模拟中国乙肝患者数据(图 3A)，可用来对中国乙肝的流行趋势进行预测分析。图 3B模拟了模型(1)预测2002年至2102年间HBV急性患者、携带者和免疫人群的流行趋势。结果发现在2014年之后，HBV急性患者和携带者人数均呈下降趋势，而免疫人群数目呈急剧上升趋势。说明GAVI乙肝疫苗合作项目和《2006-2010年全国乙型病毒性肝炎防治规划》等的实施使得大量易感人群转化为免疫人群，HBV疫情得到较好控制。
 |
| A：模型(1)与报道的HBV急性患者拟合结果；B：模型(1)预测2002-2102年间HBV急性患者、携带者和免疫人群的流行趋势，灰色区域表示95%的置信区间 图 3 模型拟合与乙型病毒性肝炎流行趋势分析 |
利用表 1和表 2的参数值可计算得出中国HBV传播的基本再生数R0=0.671 2(95%置信区间：0.659 1～0.671 8)，其中RI=0.252 5(95%置信区间：0.247 5～ 0.250 4)，RC=0.418 7(95%置信区间：0.410 9～0.422 0)。 因为R0 < 1，由Routh-Hurwitz判据可知此时模型的无病平衡点(S0,0,0,0,0)是稳定的，即当时间t充分大时，图 3B中的C(t)和I(t)都将趋于零。同文献[8]比较，我们发现HBV基本再生数从2009年的2.406大幅下降到了目前的0.671 2，说明中国HBV疫情已得到较好控制。注意到HBV携带者的基本再生数RC= 0.418 7大于HBV急性患者的基本再生数RI=0.252 5，说明目前中国乙肝传播的主要危险在于携带者。因此为了继续有效控制HBV疫情，需要重点关注对HBV携带者防控，切实降低HBV携带者对易感人群的传播风险。
2.3 敏感性分析由于在现有的社会条件和治疗水平下人口出生率μ、自然死亡率μ0、HBV潜伏期患者向急性患者转化的系数σ和HBV急性患者的转化系数γ1等参数变化很小或很难被人为调控，我们仅对模型的其余8个参数(表 3)关于基本再生数进行敏感性分析。设定抽样样本量N为1 000，参数独立，并且相应的参数基线值及其分布如表 3所示。
| 参数 | 基线值 | 分布类型 | 抽样区间 |
| μ 1 | 0.2% | 均匀分布 | (0.001 6,0.002 4) |
| β | 3.192 0 | 正态分布 | (3.126 5,3.257 5) |
| ε | 0.16 | 均匀分布 | (0.128,0.192) |
| γ 2 | 0.025 | 均匀分布 | (0.02,0.03) |
| γ 3 | 0.454 04 | 正态分布 | (0.446 56,0.461 52) |
| ω | 0.552 50 | 正态分布 | (0.545 69,0.559 31) |
| q | 0.088 5 | 均匀分布 | (0.070 8,0.106 2) |
| ψ | 0.1 | 均匀分布 | (0.08,0.12) |
图 4A、4B和4C分别为由LHS/PRCC方法计算各参数对HBV基本再生数R0、HBV急性患者基本再生数RI和HBV携带者基本再生数RC影响的PRCC结果示意图。由此我们发现HBV急性患者基本再生数RI主要对急性患者传播系数β和疫苗接种率γ3参数敏感，而HBV携带者基本再生数RC主要对携带者的康复系数γ2、携带者传播权重ε和急性患者未能清除感染的平均概率q参数敏感，并且免疫人群的失效率ψ参数对这两个基本再生数均显著敏感。根据各参数的PRCC值大小及其符号含义，如果进一步增大携带者的康复系数γ2和疫苗接种率γ3参数，并减小患者传播系数β、携带者传播权重ε、急性患者未能清除感染的平均概率q和免疫人群的失效率ψ参数，就能更好地控制HBV疫情。此外，R0与RI和RC的马尔科夫链的箱型图(图 4D)表明HBV携带者对HBV传播的影响远大于急性患者的影响，但其整体基本再生数R0始终小于1。再次说明现阶段控制HBV疫情的重点在于对HBV携带者的防控。
 |
| A：各参数对HBV基本再生数R0的PRCC值；B：各参数对HBV急性患者基本再生数RI的PRCC值；C：各参数对HBV携带者基本再生数RC的PRCC值；D：R0与RI和RC的马尔科夫链的箱型图，灰色区域表示PRCC值不显著区域 图 4 基本再生数的敏感性分析结果 |
由图 4A可知，共有6个参数对模型的基本再生数R0影响显著，图 5表示分别对它们优化20%后计算模型基本再生数的改变百分数。结果发现各参数对模型基本再生数的影响从大到小分别为HBV急性患者传播系数β、HBV免疫人群的失效率ψ、HBV疫苗接种率γ3、HBV携带者传播权重ε、HBV急性患者未能清除感染的平均概率q和HBV携带者的康复系数γ2。因此，中国防治HBV传播的最有效途径是增强全社会人群对乙型病毒性肝炎传播知识的教育以减少HBV急性患者传播系数β，以及严格实施易感人群的HBV计划免疫以降低HBV免疫人群的失效率ψ和增加HBV疫苗接种率γ3。
 |
| 图 5 敏感参数优化20%后的模型基本再生数的改变百分数结果 |
根据中国现阶段HBV传播特点，本研究构建分析了HBV传播的动力学模型并对国家卫计委公布的乙肝患者数据进行了拟合分析，结果表明本研究的模型比文献[8]能更好地拟合报道数据，可较好地模拟出自2009年开始的疫情递减趋势。同文献[6]的结论比较，我们发现我国现行HBV防控措施已取得初步成效，HBV感染基本再生数已显著下降。进一步我们发现HBV携带者的基本再生数大于HBV急性患者的基本再生数，说明目前中国乙肝传播的主要危险在于携带者。
虽然本研究发现我国的HBV疫情已得到较好控制，但由于我国HBV患者的基数十分巨大，我们现阶段仍不能对HBV传播掉以轻心，仍需国家从宏观政策方面继续加强调控。结合敏感性分析结果，我们发现进一步增强全社会人群对乙型病毒性肝炎传播知识的教育和严格实施易感人群的HBV计划免疫能更好地控制HBV疫情。
虽然本研究的模型能较好地拟合中国HBV流行动力学数据特点，但为了分析的简便，仍存在一些不足。例如，据调查，尽管我国现在10岁以下的人群HBV携带率已经下降到1%～2%，但20岁以上人群的携带率却还在8%左右，分析其原因，文献[25]认为可能和HBV与HIV的共同感染，即性传播动力学有关，需要在模型中进一步考虑相关因素的影响。此外，由于我国南方人群的HBV携带率高于北方，并且存在各区域经济发展的差异和看病难、看病贵的实际情况，因此年龄和地域分层因素，以及可能存在的未知的HBV感染者等[26, 27]，也可能对我国HBV疫情产生影响，这些都将是我们进一步考虑的重点方向。
| [1] | Dienstag J L. Hepatitis B virus infection. N Engl J Med, 2008, 359(14): 1486-1500. |
| [2] | WHO. Hepatitis B. |
| [3] | 刘耕陶, 斯崇文, 王勤环, 等. 对我国慢性病毒性肝炎防治与研究工作的建议.中华医学杂志, 2002, 82(2): 74-76. |
| [4] | 贾继东, 庄辉. 中国慢性乙型肝炎治疗进展研讨会会议纪要. 中华肝脏病杂志, 2004, 12(11): 698-699. |
| [5] | 中华人民共和国国家卫生和计划生育委员会. 2013国家法定传染病疫情情况.(2014-02-13). |
| [6] | Zou L, Zhang W, Ruan S. Modeling the transmission dynamics and control of hepatitis B virus in China. J Theor Biol, 2010, 262(2): 330-338. |
| [7] | GAVI. The Global Alliance for Vaccines and Immunization. |
| [8] | Hadler S C, Fuqiang C, Averhoff F, et al. The impact of hepatitis B vaccine in China and in the China GAVI Project. Vaccine, 2013, 31(Suppl 9): J66-J72. |
| [9] | Luo Z, Li L, Ruan B. Impact of the implementation of a vaccination strategy on hepatitis B virus infections in China over a 20-year period. Int J Infect Dis, 2012, 16(2): e82-e88. |
| [10] | 崔富强, 王富珍, 吴振华, 等. 中国2005～2010年报告乙型病毒性肝炎发病分析. 中国疫苗和免疫, 2011, 17(6): 483*486, 559. |
| [11] | Edmunds W J, Medley G F, Nokes D J. The transmission dynamics and control of hepatitis B virus in the Gambia. Stat Med, 1996, 15(20): 2215-2233. |
| [12] | Edmunds W J, Medley G F, Nokes D J. Vaccination against hepatitis B virus in highly endemic areas: waning vaccine-induced immunity and the need for booster doses. Trans R Soc Trop Med Hyg, 1996, 90(4): 436-440. |
| [13] | Hahne S, Ramsay M, Balogun K, et al. Incidence and routes of transmission of hepatitis B virus in England and Wales, 1995-2000: implications for immunisation policy. J Clin Virol, 2004, 29(4): 211-220. |
| [14] | Zhang H W, Yin J H, Li Y T, et al. Risk factors for acute hepatitis B and its progression to chronic hepatitis in Shanghai, China. Gut, 2008, 57(12): 1713-1720. |
| [15] | 中华人民共和国国家统计局. 中国统计年鉴2013. 北京: 中国统计出版社, 2013. |
| [16] | 梁晓峰, 陈园生, 王晓军, 等. 中国3岁以上人群乙型肝炎血清流行病学研究. 中华流行病学杂志, 2005, 26(9): 655-658. |
| [17] | 中国疾病控制中心. 中国有1.2亿慢性乙肝病毒携带者 3000万已发病.(2006-04-25) . |
| [18] | Diekmann O, Heesterbeek J A, Metz J A. On the definition and the computation of the basic reproduction ratio R0 in models for infectious diseases in heterogeneous populations. J Math Biol, 1990, 28(4): 365-382. |
| [19] | van-den-Driessche P, Watmough J. Reproduction numbers and sub-threshold endemic equilibria for compartmental models of disease transmission. Math Biosci, 2002, 180: 29-48. |
| [20] | Marino S, Hogue I B, Ray C J, et al. A methodology for performing global uncertainty and sensitivity analysis in systems biology. J Theor Biol, 2008, 254(1): 178-196. |
| [21] | Blower S M, Dowlatabadi H. Sensitivity and uncertainty analysis of complex models of disease transmission: an HIV model, as an example. Int Stat Rev, 1994, 62(2): 229-243. |
| [22] | Rinaldo A, Bertuzzo E, Mari L, et al. Reassessment of the 2010-2011 Haiti cholera outbreak and rainfall-driven multiseason projections. Proc Natl Acad Sci USA, 2012, 109(17): 6602-6607. |
| [23] | Gasemyr J. On an adaptive version of the Metropolis-Hastings algorithm with independent proposal distribution. Scand J Stat, 2003, 30(1): 159-173. |
| [24] | 张世强. 曲线回归的拟合优度指标的探讨. 中国卫生统计, 2002, 19(1): 9-11. |
| [25] | Zou L, Ruan S, Zhang W. On the sexual transmission dynamics of hepatitis B virus in China. J Theor Biol, 2015, 369: 1-12. |
| [26] | Lu F M, Zhuang H. Management of hepatitis B in China. Chin Med J (Engl), 2009, 122(1): 3-4. |
| [27] | 崔富强, 卢永, 王富珍 等. 中国2006年乙型肝炎报告病例构成分析. 中华流行病学杂志, 2007, 28(9): 872-874. |