带有疫苗接种的新型冠状病毒传播动力学分析

Transmission Dynamics of Corona Virus with Vaccination

摘要

本文通过数学建模研究疫苗接种和隔离措施对新型冠状病毒传播的影响机制， 通过对建立的模型进行定性分析， 利用下一代矩阵法确定了基本再生数， 并给出平衡点的存在性及其稳定性。从数值模拟可以得出， 提高疫苗有效率可以更快地控制新型冠状病毒的传播， 提高疫苗接种率可以推迟感染高峰时间和减少发病人数。除此之外， 应继续执行一些非药物干预措施， 如戴口罩、 隔离等， 可以大大减少最终感染规模和新型冠状病毒传播的时间， 最终有效控制。

Abstract

The influence mechanism of vaccination and quarantine measures on the spread of corona virus was studied through mathematical modeling. Through qualitative analysis of the established model， the basic reproductive number was determined by using the next generation matrix method， and the existence and stability of the equilibria were given. Numerical simulations show that improving vaccine effectiveness can control the spread of corona virus more quickly， and increasing vaccination rates can delay the peak time of infection and reduce the number of cases. In addition， we should continue to implement some non-drug interventions， such as wearing masks and isolation， which can greatly reduce the final scale of infection and the transmission time of corona virus， and ultimately effectively control the spread of corona virus.

Graphical abstract

关键词

疫苗接种 / 基本再生数 / 平衡点 / 稳定性

Key words

vaccination / basic reproductive number / equilibrium / stability

0 引　言

近年来， 传染病动力学研究进展迅速， 大量数学模型被用于分析各种各样的传染病问题， 以探索传染病的一般传染规律^［12］。由于新型冠状病毒对社会的危害比较严重， 许多学者在这方面做了不少工作， 提供了丰富的研究方法^［3］。

国内外学者用多种手段仿真预测了新型冠状病毒的传播， 但预测模型多是基于Kermack等^［4］提出的SIR仓室模型或SEIR仓室模型， 并通过增加仓室或改变参数等方式进行优化。针对新型冠状病毒的一些特点， 考虑了新型冠状病毒传播之间各种可能的相互作用， Herrera-Serrano等^［5］在一般流行病学假设下描述了新型冠状病毒的传播。

通过对新型冠状病毒传播的规律和特点的研究， 本文建立了新的SVAIQR模型， 并对模型的平衡点、 稳定性进行了分析， 揭示了新型冠状病毒传播的机理。从理论上验证了要实现新型冠状病毒全面快速的控制， 除了依靠疫苗接种来控制传播， 还应适当执行一些非药物干预措施。

1 数学模型

考虑到新型冠状病毒传播的特点， 将时间t时刻总人口N(t)分为易感染者S(t)、 接种疫苗者V(t)、 无症状感染者A(t)、 有症状感染者I(t)、 隔离者Q(t)、 治愈者R(t)等6类。 此模型不考虑潜伏期人口， 由于新型冠状病毒的潜伏期较短， 处于潜伏期的人口也具有传染性， 故将潜伏期并入无症状感染期， 并根据实际情况作出以下假设： 1） 只对易感染者接种疫苗； 2） 接种者也会感染新型冠状病毒并且感染后也具有传染性； 3） 无症状感染者也具有传染性。

基于以上假设的传播流程如图 1 所示。

流程图所对应的具体数学模型为

\left\{\begin{array}{l}\begin{array}{l}\begin{array}{l}\dot{S}=\Lambda -\tau S-\alpha (I+\beta A)S-\mu S,\\ \dot{V}=\tau S-\rho \alpha (I+\beta A)V-\mu V,\end{array}\\ \begin{array}{l}\dot{A}=\alpha {\lambda }_{\mathrm{1}}(I+\beta A)S+\rho \alpha {\lambda }_{\mathrm{2}}(I+\beta A)V-\\ \begin{array}{cc}& \end{array}\delta A-\eta A-\mu A,\end{array}\end{array}\\ \begin{array}{l}\dot{I}=\alpha \left(\mathrm{1}-{\lambda }_{\mathrm{1}}\right)(I+\beta A)S+\rho \alpha \left(\mathrm{1}-{\lambda }_{\mathrm{2}}\right)·\\ \begin{array}{cc}& \end{array}(I+\beta A)V+\delta A-\varphi I-\mu I,\end{array}\\ \dot{Q}=\eta A+\varphi I-\xi Q-\mu Q,\\ \dot{R}=\xi Q-\mu R,\end{array}\right.（1）

式中： \Lambda表示进入易感染者类的人口输入率； \tau为疫苗接种率； \alpha为有症状感染者的传播率； \beta为无症状感染者传播率的校正因子； \mu为所有人口的平均自然死亡率； ρ为疫苗失效率； {\lambda }_{\mathrm{1}}\left(\mathrm{0}\le {\lambda }_{\mathrm{1}}\le \mathrm{1}\right)为易感染者患病成为无症状感染者的比例； {\lambda }_{\mathrm{2}}\left(\mathrm{0}\le {\lambda }_{\mathrm{2}}\le \mathrm{1}\right)为接种疫苗者患病为有症状感染者的比例； \delta为无症状感染者转化为有症状感染者的比例； \eta为无症状感染者的隔离率； \varphi为有症状感染者的隔离率； \xi为隔离者的恢复率。

显然总人口N满足方程\dot{N}=\Lambda -\mu N,它的解为N\left(t\right)=\frac{\Lambda }{\mu }+\left(N\left(\mathrm{0}\right)-\frac{\Lambda }{\mu }\right){\mathrm{e}}^{-\mu t}。

故当t\in \left.\left[\mathrm{0},+\mathrm{\infty }\right.\right)时，

另外， 易感者S满足不等式

\dot{S}\le \Lambda -(\tau +\mu )S，

由比较定理可知， 当S(\mathrm{0})\le \frac{\Lambda }{\tau +\mu }时， 有t∈［0，+\mathrm{\infty }）时， S(t)\le \frac{\Lambda }{\tau +\mu }成立。

下面考虑系统（1）满足初始条件S(t)\ge \mathrm{0}， V(t)\ge \mathrm{0}， A(t)\ge \mathrm{0}， I(t)\ge \mathrm{0}， Q(t)\ge \mathrm{0}， R(t)\ge \mathrm{0}的解， 易知集合D=\left\{\left(S,V,A,I,Q,R\right)\begin{array}{c}\\ \end{array}\right.∈{\mathbf{R}}_{+}^{\mathrm{6}}， 0≤N≤\frac{\Lambda }{\mu }， \left.\mathrm{0}\le S\le \frac{\Lambda }{\tau +\mu }\right\}为系统 （1） 的正不变集。

由于系统（1）的前4个方程均与Q， R无关， 故系统（1）可简化为

\left\{\begin{array}{l}\dot{S}=\Lambda -\tau S-\alpha (I+\beta A)S-\mu S,\\ \dot{V}=\tau S-\rho \alpha (I+\beta A)V-\mu V,\\ \dot{A}=\alpha {\lambda }_{\mathrm{1}}(I+\beta A)S+\rho \alpha {\lambda }_{\mathrm{2}}(I+\beta A)V\\ \begin{array}{cc}& \end{array}-\delta A-\eta A-\mu A,\\ \dot{I}=\alpha \left(\mathrm{1}-{\lambda }_{\mathrm{1}}\right)(I+\beta A)S+\rho \alpha \left(\mathrm{1}-{\lambda }_{\mathrm{2}}\right)·\\ \begin{array}{cc}& \end{array}(I+\beta A)V+\delta A-\varphi I-\mu I\mathrm{。}\end{array}\right.（2）

2 基本再生数和平衡点存在性

2.1 基本再生数

令系统（2）右边全部为零， 可得方程组

\left\{\begin{array}{l}\Lambda -\tau S-\alpha (I+\beta A)S-\mu S=\mathrm{0},\\ \tau S-\rho \alpha (I+\beta A)V-\mu V=\mathrm{0},\\ \alpha {\lambda }_{\mathrm{1}}(I+\beta A)S+\rho \alpha {\lambda }_{\mathrm{2}}(I+\beta A)V-\\ \begin{array}{cc}& \end{array}\delta A-\eta A-\mu A=\mathrm{0},\\ \alpha \left(\mathrm{1}-{\lambda }_{\mathrm{1}}\right)(I+\beta A)S+\rho \alpha \left(\mathrm{1}-{\lambda }_{\mathrm{2}}\right)(I+\beta A)V+\\ \begin{array}{cc}& \end{array}\delta A-\varphi I-\mu I=\mathrm{0}\mathrm{。}\end{array}\right.（3）

令A=I=\mathrm{0}，可知系统（2）存在唯一的无病平衡点E_{\mathrm{0}}=\left(S_{\mathrm{0}},V_{\mathrm{0}},\mathrm{0,0}\right)，其中

S_{\mathrm{0}}=\frac{\mathrm{\Lambda }}{\tau +\mu }， V_{\mathrm{0}}=\frac{\tau \mathrm{\Lambda }}{\mu (\tau +\mu )}。

下面利用下一代矩阵法^［6］计算系统（2）的基本再生数， 新感染项与流出项分别为

\mathrm{ℱ}=\left[\begin{array}{c}\alpha {\lambda }_{\mathrm{1}}\left(I+\beta A\right)S+\rho \alpha {\lambda }_{\mathrm{2}}\left(I+\beta A\right)V\\ \alpha \left(\mathrm{1}-{\lambda }_{\mathrm{1}}\right)\left(I+\beta A\right)S+\rho \alpha \left(\mathrm{1}-{\lambda }_{\mathrm{2}}\right)\left(I+\beta A\right)V\end{array}\right],\mathit{𝒱}=\left(\begin{array}{c}\delta A+\eta A+\mu A\\ -\delta A+\varphi I+\mu I\end{array}\right)\mathrm{。}

可求得在E_{\mathrm{0}}处的感染矩阵\mathit{F}={\left.\frac{\partial \mathrm{ℱ}}{\partial \left(A,I\right)}\right|}_{{\mathrm{E}}_{\mathrm{0}}}和转移矩阵\mathit{V}={\left.\frac{\partial \mathit{𝒱}}{\partial \left(A,I\right)}\right|}_{{\mathrm{E}}_{\mathrm{0}}}。故系统（2）的基本再生数为

R_{\mathrm{0}}=\rho \left(\mathit{F}{\mathit{V}}^{-\mathrm{1}}\right)=\frac{\alpha {\lambda }_{\mathrm{1}}\delta \Lambda }{\left(\eta +\mu +\delta \right)\left(\varphi +\mu \right)\left(\tau +\mu \right)}+

\frac{\rho \alpha {\lambda }_{\mathrm{2}}\delta \Lambda \tau }{\mu \left(\eta +\mu +\delta \right)\left(\varphi +\mu \right)\left(\tau +\mu \right)}+

\frac{\alpha \beta {\lambda }_{\mathrm{1}}\Lambda }{\left(\eta +\mu +\delta \right)\left(\tau +\mu \right)}+\frac{\rho \alpha \beta {\lambda }_{\mathrm{2}}\Lambda \tau }{\mu \left(\eta +\mu +\delta \right)\left(\tau +\mu \right)}+\frac{\alpha \left(\mathrm{1}-{\lambda }_{\mathrm{1}}\right)\Lambda }{\left(\varphi +\mu \right)\left(\tau +\mu \right)}+\frac{\rho \alpha \left(\mathrm{1}-{\lambda }_{\mathrm{2}}\right)\Lambda \tau }{\mu \left(\varphi +\mu \right)\left(\tau +\mu \right)},

式中： \frac{\Lambda }{\mu +\tau }表示易感个体在无病人群中的比例； \frac{\alpha \beta \Lambda }{\mu +\tau }表示一个无症状个体每单位时间内感染的易感人群个体数； \frac{\mathrm{1}}{\eta +\mu +\delta }表示无症状感染者的持续时间； \frac{\alpha \beta {\lambda }_{\mathrm{1}}\Lambda }{\left(\eta +\mu +\delta \right)\left(\tau +\mu \right)}表示一个无症状感染者对易感染者的平均感染人数； \frac{\rho \alpha \beta {\lambda }_{\mathrm{2}}\Lambda \tau }{\mu \left(\eta +\mu +\delta \right)\left(\tau +\mu \right)}表示一个无症状感染者对接种疫苗者的平均感染人数； \frac{\alpha \Lambda }{\mu +\tau }表示一个有症状感染者每单位时间内感染的易感染者人群个体数； \frac{\mathrm{1}}{\varphi +\mu }表示有症状感染者的平均持续时间； \frac{\delta }{\eta +\mu +\delta }为无症状感染者发展为有症状感染者感染者的比例； \frac{\alpha \left(\mathrm{1}-{\lambda }_{\mathrm{1}}\right)\Lambda }{\left(\varphi +\mu \right)\left(\tau +\mu \right)}+\frac{\alpha {\lambda }_{\mathrm{1}}\delta \Lambda }{\left(\eta +\mu +\delta \right)\left(\varphi +\mu \right)\left(\tau +\mu \right)}表示一个有症状感染者对易感染者的平均感染人数； \frac{\rho \alpha {\lambda }_{\mathrm{2}}\delta \Lambda \tau }{\mu \left(\eta +\mu +\delta \right)\left(\varphi +\mu \right)\left(\tau +\mu \right)}+\frac{\rho \alpha \left(\mathrm{1}-{\lambda }_{\mathrm{2}}\right)\Lambda \tau }{\mu \left(\varphi +\mu \right)\left(\tau +\mu \right)}表示一个有症状感染者对接种疫苗者的平均感染人数。

2.2 平衡点的存在性

定理 1 系统（2）的平衡点如下：

1） 当R_{\mathrm{0}}\le \mathrm{1}时， 系统（2）只有无病平衡点E_{\mathrm{0}}；

2） 当R_{\mathrm{0}}>\mathrm{1}时， 系统（2）有一个无病平衡点E_{\mathrm{0}}和一个地方病平衡点E^{\mathrm{*}}。

证明 无病平衡点的存在性前面已知， 下面主要验证地方病平衡点存在的条件。

在方程组（3）中考虑满足A， I均大于零的解， 为求解简单， 令I+\beta A=P，则方程组（3）化为

\left\{\begin{array}{l}\Lambda -\tau S-\alpha PS-\mu S=\mathrm{0},\\ \tau S-\rho \alpha PV-\mu V=\mathrm{0},\\ \alpha {\lambda }_{\mathrm{1}}PS+\rho \alpha {\lambda }_{\mathrm{2}}PV-\delta A-\eta A-\mu A=\mathrm{0},\\ \alpha \left(\mathrm{1}-{\lambda }_{\mathrm{1}}\right)PS+\rho \alpha \left(\mathrm{1}-{\lambda }_{\mathrm{2}}\right)PV+\delta A-\\ \varphi I-\mu I=\mathrm{0}\mathrm{。}\end{array}\right.

可知

S=\frac{\mathrm{\Lambda }}{\tau +\alpha P+\mu },

V=\frac{\tau \mathrm{\Lambda }}{\left(\tau +\alpha P+\mu \right)\left(\rho \alpha P+\mu \right)},

A=\frac{\rho \alpha {\lambda }_{\mathrm{2}}P\tau \mathrm{\Lambda }}{\left(\tau +\alpha P+\mu \right)\left(\rho \alpha P+\mu \right)\left(\eta +\mu +\delta \right)}+\frac{\alpha {\lambda }_{\mathrm{1}}P\mathrm{\Lambda }}{\left(\tau +\alpha P+\mu \right)\left(\eta +\mu +\delta \right)},

I=\frac{\rho \alpha {\lambda }_{\mathrm{2}}P\tau \mathrm{\Lambda }\delta }{\left(\tau +\alpha P+\mu \right)\left(\rho \alpha P+\mu \right)\left(\eta +\mu +\delta \right)\left(\varphi +\mu \right)}+\frac{\alpha {\lambda }_{\mathrm{1}}P\mathrm{\Lambda }\delta }{\left(\tau +\alpha P+\mu \right)\left(\eta +\mu +\delta \right)\left(\varphi +\mu \right)}+

\frac{\rho \alpha \left(\mathrm{1}-{\lambda }_{\mathrm{2}}\right)P\tau \mathrm{\Lambda }}{\left(\tau +\alpha P+\mu \right)\left(\rho \alpha P+\mu \right)\left(\varphi +\mu \right)}+

\frac{\alpha \left(\mathrm{1}-{\lambda }_{\mathrm{1}}\right)P\mathrm{\Lambda }}{\left(\tau +\alpha P+\mu \right)\left(\varphi +\mu \right)}\mathrm{。}

从上面这些式子可以看出， 当P大于零时， S,V,A,I均大于零， 则方程组（3）的正解等价于求关于P的一元二次方程

a_{\mathrm{2}}{P}^{\mathrm{2}}+a_{\mathrm{1}}P+a_{\mathrm{0}}=\mathrm{0}（4）

的正解， 其中

a_{\mathrm{2}}=\left(\eta +\mu +\delta \right)\left(\varphi +\mu \right)\rho {\alpha }^{\mathrm{2}}>\mathrm{0},

\begin{array}{c}{a}_{\mathrm{1}}=\mu \alpha \left(\tau +\mu \right)\left(\eta +\mu +\delta \right)\left(\varphi +\mu \right)·\\ \left(\frac{\rho }{\mu }+\frac{\mathrm{1}}{\mu +\tau }-\frac{\rho \alpha \beta {\lambda }_{\mathrm{1}}\Lambda }{\mu \left(\eta +\mu +\delta \right)\left(\tau +\mu \right)}-\right.\end{array}

\frac{\rho \alpha {\lambda }_{\mathrm{1}}\delta \Lambda }{\mu \left(\eta +\mu +\delta \right)\left(\varphi +\mu \right)\left(\tau +\mu \right)}-\left.\frac{\rho \alpha \left(\mathrm{1}-{\lambda }_{\mathrm{1}}\right)\Lambda }{\mu \left(\varphi +\mu \right)\left(\tau +\mu \right)}\right)=

\left[\mu \alpha \left(\tau +\mu \right)\left(\eta +\mu +\delta \right)\left(\varphi +\mu \right)\right]\cdot \left[\frac{\rho }{\mu }\left(\mathrm{1}-R_{\mathrm{0}}\right)+\right.\frac{\mathrm{1}}{\mu +\tau }+\frac{{\rho }^{\mathrm{2}}\alpha {\lambda }_{\mathrm{2}}\delta \Lambda \tau }{{\mu }^{\mathrm{2}}\left(\eta +\mu +\delta \right)\left(\varphi +\mu \right)\left(\tau +\mu \right)}+\left.\frac{{\rho }^{\mathrm{2}}\alpha \beta {\lambda }_{\mathrm{2}}\Lambda \tau }{{\mu }^{\mathrm{2}}\left(\eta +\mu +\delta \right)\left(\tau +\mu \right)}+\frac{{\rho }^{\mathrm{2}}\alpha \left(\mathrm{1}-{\lambda }_{\mathrm{2}}\right)\Lambda \tau }{{\mu }^{\mathrm{2}}\left(\varphi +\mu \right)\left(\tau +\mu \right)}\right],

a_{\mathrm{0}}=\mu \left(\tau +\mu \right)\left(\eta +\mu +\delta \right)\left(\varphi +\mu \right)\left(\mathrm{1}-R_{\mathrm{0}}\right)\mathrm{。}

当，即R_{\mathrm{0}}>\mathrm{1}时， 方程（4）只有一个正实根P^{\mathrm{*}}，故R_{\mathrm{0}}>\mathrm{1}时系统（2）只有一个地方病平衡点E{}_{}{}^{\mathrm{*}}； 当a_{\mathrm{0}}\ge \mathrm{0}，即R_{\mathrm{0}}\le \mathrm{1}时， 易知此时a_{\mathrm{1}}>\mathrm{0}，方程（4）无正实根， 故R_{\mathrm{0}}\le \mathrm{1}时， 系统（2）无地方病平衡点。证毕。

3 平衡点的稳定性

3.1 无病平衡点的稳定性

定理 2 当时， 系统（2）的无病平衡点E_{\mathrm{0}}局部渐近稳定； 当R_{\mathrm{0}}>\mathrm{1}时， E_{\mathrm{0}}不稳定。

证明 系统（2）在无病平衡点E_{\mathrm{0}}处的雅可比矩阵\mathit{J}\left(E_{\mathrm{0}}\right)有两个负的特征值-\tau -\mu与-\mu，其余两个特征值满足方程

{\lambda }^{\mathrm{2}}+b_{\mathrm{1}}\lambda +b_{\mathrm{0}}=\mathrm{0},（5）

其中，

b_{\mathrm{0}}=\left(\eta +\delta +\mu \right)\left(\varphi +\mu \right)\left(\mathrm{1}-R_{\mathrm{0}}\right),

b_{\mathrm{1}}=\left(\eta +\delta +\mu \right)\left(\mathrm{1}-\frac{\alpha \beta {\lambda }_{\mathrm{1}}{S}_{\mathrm{0}}+\rho \alpha \beta {\lambda }_{\mathrm{2}}{V}_{\mathrm{0}}}{\eta +\mu +\delta }\right)+

\left(\varphi +\mu \right)\left(\mathrm{1}-\frac{\alpha \left(\mathrm{1}-{\lambda }_{\mathrm{1}}\right)S_{\mathrm{0}}+\rho \alpha \left(\mathrm{1}-{\lambda }_{\mathrm{2}}\right)V_{\mathrm{0}}}{\varphi +\mu }\right)>

\left(\eta +\delta +\varphi +\mathrm{2}\mu \right)\left(\mathrm{1}-R_{\mathrm{0}}\right)\mathrm{。}

当，即R_{\mathrm{0}}>\mathrm{1}时， 方程（5）有一个正实根， 故E_{\mathrm{0}}不稳定； 当b_{\mathrm{0}}>\mathrm{0}，即时， 此时b_{\mathrm{1}}>\mathrm{0}成立， 根据\mathrm{H}\mathrm{u}\mathrm{r}\mathrm{w}\mathrm{i}\mathrm{t}\mathrm{z}判据^［7］可知， 方程（5）的两个根具有负实部， 故当时， E_{\mathrm{0}}局部渐近稳定。证毕。

定理 3　当时， 系统（2）的无病平衡点E_{\mathrm{0}}在正不变集D内全局渐近稳定。

证明 根据文献［8］的方法， 令

f\left(\mathit{x},\mathit{y}\right):=\left(\mathit{F}-\mathit{V}\right)\mathit{x}-\mathrm{ℱ}\left(\mathit{x},\mathit{y}\right)+\mathit{𝒱}\left(\mathit{x},\mathit{y}\right)\mathrm{。}

系统（2）的染病仓室表示为

{\mathit{x}}^{\text{'}}=\left(\mathit{F}-\mathit{V}\right)\mathit{x}-f\left(\mathit{x},\mathit{y}\right),（6）

式中： \mathit{x}={\left(A,I\right)}^{\mathrm{T}}为染病仓室中的人群； \mathit{y}={\left(S,V\right)}^{\mathrm{T}}为非染病仓室中的人群。可求得非负矩阵{\mathit{V}}^{-\mathrm{1}}\mathit{F}对应的特征值R_{\mathrm{0}}的左特征向量{\mathit{\omega }}^{\mathrm{T}}=\left(\beta ,\mathrm{1}\right)。而

f\left(\mathit{x},\mathit{y}\right)=

\left[\begin{array}{c}\alpha {\lambda }_{\mathrm{1}}(I+\beta A)(S_{\mathrm{0}}-S)+\rho \alpha {\lambda }_{\mathrm{2}}(I+\beta A)(V_{\mathrm{0}}-V)\\ \alpha (\mathrm{1}-{\lambda }_{\mathrm{1}})(I+\beta A)(S_{\mathrm{0}}-S)+\rho \alpha (\mathrm{1}-{\lambda }_{\mathrm{2}})(I+\beta A)(V_{\mathrm{0}}-V)\end{array}\right],

易知在考虑的正不变集D内有S\le S_{\mathrm{0}}与V\le V_{\mathrm{0}}成立， 则f(\mathit{x},\mathit{y})\ge \mathrm{0}。

定义\mathrm{L}\mathrm{y}\mathrm{a}\mathrm{p}\mathrm{u}\mathrm{n}\mathrm{o}\mathrm{v}函数为

Q={\mathit{\omega }}^{\mathrm{T}}{\mathit{V}}^{-\mathrm{1}}\mathit{x}=\left[\frac{\beta }{\eta +\mu +\delta }+\right.\left.\frac{\delta }{\left(\eta +\mu +\delta \right)\left(\varphi +\mu \right)}\right]\cdot A+\frac{\mathrm{1}}{\varphi +\mu }I,

其关于系统（6）的全导数为

{\left|Q_{(\mathrm{6})}^{\text{'}}={\mathit{\omega }}^{\mathrm{T}}{\mathit{V}}^{-\mathrm{1}}{\mathit{x}}^{\mathit{\text{'}}}\right|}_{(\mathrm{6})}=(R_{\mathrm{0}}-\mathrm{1}){\mathit{\omega }}^{\mathrm{T}}\mathit{x}-{\mathit{\omega }}^{\mathrm{T}}{\mathit{V}}^{-\mathrm{1}}f(\mathit{x},\mathit{y})=(R_{\mathrm{0}}-\mathrm{1})(I+\beta A)-{\mathit{\omega }}^{\mathrm{T}}{\mathit{V}}^{-\mathrm{1}}f(\mathit{x},\mathit{y})\mathrm{。}

所以当时， 在正不变集D内Q^{\text{'}}\left|{}_{(\mathrm{6})}\le \mathrm{0}\right.,从而系统（2）的无病平衡点E_{\mathrm{0}}在正不变集D内全局渐近稳定。证毕。

3.2 地方病平衡点的稳定性

定理 4 当R_{\mathrm{0}}>\mathrm{1}时， 系统（2）唯一的地方病平衡点E^{\mathrm{*}}局部渐近稳定。

证明 系统（2）在地方病平衡点E^{\mathrm{*}}处的雅可比矩阵\mathit{J}\left(E^{\mathrm{*}}\right)的特征值满足方程

{\lambda }^{\mathrm{4}}+c_{\mathrm{1}}{\lambda }^{\mathrm{3}}+c_{\mathrm{2}}{\lambda }^{\mathrm{2}}+c_{\mathrm{3}}\lambda +c_{\mathrm{4}}=\mathrm{0},（7）

其中，

c_{\mathrm{1}}=\left(M_{\mathrm{1}}+M_{\mathrm{2}}\right)+X,

c_{\mathrm{2}}=\alpha P^{\mathrm{*}}\left(M_{\mathrm{5}}+\beta M_{\mathrm{7}}\right)+\rho \alpha P^{\mathrm{*}}\left(M_{\mathrm{6}}+\beta M_{\mathrm{8}}\right)+M_{\mathrm{1}}{M}_{\mathrm{2}}+M_{\mathrm{1}}{M}_{\mathrm{2}}X,

c_{\mathrm{3}}=\alpha P^{\mathrm{*}}{M}_{\mathrm{3}}\left(M_{\mathrm{5}}+\rho M_{\mathrm{6}}\right)+\alpha \beta P^{\mathrm{*}}{M}_{\mathrm{4}}\left(M_{\mathrm{7}}+\rho M_{\mathrm{8}}\right)+\alpha P^{\mathrm{*}}\delta \left(M_{\mathrm{7}}+\rho M_{\mathrm{8}}\right)+\alpha P^{\mathrm{*}}\left(M_{\mathrm{2}}{M}_{\mathrm{5}}+\rho M_{\mathrm{1}}{M}_{\mathrm{6}}\right)+\alpha \beta P^{\mathrm{*}}\left(M_{\mathrm{2}}{M}_{\mathrm{7}}+\rho M_{\mathrm{1}}{M}_{\mathrm{8}}\right)+M_{\mathrm{1}}{M}_{\mathrm{2}}X+\rho {\alpha }^{\mathrm{2}}\tau P^{\mathrm{*}}{S}^{\mathrm{*}},

c_{\mathrm{4}}=\alpha P^{\mathrm{*}}{M}_{\mathrm{3}}\left(M_{\mathrm{2}}{M}_{\mathrm{5}}+\rho M_{\mathrm{1}}{M}_{\mathrm{6}}+\rho \alpha \tau \left(\mathrm{1}-{\lambda }_{\mathrm{2}}\right)S^{\mathrm{*}}\right)+\alpha \beta P^{\mathrm{*}}{M}_{\mathrm{4}}\left(M_{\mathrm{2}}{M}_{\mathrm{7}}+\rho M_{\mathrm{1}}{M}_{\mathrm{8}}+\rho \alpha \tau {\lambda }_{\mathrm{2}}{S}^{\mathrm{*}}\right)+\alpha P^{\mathrm{*}}\delta \left(M_{\mathrm{2}}{M}_{\mathrm{7}}+\rho M_{\mathrm{1}}{M}_{\mathrm{8}}+\rho \alpha \tau {\lambda }_{\mathrm{2}}{S}^{\mathrm{*}}\right),

M_{\mathrm{1}}=\tau +\mu +\alpha P^{\mathrm{*}}>\mathrm{0},M_{\mathrm{2}}=\mu +\rho \alpha P^{\mathrm{*}}>\mathrm{0},

M_{\mathrm{3}}=\delta +\eta +\mu >\mathrm{0},M_{\mathrm{4}}=\varphi +\mu >\mathrm{0},

M_{\mathrm{5}}=\alpha \left(\mathrm{1}-{\lambda }_{\mathrm{1}}\right)S^{\mathrm{*}}>\mathrm{0},M_{\mathrm{6}}=\rho \alpha \left(\mathrm{1}-{\lambda }_{\mathrm{2}}\right)V^{\mathrm{*}}>\mathrm{0},M_{\mathrm{7}}=\alpha {\lambda }_{\mathrm{1}}{S}^{\mathrm{*}}>\mathrm{0},M_{\mathrm{8}}=\rho \alpha {\lambda }_{\mathrm{2}}{V}^{\mathrm{*}}>\mathrm{0},

X=\left(\varphi +\mu \right)\left(\frac{\beta A^{\mathrm{*}}}{P^{\mathrm{*}}}+\frac{\delta \left(\alpha {\lambda }_{\mathrm{1}}{S}^{\mathrm{*}}+\rho \alpha {\lambda }_{\mathrm{2}}{V}^{\mathrm{*}}\right)}{\left(\delta +\eta +\mu \right)\left(\varphi +\mu \right)}\right)+

\left(\delta +\eta +\mu \right)\frac{I^{\mathrm{*}}}{P^{\mathrm{*}}}>\mathrm{0},H_{\mathrm{1}}=c_{\mathrm{1}}>\mathrm{0},

H_{\mathrm{2}}=c_{\mathrm{1}}{c}_{\mathrm{2}}-c_{\mathrm{3}}>\rho \alpha P^{\mathrm{*}}{M}_{\mathrm{2}}\left(M_{\mathrm{6}}+\beta M_{\mathrm{8}}\right)+

\alpha P^{\mathrm{*}}\left(M_{\mathrm{5}}+\rho M_{\mathrm{6}}+\beta M_{\mathrm{7}}+\rho \beta M_{\mathrm{8}}\right)X+

\left(M_{\mathrm{1}}+M_{\mathrm{2}}\right)X^{\mathrm{2}}+{\left(M_{\mathrm{1}}+M_{\mathrm{2}}\right)}^{\mathrm{2}}X+

M_{\mathrm{1}}{M}_{\mathrm{2}}\left(M_{\mathrm{1}}+M_{\mathrm{2}}\right)>\mathrm{0},

H_{\mathrm{3}}=c_{\mathrm{3}}\left(c_{\mathrm{1}}{c}_{\mathrm{2}}-c_{\mathrm{3}}\right)-c_{\mathrm{1}}^{\mathrm{2}}{c}_{\mathrm{4}}>\left[\rho \alpha \right.P^{\mathrm{*}}{M}_{\mathrm{2}}\left(M_{\mathrm{6}}+\beta M_{\mathrm{8}}\right)+\left.\alpha P^{\mathrm{*}}\left(M_{\mathrm{5}}+\rho M_{\mathrm{6}}+\beta M_{\mathrm{7}}+\rho \beta M_{\mathrm{8}}\right)X\right]\cdot

\left[\rho {\alpha }^{\mathrm{2}}\tau P^{\mathrm{*}}\right.S^{\mathrm{*}}+\alpha P^{\mathrm{*}}{M}_{\mathrm{3}}\left(M_{\mathrm{5}}+\rho M_{\mathrm{6}}\right)+

\alpha P^{\mathrm{*}}\left(\beta M_{\mathrm{4}}+\delta \right)\left(M_{\mathrm{7}}+\rho M_{\mathrm{8}}\right)+M_{\mathrm{1}}{M}_{\mathrm{2}}X+

\alpha P^{\mathrm{*}}\left(M_{\mathrm{2}}{M}_{\mathrm{5}}+\rho M_{\mathrm{1}}{M}_{\mathrm{6}}\right)+

\left.\alpha \beta P^{\mathrm{*}}\left(M_{\mathrm{2}}{M}_{\mathrm{7}}+\rho M_{\mathrm{1}}{M}_{\mathrm{8}}\right)\right]>\mathrm{0},

H_{\mathrm{4}}=c_{\mathrm{4}}\left[c_{\mathrm{3}}\left(c_{\mathrm{1}}{c}_{\mathrm{2}}-c_{\mathrm{3}}\right)-c_{\mathrm{1}}^{\mathrm{2}}{c}_{\mathrm{4}}\right]=c_{\mathrm{4}}{H}_{\mathrm{3}}>\mathrm{0}\mathrm{。}

由Hurwitz判据可知， 地方病平衡点E^{\mathrm{*}}局部渐近稳定。证毕。

4 数值模拟

在系统（2）中选取第一组参数如表 1 第二列所示， 此时R_{\mathrm{0}}=0.272 1<1，经过模拟可看出系统（2）解的曲线最终趋向无病平衡点E_{\mathrm{0}}， 如图 2 所示。由理论分析可知， 此时E_{\mathrm{0}}全局渐近稳定， 数值模拟与理论分析一致。

在系统（2）中选取第二组参数如表 1 第4列所示， 此时R_{\mathrm{0}}=\mathrm{1.594}\mathrm{4}>\mathrm{1}，经过模拟可看出系统（2）的解的曲线最终趋向地方病平衡点E^{\mathrm{*}}，如图 3 所示。由此可看出， 当R_{\mathrm{0}}>\mathrm{1}时， E^{\mathrm{*}}稳定。

敏感性分析常用来确定模型对参数值的稳定性。用PRCC法进行敏感性分析， 如图 4 所示， 可以直观地看到各个参数与基本再生数R_{\mathrm{0}}的相关性。图 4 表明， 参数\Lambda ,\alpha ,\beta ,{\lambda }_{\mathrm{1}},{\lambda }_{\mathrm{2}},\rho与R_{\mathrm{0}}呈正相关， 而\tau ,\mu ,\eta ,\delta ,\varphi与参数R_{\mathrm{0}}呈负相关。敏感性分析说明， 在控制新型冠状病毒传播方面， 接种疫苗、 减少接触以及患病后采取一些得当的加强预防措施显得尤为重要。

根据疫苗疫苗统计数据， 取疫苗失效率\rho=0， 0.2， 0.3， 本文所建模型（2）在不同疫苗有效率下的无症状感染者和有症状感染者的解曲线， 如图 5 所示。

由图 5 中可知， 提高疫苗效率可以更快地控制COVID-19的传播， 即疫苗有效率越高（疫苗失效率越小）， 感染者达到顶峰的时间越短， 并且峰值也降得越低。

数值模拟所得不同疫苗接种率\tau=0， 0.5和0.897下感染者随时间变化的解曲线图， 如图 6 所示。由图 6 可知， 提高疫苗接种率可以让感染者的数量很快降低。

在现实中， 短时间内提高疫苗的有效率几乎不太可能， 而且仅依靠提高疫苗接种率也无法快速控制新型冠状病毒的传播。因此， 除了提高疫苗有效率和疫苗接种率之外， 还可以通过一些有效手段来减少传播率， 如佩戴口罩、 保持社交距离、 对感染者加大隔离力度， 对无症状感染者和有症状感染者提倡居家隔离或集中隔离。

取无症状感染者的隔离率\eta=0.1， 0.25和0.45，数值模拟所得感染者曲线图如图 7 所示。可以看出， \eta越大， 控制新型冠状病毒的速度越快， 当\eta=0.45时， 新型冠状病毒的传播在50 d左右基本完全控制。

类似地， 分别取有症状感染者隔离率\varphi=0.4， 0.6和0.8，如图 8 所示。可以看出， \varphi越大， 越有利于控制新型冠状病毒的传播。

5 结　论

本文为了研究新型冠状病毒的传播机理， 建立了一种具有疫苗接种和隔离措施的新型冠状病毒模型， 对模型的平衡点和稳定性进行分析， 并通过敏感性分析得到参数与R_{\mathrm{0}}的相关性， 同时通过数值模拟研究了采取一些措施对疫情的影响。仿真结果验证了疫苗接种在控制疫情传播过程中产生的影响是不可忽略的， 在提高疫苗有效率和接种率的同时要加入一些非药物干预措施， 如降低人群的接触率， 增加易感者的隔离率等措施会使控制效果更加显著。

参考文献

[1] 马知恩， 周义仓， 王稳地， 等 ．传染病动力学的数学建模与研究［M］．北京： 科学出版社， 2004.

[2] CANDIDO D， CLARO I， JESUS J， et al. Evolution and epidemic spread of SARS-CoV-2 in Brazil ［J］. Science， 2020， 369（6508）： 1255-1260.

[3] ZHOU B Q， ZHANG X H， JIANG D Q. Dynamics and density function analysis of a stochastic SVI epidemic model with half saturated incidence rate ［J］. Chaos Solitons and Fractals， 2020， 138： 109865.

[4] KERMACK W， MCKENDRICK A. A contribution to the mathematical theory of epidemics ［J］. Proceedings of Royal Society of London， Series A， 1927， 115（772）： 700-721.

[5] HERRERA-SERRANO J， MACIAS-DIAZ J， MEDINA-RAMIREZ I， et al. An efficient nonstandard computer method to solve a compartmental epidemiological model for COVID-19 with vaccination and population migration ［J］. Computer Methods and Programs in Biomedicine， 2022， 221： 106920.

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

[7] 马知恩， 周义仓. 常微分方程定性与稳定性方法［M］. 北京： 科学出版社， 2001.

[8] SHUAI Z， DRIESSCHE P. Global stability of infectious disease models using Lyapunov functions ［J］. SIAM Journal on Applied Mathematics， 2013， 73（4）： 1513-1532.

[9] SERHANI M， LABBARDI H. Mathematical modeling of COVID-19 spreading with asymptomatic infected and interacting peoples ［J］. Journal of Applied Mathematics and Computing， 2020， 66： 1-20.

基金资助

国家自然科学基金资助项目(11401541)

2016年中北大学优秀学术带头人出国进修学习资助项目