基于接种的新型冠状病毒模型及其分析：以上海为例

Novel Coronavirus Model Based on Vaccination and Its Analysis: A Case of Shanghai

摘要: 新冠病毒感染是一种由新型冠状病毒引起的急性传染性疾病，自2019年12月新冠病毒感染爆发以来，不仅危及人类的生命健康，还对全世界的经济发展造成了严重的影响，由于新冠病毒可在复制过程中不断适应宿主而产生突变，其具有更强的传播力和免疫逃逸能力，导致其对疫苗的敏感性降低，接种疫苗的人群仍旧有感染新冠病毒的风险，同时接种的疫苗在一定时间后会对人体的保护作用逐渐减弱，导致疫苗的有效性逐渐减弱，因此研究接种疫苗对新冠病毒的影响具有一定的意义，本文将接种人群单独考虑为一个仓室建立了数学模型。文中给出了模型的基本再生数以及正平衡点，通过构造Lyapunov函数证明了无病平衡点以及Rv > 1时唯一的正平衡点的稳定性。以上海2022年2月26日~5月31日的数据为例，结合当时上海所采取的政策，对模型进行数据拟合，并给出了上海在此期间的基本再生数为2.5709、5.1583、0.2982，最后通过数值模拟表明：提高每日的接种率、增加初始接种比例、能够有效降低确诊病例，从而降低新冠病毒的传染规模。但是即使初始接种比例或每日接种比例为100%，并不能完全阻断新冠病毒的传播。

Abstract: Novel coronavirus infection is an acute infectious disease caused by novel coronavirus. Since novel coronavirus infection broke out in December 2019, it not only endangers human life and health, but also has a serious impact on the economic development of the whole world. As novel coronavirus can constantly adapt to the host in the process of replication, it has a stronger transmission power and immune escape ability. Resulting in the reduction of its sensitivity to the vaccine, the vaccinat-ed population is still at risk of infection with novel coronavirus. At the same time, the protective ef-fect of the vaccinated vaccine on the human body will gradually weaken after a certain period of time, resulting in the gradual weakening of the effectiveness of the vaccine. Therefore, it is of certain significance to study the impact of vaccination on novel coronavirus. In this paper, the basic regen-eration number and the positive equilibrium point of the model are given, and the stability of the disease-free equilibrium point and the unique positive equilibrium point when Rv > 1 are proved by constructing the Lyapunov function. Taking the data of Shanghai from February 26 to May 31, 2022 as an example, and combining the policies adopted by Shanghai at that time, the model is fitted to the data. The basic regeneration numbers of Shanghai during this period are 2.5709, 5.1583 and 0.2982. Finally, numerical simulation shows that increasing the daily vaccination rate and increas-ing the initial vaccination rate can effectively reduce the confirmed cases, thereby reducing the in-fection scale of the new coronavirus. However, even if the initial vaccination rate or the daily vac-cination rate is 100%, it will not completely block the transmission of the novel coronavirus.

1. 引言

新型冠状病毒感染是一种由新型冠状病毒引起的急性传染性疾病，新冠病毒可在复制过程中不断适应宿主而产生突变。世界卫生组织(WHO)于2020年1月30日宣布将新型冠状病毒疫情列为“国际关注的突发公共卫生事件”，并且于2020年2月11日，WHO将由新型冠状病毒引发的疾病命名为2019冠状病毒病，英文名称为“Corona Virus Disease 19”(COVID-19) [1] 。同日，国际病毒命名委员会(ICTV)提议将新型冠状病毒命名为Severe Acute Respiratory Syndrome Coronavirus-2 (SARS-CoV-2) [2] 。截至 2023年3月29日，全球新型冠状病毒感染人数达761,402,282人，造成了61,887,000人的死亡 [3] 。

目前新型冠状病毒变异速度很快，以主要的流行毒株奥密克戎为例，其具有相当多的突变位点，其变异速度较快，且其下分支较多，导致其对疫苗的敏感性降低，接种疫苗人群仍旧有感染新型冠状病毒的风险 [4] [5] [6] ，随着时间推移，人体内的抗体含量持续下降，接种的疫苗在一定时间后会对人体的保护作用逐渐减弱，接种过疫苗的人体对新冠病毒的免疫能力也会下降 [7] ，因此我们考虑在接种疫苗的情况下对新型冠状病毒的影响。

近年来许多学者对新型冠状病毒的接种疫苗的分析从未停止，在 [8] 中Glenn Young等人建立了一个考虑接种的数学模型来分析新型冠状病毒的传播。研究表明由于新冠病毒疫苗接种率较低，在足够比例的人口接种疫苗之前，仍需实施非药物预防和其他控制措施来控制新型冠状病毒的传播。Cai等人 [9] 提出了具有年龄结构的新型冠状病毒随机动态模型，研究得出大多数新型冠状病毒死亡病例发生在未接种疫苗的人群中，尤其是老年群体中，得出提高老年人的接种率对减轻医疗负担有关键作用。Sun等人 [10] 考虑了医疗资源及接种疫苗等因素对新型冠状病毒的影响，建立了未接种人群与接种人群两群体分类讨论的动力学模型，分析了不同医疗资源及不同疫苗接种比例对上海新型冠状病毒的影响，得出了即使接种率达到100%，仍不能阻断新冠病毒传播的结论。Yang等人 [11] 提出了一个SVAIRS模型来分析新冠病毒是否可以通过单独接种疫苗完全消除，通过分析模型的基本再生数和平衡点的存在性，文章从理论上证明了仅仅依靠疫苗接种可能无法完全控制新型冠状病毒的传播。仍然有必要适当执行一些非药物干预措施。Tahir Khan等人 [12] 假设疫苗具有永久的免疫力，提出了SCAR数学模型，研究表明接种疫苗与新型冠状病毒的爆发有直接关系，快速的加大疫苗接种率，可以减少新型冠状病毒传播的持续时间。Mustafa等人 [13] 提出了SIRVI模型，模型中考虑接种疫苗后对于部分人群并不能提供完全有效的免疫力，导致这部分人依旧有感染新型冠状病毒的风险，通过分析表明通过疫苗接种可以有效控制新型冠状病毒的爆发。Julien Arino等人 [14] 提出了SLIARS动力学模型，进一步考虑了疫苗效能和疫苗免疫的减弱对新型冠状病毒的影响，通过动力学分析得出系统存在后向分支，并且对该模型的连续时间马尔可夫链(CTMC)进行了数值模拟，检测出后向分岔的存在性。上述学者都考虑了接种疫苗对新冠病毒传播的影响，对具有接种疫苗的传染病建模有较完善的理论研究，但是随着时间的推移，人体内中和抗体水平下降，疫苗的保护作用减弱，同时由于新型冠状病毒变异速度较快，具有较强的免疫逃逸能力，对于疫苗失效的情况的研究有所欠缺。

本文以上海新型冠状病毒的传播为例，将接种人群单独考虑为一个仓室建立了数学模型，研究了其内在机理。本文的研究结构如下：在第2节中给出了具有接种的新型冠状病毒模型及其参数的描述，第3节利用定性理论证明了正平衡点的存在性、全局稳定性。第4节运用MCMC法进行参数估计，然后根据得出的参数对上海2022年2月26~2022年5月31日的数据进行拟合，文中我们根据上海当时采取的政策将分为三个阶段处理，计算出了三个阶段的基本再生数，并且对基本再生数进行敏感性分析，然后分析不同的每日接种率和初始接种率对上海新型冠状病毒传播的影响。最后在第5节中讨论主要的结果。

2. 模型的建立

在这一部分中，我们构建 $S_1{S}_{2}EIAQR$ 模型来分析接种疫苗对新型冠状病毒的影响，把总人口分为八个仓室：未接种疫苗的易感者 $S_1\left(t\right)$ .接种疫苗的易感者 $S_2\left(t\right)$ 、暴露者 $E\left(t\right)$ 、有症状染病者 $I\left(t\right)$ 、无症状染病者 $A\left(t\right)$ 、确诊病例 $Q\left(t\right)$ 、恢复者 $R\left(t\right)$ ，以及因确诊新型冠状病毒而死亡者 $D\left(t\right)$ ，且有 $N\left(t\right)=S_1\left(t\right)+S_2\left(t\right)+E\left(t\right)+I\left(t\right)+A\left(t\right)+Q\left(t\right)+R\left(t\right)$ ，由于死亡仓室不影响模型的动力学分析，所以将死亡的确诊病例可以单独表示为 $D^{\prime }=dD$ ，对模型作如下假设：模型中采用双线性发生率；不考虑人口迁移对新型冠状病毒的影响，且新生儿不能接种疫苗；考虑到新型冠状病毒的传播方式有易感个体可以通过接触无症状和有症状的染病者而感染。可以得到流程图见图1：

Figure 1. The flow chart of novel coronavirus propagation with vaccination

图1. 具有接种的新型冠状病毒传播流程图

根据上述流程图，可以得到如下常微分方程组：

$\{\begin{array}{l}\frac{\text{d}{S}_1}{\text{d}t}=\Lambda -\beta S_1\left(I+\eta A\right)-\left(v+u\right)S_1+x{S}_2\\ \frac{\text{d}{S}_2}{\text{d}t}=v{S}_1-\theta \beta S_2\left(I+\eta A\right)-\left(u+x\right)S_2\\ \frac{\text{d}E}{\text{d}t}=\beta S_1\left(I+\eta A\right)+\theta \beta S_2\left(I+\eta A\right)-\left(\alpha +u\right)E\\ \frac{\text{d}I}{\text{d}t}=p\alpha E-\left({\lambda }_1+u\right)I\\ \frac{\text{d}A}{\text{d}t}=\left(1-p\right)\alpha E-\left({\lambda }_2+u+{\gamma }_2\right)A\\ \frac{\text{d}Q}{\text{d}t}={\lambda }_{1}I+{\lambda }_{2}A-\left({\gamma }_1+d+u\right)Q\\ \frac{\text{d}R}{\text{d}t}={\gamma }_{1}Q+{\gamma }_{2}A-uR\end{array}$ (1)

记为系统(1)，其中系统(1)中的参数及其所代表的意义如下(表1)：

Table 1. Description of parameters in system (1)

表1. 系统(1)中参数的描述

3. 动力学行为

系统(1)总存在一个无病平衡点 $E_0={\left(S_1^{\text{0}},S_2^{\text{0}},0,0,0,0,0\right)}^{\text{T}}$ 且系统(1)的可行域是正不变集

$\Omega =\left\{\left(S_1,S_2,E,I,A,Q,R\right)\in R_{+}^7:0\le S_1,S_2,E,I,A,Q,R,N\le \frac{\Lambda }{u}\right\}$ ，其中 $S_1^0=\frac{\Lambda \left(x+u\right)}{u^2+ux+uv};S_2^0=\frac{\Lambda v}{u^2+ux+uv}$ 。

3.1. 基本再生数

根据文章 [15] [16] 基本再生数以及下一代矩阵的定义，定义一组向量 $x={\left(E,I,A,Q\right)}^{\text{T}}$ ，仅包含系统(1)的致病项，由文章 [15] 的方法可得：

$F=\left(\begin{array}{cccc}0& \frac{\beta \Lambda \left(x+u+\theta v\right)}{vu+xu+u^2}& \frac{\beta \eta \Lambda \left(x+u+\theta v\right)}{vu+xu+u^2}& 0\\ 0& 0& 0& 0\\ 0& 0& 0& 0\\ 0& 0& 0& 0\end{array}\right);V=\left(\begin{array}{cccc}\alpha +u& 0& 0& 0\\ -p\alpha & {\lambda }_1+u& 0& 0\\ \left(p-1\right)\alpha & 0& {\lambda }_2+u+{\gamma }_2& 0\\ 0& -{\lambda }_1& -{\lambda }_2& d+{\gamma }_1+u\end{array}\right)$

其中 $F$ 是非负矩阵， $V$ 是非奇异矩阵

基本再生数 $R_v=\rho \left(F{V}^{-1}\right)$ ， $\rho$ 是矩阵的谱半径

$F{V}_1{}^{-1}=\left(\begin{array}{cccc}\frac{p\alpha \beta \left(S_1^{\text{0}}+\theta S_2^{\text{0}}\right)}{\left(\alpha +u\right)\left({\lambda }_1+u\right)}+\frac{\left(1-p\right)\left(\alpha +\kappa \right)\beta \eta \left(S_1^{\text{0}}+\theta S_2^{\text{0}}\right)}{\left(\alpha +u\right)\left({\lambda }_2+u+{\gamma }_2\right)}& \frac{\beta \left(S_1^{\text{0}}+\theta S_2^{\text{0}}\right)}{{\lambda }_1+u}& \frac{\beta \eta \left(S_1^{\text{0}}+\theta S_2^{\text{0}}\right)}{{\lambda }_2+{\gamma }_2+u}& 0\\ 0& 0& 0& 0\\ 0& 0& 0& 0\\ 0& 0& 0& 0\end{array}\right)$

计算可得基本再生数为：

$R_v=\phi \frac{S_1^0+\theta S_2^0}{\alpha +u}$

其中：

$\phi =\frac{p\alpha \beta }{{\lambda }_1+u}+\frac{\left(1-p\right)\alpha \beta \eta }{{\lambda }_2+u+{\gamma }_2}$

3.2. 稳定性分析

令系统(1)的右端为0，可以求得正平衡点 $E^{*}=\left(S_1{}^{*},S_2{}^{*},E^{*},I^{*},A^{*},Q^{*},R^{*}\right)$ 其中

$\{\begin{array}{l}{S}_1^{\text{*}}=\frac{\alpha +u}{\phi }-\theta S_2^{\text{*}}\\ S_2^{\text{*}}=\frac{v\left(\alpha +u\right)}{\phi \left(\theta \phi E^{*}+\theta v+x+u\right)}\\ I^{*}=\frac{p\alpha }{{\lambda }_1+u}{E}^{*}\\ A^{*}=\frac{\left(1-p\right)\alpha }{{\lambda }_2+{\gamma }_2+u}{E}^{*}\\ Q^{*}=\frac{p\alpha {\lambda }_1\left({\lambda }_2+{\gamma }_2+u\right)+\left(1-p\right)\alpha {\lambda }_2\left({\lambda }_1+u\right)}{\left({\gamma }_1+d+u\right)\left({\lambda }_1+u\right)\left({\lambda }_2+{\gamma }_2+u\right)}{E}^{*}\\ R^{*}=\frac{{\gamma }_1\left[p\alpha {\lambda }_1\left({\lambda }_2+{\gamma }_2+u\right)+\left(1-p\right)\alpha {\lambda }_2\left({\lambda }_1+u\right)\right]+{\gamma }_2\left(1-p\right)\alpha \left({\gamma }_1+d+u\right)\left({\lambda }_1+u\right)}{u\left({\gamma }_1+d+u\right)\left({\lambda }_1+u\right)\left({\lambda }_2+{\gamma }_2+u\right)}{E}^{*}\end{array}$ (2)

显然可知 $S_2^{\text{*}},E^{*},I^{*},A^{*},Q^{*},R^{*}>0$ ，现在证明 $S_1^{\text{*}}>0$ ，令 $f\left(E^{*}\right)=S_1^{\text{*}}=\frac{\alpha +u}{\phi }-\theta S_2^{\text{*}}$ ，对 $E^{*}$ 求导得：

$f^{\prime }\left(E^{*}\right)=\frac{{\theta }^{2}v\left(\alpha +u\right)}{{\left(\theta \phi E^{*}+\theta v+x+u\right)}^2}>0$

即 $S_1^{*}$ 是关于 $E^{*}$ 的单调递增函数，当 $E^{*}$ 等于0时 $S_1^{*}$ 最小，计算得 $S_1^{*}\ge \frac{\left(x+u\right)\left(\alpha +u\right)}{\phi \left(\theta v+x+u\right)}$ ，故证得 $S_1^{*}>0$ 。

接着我们将(2)代入系统 ，可得关于 $E^{*}$ 的表达式为 $f\left(E^{*}\right)=a{E}^{*}{}^2+b{E}^{*}+c=0$ ；且有

$\{\begin{array}{l}a=\theta {\phi }^2\left(\alpha +u\right)\\ b=-\theta {\phi }^2\Lambda +\left(\theta u+\theta v+x+u\right)\phi \left(\alpha +u\right)\\ c=u\left(\alpha +u\right)\left(x+v+u\right)\left(1-R_v\right)\end{array}$ (3)

其中：

$E_{1,2}^{*}=\frac{-b\pm \sqrt{b^2-4ac}}{2a},E_1^{*}+E_2^{*}=-\frac{b}{a},E_1^{*}*E_2^{*}=\frac{c}{a},a>0,sgn\left(c\right)=sgn\left(1-R_v\right)$

当 $R_v>1$ 时，有 $c<0,\frac{c}{a}<0$ ，此时 $f\left(E^{*}\right)$ 只有唯一一个正根，即系统(1)只有唯一一个正平衡点。

当 $R_v=1$ 时， $c=0$ ，根据

$R_v=\frac{\phi \Lambda \left(x+u+\theta v\right)}{\left(u^2+ux+uv\right)\left(\alpha +u\right)}=1\Rightarrow \phi \Lambda =\frac{\left(u^2+ux+uv\right)\left(\alpha +u\right)}{x+u+\theta v}$

将其代入b得

$\begin{array}{c}b=-\theta {\phi }^2\Lambda +\left(\theta u+\theta v+x+u\right)\phi \left(\alpha +u\right)\\ =\frac{\phi \left(\alpha +u\right)\left({\theta }^{2}uv+{\theta }^2{v}^2+2\theta vx+\theta uv+x^2+2xu\right)}{\theta v+x+u}>0\end{array}$

此时 $E_1^{*}=0,E_2^{*}=-\frac{b}{a}$ ， $f\left(E^{*}\right)$ 没有正根，即系统(1)只存在一个无病平衡点。

当 $R_v<1$ 时， $c>0$ ，根据

$R_v=\frac{\phi \Lambda \left(x+u+\theta v\right)}{\left(u^2+ux+uv\right)\left(\alpha +u\right)}<1\Rightarrow \phi \Lambda <\frac{\left(u^2+ux+uv\right)\left(\alpha +u\right)}{x+u+\theta v}$

将其代入b得

$\begin{array}{c}b=-\theta {\phi }^2\Lambda +\left(\theta u+\theta v+x+u\right)\phi \left(\alpha +u\right)\\ >-\theta \phi \frac{\left(u^2+ux+uv\right)\left(\alpha +u\right)}{x+u+\theta v}+\left(\theta u+\theta v+x+u\right)\phi \left(\alpha +u\right)\\ =\frac{\phi \left(\alpha +u\right)\left({\theta }^{2}uv+{\theta }^2{v}^2+2\theta vx+\theta uv+x^2+2xu\right)}{\theta v+x+u}>0\end{array}$

此时由 $E_1^{*}+E_2^{*}=-\frac{b}{a},E_1^{*}*E_2^{*}=\frac{c}{a}$ 得， $f\left(E^{*}\right)$ 没有正根，即系统(1)只存在一个无病平衡点。

综上所述，只有当 $R_v>1$ 时，系统(1)存在唯一一个正平衡点。

定理1：当 $R_v<1$ 时，系统(1)的无病平衡点在正不变集 $\Omega$ 内是局部渐进稳定的。

证明：

令 $s\left(M\right)=\max \left\{Re\lambda :\lambda 是矩阵M的特征值\right\}$ ，由参考文献的 [15] 定理2可得

$R_v>1\iff s\left(M\right)>0,R_v<1\iff s\left(M\right)<0$

$M=F-V=\left(\begin{array}{cccc}-\alpha -u& \frac{\beta \Lambda \left(x+u+\theta v\right)}{vu+xu+u^2}& \frac{\beta \eta \Lambda \left(x+u+\theta v\right)}{vu+xu+u^2}& 0\\ p\alpha & -{\lambda }_1-u& 0& 0\\ \left(1-p\right)\alpha & 0& -{\lambda }_2-{\gamma }_2-u& 0\\ 0& {\lambda }_1& {\lambda }_2& -d-{\gamma }_1-u\end{array}\right)$

在无病平衡点处得雅可比矩阵为 $J_E=\left(\begin{array}{cc}M& 0\\ -F& J_4\end{array}\right)$ ，其中

$J_4=\left(\begin{array}{ccc}-v-u& x& 0\\ v& -u-x& 0\\ 0& 0& -u\end{array}\right)$ ，

计算 $J_4$ 的特征值可得：

$s\left(J_4\right)=\max \left\{-u,-u,-u-x-v\right\}<0$ ，即 $R_v<1$ ，有 $s\left(M\right)<0$ 且 $s\left(J_4\right)<0$ ，可知系统(1)的无病平衡点在 $R_v<1$ 是局部渐进稳定的。

定理2：当 $R_v<1$ 时，系统(1)的无病平衡点在正不变集 $\Omega$ 内是全局渐进稳定的。

证明：对于无病平衡点 $E_0$ ，模型可以重新写作如下形式，记为系统(4)

$\{\begin{array}{l}\frac{\text{d}{S}_1}{\text{d}t}=S_1\left[\Lambda \left(\frac{1}{S_1}-\frac{1}{S_1^{\text{0}}}\right)+x\left(\frac{S_2}{S_1}-\frac{S_2^{\text{0}}}{S_1^{\text{0}}}\right)\beta \left(I+\eta A\right)\right]\\ \frac{\text{d}{S}_2}{\text{d}t}=S_2\left[v\left(\frac{S_1}{S_2}-\frac{S_1^{\text{0}}}{S_2^{\text{0}}}\right)-\theta \beta \left(I+\eta A\right)\right]\\ \frac{\text{d}E}{\text{d}t}=\beta \left(I+\eta A\right)\left[S_1^{\text{0}}+\theta S_2^{\text{0}}+S_1+\theta S_2-S_1^{\text{0}}-\theta S_2^{\text{0}}\right]-\left(\alpha +u\right)E\\ \frac{\text{d}I}{\text{d}t}=p\alpha E-\left({\lambda }_1+u\right)I\\ \frac{\text{d}A}{\text{d}t}=\left(1-p\right)\alpha E-\left({\lambda }_2+u+{\gamma }_2\right)A\\ \frac{\text{d}Q}{\text{d}t}={\lambda }_{1}I+{\lambda }_{2}A-\left({\gamma }_1+d+u\right)Y\end{array}$ (4)

假设 ${\left(S_1,S_2,E,I,A,Q,R\right)}^{\text{T}}$ 为系统(4)的任意正解，构造一个Lyapunov函数

$V=S_1-S_1^{\text{0}}-S_1^{\text{0}}\ln \frac{S_1}{S_1^{\text{0}}}+S_2-S_2^{\text{0}}-S_2^{\text{0}}\ln \frac{S_2}{S_2^{\text{0}}}+E+\frac{\beta \left(S_1^0+\theta S_2^0\right)}{{\lambda }_1+u}I+\frac{\beta \eta \left(S_1^0+\theta S_2^0\right)}{{\lambda }_2+{\gamma }_2+u}A$

V沿系统(4)的全导数为

$\begin{array}{c}\frac{\text{d}V}{\text{d}t}=\left(1-\frac{S_1^0}{S_1}\right)\frac{\text{d}{S}_1}{\text{d}t}+\left(1-\frac{S_2^0}{S_2}\right)\frac{\text{d}{S}_2}{\text{d}t}+\frac{\text{d}E}{\text{d}t}+\frac{\beta \left(S_1^0+\theta S_2^0\right)}{{\lambda }_1+u}\frac{\text{d}I}{\text{d}t}+\frac{\beta \eta \left(S_1^0+\theta S_2^0\right)}{{\lambda }_2+{\gamma }_2+u}\frac{\text{d}A}{\text{d}t}\\ =\Lambda \left(S_1-S_1^0\right)\left(\frac{1}{S_1}-\frac{1}{S_1^0}\right)+x\left(S_1-S_1^0\right)\left(\frac{S_2}{S_1}-\frac{S_2^0}{S_1^{\text{0}}}\right)+v\left(S_2-S_2^0\right)\left(\frac{S_1}{S_2}-\frac{S_1^0}{S_2^0}\right)\\ -\left(S_1-S_1^0\right){\beta }_1\left(I+\eta A\right)-\left(S_2-S_2^0\right)\theta \beta \left(I+\eta A\right)\\ +\beta \left(I+\eta A\right)\left(S_1^0+\theta S_2^0+S_1+\theta S_2-S_1^0-\theta S_2^0\right)-\left(\alpha +u\right)E\\ +\frac{\beta \left(S_1^0+\theta S_2^0\right)\left[p\alpha E-\left({\lambda }_1+u\right)I\right]}{{\lambda }_1+u}+\frac{\beta \eta \left(S_1^0+\theta S_2^0\right)\left[\left(1-p\right)\alpha E-\left({\lambda }_2+u+{\gamma }_2\right)A\right]}{{\lambda }_2+{\gamma }_2+u}\end{array}$

$\begin{array}{c}=\Lambda \left(2-\frac{S_1^0}{S_1}-\frac{S_1}{S_1^0}\right)+x\left(S_2-\frac{S_1{S}_2^0}{S_1^0}-\frac{S_1^0{S}_2}{S_1}+S_2^0\right)+v\left(S_1-\frac{S_2{S}_1^0}{S_2^0}-\frac{S_2^0{S}_1}{S_2}+S_1^0\right)\\ +\beta \left(I+\eta A\right)\left(S_1^0+\theta S_2^0\right)-\left(\alpha +u\right)E+\frac{\beta \left(S_1^0+\theta S_2^0\right)\left[p\alpha E-\left({\lambda }_1+u\right)I\right]}{{\lambda }_1+u}\\ +\frac{\beta \eta \left(S_1^0+\theta S_2^0\right)\left[\left(1-p\right)\alpha E-\left({\lambda }_2+u+{\gamma }_2\right)A\right]}{{\lambda }_2+{\gamma }_2+u}\end{array}$

令

$\{\begin{array}{l}\frac{\text{d}{S}_1}{\text{d}t}=S_1\left[\Lambda \left(\frac{1}{S_1}-\frac{1}{S_1^{\text{0}}}\right)+x\left(\frac{S_2}{S_1}-\frac{S_2^0}{S_1^0}\right)-\beta \left(I+\eta A\right)\right]=0\\ \frac{\text{d}{S}_2}{\text{d}t}=S_2\left[v\left(\frac{S_1}{S_2}-\frac{S_1^0}{S_2^0}\right)-\theta \beta \left(I+\eta A\right)\right]=0\end{array}$ (5)

得出

$\{\begin{array}{l}\Lambda =\left(u+v\right)S_1^0+x{S}_2^0\\ v{S}_1^0=\left(u+x\right)S_2^0\end{array}$ (6)

代入其中得

$\begin{array}{c}\frac{\text{d}V}{\text{d}t}=2\Lambda +x{S}_2^0+v{S}_1^0-\Lambda \frac{S_1^0}{S_1}-u{S}_1^{\text{0}}\frac{S_1}{S_1^0}-x{S}_2\frac{S_1^0}{S_1}-v{S}_1\frac{S_2^0}{S_2}-u{S}_2^{\text{0}}\frac{S_2}{S_2^0}-\left(\alpha +u\right)E\\ +\frac{\left[\beta \left(S_1^0+\theta S_2^0\right)\right]p\alpha E}{{\lambda }_1+u}+\frac{\left[\beta \eta \left(S_1^0+\theta S_2^0\right)\right]\left(1-p\right)\alpha E}{{\lambda }_2+{\gamma }_2+u}\\ =2\left[\left(u+v\right)S_1^0+x{S}_2^0\right]+x{S}_2^0+v{S}_1^0-\left[\left(u+v\right)S_1^0+x{S}_2^0\right]\frac{S_1^0}{S_1}\\ -u{S}_1^0\frac{S_1}{S_1^0}-x{S}_2\frac{S_1^0}{S_1}-v{S}_1\frac{S_2^0}{S_2}-u{S}_2^0\frac{S_2}{S_2^0}+\left(\alpha +u\right)E\left(R_v-1\right)\\ =u{S}_2^0\left(3-\frac{S_2}{S_2^0}-\frac{S_1^0}{S_1}-\frac{S_2^0{S}_1}{S_2{S}_1{}^0}\right)+x{S}_2^0\left(2-\frac{S_1^0{S}_2}{S_2^0{S}_1}-\frac{S_2^0{S}_1}{S_2{S}_1^0}\right)+u{S}_1^0\left(2-\frac{S_1}{S_1^0}-\frac{S_1^0}{S_1}\right)+\left(\alpha +u\right)E\left(R_v-1\right)\le 0\end{array}$

由上述可知当 $R_v<1$ 时， $\frac{\text{d}V}{\text{d}t}\le 0$ 当且仅当 $I=A=0$ ，时 $\frac{\text{d}V}{\text{d}t}=0$ ，即

$\left\{\left(S_1,S_2,E,I,A,Q,R\right)\in \Omega :\frac{\text{d}V}{\text{d}t}=0\right\}=\left\{E_0\right\}$

是系统(1)的最大不变集，由LaSalle不变集原理 [17] 可以得出无病平衡点 $E_0$ 是全局渐进稳定的。

定理3：当 $R_v>1$ 时，系统(1)的正平衡点 $E^{*}={\left(S_1^{\text{*}},S_1^{\text{*}},E^{*},I^{*},A^{*},Q^{*},R^{*}\right)}^{\text{T}}$ 在正不变集 $\Omega$ 内是全局渐进稳定的。

证明：

令 $s_1=\frac{S_1}{S_1^{*}},s_2=\frac{S_2}{S_2^{*}},e=\frac{E}{E^{*}},i=\frac{I}{I^{*}},a=\frac{A}{A^{*}},q=\frac{Q}{Q^{*}}$

模型可以重新写作如下形式：记为系统(7)

$\{\begin{array}{l}\frac{\text{d}{s}_1}{\text{d}t}=s_1\left(\frac{\Lambda }{S_1^{\text{*}}}\left(\frac{1}{s_1}-1\right)+\frac{x{S}_2^{\text{*}}}{S_1^{\text{*}}}\left(\frac{s_2}{s_1}-1\right)-\beta I^{*}\left(i-1\right)-\beta \eta A^{*}\left(a-1\right)\right)\\ \frac{\text{d}{s}_2}{\text{d}t}=s_2\left(\frac{v{S}_1^{\text{*}}}{S_2^{\text{*}}}\left(\frac{s_1}{s_2}-1\right)-\theta \beta I^{*}\left(i-1\right)-\theta \beta \eta A^{*}\left(a-1\right)\right)\\ \frac{\text{d}e}{\text{d}t}=\beta \frac{S_1^{\text{*}}{I}^{*}}{E^{*}}e\left(\frac{s_{1}i}{e}-1\right)+\beta \eta \frac{S_1^{\text{*}}{A}^{*}}{E^{*}}e\left(\frac{s_{1}a}{e}-1\right)+\theta \beta \frac{S_2^{\text{*}}{I}^{*}}{E^{*}}e\left(\frac{s_{2}i}{e}-1\right)+\theta \beta \eta \frac{S_2^{\text{*}}{a}^{*}}{E^{*}}e\left(\frac{s_{2}a}{e}-1\right)\\ \frac{\text{d}i}{\text{d}t}=p\alpha \frac{E^{*}}{I^{*}}i\left(\frac{e}{i}-1\right)\\ \frac{\text{d}a}{\text{d}t}=\left(1-p\right)\alpha \frac{E^{*}}{A^{*}}a\left(\frac{e}{a}-1\right)\\ \frac{\text{d}q}{\text{d}t}=\frac{{\lambda }_1{I}^{*}+{\lambda }_2{A}^{*}}{Q^{*}}q\left(\frac{{\lambda }_{1}I+{\lambda }_{2}A}{q\left({\lambda }_1{I}^{*}+{\lambda }_2{A}^{*}\right)}-1\right)\end{array}$ (7)

可以看出系统(7)有唯一一个正平衡点 $E^{*}\left(1,1,1,1,1\right)$ ，现证明 $E^{*}$ 的全局稳定性。

构造一个Lyapunov函数：

$\begin{array}{c}V\left(s_1,s_2,e,i,a,q\right)\\ =S_1^{*}\left(s_1-1-\ln s_1\right)+S_2^{*}\left(s_2-1-\ln s_2\right)\\ +E^{*}\left(e-1-\ln e\right)\\ +\frac{\beta \left(S_1^{*}+\theta S_2^{*}\right)I^{*}*I^{*}}{p\alpha E^{*}}\left(i-1-\ln i\right)\\ +\frac{\beta \eta \left(S_1^{*}+\theta S_2^{*}\right)A^{*}*A^{*}}{\left(1-p\right)\alpha E^{*}}\left(a-1-\ln a\right)\end{array}$

V沿系统(7)的全导数为

$\begin{array}{c}\frac{\text{d}V}{\text{d}t}=S_1^{\text{*}}\left(1-\frac{1}{s_1}\right)\frac{\text{d}{s}_1}{\text{d}t}+S_2^{\text{*}}\left(1-\frac{1}{s_2}\right)\frac{\text{d}{s}_2}{\text{d}t}+E^{*}\left(1-\frac{1}{e}\right)\frac{\text{d}e}{\text{d}t}+\frac{\beta \left(S_1^{\text{*}}+\theta S_2^{\text{*}}\right)I^{*}*I^{*}}{p\alpha E^{*}}\left(1-\frac{1}{i}\right)\frac{\text{d}i}{\text{d}t}\\ +\frac{\beta \eta \left(S_1^{\text{*}}+\theta S_2^{\text{*}}\right)A^{*}*A^{*}}{\left(1-p\right)\alpha E^{*}}\left(1-\frac{1}{a}\right)\frac{\text{d}a}{\text{d}t}\\ =\left(s_1-1\right)\left(\Lambda \left(\frac{1}{s_1}-1\right)+x{S}_2^{\text{*}}\left(\frac{s_2}{s_1}-1\right)-\beta S_1^{\text{*}}{I}^{*}\left(i-1\right)-\beta \eta S_1^{\text{*}}{A}^{*}\left(a-1\right)\right)\\ +\left(s_2-1\right)\left(v{S}_1^{\text{*}}\left(\frac{s_1}{s_2}-1\right)-\theta \beta S_2^{\text{*}}{I}^{*}\left(i-1\right)-\theta \beta \eta S_2^{\text{*}}{A}^{*}\left(a-1\right)\right)\\ +\left(e-1\right)\left(\beta S_1^{\text{*}}{I}^{*}\left(\frac{s_{1}i}{e}-1\right)+\beta \eta S_1^{\text{*}}{A}^{*}\left(\frac{s_{1}a}{e}-1\right)+\theta \beta S_2^{\text{*}}{I}^{*}\left(\frac{s_{2}i}{e}-1\right)+\theta \beta \eta S_2^{\text{*}}{a}^{*}\left(\frac{s_{2}a}{e}-1\right)\right)\end{array}$

$\begin{array}{c}+\beta \left(S_1^{\text{*}}+\theta S_2^{\text{*}}\right)I^{*}\left(e-i-\frac{e}{i}+1\right)+\beta \eta \left(S_1^{\text{*}}+\theta S_2^{\text{*}}\right)A^{*}\left(e-a-\frac{e}{a}+1\right)\\ =\Lambda \left(1-s_1-\frac{1}{s_1}+1\right)+x{S}_2^{\text{*}}\left(s_2-s_1-\frac{s_2}{s_1}+1\right)+v{S}_1^{\text{*}}\left(s_1-s_2-\frac{s_1}{s_2}+1\right)\\ -\left(s_1-1\right)\left[\beta S_1^{\text{*}}{I}^{*}i-\beta S_1^{\text{*}}{I}^{*}+\beta \eta S_1^{\text{*}}{A}^{*}a-\beta \eta S_1^{\text{*}}{A}^{*}\right]\\ -\left(s_2-1\right)\left[\theta \beta S_2^{\text{*}}{I}^{*}i-\theta \beta S_2^{\text{*}}{I}^{*}+\theta \beta \eta S_2^{\text{*}}{A}^{*}a-\theta \beta \eta S_2^{\text{*}}{A}^{*}\right]\\ +\beta S_1^{\text{*}}{I}^{*}\left(s_{1}i-e-\frac{s_{1}i}{e}+1\right)+\beta \eta S_1^{\text{*}}{A}^{*}\left(s_{1}a-e-\frac{s_{1}a}{e}+1\right)+\theta \beta S_2^{\text{*}}{I}^{*}\left(s_{2}i-e-\frac{s_{2}i}{e}+1\right)\end{array}$

$\begin{array}{c}+\theta \beta \eta S_2^{\text{*}}{a}^{*}\left(s_{2}a-e-\frac{s_{2}a}{e}+1\right)+\beta \left(S_1^{\text{*}}+\theta S_2^{\text{*}}\right)I^{*}\left(e-i-\frac{e}{i}+1\right)+\beta \eta \left(S_1^{\text{*}}+\theta S_2^{\text{*}}\right)A^{*}\left(e-a-\frac{e}{a}+1\right)\\ =2\Lambda +x{S}_2^{\text{*}}+v{S}_1^{\text{*}}-\frac{\Lambda }{s_1}-x{S}_2^{\text{*}}\frac{s_2}{s_1}-v{S}_1^{\text{*}}\frac{s_1}{s_2}-s_1\left(\Lambda +x{S}_2^{\text{*}}-v{S}_1^{\text{*}}-\beta S_1^{\text{*}}{I}^{*}-\beta \eta S_1^{\text{*}}{A}^{*}\right)\\ -s_2\left(v{S}_1^{\text{*}}-x{S}_2^{\text{*}}-\theta \beta S_2^{\text{*}}{I}^{*}-\theta \beta \eta S_2^{\text{*}}{A}^{*}\right)-\beta S_1^{\text{*}}{I}^{*}\frac{s_{1}i}{e}-\beta \eta S_1^{\text{*}}{A}^{*}\frac{s_{1}a}{e}-\theta \beta S_2^{\text{*}}{I}^{*}\frac{s_{2}i}{e}-\theta \beta \eta S_2^{\text{*}}{A}^{*}\frac{s_{2}a}{e}\\ +\beta \left(S_1^{\text{*}}+\theta S_2^{\text{*}}\right)I^{*}\left(1-\frac{e}{i}\right)+\beta \eta \left(S_1^{\text{*}}+\theta S_2^{\text{*}}\right)A^{*}\left(1-\frac{e}{a}\right)\end{array}$

令模型(1)等于0且代入 $E^{\text{*}}$ 得

$\{\begin{array}{l}\Lambda +x{S}_2^{\text{*}}=\beta S_1^{\text{*}}\left(I^{*}+\eta A^{*}\right)+\left(v+u\right)S_1^{\text{*}}\\ v{S}_1^{\text{*}}=\theta \beta S_2^{\text{*}}\left(I+\eta A\right)+\left(u+x\right)S_2^{\text{*}}\end{array}$ (8)

代入得：

$\begin{array}{l}\frac{\text{d}V}{\text{d}t}=2\Lambda +x{S}_2^{\text{*}}+v{S}_1^{\text{*}}-\frac{\Lambda }{s_1}-x{S}_2^{\text{*}}\frac{s_2}{s_1}-v{S}_1^{\text{*}}\frac{s_1}{s_2}-s_{1}u{S}_1^{\text{*}}-s_2{S}_2^{\text{*}}-\beta S_1^{\text{*}}{I}^{*}\frac{s_{1}i}{e}-\beta \eta S_1^{\text{*}}{A}^{*}\frac{s_{1}a}{e}\\ -\theta \beta S_2^{\text{*}}{I}^{*}\frac{s_{2}i}{e}-\theta \beta \eta S_2^{\text{*}}{A}^{*}\frac{s_{2}a}{e}+\beta \left(S_1^{\text{*}}+\theta S_2^{\text{*}}\right)I^{*}\left(1-\frac{e}{i}\right)+\beta \eta \left(S_1^{\text{*}}+\theta S_2^{\text{*}}\right)A^{*}\left(1-\frac{e}{a}\right)\\ =u{S}_1^{\text{*}}\left(2-\frac{1}{s_1}-s_1\right)+x{S}_2^{\text{*}}\left(2-\frac{s_2}{s_1}-\frac{s_1}{s_2}\right)+u{S}_2^{\text{*}}\left(3-\frac{s_1}{s_2}-\frac{1}{s_1}-s_2\right)\\ +\beta S_1^{\text{*}}{I}^{*}\left(3-\frac{1}{s_1}-\frac{s_{1}i}{e}-\frac{e}{i}\right)+\beta S_1^{\text{*}}\eta A^{*}\left(3-\frac{1}{s_1}-\frac{s_{1}a}{e}-\frac{e}{a}\right)\\ +\theta \beta S_2^{\text{*}}I\left(4-\frac{1}{s_1}-\frac{s_1}{s_2}-\frac{s_{2}i}{e}-\frac{e}{i}\right)+\theta \beta S_2^{\text{*}}\eta A^{*}\left(4-\frac{1}{s_1}-\frac{s_1}{s_2}-\frac{s_{2}a}{e}-\frac{e}{a}\right)\le 0\end{array}$

由上述可知 $\frac{\text{d}V}{\text{d}t}\le 0$ 当且仅当即 $S_{1,2}=S_{1,2}^{\text{*}},E=E^{*},I=I^{*},A=A^{*},Q=Q^{*},R=R^{*}$ 时 $\frac{\text{d}V}{\text{d}t}=0$ ，由LaSalle不

变集原理 [17] 可以得出唯一的正平衡点 是全局渐进稳定的。

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

4.1. 参数估计

因为上海在2月底开始爆发新型冠状病毒感染，于6月1日开始解封，故此处用上海2022年2月26日至2022年5月31日的报导的病例数据 [18] 来验证模型的精确性，模拟在此期间的累计确诊病例，用 $X\left(t\right)$ 代表累计确诊病例，所以有

$\frac{\text{d}X\left(t\right)}{\text{d}t}={\lambda }_{1}I+{\lambda }_{2}A$

此段时间内上海报告了60多万例当地的新型冠状病毒病例，根据上海政府在此期间采取的防控措施的变化，我们将此地新冠病毒的传播分为三个阶段，第一阶段是2022年2月26日至3月28日的阶段，假设它是没有任何干预的自由传播阶段。此时病例数目较少，也没有死亡病例出现，此时p = 0.04，d = 0 [18] ，随着新冠病毒的爆发，上海政府从2022年3月28日开始采取封控措施，在此期间上海基本属于封锁阶段，除了医务工作人员外，其他人几乎都待在家里，在2022年3月28日至4月16日，此时属于刚采取封控的前半个月，由于封控较晚，已经造成了社会性的传播，由于新型冠状病毒的潜伏期短，无症状染病者占绝大多数，所以即使在封控期间，会出现大量的染病者，政府在此期间采取了频繁的大规模核酸检测，导致聚集现象较多发生，因此在人群中有大量潜伏者以及无症状的感染者，易感者会与做核酸的染病者接触导致被感染，因在此期间，核酸检测频率增多，易感者被感染的机会增多，所以传染率比第一阶段增大，随着染病者的增多，医疗资源相对紧张，导致有一部人染病者不能及时被确诊。此时 ${\lambda }_1,{\lambda }_2$ 比第一阶段小，p = 0.064303379，d = 0 [18] ，第三阶段是2022年4月16日至5月31日，经过半个月的封控，每日报告的染病人数开始下降，绝大多数染病者都是在封控区被查出来，易感者与染病者接触的机会减少，以及核酸检测也有序排队，聚集现象较少发生，虽然每日报告的染病人数开始下降，但总体染病人数一直增加，所以假设染病者的确诊率还是下降的，此时p = 0.12680221，d = 0.02642364 [18] 。所以我们将分三个阶段来进行数值模拟，首先给出已知的参数，根据参考文献 [10] ，奥密克戎的平均潜伏期为1.52天，根据新型冠状病毒流行数据知，无症状人群相对有症状人群的传染性η = 0.75 [10] [19] ；由于新型冠状病毒可能逃避疫苗的免疫，中和抗体导致现有疫苗对抗奥秘克戎的有效性降低，假设θ = 0.8；假设接种疫苗的平均有效期为1年，x = 0.0028，同时考虑到疫情爆发期间，接种人数较少，所以假设持续接种率为v = 0.001，由 [20] 知，2022年上海出生人数为Λ = 288，自然死亡率为 $u=7.69{\text{e}}^{-06}$ ，根据 [10] 假设无症状染病者的平均染病时间 $1/{\gamma }_2=1/5.64$ ，确诊病例的平均染病时间 $1/{\gamma }_1=1/14$ ，上海的总人口为24871100，新型冠状病毒的平均疫苗接种率为72.6% [21] 。将初始值设置为 $S_1^0=6806881,S_2^0=18064219,$ $A^0=0,Q^0=0,R^0=0$ ；其中 $E^0,I^0$ 我们用最小二乘法进行拟合，对于未知的参数 $\beta 、{\lambda }_1、{\lambda }_2$ ，文中我们使用马尔可夫链蒙特卡洛法(MCMC) [21] 来进行估计。估计结果见表2，最后我们利用估计所得的参数计算三个阶段的基本再生数见表3。

4.2. 拟合结果

图2显示了上海2022.2.26~2022.5.31期间估计的累计确诊病例和实际确诊病例，并给出了模拟结果的95%置信区间。

Table 2. Parameter estimation with the method of MCMC (95% confidence interval)

表2. 使用MCMC方法进行参数估计(95%置信区间)

Table 3. Estimation of the value Rv in different stages (95% confidence interval)

表3. 不同阶段R_v的估计值(95%置信区间)

(a) (b)

Figure 2. (a) Fitting results of cases in Shanghai from 2022.2.26 to 2022.3.28; (b) Fitting results of cases in Shanghai from 2022.3.28 to 2022.5.31

图2. (a) 2022.2.26~2022.3.28上海累计确诊病例的拟合结果；(b) 2022.3.28~2022.5.31上海累计确诊病例的拟合结果

4.3. 敏感性分析

基本再生数 $R_v$ 是预测疾病是否会发生的阈值指数，因此确定模型(1)中对 $R_v$ 有重大影响的控制参数至关重要，以便进行相应的控制措施来消除疾病。表4给出基本再生数 $R_v$ 对每一个参数l的归一化敏感

性指数 ${\Gamma }_{R_v}^{\iota }=\frac{\partial R_v}{\partial \iota }\frac{l}{R_v}$ [22] 。

Table 4. Sensitivity index

表4. 敏感性指数

根据表4可以看出 $\beta ,\theta ,x,\eta$ 与 $R_v$ 呈正相关，即 $R_v$ 随着这些参数的增加而增加，意味着新型冠状病毒的流行更严重，而 ${\lambda }_1,{\lambda }_2,{\gamma }_2,v$ 与 $R_v$ 呈负相关。意味着增大这些参数可以减小基本再生数，从而控制新型冠状病毒的流行。因此可得结论：阻断传播、提高易感人群的接种率、提高疫苗的有效性、以及提高确诊率是控制新冠病毒的有效措施。为了更好的分析接种疫苗对新型冠状病毒传播的影响，我们考虑第一阶段，即在自由传播状态下参数对 $R_v$ 的影响。首先考虑单个参数对 $R_v$ 的影响，先考虑对 $R_v$ 影响较大的 $\eta$ 和 ${\gamma }_2$ 见图3、从图3(a)、图3(b)可以看出如果 ${\gamma }_2>0.6804,\eta <0.2667,R_v<1$ ，意味着新冠病毒将消失，从敏感性指数可知v和x对 $R_v$ 影响较小，从图3(c)、图3(d)可以看出即使每日接种率为100%，疫苗有效率为100%， $R_v$ 仍大于1，意味着新冠病毒并不会消失。即仅靠接种疫苗并不能完全阻断新型冠状病毒的传播，应该采取其他措施来控制疾病。接下来看图4所示的两个参数对 $R_v$ 的综合影响，从图4(a)中可以看到，在 ${\gamma }_2$ 和 $\eta$ 的联合作用下， $R_v$ 可以迅速减小。从图4(b)中还可以看到 ${\gamma }_2$ 和 $\eta$ 的阈值，以保持 $R_v<1$ 。

(a) (b)(c) (d)

Figure 3. The influence of parameters ${\gamma }_2$ , $\eta$ , v and x on the basic reproduction number R_v: (a) R_v in terms of different ${\gamma }_2$ ; (b) R_v in terms of different $\eta$ ; (c) R_v in terms of different v; (d) R_v in terms of different x

图3. 参数 ${\gamma }_2,\eta ,v,x$ 对基本再生数R_v的影响：(a) ${\gamma }_2$ 对R_v的影响；(b) $\eta$ 对R_v的影响；(c) v对R_v的影响；(d) x对R_v的影响

(a) (b)

Figure 4. R_v in terms of parameters $\eta$ and ${\gamma }_2$ .

图4. $\eta$ 和 ${\gamma }_2$ 对 $R_v$ 的影响

4.4. 接种率对新型冠状病毒确诊病例数的影响

为研究每日不同疫苗接种率、以及初始接种率对上海新型冠状病毒确诊病例的影响，下面给出通过变化v以及初始接种率后，所得到的相应的确诊病例数见图5，分别选取v = 0.001，v = 0.01，v = 0.05，v = 0.1，v = 0.3，v = 1，图5(a)显示随着v的增大，累计确诊病例数将减少，所以提高每日的接种率对新型冠状病毒的传播有抑制作用，但即使v = 1，上海新型冠状病毒累计确诊病例数仍有22,600。根据参考文献 [9] ，上海居民的新型冠状病毒初始接种率为72.6%，分别选取初始接种率为80%，85%，90%，95%，100%，图5(b)显示随着初始疫苗接种覆盖率的增加，上海累计确诊病例数将减少，因此我们应该尽可能扩大新型冠状病毒疫苗的推广范围以达到最大值，然而当疫苗的初始覆盖率达到100%，我们可以观察到上海累计确诊病例数将减少到27490，这表明疫苗接种可以适度缓解新型冠状病毒的传播，但是不是一项完整的控制措施，因此我们需要寻求其他策略来阻止新型冠状病毒在上海的传播，例如广泛的抗病毒药物等。

(a) (b)

Figure 5. (a) With different v values, the cumulative confirmed cases’ trajectory over time; (b) With different initial vaccination rate, the cumulative confirmed cases’ over time

图5. (a) 取不同的每日接种率v，累计确诊病例随时间变化的轨迹；(b) 取不同的初始接种率，累计确诊病例随时间变化的轨迹

5. 结论

基于2022年2月至5月上海新冠病毒传播的具体特点，我们建立了一个数学模型，模型中我们考虑具有持续的疫苗接种以及疫苗具有有效期，同时考虑疫苗并不能对人体起到100%的保护作用，即接种疫苗的人群仍有被感染新冠病毒的可能，我们建立此模型以研究接种疫苗对于上海此次新冠病毒传播的影响。首先我们根据下一代矩阵方法计算了模型的基本再生数，证明了平衡点的存在个数及其条件，通过构造Lyapunov函数证明了无病平衡点以及 $R_v>1$ 时唯一的正平衡点的稳定性。然后我们收集了上海卫生健康委员会于2022年2月26日至5月31日报告的确诊病例数据，并根据政府采取的一系列控制干预措施，将整个传播过程分为三个阶段。我们马尔可夫链蒙特卡洛法(MCMC)来估计模型的参数，并使用累积确诊病例数来验证模型的准确性，数值模拟结果如图2所示。且通过计算可得上海三个阶段的基本再生数分别为2.5709、5.1583、0.2982，意味着新型冠状病毒在上海最终会消失。根据求得的参数，对 $R_v$ 进行敏感性分析，并在图3中和图4我们分别画出了单个参数 ${\gamma }_2,\eta ,v,x$ 对 $R_v$ 的影响以及两个参数 ${\gamma }_2、\eta$ 对 $R_v$ 的影响。其次我们考虑了在不同的控制措施下对上海新冠病毒传播的最终规模的影响。根据图5我们可以得到提高每日的接种比例以及初始接种比例，能够有效降低确病例，从而降低疾病的传染规模。并且随着接种比例的增大，确诊人数减少，但是即使初始接种人数为100%，或者每日接种比例为100%，并不能完全阻断传染病的传播。所以我们应该采取其他方法来控制新型冠状病毒的传播，例如广泛的抗病毒药物等。然而本文也存在许多限制。首先我们认为所有人的传染率都是一样的。然而研究表明不同年龄结构对感染风险的影响异质性是不同的，这些留给我们在未来的工作中学习。

NOTES

^*通讯作者。