具有潜伏感染的SVEIR年龄结构传染病模型

An Age-Structured Infectious Disease Model of SVEIR with Latent Infection

摘要: 宿主的年龄是影响传染病传播与控制的重要因素，感染的风险随着年龄增长而增加。此外，疫苗接种针对不同年龄段的人群产生不同的效果。因此在本文中，在考虑疫苗接种的情况下，研究具有潜伏期感染的年龄结构SVEIR传染病模型，通过计算疾病的基本再生数R₀，进一步探究无病平衡点的全局渐近稳定性、地方病平衡点的存在性以及稳定性。这些结论对控制疾病的传播具有重要的实际意义。

Abstract: The age structure of the host population is an important factor affecting the transmission and control of infectious diseases, and the risk of infection increases with age. In addition, vaccination exerts varying effects across distinct age groups. Therefore, in this paper, considering vaccination, an infectious disease model of SVEIR age structure with latent infection was studied, and the global asymptotic stability of disease-free equilibrium point, the existence and stability of endemic equilibrium point were further explored by calculating the basic reproduction number of the disease. These conclusions have important practical significance for controlling the spread of the disease.

1. 引言

传染病是病原体通过空气、水源、接触、母婴垂直传播等多种方式引起人与人、动物与动物或人与动物之间相互传播的一类疾病[1]。它具有传染性和流行性，给人们的生命安全造成了巨大的威胁，给国家公共卫生和经济带来了不可磨灭的影响。因此有关部门通过研发疫苗，实行及时救助等多项手段在传染病预防和控制上取得显著成效。这其中包括彻底被消灭的天花、麻疹、乙肝等一系列传染病得到有效控制。但是随着经济全球一体化，人员流动加速，环境污染形势愈加严峻，导致传染病的迅速传播，甚至过去被控制的某些传染病也呈现出重来的趋势。文献[2]-[4]考虑了具有重复感染的传染病模型，文献[5]中通过研究得出，接种疫苗可以明显降低基本再生数，进一步有效控制疾病的传播；由于人群对同一种疾病的感染程度随着年龄的变化而不同，例如有一些疾病高发于儿童，如百日咳、腮腺炎等[6] [7]；一些疾病高发于成年人，例如麻疹、艾滋病等[8]-[10]。因此在建立传染病模型时，将年龄因素纳入考虑范围有着重要的理论意义。文献[11]研究了具有垂直接种和比例接种HBV年龄结构传染病模型，文献[8]考虑了年龄结构因素给麻疹传播带来的影响。

2. 模型建立

在这部分我们提出一个年龄结构的传染病模型来研究具有潜伏感染传染病模型的动力学性质，把在时刻$t$ ，年龄为$a$ 的总人群$N\left(a,t\right)$ 分为五类：易感人群$S\left(a,t\right)$ 、接种人群$V\left(a,t\right)$ 、潜伏人群$E\left(a,t\right)$ 、染病人群$I\left(a,t\right)$ 、恢复人群$R\left(a,t\right)$ 。采用感染力函数

$\lambda \left(a,t\right)=k\left(a\right){\displaystyle {\int }_0^{a^{+}}h\left(a^{\prime }\right)\frac{\epsilon E\left(a^{\prime },t\right)+I\left(a^{\prime },t\right)}{N\left(a^{\prime },t\right)}\text{d}{a}^{\prime }}$ (2.1)

其中$k\left(a\right)$ 为年龄$a$ 人群的易感性，$h\left(a^{\prime }\right)$ 为年龄$a^{\prime }$ 患病人群的感染力，$\epsilon$ 为潜伏期与感染期之比。

考虑传染病模型如下：

$\{\begin{array}{l}\frac{\partial S}{\partial a}+\frac{\partial S}{\partial t}=-\left(\psi \left(a\right)+\lambda \left(a,t\right)+\mu \left(a\right)\right)S\left(a,t\right),\\ \frac{\partial V}{\partial a}+\frac{\partial V}{\partial t}=\psi \left(a\right)S\left(a,t\right)-\delta \lambda \left(a,t\right)V\left(a,t\right)-\mu \left(a\right)V\left(a,t\right),\\ \frac{\partial E}{\partial a}+\frac{\partial E}{\partial t}=\delta \lambda \left(a,t\right)V\left(a,t\right)+\lambda \left(a,t\right)S\left(a,t\right)-\alpha \left(a\right)E\left(a,t\right)-\mu \left(a\right)E\left(a,t\right),\\ \frac{\partial I}{\partial a}+\frac{\partial I}{\partial t}=\alpha \left(a\right)E\left(a,t\right)-\gamma \left(a\right)I\left(a,t\right)-\mu \left(a\right)I\left(a,t\right),\\ \frac{\partial R}{\partial a}+\frac{\partial R}{\partial t}=\gamma \left(a\right)I\left(a,t\right)-\mu \left(a\right)R\left(a,t\right),\\ S(0,t)={\displaystyle {\int }_0^{a_{+}}b\left(a\right)N\left(a,t\right)\text{d}a,\begin{array}{cc}& V\left(0,t\right)=E\left(0,t\right)=I\left(0,t\right)=R\left(0,t\right)=0,\end{array}}\\ S\left(a,0\right)=S_0\left(a\right)\ge 0,V\left(a,0\right)=V_0\left(a\right)\ge 0,E\left(a,0\right)=E_0\left(a\right)\ge 0,\\ I\left(a,0\right)=I_0\left(a\right)\ge 0,R\left(a,0\right)=R_0\left(a\right)\ge 0.\end{array}$ (2.2)

其中，$\psi \left(a\right)$ 为疫苗有效接种率，$\mu \left(a\right)$ 为自然死亡率，${\left[\sigma \left(a\right)\right]}^{-1}$ 为平均潜伏期，$\gamma \left(a\right)$ 为恢复率，$\delta$ 为感染系数$\left(0\le \delta \le 1\right)$ 且所有参数和初始条件非负。

总人口$N\left(a,t\right)$ 满足：

$\{\begin{array}{l}\frac{\partial N}{\partial a}+\frac{\partial N}{\partial t}=-\mu \left(a\right)N\left(a,t\right),\\ N\left(0,t\right)={\displaystyle {\int }_0^{a_{+}}b\left(a\right)N\left(a,t\right)\text{d}a,}\\ N\left(a,0\right)=N_0\left(a\right),\\ N_0\left(a\right)=S_0\left(a\right)+V_0\left(a\right)+E_0\left(a\right)+I_0\left(a\right)+R_0\left(a\right).\end{array}$ (2.3)

对(2.3)沿特征线积分可得：

$N\left(a,t\right)=\{\begin{array}{l}N\left(0,t-a\right){\text{e}}^{-{\displaystyle {\int }_0^a\mu \left(\xi \right)\text{d}\zeta }},a<t\\ N_0\left(a-t\right){\text{e}}^{-{\displaystyle {\int }_0^t\mu \left(\xi \right)\text{d}\zeta }},a\ge t\end{array}$ (2.4)

将(2.4)的第一式代入(2.3)的第二式可得：

$N\left(0,t\right)=N\left(0,t\right){\displaystyle {\int }_0^{+\infty }b\left(a\right)e^{-{\displaystyle {\int }_0^a\mu \left(\xi \right)d\zeta }}\text{d}a}.$

因此，当且仅当人口增长率等于0时人口处于稳态，即

${\displaystyle {\int }_0^{+\infty }b\left(a\right)e^{-{\displaystyle {\int }_0^a\mu \left(\xi \right)d\zeta }}\text{d}a}=1$ .

对系统(2.2)标准化得：

$s\left(a,t\right)=\frac{S\left(a,t\right)}{N\left(a,t\right)},v\left(a,t\right)=\frac{V\left(a,t\right)}{N\left(a,t\right)},e\left(a,t\right)=\frac{E\left(a,t\right)}{N\left(a,t\right)},i\left(a,t\right)=\frac{I\left(a,t\right)}{N\left(a,t\right)},r\left(a,t\right)=\frac{R\left(a,t\right)}{N\left(a,t\right)}$

此时系统(2.2)变换如下形式：

$\{\begin{array}{l}\frac{\partial s}{\partial a}+\frac{\partial s}{\partial t}=-\left(\psi \left(a\right)+\lambda \left(a,t\right)\right)s\left(a,t\right),\\ \frac{\partial v}{\partial a}+\frac{\partial v}{\partial t}=\psi \left(a\right)s\left(a,t\right)-\delta \lambda \left(a,t\right)v\left(a,t\right),\\ \frac{\partial e}{\partial a}+\frac{\partial e}{\partial t}=\delta \lambda \left(a,t\right)v\left(a,t\right)+\lambda \left(a,t\right)s\left(a,t\right)-\alpha \left(a\right)e\left(a,t\right),\\ \frac{\partial i}{\partial a}+\frac{\partial i}{\partial t}=\alpha \left(a\right)e\left(a,t\right)-\gamma \left(a\right)i\left(a,t\right),\\ \frac{\partial r}{\partial a}+\frac{\partial r}{\partial t}=\gamma \left(a\right)i\left(a,t\right),\\ s\left(0,t\right)=1,v\left(0,t\right)=e\left(0,t\right)=i\left(0,t\right)=r\left(0,t\right)=0,\\ s\left(a,0\right)=s_0\left(a\right)\ge 0,v\left(a,0\right)=v_0\left(a\right)\ge 0,e\left(a,0\right)=e_0\left(a\right)\ge 0,\\ i\left(a,0\right)=i_0\left(a\right)\ge 0,r\left(a,0\right)=r_0\left(a\right)\ge 0.\end{array}$ (2.5)

其中

$\lambda \left(a,t\right)=k\left(a\right)H\left(t\right)$ (2.6)

$H\left(t\right)={\displaystyle {\int }_0^{a_{+}}h\left(a\right)\left(\epsilon e\left(a,t\right)+i\left(a,t\right)\right)da}$ (2.7)

因为前四维与$r\left(a,t\right)$ 无关，所以只研究前四维耦合系统即可。接下来的讨论中我们只讨论系统(2.8)：

$\{\begin{array}{l}\frac{\partial s}{\partial a}+\frac{\partial s}{\partial t}=-\left(\psi \left(a\right)+k\left(a\right)H\left(t\right)\right)s\left(a,t\right),\\ \frac{\partial v}{\partial a}+\frac{\partial v}{\partial t}=\psi \left(a\right)s\left(a,t\right)-\delta k\left(a\right)H\left(t\right)v\left(a,t\right),\\ \frac{\partial e}{\partial a}+\frac{\partial e}{\partial t}=\delta k\left(a\right)H\left(t\right)v\left(a,t\right)+k\left(a\right)H\left(t\right)s\left(a,t\right)-\alpha \left(a\right)e\left(a,t\right),\\ \frac{\partial i}{\partial a}+\frac{\partial i}{\partial t}=\alpha \left(a\right)e\left(a,t\right)-\gamma \left(a\right)i\left(a,t\right),\\ s\left(0,t\right)=1,v\left(0,t\right)=e\left(0,t\right)=i\left(0,t\right)=0,\\ s\left(a,0\right)=s_0\left(a\right)\ge 0,v\left(a,0\right)=v_0\left(a\right)\ge 0,\\ e\left(a,0\right)=e_0\left(a\right)\ge 0,i\left(a,0\right)=i_0\left(a\right)\ge 0.\end{array}$ (2.8)

3. 无病平衡点的存在性与稳定性

显然系统(2.8)存在无病平衡点$E^0=\left(s^0\left(a\right),v^0\left(a\right),0,0\right)$ ，且

$s^0\left(a\right)={\text{e}}^{-{\displaystyle {\int }_0^a\psi \left(\xi \right)d\xi }},v^0\left(a\right)={\displaystyle {\int }_0^a\psi \left(\eta \right)s^0\left(\eta \right)\text{d}\eta }$

接下来讨论无病平衡点$E^0$ 的稳定性，考虑以下形式的扰动：

$s\left(a,t\right)=s^0\left(a\right)+\overline{s}\left(a\right)e^{\omega t},v\left(a,t\right)=v^0\left(a\right)+\overline{v}\left(a\right)e^{\omega t},e\left(a,t\right)=\overline{x}\left(a\right)e^{\omega t},i\left(a,t\right)=\overline{i}\left(a\right)e^{\omega t}$

代入系统(2.8)并保留其线性部分得：

$\{\begin{array}{l}\frac{\text{d}\overline{s}}{\text{d}a}=-\left(\psi \left(a\right)+\omega \right)\overline{s}\left(a\right)-k\left(a\right)H_1{s}^0\left(a\right),\\ \frac{\text{d}\overline{v}}{\text{d}a}=-\omega \overline{v}\left(a\right)-\delta k\left(a\right)H_1{v}^0\left(a\right)+\psi \left(a\right)\overline{s}\left(a\right),\\ \frac{\text{d}\overline{x}}{\text{d}a}=-\left(\alpha \left(a\right)+\omega \right)\overline{x}\left(a\right)+\delta k\left(a\right)H_1{v}^0\left(a\right)+k\left(a\right)H_1{s}^0\left(a\right),\\ \frac{\text{d}\overline{i}}{\text{d}a}=-\left(\gamma \left(a\right)+\omega \right)\overline{i}\left(a\right)+\alpha \left(a\right)\overline{x}\left(a\right),\\ \overline{s}\left(0\right)=\overline{v}\left(0\right)=\overline{x}\left(0\right)=\overline{i}\left(0\right)=0,\\ H_1={\displaystyle {\int }_0^{a_{+}}h\left(a\right)\left(\epsilon \overline{x}\left(a\right)+\overline{i}\left(a\right)\right)\text{d}a.}\end{array}$ (3.1)

求解系统(3.1)得：

$\{\begin{array}{l}\overline{s}\left(a\right)=H_1{f}_{1,\omega }\left(a\right),\\ \overline{v}\left(a\right)=H_1{f}_{2,\omega }\left(a\right),\\ \overline{x}\left(a\right)=H_1{f}_{3,\omega }\left(a\right),\\ \overline{i}\left(a\right)=H_1{f}_{4,\omega }\left(a\right).\end{array}$

$f_{1,\omega }\left(a\right)=-{\displaystyle {\int }_0^{a}k\left(\eta \right)s^0\left(\eta \right){\text{e}}^{-{\displaystyle {\int }_{\eta }^a\left(\psi \left(\xi \right)+\omega \right)d\xi }}\text{d}\eta },$

$f_{2,\omega }\left(a\right)={\displaystyle {\int }_0^a\left(-\delta k\left(\eta \right)v^0\left(\eta \right)+\psi \left(\eta \right)f_{1,\omega }\left(\eta \right)\right){\text{e}}^{\omega \left(\eta -a\right)}\text{d}\eta },$

$f_{3,\omega }\left(a\right)={\displaystyle {\int }_0^a\left(\delta k\left(\eta \right)v^0\left(\eta \right)+k\left(\eta \right)s^0\left(\eta \right)\right){\text{e}}^{-{\displaystyle \underset{\eta }{\overset{a}{\int }}\left(\alpha \left(\xi \right)+\omega \right)d\xi }}\text{d}\eta },$ (3.2)

$f_{4,\omega }\left(a\right)={\displaystyle {\int }_0^a\alpha \left(\eta \right)v^0\left(\eta \right)}{f}_{3,\omega }\left(\eta \right){\text{e}}^{-{\displaystyle \underset{\eta }{\overset{a}{\int }}\left(\gamma \left(\xi \right)+\omega \right)}}\text{d}\eta .$ (3.3)

将(3.2)和(3.3)代入$H_1$ 的表达式，两边同时消去$H_1$ $\left(H_1\ne 0\right)$ 得：

$1={\displaystyle {\int }_0^{+\infty }h\left(a\right)\left(\epsilon f_{3,\omega }\left(a\right)+f_{4,\omega }\left(a\right)\right)\text{d}a}:=F\left(\omega \right)$ (3.4)

记基本再生数$R_0=F\left(0\right)$ ，即

$R_0={\displaystyle {\int }_0^{+\infty }h\left(a\right)\left(\epsilon f_{3,0}\left(a\right)+f_{4,0}\left(a\right)\right)\text{d}a}$ (3.5)

其中，

$f_{3,0}\left(a\right)={\displaystyle {\int }_0^a\left(\delta k\left(\eta \right)v^0\left(\eta \right)+k\left(\eta \right)s^0\left(\eta \right)\right){\text{e}}^{-{\displaystyle \underset{\eta }{\overset{a}{\int }}\alpha \left(\xi \right)}}\text{d}\eta },$

$f_{4,0}\left(a\right)={\displaystyle {\int }_0^a\alpha \left(\eta \right)v^0\left(\eta \right)}{f}_{3,0}\left(\eta \right){\text{e}}^{-{\displaystyle \underset{\eta }{\overset{a}{\int }}\gamma \left(\xi \right)}}\text{d}\eta .$

表示一个感染者在其平均感染周期$\left[0,a^{+}\right)$ 内所产生的所有继发感染者总人数。

从而有以下定理：

定理3.1 若$R_0<1$ ，系统(2.8)的无病平衡点$E^0=\left(s^0\left(a\right),v^0\left(a\right),0,0\right)$ 是局部渐近稳定的；如果$R_0>1$ ，

则无病平衡点$E^0$ 不稳定。

证明：若$\omega$ 是实数，则有

$F^{\prime }(\omega )<0,\underset{x\to +\infty }{\lim }F(\omega )=0,\underset{x\to -\infty }{\lim }F(\omega )=+\infty .$

又特征方程(3.4)仅存在唯一实根${\omega }^{*}$ 满足

$F\left({\omega }^{*}\right)=1.$ (3.6)

所以当$R_0<1$ 时，(3.4)存在唯一负实根${\omega }^{*}<0$ ，此时无病平衡点$E^0$ 局部渐近稳定；当$R_0>1$ 时，${\omega }^{*}>0$ ，此时$E^0$ 不稳定。

若$\omega$ 是复数，令$\omega =x+iy$ 。由(3.4)得

$\mathrm{Re}F\left(\omega \right)=1,\mathrm{Im}F\left(\omega \right)=0.$ (3.7)

又$\mathrm{Re}{\text{e}}^{\omega }\le {\text{e}}^{\mathrm{Re}\omega }$ ，根据(3.6)和(3.7)得

$F\left({\omega }^{*}\right)=1=\mathrm{Re}F\left(\omega \right)\le F\left(\mathrm{Re}\omega \right)$

易得$\mathrm{Re}\omega \le {\omega }^{*}$ 。因为$F\left(x\right)$ 是关于$x$ 的单调递减函数，所以$x\le {\omega }^{*}$ ，即${\omega }^{*}$ 是特征方程的主特征根。所以当$R_0<1$ 时，$x\le {\omega }^{*}<0$ ，无病平衡点$E^0$ 局部渐近稳定，当$R_0>1$ 时$E^0$ 不稳定。

定理3.2 若$R_0<1$ ，则系统(2.8)的无病平衡点$E^0=\left(s^0\left(a\right),v^0\left(a\right),0,0\right)$ 全局渐近稳定。

证明：对系统(2.8)，沿特征线积分可得：

$s\left(a,t\right)=\{\begin{array}{l}{s}_0\left(a-t\right){\text{e}}^{-{\displaystyle {\int }_0^t\left(\psi \left(a-t+\xi \right)+k\left(a-t+\xi \right)H\left(\xi \right)\right)\text{d}\xi }},a\ge t\\ {\text{e}}^{-{\displaystyle {\int }_0^a\left(\psi \left(\xi \right)+k\left(\xi \right)H\left(t-a+\xi \right)\right)\text{d}\xi }},a<t\end{array}$ (3.8)

$v\left(a,t\right)=\{\begin{array}{l}{v}_0\left(a-t\right){\text{e}}^{-{\displaystyle {\int }_0^t\delta k\left(a-t+\xi \right)H\left(\xi \right)\text{d}\xi }}+{\displaystyle {\int }_0^t\psi \left(a-t+\eta \right)s\left(a-t+\eta ,\eta \right){\text{e}}^{-{\displaystyle {\int }_{\eta }^t\delta k\left(a-t+\xi \right)H\left(\xi \right)\text{d}\xi }}\text{d}\eta },a\ge t\\ {\displaystyle {\int }_0^a\psi \left(\eta \right)s\left(\eta ,t-a+\eta \right){\text{e}}^{-{\displaystyle {\int }_{\eta }^a\delta k\left(\xi \right)H\left(t-a+\xi \right)\text{d}\xi }}\text{d}\eta },a<t\end{array}$

$e\left(a,t\right)=\{\begin{array}{l}{e}_0\left(a-t\right){\text{e}}^{-{\displaystyle {\int }_0^t\alpha \left(a-t+\xi \right)\text{d}\xi }}+{\displaystyle {\int }_0^{t}k\left(a-t+\eta \right)H\left(\eta \right)\left(\delta v\left(a-t+\eta ,\eta \right)+\left(s\left(a-t+\eta ,\eta \right)\right)\right){\text{e}}^{-{\displaystyle {\int }_{\eta }^t\alpha \left(a-t+\xi \right)\text{d}\xi }}\text{d}\eta },a\ge t\\ {\displaystyle {\int }_0^{a}k\left(\eta \right)H\left(t-a+\eta \right)\left(\delta v\left(\eta ,t-a+\eta \right)+s\left(\eta ,t-a+\eta \right)\right){\text{e}}^{-{\displaystyle {\int }_{\eta }^a\alpha \left(\xi \right)\text{d}\xi }}\text{d}\eta },a<t\end{array}$

$i\left(a,t\right)=\{\begin{array}{l}{i}_0\left(a-t\right){\text{e}}^{-{\displaystyle {\int }_0^t\gamma \left(\xi \right)\text{d}\xi }}+{\displaystyle {\int }_0^t\left(\alpha \left(a-t+\eta \right)e\left(a-t+\eta ,\eta \right)\right){\text{e}}^{-{\displaystyle {\int }_{\eta }^t\gamma \left(\xi \right)\text{d}\xi }}\text{d}\eta },a\ge t\\ {\displaystyle {\int }_0^a\alpha \left(\eta \right)\text{e}\left(\eta ,t-a+\eta \right){\text{e}}^{-{\displaystyle {\int }_{\eta }^a\gamma \left(\xi \right)\text{d}\xi }}\text{d}\eta },a<t\end{array}$

因为我们研究的是时间趋于无穷的解，且$\lambda \left(a,t\right)\ge 0$ ，所以对于(3.6)的第二式有：

$s\left(\eta ,t-a+\eta \right)={\text{e}}^{-{\displaystyle {\int }_0^{\eta }\left(\psi \left(\xi \right)+\lambda \left(\xi ,t-a+\xi \right)\right)\text{d}\xi }}\le {\text{e}}^{-{\displaystyle {\int }_0^{\eta }\psi \left(\xi \right)\text{d}\xi }}=s^0\left(\eta \right)$ (3.9)

同理可得

$v\left(\eta ,t-a+\eta \right)\le {\displaystyle {\int }_0^{\eta }\psi \left(\sigma \right)s^0\left(\sigma \right)\text{d}\sigma }=v^0\left(\eta \right)$ (3.10)

$e\left(a,t\right)\le {\displaystyle {\int }_0^{a}H\left(t-a+\eta \right)\left(\delta k\left(\eta \right)v^0\left(\eta \right)+k\left(\eta \right)s^0\left(\eta \right)\right){\text{e}}^{-{\displaystyle {\int }_{\eta }^a\alpha (\xi )\text{d}\xi }}\text{d}\eta }$ (3.11)$i\left(a,t\right)\le {\displaystyle {\int }_0^a\alpha \left(\eta \right){\displaystyle {\int }_0^{\eta }H\left(t-\eta +\sigma \right)\left(\delta k\left(\sigma \right)v^0\left(\sigma \right)+k\left(\sigma \right)s^0\left(\sigma \right)\right){\text{e}}^{-{\displaystyle {\int }_{\sigma }^{\eta }\alpha \left(\xi \right)\text{d}\xi }}\text{d}\sigma {\text{e}}^{-{\displaystyle {\int }_{\eta }^a\gamma \left(\xi \right)\text{d}\xi }}\text{d}\eta }}$ (3.12)

由(2.7)得

$H\left(t\right)={\displaystyle {\int }_0^{t}h\left(a\right)\left(\epsilon e\left(a,t\right)+i\left(a,t\right)\right)\text{d}a}+{\displaystyle {\int }_t^{+\infty }h\left(a\right)\left(\epsilon e\left(a,t\right)+i\left(a,t\right)\right)\text{d}a}$ (3.13)

当$t$ 足够大时，(3.13)的第二部分为0。通过(3.9)~(3.12)可得

$\begin{array}{c}H\left(t\right)={\displaystyle {\int }_0^{t}h\left(a\right)\left(\epsilon e\left(a,t\right)+i\left(a,t\right)\right)\text{d}a}\\ \le {\displaystyle {\int }_0^{t}h\left(a\right)\epsilon {\displaystyle {\int }_0^{a}H\left(t-a+\eta \right)\left(\delta k\left(\eta \right)v^0\left(\eta \right)+k\left(\eta \right)s^0\left(\eta \right)\right){\text{e}}^{-{\displaystyle {\int }_{\eta }^a\alpha \left(\xi \right)\text{d}\xi }}\text{d}\eta }}\\ +{\displaystyle {\int }_0^a\alpha \left(\eta \right){\displaystyle {\int }_0^{\eta }H\left(t-\eta +\sigma \right)\left(\delta k\left(\sigma \right)v^0\left(\sigma \right)+k\left(\sigma \right)s^0\left(\sigma \right)\right){\text{e}}^{-{\displaystyle {\int }_{\sigma }^{\eta }\alpha \left(\xi \right)\text{d}\xi }}}\text{d}\sigma }{\text{e}}^{-{\displaystyle {\int }_{\eta }^a\gamma \left(\xi \right)\text{d}\xi }}\text{d}\eta \text{d}a\end{array}$ (3.14)

记：

$C=\underset{t\to +\infty }{\lim }\sup H\left(t\right)$

对(3.14)两端同时取上极限，根据Fatou引理可得

$C\le C{\displaystyle {\int }_0^{+\infty }h\left(a\right)\left(\epsilon f_{3,0}\left(a\right)+f_{4,0}\left(a\right)\right)\text{d}a}=C{R}_0$ .

因此，若$R_0<1$ ，

$C=\underset{t\to +\infty }{\lim }\sup H\left(t\right)=0$

根据$H\left(t\right)$ 的表达式可得

$\underset{t\to +\infty }{\lim }e\left(a,t\right)=0,\underset{t\to +\infty }{\lim }i\left(a,t\right)=0$

所以，当$R_0<1$ 时，无病平衡点$E^0$ 是全局渐近稳定的。

4. 地方平衡点的存在性和稳定性

由定理3.1可得，当$R_0>1$ 时，无病平衡点$E^0$ 不稳定，此时存在地方病平衡点。

定理4.1 若$R_0>1$ ，则系统(2.8)存在地方病平衡点$E^{*}=\left(s^{*}\left(a\right),v^{*}\left(a\right),e^{*}\left(a\right),i^{*}\left(a\right)\right)$ 。

证明：若系统(2.8)存在地方病平衡点$E^{*}=\left(s^{*}\left(a\right),v^{*}\left(a\right),e^{*}\left(a\right),i^{*}\left(a\right)\right)$ ，则满足：

$\{\begin{array}{l}\frac{\text{d}{s}^{*}}{\text{d}a}=-\left(\psi \left(a\right)+k\left(a\right)H^{*}\right)s^{*}\left(a\right)\\ \frac{\text{d}{v}^{*}}{\text{d}a}=\psi \left(a\right)s^{*}\left(a\right)-\delta k\left(a\right)H^{*}{v}^{*}\left(a\right)\\ \frac{\text{d}{e}^{*}}{\text{d}a}=\delta k\left(a\right)H^{*}{v}^{*}\left(a\right)+k\left(a\right)H^{*}{s}^{*}\left(a\right)-\alpha \left(a\right)e^{*}\left(a\right)\\ \frac{\text{d}{i}^{*}}{\text{d}a}=\alpha \left(a\right)e^{*}\left(a\right)-\gamma \left(a\right)i^{*}\left(a\right)\\ s^{*}\left(0\right)=1,v^{*}\left(0\right)=e^{*}\left(0\right)=i^{*}\left(0\right)=0\\ H^{*}={\displaystyle {\int }_0^{+\infty }h\left(a\right)\left(\epsilon e^{*}\left(a\right)+i^{*}\left(a\right)\right)\text{d}a}\end{array}$ (4.1)

求解(4.1)得

$\{\begin{array}{l}{s}^{*}\left(a\right)={\text{e}}^{-{\displaystyle {\int }_0^a\left(\psi \left(\xi \right)+k\left(\xi \right)H^{*}\right)\text{d}\xi }},\\ v^{*}\left(a\right)={\displaystyle {\int }_0^a\psi \left(\eta \right)s^{*}\left(\eta \right){\text{e}}^{-{\displaystyle {\int }_{\eta }^a\delta k\left(\xi \right)H^{*}\text{d}\xi }}\text{d}\eta },\\ e^{*}\left(a\right)=H^{*}{g}_{1,H^{*}}\left(a\right),\\ i^{*}\left(a\right)=H^{*}{g}_{2,H^{*}}\left(a\right)\end{array}$ (4.2)

$g_{1,H^{*}}\left(a\right)={\displaystyle {\int }_0^a\left(\delta k\left(\eta \right)v^{*}\left(\eta \right)+k\left(\eta \right)s^{*}\left(\eta \right)\right){\text{e}}^{-{\displaystyle {\int }_{\eta }^a\alpha \left(\xi \right)\text{d}\xi }}\text{d}\eta },$ (4.3)

$g_{2,H^{*}}\left(a\right)={\displaystyle {\int }_0^a\alpha \left(\eta \right)g_{1,H^{*}}\left(\eta \right){\text{e}}^{-{\displaystyle {\int }_{\eta }^a\gamma \left(\xi \right)\text{d}\xi }}\text{d}\eta }.$ (4.4)

将(4.2)的第三式和第四式代入到$H^{*}$ 的表达式中，且两边同时消去$H^{*}\left(H^{*}\ne 0\right)$ 得

$1={\displaystyle {\int }_0^{+\infty }h\left(a\right)\left(\epsilon g_{1,H^{*}}\left(a\right)+g_{2,H^{*}}\left(a\right)\right)\text{d}a}=G\left(H^{*}\right)$ (4.5)

又因为$G\left(0\right)=R_0$ ，且$G\left(H^{*}\right)$ 是关于$H^{*}$ 的单调递减函数，所以当$R_0>1$ 时，特征方程式(4.5)存在唯一正实根${\tilde{H}}^{*}$ ，即系统(2.8)存在唯一的地方病平衡点。

接下来讨论关于地方病平衡点的稳定性，并给出以下定理：

定理4.2 若$R_0>1$ ，则系统(2.8)的地方病平衡$E^{*}=\left(s^{*}\left(a\right),v^{*}\left(a\right),e^{*}\left(a\right),i^{*}\left(a\right)\right)$ 局部渐近稳定。

证明：考虑地方病平衡点$E^{*}$ 附近的扰动变量：

$\begin{array}{l}s\left(a,t\right)=s^{*}\left(a\right)+\overline{s}\left(a\right){\text{e}}^{\omega t},v\left(a,t\right)=v^{*}\left(a\right)+\overline{v}\left(a\right){\text{e}}^{\omega t}\\ e\left(a,t\right)=x^{*}\left(a\right)+\overline{x}\left(a\right){\text{e}}^{\omega t},i\left(a,t\right)=i^{*}\left(a\right)+\overline{i}\left(a\right){\text{e}}^{\omega t}\end{array}$

代入系统(2.8)中并保留其线性系统得：

$\{\begin{array}{l}\frac{\text{d}\overline{s}}{\text{d}a}=-\left(\omega +\psi \left(a\right)+k\left(a\right)H^{*}\right)\overline{s}\left(a\right)-k\left(a\right)H_2{s}^{*}\left(a\right)\\ \frac{\text{d}\overline{v}}{\text{d}a}=-\left(\omega +\delta k\left(a\right)H^{*}\right)\overline{v}\left(a\right)+\psi \left(a\right)\overline{s}\left(a\right)-\delta k\left(a\right)H_2{v}^{*}\left(a\right)\\ \frac{\text{d}\overline{x}}{\text{d}a}=-\left(\omega +\alpha \left(a\right)\right)\overline{x}\left(a\right)+\delta k\left(a\right)H_2{v}^{*}\left(a\right)+\delta k\left(a\right)H^{*}\overline{v}\left(a\right)\\ +k\left(a\right)H_2{s}^{*}\left(a\right)+k\left(a\right)H^{*}\overline{s}\left(a\right)\\ \frac{\text{d}\overline{i}}{\text{d}a}=-\left(\omega +\gamma \left(a\right)\right)\overline{i}\left(a\right)+\alpha \left(a\right)\overline{x}\left(a\right)\\ \overline{s}\left(0\right)=\overline{v}\left(0\right)=\overline{x}\left(0\right)=\overline{i}\left(0\right)=0\\ H_2={\displaystyle {\int }_0^{a_{+}}h\left(a\right)\left(\epsilon \overline{x}\left(a\right)+\overline{i}\left(a\right)\right)\text{d}a}\end{array}$ (4.6)

求解系统(4.6)可得

$\{\begin{array}{l}\overline{s}\left(a\right)=H_2{j}_{1,\omega }\left(a\right),\\ \overline{v}\left(a\right)=H_2{j}_{2,\omega }\left(a\right),\\ \overline{x}\left(a\right)=H_2{j}_{3,\omega }\left(a\right),\\ \overline{i}\left(a\right)=H_2{j}_{4,\omega }\left(a\right).\end{array}$

$\{\begin{array}{l}{j}_{1,w}\left(a\right)=-{\displaystyle {\int }_0^{a}k\left(\eta \right)s^{*}\left(\eta \right){\text{e}}^{-{\displaystyle {\int }_{\eta }^a\left(\omega +\psi \left(\xi \right)+k\left(\xi \right)H^{*}\right)\text{d}\xi }}\text{d}\eta }<0,\\ j_{2,w}\left(a\right)={\displaystyle {\int }_0^a\left(\psi \left(\eta \right)j_{1,w}\left(\eta \right)-\delta k\left(\eta \right)v^{*}\left(\eta \right)\right){\text{e}}^{-{\displaystyle {\int }_{\eta }^a\left(\omega +\delta k\left(\xi \right)H^{*}\right)\text{d}\xi }}\text{d}\eta }<0,\\ j_{3,w}\left(a\right)={\displaystyle {\int }_0^a\left(\delta k\left(\eta \right)v^{*}\left(\eta \right)+\delta k\left(\eta \right)H^{*}{j}_{2,w}\left(\eta \right)+k\left(\eta \right)s^{*}\left(\eta \right)+k\left(\eta \right)H^{*}{j}_{1,w}\left(\eta \right)\right){\text{e}}^{-{\displaystyle {\int }_{\eta }^a\left(\omega +\alpha \left(\xi \right)\right)\text{d}\xi }}\text{d}\eta },\\ j_{4,w}\left(a\right)={\displaystyle {\int }_0^a\alpha \left(\eta \right)j_{3,w}\left(\eta \right){\text{e}}^{-{\displaystyle {\int }_{\eta }^a\left(\omega +\gamma \left(\xi \right)\right)\text{d}\xi }}\text{d}\eta }.\end{array}$ (4.7)

将(4.7)的第三式和第四式代入$H_2$ 的表达式中，两边同时除以$H_2\left(H_2\ne 0\right)$ 可得

$1={\displaystyle {\int }_0^{+\infty }h\left(a\right)\left(\epsilon j_{3,\omega }\left(a\right)+j_{4,\omega }\left(a\right)\right)}:=L\left(\omega \right)$ (4.8)

注意到

$\begin{array}{c}L\left(0\right)={\displaystyle {\int }_0^{+\infty }h\left(a\right)\left(\epsilon j_{3,0}\left(a\right)+j_{4,0}\left(a\right)\right)\text{d}a}\\ =1+{\displaystyle {\int }_0^{+\infty }h\left(a\right)\epsilon }{\displaystyle {\int }_0^a\left(\delta k\left(\eta \right)H^{*}{j}_{2,0}\left(\eta \right)+k\left(\eta \right)H^{*}{j}_{1,0}\left(\eta \right)\right){\text{e}}^{-{\displaystyle {\int }_{\eta }^a\alpha \left(\xi \right)d\xi }}\text{d}\eta }\\ +{\displaystyle {\int }_0^a\alpha \left(\eta \right){\displaystyle {\int }_0^{\eta }\left(\left(\delta k\left(\sigma \right)H^{*}{j}_{2,0}\left(\sigma \right)+k\left(\sigma \right)H^{*}{j}_{1,0}\left(\sigma \right)\right){\text{e}}^{-{\displaystyle {\int }_{\sigma }^{\eta }\alpha \left(\xi \right)d\xi }}{\text{e}}^{-{\displaystyle {\int }_{\eta }^a\gamma \left(\xi \right)d\xi }}\text{d}\sigma \text{d}\eta \right)\text{d}a}}\\ <1.\end{array}$

此外，

$\begin{array}{l}L\left(\omega \right)<{\displaystyle {\int }_0^{+\infty }h\left(a\right)\epsilon }{\displaystyle {\int }_0^a\left(\delta k\left(\eta \right)v^{*}\left(\eta \right)+k\left(\eta \right)s^{*}\left(\eta \right)\right){\text{e}}^{-{\displaystyle {\int }_{\eta }^a\left(\omega +\alpha \left(\xi \right)\right)d\xi }}}\text{d}\eta \\ +{\displaystyle {\int }_0^a\alpha \left(\eta \right)}\left(\left({\displaystyle {\int }_0^{\eta }\left(\delta k\left(\sigma \right)v^{*}\left(\sigma \right)+k\left(\sigma \right)s^{*}\left(\sigma \right)\right){\text{e}}^{-{\displaystyle {\int }_{\sigma }^{\eta }\left(\omega +\alpha \left(\xi \right)\right)\text{d}\xi }}}\right)\text{d}\sigma {\text{e}}^{-{\displaystyle {\int }_{\eta }^a\left(\omega +\gamma \left(\xi \right)\right)d\xi }}\text{d}\eta \right)\text{d}a:=\widehat{L}\left(\omega \right)\end{array}$

易得$\widehat{L}\left(\omega \right)$ 是关于$\mathrm{Re}\omega$ 的单调递减函数，故$\mathrm{Re}\omega >0$ 时，有

$L\left(\omega \right)<\widehat{L}\left(0\right)=G\left(H^{*}\right)=1$

因此特征方程(4.8)当且仅当$\mathrm{Re}\omega <0$ 时成立，即地方病平衡点$E^{*}$ 局部渐近稳定。

5. 结论

本文我们建立了年龄结构模型的基本再生数$R_0$ 的显式表达式。通过观察$R_0$ 的基本表达式可得提高疫苗接种成功率或者减小潜伏期相对感染期比值可以有效降低$R_0$ ，从而使疾病得到较好的控制。此外，本文根据$R_0$ 的显式表达式证明了无病平衡点的局部稳定性，并借助特征线法证明了无病平衡点的全局稳定性；并研究了地方病平衡点的存在性和局部渐近稳定性，但是地方病平衡点的全局稳定性还无法证明。研究地方病平衡点的全局稳定性还是一项具有挑战性的一项任务，在今后的工作中仍需继续研究解决这个问题。