0 引言

近年来，越来越多的学者利用数学模型分析宿主细胞和病毒之间的相互作用，研究过的病毒有人类免疫缺陷病毒(HIV), 乙肝病毒(HBV)，丙肝病毒(HCV), SEIR传染病毒等^[1-7].为了描述易感染细胞，感染细胞以及乙肝病毒颗粒之间的关系，建立了如下基本的病毒动力学模型：

$ \left\{ {\begin{array}{*{20}{l}} {\frac{{{\rm{d}}x}}{{{\rm{d}}t}} = \lambda - ax - \beta xz,}\\ {\frac{{{\rm{d}}y}}{{{\rm{d}}t}} = \beta xz - by,}\\ {\frac{{{\rm{d}}v}}{{{\rm{d}}t}} = cy - dv,} \end{array}} \right. $

(1)

然而要建立更精确的病毒感染数学模型，必须考虑免疫反应.免疫系统中对病毒感染有影响的两个主要反应是细胞免疫和体液免疫.细胞免疫中细胞毒性T淋巴细胞在病毒防御中起着关键的作用，因为细胞毒性T淋巴细胞可以攻击并杀死被感染细胞.而体液免疫是基于B细胞产生的抗体攻击并杀死被感染细胞.在一些病毒感染中，比如疟疾，细胞免疫比体液免疫的效果差^[8-10].Muras等^[9]已经建立了体液免疫的基本动力学模型：

$ \left\{ {\begin{array}{*{20}{l}} {\frac{{{\rm{d}}x}}{{{\rm{d}}t}} = \lambda - dx - \beta xv,}\\ {\frac{{{\rm{d}}y}}{{{\rm{d}}t}} = \beta xv - by,}\\ {\frac{{{\rm{d}}v}}{{{\rm{d}}t}} = ky - uv - rzv,}\\ {\frac{{{\rm{d}}z}}{{{\rm{d}}t}} = gzv - \mu z.} \end{array}} \right. $

(2)

其中x, y, v分别表示易感染细胞，受感染细胞以及游离病毒颗粒的数量，z表示B细胞的数量.未感染细胞以常速率λ产生，死亡速率dx，细胞被感染的速率为βxv，受感染的细胞死亡率为by.游离病毒颗粒从受感染的细胞中产生的速率为ky，死亡率为uv，同时被抗体作用移除体内的速率为rzv.B细胞被激活的速率为gzv，死亡率为μz.所有系数为正.

但是系统(2)并未考虑潜伏的感染细胞(细胞的感染率以双线性函数βxv的形式给出).2009年，Gang Huang等^[11]建立具有Bsddington-DeAngelis型功能性反应函数的病毒动力学模型，用更具一般形式的功能性反应函数$ \frac{{\beta xv}}{{1 + \gamma x + \alpha v}} $代替模型(2)中双线性函数βxv.从生物学的角度来看, 该模型更能描述病毒与正常细胞之间的动态演化作用.然而文献[10]未考虑潜伏的感染细胞和体液免疫反应.文献[12]在文献[11]的基础上，建立考虑潜伏感染细胞并具有Bsddington-DeAngelis型功能性反应函数的病毒动力学模型，但未考虑体液免疫反应。本文在文献[8-12]的启发下，建立具有Bsddington-DeAngelis型功能性反应函数传染病模型，既考虑体液免疫反应又考虑潜伏的感染细胞，模型如下：

$ \left\{ {\begin{array}{*{20}{l}} {\frac{{{\rm{d}}x}}{{{\rm{d}}t}} = \lambda - dx - \frac{{\beta xv}}{{1 + \gamma x + \alpha v}},}\\ {\frac{{{\rm{d}}x}}{{{\rm{d}}t}} = \frac{{\left( {1 - q} \right)\beta xv}}{{1 + \gamma x + \alpha v}} - \left( {e + \delta } \right)w,}\\ {\frac{{{\rm{d}}y}}{{{\rm{d}}t}} = \frac{{q\beta xv}}{{1 + \gamma x + \alpha v}} + \delta w - by,}\\ {\frac{{{\rm{d}}v}}{{{\rm{d}}t}} = ky - uv - rzv,}\\ {\frac{{{\rm{d}}z}}{{{\rm{d}}t}} = gzv - \mu z.} \end{array}} \right. $

(3)

其中γ和α是正常数，w和y分别表示潜伏的感染细胞和活动性感染细胞的浓度.潜伏感染细胞死亡的速率ew，转化为活动性感染细胞的速率δw.1－q和q(0 < q < 1)分别表示受病毒感染的细胞中变成潜伏细胞或活动细胞的概率，其他变量与系统(2)有相同的含义.

1 模型基本性质

系统(3)的所有解都非负有界.解的非负性从生物学意义上讲显然成立，只需证明其有界性.

定理1系统(3)的所有解x(t), w(t), y(t), v(t), z(t)都是非负且有界的，即存在M_i > 0, i=1, 2, 3使得0≤x(t), w(t), y(t)≤M₁，0≤v(t)≤M₂，0≤z(t)≤M₃成立.

证明设X₁(t)=x(t)+w(t)+y(t)，则有

$ \frac{{{\rm{d}}{X_1}\left( t \right)}}{{{\rm{d}}t}} = \lambda-dx\left( t \right)-ew\left( t \right)by\left( t \right) < \lambda-\\ {s_1}{X_1}\left( t \right). $

这里s₁=min {d, b, e}.因此，对充分大的t，有X₁(t) < M₁+ε(ε为任意小的正数)，其中$ {M_1} = \frac{\lambda }{{{s_1}}} $.所以0≤x(t), w(t), y(t)≤M₁.

另一方面，设$ {X_2} = v\left( t \right) + \frac{r}{g}z\left( t \right), $则有

$ \frac{{{\rm{d}}{X_2}\left( t \right)}}{{{\rm{d}}t}} = ky\left( t \right)-uv\left( t \right)-\frac{{r\mu }}{g}z\left( t \right) \le k{M_1}-\\ {s_2}\left( {v\left( t \right) + \frac{r}{g}z\left( t \right)} \right) = k{M_1} -{s_2}{X_2}\left( t \right). $

这里s₂=min {u, μ}.因此，对充分大的t，有X₂(t) < M₂+ε(ε为任意小的正数)，其中$ {M_2} = \frac{{k{M_1}}}{{{s_2}}} $.所以0≤v(t)≤M₂，0≤z(t)≤M₃，其中$ {M_3} = \frac{{g{M_2}}}{r} $.

2 模型无病平衡点的稳定性

为方便，引入函数H:(0, ∞)→[0, ∞)，定义为H(s)=s－1－lns.基本再生数R₀定义为$ {R_0} = \frac{{k\beta {x_0}\left( {eq + \delta } \right)}}{{bu\left( {e + \delta } \right)\left( {1 + \gamma {x_0}} \right)}}, $其中$ {x_0} = \frac{\lambda }{d}. $

定理2当R₀≤1时，系统(3)的无病平衡点$ {E_0} = \left( {\frac{\lambda }{d}, 0, 0, 0, 0} \right) $是全局渐近稳定的.

证明构造Lyapunov函数：

$ \begin{array}{*{20}{l}} {\;\;\;\;\;\;\;{W_0} = \frac{{{x_0}}}{{1 + \gamma {x_0}}}H\left( {\frac{x}{{{x_0}}}} \right) + \frac{\delta }{{eq + \delta }}w + \frac{{e + \delta }}{{eq + \delta }}y}\\ {\frac{{b\left( {e + \delta } \right)}}{{k\left( {eq + \delta } \right)}}v + \frac{{br\left( {e + \delta } \right)}}{{kg\left( {eq + \delta } \right)}}z.} \end{array} $

结合系统(3)有

$ \begin{array}{*{20}{l}} {\frac{{{\rm{d}}{W_0}}}{{{\rm{d}}t}} = \frac{1}{{1 + \gamma {x_0}}}\left( {1 - \frac{{{x_0}}}{x}} \right)\left( {\lambda - dx - } \right.}\\ {\left. {\frac{{\beta xv}}{{1 + \gamma x + \alpha v}}} \right) + \frac{\delta }{{eq + \delta }}\left( {\frac{{\left( {1 - q} \right)\beta xv}}{{1 + \gamma x + \alpha v}} - \left( {e + } \right.} \right.}\\ {\left. {\left. \delta \right)w} \right) + \frac{{e + \delta }}{{eq + \delta }}\left( {\frac{{q\beta xv}}{{1 + \gamma x + \alpha v}} + \delta w - by} \right) + }\\ {\frac{{b\left( {e + \delta } \right)}}{{k\left( {eq + \delta } \right)}}\left( {ky - uv - rzv} \right) + \frac{{br\left( {e + \delta } \right)}}{{kg\left( {eq + \delta } \right)}}\left( {gzv - } \right.}\\ {\left. {\mu z} \right) = - d\frac{{{{\left( {x - {x_0}} \right)}^2}}}{{x\left( {1 + \gamma {x_0}} \right)}} - \frac{{\beta xv}}{{\left( {1 + \gamma {x_0}} \right)\left( {1 + \gamma x + \alpha v} \right)}} + }\\ {\frac{{\beta {x_0}v}}{{\left( {1 + \gamma {x_0}} \right)\left( {1 + \gamma x + \alpha v} \right)}} + \frac{{\beta xv}}{{1 + \gamma x + \alpha v}} - }\\ {\frac{{bu\left( {e + \delta } \right)}}{{k\left( {eq + \delta } \right)}}v - \frac{{br\mu \left( {e + \delta } \right)}}{{kg\left( {eq + \delta } \right)}}z = - d\frac{{{{\left( {x - {x_0}} \right)}^2}}}{{x\left( {1 + \gamma {x_0}} \right)}} + }\\ {\frac{{\beta {x_0}\left( {1 + \gamma x} \right)}}{{\left( {1 + \gamma {x_0}} \right)\left( {1 + \gamma x + \alpha v} \right)}} - \frac{{bu\left( {e + \delta } \right)}}{{k\left( {eq + \delta } \right)}}v - }\\ {\frac{{br\mu \left( {e + \delta } \right)}}{{kg\left( {eq + \delta } \right)}}z = - d\frac{{{{\left( {x - {x_0}} \right)}^2}}}{{x\left( {1 + \gamma {x_0}} \right)}} - \frac{{bu\left( {e + \delta } \right)}}{{k\left( {eq + \delta } \right)}} \bullet }\\ {\frac{{\alpha {v^2}{R_0}}}{{1 + \gamma x + \alpha v}} + \frac{{bu\left( {e + \delta } \right)}}{{k\left( {eq + \delta } \right)}}\left( {{R_0} - 1} \right)v - \frac{{br\mu \left( {e + \delta } \right)}}{{kg\left( {eq + \delta } \right)}}z.} \end{array} $

(4)

因此，当R₀≤1时，对任意的x, v, z>0有$ \frac{{{\rm{d}}{W_0}}}{{{\rm{d}}t}} \le 0 $成立.又易知系统(3)的解收敛于一个集合Γ，其中Γ是$ \left\{ {\frac{{{\rm{d}}{W_0}}}{{{\rm{d}}t}} \le 0} \right\} $的最大不变子集.同时由(4)式知，$ {\frac{{{\rm{d}}{W_0}}}{{{\rm{d}}t}} = 0} $当且仅当x=x₀, v=0, z=0成立.注意到，对Γ中的任何元素，如果v=z=0，就有$ {\frac{{{\rm{d}}{W_0}}}{{{\rm{d}}t}} = 0} $.又由系统(3)的第4个方程知，$ 0 = \frac{{{\rm{d}}v}}{{{\rm{d}}t}} = ky $，因此y=0.此外，由系统(3)的第3个方程知w=0.因此，$ {\frac{{{\rm{d}}{W_0}}}{{{\rm{d}}t}} = 0} $当且仅当x=x₀, y=0, v=0, z=0成立.故由时滞微分方程的LaSalle不变原理知，无病平衡点$ {E_0} = \left( {\frac{\lambda }{d}, 0, 0, 0, 0} \right) $)是全局渐近稳定的.

3 结论

本文讨论了一类考虑潜伏感染细胞的具有体液免疫反应的传染病模型，该模型描述了未受感染的靶细胞，潜伏感染的细胞，活动性感染的细胞，自由的病毒颗粒与B细胞之间的相互作用.利用Lyapunov函数和LaSalle不变原理，证明了当基本再生数R₀≤1时，模型无病平衡点全局渐近稳定.