母婴抗体与疫苗保护的百日咳疾病模型分析

Analysis of Maternal Antibodies and Vaccine Protection in a Pertussis Disease Model

摘要

根据百日咳的传播机制，建立了一个包含母体免疫和疫苗免疫的SEIVR-M传染病模型。首先利用下一代矩阵生成法得到了基本再生数，根据线性系统的稳定性理论、比较定理以及构造Lyapunov函数证明了平衡点的稳定性。然后通过数值模拟研究了各参数对疾病传播的影响。最后引入耦合衰减项，刻画了母体抗体衰减率{\delta }_M与疫苗抗体衰减率{\delta }_V的协同作用，得到了疾病消失的临界条件。结果表明，通过精准控制母体抗体与疫苗抗体的衰减动力学，可优化群体免疫效果，最终实现疾病消除的目标。

Abstract

According to the transmission mechanism of pertussis， an SEIVR-M infectious disease model including maternal immunity and vaccine immunity was established.First， the basic reproduction number was obtained by using the next-generation matrix generation method， and the stability of the equilibrium point was proved according to the stability theory of linear systems， the comparison theorem and the Lyapunov function.Second， numerical simulations were used to study the effects of various parameters on disease transmission.Finally， the coupling attenuation term was introduced to describe the synergistic effect between the maternal antibody attenuation rate {\delta }_M and the vaccine antibody attenuation rate {\delta }_V， and the critical condition for disease extinction was obtained.The results show that by precisely controlling the attenuation kinetics of maternal antibodies and vaccine antibodies， the effect of herd immunity can be optimized， and the goal of disease elimination can be finally achieved.

Graphical abstract

关键词

母体抗体 / 疫苗保护 / 百日咳 / SEIVR-M传染病模型 / 耦合关系

Key words

maternal antibody / vaccine protection / pertussis / SEIVR-M epidemic model / coupling relationship

百日咳是由百日咳鲍特菌引起的高度传染性呼吸道疾病，其典型临床特征为阵发性痉挛咳嗽，对婴幼儿健康构成严重威胁，尤其是对于未完成基础免疫接种的婴儿，其病死率可达1%~3%^［1］。尽管全球范围内已推行了百日咳疫苗（如DTP/DTaP）的接种，但近年来，即便在疫苗接种率较高的国家，百日咳的发病率仍呈现上升趋势，形成了“百日咳再现”的现象^［2］。造成这一现象的原因有很多种，包括疫苗免疫力的衰减^［3］、病原体抗原变异引发的免疫逃逸^［4］、各国疫苗接种策略存在差异^［5］以及诊断技术改进使疾病检出率提高^［6］等。

在疫苗保护效果研究方面，Tartof等^［7］基于加州百日咳免疫接种记录与监控数据建立的多状态马尔可夫模型结果显示，完成5剂的DTaP接种虽然对重症百日咳患者的保护效果可维持在85%以上，但对轻症/无症状感染患者的保护效果在接种5年后下降至50%，从而增加潜在的传播风险。Schwartz等^［8］通过研究进一步证实，无细胞百日咳疫苗的保护效力在接种后第5年下降至50%，但适时加强免疫可使保护效力回升至82%，这为优化接种策略提供了重要依据。

值得注意的是，母婴抗体的遗传与疫苗保护的相互作用正成为新研究热点。Santilli等^［9］研究表明，通过为孕产妇接种百日咳疫苗，可使新生儿抗体滴度提升4~8倍，且不会增加不良妊娠反应。但在许多国家，孕妇的疫苗接种覆盖率仍不理想，错误信息的误导、缺乏专业医疗建议以及对接种的犹豫导致了这一接种差距。Saleem等^［10］研究证实，未接种疫苗的母亲所生婴儿发生重症百日咳的风险增加5.7倍，故孕妇若选择在怀孕期间接种疫苗，可显著降低新生儿得重症百日咳的风险。

基于上述研究并根据百日咳的传播机制，建立了带有潜伏期的疫苗接种SEIVR-M传染病模型，用来研究母体免疫和疫苗免疫的联合作用对疾病传播的影响。将模型分为6个仓室：具有母体抗体的新生儿M\left(t\right)、易感人群S\left(t\right)、潜伏期人群E\left(t\right)、感染人群I\left(t\right)、疫苗接种人群V\left(t\right)和康复人群R\left(t\right)。疾病传播情况如图1所示。

图1中，\Lambda为新生儿的数量，\beta为疾病传播过程中的有效接触率，\mu为自然死亡率，\nu为疫苗突破感染系数，{\delta }_M为母体抗体衰减率，{\delta }_V为疫苗保护衰减率，\eta为疫苗失效率，\sigma为潜伏期感染率，\gamma为感染者的恢复率。基于此，建立百日咳疾病传播SEIVR-M模型为

\left\{\begin{array}{l}\frac{\mathrm{d}M}{\mathrm{d}t}=\Lambda -{\delta }_{M}M-\mu M\\ \frac{\mathrm{d}S}{\mathrm{d}t}={\delta }_{M}M+{\delta }_{V}V-\beta SI-\mu S-\nu S\\ \frac{\mathrm{d}E}{\mathrm{d}t}=\beta SI+\eta \beta VI-\sigma E-\mu E\\ \frac{\mathrm{d}I}{\mathrm{d}t}=\sigma E-\gamma I-\mu I\\ \frac{\mathrm{d}V}{\mathrm{d}t}=\nu S-{\delta }_{V}V-\eta \beta VI-\mu V\\ \frac{\mathrm{d}R}{\mathrm{d}t}=\gamma I-\mu R\end{array}\right.（1）

1 无病平衡点与基本再生数

无病平衡点（Disease-free equilibrium，DFE）是指在传染病模型中疾病完全消失、仅存于无病人群的稳态解，此时感染者数量为0，系统处于无病传播的平衡状态。基本再生数R₀是指在易感人群中，一个感染者在整个传染期内平均引起的二代病例数，是判断疾病能否流行的关键阈值。R₀越大，疾病的总感染人数越多。根据无病平衡点的定义，令感染者和潜伏者数量全为0，代入式（1）得到

\left\{\begin{array}{l}\mathrm{0}=\Lambda -{\delta }_{M}M-\mu M\\ \mathrm{0}={\delta }_{M}M+{\delta }_{V}V-\mu S-\nu S\\ \mathrm{0}=\nu S-{\delta }_{V}V-\mu V\\ \mathrm{0}=\mu R\end{array}\right.（2）

由系统（2）中的第一个方程得到M^{\mathrm{0}}=\frac{\Lambda }{{\delta }_M+\mu },由第三个方程得到V^{\mathrm{0}}=\frac{\nu }{{\delta }_V+\mu }{S}^{\mathrm{0}}。将V^{\mathrm{0}}代入第二个方程中得到

{\delta }_M{M}^{\mathrm{0}}+{\delta }_V\left(\frac{\nu }{{\delta }_V+\mu }\right)S^{\mathrm{0}}=\left(\mu +\nu \right)S^{\mathrm{0}}，

即

{\delta }_M{M}^{\mathrm{0}}=\left(\mu +\nu -\frac{{\delta }_V\nu }{{\delta }_V+\mu }\right)S^{\mathrm{0}}，

简化系数得到

\mu +\nu -\frac{{\delta }_V\nu }{{\delta }_V+\mu }=\frac{\left(\mu +\nu \right)\left({\delta }_V+\mu \right)-{\delta }_V\nu }{{\delta }_V+\mu }=\frac{\mu {\delta }_V+{\mu }^{\mathrm{2}}+\nu \mu }{{\delta }_V+\mu }=\frac{\mu \left({\delta }_V+\mu +\nu \right)}{{\delta }_V+\mu }，

故

S^{\mathrm{0}}=\frac{{\delta }_M\left({\delta }_V+\mu \right)}{\mu \left({\delta }_V+\mu +\nu \right)}{M}^{\mathrm{0}}=\frac{\Lambda {\delta }_M\left({\delta }_V+\mu \right)}{\mu \left({\delta }_M+\mu \right)\left({\delta }_V+\mu +\nu \right)}。

综上，系统（1）的无病平衡点为

E_{\mathrm{0}}=\left(M^{\mathrm{0}},S^{\mathrm{0}},\mathrm{0,0},V^{\mathrm{0}},\mathrm{0}\right)=\left(\frac{\Lambda }{{\delta }_M+\mu },\frac{\Lambda {\delta }_M\left({\delta }_V+\mu \right)}{\mu \left({\delta }_M+\mu \right)\left({\delta }_V+\mu +\nu \right)},\mathrm{0,0},\frac{\nu \Lambda {\delta }_M}{\mu \left({\delta }_M+\mu \right)\left({\delta }_V+\mu +\nu \right)},\mathrm{0}\right)。

利用下一代矩阵法^［11］得到基本再生数R₀的求解过程：

\mathit{F}=\left[\begin{array}{cc}\mathrm{0}& \beta S^{\mathrm{0}}+\eta \beta V^{\mathrm{0}}\\ \mathrm{0}& \mathrm{0}\end{array}\right],\mathit{V}=\left[\begin{array}{cc}\sigma +\mu & \mathrm{0}\\ -\sigma & \gamma +\mu \end{array}\right],{\mathit{V}}^{-\mathrm{1}}=\frac{\mathrm{1}}{\left(\sigma +\mu \right)\left(\gamma +\mu \right)}\left[\begin{array}{cc}\gamma +\mu & \mathrm{0}\\ \sigma & \sigma +\mu \end{array}\right]，

\mathit{F}{\mathit{V}}^{-\mathrm{1}}=\left[\begin{array}{cc}\frac{\sigma \left(\beta S^{\mathrm{0}}+\eta \beta V^{\mathrm{0}}\right)}{\left(\sigma +\mu \right)\left(\gamma +\mu \right)}& \frac{\left(\beta S^{\mathrm{0}}+\eta \beta V^{\mathrm{0}}\right)\left(\sigma +\mu \right)}{\left(\sigma +\mu \right)\left(\gamma +\mu \right)}\\ \mathrm{0}& \mathrm{0}\end{array}\right]，

R_{\mathrm{0}}=\rho \left(\mathit{F}{\mathit{V}}^{-\mathrm{1}}\right)=\frac{\sigma \left(\beta S^{\mathrm{0}}+\eta \beta V^{\mathrm{0}}\right)}{\left(\sigma +\mu \right)\left(\gamma +\mu \right)}=\frac{\sigma \beta \Lambda {\delta }_M}{\left(\sigma +\mu \right)\left(\gamma +\mu \right)\left({\delta }_M+\mu \right)}\left(\frac{{\delta }_V+\mu +\eta \nu }{\mu \left({\delta }_V+\mu +\nu \right)}\right)，

其中， F 为新增感染的雅可比矩阵， V 为描述恢复、死亡、内部转移的雅可比矩阵。

2 地方病平衡点

定理1 当R_{\mathrm{0}}>\mathrm{1}时，系统（1）存在唯一的地方病平衡点。

证明 令系统（1）左端均为0，求解得到

M^{\mathrm{*}}=\frac{\Lambda }{{\delta }_M+\mu },E^{\mathrm{*}}=\frac{\gamma +\mu }{\sigma }{I}^{\mathrm{*}},V^{\mathrm{*}}=\frac{\nu S^{\mathrm{*}}}{{\delta }_V+\mu +\eta \beta I^{\mathrm{*}}},

S^{\mathrm{*}}=\frac{\left(\sigma +\mu \right)\left(\gamma +\mu \right)}{\sigma \beta }\frac{\mathrm{1}}{\mathrm{1}+\frac{\eta \beta \nu I^{\mathrm{*}}}{\left({\delta }_V+\mu +\eta \beta I^{\mathrm{*}}\right)\left(\sigma +\mu \right)}}。

迭代求解得到

I^{\mathrm{*}}=\frac{\mu \left(R_{\mathrm{0}}-\mathrm{1}\right)}{\beta \left(\mathrm{1}+\frac{\eta \nu \left(\gamma +\mu \right)}{\left({\delta }_V+\mu \right)\left(\sigma +\mu \right)}\right)}，

故当且仅当R_{\mathrm{0}}>\mathrm{1}时，I^{\mathrm{*}}>\mathrm{0},地方病平衡点存在，定理1得证。

3 无病平衡点的局部稳定性

定理2 当时，DFE局部渐近稳定。

证明 系统（1）在DFE处的雅可比矩阵为

\mathit{J}=\left[\begin{array}{cccccc}-\left({\delta }_M+\mu \right)& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}\\ {\delta }_M& -\left(\mu +\nu \right)& \mathrm{0}& -\beta S^{\mathrm{0}}& {\delta }_V& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& -\left(\sigma +\mu \right)& \beta S^{\mathrm{0}}+\eta \beta V^{\mathrm{0}}& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& \sigma & -\left(\gamma +\mu \right)& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& \nu & \mathrm{0}& -\eta \beta V^{\mathrm{0}}& -\left({\delta }_V+\mu \right)& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& \mathrm{0}& \gamma & \mathrm{0}& -\mu \end{array}\right]。

分块计算特征值：

1） M子系统：。

2） S-V子系统：特征方程为{\lambda }^{\mathrm{2}}+\left(\mathrm{2}\mu +\nu +{\delta }_V\right)\lambda +\left(\mu +\nu \right)\left({\delta }_V+\mu \right)-\nu {\delta }_V=\mathrm{0}。

由Routh-Hurwitz判据^［12］得到特征方程系数均为正且行列式{\mu }^{\mathrm{2}}+\mu {\delta }_V+\mu \nu >\mathrm{0}，因此，{\lambda }_{\mathrm{2}}\mathrm{、}{\lambda }_{\mathrm{3}}实部为负。

3） E-I子系统：特征方程为{\lambda }^{\mathrm{2}}+\left(\sigma +\mu +\gamma +\mu \right)\lambda +\left(\sigma +\mu \right)\left(\gamma +\mu \right)-\sigma \left(\beta S^{\mathrm{0}}+\eta \beta V^{\mathrm{0}}\right)=\mathrm{0}。

由于\sigma +\mathrm{2}\mu +\gamma >\mathrm{0}恒成立，且当时，\left(\sigma +\mu \right)\left(\gamma +\mu \right)-\sigma \left(\beta S^{\mathrm{0}}+\eta \beta V^{\mathrm{0}}\right)>\mathrm{0}，故方程特征根{\lambda }_{\mathrm{4}}\mathrm{、}{\lambda }_{\mathrm{5}}实部为负。

4） R子系统：。

综上，当时，所有特征值实部均为负，故DFE局部渐近稳定；当R_{\mathrm{0}}>\mathrm{1}时，E-I子系统存在正实部特征值，DFE不稳定，定理2得证。

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

定理3 当时，DFE全局渐近稳定。

证明 将系统（1）分为两部分，非疾病变量M\mathrm{、}S\mathrm{、}V\mathrm{、}R以及疾病变量E\mathrm{、}I。在无病的情况下E=I=\mathrm{0}，故系统（1）简化为

\left\{\begin{array}{l}\frac{\mathrm{d}M}{\mathrm{d}t}=\Lambda -{\delta }_{M}M-\mu M\\ \frac{\mathrm{d}S}{\mathrm{d}t}={\delta }_{M}M+{\delta }_{V}V-\mu S-\nu S\\ \frac{\mathrm{d}V}{\mathrm{d}t}=\nu S-{\delta }_{V}V-\mu V\\ \frac{\mathrm{d}R}{\mathrm{d}t}=-\mu R\end{array}\right.（3）

其中，R的方程是\frac{\mathrm{d}R}{\mathrm{d}t}=-\mu R，故其解为R\left(t\right)=R\left(\mathrm{0}\right){\mathrm{e}}^{-\mu t}，且当t\to \mathrm{\infty }时，R\left(t\right)\to \mathrm{0}。对于M，\frac{\mathrm{d}M}{\mathrm{d}t}=\Lambda -{\delta }_{M}M-\mu M是一个线性方程，其解为

M\left(t\right)=\frac{\Lambda }{{\delta }_M+\mu }+\left(M\left(\mathrm{0}\right)-\frac{\Lambda }{{\delta }_M+\mu }\right){\mathrm{e}}^{-\left({\delta }_M+\mu \right)t}，

且当t\to \mathrm{\infty }时，M\left(t\right)\to \frac{\Lambda }{{\delta }_M+\mu }=M^{\mathrm{0}}。而S和V是一个线性系统，其解可表示成矩阵的形式：

\frac{\mathrm{d}}{\mathrm{d}t}\left(\begin{array}{c}S\\ V\end{array}\right)=\mathit{A}\left(\begin{array}{c}S\\ V\end{array}\right)+\left(\begin{array}{c}{\delta }_{M}M\\ \mathrm{0}\end{array}\right),\mathit{A}=\left(\begin{array}{cc}-\left(\mu +\nu \right)& {\delta }_V\\ \nu & -\left({\delta }_V+\mu \right)\end{array}\right)（4）

因为矩阵\mathit{A}的特征值具有负实部，故系统（4）是全局渐近稳定的，且S和V会趋向于其稳定值S^{\mathrm{0}}和V^{\mathrm{0}}。

疾病变量E和I所构成的子系统是一个线性系统，且该子系统的系数矩阵依赖于S\left(t\right)和V\left(t\right)的取值。由于S\left(t\right)\to S^{\mathrm{0}}，V\left(t\right)\to V^{\mathrm{0}}，故考虑其极限系统

\left\{\begin{array}{l}\frac{\mathrm{d}E}{\mathrm{d}t}=\beta S^{\mathrm{0}}I+\eta \beta V^{\mathrm{0}}I-\left(\sigma +\mu \right)E\\ \frac{\mathrm{d}I}{\mathrm{d}t}=\sigma E-\left(\gamma +\mu \right)I\end{array}\right.（5）

系统（5）的稳定性由矩阵 J 的特征值决定， J 表示为

\mathit{J}=\left[\begin{array}{cc}-\left(\sigma +\mu \right)& \beta S^{\mathrm{0}}+\eta \beta V^{\mathrm{0}}\\ \sigma & -\left(\gamma +\mu \right)\end{array}\right]。

根据Routh-Hurwitz判据^［13］，当\left(\sigma +\mu \right)\left(\gamma +\mu \right)>\sigma \left(\beta S^{\mathrm{0}}+\eta \beta V^{\mathrm{0}}\right)时，该矩阵的特征值具有负实部。因此，当时，极限系统E和I趋向于其稳定值0。

接下来需证明对于任何初始条件，E\left(t\right)和I\left(t\right)均趋向于0。因为S\left(t\right)和V\left(t\right)趋向于S^{\mathrm{0}}和V^{\mathrm{0}}，且S\left(t\right)\le S^{\mathrm{0}}，V\left(t\right)\le V^{\mathrm{0}}，根据比较定理^［14-18］，考虑辅助系统

\left\{\begin{array}{l}\frac{\mathrm{d}E}{\mathrm{d}t}=\beta S^{\mathrm{0}}\tilde{I}+\eta \beta V^{\mathrm{0}}\tilde{I}-\left(\sigma +\mu \right)\tilde{E}\\ \frac{\mathrm{d}I}{\mathrm{d}t}=\sigma \tilde{E}-\left(\gamma +\mu \right)\tilde{I}\end{array}\right.（6）

故系统（1）的E\left(t\right)和I\left(t\right)被辅助系统（6）的解所控制。而辅助系统（6）在时趋向于0，所以系统（1）的E\left(t\right)和I\left(t\right)也趋向于0。

综上，非疾病变量M\mathrm{、}S\mathrm{、}V\mathrm{、}R在t\to \mathrm{\infty }时趋向于无病稳定值M^{\mathrm{0}}\mathrm{、}{S}^{\mathrm{0}}\mathrm{、}{V}^{\mathrm{0}}\mathrm{、}\mathrm{0}；疾病变量E\mathrm{、}I在时被一个趋向于0的线性系统所控制，因此也趋向于0，故无病平衡点在时全局渐近稳定，定理3得证。

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

定理4 当R_{\mathrm{0}}>\mathrm{1}时，地方病平衡点局部渐近稳定。

证明 在地方病平衡点处的雅可比矩阵为

\mathit{J}=\left[\begin{array}{cccccc}-\left({\delta }_M+\mu \right)& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}\\ {\delta }_M& -\left(\beta I^{\mathrm{*}}+\mu +\nu \right)& \mathrm{0}& -\beta S^{\mathrm{*}}& {\delta }_V& \mathrm{0}\\ \mathrm{0}& \beta I^{\mathrm{*}}& -\left(\sigma +\mu \right)& \beta S^{\mathrm{*}}+\eta \beta V^{\mathrm{*}}& \eta \beta I^{\mathrm{*}}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& \sigma & -\left(\gamma +\mu \right)& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& \nu & \mathrm{0}& -\eta \beta V^{\mathrm{*}}& -\left({\delta }_V+\mu +\eta \beta I^{\mathrm{*}}\right)& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& \mathrm{0}& \gamma & \mathrm{0}& -\mu \end{array}\right]，

其特征值为-\left({\delta }_M+\mu \right)和-\mu，其余4个特征根来自以下子矩阵：

\tilde{\mathit{J}}=\left[\begin{array}{cccc}-\left(\beta I^{\mathrm{*}}+\mu +\nu \right)& \mathrm{0}& -\beta S^{\mathrm{*}}& {\delta }_V\\ \beta I^{\mathrm{*}}& -\left(\sigma +\mu \right)& \beta S^{\mathrm{*}}+\eta \beta V^{\mathrm{*}}& \eta \beta I^{\mathrm{*}}\\ \mathrm{0}& \sigma & -\left(\gamma +\mu \right)& \mathrm{0}\\ \nu & \mathrm{0}& -\eta \beta V^{\mathrm{*}}& -\left({\delta }_V+\mu +\eta \beta I^{\mathrm{*}}\right)\end{array}\right]（7）

计算矩阵（7）的特征多项式得到P\left(\lambda \right)={\lambda }^{\mathrm{4}}+a_{\mathrm{1}}{\lambda }^{\mathrm{3}}+a_{\mathrm{2}}{\lambda }^{\mathrm{2}}+a_{\mathrm{3}}\lambda +a_{\mathrm{4}}。其中，

a_{\mathrm{1}}=\beta I^{\mathrm{*}}\left(\mathrm{1}+\eta \right)+{\delta }_V+\sigma +\gamma +\nu +\mathrm{4}\mu ~\psi \left(I^{\mathrm{*}}\right)，

\begin{array}{l}{a}_{\mathrm{2}}=\left(\beta I^{\mathrm{*}}+\mu +\nu \right)\left(\sigma +\mu \right)+\left(\beta I^{\mathrm{*}}+\mu +\nu \right)\left(\gamma +\mu \right)+\left(\beta I^{\mathrm{*}}+\mu +\nu \right)\left({\delta }_V+\mu +\eta \beta I^{\mathrm{*}}\right)\\ +\left(\sigma +\mu \right)\left(\gamma +\mu \right)+\left(\sigma +\mu \right)\left({\delta }_V+\mu +\eta \beta I^{\mathrm{*}}\right)+\left(\gamma +\mu \right)\left({\delta }_V+\mu +\eta \beta I^{\mathrm{*}}\right)\\ -\sigma \beta I^{\mathrm{*}}-\nu {\delta }_V~\psi \left(I^{\mathrm{*}\mathrm{2}}\right),\end{array}

\begin{array}{l}{a}_{\mathrm{3}}=\left(\beta I^{\mathrm{*}}+\mu +\nu \right)\left(\sigma +\mu \right)\left(\gamma +\mu \right)+\left(\beta I^{\mathrm{*}}+\mu +\nu \right)\left(\sigma +\mu \right)\left({\delta }_V+\mu +\eta \beta I^{\mathrm{*}}\right)\\ +\left(\beta I^{\mathrm{*}}+\mu +\nu \right)\left(\gamma +\mu \right)\left({\delta }_V+\mu +\eta \beta I^{\mathrm{*}}\right)+\left(\sigma +\mu \right)\left(\gamma +\mu \right)\left({\delta }_V+\mu +\eta \beta I^{\mathrm{*}}\right)\\ -\sigma \beta I^{\mathrm{*}}\left({\delta }_V+\mu +\eta \beta I^{\mathrm{*}}\right)-\nu {\delta }_V\left(\gamma +\mu \right)~\psi \left(I^{\mathrm{*}\mathrm{2}}\right),\end{array}

a_{\mathrm{4}}=\mu \left(\sigma +\mu \right)\left(\gamma +\mu \right)\left[\left({\delta }_V+\mu \right)\left(R_{\mathrm{0}}-\mathrm{1}\right)+\eta \beta I^{\mathrm{*}}\right]+\psi \left(I^{\mathrm{*}\mathrm{2}}\right)~\psi \left(I^{\mathrm{*}\mathrm{2}}\right)，

故a_{\mathrm{1}}{a}_{\mathrm{2}}{a}_{\mathrm{3}}~\psi \left(I^{\mathrm{*}\mathrm{5}}\right)，{a_{\mathrm{3}}}^{\mathrm{2}}~\psi \left(I^{\mathrm{*}\mathrm{4}}\right)，{a_{\mathrm{1}}}^{\mathrm{2}}{a}_{\mathrm{4}}~\psi \left(I^{\mathrm{*}\mathrm{2}}\right)\cdot \psi \left(I^{\mathrm{*}\mathrm{2}}\right)=\psi \left(I^{\mathrm{*}\mathrm{4}}\right)。当R_{\mathrm{0}}>\mathrm{1}时，I^{\mathrm{*}}足够大，各系数均大于0，且a_{\mathrm{1}}{a}_{\mathrm{2}}{a}_{\mathrm{3}}>{a_{\mathrm{3}}}^{\mathrm{2}}+{a_{\mathrm{1}}}^{\mathrm{2}}{a}_{\mathrm{4}}。由Routh-Hurwitz稳定性^［19］得到地方病平衡点局部渐近稳定，定理4得证。

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

定理5 当R_{\mathrm{0}}>\mathrm{1}时，地方病平衡点全局渐近稳定。

证明 构造Lyapunov函数：

L=S^{\mathrm{*}}\left(\frac{S}{S^{\mathrm{*}}}-\mathrm{1}-\mathrm{l}\mathrm{n}\frac{S}{S^{\mathrm{*}}}\right)+E^{\mathrm{*}}\left(\frac{E}{E^{\mathrm{*}}}-\mathrm{1}-\mathrm{l}\mathrm{n}\frac{E}{E^{\mathrm{*}}}\right)+\frac{\sigma +\mu }{\sigma }{I}^{\mathrm{*}}\left(\frac{I}{I^{\mathrm{*}}}-\mathrm{1}-\mathrm{l}\mathrm{n}\frac{I}{I^{\mathrm{*}}}\right)+V^{\mathrm{*}}\left(\frac{V}{V^{\mathrm{*}}}-\mathrm{1}-\mathrm{l}\mathrm{n}\frac{V}{V^{\mathrm{*}}}\right),

对其求导得到

\begin{array}{l}\frac{\mathrm{d}L}{\mathrm{d}\mathrm{t}}=-\left(\mu +\nu \right)\frac{{\left(S-S^{\mathrm{*}}\right)}^{\mathrm{2}}}{S}+\beta S^{\mathrm{*}}{I}^{\mathrm{*}}\left(\mathrm{2}-\frac{S^{\mathrm{*}}}{S}-\frac{SI}{S^{\mathrm{*}}{I}^{\mathrm{*}}}\right)+\left(\sigma +\mu \right)E^{\mathrm{*}}\left(\mathrm{2}-\frac{E^{\mathrm{*}}}{E}-\frac{E{I}^{\mathrm{*}}}{E^{\mathrm{*}}I}\right)\\ +\left(\sigma +\mu \right)I^{\mathrm{*}}\left(\mathrm{2}-{\frac{I}{I}}^{\mathrm{*}}-\frac{I{E}^{\mathrm{*}}}{I^{\mathrm{*}}E}\right)-{\delta }_V\frac{{\left(V-V^{\mathrm{*}}\right)}^{\mathrm{2}}}{V}+\eta \beta V^{\mathrm{*}}{I}^{\mathrm{*}}\left(\mathrm{2}-\frac{V^{\mathrm{*}}}{V}-\frac{VI}{V^{\mathrm{*}}{I}^{\mathrm{*}}}\right)\mathrm{。}\end{array}

利用几何平均不等式得到\frac{S^{\mathrm{*}}}{S}+\frac{SI}{S^{\mathrm{*}}{I}^{\mathrm{*}}}\ge \mathrm{2}\sqrt[]{\frac{S^{\mathrm{*}}}{S}+\frac{SI}{S^{\mathrm{*}}{I}^{\mathrm{*}}}}=\mathrm{2}\sqrt[]{\frac{I}{I^{\mathrm{*}}}}\ge \mathrm{2}\mathrm{。}故\mathrm{2}-\frac{S^{\mathrm{*}}}{S}-\frac{SI}{S^{\mathrm{*}}{I}^{\mathrm{*}}}\le \mathrm{0}，同理，最后得到

\frac{\mathrm{d}L}{\mathrm{d}\mathrm{t}}\le -\left(\mu +\nu \right)\frac{{\left(S-S^{\mathrm{*}}\right)}^{\mathrm{2}}}{S}-{\delta }_V\frac{{\left(V-V^{\mathrm{*}}\right)}^{\mathrm{2}}}{V}\le \mathrm{0}\mathrm{。}

当且仅当S=S^{\mathrm{*}},E=E^{\mathrm{*}},I=I^{\mathrm{*}},V=V^{\mathrm{*}}时，\frac{\mathrm{d}L}{\mathrm{d}\mathrm{t}}=\mathrm{0}，此时M\to M^{\mathrm{*}},R\to R^{\mathrm{*}}。故当R_{\mathrm{0}}>\mathrm{1}时，地方病平衡点全局渐近稳定，定理5得证。

7 数值模拟和灵敏度分析

基于Matlab软件，采用Dormand-Prince方法对模型进行分析求解，以验证上述结论。当\Lambda =\mathrm{100}，\mu =\mathrm{0.01}，{\delta }_M=\mathrm{0.002}，{\delta }_V=\mathrm{0.001}，\beta =\mathrm{0.0001}，\eta =\mathrm{0.7}，\sigma =\mathrm{0.1}，\gamma =\mathrm{0.1}，\nu =\mathrm{0.8}时，。结合图2（a）可以看出，与疾病感染程度有关的潜伏期人群仓室和染病人群仓室随着时间的推移感染人数逐渐趋于0，疾病逐渐消失。当\Lambda =\mathrm{100}，\mu =\mathrm{0.01}，{\delta }_M=\mathrm{0.05}，{\delta }_V=\mathrm{0.03}，\beta =\mathrm{0.0002}，\eta =\mathrm{0.5}，\sigma =\mathrm{0.1}，\gamma =\mathrm{0.07}，\nu =\mathrm{0.5}时，R_{\mathrm{0}}>\mathrm{1}。结合图2（b）可得，潜伏期人群仓室和染病人群仓室的感染人数随着时间的推移而上升，然后逐渐趋于一个稳定值，即染病人数不会产生太大的变化，疾病变成了一种长期存在的地方病。

接下来通过模拟\beta \mathrm{、}\eta \mathrm{、}\sigma \mathrm{、}\gamma \mathrm{、}{\delta }_M\mathrm{、}{\delta }_V不同取值，研究其对总感染人数的影响。

如图3（a）—（c）所示，当疾病传播率\beta越低时，疾病在人群中的有效传播能力越弱，单个感染者在传染期内平均能传染的人数越少。在这种情况下，基本再生数R_{\mathrm{0}}也会越小，疾病更容易被控制；当疫苗突破感染系数\eta越小时，疫苗对疾病的防护效果越强，更能有效地阻断感染和传播。此时，即使有感染者存在，接种疫苗人群被感染的概率也会显著降低，从而减少了继发病例的数量；当疾病的潜伏期感染率\sigma越小时，感染者在潜伏期内传播疾病的概率越低，感染者在出现症状前传染给他人的机会越小，这也会显著降低疾病的有效传播能力。

如图4（a）—（c）所示，当疾病的恢复率\gamma越高时，单个感染者在康复前能够传播给其他人的机会越低，从而直接减少了有效接触的数量，R_{\mathrm{0}}会相应下降；当母体抗体衰减速率{\delta }_M越慢时，新生儿通过胎盘或母乳获得的被动免疫保护作用持续时间越长，使其在生命早期更不容易感染疾病，这显著降低了易感人群的比例，从而减少疾病在人群中的传播机会；当疫苗保护衰减速率{\delta }_V越慢时，接种者维持有效免疫防护的时间越长，群体中具有抵抗力的人口比例保持较高水平，从而持续抑制疾病的传播。

8 母体抗体与疫苗抗体的衰减率耦合分析

定义耦合衰减项\kappa \left({\delta }_M,{\delta }_V\right)=\frac{\mathrm{2}{\delta }_M{\delta }_V}{{\delta }_M+{\delta }_V}，得到新的耦合方程

\left\{\begin{array}{l}\frac{\mathrm{d}M}{\mathrm{d}t}=\Lambda -\left({\delta }_M+\mu \right)M-\kappa MV\\ \frac{\mathrm{d}S}{\mathrm{d}t}={\delta }_{M}M+{\delta }_{V}V-\beta SI-\left(\mu +\nu \right)S+\kappa MV\\ \frac{\mathrm{d}V}{\mathrm{d}t}=\nu S-\left({\delta }_V+\mu \right)V-\eta \beta VI-\kappa MV\end{array}\right.（8）

计算得到系统（8）的基本再生数R₁₀为

R_{\mathrm{10}}=\frac{\beta \sigma }{\left(\sigma +\mu \right)\left(\gamma +\mu \right)}\left(S^{\mathrm{10}}+\eta V^{\mathrm{10}}\right),

其中，S^{\mathrm{10}}=\frac{\Lambda \left({\delta }_V+\mu \right)}{\left({\delta }_M+\mu \right)\left({\delta }_V+\mu \right)+\kappa \Lambda }，V^{\mathrm{10}}=\frac{\nu \Lambda }{\left({\delta }_V+\mu \right)\left(\mu +\nu \right)+\kappa \nu }。

定义临界耦合强度{\kappa }_c，该参数是当基本再生数R_{\mathrm{0}}=\mathrm{1}时的耦合强度值，即R_{\mathrm{0}}\left({\kappa }_c\right)=\mathrm{1}，此时系统将从无病稳态转变为地方病稳态。

基于临界耦合强度定义将S^{\mathrm{10}}、V^{\mathrm{10}}代入R_{\mathrm{0}}得到

\frac{\beta \sigma }{\left(\sigma +\mu \right)\left(\gamma +\mu \right)}\left(\frac{\Lambda \left({\delta }_V+\mu \right)}{\left({\delta }_M+\mu \right)\left({\delta }_V+\mu \right)+{\kappa }_c\Lambda }+\eta \frac{\nu \Lambda }{\left({\delta }_V+\mu \right)\left(\mu +\nu \right)+{\kappa }_c\nu }\right)=\mathrm{1}\mathrm{。}

假设此时疫苗接种率远大于自然死亡率（\nu \gg \mu）且\eta \approx \mathrm{1}，故可以近似忽略第二项得到

\frac{\beta \sigma }{\left(\sigma +\mu \right)\left(\gamma +\mu \right)}\frac{\Lambda \left({\delta }_V+\mu \right)}{\left({\delta }_M+\mu \right)\left({\delta }_V+\mu \right)+{\kappa }_c\Lambda }\approx \mathrm{1},

解得

{\kappa }_c\approx \frac{\beta \sigma \left({\delta }_V+\mu \right)}{\left(\sigma +\mu \right)\left(\gamma +\mu \right)}-\frac{\left({\delta }_M+\mu \right)\left({\delta }_V+\mu \right)}{\Lambda }\mathrm{。}

因此，当时，耦合作用最强，疾病灭绝；当\kappa >{\kappa }_c时，耦合不足，疾病将持续存在。

令\Psi =\frac{\beta \sigma \Lambda }{\left(\sigma +\mu \right)\left(\gamma +\mu \right)}，基于上述假设可得{\delta }_M\mathrm{、}{\delta }_V需满足

\left({\delta }_M+\mu \right)\left({\delta }_V+\mu \right)+\frac{\mathrm{2}{\delta }_M{\delta }_V\Lambda }{{\delta }_M+{\delta }_V}>\Psi \left({\delta }_V+\mu \right),

即当\left({\delta }_M+\mu \right)\left({\delta }_V+\mu \right)\left({\delta }_M+{\delta }_V\right)+\mathrm{2}{\delta }_M{\delta }_V\Lambda >\Psi \left({\delta }_V+\mu \right)\left({\delta }_M+{\delta }_V\right)时，疾病处于一个稳定的消失状态。当{\delta }_M={\delta }_V时，耦合效果最强，将上文数值模拟中的数据代入可得

\frac{R_{\mathrm{0}}-R_{\mathrm{10}}}{R_{\mathrm{0}}}\approx \frac{\kappa \Lambda }{{\left(\delta +\mu \right)}^{\mathrm{2}}+\kappa \Lambda }\approx \mathrm{20}\mathrm{\%}\mathrm{。}

综上可得，母体抗体衰减率与疫苗抗体衰减率的耦合关系对抑制疾病传播具有决定性作用。母体抗体与疫苗抗体的同步衰减可有效减少易感者积累，抑制疾病传播。当{\delta }_M={\delta }_V时，耦合效果最强，可降低R_{\mathrm{0}}约15%~25%。而临界耦合系数{\kappa }_c为稳定性阈值，在实际中需确保，进而可实现疾病的消除。

9 结论

基于百日咳的传播机制，建立了一个同时考虑母体免疫和疫苗免疫的SEIVR-M传染病模型。首先，利用下一代矩阵法推导出了模型的基本再生数R_{\mathrm{0}}，并根据线性系统的稳定性理论以及比较定理，证明了当时无病平衡点的局部和全局渐近稳定性。其次，通过构造Lyapunov函数证明了当R_{\mathrm{0}}>\mathrm{1}时地方病平衡点的全局渐近稳定性。然后，使用Matlab软件进行数值模拟，验证了上述结论。进一步，通过全局灵敏度分析方法，量化了各参数对R_{\mathrm{0}}的影响程度。最后引入了耦合衰减项，研究了母体抗体衰减率{\delta }_M与疫苗抗体衰减率{\delta }_V的协同作用，得到了疾病消失的临界条件。结果表明，通过调控母体免疫与疫苗免疫的抗体衰减动力学，可提升群体免疫效果，进而抑制疾病的传播。但此模型简化了真实的免疫过程，难以准确量化母体抗体与疫苗抗体间的复杂相互作用，此为未来模型优化的重要方向。

参考文献

[1] CHERRY J D.Pertussis：Challenges today and for the future［J］.PLoS Pathogens，2013，9（7）：e1003418.

[2] CLARK T A.Changing pertussis epidemiology：Everything old is new again［J］.Journal of Infectious Diseases，2014，209（7）：978-981.

[3] IRULAPPAN M， JACOB J J， MADHUMATHI J，et al.Pertussis in India：Vaccine-driven evolution，waning immunity，and the urgent need for Tdap boosters［J］.Indian Journal of Medical Microbiology，2025，55：100846.

[4] BART M J， HARRIS S R， ADVANI A，et al.Global population structure and evolution of Bordetella pertussis and their relationship with vaccination［J］.mBio，2014，5（2）：e01074.

[5] DOMENECH DE CELLÈS M， MAGPANTAY F M G， KING A A，et al.The pertussis enigma：Reconciling epidemiology，immunology and evolution［J］.Proceedings of the Royal Society B：Biological Sciences，2016，283（1822）：20152309.

[6] KILGORE P E， SALIM A M， ZERVOS M J，et al.Pertussis：Microbiology，disease，treatment，and prevention［J］.Clinical Microbiology Reviews，2016，29（3）：449-486.

[7] TARTOF S Y， LEWIS M， KENYON C，et al.Waning immunity to pertussis following 5 doses of DTaP［J］.Pediatrics，2013，131（4）：e1047-e1052.

[8] SCHWARTZ K L， KWONG J C， DEEKS S L，et al.Effectiveness of pertussis vaccination and duration of immunity［J］.Canadian Medical Association Journal，2016，188（16）：e399-e406.

[9] SANTILLI V， SGRULLETTI M， COSTAGLIOLA G，et al.Maternal immunization：Current evidence，progress，and challenges［J］.Vaccines，2025，13（5）：450.

[10] SALEEM A I， NOORELAHI A K， ZAKI N G，et al.Neonatal pertussis in an infant born to an immigrant unimmunized mother in Saudi Arabia［J］.Journal of Clinical Neonatology，2025，14（2）：77-79.

[11] 潘玮，靳祯.一类信息传播动力学模型与控制策略分析［J］.内蒙古大学学报（自然科学版），2017，48（6）：601-606.

[12] BATOOL H， KHAN I， LI W Y，et al.Fractional modeling and numerical investigations of COVID-19 epidemic model with non-singular fractional derivatives：A case study［J］.Scientific Reports，2025，15（1）：13256.

[13] YADAV S， KUMAR D， SINGH J，et al.Analysis and dynamics of fractional order covid-19 model with memory effect［J］.Results in Physics，2021，24：104017.

[14] GUO H， MICHAEL L， SHUAI Z.Global stability of the endemic equilibrium of multigroup SIR epidemic models［J］.Canadian Applied Mathematics Quarterly，2006，14（3）：259-284.

[15] GUO H B， LI M Y.Global dynamics of a staged-progression model with amelioration for infectious diseases［J］.Journal of Biological Dynamics，2008，2（2）：154-168.

[16] ESSAK A A， BOUKANJIME B.Global stability of an SAIRS epidemic model with vaccinations，transient immunity and treatment［J］.Nonlinear Analysis：Real World Applications，2023，73：103887.

[17] WANG Z L， YANG D S， MA T D，et al.Stability analysis for nonlinear fractional-order systems based on comparison principle［J］.Nonlinear Dynamics，2014，75（1）：387-402.

[18] GUO Y M， LI T T.Optimal control and stability analysis of an online game addiction model with two stages［J］.Mathematical Methods in the Applied Sciences，2020，43（7）：4391-4408.

[19] VAN DEN DRIESSCHE P， WATMOUGH J.Reproduction numbers and sub-threshold endemic equilibria for compartmental models of disease transmission［J］.Mathematical Biosciences，2002，180：29-48.

基金资助

国家自然科学基金项目(11971278)

山西省自然科学基金项目(202203021211086)