加权网络上一类SIRS传染病模型的动力学分析

Dynamics Analysis of an SIRS Epidemic Model on Weighted Networks

摘要: 本文主要研究基于权重网络的一类SIRS传染病模型的动力学行为，发现基本再生数R₀决定了传染疾病的传播与否。即，若R₀<1则无病平衡点E₀是全局渐近稳定性的；若R₀>1则地方病平衡点E*是存在唯一的，并且此时系统也是准持久的。最后通过数值仿真，验证了理论结果的正确性。

Abstract: In this paper, an improved SIRS epidemic model on weighted networks is investigated and we ob-tain that the basic reproduction number R₀ determines whether a disease persists, i.e., if R₀<1, then the disease-free equilibrium E₀ is globally asymptotically stable; if R₀>1, there exists a unique endemic equilibrium E*. Furthermore, the model is quasi-persistent if R₀>1. Finally, some simulations are presented to demonstrate the correctness of the theoretical results.

1. 引言

传染病一直以来都是人类生存的头号大敌，不仅危害着人们的身体健康而且影响社会的稳定。近代以来，人类在科技和医药技术方面取得了巨大的进步，但一些传染病仍然对我们的生命安全产生巨大威胁，如禽流感、甲型H1N1流感以及寨卡病毒(Zika Virus)等。因此，弄清传染病的传播机制，怎样有效地预防和控制疾病的传播成为了科学研究的焦点问题。

众所周知，数学建模是研究传染病传播机制的有效手段之一。有关这方面的研究最早可追溯到上世纪初。早在1906年，Hamer在研究麻疹时构建并分析了一个离散时间模型 [1]。随后在1926年，Kermack和McKendrick在研究黑死病以及瘟疫的传播规律时，建立了易感态–感染态–移除态(susceptible-infectious-recovered, SIR)仓室模型，之后在1932年，又构建了易感态–感染态–易感态(susceptible-infectious-susceptible, SIS)传染病模型并提出了揭示传染病传播规律的“阈值定理” [2] [3] ，这些成果至今仍然被广泛的应用并且得到了蓬勃发展。但是由于传统仓室模型过于理想化，将人群看作均匀混合的，忽略了个体之间接触的异质性。随着复杂网络理论兴起和不断发展，上述不足有了更好的解决方法。复杂网络为描述社会、经济和生物等系统的拓扑结构提供了一个强有力的框架并且异质网络的度可很好地刻画接触模式的异质性 [4]。因此将复杂网络和传染病模型结合起来，研究传染病在复杂网络上的传播动力学性质成为了近些年来的研究热点 [4] [5] [6] [7] [8]。在2001年，Pastor-Satorras和Vespiganai首次基于无标度网络提出并研究了SIS传染病模型 [5] ，结果表明在有限规模网络上，传染病的传播阈值为 $\langle k\rangle /\langle k^2\rangle$ ，而在无限规模高度异质的网络上，其传播阈值为0。Li等人在复杂网络上构建并研究了一个新的具有非线性感染率和具有出生率和死亡率的易感态–感染态–移除态–易感态(susceptible-infectious-recovered-susceptible, SIRS)模型 [6] ，Zhu等人在文献 [7] 中提出了复杂网络上的一个一般化的传染病模型，他们的结果都解释了连接异质性是如何影响传播阈值的。

然而接触模式的异质性在复杂网络中不仅表现为连边的数量不同，而且也体现在连边权重的差异上 [9]。事实上，许多常见的网络系统都是带权重的，例如交通网络、股票网络、生物网络等。在传染病传播过程中，边权可以表示两个个体之间接触的频繁程度或接触时间长短，边权越大，表示个体之间的连接越紧密，易感状态的个体被感染的几率就越大 [8] [10]。因此研究边的权重对传播过程的影响具有重要的意义。文章 [11] 利用权重来衡量个体行为对传播过程的影响，并由此估计个体的感染风险，从而采取更有效的干预措施。文章 [12] 通过边权重来表示节点之间的亲密程度，并且通过研究发现权重较大的节点更容易被感染，而对于整个网络来说，连边的权重越分散，疾病传播就会越慢。

受上述的研究成果启发，本文基于无标度网络SIRS模型，考虑感染个体的差异性，将感染个体分为低风险感染个体和高风险感染个体，并引入权重用于刻画个体之间接触的频繁程度，建立了基于权重网络的SIRS传染病模型。

2. 传染病模型

记 $S_k\left(t\right)\mathrm{,\; <mi>\; </mi>}{I}_{1k}\left(t\right)\mathrm{,\; <mi>\; </mi>}{I}_{2k}\left(t\right)\mathrm{,\; <mi>\; </mi>}{R}_k\left(t\right)$ 分别表示为t时刻度为k的易感个体的密度，低风险染病个体的密度，

高风险染病个体的密度以及康复个体的密度，其中 $k=1,2,\cdots ,n$ (n为网络的最大度)表示该节点与其他节点之间的连边数。根据异质平均场理论，可得到如下的SIRS传染病模型

$(\begin{array}{l}\frac{\text{d}{S}_k\left(t\right)}{\text{d}t}=-k{S}_k\left(t\right)\Theta \left(t\right)+\delta R_k\left(t\right)\mathrm{,}\\ \frac{\text{d}{I}_{1k}\left(t\right)}{\text{d}t}=\alpha k{S}_k\left(t\right)\Theta \left(t\right)-{\gamma }_1{I}_{1k}\left(t\right)\mathrm{,}\\ \frac{\text{d}{I}_{2k}\left(t\right)}{\text{d}t}=\left(1-\alpha \right)k{S}_k\left(t\right)\Theta \left(t\right)-{\gamma }_2{I}_{2k}\left(t\right)\mathrm{,}\\ \frac{\text{d}{R}_k\left(t\right)}{\text{d}t}={\gamma }_1{I}_{1k}\left(t\right)+{\gamma }_2{I}_{2k}\left(t\right)-\delta R_k\left(t\right)\mathrm{,}\end{array}$ (1)

其中 $k{S}_k\left(t\right)\Theta \left(t\right)$ 表示单位时间内新感染的度为k的个体密度， $\alpha$ 代表新的感染个体为低风险感染个体的比率。 ${\gamma }_1\mathrm{,\; <mi>\; </mi>}{\gamma }_2$ 分别表示低风险染病个体，高风险染病个体的治愈率。 $\delta$ 为康复个体的免疫丧失率。 $\Theta \left(t\right)={\Theta }_1\left(t\right)+{\Theta }_2\left(t\right)$ 表示为一个易感个体被一个感染的邻居节点感染的概率，并且

$\begin{array}{l}{\Theta }_1\left(t\right)={\displaystyle \underset{i=1}{\overset{n}{\sum }}}P\left(i|k\right)\frac{{\phi }_1\left(i\right)}{i}{\lambda }_{ik}{I}_{1i}\left(t\right)=\frac{1}{\langle k\rangle }{\displaystyle \underset{i=1}{\overset{n}{\sum }}}{\lambda }_{ik}{\phi }_1\left(i\right)P\left(i\right)I_{1i}\left(t\right),\\ {\Theta }_2\left(t\right)={\displaystyle \underset{i=1}{\overset{n}{\sum }}}P\left(i|k\right)\frac{{\phi }_2\left(i\right)}{i}{\lambda }_{ik}{I}_{1i}\left(t\right)=\frac{1}{\langle k\rangle }{\displaystyle \underset{i=1}{\overset{n}{\sum }}}{\lambda }_{ik}{\phi }_2\left(i\right)P\left(i\right)I_{2i}\left(t\right),\end{array}$ (2)

其中 $P\left(i|k\right)$ 是度为k的节点连向度为i的节点的概率。在度不相关的网络中 [13] ， $P\left(i|k\right)=iP\left(i\right)/\langle k\rangle$ 。 ${\phi }_1\left(k\right)\mathrm{,\; <mi>\; </mi>}{\phi }_2\left(k\right)$ 分别表示度为k的低风险感染，高风险感染节点的感染力， ${\lambda }_{ik}$ 代表度为i的节点到度为k的节点之间的感染率。 $P\left(k\right)$ 是网络中节点的度分布且 ${\displaystyle {\sum }_{k=1}^n}P\left(k\right)=1$ 。 $\langle k\rangle ={\displaystyle {\sum }_{k=1}^n}kP\left(k\right)$ 则为网络的平均度。对于一般的函数 $f\left(k\right)$ ，记 $\langle f\left(k\right)\rangle ={\displaystyle {\sum }_{k=1}^{n}f\left(k\right)P\left(k\right)}$ 。

对疾病传播而言，权重用来表示传染病传播的强弱程度，因此可以用节点的度来表示连边的权重，若度为i和j的两个节点相连，那么它们之间的权重可以用一个关于i, j的函数 $\varphi \left(i,j\right)$ 来表示 [14] [15] [16]。

其一般表达式为 $\varphi \left(i,j\right)=g\left(i\right)g\left(j\right)$ 。因为度越大，表示节点影响力越大，权重也会越大，所以 $g\left(k\right)$ 是一个关于k的增函数。函数 $\varphi \left(i,j\right)$ 表示的是连边之间的权重，类似的也可建立有关节点的权重表达式 [8] ：设 ${\Psi }_k$ 表示度为k的节点的权重，代表所有与它相连的边的权重之和，故 ${\Psi }_k=k{\displaystyle {\sum }_{i=i}^n}P\left(i|k\right)\varphi \left(i,j\right)$ ，并且在度不相关网络上，有 ${\Psi }_k=kg\left(k\right)\langle kg\left(k\right)\rangle /\langle k\rangle$ 。

设度为i的节点有一个固定的传播率 $\lambda i$ ，并且通过一条从度为i的节点到度为k的节点的连边的传播概率可由这条边的权重占节点权重的比例表示 [8] ，由此可得到感染率的表达式

${\lambda }_{ik}=\lambda i\frac{\varphi \left(i,k\right)}{{\Psi }_i}=\frac{\lambda g\left(k\right)\langle k\rangle }{\langle kg\left(k\right)\rangle }.$ (3)

将(3)代入(2)，得

${\bar{\Theta }}_1\left(t\right)=\frac{\lambda g\left(k\right)}{\langle kg\left(k\right)\rangle }{\displaystyle \underset{i=1}{\overset{n}{\sum }}}{\phi }_1\left(i\right)P\left(i\right)I_{1i}\left(t\right),\mathrm{<mi>\; </mi>\; <mi>\; </mi>}{\bar{\Theta }}_2\left(t\right)=\frac{\lambda g\left(k\right)}{\langle kg\left(k\right)\rangle }{\displaystyle \underset{i=1}{\overset{n}{\sum }}}{\phi }_2\left(i\right)P\left(i\right)I_{2i}\left(t\right).$ (4)

然后将(4)代(1)，得到如下系统

$(\begin{array}{l}\frac{\text{d}{S}_k\left(t\right)}{\text{d}t}=-\frac{\lambda k{S}_k\left(t\right)g\left(k\right)}{\langle kg\left(k\right)\rangle }\theta \left(t\right)+\delta R_k\left(t\right)\mathrm{,}\hfill \\ \frac{\text{d}{I}_{1k}\left(t\right)}{\text{d}t}=\frac{\alpha \lambda k{S}_k\left(t\right)g\left(k\right)}{\langle kg\left(k\right)\rangle }\theta \left(t\right)-{\gamma }_1{I}_{1k}\left(t\right)\mathrm{,}\hfill \\ \frac{\text{d}{I}_{2k}\left(t\right)}{\text{d}t}=\frac{\left(1-\alpha \right)\lambda k{S}_k\left(t\right)g\left(k\right)}{\langle kg\left(k\right)\rangle }\theta \left(t\right)-{\gamma }_2{I}_{2k}\left(t\right)\mathrm{,}\hfill \\ \frac{\text{d}{R}_k\left(t\right)}{\text{d}t}={\gamma }_1{I}_{1k}\left(t\right)+{\gamma }_2{I}_{2k}\left(t\right)-\delta R_k\left(t\right)\mathrm{,}\hfill \end{array}$ (5)

其中 $\theta \left(t\right)={\theta }_1\left(t\right)+{\theta }_2\left(t\right)$ ，且 ${\theta }_1\left(t\right)={\displaystyle {\sum }_{i=1}^n}{\phi }_1\left(i\right)P\left(i\right)I_{1i}\left(t\right)$ ， ${\theta }_2\left(t\right)={\displaystyle {\sum }_{i=1}^n}{\phi }_2\left(i\right)P\left(i\right)I_{2i}\left(t\right)$ 。记 $Q_k\left(t\right)=S_k\left(t\right)+I_{1k}\left(t\right)+I_{2k}\left(t\right)+R_k\left(t\right)$ ，则有 $\frac{\text{d}{Q}_k\left(t\right)}{\text{d}t}=0$ ，即总的个体数是常数，因此不妨设 $S_k\left(t\right)+I_{1k}\left(t\right)+I_{2k}\left(t\right)+R_k\left(t\right)=1$ 。且系统(5)的初始条件具有如下形式：

$\begin{array}{l}0\le S_k\left(0\right),I_{1k}\left(0\right),I_{2k}\left(0\right),R_k\left(0\right)\le 1,\\ S_k\left(0\right)+I_{1k}\left(0\right)+I_{2k}\left(0\right)+R_k\left(0\right)=1,\mathrm{<mi>\; </mi>}\theta \left(0\right)>0,\mathrm{<mi>\; </mi>}k=1,2,\cdots ,n.\end{array}$

由于 $S_k\left(t\right)=1-I_{1k}\left(t\right)-I_{2k}\left(t\right)-R_k\left(t\right)$ ，因而系统(5)等价于如下系统：

$(\begin{array}{l}\frac{\text{d}{I}_{1k}\left(t\right)}{\text{d}t}=\frac{\alpha \lambda k\left(1-I_{1k}\left(t\right)-I_{2k}\left(t\right)-R_k\left(t\right)\right)g\left(k\right)}{\langle kg\left(k\right)\rangle }\theta \left(t\right)-{\gamma }_1{I}_{1k}\left(t\right)\mathrm{,}\\ \frac{\text{d}{I}_{2k}\left(t\right)}{\text{d}t}=\frac{\left(1-\alpha \right)\lambda k\left(1-I_{1k}\left(t\right)-I_{2k}\left(t\right)-R_k\left(t\right)\right)g\left(k\right)}{\langle kg\left(k\right)\rangle }\theta \left(t\right)-{\gamma }_2{I}_{2k}\left(t\right)\mathrm{,}\\ \frac{\text{d}{R}_k\left(t\right)}{\text{d}t}={\gamma }_1{I}_{1k}\left(t\right)+{\gamma }_2{I}_{2k}\left(t\right)-\delta R_k\left(t\right)\mathrm{.}\end{array}$ (6)

因此只需在区域 $\Omega =\left\{\left(I_{1k}\left(t\right)\mathrm{,}{I}_{2k}\left(t\right)\mathrm{,}{R}_k\left(t\right)\right)|0\le I_{1k}\left(t\right)\mathrm{,}{I}_{2k}\left(t\right)\mathrm{,}{R}_k\left(t\right)\le \mathrm{1,\; <mi>\; </mi>}k=\mathrm{1,2,}\cdots \mathrm{,}n\right\}$ 上研究系统(6)即可，并且容易验证 $\Omega$ 是关于系统(6)的正不变集。

3. 动力学分析

在这一部分，我们将讨论模型(6)的一些动力学性质。

定理1 设基本再生数

${\mathcal{R}}_0=\frac{\alpha \lambda \langle k{\phi }_1\left(k\right)g\left(k\right)\rangle }{{\gamma }_1\langle kg\left(k\right)\rangle }+\frac{\left(1-\alpha \right)\lambda \langle k{\phi }_2\left(k\right)g\left(k\right)\rangle }{{\gamma }_2\langle kg\left(k\right)\rangle }\mathrm{.}$

考虑系统(6)，我们有如下结论：

(i) 系统(6)总存在一个无病平衡点 $E_0=\left(\underset{3n个}{\underset{︸}{\mathrm{0,0,0,}\cdots \mathrm{,0,0,0}}}\right)$ 。

(ii) 当 ${\mathcal{R}}_0>1$ 时，系统(6)有且仅有唯一的正平衡点 $E^{*}=\left(I_{11}^{*},I_{21}^{*},R_1^{*},\cdots ,I_{1n}^{*},I_{2n}^{*},R_n^{*}\right)$ 。并且 $E^{*}=\left(I_{11}^{*},I_{21}^{*},R_1^{*},\cdots ,I_{1n}^{*},I_{2n}^{*},R_n^{*}\right)$ 满足

$I_{1k}^{*}=\frac{\alpha \delta {\gamma }_2\lambda kg\left(k\right){\theta }^{*}}{D_k},I_{2k}^{*}=\frac{\left(1-\alpha \right)\delta {\gamma }_1\lambda kg\left(k\right){\theta }^{*}}{D_k},R_k^{*}=\frac{{\gamma }_1{\gamma }_2\lambda kg\left(k\right){\theta }^{*}}{D_k},$

其中 $D_k=\lambda kg\left(k\right){\theta }^{*}\left(\alpha \delta {\gamma }_2+\left(1-\alpha \right)\delta {\gamma }_1+{\gamma }_1{\gamma }_2\right)+\delta {\gamma }_1{\gamma }_2\langle kg\left(k\right)\rangle$ 。

证明：显然， $E_0$ 总是系统(6)的一个平衡点。我们接下来研究 $E^{\mathrm{<mo>\; *\; </mo>}}$ 的存在性。注意到平衡点 $E^{\mathrm{<mo>\; *\; </mo>}}$ 应该满足如下的方程

$(\begin{array}{l}\frac{\alpha \lambda k\left(1-I_{1k}^{\mathrm{<mo>\; *\; </mo>}}-I_{2k}^{\mathrm{<mo>\; *\; </mo>}}-R_k^{\mathrm{<mo>\; *\; </mo>}}\right)g\left(k\right)}{\langle kg\left(k\right)\rangle }{\theta }^{\mathrm{<mo>\; *\; </mo>}}-{\gamma }_1{I}_{1k}^{\mathrm{<mo>\; *\; </mo>}}=\mathrm{0,}\\ \frac{\left(1-\alpha \right)\lambda k\left(1-I_{1k}^{\mathrm{<mo>\; *\; </mo>}}-I_{2k}^{\mathrm{<mo>\; *\; </mo>}}-R_k^{\mathrm{<mo>\; *\; </mo>}}\right)g\left(k\right)}{\langle kg\left(k\right)\rangle }{\theta }^{\mathrm{<mo>\; *\; </mo>}}-{\gamma }_2{I}_{2k}^{\mathrm{<mo>\; *\; </mo>}}=\mathrm{0,}\\ {\gamma }_1{I}_{1k}^{\mathrm{<mo>\; *\; </mo>}}+{\gamma }_2{I}_{2k}^{\mathrm{<mo>\; *\; </mo>}}-\delta R_k^{\mathrm{<mo>\; *\; </mo>}}=\mathrm{0,}\end{array}$ (7)

由方程组(7)的第一个和第二个方程，我们有

$I_{2k}^{*}=\frac{{\gamma }_1\left(1-\alpha \right)}{{\gamma }_2\alpha }{I}_{1k}^{*}.$ (8)

由方程组(7)的第三个方程以及等式(8)有

$R_k^{*}=\frac{{\gamma }_1}{\mu +\delta }{I}_{1k}^{*}+\frac{{\gamma }_2}{\mu +\delta }{I}_{2k}^{*}=\frac{{\gamma }_1}{\alpha \delta }{I}_{1k}^{*}.$ (9)

将等式(8)和(9)代入(3)的第一个方程，可得

$I_{1k}^{*}=\frac{\alpha \delta {\gamma }_2\lambda kg\left(k\right){\theta }^{*}}{D_k}.$ (10)

将(8)和(10)式代入 ${\theta }^{\mathrm{<mo>\; *\; </mo>}}$ 的表达式，得到一个关于 ${\theta }^{\mathrm{<mo>\; *\; </mo>}}$ 的自洽方程

${\theta }^{*}={\displaystyle \underset{k=1}{\overset{n}{\sum }}}{\phi }_1\left(k\right)P\left(k\right)\frac{\alpha \delta {\gamma }_2\lambda kg\left(k\right){\theta }^{*}}{D_k}+{\displaystyle \underset{k=1}{\overset{n}{\sum }}}{\phi }_2\left(k\right)P\left(k\right)\frac{\left(1-\alpha \right)\delta {\gamma }_1\lambda kg\left(k\right){\theta }^{*}}{D_k}.$ (11)

令 $f\left({\theta }^{*}\right)={\displaystyle {\sum }_{k=1}^n}{\phi }_1\left(k\right)P\left(k\right)\frac{\alpha \delta {\gamma }_2\lambda kg\left(k\right)}{D_k}+{\displaystyle {\sum }_{k=1}^n}{\phi }_2\left(k\right)P\left(k\right)\frac{\left(1-\alpha \right)\delta {\gamma }_1\lambda kg\left(k\right)}{D_k}$ ，然后方程(11)可写成如下

形式

${\theta }^{*}={\theta }^{*}f\left({\theta }^{*}\right).$ (12)

由于

$(\begin{array}{l}f\left(\theta \right)<\mathrm{1,}\\ \frac{\text{d}f\left({\theta }^{\mathrm{<mo>\; *\; </mo>}}\right)}{\text{d}{\theta }^{\mathrm{<mo>\; *\; </mo>}}}=-{\displaystyle \underset{k=1}{\overset{n}{\sum }}}{\phi }_1\left(k\right)P\left(k\right)\frac{\alpha \delta {\gamma }_2\lambda kg\left(k\right)A_k}{D_k^2}-{\displaystyle \underset{k=1}{\overset{n}{\sum }}}{\phi }_2\left(k\right)P\left(k\right)\frac{\left(1-\alpha \right)\delta {\gamma }_1\lambda kg\left(k\right)A_k}{D_k^2}<0.\end{array}$

其中 $A_k=\lambda kg\left(k\right)\left(\alpha \delta {\gamma }_2+\left(1-\alpha \right)\delta {\gamma }_1+{\gamma }_1{\gamma }_2\right)$ 。若 ${\mathcal{R}}_0>1$ ，则有 $f\left(0\right)>1$ 。此时方程(12)只存在唯一一个非零解 ${\theta }^{\mathrm{<mo>\; *\; </mo>}}$ 。即当 ${\mathcal{R}}_0>1$ 时，系统(6)只存在唯一的地方病平衡点。证毕。

接下来，我们将讨论(6)的无病平衡点 $E_0$ 稳定性。

定理2 若 ${\mathcal{R}}_0<1$ ，系统(6)的无病平衡点 $E_0$ 是全局渐近稳定的。

证明：令

$\begin{array}{c}y\left(t\right)={\left(y_1\left(t\right),y_2\left(t\right),\cdots ,y_n\left(t\right),y_{n+1}\left(t\right),y_{n+2}\left(t\right),\cdots ,y_{2n},y_{2n+1}\left(t\right),y_{2n+2}\left(t\right),\cdots ,y_{3n}\right)}^{\text{T}}\\ ={\left(I_{11}\left(t\right),I_{12}\left(t\right),\cdots ,I_{1n}\left(t\right),I_{21}\left(t\right),I_{22}\left(t\right),\cdots ,I_{2n}\left(t\right),R_1\left(t\right),R_2\left(t\right),\cdots ,R_n\left(t\right)\right)}^{\text{T}},\end{array}$

将系统(6)可写成如下矩阵向量的形式：

$\frac{\text{d}y}{\text{d}t}=Ay+N\left(y\right)\mathrm{,}$ (13)

其中

$A=\left(\begin{array}{ccc}{A}_{11}& A_{12}& O\\ A_{21}& A_{22}& O\\ A_{31}& A_{32}& A_{33}\end{array}\right),$

$A_{11}=\left(\begin{array}{cccc}\alpha g\left(1\right)X_1-{\gamma }_1& \alpha g\left(1\right)X_2& \cdots & \alpha g\left(1\right)X_n\\ 2\alpha g\left(2\right)X_1& 2\alpha g\left(2\right)X_2-{\gamma }_1& \cdots & 2\alpha g\left(2\right)X_n\\ ⋮& ⋮& \ddots & ⋮\\ n\alpha g\left(n\right)X_1& n\alpha g\left(n\right)X_2& \cdots & n\alpha g\left(n\right)X_n-{\gamma }_1\end{array}\right),$

$A_{12}=\left(\begin{array}{cccc}\alpha g\left(1\right)Y_1& \alpha g\left(1\right)Y_2& \cdots & \alpha g\left(1\right)Y_n\\ 2\alpha g\left(2\right)Y_1& 2\alpha g\left(2\right)Y_2& \cdots & 2\alpha g\left(2\right)Y_n\\ ⋮& ⋮& \ddots & ⋮\\ n\alpha g\left(n\right)Y_1& n\alpha g\mathrm{<mo\; stretchy="false">\; (\; </mo>}\left(n\right)Y_2& \cdots & n\alpha g\left(n\right)Y_n\end{array}\right)\mathrm{,}$

$A_{21}=\left(\begin{array}{cccc}\left(1-\alpha \right)g\left(1\right)X_1& \left(1-\alpha \right)g\left(1\right)X_2& \cdots & \left(1-\alpha \right)g\left(1\right)X_n\\ 2\left(1-\alpha \right)g\left(2\right)X_1& 2\left(1-\alpha \right)g\left(2\right)X_2& \cdots & 2\left(1-\alpha \right)g\left(2\right)X_n\\ ⋮& ⋮& \ddots & ⋮\\ n\left(1-\alpha \right)g\left(n\right)X_1& n\left(1-\alpha \right)g\left(n\right)X_2& \cdots & n\left(1-\alpha \right)g\left(n\right)X_n\end{array}\right),$

$A_{22}=\left(\begin{array}{cccc}\left(1-\alpha \right)g\left(1\right)Y_1-{\gamma }_2& \left(1-\alpha \right)g\left(1\right)Y_2& \cdots & \left(1-\alpha \right)g\left(1\right)Y_n\\ 2\left(1-\alpha \right)g\left(2\right)Y_1& 2\left(1-\alpha \right)g\left(2\right)Y_2-{\gamma }_2& \cdots & 2\left(1-\alpha \right)g\left(2\right)Y_n\\ ⋮& ⋮& \ddots & ⋮\\ n\left(1-\alpha \right)g\left(n\right)Y_1& n\left(1-\alpha \right)g\left(n\right)Y_2& \cdots & n\left(1-\alpha \right)g\left(n\right)Y_n-{\gamma }_2\end{array}\right),$

$A_{31}={\gamma }_{1}I,\mathrm{<mi>\; </mi>}{A}_{32}={\gamma }_{2}I,\mathrm{<mi>\; </mi>}{A}_{33}=-\delta I$ 。 $X_k=\frac{\lambda {\phi }_1\left(k\right)P\left(k\right)}{\langle kg\left(k\right)\rangle },\mathrm{<mi>\; </mi>}{Y}_k=\frac{\lambda {\phi }_2\left(k\right)P\left(k\right)}{\langle kg\left(k\right)\rangle }$ 。其中O是n阶零矩阵，I是n

阶单位矩阵，而

$N\left(y\right)=-\left(\begin{array}{c}\frac{\alpha \lambda g\left(1\right)}{\langle kg\left(k\right)\rangle }\left(y_1+y_{n+1}+y_{2n+1}\right)\theta \\ ⋮\\ \frac{\alpha \lambda ng\left(n\right)}{\langle kg\left(k\right)\rangle }\left(y_n+y_{2n}+y_{3n}\right)\theta \\ \frac{\left(1-\alpha \right)\lambda g\left(1\right)}{\langle kg\left(k\right)\rangle }\left(y_1+y_{n+1}+y_{2n+1}\right)\theta \\ ⋮\\ \frac{\left(1-\alpha \right)\lambda ng\left(n\right)}{\langle kg\left(k\right)\rangle }\left(y_n+y_{2n}+y_{3n}\right)\theta \\ 0\\ ⋮\\ 0\end{array}\right).$

由于向量 $N\left(y\right)$ 每一项都是非正的，因此有 $\frac{\text{d}y}{\text{d}t}\le Ay$ ，故可得如下比较系统

$\frac{\text{d}y}{\text{d}t}=Ay\mathrm{.}$ (14)

系统(14)只有一个平衡点 $E_0\left(\mathrm{0,0,}\cdots \mathrm{,0}\right)$ ，其雅可比矩阵对应的特征方程为

${\left(x+\delta \right)}^n{\left(x+{\gamma }_1\right)}^{n-1}{\left(x+{\gamma }_2\right)}^{n-1}\left(x^2+px+q\right)=0.$

其中 $p={\gamma }_1+{\gamma }_2-\frac{\alpha \lambda \langle k{\phi }_1\left(k\right)g\left(k\right)\rangle }{\langle kg\left(k\right)\rangle }-\frac{\left(1-\alpha \right)\lambda \langle k{\phi }_2\left(k\right)g\left(k\right)\rangle }{\langle kg\left(k\right)\rangle }$ ，

$q={\gamma }_1{\gamma }_2-\frac{\alpha \lambda {\gamma }_2\langle k{\phi }_1\left(k\right)g\left(k\right)\rangle }{\langle kg\left(k\right)\rangle }-\frac{\left(1-\alpha \right)\lambda {\gamma }_1\langle k{\phi }_2\left(k\right)g\left(k\right)\rangle }{\langle kg\left(k\right)\rangle }$ 。显然 ${\mathcal{R}}_0<1\iff q>0$ ，并且 ${\mathcal{R}}_0<1$ 时， $p>0$ ，

此时A的所有特征值实部都为负的，故系统(14)的平衡点 $E_0$ 是全局渐近稳定的。又因为矩阵A的所有非对角线上元素是非负的，因此系统(14)是拟单调增的 [17] ，故由比较原理和解的正性，定理2得证。

下面，我们来讨论系统(6)的持久性。首先，我们给出如下定义

定义 若对于系统(6)的满足初值条件 $I_{1k}\left(0\right)+I_{2k}\left(0\right)+R_k\left(0\right)>0$ 的任意解，存在一个常数 $\sigma >0$ 使得

$\underset{t\to \infty }{\lim \inf }\sqrt{\left(I_{11}^2\left(t\right)+\cdots +I_{1n}^2\left(t\right)+I_{21}^2\left(t\right)+\cdots +I_{2n}^2\left(t\right)+R_1^2\left(t\right)+\cdots +I_n^2\left(t\right)\right)}\ge \sigma .$

则称该系统是准持久的。

然后，再引入如下引理

引理3 [4] 考虑如下n维的系统

$\frac{\text{d}y}{\text{d}t}=Ay+N\left(y\right)\mathrm{,}$ (15)

其中，向量 $y={\left(y_1,y_2,\cdots ,y_n\right)}^{\text{T}}$ ，A是一个 $n\times n$ 阶方阵并且向量 $N\left(y\right)$ 在区域 $D\subseteq R^n$ 内是连续可微的。假设

(i)：紧的凸集 $C\subseteq D$ 是(15)的正向不变集，并且 $0\in C$ ；

(ii)：存在某个正整数 $m\le n$ ，使得 $\underset{y\to 0}{\lim }\frac{\Vert N\left(y\right)\Vert }{{\displaystyle {\sum }_{i=1}^m}{y}_i^2}=0$ ；

(iii)：存在 $r>0$ 和 $A^{\text{T}}$ 的一个正的特征值对应的一个实特征向量 $\omega$ 使得对任意的 $y\in C$ ，都有 $\left(\omega \cdot y\right)\ge r\sqrt{{\displaystyle {\sum }_{i=1}^m}{y}_i^2}$ ；

(iv)：对所有的 $y\in C$ ，都有 $\left(\omega \cdot y\right)\le 0$ 。然后，系统(15)的任何一个初值为 $y\left(0\right)\in C-\left\{0\right\}$ 的解 $y\left(t\mathrm{<mo>\; ;\; </mo>}y\left(0\right)\right)$ 满足

$\underset{t\to \infty }{\lim \inf }\Vert y\left(t;y\left(0\right)\right)\Vert \ge \zeta ,$

其中 $\zeta >0$ 不依赖于 $y\left(0\right)$ ，另外(15) 存在一个常数解 $y=y^{*},\mathrm{<mi>\; </mi>}{y}^{*}\in C-\left\{0\right\}$ 。

因此我们可得到如下定理

定理4 若 ${\mathcal{R}}_0>1$ ，则系统(6)是准持久的。

证明：下面我们通过引理3来证明此定理

1)：紧的凸集 $\Omega$ 是(6)的正向不变集，并且 $0\in \Omega$ ，故条件(i)成立；

2)：取正整数 $m=2n$ ，易得 $\underset{y\to 0}{\lim }\frac{\Vert N\left(y\right)\Vert }{{\displaystyle {\sum }_{i=1}^m}{y}_i^2}=0$ ，因而条件(ii)成立；

3)：记 $B=\left(\begin{array}{cc}{A}_{11}& A_{12}\\ A_{21}& A_{22}\end{array}\right)$ ，则 $B^{\text{T}}$ 是一个 $2n\times 2n$ 不可约的矩阵，且当 $i\ne j$ 时， $a_{ij}\ne 0$ 。由文献 [18] 的引理4， $B^{\text{T}}$ 存在一个正的特征向量 $\bar{\omega }=\left({\vartheta }_1,{\vartheta }_2,\cdots ,{\vartheta }_{2n}\right)$ ，对应的特征值为 $S\left(B^{\text{T}}\right)=\max Re{\lambda }_i\mathrm{,\; <mi>\; </mi>}i=\mathrm{1,2,}\cdots \mathrm{,2}n$ 。令 ${\vartheta }_{2n+1}={\vartheta }_{2n+2}=\cdots ={\vartheta }_{3n}=0$ ，取 $\omega =\left({\vartheta }_1,{\vartheta }_2,\cdots ,{\vartheta }_{3n}\right)$ ，因此有 $A^{\text{T}}\omega =S\left(B^{\text{T}}\right)\omega$ ，即 $\omega$ 是矩阵 $A^{\text{T}}$ 的特征向量，对应着特征值 $S\left(B^{\text{T}}\right)$ 。记 $r=\underset{1\le k\le 2n}{\min }{\vartheta }_k$ ，故对所有的 $y\in C$ ，都有 $\left(\omega \cdot y\right)\ge r\sqrt{{\displaystyle {\sum }_{k=1}^{2n}}{y}_k^2}$ ，因此条件(iii)成立。

4)：由于 $N\left(y\right)$ 的每一项都小于等于0，而 $\omega >0$ ，故对所有的 $y\in C$ ，都有 $\left(\omega \cdot y\right)\le 0$ 。即条件(iv)也成立。因此，由引理3，存在常数 $\zeta >0$ 使得

$\underset{t\to \infty }{\lim \inf }\Vert y\left(t;y\left(0\right)\right)\Vert \ge \zeta ,$

即有

$\underset{t\to \infty }{\lim \inf }\sqrt{\left(I_{11}^2\left(t\right)+\cdots +I_{1n}^2\left(t\right)+I_{21}^2\left(t\right)+\cdots +I_{2n}^2\left(t\right)+R_1^2\left(t\right)+\cdots +I_n^2\left(t\right)\right)}\ge \zeta \mathrm{.}$

因此，定理4得证。

4. 数值模拟

这一节，我们通过两个数值仿真例子来验证相关结果的正确性。首先，根据优先连接算法我们构建了一个带有100节点的BA无标度网络。接下来我们验证 $E_0$ 的全局稳定渐近稳定性。

例1. 这里我们取 $g\left(k\right)=k^a$ ， ${\phi }_1\left(k\right)=k^b$ 以及 ${\phi }_2\left(k\right)=k^c$ 其它参数分别为 $\alpha =0.3$ ， $\lambda =0.04$ ， ${\gamma }_1=0.09$ ， ${\gamma }_2=0.08$ ， $\delta =0.05$ ， $a=1$ ， $b=0.6$ 和 $c=0.9$ 。那么， ${\mathcal{R}}_0=0.4989<1$ ，由定理2可知，如果 ${\mathcal{R}}_0=0.4989<1$ ，则无病平衡点 $E_0$ 是全局渐近稳定的。给定初值条件 $S_k\left(0\right)=0.5$ ， $I_{1k}\left(0\right)=0.15$ ， $I_{2k}\left(0\right)=0.2$ 和 $R_k\left(0\right)=0.15$ ，我们可以得到 $S_k\left(t\right)$ ， $I_{1k}\left(t\right)$ ， $I_{2k}\left(t\right)$ 以及 $R_k\left(t\right)$ 的随时间的演化图，如图1所示。从图中可以看出疾病最终都会消亡，与我们的结果相吻合。

Figure 1. The time evolution of the densities of $I_{1k}\left(t\right)$, $I_{2k}\left(t\right)$ and $R_k(t)$

图1. $I_{1k}\left(t\right)$ 、 $I_{2k}\left(t\right)$ 以及 $R_k\left(t\right)$ 的时间演化图

接下来，我们验证定理1和定理4。

例2. 这里我们分别取 ${\gamma }_1=0.2$ ， ${\gamma }_2=0.15$ 和 $\delta =0.1$ ，其它参数及初值条件分别和例1中的相同。则有 ${\mathcal{R}}_0=2.2006>1$ ，由前面的研究结果可知，如果 ${\mathcal{R}}_0>1$ ，则系统(6)存在唯一的地方病平衡点 $E^{\mathrm{<mo>\; *\; </mo>}}$ ，并且系统(6)是持久存在的 $S_k\left(t\right)$ ， $I_{1k}\left(t\right)$ ， $I_{2k}\left(t\right)$ 以及 $R_k\left(t\right)$ 的随时间的演化图，如图2所示。从图中可以看出系统最终是持久存在的，这与我们的结果是一致的。

Figure 2. The time evolution of the densities of $I_{1k}\left(t\right)$, $I_{2k}\left(t\right)$ and $R_k(t)$

图2. $I_{1k}\left(t\right)$ 、 $I_{2k}\left(t\right)$ 以及 $R_k\left(t\right)$ 的时间演化图

5. 结论

为了获得个体接触的频繁程度对疾病传播的影响，我们建立并分析了一个考虑权重的SIRS网络传染病模型，得到了基本再生数，给出了地方病平衡点的存在唯一性条件以及证明了无病平衡点的全局渐近稳定性。同时还研究了系统的持久性并通过数值模拟验证了我们的结果。由前面的研究结果

${\mathcal{R}}_0=\frac{\alpha \lambda \langle k{\phi }_1\left(k\right)g\left(k\right)\rangle }{{\gamma }_1\langle kg\left(k\right)\rangle }+\frac{\left(1-\alpha \right)\lambda \langle k{\phi }_2\left(k\right)g\left(k\right)\rangle }{{\gamma }_2\langle kg\left(k\right)\rangle }$ 可知感染者中低风险感染者所占的比率 $\alpha$ 能影响疾病传

播阈值并且还发现 ${\mathcal{R}}_0$ 随着治愈率 ${\gamma }_1$ 和 ${\gamma }_2$ 的增大而减小。因此提高疾病的治愈率是抑制疾病传播的重要措施。同时我们可以看出 ${\mathcal{R}}_0$ 对于网络拓扑结构是敏感的，网络的度和权重的改变都能导致基本再生数 ${\mathcal{R}}_0$ 的变化，进而影响疾病的传播。

基金资助

受湖南省大学生研究性学习与创新性实验计划项目(湘教通[2016] 283号(序号207))和长沙理工大学“十三五”专业综合改革试点项目(数学与应用数学专业)资助。

NOTES

^*通讯作者。