传染病动力学是对传染病进行理论性定量研究的一种重要方法^[1].在最近20年来,国内外关于传染病动力学的研究进展非常迅速^[2-3],各类微分方程模型被广泛用于研究各种各样的传染病问题,并为疾病的预防和控制提供了许多科学依据和指导.周期微分方程传染病模型由于考虑了季节更替、气候规律性变化等因素对疾病传播的影响而显得更加符合实际,在近些年随着相关理论工具的不断发展完善^[4],逐渐成为传染病动力学研究的最重要关注点之一.

在传染病动力学中,基本再生数是一个很重要的概念,其表示单个感染者在完全易感染人群中生成的第二代期望数,被用来衡量疾病的传染能力大小,可以作为判断疾病在种群中流行与消亡的阈值^[5]. Liang和Zhao等^[6-7]中引入了带有时滞的周期传染病模型基本再生数的定义方法,并证明了基本再生数与染病者变量在周期无病平衡点处线性化方程对应Poincáre映射谱半径在刻画疾病传播的"阈值"作用方面的等价性,是针对周期性传染病模型(包括具有阶段特征的种群演化周期模型)动力学研究的突破性成果.自此之后, Bai^[8]和Lou^[9]等以此为理论基础,分别研究了一类周期疟疾模型和一个寄生虫的演化模型,更多的相关研究可以参阅文献[11-12].特别指出,作为应用, Zhao等在文献[7]中研究了一个带有离散时滞的SEIR常微分方程传染病模型,定义了模型的基本再生数,并证明了基本再生数在决定疾病是否一致持久的"阈值"作用.

在离散时滞传染病模型中,时滞表示疾病的潜伏期,该潜伏期实际上是一个平均值.在现实中,潜伏期的长短受诸多随机因素的影响,特别对温度异常敏感,而温度又依赖于发病期所处的时间,因而将潜伏期视做一个依赖于时间的周期函数将更加符合实际^{[9, 13-14]}.目前国内外针对具有周期性潜伏传染病数学模型的研究成果非常少.本文将以文献[7]中的模型为基础,考虑引入周期性的潜伏期,克服由此带来的技术性困难,并利用文献[6-7]中的成果研究模型的阈值动力学.

文献[7]所考虑的时间周期时滞常微分方程模型如下:

(1.1) $ \label{i}\left\{\begin{array}{ll}\frac{{\rm d}S(t)}{{\rm d}t}=\Lambda(t)-f(t, S(t), I(t))-\mu(t)S(t)+\alpha(t)R(t), \\\frac{{\rm d}E(t)}{{\rm d}t}=f(t, S(t), I(t))-e^{-\int_{t-\tau}^{t}\mu(r){\rm d}r}f(t-\tau, S(t-\tau), I(t-\tau))-\mu(t)E(t), \\\frac{{\rm d}I(t)}{{\rm d}t}=e^{\int_{t-\tau}^{t}\mu(r){\rm d}r}f(t-\tau, S(t-\tau), I(t-\tau))-(\mu(t)+d(t)+\gamma(t))I(t), \\\frac{{\rm d}R(t)}{{\rm d}t}=\gamma(t)I(t)-\mu(t)R(t)-\alpha(t)R(t), \end{array}\right. $

其中, $S(t), E(t), I(t), R(t)$分别表示$t$时刻易感者、潜伏期内的感染者(已经染病却暂时不具备传染能力)、感染者、具有免疫力者的总量. $\mu(t), \Lambda(t)$分别表示$t$时刻种群的自然死亡率和总补充率. $d(t), \gamma(t), \alpha(t)$分别为$t$时刻的因病致死率、恢复率、和免疫丧失率.离散时滞$\tau>0$为平均的潜伏期. $ f$表示疾病的发生率. $ f(t, \cdot, \cdot)$, $\Lambda(t)$, $\mu(t)$, $d(t)$, $\alpha(t)$, $\gamma(t)$皆为$\omega$ -周期函数, $\omega>0$为常数.文献[7]对模型(1.1)中的参数做出如下假设:

(H1) $\Lambda(t)$, $\mu(t)$, $d(t)$, $\alpha(t)$, $\gamma(t)$是非负连续函数,且

$ \Lambda(t)>0, \qquad\int_0^\omega\mu(t){\rm d} t>0, \qquad \int_0^\omega\gamma(t){\rm d} t>0. $

(H2) $f(\cdot, S, I)$是$C^1$的,而且满足

(ⅰ) $f(\cdot, S, 0)\equiv0$, $f(\cdot, 0, I)\equiv0$;对所有的$S>0$, $\frac{\partial f(\cdot, S, 0)}{\partial I}>0$且关于$S$单调不减;

(ⅱ)对所有的$S, I>0$,都有$\frac{\partial f(\cdot, S, 0)}{\partial S}\geq 0$和$f(\cdot, S, I)\leq\frac{\partial f(\cdot, S, 0)}{\partial I}I$成立.

为了得到带有周期性潜伏期的SEIR模型,本文假定潜伏期是一个与时间有关的严格为正的周期函数$\tau(t)$.新模型的建模中不能通过简单地用$\tau(t)$替换系统(1.1)中的常数$\tau$来完成.此外,不失一般性,本文对于疾病的发生率$f$,采用满足(H2)的最典型的发生率—饱和发生率$f(t, S, I)=\frac{\beta(t)SI}{1+c(t)I}$,这里$c(t)$是非负连续$\omega$ -周期函数,且$c(t)\not\equiv 0, t\in[0, \omega]$.

引入染病年龄分布函数$\rho(t, a)$,这里$\rho(t, a)$表示$t$时刻染病年龄为$a$的染病者的密度.由马知恩等^[1]知$\rho(t, a)$满足

$ \frac{\partial\rho(t, a) }{\partial t}+\frac{\partial \rho(t, a) }{\partial a}=-\mu(t, a) \rho(t, a), $

其中$\mu(t, a)=\left\{\begin{array}{ll} \mu(t), ~~&\mbox{当}~a\leq\tau(a) ~\mbox{时}; \\ \mu(t)+d(t)+\gamma(t), ~~&\mbox{当}~a\geq\tau(a) ~\mbox{时}. \end{array}\right.$

令

$ E(t)=\int_{0}^{\tau(t)}\rho(t, a){\rm d}a, ~~~ I(t)=\int_{\tau(t)}^{\infty}\rho(t, a){\rm d}a. $

对$E(t)$和$I(t)$分别求导数,得到

$ \begin{eqnarray} \frac{{\rm d} E(t) }{{\rm d}t}&=& \int_{0}^{\tau(t)} \frac{\partial \rho(t, a)}{\partial t}{\rm d}a+\tau'(t)\rho(t, \tau(t))\\&=&\int_{0}^{\tau(t)} \big(-\frac{\partial \rho(t, a)}{\partial a}-\mu(t) \rho(t, a)\big){\rm d}a+\tau'(t)\rho(t, \tau(t))\\&=&-(1-\tau'(t))\rho(t, \tau(t))+\rho(t, 0)-\mu(t)E(t) \end{eqnarray} $

和

$ \begin{eqnarray} \frac{{\rm d} I(t) }{{\rm d} t}&= &\int_{\tau(t)}^{\infty} \frac{\partial \rho(t, a)}{{\rm d} t}{\rm d}a-\tau'(t)\rho(t, \tau(t))\\&=&\int_{\tau(t)}^{\infty} \big(-\frac{{\rm d} \rho(t, a)}{{\rm d} a} -(\mu(t)+d(t)+\gamma(t)) \rho(t, a)\big){\rm d}a-\tau'(t)\rho(t, \tau(t))\\&=&-\rho(t, \infty)+(1-\tau'(t))\rho(t, \tau(t))-(\mu(t)+d(t)+\gamma(t)) I(t). \end{eqnarray} $

注  根据生态学、传染病学意义^[9], $\tau(t)$满足$1-\tau'(t)>0$, $\forall t\in[0, \omega]$.

由于任何个体的生命都是有限的,因此$\rho(t, \infty)=0$,同时易见$\rho(t, 0)=\frac{\beta(t)S(t)I(t)}{1+c(t)I(t)}$.下面采用沿着特征线积分的办法计算$\rho(t, \tau(t))$.令$V^s(t)=\rho(t, t-s)$,则有

$ \left\{\begin{array}{ll}\frac{{\rm d} V^s(t)}{{\rm d} t}=-\mu(t)V^s(t), \\[3mm] V^s(s)=\rho(s, 0)=\frac{\beta(s)S(s)I(s)}{1+c(s)I(s)}.\end{array}\right. $

求解以上方程,得到

$ V^s(t)=e^{-\int_s^t\mu(r){\rm d}r}V^s(s)=e^{-\int_s^t\mu(r){\rm d}r}\frac{\beta(s)S(s)I(s)}{1+c(s)I(s)}, $

则

$ \rho(t, \tau(t))=V^{t-\tau(t)}(t)=e^{-\int_{t-\tau(t)}^t\mu(r){\rm d}r}\frac{\beta(t-\tau(t))S(t-\tau(t))I(t-\tau(t))}{1+c(t-\tau(t))I(t-\tau(t))}. $

从而得到如下带有周期性潜伏期的SEIR传染病模型

(1.2) $ \label{I4}\left\{\begin{array}{ll}\frac{{\rm d}S(t)}{{\rm d} t}=\Lambda(t)-\frac{\beta(t)S(t)I(t)}{1+c(t)I(t)}-\mu(t)S(t)+\alpha(t)R(t), \\[3mm]\frac{{\rm d}E(t)}{{\rm d} t}=\frac{\beta(t)S(t)I(t)}{1+c(t)I(t)}-(1-\tau'(t))e^{-\int_{t-\tau(t)}^t\mu(r){\rm d}r}\frac{\beta(t-\tau(t))S(t-\tau(t))I(t-\tau(t))}{1+c(t-\tau(t))I(t-\tau(t))}\\[1mm]\ \qquad\qquad-\mu(t)E(t), \\[2mm]\frac{{\rm d}I(t)}{{\rm d} t}=(1-\tau'(t))e^{-\int_{t-\tau(t)}^t\mu(r){\rm d}r}\frac{\beta(t-\tau(t))S(t-\tau(t))I(t-\tau(t))}{1+c(t-\tau(t))I(t-\tau(t))}+\gamma(t))I(t), \\[3mm]\frac{{\rm d}R(t)}{{\rm d} t}=\gamma(t)I(t)-\mu(t)R(t)-\alpha(t)R(t).\end{array}\right. $

注意到$E(t)$并没有出现在系统(1.2)其余的方程中,对其进行常数变异法,可将其化成积分方程的形式,具体过程如下:

$ \begin{eqnarray} E(t)&= &e^{ -\int_0^t\mu(r){\rm d}r }\bigg\{E(0)+\int_0^t\bigg[ \frac{\beta(s)S(s)I(s)}{1+c(s)I(s)}-(1-\tau'(s)) \\&&\times e^{-\int_{s-\tau(s)}^s\mu(r){\rm d}r}\frac{\beta(s-\tau(s))S(s-\tau(s))I(s-\tau(s))}{1+c(s-\tau(s))I(t-\tau(s))} \bigg ] e^{\int_0^s\mu(r){\rm d}r} {\rm d}s \bigg\}\\&=&e^{ -\int_0^t\mu(r){\rm d}r }\bigg(E(0)+\int_0^te^{\int_0^s\mu(r){\rm d}r} \frac{\beta(s)S(s)I(s)}{1+c(s)I(s)} {\rm d}s \\&&-\int_0^t(1-\tau'(s))e^{-\int_{s-\tau(s)}^0\mu(r){\rm d}r}\frac{\beta(s-\tau(s))S(s-\tau(s))I(s-\tau(s))}{1+c(s-\tau(s))I(t-\tau(s))} {\rm d}s \bigg)\\&=&e^{ -\int_0^t\mu(r){\rm d}r }\bigg(E(0)-\int_{-\tau(0)}^{0}e^{ \int_0^s\mu(r){\rm d}r } \frac{\beta(s)S(s)I(s)}{1+c(s)I(s)}{\rm d}s\\&&+\int_{t-\tau(t)}^t e^{\int_0^s\mu(r){\rm d}r} \frac{\beta(s)S(s)I(s)}{1+c(s)I(s)} {\rm d}s \bigg). \end{eqnarray} $

为了使系统(1.2)的解是连续的^[10],本文始终假定系统(1.2)中初值条件满足

(1.3) $ \label{E0}E(0)=\int_{-\tau(0)}^{0}e^{ \int_0^s\mu(r){\rm d}r } \frac{\beta(s)S(s)I(s)}{1+c(s)I(s)}{\rm d}s. $

由初值条件(1.3)知

(1.4) $ \label{Eint}E(t)= \int_{t-\tau(t)}^t e^{-\int_s^t\mu(r){\rm d}r} \frac{\beta(s)S(s)I(s)}{1+c(s)I(s)} {\rm d}s, $

这意味着系统(1.2)中$E(t)$可完全由$S(t)$和$I(t)$得到.

为了简化问题,将第二个方程从系统(1.2)中去掉,考虑如下系统即可

(1.5) $ \label{I}\left\{\begin{array}{ll}\frac{{\rm d}S(t)}{{\rm d} t}=\Lambda(t)-\frac{\beta(t)S(t)I(t)}{1+c(t)I(t)}-\mu(t)S(t)+\alpha(t)R(t), \\[3mm]\frac{{\rm d}I(t)}{{\rm d} t}=(1-\tau'(t))e^{-\int_{t-\tau(t)}^t\mu(r){\rm d}r}\frac{\beta(t-\tau(t))S(t-\tau(t))I(t-\tau(t))}{1+c(t-\tau(t))I(t-\tau(t))}\\[1mm]\qquad-(\mu(t)+d(t)+\gamma(t))I(t), \\[2mm]\frac{{\rm d}R(t)}{{\rm d} t}=\gamma(t)I(t)-\mu(t)R(t)-\alpha(t)R(t).\end{array}\right. $

本文将借助于耗散动力系统的相关理论结合时滞周期系统基本再生数理论研究系统(1.5) (系统(1.2))的动力学.周期性潜伏期$\tau(t)$的引入会使模型变得更加符合实际,同时会对系统动力学行为的分析过程带来新的困难,比如解空间的选取等.本文剩余部分组织如下:首先证明系统解的存在唯一性,其次定义系统的基本再生数,再证明基本再生数在决定两种不同全局动力学形态时的"阈值"作用,最后用数值方法验证所得结论,并分析结论的实际应用价值.

令$\hat{\tau}=\max\limits_{t\in[0, \omega]}\tau(t)$.再令$W:=C([-\hat{\tau}, 0], {\mathbb R}^4)$, $X:=C([-\hat{\tau}, 0], {\mathbb R}^3)$, $Y:=C([-\tau(0), 0], {\mathbb R}^2)\times C({\mathbb R})$, $W^+:=C([-\hat{\tau}, 0], {\mathbb R}^4_+)$, $X^+:=C([-\hat{\tau}, 0], {\mathbb R}^3_+)$, $Y^+:=C([-\tau(0), 0], {\mathbb R}^2_+)\times C({\mathbb R})_+$.则$(W, W^+)$, $(X, X^+)$和$(Y, Y^+)$都是序Banach空间.

显然,对于任意满足(1.3)的$\psi\in W^+$,都存在$t_\psi <\infty$,使得系统(1.2)在$[0, t_\psi)$存在唯一的解$v(t, \psi)=(S(t), E(t), I(t), R(t))$,该解满足初值条件$v_0=\psi$.令$N(t)=S(t)+E(t)+I(t)+R(t)$,则由系统(1.2)得到

$ \frac{{\rm d}N(t)}{{\rm d} t}=\Lambda(t)-\mu(t)N(t)-d(t)I(t)\leq\Lambda(t)-\mu(t)N(t). $

由于

(2.1) $ \label{d}\frac{{\rm d}w(t)}{{\rm d} t}=\Lambda(t)-\mu(t)w(t) $

有唯一的全局吸引的周期解

$ S^*(t)=\bigg(\int_0^t\Lambda(r)e^{\int_0^r\mu(s){\rm d}s} {\rm d}r+\frac{\int_0^\omega\Lambda(r)e^{\int_0^r\mu(s){\rm d}s} {\rm d}r}{e^{\int_0^\omega\mu(s){\rm d}s}-1}\bigg)e^{-\int_0^t\mu(s){\rm d}s}. $

由比较原理知对于任意满足条件(1.3)的$\psi\in W^+$,满足初值条件$w_0=\psi$的解$S(t), E(t), I(t), $$R(t)$在$[0, t_\psi)$上有界,从而由解的延拓性定理知, $t_\psi=\infty$,这意味着对于任意满足条件(1.3)的$\psi\in W^+$,系统(1.2)满足初值条件$w_0=\psi$的解是全局存在的.从而由系统(1.5)和系统(1.2)的关系得到如下结果:

引理2.1  对于任意$\psi\in X^+$,系统(1.5)在$[0, +\infty)$上存在满足$w_0=\psi$的唯一解.

引理2.2   对于任意$\phi\in Y^+$,系统(1.5)在$[0, +\infty)$上存在满足$u_0=\phi$的唯一解.

证  令$\bar{\tau}=\min\limits_{t\in[0, \omega]}\tau(t)$.对于任意$t\in[0, \bar{\tau}]$,由于$t-\tau(t)$是严格单调递增的,有

$ -\tau(0)=0-\tau(0)\leq t-\tau(t)\leq \bar{\tau}-\tau(\bar{\tau})\leq \bar{\tau}-\bar{\tau}=0, $

即$-\tau(0)\leq t-\tau(t)\leq0$,因此有

$ S(t-\tau(t))=\phi_1(t-\tau(t)), ~I(t-\tau(t))=\phi_2(t-\tau(t)). $

从而得到以下微分方程

$ \left\{\begin{array}{ll}\frac{{\rm d}S(t)}{{\rm d} t}=\Lambda(t)-\frac{\beta(t)S(t)I(t)}{1+c(t)I(t)}-\mu(t)S(t)+\alpha(t)R(t), \\\frac{{\rm d}I(t)}{{\rm d} t}=(1-\tau'(t))e^{-\int_{t-\tau(t)}^t\mu(r){\rm d}r} \frac{\beta(t-\tau(t))\phi_1(t-\tau(t))\phi_2(t-\tau(t))}{1+c(t-\tau(t))\phi_2(t-\tau(t))}\\[1mm]\qquad -(\mu(t)+d(t)+\gamma(t))I(t), \\[2mm]\frac{{\rm d}R(t)}{{\rm d} t}=\gamma(t)I(t)-\mu(t)R(t)-\alpha(t)R(t), ~~t\in[0, \bar{\tau}], \\[2mm]S(0)=\phi_1(0), I(0)=\phi_2(0), R(0)=\phi_3(0).\end{array}\right. $

对于给定$\phi\in Y^+$,以上系统的解$(S(t), I(t), R(t))$在$t\in[0, \bar{\tau}]$上是存在唯一的.换言之, $\psi_1(\theta):=S(\theta)$, $\psi_2(\theta):=I(\theta)$, $\theta\in[-\tau(0), \bar{\tau}]$是已知的, $\psi_3(\theta):=R(\theta)$, $\theta\in[0, \bar{\tau}]$是已知的.则对于$t\in[\bar{\tau}, 2\bar{\tau}]$,有

$ -\tau(0)=0-\tau(0)\leq \bar{\tau}-\tau(\bar{\tau})\leq t-\tau(t)\leq 2\bar{\tau}-2\bar{\tau}\leq2\bar{\tau}-\bar{\tau}=\bar{\tau}, $

即$-\tau(0)\leq t-\tau(t)\leq \bar{\tau}$,此时$S(t-\tau(t))=\psi_1(t-\tau(t))$, $I(t-\tau(t))=\psi_2(t-\tau(t))$.求解以下微分方程

$ \left\{\begin{array}{ll}\frac{{\rm d}S(t)}{{\rm d} t}=\Lambda(t)-\frac{\beta(t)S(t)I(t)}{1+c(t)I(t)}-\mu(t)S(t)+\alpha(t)R(t), \\[3mm]\frac{{\rm d}I(t)}{{\rm d} t}=(1-\tau'(t))e^{-\int_{t-\tau(t)}^t\mu(r){\rm d}r}\frac{\beta(t-\tau(t))\psi_1(t-\tau(t))\psi_2(t-\tau(t))}{1+c(t-\tau(t))\psi_2(t-\tau(t))}\\\qquad-(\mu(t)+d(t)+\gamma(t))I(t), \\[2mm]\frac{{\rm d}R(t)}{{\rm d} t}=\gamma(t)I(t)-\mu(t)R(t)-\alpha(t)R(t), ~~t\in[\bar{\tau}, 2\bar{\tau}], \\[2mm]S(\bar{\tau})=\psi_1(\bar{\tau}), I(\bar{\tau})=\psi_2(\bar{\tau}), R(\bar{\tau})=\psi_3(\bar{\tau}).\end{array}\right. $

从而得到以上方程的解$(S(t), I(t), R(t))$在$t\in[\bar{\tau}, 2\bar{\tau}]$上是存在唯一的.接下来,在$t\in[2\bar{\tau}, 3\bar{\tau}], [3\bar{\tau}, 4\bar{\tau}]$, $\cdots $, $[n\bar{\tau}, (n+1)\bar{\tau}], \cdots $区间上重复以上过程,可以得到得到系统(1.5)满足初值条件$u_0=\phi$的解$u(t, \phi)$对$t\in[0, +\infty)$是存在唯一的.

将系统(1.5)在$(0, S^*, 0)$处线性化,得到关于$I$的方程

(2.2) $ \label{lin}\frac{{\rm d}I(t)}{{\rm d} t}=p(t)I(t-\tau(t))-q(t)I(t), $

这里

$ p(t):=(1-\tau'(t))e^{-\int_{t-\tau(t)}^t\mu(r){\rm d}r}S^*(t-\tau(t)), ~~q(t):=\mu(t)+d(t)+\gamma(t). $

令$F(t)\phi:=p(t)\phi(-\tau(t)), V(t):=q(t)$, $\Phi(t, s):=e^{-\int_s^t q(r){\rm d}r}, \forall t\geq s, s, t\in{\mathbb R}$.用$C_\omega$表示由${\mathbb R}$映向${\mathbb R}$的$\omega$ -周期函数所组成的连续函数空间,其范数为最大模范数.假定$v\in C_\omega$为染病群体的初始分布.对于任意$s\geq0$, $F(t-s)v_{t-s}$表示$t-s$时刻的新染病者的分布,这是由在$[t-s-\hat{\tau}, t-s]$时段内的染病者所传染的.从而$\Phi(t, t-s)F(t-s)v_{t-s}$代表$t-s$时刻新产生的感染者并在$t$时刻仍然染病的个体的分布.这样,得到$t$时刻所有染病个体(由$t$时刻以前的所有可能时刻感染造成)的分布为

$ \int_0^\infty \Phi(t, t-s)F(t-s)v_{t-s}{\rm d}s. $

定义次代算子$L:C_\omega\longrightarrow C_\omega$:

$ \begin{eqnarray} L[v](t)&=&\int_0^\infty \Phi(t, t-s)F(t-s)v_{t-s}{\rm d}s \\&=&\int_0^\infty \Phi(t, t-s)F(t-s)v(t-s+\cdot){\rm d}s \\&=&\int_0^\infty \Phi(t, t-s)F(t-s)v(t-s-\tau(t-s)){\rm d}s. \end{eqnarray} $

令$t-s-\tau(t-s)=t-s_1$.考虑到$y=x-\tau(x)$是一个严格递增函数,因此其反函数是存在的,记为$x=h(y)$.因此, $t-s=h(t-s_1)$,即$s=t-h(t-s_1)$,从而${\rm d}s_1=(1-\tau'(t-s)){\rm d}s$,则${\rm d}s=\frac{1}{1-\tau'(h(t-s_1))}{\rm d}s_1$.由此得到

$ \begin{eqnarray} L[v](t)&=&\int_0^\infty \Phi(t, t-s)F(t-s)v(t-s-\tau(t-s)){\rm d}s \\&=&\int_{\tau(t)}^\infty \frac{\Phi(t, h(t-s_1))F(h(t-s_1))v(t-s_1)}{1-\tau'(h(t-s))}{\rm d}s\\&=&\int_{\tau(t)}^\infty e^{-\int_{h(t-s_1)}^t q(r){\rm d}r} e^{-\int_{h(t-s_1)-\tau(h(t-s_1))}^{h(t-s_1)}\mu(r){\rm d}r}S^*(h(t-s_1)-\tau(h(t-s_1))){\rm d}s_1\\&=&\int_{\tau(t)}^\infty e^{-\int_{h(t-s)}^t q(r){\rm d}r} e^{-\int_{h(t-s)-\tau(h(t-s))}^{h(t-s)}\mu(r){\rm d}r}S^*(h(t-s)-\tau(h(t-s)))v(t-s){\rm d}s. \end{eqnarray} $

紧接着,定义算子$L$的谱半径${\mathcal R}_0$做为模型(1.5)的基本再生数

$ {\mathcal R}_0:=r(L). $

对任意$t\geq0$,用$\hat{P}(t)$表示线性化系统(2.2)的解映射,即$\hat{P}_t(\phi)=\hat{I}(t, \phi), \forall t\geq0$,这里$\hat{I}(t, \phi)$表示方程(2.2)满足初值条件$\hat{I}_0=\phi\in C([-\hat{\tau}, 0], {\mathbb R})$的解.则$\hat{P}:=\hat{P}(\omega)$是对应于线性化系统(2.2)的Poincáre映射.令$r(\hat{P})$表示$\hat{P}$的谱半径.

对任意$t\geq0$,用$P(t)$表示线性化系统(2.2)的解映射,即$P_t(\phi)=I(t, \phi), \forall t\geq0$,这里$I(t, \phi)$为方程(2.2)满足初值条件$I_0=\phi\in C([-\tau(0), 0], {\mathbb R})$的解.则$P:=\hat{P}(\omega)$是对应于线性化系统(2.2)的Poincáre映射.令$r(P)$表示$P$的谱半径.和文献[9,引理3.8]的证明一样,可以验证$r(\hat{P})-1$与$r(P)-1$的符号相同.而根据文献[7,定理2.1]知, ${\mathcal R}_0-1$与$r(\hat{P})-1$符号相同.从而有如下结论：

引理2.3   ${\mathcal R}_0-1$与$r(P)-1$的符号相同.

引理3.1   令$\lambda=\frac{\ln r(P)}{\omega}, $则存在一个正$\omega$ -周期函数$\eta(t)$使得$e^{\lambda t}\eta(t)$是系统(2.2)的一个解.

证  令$\xi(t)$为Poincáre映射$P$对应于$r(P)$的特征函数,则$P(\xi)=r(P)\xi$且$\xi(t)\gg0$.用$w(t, \xi)$表示系统(2.2)初值为$\xi$的解,由$\xi(t)\gg0$知$w(t, \xi)\gg 0$.令$\eta(t)=e^{-\lambda t}w(t, \xi)$,则有$\eta(t)>0$, $\forall t\geq-\tau(t)$,并且

$ \frac{{\rm d} \eta(t)}{{\rm d} t}=e^{-\lambda t} \frac{{\rm d}w(t, \xi)}{{\rm d} t} -\lambda e^{-\lambda t}w(t, \xi), $

即

(3.1) $ \label{lin-mu}\frac{{\rm d}\eta(t)}{{\rm d} t}=p(t)e^{-\lambda\tau(t)}\eta(t-\tau(t))-(q(t)+\lambda)\eta(t). $

从而$\eta(t)$是方程(3.1)的一个解.由$p(t), \tau(t), q(t)$的周期性知方程(3.1)亦是一个$\omega$ -周期方程,这意味着$\eta(t+\omega)$也是方程(3.1)的一个解.

由于$\lambda=\frac{\ln r(P)}{\omega}$知$e^{-\lambda\omega}r(P)=1$,对任意$\theta\in[-\tau(0), 0]$,有

$ \begin{eqnarray} \eta( \omega+\theta)&=e^{-\lambda( \omega +\theta)}P (\xi)(\theta)=e^{-\lambda(\omega +\theta)}r(P) \xi(\theta)=e^{-\lambda \theta}w(\theta, \xi)=\eta(\theta). \end{eqnarray} $

因此, $\eta(\omega+t)=\eta(t), \forall t\in[-\tau(0), 0)$.从而由方程(3.1)解的唯一性知, $\eta(t)=\eta(t+\omega)$, $\forall t\geq 0$,即$\eta(t)$是方程(3.1)的一个$\omega$ -周期解.显然, $e^{\lambda t}\eta(t)$是方程(2.2)的一个解.

令$Y_0:=\{\phi=(\phi_1, \phi_2, \phi_3)\in Y^+: \phi_2(0)>0 \}$,现给出本文主要结果如下:

定理3.1   假定H1成立,有如下结论成立：

(ⅰ)如果${\mathcal R}_0 <1$,则系统(1.5)的无病平衡$\omega$ -周期解$(S^*(t), 0, 0)$在$Y$上是全局渐近吸引的.

(ⅱ)如果${\mathcal R}_0>1$,则系统(1.5)至少有一个正$\omega$ -周期解$(\tilde{S}(t), \tilde{I}(t), \tilde{R}(t))$,同时存在一个$\delta>0$使得对任意$\phi\in Y_0$,有

$ \liminf\limits_{t\rightarrow\infty}(S(t, \phi), I(t, \phi), R(t, \phi))>(\delta, \delta, \delta). $

证   (ⅰ)当${\mathcal R}_0 <1$时,根据引理2.3可得$r(P) <1$.对任意$\epsilon>0$,考虑方程(2.2)的带有参数$\epsilon >0$扰动的方程

(3.2) $ \label{lin-e}\frac{{\rm d}I(t)}{{\rm d} t}=p_\epsilon(t)I(t-\tau(t))-q(t)I(t), $

这里$p_\epsilon(t):=(1-\tau'(t))e^{-\int_{t-\tau(t)}^t\mu(r){\rm d}r}(S^*(t-\tau(t))+\epsilon)$.用$r(P_\epsilon)$表示对应于线性化方程(3.2)的Poincáre映射.令$r(P_\epsilon)$表示$P_\epsilon$的谱半径,由于$\lim\limits_{\epsilon\rightarrow0}r(P_\epsilon)=r(P)$.因此,可以固定一个充分小的$\epsilon_0>0$,使得$r(P_{\epsilon_0}) <1$.依引理3.1知,存在一个正$\omega$ -周期函数$\eta_{\epsilon_0}(t)$使得$e^{\lambda_{\epsilon_0} t}\eta_{\epsilon_0}(t)$是方程(3.2)的一个解(方程(3.2)中取$\epsilon=\epsilon_0$),特别地, $\eta_{\epsilon_0}(t)>0$, $\forall t>0$,这里$\lambda_{\epsilon_0}=\frac{\ln r(P_{\epsilon_0})}{\omega} <0$.

另一方面,由于方程(2.1)有唯一的全局吸引的周期解$S^*(t)$,则存在充分大的正整数$N$,使得$N\omega>\hat{\tau}$,且$S(t)\leq S^*(t)+{\epsilon_0}$, $\forall t\geq N\omega-\hat{\tau}$.由此可以得到

$ \frac{{\rm d}I(t)}{{\rm d} t}\leq p_{\epsilon_0}(t)I(t-\tau(t))-q(t)I(t), ~~~t>N\omega, $

这里$p_{\epsilon_0}(t):=(1-\tau'(t))e^{-\int_{t-\tau(t)}^t\mu(r){\rm d}r}(S^*(t-\tau(t))+{\epsilon_0})$.

存在充分大的常数$B$,使得$I(t)\leq Be^{\lambda_{\epsilon_0} t}\eta_{\epsilon_0}(t)$, $\forall t\in[N\omega-\hat{\tau}, N\omega]$.利用时滞微分方程的比较原理可以得到$I(t)\leq Be^{\lambda_{\epsilon_0} t}\eta_{\epsilon_0}(t)$, $\forall t>N\omega$.由于$\lambda_{\epsilon_0} <0$,则有$\lim\limits_{t\rightarrow\infty}I(t)=0$.最后,利用链传递集理论^[4]可以验证$\lim\limits_{t\rightarrow\infty}R(t)=0$和$\lim\limits_{t\rightarrow\infty}(S(t)-S^*(t))=0$成立.

(ⅱ)当${\mathcal R}_0>1$时,根据引理2.3可得$r(P)>1$.令

$ Y_0=\{\phi\in Y^+:\phi_2(0)\not\equiv0\}, $

$ \partial Y _0=Y^+\setminus Y _0=\{\phi\in Y^+:\phi_2(0)\not\equiv0\}, $

对于任意$t\geq0$,定义$Q_t\phi=u(t, \phi)$, $\forall \phi\in X$,这里$u(t, \phi)$系统(1.5)满足初值条件$u_0=\phi$的解.回顾周期半流的定义,容易验证, $\{Q_t\}_{t\geq0}$是$X^+$上的一个$\omega$ -周期半流.同时,由系统(1.5)解的一致有界性知, $Q_t:~X^+\rightarrow X^+$是点耗散的.令$Q:=Q_\omega$.考虑如下扰动线性方程

$ \frac{{\rm d}I(t)}{{\rm d} t}= p_{\varepsilon}(t)I(t-\tau(t))-q(t)I(t), $

这里$p_{\varepsilon}(t):=(1-\tau'(t))e^{-\int_{t-\tau(t)}^t\mu(r){\rm d}r} \big(S^*(t-\tau(t))-\varepsilon \big)$.用$r(P_\varepsilon)$表示以上方程的Poincáre映射.由于$r(P)>1$且$\lim\limits_{\varepsilon\rightarrow 0}r(P_{\varepsilon})=r(P)$.因此可以固定一个充分小的$\varepsilon_0>0$,使得$r(P_{\varepsilon_0})>1$.同时,存在一个充分小的$\sigma>0$,使得

(3.3) $ \label{*}\frac{S(t-\tau(t))-\sigma}{[1+c(t-\tau(t))I(t-\tau(t))]^2} \geq S^*(t-\tau(t))-\varepsilon , ~~~\forall I\in[0, \sigma]. $

令$E:=\{(S^*(t), 0, 0)\}$.为了证明主要结果,先验证如下论断.

论断    $E$是$ Z _0$的一个弱排斥子,即

$ \limsup\limits_{n\rightarrow\infty}\left\| Q^n(\phi)- E\right\| \geq\sigma, ~~~ \forall \phi\in Y_0. $

采用反证法,假设存在$\tilde{\phi}\in Y _0$使得

$ \limsup\limits_{n\rightarrow\infty}\big\| Q^n(\tilde{\phi})- E\big\| <\sigma. $

从而存在$n_0\in{\Bbb N}$使得$\forall \theta\in[-\tau(0), 0]$,只要$n\geq n_0, $就有

$ \vert S(n\omega+\theta, \tilde{\phi})-S^*(n\omega+\theta)\vert<\sigma, ~~\vert I(n\omega+\theta, \tilde{\phi})\vert<\sigma. $

由(3.3)式知只要$t>n_0\omega$,就有$S(t, \tilde{\phi})>S^*(t)-\sigma$以及

(3.4) $ \label{**} 0<I(t)<\sigma, $

此时, $I$满足

(3.5) $ \label{35}\frac{{\rm d}I(t)}{{\rm d} t}\geq p_{\varepsilon_0}(t)I(t-\tau(t))-q(t)I(t), $

这里$p_{\varepsilon_0}(t):=(1-\tau'(t))e^{-\int_{t-\tau(t)}^t\mu(r){\rm d}r} \big(S^*(t-\tau(t))-\varepsilon_0 \big)$.

令$\eta_{\varepsilon_0}\in C([-\tau(0), 0], {\mathbb R})$为$P_{\varepsilon_0}$的对应于$r(P_{\varepsilon_0})$的一个正特征函数.由于$I(t, \tilde{\phi})>0, $$\forall t>\tau(t)$,从而可选择充分小的$\zeta>0$,使得

$ I((n_0+1)\omega+\theta, \tilde{\phi})\geq\zeta\eta_{\varepsilon_0}(\theta), ~~\forall \theta\in[-\tau(0), 0]. $

由(3.5)式和比较原理得

$ I(t, \tilde{\phi})\geq\zeta I^{\varepsilon_0}(t-(n_0+1)\omega, \eta_{\varepsilon_0}), ~~\forall t\geq (n_0+1)\omega. $

因此, $I(n\omega, \tilde{\phi})\geq\zeta I^{\varepsilon_0}\left((n-n_0-1)\omega, \eta_{\varepsilon_0}\right)= \zeta\left(r(P_{\varepsilon_0})\right)^{(n-n_0-1)}\eta_{\varepsilon_0}(0)\rightarrow+\infty$, $n\rightarrow+\infty, $这与(3.4)式矛盾,这里$I^{\varepsilon_0}\left(t, \eta_{\varepsilon_0}\right)$是满足方程

$ \frac{{\rm d}I(t)}{{\rm d} t}= p_{\varepsilon_0}(t)I(t-\tau(t))-q(t)I(t), $

且初值为$\eta_{\varepsilon_0}$的解.至此论断得证.

由以上论断知: $E$是$Y_0$中关于$Q$的一个孤立不变集,且$W^s(E)\cap Y _0=\emptyset$, $W^s(E)$是$E$的稳定集.根据一致持久性理论, $Q:Y\rightarrow Y$关于$Y_0$是一致持久的.因此,由文献[4,定理3.1.1]知,周期半流$Q(t):Y\rightarrow Y$亦关于$Y_0$是一致持久的.因此,根据文献[16,定理4.5] (令$\rho(x)={\rm dist}(x, \partial Y_0)$即可)可知系统(1.5)有一个正的$\omega$ -周期解$\big(\tilde{S}(t), \tilde{I}(t), \tilde{R}(t)\big)\in Y _0$,且$Q: Y _0\rightarrow Y _0$有一个紧的全局吸引子${\mathcal A}$.

为了证明实际持久性,采用类似于文献[17,定理4.1]的方法.定义一个连续函数$p: Y^+\rightarrow[0, \infty), $

$ p(\phi)=\phi_2(0), \quad\forall\phi= \left(\phi_1, \phi_2, \phi_3\right)\in Y_0. $

由${\mathcal A}=Q({\mathcal A})$知对任意$\phi\in{\mathcal A}$,有$\phi_2(0)>0$.令${\mathcal B}:=\bigcup\limits_{t\in[0, \omega]}Q_t\left({\mathcal A}\right)$,则对任一$\phi\in Y _0$,都有${\mathcal B}\subset Y _0$.因为${\mathcal B}$是$ Y _0$的紧子集,从而有$\inf\limits_{\phi\in{\mathcal B}}p(\phi)=\min\limits_{\phi\in{\mathcal B}}p(\phi)>0.$同时,由于$\{Q_t\}_{t\geq0}$是$\omega$ -周期半流,则对任意$\phi\in Y_0$,都有$\lim\limits_{t\rightarrow\infty} {\rm dist} (Q_t(\phi), {\mathcal B})=0$.因此,存在一个$\delta^*>0$使得

$ \liminf\limits _{t\rightarrow\infty} I(t, \phi)=\liminf\limits _{t\rightarrow\infty}p( Q_t\phi)\geq\delta^*, ~~\forall \phi \in Y_0. $

此外,由系统(1.5)中关于$S$和$R$的方程易知,存在$\delta\in(0, \delta^*)$,使得$\liminf\limits_{t\rightarrow\infty} S(t, \phi)\geq\delta$, $\liminf\limits_{t\rightarrow\infty} R(t, \phi)\geq\delta$.至此发现,存在常数$\delta>0$,使得对任一$\phi\in Y _0$,都有

$ \liminf\limits_{t\rightarrow\infty} \big(S(t, \phi), I(t, \phi), R(t, \phi)\big)\geq(\delta, \delta, \delta), $

则$\rm(ⅱ)$中结论得证.

根据定理3.1,当${\mathcal R}_0 <1$时,由$\lim\limits_{t\rightarrow \infty} \big(S(t), I(t))-(S^*(t), 0) \big)=0 $知$\lim\limits_{t\rightarrow \infty} \frac{\beta(t)S(t)I(t)}{1+c(t)I(t)}=0 $,从而依(1.4)式可得$\lim\limits_{t\rightarrow \infty} E(t)=0 $;而当${\mathcal R}_0>1$时, $S(t), I(t)$有正的最终下界,容易验证$E(t)=\int_{t-\tau(t)}^t e^{-\int_s^t\mu(r){\rm d}r} \frac{\beta(s)S(s)I(s)}{1+c(s)I(s)} {\rm d}s$亦有正的最终下界.同时,由$\tau(t)$, $\mu(t)$, $c(t)$与$\tilde{S}(\cdot), \tilde{I}(\cdot)$的周期性可以直接验证

$ \begin{eqnarray} \tilde{E}(t):=\int_{t-\tau(t)}^t e^{-\int_s^t\mu(r){\rm d}r} \frac{\beta(s)\tilde{S}(s)\tilde{I}(s)}{1+c(s)\tilde{I}(s)} {\rm d}s, \end{eqnarray} $

亦是$\omega$ -周期的.从而,由定理3.1可得如下结论.

定理3.2  假定H1成立,有如下结论:

(ⅰ)如果${\mathcal R}_0 <1$,则系统(1.2)的无病平衡$\omega$ -周期解$(S^*(t), 0, 0, 0)$在$Y$上是全局渐近吸引的;

(ⅱ)如果${\mathcal R}_0>1$,则系统(1.5)至少有一个正$\omega$ -周期解$(\tilde{S}(t), \tilde{E}(t), \tilde{I}(t), \tilde{R}(t))$,同时存在一个$\delta'>0$使得对任意$\phi\in Y_0$,有$\liminf\limits_{t\rightarrow\infty}(S(t, \phi), E(t, \phi), I(t, \phi), R(t, \phi))>(\delta', \delta', \delta', \delta').$

本节将利用数值方法模拟、验证上节结论,并分析忽略潜伏期的周期性而取其平均值对估计疾病传播能力的造成的误差和影响.由于系统(1.5)是非自治的,无法得到基本再生数的精确值,因此采用文献[18]中提出的针对周期系统基本再生数的数值计算方法近似计算基本再生数${\mathcal R}_0$.

首先,令$\omega=\pi$, $\tau(t)=0.5(1+0.5\sin(2t))$.取$\Lambda(t)=0.5(1+0.2\cos(2t+\pi/4))$, $\beta(t)=0.1(1+ 0.5\cos(2t))$, $\mu(t)=0.13$, $d(t)=0.15$, $\alpha(t)=0.2(1+0.7\cos(2t+\pi/4))$, $\gamma(t)=0.1(1+0.5\cos(2t+\pi/4))$, $c(t)=0.01(1+0.3\cos(2t+\pi/4))$.取$S(\theta)=4$, $I(\theta)=0.4$, $R(0)=0.1$, $\theta\in[-0.5, 0]$,由(1.3)式知$E(0)=0.11$.此时, ${\mathcal R}_0=0.9278 <1$,疾病最终将灭绝(图 1).若取$\beta(t)=0.25(1+ 0.5\cos(2t))$,其余参数保持不变.取初值为$S(\theta)=0.25$, $I(\theta)=0.2$, $R(0)=0.1$, $\theta\in[-0.5, 0]$,由(1.3)式知$E(0)=0.17 $.此时, ${\mathcal R}_0=2.3209>1$,系统将一致持久,并且系统的解最终稳定到一个唯一的周期解(图 2),这一性质在理论上并不易被证实.

其次,借助数值方法检验忽略潜伏期的周期性对估计疾病传播能力的影响.固定$\Lambda(t)=0.53(1+0.2\cos(2t))$, $\beta(t)=0.035(1-0.5\cos(2t))$, $\mu(t)=0.013$, $d(t)=0.85$, $\alpha(t)=0.2(1+0.7\cos(2t))$, $\gamma(t)=0.1(1+0.5\cos(2t))$, $c(t)=0.01(1-0.3\cos(2t))$.取$\tau(t)=\tau(1+0.5\cos(2t))$,其中$\tau\geq0$为常数.注意到$\tau(t)$在$[0, \pi]$区间的平均值为

$ [\tau]:=\frac{1}{\omega}\int_0^\omega \tau(t){\rm d}t=\frac{1}{\pi}\int_0^\pi \tau(1+0.5\cos(2t)){\rm d}t= \tau. $

现分别计算$\tau(t)=0.5(1+0.5\sin(2t))$和$\tau(t)=[\tau]$时随着$\tau$变化时的${\mathcal R}_0$.结果显示(图 3),对于平均潜伏期相对较短的疾病,两者基本是吻合的,但在疾病的平均潜伏期较长时,利用平均值$[\tau]$代替周期性时滞$\tau(t)$明显低估了疾病的传播能力.

本文利用动力系统的理论和方法研究了一类具有周期性潜伏期的SEIR传染病模型,得到了阈值动力学结果：当${\mathcal R}_0 <1$时,疾病将逐渐停止传播;当${\mathcal R}_0>1$时,疾病将一致持久而形成的"地方病".数值模拟的结果进一步发现当${\mathcal R}_0>1$时,一定时间以后,疾病的传播呈现出周期性.和文献[7]中结论比较,尽管采用周期性潜伏期$\tau(t)$和采用其平均化的潜伏期$[\tau]=\frac{1}{\omega}\int_{0}^{\omega}\tau(t){\rm d}t$时阈值动力学结果是一样的,但数值模拟的结果反映出在多数情况下,尤其是平均潜伏期相对较长时,采用均值$[\tau]$时对疾病传播风险的估计存在一定低估.这一发现对于疾病预防和控制能够提供一定的指导价值.