高维传染病模型的Lyapunov函数构造

Lyapunov Functions for Higher-Dimensional Epidemiological Models

摘要: 本文主要研究几类传染病模型如SIR，SIRS，SIS和SEIR模型的Lypunov函数构造方法，从而获得传染病模型全局稳定性的结论。

Abstract: Lypunov functions for classical epidemiological models are introduced such as SIR, SIRS, SIS and SEIR. Global stability of some epidemiological models is also established.

1. 引言

动力系统的全局稳定性判断并不容易，最常用的方法就是直接构造函数 [1]。但是，构造过程中需要具有特殊性质的辅助函数，也就是说，Lypunov函数并不容易寻找，本文我们将集中探讨几类传染病模型的Lypunov函数构造。

2. SIR模型

借助古典假设 [2] [3]，我们将总人扣N分为以下几类：易感者类，染病者类和移除者类，分别用字母S，I和R表示，即 $N=S+I+R$。当个体被感染后，将从易感者类变为染病者类，随着个人的康复，从而转变为移除者类。当然也有可能因疾病而死亡。不妨假设总人口N为常数，即使出生率等于死亡率。若存在垂直传染即母婴等传染，假设p为传染水平，同时水平传染率记为 $\frac{\beta SI}{N}$。因此就产生了 $p\gamma I$，进去染病者类与此同时 $\gamma N-p\gamma I$，进入易感者类。

于是便产生了如下模型：

$\{\begin{array}{l}\stackrel{\dot{}}{S}=\gamma N-\beta \frac{SI}{N}-p\gamma I\\ \stackrel{\dot{}}{I}=\beta \frac{SI}{N}-\left(\delta -p\gamma \right)I\end{array}$ (1)

系统(1)有两个平衡点：疾病消除平衡点 $E_0=\left(\frac{\delta -p\gamma }{\beta }N,0\right)=\left(S_0,I_0\right)$ 和地方病平衡点 $E^{\ast }=\left(\frac{\delta -p\gamma }{\beta }N,\frac{\gamma }{\delta }N\right)=\left(S^{\ast },I^{\ast }\right)$。

定理1系统(1)的地方病平衡点 $E^{\ast }$ 为全局稳定的。

证明：

由微分等式可得：

$\beta \frac{S^{\ast }{I}^{\ast }}{N}=\gamma N-p\gamma I^{\ast }=\left(\delta -p\gamma \right)I^{\ast }$

和

$\beta \frac{S^{\ast }}{N}=\left(\delta -p\gamma \right)$, $I^{\ast }=\frac{\gamma }{\delta }N$.

构造Lyapunov函数如下

$V\left(S,I\right)=S^{\ast }\left(\frac{S}{S^{\ast }}-\ln \frac{S}{S^{\ast }}\right)+\frac{\delta }{\delta -p\gamma }{I}^{\ast }\left(\frac{I}{I^{\ast }}-\ln \frac{I}{I^{\ast }}\right)$

满足

$\frac{\partial V}{\partial S}=1-\frac{S^{\ast }}{S},\frac{{\partial }^{2}V}{\partial S^2}=\frac{S^{\ast }}{S^2}$,

$\frac{\partial V}{\partial I}=\frac{\delta }{\delta -p\gamma }\left(1-\frac{I^{\ast }}{I}\right),\frac{{\partial }^{2}V}{\partial I^2}=\frac{\delta }{\delta -p\gamma }\frac{I^{\ast }}{I^2}$

$\frac{{\partial }^{2}V}{\partial S\partial I}=0$

以及

$\left|\begin{array}{cc}\frac{{\partial }^{2}V}{\partial S^2}& \frac{{\partial }^{2}V}{\partial S\partial I}\\ \frac{{\partial }^{2}V}{\partial S\partial I}& \frac{{\partial }^{2}V}{\partial I^2}\end{array}\right|=\left(\frac{{\partial }^{2}V}{\partial S^2}\frac{{\partial }^{2}V}{\partial I^2}-\frac{{\partial }^{2}V}{\partial S\partial I}\right)=\frac{1}{S^2{I}^2}\frac{N\gamma }{\beta }>0$

很容易看出 $E^{\ast }$ 为唯一极值状态且为函数 $V\left(S,I\right)$ 在 $R_{+}^2$ 内的最小值。

于是

$\begin{array}{c}\stackrel{\dot{}}{V}\left(S,I\right)=\left(1-\frac{S^{\ast }}{S}\right)\left(\gamma N-\beta \frac{SI}{N}-p\gamma I\right)+\left(1-\frac{I^{\ast }}{I}\right)\left(\beta \frac{SI}{N}-\left(\delta -p\gamma \right)I\right)\\ =\gamma N-\beta \frac{SI}{N}-p\gamma I-\gamma N\frac{S^{\ast }}{S}+\beta \frac{S^{\ast }I}{N}+p\gamma I\frac{S^{\ast }}{S}+\left(I-I^{\ast }\right)\left(\beta \frac{S}{N}-\beta \frac{S^{\ast }}{N}\right)\delta \frac{N}{\beta S^{\ast }}\\ =\gamma N-\beta \frac{SI}{N}-p\gamma I-\gamma N\frac{S^{\ast }}{S}+\beta \frac{S^{\ast }I}{N}+p\gamma I\frac{S^{\ast }}{S}+\left(I-I^{\ast }\right)\left(S-S^{\ast }\right)\frac{\delta }{S^{\ast }}\\ =\gamma N\left(2-\frac{S}{S^{\ast }}-\frac{S^{\ast }}{S}\right)-p\gamma I\left(2-\frac{S}{S^{\ast }}-\frac{S^{\ast }}{S}\right)\\ =\left(2-\frac{S}{S^{\ast }}-\frac{S^{\ast }}{S}\right)\left(\gamma N-p\gamma I\right)\end{array}$

由于 $2-\left(\frac{S}{S^{\ast }}+\frac{S^{\ast }}{S}\right)\le 0$，故等式 $\stackrel{\dot{}}{V}\left(S,I\right)\le 0$, $\stackrel{\dot{}}{V}\left(S,I\right)=0$ 仅在 $S=S^{\ast }$ 时才成立。由渐近稳定定理，系统(1)的地方病平衡点 $E^{\ast }$ 为全局稳定的。

3. SEIR模型

具有非线性传染率的SEIR模型描述如下：

$\begin{array}{l}{S}^{\prime }=-\lambda I^p{S}^q+b-\mu S\\ E^{\prime }=\lambda I^p{S}^q-\left(\epsilon +\mu \right)E\\ I^{\prime }=\epsilon E-\left(\gamma +\mu \right)I\\ R^{\prime }=\gamma I-\mu R\end{array}$ (2)

其中参数 $p,q,\epsilon ,\mu ,\lambda ,\gamma ,b$ 为正数，总人口N分为以下几类：易感者类，潜伏者类，染病者类和移除者类，分别用字母S，E，I和R表示，即 $N=S+E+I+R$。非负常数 $\mu$ 表示死亡率，b为出生率(假设出生率与死亡率相等)， $\epsilon$ 表示潜伏者类到染病者类的转化率， $\gamma$ 为康复率。

当 $p=1$ 时，参数 $\sigma =\frac{\lambda \epsilon }{\left(\epsilon +\mu \right)\left(\gamma +\mu \right)}$ 表示接触数。

定理 2 [4] 如果 $0<p<1$ 或 $p=1$ 且 $\sigma >1$，地方病平衡点在区域T内是全局渐近稳定的。( $T=\left\{\left(S,E,I\right):0\le S,E,I\le 1,S+E+I\le 1\right\}$ )

当 $p>1$ 时，令 $p=2,q=1,\tau =\left(\gamma +\mu \right)t,\alpha =\frac{\mu }{\gamma +\mu },\beta =\frac{\epsilon }{\gamma +\mu },a=\frac{\lambda }{\gamma +\mu }$

系统(2)可化为

$\begin{array}{l}{S}^{\prime }=-a{I}^{2}S+\alpha -\alpha S\\ E^{\prime }=a{I}^{2}S-\left(\alpha +\beta \right)E\\ I^{\prime }=\beta E-I\end{array}$ (3)

定理3如果 $p=2,q=1$ 且 $\bar{\sigma }=\frac{a\alpha {\beta }^2}{4{\left(\alpha +\beta \right)}^2}<1$，t则只有唯一疾病消除平衡点 $P_0\left(1,0,0\right)$ 且在区域T内全局渐近稳定( $T=\left\{\left(S,E,I\right):0\le S,E,I\le 1,S+E+I\le 1\right\}$ )。

证明：构造Lyapunov函数如下

$V\left(S,E,I\right)=\left(S-S^{\ast }-S^{\ast }\ln \frac{S}{S^{\ast }}\right)+\frac{\alpha +\beta }{\beta }I\left(1+\frac{1}{p-1}{\left(\frac{I^{\ast }}{I}\right)}^p\right)+\left(E-E^{\ast }\ln E\right)$

则

$V\left(S,E,I\right)=\left(S-S^{\ast }-S^{\ast }\ln \frac{S}{S^{\ast }}\right)+\frac{\alpha +\beta }{\beta }I+E$

满足 $V\left(S^{\ast },E^{\ast },I^{\ast }\right)=0$ 当 $S,E,I>0$ 时， $\left(S^{\ast },E^{\ast },I^{\ast }\right)=\left(1,0,0\right)$

通过计算可得

$\frac{\partial V}{\partial S}=1-\frac{1}{S}$, $\frac{\partial V}{\partial E}=1$, $\frac{\partial V}{\partial I}=\frac{\alpha +\beta }{\beta }$

以及

$\begin{array}{c}{V}^{\prime }=\left(1-\frac{1}{S}\right)\left(-a{I}^{2}S+\alpha -\alpha S\right)+\left[a{I}^{2}S-\left(\alpha +\beta \right)E\right]+\frac{\alpha +\beta }{\beta }\left(\beta E-I\right)\\ =a{I}^2-\frac{\alpha +\beta }{\beta }I+\alpha \left(1-S\right)\left(1-\frac{1}{S}\right)\le 0\end{array}$

可得 $\left(1-S\right)\left(1-\frac{1}{S}\right)\le 0$ 仅当 $S=1$ 时；

令 $R_0=\frac{\lambda \epsilon }{\left(\gamma +\mu \right)\left(\epsilon +\mu \right)}=\frac{a\beta }{\alpha +\beta }$， $f\left(I\right)=a{I}^2-\frac{\alpha +\beta }{\beta }I$，可得当 $I\in \left(0,\frac{1}{R_0}\right)$ 时， $R_0<1$ 且 $f\left(I\right)=a{I}^2-\frac{\alpha +\beta }{\beta }I<0$。

由渐近稳定定理可得，疾病消除平衡点 $P_0\left(1,0,0\right)$ 在区域T内全局稳定。

定理4如果 $\bar{\sigma }=\frac{a\alpha {\beta }^2}{4{\left(\alpha +\beta \right)}^2}=1$，系统产生周期轨道。

定理5系统(3)的非常数周期解的轨迹，如果存在，则必为渐近相轨道渐近稳定。

证. 系统(2)的解 $\left(S\left(t\right),E\left(t\right),I\left(t\right)\right)$ 的线性系统为

$\begin{array}{l}{X}^{\prime }=-\left(a{I}^2+2\alpha +\beta \right)X+2aIS\left(Y+Z\right)\\ Y^{\prime }=\beta X-\left(a{I}^2+\alpha +1\right)Y\\ Z^{\prime }=a{I}^{2}Y-\left(\alpha +\beta +1\right)Z\end{array}$ (4)

想要讨论系统(4)的稳定性，需要构造下列函数

$V\left(X,Y,Z;S,E,I\right)=\sup \left\{\left|X\right|,\frac{E}{I}\left(\left|Y\right|+\left|Z\right|\right)\right\}$

假定解 $\left(S\left(t\right),E\left(t\right),I\left(t\right)\right)$ 具有最小周期 $\omega >0$，轨道 $\gamma$ 与T的边界始终有正数的距离，那么会存在常数c使得

$V\left(X,Y,Z;S,E,I\right)\ge c\sup \left\{\left|X\right|,\left|Y\right|,\left|Z\right|\right\}$

对 $\left(X,Y,Z\right)\in R_{+}^3$ 和 $\left(S,E,I\right)\in \gamma$。

通过计算可得 $V\left(t\right)$ 的右导数

$\begin{array}{c}{D}_{+}\left|X\left(t\right)\right|\le -\left(a{I}^2+2\alpha +\beta \right)\left|X\right|+2aIS\left(\left|Y\right|+\left|Z\right|\right)\\ =-\left(a{I}^2+2\alpha +\beta \right)\left|X\right|+\frac{2a{I}^{2}S}{E}\left\{\frac{E}{I}\left(\left|Y\right|+\left|Z\right|\right)\right\}\end{array}$ (5)

和

$D_{+}\left|Y\left(t\right)\right|\le \beta \left|X\right|-\left(a{I}^2+\alpha +1\right)\left|Y\right|$，

$D_{+}\left|Z\left(t\right)\right|\le a{I}^2\left|Y\right|-\left(\alpha +\beta +1\right)\left|Z\right|$，

这样便有

$\begin{array}{c}{D}_{+}\frac{E}{I}\left(\left|Y\right|+\left|Z\right|\right)=\left(\frac{E^{\prime }}{I}-\frac{I^{\prime }}{E}\right)\frac{E}{I}\left(\left|Y\right|+\left|Z\right|\right)+\frac{E}{I}{D}_{+}\left(\left|Y\right|+\left|Z\right|\right)\\ \le \left(\frac{E^{\prime }}{E}-\frac{I^{\prime }}{I}\right)\frac{E}{I}\left(\left|Y\right|+\left|Z\right|\right)+\frac{E}{I}\left(\beta \left|X\right|-\left(\alpha +1\right)\left|Y\right|+\left(\alpha +\beta +1\right)\left|Z\right|\right)\\ \le \frac{\beta E}{I}\left|X\right|+\left(\frac{E^{\prime }}{E}-\frac{I^{\prime }}{I}-\alpha -1\right)\frac{E}{I}\left(\left|Y\right|+\left|Z\right|\right)\end{array}$ (6)

联立(5) (6)可得

$D_{+}V\left(t\right)\le \sup \left\{g_1\left(t\right),g_2\left(t\right)\right\}V(\; t\; )$

其中 $g_1\left(t\right)=-\left(a{I}^2+2\alpha +\beta \right)+\frac{2a{I}^{2}S}{E}$， $g_2\left(t\right)=\frac{\beta E}{I}+\frac{E^{\prime }}{E}-\frac{I^{\prime }}{I}-\alpha -1$

结合(3)，会有

$\frac{I^{\prime }}{I}+1=\frac{\beta E}{I}$, $\frac{E^{\prime }}{E}+\alpha +\beta =\frac{a{I}^{2}S}{E}$

便可得到

$D_{+}V\left(t\right)\le \sup \left\{\frac{E^{\prime }}{E}-\alpha ,\frac{2E^{\prime }}{E}+\beta -a{I}^2\right\}V(\; t\; )$

因为 $\left(S\left(0\right),E\left(0\right),I\left(0\right)\right)\in \stackrel{0}{T}$，则存在 $\beta$，a使得 $\beta -a{I}^2<0$ 以及

${\displaystyle {\int }_0^{\omega }\sup \left\{g_1\left(t\right),g_2\left(t\right)\right\}\text{d}t}<0$

即揭示了 $V\left(t\right)\to 0$ 当 $t\to \infty$ 时，进而得到 $\left(X\left(t\right),Y\left(t\right),Z\left(t\right)\right)\to 0$ 当 $t\to \infty$ 事。因此，系统(4)渐近稳定，周期解 $\left(S\left(t\right),E\left(t\right),I\left(t\right)\right)$ 以渐近相轨道渐近稳定。

定理6如果 $\bar{\sigma }=\frac{a\alpha {\beta }^2}{4{\left(\alpha +\beta \right)}^2}>1$，系统(2.1)会产生hopf分支为稳定开关想象。

4. 结论

Lyapunov函数的构造方法同样也适用于竞争捕食系统 [5] [6]。本文考虑了具有非线性传染率的SIR和SEIR传染病模型，此外，具有非线性传染率的SIRS，SEIRS模型，也可以用同样的构造方法来解决全局稳定性。

致谢

作者对同行评阅人的意见和建议表示深深的感谢。

基金项目

本文由2020大学生创新创业训练项目(X202010452127)支持。