基于双重感染机制和Logistic增长的基孔肯雅病毒动力学模型

A Chikungunya Virus Dynamic Model Incorporating Dual Infection Mechanism and Logistic Growth

摘要: 为深入理解基孔肯雅病毒的宿主体内感染机制，本文构建并分析一个包含双重感染途径的四维动力学模型。该模型整合病毒–细胞和细胞–细胞两种传播途径、易感细胞的Logistic增长以及B细胞介导的免疫应答。本文推导基本再生数$R_0$ ，并将其作为判定病毒传播趋势的关键阈值。通过构造合适的李雅普诺夫函数，严格证明当$R_0<1$ 时，无病平衡点是全局渐近稳定的；当$R_0>1$ 时，唯一的感染平衡点是全局渐近稳定的。数值模拟结果验证理论分析的正确性。参数敏感性分析亦揭示，宿主细胞的内在动力学是对$R_0$ 影响最大的因素。本研究为理解CHIKV的复杂感染动态及评估潜在干预策略提供有力的数学工具。

Abstract: To gain a deeper insight into the in-host infection mechanism of Chikungunya virus (CHIKV), we construct and analyze a four-dimensional dynamical model incorporating dual infection pathways. The model integrates two transmission routes, virus-to-cell and cell-to-cell, along with the Logistic growth of susceptible cells and B-cell mediated immune response. We derive the basic reproduction number, $R_0$ , and establish it as a key threshold for determining the trend of viral spread. By constructing suitable Lyapunov functions, we rigorously prove that the disease-free equilibrium is globally asymptotically stable when $R_0<1$ , and the unique endemic equilibrium is globally asymptotically stable when $R_0>1$ . Numerical simulations are performed, and the results validate our theoretical analysis. Furthermore, parameter sensitivity analysis reveals that the intrinsic dynamics of host cells have the most significant impact on $R_0$ . This study provides a robust mathematical tool for understanding the complex infection dynamics of CHIKV and for evaluating potential intervention strategies.

1. 引言

基孔肯雅热(chikungunya fever)是由基孔肯雅病毒(chikungunya virus, CHIKV)引起的急性传染病，主要通过伊蚊传播，临床表现以发热、皮疹及关节疼痛为特征。自1952年坦桑尼亚首次证实其流行、1956年成功分离病毒以来，该疾病长期在非洲和东南亚地区呈地方性流行态势。进入21世纪，印度洋地区爆发大规模疫情，2005~2007年期间，法属留尼汪岛发病数高达27万例，占当地人口近40%；印度报告疑似病例超过139万例，部分地区发病率突破45%，凸显其在流行区的高强度传播能力。值得关注的是，法国、美国等非流行国家不断出现输入性病例。2025年7月，广东顺德地区基孔肯雅热病例激增，累计接近三千例，并已扩散至省外区域，引发广泛关注。

临床层面，基孔肯雅热与登革热症状高度相似，极易导致误诊。尽管其病死率较低，但在蚊媒密度较高地区易引发大规模暴发流行，对区域医疗系统造成冲击。病原学研究显示，CHIKV感染早期触发B细胞免疫应答，产生针对病毒表面糖蛋白E1和E2表位的IgM抗体，第6天中和性IgM抗体的出现与病毒血症水平降低相关；感染中后期，IgG抗体成为中和病毒的主要成分，且CHIKV感染可诱导针对其他致关节炎甲病毒的交叉中和抗体及长期记忆B细胞，相关记忆B细胞在感染后可存续长达24年[1] [2]。

在病毒传播机制方面，CHIKV除通过经典的直接感染途径扩大感染细胞范围外，还可借助病毒诱导的细胞间长延伸结构，向邻近未感染细胞传播，该结构沿延伸方向存在新生病毒粒子的活跃出芽并在尖端富集，表明CHIKV至少存在两种细胞感染途径。其发病机制涉及病毒致病因子、宿主免疫反应及细胞因子通路的复杂相互作用，而建立动力学模型可通过抽象生物学过程，简化模拟病毒在人体内的动力学行为，为深入认识疾病发生发展规律、指导临床治疗策略制定提供重要参考[3] [4]。

近年来，已有一些数学模型被用于研究病毒动力学[5] [6]。然而，现有的CHIKV宿主体内模型往往侧重于单一感染途径，或简化了宿主细胞的增殖。特别是，较少有模型同时整合CHIKV的两种传播模式(病毒–细胞和细胞–细胞)、易感细胞的Logistic自限制增长以及B细胞介导的特异性免疫应答这三个关键生物学特征。

为了填补这一空白，并更准确地刻画CHIKV在人体内的感染动态，本文构建并分析了一个新的四维动力学模型。本文的主要工作包括：(1) 推导了模型的关键阈值——基本再生数$R_0$ ；(2) 运用Lyapunov函数法，严格证明了无病平衡点和感染平衡点的全局渐近稳定性；(3) 通过参数敏感性分析，发现宿主细胞的内在动力学(死亡率$\alpha$ 和Logistic增长系数$k$ )是对$R_0$ 影响最大的因素。本研究为深入理解CHIKV的感染机制和评估潜在干预措施提供一个更完善的数学框架。

2. 基于Logistic增长、两种感染途径的CHIKV感染动力学模型

2.1. 模型建立

在经典的病毒动力学模型基础上，本研究为更准确地描述CHIKV在人体内的感染过程，引入了几个关键的生物学机制。首先，考虑到宿主细胞的自我调节与资源限制，模型采用Logistic增长来描述易感细胞的增殖，这比恒定补充率的假设更符合实际情况。其次，鉴于CHIKV不仅通过游离病毒颗粒感染健康细胞，还能通过细胞间接触直接传播，模型中包含了两种感染途径。最后，针对CHIKV感染主要激发体液免疫的特点，模型纳入了由B细胞介导的特异性免疫应答过程。

综合以上因素，我们建立如下包含Logistic增长、两种感染途径和B细胞介导的免疫应答的CHIKV感染动力学模型：

$\begin{array}{l}\frac{\text{d}S}{\text{d}t}=\mu -\alpha S+kS\left(1-\frac{S}{S_{\max }}\right)-b_{1}SV-b_{2}SI,\\ \frac{\text{d}I}{\text{d}t}=b_{1}SV+b_{2}SI-\beta I,\\ \frac{\text{d}V}{\text{d}t}=mI-rV-\rho BV,\\ \frac{\text{d}B}{\text{d}t}=\eta +cBV-\delta B,\end{array}$ (1)

其中，$S$ 、$I$ 、$V$ 与$B$ 分别代表易感人体细胞、受感染人体细胞、游离病毒颗粒和免疫B细胞的密度。易感人体细胞以$\mu$ 的速度产生，以$\alpha$ 的速度死亡，$k$ 是易感人体细胞Logistic增长系数，$S_{max}$ 细胞最大容量。在这里我们考虑两种感染途径：第一，易感人体细胞可以通过与游离病毒颗粒接触而感染，记常数$b_1$ 为它们之间的接触率；第二，易感人体细胞可以通过与受感染人体细胞接触而感染，记常数$b_2$ 为它们之间的接触率。受感染的人体细胞、游离病毒颗粒和免疫B细胞分别以$\beta$ 、$r$ 和$\delta$ 的速度死亡。病毒的增殖假设以恒定的速率$m$ 。免疫B细胞以恒定的速率$\eta$ 产生。$\rho$ 是免疫B细胞与游离病毒的接触率。最后，免疫B细胞是由于接触游离病毒颗粒受刺激而增加，比例常数为$c$ 。

2.2. 模型平衡点

首先，模型(1)有一个无病平衡点$E_0\left(S_0,0,0,B_0\right)$ ，其中

$S_0=\frac{S_{\max }\left[\left(k-\alpha \right)+\sqrt{{\left(k-\alpha \right)}^2+\frac{4k\mu }{S_{\max }}}\right]}{2k},B_0=\frac{\eta }{\delta }.$

定义基本再生数

$R_0=\frac{b_1{S}_{0}m+b_2{S}_0\left(r+\rho B_0\right)}{\beta \left(r+\rho B_0\right)}.$

当$R_0>1$ 时，模型(1)存在感染平衡点$E^{*}\left(S^{*},I^{*},V^{*},B^{*}\right)$ ，其中

$B^{*}=\frac{\eta }{\delta -c{V}^{*}},I^{*}=\frac{V^{*}\left(r+\rho B^{*}\right)}{m},S^{*}=\frac{\beta I^{*}}{b_1{V}^{*}+b_2{I}^{*}}.$

$V^{*}$ 是四次方程$A{V}^{*}{}^4+B{V}^{*}{}^3+C{V}^{*}{}^2+D{V}^{*}+E=0$ 满足$S^{*},I^{*},B^{*}$ 均大于零的正解。系数的表达式在附录中给出。

2.3. 解的存在唯一性与非负性

定理2.3 在初始条件$I\left(0\right),V\left(0\right),S\left(0\right),B\left(0\right)\ge 0$ 下，系统(1)存在唯一非负解。

证明：设$X={\left(S,I,V,B\right)}^{\text{T}}$ ，$f\left(X\right)={\left(f_1,f_2,f_3,f_4\right)}^{\text{T}}$ ，其中$f_i$ 为系统(1)各方程右侧表达式。显然，$f$ 对$S,I,V,B$ 的偏导数在$S\ge 0,I\ge 0,V\ge 0,B\ge 0$ 区域连续，故$f$ 满足局部Lipschitz条件。于是在初始条件$S\left(0\right)\ge 0,I\left(0\right)\ge 0,V\left(0\right)\ge 0,B\left(0\right)\ge 0$ 下，系统(1)在$\left[0,+\infty \right)$ 上存在唯一解[7]。

我们采用反证法来证明解的非负性。假设解$\left(S\left(t\right),I\left(t\right),V\left(t\right),B\left(t\right)\right)$ 从非负初始条件出发后，在未来某个时刻会变为负值。则存在一个时间$t^{*}>0$ ，使得$S,I,V,B$ 中的至少一个分量在该时刻恰好等于0，并即将变为负值。根据$t^{*}$ 的定义，在$t\in \left[0,t^{*}\right]$ 的时间区间内，所有四个分量都必须保持非负。我们逐一考察在$t^{*}$ 时刻的各种可能情况：若$S\left(t^{*}\right)=0$ ，则$\frac{\text{d}S}{\text{d}t}=\mu \ge 0$ 。若$B\left(t^{*}\right)=0$ ，则$\frac{\text{d}B}{\text{d}t}=\eta \ge 0$ 。若$I\left(t^{*}\right)=0$ ，则$\frac{\text{d}I}{\text{d}t}=b_{1}S\left(t^{*}\right)V\left(t^{*}\right)$ ，因为此时$S\left(t^{*}\right)\ge 0$ 且$V\left(t^{*}\right)\ge 0$ ，故$\frac{\text{d}I}{\text{d}t}\ge 0$ 。同理，若$V\left(t^{*}\right)=0$ ，则$\frac{\text{d}V}{\text{d}t}=mI\left(t^{*}\right)$ ，因为$I\left(t^{*}\right)\ge 0$ ，故$\frac{\text{d}V}{\text{d}t}\ge 0$ 。在所有情况下，即将变为负值的变量的导数都是非负的，这意味着它不可能穿过0而变为负值。这与我们最初的假设(即存在这样一个$t^{*}$ )相矛盾。因此，解$\left(S\left(t\right),I\left(t\right),V\left(t\right),B\left(t\right)\right)$ 对于所有$t\ge 0$ 必须恒保持非负。

3. 平衡点的稳定性

3.1. 无病平衡点的局部稳定性

定理3.1 当$R_0<1$ 时，系统在无病平衡点$E_0\left(S_0,0,0,B_0\right)$ 处局部渐近稳定。

证明：在无病平衡点$\left(S_0,0,0,B_0\right)$ 处，雅可比矩阵$J$ 为：

$J=\left(\begin{array}{cccc}-\alpha +k\left(1-\frac{2S_0}{S_{\max }}\right)& -b_2{S}_0& -b_1{S}_0& 0\\ 0& b_2{S}_0-\beta & b_1{S}_0& 0\\ 0& m& -r-\rho B_0& 0\\ 0& 0& c{B}_0& -\delta \end{array}\right),$

其特征方程为$\left(\lambda +\delta \right)\left(\lambda +\frac{\mu }{S_0}+\frac{k{S}_0}{S_{\max }}\right)\left({\lambda }^2+a_1\lambda +a_0\right)=0$ ，可得特征值${\lambda }_1=-\delta <0$ ，${\lambda }_2=-\left(\frac{\mu }{S_0}+\frac{k{S}_0}{S_{\max }}\right)<0$ 。其余两个特征值是从方程

${\lambda }^2+a_1\lambda +a_0=0$ (2)

解得，其中$a_1=-b_2{S}_0+\beta +r+\rho B_0$ ，$a_0=-\left(b_2{S}_0-\beta \right)\left(r+\rho B_0\right)-m{b}_1{S}_0$ 。

由$R_0=\frac{b_1{S}_{0}m+b_2{S}_0\left(r+\rho B_0\right)}{\beta \left(r+\rho B_0\right)}<1$ 得$a_0=-\left[\left(b_2{S}_0-\beta \right)\left(r+\rho B_0\right)+m{b}_1{S}_0\right]>0$ ，$a_1=-b_2{S}_0+\beta +r+\rho B_0>0$ 。

利用韦达定理，可得方程(2)的解的实部均为负数。于是$J$ 在无病平衡点的特征值的实部均为负数。因此，当$R_0<1$ 时，系统(1)在无病平衡点$E_0\left(S_0,0,0,B_0\right)$ 处局部渐近稳定。

3.2. 无病平衡点的全局渐近稳定性

定理3.2 若$R_0<1$ ，系统(1)在无病平衡点$E_0$ 是全局渐近稳定的。

证明：定义函数

$\Phi \left(x\right)=x-1-\ln x.$

容易得到$\Phi \left(x\right)\ge 0$ 对于$x\ge 0$ ，$\Phi \left(x\right)=0$ 当且仅当$x=1.$

定义Lyapunov函数：

$L_1=\Phi \left(\frac{S}{S_0}\right)+I+\frac{b_1{S}_0}{r+\rho B_0}V+\frac{b_1{S}_0\rho }{c\left(r+\rho B_0\right)}\Phi \left(\frac{B}{B_0}\right).$

沿着原系统的轨线的全导数为：

$\begin{array}{c}\frac{\text{d}{L}_1}{\text{d}t}=\left(1-\frac{S_0}{S}\right)\left[\mu -\alpha S+kS\left(1-\frac{S}{S_{\max }}\right)\right]+\beta \left(R_0-1\right)I-\frac{b_1{S}_{0}p\delta }{c\left(r+\rho B_0\right)}\frac{{\left(B-B_0\right)}^2}{B}\\ =-\frac{{\left(S-S_0\right)}^2}{S}\left(\frac{\mu }{S_0}+\frac{kS}{S_{\max }}\right)+\beta \left(R_0-1\right)I-\frac{b_1{S}_{0}p\delta }{c\left(r+\rho B_0\right)}\frac{{\left(B-B_0\right)}^2}{B}\\ \le 0\end{array}$

此外，$\frac{\text{d}{L}_1}{\text{d}t}=0$ 当且仅当系统处于无病平衡点处。结合定理3.1，进而说明原系统在无病平衡点处具有全局渐近稳定性[8]。

3.3. 感染平衡点的局部稳定性

定理3.3 当$\alpha -k+\frac{k{S}^{*}}{S_{\max }}>0,R_0>1$ 时，系统在感染平衡点$E^{*}$ 具有局部渐近稳定性。

证明：假设系统在感染平衡点处的Jacobi矩阵的特征值$\lambda$ 的实部非负，对应的特征向量为${\left[w_S,w_I,w_V,w_B\right]}^{\text{T}}$ 。记$J\left(E^{*}\right)$ 为在感染平衡点处的Jacobi矩阵。根据

$J\left(E^{*}\right){\left[w_S,w_I,w_V,w_B\right]}^{\text{T}}=\lambda {\left[w_S,w_I,w_V,w_B\right]}^{\text{T}},$

$\left(\begin{array}{cccc}-\frac{\mu }{S^{*}}-\frac{k{S}^{*}}{S_{\max }}& -b_2{S}^{*}& -b_1{S}^{*}& 0\\ \frac{\beta I^{*}}{S^{*}}& b_2{S}^{*}-\beta & b_1{S}^{*}& 0\\ 0& m& -r-\rho B^{*}& -\rho V^{*}\\ 0& 0& c{B}^{*}& c{V}^{*}-\delta \end{array}\right)\left[\begin{array}{c}{w}_S\\ w_I\\ w_V\\ w_B\end{array}\right]=\lambda \left[\begin{array}{c}{w}_S\\ w_I\\ w_V\\ w_B\end{array}\right].$

展开得到

$\left(-\frac{\mu }{S^{*}}-\frac{k{S}^{*}}{S_{\max }}\right)w_S-b_2{S}^{*}{w}_I-b_1{S}^{*}{w}_V=\lambda w_S,$

$\frac{\beta I^{*}}{S^{*}}{w}_S+\left(b_2{S}^{*}-\beta \right)w_I+b_1{S}^{*}{w}_V=\lambda w_I,$

$m{w}_I-\left(r+\rho B^{*}\right)w_V-\rho V^{*}{w}_B=\lambda w_V,$

$c{B}^{*}{w}_V+\left(c{V}^{*}-\delta \right)w_B=\lambda w_B.$

整理得到

$1-\left[\frac{b_2{S}^{*}}{\beta +\lambda }+\frac{b_1{S}^{*}m}{\left(\beta +\lambda \right)G\left(\lambda \right)}\right]=\frac{b_1{V}^{*}+b_2{I}^{*}}{k-\alpha -\frac{2k{S}^{*}}{S_{\max }}-\lambda }.$ (3)

其中

$G\left(\lambda \right)=r+\rho B^{*}+\lambda +\frac{\rho c{B}^{*}{V}^{*}}{\lambda +\delta -c{V}^{*}},\left|G\left(\lambda \right)\right|\ge r+\rho B^{*}.$

对方程(3)两边同时取实部，得到

$ℛe\left\{1-\left[\frac{b_2{S}^{*}}{\beta +\lambda }+\frac{b_1{S}^{*}m}{\left(\beta +\lambda \right)G\left(\lambda \right)}\right]\right\}=ℛe\left(\frac{b_1{V}^{*}+b_2{I}^{*}}{k-\alpha -\frac{2k{S}^{*}}{S_{\max }}-\lambda }\right).$

$\begin{array}{c}左�=1-ℛe\left[\frac{b_2{S}^{*}}{\beta +\lambda }+\frac{b_1{S}^{*}m}{\left(\beta +\lambda \right)G\left(\lambda \right)}\right]\ge 1-\left|\frac{b_2{S}^{*}}{\beta +\lambda }+\frac{b_1{S}^{*}m}{\left(\beta +\lambda \right)G\left(\lambda \right)}\right|\\ \ge 1-\left[\frac{b_2{S}^{*}}{\left|\beta +\lambda \right|}+\frac{b_1{S}^{*}m}{\left|\beta +\lambda \right|\left|G\left(\lambda \right)\right|}\right]>1-\frac{b_1{S}^{*}m+b_2{S}^{*}\left(r+\rho B^{*}\right)}{\beta \left(r+\rho B^{*}\right)}=0.\end{array}$

$\begin{array}{c}右�=ℛe\left(\frac{b_1{V}^{*}+b_2{I}^{*}}{k-\alpha -\frac{2k{S}^{*}}{S_{\max }}-\lambda }\right)=\frac{\left(b_1{V}^{*}+b_2{I}^{*}\right)\left[k-\alpha -\frac{2k{S}^{*}}{S_{\max }}-ℛe\left(\lambda \right)\right]}{{\left[k-\alpha -\frac{2k{S}^{*}}{S_{\max }}-ℛe\left(\lambda \right)\right]}^2+{\left[Im\left(\lambda \right)\right]}^2}\\ <\frac{\left(b_1{V}^{*}+b_2{I}^{*}\right)\left[k-\alpha -\frac{k{S}^{*}}{S_{\max }}-ℛe\left(\lambda \right)\right]}{{\left[k-\alpha -\frac{2k{S}^{*}}{S_{\max }}-ℛe\left(\lambda \right)\right]}^2+{\left[Im\left(\lambda \right)\right]}^2}<0.\end{array}$

这里产生了矛盾，故系统在感染平衡点处的Jacobi矩阵的特征值$\lambda$ 的实部均为负数，从而定理成立。

3.4. 感染平衡点的全局渐近稳定性

定理3.4：当$R_0>1$ 且$\alpha -k+\frac{k{S}^{*}}{S_{\max }}>0$ 时，感染平衡点是全局渐近稳定的。

证明：

定义Lyapunov函数：

$L_2=S^{*}\Phi \left(\frac{S}{S^{*}}\right)+I^{*}\Phi \left(\frac{I}{I^{*}}\right)+\frac{b_1{S}^{*}{V}^{*}{}^2}{m{I}^{*}}\Phi \left(\frac{V}{V^{*}}\right)+\frac{b_1{B}^{*}{S}^{*}{V}^{*}\rho }{cm{I}^{*}}\Phi \left(\frac{B}{B^{*}}\right).$

其中函数Φ定义于定理3.2。

沿着原系统的轨线的全导数为：

$\frac{\text{d}{L}_2}{\text{d}t}=\left(1-\frac{S^{*}}{S}\right)\frac{\text{d}S}{\text{d}t}+\left(1-\frac{I^{*}}{I}\right)\frac{\text{d}I}{\text{d}t}+\frac{b_1{S}^{*}{V}^{*}}{m{I}^{*}}\left(1-\frac{V^{*}}{V}\right)\frac{\text{d}V}{\text{d}t}+\frac{b_1{S}^{*}{V}^{*}\rho }{cm{I}^{*}}\left(1-\frac{B^{*}}{B}\right)\frac{\text{d}B}{\text{d}t}.$

由于$E^{*}$ 是原系统的感染平衡点，所以有

$\mu =b_1{S}^{*}{V}^{*}+b_2{S}^{*}{I}^{*}+\alpha S^{*}-k{S}^{*}\left(1-\frac{S^{*}}{S_{\max }}\right),$

$\beta =\frac{b_1{S}^{*}{V}^{*}}{I^{*}}+b_2{S}^{*},r=\frac{m{I}^{*}}{V^{*}}-\rho B^{*},\delta =\frac{\eta }{B^{*}}+c{V}^{*}.$

因此，有

$\begin{array}{c}\frac{\text{d}{L}_2}{\text{d}t}=\left(1-\frac{S^{*}}{S}\right)\left(b_1{S}^{*}{V}^{*}+b_2{S}^{*}{I}^{*}+\alpha S^{*}-k{S}^{*}\left(1-\frac{S^{*}}{S_{\max }}\right)-\alpha S+kS\left(1-\frac{S}{S_{\max }}\right)-b_{1}SV-b_{2}SI\right)\\ +\left(1-\frac{I^{*}}{I}\right)\left(b_{1}SV+b_{2}SI-\frac{b_1{S}^{*}{V}^{*}I}{I^{*}}+b_2{S}^{*}I\right)+\frac{b_1{S}^{*}{V}^{*}}{m{I}^{*}}\left(1-\frac{V^{*}}{V}\right)\left(mI-\frac{m{I}^{*}B}{V^{*}}-\rho B^{*}V-\rho BV\right)\\ +\frac{b_1{S}^{*}{V}^{*}\rho }{cm{I}^{*}}\left(1-\frac{B^{*}}{B}\right)\left(\eta +cBV-\frac{\eta B}{B^{*}}+c{V}^{*}B\right)\\ =-\frac{{\left(S-S^{*}\right)}^2}{S}\left[\alpha -k+\frac{k\left(S+S^{*}\right)}{S_{\max }}\right]+b_1{S}^{*}{V}^{*}+b_2{S}^{*}{I}^{*}-b_{1}SV-b_{2}SI-\frac{b_1{S}^{*}{}^2{V}^{*}}{S}-\frac{b_2{S}^{*}{}^2{I}^{*}}{S}\\ +b_{1}SV+b_{2}SI-\frac{b_1{S}^{*}{V}^{*}I}{I^{*}}+b_2{S}^{*}I-\frac{b_{1}SV{I}^{*}}{I}-b_{2}S{I}^{*}+b_1{S}^{*}{V}^{*}-b_2{S}^{*}{I}^{*}\\ +\frac{b_1{S}^{*}{V}^{*}}{m{I}^{*}}\left(mI-\frac{Vm{I}^{*}}{V^{*}}+\rho B^{*}V-\rho BV-\rho B^{*}{V}^{*}+\rho B{V}^{*}-\frac{mI{V}^{*}}{V}+m{I}^{*}\right)\\ +\frac{b_1{S}^{*}{V}^{*}\rho }{cm{I}^{*}}\left(1-\frac{B^{*}}{B}\right)\left[\eta \left(1-\frac{B}{B^{*}}\right)+cBV-CB{V}^{*}\right]\\ =-\frac{{\left(S-S^{*}\right)}^2}{S}\left[\alpha -k+\frac{k\left(S+S^{*}\right)}{S_{\max }}\right]-b_2{S}^{*}{I}^{*}\left(\frac{S}{S^{*}}+\frac{S^{*}}{S}-2\right)-\frac{b_1{S}^{*}{V}^{*}\rho \eta }{cmB{B}^{*}{I}^{*}}{\left(B-B^{*}\right)}^2\\ -b_1{S}^{*}{V}^{*}\left(\frac{S^{*}}{S}+\frac{V^{*}I}{V{I}^{*}}+\frac{I^{*}SV}{I{S}^{*}{V}^{*}}-3\right)\\ \le -\frac{{\left(S-S^{*}\right)}^2}{S}\left(\alpha -k+\frac{k{S}^{*}}{S_{\max }}\right)-b_2{S}^{*}{I}^{*}\left(\frac{S}{S^{*}}+\frac{S^{*}}{S}-2\right)-\frac{b_1{S}^{*}{V}^{*}\rho \eta }{cmB{B}^{*}{I}^{*}}{\left(B-B^{*}\right)}^2\\ -b_1{S}^{*}{V}^{*}\left(\frac{S^{*}}{S}+\frac{V^{*}I}{V{I}^{*}}+\frac{I^{*}SV}{I{S}^{*}{V}^{*}}-3\right)\\ \le 0.\end{array}$

此外，$\frac{\text{d}{L}_2}{\text{d}t}=0$ 当且仅当系统处于感染平衡点处。结合定理3.3，进而说明原系统在感染平衡点处具有全局渐近稳定性[8]。

4. 数值模拟

为验证前述理论分析的正确性，本节对模型(1)进行数值模拟，并对模型参数进行敏感性分析。所有模拟均基于MATLAB (R2023a版)平台。模型参数的含义、基准值及取值范围主要依据参考文献[9] [10]，详见表1。

Table 1. Model parameters, descriptions, units, and values

表1. 模型参数、含义、单位与取值

参数	含义	单位	取值范围	Data1	Data2
$\mu$	易感细胞产生率	Day^−^1	0.1~10	1	1
$k$	易感细胞的Logistic增长系数	—	0.01~1	0.01	0.01
$S_{\max }$	易感细胞最大容纳量	—	1000~3000	1000	1000
$b_1$	易感细胞与游离病毒接触感染的接触率	Day^−^1	0.001~1	0.001	0.001
$b_2$	易感细胞与受感染细胞接触感染的接触率	Day^−^1	0.0001~0.01	0.0001	0.0001
$\beta$	受感染人体细胞死亡率	Day^−^1	0.1~1.5	1.02	0.1
$m$	病毒增殖率	Day^−^1	0.1~10	0.5	0.5
$r$	游离病毒颗粒死亡率	Day^−^1	1~5	1	1
$\rho$	免疫细胞清除病毒率	Day^−^1	0.01~0.1	0.01	0.01
$\eta$	免疫细胞产生率	Day^−^1	0.1~1	0.1	0.1
$c$	免疫细胞受病毒刺激增加的比例常数	—	0.001~0.1	0.001	0.001
$\delta$	免疫细胞死亡率	Day^−^1	0.1~1	0.1	0.1

4.1. 无病平衡点的稳定性

为验证定理3.2，选取表1中的参数组Data1。经计算，无病平衡点为$E_0\left(316.2248,0,0,1\right)$ ，基本再生数$R_0=0.1845<1$ 。根据理论分析，此时$E_0$ 应为全局渐近稳定。我们选取四组不同的初始条件(见表2)进行模拟。

Table 2. Initial conditions for the stability simulation of the disease-free equilibrium $E_0$

表2. 无病平衡点$E_0$ 稳定性模拟的初始条件

变量初值	组1	组2	组3	组4
S(0)	50	500	150	10
I(0)	0	10	3	0.1
V(0)	0.02	5	1.5	0.01
B(0)	0.9	5	1.5	0.1

图1展示了当$R_0<1$ 时四组初始条件下系统的时间序列图。

图2展示了(S, I, V)三维相空间中的轨迹(B(0) = 10)。

模拟结果显示，尽管初始条件差异显著，所有轨迹最终均收敛至无病平衡点 $E_0\left(316.2248,0,0,1\right)$ 。即$\underset{t\to +\infty }{\lim }S\left(t\right)=316.2248$ ，$\underset{t\to +\infty }{\lim }I\left(t\right)=0$ ，$\underset{t\to +\infty }{\lim }V\left(t\right)=0$ ，$\underset{t\to +\infty }{\lim }B\left(t\right)=1$ ，此结果与理论计算的无病平衡点数据完全吻合，在数值上验证了当$R_0<1$ 时，无病平衡点$E_0$ 是全局渐近稳定的。

Figure 1. Time series plots of the system with four different sets of initial conditions for $R_0<1$

图1. $R_0<1$ 时四组初始条件下系统的时间序列图

Figure 2. 3D phase portrait of the disease-free equilibrium $E_0$ and its projections on the S-I, S-V, and I-V planes

图2. 无病平衡点$E_0$ 的三维相图及在S-I，S-V，I-V平面的投影

4.2. 感染平衡点的稳定性

为验证定理3.4，选取表1中的参数组Data2。经计算，基本再生数$R_0=1.8817>1$ ，且满足定理3.4的条件。此时，系统存在唯一的感染平衡点 $E^{*}\left(168.1042,7.1741,3.5502,1.0368\right)$ 。理论上，该平衡点是全局渐近稳定的。我们同样选取四组不同的初始条件(见表3)进行模拟。

Table 3. Initial conditions for the stability simulation of the endemic equilibrium $E^{*}$

表3. 感染平衡点$E^{*}$ 稳定性模拟的初始条件

变量初值	组1	组2	组3	组4
S(0)	50	300	150	10
I(0)	0	10	3	0
V(0)	0.002	50	1.5	0.001
B(0)	0.9	50	1.5	0.1

图3展示了当$R_0>1$ 时四组初始条件下系统的时间序列图。

Figure 3. Time series plots of the system with four different sets of initial conditions for $R_0>1$

图3. $R_0>1$ 时四组初始条件下系统的时间序列图

图4展示了(S, I, V)三维相空间中的轨迹(B(0) = 10)。

模拟结果显示，当$R_0>1$ 时，无论初始状态如何，系统最终均收敛至感染平衡点$E^{*}$ 。即$\underset{t\to +\infty }{\lim }S\left(t\right)=168.1042$ ，$\underset{t\to +\infty }{\lim }I\left(t\right)=7.1741$ ，$\underset{t\to +\infty }{\lim }V\left(t\right)=3.5502$ ，$\underset{t\to +\infty }{\lim }B\left(t\right)=1.0368$ 。相图(图4)也直观地表明$E^{*}$ 是一个稳定的焦点。此结果与理论计算的感染平衡点数据吻合，在数值上验证了定理3.4的结论。

Figure 4. 3D phase portrait of the endemic equilibrium $E^{*}$ and its projections on the S-I, S-V, and I-V planes

图4. 感染平衡点$E^{*}$ 的三维相图及在S-I，S-V，I-V平面的投影

4.3. 参数敏感性分析

为识别对病毒感染建立影响最关键的参数，本研究对基本再生数$R_0$ 进行了全局敏感性分析。我们采用拉丁超立方抽样(LHS)与偏秩相关系数(PRCC)的方法。

分析中，所有12个模型参数$\mu ,\alpha ,k,S_{\max },b_1,b_2,\beta ,m,r,\rho ,\eta ,\delta$ 以及其分布范围均以表1中Data2的参数值和取值范围为基准。LHS被用于生成$N=1000$ 组参数样本，并计算相应的PRCC值。

图5中的分析结果清晰地表明，宿主易感细胞的内在动力学是决定$R_0$ 的最关键因素。易感细胞的自然死亡率$\alpha$ (PRCC = −0.85)表现为最强的负相关参数，而其Logistic增长系数$k$ (PRCC = 0.78)则为最强的正相关参数。这一发现揭示了宿主细胞群的补充与损耗平衡是决定病毒能否成功建立感染的首要条件。

病毒生命周期中的关键参数亦表现出显著影响。病毒–细胞感染率$b_1$ (PRCC = 0.60)和病毒增殖率$m$ (PRCC = 0.60)均是$R_0$ 的强力促进因素。相反，受感染细胞的死亡率$\beta$ (PRCC = −0.50)及游离病毒的清除率$r$ (PRCC = −0.35)均起到了显著的抑制作用。值得注意的是，模型所包含的两种感染途径贡献极不均衡：与$b_1$ 相比，细胞–细胞接触感染率$b_2$ (PRCC = 0.01)对$R_0$ 的影响几乎可以忽略不计。

此外，与B细胞免疫应答相关的参数，包括$\rho$ 、$\eta$ 和$\delta$ ，其PRCC值均接近于零。这说明在$R_0$ 所衡量的感染阈值阶段，个体的基础免疫状态对阻止病毒入侵的影响甚微。综上所述，$R_0$ 阈值主要由易感细胞的内在动力学$\left(\alpha ,k\right)$ 和病毒–细胞的感染途径$\left(b_1,m,\beta \right)$ 共同决定。

Figure 5. Sensitivity analysis of $R_0$ to model parameters

图5. $R_0$ 对模型参数的敏感性分析

5. 总结

本研究成功构建并分析了一个新的基孔肯雅病毒(CHIKV)宿主体内动力学模型。该模型创新性地整合了病毒的两种感染途径(病毒–细胞和细胞–细胞)、易感细胞的Logistic增长以及B细胞介导的免疫应答，从而更贴近真实的生物学过程。在理论分析层面，本文推导了模型的关键阈值——基本再生数$R_0$ 。通过构造李雅普诺夫函数，我们严格证明了$R_0=1$ 是决定病毒在宿主体内消亡或持续的唯一拐点：当$R_0<1$ 时，无病平衡点是全局渐近稳定的，意味着病毒将被清除；$R_0>1$ ，唯一的感染平衡点在特定条件下是全局渐近稳定的，意味着感染将持续存在。

在数值分析层面，本文进行了两方面的工作。首先，通过设置多组不同的初始条件进行数值模拟，我们验证了理论分析的正确性。模拟结果显示，当$R_0<1$ 时，所有轨线均收敛至无病平衡点；当$R_0>1$ 时，所有轨线均收敛至感染平衡点，这与我们的理论证明完全吻合。其次，我们对$R_0$ 进行了参数敏感性分析，以识别关键影响因素。

敏感性分析的结果揭示了几个关键的生物学见解。分析表明，对$R_0$ 影响最显著的参数并非来自病毒本身，而是宿主易感细胞的内在动力学特性。具体而言，易感细胞的自然死亡率$\alpha$ 呈现最强的负相关性，而其Logistic增长系数$k$ 则呈现最强的正相关性。此外，病毒–细胞接触率$b_1$ 、病毒增殖速率$m$ 以及受感染细胞死亡率$\beta$ 也显示出对$R_0$ 的显著影响。

本研究的发现具有重要的理论和潜在临床意义。从理论上讲，本模型提供了一个更完善的数学框架，用于理解CHIKV复杂的宿主体内动态。从临床应用上看，敏感性分析的发现提示，除了传统的抗病毒疗法(如抑制$m$ 或$b_1$ )，靶向宿主的治疗策略(Host-Directed Therapy)可能是一种极具潜力的干预策略，例如开发保护易感细胞、调节其内在动力学(如针对$k$ 或$\alpha$ )的药物。

尽管本研究构建的模型在整合双重感染途径和免疫应答方面取得了进展，但我们必须承认模型仍存在一定的局限性，这些简化假设可能会对结果产生细微影响。首先，为了保持数学分析的可处理性，我们简化了B细胞的动力学过程(方程中假设B细胞随病毒刺激即时增长)，未详细刻画抗体产生的滞后效应、抗体类型的转换或免疫记忆的形成，这可能在一定程度上低估了病毒清除的复杂性。其次，作为常微分方程(ODE)模型，本研究基于“均匀混合”假设，忽略了细胞间直接传播的空间局限性；在实际生物组织中，细胞间接触通常发生在局部，忽略空间结构可能会高估病毒在感染早期的扩散速度。最后，模型未包含潜伏期，即忽略了细胞从被感染到具备传染性之间的时间间隔，这可能会影响对感染峰值时间和峰值浓度的瞬态预测。未来的研究可以考虑引入时滞微分方程、偏微分方程或随机模型，以更精确地捕捉病毒感染的时空特征和随机性。

综上所述，尽管存在上述简化假设，本研究构建的分析框架仍然成功捕捉了CHIKV感染过程中关键的非线性动力学特征。本文不仅在数学上严格证明了基本再生数$R_0$ 决定病毒持久性或消亡的精确阈值，更在生物学意义上揭示了宿主细胞内在动力学对感染建立的决定性作用。这些发现不仅填补了现有模型在双重感染与细胞增长机制整合方面的空白，也为未来制定更精准的、针对宿主的抗病毒干预策略提供了可靠的理论依据和定量参考。

基金项目

山东省大学生创新创业训练计划资助项目(S202510425069)。

NOTES

^*通讯作者。