0 引言

手足口病是由肠道病毒引起的传染病, 引发手足口病的肠道病毒有20多种, 其中以柯萨奇病毒A16型(Cox A16) 和肠道病毒71型(EV 71) 最为常见^[1-2].患者、隐性感染者、无症状带毒者为主要传染源^[3].人群普遍易感, 受感后可获得免疫力, 主要传染人群为5岁以下儿童.其感染途径包括消化道, 呼吸道及接触传播.大部分的患者会伴有食欲不振、恶心、呕吐、头疼等症状, 少数人会出现严重症状, 甚至导致死亡, 还有些患者则不表现任何症状(文中称这一部分为隐性患者).目前缺乏有效治疗药物.

传染病动力学是对传染病进行理论性定量研究的一种重要方法^[4-5], 对手足口病传播的理论研究, 目前已有一些相关的数学模型.文献[6]研究了一个简单的SIR模型来预测沙捞越的手足口病患病人数和疾病的持续时间.文献[7]研究了一类具有隔离且潜伏期及感染期均有传染性的手足口病SEIQR模型, 得到了无病平衡点的全局稳定性和地方病平衡点局部稳定的条件.文献[8]研究了手足口病的最优控制策略.文献[9]建立了一个带有隐性传染的手足口病模型, 证明了无病平衡点的稳定性, 疾病的一致持续性和正平衡点的存在性.文献[10-11]研究了具有周期结构的手足口病模型, 表明了隔离措施在疾病控制上有较好的作用.文献[12]建立了具有年龄结构和隔离措施的偏微分方程手足口病模型, 讨论了无病平衡态和地方病平衡态的局部渐近稳定性.文献[6-8]中均没有考虑儿童手足口病患病者中的隐性患者, 而且只给出了地方病平衡点的局部稳定性分析.文献[9]虽然考虑了隐性患者, 但是地方病平衡点的唯一性和稳定性均没有考虑.本文在文献[6-9]手足口病模型的基础上, 考虑儿童手足口病的隐性患者, 并加入隔离措施, 给出模型的全局性态分析.首先给出模型的基本再生数, 然后通过线性化方法和Lyapunov函数方法, 讨论无病平衡点和地方病平衡点的全局渐近稳定性, 所得结论推广了以往文献中的相关结论.

1 模型建立

设儿童的总人口数量为N(t), 把儿童人群分为6个仓室:易感者, 潜伏者, 显性患者, 隐性患者, 被隔离者, 恢复者, 其人口数量分别记为S(t), E(t), I₁(t), I₂(t), Q(t)和R(t).假设A为成年人群每年的出生率, βk和β分别为显性患者和隐性患者的传染率, d表示儿童的自然死亡率, σ为儿童从潜伏者到染病者的转化系数, γ₁, γ₂, γ₃分别为显性患者, 隐性患者和被隔离者的恢复率, α₁, α₂, α₃分别为显性患者, 隐性患者和被隔离者的因病死亡率, p为儿童患者中显性患者所占的比例(0 < p < 1), q为患病者的隔离率, a为从儿童群体到成人群体的转移率.

根据以上假设, 得到如下手足口病模型

$ \left\{ \begin{array}{l} \frac{{{\rm{d}}S}}{{{\rm{d}}t}} = A - \beta S\left( {k{I_1} + {I_2}} \right) - \left( {a + d} \right)S,\\ \frac{{{\rm{d}}E}}{{{\rm{d}}t}} = \beta S\left( {k{I_1} + {I_2}} \right) - \left( {\sigma + d} \right)E,\\ \frac{{{\rm{d}}{I_1}}}{{{\rm{d}}t}} = \sigma pE - \left( {q + {\gamma _1} + {\alpha _1} + d} \right){I_1},\\ \frac{{{\rm{d}}{I_2}}}{{{\rm{d}}t}} = \sigma \left( {1 - p} \right)E - \left( {{\gamma _2} + {\alpha _2} + d} \right){I_2},\\ \frac{{{\rm{d}}Q}}{{{\rm{d}}t}} = q{I_1} - \left( {{\gamma _3} + {\alpha _3} + d} \right)Q,\\ \frac{{{\rm{d}}R}}{{{\rm{d}}t}} = {\gamma _1}{I_1} + {\gamma _2}{I_2} + {\gamma _3}Q - \left( {a + d} \right)R. \end{array} \right. $

(1)

因为总人口N(t)=S(t)+E(t)+I₁(t)+I₂(t)+Q(t)+R(t), 则有

$ \frac{{{\rm{d}}N}}{{{\rm{d}}t}} = A - dN - a\left( {S + R} \right) - {\alpha _1}{I_1} - {\alpha _2}{I_2} - {\alpha _3}Q \le A - dN, $

因此$\mathop {\lim \;{\rm{sup}}}\limits_{t \to \infty } N\left( t \right) \le \frac{A}{d}$.则模型(1) 的正向不变集为

$ \mathit{\Omega = }\left\{ {\left( {S,E,{I_1},{I_2},Q,R} \right) \in {\bf{R}}_ + ^6\left| {S + E + {I_1} + {I_2} + Q + R \le \frac{A}{d}} \right.} \right\}. $

模型(1) 总有一个无病平衡点P₀=(S⁰, 0, 0, 0, 0, 0), 其中${S^0} = \frac{A}{{a + d}}$.根据文献[13]的第3部分, 计算得到模型(1) 的基本再生数为

$ {R_0} = \frac{{\beta \sigma {S^0}}}{{\sigma + d}}\left( {\frac{{kp}}{{{\omega _1}}} + \frac{{1 - p}}{{{\omega _2}}}} \right) = \frac{{\beta \sigma {S^0}\left[ {kp{\omega _2} + \left( {1 - p} \right){\omega _1}} \right]}}{{\left( {\sigma + d} \right){\omega _1}{\omega _2}}}. $

(2)

当R₀>1时, 模型(1) 还存在唯一的地方病平衡点P^*=(S^*, E^*, I₁^*, I₂^*, Q^*, R^*), 其中

$ {S^ * } = \frac{A}{{\beta \left( {kI_1^ * + I_2^ * } \right) + a + d}},{E^ * } = \frac{{{\omega _1}I_1^ * }}{{\sigma p}},I_1^ * = \frac{{{\omega _2}p\left( {a + d} \right)\left( {{R_0} - 1} \right)}}{{\beta \left[ {kp{\omega _2} + \left( {1 - p} \right){\omega _1}} \right]}},I_2^ * = \frac{{\left( {1 - p} \right){\omega _1}I_1^ * }}{{p{\omega _2}}}, $

$ {Q^ * } = \frac{{qI_1^ * }}{{{\gamma _3} + {\alpha _3} + d}},{R^ * } = \frac{{{\gamma _1}I_1^ * + {\gamma _2}I_2^ * + {\gamma _3}{Q^ * }}}{{a + d}},{\omega _1} = q + {\gamma _1} + {\alpha _1} + d,{\omega _2} = {\gamma _2} + {\alpha _2} + d. $

2 平衡点的全局渐近稳定性

本节将利用线性化方法和构造Lyapunov函数来证明平衡点的全局渐近稳定性.

定理1  当R₀ < 1时, 无病平衡点P₀=(S⁰, 0, 0, 0, 0, 0) 全局渐近稳定; 当R₀>1时, 无病平衡点P₀=(S⁰, 0, 0, 0, 0, 0) 不稳定.

证明  对模型(1) 在P₀=(S⁰, 0, 0, 0, 0, 0) 处进行线性化, 则得到P₀=(S⁰, 0, 0, 0, 0, 0) 处的特征方程为

$ \left( {\lambda + {\gamma _3} + {\alpha _3} + d} \right){\left( {\lambda + a + d} \right)^2}\left( {{\lambda ^3} + {c_1}{\lambda ^2} + {c_2}\lambda + {c_3}} \right) = 0. $

(3)

其中

$ {c_1} = \sigma + d + {\omega _1} + {\omega _2} > 0, $

$ {c_2} = \left( {{\omega _1}\left( {\sigma + d} \right) - \sigma p\beta {S^0}k} \right) + \left( {{\omega _2}\left( {\sigma + d} \right) - \sigma \beta {S^0}\left( {1 - p} \right)} \right) + {\omega _1}{\omega _2}, $

$ {c_3} = {\omega _1}{\omega _2}\left( {\sigma + d} \right)\left( {1 - {R_0}} \right). $

显然λ₁=-(γ₃+α₃+d), λ₂=λ₃=-(a+d)为特征方程(3) 的3个特征根.当R₀>1时, c₃ < 0, 因此方程(3) 有正根, 故P₀=(S⁰, 0, 0, 0, 0, 0) 不稳定.当R₀ < 1时, c₃>0, 且$\frac{{\beta \sigma {S^0}}}{{\sigma + d}}\frac{{kp}}{{{\omega _1}}} < 1$, $\frac{{\beta \sigma {S^0}}}{{\sigma + d}}\frac{{kp}}{{{\omega _2}}} < 1$, 则有c₂>ω₁ω₂, c₁c₂-c₃>(σ+d)ω₁ω₂-ω₁ω₂(σ+d)×(1-R₀)=ω₁ω₂R₀(σ+d)>0.由Routh-Hurwitz判据知方程(3) 的所有特征根均具有负实部, 则无病平衡点P₀=(S⁰, 0, 0, 0, 0, 0) 是局部渐近稳定的.

由模型(1) 得$\frac{{{\rm{d}}S}}{{{\rm{d}}t}} \le A-\left( {a + d} \right)S$, 因此对任意给定的ε>0, 存在T>0, 当t>T时, 有0≤S(t) < $\frac{A}{{a + d}} + \varepsilon = {S^0} + \varepsilon $.当R₀ < 1时, 存在充分小的ε>0, 使得

$ \frac{{\beta \left( {{S^0} + \varepsilon } \right)k\sigma p}}{{{\omega _1}}} + \frac{{\beta \left( {{S^0} + \varepsilon } \right)\sigma \left( {1 - p} \right)}}{{{\omega _2}}} - \left( {\sigma + d} \right) < 0. $

当t>T时, 构造如下Lyapunov函数

$ V\left( t \right) = E\left( t \right) + \frac{{\beta \left( {{S^0} + \varepsilon } \right)k}}{{{\omega _1}}}{I_1}\left( t \right) + \frac{{\beta \left( {{S^0} + \varepsilon } \right)}}{{{\omega _2}}}{I_2}\left( t \right), $

则有

$ \begin{array}{*{20}{c}} {V'\left( t \right)\left| {_{\left( 1 \right)}} \right. = \left[ {\frac{{\beta \left( {{S^0} + \varepsilon } \right)k\sigma p}}{{{\omega _1}}} + \frac{{\beta \left( {{S^0} + \varepsilon } \right)\sigma \left( {1 - p} \right)}}{{{\omega _2}}} - \left( {\sigma + d} \right)} \right]E\left( t \right) + }\\ {\beta \left[ {S\left( t \right) - \left( {{S^0} + \varepsilon } \right)} \right]\left( {k{I_1}\left( t \right) + {I_2}\left( t \right)} \right) \le 0.} \end{array} $

记

${M_1} = \left\{ {\left( {S,E,{I_1},{I_2},Q,R} \right) \in \mathit{\Omega }{{\left| {\mathit{V'}\left( t \right)} \right|}_{\left( 1 \right)}} = 0} \right\} = \left\{ {\left( {S,E,{I_1},{I_2},Q,R} \right) \in \mathit{\Omega }\left| {E = {I_1} = {I_2} = 0} \right.} \right\},$

则M₁中的最大不变集为{P₀}, 由LaSalle不变集原理^[14]得当R₀ < 1时P₀=(S⁰, 0, 0, 0, 0, 0) 是全局渐近稳定的.

已知函数g(x)=x-1-lnx≥0(x∈(0, +∞)), 即1-x≤-lnx, 下面利用此不等式证明地方病平衡点P^*的全局稳定性.

定理2  当R₀>1时, 地方病平衡点P^*=(S^*, E^*, I₁^*, I₂^*, Q^*, R^*)在Ω\{P₀}内全局渐近稳定.

证明  令D₁=S-S^*-${S^*}\ln \frac{S}{{{S^*}}}$+E-E^*-${E^*}\ln \frac{E}{{{E^*}}}$, D₂=I₁-I₁^*-$I_1^*\ln \frac{{{I_1}}}{{I_1^*}}$,

$ {D_3} = {I_2} - I_2^ * - I_2^ * \ln \frac{{{I_2}}}{{I_2^ * }}. $

则有

$ \begin{array}{l} {{D'}_1}\left| {_{\left( 1 \right)}} \right. = \left( {1 - \frac{{{S^ * }}}{S}} \right)\left[ {\left( {a + d} \right){S^ * } + \beta {S^ * }\left( {kI_1^ * + I_2^ * } \right) - \beta S\left( {k{I_1} + {I_2}} \right) - \left( {a + d} \right)S} \right] + \\ \;\;\;\;\;\;\;\;\;\;\;\left( {1 - \frac{{{E^ * }}}{E}} \right)\left[ {\beta S\left( {k{I_1} + {I_2}} \right) - \frac{{\beta {S^ * }\left( {kI_1^ * + I_2^ * } \right)}}{{{E^ * }}}E} \right] = \\ \;\;\;\;\;\;\;\;\;\;\; - \frac{{a + d}}{S}{\left( {S - {S^ * }} \right)^2} + \beta k{S^ * }I_1^ * \left( {2 - \frac{{{S^ * }}}{S} - \frac{{S{E^ * }{I_1}}}{{{S^ * }EI_1^ * }} + \frac{{{I_1}}}{{I_1^ * }} - \frac{E}{{{E^ * }}}} \right) + \\ \;\;\;\;\;\;\;\;\;\;\;\beta {S^ * }I_2^ * \left( {2 - \frac{{{S^ * }}}{S} - \frac{{S{E^ * }{I_2}}}{{{S^ * }EI_2^ * }} + \frac{{{I_2}}}{{I_2^ * }} - \frac{E}{{{E^ * }}}} \right) \le \\ \;\;\;\;\;\;\;\;\;\;\;\beta k{S^ * }I_1^ * \left( {\frac{{{I_1}}}{{I_1^ * }} - \frac{E}{{{E^ * }}} - \ln \frac{{{I_1}}}{{I_1^ * }} - \ln \frac{E}{{{E^ * }}}} \right) + \\ \;\;\;\;\;\;\;\;\;\;\;\beta {S^ * }I_2^ * \left( {\frac{{{I_2}}}{{I_2^ * }} - \frac{E}{{{E^ * }}} - \ln \frac{{{I_2}}}{{I_2^ * }} - \ln \frac{E}{{{E^ * }}}} \right), \end{array} $

$ \begin{array}{l} {{D'}_2}\left| {_{\left( 1 \right)}} \right. = \left( {1 - \frac{{I_1^ * }}{{{I_1}}}} \right)\left( {\sigma pE - \frac{{\sigma pE}}{{I_1^ * }}{I_1}} \right) = \sigma p{E^ * }\left( {1 + \frac{E}{{{E^ * }}} - \frac{{{I_1}}}{{I_1^ * }} - \frac{{EI_1^ * }}{{{E^ * }{I_1}}}} \right) \le \\ \;\;\;\;\;\;\;\;\;\;\;\sigma p{E^ * }\left( {\frac{E}{{{E^ * }}} - \frac{{{I_1}}}{{I_1^ * }} - \ln \frac{E}{{{E^ * }}} + \ln \frac{{{I_1}}}{{I_1^ * }}} \right), \end{array} $

$ \begin{array}{l} {{D'}_2}\left| {_{\left( 1 \right)}} \right. = \left( {1 - \frac{{I_1^ * }}{{{I_1}}}} \right)\left[ {\sigma \left( {1 - p} \right)E - \frac{{\sigma \left( {1 - p} \right){E^ * }}}{{I_2^ * }}{I_2}} \right] = \\ \;\;\;\;\;\;\;\;\;\;\;\sigma \left( {1 - p} \right){E^ * }\left( {1 + \frac{E}{{{E^ * }}} - \frac{{{I_2}}}{{I_2^ * }} - \frac{{EI_2^ * }}{{{E^ * }{I_2}}}} \right) \le \\ \;\;\;\;\;\;\;\;\;\;\;\sigma \left( {1 - p} \right){E^ * }\left( {\frac{E}{{{E^ * }}} - \frac{{{I_2}}}{{I_2^ * }}\ln \frac{E}{{{E^ * }}} + \ln \frac{{{I_2}}}{{I_2^ * }}} \right). \end{array} $

令D=a₁D₁+a₂D₂+a₃D₃, 其中a_i(i=1, 2, 3) 为待定的正常数, 得

$ \begin{array}{l} {{D'}_1}\left| {_{\left( 1 \right)}} \right. \le \left( {\frac{{{I_1}}}{{I_1^ * }} - \ln \frac{{{I_1}}}{{I_1^ * }}} \right)\left( {{a_1}\beta k{S^ * }I_1^ * - {a_2}\sigma p{E^ * }} \right) + \\ \;\;\;\;\;\;\;\;\;\;\;\left( {\frac{{{I_2}}}{{I_2^ * }} - \ln \frac{{{I_2}}}{{I_2^ * }}} \right)\left[ {{a_1}\beta {S^ * }I_2^ * - {a_3}\sigma \left( {1 - p} \right){E^ * }} \right] + \\ \;\;\;\;\;\;\;\;\;\;\;\left( {\frac{E}{{{E^ * }}} - \ln \frac{E}{{{E^ * }}}} \right)\left[ { - {a_1}\beta k{S^ * }I_1^ * + {a_2}\sigma p{E^ * } - {a_1}\beta {S^ * }I_2^ * - {a_3}\sigma \left( {1 - p} \right){E^ * }} \right]. \end{array} $

取a₁=σp(1-p)E^*, a₂=βkS^*I₁^*(1-p), a₃=βS^*I₂^*p, 得D′|₍₁₎≤0, 记

$ {M_2} = \left\{ {\left( {S,E,{I_1},{I_2},Q,R} \right) \in \mathit{\Omega }{{\left| {D'} \right|}_{\left( 1 \right)}}{\rm{ = }}0} \right\}{\rm{ = }}\left\{ {\left( {S,E,{I_1},{I_2},Q,R} \right) \in \mathit{\Omega }\left| {S = {S^ * },} \right.\frac{{{E^ * }}}{E} = \frac{{I_1^ * }}{{{I_1}}}{\rm{ = }}\frac{{I_2^ * }}{{{I_2}}}} \right\}, $

则M₂中的最大不变集为{P^*}, 由LaSalle不变集原理^[14]知当R₀>1时, P^*=(S^*, E^*, I₁^*, I₂^*, Q^*, R^*)全局渐近稳定.

当儿童患者中显性患者所占比例p=1时, 即所有患病儿童均为显性患者, 文献[7, 9, 15-17]考虑的模型均是这种情况, 此时模型(1) 变为如下形式

$ \left\{ \begin{array}{l} \frac{{{\rm{d}}S}}{{{\rm{d}}t}} = A - \beta kS{I_1} - \left( {a + d} \right)S,\\ \frac{{{\rm{d}}E}}{{{\rm{d}}t}} = \beta kS{I_1} - \left( {\sigma + d} \right)E,\\ \frac{{{\rm{d}}{I_1}}}{{{\rm{d}}t}} = \sigma E - \left( {q + {\gamma _1} + {\alpha _1} + d} \right){I_1},\\ \frac{{{\rm{d}}Q}}{{{\rm{d}}t}} = q{I_1} - \left( {{\gamma _3} + {\alpha _3} + d} \right)Q,\\ \frac{{{\rm{d}}R}}{{{\rm{d}}t}} = {\gamma _1}{I_1} + {\gamma _3}Q - \left( {a + d} \right)R. \end{array} \right. $

(4)

模型(4) 的正向不变集为

$ \mathit{\tilde \Omega = }\left\{ {\left( {S,E,{I_1},Q,R} \right) \in {\bf{R}}_ + ^5\left| {S,E,{I_1},Q,R} \right. \le \frac{A}{d}} \right\}. $

模型(4) 总有一个无病平衡点${\tilde P_0} = \left( {{S^0}, 0, 0, 0, 0} \right)$.模型(4) 的基本再生数为${\tilde R_0} = \frac{{\beta k\sigma {S^0}}}{{\left( {\sigma + d} \right){\omega _1}}}$.当R₀>1时, 模型(4) 还存在一个地方病平衡点${\tilde P^*} = \left( {{{\tilde S}^*}, {{\tilde E}^*}, {{\tilde I}^*}, {{\tilde Q}^*}, {{\tilde R}^*}} \right)$, 其中

$ {{\tilde S}^ * } = \frac{A}{{\beta k\tilde I_1^ * + a + d}},{{\tilde E}^ * } = \frac{{{\omega _1}\tilde I_1^ * }}{\sigma },\tilde I_1^ * = \frac{{\left( {a + d} \right)\left( {{{\tilde R}_0} - 1} \right)}}{{\beta k}},{{\tilde Q}^ * } = \frac{{q\tilde I_1^ * }}{{{\gamma _3} + {\alpha _3} + d}},{{\tilde R}^ * } = \frac{{{\gamma _1}\tilde I_1^ * + {\gamma _3}{{\tilde Q}^ * }}}{{a + d}}. $

当${\tilde R_0} > 1$时, 构造如下Lyapunov函数

$ L = S - {S^ * } - {S^ * }\ln \frac{S}{{{S^ * }}} + E - {E^ * } - {E^ * }\ln \frac{E}{{{E^ * }}} + \frac{{\sigma + d}}{\sigma }\left( {{I_1} - I_1^ * - I_1^ * \ln \frac{{{I_1}}}{{I_1^ * }}} \right). $

与证明定理1和定理2类似, 对模型(4) 的平衡点, 有如下全局稳定性结论.

定理3  当${\tilde R_0} < 1$时, 无病平衡点${\tilde P_0} = \left( {{S^0}, 0, 0, 0, 0} \right)$全局渐近稳定; 当${\tilde R_0} > 1$时, 无病平衡点${\tilde P_0} = \left( {{S^0}, 0, 0, 0, 0} \right)$不稳定, 地方病平衡点${\tilde P^*} = \left( {{{\tilde S}^*}, {{\tilde E}^*}, \tilde I_1^*, {{\tilde Q}^*}, {{\tilde R}^*}} \right)$在$\mathit{\tilde \Omega \backslash }\left\{ {{{\tilde P}_0}} \right\}$内全局渐近稳定.

3 结束语

本文讨论了具有隐性传染和隔离措施, 且显性患者和隐性患者都有传染性的手足口病传染病模型, 得到了模型的基本再生数, 此基本再生数完全决定了模型的动力学行为.利用特征根方法和Lyapunov函数方法得到了模型的无病平衡点和地方病平衡点的全局渐近稳定性.