1. 引言

新型冠状病毒肺炎是一种由新型冠状病毒引起的传染性疾病。在很短的时间内迅速传播并成为一种全球性流行病。世界卫生组织(WHO)于2020年1月30日宣布新冠肺炎疫情为“国际关注的突发公共卫生事件”，并于2月11日在日内瓦的全球研究与创新论坛上宣布将新型冠状病毒肺炎命名为Corona Virus Disease-19 (COVID-19) [1]。同日，国际病毒分类委员会(ICTV)提议将新型冠状病毒命名为Severe Acute Respiratory Syndrome Coronavirus-2 (SARS-CoV-2) [2]。随后于2020年3月11日世界卫生组织(WHO)宣布COVID-19成为大流行。且截至2021年3月10日，全球新冠感染人数达118,188,896人，造成了2,619,825人的死亡 [3]。目前，最有效的控制策略是接种疫苗，通过减少易感者来有效控制疾病的传播。而接种疫苗会产生一定的经济支付和疫苗副作用风险，不接种疫苗则需承担染病风险以及感染的治疗费用，所以引入博弈免疫策略，可以使得人们通过衡量自身收益来决定是否接种疫苗，从而结合经济情况做出有利于自身的决策。

近年来许多学者也已经将博弈思维结合到传染病模型中对防控措施进行评估研究。博弈是关于包含相互依存关系情况中理性行为的研究。这种相互依存是指参与者理性地采取或选择自己的策略行为，其中博弈的任何一个参与者都会受到其他参与者的影响，并且他的行为也会影响到其他参与者，在这种相互影响的关系中，每一个参与者都尽可能地通过选择行为策略来提高自己的收益 [4]。而博弈在传染病传播动力中的应用，主要是预测个体行为，并研究其优化策略。在疾病爆发期间，个体往往会通过改变自身行为策略来降低感染风险，同时会产生一定的支付，但是如果不改变自身行为策略，就要承担染病和经济损失的风险。而个体则需要权衡这两个策略之间的收益，选择最优策略。例如Zhao等人 [5] 在传统SEIR模型的基础上加入博弈论决策过程的行为模仿过程来预测武汉新冠动态，模型估计基本再生数为2.5，得出个人行为变化会影响新冠疫情，且个人模仿率增加和被感染的可能性降低都有可能缓解新冠肺炎疫情。Kabira和Tanimotoc [6] 将传染病传播模型与进化博弈论(EGT)的行为动力学相结合，发现个人对感染风险和封闭在家里的经济成本的权衡会影响封锁在家的执行率。此外，当个人封锁成本相对较低时，个人更愿意封闭在家并且这样能够有效防控疾病流行。Wei等人 [7] 将新冠初期政府和公众行为策略博弈结合到传染病模型中，根据英国和中国累计确诊病例，确定新冠在不同策略下的传播规模，结果表明政府在疫情初期采取的应急策略能够有效控制疫情的传播。以上这些研究都是根据新冠传播机理从博弈角度考虑个人和政府的行为策略对新冠传播的影响，都没有从风险和个人支付角度考虑到疫苗免疫对疾病防控的影响。对于疫苗而言，一般情况下，疫苗的接种可以大大降低许多传染病的发病率和死亡率，并且随着越来越多的个体接种疫苗，在群体免疫的作用下，其余未接种疫苗的个体被感染的可能性就会越来越小。对于疫苗覆盖率足够高的人群来说，不需要给每个人接种疫苗就能够消除疾病 [8]。但是在自愿政策下，个人决定是否接种疫苗由许多因素决定，如疫苗的价格、对疫苗安全性的考虑、疫苗的供应量等。用博弈的方法可以量化个人根据这些因素的做出决策的过程，使得模型具有现实意义。目前，新冠疫苗已经开始全球大规模接种，各个国家都实行了有计划的疫苗接种行为。由于印度人口基数较大，医疗卫生水平有限，所以，本文以印度数据为例子研究考虑加入经济因素，将经济学博弈的思维融合到传染病模型中来研究接种行为对新冠传播的影响。

出于对监控力度和医疗设施的考虑，印度的确诊病人可能没有完全被隔离起来，也具有一定的传染性，所以本文建立SEIQR模型。其次，在流行性疾病模型的基础上，考虑引入疫苗博弈免疫，将个体根据风险信息和收益选择接种行为的过程形式化、公式化，通过计算再生数和平衡点，分析其动力学形态，得出自愿接种的政策下，疫苗博弈策略对新冠疫情的影响。

2. 模型

2.1. 模型建立

由于印度人口多且密集，医疗水平和隔离治疗不够完善的特点，对于已经确诊的病人，不一定能做到全部隔离治疗，所以确诊者也具有一定的传染能力，由此我们建立SEIQR模型，其中，S表示易感者，E表示潜伏者，I表示有症状的染病者，Q表示确诊者，R表示恢复者。易感者可以通过与潜伏者，有症状的染病者和确诊者接触而感染新冠，根据新冠传播特征，模型采用标准发生率，并且假设其出生、死亡率相等。如图1为新冠传播流程图：

根据图1的传播流程图，可得下面的微分方程组：

$\{\begin{array}{l}\frac{\text{d}S}{\text{d}t}=\mu N-\frac{\beta SI+{\beta }_{1}SE+{\beta }_{2}SQ}{N}-\mu S,\\ \frac{\text{d}E}{\text{d}t}=\frac{\beta SI+{\beta }_{1}SE+{\beta }_{2}SQ}{N}-\delta E-\mu E,\\ \frac{\text{d}I}{\text{d}t}=\delta E-mI-\mu I,\\ \frac{\text{d}Q}{\text{d}t}=mI-\gamma Q-\mu Q,\\ \frac{\text{d}R}{\text{d}t}=\gamma Q-\mu R.\end{array}$ (1)

其中， $\mu$ 表示出生死亡率， $\beta$ 表示有症状的染病者对易感者的传染率， ${\beta }_{\text{1}}$ 表示潜伏者对易感者的传染率， ${\beta }_{\text{2}}$ 表示确诊者对易感者的传染率， $\text{1}/\delta$ 表示潜伏者的潜伏期， $1/m$ 表示有症状染病者从出现症状到确诊的时间， $\text{1}/\gamma$ 表示确诊者的治愈期。

2.2. 引入博弈免疫策略

疫苗的机理是刺激人体激活免疫抗原，这里假设疫苗是100%有效的，且疫苗接种针对所有年龄阶段的人，包括所有现存易感者和新生儿。假设接种疫苗的经济成本和疫苗副作用风险为 $r_v$，则对于接种者，支付函数为

$f_v=-r_v.$

其他选择不接种疫苗的人则需要承担被感染的风险和感染后治疗的经济成本，设为 $r_i$，且不接种疫苗者的支付函数与不接种疫苗感染新冠的概率和新冠的流行率有关，所以不接种疫苗者的支付函数表示为

$f_n=-r_i\eta \left(E\left(t\right)+I\left(t\right)+Q\left(t\right)\right),$

其中 $\eta$ 表示不接种疫苗的个体感染疾病的可能性。

这时，转换策略到接种的支付增加为

$\Delta E=f_v-f_n.$

假设以 $\sigma$ 的速率对个体进行抽样，则非接种者 $\sigma x$ 的速率抽样接种者，则非接种者转换策略到接种者的概率为 $\rho \Delta E$，其中 $\rho$ 是比例常数。这里用 $x$ 表示接种者的比例，用 $\left(1-x\right)$ 表示非接种者的比例。则当 $\Delta E>0$ 时，非疫苗接种者可能转换为疫苗接种者，且有

$\frac{\text{d}x}{\text{d}t}=\left(1-x\right)\cdot \sigma x\cdot \rho \Delta E,$

当 $\Delta E<0$ 时，疫苗接种者可能会转换为非疫苗接种者，此时有

$\frac{\text{d}x}{\text{d}t}=-x\cdot \sigma \left(1-x\right)\cdot \rho \left(-\Delta E\right).$

如果取 $k=\sigma \rho$ 表示个人抽样其他人和转换策略的组合模拟率，则上述两个方程可以简化为

$\frac{\text{d}x}{\text{d}t}=kx\left(1-x\right)\left(f_v-f_n\right)=kx\left(1-x\right)\left[-r_v+r_i\eta \left(E+I+Q\right)\right].$

在上述博弈策略中非疫苗接种者的支付函数与感染新冠的概率和新冠的流行率有关，即与E(t)、I(t)和Q(t)有关，而E(t)、I(t)和Q(t)可以由传染病模型来确定。所以在模型(1)的基础上加入了博弈免疫，如图2为加入博弈免疫的新冠传播流程图。其中所有易感者的接种比例和新生儿的接种比例都为x。

引入两个参数： $\kappa =k{r}_v$， $\omega =\eta \frac{r_i}{r_v}$。微分方程组为：

$\{\begin{array}{l}\frac{\text{d}S}{\text{d}t}=\mu \left(1-x\right)N-\frac{\beta SI+{\beta }_{1}SE+{\beta }_{2}SQ}{N}-xS-\mu S,\\ \frac{\text{d}E}{\text{d}t}=\frac{\beta SI+{\beta }_{1}SE+{\beta }_{2}SQ}{N}-\delta E-\mu E,\\ \frac{\text{d}I}{\text{d}t}=\delta E-mI-\mu I,\\ \frac{\text{d}Q}{\text{d}t}=mI-\gamma Q-\mu Q,\\ \frac{\text{d}R}{\text{d}t}=\mu xN+xS+\gamma Q-\mu R,\\ \frac{\text{d}x}{\text{d}t}=\kappa x\left(1-x\right)\left[-1+\omega \left(E+I+Q\right)\right].\end{array}$ (2)

3. 2.2中模型(2)的分析

3.1. 没有人接种疫苗，即x = 0

当博弈免疫模型中的免疫比例x为0时，即没有人接种疫苗，此时模型(2)与模型(1)相同，用下一代矩阵法 [9] [10] 计算基本再生数：

$F=\left[\begin{array}{ccc}\beta & {\beta }_1& {\beta }_2\\ 0& 0& 0\\ 0& 0& 0\end{array}\right],V=\left[\begin{array}{ccc}\delta +\mu & 0& 0\\ -\delta & m+\mu & 0\\ 0& -m& \gamma +\mu \end{array}\right],$

$V^{-1}=\left[\begin{array}{ccc}\frac{\text{1}}{\delta +\mu }& 0& 0\\ \frac{\delta }{\left(\delta +\mu \right)\left(m+\mu \right)}& \frac{1}{m+\mu }& 0\\ \frac{\delta m}{\left(\delta +\mu \right)\left(m+\mu \right)\left(\gamma +\mu \right)}& \frac{m}{\left(m+\mu \right)\left(\gamma +\mu \right)}& \frac{1}{\gamma +\mu }\end{array}\right],$

$R_{\text{0}}=\rho \left(F{V}^{-1}\right)=\frac{\delta \left(\gamma +\mu \right)\beta +\left(m+\mu \right)\left(\gamma +\mu \right){\beta }_1+\delta m{\beta }_2}{\left(\gamma +\mu \right)\left(m+\mu \right)(\delta +\mu )}.$

又由于在模型(1)中， $S+E+I+Q+R=N$，故R可以消掉，并对系统作归一化变换，令

$s=\frac{S}{N},e=\frac{E}{N},i=\frac{I}{N},q=\frac{Q}{N},r=\frac{R}{N},$

微分方程为：

$\{\begin{array}{l}\frac{\text{d}s}{\text{d}t}=\mu -\left(\beta si+{\beta }_{1}se+{\beta }_{2}sq\right)-\mu s,\\ \frac{\text{d}e}{\text{d}t}=\beta si+{\beta }_{1}se+{\beta }_{2}sq-\delta e-\mu e,\\ \frac{\text{d}i}{\text{d}t}=\delta e-mi-\mu i,\\ \frac{\text{d}q}{\text{d}t}=mi-\gamma q-\mu q.\end{array}$ (3)

系统(3)的可行域为

$D_1=\left\{\left(s,e,i,q\right)\in R^4|0\le s,e,i,q\le 1,0\le s+e+i+q\le 1\right\}$

且 $D_1$ 为正向不变集。

通过计算，系统有一个无病平衡点 ${\zeta }_1\left(s^0,e^0,i^0,q^0\right)=\left(1,0,0,0\right)$ 和一个地方病平衡点

${\zeta }_2\left(s^{*},e^{*},i^{*},q^{*}\right)=\left(\frac{1}{R_0},\frac{\mu \left(R_0-1\right)}{\left(\delta +\mu \right)R_0},\frac{\delta \mu \left(R_0-1\right)}{\left(\delta +\mu \right)\left(m+\mu \right)R_0},\frac{\delta m\mu \left(R_0-1\right)}{\left(\delta +\mu \right)\left(m+\mu \right)\left(\gamma +\mu \right)R_0}\right)$.

定理1 1) 当 $R_0<1$ 时，无病平衡点 ${\zeta }_1$ 在正向不变集 $D_1$ 内是全局渐近稳定的；

2) 当 $R_0>1$ 时，地方病平衡点 ${\zeta }_2$ 在正向不变集 $D_1$ 内是全局渐近稳定的。

证明：1) 令

$M=F-V=\left(\begin{array}{ccc}{\beta }_1-\delta -\mu & \beta & {\beta }_2\\ \delta & -\mu & 0\\ 0& 0& -\gamma -\mu \end{array}\right),$

令S(M) = max{Re $\lambda$ : $\lambda$ 是矩阵M的特征值}，由文献 [10] 的定理2可得

$R_0>1\iff S\left(M\right)>0,R_0<1\iff S\left(M\right)<0.$

${\zeta }_1$ 的雅可比矩阵为

$J\left({\zeta }_1\right)=\left(\begin{array}{cc}M& 0\\ J_1& J_2\end{array}\right),$

其中

$J_1=-F,J_{\text{2}}\left({\zeta }_1\right)=\left(\begin{array}{cc}-\mu & 0\\ 0& -\mu \end{array}\right).$

$J_{\text{2}}$ 的特征值为 ${\lambda }_{\text{1}}={\lambda }_2=-\mu$。当 $R_0<1$ 时， $S\left(M\right)<0$， $S\left(J_{\text{2}}\left({\zeta }_1\right)\right)<0$。所以无病平衡点 ${\zeta }_1$ 在 $R_0<1$ 时局部渐近稳定。

下证无病平衡点 ${\zeta }_1$ 的全局渐近稳定性。构造如下所示Lyapunov函数

$L_1=s-s^0-s^0\ln \frac{s^0}{s}+e+\left[\frac{\beta s^0}{m+\mu }+\frac{m{\beta }_2{s}^0}{\left(m+\mu \right)\left(\gamma +\mu \right)}\right]i+\frac{{\beta }_2{s}^0}{\gamma +\mu }q,$

$\begin{array}{c}\frac{\text{d}{L}_1}{\text{d}t}=\left(1-\frac{s^0}{s}\right)s^{\prime }+e^{\prime }+\left[\frac{\beta s^0}{m+\mu }+\frac{m{\beta }_2{s}^0}{\left(m+\mu \right)\left(\gamma +\mu \right)}\right]i^{\prime }+\frac{{\beta }_2{s}^0}{\gamma +\mu }{q}^{\prime }\\ =\left(1-\frac{s^0}{s}\right)\left[\mu -\left(\beta si+{\beta }_{1}se+{\beta }_{2}sq\right)-\mu s\right]+\left(\beta si+{\beta }_{1}se+{\beta }_{2}sq-\delta e-\mu e\right)\\ +\left[\frac{\beta s^0}{m+\mu }+\frac{m{\beta }_2{s}^0}{\left(m+\mu \right)\left(\gamma +\mu \right)}\right]\left(\delta e-mi-\mu i\right)+\frac{{\beta }_2{s}^0}{\gamma +\mu }\left(mi-\gamma q-\mu q\right)\\ =\mu s^0\left(1+\frac{1}{s}-\frac{1}{s^0}-\frac{s}{s^0}\right)+\left[{\beta }_1-\left(\delta +\mu \right)+\frac{\delta \beta }{m+\mu }+\frac{\delta m{\beta }_2}{\left(m+\mu \right)\left(\gamma +\mu \right)}\right]e\\ =\mu s^0\left(1+\frac{1}{s}-\frac{1}{s^0}-\frac{s}{s^0}\right)+\left(\delta +\mu \right)\left(R_0-1\right)e\end{array}$

由不等式 $1-a\le -\ln a$， $a>0$，有 $1+\frac{1}{s}-\frac{1}{s^0}-\frac{s}{s^0}\le \frac{1}{s}-\ln \frac{1}{s}+\ln \frac{1}{s^0}-\frac{1}{s^0}$。

因此，在基本再生数 $R_0<1$ 时， $\frac{\text{d}{L}_1}{\text{d}t}<0$ 成立。所以由LaSalle’s不变集原理 [11] 可知 $R_0<1$ 时，无病平衡点 ${\zeta }_1$ 在 $D_1$ 内是全局渐近稳定的。

2) 对于地方病平衡点 ${\zeta }_{\text{2}}\left(s^{*},e^{*},i^{*},q^{*},r^{*}\right)$，构造如下所示Lyapunov函数

$L_{\text{2}}=s-s^{\text{*}}-s^{\text{*}}\ln \frac{s^{\text{*}}}{s}+e-e^{\text{*}}-e^{\text{*}}\ln \frac{e^{\text{*}}}{e}+\frac{\beta s^{*}{i}^{*}+{\beta }_2{s}^{*}{q}^{*}}{\delta e^{*}}\left(i-i^{\text{*}}-i^{\text{*}}\ln \frac{i^{\text{*}}}{i}\right)+\frac{{\beta }_2{s}^{*}{q}^{*}}{m{i}^{*}}\left(q-q^{\text{*}}-q^{\text{*}}\ln \frac{q^{\text{*}}}{q}\right)$

$\begin{array}{c}\frac{\text{d}{L}_2}{\text{d}t}=\left(1-\frac{s^{*}}{s}\right)s^{\prime }+\left(1-\frac{e^{*}}{e}\right)e^{\prime }+\frac{\beta s^{*}{i}^{*}+{\beta }_2{s}^{*}{q}^{*}}{\delta e^{*}}\left(1-\frac{i^{*}}{i}\right)i^{\prime }+\frac{{\beta }_2{s}^{*}{q}^{*}}{m{i}^{*}}\left(1-\frac{q^{*}}{q}\right)q^{\prime }\\ =\left(1-\frac{s^{*}}{s}\right)\left[\mu -\left(\beta si+{\beta }_{1}se+{\beta }_{2}sq\right)-\mu s\right]+\left(1-\frac{e^{*}}{e}\right)\left(\beta si+{\beta }_{1}se+{\beta }_{2}sq-\delta e-\mu e\right)\\ +\frac{\beta s^{*}{i}^{*}+{\beta }_2{s}^{*}{q}^{*}}{\delta e^{*}}\left(1-\frac{i^{*}}{i}\right)\left(\delta e-mi-\mu i\right)+\frac{{\beta }_2{s}^{*}{q}^{*}}{m{i}^{*}}\left(1-\frac{q^{*}}{q}\right)\left(mi-\gamma q-\mu q\right)\\ =\beta s^{*}{i}^{*}\left[\left(1-\frac{s^{*}}{s}\right)\left(1-\frac{si}{s^{*}{i}^{*}}\right)+\left(1-\frac{e^{*}}{e}\right)\left(\frac{si}{s^{*}{i}^{*}}-\frac{e}{e^{*}}\right)\right]+\frac{\beta s^{*}{i}^{*}+{\beta }_2{s}^{*}{q}^{*}}{\delta e^{*}}\delta e^{*}\left(1-\frac{i^{*}}{i}\right)\left(\frac{e}{e^{*}}-\frac{i}{i^{*}}\right)\end{array}$

$\begin{array}{l}+\beta s^{*}{e}^{*}\left[\left(1-\frac{s^{*}}{s}\right)\left(1-\frac{se}{s^{*}{e}^{*}}\right)+\left(1-\frac{e^{*}}{e}\right)\left(\frac{se}{s^{*}{e}^{*}}-\frac{e}{e^{*}}\right)\right]+\mu s^{*}\left(1-\frac{s^{*}}{s}\right)\left(1-\frac{s}{s^{*}}\right)\\ +\beta s^{*}{q}^{*}\left[\left(1-\frac{s^{*}}{s}\right)\left(1-\frac{sq}{s^{*}{q}^{*}}\right)+\left(1-\frac{e^{*}}{e}\right)\left(\frac{sq}{s^{*}{q}^{*}}-\frac{e}{e^{*}}\right)\right]+\frac{{\beta }_2{s}^{*}{q}^{*}}{m{i}^{*}}m{i}^{*}\left(1-\frac{q^{*}}{q}\right)\left(\frac{i}{i^{*}}-\frac{q}{q^{*}}\right)\\ =\beta s^{*}{i}^{*}\left(2-\frac{s^{*}}{s}+\frac{i}{i^{*}}-\frac{e^{*}}{e}-\frac{si{e}^{*}}{s^{*}{i}^{*}e}\right)+\frac{\beta s^{*}{i}^{*}+{\beta }_2{s}^{*}{q}^{*}}{\delta e^{*}}\delta e^{*}\left(1+\frac{e}{e^{*}}-\frac{i^{*}}{i}-\frac{i^{*}e}{i{e}^{*}}\right)\\ +{\beta }_2{s}^{*}{q}^{*}\left(2-\frac{s^{*}}{s}+\frac{q}{q^{*}}-\frac{e^{*}}{e}-\frac{sq{e}^{*}}{s^{*}{q}^{*}e}\right)+\frac{{\beta }_2{s}^{*}{q}^{*}}{m{i}^{*}}m{i}^{*}\left(1+\frac{i}{i^{*}}-\frac{q^{*}}{q}-\frac{i{q}^{*}}{i^{*}q}\right)+\left({\beta }_1{s}^{*}{e}^{*}+\mu s^{*}\right)\left(2-\frac{s^{*}}{s}-\frac{s}{s^{*}}\right)\end{array}$

由于 $2-\frac{s^{*}}{s}-\frac{s}{s^{*}}\le \text{0}$，且用到不等式 $1-a\le -\ln a$， $a>0$，

$\begin{array}{l}\text{1+}\frac{e}{e^{*}}-\frac{i}{i^{*}}-\frac{i^{*}e}{i{e}^{*}}\le \frac{e}{e^{*}}-\ln \frac{e}{e^{*}}+\ln \frac{i}{i^{*}}-\frac{i}{i^{*}},\\ \text{1+}\frac{i}{i^{*}}-\frac{q}{q^{*}}-\frac{i{q}^{*}}{i^{*}q}\le \frac{i}{i^{*}}-\ln \frac{i}{i^{*}}+\ln \frac{q}{q^{*}}-\frac{q}{q^{*}},\end{array}$

$\begin{array}{l}2-\frac{s^{*}}{s}+\frac{i}{i^{*}}-\frac{e^{*}}{e}-\frac{si{e}^{*}}{s^{*}{i}^{*}e}\le -\ln \frac{s}{s^{*}}+\frac{i}{i^{*}}-\frac{e}{e^{*}}-\ln \frac{si{e}^{*}}{s^{*}{i}^{*}e}=\frac{i}{i^{*}}-\ln \frac{i}{i^{*}}+\ln \frac{e}{e^{*}}-\frac{e}{e^{*}}\\ 2-\frac{s^{*}}{s}+\frac{q}{q^{*}}-\frac{e^{*}}{e}-\frac{sq{e}^{*}}{s^{*}{q}^{*}e}\le -\ln \frac{s}{s^{*}}+\frac{q}{q^{*}}-\frac{e}{e^{*}}-\ln \frac{sq{e}^{*}}{s^{*}{q}^{*}e}=\frac{q}{q^{*}}-\ln \frac{q}{q^{*}}+\ln \frac{e}{e^{*}}-\frac{e}{e^{*}}\end{array}$

因此，

$\begin{array}{l}\frac{\text{d}{L}_{\text{2}}}{\text{d}t}\le \left(\beta s^{*}{i}^{*}+{\beta }_2{s}^{*}{q}^{*}\right)\left(\frac{e}{e^{*}}-\ln \frac{e}{e^{*}}-\frac{e}{e^{*}}+\ln \frac{e}{e^{*}}\right)+\left(\beta s^{*}{i}^{*}+{\beta }_2{s}^{*}{q}^{*}\right)\left(\frac{i}{i^{*}}-\ln \frac{i}{i^{*}}-\frac{i}{i^{*}}+\ln \frac{i}{i^{*}}\right)\\ +{\beta }_2{s}^{*}{q}^{*}\left(\frac{q}{q^{*}}-\ln \frac{q}{q^{*}}-\frac{q}{q^{*}}+\ln \frac{q}{q^{*}}\right)+\left({\beta }_1{s}^{*}{e}^{*}+\mu s^{*}\right)\left(2-\frac{s^{*}}{s}-\frac{s}{s^{\text{*}}}\right)\le \text{0}\text{.}\end{array}$

所以，根据 [12] 中定理3.1的定义，函数 $L_2$ 是系统(3)的Lyapunov函数，当且仅当 $s=s^{*},e=e^{*},i=i^{*},q=q^{*}$ 时， ${L^{\prime }}_2=0$ 成立，且对于可行域 $D_1$ 中的所有点，都有 ${L^{\prime }}_2\le 0$。根据LaSalle’s不变原理 [11]，当 $R_0>1$ 时，地方病平衡点 ${\zeta }_2$ 在 $D_1$ 内是全局渐近稳定的。

3.2. 疫苗接种比例为常数

当博弈免疫模型中的免疫比例x为常数时，则每天按固定比例对易感者和新生儿接种，且设这个固定比例为v，此时模型为：

$\{\begin{array}{l}\frac{\text{d}S}{\text{d}t}=\mu \left(1-v\right)N-\frac{\beta SI+{\beta }_{1}SE+{\beta }_{2}SQ}{N}-vS-\mu S,\\ \frac{\text{d}E}{\text{d}t}=\frac{\beta SI+{\beta }_{1}SE+{\beta }_{2}SQ}{N}-\delta E-\mu E,\\ \frac{\text{d}I}{\text{d}t}=\delta E-mI-\mu I,\\ \frac{\text{d}Q}{\text{d}t}=mI-\gamma Q-\mu Q,\\ \frac{\text{d}R}{\text{d}t}=\mu vN+vS+\gamma Q-\mu R.\end{array}$ (4)

与3.1相同，对于模型(4)消除R，并作归一化变换，微分方程为：

$\{\begin{array}{l}\frac{\text{d}s}{\text{d}t}=\mu \left(1-v\right)-\left(\beta si+{\beta }_{1}se+{\beta }_{2}sq\right)-vs-\mu s,\\ \frac{\text{d}e}{\text{d}t}=\beta si+{\beta }_{1}se+{\beta }_{2}sq-\delta e-\mu e,\\ \frac{\text{d}i}{\text{d}t}=\delta e-mi-\mu i,\\ \frac{\text{d}q}{\text{d}t}=mi-\gamma q-\mu q.\end{array}$ (5)

系统(5)的可行域为

$D_{\text{2}}=\left\{\left(s,e,i,q\right)\in R^4|0\le s,e,i,q\le 1,0\le s+e+i+q\le 1\right\}$

且 $D_{\text{2}}$ 为正向不变集。

通过计算，模型存在一个无病平衡点 ${\xi }_{\text{1}}\left(\overline{s^0},\overline{e^0},\overline{i^0},\overline{q^0}\right)=\left(\frac{\mu \left(1-v\right)}{v+\mu },0,0,0\right)$ 和一个地方病平衡点 ${\xi }_{\text{2}}\left(\overline{s^{\text{*}}},\overline{e^{\text{*}}},\overline{i^{\text{*}}},\overline{q^{\text{*}}}\right)$，其中

$\overline{s^{*}}=\frac{1}{R_0},\overline{e^{*}}=\frac{\mu \left(1-v\right)R_0-\mu -v}{\left(\delta +\mu \right)R_0},\overline{i^{*}}=\frac{\delta \left(\mu \left(1-v\right)R_0-\mu -v\right)}{\left(\delta +\mu \right)\left(m+\mu \right)R_0},\overline{q^{*}}=\frac{m\delta \left(\mu \left(1-v\right)R_0-\mu -v\right)}{\left(\delta +\mu \right)\left(\gamma +\mu \right)\left(m+\mu \right)R_0}.$

由于 ${\xi }_{\text{2}}$ 的存在性，可得有效再生数为

$R_v=\left(1-v_0\right)R{}_0,v_0=\frac{v\left(1+\mu \right)}{\mu +v}.$

定理2 1) 当 $R_v<1$ 时，无病平衡点 ${\xi }_1$ 在正向不变集 $D_{\text{2}}$ 内是全局渐近稳定的；

2) 当 $R_v>1$ 时，地方病平衡点 ${\xi }_2$ 在正向不变集 $D_{\text{2}}$ 内是全局渐近稳定的。

证明同定理1类似。

其中Lyapunov函数分别为

$\begin{array}{l}{L}_3=s-\overline{s^0}-\overline{s^0}\ln \frac{\overline{s^0}}{s}+e+\left[\frac{\beta \overline{s^0}}{m+\mu }+\frac{m{\beta }_2\overline{s^0}}{\left(m+\mu \right)\left(\gamma +\mu \right)}\right]i+\frac{{\beta }_2\overline{s^0}}{\gamma +\mu }q\\ L_4=s-\overline{s^{*}}-\overline{s^{*}}\ln \frac{\overline{s^{*}}}{s}+e-\overline{e^{*}}-\overline{e^{*}}\ln \frac{\overline{e^{*}}}{e}+\frac{\beta \overline{s^{*}}\overline{i^{*}}+{\beta }_2\overline{s^{*}}\overline{q^{*}}}{\delta \overline{e^{*}}}\left(i-\overline{i^{*}}-\overline{i^{*}}\ln \frac{\overline{i^{*}}}{i}\right)\\ +\frac{{\beta }_2\overline{s^{*}}\overline{q^{*}}}{m\overline{i^{*}}}\left(q-\overline{q^{*}}-\overline{q^{*}}\ln \frac{\overline{q^{*}}}{q}\right)\end{array}$

3.3. 疫苗接种比例随染病人数变化

当博弈免疫模型中的免疫比例x为随时间按变化时，就是加入博弈免疫的模型(2)，同3.1相同对于模型(2)消除R，微分方程为：

$\{\begin{array}{l}\frac{\text{d}S}{\text{d}t}=\mu \left(1-x\right)N-\frac{\beta SI+{\beta }_{1}SE+{\beta }_{2}SQ}{N}-xS-\mu S,\\ \frac{\text{d}E}{\text{d}t}=\frac{\beta SI+{\beta }_{1}SE+{\beta }_{2}SQ}{N}-\delta E-\mu E,\\ \frac{\text{d}I}{\text{d}t}=\delta E-mI-\mu I,\\ \frac{\text{d}Q}{\text{d}t}=mI-\gamma Q-\mu Q,\\ \frac{\text{d}x}{\text{d}t}=\kappa x\left(1-x\right)\left[-1+\omega \left(E+I+Q\right)\right].\end{array}$ (6)

系统(6)的可行域为

$D_{\text{3}}=\left\{\left(S,E,I,Q,x\right)\in R^5|0\le S,E,I,Q\le N,0\le x\le 1,0\le S+E+I+Q\le N\right\}$

且 $D_{\text{3}}$ 为正向不变集。

当 $x\ne 0$ 时，系统有2个平衡点， $x=1$ 时，有平衡点 ${\epsilon }_{\text{1}}\left(\overline{S_1},\overline{E_1},\overline{I_1},\overline{Q_1},\overline{x_1}\right)=\left(0,0,0,0,1\right)$，此时所有人全部免疫，且 ${\epsilon }_{\text{1}}$ 是全局稳定的。 $0<x<1$ 时，有一部分人选择接种疫苗，且当 $R_{\omega }>1$ 时，系统有一个地方病平衡点 ${\epsilon }_{\text{2}}\left(\overline{S_{\text{2}}},\overline{E_{\text{2}}},\overline{I_{\text{2}}},\overline{Q_{\text{2}}},\overline{x_{\text{2}}}\right)$，其中

$\begin{array}{l}\overline{S_2}=\frac{N}{R_0},\\ \overline{E_2}=\frac{\left(\gamma +\mu \right)\left(m+\mu \right)}{\left[\left(\gamma +\mu \right)\left(m+\mu \right)+\delta \left(\gamma +\mu \right)+\delta m\right]\omega },\\ \overline{I_2}=\frac{\delta \left(\gamma +\mu \right)}{\left[\left(\gamma +\mu \right)\left(m+\mu \right)+\delta \left(\gamma +\mu \right)+\delta m\right]\omega },\\ \overline{Q_2}=\frac{\delta m}{\left[\left(\gamma +\mu \right)\left(m+\mu \right)+\delta \left(\gamma +\mu \right)+\delta m\right]\omega },\\ \overline{x_2}=\frac{\mu \left(R_0-1\right)}{\mu R_0+1}-\frac{\left(\delta +\mu \right)\left(\gamma +\mu \right)\left(m+\mu \right)}{\omega N\left(\mu R_0+1\right)\left[\left(\gamma +\mu \right)\left(m+\mu \right)+\delta \left(\gamma +\mu \right)+\delta m\right]},\end{array}$

模型的有效再生数 $R_{\omega }$ 由 ${\epsilon }_{\text{2}}$ 的存在性得出

$R_{\omega }=\left(1-{\omega }_0\right)R_0,{\omega }_0=\frac{\left(\delta +\mu \right)\left(\gamma +\mu \right)\left(m+\mu \right)}{\mu \omega N\left(\mu R_0+1\right)\left[\left(\gamma +\mu \right)\left(m+\mu \right)+\delta \left(\gamma +\mu \right)+\delta m\right]}.$

4. 数值模拟

4.1. 数据拟合

我们收集了2020年4月1日至2020年7月29日这120天的新增病例和累计病例数据，首先用模型(1)模拟印度2020年4月1日至2020年7月29日期间的病例数。其中，Z(t)代表累积的病例数，Z(t)表示为：

$\frac{\text{d}Z}{\text{d}t}=mI\left(t\right).$

系统(1)的初始条件和其他参数估计如表1和表2所示。用最小二乘法分别对120天印度新冠肺炎的每日新增病例和累计病例进行数据模拟，并与实际数据进行对比如图3所示。此外，通过最小二乘拟合得出参数 $\beta$ 和m的值：

$\beta =0.1439,m=\mathrm{0.2608.}$

4.2. 免疫比例x为常数

将最小二乘拟合得到 $\beta =0.1439,m=0.2608$ 代入到模型(4)，当免疫比例x为常数v时，即每天给一定比例的人接种疫苗，其中图4(a)表示新增病例的变化，图4(b)表示累计病例的变化，初值取最小二乘拟合的末值，三条线颜色不同的线表示每天的疫苗接种比例不同对疾病传播的影响也是不同的。我们可以清楚的看到，按比例接种后，每日新增病例趋势是先增后减，随后减小到0，且保持在0增长；累计病例也会由增加，随后保持在固定值。

当v取不同值时，即每天对印度不同比例的人口接种疫苗，从图4可以看到当v的值越大时，新增病例最大值会越低，且到达最大值的时间越早。同样地，v取值越大时，累计病例最终稳定值会越小，而且到达稳定值的时间也越早。所以，在按比例接种疫苗时，增大疫苗接种比例能有效减少染病人数，且会减少疫情持续时间。

4.3. 疫苗博弈的免疫效果模拟

将最小二乘拟合得到的参数 $\beta =0.1439,m=0.2608$ 代入到模型(2)中，由于在模型(6)中要x的取值为 $0<x<1$，所以取 $\omega =0.000001,k=0.001$，可得如图5所示在博弈免疫策略下接种疫苗后的染病者和疫苗

接种者的变化趋势及其相互关系。同样初值取最小二乘拟合的末值，而免疫比例x的初值设为0.001。从图5可以看到，在开始接种疫苗的一段时间中，染病者的人数持续增加，同时疫苗接种量很低，当染病者增长到某个数量，疫苗接种开始快速增长，而随着疫苗接种的增长，染病者数量开始下降，并且很快下降到0，然后持续在0染病状态，同时，在染病人数降低到一定程度后，疫苗接种量也开始就下降，最终可以下降到0。

在图6、图7中，我们对疫苗博弈策略之后一段时间的每日新增病例数和累计病例数进行了短期预测。首先，从图6中可以清楚地看出，在博弈免疫策略下每日新增病例数量在先增加后减小，然后趋于0。然后，从图7中可以清楚地看出，累计病例接种初期还是呈现增长趋势，但是在一段时间后，累计病例不再增加而是稳定状态。因此，在博弈疫苗策略下，能有效减少染病人数，同时大大缩短疫情的时间。

对于博弈模型中的两个参数 $\omega$ 和 $\kappa$， $\kappa =k{r}_v$， $\omega =\eta \left({r_i/r}_v\right)$，其中k表示组合模仿率，即学习模仿他人的能力， $\eta$ 表示不接种疫苗的个体感染疾病的可能性，k和 $\eta$ 的值会影响疾病的传播，即 $\omega$ 和 $\kappa$ 的值会影响疾病的传播。在图6(a)、图6(b)和图7(a)、图7(b)中，我们对不同 $\omega$ 和 $\kappa$ 值下的累计病例数和每日新增病例数进行了短期预测。从图6(a)和图7(a)可以清楚地看出，当 $\omega$ 值固定为0.000001时，不同值的 $\kappa$ 在博弈免疫初期对新增病例和累计病例的影响不大，而在一段时间的接种后，随着 $\kappa$ 值的增大，新增病例最大值会减小，达到新增病例最大值的时间会越来越早，且达到0增长的时间也会越来越早；同时，随着 $\kappa$ 值的增大，累计病例的最终稳定值会减小，达到稳定值的时间也会越早。所以，人们的学习能力越强(k越大)时，对疾病的控制效果越好。

同样地，从图6(b)和图7(b)可以清楚地看出，当 $\kappa$ 值固定为0.005时，不同值的 $\omega$ 会对博弈免疫产生相同的影响，且不同值的 $\omega$ 在博弈免疫初期对新增病例和累计病例的影响不大，而随着 $\omega$ 值的增大，新增病例最大值会减小，达到最大值的时间会提前，达到0增长的时间也会提前；同时，随着 $\omega$ 值的增大，累计病例的最终稳定值会减小，达到稳定值的时间也会提前。所以，在不接种疫苗时，个体感染疾病的可能性越大( $\eta$ 越大)时，对疾病的控制效果越好。因此，在人们学习能力强(k较大)或者不接种疫苗的个体感染疾病的可能性更大 $\eta$ 较大)时，会促使人们接种疫苗，更快控制疾病且达到消灭疾病。

4.4. 基本再生数R₀和有效再生数R_v, R_ω的比较

通过计算和动力学分析，可以得到在不考虑疫苗时，模型的基本再生数为1.5641，且有效再生数 $R_v,R_{\omega }$ 的表达式如下表3，其中

$v_0=\frac{v\left(1+\mu \right)}{\mu +v},{\omega }_0=\frac{\left(\delta +\mu \right)\left(\gamma +\mu \right)\left(m+\mu \right)}{\mu \omega N\left(\mu R_0+1\right)\left[\left(\gamma +\mu \right)\left(m+\mu \right)+\delta \left(\gamma +\mu \right)+\delta m\right]}.$

通过分析可知，在不考虑疫苗时，模型的基本再生数大于1，表明此时如果不采取控制措施，疾病将会继续传播。而对于免疫比例为常数v时，由有效再生数 $R_v=1$ 可得此时出生率和免疫率的比例，即 $\mu /v=1.7729$，所以免疫率v与出生率 $\mu$ 有关。当出生率和免疫率的比例大于某个固定值时，有效再生数 $R_v$ 会小于1，此时疾病传播将会得到有效控制，疾病不会继续传播；而当出生率和免疫率的比例小于某个固定值时，有效再生数 $R_v$ 大于1，此时疾病将继续传播。且当免疫率v越大时，染病人数会大规模减少，且疫情持续时间也会缩短很多。如图4所示，免疫比例v从0.005增大到0.01，新增病例的最大值几乎减小一半，且达到最大值的时间也提前将近50天，新增病例减小到0的时间提前将近100天，而累计病例的最终规模减小了1.5倍，到达最终稳定值的时间也提前将近100天。这表明了按比例接种疫苗接种对控制疾病传播的有效性，当出生率和免疫率的比例大于某个固定值，能够有效减少染病人数，且缩短疫情持续时间。

而对于博弈免疫的有效再生数，由 $R_{\omega }=1$ 可得 $\omega$ 的阈值为0.000001517，且当 $\omega$ 大于这个阈值时，有 $R_{\omega }$ 大于1；而当 $\omega$ 小于这个阈值时，有效再生数 $R_{\omega }$ 小于1。且 $\omega$ 越大，有效再生数 $R_{\omega }$ 越大，即不接种疫苗的个体感染疾病的可能性更大时( $\eta$ 较大)，有效再生数 $R_{\omega }$ 越大。而疾病的控制不止与不接种疫苗的个体感染疾病的可能性( $\eta$ )有关，还与人们学习能力(k)有关，从图6、图7可以看出，在人们学习能力强(k较大)或者不接种疫苗的个体感染疾病的可能性更大( $\eta$ 较大)时，新增病例最大值和累计病例最终规模都会减少，染病人数大量减少，且疫情持续时间大大缩短。

5. 结论

本文针对印度实际情况，建立SEIQR模型，在此基础上加入博弈免疫策略，并建立相应模型，计算其基本再生数 $R_{\text{0}}$ 、按比例接种疫苗的有效再生数 $R_v$ 和疫苗博弈策略的有效再生数 $R_{\omega }$。结论表明印度新冠在2020年4月1日到7月29日的基本再生数为1.5941，考虑加入免疫，若按比例接种疫苗，有效再生数 $R_v$ 与出生率和免疫率的比例有关，比例大于某一固定值时，有效再生数 $R_v$ 大于1，此时疾病将继续传播；比例小于某一固定值时，有效再生数 $R_v$ 小于1，此时疾病传播将会得到有效控制，疾病将不会继续传播。且当免疫比例v越大时，染病人数减少越多，且疫情持续时间也会缩短越多。若按照博弈免疫策略进行疫苗接种，即在自愿政策下进行免疫，这时有效再生数 $R_{\omega }$ 与参数 $\omega$ 有关，有效再生数 $R_{\omega }=1$ 可得 $\omega$ 的阈值，当 $\omega$ 大于这个阈值时，有效再生数 $R_{\omega }$ 大于1，此时地方病平衡点存在，而当 $\omega$ 小于这个阈值时，有效再生数 $R_{\omega }$ 小于1。且 $\omega$ 越大，有效再生数 $R_{\omega }$ 越大，即不接种疫苗的个体感染疾病的可能性更大时( $\eta$ 较大)，有效再生数 $R_{\omega }$ 越大。根据数值模拟的结果，在人们学习能力强(k较大)或者不接种疫苗的个体感染疾病的可能性更大( $\eta$ 较大)时，对于个体，此时接种疫苗的收益比不接种疫苗的收益更大，这会促使人们寻求疫苗接种，从而能够更快地减少染病人数，同时大大缩短疫情的持续时间，更好控制疾病传播且达到消灭疾病的目的。所以在自愿接种的政策下，考虑经济因素的疫苗博弈策略，既能够有效控制疾病传播，又能保障人民自身利益。

致谢

本文作者衷心感谢审稿人的意见和建议。

基金项目

本文受到国家自然科学基金(批准号：11801398)和山西省应用基础研究面上青年项目(批准号：201801D221024)的资助。

参考文献