受媒体报道影响的SEIAQRS传染病模型分析

Analysis of SEIAQRS Epidemic Model Effected by Media Coverage

摘要

为深入了解新型冠状病毒肺炎的传播机制，建立了具有媒体报道和隔离的SEIAQRS传染病模型。利用基本再生数R_c讨论了平衡点的存在性与稳定性。当时，无病平衡点是全局渐近稳定的，当R_c>\mathrm{1}时，存在唯一的局部渐近稳定的地方病平衡点。分析了模型参数\alpha（有症状感染者的隔离率）、p（有症状感染者的恢复率）和m_{\mathrm{0}}（媒体报道系数）对基本再生数的影响，发现隔离有症状感染者、提高有症状感染者的恢复率和加大常态化宣传传染病知识的媒体报道，能降低基本再生数R_c。最后，通过数值模拟验证了媒体报道和有效隔离能减少感染人数，进而有效控制新型冠状病毒肺炎的传播。

Abstract

A SEIAQRS infectious disease model with media coverage and isolation is established to understand the transmission mechanism of COVID-19. The existence and stability of equilibrium point are discussed by using the basic regeneration number R_c. When ， the disease-free equilibrium is globally asymptotically stable， and when R_c>\mathrm{1}， there is a unique locally asymptotically stable endemic disease equilibrium. The effects of model parameters \alpha（isolation rate of symptomatic infected persons）， p（recovery rate of symptomatic infected persons） and m_{\mathrm{0}}（media coverage coefficient） on the basic regeneration number are analyzed. It is obtained that isolating symptomatic infected persons， improving the recovery rate of symptomatic infected persons and increasing the media coverage of regular dissemination of infectious disease knowledge can reduce the basic regeneration number R_c. Finally， through numerical simulation， media coverage and effective isolation can greatly reduce the number of infected people， and thus effectively control the spread of COVID-19.

Graphical abstract

关键词

媒体报道 / 传播率 / 隔离 / SEIAQRS传染病模型

Key words

media coverage / transmission rate / isolation / SEIAQRS epidemic model

0 引 言

当传染病在人群中暴发时，人类对其反应取决于他们对风险的认知，而这种认知受媒体报道的影响^［1‐6］。媒体报道有利于人们对风险程度和风险地区预防需求的了解，使得带口罩人数增加、公共场所人员密集度降低等，进而助力于疫情管控。广泛的媒体报道和信息的快速流动会对公众产生深刻的心理影响^［7‐8］，比如易感个体产生警觉，避免与受感染个体发生不必要的接触，感染个体自觉隔离、取消行程等。因此，媒体报道可以降低人类接触率和传播概率，许多学者使用动力学模型来描述传染病传播过程中媒体报道发挥的作用^［9‐10］。

在传染病传播动力学建模的研究中，发病率函数的选取极其重要，因为它决定了传染病^［11］的流行趋势。在许多传染病模型中，双线性发病率\beta SI和标准发病率\beta SI/N常被用来表征传染病的发病率，其中，\beta表示传染病的传播率。然而，这些发病率函数没有考虑媒体报道对传染病传播的影响^［12］。考虑到媒体报道，Cui等^［9，13］采用函数\beta {\mathrm{e}}^{-mI}来描述媒体影响因素，提出了一种新的传播速率函数。然而，当I\to \mathrm{\infty }时，\beta {\mathrm{e}}^{-mI}\to \mathrm{0}，其值与m无关，这是不合理的。Liu等^［4］建立了透射系数为\beta {\mathrm{e}}^{-({\alpha }_{\mathrm{1}}E+{\alpha }_{\mathrm{2}}I+{\alpha }_{\mathrm{3}}H)}的EIH模型，其中H表示住院患者，观察到感染可能多次爆发和持续性的周期波动，而描述这种指数递减因子的影响是很难实现的。Cui等^［13］提出一种接触率\beta (I)=c_{\mathrm{1}}-c_{\mathrm{2}}f(I)，以反应媒体报道对传染病传播的影响，其中，c_{\mathrm{1}}是不考虑媒体报道的接触率。Huo等^［14］通过将信息集设置成一个单独的仓室M，并使用指数函数\beta SI{e}^{-\alpha M}作为传播率，其中，\alpha是媒体报道影响疾病传播程度的因子。Sahu等^［15］建立了流行模型，以\beta {\mathrm{e}}^{-\frac{m_{\mathrm{1}}I+m_{\mathrm{2}}H}{N}}\frac{S(I+\eta H)}{N}作为疾病的传播期，提出了一种媒体感知系数m的估计方法。

本文建立了一种传播率为\beta {\mathrm{e}}^{-\frac{m_{\mathrm{1}}I+m_{\mathrm{2}}A+m_{\mathrm{0}}N}{N}}\frac{S(I+\theta A)}{N}的SEIAQRS模型，其中，m_{\mathrm{0}}表示常态化宣传传染病知识的媒体报道系数。该模型相比传统的SEIR模型增加了无症状感染者（A）和隔离者（Q），同时，本文考虑了疾病的复发，即移出者（R）会再一次变为易感者（S）。本模型的隔离对象是感染者而非易感者，相比Shi等^［16］提出的隔离易感者，防疫更具有针对性，该模型更符合新型冠状病毒肺炎传播的实际情况，对于有效预防和控制新型冠状病毒肺炎的传播有着重要的作用。

1 SEIAQRS模型构建

将传染病传播区域的总人数分为易感者（S）、潜伏者（E）、有症状感染者（I）、无症状感染者（A）、隔离者（Q）、移出者（R），其中N=S+E+I+A+Q+R。根据病毒的传播机理，构建了具有媒体报道影响的SEIAQRS传染病传播动力学模型，即

\left\{\begin{array}{l}\frac{\mathrm{d}\tilde{S}}{\mathrm{d}\tilde{t}}=\mathrm{\Lambda }-\tilde{\mu }\tilde{S}-\tilde{\beta }{\mathrm{e}}^{-\frac{m_{\mathrm{1}}\tilde{I}+m_{\mathrm{2}}\tilde{A}+m_{\mathrm{0}}\tilde{N}}{\tilde{N}}}\frac{\tilde{S}(\tilde{I}+\eta \tilde{A})}{\tilde{N}}+{\tilde{\theta }}_{\mathrm{2}}\tilde{R}，\\ \frac{\mathrm{d}\tilde{E}}{\mathrm{d}\tilde{t}}=\tilde{\beta }{\mathrm{e}}^{-\frac{m_{\mathrm{1}}\tilde{I}+m_{\mathrm{2}}\tilde{A}+m_{\mathrm{0}}\tilde{N}}{\tilde{N}}}\frac{\tilde{S}(\tilde{I}+\eta \tilde{A})}{\tilde{N}}-(\tilde{\mu }+{\tilde{\theta }}_{\mathrm{1}}+\tilde{\sigma })\tilde{E}，\\ \frac{\mathrm{d}\tilde{I}}{\mathrm{d}\tilde{t}}=\tilde{\sigma }(\mathrm{1}-\tilde{\rho })\tilde{E}-({\tilde{\mu }}_{\mathrm{1}}+\tilde{\mu }+\tilde{\alpha }+\tilde{p})\tilde{I}，\\ \frac{\mathrm{d}\tilde{A}}{\mathrm{d}\tilde{t}}=\tilde{\sigma }\tilde{\rho }\tilde{E}-(\tilde{\mu }+{\tilde{\omega }}_{\mathrm{1}}+\tilde{\lambda })\tilde{A}，\\ \frac{\mathrm{d}\tilde{Q}}{\mathrm{d}\tilde{t}}=\tilde{\alpha }\tilde{I}+\tilde{\lambda }\tilde{A}-(\tilde{\mu }+{\tilde{\omega }}_{\mathrm{2}})\tilde{Q}，\\ \frac{\mathrm{d}\tilde{R}}{\mathrm{d}\tilde{t}}={\tilde{\omega }}_{\mathrm{1}}\tilde{A}+{\tilde{\omega }}_{\mathrm{2}}\tilde{Q}+\tilde{p}\tilde{I}+{\tilde{\theta }}_{\mathrm{1}}\tilde{E}-(\tilde{\mu }+{\tilde{\theta }}_{\mathrm{2}})\tilde{R}，\end{array}\right.（1）

式中：\mathrm{\Lambda }表示易感者的常数输入率；\tilde{\mu }表示普通人的自然死亡率；{\tilde{\mu }}_{\mathrm{1}}表示有症状感染者的因病死亡率；\tilde{\beta }表示易感者与感染者之间接触的感染率；\eta是降低无症状感染者的修正参数；m_{\mathrm{1}}和m_{\mathrm{2}}分别对应于\tilde{I}和\tilde{A}的媒体覆盖系数；m_{\mathrm{0}}表示常态化宣传传染病知识的媒体报道系数；\tilde{\sigma }表示从暴露到感染级的进展率；\tilde{\rho }表示无症状感染者的比例；{\tilde{\omega }}_{\mathrm{1}}、{\tilde{\omega }}_{\mathrm{2}}、{\tilde{\theta }}_{\mathrm{1}}、\tilde{p}分别表示无症状感染者、隔离者、潜伏者和有症状感染者的恢复率；{\tilde{\theta }}_{\mathrm{2}}表示移出者进入易感级的比率；\tilde{\lambda }、\tilde{\alpha }分别表示无症状感染者和有症状感染者的隔离率。

假设模型(\mathrm{1})满足的原始条件为\tilde{S}(\mathrm{0})={\tilde{S}}_{\mathrm{0}}>\mathrm{0}，\tilde{E}(\mathrm{0})={\tilde{E}}_{\mathrm{0}}>\mathrm{0}，\tilde{I}(\mathrm{0})={\tilde{I}}_{\mathrm{0}}>\mathrm{0}，\tilde{A}(\mathrm{0})={\tilde{A}}_{\mathrm{0}}>\mathrm{0}，\tilde{Q}(\mathrm{0})={\tilde{Q}}_{\mathrm{0}}>\mathrm{0}，\tilde{R}(\mathrm{0})={\tilde{R}}_{\mathrm{0}}>\mathrm{0}，其中，\tilde{N}=\tilde{S}+\tilde{E}+\tilde{I}+\tilde{A}+\tilde{Q}+\tilde{R}。由于\frac{\mathrm{d}\tilde{N}}{\mathrm{d}\tilde{t}}=\mathrm{\Lambda }-\tilde{\mu }\tilde{N}-{\tilde{\mu }}_{\mathrm{1}}\tilde{I}，则模型(\mathrm{1})的约束区域为

\begin{array}{c}\mathrm{\Gamma }=\left\{(\tilde{S},\tilde{E},\tilde{I},\tilde{A},\tilde{Q},\tilde{R})\in \right.\\ {\mathbf{R}}_{+}^{\mathrm{6}}:\tilde{S},\tilde{E},\tilde{I},\tilde{A},\tilde{Q},\tilde{R}\ge \mathrm{0},\end{array}

\left.\tilde{S}+\tilde{E}+\tilde{I}+\tilde{A}+\tilde{Q}+\tilde{R}\le \frac{\mathrm{\Lambda }}{\tilde{\mu }}\right\}\mathrm{。}

考虑无量纲变换

S=\frac{\tilde{S}}{\tilde{N}}，E=\frac{\tilde{E}}{\tilde{N}}，I=\frac{\tilde{I}}{\tilde{N}}，A=\frac{\tilde{A}}{\tilde{N}}，Q=\frac{\tilde{Q}}{\tilde{N}}，R=\frac{\tilde{R}}{\tilde{N}}，N=\frac{\tilde{N}}{{\tilde{N}}^{\mathrm{0}}}，{\tilde{N}}^{\mathrm{0}}=\frac{\mathrm{\Lambda }}{\tilde{\mu }}，t=\tilde{\mu }\tilde{t}。

因为S=\mathrm{1}-E-I-A-Q-R，那么

\begin{array}{c}\frac{\mathrm{d}S}{\mathrm{d}t}=\frac{\mathrm{1}}{N}-\beta {\mathrm{e}}^{-m_{\mathrm{1}}I-m_{\mathrm{2}}A-m_{\mathrm{0}}}S(I+\eta A)+\\ {\theta }_{\mathrm{2}}R-S-\frac{S}{N}\frac{\mathrm{d}N}{\mathrm{d}t}\mathrm{。}\end{array}

因此，可将模型(\mathrm{1})简化为模型(\mathrm{2})，其中模型(\mathrm{2})符合的初始条件为E(\mathrm{0})=E_{\mathrm{0}}>\mathrm{0}，I(\mathrm{0})=I_{\mathrm{0}}>\mathrm{0}，A(\mathrm{0})=A_{\mathrm{0}}>\mathrm{0}，Q(\mathrm{0})=Q_{\mathrm{0}}>\mathrm{0}，R(\mathrm{0})=R_{\mathrm{0}}>\mathrm{0}, N(\mathrm{0})=N_{\mathrm{0}}>\mathrm{0}。其中, \beta =\frac{\tilde{\beta }}{\tilde{\mu }}, {\mu }_{\mathrm{1}}=\frac{{\tilde{\mu }}_{\mathrm{1}}}{\tilde{\mu }},\sigma =\frac{\tilde{\sigma }}{\tilde{\mu }},\rho =\frac{\tilde{\rho }}{\tilde{\mu }}, {\omega }_{\mathrm{1}}=\frac{{\tilde{\omega }}_{\mathrm{1}}}{\tilde{\mu }}, {\omega }_{\mathrm{2}}=\frac{{\tilde{\omega }}_{\mathrm{2}}}{\tilde{\mu }}, {\theta }_{\mathrm{1}}=\frac{{\tilde{\theta }}_{\mathrm{1}}}{\tilde{\mu }}, {\theta }_{\mathrm{2}}=\frac{{\tilde{\theta }}_{\mathrm{2}}}{\tilde{\mu }},p=\frac{\tilde{p}}{\tilde{\mu }}, \lambda =\frac{\lambda }{\tilde{\mu }}, \alpha =\frac{\tilde{\alpha }}{\tilde{\mu }}。

\left\{\begin{array}{l}\frac{\mathrm{d}E}{\mathrm{d}t}=\beta {\mathrm{e}}^{-m_{\mathrm{1}}I-m_{\mathrm{2}}A-m_{\mathrm{0}}}S(I+\eta A)-(\mathrm{1}+{\theta }_{\mathrm{1}}+\sigma )E-\frac{E}{N}\frac{\mathrm{d}N}{\mathrm{d}t}≔f_{\mathrm{1}},\\ \frac{\mathrm{d}I}{\mathrm{d}t}=\sigma (\mathrm{1}-\rho )E-(\mathrm{1}+{\mu }_{\mathrm{1}}+\alpha +p)I-\frac{I}{N}\frac{\mathrm{d}N}{\mathrm{d}t}≔f_{\mathrm{2}},\\ \frac{\mathrm{d}A}{\mathrm{d}t}=\sigma \rho E-(\mathrm{1}+{\omega }_{\mathrm{1}}+\lambda )A-\frac{A}{N}\frac{\mathrm{d}N}{\mathrm{d}t}≔f_{\mathrm{3}},\\ \frac{\mathrm{d}Q}{\mathrm{d}t}=\alpha I+\lambda A-(\mathrm{1}+{\omega }_{\mathrm{2}})Q-\frac{Q}{N}\frac{\mathrm{d}N}{\mathrm{d}t}≔f_{\mathrm{4}},\\ \frac{\mathrm{d}R}{\mathrm{d}t}={\omega }_{\mathrm{1}}A+{\omega }_{\mathrm{2}}Q+pI+{\theta }_{\mathrm{1}}E-{\theta }_{\mathrm{2}}R-R-\frac{R}{N}\frac{\mathrm{d}N}{\mathrm{d}t}≔f_{\mathrm{5}},\\ \frac{\mathrm{d}N}{\mathrm{d}t}=\mathrm{1}-(\mathrm{1}+{\mu }_{\mathrm{1}}I)N≔f_{\mathrm{6}},\end{array}\right.（2）

故模型(\mathrm{2})的可行域为集合

\mathrm{\Omega }=\left\{(E,I,A,Q,R,N):\mathrm{0}\le E,I,A,Q,R,N\le \mathrm{1}\right\}。

2 基本再生数和平衡点的存在性

易知模型(\mathrm{2})的一个无病平衡点E^{\mathrm{0}}=(\mathrm{0,0},\mathrm{0,0},\mathrm{0,1})。利用下一代矩阵^［17］推导出基本再生数（R_c）的表达式为

\mathit{F}=\left[\begin{array}{llll}\mathrm{0}& \beta {\mathrm{e}}^{-m_{\mathrm{0}}}& \beta {\mathrm{e}}^{-m_{\mathrm{0}}}\eta & \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],

\mathit{V}=\left(\begin{array}{cccc}\mathrm{1}+{\theta }_{\mathrm{1}}+\sigma & \mathrm{0}& \mathrm{0}& \mathrm{0}\\ -\sigma (\mathrm{1}-\rho )& \mathrm{1}+{\mu }_{\mathrm{1}}+\alpha +p& \mathrm{0}& \mathrm{0}\\ -\sigma \rho & \mathrm{0}& \mathrm{1}+{\omega }_{\mathrm{1}}+\lambda & \mathrm{0}\\ \mathrm{0}& -\alpha & -\lambda & \mathrm{1}+{\omega }_{\mathrm{2}}\end{array}\right)，

\mathit{F}{\mathit{V}}^{-\mathrm{1}}=\left(\begin{array}{cccc}\frac{\beta {\mathrm{e}}^{-m_{\mathrm{0}}}\sigma (\mathrm{1}-\rho )(\mathrm{1}+{\omega }_{\mathrm{1}}+\lambda )+\beta {\mathrm{e}}^{-m_{\mathrm{0}}}\eta \sigma \rho (\mathrm{1}+{\mu }_{\mathrm{1}}+\alpha +p)}{(\mathrm{1}+{\theta }_{\mathrm{1}}+\sigma )(\mathrm{1}+{\mu }_{\mathrm{1}}+\alpha +p)(\mathrm{1}+{\omega }_{\mathrm{1}}+\lambda )}& \frac{\beta {\mathrm{e}}^{-m_{\mathrm{0}}}}{\mathrm{1}+{\mu }_{\mathrm{1}}+\alpha +p}& \frac{\beta {\mathrm{e}}^{-m_{\mathrm{0}}}\eta }{\mathrm{1}+{\omega }_{\mathrm{1}}+\lambda }& \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{。}

经过推导，基本再生数R_c为

\begin{array}{c}{R}_c=\\ \frac{\beta \sigma (\mathrm{1}-\rho )(\mathrm{1}+{\omega }_{\mathrm{1}}+\lambda )+\beta \eta \sigma \rho (\mathrm{1}+{\mu }_{\mathrm{1}}+\alpha +p)}{(\mathrm{1}+{\theta }_{\mathrm{1}}+\sigma )(\mathrm{1}+{\mu }_{\mathrm{1}}+\alpha +p)(\mathrm{1}+{\omega }_{\mathrm{1}}+\lambda )}\cdot \\ {\mathrm{e}}^{-m_{\mathrm{0}}}=\frac{\beta {\mathrm{e}}^{-m_{\mathrm{0}}}(\mathrm{1}+b_{\mathrm{2}}\eta )}{b_{\mathrm{1}}(\mathrm{1}+{\theta }_{\mathrm{1}}+\sigma )},\end{array}

其中，

b_{\mathrm{1}}=\frac{\mathrm{1}+{\mu }_{\mathrm{1}}+\alpha +p}{\sigma (\mathrm{1}-\rho )}，

b_{\mathrm{2}}=\frac{\rho (\mathrm{1}+{\mu }_{\mathrm{1}}+\alpha +p)}{(\mathrm{1}-\rho )(\mathrm{1}+{\omega }_{\mathrm{1}}+\lambda )}，

b_{\mathrm{3}}=\frac{\alpha +\lambda b_{\mathrm{2}}}{\mathrm{1}+{\omega }_{\mathrm{2}}}。

定理1 如果R_c>\mathrm{1}，则模型(\mathrm{2})存在唯一的地方病平衡点\overline{E}=(E^{\mathrm{*}},I^{\mathrm{*}},A^{\mathrm{*}},Q^{\mathrm{*}},R^{\mathrm{*}},N^{\mathrm{*}})。

证明 令模型(\mathrm{2})中等号右面等于0，将地方病平衡点\overline{E}=(E^{\mathrm{*}},I^{\mathrm{*}},A^{\mathrm{*}},Q^{\mathrm{*}},R^{\mathrm{*}},N^{\mathrm{*}})代入，可得

E^{\mathrm{*}}=b_{\mathrm{1}}{I}^{\mathrm{*}}，A^{\mathrm{*}}=b_{\mathrm{2}}{I}^{\mathrm{*}}，Q^{\mathrm{*}}=b_{\mathrm{3}}{I}^{\mathrm{*}}N^{\mathrm{*}}=\frac{\mathrm{1}}{\mathrm{1}+{\mu }_{\mathrm{1}}{I}^{\mathrm{*}}}，R^{\mathrm{*}}=b_{\mathrm{4}}{I}^{\mathrm{*}}。

其中，b_{\mathrm{4}}=\frac{{\theta }_{\mathrm{1}}{b}_{\mathrm{1}}+{\omega }_{\mathrm{1}}{b}_{\mathrm{2}}+{\omega }_{\mathrm{2}}{b}_{\mathrm{3}}+p}{{\theta }_{\mathrm{2}}+\mathrm{1}}，b₅=1+b₁+b₂+b₃+b₄，m=m_{\mathrm{1}}+m_{\mathrm{2}}{b}_{\mathrm{2}}。

I^{\mathrm{*}}的值可由方程(\mathrm{3})的解给出。

\mathrm{1}-b_{\mathrm{5}}{I}^{\mathrm{*}}=\frac{{\mathrm{e}}^{m{I}^{\mathrm{*}}}}{R_c}\mathrm{。}（3）

假设系数m_{\mathrm{1}}=\mathrm{0}，m_{\mathrm{2}}=\mathrm{0}，即m=\mathrm{0}，可得I^*=\frac{\mathrm{1}}{b_{\mathrm{5}}}\left(\mathrm{1}-\frac{\mathrm{1}}{R_c}\right)。因此，当且仅当R_c>\mathrm{1}时，I^{\mathrm{*}}存在正值，可以知道它是唯一的地方病平衡点，否则I^{\mathrm{*}}的值由方程\mathrm{1}-b_{\mathrm{5}}{I}^{\mathrm{*}}=\frac{{\mathrm{e}}^{m{I}^{\mathrm{*}}}}{R_c}得出。

通过图示法来证明R_c>\mathrm{1}时，I^{\mathrm{*}}的存在性。在图 1 中，画出了曲线\frac{{\mathrm{e}}^{m{I}^{\mathrm{*}}}}{R_c}和直线\mathrm{1}-b_{\mathrm{5}}{I}^{\mathrm{*}}关于I^{\mathrm{*}}在区间[\mathrm{0,1}]上的图形。可知，当时，模型(\mathrm{2})没有地方病平衡点（见图 \mathrm{1}（a））；当R_c>\mathrm{1}时，I^{\mathrm{*}}在正水平方向上存在且是唯一的，即模型(\mathrm{2})中的地方病平衡点是唯一存在的（见图 1（b））。

3 平衡点的稳定性

定理2 如果，模型(\mathrm{2})的E^{\mathrm{0}}在\mathrm{\Omega }内是局部渐近稳定的。如果R_c>\mathrm{1}，E^{\mathrm{0}}是不稳定的。

证明 模型(\mathrm{2})在E^{\mathrm{0}}点处的Jacobian矩阵为

\mathit{J}(E^{\mathrm{0}})=\left[\begin{array}{cccccc}-(\mathrm{1}+{\theta }_{\mathrm{1}}+\sigma )& \beta e^{-m_{\mathrm{0}}}& \beta e^{-m_{\mathrm{0}}}\eta & \mathrm{0}& \mathrm{0}& \mathrm{0}\\ \sigma (\mathrm{1}-\rho )& -(\mathrm{1}+{\mu }_{\mathrm{1}}+\alpha +p)& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}\\ \sigma \rho & \mathrm{0}& -(\mathrm{1}+{\omega }_{\mathrm{1}}+\lambda )& \mathrm{0}& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& \alpha & \lambda & -(\mathrm{1}+{\omega }_{\mathrm{2}})& \mathrm{0}& \mathrm{0}\\ {\theta }_{\mathrm{1}}& p& {\omega }_{\mathrm{1}}& {\omega }_{\mathrm{2}}& -({\theta }_{\mathrm{2}}+\mathrm{1})& \mathrm{0}\\ \mathrm{0}& -{\mu }_{\mathrm{1}}& \mathrm{0}& \mathrm{0}& \mathrm{0}& -\mathrm{1}\end{array}\right]\mathrm{。}

令A^{\text{'}}=\mathrm{1}+{\theta }_{\mathrm{1}}+\sigma，B^{\text{'}}=\mathrm{1}+{\mu }_{\mathrm{1}}+\alpha +p，C^{\text{'}}=\mathrm{1}+{\omega }_{\mathrm{1}}+\lambda，则\mathit{J}(E^{\mathrm{0}})对应的特征方程为

\begin{array}{c}(\lambda +\mathrm{1})[\lambda +({\theta }_{\mathrm{2}}+\mathrm{1})][\lambda +(\mathrm{1}+{\omega }_{\mathrm{2}})]·\\ \left\{{\lambda }^{\mathrm{3}}+{\lambda }^{\mathrm{2}}(A^{\text{'}}+B^{\text{'}}+C^{\text{'}})\right.\lambda \left[A^{\text{'}}{B}^{\text{'}}+B^{\text{'}}{C}^{\text{'}}+A^{\text{'}}{C}^{\text{'}}-\beta {\mathrm{e}}^{-m_{\mathrm{0}}}\eta \sigma \rho -\beta {\mathrm{e}}^{-m_{\mathrm{0}}}\sigma (\mathrm{1}-\rho )\right]-\\ \left.\beta {\mathrm{e}}^{-m_{\mathrm{0}}}\eta \sigma \rho B^{\text{'}}-\beta {\mathrm{e}}^{-m_{\mathrm{0}}}\sigma (\mathrm{1}-\rho )C^{\text{'}}+A^{\text{'}}{B}^{\text{'}}{C}^{\text{'}}\right\}=\mathrm{0}\mathrm{。}\end{array}（4）

通过求解方程(\mathrm{4})，可得\mathit{J}(E^{\mathrm{0}})的特征值，即{\lambda }_{\mathrm{1}}=-\mathrm{1}，{\lambda }_{\mathrm{2}}=-(\mathrm{1}+{\theta }_{\mathrm{2}})，{\lambda }_{\mathrm{3}}=-(\mathrm{1}+{\omega }_{\mathrm{2}})。其他根的特征方程可简记为{\lambda }^{\mathrm{3}}+n_{\mathrm{1}}{\lambda }^{\mathrm{2}}+n_{\mathrm{2}}\lambda +n_{\mathrm{3}}=\mathrm{0}的形式， 其中， n_{\mathrm{1}}=A^{\text{'}}+B^{\text{'}}+C^{\text{'}}， n_{\mathrm{2}}=A^{\text{'}}{B}^{\text{'}}+B^{\text{'}}{C}^{\text{'}}+A^{\text{'}}{C}^{\text{'}}-\beta {\mathrm{e}}^{-m_{\mathrm{0}}}\eta \sigma \rho -\beta {\mathrm{e}}^{-m_{\mathrm{0}}}\sigma (\mathrm{1}-\rho )， n_{\mathrm{3}}=-\beta {\mathrm{e}}^{-m_{\mathrm{0}}}\eta \sigma \rho B^{\text{'}}－\beta {\mathrm{e}}^{-m_{\mathrm{0}}}\sigma (\mathrm{1}-\rho )C^{\text{'}}+A^{\text{'}}{B}^{\text{'}}{C}^{\text{'}}。

当时，有\beta {\mathrm{e}}^{-m_{\mathrm{0}}}\eta \sigma \rho B^{\text{'}}+\beta {\mathrm{e}}^{-m_{\mathrm{0}}}\sigma (\mathrm{1}-\rho )·C^{\text{'}}<A^{\text{'}}{B}^{\text{'}}{C}^{\text{'}}，，A^{\text{'}}{C}^{\text{'}}，所以n_{\mathrm{1}}>\mathrm{0}，n_{\mathrm{2}}>\mathrm{0}，n_{\mathrm{3}}>\mathrm{0}\mathrm{。}经计算可知，

\begin{array}{c}{n}_{\mathrm{1}}{n}_{\mathrm{2}}-n_{\mathrm{3}}={A^{\text{'}}}^{\mathrm{2}}{B}^{\text{'}}+\mathrm{2}{A}^{\text{'}}{B}^{\text{'}}{C}^{\text{'}}+{A^{\text{'}}}^{\mathrm{2}}{C}^{\text{'}}+\\ A^{\text{'}}{B^{\text{'}}}^{\mathrm{2}}+{B^{\text{'}}}^{\mathrm{2}}{C}^{\text{'}}+B^{\text{'}}{C^{\text{'}}}^{\mathrm{2}}+A^{\text{'}}{C^{\text{'}}}^{\mathrm{2}}-\end{array}

\begin{array}{c}{A}^{\text{'}}\beta e^{-m_{\mathrm{0}}}\eta \sigma \rho -A^{\text{'}}\beta e^{-m_{\mathrm{0}}}\sigma (\mathrm{1}-\rho )-\\ B^{\text{'}}\beta e^{-m_{\mathrm{0}}}\sigma (\mathrm{1}-\rho )-C^{\text{'}}\beta e^{-m_{\mathrm{0}}}\eta \sigma \rho >\mathrm{0}\mathrm{。}\end{array}

根据Routh-Hurwitz判据^［18］可知，当时，模型(\mathrm{2})的无病平衡点E^{\mathrm{0}}是局部渐近稳定的。当R_c>\mathrm{1}时，无病平衡点E^{\mathrm{0}}是不稳定的。

定理3 如果R_c>\mathrm{1}，模型(\mathrm{2})的地方病平衡点\overline{E}是局部渐近稳定的。

由于通过Routh-Hurwitz判据难以证明地方病平衡点\overline{E}是局部渐近稳定的，本文结合R_c=\mathrm{1}时有跨临界分支现象发生来证明。

证明 模型(\mathrm{2})在E^{\mathrm{0}}处的Jacobian矩阵为\mathit{J}(E^{\mathrm{0}})，基于中心流行理论，以\beta {\mathrm{e}}^{-m_{\mathrm{0}}}为分叉参数^［19］创建了地方病平衡点的局部稳定性。当R_c=\mathrm{1}时，分叉参数\beta {\mathrm{e}}^{-m_{\mathrm{0}}}的临界值为

{\beta }_c{\mathrm{e}}^{-m_{\mathrm{0}}}=\frac{(\mathrm{1}+{\theta }_{\mathrm{1}}+\sigma )(\mathrm{1}+{\mu }_{\mathrm{1}}+\alpha +p)(\mathrm{1}+{\omega }_{\mathrm{1}}+\lambda )}{\sigma (\mathrm{1}-\rho )(\mathrm{1}+{\omega }_{\mathrm{1}}+\lambda )+\eta \sigma \rho (\mathrm{1}+{\mu }_{\mathrm{1}}+\alpha +p)}。

可以验证在\beta {\mathrm{e}}^{-m_{\mathrm{0}}}={\beta }_c{\mathrm{e}}^{-m_{\mathrm{0}}}处的Jacobian矩阵\mathit{J}(E^{\mathrm{0}})的一个右特征向量为

\mathit{Y}=(y_{\mathrm{1}},y_{\mathrm{2}},y_{\mathrm{3}},y_{\mathrm{4}},y_{\mathrm{5}},y_{\mathrm{6}}{)}^{\mathrm{T}}，

其中， y_{\mathrm{1}}=\mathrm{1}+{\omega }_{\mathrm{1}}+\lambda，

y_{\mathrm{2}}=\frac{\sigma (\mathrm{1}-\rho )(\mathrm{1}+{\omega }_{\mathrm{1}}+\lambda )}{\mathrm{1}+{\mu }_{\mathrm{1}}+\alpha +p}， y_{\mathrm{3}}=\sigma \rho，

y_{\mathrm{4}}=\frac{\alpha \sigma (\mathrm{1}-\rho )(\mathrm{1}+{\omega }_{\mathrm{1}}+\lambda )+\lambda \sigma \rho (\mathrm{1}+{\mu }_{\mathrm{1}}+\alpha +p)}{(\mathrm{1}+{\mu }_{\mathrm{1}}+\alpha +p)(\mathrm{1}+{\omega }_{\mathrm{2}})},

\begin{array}{c}{y}_{\mathrm{5}}=\frac{{\theta }_{\mathrm{1}}(\mathrm{1}+{\omega }_{\mathrm{1}}+\lambda )(\mathrm{1}+{\mu }_{\mathrm{1}}+\alpha +p)}{(\mathrm{1}+{\mu }_{\mathrm{1}}+\alpha +p)(\mathrm{1}+{\theta }_{\mathrm{2}})}+\\ \frac{p\sigma (\mathrm{1}-\rho )(\mathrm{1}+{\omega }_{\mathrm{1}}+\lambda )}{(\mathrm{1}+{\mu }_{\mathrm{1}}+\alpha +p)(\mathrm{1}+{\theta }_{\mathrm{2}})}+\\ \frac{{\omega }_{\mathrm{1}}(\mathrm{1}+{\mu }_{\mathrm{1}}+\alpha +p)\sigma \rho }{(\mathrm{1}+{\mu }_{\mathrm{1}}+\alpha +p)(\mathrm{1}+{\theta }_{\mathrm{2}})}+\\ \frac{{\omega }_{\mathrm{2}}(\mathrm{1}+{\mu }_{\mathrm{1}}+\alpha +p)y_{\mathrm{4}}}{(\mathrm{1}+{\mu }_{\mathrm{1}}+\alpha +p)(\mathrm{1}+{\theta }_{\mathrm{2}})},\end{array}

y_{\mathrm{6}}=\frac{-{\mu }_{\mathrm{1}}\sigma (\mathrm{1}-\rho )(\mathrm{1}+{\omega }_{\mathrm{1}}+\lambda )}{\mathrm{1}+{\mu }_{\mathrm{1}}+\alpha +p}。

此外，对应于零特征值，左特征向量的分量\mathit{V}=(v_{\mathrm{1}},v_{\mathrm{2}},v_{\mathrm{3}},v_{\mathrm{4}},v_{\mathrm{5}},v_{\mathrm{6}})必须满足等式\mathit{V}\cdot \mathit{J}(E^{\mathrm{0}})=\mathrm{0}和\mathit{V}\cdot \mathit{Y}=\mathrm{1}。因此，v_{\mathrm{4}}=v_{\mathrm{5}}=v_{\mathrm{6}}=\mathrm{0}，

\begin{array}{c}{v}_{\mathrm{1}}=\\ \frac{v_{\mathrm{3}}[\sigma (\mathrm{1}-\rho )(\mathrm{1}+{\omega }_{\mathrm{1}}+\lambda )+\eta \sigma \rho (\mathrm{1}+{\mu }_{\mathrm{1}}+\alpha +p)]}{\eta (\mathrm{1}+{\theta }_{\mathrm{1}}+\sigma )(\mathrm{1}+{\mu }_{\mathrm{1}}+\alpha +p)},\end{array}

v_{\mathrm{2}}=\frac{v_{\mathrm{3}}(\mathrm{1}+{\omega }_{\mathrm{1}}+\lambda )}{\eta (\mathrm{1}+{\mu }_{\mathrm{1}}+\alpha +p)}， v_{\mathrm{3}}=\frac{\mathrm{1}}{y_{\mathrm{3}}+B_{\mathrm{1}}+B_{\mathrm{2}}}。

其中，B_{\mathrm{1}}=\frac{\sigma (\mathrm{1}-\rho )(\mathrm{1}+{\omega }_{\mathrm{1}}{+\lambda )}^{\mathrm{2}}}{\eta (\mathrm{1}+{\mu }_{\mathrm{1}}{+\alpha +p)}^{\mathrm{2}}}，

B_{\mathrm{2}}=\frac{(\mathrm{1}+{\omega }_{\mathrm{1}}+\lambda )[\sigma (\mathrm{1}-\rho )(\mathrm{1}+{\omega }_{\mathrm{1}}+\lambda )}{\eta (\mathrm{1}+{\theta }_{\mathrm{1}}+\sigma )(\mathrm{1}+{\mu }_{\mathrm{1}}+\alpha +p)}+

\frac{\eta \sigma \rho (\mathrm{1}+{\mu }_{\mathrm{1}}+\alpha +p)]}{\eta (\mathrm{1}+{\theta }_{\mathrm{1}}+\sigma )(\mathrm{1}+{\mu }_{\mathrm{1}}+\alpha +p)}\mathrm{。}

使用符号x_{\mathrm{1}}\equiv E，x_{\mathrm{2}}\equiv I，x_{\mathrm{3}}\equiv A，x_{\mathrm{4}}\equiv Q，x_{\mathrm{5}}\equiv R，x_{\mathrm{6}}\equiv N，得

a={\displaystyle \sum _{k,i,j=\mathrm{1}}^{\mathrm{6}}}{v}_k{y}_i{y}_j\frac{{\partial }^{\mathrm{2}}{f}_k(\mathrm{0,0})}{\partial x_i\partial x_j},

b={\displaystyle \sum _{k,i,j=\mathrm{1}}^{\mathrm{6}}}{v}_k{y}_i\frac{{\partial }^{\mathrm{2}}{f}_k(\mathrm{0,0})}{\partial x_i\partial \beta {\mathrm{e}}^{-m_{\mathrm{0}}}}。

当\beta {\mathrm{e}}^{-m_{\mathrm{0}}}={\beta }_c{\mathrm{e}}^{-m_{\mathrm{0}}}时，将无病平衡点E^{\mathrm{0}}处计算的所有二阶导数的值代入等式a和b，得

\begin{array}{c}a=v_{\mathrm{1}}[\mathrm{2}{y}_{\mathrm{1}}{y}_{\mathrm{6}}-\mathrm{2}{y}_{\mathrm{2}}{y}_{\mathrm{4}}{\beta }_c{\mathrm{e}}^{-m_{\mathrm{0}}}-\mathrm{2}{y}_{\mathrm{2}}{y}_{\mathrm{5}}{\beta }_c{\mathrm{e}}^{-m_{\mathrm{0}}}-\\ \mathrm{2}{y}_{\mathrm{3}}{y}_{\mathrm{4}}\eta {\beta }_c{\mathrm{e}}^{-m_{\mathrm{0}}}-\mathrm{2}{y}_{\mathrm{3}}{y}_{\mathrm{5}}\eta {\beta }_c{\mathrm{e}}^{-m_{\mathrm{0}}}+\mathrm{2}{y}_{\mathrm{2}}{y}_{\mathrm{3}}·\\ (-{\beta }_c{\mathrm{e}}^{-m_{\mathrm{0}}}-m_{\mathrm{2}}{\beta }_c{\mathrm{e}}^{-m_{\mathrm{0}}}-\eta {\beta }_c{\mathrm{e}}^{-m_{\mathrm{0}}}-m_{\mathrm{1}}\eta {\beta }_c{\mathrm{e}}^{-m_{\mathrm{0}}})+\end{array}

\begin{array}{c}{y}_{\mathrm{2}}^{\mathrm{2}}(-\mathrm{2}{\beta }_c{e}^{-m_{\mathrm{0}}}-\mathrm{2}{m}_{\mathrm{1}}{\beta }_c{e}^{-m_{\mathrm{0}}})+\mathrm{2}{y}_{\mathrm{1}}{y}_{\mathrm{2}}(-{\beta }_c{\mathrm{e}}^{-m_{\mathrm{0}}}+\\ {\mu }_{\mathrm{1}})+\mathrm{2}{y}_{\mathrm{1}}{y}_{\mathrm{3}}(-\eta {\beta }_c{\mathrm{e}}^{-m_{\mathrm{0}}})+y_{\mathrm{3}}^{\mathrm{2}}(-\mathrm{2}\eta {\beta }_c{\mathrm{e}}^{-m_{\mathrm{0}}}-\\ \mathrm{2}{m}_{\mathrm{2}}\eta {\beta }_c{\mathrm{e}}^{-m_{\mathrm{0}}})]+v_{\mathrm{2}}(\mathrm{2}{y}_{\mathrm{2}}{y}_{\mathrm{6}}+\mathrm{2}{y}_{\mathrm{2}}^{\mathrm{2}}{\mu }_{\mathrm{1}})+\\ v_{\mathrm{3}}(\mathrm{2}{y}_{\mathrm{3}}{y}_{\mathrm{6}}+\mathrm{2}{y}_{\mathrm{2}}{y}_{\mathrm{3}}{\mu }_{\mathrm{1}})+v_{\mathrm{4}}(\mathrm{2}{y}_{\mathrm{4}}{y}_{\mathrm{6}}+\mathrm{2}{y}_{\mathrm{2}}{y}_{\mathrm{4}}{\mu }_{\mathrm{1}})+\\ v_{\mathrm{5}}(\mathrm{2}{y}_{\mathrm{5}}{y}_{\mathrm{6}}+\mathrm{2}{y}_{\mathrm{2}}{y}_{\mathrm{5}}{\mu }_{\mathrm{1}})+v_{\mathrm{6}}(-\mathrm{2}{y}_{\mathrm{3}}{y}_{\mathrm{6}}{\mu }_{\mathrm{1}}),\end{array}

b=v_{\mathrm{1}}(y_{\mathrm{2}}+\eta y_{\mathrm{3}})。

最后，代入\mathit{V}和\mathit{Y}的值，得

a=B_{\mathrm{3}}+B_{\mathrm{4}}+B_{\mathrm{5}}， b=\frac{v_{\mathrm{3}}{B}_{\mathrm{2}}^{\mathrm{2}}\eta (\mathrm{1}+{\theta }_{\mathrm{1}}+\sigma )}{(\mathrm{1}+{\omega }_{\mathrm{2}}{+\lambda )}^{\mathrm{2}}}，

其中，

\begin{array}{c}{B}_{\mathrm{4}}=v_{\mathrm{1}}[-\mathrm{2}{y}_{\mathrm{2}}{y}_{\mathrm{4}}{\beta }_c{\mathrm{e}}^{-m_{\mathrm{0}}}-\mathrm{2}{y}_{\mathrm{2}}{y}_{\mathrm{5}}{\beta }_c{\mathrm{e}}^{-m_{\mathrm{0}}}-\\ \mathrm{2}{y}_{\mathrm{3}}{y}_{\mathrm{4}}\eta {\beta }_c{\mathrm{e}}^{-m_{\mathrm{0}}}-\mathrm{2}{y}_{\mathrm{3}}{y}_{\mathrm{5}}\eta {\beta }_c{\mathrm{e}}^{-m_{\mathrm{0}}}+\mathrm{2}{y}_{\mathrm{2}}{y}_{\mathrm{3}}(-{\beta }_c{\mathrm{e}}^{-m_{\mathrm{0}}}-\\ m_{\mathrm{2}}{\beta }_c{\mathrm{e}}^{-m_{\mathrm{0}}}-\eta {\beta }_c{\mathrm{e}}^{-m_{\mathrm{0}}}-m_{\mathrm{1}}\eta {\beta }_c{\mathrm{e}}^{-m_{\mathrm{0}}})+\\ y_{\mathrm{2}}^{\mathrm{2}}(-\mathrm{2}{\beta }_c{\mathrm{e}}^{-m_{\mathrm{0}}}-\mathrm{2}{m}_{\mathrm{1}}{\beta }_c{\mathrm{e}}^{-m_{\mathrm{0}}})+\mathrm{2}{y}_{\mathrm{1}}{y}_{\mathrm{3}}(-\eta {\beta }_c{\mathrm{e}}^{-m_{\mathrm{0}}})+\end{array}

\begin{array}{c}{B}_{\mathrm{5}}=v_{\mathrm{2}}(\mathrm{2}{y}_{\mathrm{2}}{y}_{\mathrm{6}}+\mathrm{2}{y}_{\mathrm{2}}^{\mathrm{2}}{\mu }_{\mathrm{1}})+v_{\mathrm{3}}(\mathrm{2}{y}_{\mathrm{3}}{y}_{\mathrm{6}}+\mathrm{2}{y}_{\mathrm{2}}{y}_{\mathrm{3}}{\mu }_{\mathrm{1}})+\\ v_{\mathrm{4}}(\mathrm{2}{y}_{\mathrm{4}}{y}_{\mathrm{6}}+\mathrm{2}{y}_{\mathrm{2}}{y}_{\mathrm{4}}{\mu }_{\mathrm{1}})+v_{\mathrm{5}}(\mathrm{2}{y}_{\mathrm{5}}{y}_{\mathrm{6}}+\mathrm{2}{y}_{\mathrm{2}}{y}_{\mathrm{5}}{\mu }_{\mathrm{1}})+\\ v_{\mathrm{6}}(-\mathrm{2}{y}_{\mathrm{3}}{y}_{\mathrm{6}}{\mu }_{\mathrm{1}}）=\mathrm{0}\mathrm{。}\end{array}

由于\beta {\mathrm{e}}^{-m_{\mathrm{0}}}={\beta }_c{\mathrm{e}}^{-m_{\mathrm{0}}}时，且b>\mathrm{0}，则在R_c=\mathrm{1}时有跨临界分支现象发生，并且对于R_c>\mathrm{1}，唯一地方病平衡点\overline{E}是局部渐近稳定的。

定理4 如果且{\mu }_{\mathrm{1}}=\mathrm{0}，则模型(\mathrm{2})的无病平衡点E^{\mathrm{0}}是全局渐近稳定的。

证明 当{\mu }_{\mathrm{1}}=\mathrm{0}时，得到\frac{\mathrm{d}N}{\mathrm{d}t}=\mathrm{1}-N。当t\to \mathrm{\infty }时，N\to \mathrm{1}。可在极限条件下得N=\mathrm{1}，则模型(\mathrm{2})可被改写为

\left\{\begin{array}{l}\frac{\mathrm{d}E}{\mathrm{d}t}=\beta {\mathrm{e}}^{-m_{\mathrm{1}}I-m_{\mathrm{2}}A-m_{\mathrm{0}}}S(I+\eta A)-(\mathrm{1}+{\theta }_{\mathrm{1}}+\sigma )E,\\ \frac{\mathrm{d}I}{\mathrm{d}t}=\sigma (\mathrm{1}-\rho )E-(\mathrm{1}+\alpha +p)I,\\ \frac{\mathrm{d}A}{\mathrm{d}t}=\sigma \rho E-(\mathrm{1}+{\omega }_{\mathrm{1}}+\lambda )A,\\ \frac{\mathrm{d}Q}{\mathrm{d}t}=\alpha I+\lambda A-(\mathrm{1}+{\omega }_{\mathrm{2}})Q,\\ \frac{\mathrm{d}R}{\mathrm{d}t}={\omega }_{\mathrm{1}}A+{\omega }_{\mathrm{2}}Q+pI+{\theta }_{\mathrm{1}}E-{\theta }_{\mathrm{2}}R-R\mathrm{。}\end{array}\right.（5）

设\mathit{X}=(R)，\mathit{Z}=(E,I,A,Q)，这里{\mathit{U}}_{\mathrm{0}}=（X^{\mathrm{0}}，Z^{\mathrm{0}}），其中{\mathit{X}}^{\mathrm{0}}=(\mathrm{0})，{\mathit{Z}}^{\mathrm{0}}=(\mathrm{0,0},\mathrm{0},\mathrm{0})， 得

\begin{array}{c}\frac{\mathrm{d}X}{\mathrm{d}t}=\mathit{F}(X,Z)={\omega }_{\mathrm{1}}A+{\omega }_{\mathrm{2}}Q+pI+\\ {\theta }_{\mathrm{1}}E-({\theta }_{\mathrm{2}}+\mathrm{1})R\mathrm{。}\end{array}

在\mathit{Z}={\mathit{Z}}^{\mathrm{0}}时，\mathit{G}(X,\mathrm{0})=\mathrm{0}。等式\frac{\mathrm{d}X}{\mathrm{d}t}=\mathit{F}(X,\mathrm{0})=-(\mathrm{1}+{\theta }_{\mathrm{2}})X，当t\to \mathrm{\infty }时，X\to X^{\mathrm{0}}。所以可知，X=X^{\mathrm{0}}(=R^{\mathrm{0}}=\mathrm{0})是全局渐近稳定的。从模型(\mathrm{5})可得

\frac{\mathrm{d}Z}{\mathrm{d}t}=\mathit{G}(X,Z)=\mathit{B}\mathit{Z}-\widehat{\mathit{G}}(X,Z)，

其中，

\mathit{B}=\left[\begin{array}{cccc}-\mathrm{1}-{\theta }_{\mathrm{1}}-\sigma & \beta {\mathrm{e}}^{-m_{\mathrm{0}}}& \beta {\mathrm{e}}^{-m_{\mathrm{0}}}\eta & \beta {\mathrm{e}}^{-m_{\mathrm{0}}}\\ \sigma (\mathrm{1}-\rho )& -\mathrm{1}-\alpha -p& \mathrm{0}& \mathrm{0}\\ \sigma \rho & \mathrm{0}& -\mathrm{1}-{\omega }_{\mathrm{1}}-\lambda & \mathrm{0}\\ \mathrm{0}& \alpha & \lambda & -\mathrm{1}-{\omega }_{\mathrm{2}}\end{array}\right],

\begin{array}{c}\widehat{\mathit{G}}(X,Z)=\\ \left[\begin{array}{c}\beta e^{-m_{\mathrm{0}}}(I+\eta A)(\mathrm{1}-{\mathrm{e}}^{-m_{\mathrm{1}}I-m_{\mathrm{2}}A}S)+\beta {\mathrm{e}}^{-m_{\mathrm{0}}}Q\\ \mathrm{0}\\ \mathrm{0}\\ \mathrm{0}\end{array}\right]\end{array}。

显然，矩阵\mathit{B}为 M 矩阵，且不可约。当I\ge \mathrm{0}，A\ge \mathrm{0}时，有。当\mathrm{0}\le S\le \mathrm{1}时，有\widehat{\mathit{G}}(X,Z)\ge \mathrm{0}，满足条件(H_{\mathrm{1}})和(H_{\mathrm{2}})。因此，如果，无病平衡点E^{\mathrm{0}}是全局渐近稳定的^［20］。

4 阈值分析

通过计算基本再生数R_c相对于参数\alpha（有症状感染者的隔离率）、参数p（有症状感染者的恢复率）和参数m_{\mathrm{0}}（媒体报道系数）的偏导数，定性测量了隔离、治疗和媒体报道对疾病传播动力学的影响。对参数\alpha、参数p和参数m_{\mathrm{0}}进行阈值分析，可得

通过讨论，可知隔离参数\alpha越大，隔离效果越好。隔离参数\alpha增大，基本再生数R_c反而减小，即对传染性个体隔离将降低人类接触率和传播概率。有症状感染者的恢复率p增大，基本再生数R_c减小。媒体报道系数m_{\mathrm{0}}增大，基本再生数R_c减小。

5 数值模拟

针对模型(\mathrm{2})的无病平衡点、地方病平衡点的稳定性问题，分别进行了数值模拟。选择初始数值{\tilde{S}}_{\mathrm{0}}=400， {\tilde{E}}_{\mathrm{0}}=60， {\tilde{I}}_{\mathrm{0}}=250，{\tilde{A}}_{\mathrm{0}}=100， {\tilde{Q}}_{\mathrm{0}}=40， {\tilde{R}}_{\mathrm{0}}=150。无量纲形式的相应初始数值为S_{\mathrm{0}}=0.04， E_{\mathrm{0}}=0.06， I_{\mathrm{0}}=0.25， A_{\mathrm{0}}=0.10， Q_{\mathrm{0}}=0.04， R_{\mathrm{0}}=0.15。以下部分参数选取自参考文献［15］。令\mathrm{\Lambda }=136， w_{\mathrm{1}}=\frac{\mathrm{1}}{\mathrm{6}}， w_{\mathrm{2}}=\frac{\mathrm{1}}{\mathrm{14}}，m_{\mathrm{1}}=0.2， p=\frac{\mathrm{1}}{\mathrm{100}}， \mu=0.075 1， m_{\mathrm{2}}=0.2， {\mu }_{\mathrm{1}}=0.035 21， \eta=0.6， \sigma=\frac{\mathrm{1}}{\mathrm{5.2}}， {\theta }_{\mathrm{1}}=0.5， \lambda=0.95， \alpha=\frac{\mathrm{1}}{\mathrm{14}}， \rho=0.74， {\theta }_{\mathrm{2}}=0.5， \beta {\mathrm{e}}^{-m_{\mathrm{0}}}=［0.25，0.95］。

对于\beta {\mathrm{e}}^{-m_{\mathrm{0}}}=\mathrm{0.25}，R_c=\mathrm{0.622}\mathrm{83}，无病平衡点E^{\mathrm{0}}在\mathrm{\Omega }内是全局渐近稳定的，如图 2（a）。对于\beta {\mathrm{e}}^{-m_{\mathrm{0}}}=\mathrm{0.95}，R_c=\mathrm{2.507}\mathrm{20}，地方病平衡点\overline{E}是局部渐近稳定的，如图 2（b）。图 3 分别显示了参数m_{\mathrm{1}}和m_{\mathrm{2}}对受感染个体的影响，这表明媒体报道降低了感染率。当和R_c=\mathrm{2.507}\mathrm{20}>\mathrm{1}的初始值不同时，有症状感染人数占总人数的比例随时间的函数关系分别如图 4（a） 和 4（b） 所示。可以得出，当R_c为确定值时，有症状感染个体的趋势不会随着有症状感染初始值的不同而改变。从图 5 可以看出，基本再生数随着隔离参数\alpha、有症状感染者的恢复率参数p和媒体报道系数m_{\mathrm{0}}的增大而减少。参数\sigma对感染人数的影响如图 6（a） 所示，当时，感染人数减少并趋于0，而\sigma >\mathrm{0.21}时，感染人数增加。参数\alpha对感染人数的影响如图 6（b） 所示，随着隔离参数的增大，感染人数将减少。

6 结 论

考虑到媒体报道以及隔离对传染病传播的影响，本文构建并分析了一种含媒体报道和隔离的SEIAQRS传染病传播动力学模型。通过推导得出了传染病的传播阈值R_c，发现当时，无病平衡点E^{\mathrm{0}}全局渐近稳定。当R_c>\mathrm{1}时，地方病平衡点\overline{E}局部渐近稳定。数值模拟验证了理论分析的正确性。基本再生数R_c受媒体报道系数m_{\mathrm{0}}的影响，当m_{\mathrm{0}}增大时，基本再生数R_c减小。当疾病不存在时，I=\mathrm{0},A=\mathrm{0}，媒体报道系数m_{\mathrm{1}}和m_{\mathrm{2}}不影响基本再生数R_c。然而，如果疾病发生大规模流行，其将有助于减轻疾病负担，因为媒体报道系数m_{\mathrm{1}}和m_{\mathrm{2}}降低了地方病稳定状态下的感染人数。该模型适用于新型冠状病毒肺炎的传播，分析模型参数对基本再生数和感染人数的影响，可以从媒体报道、检疫、治疗和恢复等方面为新型冠状病毒肺炎的预防和治疗提供建议。

参考文献

[1] GUAN J X， WEI Y Y， ZHAO Y，et al.Modeling the transmission dynamics of COVID-19 epidemic：A systematic review［J］.Journal of Biomedical Research，2020，34（6）：422-430.

[2] TANG B， WANG X， LI Q，et al.Estimation of the transmission risk of the 2019-nCoV and its implication for public health interventions［J］.Journal of Clinical Medicine，2020，9（2）：462.

[3] DIN R， SEADAWY A R， SHAH K，et al.Study of global dynamics of COVID-19 via a new mathematical model［J］.Results in Physics，2020，19：103468.

[4] LIU R， WU J， ZHU H.Media/psychological impact on multiple outbreaks of emerging infectious diseases［J］.Computational and Mathematical Methods in Medicine，2007，8（3）：153-164.

[5] LI Y F， CUI J A. The effect constant and pulse vaccination on SIS epidemic models incorporate media coverage［J］. Communications in Nonlinear Science and Numerical Simulation，2009，14（5）：2353-2365.

[6] TCHUENCHE J M， DUBE N， BHUNU C P，et al.The impact of media coverage on the transmission dynamics of human influenza［J］.BMC Public Health，2011，11（1）：S5.

[7] RONG X M， YANG L， CHU H D，et al.Effect of delay in diagnosis on transmission of COVID-19［J］.Mathematical Biosciences and Engineering，2020，17（3）：2725-2740.

[8] WANG A， XIAO Y.A Filippov system describing media effects on the spread of infectious diseases［J］.Nonlinear Analysis：Hybrid Systems，2014，11：84-97.

[9] CUI J， SUN Y， ZHU H.The impact of media on the control of infectious diseases［J］.Journal of Dynamics and Differential Equations，2007，20：31-53.

[10] LIU Y P， CUI J A，The impact of media coverage on the dynamics of infectious disease［J］.International Journal of Biomathematics，2008，1（1）：65-74.

[11] SAHU G P， DHAR J，Analysis of an SVEIS epidemic model with partial temporary immuneity and saturation incidence rate［J］.Applied Mathematical Modelling，2012，36（3）：908-923.

[12] 武婧媛，两类流感疾病模型的研究［D］.临汾：山西师范大学，2016.

[13] CUI J， TAO X， ZHU H P.An SIS infection on model in corporating media coverage［J］.Rocky Moutain Journal of Mathematics，2008，38 （5）：1323-1334.

[14] HUO H F， YANG P， XIANG H.Stability and bifurcation for an SEIS epidemic model with the impact of media［J］.Physica A：Statistical Mechanics and Its Applications，2018，490：702-720.

[15] SAHU G P， DHAR J.Dynamics of an SEQIHRS epidemic model with media coverage，quarantine and isolation in a community with pre-existing immunity［J］.Journal of Mathematical Analysis and Applications，2015，421（2）：1651-1672.

[16] SHI X， GAO X， ZHOU X，et al.Analysis of an SQEIAR epidemic model with media coverage and asymptomatic infection［J］.AIMS Mathematics，2021，6（11）：12298-12320.

[17] DRIESSCHE P， WATMOUGH J.Reproduction numbers and sub-threshold endemic equilibria for compartmental models of disease transmission［J］.Mathematical Biosciences，2002，180 （1/2）：29-48.

[18] 王琪，薛亚奎.具有心理效应的媒介传染病模型的研究［J］.重庆理工大学学报（自然科学版），2021，35（2）：251-257.

[19] WANG Qi， XUE Yakui.Study on the vector-borne infectious disease model with psychological effect［J］.Journal of Chongqing University of Technology（Natural Science），2021，35（2）：251-257. （in Chinese）

[20] CASTILLO-CHAVEZ C， SONG B.Dynamical models of tuberculosis and their applications［J］.Mathematical Biosciences Engineering，2004，1（2）：361-404.

[21] CASTILLO-CHAVEZ C， FENG Z L， HUANG W Z.Mathematical approaches for emerging and reemerging infectious diseases：An introduction［M］.New York： Springer-Verlag，2002.

基金资助

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

山西省自然科学青年基金资助项目(201801D221040)