带有疫苗接种的分数阶SEIQR传染病斑块模型分析

Analysis of the Fractional-Order SEIQR Patch Model with Vaccination for Infections Diseases

摘要

根据麻疹病毒的传播机制，建立了一类分数阶SEIQR传染病斑块模型，研究了疫苗接种和迁移运动对传染病控制的影响。利用下一代矩阵法得到了基本再生数R_{\mathrm{0}}，根据Mittag-Leffler函数和Laplace变换得出了分数阶模型的正向不变集以及阈值条件下平衡点的稳定性。最后针对两个斑块进行数值模拟，结果表明迁移运动可以使疾病消失，同时增加疫苗接种率能够减小基本再生数，以及疾病趋于稳定的时间会随着分数阶参数α的减小而延长。

Abstract

According to the transmission mechanism of measles virus， a fractional order SEIQR patch epidemic model was established to study the effects of vaccination and migration on the control of infectious diseases. The basic reproducing number R_{\mathrm{0}} was obtained by using the next generation matrix method， and the forward invariant set of the fractional order model and the stability of the equilibrium point under threshold conditions were obtained based on the Mittag-Leffler function and Laplace transform. Finally， by selecting two patches for numerical simulation， it is found that the migration movement can make the disease disappear， at the same time， increasing the vaccination rate can reduce the basic reproduction number， and it is found that the time for the disease to stabilize increases as α decreases.

Graphical abstract

关键词

疫苗接种 / SEIQR传染病斑块模型 / 基本再生数 / 分数阶 / Caputo 导数

Key words

vaccination / SEIQR patch model for infections disease / basic reproductive number / fractional order / Caputo derivative

麻疹，也被称为红疹或莫比利，是一种由副黏病毒科的莫比利病毒引起的传染病^［1］。在现代医学的发展下，麻疹可以通过接种疫苗预防，但在一些疫苗接种率较低的地区，麻疹依然对公共卫生安全造成了严重威胁。目前已经有许多学者研究了麻疹病毒，Aldila等 ^［2］为了控制麻疹在人群中的传播，实施了两次疫苗接种的干预措施，建立了麻疹感染的改良\mathrm{S}\mathrm{V}\mathrm{I}\mathrm{Q}\mathrm{R}（易感-接种疫苗-感染-隔离-康复）确定性模型。Berhe等^［3］提出了一个用微分方程稳定性理论进行分析的麻疹传播动力学的模型，描述了麻疹疫情在人群中的传播动态。Memon等^［4］建立\mathrm{S}\mathrm{V}\mathrm{E}\mathrm{I}\mathrm{R}模型，使用2019年1—10月报告的巴基斯坦的真实麻疹病例进行模拟，得出提高疫苗效力和覆盖率可以大幅降低麻疹流行风险的结论。

分数阶微分方程作为传统整数阶微分方程的自然延伸，已经引起了广泛的关注，因为它们能够更准确地描述许多非局部和非线性现象。分数阶微分方程在许多领域都得到了应用，比如量子力学、信号与图像处理、经济学等^［5-8］领域。分数阶导数的定义是以积分的形式表示，分数阶微分方程不仅考虑当前时刻状态还考虑了之前时刻状态，引入了记忆和遗传特性，可以更好地描述传染病在时间上的动力学行为，尤其是具有长时间记忆效应的疾病传播，因此也被应用于传染病动力学分析^［9-10］。相比于整数阶传染病模型，分数阶传染病模型可以更好地刻画生物的遗传性和记忆性特征，并且克服了整数阶传染病模型无法完全拟合数据的缺点，因此，分数阶传染病模型的研究具有重要的现实意义。

许多研究人员意识到迁移运动会快速将传染病从一个地方传播到另一个地方，因此他们对斑块传染病进行了研究。Driessche等^［11］基于具有斑块的空间异质环境，建立了疾病传播模型，结果表明，感染者在斑块环境中的移动可能对疾病传播有重要影响。Cui等^［12］建立了一类具有疫苗接种的\mathrm{S}\mathrm{I}\mathrm{R}斑块模型，研究了人口流动性增加和疫苗接种对疾病传播的影响。

在预防和控制疾病传播的过程中，接种疫苗和隔离是经常采用的两种手段。基于上述考虑以及病毒传播的记忆性特征，使用具有记忆效应的分数阶微分系统刻画传染病模型比整数阶微分系统更为精确。因此，本文根据麻疹病毒的传播机制，建立了一类带有疫苗接种的分数阶\mathrm{S}\mathrm{E}\mathrm{I}\mathrm{Q}\mathrm{R}传染病斑块模型，研究在不同的分数阶参数\alpha下，迁移运动和疫苗接种对控制疾病传播的影响。

1 分数阶SEIQR传染病斑块模型的建立

根据麻疹病毒的传播机制建立了一类基于Caputo分数阶导数的\mathrm{S}\mathrm{E}\mathrm{I}\mathrm{Q}\mathrm{R}传染病斑块模型，将每个斑块的人群划分为5个仓室：易感者类S_i\left(t\right)，潜伏者类E_i\left(t\right)，感染者类I_i\left(t\right)，隔离者类Q_i\left(t\right)和恢复者类R_i\left(t\right)。假设隔离者不能在斑块间流动，其他子群体都能够在任意两个斑块间流动，斑块i中的疾病传播情况如图1所示。

考虑病毒传播时记忆和遗传效应的影响，疾病传播的流程图所对应的分数阶\mathrm{S}\mathrm{E}\mathrm{I}\mathrm{Q}\mathrm{R}传染病斑块的数学模型为

\left\{\begin{array}{c}{}_{\mathrm{0}}{}^c{D}_t^{\alpha }{S}_i={\Lambda }_i-{\beta }_i{S}_i{I}_i-(p_i{e}_i+{\mu }_i)S_i+{\displaystyle \sum _{j\ne i}^r}\left(a_{ij}{S}_j-a_{ji}{S}_i\right)\\ {}_{\mathrm{0}}{}^c{D}_t^{\alpha }{E}_i={\beta }_i{S}_i{I}_i-\left({\sigma }_i+{\mu }_i\right)E_i+{\displaystyle \sum _{j\ne i}^r}\left(b_{ij}{E}_j-b_{ji}{E}_i\right)\\ {}_{\mathrm{0}}{}^c{D}_t^{\alpha }{I}_i={\sigma }_i{E}_i-({\gamma }_i+{\epsilon }_i+{\mu }_i)I_i+{\displaystyle \sum _{j\ne i}^r}\left(c_{ij}{I}_j-c_{ji}{I}_i\right)\\ {}_{\mathrm{0}}{}^c{D}_t^{\alpha }{Q}_i={\epsilon }_i{I}_i-\left({\mu }_i+h_i\right)Q_i\\ {}_{\mathrm{0}}{}^c{D}_t^{\alpha }{R}_i={\gamma }_i+p_i{e}_i{S}_i+h_i{Q}_i-{\mu }_i{R}_i+{\displaystyle \sum _{j\ne i}^r}\left(d_{ij}{R}_j-d_{ji}{R}_i\right)\end{array}\right.（1）

式中，{\Lambda }_i为第i个斑块中新生儿的数量，i=\mathrm{1,2},\cdots ,r；{\beta }_i为第i个斑块中有效接触率；{\mu }_i为第i个斑块中自然死亡率；p_i为第i个斑块中疫苗接种率；e_i为第i个斑块中疫苗接种有效率；{\sigma }_i为第i个斑块中潜伏者到感染者的转换率；{\gamma }_i为第i个斑块中感染者的恢复率；{\epsilon }_i为第i个斑块中感染者的隔离率；h_i为第i个斑块中隔离者的恢复率；a_{ij}、b_{ij}、c_{ij}和d_{ij}分别为第j个斑块中易感者、潜伏者、感染者和恢复者到第i个斑块中的迁移率；\alpha \in (\mathrm{0,1}]对应Caputo分数阶阶次。以上所有参数都为正数。

由于R_i只出现在系统（1）的最后一个方程中，所以只考虑分数阶\mathrm{S}\mathrm{E}\mathrm{I}\mathrm{Q}\mathrm{R}模型的简化系统：

\left\{\begin{array}{l}{}_{\mathrm{0}}{}^c{D}_t^{\alpha }{S}_i={\Lambda }_i-{\beta }_i{S}_i{I}_i-(p_i{e}_i+{\mu }_i)S_i+{\displaystyle \sum _{j\ne i}^r}\left(a_{ij}{S}_j-a_{ji}{S}_i\right)\\ {}_{\mathrm{0}}{}^c{D}_t^{\alpha }{E}_i={\beta }_i{S}_i{I}_i-\left({\sigma }_i+{\mu }_i\right)E_i+{\displaystyle \sum _{j\ne i}^r}\left(b_{ij}{E}_j-b_{ji}{E}_i\right)\\ {}_{\mathrm{0}}{}^c{D}_t^{\alpha }{I}_i={\sigma }_i{E}_i-({\gamma }_i+{\epsilon }_i+{\mu }_i)I_i+{\displaystyle \sum _{j\ne i}^r}\left(c_{ij}{I}_j-c_{ji}{I}_i\right)\\ {}_{\mathrm{0}}{}^c{D}_t^{\alpha }{Q}_i={\epsilon }_i{I}_i-\left({\mu }_i+h_i\right)Q_i\end{array}\right.（2）

2 平衡点的存在性

定理1 系统（2）有唯一的无病平衡点。

证明 根据无病平衡点的定义，将无病状态的情形I_i=\mathrm{0}\left(i=\mathrm{1,2},\cdots ,r\right)代入系统（2）中有

\left\{\begin{array}{l}{\Lambda }_i-(p_i{e}_i+{\mu }_i)S_i+{\displaystyle \sum _{j\ne i}^r}\left(a_{ij}{S}_j-a_{ji}{S}_i\right)=\mathrm{0}\\ -\left({\sigma }_i+{\mu }_i\right)E_i+{\displaystyle \sum _{j\ne i}^r}\left(b_{ij}{E}_j-b_{ji}{E}_i\right)=\mathrm{0}\\ {\sigma }_i{E}_i=\mathrm{0}\\ -\left({\mu }_i+h_i\right)Q_i=\mathrm{0}\end{array}\right.（3）

那么，可以得到E_i=Q_i=\mathrm{0}\left(i=\mathrm{1,2},\cdots ,r\right)且

{\Lambda }_i-(p_i{e}_i+{\mu }_i)S_i+{\displaystyle \sum _{j\ne i}^r}\left(a_{ij}{S}_j-a_{ji}{S}_i\right)=\mathrm{0}（4）

将式（4）化成矩阵的形式

\mathit{J}\mathit{S}=\mathit{\Lambda }（5）

其中

\mathit{J}=\left[\begin{array}{cccc}{\mu }_{\mathrm{1}}+p_{\mathrm{1}}{e}_{\mathrm{1}}+{\displaystyle \sum _{j\ne \mathrm{1}}^r}{a}_{j\mathrm{1}}& -a_{\mathrm{12}}& \cdots & -a_{\mathrm{1}r}\\ -a_{\mathrm{21}}& {\mu }_{\mathrm{2}}+p_{\mathrm{2}}{e}_{\mathrm{2}}+{\displaystyle \sum _{j\ne \mathrm{2}}^r}{a}_{j\mathrm{2}}& \cdots & -a_{\mathrm{2}r}\\ ⋮& ⋮& & ⋮\\ -a_{r\mathrm{1}}& -a_{r\mathrm{2}}& \cdots & {\mu }_r+p_r{e}_r+{\displaystyle \sum _{j\ne r}^r}{a}_{jr}\end{array}\right]，

\mathit{\Lambda }={\left({\Lambda }_{\mathrm{1}},{\Lambda }_{\mathrm{2}},\cdots ,{\Lambda }_r\right)}^{\mathrm{T}}，\mathit{S}={\left(S_{\mathrm{1}},S_{\mathrm{2}},\cdots ,S_r\right)}^{\mathrm{T}}，

由式（5）可得

\mathit{S}={\mathit{J}}^{-\mathrm{1}}\mathit{\Lambda }。

由于\mathit{J}的所有非对角线元素都非正且\mathit{J}的每一列的和都为正数，那么根据文献［13］可得\mathit{J}是一个非奇异\mathit{M}矩阵，所以{\mathit{J}}^{-\mathrm{1}}\ge \mathrm{0}。因此，\mathit{J}可以写成\mathit{J}=K\mathit{I}-\mathit{B}，其中\mathit{I}是单位矩阵，\mathit{B}是非负矩阵，且K>\rho \left(\mathit{B}\right)，\rho \left(\mathit{B}\right)是\mathit{B}的谱半径，所以系统（2）有唯一的正解{\mathit{S}}^{\mathrm{0}}=\left(S_{\mathrm{1}}^{\mathrm{0}},S_{\mathrm{2}}^{\mathrm{0}},\cdots ,S_r^{\mathrm{0}}\right)={\mathit{J}}^{-\mathrm{1}}\mathit{\Lambda }>\mathrm{0}，以及唯一的无病平衡点E_{\mathrm{0}}=\left(S_{\mathrm{1}}^{\mathrm{0}},\mathrm{0,0},\mathrm{0},S_{\mathrm{2}}^{\mathrm{0}},\mathrm{0,0},\mathrm{0},\cdots ,S_r^{\mathrm{0}},\mathrm{0,0},\mathrm{0}\right)。

3 正向不变集

定理2 系统（2）的一个正向不变集为

，

其中N={\displaystyle \sum _{i=\mathrm{1}}^r}\left(S_i+E_i+I_i+Q_i\right)，\mu =\mathrm{m}\mathrm{i}\mathrm{n}\left\{{\mu }_{\mathrm{1}},{\mu }_{\mathrm{2}},\cdots ,{\mu }_r\right\}，\overline{\Lambda }={\displaystyle \sum _{i=\mathrm{1}}^r}{\Lambda }_i，\overline{{\rm N}}=\mathrm{m}\mathrm{a}\mathrm{x}\left\{N_{\mathrm{0}},\frac{\overline{\Lambda }}{\mu }\right\}，N_{\mathrm{0}}=N\left(\mathrm{0}\right)。

证明 由文献［14］可得系统（2）在\left(\mathrm{0},+\mathrm{\infty }\right)上解的存在性和唯一性。首先证明S_i>\mathrm{0},N>\mathrm{0},i=\mathrm{1,2},\cdots ,r。

令{\Omega }_{+}=\left\{\left(S_{\mathrm{1}},I_{\mathrm{1}},E_{\mathrm{1}},Q_{\mathrm{1}},\cdots ,S_r,I_r,E_r,Q_r\right)\left|S_i\right.>\mathrm{0},I_i\ge \mathrm{0},Q_i\ge \mathrm{0},E_i\ge \mathrm{0},i=\mathrm{1,2},\cdots ,r\right\}，在{\Omega }_{+}的超平面上有

\left\{\begin{array}{l}{}_{\mathrm{0}}{}^c{D}_t^{\alpha }{S}_i\left|{}_{S_i=\mathrm{0}}\right.={\Lambda }_i+{\displaystyle \sum _{j\ne i}^r}{a}_{ij}{S}_j\ge \mathrm{0}\\ {}_{\mathrm{0}}{}^c{D}_t^{\alpha }{E}_i\left|{}_{E_i=\mathrm{0}}\right.={\beta }_i{S}_i{I}_i+{\displaystyle \sum _{j\ne i}^r}{b}_{ij}{E}_j\ge \mathrm{0}\\ {}_{\mathrm{0}}{}^c{D}_t^{\alpha }{I}_i\left|{}_{I_i=\mathrm{0}}\right.={\sigma }_i{E}_i+{\displaystyle \sum _{j\ne i}^r}{c}_{ij}{I}_j\ge \mathrm{0}\\ {}_{\mathrm{0}}{}^c{D}_t^{\alpha }{Q}_i\left|{}_{Q_i=\mathrm{0}}\right.={\epsilon }_i{I}_i\ge \mathrm{0}\end{array}\right.（6）

由式（6）及文献［15］得，若\left\{S_{\mathrm{1}}\left(\mathrm{0}\right),E_{\mathrm{1}}\left(\mathrm{0}\right),I_{\mathrm{1}}\left(\mathrm{0}\right),Q_{\mathrm{1}}\left(\mathrm{0}\right),\cdots ,S_r\left(\mathrm{0}\right),E_r\left(\mathrm{0}\right),I_r\left(\mathrm{0}\right),Q_r\left(\mathrm{0}\right)\right\}\in {\Omega }_{+}，那么系统（2）的解\left(S_{\mathrm{1}},E_{\mathrm{1}},I_{\mathrm{1}},Q_{\mathrm{1}},\cdots ,S_r,E_r,I_r,Q_r\right)不会脱离S_i=\mathrm{0},E_i=\mathrm{0},I_i=\mathrm{0},Q_i=\mathrm{0}的超平面，其中每个超平面上的向量场与该平面相切且指向区域{\Omega }_{+}的内部，即系统（2）的解总保持在区域{\Omega }_{+}中，因此得证\mathit{S}>\mathrm{0},N>\mathrm{0},i=\mathrm{1,2},\cdots ,r。

接下来证明\mathit{S}\le {\mathit{S}}^{\mathrm{0}},N\le \overline{N},i=\mathrm{1,2},\cdots ,r。系统（2）的4个式子相加得到

{}_{\mathrm{0}}{}^C{}_t^{\alpha }N\left(t\right)\le \overline{\Lambda }-\mu N（7）

对式（7）进行Caputo分数阶导数的Laplace变换^［16］，可得

s^{\alpha }L\left(N\right)-s^{\alpha -\mathrm{1}}N\left(\mathrm{0}\right)\le \frac{\overline{\Lambda }}{s}-\mu L\left(N\right)，

L\left(N\right)\le \frac{s^{-\mathrm{1}}}{s^{\alpha }+\mu }\overline{\Lambda }+\frac{s^{\alpha -\mathrm{1}}}{s^{\alpha }+\mu }N\left(\mathrm{0}\right)（8）

若\left\{S_{\mathrm{1}}\left(\mathrm{0}\right),E_{\mathrm{1}}\left(\mathrm{0}\right),I_{\mathrm{1}}\left(\mathrm{0}\right),Q_{\mathrm{1}}\left(\mathrm{0}\right),\cdots ,S_r\left(\mathrm{0}\right),E_r\left(\mathrm{0}\right),I_r\left(\mathrm{0}\right),Q_r\left(\mathrm{0}\right)\right\}\in \Omega，那么根据文献［17］中的Mittag-Leffler函数及其Laplace变换，式（8）可化成

\begin{array}{l}N\left(t\right)\le \overline{\Lambda }{t}^{\alpha }{E}_{\alpha ,\alpha +\mathrm{1}}\left(-\mu t^{\alpha }\right)+E_{\alpha ,\mathrm{1}}\left(-\mu t^{\alpha }\right)N\left(\mathrm{0}\right)\\ \le \overline{N}\left(\mu E_{\alpha ,\alpha +\mathrm{1}}\left(-\mu t^{\alpha }\right)+E_{\alpha ,\mathrm{1}}\left(-\mu t^{\alpha }\right)\right)=\overline{N}\frac{\mathrm{1}}{\mathrm{\Gamma }\left(\mathrm{1}\right)}=\overline{N}\end{array}（9）

由系统（2）的第1个式子得

{}_{\mathrm{0}}{}^C{}_t^{\alpha }\mathit{S}\le \mathit{\Lambda }-\mathit{J}\mathit{S}=\mathit{J}{\mathit{S}}^{\mathrm{0}}-\mathit{J}\mathit{S}（10）

对式（10）进行Caputo分数阶导数的Laplace变换^［16］，可得

s^{\alpha }L\left(\mathit{S}\right)-s^{\alpha -\mathrm{1}}\mathit{S}\left(\mathrm{0}\right)\le \mathit{J}{\mathit{S}}^{\mathrm{0}}{s}^{-\mathrm{1}}-\mathit{J}L\left(\mathit{S}\right)，

L(\mathit{S})\le s^{-\mathrm{1}}(s^{\alpha }{\mathit{I}+\mathit{J})}^{-\mathrm{1}}\mathit{J}{\mathit{S}}^{\mathrm{0}}+s^{\alpha -\mathrm{1}}(s^{\alpha }{\mathit{I}+\mathit{J})}^{-\mathrm{1}}\mathit{S}(\mathrm{0})（11）

其中\mathit{I}是单位矩阵。若\left\{S_{\mathrm{1}}\left(\mathrm{0}\right),E_{\mathrm{1}}\left(\mathrm{0}\right),I_{\mathrm{1}}\left(\mathrm{0}\right),Q_{\mathrm{1}}\left(\mathrm{0}\right),\cdots ,S_r\left(\mathrm{0}\right),E_r\left(\mathrm{0}\right),I_r\left(\mathrm{0}\right),Q_r\left(\mathrm{0}\right)\right\}\in \Omega，那么由文献［17］中的Mittag-Leffler函数及其Laplace变换，式（11）可化成

\begin{array}{l}\mathit{S}\left(t\right)\le t^{\alpha }{E}_{\alpha ,\alpha +\mathrm{1}}\left(-\mathit{J}{t}^{\alpha }\right)\mathit{J}{\mathit{S}}^{\mathrm{0}}+E_{\alpha ,\mathrm{1}}\left(-\mathit{J}{t}^{\alpha }\right)\mathit{S}\left(\mathrm{0}\right)\\ \le \left(t^{\alpha }{E}_{\alpha ,\alpha +\mathrm{1}}\left(-\mathit{J}{t}^{\alpha }\right)+E_{\alpha ,\mathrm{1}}\left(-\mathit{J}{t}^{\alpha }\right)\right)={\mathit{S}}^{\mathrm{0}}\frac{\mathrm{1}}{\mathrm{\Gamma }\left(\mathrm{1}\right)}={\mathit{S}}^{\mathrm{0}}\end{array}（12）

由式（9）和式（12）得\mathit{S}\le {\mathit{S}}^{\mathrm{0}},N\le \overline{N},i=\mathrm{1,2},\cdots ,r。

综上所述，\left\{S_{\mathrm{1}}\left(t\right),E_{\mathrm{1}}\left(t\right),I_{\mathrm{1}}\left(t\right),Q_{\mathrm{1}}\left(t\right),\cdots ,S_r\left(t\right),E_r\left(t\right),I_r\left(t\right),Q_r\left(t\right)\right\}\in \Omega。因此，\Omega是一个正不变集。

4 基本再生数

本节计算系统（2）的基本再生数。通过下一代矩阵法^［18］可得

\mathit{F}\left|{}_{E_{\mathrm{0}}}\right.=\left[\begin{array}{ccc}\mathrm{0}& {\mathit{F}}_{\mathrm{12}}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& \mathrm{0}\end{array}\right]，\mathit{V}\left|{}_{E_{\mathrm{0}}}\right.=\left[\begin{array}{ccc}{\mathit{V}}_{\mathrm{11}}& \mathrm{0}& \mathrm{0}\\ {\mathit{V}}_{\mathrm{21}}& {\mathit{V}}_{\mathrm{22}}& \mathrm{0}\\ \mathrm{0}& {\mathit{V}}_{\mathrm{32}}& {\mathit{V}}_{\mathrm{33}}\end{array}\right]，

其中

{\mathit{F}}_{\mathrm{12}}=\mathrm{d}\mathrm{i}\mathrm{a}\mathrm{g}\left({\beta }_{\mathrm{1}}{S}_{\mathrm{1}}^{\mathrm{0}},{\beta }_{\mathrm{2}}{S}_{\mathrm{2}}^{\mathrm{0}},\cdots ,{\beta }_r{S}_r^{\mathrm{0}}\right)，{\mathit{V}}_{\mathrm{21}}=\mathrm{d}\mathrm{i}\mathrm{a}\mathrm{g}\left(-{\sigma }_{\mathrm{1}},-{\sigma }_{\mathrm{2}},\cdots ,-{\sigma }_r\right)，

{\mathit{V}}_{\mathrm{32}}=\mathrm{d}\mathrm{i}\mathrm{a}\mathrm{g}\left(-{\epsilon }_{\mathrm{1}},-{\epsilon }_{\mathrm{2}},\cdots ,-{\epsilon }_r\right)，{\mathit{V}}_{\mathrm{33}}=\mathrm{d}\mathrm{i}\mathrm{a}\mathrm{g}\left({\mu }_{\mathrm{1}}+h_{\mathrm{1}},{\mu }_{\mathrm{2}}+h_{\mathrm{2}},\cdots ,{\mu }_r+h_r\right)，

{\mathit{V}}_{\mathrm{11}}=\left[\begin{array}{cccc}{\sigma }_{\mathrm{1}}+{\mu }_{\mathrm{1}}+{\displaystyle \sum _{j\ne \mathrm{1}}^r}{b}_{j\mathrm{1}}& -b_{\mathrm{12}}& \cdots & -b_{\mathrm{1}r}\\ -b_{\mathrm{21}}& {\sigma }_{\mathrm{2}}+{\mu }_{\mathrm{2}}+{\displaystyle \sum _{j\ne \mathrm{2}}^r}{b}_{j\mathrm{2}}& \cdots & -b_{\mathrm{2}r}\\ ⋮& ⋮& & ⋮\\ -b_{r\mathrm{1}}& -b_{r\mathrm{2}}& \cdots & {\sigma }_r+{\mu }_r+{\displaystyle \sum _{j\ne r}^r}{b}_{jr}\end{array}\right]，

{\mathit{V}}_{\mathrm{22}}=\left[\begin{array}{cccc}{\gamma }_{\mathrm{1}}+{\epsilon }_{\mathrm{1}}+{\mu }_{\mathrm{1}}+{\displaystyle \sum _{j\ne \mathrm{1}}^r}{c}_{j\mathrm{1}}& -c_{\mathrm{12}}& \cdots & -c_{\mathrm{1}r}\\ -c_{\mathrm{21}}& {\gamma }_{\mathrm{2}}+{\epsilon }_{\mathrm{2}}+{\mu }_{\mathrm{2}}+{\displaystyle \sum _{j\ne \mathrm{2}}^r}{c}_{j\mathrm{2}}& \cdots & -c_{\mathrm{2}r}\\ ⋮& ⋮& & ⋮\\ -c_{r\mathrm{1}}& -c_{r\mathrm{2}}& \cdots & {\gamma }_r+{\epsilon }_r+{\mu }_r+{\displaystyle \sum _{j\ne r}^r}{c}_{jr}\end{array}\right]，

那么

\mathit{F}{\mathit{V}}^{-\mathrm{1}}=\left[\begin{array}{ccc}-{\mathit{F}}_{\mathrm{12}}{\mathit{V}}_{\mathrm{11}}^{-\mathrm{1}}{\mathit{V}}_{\mathrm{21}}{\mathit{V}}_{\mathrm{22}}^{-\mathrm{1}}& {\mathit{F}}_{\mathrm{12}}{\mathit{V}}^{-\mathrm{1}}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& \mathrm{0}\end{array}\right]。

因此

R_{\mathrm{0}}=\rho \left(\mathit{F}{\mathit{V}}^{-\mathrm{1}}\right)=\rho \left(-{\mathit{F}}_{\mathrm{12}}{\mathit{V}}_{\mathrm{11}}^{-\mathrm{1}}{\mathit{V}}_{\mathrm{21}}{\mathit{V}}_{\mathrm{22}}^{-\mathrm{1}}\right)。

所以，在斑块与斑块之间无迁移运动时的基本再生数是

R_{\mathrm{0}}^{\left(i\right)}=\frac{{\beta }_i{S}_i^{\mathrm{0}}{\sigma }_i}{\left({\sigma }_i+{\mu }_i\right)\left({\gamma }_i+{\mu }_i+{\epsilon }_i\right)}。

5 无病平衡点的稳定性分析

5.1 无病平衡点的局部稳定性

定理3 当时，无病平衡点E_{\mathrm{0}}在\Omega内局部渐近稳定；当R_{\mathrm{0}}>\mathrm{1}时，无病平衡点E_{\mathrm{0}}不稳定。

证明 系统（2）的Jacobian矩阵为

\mathit{A}\left|{}_{E_{\mathrm{0}}}\right.=\left[\begin{array}{cc}-\mathit{J}& \left[\begin{array}{cc}{\mathit{F}}_{\mathrm{11}}& {\mathit{F}}_{\mathrm{12}}\end{array}\right]\\ \mathrm{0}& \mathit{F}-\mathit{V}\end{array}\right]，

则\mathit{A}\left|{}_{E_{\mathrm{0}}}\right.的特征值是矩阵-\mathit{J}和\mathit{F}-\mathit{V}的特征值。又因为-\mathit{J}是一个非奇异 M -矩阵，所以-\mathit{J}的所有特征值具有负实部，因此，E_{\mathrm{0}}的局部稳定性取决于矩阵\mathit{F}-\mathit{V}的特征值。由于\mathit{V}的所有非对角线元素都非正，且\mathit{V}的每一列的元素之和都非负，\mathit{F}是一个非负矩阵，所以\mathit{F}-\mathit{V}的谱半径，即\mathit{F}-\mathit{V}的所有特征值都是负的，根据文献［17］得到，在时，系统（2）在E_{\mathrm{0}}处是局部稳定的。

若R_{\mathrm{0}}>\mathrm{1}，则s\left(\mathit{F}-\mathit{V}\right)>\mathrm{0}，因此至少有一个\mathit{F}-\mathit{V}的特征值有正实部，系统（2）在E_{\mathrm{0}}处是不稳定的。

5.2 无病平衡点的全局稳定性

定理4 设\mathit{B}=\left(b_{ij}\right)和\mathit{C}=\left(c_{ij}\right)是不可约矩阵。若时，无病平衡点E_{\mathrm{0}}在\Omega内全局渐近稳定；若R_{\mathrm{0}}>\mathrm{1}时，无病平衡点E_{\mathrm{0}}不稳定。

证明 设条件为\left\{S_{\mathrm{1}}\left(\mathrm{0}\right),I_{\mathrm{1}}\left(\mathrm{0}\right),E_{\mathrm{1}}\left(\mathrm{0}\right),Q_{\mathrm{1}}\left(\mathrm{0}\right),\cdots ,\right.\left.S_r\left(\mathrm{0}\right),I_r\left(\mathrm{0}\right),E_r\left(\mathrm{0}\right),Q_r\left(\mathrm{0}\right)\right\}的分数阶系统（2）的解为\left(S_{\mathrm{1}},I_{\mathrm{1}},E_{\mathrm{1}},Q_{\mathrm{1}},\cdots ,S_r,I_r,E_r,Q_r\right)，由文献［17］可得系统（2）在\Omega边界上存在唯一的平衡点E_{\mathrm{0}}，因此，为了得到定理4的证明，只需要证明当t\to \mathrm{\infty }时，系统（2）的正解趋向E_{\mathrm{0}}。

因为S_i\le S_i^{\mathrm{0}}，根据系统（2）的第2个式子可得

{}_{\mathrm{0}}{}^c{}_t^{\alpha }{E}_i\le {\beta }_i{S}_i^{\mathrm{0}}{I}_i-\left({\sigma }_i+{\mu }_i\right)E_i+{\displaystyle \sum _{j\ne i}^r}\left(b_{ij}{E}_j-b_{ji}{E}_i\right)（13）

根据不等式（13）的性质和系统（2）的第3、4个式子定义一个辅助线性系统

\left\{\begin{array}{l}{}_{\mathrm{0}}{}^c{D}_t^{\alpha }{E}_i={\beta }_i{S}_i^{\mathrm{0}}{I}_i-\left({\sigma }_i+{\mu }_i\right)E_i+{\displaystyle \sum _{j\ne i}^r}\left(b_{ij}{E}_j-b_{ji}{E}_i\right)\\ {}_{\mathrm{0}}{}^c{D}_t^{\alpha }{I}_i={\sigma }_i{E}_i-({\gamma }_i+{\epsilon }_i+{\mu }_i)I_i+{\displaystyle \sum _{j\ne i}^r}\left(c_{ij}{I}_j-c_{ji}{I}_i\right)\\ {}_{\mathrm{0}}{}^c{D}_t^{\alpha }{Q}_i={\epsilon }_i{I}_i-\left({\mu }_i+h_i\right)Q_i\end{array}\right.（14）

可见系统（14）的右侧是\mathit{F}-\mathit{V}的系数矩阵。当时，\mathit{F}-\mathit{V}的每个特征值都具有负实部。根据文献［17］中的定理1得到系统（14）的每个正解满足\underset{t\to \mathrm{\infty }}{\mathrm{l}\mathrm{i}\mathrm{m}}{E}_i\left(t\right)=\mathrm{0}，\underset{t\to \mathrm{\infty }}{\mathrm{l}\mathrm{i}\mathrm{m}}{I}_i\left(t\right)=\mathrm{0}，\underset{t\to \mathrm{\infty }}{\mathrm{l}\mathrm{i}\mathrm{m}}{Q}_i\left(t\right)=\mathrm{0}，因为系统（14）是一个线性系统，所以系统（14）的无病平衡点是全局渐进稳定的。又因为系统（2）的所有变量都是非负，则根据分数阶系统（FDE）的比较理论^［19］得

\underset{t\to \mathrm{\infty }}{\mathrm{l}\mathrm{i}\mathrm{m}}{E}_i\left(t\right)=\mathrm{0}，\underset{t\to \mathrm{\infty }}{\mathrm{l}\mathrm{i}\mathrm{m}}{I}_i\left(t\right)=\mathrm{0}，\underset{t\to \mathrm{\infty }}{\mathrm{l}\mathrm{i}\mathrm{m}}{Q}_i\left(t\right)=\mathrm{0}。

根据系统（2）的第1个式子，当t\to \mathrm{\infty }时，E_i,I_i,Q_i\to \mathrm{0}，有

{}_{\mathrm{0}}{}^c{}_t^{\alpha }{\mathit{S}}_i=\mathit{\Lambda }-\mathit{J}\mathit{S}=\mathit{J}{\mathit{S}}^{\mathrm{0}}-\mathit{J}\mathit{S}（15）

运用类似于上述定理2的证明方法，式（15）的解可以推导为

\begin{array}{l}\mathit{S}(t)=t^{\alpha }{E}_{\alpha ,\alpha +\mathrm{1}}(-\mathit{J}{t}^{\alpha })\mathit{J}{\mathit{S}}^{\mathrm{0}}+E_{\alpha ,\mathrm{1}}(-\mathit{J}{t}^{\alpha })\mathit{S}(\mathrm{0})\\ =(\mathrm{1}-E_{\alpha ,\mathrm{1}}(-\mathit{J}{t}^{\alpha })){\mathit{S}}^{\mathrm{0}}+E_{\alpha ,\mathrm{1}}(-\mathit{J}{t}^{\alpha })\mathit{S}(\mathrm{0})\\ ={\mathit{S}}^{\mathrm{0}}+E_{\alpha ,\mathrm{1}}(-\mathit{J}{t}^{\alpha })(\mathit{S}(\mathrm{0})-{\mathit{S}}^{\mathrm{0}})\mathrm{。}\end{array}

根据文献［17］得知Mittag-Leffler具有渐近行为，那么。

由-\mathit{J}具有所有实部为负的特征值可得，当t\to \mathrm{\infty }时，\mathit{S}\left(t\right)\to {\mathit{S}}^{\mathrm{0}}，得证无病平衡点E_{\mathrm{0}}是全局渐近稳定的。根据定理3可知，当R_{\mathrm{0}}>\mathrm{1}时，无病平衡点E_{\mathrm{0}}是不稳定的。

6 数值模拟

利用Matlab进行数值模拟，验证了无病平衡点E_{\mathrm{0}}的稳定性。以两个斑块为例，选择初值S_{\mathrm{1}}\left(\mathrm{0}\right)=\mathrm{1000}，E_{\mathrm{1}}\left(\mathrm{0}\right)=\mathrm{120}，I_{\mathrm{1}}\left(\mathrm{0}\right)=\mathrm{100}，Q_{\mathrm{1}}\left(\mathrm{0}\right)=\mathrm{80}，R_{\mathrm{1}}\left(\mathrm{0}\right)=\mathrm{250}，S_{\mathrm{2}}\left(\mathrm{0}\right)=\mathrm{800}，E_{\mathrm{2}}\left(\mathrm{0}\right)=\mathrm{120}，I_{\mathrm{2}}\left(\mathrm{0}\right)=\mathrm{80}，Q_{\mathrm{2}}\left(\mathrm{0}\right)=\mathrm{60}，R_{\mathrm{2}}\left(\mathrm{0}\right)=\mathrm{200}，假设a_{ij}=b_{ij}=c_{ij}=d_{ij}=m_{ij}，其中i,j=\mathrm{1,2}且i\ne j。

6.1 迁移运动对传染病的影响

图2取参数{\Lambda }_{\mathrm{1}}=\mathrm{100}，{\Lambda }_{\mathrm{2}}=\mathrm{300}，{\beta }_{\mathrm{1}}=\mathrm{0.0004}，{\beta }_{\mathrm{2}}=\mathrm{0.0008}，p_{\mathrm{1}}=\mathrm{0.6}，p_{\mathrm{2}}=\mathrm{0.7}，e_{\mathrm{1}}=\mathrm{0.8}，e_{\mathrm{2}}=\mathrm{0.9}，{\mu }_{\mathrm{1}}=\mathrm{0.02}，{\mu }_{\mathrm{2}}=\mathrm{0.01}，{\sigma }_{\mathrm{1}}=\mathrm{0.14}，{\sigma }_{\mathrm{2}}=\mathrm{0.18}，{\gamma }_{\mathrm{1}}=\mathrm{0.15}，{\gamma }_{\mathrm{2}}=\mathrm{0.18}，{\epsilon }_{\mathrm{1}}=\mathrm{0.4}，{\epsilon }_{\mathrm{2}}=\mathrm{0.3}，h_{\mathrm{1}}=\mathrm{0.1}，h_{\mathrm{2}}=\mathrm{0.15}，m_{\mathrm{12}}=\mathrm{0.04}，m_{\mathrm{21}}=\mathrm{0.02}。此时R_{\mathrm{0}}^{\mathrm{1}}=\mathrm{0.6140}，R_{\mathrm{0}}^{\mathrm{2}}=\mathrm{1.2374}且，其中R_{\mathrm{0}}^i表示没有迁移运动时的第i个斑块的基本再生数。

为了分析处于平衡状态的各仓室人群的动态，设置\alpha值分别为\alpha =\mathrm{0.7}，\alpha =\mathrm{0.8}，\alpha =\mathrm{0.9}和\alpha =\mathrm{1.0}（整数阶）。

图2（a）（c）（e）（g）和（i）以及图2（b）（d）（f）（h）和（j）中的每一条轨迹分别反映了斑块1、斑块2在迁移运动的情形下，在一组特定的参数值中，当α取不同的值时，易感者、潜伏者、感染者、隔离者以及恢复者达到稳定状态的变化趋势。可以发现随着\alpha的降低，感染者数量随时间的变化趋势变得更加平缓，这表明系统的响应变得更加缓慢，传染病传播和康复过程变得更加缓和，系统达到稳定状态的速度减慢，时间延长。因此，分数阶参数\alpha在描述各仓室的疾病动力学中起着重要作用，同时可以得到无病平衡点E_{\mathrm{0}}是全局渐进稳定的。此外，如果不考虑迁移运动，R_{\mathrm{0}}^{\mathrm{2}}=\mathrm{1.2374}>\mathrm{1}，那么疾病将会在斑块2中传播。但是，在考虑迁移运动后，，疾病会在斑块2中灭绝。

6.2 接种疫苗对传染病的影响

假设p_{\mathrm{1}}=p_{\mathrm{2}}，其他参数保持不变。在p_{\mathrm{1}}=p_{\mathrm{2}}=\mathrm{0.01}时，基本再生数R_{\mathrm{0}}=\mathrm{11.1563}，斑块i中的染病者数量随时间变化如图3（a）和（d）所示；在p_{\mathrm{1}}=p_{\mathrm{2}}=\mathrm{0.30}时，基本再生数R_{\mathrm{0}}=\mathrm{1.2140}，斑块i中的染病者数量随时间变化如图3（b）和（e）所示；在p_{\mathrm{1}}=p_{\mathrm{2}}=\mathrm{0.8}0时，基本再生数R_{\mathrm{0}}=\mathrm{0.4870}，斑块i中的染病者数量随时间变化如图3（c）和（f）所示。可见，增加疫苗接种率使得基本再生数R_{\mathrm{0}}减小，从而有效控制疾病的传播。

从图3可以观察到，随着疫苗接种率的增加，感染人口的数量减少，表明了疫苗接种策略在预防疾病传播中起着重要作用。此外，根据分数阶参数\alpha，每个图都有特定的轨迹和收敛速度。如果\alpha值较小，则收敛速度较慢，反之亦然。

7 结论

本文根据麻疹病毒的传播机制，建立了一类带有疫苗接种的分数阶\mathrm{S}\mathrm{E}\mathrm{I}\mathrm{Q}\mathrm{R}传染病斑块模型。利用下一代矩阵法计算了模型的基本再生数R_{\mathrm{0}}，根据Mittag-Leffler函数和Laplace变换得到了无病平衡点E_{\mathrm{0}}在时是局部和全局渐近稳定的，在R_{\mathrm{0}}>\mathrm{1}时是不稳定的。使用Matlab进行数值模拟，观察到模型解的收敛速度随分数阶次的降低而减慢，表明流行病患者全部康复所需时间随分数阶次的降低而变长，即分数阶传染病模型中的分数阶次\alpha对系统的稳定性条件具有重要影响，且随着\alpha的降低，系统的稳定性条件变得更加宽松。斑块之间进行迁移运动可能使得传染病在两个斑块中灭绝，同时提高疫苗的接种率，可以减少基本再生数的大小，从而控制疾病的传播。

参考文献

[1] GRIFFIN D E.The immune response in measles：Virus control，clearance and protective immunity［J］.Viruses，2016，8（10）：282.

[2] ALDILA D， ASRIANTI D.A deterministic model of measles with imperfect vaccination and quarantine intervention［J］.Journal of Physics：Conference Series，2019，1218（1）：012044.

[3] BERHE H W， MAKINDE O D.Computational modelling and optimal control of measles epidemic in human population［J］.Biosystems，2020，190：104102.

[4] MEMON Z， QURESHI S， MEMON B R.Mathematical analysis for a new nonlinear measles epidemiological system using real incidence data from Pakistan［J］.European Physical Journal Plus，2020，135（4）：378.

[5] FELLAH Z E A， FELLAH M， OGAM E，et al.Reflection and transmission of transient ultrasonic wave in fractal porous material：Application of fractional calculus［J］.Wave Motion，2021，106：102804.

[6] YU M， YU K S， HAN T Z，et al.Research on application of fractional calculus in signal analysis and processing of stock market［J］.Chaos Solitons & Fractals，2020，131：109468.

[7] 杨贵宝，高霞.基于分数阶logistic映射与随机变换的双图像加密算法［J］.内蒙古大学学报（自然科学版），2017，48（2）：189-195.

[8] ANNABY M H， AYAD H A， RUSHDI M A，et al.Difference operators and generalized discrete fractional transforms in signal and image processing［J］.Signal Processing，2018，151：1-18.

[9] KHAN A A， AMIN R， ULLAH S，et al.Numerical simulation of a Caputo fractional epidemic model for the novel coronavirus with the impact of environmental transmission［J］.Alexandria Engineering Journal，2022，61（7）：5083-5095.

[10] NISAR K S， AHMAD S， ULLAH A，et al.Mathematical analysis of SIRD model of COVID-19 with Caputo fractional derivative based on real data［J］.Results in Physics，2021，21：103772.

[11] DRIESSCHE P V D， SALMANI M.A model for disease transmission in a patchy environment［J］.Discrete and Continuous Dynamical Systems-B，2017，6（1）：185-202.

[12] CUI Q Q， QIU Z P， DING L.An SIR epidemic model with vaccination in a patchy environment［J］.Mathematical Biosciences and Engineering，2017，14（5/6）：1141-1157.

[13] BERMAN A.Nonnegative matrices in mathematical sciences［M］.New York：Academic Press，1979.

[14] LIN W.Global existence theory and chaos control of fractional differential equations［J］.Journal of Mathematical Analysis and Applications，2007，332（1）：709-726.

[15] ODIBAT Z M， SHAWAGFEH N T.Generalized Taylor's formula［J］.Applied Mathematics and Computation，2007，186（1）：286-293.

[16] SALEEM S， RAFIQ M， AHMED N，et al.Fractional epidemic model of coronavirus disease with vaccination and crowding effects［J］.Scientific Reports，2024，14（1）：8157.

[17] KHEIRI H， JAFARI M.Stability analysis of a fractional order model for the HIV/AIDS epidemic in a patchy environment［J］.Journal of Computational and Applied Mathematics，2019，346：323-339.

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

[19] WANG Z L， YANG D S， MA T D，et al.Stability analysis for nonlinear fractional-order systems based on comparison principle［J］.Nonlinear Dynamics，2014，75（1）：387-402.

基金资助

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

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