基于免疫力变化的百日咳疫苗有效性模型分析

Analysis of pertussis vaccine efficacy model based on immune system changes

摘要

百日咳是一种常见的呼吸道传染性疾病，可通过飞沫传播，并存在二次感染的可能性.为探究疫苗有效性对疾病传播的影响，本文基于免疫力动态变化特征，将免疫人群细分为高度免疫个体以及免疫减弱个体，在此基础上建立了百日咳{\mathrm{S}}_{\mathrm{1}}\mathrm{E}\mathrm{I}\mathrm{V}{\mathrm{R}}_{\mathrm{2}}{\mathrm{R}}_{\mathrm{1}}{\mathrm{S}}_{\mathrm{2}}传播动力学模型，分析并计算得出疾病流行的阈值，研究系统平衡点的稳定性.进一步采用参数耦合分析方法，揭示疫苗有效性与免疫力衰退率对百日咳传播动力学的联合作用机制.研究结论表明，通过优化疫苗效力与免疫力持续时间的协同关系，可以确定一个最优平衡点，从而显著降低百日咳的流行规模；将非药物干预措施与提高人群免疫水平相结合，可形成更有效的综合防控策略.

Abstract

Pertussis is a common respiratory infectious disease transmitted through droplets， with the possibility of reinfection. To investigate the impact of vaccine efficacy on disease transmission， we developed a pertussis {\mathrm{S}}_{\mathrm{1}}\mathrm{E}\mathrm{I}\mathrm{V}{\mathrm{R}}_{\mathrm{2}}{\mathrm{R}}_{\mathrm{1}}{\mathrm{S}}_{\mathrm{2}} transmission dynamics model based on the dynamic characteristics of immunity， which categorizes the immune population into highly immune individuals and waning immunity individuals. The epidemic threshold of the disease was analytically derived， and the stability of the system’s equilibrium points was examined. Furthermore， a parameter coupling analysis method was employed to elucidate the joint mechanism of vaccine efficacy and immunity waning rate on pertussis transmission dynamics. The findings demonstrate that by optimizing the synergistic relationship between vaccine efficacy and immunity duration， an optimal equilibrium point can be identified， significantly reducing the scale of pertussis outbreaks. Combining non-pharmaceutical interventions with enhanced population immunity levels can form a more effective comprehensive prevention and control strategy.

Graphical abstract

关键词

传染病模型 / 免疫力变化 / 疫苗接种 / 全局稳定性 / 参数耦合

Key words

infectious disease model / immunity dynamics / vaccination / global stability / parameter coupling

百日咳（pertussis）是一种由百日咳杆菌引起的急性呼吸道传染病，主要通过飞沫传播.它具有高度的传染性、致命性，且在全球范围内广泛传播，尤其在未接种疫苗的婴幼儿中发病率较高.百日咳疫苗的接种已成为控制该疾病的最有效手段.根据相关研究，疫苗接种能够显著减少百日咳的发病率和致死率，尤其是在高风险人群中^［1］.然而，百日咳的流行现象依然存在，美国等发达国家仍然会出现周期性的百日咳疫情.

虽然疫苗接种是目前最为有效的传染病控制手段，但研究发现，接种疫苗后并不会产生永久的免疫力，且免疫力会随着时间的推移而减弱^［2-4］.通过疫苗接种和自然感染获得的免疫力在体内持续的时间不同，前者持续时间通常短于后者.自然感染的患者抗体浓度要高于疫苗接种者^［5-6］.例如，Broutin等^［7］在一项针对塞内加尔人口的研究中发现，未接种疫苗的儿童两次百日咳发作的平均间隔为7.1年，而接种过疫苗的儿童为5.1年.二次感染能够增强减弱的免疫力，虽然这种免疫力仍然比初次感染时弱^［8-9］.

近年来许多研究者已尝试量化免疫动态对百日咳传播的影响，Carlsson等^［10］提出一种涉及离散免疫类别、但具有年龄和时间连续变化的偏微分方程\mathrm{S}\mathrm{I}\mathrm{S}模型.Ehrhardt等^［11］通过推导\mathrm{S}\mathrm{I}\mathrm{R}模型来描述疫苗接种和免疫力减弱的过程.Elbasha等^［12］通过构建不完美疫苗接种的仓室模型，研究了自然获得免疫力与通过接种疫苗获得免疫力的差异及其减弱的情况.Opoku-Sarkodie等^［13］则扩展了基本的\mathrm{S}\mathrm{I}\mathrm{R}\mathrm{S}模型，考虑了在没有任何疫苗接种方案的情况下，反复接触感染会导致免疫力的减弱和增强.然而，现有模型多假设疫苗接种后个体直接获得完全免疫，未充分考虑疫苗失效的现实情况.

疫苗效力是决定疾病传播控制效果的关键参数，但现实中，疫苗并不能100%地发挥作用^［14］.百日咳疫苗也不例外^［15-18］.考虑到这一实际因素，本文将免疫人群细分为高度免疫个体（近期接种疫苗或感染后获得强免疫保护者）和免疫减弱个体（疫苗保护效果随年限下降或抗体水平衰减者），研究在免疫力变化的前提下，疫苗效力对百日咳疾病传播的影响.通过此研究，能够更好地理解疫苗接种对百日咳传播的控制效应，并为政策制定者提供科学依据，从而有效控制百日咳的传播，保护公共卫生安全.

1 模型的建立

将人群分为原发易感者S_{\mathrm{1}}，潜伏期患者E，染病者I，疫苗接种个体V，高度免疫个体R_{\mathrm{2}}，免疫力减弱个体R_{\mathrm{1}}以及继发易感者S_{\mathrm{2}}.该模型的传播流程图如图1所示.

根据流程图建立以下模型：

\left\{\begin{array}{l}\frac{\mathrm{d}{S}_{\mathrm{1}}}{\mathrm{d}t}=\left(\mathrm{1}-p\right)\mathrm{\Lambda }-{\beta }_{\mathrm{1}}{S}_{\mathrm{1}}E-{\beta }_{\mathrm{2}}{S}_{\mathrm{1}}I-\epsilon S_{\mathrm{1}}-\mu S_{\mathrm{1}},\\ \frac{\mathrm{d}E}{\mathrm{d}t}={\beta }_{\mathrm{1}}{S}_{\mathrm{1}}E+{\beta }_{\mathrm{2}}{S}_{\mathrm{1}}I+w{\beta }_{\mathrm{1}}VE+w{\beta }_{\mathrm{2}}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}=\epsilon S_{\mathrm{1}}+p\mathrm{\Lambda }-w{\beta }_{\mathrm{1}}VE-w{\beta }_{\mathrm{2}}VI-\left(\mathrm{1}-w\right)b{\alpha }_{\mathrm{1}}V-\mu V,\\ \frac{\mathrm{d}{R}_{\mathrm{2}}}{\mathrm{d}t}=\gamma I+g{\beta }_{\mathrm{1}}{R}_{\mathrm{1}}E+g{\beta }_{\mathrm{2}}{R}_{\mathrm{1}}I-{\alpha }_{\mathrm{1}}{R}_{\mathrm{2}}-\mu R_{\mathrm{2}},\\ \frac{\mathrm{d}{R}_{\mathrm{1}}}{\mathrm{d}t}=\left(\mathrm{1}-w\right)b{\alpha }_{\mathrm{1}}V+{\alpha }_{\mathrm{1}}{R}_{\mathrm{2}}-g{\beta }_{\mathrm{1}}{R}_{\mathrm{1}}E-g{\beta }_{\mathrm{2}}{R}_{\mathrm{1}}I-{\alpha }_{\mathrm{2}}{R}_{\mathrm{1}}-\mu R_{\mathrm{1}},\\ \frac{\mathrm{d}{S}_{\mathrm{2}}}{\mathrm{d}t}={\alpha }_{\mathrm{2}}{R}_{\mathrm{1}}-\mu S_{\mathrm{2}}.\end{array}\right.（1）

假设个体出生时都为原发易感者，并以\mu的比率自然死亡，其中比例为p的新生儿在出生时接种了疫苗.原发易感者接种疫苗的比例为\epsilon.由于与潜伏者的成功接触，原发易感者以{\beta }_{\mathrm{1}}E的感染率获得感染.与染病者成功接触，以{\beta }_{\mathrm{2}}I的感染率获得感染.\sigma是由E到I的转化率.染病者以\gamma的速率从感染中恢复到R_{\mathrm{2}}，但其自然获得的免疫力会下降，并以{\alpha }_{\mathrm{1}}的速率减弱.接种疫苗的个体疫苗失效的可能性为w且获得性免疫力以b{\alpha }_{\mathrm{1}}的速度下降，其中，b是一个重新调整参数，说明通过疫苗接种获得的免疫力相对于自然获得的免疫力的损失.减弱个体R_{\mathrm{1}}的免疫力要么与潜伏个体或染病个体接触以g{\beta }_{\mathrm{1}}E、g{\beta }_{\mathrm{2}}I的速度上升，要么以{\alpha }_{\mathrm{2}}的速度继续下降，变成继发易感者S_{\mathrm{2}}.其中，g是一个重新调整参数.

设N为人口的总数，它满足N=S_{\mathrm{1}}+E+I+V+R_{\mathrm{2}}+R_{\mathrm{1}}+S_{\mathrm{2}}，将系统（1）中的方程相加得到

\frac{\mathrm{d}N}{\mathrm{d}t}=\mathrm{\Lambda }-\mu N，

则系统（1）的正向不变集为

.

2 平衡点和基本再生数

通过计算可得系统（1）的无病平衡点：

\begin{array}{l}{E}^{\mathrm{0}}=\left(S_{\mathrm{10}},E_{\mathrm{0}},I_{\mathrm{0}},V_{\mathrm{0}},R_{\mathrm{20}},R_{\mathrm{10}},S_{\mathrm{20}}\right)=\\ \left(\frac{\mathrm{\Lambda }\left(\mathrm{1}-p\right)}{\epsilon +\mu }\right.,\mathrm{0,0},\frac{\mathrm{\Lambda }\left(\epsilon +p\mu \right)}{\left(\epsilon +\mu \right)\left[b{\alpha }_{\mathrm{1}}\left(\mathrm{1}-w\right)+\mu \right]},\mathrm{0},\\ \frac{b{\alpha }_{\mathrm{1}}\mathrm{\Lambda }\left(\mathrm{1}-w\right)\left(\epsilon +p\mu \right)}{\left(\epsilon +\mu \right)\left({\alpha }_{\mathrm{2}}+\mu \right)\left[b{\alpha }_{\mathrm{1}}\left(\mathrm{1}-w\right)+\mu \right]},\\ \left.\frac{b{\alpha }_{\mathrm{1}}{\alpha }_{\mathrm{2}}\mathrm{\Lambda }\left(\mathrm{1}-w\right)\left(\epsilon +p\mu \right)}{\mu \left(\epsilon +\mu \right)\left({\alpha }_{\mathrm{2}}+\mu \right)\left[b{\alpha }_{\mathrm{1}}\left(\mathrm{1}-w\right)+\mu \right]}\right).\end{array}

基本再生数R_{\mathrm{0}}的计算过程如下：

F=\left[\begin{array}{cc}{\beta }_{\mathrm{1}}{S}_{\mathrm{10}}+w{\beta }_{\mathrm{1}}{V}_{\mathrm{0}}& {\beta }_{\mathrm{2}}{S}_{\mathrm{10}}+w{\beta }_{\mathrm{2}}{V}_{\mathrm{0}}\\ \mathrm{0}& \mathrm{0}\end{array}\right],

V=\left[\begin{array}{cc}\sigma +\mu & \mathrm{0}\\ -\sigma & \gamma +\mu \end{array}\right].

\begin{array}{l}{R}_{\mathrm{0}}=\rho \left(F{V}^{-\mathrm{1}}\right)=\\ \frac{{\beta }_{\mathrm{1}}{S}_{\mathrm{10}}+w{\beta }_{\mathrm{1}}{V}_{\mathrm{0}}}{\sigma +\mu }+\frac{\sigma \left({\beta }_{\mathrm{2}}{S}_{\mathrm{10}}+w{\beta }_{\mathrm{2}}{V}_{\mathrm{0}}\right)}{\left(\sigma +\mu \right)\left(\gamma +\mu \right)}=\\ \left(S_{\mathrm{10}}+w{V}_{\mathrm{0}}\right)\left(\frac{{\beta }_{\mathrm{1}}}{\sigma +\mu }+\frac{{\beta }_{\mathrm{2}}\sigma }{\left(\sigma +\mu \right)\left(\gamma +\mu \right)}\right).\end{array}

定理1 当R_{\mathrm{0}}>\mathrm{1}时，系统（1）在\mathrm{\Omega }内存在唯一的地方病平衡点E_{\mathrm{*}}=\left(S_{\mathrm{1}}^{\mathrm{*}},E^{\mathrm{*}},I^{\mathrm{*}},V^{\mathrm{*}},R_{\mathrm{2}}^{\mathrm{*}},R_{\mathrm{1}}^{\mathrm{*}},S_{\mathrm{2}}^{\mathrm{*}}\right).

证明 令系统（1）的方程右端均等于0，可得

I^{\mathrm{*}}=\frac{\sigma }{a}{E}^{\mathrm{*}},

{S_{\mathrm{1}}}^{\mathrm{*}}=\frac{\mathrm{\Lambda }\left(\mathrm{1}-p\right)}{A{E}^{\mathrm{*}}+\epsilon +\mu },

V^{\mathrm{*}}=\frac{\epsilon \mathrm{\Lambda }\left(\mathrm{1}-p\right)+p\mathrm{\Lambda }\left(\epsilon +\mu +A{E}^{\mathrm{*}}\right)}{\left[wA{E}^{\mathrm{*}}+\left(\mathrm{1}-w\right)b{\alpha }_{\mathrm{1}}+\mu \right]\left(A{E}^{\mathrm{*}}+\epsilon +\mu \right)},

R_{\mathrm{1}}^{\mathrm{*}}=\frac{b{\alpha }_{\mathrm{1}}\mu V^{\mathrm{*}}\left(\mathrm{1}-w\right)+{\alpha }_{\mathrm{1}}\gamma I^{\mathrm{*}}+b{{\alpha }_{\mathrm{1}}}^{\mathrm{2}}{V}^{\mathrm{*}}\left(\mathrm{1}-w\right)}{{\alpha }_{\mathrm{1}}{\alpha }_{\mathrm{2}}+{\alpha }_{\mathrm{1}}\mu +\mu gA{E}^{\mathrm{*}}+\mu {\alpha }_{\mathrm{2}}+{\mu }^{\mathrm{2}}},

R_{\mathrm{2}}^{\mathrm{*}}=\frac{\gamma I^{\mathrm{*}}+\left(\mathrm{1}-w\right)b{\alpha }_{\mathrm{1}}{V}^{\mathrm{*}}-{\alpha }_{\mathrm{2}}{R}_{\mathrm{1}}^{\mathrm{*}}-\mu R_{\mathrm{1}}^{\mathrm{*}}}{\mu },

{S_{\mathrm{2}}}^{\mathrm{*}}=\frac{{\alpha }_{\mathrm{2}}}{\mu }{R_{\mathrm{1}}}^{\mathrm{*}}.

其中，A{E}^{\mathrm{*}}={\beta }_{\mathrm{1}}{E}^{\mathrm{*}}+{\beta }_{\mathrm{2}}{I}^{\mathrm{*}}，a=\gamma +\mu，c=\sigma +\mu.

化简得到关于E^{\mathrm{*}}的一元二次方程a_{\mathrm{1}}{E^{\mathrm{*}}}^{\mathrm{2}}+a_{\mathrm{2}}{E}^{\mathrm{*}}+a_{\mathrm{3}}=\mathrm{0}，其中，

a_{\mathrm{1}}=-\left(\sigma +\mu \right)w{A}^{\mathrm{2}},

a_{\mathrm{2}}=A^{\mathrm{2}}\mathrm{\Lambda }w-\left[Acw\left(\epsilon +\mu \right)+Ab{\alpha }_{\mathrm{1}}c\left(\mathrm{1}-w\right)+Ac\mu \right],

\begin{array}{l}{a}_{\mathrm{3}}=A\mathrm{\Lambda }b{\alpha }_{\mathrm{1}}\left(\mathrm{1}-w\right)\left(\mathrm{1}-p\right)+\\ A\mathrm{\Lambda }\mu \left(\mathrm{1}-p\right)+A\mathrm{\Lambda }w\left(\epsilon +p\mu \right)-\\ \left[cb{\alpha }_{\mathrm{1}}\left(\epsilon +\mu \right)\left(\mathrm{1}-w\right)+\mu c\left(\epsilon +\mu \right)\right].\end{array}

当R_{\mathrm{0}}>\mathrm{1}时，有

a_{\mathrm{3}}>\mathrm{0},

故方程只有一个正根，所以系统（1）只存在一个地方病平衡点E_{\mathrm{*}}，定理得证.

3 稳定性分析

定理2 当时，系统（1）的无病平衡点E^{\mathrm{0}}在\mathrm{\Omega }内局部渐近稳定；当R_{\mathrm{0}}>\mathrm{1}时，无病平衡点E^{\mathrm{0}}在\mathrm{\Omega }内不稳定.

证明 系统（1）在E^{\mathrm{0}}处的Jacobian矩阵为：

\left[\begin{array}{ccccccc}-\epsilon -\mu & -{\beta }_{\mathrm{1}}{S}_{\mathrm{10}}& -{\beta }_{\mathrm{2}}{S}_{\mathrm{10}}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& {\beta }_{\mathrm{1}}{S}_{\mathrm{10}}+w{\beta }_{\mathrm{1}}{V}_{\mathrm{0}}-\sigma -\mu & {\beta }_{\mathrm{2}}{S}_{\mathrm{10}}+w{\beta }_{\mathrm{2}}{V}_{\mathrm{0}}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& \sigma & -\left(\gamma +\mu \right)& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}\\ \epsilon & -w{\beta }_{\mathrm{1}}{V}_{\mathrm{0}}& -w{\beta }_{\mathrm{2}}{V}_{\mathrm{0}}& -\left(\mathrm{1}-w\right)b{\alpha }_{\mathrm{1}}-\mu & \mathrm{0}& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& g{\beta }_{\mathrm{1}}{R}_{\mathrm{10}}& \gamma +g{\beta }_{\mathrm{2}}{R}_{\mathrm{10}}& \mathrm{0}& -\left({\alpha }_{\mathrm{1}}+\mu \right)& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& -g{\beta }_{\mathrm{1}}{R}_{\mathrm{10}}& -g{\beta }_{\mathrm{2}}{R}_{\mathrm{10}}& \left(\mathrm{1}-w\right)b{\alpha }_{\mathrm{1}}& {\alpha }_{\mathrm{1}}& -\left({\alpha }_{\mathrm{2}}+\mu \right)& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& {\alpha }_{\mathrm{2}}& -\mu \end{array}\right].

显然有-\left(\epsilon +\mu \right)、-\left(\mathrm{1}-w\right)b{\alpha }_{\mathrm{1}}-\mu、-\left({\alpha }_{\mathrm{1}}+\mu \right)、-\left({\alpha }_{\mathrm{2}}+\mu \right)、-\mu这5个负的特征值，另外的2个特征值由如下特征多项式确定：

\begin{array}{l}\left(\lambda -{\beta }_{\mathrm{1}}{S}_{\mathrm{10}}-w{\beta }_{\mathrm{1}}{V}_{\mathrm{0}}+\sigma +\mu \right)\left(\lambda +\gamma +\mu \right)-\\ \sigma \left({\beta }_{\mathrm{2}}{S}_{\mathrm{10}}+w{\beta }_{\mathrm{2}}{V}_{\mathrm{0}}\right)=\mathrm{0}.\end{array}（2）

由求根公式可得

{\lambda }_{\mathrm{1}}+{\lambda }_{\mathrm{2}}=-\left(\gamma +\mu -{\beta }_{\mathrm{1}}{S}_{\mathrm{10}}-w{\beta }_{\mathrm{1}}{V}_{\mathrm{0}}+\sigma +\mu \right),

\begin{array}{l}{\lambda }_{\mathrm{1}}{\lambda }_{\mathrm{2}}=\left(-{\beta }_{\mathrm{1}}{S}_{\mathrm{10}}-w{\beta }_{\mathrm{1}}{V}_{\mathrm{0}}+\sigma +\mu \right)\left(\gamma +\mu \right)-\\ \sigma \left({\beta }_{\mathrm{2}}{S}_{\mathrm{10}}+w{\beta }_{\mathrm{2}}{V}_{\mathrm{0}}\right).\end{array}

当时，{\lambda }_{\mathrm{1}}{\lambda }_{\mathrm{2}}>\mathrm{0}，，方程（2）所有根的实部均为负，无病平衡点E^{\mathrm{0}}局部渐近稳定.当R_{\mathrm{0}}>\mathrm{1}时，，即方程（2）有根具有正实部，平衡点E^{\mathrm{0}}不稳定，定理得证.

定理3 当时，系统（1）的无病平衡点E^{\mathrm{0}}在\mathrm{\Omega }内全局渐近稳定.

证明 由系统（1）可知\underset{t\to +\mathrm{\infty }}{\mathrm{l}\mathrm{i}\mathrm{m}}{S}_{\mathrm{1}}\left(t\right)=S_{\mathrm{10}}，\underset{t\to +\mathrm{\infty }}{\mathrm{l}\mathrm{i}\mathrm{m}}V\left(t\right)=V_{\mathrm{0}}.

构造Lyapunov函数

V_{\mathrm{1}}=\sigma \left(\gamma +\mu \right)E+\left(\sigma +\mu \right)I,

则V_{\mathrm{1}}沿着系统（1）的全导数为

\begin{array}{l}\frac{\mathrm{d}{V}_{\mathrm{1}}}{\mathrm{d}t}=\sigma \left(\sigma +\mu \right)\left(\gamma +\mu \right)\left(R_{\mathrm{0}}-\mathrm{1}\right)E+\\ \left[\sigma \left(\sigma +\mu \right)-{\sigma }^{\mathrm{2}}\left({\beta }_{\mathrm{2}}{S}_{\mathrm{10}}+w{\beta }_{\mathrm{2}}V\right)\right]E+\\ \left(\gamma +\mu \right)\left(\sigma +\mu \right)\left(\gamma +\mu \right)\left(R_{\mathrm{0}}-\mathrm{1}\right)I+\\ \left(\gamma +\mu \right)\left(\sigma +\mu \right)\left(\gamma +\mu \right)I+\end{array}

\begin{array}{l}\left[-\left({\beta }_{\mathrm{1}}{S}_{\mathrm{10}}+w{\beta }_{\mathrm{1}}{V}_{\mathrm{0}}\right)\left(\gamma +\mu \right)-\left(\sigma +\mu \right)\right]I\le \\ \sigma \left[\left(\sigma +\mu \right)\left(\gamma +\mu \right)\left(R_{\mathrm{0}}-\mathrm{1}\right)+\left(\sigma +\mu \right)\right]E+\\ \left(\gamma +\mu \right)\left(\sigma +\mu \right)\left(\gamma +\mu \right)\left(R_{\mathrm{0}}-\mathrm{1}\right)I+\\ \left(\sigma +\mu \right)\left(\gamma +\mu \right)I.\end{array}

当R_{\mathrm{0}}\le \mathrm{1}时，\frac{\mathrm{d}{V}_{\mathrm{1}}}{\mathrm{d}t}\le \mathrm{0}.由LaSalle不变集原理^［19］可知，方程（1）在\mathrm{\Omega }上的全部轨线趋向于集合\left\{\left(S_{\mathrm{1}},V,R_{\mathrm{1}},S_{\mathrm{2}}\right)\in \mathrm{\Omega }|E,I=\mathrm{0}\right\}的最大不变集，即无病平衡点全局稳定.定理得证.

定理4 当R_{\mathrm{0}}>\mathrm{1}时，系统（1）的地方病平衡点E_{\mathrm{*}}在\mathrm{\Omega }内全局渐近稳定.

证明 因为（1）中前4个方程与R_{\mathrm{2}}、R_{\mathrm{1}}、S_{\mathrm{2}}无关，所以只需要考虑其子方程的稳定性即可，则有

\left\{\begin{array}{l}\frac{\mathrm{d}{S}_{\mathrm{1}}}{\mathrm{d}t}=S_{\mathrm{1}}\left[\mathrm{\Lambda }\left(\mathrm{1}-p\right)\left(\frac{\mathrm{1}}{S_{\mathrm{1}}}-\frac{\mathrm{1}}{{S_{\mathrm{1}}}^{\mathrm{*}}}\right)-{\beta }_{\mathrm{1}}\left(E-E^{\mathrm{*}}\right)-{\beta }_{\mathrm{2}}\left(I-I^{\mathrm{*}}\right)\right],\\ \frac{\mathrm{d}E}{\mathrm{d}t}=E\left[{\beta }_{\mathrm{1}}\left(S_{\mathrm{1}}-{S_{\mathrm{1}}}^{\mathrm{*}}\right)+{\beta }_{\mathrm{2}}\left(\frac{S_{\mathrm{1}}I}{E}-\frac{{S_{\mathrm{1}}}^{\mathrm{*}}{I}^{\mathrm{*}}}{E^{\mathrm{*}}}\right)+w{\beta }_{\mathrm{1}}\left(V-V^{\mathrm{*}}\right)+w{\beta }_{\mathrm{2}}\left(\frac{VI}{E}-\frac{V^{\mathrm{*}}{I}^{\mathrm{*}}}{E^{\mathrm{*}}}\right)\right],\\ \frac{\mathrm{d}I}{\mathrm{d}t}=I\left(\sigma \frac{E}{I}-\sigma \frac{E^{\mathrm{*}}}{I^{\mathrm{*}}}\right),\\ \frac{\mathrm{d}V}{\mathrm{d}t}=V\left[p\mathrm{\Lambda }\left(\frac{\mathrm{1}}{V}-\frac{\mathrm{1}}{V^{\mathrm{*}}}\right)+\epsilon \left(\frac{S_{\mathrm{1}}}{V}-\frac{{S_{\mathrm{1}}}^{\mathrm{*}}}{V^{\mathrm{*}}}\right)-w{\beta }_{\mathrm{1}}\left(E-E^{\mathrm{*}}\right)-w{\beta }_{\mathrm{2}}\left(I-I^{\mathrm{*}}\right)\right].\end{array}\right.（3）

令S_{\mathrm{1}}=x{S_{\mathrm{1}}}^{\mathrm{*}}，E=y{E}^{\mathrm{*}}，I=z{I}^{\mathrm{*}}，V=u{V}^{\mathrm{*}}，代入（3）得

\left\{\begin{array}{l}\frac{\mathrm{d}{S}_{\mathrm{1}}}{\mathrm{d}t}=\mathrm{\Lambda }\left(\mathrm{1}-p\right)\left(\mathrm{1}-x\right)-{\beta }_{\mathrm{1}}{S_{\mathrm{1}}}^{\mathrm{*}}{E}^{\mathrm{*}}x\left(y-\mathrm{1}\right)-\\ {\beta }_{\mathrm{2}}{S_{\mathrm{1}}}^{\mathrm{*}}{I}^{\mathrm{*}}x\left(z-\mathrm{1}\right),\\ \frac{\mathrm{d}E}{\mathrm{d}t}={\beta }_{\mathrm{1}}{S_{\mathrm{1}}}^{\mathrm{*}}{E}^{\mathrm{*}}y\left(x-\mathrm{1}\right)+{\beta }_{\mathrm{2}}{S_{\mathrm{1}}}^{\mathrm{*}}{I}^{\mathrm{*}}\left(xz-y\right)+\\ w{\beta }_{\mathrm{1}}{E}^{\mathrm{*}}{V}^{\mathrm{*}}y\left(u-\mathrm{1}\right)+w{\beta }_{\mathrm{2}}{V}^{\mathrm{*}}{I}^{\mathrm{*}}\left(uz-y\right),\\ \frac{\mathrm{d}I}{\mathrm{d}t}=\sigma E^{\mathrm{*}}\left(y-z\right),\\ \frac{\mathrm{d}V}{\mathrm{d}t}=p\mathrm{\Lambda }\left(\mathrm{1}-u\right)+\epsilon {S_{\mathrm{1}}}^{\mathrm{*}}\left(x-u\right)-\\ w{\beta }_{\mathrm{1}}{V}^{\mathrm{*}}{E}^{\mathrm{*}}u\left(y-\mathrm{1}\right)-w{\beta }_{\mathrm{2}}{V}^{\mathrm{*}}{I}^{\mathrm{*}}u\left(z-\mathrm{1}\right).\end{array}\right.

构造Lyapunov函数

\begin{array}{l}{V}_{\mathrm{2}}=S_{\mathrm{1}}^{\mathrm{*}}\left(x-\mathrm{1}-\mathrm{l}\mathrm{n}x\right)+E^{\mathrm{*}}\left(y-\mathrm{1}-\mathrm{l}\mathrm{n}y\right)+\\ V^{\mathrm{*}}\left(u-\mathrm{1}-\mathrm{l}\mathrm{n}u\right)+t{I}^{\mathrm{*}}\left(z-\mathrm{1}-\mathrm{l}\mathrm{n}z\right).\end{array}

计算V_{\mathrm{2}}函数的导数得

\begin{array}{l}\frac{\mathrm{d}{V}_{\mathrm{2}}}{\mathrm{d}t}=\mathrm{\Lambda }\left(\mathrm{1}-p\right)\left(\mathrm{1}-x\right)\left(\mathrm{1}-\frac{\mathrm{1}}{x}\right)+\\ p\mathrm{\Lambda }\left(\mathrm{1}-u\right)\left(\mathrm{1}-\frac{\mathrm{1}}{u}\right)+\\ \epsilon {S_{\mathrm{1}}}^{\mathrm{*}}\left(x-u\right)\left(\mathrm{1}-\frac{\mathrm{1}}{u}\right)+t\sigma E^{\mathrm{*}}+x{\beta }_{\mathrm{2}}{S_{\mathrm{1}}}^{\mathrm{*}}{I}^{\mathrm{*}}-\\ y\left({\beta }_{\mathrm{2}}{S_{\mathrm{1}}}^{\mathrm{*}}{I}^{\mathrm{*}}+w{\beta }_{\mathrm{2}}{V}^{\mathrm{*}}{I}^{\mathrm{*}}-t\sigma E^{\mathrm{*}}\right)+\\ z\left({\beta }_{\mathrm{2}}{S_{\mathrm{1}}}^{\mathrm{*}}{I}^{\mathrm{*}}+w{\beta }_{\mathrm{2}}{V}^{\mathrm{*}}{I}^{\mathrm{*}}-t\sigma E^{\mathrm{*}}\right)+\\ uw{\beta }_{\mathrm{2}}{V}^{\mathrm{*}}{I}^{\mathrm{*}}-{\beta }_{\mathrm{2}}{S_{\mathrm{1}}}^{\mathrm{*}}{I}^{\mathrm{*}}\frac{xz}{y}-w{\beta }_{\mathrm{2}}{V}^{\mathrm{*}}{I}^{\mathrm{*}}\frac{uz}{y}-t\sigma E^{\mathrm{*}}\frac{y}{z}.\end{array}

令

t=\frac{{\beta }_{\mathrm{2}}{S_{\mathrm{1}}}^{\mathrm{*}}{I}^{\mathrm{*}}+w{\beta }_{\mathrm{2}}{V}^{\mathrm{*}}{I}^{\mathrm{*}}}{\sigma E^{\mathrm{*}}},

故

\begin{array}{l}\frac{\mathrm{d}{V}_{\mathrm{2}}}{\mathrm{d}t}=\mathrm{\Lambda }\left(\mathrm{1}-p\right)\left(\mathrm{1}-x\right)\left(\mathrm{1}-\frac{\mathrm{1}}{x}\right)+\\ p\mathrm{\Lambda }\left(\mathrm{1}-u\right)\left(\mathrm{1}-\frac{\mathrm{1}}{u}\right)+\\ {\beta }_{\mathrm{2}}{S_{\mathrm{1}}}^{\mathrm{*}}{I}^{\mathrm{*}}+w{\beta }_{\mathrm{2}}{V}^{\mathrm{*}}{I}^{\mathrm{*}}+\epsilon {S_{\mathrm{1}}}^{\mathrm{*}}\left(x-u\right)\left(\mathrm{1}-\frac{\mathrm{1}}{u}\right)+\\ x{\beta }_{\mathrm{2}}{S_{\mathrm{1}}}^{\mathrm{*}}{I}^{\mathrm{*}}+uw{\beta }_{\mathrm{2}}{V}^{\mathrm{*}}{I}^{\mathrm{*}}-{\beta }_{\mathrm{2}}{S_{\mathrm{1}}}^{\mathrm{*}}{I}^{\mathrm{*}}\frac{xz}{y}-\\ w{\beta }_{\mathrm{2}}{V}^{\mathrm{*}}{I}^{\mathrm{*}}\frac{uz}{y}-\left({\beta }_{\mathrm{2}}{S_{\mathrm{1}}}^{\mathrm{*}}{I}^{\mathrm{*}}+w{\beta }_{\mathrm{2}}{V}^{\mathrm{*}}{I}^{\mathrm{*}}\right)\frac{y}{z}.\end{array}

根据文献［20］中的方法，令

\left\{\begin{array}{l}{B}_{\mathrm{1}}=b_{\mathrm{1}}\left(\mathrm{2}-x-\frac{\mathrm{1}}{x}\right),\\ B_{\mathrm{2}}=b_{\mathrm{2}}\left(\mathrm{2}-u-\frac{\mathrm{1}}{u}\right),\\ B_{\mathrm{3}}=b_{\mathrm{3}}\left(\mathrm{3}-\frac{\mathrm{1}}{x}-\frac{y}{z}-\frac{xz}{y}\right),\\ B_{\mathrm{4}}=b_{\mathrm{4}}\left(\mathrm{3}-\frac{\mathrm{1}}{u}-\frac{y}{z}-\frac{uz}{y}\right),\\ B_{\mathrm{5}}=b_{\mathrm{5}}\left(\mathrm{3}-u-\frac{\mathrm{1}}{x}-\frac{x}{u}\right),\\ B_{\mathrm{6}}=b_{\mathrm{6}}\left(\mathrm{4}-\frac{\mathrm{1}}{x}-\frac{x}{u}-\frac{uz}{y}-\frac{y}{z}\right).\end{array}\right.

根据对应系数相等可得

\left\{\begin{array}{l}{b}_{\mathrm{1}}=\mathrm{\Lambda }\left(\mathrm{1}-p\right)-{\beta }_{\mathrm{2}}{S_{\mathrm{1}}}^{\mathrm{*}}{I}^{\mathrm{*}}-\epsilon {S_{\mathrm{1}}}^{\mathrm{*}}={\beta }_{\mathrm{1}}{S_{\mathrm{1}}}^{\mathrm{*}}{E}^{\mathrm{*}}+\mu {S_{\mathrm{1}}}^{\mathrm{*}},\\ b_{\mathrm{2}}=p\mathrm{\Lambda }-w{\beta }_{\mathrm{2}}{V}^{\mathrm{*}}{I}^{\mathrm{*}}+b_{\mathrm{6}},\\ b_{\mathrm{3}}={\beta }_{\mathrm{2}}{S_{\mathrm{1}}}^{\mathrm{*}}{I}^{\mathrm{*}},\\ b_{\mathrm{4}}=w{\beta }_{\mathrm{2}}{V}^{\mathrm{*}}{I}^{\mathrm{*}}-b_{\mathrm{6}},\\ b_{\mathrm{5}}=\epsilon {S_{\mathrm{1}}}^{\mathrm{*}}-b_{\mathrm{6}}.\end{array}\right.

当系数非负的时候，要求

-p\mathrm{\Lambda }+w{\beta }_{\mathrm{2}}{V}^{\mathrm{*}}{I}^{\mathrm{*}}\le b_{\mathrm{6}}\le \mathrm{m}\mathrm{i}\mathrm{n}\left\{w{\beta }_{\mathrm{2}}{V}^{\mathrm{*}}{I}^{\mathrm{*}},\epsilon S_{\mathrm{1}}^{\mathrm{*}}\right\}.

1） 当p\mathrm{\Lambda }\ge w{\beta }_{\mathrm{2}}{V}^{*}{I}^{*}时，取b_{\mathrm{6}}=\mathrm{0}，则

\begin{array}{l}\frac{\mathrm{d}{V}_{\mathrm{2}}}{\mathrm{d}t}=\left({\beta }_{\mathrm{1}}{S_{\mathrm{1}}}^{\mathrm{*}}{E}^{\mathrm{*}}+\mu {S_{\mathrm{1}}}^{\mathrm{*}}\right)\left(\mathrm{2}-x-\frac{\mathrm{1}}{x}\right)+\\ \left(p\mathrm{\Lambda }-w{\beta }_{\mathrm{2}}{V}^{\mathrm{*}}{I}^{\mathrm{*}}\right)\left(\mathrm{2}-u-\frac{\mathrm{1}}{u}\right)+\\ {\beta }_{\mathrm{2}}{S_{\mathrm{1}}}^{\mathrm{*}}{I}^{\mathrm{*}}\left(\mathrm{3}-\frac{\mathrm{1}}{x}-\frac{y}{z}-\frac{xz}{y}\right)+\\ w{\beta }_{\mathrm{2}}{V}^{\mathrm{*}}{I}^{\mathrm{*}}\left(\mathrm{3}-\frac{\mathrm{1}}{u}-\frac{y}{z}-\frac{uz}{y}\right)+\\ \epsilon {S_{\mathrm{1}}}^{\mathrm{*}}\left(\mathrm{3}-u-\frac{\mathrm{1}}{x}-\frac{x}{u}\right)\le \mathrm{0}.\end{array}

2） 当p\mathrm{\Lambda }\le w{\beta }_{\mathrm{2}}{V}^{*}{I}^{*}时，取b_{\mathrm{6}}=w{\beta }_{\mathrm{2}}{V}^{*}{I}^{*}-p\mathrm{\Lambda }，则

\begin{array}{l}\frac{\mathrm{d}{V}_{\mathrm{2}}}{\mathrm{d}t}=\left({\beta }_{\mathrm{1}}{S_{\mathrm{1}}}^{\mathrm{*}}{E}^{\mathrm{*}}+\mu {S_{\mathrm{1}}}^{\mathrm{*}}\right)\left(\mathrm{2}-x-\frac{\mathrm{1}}{x}\right)+\\ {\beta }_{\mathrm{2}}{S_{\mathrm{1}}}^{\mathrm{*}}{I}^{\mathrm{*}}\left(\mathrm{3}-\frac{\mathrm{1}}{x}-\frac{y}{z}-\frac{xz}{y}\right)+\\ p\mathrm{\Lambda }\left(\mathrm{3}-\frac{\mathrm{1}}{u}-\frac{y}{z}-\frac{uz}{y}\right)+\\ \left(\epsilon {S_{\mathrm{1}}}^{\mathrm{*}}-w{\beta }_{\mathrm{2}}{V}^{\mathrm{*}}{I}^{\mathrm{*}}+p\mathrm{\Lambda }\right)\left(\mathrm{3}-u-\frac{\mathrm{1}}{x}-\frac{x}{u}\right)+\end{array}

\left(w{\beta }_{\mathrm{2}}{V}^{\mathrm{*}}{I}^{\mathrm{*}}-p\mathrm{\Lambda }\right)\left(\mathrm{4}-\frac{\mathrm{1}}{x}-\frac{x}{u}-\frac{uz}{y}-\frac{y}{z}\right)\le \mathrm{0}.

当\frac{\mathrm{d}{V}_{\mathrm{2}}}{\mathrm{d}t}=\mathrm{0}时，集合上的不变集是地方病平衡点，所以由LaSalle不变集原理可知，地方病平衡点E_{\mathrm{*}}在\mathrm{\Omega }中全局渐近稳定.定理得证.

4 数值模拟

为了证明平衡点的稳定性，选择参数值为\mathrm{\Lambda }=\mathrm{10}，{\beta }_{\mathrm{1}}=\mathrm{0.003}，{\beta }_{\mathrm{2}}=\mathrm{0.002}，p=\mathrm{0.85}^［21］，w=\mathrm{0.2}，\mu =\mathrm{0.014}，\epsilon =\mathrm{0.85}，b=\mathrm{2}^［21］，g=\mathrm{0.5}，{\alpha }_{\mathrm{1}}=\mathrm{0.09}^［22］，\sigma =\mathrm{0.2}，{\alpha }_{\mathrm{2}}=\mathrm{0.16}，\gamma =\mathrm{0.06}，通过计算可得，从图2a中可以看出E\left(t\right)\mathrm{、}I\left(t\right)的数量趋近于0，则疾病最终会灭绝，且参数值满足定理3的条件，说明E^{\mathrm{0}}全局渐近稳定.

选取参数值为\mathrm{\Lambda }=\mathrm{10}，{\beta }_{\mathrm{1}}=\mathrm{0.008}，{\beta }_{\mathrm{2}}=\mathrm{0.002}，p=\mathrm{0.85}，w=\mathrm{0.54}，\mu =\mathrm{0.014}，\epsilon =\mathrm{0.85}，b=\mathrm{2}，g=\mathrm{0.5}，{\alpha }_{\mathrm{1}}=\mathrm{0.05}，\sigma =\mathrm{0.2}，{\alpha }_{\mathrm{2}}=\mathrm{0.16}，\gamma =\mathrm{0.01}，通过计算可得R_{\mathrm{0}}>\mathrm{1}.从图2b中可以看出E\left(t\right)\mathrm{、}I\left(t\right)最后的数量趋于一个稳定值，意味着疾病将一直持续，且参数值满足定理4的条件，则E_{\mathrm{*}}全局渐近稳定.

5 敏感性分析

本节采用阈值进行全局敏感性分析来衡量模型主导因素变化的效果，并指出模型中对疾病流行产生最大影响的参数^［23-24］.敏感性分析的结果为制定有效的控制策略提供了依据.变量R_{\mathrm{0}}对于参数是可微的，其灵敏度指数\mathrm{\Pi }定义为^［25］：

{\mathrm{\Pi }}_{\mathrm{\Phi }}^{R_{\mathrm{0}}}=\frac{\partial R_{\mathrm{0}}}{\partial \mathrm{\Phi }}\times \frac{\mathrm{\Phi }}{R_{\mathrm{0}}}，

式中，\mathrm{\Phi }表示所有基本参数.

例如，R_{\mathrm{0}}相对于基本参数{\beta }_{\mathrm{1}}的灵敏度指数为：

\begin{array}{l}{\mathrm{\Pi }}_{{\beta }_{\mathrm{1}}}^{R_{\mathrm{0}}}=\frac{\partial R_{\mathrm{0}}}{\partial {\beta }_{\mathrm{1}}}\times \frac{{\beta }_{\mathrm{1}}}{R_{\mathrm{0}}}=\\ \frac{\mathrm{\Lambda }\left[b{\alpha }_{\mathrm{1}}\left(\mathrm{1}-w\right)\left(\mathrm{1}-p\right)+\mu \left(\mathrm{1}-p\right)+w\left(\epsilon +p\mu \right)\right]}{\left(\sigma +\mu \right)\left(\epsilon +\mu \right)\left[b{\alpha }_{\mathrm{1}}\left(\mathrm{1}-w\right)+\mu \right]}>\\ \mathrm{0}.\end{array}

以相同的方法得到{\mathrm{\Pi }}_{{\beta }_{\mathrm{2}}}^{R_{\mathrm{0}}}、{\mathrm{\Pi }}_p^{R_{\mathrm{0}}}、{\mathrm{\Pi }}_w^{R_{\mathrm{0}}}、{\mathrm{\Pi }}_{\mu }^{R_{\mathrm{0}}}、{\mathrm{\Pi }}_{\epsilon }^{R_{\mathrm{0}}}、{\mathrm{\Pi }}_{{\alpha }_{\mathrm{1}}}^{R_{\mathrm{0}}}、{\mathrm{\Pi }}_{\sigma }^{R_{\mathrm{0}}}和{\mathrm{\Pi }}_{\gamma }^{R_{\mathrm{0}}}，其灵敏度指标如表1所示.

结合图3可明显看出，{\beta }_{\mathrm{2}}和w对R_{\mathrm{0}}有正向影响.这表明，假设这些参数增加10%，基本再生数R_{\mathrm{0}}将分别增加5.4%和10.8%.\gamma和{\alpha }_{\mathrm{1}}每增加10%，基本再生数R_{\mathrm{0}}的值将分别下降6.7%和8%，对R_{\mathrm{0}}有负向影响.

图4a模拟了有效接触率{\beta }_{\mathrm{2}}取不同值时，总感染人数E\left(t\right)+I\left(t\right)随时间变化的曲线.取与图2相同的参数，通过分析与染病者不同程度的接触对疾病传播的影响，很明显，当接触率降低时，可以显著降低总感染人数.

图4b模拟了染病者的恢复率\gamma对疾病感染人数的影响，即患者的恢复速度越快，基本再生数R_{\mathrm{0}}就会越低.当到达某一个特定值以后，E\left(t\right)+I\left(t\right)的人数将会趋于0，疾病最终会得到控制.

6 参数耦合关系分析

因为V_{\mathrm{0}}同时依赖w和{\alpha }_{\mathrm{1}}，

V_{\mathrm{0}}\left(w,{\alpha }_{\mathrm{1}}\right)=\frac{C}{b{\alpha }_{\mathrm{1}}\left(\mathrm{1}-w\right)+\mu }，

其中，C=\epsilon S_{\mathrm{10}}+p\mathrm{\Lambda }，故w和{\alpha }_{\mathrm{1}}对R_{\mathrm{0}}的影响应联合分析，

\begin{array}{l}{R}_{\mathrm{0}}\left(w,{\alpha }_{\mathrm{1}}\right)=\\ \left(S_{\mathrm{10}}+w{V}_{\mathrm{0}}\right)\left(\frac{{\beta }_{\mathrm{1}}}{\sigma +\mu }+\frac{{\beta }_{\mathrm{2}}\sigma }{\left(\sigma +\mu \right)\left(\gamma +\mu \right)}\right)=\\ \left[\frac{\left(\mathrm{1}-p\right)\mathrm{\Lambda }}{\epsilon +\mu }+w\frac{\frac{\left(\mathrm{1}-p\right)\mathrm{\Lambda }}{\epsilon +\mu }+p\mathrm{\Lambda }}{b{\alpha }_{\mathrm{1}}\left(\mathrm{1}-w\right)+\mu }\right]K,\end{array}

其中，K=\frac{{\beta }_{\mathrm{1}}}{\sigma +\mu }+\frac{{\beta }_{\mathrm{2}}\sigma }{\left(\sigma +\mu \right)\left(\gamma +\mu \right)}.

联合影响体现在乘积项w{V}_{\mathrm{0}}\left(w,{\alpha }_{\mathrm{1}}\right)中.

1） 对w求偏导：

\frac{\partial V_{\mathrm{0}}}{\partial w}=\frac{b{\alpha }_{\mathrm{1}}C}{{\left[b{\alpha }_{\mathrm{1}}\left(\mathrm{1}-w\right)+\mu \right]}^{\mathrm{2}}}>\mathrm{0}，

\frac{\partial R_{\mathrm{0}}}{\partial w}=\left[V_{\mathrm{0}}+w\frac{b{\alpha }_{\mathrm{1}}C}{{\left[b{\alpha }_{\mathrm{1}}\left(\mathrm{1}-w\right)+\mu \right]}^{\mathrm{2}}}\right]K>\mathrm{0}.

w增大，分母b{\alpha }_{\mathrm{1}}\left(\mathrm{1}-w\right)+\mu减小，V_{\mathrm{0}}增大.即疫苗失效率提高会增加疫苗接种者数量，从而增大R_{\mathrm{0}}.

由图5可得疫苗失效性越低，疾病的峰值越低，同时疾病达到峰值的时间也在逐渐推后.当疫苗失效率较高时，更多的人会被感染，从而导致疾病扩散.通过加大疫苗研发与改进，降低失效率，能够显著减少疫苗无效的情况，从而减少病毒传播的风险，达到控制疾病蔓延的效果.

2） 对{\alpha }_{\mathrm{1}}求偏导：

，

.

在模型中疫苗接种者V的流失率包含两项：免疫衰弱流失\left(\mathrm{1}-w\right)b{\alpha }_{\mathrm{1}}V和自然流失\mu V.当{\alpha }_{\mathrm{1}}减低时，\left(\mathrm{1}-w\right)b{\alpha }_{\mathrm{1}}{V}_{\mathrm{0}}减小，V_{\mathrm{0}}的留存率提高，接种者数量V_{\mathrm{0}}增加.在疫苗保护不完全时，V_{\mathrm{0}}增加会通过w{V}_{\mathrm{0}}提高R_{\mathrm{0}}.结果如图6所示.

若免疫力衰弱适度，接种者转入免疫状态R_{\mathrm{1}}，减少V中存量.虽然R_{\mathrm{1}}会逐渐失去免疫力，但短期内降低了可感染人群，从而降低了传播风险.故单研究疫苗的有效性不足以控制疫情，必须配合疫苗免疫持续时间.

3） w与{\alpha }_{\mathrm{1}}的联合效应

基于图2的数据，引入拉格朗日函数

\begin{array}{l}L\left(w,{\alpha }_{\mathrm{1}},\lambda \right)=\mathrm{0.83}+\mathrm{50}(\mathrm{1}-w)+\mathrm{100}{\alpha }_{\mathrm{1}}+\\ \lambda \left(w+\mathrm{1.8}{\alpha }_{\mathrm{1}}-\mathrm{0.35}\right),\end{array}

得到最优解

\left\{\begin{array}{c}{w}^{\mathrm{*}}=\mathrm{0.15},\\ {{\alpha }_{\mathrm{1}}}^{\mathrm{*}}=\mathrm{0.11},\\ R_{\mathrm{0}}^{\mathrm{*}}=\mathrm{0.92}.\end{array}\right.

通过图7中的对比可得，在传染病防控中，疫苗失效率与免疫力衰退率的协同配合展现出显著的双重抑制效果.当疫苗有效率提升至85%，可直接将接种人群的感染风险降低约72%，同时通过群体免疫效应间接减少病毒传播链；而将衰退率提高11%则能降低V_{\mathrm{0}}中人数，使平均传染期缩短.二者的最优组合可形成“预防—治疗”闭环.

综合上述过程，当w^{\mathrm{*}}=\mathrm{0.15}与{{\alpha }_{\mathrm{1}}}^{\mathrm{*}}=\mathrm{0.11}配合时，可将疫情传播周期缩短58%，且二次感染风险降低72%，这种协同作用比单独优化任何一个参数的效果提升2.3倍.值得注意的是，过低的免疫力衰退率，反而会导致免疫系统刺激不足，这与现代免疫学中的“适度抗原刺激”理论相吻合.该参数组合为制定差异化的疫苗接种策略和免疫加强方案提供了量化依据.

7 结论

疫苗接种虽已成为防控疾病的最重要手段，但是疫苗的有效性仍是人们需要攻克的难关.由此，本文建立了一个{\mathrm{S}}_{\mathrm{1}}\mathrm{E}\mathrm{I}\mathrm{V}{\mathrm{R}}_{\mathrm{2}}{\mathrm{R}}_{\mathrm{1}}{\mathrm{S}}_{\mathrm{2}}动力学模型，考虑了在免疫力变化的情况下，疫苗有效性对疾病传播的影响.

本文首先对模型进行定性分析，求解了基本再生数R_{\mathrm{0}}，并进一步分析了平衡点的稳定性.接着通过数值模拟和分析，重点考察了疫苗失效率与免疫力衰退率对疾病传播的协同影响.研究结果表明，这两个关键参数的耦合作用会显著改变疾病的传播动力学特征.

参考文献

[1] SAFAN M， KRETZSCHMAR M， HADELER K P. Vaccination based control of infections in SIRS models with reinfection： special reference to pertussis［J］. Journal of Mathematical Biology，2013，67（5）：1083-1110.

[2] GU X X， PLOTKIN S A， EDWARDS K M， et al. Waning immunity and microbial vaccines-workshop of the National Institute of Allergy and Infectious Diseases［J/OL］. Clinical and Vaccine Immunology， 2017，24（7） ［2025-09-16］.

[3] MOOI F R， VAN D M N A T， DE MELKER H E. Pertussis resurgence： waning immunity and pathogen adaptation - two sides of the same coin［J］. Epidemiology & Infection， 2014， 142（4）：685-694.

[4] WENDELBOE A M， RIE A V， SALMASO S， et al. Duration of immunity against pertussis after natural infection or vaccination［J］. Pediatric Infectious Disease Journal， 2005， 24（5）：58-61.

[5] CHAN K H， TO K K W， HUNG I F N， et al. Differences in antibody responses of individuals with natural infection and those vaccinated against pandemic H1N1 2009 influenza［J］. Clinical and Vaccine Immunology， 2011， 18（5）：867-873.

[6] HILTON C， WHITTLE. Effect of subclinical infection on maintaining immunity against measles in vaccinated children in West Africa［J］. Lancet，1999，353（9147）：98-102.

[7] HELENE B， SIMONDON F， ROHANI P， et al. Loss of immunity to pertussis in a rural community in Senegal［J］. Vaccine， 2004， 22（5-6）：594-596.

[8] HEFFERNAN J M， KEELING M J. Implications of vaccination and waning immunity［J］. Biological Sciences， 2009， 276（1664）： 2071-2080.

[9] LAVINE J S， KING A， BJORNSTAD O N. Natural immune boosting in pertussis dynamics and the potential for long-term vaccine failure［J］. Proceedings of the National Academy of Sciences of the United States of America， 2011，108（17）：7259-7264.

[10] CARLSSON R M， CHILDS L M， FENG Z， et al. Modeling the waning and boosting of immunity from infection or vaccination［J/OL］. Journal of Theoretical Biology，2020，497［2025-09-16］.

[11] GASPERJAN， SONA KILIANOVA， EHRHARDT M. SIR-based mathematical modeling of infectious diseases with vaccination and waning immunity［J/OL］. Journal of Computational Science，2019，37［2025-09-16］.

[12] ELBASHA E H， PODDER C N， GUMEL A B. Analyzing the dynamics of an SIRS vaccination model with waning natural and vaccine-induced immunity［J］. Nonlinear Analysis Real World Applications，2011，12（5）：2692-2705.

[13] OPOKU-SARKODIE R， BARTHA F A， POLNER M， et al. Dynamics of an SIRWS model with waning of immunity and varying immune boosting period［J］. Journal of Biological Dynamics， 2022，16（1）：596-618.

[14] 高杰，李毅红，罗晓峰. 考虑疫苗有效性的传染病渗流模型［J］.中北大学学报（自然科学版），2025，46（2）：254-260.

[15] GAO J， LI Y H， LUO X F. Infectious disease percolation model considering vaccine efficacy［J］. Journal of North University of China （Natural Science Edition）， 2025， 46（2）： 254-260. （Ch）.

[16] CALVO E A， URRUTIA T G A， ZAMBRANO V C J， et al. Pertussis vaccine effectiveness following country-wide implementation of a hexavalent acellular pertussis immunization schedule in infants and children in Panama［J/OL］. Human Vaccines，2024， 20 （1）［2025-09-16］.

[17] TAN C， XIAO Y， CHEN S， et al. Bibliometrics analysis and knowledge mapping of pertussis vaccine research： trends from 1994 to 2023［J］. Infection. 2025， 53（3）：1001-1012.

[18] BLANCHARD E， CHAVADE D， WAZIERES D B， et al. Pertussis vaccination in adults in France： overview and suggestions for improvement［J/OL］. Infectious Diseases Now， 2024， 54 （6） ［2025-09-16］.

[19] PHAM N， LE N， OLSON L， et al. Epidemiology and effective pertussis vaccination among hospitalized children in Vietnam National Children’s Hospital， 2015–2018［J/OL］. International Journal of Infectious Diseases， 2020， 101［2025-09-16］ .

[20] 廖晓昕.稳定性的理论，方法和应用［M］.武汉：华中理工大学出版社，1999：50-78.

[21] LIAO X X. Theory， methods， and applications of stability［M］. Wuhan： Huazhong University of Science and Technology Press， 1999：50-78. （Ch）.

[22] LI J， YANG Y， ZHOU Y. Global stability of an epidemic model with latent stage and vaccination［J］. Nonlinear Analysis Real World Applications， 2011，12（4）：2163-2173.

[23] RACHLIN A， DANOVARO-HOLLIDAY M C， WALLACE S A S. Routine vaccination coverage-worldwide， 2021［J］. Morbidity and Mortality Weekly Report： CDC Surveillance Summaries， 2022，71（44）：1396-1400.

[24] WEARING H J， ROHANI P. Estimating the duration of pertussis immunity using epidemiological signatures［J/OL］. PLOS Pathogens， 2009， 5（10） ［2025-09-16］.

[25] YOUSSEF H M， NAJAT A， MAGDY A E， et al. A proposed modified SEIQR epidemic model to analyze the COVID-19 spreading in Saudi Arabia［J］. Alexandria Engineering Journal， 2021，61（3）：2456-2470.

[26] 张轶菲，薛亚奎.一类布鲁氏菌病的多仓室传播动力学模型分析和应用研究［J］.华中师范大学学报（自然科学版），2025，59（2）：179-187.

[27] ZHANG Y F， XUE Y K. Analysis and application of a multi-compartment transmission dynamics model for brucellosis［J］. Journal of Central China Normal University （Natural Sciences）， 2025， 59（2）： 179-187. （Ch）.

[28] HABENOM H， AYCHLUH M， SUTHAR D L， et al. Modeling and analysis on the transmission of COVID-19 Pandemic in Ethiopia［J］. Alexandria Engineering Journal， 2021， 61：5323-5342.

基金资助

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

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