手足口病是一种常见的传染病，主要由多种肠道病毒引起。感染通常发生于5岁以下、尤其是3岁以下的儿童^[1]。手足口病一年四季均可发生，但夏季和秋季发病率较高^[2]。导致手足口病的肠道病毒有20多种，其中以肠道病毒71型（EV-A71）和柯萨奇病毒A16型（CV-A16）最为常见。中国大陆70%的重症手足口病病例和90%的手足口病相关死亡病例由EV-A71型所引起。

手足口病已成为我国严重的公共卫生问题。自2007年以来，我国手足口病发病率不断上升并伴有严重疫情^[3]，其中上海^[4]和浙江^[5]的疫情也比较严重。鉴于手足口病的广泛传播和潜在严重性，中国和其他国家通过动态模型来解释手足口病传播动态并进行公共卫生干预措施。此前，Chuo等^[6]建立了一个简单的确定性易感-感染-恢复(SIR)模型来预测马来西亚沙捞越发生手足口病期间的感染人数。Yang等^[7]使用SEIR模型，根据手足口病的潜伏期增加一个仓室表示暴露人群(E)。通过在SEIR模型中加入隔离室(Q)，建立SEIQR模型模拟手足口病传播，发现隔离是手足口病防护的最佳策略^[8]。但是这些模型未考虑接种疫苗的因素，目前也没有针对上海市和浙江省手足口病动力学模型的研究。

手足口病临床上无特效药物，主要防治措施为隔离传染源和切断传播途径。预防EV-A71所致手足口病最好的方法就是接种疫苗。目前中国大陆的三家生产商已完成EV-A71疫苗的三期临床试验，并于2015年12月获得了中国食品药品监督管理局的批准。共有3种单价EV-A71疫苗获得了许可，而两价EV-A71/CV-A16疫苗也在开发中。目前还需要更多的研究来评估这些不同疫苗在中国和其他地区预防手足口病暴发的安全性和有效性。

本研究通过建立SEIR动力学模型，了解手足口病在上海市和浙江省的传播情况，并预测不同接种方案下手足口病的负担变化。利用现有的EV-A71疫苗，模拟不同的覆盖率场景下对手足口病发病率产生的不同影响。此外我们还进行了敏感性分析，以发现哪些关键因素可能与手足口病的传播有关。

资料和方法

数据来源   2015—2017年上海市和浙江省每月手足口病例数来自上海市^[9]和浙江省^[10]疾病预防控制中心传染病报告系统，如表 1所示；人口学数据包括5岁以下儿童人数、出生率和5岁以下儿童死亡率, 来自国家统计局^[11]；传播率、疾病持续时间、血清型分布和疫苗效力等参数来自于以往文献报道^[12-13]。

模型描述   根据手足口病的自然病史和接种疫苗的干预，将人群分为8个仓室，建立了扩展SEIR模型来评估上海市和浙江省手足口病的传播和发病率。总人数（N）分为4类：手足口病易感人群（S）、手足口病的相关病原体暴露人群（E），感染人群（I），康复人群（R）。除了康复人群这一类别之外，其余3类按照疫苗接种、病原体类型和住院情况再分为2~3个组别：接种人群（V）、未接种疫苗的暴露人群（E₁）、接种EV-A71疫苗的暴露人群（E₂），未去医院就诊的感染人数（I₁），感染非EV-A71所致的手足口病（I₂）和感染EV-A71所致的手足口病（I₃）。因此，总人数用N=S+V+E₁+E₂+I₁+I₂+I₃+R表示。模型中涉及的参数和描述如图 1和表 2所示。

SEIR模型传播动力学微分方程描述如下:

$ \left\{\begin{array}{l} \frac{d \mathrm{S}}{d \mathrm{t}}=\pi N+\gamma R-\frac{\beta S \sum I}{N}-v-\mu S \\ \frac{d V}{d t}=v-\frac{\beta V \sum I}{N}-\mu V \\ \frac{d E_{1}}{d t}=\frac{\beta S \sum I}{N}-\alpha P_{j} E_{1}-\mu E_{1}(\text { where } j=1,2,3) \\ \frac{d E_{2}}{d t}=\frac{\beta V \sum I}{N}-\alpha \times\left[P_{1}+P_{2}+(1-V E) \times P_{3}\right] E_{2}-\mu E_{2} \\ \frac{d I_{1}}{d t}=\alpha P_{1} E_{i}-\gamma_{1} I_{1}-\mu I_{1}(\text { where } i=1,2) \\ \frac{d I_{2}}{d t}=\alpha P_{2} E_{i}-\gamma_{2} I_{2}-\mu I_{2} \\ \frac{d I_{3}}{d t}=\alpha P_{3} E_{1}+(1-V E) \alpha P_{3} E_{2}-\gamma_{3} I_{3}-\mu I_{3} \\ \left.\frac{d R}{d t}=\gamma_{1} I_{N}+\gamma_{2} I_{2}+\gamma_{3} I_{3}-\gamma R-\mu R\right] \end{array}\right. $

(1)

模型假设   该模型基于以下4个假设:

(1) 5岁以下儿童数初值为S(0)，所有新生儿均为易感人群。将1月份报告的病例数除以30，得到I₂(0)+I₃(0)的初值。

(2) 假设接种疫苗的人在接种EV-A71疫苗后立即受到保护，不存在对其他病毒株的交叉保护。未考虑疫苗接种后手足口病原体变化的可能性。

(3) 据中国CDC手足口病报告可知手足口病死亡人数极少，约为万分之一。因此, 在本模型中没有考虑手足口病死亡率，并将死亡率设置为μ=0。

(4) 自从EV-A71疫苗于2016年12月在中国推出以来，疫苗的使用量仍然非常低。我们假设没有明显影响疾病的传播率β, 在2016年之前和研究期间每个人都易感。

参数估计  使用MATLAB R2018b (version 9.5)中的fminsearch工具^[12]计算最小平方和^[16](minimun sum of square, MSS)来估计参数，该工具用于求解无约束非线性规划问题的最小值。只有当fminsearch的结果收敛后，才能得到所有的最优参数值。

$ \begin{aligned} M S S=& \sum\left( log _{2}(1+\text { data per month })-log _{2}\right.\\ &\left(1+I_{2} \text { and } I_{3} \text { on day } 1+\right.\text { on day } 2+\ldots+\\ &\text { on day } 28, \text { or } 29,30,31))^{2} \end{aligned} $

(2)

其中I₂表示非EV-A71所致手足口病的患者数量，I₃表示EV-A71手足口病的患者数量。

基本再生数（R₀）的推导   基本再生数（R₀）在流行病学研究中具有重要意义。它是在引入一个典型的感染个体后，在完全易感人群中产生的继发性病例的预期数量。本研究中，我们使用van den Driessche方法计算了2015年和2016年的基本再生数（R₀）^[17]。R₀与传输速率（β），进展速度（α）的比例I₁，I₂，I₃（P₁，P₂，P₃）, 康复率（γ₁，γ₂，γ₃）和死亡率（μ）有关，方程表达式如下：

$ \begin{aligned} R_{0}=& \frac{ \mathit{β}}{ \mathit{α}+\mu}+\frac{\beta \alpha P_{1}}{( \mathit{α}+\mu)\left(\gamma_{1}+\mu\right)}+\\ & \frac{0.3 \times \beta \alpha\left(1-P_{1}\right)}{(\alpha+\mu)\left(\gamma_{2}+\mu\right)}+\frac{0.7 \times \beta \alpha\left(1-P_{1}\right)}{(\alpha+\mu)\left(\gamma_{3}+\mu\right)} \end{aligned} $

(3)

结果

模型拟合   将拟合好的参数取值代入模型，计算出2015到2017年上海市和浙江省手足口病发病数拟合值。结果如图 2所示，手足口病发病呈现明显的季节性，2015—2016年上海市和2015—2017年浙江省模型拟合数据与真实数据近似，一年内都存在两个发病高峰，分别在6~7月和11~12月出现；但2017年上海市发病高峰在10~11月，小高峰在6~7月。2016年发病人数及发病高峰人数高于2015年，2017年的发病人数及发病高峰人数呈现明显下降趋势。通过模型方程式和van den Driessche的计算方法，推导出本次研究模型的R₀值，2015—2017年上海市手足口病基本再生数R₀的值分别为1.578、1.602和1.931；2015—2017年浙江省手足口病基本再生数R₀的值分别为1.599、1.724和1.657，R₀均大于1，说明手足口病将会在人群中流行甚至暴发。

为了检验模型拟合结果与现实报告病例数之间的拟合优度，我们使用χ²检验并考虑以下假设。H0：拟合结果（HFMD病例数）等于表 1中所示的报告数据；H1：拟合结果不等于表 1中所示的报告数据。χ²值和P值显示在表 3中。所有P值均大于0.05，因此，不能拒绝5%显著性水平的原假设，拟合结果是对真实值的良好估计。

不同疫苗覆盖场景下上海市和浙江省疾病负担预测   为了解不同疫苗覆盖场景下手足口发病率的差异，引入如下疫苗接种比例：（1）0%，全年不接种EV-A71疫苗；（2）10% EV-A71疫苗覆盖率；（3）30% EV-A71疫苗覆盖率；（4）50% EV-A71疫苗覆盖率；（5）80% EV-A71疫苗覆盖率；（6）100% EV-A71疫苗覆盖率。使用模型来模拟不同疫苗覆盖情景下上海市和浙江省EV-A71手足口病的发病率。结果如图 3所示，在上海市和浙江省随着疫苗接种率的上升，EV-A71手足口病的发病率呈现明显下降趋势，接种率达到80%，EV-A71所致手足口病病例数将减少1/3左右，发病仍然存在季节性趋势，发病高峰在6~7月。

敏感性分析   考虑模型中对R₀值产生影响的参数，R₀的表达式中与疾病属性相关的参数β、α、γ1、γ2、γ3可看作模型中的独立变量，而其他与人口统计学数值相关的参数可看作为常量。估计5个参数在拟合值上下5%区间内的变化对EV-A71手足口病病例数的影响，并根据敏感性分析得出影响范围，结果如图 4和5所示。敏感性分析曲线表明，上海市和浙江省的疾病传播率（β）、未去医院就诊比例（P1）、未去医院就诊的康复率（γ1）、非EV-A71感染患者的康复率（γ2）、EV-A71感染患者康复率（γ3）对EV-A71手足口病的发病率影响方向和程度存在不同。随着参数β的增加，手足口病例数也随之增加，其他4个参数呈现相反趋势。同时，可以清楚地看到P₁、β和γ1对病例数的影响最大，γ2和γ3对病例数的影响相对较轻微。

讨论

目前有很多研究从公共卫生或统计模型的角度来控制和预防手足口病，然而基于数学模型来模拟手足口病的发病动态，预测不同疫苗覆盖情况下手足口病疾病负担的文章却很少。本研究根据手足口病流行特点，考虑疾病的潜伏期和EV-A71疫苗接种情况，利用上海市和浙江省的手足口病发病例数构建了含有8个参数的SEIR模型。2015—2017年上海市和的浙江省模型拟合数据与真实数据近似，并且发现手足口病发病呈现明显的季节性，这与国内其他地区^[18]的发病情况和模型拟合结果一致。

R₀是疾病暴发的重要指标。它的值表示一定时间内一名感染者在易感人群中传染而发生新感染的人数。如果R₀大于1，即平均每一个处于传染期的感染个体在易感人群中传染大于1人，疾病就会在人群中流行甚至暴发。如果R₀小于1，平均每一个处于传染期的感染个体在易感人群中传染少于1人，传染规模不会扩大，一般认为疾病不会流行，感染就会消失。在我们的研究中，R₀的公式由4部分组成。这与中国台湾学者^[19]的一项研究相一致。在本研究中，2015—2017年上海市和浙江省的R₀值都大于1。从理论上讲，如果没有其他干预措施，手足口病在人群中仍会流行甚至有暴发的风险。

除了建立SEIR模型来了解报道的手足口病数据外，本研究还引入了疫苗接种来预测手足口病发病率如何随着疫苗覆盖率的不同而变化。结果显示在上海市和浙江省，随着疫苗接种率的上升，EV-A71手足口病的发病率呈现明显下降趋势，发病仍然存在季节性趋势，发病高峰在6~7月。因此可以看出EV-A71疫苗接种对于控制EV-A71所致手足口病来说有重要影响。虽然一般认为EV-A71疫苗可以预防EV-A71，但不能预防CV-A16感染，其他血清型别所致的手足口病可能无法因接种EV-A71疫苗而得到很好控制。

在本研究中，我们还进行了敏感性分析，研究了不同参数对手足口病发病率的影响。发现这5个参数对EV-A71手足口病的发病率影响方向和程度存在不同。随着参数β的增加，手足口病例数也随之增加，其他4个参数呈现相反趋势。同时，可以清楚地看到P₁、β和γ₁对病例数的影响最大，γ₂和γ₃对病例数的影响相对较轻微。从生物学意义角度来考虑，β也称为传输速率，由两部分组成：接触频率×感染效率，感染效率是由不同病原体自身的特性决定的，很难进行干预，因此可以通过其他方式改变接触频率的值。γ₁表示感染但没有住院的人的康复率，等于恢复时间的倒数。由于大多数手足口病是自限性疾病，感染但未住院的患者占所有手足口病患者的绝大多数，因此家庭护理和抗病毒药物自购需求是影响疾病传播的关键因素。总的来说，可以从以下几个方面来控制手足口病：包括在传播的早期隔离患者，更规律地洗手和良好通风^[20]，以及为患者提供更好的治疗条件，从而增加γ₁、γ₂和γ₃的值。

本研究尚存在以下不足：第一，由于手足口病可由多种血清型肠病毒引起，现有的EV-A71疫苗对其他肠病毒亚型没有交叉保护，虽然EV-A71疫苗为控制由EV-A71引起的手足口病提供了更具体、更有效的方法，但其对手足口病整体流行的预防作用非常有限。第二，研究所采用的手足口病疫情资料按月统计，无法获得每日的病例数，这可能会对模型的准确性产生一定影响。第三，EV-A71疫苗于2016年底在中国推出，但是在调查中，我们发现人们对这种疫苗的接受程度普遍较低，因为它的有效性有限且价格昂贵。因此，本研究中疫苗覆盖率只是基于我们的假设，考虑到成本效益，可能在未来无法实现较高的覆盖率。