无标度网络上带有隔离项的时滞SIQRS模型稳定性分析

Stability Analysis of Delayed SIQRS Model with Isolation Item on Scale-Free Networks

摘要

考虑了无标度网络上一类带有隔离项的时滞SIQRS传染病模型的动力学行为。首先，利用地方病平衡态的存在性得到了基本再生数R_{\mathrm{0}}的表达式；接着，通过构造适当的李雅普诺夫函数证明了平衡态的全局稳定性，结果表明：当R_{\mathrm{0}}\le \mathrm{1}时，无病平衡态全局渐近稳定，隔离措施对控制疾病传播没有影响；否则地方病平衡态全局渐近稳定，隔离措施在一定程度上降低了最终感染规模。最后，利用MATLAB软件对分析结果进行了数值模拟，发现当时，时滞减缓了疾病的传播；当R_{\mathrm{0}}>\mathrm{1}时，时滞促进了疾病的传播。

Abstract

A scale-free network based SIQRS epidemic model with isolation item and time delay is considered in the present paper. At first, the formula of the basic reproduction number is given by utilizing the existence of the endemic equilibrium. Next, the global stability of the equilibrium is discussed by constructing appropriate Lyapunov functions, and the obtained results manifest that the disease-free equilibrium is globally asymptotically stable when R_{\mathrm{0}}\le \mathrm{1}, and isolate measures have no effect on controlling the spread of the disease. Otherwise the endemic equilibrium is globally asymptotically stable, and isolate measures reduce the final infection scale to a certain extent. Furthermore, the results are simulated by MATLAB. It is found that time delay slows down the spread of disease when and time delay promotes the spread of disease when R_{\mathrm{0}}>\mathrm{1}.

Graphical abstract

关键词

时滞 / 隔离措施 / 无标度网络 / 全局稳定性

Key words

time delay / isolate measures / scale-free network / global stability

0　 引 言

随着复杂网络这门新兴学科的迅速兴起及无标度网络的提出^［1］，人们发现生活中大量的复杂系统均表现出无标度性（度分布为），如万维网、交通网、接触网等^［2，3］。于是利用复杂网络研究传染病的动力学行为引起了许多学者的广泛关注^［4~7］，其中最为经典的是文献［6］，它在研究无标度网络上的SIS模型时发现：当网络规模趋于无穷大时，传播阈值趋于零，这很好地解释了计算机病毒长期存在而不能绝迹的原因。隔离措施在控制传染病的传播中发挥着举足轻重的作用，为探究隔离措施对控制疾病传播的影响，文献［7］提出了复杂网络上带有隔离项及疫苗接种的SIQRS传染病模型

\left\{\begin{array}{l}\frac{\mathrm{d}{S}_k\left(t\right)}{dt}=-\beta k{S}_k\left(t\right)\mathrm{\Theta }\left(t\right)-\alpha S_k\left(t\right)+\delta R_k\left(t\right)\\ \frac{\mathrm{d}{I}_k\left(t\right)}{\mathrm{d}t}=\beta k{S}_k\left(t\right)\Theta \left(t\right)-\sigma I_k\left(t\right)-\gamma I_k\left(t\right)\\ \frac{\mathrm{d}{Q}_k\left(t\right)}{dt}=\sigma I_k\left(t\right)-\epsilon Q_k\left(t\right)\\ \frac{\mathrm{d}{R}_k\left(t\right)}{\mathrm{d}t}=\alpha S_k\left(t\right)+\gamma I_k\left(t\right)+\epsilon Q_k\left(t\right)-\delta R_k\left(t\right)\end{array}\right.（1）

其中，S_k\left(t\right),I_k\left(t\right),Q_k\left(t\right),R_k\left(t\right)分别表示t时刻度为k\left(=\mathrm{1,2},\cdots ,n\right)的易感者、感染者、隔离者、恢复者的密度，\Theta \left(t\right)=\frac{\mathrm{1}}{〈k〉}{\displaystyle \sum _{k=\mathrm{1}}^n}kP\left(k\right)I_k\left(t\right)表示任意一条给定的边与一个染病者相连的概率，k\left(=\mathrm{1,2},\cdots ,n\right)表示节点的度（即一个人在单位时间内接触的人数），\beta表示传播率系数，\alpha是易感者的疫苗接种率，\delta是恢复者的失去免疫率，\sigma是感染者的隔离率，\gamma是感染者的治愈率，\epsilon是隔离者的治愈率。

事实上，许多传染病在传播过程均表现出时滞性，主要有潜伏期时滞、免疫期时滞、感染期时滞^［8，9］等。文献［10］分析了时滞、非线性发病率等因素对无标度网络上的SIS模型动力学行为的影响。文献［11］研究了无标度网络上具有时滞的SEIRS传染病模型，表明时滞对基本再生数无影响。文献［12］通过选取时滞作为分支参数，得到了具有非线性发病率的时滞SEIS模型Hopf分支存在的条件。文献［13］研究了预警意识、延迟提供预警资金对经典SEIR模型动力学行为的影响。这些研究在建立模型时均假设时滞为单值，然而许多病原体的感染期和潜伏期等会随着环境的变化而发生变化，不同的感染者由于自身健康状况的差异性，恢复期也会有所不同，为更好地揭示这种时滞现象，文献［14］通过将t时刻的发病率假设为\beta k{S}_k\left(t\right){\int }_{\mathrm{0}}^{\mathrm{\infty }}f\left(\tau \right){\mathrm{e}}^{-\left(\mu +\lambda \right)\tau }\Theta \left(t-\tau \right)\mathrm{d}\tau，提出了复杂网络上具有分布时滞的SIR模型，其中{\mathrm{e}}^{-\left(\mu +\lambda \right)\tau }表示感染个体在潜伏期的存活概率，可积函数f\left(\tau \right)用来描述感染性随时间变化的情况，满足f\left(\tau \right)\ge \mathrm{0},\tau \in \left[\mathrm{0},\mathrm{\infty }\right)且{\int }_{\mathrm{0}}^{\mathrm{\infty }}f\left(\tau \right)\mathrm{d}\tau =\mathrm{1}。

由于复杂网络上的传染病模型维数较高，分析其动力学行为较为困难，所以文献［7］未讨论地方病平衡态的稳定性；文献［14］所构建的SIR模型未考虑疫苗注射、隔离措施等因素的影响，所以模型耦合度较低，故作者将模型解耦后利用基尔霍夫矩阵树定理分析了模型的全局动力学行为。本文综合考虑了文献［7］中的隔离措施、疫苗注射及文献［14］中的时滞因素，提出了无标度网络上带有隔离项的时滞SIQRS模型

\left\{\begin{array}{l}\frac{\mathrm{d}{S}_k\left(t\right)}{\mathrm{d}t}=-\beta k{S}_k\left(t\right){\int }_{\mathrm{0}}^{\mathrm{\infty }}f\left(\tau \right){\mathrm{e}}^{-\left(\mu +\gamma +\sigma \right)\tau }\mathrm{\Theta }\left(t-\tau \right)\mathrm{d}\tau -\alpha S_k\left(t\right)-\mu S_k\left(t\right)+\delta R_k\left(t\right)+\mu \\ \frac{\mathrm{d}{I}_k\left(t\right)}{\mathrm{d}t}=\beta k{S}_k\left(t\right){\int }_{\mathrm{0}}^{\mathrm{\infty }}f\left(\tau \right){\mathrm{e}}^{-\left(\mu +\gamma +\sigma \right)\tau }\mathrm{\Theta }\left(t-\tau \right)\mathrm{d}\tau -\sigma I_k\left(t\right)-\gamma I_k\left(t\right)-\mu I_k\left(t\right)\\ \frac{\mathrm{d}{Q}_k\left(t\right)}{\mathrm{d}t}=\sigma I_k\left(t\right)-\epsilon Q_k\left(t\right)-\mu Q_k\left(t\right)\\ \frac{\mathrm{d}{R}_k\left(t\right)}{\mathrm{d}t}=\alpha S_k\left(t\right)+\gamma I_k\left(t\right)+\epsilon Q_k\left(t\right)-\mu R_k\left(t\right)-\delta R_k\left(t\right)\end{array}\right.（2）

其中，\beta k{S}_k\left(t\right){\int }_{\mathrm{0}}^{\mathrm{\infty }}f\left(\tau \right){\mathrm{e}}^{-\left(\mu +\gamma +\sigma \right)\tau }\Theta \left(t-\tau \right)\mathrm{d}\tau表示t时刻的发病率，满足{\int }_{\mathrm{0}}^{\mathrm{\infty }}f\left(\tau \right)\mathrm{d}\tau =\mathrm{1}的非负可积函数f\left(\tau \right)表示感染性随时间变化的情况，{\mathrm{e}}^{-\left(\mu +\gamma +\sigma \right)\tau }表示感染个体在潜伏期内的存活概率，\mu表示个体的常数输入及自然死亡率，其他符号含义如（1）式。

1　 平衡态的存在性

显然，模型（2）的初值满足

\Theta \left(\mathrm{0}\right)>\mathrm{0},S_k\left(\mathrm{0}\right)+I_k\left(\mathrm{0}\right)+Q_k\left(\mathrm{0}\right)+R_k\left(\mathrm{0}\right)=\mathrm{1}\left(k=\mathrm{1,2},\cdots ,n\right)（3）

定理1 令\left(S_{\mathrm{1}}\left(t\right),I_{\mathrm{1}}\left(t\right),Q_{\mathrm{1}}\left(t\right),R_{\mathrm{1}}\left(t\right)\cdots ,S_n\left(t\right),I_n\left(t\right),Q_n\left(t\right),R_n\left(t\right)\right)为模型（2）满足初值条件（3）的解，则当t>\mathrm{0}时，有且。

证 首先证明当t>\mathrm{0}时，有。因为I_k\left(t\right)连续，所以存在一个\xi >\mathrm{0}，使得当t\in \left(\mathrm{0},\xi \right)有。假设当t>\mathrm{0}时，I_k\left(t\right)=\mathrm{1}\left(k=\mathrm{1,2},\cdots ,n\right)，则存在一个t_{\mathrm{0}}\ge \xi >\mathrm{0}，使得I_k\left(t_{\mathrm{0}}\right)=\mathrm{1}，且当t\in \left(\mathrm{0},t_{\mathrm{0}}\right)有。

由模型（2）的第二个式子可得，故存在\tilde{t}\in \left(\mathrm{0},t_{\mathrm{0}}\right)，使得I_k\left(\tilde{t}\right)>\mathrm{1}，与假设矛盾。所以当t>\mathrm{0}时，有。

接着证明当t>\mathrm{0}时，有I_k\left(t\right)>\mathrm{0}\left(k=\mathrm{1,2},\cdots ,n\right)。对模型（2）的第二个式子两边从0到t积分可得

I_k\left(t\right)=I_k\left(\mathrm{0}\right){\mathrm{e}}^{-\left(\sigma +\gamma +\mu \right)t}+{\int }_{\mathrm{0}}^t{\mathrm{e}}^{-\left(\sigma +\gamma +\mu \right)\left(t-u\right)}\beta k{S}_k\left(u\right){\int }_{\mathrm{0}}^{\mathrm{\infty }}f\left(\tau \right){\mathrm{e}}^{-\left(\sigma +\gamma +\mu \right)\tau }\Theta \left(u-\tau \right)\mathrm{d}\tau \mathrm{d}u（4）

假设当t>\mathrm{0}时，有I_k\left(t\right)=\mathrm{0}\left(k=\mathrm{1,2},\cdots ,n\right)，则存在一个t_{\mathrm{1}}>\mathrm{0}，使得I_k\left(t_{\mathrm{1}}\right)=\mathrm{0}。而由（4）式可知

I_k\left(t_{\mathrm{1}}\right)=I_k\left(\mathrm{0}\right){\mathrm{e}}^{-\left(\sigma +\gamma +\mu \right)t_{\mathrm{1}}}+{\int }_{\mathrm{0}}^{t_{\mathrm{1}}}{\mathrm{e}}^{-\left(\sigma +\gamma +\mu \right)\left(t_{\mathrm{1}}-u\right)}\beta k{S}_k\left(u\right){\int }_{\mathrm{0}}^{\mathrm{\infty }}f\left(\tau \right){\mathrm{e}}^{-\left(\sigma +\gamma +\mu \right)\tau }\mathrm{\Theta }\left(u-\tau \right)\mathrm{d}\tau \mathrm{d}u>\mathrm{0}

与假设矛盾，所以当t>\mathrm{0}时，有I_k\left(t\right)>\mathrm{0}\left(k=\mathrm{1,2},\cdots ,n\right)。

综上可得，当t>\mathrm{0}时，有。

显然，当t>\mathrm{0}时，有。同理可证，当t>\mathrm{0}时，有\left(k=\mathrm{1,2},\cdots ,n\right)。

定理2 定义R_{\mathrm{0}}=\frac{〈k^{\mathrm{2}}〉}{〈k〉}\frac{\left(\mu +\delta \right)\beta \tilde{f}}{\left(\sigma +\gamma +\mu \right)\left(\mu +\alpha +\delta \right)}，则下列命题成立

1） 模型（2）总存在无病平衡态E_k^{\mathrm{0}}={\left\{\frac{\mu +\delta }{\mu +\alpha +\delta },\mathrm{0,0},\frac{\alpha }{\mu +\alpha +\delta }\right\}}_{k=\mathrm{1}}^n。

2） 当R_{\mathrm{0}}>\mathrm{1}时，模型（2）存在地方病平衡态E_k^{\mathrm{*}}={\left\{S_k^{\mathrm{*}},I_k^{\mathrm{*}},Q_k^{\mathrm{*}},R_k^{\mathrm{*}}\right\}}_{k=\mathrm{1}}^n，其中

S_k^{\mathrm{*}}=\frac{\sigma +\gamma +\mu }{\beta k\tilde{f}{\mathrm{\Theta }}^{\mathrm{*}}}{I}_k^{\mathrm{*}},Q_k^{\mathrm{*}}=\frac{\sigma }{\mu +\epsilon }{I}_k^{\mathrm{*}},R_k^{\mathrm{*}}=\frac{\left(\sigma +\gamma +\mu \right)\left(\beta k\tilde{f}{\mathrm{\Theta }}^{\mathrm{*}}+\alpha +\mu \right)}{\delta \beta k\tilde{f}{\mathrm{\Theta }}^{\mathrm{*}}}{I}_k^{\mathrm{*}}-\frac{\mu }{\delta }

I_k^{\mathrm{*}}=\frac{\left(\mu +\delta \right)\left(\mu +\epsilon \right)\beta k\tilde{f}{\Theta }^{\mathrm{*}}}{\left[\sigma \delta +\left(\mu +\epsilon \right)\left(\delta +\sigma +\gamma +\mu \right)\right]\beta k\tilde{f}{\Theta }^{\mathrm{*}}+\left(\mu +\epsilon \right)\left(\sigma +\gamma +\mu \right)\left(\delta +\alpha +\mu \right)}

证 显然，模型（2）总存在无病平衡态E_k^{\mathrm{0}}。令\tilde{f}={\int }_{\mathrm{0}}^{\mathrm{\infty }}f\left(\tau \right){\mathrm{e}}^{-\left(\mu +\gamma +\sigma \right)\tau }\mathrm{d}\tau，在平衡态E_k^{\mathrm{*}}有

\left\{\begin{array}{l}-\beta k{S}_k^{\mathrm{*}}\tilde{f}{\mathrm{\Theta }}^{\mathrm{*}}-\alpha S_k^{\mathrm{*}}-\mu S_k^{\mathrm{*}}+\delta R_k^{\mathrm{*}}+\mu =\mathrm{0}\\ \beta k{S}_k^{\mathrm{*}}\tilde{f}{\mathrm{\Theta }}^{\mathrm{*}}-\sigma I_k^{\mathrm{*}}-\gamma I_k^{\mathrm{*}}-\mu I_k^{\mathrm{*}}=\mathrm{0}\\ \sigma I_k^{\mathrm{*}}-\epsilon Q_k^{\mathrm{*}}-\mu Q_k^{\mathrm{*}}=\mathrm{0}\\ \alpha S_k^{\mathrm{*}}+\gamma I_k^{\mathrm{*}}+\epsilon Q_k^{\mathrm{*}}-\mu R_k^{\mathrm{*}}-\delta R_k^{\mathrm{*}}=\mathrm{0}\end{array}\right.

其中{\Theta }^{\mathrm{*}}=\frac{\mathrm{1}}{〈k〉}{\displaystyle \sum _{k=\mathrm{1}}^n}kP\left(k\right)I_k^{\mathrm{*}}，通过计算可得

\left\{\begin{array}{l}{S}_k^{\mathrm{*}}=\frac{\sigma +\gamma +\mu }{\beta k\tilde{f}{\Theta }^{\mathrm{*}}}{I}_k^{\mathrm{*}}\\ Q_k^{\mathrm{*}}=\frac{\sigma }{\mu +\epsilon }{I}_k^{\mathrm{*}}\\ R_k^{\mathrm{*}}=\frac{\left(\sigma +\gamma +\mu \right)\left(\beta k\tilde{f}{\Theta }^{\mathrm{*}}+\alpha +\mu \right)}{\delta \beta k\tilde{f}{\Theta }^{\mathrm{*}}}{I}_k^{\mathrm{*}}-\frac{\mu }{\delta }\end{array}\right.

使用归一化条件，有

I_k^{\mathrm{*}}=\frac{\left(\mu +\delta \right)\left(\mu +\epsilon \right)\beta k\tilde{f}{\Theta }^{\mathrm{*}}}{\left[\sigma \delta +\left(\mu +\epsilon \right)\left(\delta +\sigma +\gamma +\mu \right)\right]\beta k\tilde{f}{\Theta }^{\mathrm{*}}+\left(\mu +\epsilon \right)\left(\sigma +\gamma +\mu \right)\left(\delta +\alpha +\mu \right)}（5）

将（5）式代入\Theta得

{\Theta }^{\mathrm{*}}=\frac{{\displaystyle \sum _{k=\mathrm{1}}^n}kP\left(k\right)}{〈k〉}\frac{\left(\mu +\delta \right)\left(\mu +\epsilon \right)\beta k\tilde{f}{\Theta }^{\mathrm{*}}}{\left[\sigma \delta +\left(\mu +\epsilon \right)\left(\delta +\sigma +\gamma +\mu \right)\right]\beta k\tilde{f}{\Theta }^{\mathrm{*}}+\left(\mu +\epsilon \right)\left(\sigma +\gamma +\mu \right)\left(\delta +\alpha +\mu \right)}≕f\left({\Theta }^{\mathrm{*}}\right)（6）

显然，{\Theta }^{\mathrm{*}}=\mathrm{0}是方程（6）的一个平凡解，其对应着模型（2）的无病平衡态。为使方程（6）存在非平凡解，则需满足{\left.\frac{\mathrm{d}f\left({\Theta }^{\mathrm{*}}\right)}{\mathrm{d}{\Theta }^{\mathrm{*}}}\right|}_{{\Theta }^{\mathrm{*}}=\mathrm{0}}>\mathrm{1}且f\left(\mathrm{1}\right)\le \mathrm{1}。即

\left\{\begin{array}{l}\frac{〈k^{\mathrm{2}}〉}{〈k〉}\frac{\left(\delta +\mu \right)\beta \tilde{f}}{\left(\sigma +\gamma +\mu \right)\left(\delta +\alpha +\mu \right)}>\mathrm{1}\\ \frac{〈k^{\mathrm{2}}〉}{〈k〉}\frac{\left(\mu +\delta \right)\left(\mu +\epsilon \right)\beta \tilde{f}}{\left[\sigma \delta +\left(\mu +\epsilon \right)\left(\delta +\sigma +\gamma +\mu \right)\right]\beta k\tilde{f}+\left(\mu +\epsilon \right)\left(\sigma +\gamma +\mu \right)\left(\delta +\alpha +\mu \right)}\le \mathrm{1}\end{array}\right.（7）

要使（7）式成立，则不等式\frac{〈k^{\mathrm{2}}〉}{〈k〉}\frac{\left(\delta +\mu \right)\beta \tilde{f}}{\left(\sigma +\gamma +\mu \right)\left(\delta +\alpha +\mu \right)}>\mathrm{1}需成立。故定义R_{\mathrm{0}}=\frac{〈k^{\mathrm{2}}〉}{〈k〉}\cdot\frac{\left(\delta +\mu \right)\beta \tilde{f}}{\left(\sigma +\gamma +\mu \right)\left(\delta +\alpha +\mu \right)}，当且仅当R_{\mathrm{0}}>\mathrm{1}时，（6）式存在非平凡解，即模型（2）的地方病平衡态存在。

2　 平衡态的全局稳定性

2.1　无病平衡态的稳定性

定理3 当，模型（2）的无病平衡态E_k^{\mathrm{0}}局部渐近稳定；当R_{\mathrm{0}}>\mathrm{1}，无病平衡态E_k^{\mathrm{0}}不稳定。

证 不失一般性，令{\tilde{S}}_k=S_k-S_k^{\mathrm{0}},{\tilde{I}}_k=I_k-I_k^{\mathrm{0}},{\tilde{Q}}_k=Q_k-Q_k^{\mathrm{0}},{\tilde{R}}_k=R_k-R_k^{\mathrm{0}}。则模型（2）可被线性化为

\left\{\begin{array}{l}\frac{\mathrm{d}{S}_k\left(t\right)}{\mathrm{d}t}=-\beta k{S}_k^{\mathrm{0}}{\int }_{\mathrm{0}}^{\mathrm{\infty }}f\left(\tau \right){\mathrm{e}}^{-\left(\mu +\gamma +\sigma \right)\tau }\Theta \left(t-\tau \right)\mathrm{d}\tau -\alpha S_k\left(t\right)-\mu S_k\left(t\right)+\delta R_k\left(t\right)\\ \frac{\mathrm{d}{I}_k\left(t\right)}{\mathrm{d}t}=\beta k{S}_k^{\mathrm{0}}{\int }_{\mathrm{0}}^{\mathrm{\infty }}f\left(\tau \right){\mathrm{e}}^{-\left(\mu +\gamma +\sigma \right)\tau }\Theta \left(t-\tau \right)\mathrm{d}\tau -\sigma I_k\left(t\right)-\gamma I_k\left(t\right)-\mu I_k\left(t\right)\\ \frac{\mathrm{d}{Q}_k\left(t\right)}{\mathrm{d}t}=\sigma I_k\left(t\right)-\epsilon Q_k\left(t\right)-\mu Q_k\left(t\right)\\ \frac{\mathrm{d}{R}_k\left(t\right)}{\mathrm{d}t}=\alpha S_k\left(t\right)+\gamma I_k\left(t\right)+\epsilon Q_k\left(t\right)-\mu R_k\left(t\right)-\delta R_k\left(t\right)\end{array}\right.（8）

为书写方便，这里将“~”省略。接着探究线性模型（8）的稳定性，不妨假设其有如下形式的指数解

S_k\left(t\right)=S_{k\mathrm{0}}{\mathrm{e}}^{\rho t},I_k\left(t\right)=I_{k\mathrm{0}}{\mathrm{e}}^{\rho t},Q_k\left(t\right)=Q_{k\mathrm{0}}{\mathrm{e}}^{\rho t},R_k\left(t\right)=R_{k\mathrm{0}}{\mathrm{e}}^{\rho t},k=\mathrm{1,2},\cdots ,n（9）

其中S_{k\mathrm{0}},I_{k\mathrm{0}},Q_{k\mathrm{0}},R_{k\mathrm{0}}\left(k=\mathrm{1,2},\cdots ,n\right)表示模型（8）的初值，将（9）式代入模型（8）得

\rho {\left(S_{\mathrm{10}},\cdots ,S_{n\mathrm{0}},I_{\mathrm{10}},\cdots ,I_{n\mathrm{0}},Q_{\mathrm{10}},\cdots ,Q_{n\mathrm{0}},R_{\mathrm{10}},\cdots ,R_{n\mathrm{0}}\right)}^{\mathrm{T}}=\mathit{A}{\left(S_{\mathrm{10}},\cdots ,S_{n\mathrm{0}},I_{\mathrm{10}},\cdots ,I_{n\mathrm{0}},Q_{\mathrm{10}},\cdots ,Q_{n\mathrm{0}},R_{\mathrm{10}},\cdots ,R_{n\mathrm{0}}\right)}^{\mathrm{T}}

其中，\rho为模型（8）的特征值。令

\Phi ={\int }_{\mathrm{0}}^{\mathrm{\infty }}f\left(\tau \right){\mathrm{e}}^{-\left(\rho +\mu +\gamma +\sigma \right)\tau }\mathrm{d}\tau

a_k=\frac{\mathrm{1}}{〈k〉}kP\left(k\right)

\mathit{A}={\left(\begin{array}{cc}{\mathit{A}}_{\mathrm{1}}& {\mathit{A}}_{\mathrm{2}}\\ {\mathit{A}}_{\mathrm{3}}& {\mathit{A}}_{\mathrm{4}}\end{array}\right)}_{{}_{\mathrm{4}n\times \mathrm{4}n}}

{\mathit{A}}_{\mathrm{1}}={\left(\begin{array}{cccccc}-\left(\alpha +\mu \right)& \cdots & \mathrm{0}& -\beta S_{\mathrm{1}}^{\mathrm{0}}{a}_{\mathrm{1}}\Phi & \cdots & -\beta S_{\mathrm{1}}^{\mathrm{0}}{a}_n\Phi \\ ⋮& \ddots & ⋮& ⋮& \ddots & ⋮\\ \mathrm{0}& \cdots & -\left(\alpha +\mu \right)& -\beta n{S}_n^{\mathrm{0}}{a}_{\mathrm{1}}\Phi & \cdots & -\beta n{S}_n^{\mathrm{0}}{a}_n\Phi \\ \mathrm{0}& \cdots & \mathrm{0}& -\left(\sigma +\gamma +\mu \right)+\beta S_{\mathrm{1}}^{\mathrm{0}}{a}_{\mathrm{1}}\Phi & \cdots & \beta S_{\mathrm{1}}^{\mathrm{0}}{a}_n\Phi \\ ⋮& \ddots & ⋮& ⋮& \ddots & ⋮\\ \mathrm{0}& \cdots & \mathrm{0}& \beta n{S}_n^{\mathrm{0}}{a}_{\mathrm{1}}\Phi & \cdots & -\left(\sigma +\gamma +\mu \right)+\beta n{S}_n^{\mathrm{0}}{a}_n\Phi \end{array}\right)}_{{}_{\mathrm{2}n\times \mathrm{2}n}},

{\mathit{A}}_{\mathrm{2}}={\left(\begin{array}{cccccc}\mathrm{0}& \cdots & \mathrm{0}& \delta & \cdots & \mathrm{0}\\ ⋮& \ddots & ⋮& ⋮& \ddots & ⋮\\ \mathrm{0}& \cdots & \mathrm{0}& \mathrm{0}& \cdots & \delta \\ \mathrm{0}& \cdots & \mathrm{0}& \mathrm{0}& \cdots & \mathrm{0}\\ ⋮& \ddots & ⋮& ⋮& \ddots & ⋮\\ \mathrm{0}& \cdots & \mathrm{0}& \mathrm{0}& \cdots & \mathrm{0}\end{array}\right)}_{{}_{\mathrm{2}n\times \mathrm{2}n}},{\mathit{A}}_{\mathrm{3}}={\left(\begin{array}{cccccc}\mathrm{0}& \cdots & \mathrm{0}& \sigma & \cdots & \mathrm{0}\\ ⋮& \ddots & ⋮& ⋮& \ddots & ⋮\\ \mathrm{0}& \cdots & \mathrm{0}& \mathrm{0}& \cdots & \sigma \\ \alpha & \cdots & \mathrm{0}& \gamma & \cdots & \mathrm{0}\\ ⋮& \ddots & ⋮& ⋮& \ddots & ⋮\\ \mathrm{0}& \cdots & \alpha & \mathrm{0}& \cdots & \gamma \end{array}\right)}_{{}_{\mathrm{2}n\times \mathrm{2}n}},

{\mathit{A}}_{\mathrm{4}}={\left(\begin{array}{cccccc}-\left(\epsilon +\mu \right)& \cdots & \mathrm{0}& \mathrm{0}& \cdots & \mathrm{0}\\ ⋮& \ddots & ⋮& ⋮& \ddots & ⋮\\ \mathrm{0}& \cdots & -\left(\epsilon +\mu \right)& \mathrm{0}& \cdots & \mathrm{0}\\ \epsilon & \cdots & \mathrm{0}& -\left(\mu +\delta \right)& \cdots & \mathrm{0}\\ ⋮& \ddots & ⋮& ⋮& \ddots & ⋮\\ \mathrm{0}& \cdots & \epsilon & \mathrm{0}& \cdots & -\left(\mu +\delta \right)\end{array}\right)}_{\mathrm{2}n\times \mathrm{2}n}

通过计算得：矩阵\mathit{A}有n重特征值-\mu，n重特征值-\left(\mu +\delta +\alpha \right)及n重特征值-\left(\mu +\epsilon \right)，为得到矩阵\mathit{A}的其余特征值，定义

\mathit{\Lambda }={\left(\begin{array}{cccc}-\left(\sigma +\gamma +\mu \right)+\beta S_{\mathrm{1}}^{\mathrm{0}}{a}_{\mathrm{1}}\Phi & \beta S_{\mathrm{1}}^{\mathrm{0}}{a}_{\mathrm{2}}\Phi & \cdots & \beta S_{\mathrm{1}}^{\mathrm{0}}{a}_n\Phi \\ \beta \mathrm{2}{S}_{\mathrm{2}}^{\mathrm{0}}{a}_{\mathrm{1}}\Phi & -\left(\sigma +\gamma +\mu \right)+\beta \mathrm{2}{S}_{\mathrm{2}}^{\mathrm{0}}{a}_{\mathrm{2}}\Phi & \cdots & \beta \mathrm{2}{S}_{\mathrm{2}}^{\mathrm{0}}{a}_n\Phi \\ \beta \mathrm{3}{S}_{\mathrm{3}}^{\mathrm{0}}{a}_{\mathrm{1}}\Phi & \beta \mathrm{3}{S}_{\mathrm{3}}^{\mathrm{0}}{a}_{\mathrm{2}}\Phi & \cdots & \beta \mathrm{3}{S}_{\mathrm{3}}^{\mathrm{0}}{a}_n\Phi \\ ⋮& ⋮& \ddots & ⋮\\ \beta n{S}_n^{\mathrm{0}}{a}_{\mathrm{1}}\Phi & \beta n{S}_n^{\mathrm{0}}{a}_{\mathrm{2}}\Phi & \cdots & -\left(\sigma +\gamma +\mu \right)+\beta n{S}_n^{\mathrm{0}}{a}_n\Phi \end{array}\right)}_{n\times n}

通过相似变换，可得矩阵\mathit{\Lambda }的相似矩阵

{\mathit{\Lambda }}^{\mathbf{*}}={\left(\begin{array}{cccc}-\left(\sigma +\gamma +\mu \right)& \mathrm{0}& \cdots & \beta \Phi S_{\mathrm{1}}^{\mathrm{0}}{a}_n\\ \mathrm{0}& -\left(\sigma +\gamma +\mu \right)& \cdots & \beta \Phi \left(S_{\mathrm{1}}^{\mathrm{0}}{a}_{\mathrm{1}}+\mathrm{2}{S}_{\mathrm{2}}^{\mathrm{0}}{a}_{\mathrm{2}}\right)\frac{a_n}{a_{\mathrm{2}}}\\ \mathrm{0}& \mathrm{0}& \cdots & \beta \Phi \left(S_{\mathrm{1}}^{\mathrm{0}}{a}_{\mathrm{1}}+\mathrm{2}{S}_{\mathrm{2}}^{\mathrm{0}}{a}_{\mathrm{2}}+\mathrm{3}{S}_{\mathrm{3}}^{\mathrm{0}}{a}_{\mathrm{3}}\right)\frac{a_n}{a_{\mathrm{3}}}\\ ⋮& ⋮& \ddots & ⋮\\ \mathrm{0}& \mathrm{0}& \cdots & -\left(\sigma +\gamma +\mu \right)+\beta \Phi {\displaystyle \sum _{k=\mathrm{1}}^n}k{S}_k^{\mathrm{0}}{a}_k\end{array}\right)}_{n\times n}

显然，矩阵{\mathit{\Lambda }}^{\mathrm{*}}有n-\mathrm{1}重特征值-\left(\sigma +\gamma +\mu \right)及特征值

\eta =-\left(\sigma +\gamma +\mu \right)+\beta \Phi {\displaystyle \sum _{k=\mathrm{1}}^n}k{S}_k^{\mathrm{0}}{a}_k=

\left(\sigma +\gamma +\mu \right)\left[\frac{\beta \Phi {\displaystyle \sum _{k=\mathrm{1}}^n}k\mu \left(\mu +\delta \right)kP\left(k\right)}{\left(\sigma +\gamma +\mu \right)\left(\delta +\alpha +\mu \right)〈k〉}-\mathrm{1}\right]=

\left(\sigma +\gamma +\mu \right)\left[\frac{〈k^{\mathrm{2}}〉\left(\delta +\mu \right)}{〈k〉\left(\sigma +\gamma +\mu \right)\left(\delta +\alpha +\mu \right)}\cdot {\int }_{\mathrm{0}}^{\mathrm{\infty }}f\left(\tau \right){\mathrm{e}}^{-\left(\eta +\mu +\gamma +\sigma \right)\tau }\mathrm{d}\tau -\mathrm{1}\right]（10）

假设\eta =a+b\mathrm{i}，则有

a+b\mathrm{i}=\left(\sigma +\gamma +\mu \right)R_{\mathrm{0}}\frac{{\int }_{\mathrm{0}}^{\mathrm{\infty }}f\left(\tau \right){\mathrm{e}}^{-\left(a+\mu +\gamma +\sigma \right)\tau }\left(\mathrm{c}\mathrm{o}\mathrm{s}b\tau -\mathrm{i}\mathrm{s}\mathrm{i}\mathrm{n}b\tau \right)\mathrm{d}\tau }{\tilde{f}}-\left(\sigma +\gamma +\mu \right)

分离实部和虚部得

\begin{array}{l}a=\left(\sigma +\gamma +\mu \right)\left[R_{\mathrm{0}}\frac{{\int }_{\mathrm{0}}^{\mathrm{\infty }}f\left(\tau \right){\mathrm{e}}^{-\left(a+\mu +\gamma +\sigma \right)\tau }\mathrm{c}\mathrm{o}\mathrm{s}b\tau \mathrm{d}\tau }{\tilde{f}}-\mathrm{1}\right]\\ b=-\left(\sigma +\gamma +\mu \right)R_{\mathrm{0}}\frac{{\int }_{\mathrm{0}}^{\mathrm{\infty }}f\left(\tau \right){\mathrm{e}}^{-\left(a+\mu +\gamma +\sigma \right)\tau }\mathrm{s}\mathrm{i}\mathrm{n}b\tau \mathrm{d}\tau }{\tilde{f}}\end{array}

令

V\left(x\right)=x-\left(\sigma +\gamma +\mu \right)\left[R_{\mathrm{0}}\frac{{\int }_{\mathrm{0}}^{\mathrm{\infty }}f\left(\tau \right){\mathrm{e}}^{-\left(x+\mu +\gamma +\sigma \right)\tau }\mathrm{d}\tau }{\tilde{f}}-\mathrm{1}\right]

则

\frac{\mathrm{d}V\left(x\right)}{\mathrm{d}x}=\mathrm{1}+\left(\sigma +\gamma +\mu \right)R_{\mathrm{0}}\frac{{\int }_{\mathrm{0}}^{\mathrm{\infty }}f\left(\tau \right)\tau {\mathrm{e}}^{-\left(x+\mu +\gamma +\sigma \right)\tau }\mathrm{d}\tau }{\tilde{f}}>\mathrm{1}

当R_{\mathrm{0}}>\mathrm{1}时，有，则V\left(x\right)=\mathrm{0}有正实根，故无病平衡态E_k^{\mathrm{0}}不稳定；当时，有V\left(\mathrm{0}\right)>\mathrm{0}，则V\left(x\right)=\mathrm{0}有负实根，故无病平衡态E_k^{\mathrm{0}}局部渐近稳定。

接着，考虑无病平衡态的全局稳定性。

引理1^［15］ 如果a>\mathrm{0},b>\mathrm{0}且\frac{\mathrm{d}x}{\mathrm{d}t}\ge b-ax，当t\ge \mathrm{0}且x\left(\mathrm{0}\right)\ge \mathrm{0}，有\underset{t\to +\mathrm{\infty }}{\mathrm{l}\mathrm{i}\mathrm{m}}\mathrm{i}\mathrm{n}\mathrm{f}x\left(t\right)\ge \frac{b}{a}。

如果a>\mathrm{0},b>\mathrm{0}且\frac{\mathrm{d}x}{\mathrm{d}t}\le b-ax，当t\ge \mathrm{0}且x\left(\mathrm{0}\right)\ge \mathrm{0}，有\underset{t\to +\mathrm{\infty }}{\mathrm{l}\mathrm{i}\mathrm{m}}\mathrm{s}\mathrm{u}\mathrm{p}x\left(t\right)\le \frac{b}{a}。

定理4 当R_{\mathrm{0}}\le \mathrm{1}，模型（2）的无病平衡态E_k^{\mathrm{0}}全局渐近稳定。

证 利用归一化条件，模型（2）可被重新写为

\left\{\begin{array}{l}\frac{\mathrm{d}{S}_k\left(t\right)}{\mathrm{d}t}=-\beta k{S}_k\left(t\right){\int }_{\mathrm{0}}^{\mathrm{\infty }}f\left(\tau \right){\mathrm{e}}^{-\left(\mu +\gamma +\sigma \right)\tau }\Theta \left(t-\tau \right)\mathrm{d}\tau -\left(\alpha +\mu +\delta \right)S_k\left(t\right)+\delta \left(\mathrm{1}-I_k\left(t\right)-Q_k\left(t\right)\right)+\mu \\ \frac{\mathrm{d}{I}_k\left(t\right)}{\mathrm{d}t}=\beta k{S}_k\left(t\right){\int }_{\mathrm{0}}^{\mathrm{\infty }}f\left(\tau \right){\mathrm{e}}^{-\left(\mu +\gamma +\sigma \right)\tau }\Theta \left(t-\tau \right)\mathrm{d}\tau -\left(\sigma +\gamma +\mu \right)I_k\left(t\right)\\ \frac{\mathrm{d}{Q}_k\left(t\right)}{\mathrm{d}t}=\sigma I_k\left(t\right)-\epsilon Q_k\left(t\right)-\mu Q_k\left(t\right)\end{array}\right.（11）

首先，证明\underset{t\to +\mathrm{\infty }}{\mathrm{l}\mathrm{i}\mathrm{m}}{I}_k\left(t\right)=\mathrm{0}。由（11）式的第一个式子得： \frac{\mathrm{d}{S}_k\left(t\right)}{\mathrm{d}t}\le \left(\mu +\delta \right)-\left(\mu +\alpha +\delta \right)S_k\left(t\right)。由引理2，有

\underset{t\to +\mathrm{\infty }}{\mathrm{l}\mathrm{i}\mathrm{m}}\mathrm{s}\mathrm{u}\mathrm{p}{S}_k\left(t\right)\le \frac{\mu +\delta }{\mu +\alpha +\delta }≕{\tilde{S}}_k^{\mathrm{0}}（12）

故对任意充分小的{\epsilon }_{\mathrm{1}}>\mathrm{0}，存在t_{\mathrm{1}}>\mathrm{0}使得S_k\left(t\right)\le {\tilde{S}}_k^{\mathrm{0}}+{\epsilon }_{\mathrm{1}}对所有的t>t_{\mathrm{1}}成立。因此，当t>t_{\mathrm{1}}时，有

\frac{\mathrm{d}{I}_k\left(t\right)}{\mathrm{d}t}\le \beta k\left({\tilde{S}}_k^{\mathrm{0}}+{\epsilon }_{\mathrm{1}}\right){\int }_{\mathrm{0}}^{\mathrm{\infty }}f\left(\tau \right){\mathrm{e}}^{-\left(\mu +\gamma +\sigma \right)\tau }\Theta \left(t-\tau \right)\mathrm{d}\tau -\left(\sigma +\gamma +\mu \right)I_k\left(t\right)

接着，构造一个满足M_k\left(\mathrm{0}\right)=I_k\left(\mathrm{0}\right)\ge \mathrm{0}的等价系统如下

\frac{\mathrm{d}{M}_k\left(t\right)}{\mathrm{d}t}=\beta k\left({\tilde{S}}_k^{\mathrm{0}}+{\epsilon }_{\mathrm{1}}\right){\int }_{\mathrm{0}}^{\mathrm{\infty }}f\left(\tau \right){\mathrm{e}}^{-\left(\mu +\gamma +\sigma \right)\tau }\tilde{\Theta }\left(t-\tau \right)\mathrm{d}\tau -\left(\sigma +\gamma +\mu \right)M_k\left(t\right)(13)

其中\tilde{\Theta }\left(t-\tau \right)=\frac{\mathrm{1}}{〈k〉}{\displaystyle \sum _{k=\mathrm{1}}^n}kP\left(k\right)M_k\left(t\right)。

构造李雅普诺夫函数N\left(t\right)={\displaystyle \sum _{k=\mathrm{1}}^n}\frac{kP\left(k\right)}{〈k〉\left(\sigma +\gamma +\mu \right)}{M}_k\left(t\right)，则

\frac{\mathrm{d}N\left(t\right)}{\mathrm{d}t}\left|{}_{\left(\mathrm{13}\right)}\right.={\displaystyle \sum _{k=\mathrm{1}}^n}\left[\frac{\beta k^{\mathrm{2}}P\left(k\right)}{〈k〉\left(\sigma +\gamma +\mu \right)}\left({\tilde{S}}_k^{\mathrm{0}}+{\epsilon }_{\mathrm{1}}\right){\int }_{\mathrm{0}}^{\mathrm{\infty }}f\left(\tau \right){\mathrm{e}}^{-\left(\mu +\gamma +\sigma \right)\tau }\tilde{\Theta }\left(t-\tau \right)\mathrm{d}\tau -\frac{kP\left(k\right)}{〈k〉}{M}_k\left(t\right)\right]=

\tilde{\Theta }\left(t-\tau \right)\left[R_{\mathrm{0}}+\frac{\beta 〈k^{\mathrm{2}}〉}{\left(\sigma +\gamma +\mu \right)〈k〉}{\epsilon }_{\mathrm{1}}{\int }_{\mathrm{0}}^{\infty }f\left(\tau \right)e^{-\left(\mu +\gamma +\sigma \right)\tau }d\tau -\mathrm{1}\right]

由于R_{\mathrm{0}}\le \mathrm{1}，故存在一个充分小的{\epsilon }_{\mathrm{1}}>\mathrm{0}使得R_{\mathrm{0}}+\frac{\beta 〈k^{\mathrm{2}}〉}{\left(\sigma +\gamma +\mu \right)〈k〉}{\epsilon }_{\mathrm{1}}{\int }_{\mathrm{0}}^{\infty }f\left(\tau \right)e^{-\left(\mu +\gamma +\sigma \right)\tau }d\tau \le \mathrm{1}，即\frac{\mathrm{d}N\left(t\right)}{\mathrm{d}t}\le \mathrm{0}对所有的M_k\left(t\right)\ge \mathrm{0}成立，并且\frac{\mathrm{d}N\left(t\right)}{\mathrm{d}t}=\mathrm{0}当且仅当\tilde{\Theta }\left(t-\tau \right)=\mathrm{0}，也就是M_k\left(t\right)=\mathrm{0}，因此\underset{t\to +\mathrm{\infty }}{\mathrm{l}\mathrm{i}\mathrm{m}}{M}_k\left(t\right)=\mathrm{0}。又因为\mathrm{0}\le I_k\left(t\right)\le M_k\left(t\right)对所有的t>\mathrm{0}成立，由比较原则可得\underset{t\to +\mathrm{\infty }}{\mathrm{l}\mathrm{i}\mathrm{m}}{I}_k\left(t\right)=\mathrm{0}。由（11）式的第三个式子易得\underset{t\to +\mathrm{\infty }}{\mathrm{l}\mathrm{i}\mathrm{m}}{Q}_k\left(t\right)=\mathrm{0}。

接着，证明\underset{t\to +\mathrm{\infty }}{\mathrm{l}\mathrm{i}\mathrm{m}}{S}_k\left(t\right)={\tilde{S}}_k^{\mathrm{0}}。

因为\underset{t\to +\mathrm{\infty }}{\mathrm{l}\mathrm{i}\mathrm{m}}{I}_k\left(t\right)=\mathrm{0}且\underset{t\to +\mathrm{\infty }}{\mathrm{l}\mathrm{i}\mathrm{m}}{Q}_k\left(t\right)=\mathrm{0}，则对任意充分小的{\epsilon }_{\mathrm{2}}>\mathrm{0}，存在t_{\mathrm{2}}>\mathrm{0}使得\mathrm{0}\le I_k\left(t\right)\le {\epsilon }_{\mathrm{2}}，\mathrm{0}\le Q_k\left(t\right)\le {\epsilon }_{\mathrm{2}}对t>t_{\mathrm{2}}成立。因此

\frac{\mathrm{d}{S}_k\left(t\right)}{\mathrm{d}t}\le \left(\mu +\delta \right)-\left[\mu +\alpha +\delta +\beta k{\epsilon }_{\mathrm{2}}{\int }_{\mathrm{0}}^{\mathrm{\infty }}f\left(\tau \right){\mathrm{e}}^{-\left(\mu +\gamma +\sigma \right)\tau }\frac{\mathrm{1}}{〈k〉}{\displaystyle \sum _{k=\mathrm{1}}^n}kP\left(k\right)\mathrm{d}\tau \right]S_k\left(t\right)-\mathrm{2}\delta {\epsilon }_{\mathrm{2}}

由比较原则得

\underset{t\to +\mathrm{\infty }}{\mathrm{l}\mathrm{i}\mathrm{m}}\mathrm{i}\mathrm{n}\mathrm{f}{S}_k\left(t\right)\ge \frac{\left(\mu +\delta \right)-\mathrm{2}\delta {\epsilon }_{\mathrm{2}}}{\mu +\alpha +\delta +\beta k{\epsilon }_{\mathrm{2}}{\int }_{\mathrm{0}}^{\mathrm{\infty }}f\left(\tau \right){\mathrm{e}}^{-\left(\mu +\gamma +\sigma \right)\tau }\frac{\mathrm{1}}{〈k〉}{\displaystyle \sum _{k=\mathrm{1}}^n}kP\left(k\right)\mathrm{d}\tau }

令{\epsilon }_{\mathrm{2}}\to \mathrm{0}有

\underset{t\to +\mathrm{\infty }}{\mathrm{l}\mathrm{i}\mathrm{m}}\mathrm{i}\mathrm{n}\mathrm{f}{S}_k\left(t\right)\ge \frac{\left(\mu +\delta \right)}{\mu +\alpha +\delta }={\tilde{S}}_k^{\mathrm{0}}（14）

结合（12）式和（14）式，有\underset{t\to +\mathrm{\infty }}{\mathrm{l}\mathrm{i}\mathrm{m}}{S}_k\left(t\right)=\frac{\mu +\delta }{\mu +\alpha +\delta }。

最后，由归一化条件S_k\left(t\right)+I_k\left(t\right)+Q_k\left(t\right)+R_k\left(t\right)=\mathrm{1}易得\underset{t\to +\mathrm{\infty }}{\mathrm{l}\mathrm{i}\mathrm{m}}{R}_k\left(t\right)=\frac{\alpha }{\mu +\alpha +\delta }。

综上，可得模型（2）的无病平衡态E_k^{\mathrm{0}}在R_{\mathrm{0}}\le \mathrm{1}时全局吸引。结合定理2易得无病平衡态E_k^{\mathrm{0}}在R_{\mathrm{0}}\le \mathrm{1}时全局渐近稳定。

2.2　地方病平衡态的稳定性

定理5 定义

当R_{\mathrm{0}}>\mathrm{1}且\frac{\epsilon +\mu }{\sigma }-\frac{\epsilon -\alpha }{\mathrm{2}\left(\gamma -\alpha \right)}>\mathrm{0}，\frac{\mu +\delta +\alpha }{\gamma }-\frac{\epsilon -\alpha }{\mathrm{2}\left(\gamma -\alpha \right)}>\mathrm{0}时，模型（2）的地方病平衡态E_k^{\mathrm{*}}在D_k内全局渐近稳定。

证 根据归一化条件S_k+I_k+Q_k+R_k=\mathrm{1}，得模型（2）的等价模型如下

\left\{\begin{array}{l}\frac{\mathrm{d}{I}_k}{\mathrm{d}t}=\beta k\left(\mathrm{1}-I_k-Q_k-R_k\right){\int }_{\mathrm{0}}^{\mathrm{\infty }}f\left(\tau \right){\mathrm{e}}^{-\left(\mu +\gamma +\sigma \right)\tau }\Theta \left(t-\tau \right)\mathrm{d}\tau -\left(\sigma +\gamma +\mu \right)I_k\\ \frac{\mathrm{d}{Q}_k}{\mathrm{d}t}=\sigma I_k-\left(\epsilon +\mu \right)Q_k\\ \frac{\mathrm{d}{R}_k}{\mathrm{d}t}=\alpha \left(\mathrm{1}-I_k-Q_k-R_k\right)+\gamma I_k+\epsilon Q_k-\left(\mu +\delta \right)R_k\end{array}\right.(15)

令\tilde{f}={\int }_{\mathrm{0}}^{\mathrm{\infty }}f\left(\tau \right){\mathrm{e}}^{-\left(\mu +\gamma +\sigma \right)\tau }\mathrm{d}\tau，对模型（15）作变换x_k=I_k-I_k^{\mathrm{*}},y_k=Q_k-Q_k^{\mathrm{*}},z_k=R_k-R_k^{\mathrm{*}}可得

\left\{\begin{array}{l}\frac{\mathrm{d}{x}_k}{\mathrm{d}t}=\beta k\left(\mathrm{1}-x_k-y_k-z_k+X_k^{\mathrm{*}}\right){\int }_{\mathrm{0}}^{\mathrm{\infty }}f\left(\tau \right){\mathrm{e}}^{-\left(\mu +\gamma +\sigma \right)\tau }\frac{\mathrm{1}}{〈k〉}{\displaystyle \sum _{k=\mathrm{1}}^n}kP\left(k\right)x_k\left(t-\tau \right)\mathrm{d}\tau \\ -\left(\beta k\frac{\mathrm{1}}{〈k〉}{\displaystyle \sum _{k=\mathrm{1}}^n}kP\left(k\right)I_k^{\mathrm{*}}\tilde{f}+\sigma +\gamma +\mu \right)x_k-\beta k\frac{\mathrm{1}}{〈k〉}{\displaystyle \sum _{k=\mathrm{1}}^n}kP\left(k\right)I_k^{\mathrm{*}}\tilde{f}\left(y_k+z_k\right)\\ \frac{\mathrm{d}{y}_k}{\mathrm{d}t}=\sigma x_k-\left(\epsilon +\mu \right)y_k\\ \frac{\mathrm{d}{z}_k}{\mathrm{d}t}=\left(\gamma -\alpha \right)x_k+\left(\epsilon -\alpha \right)y_k-\left(\mu +\delta +\alpha \right)z_k\end{array}\right.(16)

其中X_k^{\mathrm{*}}=-I_k^{\mathrm{*}}-Q_k^{\mathrm{*}}-R_k^{\mathrm{*}}。

定义李雅普诺夫函数

V_k=\frac{\mathrm{1}}{\mathrm{2}}{x}_k^{\mathrm{2}}+\frac{\beta k\frac{\mathrm{1}}{〈k〉}{\displaystyle \sum _{k=\mathrm{1}}^n}kP\left(k\right)I_k^{\mathrm{*}}\tilde{f}}{\mathrm{2}\sigma }{y}_k^{\mathrm{2}}+\frac{\beta k\frac{\mathrm{1}}{〈k〉}{\displaystyle \sum _{k=\mathrm{1}}^n}kP\left(k\right)I_k^{\mathrm{*}}\tilde{f}}{\mathrm{2}\left(\gamma -\alpha \right)}{z}_k^{\mathrm{2}}+

\frac{\beta k\frac{\mathrm{1}}{〈k〉}{\displaystyle \sum _{k=\mathrm{1}}^n}kP\left(k\right)I_k^{\mathrm{*}}\tilde{f}+\sigma +\gamma +\mu }{\mathrm{2}}{\int }_{t-\tau }^t{x}_k^{\mathrm{2}}\left(v\right)\mathrm{d}v

则

\mathrm{d}{V}_k=\beta k\left(\mathrm{1}-x_k-y_k-z_k+X_k^{\mathrm{*}}\right)\frac{\mathrm{1}}{〈k〉}{\displaystyle \sum _{k=\mathrm{1}}^n}kP\left(k\right)\tilde{f}{x}_k\left(t\right)x_k\left(t-\tau \right)-

x_k^{\mathrm{2}}\left[\beta k\frac{\mathrm{1}}{〈k〉}{\displaystyle \sum _{k=\mathrm{1}}^n}kP\left(k\right)I_k^{\mathrm{*}}\tilde{f}+\left(\sigma +\gamma +\mu \right)\right]-\beta k\frac{\mathrm{1}}{〈k〉}{\displaystyle \sum _{k=\mathrm{1}}^n}kP\left(k\right)I_k^{\mathrm{*}}\tilde{f}{x}_k\left(y_k+z_k\right)

+\frac{\beta k\frac{\mathrm{1}}{〈k〉}{\displaystyle \sum _{k=\mathrm{1}}^n}kP\left(k\right)I_k^{\mathrm{*}}\tilde{f}}{\sigma }{y}_k\left[\sigma x_k-\left(\epsilon +\mu \right)y_k\right]+

\frac{\beta k\frac{\mathrm{1}}{〈k〉}{\displaystyle \sum _{k=\mathrm{1}}^n}kP\left(k\right)I_k^{\mathrm{*}}\tilde{f}}{\gamma -\alpha }{z}_k\left[\left(\gamma -\alpha \right)x_k+\left(\epsilon -\alpha \right)y_k-\left(\mu +\delta +\alpha \right)z_k\right]+

\frac{\beta k\frac{\mathrm{1}}{〈k〉}{\displaystyle \sum _{k=\mathrm{1}}^n}kP\left(k\right)I_k^{\mathrm{*}}\tilde{f}+\left(\sigma +\gamma +\mu \right)}{\mathrm{2}}{x}_k^{\mathrm{2}}-\frac{\beta k\frac{\mathrm{1}}{〈k〉}{\displaystyle \sum _{k=\mathrm{1}}^n}kP\left(k\right)I_k^{\mathrm{*}}\tilde{f}+\left(\sigma +\gamma +\mu \right)}{\mathrm{2}}{x}_k^{\mathrm{2}}\left(t-\tau \right)=

-\frac{\beta k\frac{\mathrm{1}}{〈k〉}{\displaystyle \sum _{k=\mathrm{1}}^n}kP\left(k\right)I_k^{\mathrm{*}}\tilde{f}+\left(\sigma +\gamma +\mu \right)}{\mathrm{2}}{x}_k^{\mathrm{2}}+\beta k\left(\mathrm{1}-x_k-y_k-z_k+X_k^{\mathrm{*}}\right)\frac{\mathrm{1}}{〈k〉}{\displaystyle \sum _{k=\mathrm{1}}^n}kP\left(k\right)\tilde{f}{x}_k\left(t\right)x_k\left(t-\tau \right)-

\frac{\beta k\frac{\mathrm{1}}{〈k〉}{\displaystyle \sum _{k=\mathrm{1}}^n}kP\left(k\right)I_k^{\mathrm{*}}\tilde{f}+\left(\sigma +\gamma +\mu \right)}{\mathrm{2}}{x}_k^{\mathrm{2}}\left(t-\tau \right)-\frac{\left(\epsilon +\mu \right)\beta k\frac{\mathrm{1}}{〈k〉}{\displaystyle \sum _{k=\mathrm{1}}^n}kP\left(k\right)I_k^{\mathrm{*}}\tilde{f}}{\sigma }{y}_k^{\mathrm{2}}-

\frac{\left(\mu +\delta +\alpha \right)\beta k\frac{\mathrm{1}}{〈k〉}{\displaystyle \sum _{k=\mathrm{1}}^n}kP\left(k\right)I_k^{\mathrm{*}}\tilde{f}}{\gamma -\alpha }{z}_k^{\mathrm{2}}+\frac{\left(\epsilon -\alpha \right)\beta k\frac{\mathrm{1}}{〈k〉}{\displaystyle \sum _{k=\mathrm{1}}^n}kP\left(k\right)I_k^{\mathrm{*}}\tilde{f}}{\gamma -\alpha }{y}_k{z}_k\le

-\frac{\beta k\frac{\mathrm{1}}{〈k〉}{\displaystyle \sum _{k=\mathrm{1}}^n}kP\left(k\right)I_k^{\mathrm{*}}\tilde{f}+\left(\sigma +\gamma +\mu \right)}{\mathrm{2}}{x}_k^{\mathrm{2}}+\beta k\left(\mathrm{1}-x_k-y_k-z_k+X_k^{\mathrm{*}}\right)\frac{\mathrm{1}}{〈k〉}{\displaystyle \sum _{k=\mathrm{1}}^n}kP\left(k\right)\tilde{f}{x}_k\left(t\right)x_k\left(t-\tau \right)-

\frac{\beta k\frac{\mathrm{1}}{〈k〉}{\displaystyle \sum _{k=\mathrm{1}}^n}kP\left(k\right)I_k^{\mathrm{*}}\tilde{f}+\left(\sigma +\gamma +\mu \right)}{\mathrm{2}}{x}_k^{\mathrm{2}}\left(t-\tau \right)-\beta k\frac{\mathrm{1}}{〈k〉}{\displaystyle \sum _{k=\mathrm{1}}^n}kP\left(k\right)I_k^{\mathrm{*}}\tilde{f}\left[\frac{\epsilon +\mu }{\sigma }-\frac{\epsilon -\alpha }{\mathrm{2}\left(\gamma -\alpha \right)}\right]y_k^{\mathrm{2}}-

\beta k\frac{\mathrm{1}}{〈k〉}{\displaystyle \sum _{k=\mathrm{1}}^n}kP\left(k\right)I_k^{\mathrm{*}}\tilde{f}\left[\frac{\mu +\delta +\alpha }{\gamma -\alpha }-\frac{\epsilon -\alpha }{\mathrm{2}\left(\gamma -\alpha \right)}\right]z_k^{\mathrm{2}}=

-\frac{\mathrm{1}}{\mathrm{2}}\left[{\mathit{w}}_k^{\mathrm{T}}\left(t,\tau \right){\mathit{A}}_k\left(t\right){\mathit{w}}_k\left(t,\tau \right)\right]-\beta k\frac{\mathrm{1}}{〈k〉}{\displaystyle \sum _{k=\mathrm{1}}^n}kP\left(k\right)I_k^{\mathrm{*}}\tilde{f}\left[\frac{\epsilon +\mu }{\sigma }-\frac{\epsilon -\alpha }{\mathrm{2}\left(\gamma -\alpha \right)}\right]y_k^{\mathrm{2}}

-\beta k\frac{\mathrm{1}}{〈k〉}{\displaystyle \sum _{k=\mathrm{1}}^n}kP\left(k\right)I_k^{\mathrm{*}}\tilde{f}\left[\frac{\mu +\delta +\alpha }{\gamma -\alpha }-\frac{\epsilon -\alpha }{\mathrm{2}\left(\gamma -\alpha \right)}\right]z_k^{\mathrm{2}}

其中，

{\mathit{w}}_k^{\mathrm{T}}=\mathrm{c}\mathrm{o}\mathrm{l}\left(x_k\left(t\right),x_k\left(t-\tau \right)\right)

{\mathit{A}}_k=\left[\begin{array}{cc}\beta k\frac{\mathrm{1}}{〈k〉}{\displaystyle \sum _{k=\mathrm{1}}^n}kP\left(k\right)I_k^{\mathrm{*}}\tilde{f}+\sigma +\gamma +\mu & -\beta k\left(\mathrm{1}-x_k-y_k-z_k+X_k^{\mathrm{*}}\right)\tilde{f}\frac{\mathrm{1}}{〈k〉}{\displaystyle \sum _{k=\mathrm{1}}^n}kP\left(k\right)\\ -\beta k\left(\mathrm{1}-x_k-y_k-z_k+X_k^{\mathrm{*}}\right)\tilde{f}\frac{\mathrm{1}}{〈k〉}{\displaystyle \sum _{k=\mathrm{1}}^n}kP\left(k\right)& \beta k\frac{\mathrm{1}}{〈k〉}{\displaystyle \sum _{k=\mathrm{1}}^n}kP\left(k\right)I_k^{\mathrm{*}}\tilde{f}+\sigma +\gamma +\mu \end{array}\right]

显然，当时，矩阵{\mathit{A}}_k是正定的。对任意给定的\epsilon >\mathrm{0}，定义

对区域D_k^{\epsilon }上的任意点\left(I_k,Q_k,R_k\right)，用{\lambda }_k表示矩阵{\mathit{A}}_k的最小特征值，则{\lambda }_k>\mathrm{0}，且有{\mathit{w}}_k^{\mathrm{T}}\left(t,\tau \right){\mathit{A}}_k\left(t\right){\mathit{w}}_k\left(t,\tau \right)\ge {\lambda }_k\left(x_k^{\mathrm{2}}\left(t\right)+x_k^{\mathrm{2}}\left(t-\tau \right)\right)，又由于\frac{\epsilon +\mu }{\sigma }-\frac{\epsilon -\alpha }{\mathrm{2}\left(\gamma -\alpha \right)}>\mathrm{0}，\frac{\mu +\delta +\alpha }{\gamma }-\frac{\epsilon -\alpha }{\mathrm{2}\left(\gamma -\alpha \right)}>\mathrm{0}，从而

\mathrm{d}{V}_k\le -\frac{\mathrm{1}}{\mathrm{2}}{\lambda }_k\left(x_k^{\mathrm{2}}\left(t\right)+x_k^{\mathrm{2}}\left(t-\tau \right)\right)-\beta k\frac{\mathrm{1}}{〈k〉}{\displaystyle \sum _{k=\mathrm{1}}^n}kP\left(k\right)I_k^{\mathrm{*}}\tilde{f}\left[\frac{\epsilon +\mu }{\sigma }-\frac{\epsilon -\alpha }{\mathrm{2}\left(\gamma -\alpha \right)}\right]y_k^{\mathrm{2}}

-\beta k\frac{\mathrm{1}}{〈k〉}{\displaystyle \sum _{k=\mathrm{1}}^n}kP\left(k\right)I_k^{\mathrm{*}}\tilde{f}\left[\frac{\mu +\delta +\alpha }{\gamma -\alpha }-\frac{\epsilon -\alpha }{\mathrm{2}\left(\gamma -\alpha \right)}\right]z_k^{\mathrm{2}}\le \mathrm{0}

故模型（16）在区域D_k^{\epsilon }内全局渐近稳定，由\epsilon的任意性可知，模型（16）在区域D_k内全局渐近稳定，因此模型（2）的地方病平衡态E_k^{\mathrm{*}}在区域D_k内全局渐近稳定。

注1 由

I_k^{\mathrm{*}}=\frac{\left(\mu +\delta \right)\left(\mu +\epsilon \right)\beta k\tilde{f}{\Theta }^{\mathrm{*}}}{\left[\sigma \delta +\left(\mu +\epsilon \right)\left(\delta +\sigma +\gamma +\mu \right)\right]\beta k\tilde{f}{\Theta }^{\mathrm{*}}+\left(\mu +\epsilon \right)\left(\sigma +\gamma +\mu \right)\left(\delta +\alpha +\mu \right)}

其中，\tilde{f}={\int }_{\mathrm{0}}^{\mathrm{\infty }}f\left(\tau \right){\mathrm{e}}^{-\left(\mu +\gamma +\sigma \right)\tau }\mathrm{d}\tau，可知当R_{\mathrm{0}}>\mathrm{1}时，时滞促进了疾病的传播，增加了疾病的最终感染规模。

3　 数值模拟

首先，利用BA算法^［1］生成一个含有800个节点的无标度网络Q，生成参数为m_{\mathrm{0}}=\mathrm{5},m=\mathrm{1},下面的数值模拟基于网络Q进行。

例1 选取参数\mu =\mathrm{0.5},\delta =\mathrm{0.7},\sigma =\mathrm{0.8},\gamma =\mathrm{0.5},\alpha =\mathrm{0.9},\beta =\mathrm{0.6},\epsilon =\mathrm{0.9},f\left(\tau \right)=\mathrm{0.5}{\mathrm{e}}^{-\mathrm{0.5}\tau },\mathrm{0}\le \tau \le \mathrm{\infty }，得。易感者、感染者、隔离者、恢复者随时间变化的轨迹如图1。从图1可以看出无病平衡态全局渐近稳定，这说明了定理3的有效性。当\tau分别取3、4、5，感染者的密度随时间的变化轨迹如图2，疾病最终消失。

例2 选取参\beta =\mathrm{0.6},\mu =\mathrm{0.05},\delta =\mathrm{0.08},\sigma =\mathrm{0.06},\gamma =\mathrm{0.15},\alpha =\mathrm{0.1},\epsilon =\mathrm{0.1},f\left(\tau \right)=\mathrm{0.5}{\mathrm{e}}^{-\mathrm{0.5}\tau },\mathrm{0}\le \tau \le \mathrm{\infty }，得R_{\mathrm{0}}=\frac{〈k^{\mathrm{2}}〉}{〈k〉}\frac{\left(\mu +\delta \right)\beta \tilde{f}}{\left(\sigma +\gamma +\mu \right)\left(\mu +\alpha +\delta \right)}\approx \mathrm{5.89}>\mathrm{1}。易感者、感染者、隔离者、恢复者随时间变化的轨迹如图3。从图3可看出地方病平衡态全局渐近稳定，这说明了定理4的有效性。当\tau分别取3、4、5，感染者的密度随时间的变化轨迹如图4，疾病永久存在。

4　 结 语

本文综合考虑了接触异质性、隔离措施、疫苗注射、时滞等因素，建立了更加符合实际的无标度网络上带有隔离项的时滞SIQRS模型。分析了网络的拓扑结构、时滞、隔离措施等因素对基本再生数R_{\mathrm{0}}大小的影响，发现基本再生数R_{\mathrm{0}}与网络的异质性〈k^{\mathrm{2}}〉/〈k〉呈正相关、与隔离率\sigma呈负相关、与疾病的感染性f\left(\tau \right)呈正相关。通过构造合适的李雅普诺夫函数证明了平衡态的全局稳定性，结果表明：当R_{\mathrm{0}}\le \mathrm{1}时，无病平衡态全局渐近稳定，即疾病最终消失；否则地方病平衡态全局渐近稳定，即疾病永久存在。利用MATLAB软件对分析结果进行了数值模拟，发现当时，时滞减缓了疾病的传播；当R_{\mathrm{0}}>\mathrm{1}时，时滞促进了疾病的传播。

虽然隔离措施可以降低基本再生数，但由于隔离措施会对经济发展产生消极的影响，且当基本再生数R_{\mathrm{0}}\le \mathrm{1}时，隔离措施无法起到控制疾病传播的作用，只有当基本再生数R_{\mathrm{0}}>\mathrm{1}时，隔离措施可以在一定程度控制疾病的传播。因此，选择恰当的时机实施隔离措施既能有效地控制疾病的传播，也不会对经济发展产生太大的影响。

由于本文所建模型耦合度较高，所以如何构造合适的李雅普诺夫函数来分析其动力学行为有一定的挑战性，故本文利用归一化条件和线性变换将模型进行降维及线性化处理。其次，由于复杂网络上的传染病模型维度较高，因此如何从理论上给出复杂网络上的时滞传染病模型Hopf分支存在的条件是一个值得研究的课题。事实上，由于生物系统中环境的变化，疾病在传播过程中或多或少都会受到随机噪声的影响，所以将随机干扰引入到模型中是随后的一个研究方向。

参考文献

[1] BARABÁSI A L， ALBERT R. Emergence of scaling in random networks ［J］. Science， 1999， 286（5439）： 509-512. DOI： 10.1126/science.286.5439.509 .

[2] ALBERT R， JEONG H， BARABÁSI A L. Diameter of the World Wide Web ［J］. Nature， 1999， 401（6）： 130-131. DOI： 10.1038/43601 .

[3] FALOUTSOS M， FALOUTSOS P， FALOUTSOS C. On power-law relationships of the internet topology ［J］. Computer Communication Review， 1999， 29（4）： 251-262. DOI： 10.1145/316194.316229 .

[4] LIU X D， LI T， WANG Y M， et al. An SIS epidemic model with infective medium and feedback mechanism on scale-free networks ［J］. Open Access Library Journal， 2017， 4（5）： 1-9. DOI：10.4236/oalib.1103598 .

[5] SAMIRA M， NAHID M J. Epidemic spreading phenomena on a scale-free network with time varying transmission rate due to social responses ［J］. International Journal of Modern Physics C， 2020， 31： 2050148. DOI： 10.1142/s012918312050148x .

[6] PASTOR-SATORRAS R， VESPIGNANI A. Epidemic spreading in scale-free networks ［J］. Physical Review Letters， 2001， 86（14）： 3200-3203. DOI：10.1103/PhysRevLett.86.3200 .

[7] LI T， WANG Y M， GUAN Z H. Spreading dynamics of a SIQRS epidemic model on scale-free networks ［J］. Communications in Nonlinear Science and Numerical Simulation， 2014， 19（3）： 686-692. DOI：10.1016/j.cnsns.2013.07.010 .

[8] LAHROUZ A. Dynamics of a delayed epidemic model with varying immunity period and nonlinear transmission ［J］. International Journal of Biomathematics， 2015， 8（2）： 16. DOI：10.1142/S1793524515500278 .

[9] SCHEINER S， UKAJ N， HELLMICH C. Mathematical modeling of COVID-19 fatality trends： Death kinetics law versus infection-to-death delay rule ［J］. Chaos， Solitons & Fractals， 2020， 136： 109891. DOI：10.1016/j.chaos.2020.109891 .

[10] SUN M F， ZHANG H F， KANG H Y， et al. Epidemic spreading on adaptively weighted scale-free networks ［J］. Journal of Mathematical Biology， 2017， 74（5）： 1263-1298. DOI：10.1007/s00285-016-1057-6 .

[11] LIU J. Bifurcation of a delayed SEIS epidemic model with a changing delitescence and nonlinear incidence rate ［J］. Discrete Dynamics in Nature and Society， 2017， 2017： 1-9. DOI：10.1155/2017/2340549 .

[12] YANG P， WANG Y S. Dynamics for an SEIRS epidemic model with time delay on a scale-free network ［J］. Physica A： Statistical Mechanics and Its Applications， 2019， 527： 121290. DOI：10.1016/j.physa.2019.121290 .

[13] HAJJAMI R， JARROUDI M E， LAHROUZ A， et al. Dynamic analysis of an SEIR epidemic model with a time lag in awareness allocated funds ［J］. Discrete and Continuous Dynamical Systems-Series B， 2021，26（8）： 4191-4225. DOI：10.3934/dcdsb.2020285 .

[14] HUANG C X， CAO J， WEN F H， et al. Stability analysis of SIR model with distributed delay on complex networks ［J］. PLoS One， 2016， 11（8）： e0158813. DOI：10.1371/journal.pone.0158813 .

[15] CHEN F D. On a nonlinear nonautonomous predator-prey model with diffusion and distributed delay ［J］. Journal of Computational and Applied Mathematics， 2005， 180（1）： 33-49. DOI：10.1016/j.cam.2004.10.001 .

基金资助

国家自然科学基金(11361044)

宁夏自然科学基金(2021AAC03030)

宁夏大学研究生创新项目(GIP2020060)