有限医疗资源影响下的多尺度传染病模型研究

Analysis of Multi-scale Infectious Disease Model Under the Influence of Limited Medical Resources

摘要

考虑到传染病传播过程中宿主体内病毒载量和医疗资源有限对疾病传播的影响，提出一类具有饱和治疗率的宿主体内和宿主之间耦合的多尺度传染病模型。首先，对于宿主之间疾病传播模型，通过计算得到基本再生数，并给出刻画其无病平衡点和地方病平衡点的存在性与稳定性的判别准则。特别地，当基本再生数小于1时，模型存在两个地方病平衡点，其中一个局部渐近稳定，另一个不稳定。当基本再生数大于1时，模型存在唯一的地方病平衡点，它在特定条件下是全局渐近稳定的。其次，探讨医疗资源有限导致的前向/后向分支的存在性，这意味着疾病的消除或流行不再取决于基本再生数，而是与感染者的初始状态以及医疗资源的供应密切相关。最后，通过数值模拟解释了本文的主要理论结果。

关键词

多尺度传染病模型 / 基本再生数 / 平衡点 / 渐近稳定性 / 前向/后向分支

Key words

multi-scale infectious disease model / basic reproduction number / equilibrium / asymptotic stability / forward/backward bifurcation

0 引言

传染病严重危害着人类的生命健康，并给人们的工作和生活带来诸多不便。近年来，一些由病毒引起的传染病，例如埃博拉出血热（Ebola Hemorrhagic Fever，EBHF）^［1］，艾滋病（Acquired Immune Deficiency Syndrome，AIDS）^［2］，新型冠状病毒感染（Coronavirus Disease 2019，COVID-19）^［3-4］等，仍然是众多学者的研究热点。通常情况下，病毒性传染病的传播至少涉及两个关键过程，分别是宿主体内病毒感染过程和宿主之间疾病传播过程^［5］，将两个过程结合起来考虑对预防和控制病毒性传染病具有重要意义。为此，一些国内外学者建立了宿主体内和宿主之间的疾病传播耦合模型^［6-11］，讨论了宿主体内病毒浓度对疾病在宿主之间传播的影响。特别地，Gilchrist等^［9］将宿主体内病毒感染模型嵌套在宿主之间疾病传播模型中，推导出了疾病传播与毒力之间的函数关系。Feng等^［10］在文献［9］的基础上引入依赖于宿主体内病毒浓度的传播率，使得宿主体内模型和宿主之间模型存在明显的耦合关系，结果表明了耦合模型更符合传染病的传播规律。

2019年末，新型冠状病毒感染疫情的突然暴发和短期内病例的迅速增加，使得各国的医疗资源如防护服、口罩、医院床位等十分紧缺。在面对重大公共卫生事件时，医疗资源短缺无疑会增加疫情防控难度，甚至会耽误对患者的治疗。有限的医疗资源对传染病传播和控制的影响引起了国内外学者的广泛关注，例如，Yadav等^［12］提出了具有饱和治疗率的SIRS（Susceptible-infected-recovered-susceptible）模型，分析了由于医疗资源不足而导致的治疗率饱和对疾病动态的影响。Hu等^［13］建立了两斑块数学模型来模拟SARS-CoV-2（Severe Acute Respiratory Syndrome Coronavirus 2）的传播，并引入半饱和检测率函数来描述医疗资源有限，结果表明医疗资源的供应对控制疾病传播具有重要作用。

基于上述讨论，考虑到病毒性传染病传播过程中的复杂性和医疗资源的有限性，本文将建立一类具有饱和治疗率的多尺度疾病传播耦合模型，讨论宿主体内病毒载量和有限医疗资源如何影响疾病在宿主之间的传播。

1 模型的建立

将某地区总人口分为三类：易感者类、感染者类和恢复者类，其在 t时刻的数量分别为 S \left(t\right)， I \left(t\right)和 R \left(t\right)。考虑到病毒载量对疾病传播的影响^［14］，即疾病传播率与感染宿主体内的病毒浓度 V有关，不妨设传播率为 \beta \left(V\right) ^［10］。此外，假设感染力随着感染者数量的增加而饱和， \alpha决定了感染力饱和的水平。因此，病毒在宿主之间的传播可由下面模型（1）的微分方程组描述

\left\{\begin{array}{l}\frac{\mathrm{d}S\left(t\right)}{\mathrm{d}t}={\Lambda }_h-\frac{\beta \left(V\right)S\left(t\right)I\left(t\right)}{\mathrm{1}+\alpha I\left(t\right)}-\mu S\left(t\right)+\kappa R\left(t\right),\\ \frac{\mathrm{d}I\left(t\right)}{\mathrm{d}t}=\frac{\beta \left(V\right)S\left(t\right)I\left(t\right)}{\mathrm{1}+\alpha I\left(t\right)}-\frac{\gamma I\left(t\right)}{\mathrm{1}+\theta I\left(t\right)}-\mu I\left(t\right),\\ \frac{\mathrm{d}R\left(t\right)}{\mathrm{d}t}=\frac{\gamma I\left(t\right)}{\mathrm{1}+\theta I\left(t\right)}-\left(\mu +\kappa \right)R\left(t\right),\end{array}\right. （1）

其中 {\Lambda }_h表示易感者的补充率， \mu表示宿主（包括易感者，感染者和恢复者）的自然死亡率， \kappa表示恢复者丧失免疫力重新变得易感的速率。饱和治疗函数 \gamma I / ( \mathrm{1} + \theta I )描述了医疗资源的有限性，这里 \gamma表示治疗率， \theta是饱和系数， \theta = \mathrm{0}表示所有感染者都可以接受治疗， \theta越大表示医疗资源越有限。

在上述模型中，假设 \beta \left(V\right)满足 \beta \left(\mathrm{0}\right) = \mathrm{0} , \beta \left(V\right) \ge \mathrm{0} , \beta \mathrm{\text{'}} \left(V\right) > \mathrm{0} , \beta \mathrm{\text{'}} \mathrm{\text{'}} \left(V\right) \le \mathrm{0}。 \beta \left(V\right)的最简单形式是 \beta \left(V\right) = a V，其他形式包括 \beta \left(V\right) = a V / \left(\mathrm{1}+bV\right)和 ，其中a和b是正常数。

为了得到V的具体形式，需要进一步研究病毒在宿主体内的感染过程。基于文献［10］，病毒在宿主体内的演变规律由下面模型（2）的微分方程组描述

\left\{\begin{array}{l}\frac{\mathrm{d}T\left(s\right)}{\mathrm{d}s}={\Lambda }_c-{\kappa }_{\mathrm{1}}T\left(s\right)V\left(s\right)-{\kappa }_{\mathrm{2}}T\left(s\right)T^{\mathrm{*}}\left(s\right)-{\mu }_{c}T\left(s\right),\\ \frac{\mathrm{d}{T}^{\mathrm{*}}\left(s\right)}{\mathrm{d}s}={\kappa }_{\mathrm{1}}T\left(s\right)V\left(s\right)+{\kappa }_{\mathrm{2}}T\left(s\right)T^{\mathrm{*}}\left(s\right)-\left({\mu }_c+{\delta }_c\right)T^{\mathrm{*}}\left(s\right),\\ \frac{\mathrm{d}V\left(s\right)}{\mathrm{d}s}=p{T}^{\mathrm{*}}\left(s\right)-{\mu }_{v}V\left(s\right),\end{array}\right. （2）

其中 T \left(s\right) , T^{\mathrm{*}} \left(s\right)和 V \left(s\right)分别表示健康体细胞、感染体细胞和游离病毒在时刻s的浓度。 {\Lambda }_c表示健康体细胞的补充率， {\kappa }_{\mathrm{1}}表示病毒的入侵感染率， {\kappa }_{\mathrm{2}}表示感染体细胞的接触感染率， {\mu }_c表示细胞（健康体细胞，感染体细胞）的自然死亡率， {\delta }_c表示感染体细胞的因病死亡率， p表示感染体细胞的病毒释放率， {\mu }_v表示体内病毒的清除率。

在模型（1）中，t表示病毒在宿主之间传播的时间尺度，在模型（2）中，s表示体细胞，病毒在宿主体内演化的时间尺度。通常情况下，病毒在宿主体内的演化速率要高于其在宿主之间的传播速率，因此t是一个慢时间尺度，而s是一个快时间尺度^［11］。不妨假设 t = ϵ s，其中 ϵ是一个非常小的正常数。

2 宿主体内病毒感染模型

根据模型（2）的生物学意义，假设它的任意解 \left(T\left(s\right),T^{\mathrm{*}}\left(s\right),V\left(s\right)\right)满足初始条件

T \left(\mathrm{0}\right) > \mathrm{0} , T^{\mathrm{*}} \left(\mathrm{0}\right) > \mathrm{0} , V \left(\mathrm{0}\right) > \mathrm{0} \mathrm{。} （3）

关于模型（2）解的正性和有界性，下面的定理显然成立。

定理1 模型（2）满足初始条件（3）的解 \left(T\left(s\right),T^{\mathrm{*}}\left(s\right),V\left(s\right)\right)对任何 s \ge \mathrm{0}是正的和最终有界的。

定义宿主体内模型（2）的基本再生数为

{\mathrm{ℛ}}_{\mathrm{0}}^w = \frac{\left({\kappa }_{\mathrm{1}}p+{\kappa }_{\mathrm{2}}{\mu }_v\right)T_{\mathrm{0}}}{{\mu }_v\left({\mu }_c+{\delta }_c\right)} ,

其中 T_{\mathrm{0}} = {\Lambda }_c / {\mu }_c。关于模型（2）的动力学性质，在文献［15］中已有相应结论，现总结如下。

关于模型（2）平衡点的存在性，有下面的结论。

定理2 若 {\mathrm{ℛ}}_{\mathrm{0}}^w \le \mathrm{1}，则模型（2）仅存在无感染平衡点 U_{\mathrm{0}} = \left(T_{\mathrm{0}},\mathrm{0},\mathrm{0}\right)；若 {\mathrm{ℛ}}_{\mathrm{0}}^w > \mathrm{1}，则模型（2）除 U_{\mathrm{0}}外还存在唯一的感染平衡点 \tilde{U} = \left(\tilde{T},{\tilde{T}}^{\mathrm{*}},\tilde{V}\right)，这里 \tilde{T} = \frac{T_{\mathrm{0}}}{\mathrm{ℛ}{}_{\mathrm{0}}^w} , {\tilde{T}}^{\mathrm{*}} = \frac{{\mu }_v{\mu }_c\left(\mathrm{ℛ}{}_{\mathrm{0}}^w-\mathrm{1}\right)}{{\kappa }_{\mathrm{1}}p+{\kappa }_{\mathrm{2}}{\mu }_v} , \tilde{V} = \frac{p{\mu }_c\left(\mathrm{ℛ}{}_{\mathrm{0}}^w-\mathrm{1}\right)}{{\kappa }_{\mathrm{1}}p+{\kappa }_{\mathrm{2}}{\mu }_v} \mathrm{。}

关于 U_{\mathrm{0}}和 \tilde{U}的全局稳定性，有下面的结论。

定理3 若 {\mathrm{ℛ}}_{\mathrm{0}}^w \le \mathrm{1}，则无感染平衡点 U_{\mathrm{0}}全局渐近稳定；若 {\mathrm{ℛ}}_{\mathrm{0}}^w > \mathrm{1}，则感染平衡点 \tilde{U}全局渐近稳定。

3 宿主之间SIRS耦合模型

由于病毒在宿主体内的感染变化速度高于病毒在宿主之间的传播速度，因此，假设宿主之间的疾病传播过程还没有发生显著变化时，宿主体内的病毒载量 V已经达到其稳态 \widehat{V}。即当 {\mathrm{ℛ}}_{\mathrm{0}}^w \le \mathrm{1}时， \widehat{V} = \mathrm{0}；当 {\mathrm{ℛ}}_{\mathrm{0}}^w > \mathrm{1}时， \widehat{V} = \tilde{V}。当 {\mathrm{ℛ}}_{\mathrm{0}}^w \le \mathrm{1}时，感染者体内的病毒数量最终被清除，则由病毒引起的传染病最终不流行。因此，在接下来的讨论中，我们总是假设 {\mathrm{ℛ}}_{\mathrm{0}}^w > \mathrm{1}。如果 {\mathrm{ℛ}}_{\mathrm{0}}^w > \mathrm{1}，则 \beta \left(\widehat{V}\right) = \beta \left(\tilde{V}\right)，从而模型（1）可改写为

\left\{\begin{array}{l}\frac{\mathrm{d}S\left(t\right)}{\mathrm{d}t}={\Lambda }_h-\frac{\beta \left(\tilde{V}\right)S\left(t\right)I\left(t\right)}{\mathrm{1}+\alpha I\left(t\right)}-\mu S\left(t\right)+\kappa R\left(t\right),\\ \frac{\mathrm{d}I\left(t\right)}{\mathrm{d}t}=\frac{\beta \left(\tilde{V}\right)S\left(t\right)I\left(t\right)}{\mathrm{1}+\alpha I\left(t\right)}-\frac{\gamma I\left(t\right)}{\mathrm{1}+\theta I\left(t\right)}-\mu I\left(t\right),\\ \frac{\mathrm{d}R\left(t\right)}{\mathrm{d}t}=\frac{\gamma I\left(t\right)}{\mathrm{1}+\theta I\left(t\right)}-\left(\mu +\kappa \right)R\left(t\right)\mathrm{。}\end{array}\right. （4）

令 N \left(t\right) = S \left(t\right) + I \left(t\right) + R \left(t\right)，由模型（4）可得 \mathrm{d} N \left(t\right) / \mathrm{d} t = {\Lambda }_h - \mu N \left(t\right)。容易计算， \underset{t\to \mathrm{\infty }}{\mathrm{l}\mathrm{i}\mathrm{m}} N \left(t\right) = {\Lambda }_h / \mu。 因此，考虑到模型（4）的生物学背景，仅在如下区域讨论其动力学行为

\mathrm{\Omega } = \left\{\left(S\left(t\right),I\left(t\right),R\left(t\right)\right)\in {\mathbb{R}}_{+}^{\mathrm{3}}:\mathrm{0}\le S\left(t\right)+I\left(t\right)+R\left(t\right)\le \frac{{\Lambda }_h}{\mu }\right\} ,

其中 {\mathbb{R}}_{+}^{\mathrm{3}} ≔ \left\{\left(x_{\mathrm{1}},x_{\mathrm{2}},x_{\mathrm{3}}\right):x_i\ge \mathrm{0},i=\mathrm{1,2},\mathrm{3}\right\}。此外，定义模型（4）的基本再生数为

{\mathrm{ℛ}}_{\mathrm{0}}^h = \frac{S_{\mathrm{0}}\beta \left(\tilde{V}\right)}{\mu +\gamma } \mathrm{。}

3.1　平衡点的存在性

模型（4）始终存在无病平衡点 W_{\mathrm{0}} = \left(S_{\mathrm{0}},\mathrm{0},\mathrm{0}\right)。若模型（4）存在地方病平衡点 \tilde{W} = \left(\tilde{S},\tilde{I},\tilde{R}\right)，则其满足下面的方程组

\left\{\begin{array}{l}{\Lambda }_h-\frac{\beta \left(\tilde{V}\right)\tilde{S}\tilde{I}}{\mathrm{1}+\alpha \tilde{I}}-\mu \tilde{S}+\kappa \tilde{R}=\mathrm{0},\\ \frac{\beta \left(\tilde{V}\right)\tilde{S}\tilde{I}}{\mathrm{1}+\alpha \tilde{I}}-\frac{\gamma \tilde{I}}{\mathrm{1}+\theta \tilde{I}}-\mu \tilde{I}=\mathrm{0},\\ \frac{\gamma \tilde{I}}{\mathrm{1}+\theta \tilde{I}}-\left(\mu +\kappa \right)\tilde{R}=\mathrm{0}\mathrm{。}\end{array}\right. （5）

由式（5）的第2个和第3个方程可以得到

\tilde{S} = \frac{\mathrm{1}}{\beta \left(\tilde{V}\right)} \left(\frac{\gamma }{\mathrm{1}+\theta \tilde{I}}+\mu \right) \left(\mathrm{1}+\alpha \tilde{I}\right) , \tilde{R} = \frac{\gamma \tilde{I}}{\left(\mathrm{1}+\theta \tilde{I}\right)\left(\mu +\kappa \right)} \mathrm{。} （6）

将式（6）代入式（5）中的第1个方程，得到 \tilde{I}满足

a {\tilde{I}}^{\mathrm{2}} + b \tilde{I} + c = \mathrm{0} , （7）

其中

\begin{array}{c}a=\mu \theta \left(\frac{\mu \alpha }{\beta \left(\tilde{V}\right)}+\mathrm{1}\right),b=\frac{\mu }{\beta \left(\tilde{V}\right)}\left(\left(\mu +\gamma \right)\left(\frac{\beta \left(\tilde{V}\right)}{\mu }+\alpha \right)+\mu \theta \right)-\mu S_{\mathrm{0}}\theta -\frac{\kappa \gamma }{\mu +\kappa },\\ c=\frac{\mathrm{1}}{\beta \left(\tilde{V}\right)}\mu \left(\mu +\gamma \right)-\mu S_{\mathrm{0}}=\frac{\mathrm{1}}{\beta \left(\tilde{V}\right)}\mu \left(\mu +\gamma \right)\left(\mathrm{1}-{\mathrm{ℛ}}_{\mathrm{0}}^h\right)\mathrm{。}\end{array}

因此， \mathrm{\Delta } = b^{\mathrm{2}} - \mathrm{4} a c = b^{\mathrm{2}} - \mathrm{4} a \mu \left(\mu +\gamma \right) \left(\mathrm{1}-{\mathrm{ℛ}}_{\mathrm{0}}^h\right) / \beta \left(\tilde{V}\right)，且当 \mathrm{\Delta } \ge \mathrm{0}时，式（7）存在实根。现令

{\widetilde{\mathrm{ℛ}}}_{\mathrm{0}}^h = \mathrm{1} - \frac{b^{\mathrm{2}}\beta \left(\tilde{V}\right)}{\mathrm{4}a\mu \left(\mu +\gamma \right)} ,

显然可以得到

因此，当 {\mathrm{ℛ}}_{\mathrm{0}}^w > \mathrm{1}时，关于地方病平衡点的存在性，有下面的结论。

定理4 假设 {\mathrm{ℛ}}_{\mathrm{0}}^w > \mathrm{1}，则

1） 模型（4）不存在地方病平衡点当且仅当下列条件之一成立。

（i） {\widetilde{\mathrm{ℛ}}}_{\mathrm{0}}^h > {\mathrm{ℛ}}_{\mathrm{0}}^h；（ii） {\widetilde{\mathrm{ℛ}}}_{\mathrm{0}}^h = {\mathrm{ℛ}}_{\mathrm{0}}^h并且 b \ge \mathrm{0}；（iii） 并且 b > \mathrm{0}。

2） 模型（4）存在唯一的地方病平衡点 \overline{W}当且仅当下列条件之一成立。

（iv） 并且 ；（v） {\mathrm{ℛ}}_{\mathrm{0}}^h > \mathrm{1}；（vi） 并且 。

3） 模型（4）存在两个地方病平衡点 W_{\mathrm{*}} = \left(S_{\mathrm{*}},I_{\mathrm{*}},R_{\mathrm{*}}\right)和 W^{\mathrm{*}} = \left(S^{\mathrm{*}},I^{\mathrm{*}},R^{\mathrm{*}}\right)当且仅当下面的条件成立。

（vii） 并且 。

这里， I_{\mathrm{*}} = \left(-b-\sqrt[]{\mathrm{\Delta }}\right) / \mathrm{2} a和 I^{\mathrm{*}} = \left(-b+\sqrt[]{\mathrm{\Delta }}\right) / \mathrm{2} a。

3.2　平衡点的稳定性

首先讨论 W_{\mathrm{0}}的局部稳定性。

定理5 当 时， W_{\mathrm{0}}是局部渐近稳定的；当 {\mathrm{ℛ}}_{\mathrm{0}}^h > \mathrm{1}时， W_{\mathrm{0}}不稳定。

证明 模型（4）在 W_{\mathrm{0}}处的Jacobi矩阵为

\mathit{J} \left(W_{\mathrm{0}}\right) = \left(\begin{array}{ccc}-\mu & -\beta \left(\tilde{V}\right)S_{\mathrm{0}}& \kappa \\ \mathrm{0}& \beta \left(\tilde{V}\right)S_{\mathrm{0}}-\gamma -\mu & \mathrm{0}\\ \mathrm{0}& \gamma & -\mu -\kappa \end{array}\right) , （8）

因此，特征方程为

\left(\lambda +\mu \right) \left(\lambda +\mu +\kappa \right) \left(\lambda -\beta \left(\tilde{V}\right)S_{\mathrm{0}}+\gamma +\mu \right) = \mathrm{0} ,

此特征方程有两个负特征值 - \mu和 - \mu - \kappa，还有一个特征值

\lambda = \beta \left(\tilde{V}\right) S_{\mathrm{0}} - \gamma - \mu = \left(\gamma +\mu \right) \left(\mathrm{ℛ}{}_{\mathrm{0}}^h-\mathrm{1}\right) \mathrm{。}

当 时，所有的特征值都具有负实部，因此 W_{\mathrm{0}}是局部渐近稳定的；而当 {\mathrm{ℛ}}_{\mathrm{0}}^h > \mathrm{1}时，存在一个特征值大于零，即 W_{\mathrm{0}}是不稳定的。证毕。

下面讨论 W_{\mathrm{*}}和 W^{\mathrm{*}}的局部稳定性。模型（4）在任意地方病平衡点处的Jacobi矩阵为

\mathit{J} = \left(\begin{array}{ccc}-\frac{\beta \left(\tilde{V}\right)I}{\mathrm{1}+\alpha I}-\mu & -\frac{\mathrm{1}}{\mathrm{1}+\alpha I}\left(\frac{\gamma }{\mathrm{1}+\theta I}+\mu \right)& \kappa \\ \frac{\beta \left(\tilde{V}\right)I}{\mathrm{1}+\alpha I}& \frac{\mathrm{1}}{\mathrm{1}+\alpha I}\left(\frac{\gamma }{\mathrm{1}+\theta I}+\mu \right)-\frac{\gamma }{{(\mathrm{1}+\theta I)}^{\mathrm{2}}}-\mu & \mathrm{0}\\ \mathrm{0}& \frac{\gamma }{{(\mathrm{1}+\theta I)}^{\mathrm{2}}}& -\mu -\kappa \end{array}\right) ,

因此，特征方程为

\begin{array}{c}\left(\lambda +\mu \right)\left({\lambda }^{\mathrm{2}}-\mathrm{t}\mathrm{r}{J}_{\mathrm{1}}\left(I\right)\lambda +\mathrm{d}\mathrm{e}\mathrm{t}{J}_{\mathrm{1}}\left(I\right)\right)=\mathrm{0}\end{array} , （9）

其中

\begin{array}{c}\mathrm{d}\mathrm{e}\mathrm{t}{J}_{\mathrm{1}}\left(I\right)=\frac{I}{\left(\mathrm{1}+\alpha I\right){(\mathrm{1}+\theta I)}^{\mathrm{2}}}F\left(I\right),\\ F\left(I\right)=\left(\mu +\kappa \right)\left(\alpha \mu +\beta \left(\tilde{V}\right)\right){\theta }^{\mathrm{2}}{I}^{\mathrm{2}}+\mathrm{2}\theta \left(\mu +\kappa \right)\left(\alpha \mu +\beta \left(\tilde{V}\right)\right)I+\\ \beta \left(\tilde{V}\right)\left(\mu +\kappa +\gamma \right)+\alpha \left(\mu +\gamma \right)\left(\mu +\kappa \right)-\gamma \theta \left(\mu +\kappa \right);\\ \mathrm{t}\mathrm{r}{J}_{\mathrm{1}}\left(I\right)=\frac{\mathrm{1}}{\left(\mathrm{1}+\alpha I\right){(\mathrm{1}+\theta I)}^{\mathrm{2}}}G\left(I\right),G\left(I\right)=A{I}^{\mathrm{3}}+B{I}^{\mathrm{2}}+CI+D,\\ A=-{\theta }^{\mathrm{2}}\left(\beta \left(\tilde{V}\right)+\alpha \left(\mathrm{2}\mu +\kappa \right)\right),B=-\left({\theta }^{\mathrm{2}}\left(\mu +\kappa \right)+\mathrm{2}\theta \left(\beta \left(\tilde{V}\right)+\alpha \left(\mathrm{2}\mu +\kappa \right)\right)\right),\\ C=\gamma \theta -\mathrm{2}\theta \left(\mu +\kappa \right)-\alpha \left(\mathrm{2}\mu +\kappa +\gamma \right)-\beta \left(\tilde{V}\right),D=-\left(\mu +\kappa \right)\mathrm{。}\end{array}

定理6 如果 且 ，模型（4）的地方病平衡点 W_{\mathrm{*}}是鞍点。

证明 特征方程（9）显然有一个特征值 - \mu，其余特征值由 {\lambda }^{\mathrm{2}} - \mathrm{t} \mathrm{r} J_{\mathrm{1}} \left(I\right) \lambda + \mathrm{d} \mathrm{e} \mathrm{t} J_{\mathrm{1}} \left(I\right) = \mathrm{0}确定。为此，将 I_{\mathrm{*}} = \left(-b-\sqrt[]{\mathrm{\Delta }}\right) / \left(\mathrm{2}a\right)代入 F \left(I\right)得

F \left(I_{\mathrm{*}}\right) = \left(\mu +\kappa \right) \left(\alpha \mu +\beta \left(\tilde{V}\right)\right) {\theta }^{\mathrm{2}} {\left(\frac{-b-\sqrt[]{\mathrm{\Delta }}}{\mathrm{2}a}\right)}^{\mathrm{2}} - \gamma \theta \left(\mu +\kappa \right) + \mathrm{2} \theta \left(\mu +\kappa \right) \left(\alpha \mu +\beta \left(\tilde{V}\right)\right) \frac{-b-\sqrt[]{\mathrm{\Delta }}}{\mathrm{2}a} +

即 ，则当 W_{\mathrm{*}}存在时，它是鞍点。证毕。

定理7 当 并且 时，

1） 若 ，模型（4）的地方病平衡点 W^{\mathrm{*}}是局部渐近稳定的；

2） 若 C > \mathrm{0}，模型（4）的地方病平衡点 W^{\mathrm{*}}是局部渐近稳定的当且仅当下列条件之一成立，

（i） ；（ii） G \left(I_{\mathrm{3}}\right) > \mathrm{0}且 I^{\mathrm{*}} \in \left(\mathrm{0},I_{\mathrm{2}}\right) \bigcup \left(I_{\mathrm{4}},+\mathrm{\infty }\right)，

这里， I_{\mathrm{1}}和 I_{\mathrm{3}}分别为 G \left(I\right)的极小值点和极大值点， I_{\mathrm{2}}和 I_{\mathrm{4}}为 G \left(I\right) = \mathrm{0}的两个正实根。

证明 类似于定理6的证明，将 I^{\mathrm{*}} = \left(-b+\sqrt[]{\mathrm{\Delta }}\right) / \left(\mathrm{2}a\right)代入 F \left(I\right)得

F \left(I^{\mathrm{*}}\right) = \frac{\theta \beta \left(\tilde{V}\right)\left(\mu +\kappa \right)}{\mu } \sqrt[]{\mathrm{\Delta }} \left(\frac{\mathrm{1}}{\theta }+\frac{-b+\sqrt[]{\mathrm{\Delta }}}{\mathrm{2}a}\right) > \mathrm{0} ,

即 \mathrm{d} \mathrm{e} \mathrm{t} J_{\mathrm{1}} \left(I^{\mathrm{*}}\right) > \mathrm{0}。接下来考虑 \mathrm{t} \mathrm{r} J_{\mathrm{1}} \left(I^{\mathrm{*}}\right)的符号：

1） 当 时，由 G \left(I\right)的表达式可知， ，显然有 ，则有 ；

2） 当 C > \mathrm{0}时，易知 I_{\mathrm{1}} = \left(-\mathrm{2}B+\sqrt[]{\mathrm{4}{B}^{\mathrm{2}}-\mathrm{12}AC}\right) / \left(\mathrm{6}A\right)和 I_{\mathrm{3}} = \left(-\mathrm{2}B-\sqrt[]{\mathrm{4}{B}^{\mathrm{2}}-\mathrm{12}AC}\right) / \left(\mathrm{6}A\right)分别是 G \left(I\right)的极小值点和极大值点。（i）当 时， ，则有 ；（ii）当 G \left(I_{\mathrm{3}}\right) > \mathrm{0}时，若 I^{\mathrm{*}} \in \left(\mathrm{0},I_{\mathrm{2}}\right) \bigcup \left(I_{\mathrm{4}},+\mathrm{\infty }\right)，则有 。

综上所述可得 W^{\mathrm{*}}是局部渐近稳定的。证毕。

最后讨论当 {\mathrm{ℛ}}_{\mathrm{0}}^h > \mathrm{1}时，地方病平衡点 \overline{W} = \left(\overline{S},\overline{I},\overline{R}\right)的全局稳定性。

定理8 若 {\mathrm{ℛ}}_{\mathrm{0}}^h > \mathrm{1}，则

1） 当 \alpha = \theta = \mathrm{0}时，地方病平衡点是全局渐近稳定的；

2） 当 \alpha = \mathrm{0} , \theta \ne \mathrm{0}时，地方病平衡点是全局渐近稳定的当且仅当下面的条件成立，

\mathrm{1} - \frac{\mu }{\left(\mu +\theta {\Lambda }_h\right)\left(\mathrm{1}+\theta \overline{I}\right)} \le \frac{\mathrm{2}(\mu +\kappa )}{\gamma } \mathrm{。} （10）

证明 对地方病平衡点 \overline{W}，构造Lyapunov函数如下

L_{\mathrm{2}} = \frac{\mathrm{1}}{\mathrm{2}} {(S-\overline{S}+I-\overline{I}+R-\overline{R})}^{\mathrm{2}} + \frac{\mathrm{2}\mu }{\beta } \left(I-\overline{I}-\overline{I}\mathrm{l}\mathrm{n}\frac{I}{\overline{I}}\right) + \frac{\mu }{\gamma } {(R-\overline{R})}^{\mathrm{2}} \mathrm{。}

沿着模型（4）计算 L_{\mathrm{2}}的全导数，可以得到

\begin{array}{c}\frac{\mathrm{d}{L}_{\mathrm{2}}}{\mathrm{d}t}=\left(S-\overline{S}+I-\overline{I}+R-\overline{R}\right)\left(\frac{\mathrm{d}S}{\mathrm{d}t}+\frac{\mathrm{d}I}{\mathrm{d}t}+\frac{\mathrm{d}R}{\mathrm{d}t}\right)+\frac{\mathrm{2}\mu }{\beta }\left(\mathrm{1}-\frac{\overline{I}}{I}\right)\frac{\mathrm{d}I}{\mathrm{d}t}+\frac{\mathrm{2}\mu }{\gamma }\left(R-\overline{R}\right)\frac{\mathrm{d}R}{\mathrm{d}t}=\\ -\mu {(S-\overline{S})}^{\mathrm{2}}-\left(\mu +\frac{\mathrm{2}\mu \gamma \theta }{\beta (\mathrm{1}+\theta I)(\mathrm{1}+\theta \overline{I})}\right){(I-\overline{I})}^{\mathrm{2}}-\left(\mu +\frac{\mathrm{2}\mu (\mu +\kappa )}{\gamma }\right){(R-\overline{R})}^{\mathrm{2}}-\mathrm{2}\mu \left(S-\overline{S}\right)\left(I-\overline{I}\right)-\\ \mathrm{2}\mu \left(S-\overline{S}\right)\left(R-\overline{R}\right)-\mathrm{2}\mu \left(I-\overline{I}\right)\left(R-\overline{R}\right)+\frac{\mathrm{2}\mu \left(I-\overline{I}\right)\left(R-\overline{R}\right)}{(\mathrm{1}+\theta I)(\mathrm{1}+\theta \overline{I})}+\frac{\mathrm{2}\mu }{\beta }\left(I-\overline{I}\right)\left(\frac{\beta S}{\mathrm{1}+\alpha I}-\frac{\beta \overline{S}}{\mathrm{1}+\alpha \overline{I}}\right)\mathrm{。}\end{array}

当 \alpha = \theta = \mathrm{0}时，

\frac{\mathrm{d}{L}_{\mathrm{2}}}{\mathrm{d}t} = - \mu {(S-\overline{S})}^{\mathrm{2}} - \mu {(I-\overline{I})}^{\mathrm{2}} - \left(\mu +\frac{\mathrm{2}\mu (\mu +\kappa )}{\gamma }\right) {(R-\overline{R})}^{\mathrm{2}} - \mathrm{2} \mu \left(S-\overline{S}\right) \left(R-\overline{R}\right) \le \mathrm{0} \mathrm{。}

当 \alpha = \mathrm{0} , \theta \ne \mathrm{0}时，

\begin{array}{c}\frac{\mathrm{d}{L}_{\mathrm{2}}}{\mathrm{d}t}=-\mu {(S-\overline{S})}^{\mathrm{2}}-\left(\mu +\frac{\mathrm{2}\mu \gamma \theta }{\beta \left(\mathrm{1}+\theta I\right)\left(\mathrm{1}+\theta \overline{I}\right)}\right){(I-\overline{I})}^{\mathrm{2}}-\left(\mu +\frac{\mathrm{2}\mu (\mu +\kappa )}{\gamma }\right){(R-\overline{R})}^{\mathrm{2}}-\\ \mathrm{2}\mu (S-\overline{S})(R-\overline{R})-\mathrm{2}\mu \left(\mathrm{1}-\frac{\mathrm{1}}{\left(\mathrm{1}+\theta I\right)\left(\mathrm{1}+\theta \overline{I}\right)}\right)(I-\overline{I})(R-\overline{R})\le -\left(\mu +\frac{\mathrm{2}\mu \gamma \theta }{\beta \left(\mathrm{1}+\theta I\right)\left(\mathrm{1}+\theta \overline{I}\right)}\right)\\ {(I-\overline{I})}^{\mathrm{2}}-\frac{\mathrm{2}\mu (\mu +\kappa )}{\gamma }{(R-\overline{R})}^{\mathrm{2}}+\mathrm{2}\mu \left(\mathrm{1}-\frac{\mu }{(\mu +\theta {\mathrm{\Lambda }}_h)\left(\mathrm{1}+\theta \overline{I}\right)}\right)\left|I-\overline{I}\right|\left|R-\overline{R}\right|,\end{array}

当 {\mathrm{ℛ}}_{\mathrm{0}}^h > \mathrm{1}时，若条件（10）成立，则对任意的 \left(S\left(t\right),I\left(t\right),R\left(t\right)\right)有 。当且仅当 S = \overline{S} , I = \overline{I} , R = \overline{R}时，有 \mathrm{d} L_{\mathrm{2}} / \mathrm{d} t = \mathrm{0}。因而， \left\{\overline{W}\right\}是 \left\{\left(S\left(t\right),I\left(t\right),R\left(t\right)\right)\in \mathrm{\Omega }|\mathrm{d}{L}_{\mathrm{2}}\left(t\right)/\mathrm{d}t=\mathrm{0}\right\}中的最大不变集，由Lyapunov定理^［16］和LaSalle不变原理^［17］可知 \overline{W}是全局渐近稳定的。证毕。

注 在定理8中，只建立了当 {\mathrm{ℛ}}_{\mathrm{0}}^h > \mathrm{1}时， \alpha = \theta = \mathrm{0}和 \alpha = \mathrm{0} , \theta \ne \mathrm{0}情况下的地方病平衡点的全局稳定性。是否可以通过构造合适的Lyapunov函数或给出不含变量的条件去证明当 \alpha \ne \mathrm{0} , \theta \ne \mathrm{0}时地方病平衡点的全局稳定性是一个有待解决的问题，我们将在数值模拟部分验证当 \alpha \ne \mathrm{0} , \theta \ne \mathrm{0}时地方病平衡点的全局稳定性。

3.3　前向/后向分支的存在性

关于模型（4）前向/后向分支的存在性，有下面的结论。

定理9 若 \theta > \widehat{\theta }，则模型（4）在 {\mathrm{ℛ}}_{\mathrm{0}}^h = \mathrm{1}处出现后向分支；若 ，则模型（4）在 {\mathrm{ℛ}}_{\mathrm{0}}^h = \mathrm{1}处出现前向分支，其中

\widehat{\theta } = \frac{\left(\mu +\gamma \right)\left(\mathrm{2}\alpha \left(\mu +\kappa \right)S_{\mathrm{0}}+\mu +\kappa +\gamma \right)}{\mathrm{2}\gamma \left(\mu +\kappa \right)S_{\mathrm{0}}} \mathrm{。}

证明 令 S = x_{\mathrm{1}} , I = x_{\mathrm{2}} , R = x_{\mathrm{3}}，则模型（4）可改写为

\left\{\begin{array}{l}\frac{\mathrm{d}{x}_{\mathrm{1}}\left(t\right)}{\mathrm{d}t}={\Lambda }_h-\frac{\beta \left(\tilde{V}\right)x_{\mathrm{1}}\left(t\right)x_{\mathrm{2}}\left(t\right)}{\mathrm{1}+\alpha x_{\mathrm{2}}\left(t\right)}-\mu x_{\mathrm{1}}\left(t\right)+\kappa x_{\mathrm{3}}\left(t\right)≔f_{\mathrm{1}},\\ \frac{\mathrm{d}{x}_{\mathrm{2}}\left(t\right)}{\mathrm{d}t}=\frac{\beta \left(\tilde{V}\right)x_{\mathrm{1}}\left(t\right)x_{\mathrm{2}}\left(t\right)}{\mathrm{1}+\alpha x_{\mathrm{2}}\left(t\right)}-\frac{\gamma x_{\mathrm{2}}\left(t\right)}{\mathrm{1}+\theta x_{\mathrm{2}}\left(t\right)}-\mu x_{\mathrm{2}}\left(t\right)≔f_{\mathrm{2}},\\ \frac{\mathrm{d}{x}_{\mathrm{3}}\left(t\right)}{\mathrm{d}t}=\frac{\gamma x_{\mathrm{2}}\left(t\right)}{\mathrm{1}+\theta x_{\mathrm{2}}\left(t\right)}-\mu x_{\mathrm{3}}\left(t\right)-\kappa x_{\mathrm{3}}\left(t\right)≔f_{\mathrm{3}}\mathrm{。}\end{array}\right.

当 {\mathrm{ℛ}}_{\mathrm{0}}^h = \mathrm{1}时，由（8）可知 \mathit{J} \left(W_{\mathrm{0}}\right)对应的特征方程为

\lambda \left(\lambda +\mu \right) \left(\lambda +\mu +\kappa \right) = \mathrm{0} \mathrm{。}

因此， \mathit{J} \left(W_{\mathrm{0}}\right)的特征值分别为 {\lambda }_{\mathrm{1}} = \mathrm{0} , {\lambda }_{\mathrm{2}} = - \mu , {\lambda }_{\mathrm{3}} = - \mu - \kappa。显然 {\lambda }_{\mathrm{1}} = \mathrm{0}是一个零特征值， \mathit{J} \left(W_{\mathrm{0}}\right)的其他特征值都有负实部。因此，可以利用文献［18］中的中心流形理论来讨论分支的存在性。

选择病毒传播率 \beta \left(\tilde{V}\right)作为分支参数，则 {\mathrm{ℛ}}_{\mathrm{0}}^h = \mathrm{1}等价于

\beta \left(\tilde{V}\right) = \frac{\mu +\gamma }{S_{\mathrm{0}}} ≔ {\beta }^{\mathrm{*}} \mathrm{。}

因此，当 \beta \left(\tilde{V}\right) = {\beta }^{\mathrm{*}}时，由文献［19］中的定理4.1可知无病平衡点 W_{\mathrm{0}}是一个非双曲平衡点。

令 \mathit{w} = ( {\omega }_{\mathrm{1}} , {\omega }_{\mathrm{2}} , {\omega }_{\mathrm{3}} {)}^T表示与零特征值 {\lambda }_{\mathrm{1}} = \mathrm{0}对应的右特征向量，其满足 \mathit{J} \left(W_{\mathrm{0}},{\beta }^{\mathrm{*}}\right) \mathit{w} = \mathrm{0}，即

\left\{\begin{array}{l}-\mu {\omega }_{\mathrm{1}}-{\beta }^{\mathrm{*}}{S}_{\mathrm{0}}{\omega }_{\mathrm{2}}+\kappa {\omega }_{\mathrm{3}}=\mathrm{0},\\ \left({\beta }^{\mathrm{*}}{S}_{\mathrm{0}}-\mu -\gamma \right){\omega }_{\mathrm{2}}=\mathrm{0},\\ \gamma {\omega }_{\mathrm{2}}-\left(\mu +\kappa \right){\omega }_{\mathrm{3}}=\mathrm{0}\mathrm{。}\end{array}\right.

由此得

\begin{array}{c}{\omega }_{\mathrm{1}}=-\left(\mu +\kappa +\gamma \right),{\omega }_{\mathrm{2}}=\mu +\kappa ,{\omega }_{\mathrm{3}}=\gamma \end{array} \mathrm{。} （11）

此外，满足 \mathit{v} \cdot \mathit{w} = \mathrm{1}的与零特征值对应的左特征向量 \mathit{v} = \left(v_{\mathrm{1}},v_{\mathrm{2}},v_{\mathrm{3}}\right)由下面的方程组给出

\left\{\begin{array}{l}-\mu v_{\mathrm{1}}=\mathrm{0},\\ -{\beta }^{\mathrm{*}}{S}_{\mathrm{0}}{v}_{\mathrm{1}}+\left({\beta }^{\mathrm{*}}{S}_{\mathrm{0}}-\mu -\gamma \right)v_{\mathrm{2}}+\gamma v_{\mathrm{3}}=\mathrm{0},\\ \kappa v_{\mathrm{1}}-\left(\mu +\kappa \right)v_{\mathrm{3}}=\mathrm{0}\mathrm{。}\end{array}\right.

由此得

v_{\mathrm{1}} = \mathrm{0} , v_{\mathrm{2}} = \frac{\mathrm{1}}{\mu +\kappa } , v_{\mathrm{3}} = \mathrm{0} \mathrm{。} （12）

计算在无病平衡点 W_{\mathrm{0}}处的二阶偏导数得

\begin{array}{c}\frac{{\partial }^{\mathrm{2}}{f}_{\mathrm{1}}}{\partial x_{\mathrm{1}}\partial x_{\mathrm{2}}}=\frac{{\partial }^{\mathrm{2}}{f}_{\mathrm{1}}}{\partial x_{\mathrm{2}}\partial x_{\mathrm{1}}}=-\beta \left(\tilde{V}\right),\frac{{\partial }^{\mathrm{2}}{f}_{\mathrm{1}}}{\partial x_{\mathrm{2}}\partial \beta \left(\tilde{V}\right)}=-S_{\mathrm{0}},\frac{{\partial }^{\mathrm{2}}{f}_{\mathrm{1}}}{\partial x_{\mathrm{2}}^{\mathrm{2}}}=\mathrm{2}\alpha \beta \left(\tilde{V}\right)S_{\mathrm{0}},\frac{{\partial }^{\mathrm{2}}{f}_{\mathrm{3}}}{\partial x_{\mathrm{2}}^{\mathrm{2}}}=-\mathrm{2}\theta \gamma ,\\ \frac{{\partial }^{\mathrm{2}}{f}_{\mathrm{2}}}{\partial x_{\mathrm{1}}\partial x_{\mathrm{2}}}=\frac{{\partial }^{\mathrm{2}}{f}_{\mathrm{2}}}{\partial x_{\mathrm{2}}\partial x_{\mathrm{1}}}=\beta \left(\tilde{V}\right),\frac{{\partial }^{\mathrm{2}}{f}_{\mathrm{2}}}{\partial x_{\mathrm{2}}\partial \beta \left(\tilde{V}\right)}=S_{\mathrm{0}},\frac{{\partial }^{\mathrm{2}}{f}_{\mathrm{2}}}{\partial x_{\mathrm{2}}^{\mathrm{2}}}=-\mathrm{2}\alpha \beta \left(\tilde{V}\right)S_{\mathrm{0}}+\mathrm{2}\theta \gamma ,\end{array} （13）

且其他二阶偏导数都等于零。因此由式（11），式（12）和式（13），可以计算文献［19］定理4.1中定义的系数 a和 b如下，

a = {\displaystyle \sum _{k,i,j=\mathrm{1}}^n} v_{\kappa } {\omega }_i {\omega }_j \frac{{\partial }^{\mathrm{2}}{f}_{\kappa }}{\partial x_i\partial x_j} \left(W_{\mathrm{0}}\right) = - \left(\mu +\kappa +\gamma \right) {\beta }^{\mathrm{*}} - \mathrm{2} \alpha {\beta }^{\mathrm{*}} S_{\mathrm{0}} \left(\mu +\kappa \right) + \mathrm{2} \theta \gamma \left(\mu +\kappa \right) \mathrm{。}

b = {\displaystyle \sum _{k,i=\mathrm{1}}^n} v_{\kappa } {\omega }_i \frac{{\partial }^{\mathrm{2}}{f}_{\kappa }}{\partial x_i\partial \beta \left(\tilde{V}\right)} \left(W_{\mathrm{0}}\right) = v_{\mathrm{2}} {\omega }_{\mathrm{2}} \frac{{\partial }^{\mathrm{2}}{f}_{\mathrm{2}}}{\partial x_{\mathrm{2}}\partial \beta \left(\tilde{V}\right)} \left(W_{\mathrm{0}}\right) = S_{\mathrm{0}} > \mathrm{0} \mathrm{。}

由于系数 b总为正，因此，若 \theta > \widehat{\theta }，即系数 a为正，则模型（4）存在后向分支；若 ，即系数 a为负，则模型（4）存在前向分支。证毕。

4 数值模拟

本节使用MATLAB软件进行数值模拟。首先，解释当宿主体内模型（2）的基本再生数 {\mathrm{ℛ}}_{\mathrm{0}}^w > \mathrm{1}时，宿主体内病毒浓度对宿主间疾病传播的影响。以函数 \beta \left(V\right) = a V为例，选取 a = \mathrm{0.77} \times {\mathrm{10}}^{-\mathrm{11}}。关于宿主间模型（4），选取 {\Lambda }_h = \mathrm{750}， \mu = \mathrm{0.000} \mathrm{182}， \kappa = \mathrm{0.002}， \gamma = \mathrm{0.11}， \theta = \mathrm{1.85}， \alpha = \mathrm{0.20}，参数单位为d，有下列三种情形。

（1）当 {\Lambda }_c = \mathrm{1} \mathrm{500}， p = \mathrm{23}， {\kappa }_{\mathrm{1}} = \mathrm{18} \times {\mathrm{10}}^{-\mathrm{7}}， {\kappa }_{\mathrm{2}} = \mathrm{15} \times {\mathrm{10}}^{-\mathrm{7}}， {\mu }_c = \mathrm{0.011} \mathrm{05}， {\delta }_c = \mathrm{0.68}， {\mu }_v = \mathrm{10.258}（参数单位为d）时，计算得到 {\mathrm{ℛ}}_{\mathrm{0}}^w = \mathrm{1.087} \mathrm{4} > \mathrm{1}，从而模拟得到体内病毒浓度随着时间 s的变化趋于稳定值 \tilde{V} = \mathrm{391.342} \mathrm{3}，如图1（a）所示。此时，宿主间模型（4）的基本再生数 ，在这种情况下，无病平衡点 W_{\mathrm{0}}是全局渐近稳定的，意味着疾病灭绝，如图1（b）所示。

（2）当 {\Lambda }_c = \mathrm{1} \mathrm{290}， p = \mathrm{30}， {\kappa }_{\mathrm{1}} = \mathrm{55} \times {\mathrm{10}}^{-\mathrm{7}}， {\kappa }_{\mathrm{2}} = \mathrm{50} \times {\mathrm{10}}^{-\mathrm{7}}， {\mu }_c = \mathrm{0.011} \mathrm{2}， {\delta }_c = \mathrm{0.58}， {\mu }_v = \mathrm{12}（参数单位为d）时，计算得到 {\mathrm{ℛ}}_{\mathrm{0}}^w = \mathrm{3.652} \mathrm{9} > \mathrm{1}，从而模拟得到体内病毒浓度随着时间 s的变化趋于稳定值 \tilde{V} = \mathrm{3} \mathrm{961.7}，如图2（a）所示。此时，宿主间模型（4）的基本再生数 {\mathrm{ℛ}}_{\mathrm{0}}^h = \mathrm{1.140} \mathrm{9} > \mathrm{1}，在这种情况下，存在唯一的地方病平衡点 \overline{W}，且其是全局渐近稳定的，意味着疾病持续存在，如图2（b）所示。

（3）当 {\Lambda }_c = \mathrm{1} \mathrm{500}， p = \mathrm{23}， {\kappa }_{\mathrm{1}} = \mathrm{18} \times {\mathrm{10}}^{-\mathrm{7}}， {\kappa }_{\mathrm{2}} = \mathrm{15} \times {\mathrm{10}}^{-\mathrm{7}}， {\mu }_c = \mathrm{0.011} \mathrm{05}， {\delta }_c = \mathrm{0.68}， {\mu }_v = \mathrm{9.248}（参数单位为d）时，计算得到 {\mathrm{ℛ}}_{\mathrm{0}}^w = \mathrm{1.174} \mathrm{0} > \mathrm{1}，从而模拟得到体内病毒浓度随着时间 s的变化趋于稳定值 \tilde{V} = \mathrm{800}，如图3（a）所示。此时，宿主间模型（4）的基本再生数 ，并且存在双稳态，这意味着基本再生数不再是确定疾病流行的阈值条件，疾病的持续或灭绝更多地取决于感染者的初始数量，如图3（b）所示。

接下来，选取参数值 {\Lambda }_h = \mathrm{750}， \mu = \mathrm{0.000} \mathrm{182}， \kappa = \mathrm{0.002}， \gamma = \mathrm{0.11}， \alpha = \mathrm{0.20}， \tilde{V} = \mathrm{800}， \theta = \mathrm{8.338}，计算得到 {\widetilde{\mathrm{ℛ}}}_{\mathrm{0}}^h = \mathrm{0.231}， 。注意到，当 并且 时，模型（4）存在两个地方病平衡点，一个不稳定，另一个局部渐近稳定；当 {\mathrm{ℛ}}_{\mathrm{0}}^h > \mathrm{1}时，模型（4）存在唯一的地方病平衡点。这表示控制疾病传播的阈值是 {\mathrm{ℛ}}_{\mathrm{0}}^h = {\widetilde{\mathrm{ℛ}}}_{\mathrm{0}}^h = \mathrm{0.231}。换句话说，传染病只能在 之前得到控制，如图4（a）所示。另外，改变 \theta的取值，可以得到当 \theta > \widehat{\theta } = \mathrm{8.338}时，模型（4）存在两个地方病平衡点，这意味着， \theta越大，即医疗资源越有限，传染病越难得到控制，如图4（b）所示。

5 结论

基于传染病传播过程中的复杂性和医疗资源的有限性，本文提出了一类具有饱和治疗率的多尺度传染病模型。对于宿主体内病毒感染模型（2），得到了基本再生数 {\mathrm{ℛ}}_{\mathrm{0}}^w的表达式，建立了无感染平衡点和感染平衡点的全局渐近稳定性的判别准则。对于宿主之间疾病传播模型（4），得到了基本再生数 {\mathrm{ℛ}}_{\mathrm{0}}^h的表达式，建立了各类平衡点存在性和稳定性的判别准则。特别地，当 并且 时，模型（4）存在两个地方病平衡点 W_{\mathrm{*}}和 W^{\mathrm{*}}，其中 W_{\mathrm{*}}不稳定， W^{\mathrm{*}}局部渐近稳定；而当 {\mathrm{ℛ}}_{\mathrm{0}}^h > \mathrm{1}时，模型（4）存在唯一的地方病平衡点 \overline{W}，它在特定条件下是全局渐近稳定的。

注意到当 \theta > \widehat{\theta }时，模型（4）在 {\mathrm{ℛ}}_{\mathrm{0}}^h = \mathrm{1}处存在后向分支，由此可知饱和治疗可能导致多个地方病平衡点的存在，这使得对疾病的控制变得更加困难。此外，由 {\mathrm{ℛ}}_{\mathrm{0}}^h的表达式可以看出， \beta \left(\tilde{V}\right)可以决定 {\mathrm{ℛ}}_{\mathrm{0}}^h > \mathrm{1}或 。若初始感染者数量多，感染宿主体内病毒浓度越高，由 \beta \left(V\right)的性质可知 \beta \left(\tilde{V}\right)越大，直至使得 ，此时模型（4）存在两个地方病平衡点，这表明即使基本再生数 ，宿主之间疾病仍可能持续存在。因此，增加医疗资源并降低宿主体内病毒浓度对控制疾病传播甚至根除疾病具有不可或缺的作用。病毒在宿主体内会存在潜伏期，所以在模型中考虑潜伏期时滞等符合实际的因素，是进一步值得研究的问题。

参考文献

[1] RHOUBARI Z EL, HATTAF K, YOUSFI N. A Class of Ebola Virus Disease Models with Post-death Transmission and Environmental Contamination[M]//Studies in Systems, Decision and Control. Cham: Springer International Publishing, 2020, 302: 295-321. DOI: 10.1007/978-3-030-49896-2_11 .

[2] XU R. Global Stability of an HIV-1 Infection Model with Saturation Infection and Intracellular Delay[J]. J Math Anal Appl, 2011, 375(1): 75-81. DOI: 10.1016/j.jmaa.2010.08.055 .

[3] 李淑萍, 苗慧. 数据驱动下带环境病毒的新冠肺炎传播模型[J]. 山西大学学报(自然科学版), 2023, 46(1): 131-140. DOI: 10.13451/j.sxu.ns.2021178 .

[4] LI S P, MIAO H. A COVID-19 Spreading Model with Environmental Viruses Based on Data-driven[J]. J Shanxi Univ Nat Sci Ed, 2023, 46(1): 131-140. DOI: 10.13451/j.sxu.ns.2021178 .

[5] 李晓伟, 李桂花. 考虑环境病毒影响的COVID-19模型的动力学性态研究[J]. 山东大学学报(理学版), 2023, 58(1): 10-15. DOI: 10.6040/j.issn.1671-9352.0.2021.460 .

[6] LI X W, LI G H. Dynamic Behaviors Analysis of COVID-19 Model with Environmental Virus Effects[J]. J Shandong Univ Nat Sci, 2023, 58(1): 10-15. DOI: 10.6040/j.issn.1671-9352.0.2021.460 .

[7] AILI A, TENG Z D, ZHANG L. Dynamics in a Disease Transmission Model Coupled Virus Infection in Host with Incubation Delay and Environmental Effects[J].J Appl Math Comput, 2022, 68(6): 4331-4359. DOI: 10.1007/s12190-022-01709-y .

[8] ALMOCERA A E S, NGUYEN V K, HERNANDEZ-VARGAS E A. Multiscale Model Within-host and Between-host for Viral Infectious Diseases[J]. J Math Biol, 2018, 77(4): 1035-1057. DOI: 10.1007/s00285-018-1241-y .

[9] FENG Z L, CEN X L, ZHAO Y L, et al. Coupled Within-host and Between-host Dynamics and Evolution of Virulence[J]. Math Biosci, 2015, 270: 204-212. DOI: 10.1016/j.mbs.2015.02.012 .

[10] 唐思甜, 滕志东. 一类具有体液免疫的宿主内部和宿主之间的疾病传播耦合模型[J]. 数学的实践与认识, 2019, 49(7): 276-287. DOI: 10.1016/j.mbs.2012.09.004 .

[11] TANG S T, TENG Z D. Coupled Within-host and Between-host Dynamics of Infectious Disease Model with Humoral Immunity[J]. Math Pract Theory, 2019, 49(7): 276-287. DOI: 10.1016/j.mbs.2012.09.004 .

[12] GILCHRIST M A, COOMBS D. Evolution of Virulence: Interdependence, Constraints, and Selection Using Nested Models[J]. Theor Popul Biol, 2006, 69(2): 145-153. DOI: 10.1016/j.tpb.2005.07.002 .

[13] FENG Z L, VELASCO-HERNANDEZ J, TAPIA-SANTOS B, et al. A Model for Coupling Within-host and Between-host Dynamics in an Infectious Disease[J]. Nonlinear Dyn, 2012, 68(3): 401-411. DOI: 10.1007/s11071-011-0291-0 .

[14] FENG Z L, VELASCO-HERNANDEZ J, TAPIA-SANTOS B. A Mathematical Model for Coupling Within-host and Between-host Dynamics in an Environmentally-driven Infectious Disease[J]. Math Biosci, 2013, 241(1): 49-55. DOI: 10.1016/j.mbs.2012.09.004 .

[15] YADAV A, SRIVASTAVA P K. The Impact of Information and Saturated Treatment with Time Delay in an Infectious Disease Model[J].J Appl Math Comput, 2021, 66(1/2): 277-305. DOI: 10.1007/s12190-020-01436-2 .

[16] HU L, WANG S F, ZHENG T T, et al. The Effects of Migration and Limited Medical Resources of the Transmission of SARS-COV-2 Model with Two Patches[J].Bull Math Biol, 2022, 84(5): 1-25. DOI: 10.1007/s11538-022-01010-w .

[17] EDENBOROUGH K M, GILBERTSON B P, BROWN L E. A Mouse Model for the Study of Contact-dependent Transmission of Influenza a Virus and the Factors that Govern Transmissibility[J]. J Virol, 2012, 86(23): 12544-12551. DOI: 10.1128/JVI.00859-12 .

[18] LIN J Z, XU R, TIAN X H. Threshold Dynamics of an HIV-1 Virus Model with both Virus-to-cell and Cell-to-cell Transmissions, Intracellular Delay, and Humoral Immunity[J]. Appl Math Comput, 2017, 315: 516-530. DOI: 10.1016/j.amc.2017.08.004 .

[19] DUMORTIER F, LLIBRE J, ARTÉS J C. Qualitative theory of planar differential systems[M]. Berlin: Springer, 2006.

[20] LASALLE J P. The stability of dynamical systems[M]. Philadelphia: Society for Industrial and Applied Mathematics, 1976.

[21] GUCKENHEIMER J, HOLMES P. Nonlinear Oscillations, Dynamical Systems, and Bifurcations of Vector Fields[M]. New York, NY: Springer New York, 1983. DOI: 10.1007/978-1-4612-1140-2 .

[22] CASTILLO-CHAVEZ C, SONG B J. Dynamical Models of Tuberculosis and Their Applications[J]. Math Biosci Eng, 2004, 1(2): 361-404. DOI: 10.3934/mbe.2004.1.361 .

基金资助

国家自然科学基金(11961066)

新疆维吾尔自治区自然科学基金(2021D01E12)