具有有限医疗资源的手足口病模型分析

Analysis of Hand-Foot-Mouth Disease Models with Limited Medical Resources

摘要: 本文建立了一个具有有限医疗资源的手足口病传播模型，首先计算了模型的基本再生数，讨论了地方病平衡点的存在性，并分析了无病平衡点和地方病平衡点的稳定性，如果假定在治疗能力之内时治疗函数满足一个非线性函数，当染病者数量超过了治疗能力承受的界限时，治疗函数为某一常数形式。模拟发现，若治疗函数较小时，会出现后向分支，这表明基本再生数小于1不一定能够使疾病消除。最后结合2023年中国大陆手足口病月报告病例数据进行拟合，通过数值模拟来展示有限医疗资源对手足口病传播的影响。

Abstract: In this paper, we develop a model for the spread of HFMD with limited medical resources, first calculate the basic regeneration number of the model, discuss the existence of the endemic equilibrium point, and analyze the stability of the disease-free equilibrium point and the endemic equilibrium point, if it is assumed that the treatment function satisfies a nonlinear function when it is within the treatment capacity, and when the number of infected people exceeds the limit of the treatment capacity to withstand the treatment, the treatment function is a certain constant form. The simulation finds that backward branching occurs if the treatment function is small, which suggests that a basic regeneration number less than 1 does not necessarily lead to disease elimination. Finally, numerical simulations are fitted with the data of monthly reported cases of HFMD in mainland China in 2023 to demonstrate the impact of limited medical resources on the spread of HFMD.

1. 引言

手足口病(Hand, Foot, and Mouth Disease, HFMD)是由多种肠道病毒引起的儿童传染病[1]，又名发疹性水疱性口腔炎。常见的病原体包括柯萨奇病毒A16型(CoX A16)和肠道病毒71型(EV71)，其中，EV71型病毒具有较强的毒性，容易引发神经系统并发症，是造成手足口病的主要病原体[2]。通常表现为口腔疼痛、低烧，以及手、足、口腔等部位出现疱疹或溃疡，严重时可致命。该病具有高度传染性，主要通过密切接触传播，途径包括粪–口或呼吸道。

目前，手足口病仍然对儿童的身心健康造成威胁。2003年，中国大规模爆发了SARS传染病，中国卫生部开始每月公布国家法定传染病报告病例，2008年3月安徽阜阳爆发EV71型疫情，累计报告超9000例，并导致20余名儿童死亡，这一重大公共卫生事件直接促使国家将其纳入丙类传染病管理，当年全国报告病例数达49万，从中国统计年鉴查到2004到2008年，0至5岁以下儿童有8000万人，发病率达6.13‰，到了2014年，患病人数达277.9万例，这个病例数远远高于2008年患病人数(见图1和图2)，因此手足口病已严重影响5岁以下儿童身心健康。幸运的是，2016年中国研制出针对手足口病EV71型病毒引发的疾病疫苗，投入市场后手足口病患病人数没有明显下降，造成这一原因是目前使用的疫苗仅针对肠道病毒71型，而对柯萨奇病毒A16型等其他病原体缺乏有效防护。此外，不同地区的医疗资源配置不平衡(如临床医生、护士、病床数，先进设备、药品以及能否得到及时治疗等)。同时，公众对手足口病的认知水平参差不齐，特别是在传播途径和预防措施方面存在显著差异，进一步加剧了疾病的传播风险，且在预防和控制疾病上没有做出相应的措施。此外，感染手足口病之后，多数患者在1~2周内可自愈。然而，手足口病在我国呈现明显的季节性高峰(5~7月为主高峰，9~11月为次高峰)，尤其对于EV71型病毒易引发重症的病毒株感染患者，如果延迟治疗(如未能及时就医或忽视早期症状)可能导致严重后果，容易导致疾病持续存在。

关于手足口病模型的研究目前非常多，Ma等人[3]基于山东省手足口病的病例数据，提出了一个具有季节性传播的SEIIeQR模型，通过对一些关键参数的敏感性分析，发现隐性感染群体在疾病传播中起着重要作用，隔离人群对疾病的控制作用较小，即使隔离率值很大，基本再生数仍大于1，因此只隔离感染个体并不能有效控制手足口病的传播。Wang和Xiao [4]提出了一个SEIIeRW的数学模型来研究手足口病在中国大陆的传播动态。数值模拟表明，提高病毒清除率和降低传播率可以延缓疫情爆发，减轻手足口病的严重程度，敏感性分析表明，基本再生数对无症状感染者引起的传播率以及与污染环境相关的参数(如间接传播率、清除率和病毒脱落率)敏感。Wang等人[5]将人群划分为7岁以下儿童和成人两个群体，通过拟合中国大陆2011年1月至2018年12月的手足口病数据，得出基本再生数为0.9599，敏感性分析表明降低7岁以下儿童群体的感染风险，同时加强对高龄成人群体的隔离管理，能够合理的控制手足口病的传播。

在手足口病等突发性疫情中，当感染人数激增时，医院的接受效率可能因资源饱和而显著下降，为更真实地反映医疗资源受限时的治疗能力，在构建模型时根据不同的情况采用不同的治疗函数，2004年Wang和Ruan [6]提出了一个恒定恢复率函数来分析传染病，表达式如下

$T_1=\{\begin{array}{l}0,I=0,\\ k,I>0.\end{array}$

其中，$k>0$ 为常数表示感染者的治疗能力，$I$ 表示感染个体的数量，说明为了让疾病完全消失，可以采取最大治疗能力，但这种恒定恢复率只适用于感染人数较少和医疗资源充足时。随后，Wang [7]又提出了如下分段线性治疗函数

$T_2=\{\begin{array}{l}rI,0\le I\le I_0,\\ k,I>I_0.\end{array}$

其中，$k=r{I}_0$ ，这里$I_0$ 代表最大治疗能力，即在医疗资源达到饱和之前，治疗率与感染个体数量成正比，当医疗资源达到饱和之后，治疗率保持在$k$ ，然而上述的治疗函数并不是处处可微的，无法解释当感染个体数量趋于无穷时，治疗率无限接近某一值。因此，Zhang和Liu [8]提出了一个饱和治疗函数为

$h\left(I\right)=\frac{cI}{1+bI}.$

其中，$b\left(\ge 0\right)$ 表示衡量被感染者延迟治疗效果的饱和因子，$c$ 表示治疗率，当$I$ 很小时，$h\left(I\right)$ 可以用$cI$ 来表示，当$I$ 较大时，治疗率达到$c/b$ ，当$b=0$ 时，$h\left(I\right)$ 变成线性治疗率，所以函数$h\left(I\right)$ 表示有限的医疗资源是可行的，更符合实际情况。

本文建立一个具有有限医疗资源的手足口病模型来研究中国大陆手足口病模型传播。首先计算了模型的基本再生数，讨论了地方病平衡点的存在性，并分析了无病平衡点的稳定性和地方病平衡点的稳定性。此外还分析了模型出现的后向分支。最后利用中国大陆2023年1月至12月的手足口病月数据进行拟合，通过数值模拟来展示不同控制措施对手足口病传播的影响。

Figure 1. Number of reported cases of hand, foot and mouth disease from 2008 to 2024

图1. 2008~2024年报道的手足口病病例数

Figure 2. Reported deaths from hand, foot, and mouth disease from 2008 to 2024

图2. 2008~2024年报道的手足口病死亡数

2. 模型建立

本文将考虑手足口病5岁以下儿童，将人群分为：易感者(S)，潜伏者(E)，显性感染者(I)，隐性感染者(I_e)，恢复者(R)，在模型中，易感者通过与显性感染者和隐性感染者接触后先进入潜伏期，在潜伏期之后表现出症状的为显性感染者，不表现出任何症状的为隐性感染者，其中选用饱和治疗函数来衡量有限医疗资源，并考虑了疫苗接种对手足口病的影响，建立的模型由以下微分方程表示。

$\{\begin{array}{l}\frac{\text{d}S}{\text{d}t}=\Lambda -{\beta }_{1}SI-{\beta }_{2}S{I}_e-dS-vS，\\ \frac{\text{d}E}{\text{d}t}={\beta }_{1}SI+{\beta }_{2}S{I}_e-\eta E-dE，\\ \frac{\text{d}I}{\text{d}t}=\rho \eta E-dI-\delta I-{\gamma }_{1}I-\frac{cI}{1+bI},\\ \frac{\text{d}{I}_e}{\text{d}t}=\left(1-\rho \right)\eta E-d{I}_e-{\gamma }_2{I}_e,\\ \frac{\text{d}R}{\text{d}t}=\frac{cI}{1+bI}+{\gamma }_{1}I+{\gamma }_2{I}_e+vS-dR.\end{array}$ (1)

其中所有参数值均为正值，$\Lambda$ 为人口出生率，$d$ 为自然死亡率，${\beta }_1$ 为易感者与显性感染者接触的传染率，${\beta }_2$ 为易感者与隐性感染者接触的传染率，$v$ 为疫苗接种率，$\rho$ 为显性感染者的比例，$1/\eta$ 为疾病平均潜伏期，${\gamma }_1$ 为显性感染者的恢复率，${\gamma }_2$ 为隐性感染者的恢复率，$\delta$ 为显性感染者的因病死亡率，$c$ 为治疗率，$b$ 为衡量延迟治疗对感染个体的影响因子。

系统(1)满足初始条件：$S\left(0\right)\ge 0,E\left(0\right)\ge 0,I\left(0\right)\ge 0,I_e\left(0\right)\ge 0,R\left(0\right)\ge 0$ ，设$N\left(t\right)$ 为$t$ 时刻种群中个体的总数，有$N\left(t\right)=S\left(t\right)+E\left(t\right)+I\left(t\right)+I_e\left(t\right)+R\left(t\right)$ ，通过系统(1)，可以得到

$\begin{array}{c}\frac{\text{d}N\left(t\right)}{\text{d}t}=\Lambda -dS\left(t\right)-dE\left(t\right)-dI\left(t\right)-d{I}_e\left(t\right)-dR\left(t\right)-\delta I\\ =\Lambda -dN\left(t\right)-\delta I\\ \le \Lambda -dN\left(t\right)\end{array}$

计算得

$\underset{t\to +\infty }{\lim }\sup \left(N\left(t\right)\right)\le \frac{\Lambda }{d}$ ，

则

$X=\left\{\left(S,E,I,I_e,R\right)\in {\Re }_{+}^5:S+E+I+I_e+R\le \frac{\Lambda }{d}\right\}$

是系统(1)的正不变集，由于前四个方程不包含变量$R$ ，因此只需研究以下模型

$\{\begin{array}{l}\frac{\text{d}S}{\text{d}t}=\Lambda -{\beta }_{1}SI-{\beta }_{2}S{I}_e-dS-vS,\\ \frac{\text{d}E}{\text{d}t}={\beta }_{1}SI+{\beta }_{2}S{I}_e-\eta E-dE,\\ \frac{\text{d}I}{\text{d}t}=\rho \eta E-dI-\delta I-{\gamma }_{1}I-\frac{cI}{1+bI},\\ \frac{\text{d}I{}_e}{\text{d}t}=\left(1-\rho \right)\eta E-d{I}_e-{\gamma }_2{I}_e.\end{array}$ (2)

通过求解可知，系统(2)存在无病平衡点$P_0=\left(S_0,0,0,0\right)$ ，其中$S_0=\Lambda /\left(d+v\right)$ ，下面通过Driessche and Watmough的下一代矩阵演绎法[9]，计算系统(2)的基本再生数为

$\begin{array}{c}{R}_0=\frac{S_0{\beta }_1\eta \rho }{\left(\eta +d\right)\left(d+\delta +{\gamma }_1+c\right)}+\frac{S_0{\beta }_1\eta \left(1-\rho \right)}{\left(\eta +d\right)\left(d+{\gamma }_2\right)}\\ =R_1+R_2,\end{array}$

其中$R_1$ 表示显性感染者对易感个体的传播能力，$R_2$ 表示隐性感染者对易感个体的传播能力。

3. 地方病平衡点的存在性

本节中，将讨论系统(2)的地方病平衡点的存在性，令系统(2)的左端等于0，得到

$\{\begin{array}{l}\Lambda -{\beta }_{1}SI-{\beta }_{2}S{I}_e-dS-vS=0,\\ {\beta }_{1}SI+{\beta }_{2}S{I}_e-\eta E-dE=0,\\ \rho \eta E-dI-\delta I-{\gamma }_{1}I-\frac{cI}{1+bI}=0,\\ \left(1-\rho \right)\eta E-d{I}_e-{\gamma }_2{I}_e=0.\end{array}$ (3)

由方程(3)，通过计算可得到关系式

$S^{*}=\frac{\Lambda }{{\beta }_1{I}^{*}+{\beta }_2{I}_e^{*}+M_3},E^{*}=\frac{M_1+M_{4}b{I}^{*}}{\rho \eta \left(1+b{I}^{*}\right)}{I}^{*},I_e^{*}=\frac{M_1+M_{4}b{I}^{*}}{1+b{I}^{*}}A{I}^{*}.$

其中，

$M_1=d+\delta +{\gamma }_1+c,M_2=d+{\gamma }_2,M_3=d+v,$

$M_4=d+\delta +{\gamma }_1,A=\frac{1-\rho }{M_2\rho }.$

将上述$S^{*},E^{*},I_e^{*}$ 表达式代入系统(2)的第二个等式中，得到一个关于$I^{*}$ 的一元三次方程

$f\left(I\right)=a_0{I}^3+a_1{I}^2+a_{2}I+a_3=0$ ，

其中$a_0,a_1,a_2,a_3,B$ 分别为

$a_0={\beta }_{1}B{M}_4{b}^2+{\beta }_{2}BA{M}_4^2{b}^2>0$ ，

$a_1={\beta }_{1}B{M}_{1}b+{\beta }_{2}BA{M}_1{M}_{4}b+{\beta }_{1}B{M}_{4}b+{\beta }_{2}BA{M}_1{M}_{4}b+B{M}_3{M}_4{b}^2-{\beta }_1\Lambda b^2-{\beta }_2\Lambda A{M}_4{b}^2$ ，

$a_2={\beta }_{1}B{M}_1+{\beta }_{2}B{M}_1^{2}A+B{M}_1{M}_{3}b+B{M}_3{M}_{4}b-2{\beta }_1\Lambda b-{\beta }_2\Lambda A{M}_{1}b-{\beta }_2\Lambda A{M}_{4}b$ ，

$a_3=B{M}_1{M}_3-{\beta }_1\Lambda -{\beta }_2\Lambda A{M}_1=B{M}_1{M}_3\left(1-R_0\right)$ ，

$B=\frac{\eta +d}{\rho \eta }$ .

令

$A_0=a_1^2-3a_0{a}_2,B_0=a_1{a}_2-9a_0{a}_3,C_0=a_2^2-3a_1{a}_3.$

定义${\Delta }_0$ 是$f\left(I\right)$ 关于$I$ 的判别式，则

${\Delta }_0=B_0^2-4A_{0}C$ ，

$f\left(I\right)$ 关于$I$ 的求导为

$f^{\prime }\left(I\right)=3a_0{I}^2+2a_{1}I+a_2$ ，

定义$\Delta$ 为$f^{\prime }\left(I\right)=0$ 关于$I$ 的判别式，其中$\Delta ={\left(2a_1\right)}^2-4\left(3a_0\right)a_2$ 。

注意到$f\left(0\right)=B{M}_1{M}_3\left(1-R_0\right)$ ，$f^{\prime }\left(0\right)=a_2$ ，则$f\left(I\right)$ 的正根有以下几种情况：

一、$R_0>1$ ，在这种情况下，$a_3<0$ ，

1. $a_1>0$ ，无论$a_2$ 为正还是负，$f\left(I\right)=0$ 均只有一个正根。

2. $a_1<0$ ，

(1) $a_2>0$ ，如果$\Delta \le 0$ ，方程$f\left(I\right)=0$ 只有一个正根；如果$\Delta >0$ ，方程$f\left(I\right)=0$ 有三个正根。

(2) $a_2<0$ ，方程$f\left(I\right)=0$ 只有一个正根。

二、$R_0<1$ ，在这种情况下，$a_3>0$ ，

1. $a_1>0$ ，

(1) $a_2>0$ ，方程$f\left(I\right)=0$ 没有正根，

(2) $a_2<0$ ，类似于Huang和Ruan [10]在引理2.1中所描述的，如果${\Delta }_0>0$ ，方程$f\left(I\right)=0$ 没有正根；如果${\Delta }_0=0$ ，方程$f\left(I\right)=0$ 存在唯一的正根；如果${\Delta }_0<0$ ，方程$f\left(I\right)=0$ 有两个正根。

2. $a_1<0$ ，

(1) $a_2>0$ ，如果$\Delta \le 0$ ，方程$f\left(I\right)=0$ 没有正根；如果$\Delta >0$ ，方程$f\left(I\right)=0$ 有两个正根。

(2) $a_2<0$ ，同情况二、1.(2)。

三、$R_0=1$ ，方程$f\left(I\right)=0$ 变为$f\left(I\right)=a_0{I}^3+a_1{I}^2+a_{2}I=0$ ，

1. $a_1>0$ ，

(1) $a_2>0$ ，方程$f\left(I\right)=0$ 没有正根，

(2) $a_2<0$ ，方程$f\left(I\right)=0$ 只有一个正根。

2. $a_1<0$ ，

(1) $a_2>0$ ，如果${\Delta }_0>0$ ，方程$f\left(I\right)=0$ 没有正根；如果${\Delta }_0=0$ ，方程$f\left(I\right)=0$ 存在唯一的正根；如果${\Delta }_0<0$ ，方程$f\left(I\right)=0$ 有两个正根。

(2) $a_2<0$ ，同情况三、1.(2)。

定理1 系统(2)的地方病平衡点的存在性有以下结果成立。

1. 当$R_0>1$ ，

(1) $a_1>0$ ，系统(2)存在唯一的地方病平衡点；

(2) $a_1<0$ ，$a_2>0$ ，$\Delta \le 0$ ，系统(2)存在唯一的地方病平衡点；

(3) $a_1<0$ ，$a_2>0$ ，$\Delta >0$ ，系统(2)存在三个地方病平衡点；

(4) $a_1<0$ ，$a_2<0$ ，系统(2)存在唯一的地方病平衡点。

2. 当$R_0<1$ ，

(1) $a_1>0$ ，$a_2<0$ ，${\Delta }_0=0$ ，系统(2)存在唯一的地方病平衡点；

(2) $a_1>0$ ，$a_2<0$ ，${\Delta }_0<0$ ，系统(2)存在两个地方病平衡点；

(3) $a_1<0$ ，$a_2>0$ ，$\Delta >0$ ，系统(2)存在两个地方病平衡点；

(4) $a_1<0$ ，$a_2<0$ ，${\Delta }_0=0$ ，系统(2)存在唯一的地方病平衡点；

(5) $a_1<0$ ，$a_2<0$ ，${\Delta }_0<0$ ，系统(2)存在两个地方病平衡点。

3. 当$R_0=1$ ，

(1) $a_1>0$ ，$a_2<0$ ，系统(2)存在唯一的地方病平衡点；

(2) $a_1<0$ ，$a_2>0$ ，${\Delta }_0=0$ ，系统(2)存在唯一的地方病平衡点；

(3) $a_1<0$ ，$a_2>0$ ，${\Delta }_0<0$ ，系统(2)存在两个地方病平衡点；

(4) $a_1<0$ ，$a_2<0$ ，系统(2)存在唯一的地方病平衡点。

4. 平衡点的稳定性

4.1. 无病平衡点的局部稳定性

定理2 对于系统(2)，如果$R_0<1$ ，那么无病平衡点$P_0=\left(S_0,0,0,0\right)$ 是局部渐近稳定的。

证明 系统(2)在$P_0$ 处的雅可比矩阵是

$J\left(P_0\right)=\left(\begin{array}{cccc}-\left(d+v\right)& 0& -{\beta }_1{S}_0& -{\beta }_2{S}_0\\ 0& -\left(\eta +d\right)& {\beta }_1{S}_0& {\beta }_2{S}_0\\ 0& \rho \eta & -\left(d+\delta +{\gamma }_1+c\right)& 0\\ 0& \left(1-\rho \right)\eta & 0& -\left(d+{\gamma }_2\right)\end{array}\right),$

计算出一个特征根为$-\left(d+v\right)$ ，其他特征根由以下方程给出

${\lambda }^3+A_1{\lambda }^2+A_2\lambda +A_3=0,$

其中，

$A_1=\left(\eta +d\right)+\left(d+{\gamma }_2\right)+\left(d+\delta +{\gamma }_1+c\right)$ ，

$A_2=\left(d+\delta +{\gamma }_1+c\right)\left(d+{\gamma }_2\right)+\left(\eta +d\right)\left(d+{\gamma }_2\right)+\left(\eta +d\right)\left(d+\delta +{\gamma }_1+c\right)-S_0\left[{\beta }_1\eta \rho +{\beta }_2\eta \left(1-\rho \right)\right]$ ，

$\begin{array}{c}{A}_3=\left(\eta +d\right)\left(d+{\gamma }_2\right)\left(d+\delta +{\gamma }_1+c\right)-{\beta }_1{S}_0\eta \rho \left(d+{\gamma }_2\right)-{\beta }_2{S}_0\eta \left(1-\rho \right)\left(d+\delta +{\gamma }_1+c\right)\\ =\left(\eta +d\right)\left(d+{\gamma }_2\right)\left(d+\delta +{\gamma }_1+c\right)\left(1-R_0\right).\end{array}$

注意到$A_1>0$ ，若$R_0<1$ ，则$A_3>0$ ，并进一步得到

$\begin{array}{c}{A}_1{A}_2-A_3=\left(\eta +d\right)\left(d+\delta +{\gamma }_1+c\right)\left(\eta +2d+\delta +{\gamma }_1+c\right)\left(1-R_1\right)\\ +\left(\eta +d\right)\left(d+{\gamma }_2\right)\left(\eta +2d+{\gamma }_2\right)\left(1-R_2\right)\\ +\left(d+\delta +{\gamma }_1+c\right)\left(d+{\gamma }_2\right)\left(2\eta +4d+{\gamma }_2+\delta +{\gamma }_1+c\right).\end{array}$

当$R_0<1$ 时，那么$R_1<1$ ，$R_2<1$ ，意味着如果$R_0<1$ ，那么$A_1{A}_2-A_3>0$ ，因此由Routh-Hurwitz准则[11]可知，$P_0\left(\lambda \right)=0$ 的所有根均具有负实部，即当$R_0<1$ 时，无病平衡点$P_0$ 是局部渐近稳定的。

4.2. 无病平衡点的全局稳定性

定理3 如果$R_0<1$ ，系统(2)的无病平衡点$P_0$ 是全局渐近稳定的。

证明 构造Lyapunov函数

$V\left(t\right)=\left(d+{\gamma }_2\right)\left(d+\delta +{\gamma }_1+c\right)E\left(t\right)+\left(d+{\gamma }_2\right){\beta }_1{S}_{0}I\left(t\right)+\left(d+\delta +{\gamma }_1+c\right){\beta }_2{S}_0{I}_e\left(t\right)$ ，

则函数$V\left(t\right)$ 沿系统(2)的全导数为

$\begin{array}{c}\frac{\text{d}V\left(t\right)}{\text{d}t}=\left(d+{\gamma }_2\right)\left(d+\delta +{\gamma }_1+c\right){\beta }_{1}SI+\left(d+{\gamma }_2\right)\left(d+\delta +{\gamma }_1+c\right){\beta }_{2}S{I}_e\\ -\left(d+{\gamma }_2\right)\left(d+\delta +{\gamma }_1+c\right)\left(\eta +d\right)E+\left(d+{\gamma }_2\right){\beta }_1{S}_0\rho \eta E\\ -\left(d+{\gamma }_2\right){\beta }_1{S}_0\left(d+\delta +{\gamma }_1+c/1+bI\right)I\\ +\left(d+\delta +{\gamma }_1+c\right)\left(d+{\gamma }_2\right){\beta }_2{S}_0\left(1-\rho \right)\eta E-\left(d+\delta +{\gamma }_1+c\right)\left(d+{\gamma }_2\right){\beta }_2{S}_0{I}_e\\ \le \left(d+{\gamma }_2\right)\left(d+\delta +{\gamma }_1+c\right){\beta }_1{S}_{0}I+\left(d+{\gamma }_2\right)\left(d+\delta +{\gamma }_1+c\right){\beta }_2{S}_0{I}_e\\ -\left(d+{\gamma }_2\right)\left(d+\delta +{\gamma }_1+c\right)\left(\eta +d\right)E+\left(d+{\gamma }_2\right){\beta }_1{S}_0\rho \eta E\\ -{\beta }_1{S}_0\left(d+{\gamma }_2\right)\left(d+\delta +{\gamma }_1+c\right)I+{\beta }_2{S}_0\left(d+\delta +{\gamma }_1+c\right)\left(1-\rho \right)\eta E\\ -{\beta }_2{S}_0\left(d+\delta +{\gamma }_1+c\right)\left(d+{\gamma }_2\right)I_e\\ ={\beta }_1{S}_0\rho \eta \left(d+{\gamma }_2\right)+{\beta }_2{S}_0\left(d+\delta +{\gamma }_1+c\right)-\left(d+{\gamma }_2\right)\left(d+\delta +{\gamma }_1+c\right)\left(\eta +d\right)E\\ =\left(R_0-1\right)E,\end{array}$

由于所有参数都是非负的，因此当$R_0<1$ 时，$\frac{\text{d}V}{\text{d}t}\le 0$ ，$\frac{\text{d}V}{\text{d}t}=0$ 且仅在$S=S_0$ ，$E=0$ ，$b=0$ 时成立，由LaSalle不变集原理[12]知，当$R_0<1$ ，无病平衡点是全局渐近稳定的。

4.3. 地方病平衡点的局部稳定性

定理4 如果$R_0>1$ ，则定理1中存在唯一的地方病平衡点$P^{*}$ 是局部渐近稳定的。

证明 系统(2)在$P^{*}$ 处的雅可比矩阵是

$J\left(P^{*}\right)=\left(\begin{array}{cccc}-{\beta }_1{I}^{*}-{\beta }_2{I}_e^{*}-\left(d+v\right)& 0& -{\beta }_1{S}^{*}& -{\beta }_2{S}^{*}\\ {\beta }_1{I}^{*}+{\beta }_2{I}_e^{*}& -\left(\eta +d\right)& {\beta }_1{S}^{*}& {\beta }_2{S}^{*}\\ 0& \rho \eta & -\left[d+\delta +{\gamma }_1+\frac{c}{{\left(1+b{I}^{*}\right)}^2}\right]& 0\\ 0& \left(1-\rho \right)\eta & 0& -\left(d+{\gamma }_2\right)\end{array}\right)$ ，

特征方程为

${\lambda }^4+B_1{\lambda }^3+B_2{\lambda }^2+B_3\lambda +B_4=0.$ (4)

令

$b_{11}=-{\beta }_1{I}^{*}-{\beta }_2{I}_e^{*}-\left(d+v\right),b_{13}=-{\beta }_1{S}^{*},b_{14}=-{\beta }_2{S}^{*},b_{21}={\beta }_1{I}^{*}+{\beta }_2{I}_e^{*},b_{22}=-\left(\eta +d\right),$

$b_{23}={\beta }_1{S}^{*},b_{24}={\beta }_2{S}^{*},b_{32}=\rho \eta ,b_{33}=-\left[d+\delta +{\gamma }_1+\frac{c}{{\left(1+b{I}^{*}\right)}^2}\right],b_{42}=\left(1-\rho \right)\eta ,b_{44}=-\left(d+{\gamma }_2\right).$

通过计算得

$B_1=-b_{11}-b_{22}-b_{33}-b_{44},$

$B_2=b_{11}{b}_{22}+b_{11}{b}_{33}+b_{11}{b}_{44}+b_{22}{b}_{33}+b_{22}{b}_{44}+b_{33}{b}_{44}-b_{23}{b}_{32}-b_{24}{b}_{42},$

$\begin{array}{c}{B}_3=b_{11}{b}_{24}{b}_{42}+b_{11}{b}_{23}{b}_{32}+b_{23}{b}_{32}{b}_{44}+b_{24}{b}_{33}{b}_{42}-b_{11}{b}_{22}{b}_{33}-b_{11}{b}_{22}{b}_{44}\\ -b_{11}{b}_{33}{b}_{44}-b_{22}{b}_{33}{b}_{44}-b_{14}{b}_{21}{b}_{42}-b_{13}{b}_{21}{b}_{32},\end{array}$

$B_4=b_{11}{b}_{22}{b}_{33}{b}_{44}+b_{14}{b}_{21}{b}_{33}{b}_{42}+b_{13}{b}_{21}{b}_{32}{b}_{44}-b_{11}{b}_{24}{b}_{33}{b}_{42}-b_{11}{b}_{23}{b}_{32}{b}_{44}.$

由Routh-Hurwitz准则[11]，特征方程(4)的所有特征根均具有负实部 当且仅当$B_1,B_2,B_3,B_4$ ，

$\left|\begin{array}{cc}{B}_1& B_3\\ 1& B_2\end{array}\right|>0$ ，$\left|\begin{array}{ccc}{B}_1& B_3& 0\\ 1& B_2& B_4\\ 0& B_1& B_3\end{array}\right|>0$ ，$\left|\begin{array}{cccc}{B}_1& B_3& 0& 0\\ 1& B_2& B_4& 0\\ 0& B_1& B_3& 0\\ 0& 1& B_2& B_4\end{array}\right|>0$ ，

因此当$R_0>1$ 时，唯一的地方病平衡点$P^{*}$ 为局部渐近稳定。

5. 后向分支

在传染病动力学研究中，后向分支是分支流行病学模型中的一个重要现象，这种现象表明，当基本再生数$R_0$ 降至临界值1以下时，疾病可能仍然持续存在。因此通过利用分支理论[13]中的方法计算发生后向分支的条件。

定理5 如果$b>b^{*}$ 时，系统(2)在$R_0=1$ 处会发生后向分支，其中，

$b^{*}=\frac{d{\left(\eta +d\right)}^2\left(d+\delta +{\gamma }_1+c\right)\left(d+\delta +{\gamma }_1+2c\right)}{{\rho }^2{\eta }^2\Lambda \left(d+v\right)c}.$

证明 令$S=x_1$ ，$E=x_2$ ，$I=x_3$ ，$I_e=x_4$ ，则系统(2)化为

$\{\begin{array}{l}\frac{\text{d}{x}_1}{\text{d}t}=f_1=\Lambda -{\beta }_1{x}_1{x}_3-{\beta }_2{x}_1{x}_4-d{x}_1-v{x}_1,\\ \frac{\text{d}{x}_2}{\text{d}t}=f_2={\beta }_1{x}_1{x}_3+{\beta }_2{x}_1{x}_4-\eta x_2-d{x}_2,\\ \frac{\text{d}{x}_3}{\text{d}t}=f_3=\rho \eta x_2-d{x}_3-\delta x_3-{\gamma }_1{x}_3-\frac{c{x}_3}{1+b{x}_3},\\ \frac{\text{d}{x}_4}{\text{d}t}=f_4=\left(1-\rho \right)\eta x_2-d{x}_4-{\gamma }_2{x}_4.\end{array}$ (5)

当$R_0=1$ 时，选择${\beta }_1$ 为分支参数，因此

${\beta }_1={\beta }_1^{*}=\frac{\left(d+\delta +{\gamma }_1+c\right)\left[\left(\eta +d\right)\left(d+{\gamma }_2\right)-{\beta }_2{S}_0\left(1-\rho \right)\eta \right]}{S_0\rho \eta \left(d+{\gamma }_2\right)}.$

当${\beta }_1={\beta }_1^{*}$ 时，系统(5)在$P_0=\left(S_0,0,0,0\right)$ 处的雅可比矩阵为

$J\left(P_0,{\beta }_1^{*}\right)=\left(\begin{array}{cccc}-\left(d+v\right)& 0& -{\beta }_1^{*}{S}_0& -{\beta }_2{S}_0\\ 0& -\left(\eta +d\right)& {\beta }_1^{*}{S}_0& {\beta }_2{S}_0\\ 0& \rho \eta & -\left(d+\delta +{\gamma }_1+c\right)& 0\\ 0& \left(1-\rho \right)\eta & 0& -\left(d+{\gamma }_2\right)\end{array}\right)$

其中，$\lambda =0$ 是$J\left(P_0,{\beta }_1^{*}\right)$ 的一个特征根，并且其他特征根均具有负实部，$J\left(P_0,{\beta }_1^{*}\right)$ 的零特征根对应的右特征向量设$w={\left(w_1,w_2,w_3,w_4\right)}^{\text{T}}$ ，则

$\left(\begin{array}{cccc}-\left(d+v\right)& 0& -{\beta }_1^{*}{S}_0& -{\beta }_2{S}_0\\ 0& -\left(\eta +d\right)& {\beta }_1^{*}{S}_0& {\beta }_2{S}_0\\ 0& \rho \eta & -\left(d+\delta +{\gamma }_1+c\right)& 0\\ 0& \left(1-\rho \right)\eta & 0& -\left(d+{\gamma }_2\right)\end{array}\right)\left(\begin{array}{c}{w}_1\\ w_2\\ w_3\\ w_4\end{array}\right)=0.$

计算可得

$w={\left(w_1,w_2,w_3,w_4\right)}^{\text{T}}={\left(-\frac{\eta +d}{d+v},1,\frac{\rho \eta }{d+\delta +{\gamma }_1+c},\frac{\left(1-\rho \right)\eta }{d+{\gamma }_2}\right)}^{\text{T}}.$

同理计算$J\left(P_0,{\beta }_1^{*}\right)$ 的零特征根对应的左特征向量，设$u=\left(u_1,u_2,u_3,u_4\right)$ ，且$u\times w=1$ ，因此

$\begin{array}{c}u=\left(u_1,u_2,u_3,u_4\right)\\ =\left(0,\frac{d+\delta +{\gamma }_1+2c}{d+\delta +{\gamma }_1+2c+\rho \eta },\frac{d+\delta +{\gamma }_1+c}{d+\delta +{\gamma }_1+2c+\rho \eta },0\right).\end{array}$

下面求得系统(5)的函数在$\left(P_0,{\beta }_1^{*}\right)$ 处的偏导数为

$\frac{{\partial }^2{f}_1}{\partial x_1\partial x_3}\left(P_0,{\beta }_1^{*}\right)=\frac{{\partial }^2{f}_1}{\partial x_3\partial x_1}\left(P_0,{\beta }_1^{*}\right)=-{\beta }_1^{*},\frac{{\partial }^2{f}_1}{\partial x_1\partial x_4}\left(P_0,{\beta }_1^{*}\right)=\frac{{\partial }^2{f}_1}{\partial x_4\partial x_1}\left(P_0,{\beta }_1^{*}\right)=-{\beta }_2,$

$\frac{{\partial }^2{f}_2}{\partial x_1\partial x_3}\left(P_0,{\beta }_1^{*}\right)=\frac{{\partial }^2{f}_2}{\partial x_3\partial x_1}\left(P_0,{\beta }_1^{*}\right)={\beta }_1^{*},\frac{{\partial }^2{f}_2}{\partial x_1\partial x_4}\left(P_0,{\beta }_1^{*}\right)=\frac{{\partial }^2{f}_2}{\partial x_4\partial x_1}\left(P_0,{\beta }_1^{*}\right)={\beta }_2,$

$\frac{{\partial }^2{f}_1}{\partial x_3\partial {\beta }_1}\left(P_0,{\beta }_1^{*}\right)=-\frac{\Lambda }{d},\frac{{\partial }^2{f}_2}{\partial x_3\partial {\beta }_1}\left(P_0,{\beta }_1^{*}\right)=\frac{\Lambda }{d},\frac{{\partial }^2{f}_3}{\partial x_3^2}\left(P_0,{\beta }_1^{*}\right)=2bc.$

其余偏导数均为0，因此

$\begin{array}{c}a={\displaystyle \sum _{k,i,j}^4{u}_k{w}_i{w}_j\frac{{\partial }^2{f}_k}{\partial x_i\partial x_j}\left(P_0,{\beta }_1^{*}\right)}\\ =2u_1{w}_1{w}_3\frac{{\partial }^2{f}_1}{\partial x_1\partial x_3}\left(P_0,{\beta }_1^{*}\right)+2u_1{w}_1{w}_4\frac{{\partial }^2{f}_1}{\partial x_1\partial x_4}\left(P_0,{\beta }_1^{*}\right)+2u_2{w}_1{w}_3\frac{{\partial }^2{f}_2}{\partial x_1\partial x_3}\left(P_0,{\beta }_1^{*}\right)\\ +2u_2{w}_1{w}_4\frac{{\partial }^2{f}_2}{\partial x_1\partial x_4}\left(P_0,{\beta }_1^{*}\right)+u_3{w}_3{w}_3\frac{{\partial }^2{f}_3}{\partial x_3^2}\left(P_0,{\beta }_1^{*}\right)\\ =2\frac{d+\delta +{\gamma }_1+2c}{\rho \eta +d+\delta +{\gamma }_1+2c}\left(-\frac{\eta +d}{d+v}\right)\frac{\rho \eta }{d+\delta +{\gamma }_1+c}{\beta }_1^{*}+2\frac{d+\delta +{\gamma }_1+2c}{\rho \eta +d+\delta +{\gamma }_1+2c}\left(-\frac{\eta +d}{d+v}\right)\frac{\left(1-\rho \right)\eta }{d+{\gamma }_2}{\beta }_2\\ +\frac{d+\delta +{\gamma }_1+c}{\rho \eta +d+\delta +{\gamma }_1+2c}\frac{\rho \eta }{d+\delta +{\gamma }_1+c}\frac{\rho \eta }{d+\delta +{\gamma }_1+c}2bc,\end{array}$

$\begin{array}{c}b={\displaystyle \sum _{k,i}^4{u}_k{w}_i\frac{{\partial }^2{f}_k}{\partial x_i\partial {\beta }_1}\left(P_0,{\beta }_1^{*}\right)}\\ =u_2{w}_3\frac{{\partial }^2{f}_2}{\partial x_3\partial {\beta }_1}\left(P_0,{\beta }_1^{*}\right)\\ =\frac{d+\delta +{\gamma }_1+2c}{d+\delta +{\gamma }_1+2c+\rho \eta }\frac{\rho \eta }{d+\delta +{\gamma }_1+c}\frac{\Lambda }{d}>0.\end{array}$

根据产生后向分支的条件可知$a>0\iff b>b^{*}$ ，其中，

$b^{*}=\frac{d{\left(\eta +d\right)}^2\left(d+\delta +{\gamma }_1+c\right)\left(d+\delta +{\gamma }_1+2c\right)}{{\rho }^2{\eta }^2\Lambda \left(d+v\right)c}.$ (6)

下面取$\Lambda =120$ ，$d=0.16$ ，${\beta }_2=0.00009$ ，$v=0.6$ ，$\rho =0.1$ ，$\eta =0.7$ ，${\gamma }_1=0.05$ ，${\gamma }_2=0.05$ ，$\delta =0.9$ ，$c=0.07$ ，由$R_0=1$ 可以解出${\beta }_1=0.087264$ ，因此通过式(6)可以求出$b^{*}=5.5798$ 时，取$b=2,6,10$ ，如图3所示。比如2008年安徽阜阳手足口病疫情为后向分支现象提供了典型案例。尽管当地卫生部门迅速采取了控制措施，但由于医疗资源严重不足，疫情仍持续了较长时间。这表明，即使将基本再生数降到1以下，疾病仍然存在，因此要消除疾病，就要减小$b$ 的值，即确保患者得到及时治疗。

Figure 3. Backward branching graph

图3. 后向分支图

6. 数值模拟

本小节从中国疾病预防控制中心官网获取了2023年1月至12月的手足口病月病例数，并将该数据与系统(2)进行拟合。通过查阅中国统计年鉴，2023年出生人口为902万人，平均每月新增人口数大约为$902/12=75.2$ 万人，即$\Lambda =752000$ ；自然死亡率用平均寿命的倒数来表示，即$d$ 约为0.0128；该年手足口病死亡1人，患病总人数为168万人，计算因病死亡率$\delta =6\times {10}^{-7}$ ，手足口病的平均潜伏期约为4天，计算$1/\eta =4/30$ ；患者在感染后平均恢复时间为8.5天，得到恢复率${\gamma }_1$ 为$30/8.5=3.5293$ ，并假设${\gamma }_1={\gamma }_2$ ，初始值$S\left(0\right)$ 可以通过中国统计年鉴来得到；$I\left(0\right)$ 可以由中国疾病预防控制中心来获取，其他参数值及初始值无法直接获得，因此通过使用最小二乘法拟合得到，如表1和表2。系统(1)对2023年中国大陆报告手足口病月病例数拟合程度见图4。

由图5可知，疫苗接种率$v$ 与基本再生数$R_0$ 之间存在显著的负相关关系，为更直观展示疫苗接种对手足口病的效果，我们进一步模拟了不同接种率下各仓室的变化情况，如图6和图7，疫苗接种率的提升显著降低了易感者和感染者的数量。因此在手足口病高峰时期，应当加大对儿童的疫苗接种。

当固定其他值，选取接种率$v=0.01$ ，此时计算基本再生数为$R_0=0.71738<1$ ，得到了五个仓室随时间的变化图，如图8，其中，易感者和恢复者数量最终趋于稳定，而潜伏者，感染者以及隐性感染者数量随时间逐渐趋近于零，因此手足口病在一定时间内会消亡。当选取接种率$v=0.0006$ ，此时计算基本再生数为$R_0=1.2206>1$ ，潜伏者、感染者以及隐性感染者的数量不为零，如图9，手足口病呈现持续流行态势。因此当疫苗接种达到一定的阈值时，才能很好的控制手足口病疫情。

Figure 4. The fitting of System (1) to the monthly reported data from Mainland China in 2022

图4. 系统(1)与中国大陆2022年月报道数据的拟合

Figure 5. Plot $R_0$ with inoculation rate $v$

图5. $R_0$ 与接种率$v$ 的关系图

Table 1. Parameter values and sources

表1. 参数取值及来源

参数	取值	来源	参数	取值	来源
$\Lambda$	752,000	[14]	$d$	0.0128	[14]
${\beta }_1$	2 × 10^−^7	拟合	${\beta }_2$	9.31 × 10^−^9	拟合
$v$	3.6 × 10^−^3	拟合	$\rho$	0.36	拟合
$\eta$	7.5	[16]	${\gamma }_1$	3.5293	[16]
${\gamma }_2$	3.5293	[16]	$\delta$	6 × 10^−^7	[15]
$b$	9 × 10^−^3	拟合	$c$	8 × 10^−^2	拟合

Table 2. Initial values and sources

表2. 初始值及来源

初值	取值	来源	初值	取值	来源
$S\left(0\right)$	6 × 10^7	[14]	$E\left(0\right)$	1.1 × 10^5	拟合
$I\left(0\right)$	2484	[15]	$I_e$	7.4 × 10^4	拟合
$R\left(0\right)$	2 × 10^4	拟合

Figure 6. Impact of vaccination on susceptible persons

图6. 疫苗接种对易感者的影响

Figure 7. Impact of vaccination on infected persons

图7. 疫苗接种对感染者的影响

Figure 8. Plots of individual bins over time for $R_0<1$

图8. $R_0<1$ 时各个仓室随时间变化图

Figure 9. Plot of individual bins with time for $R_0>1$

图9. $R_0>1$ 时各个仓室随时间变化图

7. 总结

本文建立了一个具有有限医疗资源手足口病模型，并分析了疫苗接种对手足口病的影响。通过计算得到了模型的基本再生数和无病平衡点，讨论了地方病平衡点的存在性条件，证明了无病平衡点和地方病平衡点的稳定性，此外模型在一定条件下会发生后向分支。最后结合中国大陆2023年手足口病月报告病例数进行拟合，通过数值模拟来展示有限医疗资源对手足口病传播的影响。因此，建议在手足口病高发季节，加大对儿童的疫苗接种力度，从而有效控制手足口病的传播。

基金项目

重庆市教育委员会科学技术研究项目(KJQN202100709)。

重庆市研究生院级科研创新项目(2024ST003)。

NOTES

^*通讯作者。