一类具有两菌株突变的传染病模型

A Class of Infectious Disease Models with Mutations in Two Strains

摘要: 病原体在传播过程中所引起的变异严重影响对于疾病传播动力学的分析。为了更好地了解病原体突变对疾病传播的影响，本文建立了一个具有两菌株在传播过程中发生变异的传染病模型。我们主要用常微分方程描述了宿主内两种菌株变异的传染病模型，采用在模型平衡点处线性化的方法，分析了该模型的平衡点以及平衡点附近的稳定性。我们得出在模型在无病平衡点处是全局渐近稳定的，模型在边界平衡点处是局部渐近稳定的，最后我们给出了模型共存平衡点存在的条件。

Abstract: The variation caused by the pathogen during transmission has a serious impact on the analysis of disease transmission dynamics. In order to better understand the effects of pathogens mutation on the spread of disease, the article established a mutation in the process of two strains in the spread of infectious disease model. We describe an infectious disease model with variation of two strains in the host mainly by ordinary differential equations, and analyze the stability of the model at and near the equilibrium point by linearizing the equilibrium point of the model. We know that the model is globally asymptotically stable at the disease-free equilibrium point and locally asymp-totically stable at the boundary equilibrium point. Finally, we give the conditions for the existence of coexistence equilibrium.

1. 引言

最近，对于流行病的研究引起了许多学者的广泛关注，理解物种多样化的机制仍然是进化生态学的核心问题，为了制定合理的控制措施，数学模型被提出。数学模型对于研究病原体的变异以及分析疾病传播的动力学具有重要意义。并且数学模型被用于研究由RNA病毒引起的疾病，特别是HIV病毒，已有多年的历史 [1]。突变是微生物DNA或RNA序列的变化，是由DNA或RNA复制过程中的错误引起的。使用RNA作为基因物质的微生物(例如病毒)具有快速变异率，这对于病毒的生存可能是一个优势，因为这些病原体不断迅速进化，形成不同的抗原特性，从而能够逃避人类免疫系统的防御反应。

由于寄生虫的相互作用，人类、动物和植物可以被同一种寄生虫的不同基因型多重感染，在传染病模型中，突变可以通过一个术语来解释，也就是将感染其中一种菌株的个体转变为感染另一种菌株的个体 [2]。一般情况下，病毒RNA聚合酶缺乏 [3] 的验证能力，因此RNA病毒的突变率高于DNA病毒。在自然选择条件下，RNA病毒较短的生成时间和较高的突变率可以帮助其快速适应宿主环境的变化。

然而，突变可能会改变之前的竞争平衡。其中文献 [4] [5] 普遍认为，免疫上调将使宿主付出代价，主要包括能量和合成代谢过程，而这反过来又可能降低宿主的生存率和繁殖率。不同宿主的免疫上调程度可能不同，这可能导致病毒在宿主内发生变异。有关文献 [6] [7] [8] [9] 还表明，在感染过程中控制传播模式和毒力的宿主特性可以极大地影响病原体进化的毒力。这些因素可能是与免疫相关的病理，其动机是受一种特定系统的驱动，即宿主感染后的传播非常迅速，但病理发展较慢，许多新出现的人类和野生动物疾病都具有这种模式 [10] [11] [12] [13]。因此研究菌株变异模型可以更好地了解菌株进化后的毒性强弱等问题。

病原体突变主要为毒性的减弱或者增强，体现在 [3] [14] [15] [16] 祖老师等提出的捕食食饵模型当中，捕食者的捕食攻击能力可以适应进化，被捕食者的躲避能力也会随捕食者的捕食能力适应性进化，即可能会共同进化，也可能会产生进化分支，并且发现当捕食者物种表现出种内竞争(即干扰竞争)时，攻击能力的适应性进化可以促进捕食者物种的多样性。

本文其它部分构成如下：在第二节中，我们建立一个具有两种病毒变异的免疫流行病学模型，第三节中我们得到了模型无病平衡点的全局稳定性、唯一边界平衡点的存在性和稳定性，以及两种菌株共存平衡点的存在性。最后，我们对我们的成果进行了总结。

2. 两菌株变异模型

2.1. 模型提出

为了更好地解释菌株变异过程，我们建立如下传染病模型 [17]。

$\{\begin{array}{l}\frac{\text{d}S}{\text{d}t}=\Lambda -\mu S-{\beta }_{1}SI-{\beta }_{2}SJ\\ \frac{\text{d}I}{\text{d}t}={\beta }_{1}SI-\left(\mu +{\alpha }_1+m\right)I\\ \frac{\text{d}J}{\text{d}t}={\beta }_{2}SJ-\left(\mu +{\alpha }_2\right)J+mI\end{array}$ (1)

在系统(1)中，S是易感者， $\Lambda$ 表示种群出生率， $\mu$ 表示自然死亡率，I和J分别表示两种菌株， ${\alpha }_1$ 表示感染菌株I的因病死亡率或者表示菌株的毒性， ${\alpha }_2$ 表示感染菌株J的因病死亡率， ${\beta }_1$ 是菌株I的感染率， ${\beta }_2$ 是菌株J的感染率，当 ${\alpha }_1$ 增加时，m表示菌株I突变成为菌株J的变异率，具有不同毒性的病原体可以对相同的宿主S进行感染，因此两种不同毒性的病原体对于易感者的争夺会形成竞争机制。

此外，我们假设模型的所有参数都是非负的，并且 $\mu >0$，总人口由模型中所有易感者与感染者之和给出：

$N=S+I+J$

将系统(2.1)中的方程相加，我们得到以下方程式：

$\frac{\text{d}N}{\text{d}t}=\Lambda -\mu S-\left(\mu +{\alpha }_1\right)I-\left(\mu +{\alpha }_2\right)J\le \Lambda -\mu N(\; t\; )$

因此，模型总的人口数

$\underset{t\to \infty }{\lim }N\left(t\right)\le \frac{\Lambda }{\mu }$

由此可以得出这样的结论

$N\left(t\right)\le N_0{\text{e}}^{-\mu t}+\frac{\Lambda }{\mu }\left(1-{\text{e}}^{-\mu t}\right)$

2.2. 模型的平衡点存在及稳定性

可以看出系统(1)显然存在无病平衡点 ${\epsilon }_0^{*}=\left(S_0^{*},0,0\right)$，其中 $S_0^{*}=\frac{\Lambda }{\mu }$。我们用 ${\epsilon }_1^{*}=\left(S_1^{*},0,J_1^{*}\right)$ 表示感染第二种菌株J相对应的边界平衡点，用 ${\epsilon }^{*}=\left(S^{*},I^{*},J^{*}\right)$ 表示菌株变异后I和J两种菌株的共存平衡点。基于生物意义，首先我们定义了基本再生数 $R_1=\frac{{\beta }_1}{\mu +{\alpha }_1+m}$，基本再生数 $R_2=\frac{{\beta }_2}{\mu +{\alpha }_2}$ 它表示菌株J的进化平衡。

定理2.1如果 $R_1<1$ 且 $R_2<1$，那么无病平衡点 ${\epsilon }_0^{*}=\left(S_0^{*},0,0\right)$， $S_0^{*}=N_0^{*}=\frac{\Lambda }{\mu }$ 是局部渐近稳定的，否则 $R_1>1$ 或 $R_2>1$，那么无病平衡点 ${\epsilon }_0^{*}$ 是不稳定的。

证明：让 $S=S_0^{*}+x\left(t\right),I=y\left(t\right),J=z\left(t\right)$，通过将系统(2.1)在无病稳态周围线性化，那么我们可以得到：

$\{\begin{array}{l}\frac{\text{d}x\left(t\right)}{\text{d}t}=-\mu x\left(t\right)-{\beta }_1{S}_0^{*}y\left(t\right)-{\beta }_2{S}_0^{*}z\left(t\right)\\ \frac{\text{d}y\left(t\right)}{\text{d}t}={\beta }_1{S}_0^{*}y\left(t\right)-\left(\mu +{\alpha }_1+m\right)y\left(t\right)\\ \frac{\text{d}z\left(t\right)}{\text{d}t}={\beta }_2{S}_0^{*}z\left(t\right)-\left(\mu +{\alpha }_2\right)z\left(t\right)+my\left(t\right)\end{array}$ (2)

为了分析 ${\epsilon }_0^{*}$ 附近的渐近行为，令 $x\left(t\right)=x{\text{e}}^{\lambda t},y\left(t\right)=y\left(t\right){\text{e}}^{\lambda t},z\left(t\right)=z{\text{e}}^{\lambda t}$

那么，系统(2)变为

$\{\begin{array}{l}\lambda x=-\mu x-{\beta }_1{S}_0^{*}y-{\beta }_2{S}_0^{*}z\\ \lambda y={\beta }_1{S}_0^{*}y-\left(\mu +{\alpha }_1+m\right)y\\ \lambda z={\beta }_2{S}_0^{*}z-\left(\mu +{\alpha }_2\right)z+my\end{array}$ (3)

得到系统(3)的雅可比矩阵

$D=\left|\begin{array}{ccc}-\mu & -{\beta }_1{S}_0^{*}& -{\beta }_2{S}_0^{*}\\ 0& {\beta }_1{S}_0^{*}-\left(\mu +{\alpha }_1+m\right)& 0\\ 0& m& {\beta }_2{S}_0^{*}-\left(\mu +{\alpha }_2\right)\end{array}\right|$ (4)

因此我们可以得到特征方程

$\begin{array}{c}\left|\lambda E-D\right|=\left|\begin{array}{ccc}\lambda +\mu & {\beta }_1{S}_0^{*}& {\beta }_{2}S\\ 0& \lambda -{\beta }_1{S}_0^{*}+\left(\mu +{\alpha }_1+m\right)& 0\\ 0& -m& \lambda -{\beta }_2{S}_0^{*}+\left(\mu +{\alpha }_2\right)\end{array}\right|\\ =\left(\lambda +\mu \right)\left(\lambda +\mu +{\alpha }_1+m-{\beta }_1\right)\left(\lambda +\mu +{\alpha }_2-{\beta }_2\right)=0\end{array}$ (5)

显然，由上式特征方程得出 ${\lambda }_1=-\mu$， ${\lambda }_2={\beta }_1-\left(\mu +{\alpha }_1+m\right)$， ${\lambda }_3={\beta }_2-\left(\mu +{\alpha }_2\right)$，那么它的稳定性取决于 ${\lambda }_2$ 和 ${\lambda }_3$，可以看出当 $R_1<1$ 时 ${\lambda }_2<0$，且 $R_2<1$ 时 ${\lambda }_3<0$ 时，由Routh-Hurwitz准则可以得出，当 $R_1<1$ 且 $R_2<1$ 时，无病平衡点 ${\epsilon }_0^{*}$ 是局部渐近稳定的。

定理2.2如果 $R_2<1$ 且 $R_3=\frac{{\beta }_1}{\mu +{\alpha }_1}<1$，那么无病平衡点 ${\epsilon }_0^{*}$ 是全局渐近稳定的，其中 $S_0^{*}=N_0^{*}=\frac{\Lambda }{\mu }$。

证明：下面我们构造Lyapunov函数验证无病平衡点的全局稳定性，令

$V\left(t\right)=I\left(t\right)+J\left(t\right)$ (6)

然后计算函数 $V\left(t\right)$ 沿系统(1)的时间导数，我们得到

$\begin{array}{c}\frac{\text{d}V\left(t\right)}{\text{d}t}=\frac{\text{d}I\left(t\right)}{\text{d}t}+\frac{\text{d}J\left(t\right)}{\text{d}t}\\ ={\beta }_{1}SI-\left(\mu +{\alpha }_1+m\right)I+{\beta }_{2}SJ-\left(\mu +{\alpha }_2\right)J+mI\\ =\left({\beta }_1-\mu -{\alpha }_1\right)I+\left({\beta }_2-\mu -{\alpha }_2\right)J\end{array}$ (7)

由此看出满足 $\frac{\text{d}V}{\text{d}t}<0$， $R_2<1$ 且 $R_3<1$，所以等式 $\frac{\text{d}V}{\text{d}t}=0$ 成立的充要条件是： $I\left(t\right)=0,J\left(t\right)=0$。因此无病平衡态是全局渐近稳定的，即在人群中无论初始感染者有多少，疾病都不会再流行而是逐渐消失。

2.3. 模型的边界平衡点存在及稳定性

定理2.3当 $R_2>\frac{\mu }{\Lambda }$ 时候，边界平衡点 ${\epsilon }_1^{*}=\left(S_1^{*},0,J_1^{*}\right)$ 存在。

证明：系统(1)变为

$\{\begin{array}{l}0=\Lambda -\mu S_1^{*}-{\beta }_2{S}_1^{*}{J}_1^{*}\\ 0={\beta }_2{S}_1^{*}{J}_1^{*}-\left(\mu +{\alpha }_2\right)J_1^{*}\end{array}$ (8)

同样的，由(8)第二个式子可以得到

$s_1^{*}=\frac{\mu +{\alpha }_2}{{\beta }_2}=\frac{1}{R_2}$ (9)

将(9)带入(8)的第一个式子得到：

$j_1^{*}=\frac{\Lambda R_2-\mu }{{\beta }_2}$ (10)

因此，如果 $R_2>\frac{\mu }{\Lambda }$，可以得出存在边界平衡态 ${\phi }_1^{*}=\left(s_1^{*},0,j_1^{*}\right)$ 存在，也就是边界平衡点 ${\epsilon }_1^{*}=\left(S_1^{*},0,J_1^{*}\right)$ 存在。

定理2.4当 $R_2>\max \left\{R_1,\frac{\mu }{\Lambda }\right\}$ 时，如系统恰好存在一个简单的地方病平衡，那么边界平衡点 ${\epsilon }_1^{*}=\left(S_1^{*},0,J_1^{*}\right)$ 是局部渐近稳定的，否则是不稳定的。

证明：用类似的方法将系统(1)在边界平衡点 ${\epsilon }_1^{*}=\left(S_1^{*},0,J_1^{*}\right)$ 的处线性化，那么我们得到以下雅可比矩阵：

$T=\left|\begin{array}{ccc}-\mu -{\beta }_2{j}_1^{*}& -{\beta }_1{s}_1^{*}& -{\beta }_2{s}_1^{*}\\ 0& {\beta }_1{s}_1^{*}-\left(\mu +{\alpha }_1+m\right)& 0\\ {\beta }_2{j}_1^{*}& m& {\beta }_2{s}_1^{*}-\left(\mu +{\alpha }_2\right)\end{array}\right|$ (11)

因此可以得到如下特征方程：

$\begin{array}{c}\left|\lambda E-T\right|=\left|\begin{array}{ccc}\lambda +\mu +{\beta }_2{j}_1^{*}& {\beta }_1{s}_2^{*}& {\beta }_2{s}_2^{*}\\ 0& \lambda -{\beta }_1{s}_2^{*}+\left(\mu +{\alpha }_1+m\right)& 0\\ -{\beta }_2{j}_1^{*}& -m& \lambda -{\beta }_2{s}_2^{*}+\left(\mu +{\alpha }_2\right)\end{array}\right|\\ =\left(\lambda -{\beta }_1{s}_2^{*}+\mu +{\alpha }_1+m\right)\left|\begin{array}{cc}\lambda +\mu +{\beta }_2{j}_1^{*}& {\beta }_2{s}_2^{*}\\ -{\beta }_2{j}_1^{*}& \lambda -{\beta }_2{s}_2^{*}+\left(\mu +{\alpha }_2\right)\end{array}\right|\\ =\left(\lambda -{\beta }_1{s}_2^{*}+\mu +{\alpha }_1+m\right)T\left(\lambda \right)\end{array}$ (12)

因此我们可以得到关于 $\lambda$ 的函数 $T(\; \lambda \; )$

$T\left(\lambda \right)={\lambda }^2+\left(2\mu +{\alpha }_2+{\beta }_2{j}_1^{*}-{\beta }_2{s}_1^{*}\right)\lambda +\left(\mu +{\beta }_2{j}_1^{*}\right)\left(\mu +{\alpha }_2-{\beta }_2{s}_1^{*}\right)+{\beta }_{{}_2}^2{j}_1^{*}{s}_1^{*}$

因此我们可以得到特征根

${\lambda }_1={\beta }_1{s}_1^{*}-\mu -{\alpha }_1-m$

如果 ${\lambda }_1<0$，也就是

$\begin{array}{l}\frac{{\beta }_1}{R_2}-\mu -{\alpha }_1-m<0\\ -\mu -{\alpha }_1-m<-\frac{{\beta }_1}{R_2}\\ \frac{{\beta }_1}{\mu +{\alpha }_1+m}<R_2\end{array}$

也就是

$R_1<R_2$

此时在 $R_1<R_2$ 条件下，边界平衡点的稳定性取决于 $T\left(\lambda \right)$ 的根实部的符号，也就是

$a_1=2\mu +{\alpha }_2+{\beta }_2{J}^{*}-{\beta }_2{S}^{*}$

以及

$a_2={\mu }^2+\mu {\alpha }_2+\mu {\beta }_2{J}^{*}+{\alpha }_2{\beta }_2{J}^{*}-\mu {\beta }_2{S}^{*}-{\beta }_2{J}^{*}{S}^{*}+{\beta }_2^2{J}^{*}{S}^{*}$

代入可以得到

$a_1=\Lambda {\Re }_2$ (13)

$a_2=\left(\Lambda {\Re }_2-\mu \right)\left(\mu +{\alpha }_2\right)$ (14)

显然 $a_1>0$，因此当 $R_2>\frac{\mu }{\Lambda }$ 时， $a_2>0$，所以由Routh-Hurwitz判据准则知，当 $R_2>\frac{\mu }{\Lambda }$ 时，边界平衡点 ${\epsilon }_1^{*}=\left(S_1^{*},0,J_1^{*}\right)$ 显然是稳定的，否则是不稳定的。

3. 模型共存平衡点的存在性

定理2.5当 $R_1>\max \left\{R_2,\frac{\mu }{\Lambda },\frac{\mu }{{\beta }_1}\right\}$ 时候，共存平衡点 ${\epsilon }^{*}=\left(S^{*},I^{*},J^{*}\right)$ 存在。

证明：同样的，要想证明共存平衡点的存在性，可由系统(1)得到

$\{\begin{array}{l}0=\Lambda -\mu S^{*}-{\beta }_1{S}^{*}{I}^{*}-{\beta }_2{S}^{*}{J}^{*}\\ 0={\beta }_1{S}^{*}{I}^{*}-\left(\mu +{\alpha }_1+m\right)I^{*}\\ 0={\beta }_2{S}^{*}{J}^{*}-\left(\mu +{\alpha }_2\right)J^{*}+m{I}^{*}\end{array}$ (15)

解系统(15)的第二个方程可得

$S^{*}=\frac{\mu +{\alpha }_1+m}{{\beta }_1}=\frac{1}{R_1}$

将带入的第二、三个方程得到

$I^{*}=\frac{{\beta }_2\left(\frac{1}{R_2}-\frac{1}{R_1}\right)\left(\Lambda R_1-\mu \right)}{{\beta }_1\left(\mu +{\alpha }_2\right)-{\beta }_2\left(\mu +{\alpha }_1\right)}$, $J^{*}=\frac{m\left(\Lambda R_1-\mu \right)}{{\beta }_1\left(\mu +{\alpha }_2\right)-{\beta }_2\left(\mu +{\alpha }_1\right)}$

为确保参数的生物学意义，我们假设所有参数大于零。

我们注意到当下式成立时

$R_1>\max \left\{R_2,\frac{\mu }{\Lambda },\frac{\mu }{{\beta }_1}\right\}$ (16)

此时 $I^{*}$ 和 $J^{*}$ 大于零，共存平衡点 ${\epsilon }^{*}=\left(S^{*},I^{*},J^{*}\right)$ 存在。

4. 总结

在这篇文章中，我们主要建立了一个具有菌株进化的传染病模型，随着时间t的增加，菌株I通过一定的变异率进化为菌株J，进化的菌株和原来的都可以继续感染原始宿主进而我们证明了模型无病平衡点 ${\epsilon }_0^{*}$ 的全局渐近稳定性、边界平衡点 ${\epsilon }_1^{*}$ 的存在性和局部稳定性，以及共存平衡点 ${\epsilon }^{*}$ 的存在性。当病原体入侵一个新的宿主种群时，初始毒力是宿主和病原体特性之间的复杂相互作用，可能远远不是共存稳定。因此，病原体在新的宿主群体中出现后，毒力可能增加或减少。当两个菌株个体的繁殖数量都小于1时，两株病毒被证明灭绝；否则，它们的进化将遵循竞争排斥原则，即强者最终生存 [18] [19]。除此之外，对于病原体突变的进一步研究是宿主病与原体的共同进化，一般文献 [20] [21] [22] [23] 会采用自适应动力学的方法来表述共同进化传染过程。