具有体液免疫和凋亡率的病毒模型分析

Analysis of Viral Infection Models with Humoral Immunity and Apoptosis Rates

摘要

为了研究体液免疫和感染细胞诱导未感染细胞凋亡等因素对病毒传播的影响，本文提出一类具有体液免疫，细胞间传播和凋亡率的病毒感染动力学模型。首先，给出模型解的非负性和有界性，通过简单计算，得到模型的基本再生数R₀和体液免疫激活基本再生数R_H。其次，利用Routh-Hurwitz判据和Lyapunov-LaSalle定理，证明了无病平衡点和免疫未激活平衡点的局部和全局渐近稳定性。最后，通过数值模拟验证了理论结果的正确性，发现凋亡率的增加对病毒感染的控制具有积极作用。

Abstract

In order to investigate the effects of humoral immunity and induction of apoptosis in uninfected cells by infected cells on viral transmission, this paper presents a dynamic model of viral infection with humoral immunity, cell-to-cell transmission and apoptosis rate. Firstly, the non-negativity and boundedness of the model solutions are given, and the basic reproduction number R₀ and the humoral immune activation basic reproduction number R_H of the model are obtained by simple calculations. Secondly, the local and global asymptotic stability of the disease-free equilibrium and the immune response free equilibrium are proved by using the Routh-Hurwitz criterion and Lyapunov-LaSalle theorem. Finally, by numerical simulation verified the correctness of the theoretical results, and the increase in apoptosis rate was found to have a positive effect on the control of viral infection.

Graphical abstract

关键词

病毒感染模型 / 细胞间传播 / 细胞诱导凋亡 / 体液免疫

Key words

viral infection model / cell-to-cell transmission / induced apoptosis / humoral immunity

0 引言

众所周知，病毒与疾病密切相关，大多数的传染病都是由病毒引起的。例如人类免疫缺陷病毒（Human Immunodeficiency Virus， HIV）、乙型肝炎病毒（Hepatitis B Virus， HBV）、丙型肝炎病毒（Hepatitis C Virus， HCV）、埃博拉病毒（Ebola）、COVID-19等。病毒传播主要有病毒对细胞感染^［1］和细胞间传播^［2］两种途径。因此，在建立数学模型时，不仅要考虑病毒对细胞的感染，也要考虑细胞间的传播。例如Guo等^［3］建立了一个细胞间传播和病毒对细胞感染的HIV感染动力学模型，发现细胞间传播明显增加了感染细胞浓度。Pan等^［4］讨论了一个包括感染和传播途径的HCV感染模型，并对四种不同的HCV模型进行了数值解释。Wang等^［5］和Luo等^［6］对细胞间传播和病毒对细胞感染进行了研究，发现随着细胞传播率的增加，感染细胞和病毒粒子浓度增加，健康细胞浓度减少。

当病毒进入人体后，免疫系统发挥重要作用，它可以有效地控制病毒在宿主体内的复制和传播^［7-8］，其中免疫系统由细胞免疫和体液免疫两部分组成。近年来，对于病毒动力学与细胞免疫的数学模型，研究者们取得了较多成果^［9-10］。但是，在一些感染中，体液免疫比细胞免疫更加有效^［11-13］。例如Wang等^［12］指出，当病毒在宿主内复制和传播时，体液免疫更为重要。

此外，研究者们发现可通过诱导细胞凋亡的途径来杀死被病毒入侵的细胞。关于考虑细胞凋亡的病毒模型的研究已有一些成果^［14-19］。例如Fan等^［15］在模型中引入细胞凋亡，发现细胞凋亡可以促进或缓解HIV感染的长期演变。Guo等^［17］在HIV感染时滞模型中考虑了细胞凋亡这一因素，并分析了模型的全局动力学。然而，上述模型大都是单一考虑细胞间传播、细胞凋亡或体液免疫的影响，为了综合考虑上述因素对病毒感染动力学的影响，本文提出一类具有体液免疫、细胞间传播和凋亡率的病毒感染动力学模型，进行研究讨论。

本文结构如下。在第1节中介绍了模型的建立，并给出了解的正性和有界性。第2节中，通过计算得到模型的基本再生数R_{\mathrm{0}}和体液免疫激活基本再生数R_H，并利用Routh-Hurwitz判据给出了无病平衡点和免疫未激活平衡点的局部渐近稳定性。无病平衡点和免疫未激活平衡点的全局渐近稳定性在第3节中给出。第4节中，通过数值模拟对所得的理论结果进行验证。在第5节中给出了结论。

1 建立模型

考虑体液免疫，细胞间传播以及细胞诱导凋亡等因素，利用图1来说明各仓室之间的传播动力学。

建立对应的微分方程模型如下：

\left\{\begin{array}{l}\frac{\mathrm{d}T}{\mathrm{d}t}=\lambda -{\beta }_{\mathrm{1}}TV-{\beta }_{\mathrm{2}}TI-\alpha TI-d_{\mathrm{1}}T,\\ \frac{\mathrm{d}I}{\mathrm{d}t}={\beta }_{\mathrm{1}}TV+{\beta }_{\mathrm{2}}TI-d_{\mathrm{2}}I,\\ \frac{\mathrm{d}V}{\mathrm{d}t}=kI-d_{\mathrm{3}}V-pVZ,\\ \frac{\mathrm{d}Z}{\mathrm{d}t}=cVZ-d_{\mathrm{4}}Z,\end{array}\right.（1）

其中T，I和V分别代表t时刻未感染细胞、感染细胞和游离病毒粒子的浓度，Z代表体液免疫反应（B细胞或抗体）。未感染细胞以速率\lambda从宿主内产生，d_{\mathrm{1}},d_{\mathrm{2}},d_{\mathrm{3}},d_{\mathrm{4}}分别表示未感染细胞、感染细胞、病毒和B细胞的自然死亡率，病毒感染未感染细胞的速率为{\beta }_{\mathrm{1}}，感染细胞感染未感染细胞的速率为{\beta }_{\mathrm{2}}，\alpha表示感染细胞诱导未感染细胞的凋亡率，感染细胞产生游离病毒粒子速率为k，病毒粒子进入后，体液免疫反应通过免疫细胞以c的速率产生而被触发，游离病毒粒子以速率p被B细胞中和。出于生物学上的考虑，假设模型的所有参数都是正的，且系统（1）的初始条件为(T(\mathrm{0}),I(\mathrm{0}),V(\mathrm{0}),Z(\mathrm{0}))\in R_{+}^{\mathrm{4}}，并有以下引理成立。

引理1^［4］ 系统（1）具有非负初始条件(T(\mathrm{0}),I(\mathrm{0}),V(\mathrm{0}),Z(\mathrm{0}))\in R_{+}^{\mathrm{4}}的解是正的，并且是最终有界的。

2 平衡点的存在性和局部稳定性

通过简单计算，系统（1）在任意平衡点E_i(T_i,I_i,V_i,Z_i)处的Jacobi矩阵为：

\mathit{J}(E_i)=\left[\begin{array}{cccc}-{\beta }_{\mathrm{1}}{V}_i-{\beta }_{\mathrm{2}}{I}_i-\alpha I_i-d_{\mathrm{1}}& -{\beta }_{\mathrm{2}}{T}_i-\alpha T_i& -{\beta }_{\mathrm{1}}{T}_i& \mathrm{0}\\ {\beta }_{\mathrm{1}}{V}_i+{\beta }_{\mathrm{2}}{I}_i& {\beta }_{\mathrm{2}}{T}_i-d_{\mathrm{2}}& {\beta }_{\mathrm{1}}{T}_i& \mathrm{0}\\ \mathrm{0}& k& -d_{\mathrm{3}}-p{Z}_i& -p{V}_i\\ \mathrm{0}& \mathrm{0}& \mathrm{0}& c{V}_i-d_{\mathrm{4}}\end{array}\right]。（2）

显然，系统（1）存在无病平衡点E_{\mathrm{0}}=(T_{\mathrm{0}},I_{\mathrm{0}},V_{\mathrm{0}},Z_{\mathrm{0}})=\left(\frac{\lambda }{d_{\mathrm{1}}},\mathrm{0,0},\mathrm{0}\right)。

利用下一代矩阵法^［20］计算得系统（1）的基本再生数R_{\mathrm{0}}：

R_{\mathrm{0}}=\rho (\mathrm{ℱ}{\mathrm{𝒱}}^{-\mathrm{1}})=\frac{\lambda k{\beta }_{\mathrm{1}}+\lambda {\beta }_{\mathrm{2}}{d}_{\mathrm{3}}}{d_{\mathrm{1}}{d}_{\mathrm{2}}{d}_{\mathrm{3}}}=R_{\mathrm{01}}+R_{\mathrm{02}}，

其中R_{\mathrm{01}}=\frac{k{\beta }_{\mathrm{1}}{T}_{\mathrm{0}}}{d_{\mathrm{2}}{d}_{\mathrm{3}}}，R_{\mathrm{02}}=\frac{{\beta }_{\mathrm{2}}{T}_{\mathrm{0}}}{d_{\mathrm{2}}}。

注1 在生物学上R_{\mathrm{01}}和R_{\mathrm{02}}分别表示病毒对细胞和细胞对细胞的感染相对应的基本再生数。

定理1 对于系统（1），当时，无病平衡点E_{\mathrm{0}}是局部渐近稳定的。

证明 将无病平衡点E_{\mathrm{0}}=(T_{\mathrm{0}},I_{\mathrm{0}},V_{\mathrm{0}},Z_{\mathrm{0}})代入（2）式中，经过简单计算可得E_{\mathrm{0}}处的特征方程为：

(r+d_{\mathrm{1}})(r+d_{\mathrm{4}})(r^{\mathrm{2}}+A_{\mathrm{1}}r+A_{\mathrm{2}})=\mathrm{0}，（3）

其中

A_{\mathrm{1}}=d_{\mathrm{2}}+d_{\mathrm{3}}-\frac{\lambda {\beta }_{\mathrm{2}}}{d_{\mathrm{1}}}=d_{\mathrm{3}}+d_{\mathrm{2}}(\mathrm{1}-R_{\mathrm{0}}+R_{\mathrm{01}})，

A_{\mathrm{2}}=d_{\mathrm{2}}{d}_{\mathrm{3}}-\frac{\lambda }{d_{\mathrm{1}}}(k{\beta }_{\mathrm{1}}+{\beta }_{\mathrm{2}}{d}_{\mathrm{3}})=d_{\mathrm{3}}{d}_{\mathrm{2}}(\mathrm{1}-R_{\mathrm{0}})，

显然，特征方程（3）的两个特征值分别为r_{\mathrm{1}}=-d_{\mathrm{1}}，r_{\mathrm{2}}=-d_{\mathrm{4}}，均为负数。当时，A_{\mathrm{1}}>\mathrm{0}，A_{\mathrm{2}}>\mathrm{0}。因此，根据Routh-Hurwitz判据^［4］得，当时，无病平衡点E_{\mathrm{0}}是局部渐近稳定的，R_{\mathrm{0}}>\mathrm{1}时E_{\mathrm{0}}是不稳定的。

\left.\mathrm{i}\right) 当R_{\mathrm{0}}>\mathrm{1}时，系统（1）存在免疫未激活平衡点E_{\mathrm{1}}：

E_{\mathrm{1}}=(T_{\mathrm{1}},I_{\mathrm{1}},V_{\mathrm{1}},Z_{\mathrm{1}})=\left(\frac{d_{\mathrm{2}}{d}_{\mathrm{3}}}{{\beta }_{\mathrm{1}}k+{\beta }_{\mathrm{2}}{d}_{\mathrm{3}}},\frac{d_{\mathrm{1}}{T}_{\mathrm{1}}}{\alpha T_{\mathrm{1}}+d_{\mathrm{2}}}(R_{\mathrm{0}}-\mathrm{1}),\frac{k}{d_{\mathrm{3}}}{I}_{\mathrm{1}},\mathrm{0}\right)。

定义参数R_{\mathrm{1}}=\frac{ck}{d_{\mathrm{3}}{d}_{\mathrm{4}}}{I}^{\mathrm{*}}，R_{\mathrm{1}}表示在感染的慢性阶段，未发生体液免疫反应的病毒粒子数量。

\left.\mathrm{i}\mathrm{i}\right)当R_{\mathrm{1}}>\mathrm{1}时，系统（1）存在免疫激活平衡点E^{\mathrm{*}}：

E^{\mathrm{*}}=(T^{\mathrm{*}},I^{\mathrm{*}},V^{\mathrm{*}},Z^{\mathrm{*}})=\left(\frac{d_{\mathrm{2}}{I}^{\mathrm{*}}}{{\beta }_{\mathrm{1}}{V}^{\mathrm{*}}+{\beta }_{\mathrm{2}}{I}^{\mathrm{*}}},\frac{-m_{\mathrm{2}}+\sqrt[]{{m_{\mathrm{2}}}^{\mathrm{2}}+\mathrm{4}{m}_{\mathrm{1}}{m}_{\mathrm{3}}}}{\mathrm{2}{m}_{\mathrm{1}}},\frac{d_{\mathrm{4}}}{c},\frac{d_{\mathrm{3}}}{p}(R_{\mathrm{1}}-\mathrm{1})\right)，

其中m_{\mathrm{1}}=c{\beta }_{\mathrm{2}}{d}_{\mathrm{2}}+c\alpha d_{\mathrm{2}}，m_{\mathrm{2}}={\beta }_{\mathrm{1}}{d}_{\mathrm{2}}{d}_{\mathrm{4}}+c{d}_{\mathrm{1}}{d}_{\mathrm{2}}-\lambda c{\beta }_{\mathrm{2}}，m_{\mathrm{3}}=\lambda {\beta }_{\mathrm{1}}{d}_{\mathrm{4}}。

注2 定义R_H为体液免疫激活基本再生数，其中

R_H=\frac{\lambda ck(k{\beta }_{\mathrm{1}}+{\beta }_{\mathrm{2}}{d}_{\mathrm{3}})}{ck{d}_{\mathrm{1}}{d}_{\mathrm{2}}{d}_{\mathrm{3}}+d_{\mathrm{2}}{d}_{\mathrm{3}}{d}_{\mathrm{4}}[k{\beta }_{\mathrm{1}}+(\alpha +{\beta }_{\mathrm{2}})d_{\mathrm{3}}]}，

表示在存在体液免疫反应的情况下产生的二次感染细胞的平均数量。

引理2(\mathrm{i})R_{\mathrm{1}}>\mathrm{1}\iff R_H>\mathrm{1}，(\mathrm{i}\mathrm{i})R_{\mathrm{1}}=\mathrm{1}\iff R_H=\mathrm{1}，(\mathrm{i}\mathrm{i}\mathrm{i})。

证明 (\mathrm{i})因为R_{\mathrm{1}}>\mathrm{1}\iff I^{\mathrm{*}}>\frac{d_{\mathrm{3}}{d}_{\mathrm{4}}}{ck}\iff \frac{-m_{\mathrm{2}}+\sqrt[]{{m_{\mathrm{2}}}^{\mathrm{2}}+\mathrm{4}{m}_{\mathrm{1}}{m}_{\mathrm{3}}}}{\mathrm{2}{m}_{\mathrm{1}}}>\frac{d_{\mathrm{3}}{d}_{\mathrm{4}}}{ck}，即

R_{\mathrm{1}}>\mathrm{1}\iff ({m_{\mathrm{2}}}^{\mathrm{2}}+\mathrm{4}{m}_{\mathrm{1}}{m}_{\mathrm{3}})-{\left(\frac{\mathrm{2}{m}_{\mathrm{1}}{d}_{\mathrm{3}}{d}_{\mathrm{4}}}{ck}+m_{\mathrm{2}}\right)}^{\mathrm{2}}>\mathrm{0}，

化简不等式得到

R_{\mathrm{1}}>\mathrm{1}\iff \frac{\mathrm{4}({\beta }_{\mathrm{2}}+\alpha )d_{\mathrm{2}}{d}_{\mathrm{4}}}{k^{\mathrm{2}}}[ck{d}_{\mathrm{1}}{d}_{\mathrm{2}}{d}_{\mathrm{3}}+d_{\mathrm{2}}{d}_{\mathrm{3}}{d}_{\mathrm{4}}({\beta }_{\mathrm{1}}k+({\beta }_{\mathrm{2}}+\alpha )d_{\mathrm{3}})](R_H-\mathrm{1})>\mathrm{0}，

因此，R_{\mathrm{1}}>\mathrm{1}\iff R_H>\mathrm{1}。 同理可证(\mathrm{i})和(\mathrm{i}\mathrm{i})。 因为，所以由引理2可得，当时，当R_{\mathrm{1}}>\mathrm{1}时R_{\mathrm{0}}>\mathrm{1}。

定理2 对于系统（1），当时，免疫未激活平衡点E_{\mathrm{1}}是局部渐近稳定的。

证明 将免疫未激活平衡点E_{\mathrm{1}}=(T_{\mathrm{1}},I_{\mathrm{1}},V_{\mathrm{1}},Z_{\mathrm{1}})代入（2）式中，经过简单计算可得E_{\mathrm{1}}处的特征方程为：

(r+d_{\mathrm{4}}-c{V}_{\mathrm{1}})(r^{\mathrm{3}}+B_{\mathrm{1}}{r}^{\mathrm{2}}+B_{\mathrm{2}}r+B_{\mathrm{3}})=\mathrm{0}，（4）

其中

\begin{array}{c}{B}_{\mathrm{1}}=d_{\mathrm{1}}+d_{\mathrm{3}}+d_{\mathrm{2}}\frac{R_{\mathrm{01}}}{R_{\mathrm{0}}}+d_{\mathrm{1}}(R_{\mathrm{0}}-\mathrm{1}),\\ B_{\mathrm{2}}=d_{\mathrm{1}}{d}_{\mathrm{3}}+d_{\mathrm{1}}{d}_{\mathrm{2}}\frac{R_{\mathrm{01}}}{R_{\mathrm{0}}}+d_{\mathrm{1}}\left(d_{\mathrm{2}}+d_{\mathrm{3}}+\frac{d_{\mathrm{2}}\alpha T_{\mathrm{1}}}{d_{\mathrm{2}}+\alpha T_{\mathrm{1}}}\frac{R_{\mathrm{01}}}{R_{\mathrm{0}}}\right)(R_{\mathrm{0}}-\mathrm{1}),\\ B_{\mathrm{3}}=d_{\mathrm{1}}{d}_{\mathrm{2}}{d}_{\mathrm{3}}(R_{\mathrm{0}}-\mathrm{1}),\end{array}

显然，特征方程（4）的一个特征值为c{V}_{\mathrm{1}}-d_{\mathrm{4}}=\frac{\lambda ck}{d_{\mathrm{2}}{d}_{\mathrm{3}}{R}_H}(R_H-\mathrm{1})，当即（引理2），该特征值为负。此外，当R_{\mathrm{0}}>\mathrm{1}时，B_{\mathrm{1}},B_{\mathrm{2}},B_{\mathrm{3}}>\mathrm{0}且

\begin{array}{c}{B}_{\mathrm{1}}{B}_{\mathrm{2}}-B_{\mathrm{3}}=d_{\mathrm{1}}{d}_{\mathrm{3}}(d_{\mathrm{1}}+d_{\mathrm{3}})+d_{\mathrm{1}}{d}_{\mathrm{2}}(d_{\mathrm{1}}+d_{\mathrm{3}})\frac{R_{\mathrm{01}}}{R_{\mathrm{0}}}+d_{\mathrm{1}}\frac{R_{\mathrm{01}}}{R_{\mathrm{0}}}\left(d_{\mathrm{1}}+\frac{\alpha T_{\mathrm{1}}{d}_{\mathrm{2}}}{\alpha T_{\mathrm{1}}+d_{\mathrm{2}}}\right)(R_{\mathrm{0}}-\mathrm{1})+\\ {d_{\mathrm{1}}}^{\mathrm{2}}(d_{\mathrm{2}}+d_{\mathrm{3}})(R_{\mathrm{0}}{-\mathrm{1})}^{\mathrm{2}}+d_{\mathrm{1}}(d_{\mathrm{1}}{d}_{\mathrm{2}}+\mathrm{2}{d}_{\mathrm{1}}{d}_{\mathrm{3}}+{d_{\mathrm{3}}}^{\mathrm{2}})(R_{\mathrm{0}}-\mathrm{1})>\mathrm{0}\mathrm{。}\end{array}

综上所述，当时，特征方程（4）的所有特征值均具有负实部。根据Routh-Hurwitz判据^［4］知，当时，免疫未激活平衡点E_{\mathrm{1}}是局部渐近稳定的。此外，当R_H>\mathrm{1}时，即R_{\mathrm{1}}>\mathrm{1}时（引理2），特征值c{V}_{\mathrm{1}}-d_{\mathrm{4}}>\mathrm{0}。所以，R_{\mathrm{1}}>\mathrm{1}时免疫未激活平衡点E_{\mathrm{1}}是不稳定的。

3 平衡点的全局稳定性

首先，考虑函数g(x)=x-\mathrm{1}-\mathrm{l}\mathrm{n}x，x>\mathrm{0}，注意到，g(x)\ge \mathrm{0}对任意x>\mathrm{0}成立，且g(x)=\mathrm{0}当且仅当x=\mathrm{1}。

定理3 当R_{\mathrm{0}}\le \mathrm{1}时，系统（1）的无病平衡点E_{\mathrm{0}}是全局渐近稳定的。

证明 由系统（1）的第一个方程知有\frac{\mathrm{d}T}{\mathrm{d}t}\le \lambda -d_{\mathrm{1}}T，通过比较原理知\underset{t\to \mathrm{\infty }}{\mathrm{l}\mathrm{i}\mathrm{m}}\mathrm{s}\mathrm{u}\mathrm{p}T(t)\le T_{\mathrm{0}}，其中T_{\mathrm{0}}=\frac{\lambda }{d_{\mathrm{1}}}。 假设T(t)\le T_{\mathrm{0}}，对任意的t\ge \mathrm{0}，显然有以下不等式组成立：

\left\{\begin{array}{l}\frac{\mathrm{d}I}{\mathrm{d}t}\le {\beta }_{\mathrm{1}}{T}_{\mathrm{0}}V+{\beta }_{\mathrm{2}}{T}_{\mathrm{0}}I-d_{\mathrm{2}}I,\\ \frac{\mathrm{d}V}{\mathrm{d}t}\le kI-d_{\mathrm{3}}V,\end{array}\right.（5）

对应的比较系统为：

\left\{\begin{array}{l}\frac{\mathrm{d}{I}_{\mathrm{1}}}{\mathrm{d}t}={\beta }_{\mathrm{1}}{T}_{\mathrm{0}}{V}_{\mathrm{1}}+{\beta }_{\mathrm{2}}{T}_{\mathrm{0}}{I}_{\mathrm{1}}-d_{\mathrm{2}}{I}_{\mathrm{1}},\\ \frac{\mathrm{d}{V}_{\mathrm{1}}}{\mathrm{d}t}=k{I}_{\mathrm{1}}-d_{\mathrm{3}}{V}_{\mathrm{1}},\end{array}\right.（6）

进一步得到矩阵

\mathit{A}=\left(\begin{array}{cc}{\beta }_{\mathrm{2}}{T}_{\mathrm{0}}-d_{\mathrm{2}}& {\beta }_{\mathrm{1}}{T}_{\mathrm{0}}\\ k& -d_{\mathrm{3}}\end{array}\right)，

从而矩阵 A 特征方程的根为

r=\frac{-\left(d_{\mathrm{2}}+d_{\mathrm{3}}-{\beta }_{\mathrm{1}}{T}_{\mathrm{0}}\right)\pm \sqrt[]{{\left(d_{\mathrm{2}}+d_{\mathrm{3}}-{\beta }_{\mathrm{1}}{T}_{\mathrm{0}}\right)}^{\mathrm{2}}+\mathrm{4}{d}_{\mathrm{2}}{d}_{\mathrm{3}}\left(R_{\mathrm{0}}-\mathrm{1}\right)}}{\mathrm{2}}，

若，那么。当t\to \mathrm{\infty }时，方程组（6）的解(I_{\mathrm{1}}(t),V_{\mathrm{1}}(t))\to (\mathrm{0,0})。由比较原理知，当t\ge \mathrm{0}时，(I(t),V(t))\le(I_{\mathrm{1}}(t),V_{\mathrm{1}}(t))。所以当t\to \mathrm{\infty }时，(I(t),V(t))\to (\mathrm{0,0})。进一步，由系统（1）的第一个方程和最后一个方程得如下极限方程：

\left\{\begin{array}{l}\frac{\mathrm{d}T}{\mathrm{d}t}=\lambda -d_{\mathrm{1}}T,\\ \frac{\mathrm{d}Z}{\mathrm{d}t}=-d_{\mathrm{4}}Z\mathrm{。}\end{array}\right.

由渐近自治半流理论^［6］得，t\to \mathrm{\infty }时，有Z(t)\to \mathrm{0}，T(t)\to T_{\mathrm{0}}。因此无病平衡点E_{\mathrm{0}}是全局渐近稳定的。

注3 在证明免疫未激活平衡点E_{\mathrm{1}}的全局渐近稳定性之前，我们首先给出如下假设：

（H1） \mathit{B}=\left(\begin{array}{cc}{B}_{\mathrm{11}}& -B_{\mathrm{12}}\\ -B_{\mathrm{12}}& B_{\mathrm{22}}\end{array}\right)为正定矩阵，其中

B_{\mathrm{11}}=\frac{\mathrm{1}}{T_m}(d_{\mathrm{1}}+\alpha I_{\mathrm{1}})+d_{\mathrm{1}}，B_{\mathrm{22}}=d_{\mathrm{2}}+\alpha T_{\mathrm{1}}，\mathrm{2}{B}_{\mathrm{12}}=d_{\mathrm{2}}+d_{\mathrm{1}}+\alpha T_{\mathrm{1}}+\alpha +\alpha I_m，

T_m，I_m是T(t)，I(t)的上界。

定理4 当且（H1）成立时，免疫未激活平衡点E_{\mathrm{1}}是全局渐近稳定的。

证明 定义Lyapunov函数U_{\mathrm{1}}(T,I,V,Z)：

U_{\mathrm{1}}=\left(T-T_{\mathrm{1}}-T_{\mathrm{1}}\mathrm{l}\mathrm{n}\frac{T}{T_{\mathrm{1}}}\right)+\left(I-I_{\mathrm{1}}-I_{\mathrm{1}}\mathrm{l}\mathrm{n}\frac{I}{I_{\mathrm{1}}}\right)+\frac{{\beta }_{\mathrm{1}}{T}_{\mathrm{1}}{V}_{\mathrm{1}}}{k{I}_{\mathrm{1}}}\left(V-V_{\mathrm{1}}-V_{\mathrm{1}}\mathrm{l}\mathrm{n}\frac{V}{V_{\mathrm{1}}}\right)+\frac{p{\beta }_{\mathrm{1}}{T}_{\mathrm{1}}{V}_{\mathrm{1}}}{ck{I}_{\mathrm{1}}}Z，

则

\begin{array}{c}\frac{\mathrm{d}{U}_{\mathrm{1}}}{\mathrm{d}t}=-\frac{(T-T_{\mathrm{1}}{)}^{\mathrm{2}}}{T}(d_{\mathrm{1}}+\alpha I_{\mathrm{1}})+\frac{p{\beta }_{\mathrm{1}}{T}_{\mathrm{1}}}{d_{\mathrm{3}}}\frac{\lambda Z}{d_{\mathrm{2}}{d}_{\mathrm{3}}{R}_H}(R_H-\mathrm{1})-{\beta }_{\mathrm{1}}{T}_{\mathrm{1}}{V}_{\mathrm{1}}\left[\frac{T_{\mathrm{1}}}{T}+\frac{I{V}_{\mathrm{1}}}{I_{\mathrm{1}}V}+\frac{TV{I}_{\mathrm{1}}}{T_{\mathrm{1}}{V}_{\mathrm{1}}I}-\mathrm{3}\right]-\\ {\beta }_{\mathrm{2}}{T}_{\mathrm{1}}{I}_{\mathrm{1}}\left[\frac{T_{\mathrm{1}}}{T}+\frac{T}{T_{\mathrm{1}}}-\mathrm{2}\right]-\alpha (T-T_{\mathrm{1}})(I-I_{\mathrm{1}})\mathrm{。}\end{array}

定义Lyapunov函数U_{\mathrm{2}}(T,I,V,Z)：

U_{\mathrm{2}}=\frac{\mathrm{1}}{\mathrm{2}}[(T-T_{\mathrm{1}})+(I-I_{\mathrm{1}}{)]}^{\mathrm{2}}，

则

\begin{array}{c}\frac{\mathrm{d}{U}_{\mathrm{2}}}{\mathrm{d}t}=[(T-T_{\mathrm{1}})+(I-I_{\mathrm{1}})]\left(\frac{\mathrm{d}T}{\mathrm{d}t}+\frac{\mathrm{d}I}{\mathrm{d}t}\right)=\\ -d_{\mathrm{1}}(T-T_{\mathrm{1}}{)}^{\mathrm{2}}-d_{\mathrm{1}}(T-T_{\mathrm{1}})(I-I_{\mathrm{1}})-d_{\mathrm{2}}(T-T_{\mathrm{1}})(I-I_{\mathrm{1}})-d_{\mathrm{2}}(I-I_{\mathrm{1}}{)}^{\mathrm{2}}-\\ \alpha T_{\mathrm{1}}(T-T_{\mathrm{1}})(I-I_{\mathrm{1}})-\alpha T_{\mathrm{1}}(I-I_{\mathrm{1}}{)}^{\mathrm{2}}-\alpha I(T-T_{\mathrm{1}}{)}^{\mathrm{2}}-\alpha I(T-T_{\mathrm{1}})(I-I_{\mathrm{1}})\mathrm{。}\end{array}

定义Lyapunov函数U(T,I,V,Z)：

U(T,I,V,Z)=U_{\mathrm{1}}(T,I,V,Z)+U_{\mathrm{2}}(T,I,V,Z)，

则

\begin{array}{c}\frac{\mathrm{d}{U}_{\mathrm{1}}}{\mathrm{d}t}+\frac{\mathrm{d}{U}_{\mathrm{2}}}{\mathrm{d}t}=-\frac{(T-T_{\mathrm{1}}{)}^{\mathrm{2}}}{T}(d_{\mathrm{1}}+\alpha I_{\mathrm{1}})+\frac{p{\beta }_{\mathrm{1}}{T}_{\mathrm{1}}}{d_{\mathrm{3}}}\frac{\lambda z}{d_{\mathrm{2}}{d}_{\mathrm{3}}{R}_H}(R_H-\mathrm{1})-{\beta }_{\mathrm{1}}{T}_{\mathrm{1}}{V}_{\mathrm{1}}\left[\frac{T_{\mathrm{1}}}{T}+\frac{I{V}_{\mathrm{1}}}{I_{\mathrm{1}}V}+\frac{TV{I}_{\mathrm{1}}}{T_{\mathrm{1}}{V}_{\mathrm{1}}I}-\mathrm{3}\right]-\\ {\beta }_{\mathrm{2}}{T}_{\mathrm{1}}{I}_{\mathrm{1}}\left[\frac{T_{\mathrm{1}}}{T}+\frac{T}{T_{\mathrm{1}}}-\mathrm{2}\right]-\alpha (T-T_{\mathrm{1}})(I-I_{\mathrm{1}})-d_{\mathrm{1}}(T-T_{\mathrm{1}}{)}^{\mathrm{2}}-d_{\mathrm{1}}(T-T_{\mathrm{1}})(I-I_{\mathrm{1}})-\\ d_{\mathrm{2}}(T-T_{\mathrm{1}})(I-I_{\mathrm{1}})-d_{\mathrm{2}}(I-I_{\mathrm{1}}{)}^{\mathrm{2}}-\alpha T_{\mathrm{1}}(T-T_{\mathrm{1}})(I-I_{\mathrm{1}})-\alpha T_{\mathrm{1}}(I-I_{\mathrm{1}}{)}^{\mathrm{2}}-\\ \alpha I(T-T_{\mathrm{1}}{)}^{\mathrm{2}}-\alpha I(T-T_{\mathrm{1}})(I-I_{\mathrm{1}})\le \\ -B_{\mathrm{11}}(T-T_{\mathrm{1}}{)}^{\mathrm{2}}-B_{\mathrm{22}}(I-I_{\mathrm{1}}{)}^{\mathrm{2}}+\mathrm{2}{B}_{\mathrm{12}}\left|(T-T_{\mathrm{1}})(I-I_{\mathrm{1}})\right|+\frac{p{\beta }_{\mathrm{1}}{T}_{\mathrm{1}}}{d_{\mathrm{3}}}\frac{\lambda z}{d_{\mathrm{2}}{d}_{\mathrm{3}}{R}_H}(R_H-\mathrm{1})-\\ {\beta }_{\mathrm{1}}{T}_{\mathrm{1}}{V}_{\mathrm{1}}\left[\frac{T_{\mathrm{1}}}{T}+\frac{I{V}_{\mathrm{1}}}{I_{\mathrm{1}}V}+\frac{TV{I}_{\mathrm{1}}}{T_{\mathrm{1}}{V}_{\mathrm{1}}I}-\mathrm{3}\right]-{\beta }_{\mathrm{2}}{T}_{\mathrm{1}}{I}_{\mathrm{1}}\left[\frac{T_{\mathrm{1}}}{T}+\frac{T}{T_{\mathrm{1}}}-\mathrm{2}\right]=\\ \left|T-T_{\mathrm{1}}\right|,\left|I-I_{\mathrm{1}}\right|)B(\left|T-T_{\mathrm{1}}\right|,\left|I-I_{\mathrm{1}}\right|)+\frac{p{\beta }_{\mathrm{1}}{T}_{\mathrm{1}}}{d_{\mathrm{3}}}\frac{\lambda z}{d_{\mathrm{2}}{d}_{\mathrm{3}}{R}_H}(R_H-\mathrm{1})-\\ (-{\beta }_{\mathrm{1}}{T}_{\mathrm{1}}{V}_{\mathrm{1}}[T_{\mathrm{1}}T+I{V}_{\mathrm{1}}{I}_{\mathrm{1}}V+TV{I}_{\mathrm{1}}{T}_{\mathrm{1}}{V}_{\mathrm{1}}I-\mathrm{3}]-{\beta }_{\mathrm{2}}{T}_{\mathrm{1}}{I}_{\mathrm{1}}[T_{\mathrm{1}}T+T{T}_{\mathrm{1}}-\mathrm{2}]\mathrm{。}\end{array}

显然，当，即（引理2）且（H1）成立时，由几何-平均不等式，可得成立。 当且仅当(T,I,V,Z)=(T_{\mathrm{1}},I_{\mathrm{1}},V_{\mathrm{1}},Z_{\mathrm{1}})时，\frac{\mathrm{d}U}{\mathrm{d}t}=\mathrm{0}严格成立。由LaSalle不变集原理^［21］知，当且（H1）成立时，免疫未激活平衡点E_{\mathrm{1}}是全局渐近稳定的。

4 数值模拟

在本节中，分别给定初始值为i{c}_{\mathrm{1}}≔(\mathrm{70,2},\mathrm{1,200})，i{c}_{\mathrm{2}}≔(\mathrm{50,5},\mathrm{5,300})和i{c}_{\mathrm{3}}≔(\mathrm{30,10,2},\mathrm{500})。绘制未感染细胞（图2（a）、图3（a））、感染细胞（图2（b）、图3（b））、病毒粒子（图2（c）、图3（c））和B细胞（图2（d）、图3（d））在不同参数值情况下随时间的演化行为。

数值算例1 取参数值为\lambda =\mathrm{0.75},{\beta }_{\mathrm{1}}=\mathrm{0.01},{\beta }_{\mathrm{2}}=\mathrm{0.01},d_{\mathrm{1}}=\mathrm{0.01},d_{\mathrm{2}}=\mathrm{1},d_{\mathrm{3}}=\mathrm{6},d_{\mathrm{4}}=\mathrm{0.3},k=\mathrm{1.9},p=\mathrm{0.006},\alpha =\mathrm{0}.\mathrm{01},c=\mathrm{0.1}。经过计算得，此时系统（1）的解轨线如图2所示。无病平衡点E_{\mathrm{0}}=(\mathrm{75,0},\mathrm{0,0})是全局渐近稳定的，满足定理3。

数值算例2 取参数值为\lambda =\mathrm{1},{\beta }_{\mathrm{1}}=\mathrm{1.06},{\beta }_{\mathrm{2}}=\mathrm{1},d_{\mathrm{1}}=\mathrm{1},d_{\mathrm{2}}=\mathrm{1.2},d_{\mathrm{3}}=\mathrm{1},d_{\mathrm{4}}=\mathrm{1},k=\mathrm{1},p=\mathrm{0.006},\alpha =\mathrm{0.01},c=\mathrm{1}，经过计算得R_{\mathrm{0}}=\mathrm{1.72}>\mathrm{1}和，同时B是正定的，此时系统（1）的解轨线如图3所示。免疫未激活平衡点E_{\mathrm{1}}=(\mathrm{0.583,0.346,0.346,0})是全局渐近稳定的，满足定理4。特别地，取初始值(\mathrm{30,10,2},\mathrm{500})，参数值\lambda =\mathrm{10},{\beta }_{\mathrm{1}}=\mathrm{0.01},{\beta }_{\mathrm{2}}=\mathrm{0.01},d_{\mathrm{1}}=\mathrm{0.01},d_{\mathrm{2}}=\mathrm{1},d_{\mathrm{3}}=\mathrm{6},d_{\mathrm{4}}=\mathrm{0.3}，k=\mathrm{2.9},p=\mathrm{0.006}，\alpha =\mathrm{0.001},c=\mathrm{0.1}。由图4可得，随着细胞凋亡率\alpha的增加，体液免疫激活再生数，系统平衡点从免疫激活平衡点转向免疫未激活平衡点。所以，细胞凋亡率的增加对病毒感染的控制有积极作用。

5 结论

本文建立了一类具有体液免疫、细胞间传播和凋亡率的病毒感染动力学模型，并对模型的动力学进行讨论。结果表明，时，无病平衡点E_{\mathrm{0}}是全局渐近稳定的。且（H1）成立时，免疫未激活平衡点E_{\mathrm{1}}是全局渐近稳定的。由于模型同时考虑体液免疫、细胞间传播和感染细胞诱导未感染细胞凋亡多种因素，增加了该病毒模型的耦合程度，因此对免疫激活平衡点的局部和全局稳定性的讨论变得较为困难，仍然需要我们继续努力。此外，本文提出的具有体液免疫和凋亡率的病毒模型分析具有非常重要的现实意义。主要体现为以下两个方面，一方面，病毒进入人体后，会激活免疫系统，阻止病毒进一步感染。另一方面，HIV基因表达产物可以产生有毒物质，直接或间接诱导未感染的CD4+T细胞凋亡^［17］，因此，凋亡率在病毒感染模型建模过程中是不可忽视的。例如，本文研究发现体液免疫和凋亡率有助于控制病毒感染，降低病毒对人体的伤害。另外，本文建立的模型及方法可类似地用于分析HIV、HCV等病毒感染动力学模型。

参考文献

[1] NOWAK M A, BANGHAM C R. Population Dynamics of Immune Responses to Persistent Viruses[J]. Science, 1996, 272(5258): 74-79. DOI: 10.1126/science.272.5258.74 .

[2] MOTHES W, SHERER N M, JIN J, et al. Virus Cell-to-cell Transmission[J]. J Virol, 2010, 84(17): 8360-8368. DOI: 10.1128/jvi.00443-10 .

[3] GUO T, QIU Z P, RONG L B. Analysis of an HIV Model with Immune Responses and Cell-to-cell Transmission[J]. Bull Malays Math Sci Soc, 2020, 43(1): 581-607. DOI: 10.1007/s40840-018-0699-5 .

[4] PAN S, CHAKRABARTY S P. Threshold Dynamics of HCV Model with Cell-to-cell Transmission and a Non-cytolytic Cure in the Presence of Humoral Immunity[J]. Commun Nonlinear Sci Numer Simul, 2018, 61: 180-197. DOI: 10.1016/j.cnsns.2018.02.010 .

[5] WANG Y, LU M M, JIANG D Q. Dynamic Behavior of a General Stochastic HIV Model with Virus-to-cell Infection, Cell-to-cell Transmission, Immune Response and Distributed Delays[J]. J Nonlinear Sci, 2023, 33(5): 97. DOI: 10.1007/s00332-023-09955-5 .

[6] LUO Y T, ZHANG L, ZHENG T T, et al. Analysis of a Diffusive Virus Infection Model with Humoral Immunity, Cell-to-cell Transmission and Nonlinear Incidence[J]. Phys A Stat Mech Appl, 2019, 535: 122415. DOI: 10.1016/j.physa.2019.122415 .

[7] TANG S T, TENG Z D, MIAO H. Global Dynamics of a Reaction-Diffusion Virus Infection Model with Humoral Immunity and Nonlinear Incidence[J]. Comput Math Appl, 2019, 78(3): 786-806. DOI: 10.1016/j.camwa.2019.03.004 .

[8] AHMAD NAIK P, ZU J, OWOLABI K M. Modeling the Mechanics of Viral Kinetics under Immune Control during Primary Infection of HIV-1 with Treatment in Fractional Order[J]. Phys A Stat Mech Appl, 2020, 545: 123816. DOI: 10.1016/j.physa.2019.123816 .

[9] HU Z X, ZHANG J J, WANG H, et al. Dynamics Analysis of a Delayed Viral Infection Model with Logistic Growth and Immune Impairment[J]. Appl Math Model, 2014, 38(2): 524-534. DOI: 10.1016/j.apm.2013.06.041 .

[10] WANG Y, ZHOU Y C, BRAUER F, et al. Viral Dynamics Model with CTL Immune Response Incorporating Antiretroviral Therapy[J]. J Math Biol, 2013, 67(4): 901-934. DOI: 10.1007/s00285-012-0580-3 .

[11] DEANS J A, COHEN S. Immunology of Malaria[J]. Annu Rev Microbiol, 1983, 37: 25-49. DOI: 10.1146/annurev.mi.37.100183.000325 .

[12] WANG S F, ZOU D Y. Global Stability of In-host Viral Models with Humoral Immunity and Intracellular Delays[J]. Appl Math Model, 2012, 36(3): 1313-1322. DOI: 10.1016/j.apm.2011.07.086 .

[13] MIAO H, LIU R, JIAO M Y. Global Dynamics of a Delayed Latent Virus Model with both Virus-to-cell and Cell-to-cell Transmissions and Humoral Immunity[J]. J Inequal Appl, 2021, 2021(1): 156. DOI: 10.1186/s13660-021-02691-y .

[14] LI M, GUO K, MA W B. Uniform Persistence and Global Attractivity in a Delayed Virus Dynamic Model with Apoptosis and both Virus-to-cell and Cell-to-cell Infections[J]. Mathematics, 2022, 10(6): 975. DOI: 10.3390/math10060975 .

[15] FAN R L, DONG Y P, HUANG G, et al. Apoptosis in Virus Infection Dynamics Models[J]. J Biol Dyn, 2014, 8(1): 20-41. DOI: 10.1080/17513758.2014.895433 .

[16] ZHANG T Q, WANG J L, SONG Y, et al. Dynamical Analysis of a Delayed HIV Virus Dynamic Model with Cell-to-cell Transmission and Apoptosis of Bystander Cells[J]. Complexity, 2020, 2020: 2313102. DOI: 10.1155/2020/2313102 .

[17] GUO S B, MA W B. Global Behavior of Delay Differential Equations Model of HIV Infection with Apoptosis[J]. Discrete Contin Dyn Syst Ser B, 2016, 21(1): 103-119. DOI: 10.3934/dcdsb.2016.21.103 .

[18] SELLIAH N, FINKEL T H. Biochemical Mechanisms of HIV Induced T Cell Apoptosis[J]. Cell Death Differ, 2001, 8(2): 127-136. DOI: 10.1038/sj.cdd.4400822 .

[19] MENG R, LUO Y T, ZHENG T T. Stability Analysis for a HIV Model with Cell-to-cell Transmission, Two Immune Responses and Induced Apoptosis[J]. AIMS Math, 2024, 9(6): 14786-14806. DOI: 10.3934/math.2024719 .

[20] 马怡婷, 张太雷, 邓金超, 等. 一类潜伏期传染且具有病毒变异的传染病模型[J]. 山西大学学报(自然科学版), 2023, 46(3): 574-584. DOI: 10.13451/j.sxu.ns.2022083 .

[21] MA Y T, ZHANG T L, DENG J C, et al. An Infectious Disease Model with Virus Variation and Infectivity in the Incubation Period[J]. J Shanxi Univ Nat Sci Ed, 2023, 46(3): 574-584. DOI: 10.13451/j.sxu.ns.2022083 .

[22] HALE J K, LUNEL S M V. Introduction to Functional Differential Equations[M]. Heidelberg: Springer Science & Business Media, 2013.

基金资助

新疆维吾尔自治区自然科学基金(2022D01C64)

国家自然科学基金(12201540)