一类具有疫苗接种行为决策函数的传染病模型的敏感性分析

Sensitivity Analysis of Infectious Disease Models with Vaccination Behavior Decision Functions

摘要: 文章构建了一个融合了疫苗决策功能的SEIRV动力学模型。考虑基本的疾病传播过程的同时，引入疫苗接种决策变量，计算了模型边界平衡点和内部平衡点以及基本再生数和控制再生数，验证了不同条件下边界平衡点的稳定性。利用拉丁超立方体抽样(LHS)与偏秩相关系数(PRCC)方法分析模型部分参数对疫情规模的影响程度，观察到提高疫苗有效性、延长疫苗保护时间以及采取更严格的措施可以有效抑制病毒的传播，并显著缩短其持续时间。

Abstract: This paper presents the construction of a SEIRV dynamical model that incorporates a vaccine decision-making functionality. In addition to considering the fundamental processes of disease transmission, the model introduces a vaccination decision variable. We compute the model boundary, internal equilibrium points, basic reproduction numbers, and control reproduction numbers. The stability of the boundary equilibrium points under various conditions is rigorously verified. Utilizing Latin Hypercube Sampling (LHS) in conjunction with the Partial Rank Correlation Coefficient (PRCC) method, we analyze the impact of certain model parameters on the scale of the epidemic. Our findings reveal that increasing vaccine efficacy, extending vaccine protection duration, and implementing stricter measures can effectively mitigate viral transmission and substantially reduce its duration.

1. 绪论

自1918年的西班牙流感大流行，至2003年的严重急性呼吸综合征(SARS)疫情、2009年的甲型H1N1流感疫情，再到2019年的COVID-19疫情，我们接连遭遇了多种新发流行病。这一系列事件凸显了一个事实，即全球大流行病始终构成一种持续威胁，危及现代社会体系的安全与稳定。在此背景下，疫苗接种已成为预防和控制传染病最具成本效益和效率的公共卫生干预手段。对于家庭而言，疫苗接种是减少家庭成员患病概率及后续医疗费用支出的关键工具。

依据Kermack与McKendrick所提出的流行病学理论框架[1]，传染病的传播机制可通过构建隔离模型来进行数学化描述。因此，众多学者建议对SEIR模型进行拓展，纳入隔离、死亡以及潜在感染隔离等额外模块[2]。2006年，加拿大研究团队在探索SARS病毒疫苗的过程中，首次引入了包含疫苗仓室V的传染病模型。随后，该团队计算了模型在可行域上的平衡点，并给出了基本再生数以探讨其相关特性[3]。

在构建疾病传播模型的过程中，一个不可忽视的环节是不仅要精确模拟疾病的动态传播过程，还要融入人类决策行为的社会动态因素，包括保持社交距离、佩戴口罩、自我隔离在内的多种干预措施。博弈论作为一种重要的理论工具，在疫苗接种决策的制定过程中得到了广泛应用。在疫苗博弈的建模过程中，采用非线性响应函数来描述疫苗接种决策过程更为符合实际情况。这是因为人类的决策行为往往受到多种因素的影响，这些因素之间相互作用、相互制约，形成了一个复杂的非线性系统[4] [5]。因此，在本文中，我们构建了一个包含非线性响应函数的疫苗博弈数学模型，以期更准确地刻画个体在疫苗接种决策中的行为特征。

2. 传染病数学模型的建立

2.1. 带有疫苗接种行为的传染病模型

我们的模型基于SEIRV过程，适用于描述麻疹、流感及COVID-19等传染病的传播动力学。在此模型中，总人口被划分为五个不同的群组：易感人群(S)、暴露人群(E)、已接种人群(V)、感染人群(I)和具有自然免疫力的康复人群(R)。易感个体在与感染者有效接触后，会转变为暴露人群。暴露个体在病毒潜伏期结束后会转变为感染人群。感染个体最终会因自然康复、隔离或医疗干预而进入康复类别。值得注意的是，对于COVID-19和流感等疾病，近期接种疫苗或自然感染后的个体在一段时间内通常对再次感染具有免疫力。然而，由于疫苗效力的减弱，已接种个体在一段时间后可能会重新转变为易感人群。其中状态变量的动态变化通过以下微分方程组进行数学描述：

$\begin{array}{l}{S}^{\prime }\left(t\right)=\mu N-\left(\beta SI/N\right)-\sigma xS+\phi V-\mu S\\ E^{\prime }\left(t\right)=\left(\beta SI/N\right)-ϵE-\mu E\\ I^{\prime }\left(t\right)=ϵE-\gamma I-\mu I\\ R^{\prime }\left(t\right)=\gamma I-\mu R\\ V^{\prime }\left(t\right)=\sigma xS-\phi V-\mu V\end{array}$ (1)

该方程组具有非负初始条件$\left(S\left(0\right),E\left(0\right),I\left(0\right),R\left(0\right),V\left(0\right)\right)\in {ℝ}_{+}^5$ ，将模型(1)中的所有方程相加，得到

$\frac{\text{d}}{\text{d}t}\left(S+E+I+R+V\right)\le \mu N-\mu \left(S+E+I+R+V\right).$ (2)

这意味着

$\underset{t\to \infty }{\lim \sup }\left(S+E+I+R+V\right)\le N.$

我们将在如下的后续可行域内对模型(1)进行分析：

$\Gamma =\left\{\left(S,E,I,R,V\right)\in {ℝ}_{+}^5|0\le S+E+I+R+V\le N\right\}.$ (3)

模型中涉及的参数具体解释如下：$\mu$ 代表人均自然死亡率；$\varphi$ 代表疫苗效力随时间衰减的比率；$ϵ$ 代表个体从暴露状态向感染状态转化的速率；$\sigma$ 代表疫苗的有效率，其取值范围限定在$\left[0,1\right]$ 之间；$\beta$ 代表有症状病例向易感人群的传播速率；$\gamma$ 代表有症状病例的恢复速率；x代表接种疫苗的比例，其取值范围同样限定在$\left[0,1\right]$ 之间。值得注意的是，x会随着接种行为的变化而动态调整，这构成了一个复杂的动态决策博弈过程。在给定模型中，我们假设人口的出生率和死亡率保持平衡，且未将疾病导致的死亡率纳入考量。

因此，我们设定人口规模N为一个恒定值，且满足$N=S+E+V+I+R$ 。为了简化分析，我们引入以下比例变量：

$s=\frac{S}{N},e=\frac{E}{N},i=\frac{I}{N},r=\frac{R}{N},v=\frac{V}{N}$ .

显然，这些比例变量的和满足$1=s+e+v+i+r$ 。基于这些比例变量，我们可以将原微分方程系统(1)转化为以下形式：

$\begin{array}{l}{s}^{\prime }\left(t\right)=\mu -\beta si-\sigma xs-\mu s+\varphi v,\\ e^{\prime }\left(t\right)=\beta si-ϵe-\mu e,\\ i^{\prime }\left(t\right)=ϵe-\mu i-\gamma i,\\ v^{\prime }\left(t\right)=\sigma xs-\varphi v-\mu v.\end{array}$ (4)

2.2. 行为决策函数

接下来，我们可以对关于x的博弈过程进行描述。首先，疫苗接种的成本可以表示为$c_v$ ，而感染后的成本则定义为$c_i$ 。因此，我们可以将$h_v$ 定义为疫苗接种的收益，将$h_i$ 定义为不接种疫苗的收益。其中，根据传染病动力学的知识，我们可以得知未接种疫苗个体被感染的概率是$\beta i$ ，而疫苗失效后接种疫苗个体的感染风险为$\varphi \beta i$ 。请注意，$h_v$ 和$h_i$ 被定义为收益，因此它们的表达式为：

$\begin{array}{l}{h}_v=-\varphi \beta i{c}_i-c_v,\\ h_i=-\beta i{c}_i.\end{array}$

由此，我们可以得到一个信息函数$g_1\left(t\right)=h_v-h_i$ ，用于描述疫苗接种的成本。同时，对于疫情的传播，新增感染病例数将作为第二个信息函数$g_2\left(t\right)=\beta siN$ 。该函数展示了流行率变化对疫苗接种决策的影响。然后，我们使用m作为权重系数来连接这两个函数，其中$0<m<1$ 。若$m=1$ ，则决策完全取决于疫苗接种的成本；若$m=0$ ，则仅考虑新增感染的速度。此时，我们可以假设我们的净信息函数为

$g\left(t\right)=m{g}_1\left(t\right)+\left(1-m\right)g_2\left(t\right).$

最后，我们设定人群对信息的敏感度为$\rho$ ，且$\rho >0$ 。因此，行为响应方程为$M\left(t\right)=\rho g\left(t\right)$ 。我们可以得到最终的反应方程为

$\text{Response}\left(M\right)=\frac{1}{1+{\text{e}}^{-bM\left(t\right)}}-\frac{1}{2},$

其中$b>0$ 为信息响应强度。参考模仿行为的动力学，写出$x^{\prime }\left(t\right)$ 的表达式，并将其与我们微分方程组中的x的表达式进行对照：

$\begin{array}{l}{s}^{\prime }\left(t\right)=\mu -\beta si-\sigma xs-\mu s+\varphi v,\\ e^{\prime }\left(t\right)=\beta si-ϵe-\mu e,\\ i^{\prime }\left(t\right)=ϵe-\mu i-\gamma i,\\ v^{\prime }\left(t\right)=\sigma xs-\varphi v-\mu v,\\ x^{\prime }\left(t\right)=\nu x\left(1-x\right)\left(\frac{1}{1+{\text{e}}^{-bM\left(t\right)}}-\frac{1}{2}\right).\end{array}$ (5)

其中$\nu >0$ 为模仿强度的系数。

3. 平衡点的存在性和稳定性

3.1. 平衡点的计算和基本再生数

在本节中，我们将进行稳定性分析，以求解无病平衡点和地方病平衡点。为了找到这些点，我们将使方程组(5)中的每个方程等于零。平衡点被设定为$E_n=\left(s_n,e_n,i_n,v_n,x_n\right)$ 。

首先，存在一个无病平衡点。当没有接种意愿($x=0$ )时，可以求得平衡点$E_1$ 为

$E_1=\left(1,0,0,0,0\right).$

相反，如果每个人都自愿接种疫苗($x=1$ )，则可以得到另一个无病平衡点$E_2$ 为

$E_2=\left(\frac{\mu +\varphi }{\mu +\varphi +\sigma },0,0,\frac{\sigma }{\mu +\varphi +\sigma },1\right).$

感染个体引起的平均感染数由基本再生数$R_0$ 给出。基本再生数可用于估计病毒的最终规模。我们采用下一代矩阵方法[6]来确定基本再生数$R_0$ 。利用方程组(3)，无病平衡点$E_1$ 的雅可比矩阵如下：

$\begin{array}{l}F=\left[\begin{array}{ll}0\hfill & \beta S_0\hfill \\ 0\hfill & 0\hfill \end{array}\right],V=\left[\begin{array}{cc}\mu +ϵ& 0\\ -ϵ& \gamma +\mu \end{array}\right],\\ V^{-1}=\frac{1}{\left(\mu +ϵ\right)\left(\gamma +\mu \right)}\left[\begin{array}{cc}\gamma +\mu & 0\\ ϵ& \mu +ϵ\end{array}\right],\\ F{V}^{-1}=\frac{1}{\left(\mu +ϵ\right)\left(\gamma +\mu \right)}\left[\begin{array}{cc}\beta ϵS_0& \beta S_0\left(\mu +ϵ\right)\\ 0& 0\end{array}\right].\end{array}$

根据下一代矩阵方法，我们可以推导出：

$R_0=\rho \left(F{V}^{-1}\right)=\frac{ϵ\beta S_0}{\left(\mu +ϵ\right)\left(\gamma +\mu \right)}.$

将无病平衡点

$\begin{array}{l}{E}_1=\left(1,0,0,0,0\right)\\ E_2=\left(\frac{\mu +\varphi }{\mu +\varphi +\sigma },0,0,\frac{\sigma }{\mu +\varphi +\sigma },1\right)\end{array}$

代入相关方程，我们得到以下两个表达式：

$\begin{array}{l}{R}_{\nu }^0=R_0=\frac{ϵ\beta }{\left(\mu +ϵ\right)\left(\mu +\gamma \right)},\\ R_{\nu }^1=R_0\frac{\mu +\varphi }{\mu +\varphi +\sigma }.\end{array}$

接下来，若$i\ne 0,R_{\nu }^0>1$ ，且$x=0$ ，则在边界上存在一个地方性平衡点$E_3$ ：

$E_3=\left(\frac{\left(\mu +ϵ\right)\left(\mu +\gamma \right)}{ϵ\beta },\frac{\mu \left(R_{\nu }^0-1\right)\left(\mu +ϵ\right)}{R_{\nu }^0},\frac{\mu \left(\frac{\sigma }{\mu +\varphi +\sigma }-1\right)}{\beta },0,0\right).$

此外，若$i\ne 0,R_{\nu }^1>1$ ，且$x=1$ ，则在边界上存在另一个地方性平衡点$E_4$ ：

$E_4=\left(\frac{\left(\mu +ϵ\right)\left(\gamma +\mu \right)}{ϵ\beta },\left(\mu +\gamma \right)\frac{\mu \left(R_{\nu }^0-1\right)\left(\mu +\varphi \right)-\mu \sigma }{ϵ\beta \left(\mu +\varphi \right)},\frac{\mu \left(R_{\nu }^0-1\right)\left(\mu +\varphi \right)-\mu \sigma }{\beta \left(\mu +\varphi \right)},\frac{\sigma }{R_{\nu }^0\left(\mu +\varphi \right)},1\right).$

上述四个平衡点均为边界平衡点。接下来，我们将进一步计算系统的内部平衡点。

为此，我们将系统方程(5)的左侧设为零，得到以下方程组：

$\begin{array}{l}0=\mu -\beta s_5{i}_5-\sigma x_5{s}_5-\mu s_5+\varphi v_5,\\ 0=\beta s_5{i}_5-ϵe_5-\mu e_5,\\ 0=ϵe_5-\mu i_5-\gamma i_5,\\ 0=\sigma x_5{s}_5-\varphi v_5-\mu v_5.\end{array}$

通过解这个方程组，我们可以得出$s_5=1/R_{\nu }^0$ 。并且，当$x\left(t\right)=0$ 时，我们可以得出$M\left(t\right)=0$ 和$g\left(t\right)=0$ 。进一步推导，我们可以得到$i_5$ 的表达式：

$i_5=\frac{m{c}_v}{c_{i}m\beta \left(1-\varphi \right)+N\left(1-m\right)\beta s_5}=\frac{R_{\nu }^{0}m{c}_v}{R_{\nu }^0{c}_{i}m\beta \left(1-\varphi \right)+N\left(1-m\right)\beta }.$

然后，我们可以从方程(5)中推导出$e_5$ 、$v_5$ 和$x_5$ 的表达式：

$\begin{array}{l}{e}_5=\frac{\beta i_5}{R_{\nu }^0\left(ϵ+\mu \right)},\\ v_5=\frac{\mu \left(R_{\nu }^0-1\right)-\beta i_5}{\mu R_{\nu }^0},\\ x_5=\frac{\left(\varphi +\mu \right)\left[\mu \left(R_{\nu }^0-1\right)-\beta i_5\right]}{\sigma \mu }.\end{array}$

接下来，我们探讨平衡点$E_5=\left(s_5,e_5,i_5,v_5,x_5\right)$ 的存在条件，即该平衡点应满足以下条件：

$E_5\in \Gamma =\left\{\left(s_n,e_n,i_n,v_n,x_n\right)|s_n\ge 0,e_n\ge 0,i_n\ge 0,v_n\ge 0,0\le x_n\le 1\right\}\subset {ℝ}^5.$

当基本再生数$R_{\nu }^0>1$ 时，我们发现易感者比例$s_5$ 的取值范围为$0<s_5<1$ 。

由于参数m的取值范围在0到1之间，因此感染者比例$i_5$ 必然大于0。为了更深入地探讨恢复者比例$v_5$ 的取值范围，我们引入了一个新的常数$R_c$ 。为了简化计算，我们将$R_c$ 定义为：

$R_c=\frac{\frac{\mu ϵ}{\left(\mu +\gamma \right)\left(\mu +ϵ\right)}}{\frac{\mu ϵ}{\left(\mu +\gamma \right)\left(\mu +ϵ\right)}-i_5}.$

通过计算，我们发现当$R_c>R_{\nu }^0$ 时，恢复者比例$v_5$ 大于0。同时，根据$v_5$ 的表达式，我们可以推断出$v_5$ 必然小于1。

显然，接种疫苗者比例$x_5$ 大于0。此外，我们还推导出了$x_5<1$ 和$i_5<1$ 的条件，即当$R_c>R_{\nu }^0>1$ 时，这两个不等式同时成立。这表示在疾病传播过程中，感染者和接种疫苗者的比例都不会超过1。

3.2. 平衡点的稳定性

在此部分，我们着手探讨上述平衡点$E_n=\left(s_n,e_n,i_n,v_n,x_n\right)$ (其中$n=1,2,3,4,5$ )的稳定性(见表1)。

在此之前，我们首先在平衡点$E_n$ 处对原系统(5)重新表述为：

$\begin{array}{l}{s}^{\prime }\left(t\right)=\mu -\beta s_n{i}_n-\sigma x_n{s}_n-\mu s_n+\varphi v_n,\\ e^{\prime }\left(t\right)=\beta s_n{i}_n-ϵe_n-\mu e_n,\\ i^{\prime }\left(t\right)=ϵe_n-\mu i_n-\gamma i_n,\\ v^{\prime }\left(t\right)=\sigma x_n{s}_n-\varphi v_n-\mu v_n,\\ x^{\prime }\left(t\right)=M_1\left(n\right)s_n+M_2\left(n\right)i_n+M_3\left(n\right)x_n.\end{array}$ (6)

其中，我们设定

$\begin{array}{l}{M}_0\left(n\right)=\rho \left(m\left(\left(1-\varphi \right)\beta i_n{c}_i-c_v\right)+\left(1-m\right)\beta N{s}_n{i}_n\right),\\ M_1\left(n\right)=\nu x_n\left(1-x_n\right)\cdot \left(1-m\right)\rho \beta N\frac{b{\text{e}}^{-b{M}_0\left(n\right)}}{{\left(1+{\text{e}}^{-b{M}_0\left(n\right)}\right)}^2}{i}_n,\\ M_2\left(n\right)=\nu x_n\left(1-x_n\right)\cdot \frac{b{\text{e}}^{-b{M}_0\left(n\right)}}{{\left(1+{\text{e}}^{-b{M}_0\left(n\right)}\right)}^2}\left[m\left(1-\varphi \right)\beta c_i\rho +\left(1-m\right)\beta N{s}_n\rho \right],\\ M_3\left(n\right)=\nu \left(1-2x_n\right)\left(\frac{1}{1+{\text{e}}^{-b{M}_0\left(n\right)}}-\frac{1}{2}\right).\end{array}$

系统(5)的雅可比矩阵为

$\left|\begin{array}{ccccc}-\left(\beta i_n+\mu +\sigma x_n\right)& 0& -\beta s_n& \varphi & -\sigma s_n\\ \beta i_n& -\left(\mu +ϵ\right)& \beta s_n& 0& 0\\ 0& ϵ& -\left(\mu +\gamma \right)& 0& 0\\ \sigma x_n& 0& 0& -\left(\varphi +\mu \right)& \sigma s_n\\ M_1\left(n\right)& 0& M_2\left(n\right)& 0& M_3\left(n\right)\end{array}\right|=J.$ (7)

定理1：$R_{\nu }^0<1$ 时，则平衡点$E_1$ 是局部渐近稳定的。

证明：无病平衡点$E_1$ 处的雅可比矩阵为

$\left|\begin{array}{ccccc}-\left(\mu +\lambda \right)& 0& -\beta & \varphi & -\sigma \\ 0& -\left(\mu +ϵ+\lambda \right)& \beta & 0& 0\\ 0& ϵ& -\left(\mu +\gamma +\lambda \right)& 0& 0\\ 0& 0& 0& -\left(\varphi +\mu +\lambda \right)& \sigma \\ 0& 0& 0& 0& M_3\left(1\right)-\lambda \end{array}\right|=0.$ (8)

容易得出该方程关于${\lambda }_n$ (其中$n=1,2,3,4,5$ )的三个解。其中，${\lambda }_1=-\mu$ ，${\lambda }_2=M_3\left(1\right)<0$ ，${\lambda }_3=-\left(\mu +\varphi \right)$ 。原方程可简化为

$\left|\begin{array}{cc}-\left(\mu +ϵ+\lambda \right)& \beta \\ ϵ& -\left(\mu +\gamma +\lambda \right)\end{array}\right|={\lambda }^2+\left(2\mu +ϵ+\gamma \right)\lambda +\left(\mu +ϵ\right)\left(\mu +\gamma \right)-ϵ\beta =0.$

当$R_{\nu }^0<1$ 时，容易证明剩余的两个解${\lambda }_4,{\lambda }_5<0$ 。证毕。

定理2：$E_2$ 总是不稳定的。

证明：无病平衡点$E_2$ 处的雅可比矩阵为

$\left|\begin{array}{ccccc}-\left(\mu +\sigma \right)-\lambda & 0& -\beta \frac{\mu +\varphi }{\mu +\varphi +\sigma }& \varphi & -\sigma \frac{\mu +\varphi }{\mu +\varphi +\sigma }\\ 0& -\left(\mu +ϵ\right)-\lambda & \beta \frac{\mu +\varphi }{\mu +\varphi +\sigma }& 0& 0\\ 0& ϵ& -\left(\mu +\gamma \right)-\lambda & 0& 0\\ \sigma & 0& 0& -\left(\varphi +\mu \right)-\lambda & \sigma \frac{\mu +\varphi }{\mu +\varphi +\sigma }\\ 0& 0& 0& 0& M_3\left(2\right)-\lambda \end{array}\right|=0.$ (9)

由于

${\lambda }_1=M_3\left(2\right)=\nu \left(1-2x_1\right)\left(\frac{1}{1+{\text{e}}^{-b{M}_0\left(2\right)}}-\frac{1}{2}\right)$ ,

且$M_0\left(2\right)<0$ ，我们得知${\lambda }_1=M_3\left(2\right)>0$ ，因此$E_2$ 总是不稳定的。证毕。

定理3：若$R_c>R_{\nu }^1>1$ ，则$E_3$ 是局部渐近稳定的。

证明：地方病平衡点$E_3$ 的雅可比矩阵为

$\left|\begin{array}{ccccc}-\mu R_{\nu }^0-\lambda & 0& -\frac{\beta }{R_{\nu }^0}& \varphi & -\frac{\sigma }{R_{\nu }^0}\\ \mu \left(R_{\nu }^0-1\right)& -\left(\mu +ϵ\right)-\lambda & \frac{\beta }{R_{\nu }^0}& 0& 0\\ 0& ϵ& -\left(\mu +\gamma \right)-\lambda & 0& 0\\ 0& 0& 0& -\left(\varphi +\mu \right)-\lambda & \frac{\sigma }{R_{\nu }^0}\\ 0& 0& 0& 0& M_3\left(3\right)-\lambda \end{array}\right|=0.$ (10)

我们可以推导出，当基本再生数$R_{\nu }^0>1$ 且${\lambda }_2=-\left(\mu +\varphi \right)<0$ 时，有$M_0\left(3\right)<0$ 及

${\lambda }_1=M_3\left(3\right)=\nu \left(1-2x_3\right)\left(\frac{1}{1+{\text{e}}^{-b{M}_0\left(2\right)}}-\frac{1}{2}\right)<0$ ,

则原方程可以简化为：

$\left|\begin{array}{ccc}-\mu R_{\nu }^0-\lambda & 0& -\frac{\beta }{R_{\nu }^0}\\ \mu \left(R_{\nu }^0-1\right)& -\left(\mu +ϵ\right)-\lambda & \frac{\beta }{R_{\nu }^0}\\ 0& ϵ& -\left(\mu +\gamma \right)-\lambda \end{array}\right|={\lambda }^3+a_{31}{\lambda }^2+a_{32}\lambda +a_{33}=0.$

其中，系数$a_{31}=\mu R_{\nu }^0+2\mu +ϵ+\gamma >0$ ，$a_{32}=\mu R_{\nu }^0\left(\mu +ϵ\right)+\mu R_{\nu }^0\left(\mu +\gamma \right)>0$ ，$a_{33}=\mu \left(R_{\nu }^0-1\right)\left(\mu +ϵ\right)\left(\mu +\gamma \right)>0$ 。容易验证$a_{31}{a}_{32}-a_{33}>0$ 。根据赫尔维茨判据，平衡点$E_3$ 是渐近稳定的。证毕。

定理4：若$R_{\nu }^1>R_c$ ，则平衡点$E_4$ 是局部渐近稳定的。其中，临界值$R_c$ 定义为：

$R_c=\frac{\frac{\mu ϵ}{\left(\gamma +\mu \right)\left(ϵ+\mu \right)}}{\frac{\mu ϵ}{\left(\gamma +\mu \right)\left(ϵ+\mu \right)}-i_5}.$

证明：地方病平衡点$E_4$ 的表达式为

$\left(\frac{\left(\mu +ϵ\right)\left(\gamma +\mu \right)}{ϵ\beta },\left(\mu +\gamma \right)\frac{\mu \left(R_{\nu }^0-1\right)\left(\mu +\varphi \right)-\mu \sigma }{ϵ\beta \left(\mu +\varphi \right)},\frac{\mu \left(R_{\nu }^0-1\right)\left(\mu +\varphi \right)-\mu \sigma }{\beta \left(\mu +\varphi \right)},\frac{\sigma }{R_{\nu }^0\left(\mu +\varphi \right)},1\right).$

其特征方程为：

$\left|\begin{array}{ccccc}-\mu -\beta i_4-\sigma -\lambda & 0& -\frac{\beta }{R_{\nu }^0}& \varphi & -\frac{\sigma }{R_{\nu }^0}\\ \beta i_4& -\left(\mu +ϵ\right)-\lambda & \frac{\beta }{R_{\nu }^0}& 0& 0\\ 0& ϵ& -\left(\mu +\gamma \right)-\lambda & 0& 0\\ \sigma & 0& 0& -\left(\varphi +\mu \right)-\lambda & \frac{\sigma }{R_{\nu }^0}\\ 0& 0& 0& 0& M_3\left(4\right)-\lambda \end{array}\right|=0.$ (11)

我们可以得到$M_0\left(4\right)>0$ 及

${\lambda }_1=M_3\left(4\right)=\nu \left(1-2x_4\right)\left(\frac{1}{1+{\text{e}}^{-b{M}_0\left(4\right)}}-\frac{1}{2}\right)<0$ ,

因此，特征方程可以化简为：

${\lambda }^4+a_{41}{\lambda }^3+a_{42}{\lambda }^2+a_{43}\lambda +a_{44}=0.$

这里，我们令$A_1=\mu +ϵ,A_2=\mu +\gamma ,A_3=\mu +\varphi ,A_4=\mu +\sigma ,B=A_1+A_2+A_3+A_4$ ，因此

$\begin{array}{l}{a}_{41}=B+\beta i_4,\\ a_{42}=-\varphi \sigma +\left(A_1+A_2\right)A_3+\left(A_4+\beta i_4\right)\left(B-A_4\right),\\ a_{43}=\left(A_1+A_2\right)\left(A_3{A}_4-\varphi \sigma \right)+\beta i_4\left(A_1{A}_2+A_1{A}_3+A_2{A}_3\right),\\ a_{44}=\beta A_1{A}_2{A}_3{i}_4.\end{array}$

通过计算我们得到：

$a_{41}>0,a_{42}>0,a_{43}>0,a_{44}>0,$

$\left|\begin{array}{cc}{a}_{41}& a_{43}\\ 1& a_{42}\end{array}\right|=a_{41}{a}_{42}-a_{43}>0,$

$\left|\begin{array}{ccc}{a}_{41}& a_{43}& 1\\ 1& a_{42}& a_{44}\\ 0& a_{41}& a_{43}\end{array}\right|=a_{43}\left(a_{41}{a}_{42}-a_{43}\right)-a_{41}^2{a}_{44}>0.$

根据Routh-Hurwitz判据，平衡点$E_4$ 是渐近稳定的。证毕。

Table 1. The stability of boundary equilibrium points

表1. 边界平衡点的稳定性

平衡点	稳定性	稳定性条件
$E_1$	局部渐近稳定	$R_{\nu }^1<1$
$E_2$	不稳定
$E_3$	局部渐近稳定	$R_c>R_{\nu }^1>1$
$E_4$	局部渐近稳定	$R_{\nu }^1>R_c$

4. 数值模拟和敏感性分析

在本节中运用MATLAB对已建立的模型进行了一系列数值模拟。首先，我们展示了一些模型校准的结果，这些结果证实了系统(5)部分平衡点的全局稳定性。

4.1. 模型校准

首先我们设定以下参数值：$N=1000$ ，$\mu =0.001$ ，$\beta =0.121$ ，$\gamma =0.1$ ，$c_v=0.1$ ，$ϵ=0.1$ ，$\varphi =0.055$ ，$c_i=5000$ ，$m=0.4$ ，$b=0.1$ ，$\nu =0.002$ ，$\sigma =0.0008$ ，$\rho =0.01$ 。不难计算出对应的$R_{\nu }^0=1.0001$ ，$R_{\nu }^1=0.9878$ ，$R_c=1.0016$ 。依据第三章内容，此时$R_c>R_{\nu }^1>1$ ，对应系统内部平衡点$E_5$ 。在图1中，我们可以观察到在$t\to \infty$ 时所提系统其收敛于平衡点，确认此时其全局稳定性。

Figure 1. The dynamic curves of $S\left(t\right),E\left(t\right),I\left(t\right),R\left(t\right),V\left(t\right)$ and $x\left(t\right)$ at the equilibrium E5

图1. E5点处$S\left(t\right)，E\left(t\right)，I\left(t\right)，R\left(t\right)，V\left(t\right)$ 和$x\left(t\right)$ 的动力学曲线图

(a)

(b) (c)

Figure 2. The dynamic curves of $S\left(t\right),E\left(t\right),I\left(t\right),R\left(t\right),V\left(t\right)$ and $x\left(t\right)$ at the equilibrium E4. (a) when $x\left(0\right)\in \left(0,1\right)$ ; (b) when $x\left(0\right)=0$ ; (c) when $x\left(0\right)=1$

图2. E4点处$S\left(t\right)，E\left(t\right)，I\left(t\right)，R\left(t\right)，V\left(t\right)$ 和$x\left(t\right)$ 的动力学曲线图。(a) 当$x\left(0\right)\in \left(0,1\right)$ 时；(b) 当$x\left(0\right)=0$ 时；(c) 当$x\left(0\right)=1$ 时

接下来我们设定以下参数值：$N=1000$ ，$\mu =0.001$ ，$\beta =3.0$ ，$\gamma =0.14$ ，$c_v=500$ ，$ϵ=0.5$ ，$\varphi =0.055$ ，$c_i=5000$ ，$m=0.4$ ，$b=0.1$ ，$\nu =0.002$ ，$\sigma =0.8$ ，$\rho =0.01$ 。不难计算出对应的$R_{\nu }^0=21.2341$ ，$R_{\nu }^1=1.3891$ ，$R_c=-0.2557$ ，对应第三章内容，此时$R_{\nu }^1>R_c$ ，对应平衡点$E_4$ 。在图2中，我们分别设定了$x\left(0\right)=0$ 、$x\left(0\right)=1$ 和$x\left(0\right)\in \left(0,1\right)$ 三种情况，验证了所提系统其收敛于平衡点，确认其全局稳定性。

在疫苗接种成本($C_v$ )较高或感染成本极低的情况下，个体接种疫苗的意愿预计将出现明显下降。这一变化进而导致人群中缺乏足够免疫保护的比例上升，从而增加了感染的风险。图3(a)中的数值模拟结果有力地证实了这一现象，该结果显示，在平衡状态$E_4$ 下，感染个体的比例与疫苗接种成本($C_v$ )呈正相关关系。此外，当权重系数(m)倾向于强调疫苗接种成本时，从图3(b)的分析中可以明显看出，平衡状态$E_4$ 下的感染比例随m的增大而呈现出显著的上升趋势。这凸显了疫苗接种成本、个体决策过程以及它们对特定人群中疾病流行率的后续影响之间复杂的动态关系。

(a) (b)

Figure 3. The impact of vaccination cost $C_v$ and strategy weight coefficient m on the infection proportion $i_4$ under equilibrium state $E_4$

图3. 疫苗接种成本$C_v$ 和策略权重系数m对平衡状态$E_4$ 下感染比例$i_4$ 的影响

4.2. 敏感性分析

在传染病建模领域，偏秩相关系数(PRCC)敏感性分析作为一种关键的统计方法，对于定量评估不同参数对模型结果影响的大小具有重要意义。

Table 2. Specific descriptions of model variables and parameters

表2. 模型变量和参数的具体描述

符号	符号描述	数值	单位	来源
仓室变量
S(t)	易感人群
V(t)	已接种疫苗人群
E(t)	暴露者人群
I(t)	受感染人群
R(t)	恢复者人群
x(t)	行为动力学方程
参数设置
N	区域总人口数	1.09 * 106	人	[8]
$\mu$	人口自然衰减率	1/(78.87 * 365)	1/天	[9]
m	权重系数	0.1		[5]
$\beta$	接触率	0.45	1/(人*天)	[5]
$ϵ$	传递速率	0.5	1/天	[10]
$\sigma$	疫苗有效率	0.8	1/天	[11]
$\varphi$	疫苗防护衰退速率	0.0055	1/天	[12]
$\gamma$	感染恢复率	0.14	1/天	[13]
Ci	感染风险成本	5000
Cv	疫苗接种风险成本	500
$\rho$	疫苗注射比例	0.01	1/天
b	信息响应强度	0.1
$\nu$	模仿率	0.002	1/天

拉丁超立方抽样(LHS)是一种用于从多元参数分布中生成近似随机样本的方法[7]。

相关系数(通常用r表示)是变量之间线性关联程度的定量度量。其取值范围在−1到+1之间，正相关系数表示一个变量随另一个变量的增加而直接正比例增加，负相关系数则表示一个变量随另一个变量的增加而反比例减少。系数的绝对值接近0，表明两个变量之间的线性关系较弱。

我们以表2所提供的美国罗德岛州COVID-19疫情相关研究的参数[5] [8]-[13]为基础，将本模型中不同参数的影响可视化。图4展示了${ℛ}_{\nu }^0$ 、$R_{\nu }^1$ 、$R_c$ 与其对应的模型参数：$\beta$ 、$ϵ$ 、$\mu$ 、$\gamma$ 、$m$ 、$\sigma$ 、$\varphi$ 、$C_i$ 和$C_v$ 之间的相关性。在本文的传染病模型中，${ℛ}_{\nu }^0$ 、$R_{\nu }^1$ 和$R_c$ 共同决定了疫情的最终规模和趋势。通过观察下列柱状图，我们可以得出以下结论：

(a) (b)

(c) (d)

(e) (f)

Figure 4. Histogram of the distribution of the reproduction number ${ℛ}_{\nu }^0,R_{\nu }^1,R_c$ and sensitivity analysis of model parameters 图4. 再生数${ℛ}_{\nu }^0,R_{\nu }^1,R_c$ 的分布直方图及对模型参数的敏感性分析

首先，参数$\beta$ 在图4(b)中显示出最高的正偏秩相关系数(PRCC)，表明该因素显著加速了控制再生数$R_{\nu }^0$ 的增长；相反，$\gamma$ 的负PRCC值表明这些参数的增加导致了$R_{\nu }^0$ 的降低。

其次，在图4(d)中，参数$\varphi$ 展现出最显著的正偏秩相关系数(PRCC)，清晰表明其对加速控制再生数$R_{\nu }^1$ 的巨大贡献。相反，与$\sigma$ 相关的负PRCC值表明该参数的增加导致了$R_{\nu }^1$ 值的减小。

此外，在图4(e)中，参数$\varphi$ 和$m$ 显示出最高的正偏秩相关系数(PRCC)，表明它们在加速控制再生数$R_{\nu }^c$ 中起关键作用。相反，$\mu$ 和$C_i$ 的负PRCC值表明这些参数的增加导致了有效再生数$R_c$ 的降低，凸显了对疾病传播的抑制作用。

因此，显然这些参数在调节繁殖频率方面发挥着至关重要的作用，凸显了在控制疾病传播中的重要性。

4.3. 控制策略

我们可以观察到${ℛ}_{\nu }^0$ 、$R_{\nu }^1$ 、$R_c$ 与这些参数之间可能存在的正相关和负相关关系，并据此确定采取控制措施的必要性。为了有效实施针对病毒传播的必要防控措施，我们特别关注影响疫情的四个关键参数：传播率$\beta$ ；疫苗有效期$\varphi$ ；疫苗效力$\sigma$ ；感染风险成本$C_i$ 。为了阐明这些参数对暴露者和感染者数量的具体影响，我们利用表2中提供的参数值，进行了具体分析，如图5所示。

我们发现，在保持其他参数不变的情况下，增加感染患者的治疗成本($C_i$ )和疫苗效力($\sigma$ )会显著降低暴露者和感染者中疫情的规模和峰值。相反，降低传播率($\beta$ )和疫苗有效期($\varphi$ )则伴随着疫情规模和严重程度急剧上升。这一观察结果与我们的常识以及从疫情防控工作中获得的经验知识相吻合。例如，严格的疫情控制政策可以人为降低病毒传播率，使疫情的最终规模控制在可接受的范围内。快速、高效且长期有效的疫苗也能达到上述效果。同时，降低疫苗价格和提高病毒感染后的治疗成本，可以引导更多人及时接种疫苗，实现群体免疫。

综上所述，本研究为疫情防控提供了新的见解和策略，强调了提高疫苗接种率和有效性以及非药物干预措施在减少疾病传播方面的重要性。这些发现对指导实际疫情防控工作具有重要意义。

5. 结果与讨论

通过证明，我们确立了当${ℛ}_{\nu }^0\le 1$ 时，系统存在无病平衡点且局部稳定，为疾病防控工作提供了重要的理论支持。此外，我们还证明了当${ℛ}_{\nu }^0>1$ 时，某些边界地方性平衡点也存在且局部稳定，具体理论结果可见表1。这一发现对于理解不同强度下传染病传播的动态变化具有重要意义。然而，由于模型较为复杂，我们仅给出了内部地方性平衡点的存在，对于其稳定性还需进一步研究。

(a) (b)

(c) (d)

Figure 5. Graph of the relationship between some key parameters and the total infected population ($E+I$ ). (a) Epidemic transmission rate $\beta$ and proportion of total infected population ($E+I$ ); (b) Graph of vaccine failure time $\varphi$ versus proportion of total infected population ($E+I$ ); (c) Graph of vaccine efficacy $\sigma$ versus the proportion of total infected population ($E+I$ ); (d) Graph of infection risk cost $C_i$ versus proportion of total infected population ($E+I$ )

图5. 部分参数与感染总人口($E+I$ )关系图。(a) 疫情传播率$\beta$ 与感染总人口($E+I$ )比例图；(b) 疫苗失效时间$\varphi$ 与感染总人口($E+I$ )比例图；(c) 疫苗效力$\sigma$ 与感染总人口($E+I$ )比例图；(d) 感染风险成本$C_i$ 与感染总人口($E+I$ )比例图

我们对某些重要参数实施了控制策略，并观察到提高疫苗有效性、延长保护时间以及采取更严格的措施可以有效抑制COVID-19的传播，并显著缩短其持续时间。此外，降低居民接种成本并强调治疗疾病的高成本也可以带来这些效果。这些发现为这些措施的有效性提供了理论验证。

基金项目

吉林省自然科学基金面上项目20230101291JC。

NOTES

^*通讯作者。