带有时滞的SEAIQRS网络传染病模型动力学分析

Dynamic Analysis of SEAIQRS Network Epidemic Model with Time Delay

摘要

大多数经典的传染病模型并没有考虑到人口的异质性， 为了更准确地描述新型冠状病毒的传播机制， 考虑了基于无标度网络的带有时滞的SEAIQRS传染病模型。通过下一代矩阵法求得其基本再生数， 并证明了正平衡点的存在唯一性， 然后证明了模型在无病平衡点处Hopf分岔的存在性， 最后利用MATLAB模拟了本文模型， 通过在节点数为1 000的无标度网络上运行30次取平均值得到关于各类人口的时间演化曲线， 并利用PRCC敏感性分析研究了敏感性更强的参数对疾病传播的影响。结果表明， 当R₀ >1时， 疾病存在唯一正平衡点且是局部渐近稳定的， 并且发现染病率和治愈率是影响疾病传播的关键参数。文中考虑了无症状染病者变为有症状染病者的时间滞后\tau， 这可能会导致传染病在传播过程中出现周期性传播的现象。为了达到控制传染病暴发的目的， 通过数值模拟绘制了平面投影图来表明R_{\mathrm{0}}， {\beta }_{\mathrm{1}}和{\gamma }_{\mathrm{1}}三者之间的关系。

Abstract

In most classical epidemic models， the spatial heterogeneity of population is not considered. To describe the transmission mechanism of corona virus more appropriately， an SEAIQRS epidemic model with time delay based on scale-free network was established in this paper. The next-generation matrix method was used to calculate the basic reproduction number. Next， the uniqueness and the existence of the endemic equilibrium were demonstrated， and the existence of the Hopf bifurcation at the disease-free equilibrium was also proved. Finally， MATLAB was used to simulate the model of this paper. A scale-free network with 1 000 nodes was averaged thirty times to obtain time evolution curves about each type of population， and the PRCC sensitivity analysis was performed to analyze the impact of more sensitive parameters on disease transmission. It is shown that when R_{\mathrm{0}}>1， the disease has a unique endemic equilibrium， which is locally asymptotically stable， and the infection rate and cure rate are the key parameters affecting the spread of the disease. The time lag \tau from asymptomatic infected people to symptomatic infected people is considered in this paper， which may lead to periodic transmission of diseases. Moreover， projection diagram planes are presented in the numerical simulation section to reflect the relationship between R_{\mathrm{0}}， {\beta }_{\mathrm{1}} and {\gamma }_{\mathrm{1}} in order to achieve the purpose of control the outbreaks.

Graphical abstract

关键词

传染病模型 / 稳定定性分析 / 无标度网络 / 时滞 / Hopf分岔

Key words

epidemic model / stability analysis / scale-free network / time delay / Hopf bifurcation

0 引　言

传染病不仅会对个体健康造成危害， 而且会给全球的公共卫生和社会经济发展带来巨大的挑战^［1］。因此， 为了社会秩序的正常运行和经济的高速发展， 建立传染病传播的数学模型已经成为研究人员的重要任务之一。在传染病传播的数学模型中， 网络科学已经成为一个热门的研究领域^［2-3］。网络科学的基本思想是将一个复杂的系统看作是由多个相互作用的单元构成的网络， 通过分析网络的结构和特征来更好地理解和模拟系统的行为和演化。在传染病传播的数学模型中， 网络结构可以用来描述人与人之间的接触和传播方式。传染病模型可以分为随机网络模型和确定性网络模型两种类型， 其中， 随机网络模型通常采用均匀随机网络、 ER网络等， 而确定性网络模型通常采用规则网络、 小世界网络等^［4‐8］。然而， 实际中的传染病传播往往不是完全随机的， 而是具有一定的规律性和结构性。无标度网络是一种具有高聚集度和短路径长度的网络结构， 能够更好地模拟真实世界中的复杂网络。许多传染病的传播模式都具有无标度网络的特征， 因此， 使用无标度网络模型对传染病的传播进行建模具有重要的意义^［9‐12］。

本文将基于无标度网络构建一种新的SEAIQRS时滞传染病模型， 该模型考虑了潜伏期、 染病态（无症状与有症状）和恢复态的传染病传播过程。具体是将网络节点视为人群中的个体， 用节点之间的连边表示两个个体之间的接触， 然后模拟传染病在节点之间的传播过程。本文模型将引入时滞效应， 以更好地模拟实际中传染病的传播过程； 然后证明了模型在无病平衡点处Hopf分岔的存在性； 最后利用MATLAB对模型进行模拟， 分析了模型参数对疾病的影响。研究结果将有助于我们更好地理解和预测传染病传播的行为和演化， 为公共卫生和政策制定提供科学依据。

1 模型构造

假设总人群的接触网络基于无标度网络， 每个节点都代表一个人， 网络上的连边表示人与人之间的接触， 假设：

1） 网络上的人分为6类： 易感者S（Susceptible）， 潜伏者E（Exposed）， 无症状染病者A（Asymptomatic infected）， 有症状染病者I（Symptomatic infected）， 隔离者Q（Quarantine）， 恢复者R（Recovered）， 个体总数为N。S_k(t)， E_k(t)， A_k(t)， I_k(t)， Q_k(t)， R_k(t)分别表示在t时刻度为k（k=\mathrm{1}，\cdots ，n）的易感者、 潜伏者、 无症状染病者、 有症状染病者、 隔离者、 恢复者的密度。

2） 易感个体S通过接触染病者（包括无症状染病者与有症状染病者）变成潜伏者E， 潜伏者E又以概率\alpha变成染病者， 无症状染病者A以概率q变为有症状染病者I， 以概率{\gamma }_{\mathrm{1}}恢复健康， 变成恢复者R， 有症状染病者A以概率m被隔离， 成为隔离者Q， 以概率{\gamma }_{\mathrm{2}}恢复健康， 成为恢复者R， 隔离者Q以概率{\gamma }_{\mathrm{3}}恢复健康， 成为恢复者R， 恢复者R又以概率\delta重新变为易感者S。假设无症状染病者变为有症状染病者之后不能立刻被发现， 而是要经过时间\tau才能被发现。

基于以上假设， 本文考虑了基于无标度网络的时滞传染病模型

\left\{\begin{array}{l}\begin{array}{l}\frac{\mathrm{d}{S}_k(t)}{\mathrm{d}t}=-{\beta }_{\mathrm{1}}k{S}_k(t){\theta }_{\mathrm{1}}(t)-{\beta }_{\mathrm{2}}k{S}_k(t){\theta }_{\mathrm{2}}(t)+\\ \delta R_k(t),\end{array}\\ \begin{array}{l}\frac{\mathrm{d}{E}_k(t)}{\mathrm{d}t}={\beta }_{\mathrm{1}}k{S}_k(t){\theta }_{\mathrm{1}}(t)+{\beta }_{\mathrm{2}}k{S}_k(t){\theta }_{\mathrm{2}}(t)-\\ \alpha E_k(t),\end{array}\\ \frac{\mathrm{d}{A}_k(t)}{\mathrm{d}t}=p\alpha E_k(t)-q{A}_k(t-\tau )-{\gamma }_{\mathrm{1}}{A}_k(t),\\ \begin{array}{l}\frac{\mathrm{d}{I}_k(t)}{\mathrm{d}t}=(\mathrm{1}-p)\alpha E_k(t)+q{A}_k(t-\tau )-\\ m{I}_k(t)-{\gamma }_{\mathrm{2}}{I}_k(t),\end{array}\\ \frac{\mathrm{d}{Q}_k(t)}{\mathrm{d}t}=m{I}_k(t)-{\gamma }_{\mathrm{3}}{Q}_k(t),\\ \begin{array}{l}\frac{\mathrm{d}{R}_k(t)}{\mathrm{d}t}={\gamma }_{\mathrm{1}}{A}_k(t)+{\gamma }_{\mathrm{2}}{I}_k(t)+{\gamma }_{\mathrm{3}}{Q}_k(t)-\\ \delta R_k(t)\mathrm{。}\end{array}\end{array}\right.（1）

模型（1）的所有参数{\beta }_{\mathrm{1}}， {\beta }_{\mathrm{2}}， \delta， \alpha， p， q， m， {\gamma }_{\mathrm{1}}， {\gamma }_{\mathrm{2}}， {\gamma }_{\mathrm{3}}都是非负的， 并且都满足0\le{\beta }_{\mathrm{1}}， {\beta }_{\mathrm{2}}， \delta， \alpha， p， q， m， {\gamma }_{\mathrm{1}}， {\gamma }_{\mathrm{2}}， {\gamma }_{\mathrm{3}}\le1； k表示每个节点的度； {\theta }_{\mathrm{1}}(t)和{\theta }_{\mathrm{2}}(t)分别表示度为k的易感者接触到无症状染病者和有症状染病者的概率， 并且

{\theta }_{\mathrm{1}}(t)={\displaystyle \sum _{k=\mathrm{1}}^n}\frac{kp(k)A_k(t)}{〈k〉}，

{\theta }_{\mathrm{2}}(t)={\displaystyle \sum _{k=\mathrm{1}}^n}\frac{kp(k)I_k(t)}{〈k〉},

式中： 〈k〉={\displaystyle \sum _{k=\mathrm{1}}^n}kp(k)表示无标度网络的平均度； p（k）表示网络的度分布。

2 动力学分析

通过之前对模型的一些说明， 易知S_k(t)+E_k(t)+A_k(t)+I_k(t)+Q_k(t)+R_k(t)=\mathrm{1}， 并且集合\Omega=\left\{S_k(t)\right.>0， E_k(t)≥0， A_k(t)≥0， I_k(t)≥0， Q_k(t)≥0， R_k(t)≥0， S_k(t)+E_k(t)+A_k(t)+I_k(t)+Q_k(t)+\left.R_k(t)=\mathrm{1}\right\}是模型（1）的正不变集。

下面计算模型（1）的基本再生数R₀。显然， 模型（1）具有无病平衡E_{\mathrm{0}}=(\underset{n}{\underbrace{\mathrm{1},\cdots ,\mathrm{1}}},\underset{\mathrm{5}n}{\underbrace{\mathrm{0},\cdots ,\mathrm{0}}})， 根据下一代矩阵法， 结合模型（1）中的E_k， A_k， I_k， Q_k仓室， 可以得到 F （新增感染率）和 V （仓室的移出率）

\mathit{F}={\left[\begin{array}{cccc}\mathrm{0}& {\beta }_{\mathrm{1}}\mathrm{*}P& {\beta }_{\mathrm{2}}\mathrm{*}P& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}\end{array}\right]}_{\mathrm{4}n\times \mathrm{4}n},

\mathit{V}={\left[\begin{array}{cccc}\alpha {\mathit{I}}_n& \mathrm{0}& \mathrm{0}& \mathrm{0}\\ -p\alpha {\mathit{I}}_n& (q+{\gamma }_{\mathrm{1}}){\mathit{I}}_n& \mathrm{0}& \mathrm{0}\\ -(\mathrm{1}-p)\alpha {\mathit{I}}_n& -q{\mathit{I}}_n& (m+{\gamma }_{\mathrm{2}}){\mathit{I}}_n& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& -m{\mathit{I}}_n& {\gamma }_{\mathrm{3}}{\mathit{I}}_n\end{array}\right]}_{\mathrm{4}n\times \mathrm{4}n},

其中， {\mathit{I}}_n为n阶单位矩阵。

\mathit{P}=\frac{\mathrm{1}}{〈k〉}{\left[\begin{array}{cccc}{k}_{\mathrm{1}}^{\mathrm{2}}p(k_{\mathrm{1}})& k_{\mathrm{1}}^{}{k}_{\mathrm{2}}p(k_{\mathrm{2}})& \cdots & k_{\mathrm{1}}^{}{k}_{n}p(k_n)\\ k_{\mathrm{2}}{k}_{\mathrm{1}}p(k_{\mathrm{1}})& k_{\mathrm{2}}^{\mathrm{2}}p(k_{\mathrm{2}})& \cdots & k_{\mathrm{2}}{k}_{n}p(k_n)\\ ⋮& ⋮& \ddots & ⋮\\ k_n{k}_{\mathrm{1}}p(k_{\mathrm{1}})& k_n{k}_{\mathrm{2}}p(k_{\mathrm{2}})& \cdots & k_n^{\mathrm{2}}p(k_n)\end{array}\right]}_{n\times n},

\mathit{F}{\mathit{V}}^{-\mathrm{1}}={\left[\begin{array}{cccc}\mathrm{0}& {\beta }_{\mathrm{1}}\mathrm{*}P& {\beta }_{\mathrm{2}}\mathrm{*}P& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}\end{array}\right]}_{\mathrm{4}n\times \mathrm{4}n}\cdot

{\left[\begin{array}{cccc}\frac{\mathrm{1}}{\alpha }{\mathit{I}}_n& \mathrm{0}& \mathrm{0}& \mathrm{0}\\ \frac{p}{q+{\gamma }_{\mathrm{1}}}{\mathit{I}}_n& \frac{\mathrm{1}}{q+{\gamma }_{\mathrm{1}}}{\mathit{I}}_n& \mathrm{0}& \mathrm{0}\\ \frac{(\mathrm{1}-p){\gamma }_{\mathrm{1}}+q}{(q+{\gamma }_{\mathrm{1}})(m+{\gamma }_{\mathrm{2}})}{\mathit{I}}_n& \frac{q}{(q+{\gamma }_{\mathrm{1}})(m+{\gamma }_{\mathrm{2}})}{\mathit{I}}_n& \frac{\mathrm{1}}{m+{\gamma }_{\mathrm{2}}}{\mathit{I}}_n& \mathrm{0}\\ \frac{m((\mathrm{1}-p){\gamma }_{\mathrm{1}}+q)}{(q+{\gamma }_{\mathrm{1}})(m+{\gamma }_{\mathrm{2}}){\gamma }_{\mathrm{3}}}{\mathit{I}}_n& \frac{mq}{(q+{\gamma }_{\mathrm{1}})(m+{\gamma }_{\mathrm{2}}){\gamma }_{\mathrm{3}}}{\mathit{I}}_n& \frac{m}{(m+{\gamma }_{\mathrm{2}}){\gamma }_{\mathrm{3}}}{\mathit{I}}_n& \frac{\mathrm{1}}{{\gamma }_{\mathrm{3}}}{\mathit{I}}_n\end{array}\right]}_{\mathrm{4}n\times \mathrm{4}n}=

{\left[\begin{array}{cccc}\left(\frac{{\beta }_{\mathrm{1}}p}{q+{\gamma }_{\mathrm{1}}}+\frac{{\beta }_{\mathrm{2}}((\mathrm{1}-P){\gamma }_{\mathrm{1}}+q)}{(q+{\gamma }_{\mathrm{1}})(m+{\gamma }_{\mathrm{2}})}\right)P& (\frac{{\beta }_{\mathrm{1}}}{q+{\gamma }_{\mathrm{1}}}+\frac{{\beta }_{\mathrm{2}}q}{(q+{\gamma }_{\mathrm{1}})(m+{\gamma }_{\mathrm{2}})})P& \frac{{\beta }_{\mathrm{2}}}{m+{\gamma }_{\mathrm{2}}}P& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}\end{array}\right]}_{\mathrm{4}n\times \mathrm{4}n}\mathrm{。}

矩阵的谱半径为

\rho (\mathit{F}{\mathit{V}}^{-\mathrm{1}})=

\mathrm{t}\mathrm{r}\mathrm{a}\mathrm{c}\mathrm{e}\left(\left(\frac{{\beta }_{\mathrm{1}}p}{q+{\gamma }_{\mathrm{1}}}+\frac{{\beta }_{\mathrm{2}}((\mathrm{1}-p){\gamma }_{\mathrm{1}}+q)}{(q+{\gamma }_{\mathrm{1}})(m+{\gamma }_{\mathrm{2}})}\right)P\right)=

\left(\frac{{\beta }_{\mathrm{1}}p}{q+{\gamma }_{\mathrm{1}}}+\frac{{\beta }_{\mathrm{2}}((\mathrm{1}-p){\gamma }_{\mathrm{1}}+q)}{(q+{\gamma }_{\mathrm{1}})(m+{\gamma }_{\mathrm{2}})}\right)\mathrm{t}\mathrm{r}\mathrm{a}\mathrm{c}\mathrm{e}P=

\left(\frac{{\beta }_{\mathrm{1}}p}{q+{\gamma }_{\mathrm{1}}}+\frac{{\beta }_{\mathrm{2}}((\mathrm{1}-p){\gamma }_{\mathrm{1}}+q)}{(q+{\gamma }_{\mathrm{1}})(m+{\gamma }_{\mathrm{2}})}\right)\frac{〈k^{\mathrm{2}}〉}{〈k〉}\mathrm{。}

因此， 可以得到基本再生数

R_{\mathrm{0}}=\left(\frac{{\beta }_{\mathrm{1}}p}{q+{\gamma }_{\mathrm{1}}}+\frac{{\beta }_{\mathrm{2}}((\mathrm{1}-p){\gamma }_{\mathrm{1}}+q)}{(q+{\gamma }_{\mathrm{1}})(m+{\gamma }_{\mathrm{2}})}\right)\frac{〈k^{\mathrm{2}}〉}{〈k〉}\mathrm{。}

定理 1　当且仅当R₀ > 1， 模型（1）存在唯一的正平衡点。

证明　令\theta ={\beta }_{\mathrm{1}}{\theta }_{\mathrm{1}}+{\beta }_{\mathrm{2}}{\theta }_{\mathrm{2}}， 根据模型（1）， 可以计算出

E_k^{\mathrm{*}}=\frac{k({\beta }_{\mathrm{1}}{\theta }_{\mathrm{1}}+{\beta }_{\mathrm{2}}{\theta }_{\mathrm{2}})}{\alpha }{S}_k^{\mathrm{*}},（2）

{A_k}^{\mathrm{*}}=\frac{pk({\beta }_{\mathrm{1}}{\theta }_{\mathrm{1}}+{\beta }_{\mathrm{2}}{\theta }_{\mathrm{2}})}{q+{\gamma }_{\mathrm{1}}}{S}_k^{\mathrm{*}},（3）

I_k^{\mathrm{*}}=\frac{k((\mathrm{1}-p){\gamma }_{\mathrm{1}}+q)({\beta }_{\mathrm{1}}{\theta }_{\mathrm{1}}+{\beta }_{\mathrm{2}}{\theta }_{\mathrm{2}})}{(q+{\gamma }_{\mathrm{1}})(m+{\gamma }_{\mathrm{2}})}{S}_k^{\mathrm{*}},（4）

Q_k^{\mathrm{*}}=\frac{mk((\mathrm{1}-p){\gamma }_{\mathrm{1}}+q)({\beta }_{\mathrm{1}}{\theta }_{\mathrm{1}}+{\beta }_{\mathrm{2}}{\theta }_{\mathrm{2}})}{(q+{\gamma }_{\mathrm{1}})(m+{\gamma }_{\mathrm{2}}){\gamma }_{\mathrm{3}}}{S}_k^{\mathrm{*}},（5）

R_k^{\mathrm{*}}=\frac{k({\beta }_{\mathrm{1}}{\theta }_{\mathrm{1}}+{\beta }_{\mathrm{2}}{\theta }_{\mathrm{2}})}{\delta }{S}_k^{\mathrm{*}}\mathrm{。}（6）

可求得

S_k^{\mathrm{*}}=\frac{A_{\mathrm{1}}}{A_{\mathrm{1}}+k{A}_{\mathrm{2}}},（7）

其中， A_{\mathrm{1}}=\alpha \delta (q+{\gamma }_{\mathrm{1}})(m+{\gamma }_{\mathrm{2}}){\gamma }_{\mathrm{3}}，A_{\mathrm{2}}=p\alpha \delta (m+{\gamma }_{\mathrm{2}}){\gamma }_{\mathrm{3}}+(\delta +\alpha )(q+{\gamma }_{\mathrm{1}})(m+{\gamma }_{\mathrm{2}}){\gamma }_{\mathrm{3}}-(p{\gamma }_{\mathrm{1}}-q-{\gamma }_{\mathrm{1}})({\gamma }_{\mathrm{3}}+m)\alpha \delta \mathrm{。}

将式（3）和式（4）代入\theta ={\beta }_{\mathrm{1}}{\theta }_{\mathrm{1}}+{\beta }_{\mathrm{2}}{\theta }_{\mathrm{2}}中可得

\theta =\frac{A_{\mathrm{1}}\theta }{A_{\mathrm{1}}+k\theta A_{\mathrm{2}}}\left(\frac{{\beta }_{\mathrm{1}}p}{q+{\gamma }_{\mathrm{1}}}+\frac{{\beta }_{\mathrm{2}}((\mathrm{1}-p){\gamma }_{\mathrm{1}}+q)}{(q+{\gamma }_{\mathrm{1}})(m+{\gamma }_{\mathrm{2}})}\right)\frac{〈k^{\mathrm{2}}〉}{〈k〉}\mathrm{。}

令

\begin{array}{c}F(\theta )=\theta -\frac{A_{\mathrm{1}}\theta }{A_{\mathrm{1}}+k\theta A_{\mathrm{2}}}\left(\frac{{\beta }_{\mathrm{1}}p}{q+{\gamma }_{\mathrm{1}}}+\frac{{\beta }_{\mathrm{2}}((\mathrm{1}-p){\gamma }_{\mathrm{1}}+q)}{(q+{\gamma }_{\mathrm{1}})(m+{\gamma }_{\mathrm{2}})}\right)\frac{〈k^{\mathrm{2}}〉}{〈k〉}=\\ \theta \left(\mathrm{1}-\frac{A_{\mathrm{1}}}{A_{\mathrm{1}}+k\theta A_{\mathrm{2}}}\left(\frac{{\beta }_{\mathrm{1}}p}{q+{\gamma }_{\mathrm{1}}}+\frac{{\beta }_{\mathrm{2}}((\mathrm{1}-p){\gamma }_{\mathrm{1}}+q)}{(q+{\gamma }_{\mathrm{1}})(m+{\gamma }_{\mathrm{2}})}\right)\frac{〈k^{\mathrm{2}}〉}{〈k〉}\right)=\theta f(\theta )\mathrm{。}\end{array}

可知， 0是F(\theta )的一个零解。接下来证F(\theta )有一个正解。当R_{\mathrm{0}}>\mathrm{1}时，

\frac{\mathrm{d}f(\theta )}{\mathrm{d}\theta }=\frac{A_{\mathrm{1}}{A}_{\mathrm{2}}k}{(A_{\mathrm{1}}+k\theta A_{\mathrm{2}}{)}^{\mathrm{2}}}\left(\frac{{\beta }_{\mathrm{1}}p}{q+{\gamma }_{\mathrm{1}}}+\frac{{\beta }_{\mathrm{2}}((\mathrm{1}-p){\gamma }_{\mathrm{1}}+q)}{(q+{\gamma }_{\mathrm{1}})(m+{\gamma }_{\mathrm{2}})}\right)\frac{〈k^{\mathrm{2}}〉}{〈k〉}>\mathrm{0},\begin{array}{l}f(\theta )=\mathrm{1}-\frac{A_{\mathrm{1}}}{A_{\mathrm{1}}+k\theta A_{\mathrm{2}}}\left(\frac{{\beta }_{\mathrm{1}}p}{q+{\gamma }_{\mathrm{1}}}+\frac{{\beta }_{\mathrm{2}}((\mathrm{1}-p){\gamma }_{\mathrm{1}}+q)}{(q+{\gamma }_{\mathrm{1}})(m+{\gamma }_{\mathrm{2}})}\right)\frac{〈k^{\mathrm{2}}〉}{〈k〉}>\mathrm{1}-\\ \frac{A_{\mathrm{1}}}{k\theta A_{\mathrm{2}}}\left(\frac{{\beta }_{\mathrm{1}}p}{q+{\gamma }_{\mathrm{1}}}+\frac{{\beta }_{\mathrm{2}}((\mathrm{1}-p){\gamma }_{\mathrm{1}}+q)}{(q+{\gamma }_{\mathrm{1}})(m+{\gamma }_{\mathrm{2}})}\right)\frac{〈k^{\mathrm{2}}〉}{〈k〉}>\mathrm{1}-\frac{A_{\mathrm{1}}}{\theta A_{\mathrm{2}}}\left(\frac{{\beta }_{\mathrm{1}}p}{q+{\gamma }_{\mathrm{1}}}+\frac{{\beta }_{\mathrm{2}}((\mathrm{1}-p){\gamma }_{\mathrm{1}}+q)}{(q+{\gamma }_{\mathrm{1}})(m+{\gamma }_{\mathrm{2}})}\right)\mathrm{。}\end{array}

为使f(\theta )>\mathrm{0}(\theta >\mathrm{0})， 取{\theta }^{\mathrm{*}}=\mathrm{2}\frac{A_{\mathrm{1}}}{A_{\mathrm{2}}}\left(\frac{{\beta }_{\mathrm{1}}p}{q+{\gamma }_{\mathrm{1}}}\right.+\left.\frac{{\beta }_{\mathrm{2}}((\mathrm{1}-p){\gamma }_{\mathrm{1}}+q)}{(q+{\gamma }_{\mathrm{1}})(m+{\gamma }_{\mathrm{2}})}\right)， 可使f({\theta }^{\mathrm{*}})>\mathrm{0}。因此， f(\theta )存在正解， 即F(\theta )存在一个正解， 使得模型（1）存在唯一的正平衡点。定理1得证。

下面对模型在无病平衡点处的稳定性进行分析。

在无病平衡点E₀处将模型（1）线性化可以得到系统

\begin{array}{c}\frac{\mathrm{d}X(t)}{\mathrm{d}t}={\left[\begin{array}{cccccc}\mathrm{0}& \mathrm{0}& -{\beta }_{\mathrm{1}}\mathit{P}& -{\beta }_{\mathrm{2}}\mathit{P}& \mathrm{0}& \delta {\mathit{I}}_n\\ \mathrm{0}& -\alpha {\mathit{I}}_n& {\beta }_{\mathrm{1}}\mathit{P}& {\beta }_{\mathrm{2}}\mathit{P}& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& p\alpha {\mathit{I}}_n& -{\gamma }_{\mathrm{1}}{\mathit{I}}_n& \mathrm{0}& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& (\mathrm{1}-p)\alpha {\mathit{I}}_n& \mathrm{0}& -(m+{\gamma }_{\mathrm{2}}){\mathit{I}}_n& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& \mathrm{0}& m{\mathit{I}}_n& -{\gamma }_{\mathrm{3}}{\mathit{I}}_n& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& {\gamma }_{\mathrm{1}}{\mathit{I}}_n& {\gamma }_{\mathrm{2}}{\mathit{I}}_n& {\gamma }_{\mathrm{3}}{\mathit{I}}_n& -\delta {\mathit{I}}_n\end{array}\right]}_{\mathrm{6}n\times \mathrm{6}n}·X(t)+\\ {\left[\begin{array}{cccccc}\mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& -q{\mathit{I}}_n& \mathrm{0}& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& q{\mathit{I}}_n& \mathrm{0}& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}\end{array}\right]}_{\mathrm{6}n\times \mathrm{6}n}·X(t-\tau )\mathrm{。}\end{array}（8）

系统（8）的系数矩阵的特征多项式为

f(\lambda )={\left|\begin{array}{cccccc}\lambda {\mathit{I}}_n& \mathrm{0}& {\beta }_{\mathrm{1}}\mathit{P}& {\beta }_{\mathrm{2}}\mathit{P}& \mathrm{0}& -\delta {\mathit{I}}_n\\ \mathrm{0}& (\lambda +\alpha ){\mathit{I}}_n& -{\beta }_{\mathrm{1}}\mathit{P}& -{\beta }_{\mathrm{2}}\mathit{P}& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& -p\alpha {\mathit{I}}_n& (\lambda +{\gamma }_{\mathrm{1}}+q{e}^{-\lambda \tau }){\mathit{I}}_n& \mathrm{0}& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& -(\mathrm{1}-p)\alpha {\mathit{I}}_n& -q{\mathrm{e}}^{-\lambda \tau }{\mathit{I}}_n& (\lambda +{\gamma }_{\mathrm{2}}+m){\mathit{I}}_n& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& \mathrm{0}& -m{\mathit{I}}_n& (\lambda +{\gamma }_{\mathrm{3}}){\mathit{I}}_n& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& -{\gamma }_{\mathrm{1}}{\mathit{I}}_n& -{\gamma }_{\mathrm{2}}{\mathit{I}}_n& -{\gamma }_{\mathrm{3}}{\mathit{I}}_n& (\lambda +\delta ){\mathit{I}}_n\end{array}\right|}_{\mathrm{6}n\times \mathrm{6}n}=

\begin{array}{c}\left|\lambda {\mathit{I}}_n\right|\left|(\lambda +\delta ){\mathit{I}}_n\right|\left|(\lambda +{\gamma }_{\mathrm{3}}){\mathit{I}}_n\right|\left|(\lambda +\alpha )(\lambda +{\gamma }_{\mathrm{1}}+q{\mathrm{e}}^{-\lambda \tau })(\lambda +{\gamma }_{\mathrm{2}}+m){\mathit{I}}_n-\right.\\ {\beta }_{\mathrm{1}}p\alpha (\lambda +{\gamma }_{\mathrm{2}}+m)\mathit{P}-{\beta }_{\mathrm{2}}[(\mathrm{1}-p)\alpha \lambda +(\mathrm{1}-p)\alpha {\gamma }_{\mathrm{1}}+\alpha q{\mathrm{e}}^{-\lambda \tau }]\mathit{P}\left|=\right.\\ \left|\lambda {\mathit{I}}_n\right|\left|(\lambda +\delta ){\mathit{I}}_n\right|\left|(\lambda +{\gamma }_{\mathrm{3}}){\mathit{I}}_n\right|\left|{\lambda }^{\mathrm{3}}+({\gamma }_{\mathrm{1}}+{\gamma }_{\mathrm{2}}+\alpha +m){\lambda }^{\mathrm{2}}\right|+\\ [{\gamma }_{\mathrm{1}}{\gamma }_{\mathrm{2}}+{\gamma }_{\mathrm{1}}m+\alpha ({\gamma }_{\mathrm{1}}+{\gamma }_{\mathrm{2}}+m)-(p\alpha {\beta }_{\mathrm{1}}+(\mathrm{1}-p)\alpha {\beta }_{\mathrm{2}})\mathit{P}]\lambda +\\ [\alpha {\gamma }_{\mathrm{1}}{\gamma }_{\mathrm{2}}+\alpha {\gamma }_{\mathrm{1}}m-(p\alpha ({\gamma }_{\mathrm{2}}+m){\beta }_{\mathrm{1}}+(\mathrm{1}-p)\alpha {\gamma }_{\mathrm{1}}{\beta }_{\mathrm{2}})\mathit{P}]+\\ +[{\lambda }^{\mathrm{2}}+(\alpha +{\gamma }_{\mathrm{2}}+m)\lambda +\alpha ({\gamma }_{\mathrm{2}}+m)-\alpha {\beta }_{\mathrm{2}}P]q{\mathrm{e}}^{-\lambda \tau }\left|\mathrm{。}\right.\end{array}

根据舒尔引理， 存在一个酉矩阵 U，使得\mathit{P}={\mathit{U}}^{-\mathrm{1}}\overline{\mathit{P}}\mathit{U}， 可以得到公式

\begin{array}{c}f(\lambda )=\left|\lambda {\mathit{I}}_n\right|\left|(\lambda +\delta ){\mathit{I}}_n\right|\left|(\lambda +{\gamma }_{\mathrm{3}}){\mathit{I}}_n\right|\left|{\lambda }^{\mathrm{3}}+({\gamma }_{\mathrm{1}}+{\gamma }_{\mathrm{2}}+\alpha +m){\lambda }^{\mathrm{2}}+\right.\\ [{\gamma }_{\mathrm{1}}{\gamma }_{\mathrm{2}}+{\gamma }_{\mathrm{1}}m+\alpha ({\gamma }_{\mathrm{1}}+{\gamma }_{\mathrm{2}}+m)-(p\alpha {\beta }_{\mathrm{1}}+(\mathrm{1}-p)\alpha {\beta }_{\mathrm{2}})\stackrel{-}{\mathit{P}}]\lambda +\\ [\alpha {\gamma }_{\mathrm{1}}{\gamma }_{\mathrm{2}}+\alpha {\gamma }_{\mathrm{1}}m-(p\alpha ({\gamma }_{\mathrm{2}}+m){\beta }_{\mathrm{1}}+(\mathrm{1}-p)\alpha {\gamma }_{\mathrm{1}}{\beta }_{\mathrm{2}})\stackrel{-}{\mathit{P}}]+\\ [{\lambda }^{\mathrm{2}}+(\alpha +{\gamma }_{\mathrm{2}}+m)\lambda +\alpha ({\gamma }_{\mathrm{2}}+m)-\alpha {\beta }_{\mathrm{2}}\stackrel{-}{\mathit{P}}]q{\mathrm{e}}^{-\lambda \tau }\left|=\right.\\ {\lambda }^n{(\lambda +\delta )}^n(\lambda +{\gamma }_{\mathrm{3}}{)}^n{(\lambda +\alpha )}^{n-\mathrm{1}}(\lambda +m+{\gamma }_{\mathrm{2}}{)}^{n-\mathrm{1}}(\lambda +{\gamma }_{\mathrm{1}}+q{\mathrm{e}}^{-\lambda \tau }{)}^{n-\mathrm{1}}·\\ （{\lambda }^{\mathrm{3}}+q_{\mathrm{1}}{\lambda }^{\mathrm{2}}+q_{\mathrm{2}}\lambda +q_{\mathrm{3}}+({\lambda }^{\mathrm{2}}+q_{\mathrm{4}}\lambda +q_{\mathrm{5}})q{\mathrm{e}}^{-\lambda \tau })，\end{array}（9）

其中，

\stackrel{-}{\mathit{P}}=\left[\begin{array}{cccc}K& \frac{k_{\mathrm{1}}{k}_{\mathrm{2}}p(k_{\mathrm{2}})}{〈k〉}& \cdots & \frac{k_{\mathrm{1}}{k}_{n}p(k_n)}{〈k〉}\\ \mathrm{0}& \mathrm{0}& \cdots & \mathrm{0}\\ ⋮& ⋮& ⋮& ⋮\\ \mathrm{0}& \mathrm{0}& \cdots & \mathrm{0}\end{array}\right]，

K=\frac{{\displaystyle \sum _{k=\mathrm{1}}^n}{k}^{\mathrm{2}}p(k)}{〈k〉}=\frac{〈k^{\mathrm{2}}〉}{〈k〉}，

q_{\mathrm{1}}=\alpha +{\gamma }_{\mathrm{1}}+{\gamma }_{\mathrm{2}}+m，

q_{\mathrm{2}}=\alpha {\gamma }_{\mathrm{1}}+\alpha {\gamma }_{\mathrm{2}}+\alpha m+{\gamma }_{\mathrm{1}}{\gamma }_{\mathrm{2}}+{\gamma }_{\mathrm{1}}m-[(\mathrm{1}-p){\beta }_{\mathrm{2}}+p{\beta }_{\mathrm{1}}]\alpha K，

\begin{array}{c}{q}_{\mathrm{3}}=\alpha {\gamma }_{\mathrm{1}}{\gamma }_{\mathrm{2}}+\alpha {\gamma }_{\mathrm{1}}m-[(\mathrm{1}-p){\gamma }_{\mathrm{1}}{\beta }_{\mathrm{2}}+\\ p(m-{\gamma }_{\mathrm{2}}){\beta }_{\mathrm{1}}]\alpha K\end{array}，

q_{\mathrm{4}}=\alpha +m+{\gamma }_{\mathrm{2}}，

q_{\mathrm{5}}=\alpha {\gamma }_{\mathrm{2}}+\alpha m-\alpha {\beta }_{\mathrm{2}}K。

为便于后续应用， 假设以下条件成立， 即

由式（9）可知， 对应有n重特征根\lambda=0， -\delta， -{\gamma }_{\mathrm{3}}， n-\mathrm{1}重特征根\lambda=-\alpha， -(m+{\gamma }_{\mathrm{2}})， 其余特征根满足以下两个方程

f_{\mathrm{1}}(\lambda )=\lambda +{\gamma }_{\mathrm{1}}+q{\mathrm{e}}^{-\lambda \tau }，

f_{\mathrm{2}}(\lambda )={\lambda }^{\mathrm{3}}+q_{\mathrm{1}}{\lambda }^{\mathrm{2}}+q_{\mathrm{2}}\lambda +q_{\mathrm{3}}+({\lambda }^{\mathrm{2}}+q_{\mathrm{4}}\lambda +q_{\mathrm{5}})q{\mathrm{e}}^{-\lambda \tau }。

下面进行分析：

1） \tau =\mathrm{0}

f(\lambda )={\lambda }^n{(\lambda +\delta )}^n(\lambda +{\gamma }_{\mathrm{3}}{)}^n{(\lambda +\alpha )}^{n-\mathrm{1}}·

(\lambda +m+{\gamma }_{\mathrm{2}}{)}^{n-\mathrm{1}}(\lambda +{\gamma }_{\mathrm{1}}{+q)}^{n-\mathrm{1}}({\lambda }^{\mathrm{3}}+(q_{\mathrm{1}}+q){\lambda }^{\mathrm{2}}+

(q_{\mathrm{2}}+q_{\mathrm{4}}q)\lambda +q_{\mathrm{3}}+q_{\mathrm{5}}q)\mathrm{。}

当q_{\mathrm{2}}+q_{\mathrm{4}}q>\mathrm{0}并且q_{\mathrm{3}}+q_{\mathrm{5}}q>\mathrm{0}时， f(\lambda )的特征值的实部小于等于0， 即

\left\{\begin{array}{l}\alpha {\gamma }_{\mathrm{1}}+\alpha {\gamma }_{\mathrm{2}}+\alpha m+{\gamma }_{\mathrm{1}}{\gamma }_{\mathrm{2}}+{\gamma }_{\mathrm{1}}m+\alpha q+mq+{\gamma }_{\mathrm{2}}q-((\mathrm{1}-p){\beta }_{\mathrm{2}}+p{\beta }_{\mathrm{1}})\alpha K=\\ \alpha {\gamma }_{\mathrm{1}}+\alpha {\gamma }_{\mathrm{2}}+\alpha m+{\gamma }_{\mathrm{1}}{\gamma }_{\mathrm{2}}+{\gamma }_{\mathrm{1}}m+\alpha q+mq+{\gamma }_{\mathrm{2}}q(\mathrm{1}-R_{\mathrm{1}}),\\ \alpha ({\gamma }_{\mathrm{1}}{\gamma }_{\mathrm{2}}+{\gamma }_{\mathrm{1}}m+{\gamma }_{\mathrm{2}}q+mq)-[(\mathrm{1}-p){\gamma }_{\mathrm{1}}{\beta }_{\mathrm{2}}+p(m+{\gamma }_{\mathrm{2}}){\beta }_{\mathrm{1}}]\alpha K-\alpha {\beta }_{\mathrm{2}}qK=\\ \alpha (q+{\lambda }_{\mathrm{1}})(m+{\gamma }_{\mathrm{2}})(\mathrm{1}-R_{\mathrm{0}}),\end{array}\right.

其中，

\begin{array}{c}{R}_{\mathrm{1}}=\\ \frac{((\mathrm{1}-p){\beta }_{\mathrm{2}}+p{\beta }_{\mathrm{1}})\alpha K}{\alpha {\gamma }_{\mathrm{1}}+\alpha {\gamma }_{\mathrm{2}}+\alpha m+{\gamma }_{\mathrm{1}}{\gamma }_{\mathrm{2}}+{\gamma }_{\mathrm{1}}m+\alpha q+mq+{\gamma }_{\mathrm{2}}q}\end{array}。

因此， 当， 时， q_{\mathrm{2}}+q_{\mathrm{4}}q>\mathrm{0}且q_{\mathrm{3}}+q_{\mathrm{5}}q>\mathrm{0}， f(\lambda )的特征值的实部小于等于0。

2） \tau >\mathrm{0}

假设系统（9）有一对纯虚特征值\pm \omega \mathrm{i}（\omega \ne \mathrm{0}）， 显而易见， \lambda满足f_{\mathrm{1}}(\lambda )=\lambda +{\gamma }_{\mathrm{1}}+q{\mathrm{e}}^{-\lambda \tau }=\mathrm{0}或者f_{\mathrm{2}}(\lambda )={\lambda }^{\mathrm{3}}+q_{\mathrm{1}}{\lambda }^{\mathrm{2}}+q_{\mathrm{2}}\lambda +q_{\mathrm{3}}+({\lambda }^{\mathrm{2}}+q_{\mathrm{4}}\lambda +q_{\mathrm{5}})q{\mathrm{e}}^{-\lambda \tau }=\mathrm{0}。令\lambda =\omega \mathrm{i}， 讨论两种情况。

情况 1： \lambda满足f_{\mathrm{1}}(\lambda )=\lambda +{\gamma }_{\mathrm{1}}+q{\mathrm{e}}^{-\lambda \tau }=\mathrm{0}。

将\lambda =\omega \mathrm{i}代入f_{\mathrm{1}}(\lambda )=\lambda +{\gamma }_{\mathrm{1}}+q{\mathrm{e}}^{-\lambda \tau }=\mathrm{0}可得 \omega \mathrm{i}+{\gamma }_{\mathrm{1}}+q(\mathrm{c}\mathrm{o}\mathrm{s}(\omega \tau )-\mathrm{i}\mathrm{s}\mathrm{i}\mathrm{n}(\omega \tau ))=\mathrm{0}。

分离实部和虚部可得

\left\{\begin{array}{l}{\gamma }_{\mathrm{1}}+q\mathrm{c}\mathrm{o}\mathrm{s}(\omega \tau )=\mathrm{0}，\\ \omega -q\mathrm{s}\mathrm{i}\mathrm{n}(\omega \tau )=\mathrm{0}\mathrm{。}\end{array}\right.

因为\mathrm{s}\mathrm{i}\mathrm{n}{(\omega \tau )}^{\mathrm{2}}+\mathrm{c}\mathrm{o}\mathrm{s}{(\omega \tau )}^{\mathrm{2}}=\mathrm{1}， 所以可以得到{{\omega }_{\mathrm{1}}}^{\mathrm{2}}=q^{\mathrm{2}}-{{\gamma }_{\mathrm{1}}}^{\mathrm{2}}， 因此， 只有当q>{\gamma }_{\mathrm{1}}时， 存在{\omega }_{\mathrm{1}}=\sqrt[]{q^{\mathrm{2}}-{\gamma }_{\mathrm{1}}^{\mathrm{2}}}， 对应的{\tau }_{\mathrm{1}}^j=\frac{\mathrm{1}}{\sqrt[]{q^{\mathrm{2}}-{\gamma }_{\mathrm{1}}^{\mathrm{2}}}}·\left(\mathrm{a}\mathrm{r}\mathrm{c}\mathrm{s}\mathrm{i}\mathrm{n}(\frac{\sqrt[]{q^{\mathrm{2}}-{\gamma }_{\mathrm{1}}^{\mathrm{2}}}}{q})+\mathrm{2}j\mathrm{\pi }\right)， j=\mathrm{0}，\mathrm{1}，\cdots。当时， 不存在实的{\omega }_{\mathrm{1}}。

情况 2： \lambda满足f_{\mathrm{2}}(\lambda )={\lambda }^{\mathrm{3}}+q_{\mathrm{1}}{\lambda }^{\mathrm{2}}+q_{\mathrm{2}}\lambda +q_{\mathrm{3}}+({\lambda }^{\mathrm{2}}+q_{\mathrm{4}}\lambda +q_{\mathrm{5}})q{\mathrm{e}}^{-\lambda \tau }=\mathrm{0}。

类似地， 可以得到

\left\{\begin{array}{l}-q_{\mathrm{1}}{\omega }^{\mathrm{2}}+q_{\mathrm{3}}+(q_{\mathrm{5}}q-w^{\mathrm{2}}q)\mathrm{c}\mathrm{o}\mathrm{s}(\omega \tau )+q_{\mathrm{4}}\omega q\mathrm{s}\mathrm{i}\mathrm{n}(\omega \tau )=\mathrm{0}，\\ -{\omega }^{\mathrm{3}}+q_{\mathrm{2}}\omega +q_{\mathrm{4}}\omega q\mathrm{c}\mathrm{o}\mathrm{s}(\omega \tau )-(q_{\mathrm{5}}q-w^{\mathrm{2}}q)\mathrm{s}\mathrm{i}\mathrm{n}(\omega \tau )=\mathrm{0},\end{array}\right.

则{\omega }_{\mathrm{2}}满足

({\omega }^{\mathrm{2}}{)}^{\mathrm{3}}+(q_{\mathrm{1}}^{\mathrm{2}}+\mathrm{2}{q}_{\mathrm{2}}-q^{\mathrm{2}})({\omega }^{\mathrm{2}}{)}^{\mathrm{2}}+({q_{\mathrm{2}}}^{\mathrm{2}}-\mathrm{2}{q}_{\mathrm{1}}{q}_{\mathrm{3}}+\mathrm{2}{q}_{\mathrm{5}}{q}^{\mathrm{2}}-q_{\mathrm{4}}^{\mathrm{2}}{q}^{\mathrm{2}}){\omega }^{\mathrm{2}}+q_{\mathrm{3}}^{\mathrm{2}}-q_{\mathrm{5}}^{\mathrm{2}}{q}^{\mathrm{2}}=\mathrm{0}。

因此， 当假设H_{\mathrm{1}}和H_{\mathrm{2}}或者H_{\mathrm{1}}和H_{\mathrm{3}}成立时， 存在正根{\omega }^{\mathrm{2}}， 可以得到{\omega }_{\mathrm{2}}和

{\tau }_{\mathrm{2}}^j=\frac{\mathrm{1}}{{\omega }_{\mathrm{2}}}\left(\mathrm{a}\mathrm{r}\mathrm{c}\mathrm{s}\mathrm{i}\mathrm{n}\left(\frac{(q_{\mathrm{5}}-{\omega }^{\mathrm{2}})(-{\omega }^{\mathrm{3}}+q_{\mathrm{2}}\omega )+q_{\mathrm{4}}\omega (q_{\mathrm{1}}{\omega }^{\mathrm{2}}-q_{\mathrm{3}})}{q_{\mathrm{4}}\omega q+(q_{\mathrm{5}}-{\omega }^{\mathrm{2}}{)}^{\mathrm{2}}q})+\mathrm{2}\mathrm{i}\mathrm{\pi }\right)\right)，j=\mathrm{0}，\mathrm{1}，\cdots。

下面回到分支分析， 用\tau作为分岔参数， f_{\mathrm{1}}(\lambda )和f_{\mathrm{2}}(\lambda )都是关于\tau的函数， 令\lambda (\tau )=\alpha (\tau )+\omega (\tau )\mathrm{i}， 分岔参数的初值{\tau }_{\mathrm{0}}有\alpha ({\tau }_{\mathrm{0}})=\mathrm{0}， \omega ({\tau }_{\mathrm{0}})={\omega }_{\mathrm{0}}（{\omega }_{\mathrm{0}}>\mathrm{0}）， 为了定义Hopf分岔， 有引理 1。

引理 1　1） 如果q>{\gamma }_{\mathrm{1}}， 则f_{\mathrm{1}}(\lambda )=\mathrm{0}的横截性条件为

{\left.{\left(\frac{\mathrm{d}(\mathrm{R}\mathrm{e}\lambda )}{\mathrm{d}\tau }\right)}^{-\mathrm{1}}\right|}_{{}_{\tau ={\tau }_{\mathrm{1}}^j}}>\mathrm{0}， j=\mathrm{0}，\mathrm{1}，\dots \mathrm{。}

2） 当假设H_{\mathrm{1}}和H_{\mathrm{2}}或者H_{\mathrm{1}}和H_{\mathrm{3}}成立时， f_{\mathrm{2}}(\lambda )=\mathrm{0}的横截性条件为

{\left.{\left(\frac{\mathrm{d}(\mathrm{R}\mathrm{e}\lambda )}{\mathrm{d}\tau }\right)}^{-\mathrm{1}}\right|}_{{}_{\tau ={\tau }_{\mathrm{2}}^j}}>\mathrm{0}， j=\mathrm{0}，\mathrm{1}，\dots \mathrm{。}

证明 1） 对f_{\mathrm{1}}(\lambda )关于\tau求导， 可得(\mathrm{1}-q\tau {\mathrm{e}}^{-\lambda \tau })\frac{\mathrm{d}\lambda }{\mathrm{d}\tau }-q\lambda {\mathrm{e}}^{-\lambda \tau }=\mathrm{0}，也就是说， {\left(\frac{\mathrm{d}\lambda }{\mathrm{d}t}\right)}^{-\mathrm{1}}=\frac{\mathrm{1}-q\tau {\mathrm{e}}^{-\lambda \tau }}{q\lambda {\mathrm{e}}^{-\lambda \tau }}=-\frac{\mathrm{1}}{\lambda (\lambda +{\gamma }_{\mathrm{1}})}-\frac{\tau }{\lambda }。

故

\begin{array}{c}sign{\left\{\mathrm{R}\mathrm{e}\frac{\mathrm{d}\lambda }{\mathrm{d}\tau }\right\}}_{\tau ={\tau }_{\mathrm{1}}^j}=sign{\left\{\mathrm{R}\mathrm{e}\frac{\mathrm{d}\tau }{\mathrm{d}\lambda }\right\}}_{\lambda ={\omega }_{\mathrm{1}}\mathrm{i}}=\\ \mathrm{s}ign{\left\{\mathrm{R}\mathrm{e}\left(-\frac{\mathrm{1}}{\lambda (\lambda +{\gamma }_{\mathrm{1}})}-\frac{\tau }{\lambda }\right)\right\}}_{\lambda ={\omega }_{\mathrm{1}}\mathrm{i}}=\\ \mathrm{s}ign\left\{\mathrm{R}\mathrm{e}\left(-\frac{\mathrm{1}}{-{\omega }^{\mathrm{2}}+{\gamma }_{\mathrm{1}}\omega \mathrm{i}}-\frac{\tau }{\omega \mathrm{i}}\right)\right\}=sign\left\{\frac{\mathrm{1}}{{\omega }^{\mathrm{2}}+{\gamma }_{\mathrm{1}}^{\mathrm{2}}}\right\}>\mathrm{0},\end{array}

因此， {\left.{\left(\frac{\mathrm{d}(\mathrm{R}\mathrm{e}\lambda )}{\mathrm{d}\tau }\right)}^{-\mathrm{1}}\right|}_{{}_{\tau ={\tau }_{\mathrm{1}}^j}}>\mathrm{0}。

2） 利用相似的方法， 可以得到{\left.{\left(\frac{\mathrm{d}(\mathrm{R}\mathrm{e}\lambda )}{\mathrm{d}\tau }\right)}^{-\mathrm{1}}\right|}_{{}_{\tau ={\tau }_{\mathrm{2}}^j}}>\mathrm{0}， j=\mathrm{0}，\mathrm{1}，\dots。

引理1得证。

定理 2　1） 如果R_{\mathrm{0}}>\mathrm{1}， 对于任意\tau \ge \mathrm{0}， 无病平衡点E_{\mathrm{0}}是不稳定的。

2） 如果， ， 并且假设H_{\mathrm{1}}不成立， 则对于任意\tau \ge \mathrm{0}， 模型（1）在无病平衡点E_{\mathrm{0}}处稳定。

3） 如果， ， q>{\gamma }_{\mathrm{1}}并且假设H_{\mathrm{1}}不成立。如果\tau \in [\mathrm{0},{\tau }_{\mathrm{1}}^{\mathrm{0}})， 无病平衡点E_{\mathrm{0}}是稳定的； 如果\tau >{\tau }_{\mathrm{1}}^{\mathrm{0}}， 无病平衡点E_{\mathrm{0}}是不稳定的。在\tau ={\tau }_{\mathrm{1}}^j， j=\mathrm{1}，\mathrm{2}，\dots， 发生Hopf分岔。

4） 如果， ， 并且假设H_{\mathrm{1}}和假设H_{\mathrm{2}}或者假设H_{\mathrm{1}}和假设H_{\mathrm{3}}成立。如果\tau \in [\mathrm{0},{\tau }_{\mathrm{2}}^{\mathrm{0}})， 无病平衡点E_{\mathrm{0}}是稳定的； 如果\tau >{\tau }_{\mathrm{2}}^{\mathrm{0}}， 无病平衡点E_{\mathrm{0}}是不稳定的。在\tau ={\tau }_{\mathrm{2}}^j， j=1，2，…， 发生Hopf分岔。

5） 如果， ， q>{\gamma }_{\mathrm{1}}并且假设H_{\mathrm{1}}和假设H_{\mathrm{2}}或者假设H_{\mathrm{1}}和假设H_{\mathrm{3}}成立。如果\tau \in [\mathrm{0},\mathrm{m}\mathrm{i}\mathrm{n}({\tau }_{\mathrm{1}}^{\mathrm{0}},{\tau }_{\mathrm{2}}^{\mathrm{0}}))， 无病平衡点E_{\mathrm{0}}是稳定的； 如果\tau >\mathrm{m}\mathrm{i}\mathrm{n}({\tau }_{\mathrm{1}}^{\mathrm{0}},{\tau }_{\mathrm{2}}^{\mathrm{0}})， 无病平衡点E_{\mathrm{0}}是不稳定的。在\tau ={\tau }_{\mathrm{1}}^j({\tau }_{\mathrm{2}}^j)时， 发生Hopf分岔。

3 数值模拟

本文基于一个节点数为1 000的无标度网络进行数值模拟， 具体来说就是从初始节点m₀=5开始， 每增加一个节点就新增三条连边， 连边连接到该节点的概率与该节点本身已有的边的数目成正比。由于网络的生成具有随机性， 所以在实验中进行30次数值模拟， 然后取其平均值， 得到关于S(t)， E(t)， A(t)， I(t)， Q(t)， R(t)的时间演化曲线。

首先， 研究感染率{\beta }_{\mathrm{1}}与{\beta }_{\mathrm{2}}对疾病传播的影响。由图 1 可知， 当R₀<1时， 随着两种感染率的降低， 住院治疗的人群显著减少， 易感者也从先减再增变为直接增加。

当R₀>1时， 随着时间变化， 传染病趋于稳定， 并且随着两种感染率的降低， 易感者比例逐渐增大， 说明降低染病率可以有效抑制疾病的传播， 有助于疫情防控。

其次， 对影响R_{\mathrm{0}}的各个参数进行PRCC敏感性分析 ，结果如图 2 所示。由图 2 可以看出， p， {\beta }_{\mathrm{1}}和{\beta }_{\mathrm{2}}对R_{\mathrm{0}}有正影响， q， m， p， {\gamma }_{\mathrm{1}}和{\gamma }_{\mathrm{2}}对R_{\mathrm{0}}有负影响。挑选其中对R_{\mathrm{0}}影响较大的两个参数{\beta }_{\mathrm{1}}， {\gamma }_{\mathrm{1}}， 分析其对R_{\mathrm{0}}的影响。其中， 参数{\beta }_{\mathrm{1}}， {\beta }_{\mathrm{2}}， q， m， p， {\gamma }_{\mathrm{1}}， {\gamma }_{\mathrm{2}}服从正态分布。

由图 3 可以看出， R_{\mathrm{0}}随着{\gamma }_{\mathrm{1}}的增大而减小， 随着{\beta }_{\mathrm{1}}的增大而增大， 并且可以观察到{\gamma }_{\mathrm{1}}对R_{\mathrm{0}}的影响更为显著， 说明在传染病的控制过程中可以通过提高治疗率来更加有效地抑制疾病的传播。

由图 4 可以看出： 在区域Ⅰ中， 说明疾病会灭绝； 在区域Ⅱ中R_{\mathrm{0}}>\mathrm{1}， 说明疾病将会持续存在； 当时， R_{\mathrm{0}}始终处于区域Ⅰ中， 传染病将会灭绝； 当{\beta }_{\mathrm{1}}>\mathrm{0.22}时， R_{\mathrm{0}}始终处于区域Ⅱ中， 此时无论如何调控{\gamma }_{\mathrm{1}}， 疾病终将持续爆发； 当时， 区域Ⅰ和区域Ⅱ的边界是函数{\beta }_{\mathrm{1}}=f({\gamma }_{\mathrm{1}})， 每取定一个{{\beta }_{\mathrm{1}}}^{\mathrm{*}}， 都有唯一{\gamma }_{\mathrm{1}}^{\mathrm{*}}， 使得点（{\gamma }_{\mathrm{1}}^{\mathrm{*}}， {\beta }_{\mathrm{1}}^{\mathrm{*}}）在边界上， 这意味着， 传染病爆发后， 染病率\beta在一定范围内时， 可以通过精确控制治愈率{\gamma }_{\mathrm{1}}的值， 使疾病趋于灭绝， 这对于抑制传染病传播和节约医疗资源有一定的意义。

4 结　论

本文考虑了带有时滞的SEAIQRS传染病模型。首先， 利用下一代矩阵法求得模型的基本再生数并且证明了正平衡点的存在唯一性。随后， 对模型进行了动力学分析， 证明模型在无病平衡点处Hopf分岔的存在性。时滞\tau的加入可能会使得正常情况下稳定的无病平衡点出现周期性的爆发， 但由于没有合适的数值模拟方法， 所以这里并不能直观地表示出来。最后， 进行了相关数值模拟， 通过对模型参数的分析可以发现， 提高治疗率和降低染病率（佩戴口罩， 保持社交距离）可以有效抑制疾病传播。本文没有考虑出生与死亡等因素， 也没有考虑其他因素对模型的扰动， 较为理想。后续可以改善模型， 使模型更加贴切实际， 如考虑动态网络等。

参考文献

[1] WEI G W， WEI Y， LI X F， et al. Infectious disease pandemic and permanent volatility of international stock markets： A long-term perspective［J］. Finance Research Letters， 2021， 40： 101709.

[2] 顾梓玉， 党元辰， 何沛妍， 等. 复杂网络中SAIQRS传染病模型动力学分析［J］. 河北师范大学学报（自然科学版）， 2023， 47（1）： 41-46.

[3] GU Ziyu， DANG Yuanchen， HE Peiyan， et al. Dynamics analysis of SAIQRS infectious disease model in complex networks［J］. Journal of Hebei Normal University（Natural Science）， 2023， 47（1）： 41-46. （in Chinese）

[4] SOTTILE S， LIU X. Time-varying epidemic transmission in heterogeneous networks and applications to measles［J］. Journal of Biological Systems， 2020， 28（4）： 901-926.

[5] MILLER， JOEL C. Percolation and epidemics in random clustered networks［J］. Physical Review E Statistical Nonlinear & Soft Matter Physics， 2009， 80（2）： 020901.

[6] ZHANG R X， JJN Z， LI S P. Epidemic spreading with time delay on complex networks［J］. Acta Mathematicae Applicatae Sinica， 2016， 32（2）： 319-326.

[7] PIET V M， CLAUDIO C， ROMUALDO P S， et al. Epidemic processes in complex networks［J］. Reviews of Modern Physics， 2015， 87（3）： 925.

[8] CHRISTLEY R M， PINCHBECK G L， BOWERS R G， et al. Infection in social networks： Using network analysis to identify high-risk individuals［J］. American Journal of Epidemiology， 2005（10）： 10.

[9] SALATHÉ M， KAZANDJIEVA M， LEE J W， et al. A high-resolution human contact network for infectious disease transmission［J］. Proceedings of the National Academy of Sciences of the United States of America， 2010， 107（51）： 22020-22025.

[10] 徐云程. 无标度网络上带有非单调发病率的随机SIQR传染病模型研究［D］. 银川： 宁夏大学， 2022.

[11] 赵润东， 孙梅慈， 刘启明.无标度网络上一类具有隔离项的时滞传染病模型研究［J］. 西南师范大学学报（自然科学版）， 2022， 47（1）： 28-35.

[12] ZHAO Rundong， SUN Meici， LIU Qiming. A model study of time-delayed infectious diseases with isolation terms on a scale-free network［J］. Journal of Southwest China Normal University（Natural Science Edition）， 2022， 47（1）： 28-35. （in Chinese）

[13] 柯青.无标度网络上具有不完美接种的SIVS模型的动力学分析与最优控制［D］. 上海： 上海大学， 2020.

[14] LYU W， KE Q， LI K. Dynamical analysis and control strategies of an SIVS epidemic model with imperfect vaccination on scale-free networks［J］. Nonlinear Dynamics， 2020， 99（2）： 1507-1523.

基金资助

国家自然科学基金资助项目(12071445)

国家自然科学基金资助项目(12001501)

山西省自然科学基金资助项目(201901D211275)