溶瘤病毒-CAR-T细胞联合治疗的随机脉冲动力学模型及其平稳分布

Stochastic Impulsive Dynamical Model and Stationary Distribution of Oncolytic Virus-CAR-T Cell Combination Therapy

摘要: 近年来，溶瘤病毒与CAR-T细胞联合治疗在癌症免疫治疗领域展现出广阔的应用前景。本研究构建了具有脉冲输注和随机扰动的肿瘤免疫动力学模型，探讨CAR-T细胞在溶瘤病毒联合治疗中的作用机制。通过构造辅助函数，建立了系统存在唯一遍历平稳分布的条件，并推导出CAR-T细胞脉冲输注次数的约束条件。结果表明，输注频率过高可能降低疗效甚至产生负面影响，强调合理控制输注频率的重要性。该研究为优化CAR-T细胞治疗方案提供数学依据，有助于提高联合治疗的临床效果。

Abstract: In recent years, the combination of oncolytic viruses (OVs) and chimeric antigen receptor T (CAR-T) cells has shown a broad application prospect in cancer immunotherapy. This study develops a tumor immunodynamic model incorporating pulse infusion and stochastic perturbations to investigate the mechanistic role of CAR-T cells in OV-based combination therapy. By constructing an auxiliary function, we establish conditions ensuring the existence of a unique ergodic stationary distribution and derive constraints on the frequency of CAR-T cell pulse infusions. The results indicate that excessively frequent infusions may diminish therapeutic efficacy or even induce adverse effects, highlighting the importance of optimizing infusion schedules. This study provides a mathematical foundation for refining CAR-T cell treatment strategies, contributing to the enhancement of clinical outcomes in combination therapy.

1. 引言

嵌合抗原受体T (CAR-T)细胞疗法是一种新兴的过继细胞治疗策略，通过基因工程改造患者T细胞，使其表达嵌合抗原受体CAR，从而特异性识别并清除肿瘤细胞[1]。尽管CAR-T疗法在血液肿瘤中取得突破，但在实体瘤治疗中仍面临挑战，如肿瘤微环境的免疫抑制作用。

溶瘤病毒(OVs)是一类能够选择性感染并裂解肿瘤细胞的病毒，同时能激活宿主免疫系统，促进抗肿瘤免疫反应[2]-[5]。研究表明，OVs可以克服CAR-T疗法的部分局限性，例如通过感染肿瘤细胞，使其表达CAR-T可识别的抗原，并改善肿瘤微环境以增强CAR-T细胞的疗效[6]。此外，某些OVs还能促进CAR-T细胞募集，以提高肿瘤杀伤能力[7]。研究还发现，CD19-CAR-T细胞可促使死亡的肿瘤细胞释放更多完整的病毒颗粒，这进一步增强了病毒在肿瘤中的扩散，并形成CAR-T细胞和OVs之间的协同作用，以提高抗肿瘤疗效[8]。

近年来，数学建模已成为研究CAR-T细胞和OVs相互作用的重要工具。2016年Walker等人[9]提出了CAR-T细胞与OVs协同作用的数学模型，随后Mahasa等人[10]简化了该模型，并探讨了单剂量和多剂量给药方式对疗效的影响。然而，生物系统常受到环境因素的影响，表现出随机性和不确定性[11] [12]。此外，临床上治疗的实施通常是脉冲式注射药物，我们考虑CAR-T细胞的输注采用脉冲式策略，以维持其在肿瘤微环境中的杀伤效能并缓解T细胞衰竭，多次脉冲输注策略不仅可以增强CAR-T细胞治疗的持久性，还可能缓解T细胞衰竭问题，从而对肿瘤细胞施加持续免疫压力[13] [14]。

基于此，我们采用脉冲随机微分方程构建了一个动力学模型，描述未感染肿瘤细胞、受感染肿瘤细胞、溶瘤病毒及CAR-T细胞之间的相互作用，如下所示：

$\{\begin{array}{l}\text{d}U\left(t\right)=\left[r\left(t\right)U\left(1-\frac{U+I}{K\left(t\right)}\right)-\beta \left(t\right)UV-\alpha \left(t\right)CU-{\mu }_u\left(t\right)U\right]\text{d}t+{\sigma }_1\left(t\right)Ud{B}_1\left(t\right),t\ne t_k,k\in N\\ \text{d}I\left(t\right)=\left[\beta \left(t\right)UV-\delta \left(t\right)I-\alpha \left(t\right)CI\right]\text{d}t+{\sigma }_2\left(t\right)I\text{d}{B}_2\left(t\right),t\ne t_k,k\in N\\ \text{d}V\left(t\right)=\left[\eta \left(t\right)\left(1+\omega \left(t\right)C\right)I-\theta \left(t\right)CV-{\mu }_v\left(t\right)V\right]\text{d}t+{\sigma }_3\left(t\right)V\text{d}{B}_{\text{3}}\left(t\right),t\ne t_k,k\in N\\ \text{d}C\left(t\right)=\left[\lambda \left(t\right)+\xi \left(t\right)CV+\rho \left(t\right)\left(U+I\right)C-{\mu }_c\left(t\right)C\right]\text{d}t+{\sigma }_4\left(t\right)C\text{d}{B}_{\text{4}}\left(t\right),t\ne t_k,k\in N\\ U\left(t_k^{+}\right)=\left(1+R_{1k}\right)U\left(t\right),t=t_k,k\in N\\ I\left(t_k^{+}\right)=\left(1+R_{2k}\right)I\left(t\right),t=t_k,k\in N\end{array}$ (1)

其中，$r$ 为未感染肿瘤细胞的增殖率，$K$ 表示肿瘤环境承载量，$\beta$ 描述游离病毒对未感染肿瘤细胞的感染率，$\alpha$ 表示CAR-T细胞对肿瘤细胞的杀伤率，$\delta$ 表示被感染肿瘤细胞的裂解率，$\eta$ 表示裂解细胞释放的病毒载量，$\omega$ 表示CAR-T细胞诱导的病毒释放增强率，$\theta$ 表示CAR-T细胞对病毒的清除速率，$\lambda$ 表示CAR-T细胞疗法，$\xi$ 表示病毒诱导的CAR-T细胞募集速率，$\rho$ 表示CAR-T细胞的增殖率，${\mu }_u$ 、${\mu }_v$ 以及${\mu }_c$ 分别表示未感染肿瘤死亡率、病毒清除率、CAR-T细胞清除率及其与肿瘤细胞相互作用的失活率，$R_{1k}<0$ 是指在$t=t_k$ 时注射CAR-T细胞导致的未感染的肿瘤细胞的死亡，$R_{2k}<0$ 是指在$t=t_k$ 时注射CAR-T细胞导致的受感染的肿瘤细胞的死亡，当$R_{ik}\left(i=1,2\right)>0$ 是指肿瘤细胞从其他组织转移到所考虑的组织。我们必须指出，所有参数都是正的。这里$B_i\left(t\right)\left(i=1,2,3,4\right)$ 表示在完备概率空间$\left(\Omega ,ℱ,{\left\{{ℱ}_t\right\}}_{t\ge 0},P\right)$ 上的标准布朗运动，其中$\Omega$ 是样本空间，$P$ 是概率测度，${\left\{{ℱ}_t\right\}}_{t\ge 0}$ 是满足通常条件的滤波，${\sigma }_i\left(i=1,2,3,4\right)$ 表示白噪声的强度。

本文的组织结构如下：第2部分证明了唯一遍历平稳分布存在的充分条件；第3部分提供数值实例验证结论的正确性以及给出优化肿瘤治疗的策略；最后，我们在第4部分进行了总结和展望。

2. 唯一遍历平稳分布

我们主要研究模型(1)的遍历平稳分布。首先，我们陈述其全局正解的存在性，这是随机动力系统研究中的关键步骤。它确保了模型的生物学合理性、数学可解释性，并为进一步研究系统的稳定性、遍历性和平稳分布提供理论基础。

定理2.1. 全局正解的存在性

对于任意的初值$\left(U\left(0\right),I\left(0\right),V\left(0\right),C\left(0\right)\right)\in R_{+}^4$ ，对所有$t\ge 0$ ，模型(1)以概率为1的存在唯一的全局正解$\left(U\left(t\right),I\left(t\right),V\left(t\right),C\left(t\right)\right)$ 。

定理2.1的证明是标准推导，因此我们省略了详细过程，读者可以参考文献[15]。

定理2.2. 唯一遍历平稳分布

在本节中，将使用文献[16]提出的方法研究模型(1)唯一遍历平稳分布的存在性。首先，考虑以下不带脉冲的随机微分方程：

$\{\begin{array}{l}\text{d}{x}_1\left(t\right)=\left[r\left(t\right)x_1\left(1-\frac{x_1{D}_1\left(t\right)+x_2{D}_2\left(t\right)}{K\left(t\right)}\right)-\beta \left(t\right)x_1{x}_3-\alpha \left(t\right)x_1{x}_4-{\mu }_u\left(t\right)x_1\right]\text{d}t+{\sigma }_1\left(t\right)x_1\text{d}{B}_1\left(t\right)\\ \text{d}{x}_2\left(t\right)=\left[\frac{\beta \left(t\right)D_1\left(t\right)}{D_{\text{2}}\left(t\right)}{x}_1{x}_{\text{3}}-\delta \left(t\right)x_2-\alpha \left(t\right)x_2{x}_{\text{4}}\right]\text{d}t+{\sigma }_2\left(t\right)x_2\text{d}{B}_2\left(t\right)\\ \text{d}{x}_{\text{3}}\left(t\right)=\left[\eta \left(t\right)D_{\text{2}}\left(t\right)\left(1+\omega \left(t\right)x_{\text{4}}\right)x_{\text{2}}-\theta \left(t\right)x_{\text{3}}{x}_{\text{4}}-{\mu }_v\left(t\right)x_{\text{3}}\right]\text{d}t+{\sigma }_3\left(t\right)x_{\text{3}}\text{d}{B}_{\text{3}}\left(t\right)\\ \text{d}{x}_{\text{4}}\left(t\right)=\left[\lambda \left(t\right)+\xi \left(t\right)x_{\text{3}}{x}_{\text{4}}+\rho \left(t\right)\left(x_1{D}_1\left(t\right)+x_2{D}_2\left(t\right)\right)x_{\text{4}}-{\mu }_c\left(t\right)x_{\text{4}}\right]\text{d}t+{\sigma }_4(t)x_{\text{4}}\text{d}{B}_{\text{4}}\left(t\right)\end{array}$ (2)

其中$U\left(t\right)=x_1\left(t\right)D_1\left(t\right)$ ，$I\left(t\right)=x_2\left(t\right)D_2\left(t\right)$ ，$V\left(t\right)=x_3\left(t\right)$ ，$C\left(t\right)=x_4\left(t\right)$ 以及$D_i\left(t\right)={\displaystyle \prod _{0<t_k<t}\left(1+R_{ik}\right)},i=1,2$ ，模型(2)的初值为：$X\left(0\right)=\left(x_1\left(0\right),x_2\left(0\right),x_3\left(0\right),x_4\left(0\right)\right)=\left(U\left(0\right),I\left(0\right),V\left(0\right),C\left(0\right)\right)=Z\left(0\right)\in R_{\text{+}}^{\text{4}}$ ，我们有以下假设：如果

$ℛ:=\underset{\left(x_1,x_2,x_3,x_4\right)\in {ℝ}_{+}^4}{\sup }\left\{\frac{r\left(t\right)\beta \left(t\right)\eta \left(t\right)\omega \left(t\right)\lambda \left(t\right)D_1\left(t\right)}{\left({\mu }_u\left(t\right)+\frac{{\sigma }_1^2\left(t\right)}{2}\right)\left(\delta \left(t\right)+\frac{{\sigma }_2^2\left(t\right)}{2}\right)\left({\mu }_v\left(t\right)+\frac{{\sigma }_3^2\left(t\right)}{2}\right)\left({\mu }_c\left(t\right)+\frac{{\sigma }_4^2\left(t\right)}{2}\right)}\right\}<1$ ，(3)

则模型(2)有唯一遍历性平稳分布。

证明：为了简化，我们把$x_1\left(t\right),x_2\left(t\right),x_3\left(t\right),x_4\left(t\right)$ 写成$x_1,x_2,x_3,x_4$ ，首先给出模型(2)的扩散矩阵$b\left(X\right)=diag\left({\sigma }_1^2\left(t\right)x_1^2,{\sigma }_2^2\left(t\right)x_2^2,{\sigma }_3^2\left(t\right)x_3^2,{\sigma }_4^2\left(t\right)x_4^2\right)$ ，选取

$G=\underset{{}_{\left(x_1,x_2,x_3,x_4\right)\in {\mathcal{U}}_{\epsilon }}}{\min }\left\{{\sigma }_1^2\left(t\right)x_1^2,{\sigma }_2^2\left(t\right)x_2^2,{\sigma }_3^2\left(t\right)x_3^2,{\sigma }_4^2\left(t\right)x_4^2\right\}$ ，则：

${\displaystyle \sum _{i,j=1}^4{b}_{i,j}{\xi }_i{\xi }_j}={\sigma }_1^2\left(t\right)x_1^2{\xi }_1^2+{\sigma }_2^2\left(t\right)x_2^2{\xi }_2^2+{\sigma }_3^2\left(t\right)x_3^2{\xi }_3^2+{\sigma }_4^2\left(t\right)x_4^2{\xi }_4^2\ge G{\Vert \xi \Vert }^2$ ，

其中$X=\left(x_1,x_2,x_3,x_4\right)\in {\mathcal{U}}_{\epsilon }$ ，$\xi =\left({\xi }_1,{\xi }_2,{\xi }_3,{\xi }_4\right)\in R_{+}^4$ ，${\mathcal{U}}_{\epsilon }$ 的定义将在后文给出，可知$b\left(X\right)$ 相对于${\mathcal{U}}_{\epsilon }$ 是非退化的。定义Lyapunov函数，令：

${\mathcal{V}}_1=\frac{1}{x_1}-\ln x_1-c_1\left(t\right)\ln x_2-c_2\left(t\right)\ln x_3-c_3\left(t\right)\ln x_4$ ，

其中$c_1\left(t\right)$ ，$c_2\left(t\right)$ 和$c_3\left(t\right)$ 是稍后确定的参数。应用Itô’s公式得：

$\text{d}{\mathcal{V}}_1=ℒ{\mathcal{V}}_1\text{d}t-{\sigma }_1\left(t\right)\left(1+\frac{1}{x_1}\right)\text{d}{B}_1\left(t\right)-c_1\left(t\right){\sigma }_2\left(t\right)\text{d}{B}_2\left(t\right)-c_2\left(t\right){\sigma }_3\left(t\right)\text{d}{B}_3\left(t\right)-c_3\left(t\right){\sigma }_4\left(t\right)\text{d}{B}_4\left(t\right),$

其中：

$\begin{array}{c}ℒ{\mathcal{V}}_1=-\frac{r\left(t\right)}{x_1}+\frac{r\left(t\right)D_1\left(t\right)}{K\left(t\right)}+\frac{1}{x_1}\left(\frac{r\left(t\right)D_2\left(t\right)}{K\left(t\right)}{x}_2+\beta \left(t\right)x_3+\alpha \left(t\right)x_4+{\mu }_u\left(t\right)+{\sigma }_1^2\left(t\right)\right)-r\left(t\right)+\frac{r\left(t\right)D_1\left(t\right)}{K\left(t\right)}{x}_1\\ +\frac{r\left(t\right)D_2\left(t\right)}{K\left(t\right)}{x}_2+\beta \left(t\right)x_3+\alpha \left(t\right)x_4+{\mu }_u\left(t\right)+\frac{{\sigma }_1^2\left(t\right)}{2}-\frac{c_1\left(t\right)\beta \left(t\right)D_1\left(t\right)x_1{x}_3}{D_2\left(t\right)x_2}+c_1\left(t\right)\left(\delta \left(t\right)+\frac{{\sigma }_2^2\left(t\right)}{2}\right)\\ +c_1\left(t\right)\alpha \left(t\right)x_4-\frac{c_2\left(t\right)\eta \left(t\right)\omega \left(t\right)D_2\left(t\right)x_2{x}_4}{x_3}-\frac{c_2\left(t\right)\eta \left(t\right)D_2\left(t\right)x_2}{x_3}+c_2\left(t\right)\left({\mu }_v\left(t\right)+\frac{{\sigma }_3^2\left(t\right)}{2}\right)\\ +c_2\left(t\right)\theta \left(t\right)x_4-\frac{c_3\left(t\right)\lambda \left(t\right)}{x_4}-c_3\left(t\right)\xi \left(t\right)x_3-c_3\left(t\right)\rho \left(t\right)\left(x_1{D}_1\left(t\right)+x_2{D}_2\left(t\right)\right)+c_3\left(t\right)\left({\mu }_c\left(t\right)+\frac{{\sigma }_4^2\left(t\right)}{2}\right)\\ \le -4\sqrt[4]{r\left(t\right)\beta \left(t\right)\eta \left(t\right)\omega \left(t\right)\lambda \left(t\right)D_1\left(t\right)c_1\left(t\right)c_2\left(t\right)c_3\left(t\right)}\\ +\frac{r\left(t\right)D_1\left(t\right)}{K\left(t\right)}+\frac{1}{x_1}\left(\frac{r\left(t\right)D_2\left(t\right)}{K\left(t\right)}{x}_2+\beta \left(t\right)x_3+\alpha \left(t\right)x_4+{\mu }_u\left(t\right)+{\sigma }_1^2\left(t\right)\right)+\frac{r\left(t\right)D_1\left(t\right)}{K\left(t\right)}{x}_1+\frac{r\left(t\right)D_2\left(t\right)}{K\left(t\right)}{x}_2\\ +\beta \left(t\right)x_3+\left(\alpha \left(t\right)+c_1\left(t\right)\alpha \left(t\right)+c_2\left(t\right)\theta \left(t\right)\right)x_4+{\mu }_u\left(t\right)+\frac{{\sigma }_1^2\left(t\right)}{2}\\ +c_1\left(t\right)\left(\delta \left(t\right)+\frac{{\sigma }_2^2\left(t\right)}{2}\right)+c_2\left(t\right)\left({\mu }_v\left(t\right)+\frac{{\sigma }_3^2\left(t\right)}{2}\right)+c_3\left(t\right)\left({\mu }_c\left(t\right)+\frac{{\sigma }_4^2\left(t\right)}{2}\right).\end{array}$

令：$\begin{array}{c}{c}_1\left(t\right)\left(\delta \left(t\right)+\frac{{\sigma }_2^2\left(t\right)}{2}\right)=c_2\left(t\right)\left({\mu }_v\left(t\right)+\frac{{\sigma }_3^2\left(t\right)}{2}\right)=c_3\left(t\right)\left({\mu }_c\left(t\right)+\frac{{\sigma }_4^2\left(t\right)}{2}\right)\\ =\frac{r\left(t\right)\beta \left(t\right)\eta \left(t\right)\omega \left(t\right)\lambda \left(t\right)D_1\left(t\right)}{\left(\delta \left(t\right)+\frac{{\sigma }_2^2\left(t\right)}{2}\right)\left({\mu }_v\left(t\right)+\frac{{\sigma }_3^2\left(t\right)}{2}\right)\left({\mu }_c\left(t\right)+\frac{{\sigma }_4^2\left(t\right)}{2}\right)}\end{array}$ ，

因此，

$ℒ{\mathcal{V}}_1\le -\frac{r\left(t\right)\beta \left(t\right)\eta \left(t\right)\omega \left(t\right)\lambda \left(t\right)D_1\left(t\right)}{\left(\delta \left(t\right)+\frac{{\sigma }_2^2\left(t\right)}{2}\right)\left({\mu }_v\left(t\right)+\frac{{\sigma }_3^2\left(t\right)}{2}\right)\left({\mu }_c\left(t\right)+\frac{{\sigma }_4^2\left(t\right)}{2}\right)}+\frac{r\left(t\right)D_1\left(t\right)}{K\left(t\right)}$

$\begin{array}{c}+\frac{1}{x_1}\left(\frac{r\left(t\right)D_2\left(t\right)}{K\left(t\right)}{x}_2+\beta \left(t\right)x_3+\alpha \left(t\right)x_4+{\mu }_u\left(t\right)+{\sigma }_1^2\left(t\right)\right)+\frac{r\left(t\right)D_1\left(t\right)}{K\left(t\right)}{x}_1\\ +\frac{r\left(t\right)D_2\left(t\right)}{K\left(t\right)}{x}_2+\beta \left(t\right)x_3+\left(\alpha \left(t\right)+c_1\left(t\right)\alpha \left(t\right)+c_2\left(t\right)\theta \left(t\right)\right)x_4+\left({\mu }_u\left(t\right)+\frac{{\sigma }_1^2\left(t\right)}{2}\right)\\ \le \left({\mu }_u\left(t\right)+\frac{{\sigma }_1^2\left(t\right)}{2}\right)\left(1-\frac{r\left(t\right)\beta \left(t\right)\eta \left(t\right)\omega \left(t\right)\lambda \left(t\right)D_1\left(t\right)}{\left({\mu }_u\left(t\right)+\frac{{\sigma }_1^2\left(t\right)}{2}\right)\left(\delta \left(t\right)+\frac{{\sigma }_2^2\left(t\right)}{2}\right)\left({\mu }_v\left(t\right)+\frac{{\sigma }_3^2\left(t\right)}{2}\right)\left({\mu }_c\left(t\right)+\frac{{\sigma }_4^2\left(t\right)}{2}\right)}\right)\\ +\frac{r\left(t\right)D_1\left(t\right)}{K\left(t\right)}{x}_1+T\\ \le -\gamma +\frac{r\left(t\right)D_1\left(t\right)}{K\left(t\right)}{x}_1+T,\end{array}$

其中：

$\begin{array}{l}\gamma =\underset{\left(x_1,x_2,x_3,x_4\right)\in {ℝ}_{+}^4}{\sup }\left\{\left({\mu }_u\left(t\right)+\frac{{\sigma }_1^2\left(t\right)}{2}\right)\left(\frac{r\left(t\right)\beta \left(t\right)\eta \left(t\right)\omega \left(t\right)\lambda \left(t\right)D_1\left(t\right)}{\left({\mu }_u\left(t\right)+\frac{{\sigma }_1^2\left(t\right)}{2}\right)\left(\delta \left(t\right)+\frac{{\sigma }_2^2\left(t\right)}{2}\right)\left({\mu }_v\left(t\right)+\frac{{\sigma }_3^2\left(t\right)}{2}\right)\left({\mu }_c\left(t\right)+\frac{{\sigma }_4^2\left(t\right)}{2}\right)}-1\right)\right\},\\ T=\underset{\left(x_1,x_2,x_3,x_4\right)\in {ℝ}_{+}^4}{\sup }\{\frac{1}{x_1}\left(\frac{r\left(t\right)D_2\left(t\right)}{K\left(t\right)}{x}_2+\beta \left(t\right)x_3+\alpha \left(t\right)x_4+{\mu }_u\left(t\right)+{\sigma }_1^2\left(t\right)\right)+\frac{r\left(t\right)D_2\left(t\right)}{K\left(t\right)}{x}_2+\beta \left(t\right)x_3\\ +\left(\alpha \left(t\right)+c_1\left(t\right)\alpha \left(t\right)+c_2\left(t\right)\theta \left(t\right)\right)x_4+\frac{r\left(t\right)D_1\left(t\right)}{K\left(t\right)}\},\end{array}$

根据(3)式可知$\gamma >0$ 。构造$C^2$ -函数$\overline{\mathcal{V}}:R_{+}^4\to R$ ，定义如下：

$\begin{array}{c}\overline{\mathcal{V}}\left(x_1,x_2,x_3,x_4\right)=M{\mathcal{V}}_1-\ln x_2-\ln x_3-\ln x_4+\frac{1}{s+1}{\left(p\left(t\right)x_1+\frac{p\left(t\right)D_2\left(t\right)}{D_1\left(t\right)}{x}_2+g\left(t\right)\frac{\delta \left(t\right)}{2\eta \left(t\right)}{x}_3+q\left(t\right)x_4\right)}^{s+1}\\ :=M{\mathcal{V}}_1+{\mathcal{V}}_2+{\mathcal{V}}_3+{\mathcal{V}}_4+{\mathcal{V}}_5\end{array}$

其中$s\in \left(0,1\right)$ ，$p\left(t\right)$ ，$g\left(t\right)$ 和$q\left(t\right)$ 是稍后确定的参数，且$M>0$ 满足以下条件：

$-M\gamma +A_1\le -2$ (4)

这里：

$\begin{array}{c}{A}_1=\underset{\left(x_1,x_2,x_3,x_4\right)\in R_{+}^4}{\sup }\{-\frac{p^{s+1}\left(t\right)r\left(t\right)D_1\left(t\right)}{2K\left(t\right)}{x}_1^{s+2}-\frac{\delta \left(t\right)p^s\left(t\right)D_2^{s+1}\left(t\right)}{4D_1^s\left(t\right)}\left(\frac{2p\left(t\right)}{D_1\left(t\right)}-g\left(t\right)\right)x_2^{s+1}\\ -{\left(\frac{g\left(t\right)\delta \left(t\right)}{2\eta \left(t\right)}\right)}^{s+1}\frac{{\mu }_v\left(t\right)}{2}{x}_3^{s+1}-q^{s+1}\left(t\right)\frac{{\mu }_c\left(t\right)}{2}{x}_4^{s+1}+Y\},\end{array}$

$Y=\underset{\left(x_1,x_2,x_3,x_4\right)\in R_{+}^4}{\sup }\left\{B+MT+\left(\alpha \left(t\right)+\theta \left(t\right)\right)x_4+\delta \left(t\right)+{\mu }_v\left(t\right)+{\mu }_c\left(t\right)+\frac{{\sigma }_2^2\left(t\right)}{\text{2}}+\frac{{\sigma }_3^2\left(t\right)}{\text{2}}+\frac{{\sigma }_4^2\left(t\right)}{\text{2}}\right\},$

$B=\underset{\left(x_1,x_2,x_3,x_4\right)\in R_{+}^4}{\sup }\{-\frac{p^{s+1}\left(t\right)r\left(t\right)D_1\left(t\right)}{2K\left(t\right)}{x}_1^{s+2}-\frac{\delta \left(t\right)p^s\left(t\right)D_2^{s+1}\left(t\right)}{4D_1^s\left(t\right)}\left(\frac{2p\left(t\right)}{D_1\left(t\right)}-g\left(t\right)\right)x_2^{s+1}-{\left(\frac{g\left(t\right)\delta \left(t\right)}{2\eta \left(t\right)}\right)}^{s+1}\frac{{\mu }_v\left(t\right)}{2}{x}_3^{s+1}$

$\begin{array}{c}-q^{s+1}\left(t\right)\frac{{\mu }_c\left(t\right)}{2}{x}_4^{s+1}+p\left(t\right)r\left(t\right)x_1(p\left(t\right)x_1+\frac{p\left(t\right)D_2\left(t\right)}{D_1\left(t\right)}{x}_2{+\frac{g\left(t\right)\delta \left(t\right)}{2\eta \left(t\right)}{x}_3+q\left(t\right)x_4)}^s\\ +q\left(t\right)\lambda \left(t\right){\left(p\left(t\right)x_1+\frac{p\left(t\right)D_2\left(t\right)}{D_1\left(t\right)}{x}_2+\frac{g\left(t\right)\delta \left(t\right)}{2\eta \left(t\right)}{x}_3+q\left(t\right)x_4\right)}^s\\ +\frac{s}{2}\left[p^{s+1}\left(t\right){\sigma }_1^2\left(t\right)x_1^{s+1}+{\left(\frac{p\left(t\right)D_2\left(t\right)}{D_1\left(t\right)}\right)}^{s+1}{\sigma }_2^2\left(t\right)x_2^{s+1}+{\left(\frac{g\left(t\right)\delta \left(t\right)}{2\eta \left(t\right)}\right)}^{s+1}{\sigma }_3^2\left(t\right)x_3^{s+1}+q^{s+1}\left(t\right){\sigma }_4^2\left(t\right)x_4^{s+1}\right]\}\\ <\infty .\end{array}$

因此得出：

$\underset{n\to \infty ,\left(x_1,x_2,x_3,x_4\right)\in R_{+}^4\backslash {\mathcal{U}}_n}{\lim }\overline{\mathcal{V}}\left(x_1,x_2,x_3,x_4\right)=+\infty$ ，

其中${\mathcal{U}}_n=\left(\frac{1}{n},n\right)\times \left(\frac{1}{n},n\right)\times \left(\frac{1}{n},n\right)\times \left(\frac{1}{n},n\right)$ 。此外，$\overline{\mathcal{V}}\left(x_1,x_2,x_3,x_4\right)$ 是连续函数并且存在最小下界，假设在$R_{+}^4$ 中有一点$\left({\overline{x}}_1,{\overline{x}}_2,{\overline{x}}_3,{\overline{x}}_4\right)$ ，函数$\overline{\mathcal{V}}$ 在该点处达到最小值。定义非负的$C^2$ -函数$\widehat{\mathcal{V}}:R_{+}^4\to R$ ，如下：

$\widehat{\mathcal{V}}\left(x_1,x_2,x_3,x_4\right)=\overline{\mathcal{V}}\left(x_1,x_2,x_3,x_4\right)-\overline{\mathcal{V}}\left({\overline{x}}_1,{\overline{x}}_2,{\overline{x}}_3,{\overline{x}}_4\right).$

同理，分别由作用在函数${\mathcal{V}}_2,{\mathcal{V}}_3,{\mathcal{V}}_4,{\mathcal{V}}_5$ 上的微分算子可得：

$ℒ{\mathcal{V}}_2=-\frac{\beta \left(t\right)D_1\left(t\right)x_1{x}_3}{D_2\left(t\right)x_2}+\delta \left(t\right)+\alpha \left(t\right)x_4+\frac{{\sigma }_2^2\left(t\right)}{2},$

$ℒ{\mathcal{V}}_3=-\frac{\eta \left(t\right)\omega \left(t\right)D_2\left(t\right)x_2{x}_4}{x_3}-\frac{\eta \left(t\right)D_2\left(t\right)x_2}{x_3}+\theta \left(t\right)x_4+{\mu }_v\left(t\right)+\frac{{\sigma }_3^2\left(t\right)}{2},$

$ℒ{\mathcal{V}}_4=-\frac{\lambda \left(t\right)}{x_4}-\xi \left(t\right)x_3-\rho \left(t\right)\left(x_1{D}_1\left(t\right)+x_2{D}_2\left(t\right)\right)+{\mu }_c\left(t\right)+\frac{{\sigma }_4^2\left(t\right)}{2},$

$\begin{array}{c}ℒ{\mathcal{V}}_5={\left(p\left(t\right)x_1+\frac{p\left(t\right)D_2\left(t\right)}{D_1\left(t\right)}{x}_2+\frac{g\left(t\right)\delta \left(t\right)}{2\eta \left(t\right)}{x}_3+q\left(t\right)x_4\right)}^{s}d\left(p\left(t\right)x_1+\frac{p\left(t\right)D_2\left(t\right)}{D_1\left(t\right)}{x}_2+\frac{g\left(t\right)\delta \left(t\right)}{2\eta \left(t\right)}{x}_3+q\left(t\right)x_4\right)\\ +\frac{s}{2}{\left(p\left(t\right)x_1+\frac{p\left(t\right)D_2\left(t\right)}{D_1\left(t\right)}{x}_2+\frac{g\left(t\right)\delta \left(t\right)}{2\eta \left(t\right)}{x}_3+q\left(t\right)x_4\right)}^{s-1}[p^2\left(t\right){\sigma }_1^2\left(t\right)x_1^2+{\left(\frac{p\left(t\right)D_2\left(t\right)}{D_1\left(t\right)}\right)}^2{\sigma }_2^2\left(t\right)x_2^2\\ +{\left(\frac{g\left(t\right)\delta \left(t\right)}{2\eta \left(t\right)}\right)}^2{\sigma }_3^2\left(t\right)x_3^2+q^2\left(t\right){\sigma }_4^2\left(t\right)x_4^2]\\ \le \left(-\frac{p\left(t\right)r\left(t\right)D_1\left(t\right)}{K\left(t\right)}{x}_1^2-\frac{\delta \left(t\right)D_2\left(t\right)}{2}\left(\frac{2p\left(t\right)}{D_1\left(t\right)}-g\left(t\right)\right)x_2-\frac{g\left(t\right)\delta \left(t\right){\mu }_v\left(t\right)}{2\eta \left(t\right)}{x}_3-q\left(t\right){\mu }_c\left(t\right)x_4\right)\\ \times (p\left(t\right)x_1+\frac{p\left(t\right)D_2\left(t\right)}{D_1\left(t\right)}{x}_2{+\frac{g\left(t\right)\delta \left(t\right)}{2\eta \left(t\right)}{x}_3+q\left(t\right)x_4)}^s+p\left(t\right)r\left(t\right)x_1(p\left(t\right)x_1+\frac{p\left(t\right)D_2\left(t\right)}{D_1\left(t\right)}{x}_2\\ +{\frac{g\left(t\right)\delta \left(t\right)}{2\eta \left(t\right)}{x}_3+q\left(t\right)x_4)}^s+q\left(t\right)\lambda \left(t\right){\left(p\left(t\right)x_1+\frac{p\left(t\right)D_2\left(t\right)}{D_1\left(t\right)}{x}_2+\frac{g\left(t\right)\delta \left(t\right)}{2\eta \left(t\right)}{x}_3+q\left(t\right)x_4\right)}^s\\ +D_1\left(t\right)\left(q\left(t\right)\rho \left(t\right)-\frac{p\left(t\right)\alpha \left(t\right)}{D_1\left(t\right)}\right)x_1{x}_4{\left(p\left(t\right)x_1+\frac{p\left(t\right)D_2\left(t\right)}{D_1\left(t\right)}{x}_2+\frac{g\left(t\right)\delta \left(t\right)}{2\eta \left(t\right)}{x}_3+q\left(t\right)x_4\right)}^s\end{array}$

$\begin{array}{c}+D_2\left(t\right)\left(\frac{g\left(t\right)\delta \left(t\right)\omega \left(t\right)}{2}+q\left(t\right)\rho \left(t\right)-\frac{p\left(t\right)\alpha \left(t\right)}{D_1\left(t\right)}\right)x_2{x}_4(p\left(t\right)x_1\left(t\right)+\frac{p\left(t\right)D_2\left(t\right)}{D_1\left(t\right)}{x}_2\\ +{\frac{g\left(t\right)\delta \left(t\right)}{2\eta \left(t\right)}{x}_3+q\left(t\right)x_4)}^s+\left(q\left(t\right)\xi \left(t\right)-\frac{g\left(t\right)\delta \left(t\right)\theta \left(t\right)}{2\eta \left(t\right)}\right)x_3{x}_4(p\left(t\right)x_1+\frac{p\left(t\right)D_2\left(t\right)}{D_1\left(t\right)}{x}_2\\ +{\frac{g\left(t\right)\delta \left(t\right)}{2\eta \left(t\right)}{x}_3+q\left(t\right)x_4)}^s+\frac{s}{2}(p\left(t\right)x_1+\frac{p\left(t\right)D_2\left(t\right)}{D_1\left(t\right)}{x}_2{+\frac{g\left(t\right)\delta \left(t\right)}{2\eta \left(t\right)}{x}_3+q\left(t\right)x_4)}^{s-1}\\ \times \left[p^2\left(t\right){\sigma }_1^2\left(t\right)x_1^2+{\left(\frac{p\left(t\right)D_2\left(t\right)}{D_1\left(t\right)}\right)}^2{\sigma }_2^2\left(t\right)x_2^2+{\left(\frac{g\left(t\right)\delta \left(t\right)}{2\eta \left(t\right)}\right)}^2{\sigma }_3^2\left(t\right)x_3^2+q^2\left(t\right){\sigma }_4^2\left(t\right)x_4^2\right].\end{array}$

存在$p\left(t\right),g\left(t\right)$ 满足$\frac{2p\left(t\right)}{D_1\left(t\right)}-g\left(t\right)\ge 0$ 且$q\left(t\right)\in \left(0,\frac{p\left(t\right)\alpha \left(t\right)}{\rho \left(t\right)D_1\left(t\right)}-\frac{g\left(t\right)\delta \left(t\right)\omega \left(t\right)}{2\rho \left(t\right)}\right]\cap \left(0,\frac{g\left(t\right)\delta \left(t\right)\theta \left(t\right)}{2\eta \left(t\right)\xi \left(t\right)}\right]$ ，因此，

$\begin{array}{c}ℒ{\mathcal{V}}_5\le \left(-\frac{p\left(t\right)r\left(t\right)D_1\left(t\right)}{K\left(t\right)}{x}_1^2-\frac{\delta \left(t\right)D_2\left(t\right)}{2}\left(\frac{2p\left(t\right)}{D_1\left(t\right)}-g\left(t\right)\right)x_2-\frac{g\left(t\right)\delta \left(t\right){\mu }_v\left(t\right)}{2\eta \left(t\right)}{x}_3-q\left(t\right){\mu }_c\left(t\right)x_4\right)\\ \times (p\left(t\right)x_1+\frac{p\left(t\right)D_2\left(t\right)}{D_1\left(t\right)}{x}_2{+\frac{g\left(t\right)\delta \left(t\right)}{2\eta \left(t\right)}{x}_3+q\left(t\right)x_4)}^s+p\left(t\right)r\left(t\right)x_1(p\left(t\right)x_1+\frac{p\left(t\right)D_2\left(t\right)}{D_1\left(t\right)}{x}_2\\ +{\frac{g\left(t\right)\delta \left(t\right)}{2\eta \left(t\right)}{x}_3+q\left(t\right)x_4)}^s+q\left(t\right)\lambda (p\left(t\right)x_1+\frac{p\left(t\right)D_2\left(t\right)}{D_1\left(t\right)}{x}_2{+\frac{g\left(t\right)\delta \left(t\right)}{2\eta \left(t\right)}{x}_3+q\left(t\right)x_4)}^s\\ +\frac{s}{2}\left[p^{s+1}\left(t\right){\sigma }_1^2\left(t\right)x_1^{s+1}+{\left(\frac{p\left(t\right)D_2\left(t\right)}{D_1\left(t\right)}\right)}^{s+1}{\sigma }_2^2\left(t\right)x_2^{s+1}+{\left(\frac{g\left(t\right)\delta \left(t\right)}{2\eta \left(t\right)}\right)}^{s+1}{\sigma }_3^2\left(t\right)x_3^{s+1}+q^{s+1}{\sigma }_4^2\left(t\right)x_4^{s+1}\right]\\ \le -\frac{p^{s+1}\left(t\right)r\left(t\right)D_1\left(t\right)}{2K\left(t\right)}{x}_1^{s+2}-\frac{\delta \left(t\right)p^s\left(t\right)D_2^{s+1}\left(t\right)}{4D_1^s\left(t\right)}\left(\frac{2p\left(t\right)}{D_1\left(t\right)}-g\left(t\right)\right)x_2^{s+1}-{\left(\frac{g\left(t\right)\delta \left(t\right)}{2\eta \left(t\right)}\right)}^{s+1}\frac{{\mu }_v\left(t\right)}{2}{x}_3^{s+1}\\ -q^{s+1}\left(t\right)\frac{{\mu }_c\left(t\right)}{2}{x}_4^{s+1}+B.\end{array}$

综上可得：

$\begin{array}{c}ℒ\widehat{\mathcal{V}}\le -M\gamma +\frac{Mr\left(t\right)D_1\left(t\right)}{K\left(t\right)}{x}_1-\frac{\beta \left(t\right)D_1\left(t\right)x_1{x}_3}{D_2\left(t\right)x_2}-\frac{\eta \left(t\right)D_2\left(t\right)x_2}{x_3}-\frac{\lambda \left(t\right)}{x_4}-\frac{p^{s+1}\left(t\right)r\left(t\right)D_1\left(t\right)}{2K\left(t\right)}{x}_1^{s+2}\\ -\frac{\delta \left(t\right)p^s\left(t\right)D_2^{s+1}\left(t\right)}{4D_1^s\left(t\right)}\left(\frac{2p\left(t\right)}{D_1\left(t\right)}-g\left(t\right)\right)x_2^{s+1}-{\left(\frac{g\left(t\right)\delta \left(t\right)}{2\eta \left(t\right)}\right)}^{s+1}\frac{{\mu }_v\left(t\right)}{2}{x}_3^{s+1}-q^{s+1}\left(t\right)\frac{{\mu }_c\left(t\right)}{2}{x}_4^{s+1}+Y.\end{array}$

现在我们可以构造一个紧子集${\mathcal{U}}_{\epsilon }$ ，定义以下有界闭集：

${\mathcal{U}}_{\epsilon }=\left\{\epsilon \le x_1\le \frac{1}{\epsilon },{\epsilon }^3\le x_2\le \frac{1}{{\epsilon }^3},{\epsilon }^2\le x_3\le \frac{1}{{\epsilon }^2},\epsilon \le x_4\le \frac{1}{\epsilon }\right\},$

其中$\epsilon >0$ 是一个足够小的正数。在集合$R_{+}^4\backslash {\mathcal{U}}_{\epsilon }$ 中，我们选择足够小的$\epsilon$ ，使得以下条件成立：

$-M\gamma +\frac{Mr{D}_1}{K}\epsilon +A_1\le -1,$ (5)

$-\frac{\beta D_1}{D_2\epsilon }+A_2\le -1,$ (6)

$-\frac{\eta D_2}{\epsilon }+A_2\le -1,$ (7)

$-\frac{\lambda }{\epsilon }+A_2\le -1,$ (8)

$-\frac{p^{s+1}r{D}_1}{2K}\frac{1}{{\epsilon }^{s+2}}+A_3\le -1,$ (9)

$-\frac{\delta p^s{D}_2^{s+1}}{4D_1^s}\left(\frac{2p}{D_1}-g\right)\frac{1}{{\epsilon }^{3\left(s+1\right)}}+A_4\le -1,$ (10)

$-{\left(\frac{g\delta }{2\eta }\right)}^{s+1}\frac{{\mu }_v}{2}\frac{1}{{\epsilon }^{2\left(s+1\right)}}+A_5\le -1,$ (11)

$-q^{s+1}\frac{{\mu }_c}{2}\frac{1}{{\epsilon }^{s+1}}+A_6\le -1.$ (12)

其中$A_1,A_2,A_3,\cdots ,A_6$ 是正常数，由下面的不等式(13) (14)，(17)~(20)给出。为了方便起见，我们将$R_{+}^4\backslash {\mathcal{U}}_{\epsilon }$ 划分为以下八个区域：

$\begin{array}{l}{\mathcal{U}}_1=\left\{0<x_1<\epsilon \right\},\\ {\mathcal{U}}_3=\left\{0<x_3<{\epsilon }^2,x_2\ge \epsilon \right\},\\ {\mathcal{U}}_5=\left\{x_1>\frac{1}{\epsilon }\right\},\\ {\mathcal{U}}_7=\left\{x_3>\frac{1}{{\epsilon }^2}\right\},\end{array}$ $\begin{array}{l}{\mathcal{U}}_2=\left\{0<x_2<{\epsilon }^3,x_1\ge \epsilon ,x_3\ge \epsilon \right\},\\ {\mathcal{U}}_4=\left\{0<x_4<\epsilon \right\},\\ {\mathcal{U}}_6=\left\{x_2>\frac{1}{{\epsilon }^3}\right\},\\ {\mathcal{U}}_8=\left\{x_4>\frac{1}{\epsilon }\right\}.\end{array}$

显然，${\mathcal{U}}_{\epsilon }={\mathcal{U}}_1\cup {\mathcal{U}}_2\cup {\mathcal{U}}_3\cup {\mathcal{U}}_4\cup {\mathcal{U}}_5\cup {\mathcal{U}}_6\cup {\mathcal{U}}_7\cup {\mathcal{U}}_8$ 。接下来，我们证明${\mathcal{U}}_{\epsilon }$ 上$ℒ\widehat{\mathcal{V}}\left(x_1,x_2,x_3,x_4\right)\le -1$ ，即等价于在上述八个定义域上证明。

情形1，若$\left(x_1,x_2,x_3,x_4\right)\in {\mathcal{U}}_1$ ，则：

$\begin{array}{c}ℒ\widehat{\mathcal{V}}\le -M\gamma +\frac{Mr\left(t\right)D_1\left(t\right)}{K\left(t\right)}{x}_1-\frac{p^{s+1}\left(t\right)r\left(t\right)D_1\left(t\right)}{2K\left(t\right)}{x}_1^{s+2}-\frac{\delta \left(t\right)p^s\left(t\right)D_2^{s+1}\left(t\right)}{4D_1^s\left(t\right)}\left(\frac{2p\left(t\right)}{D_1\left(t\right)}-g\left(t\right)\right)x_2^{s+1}\\ -{\left(\frac{g\left(t\right)\delta \left(t\right)}{2\eta \left(t\right)}\right)}^{s+1}\frac{{\mu }_v\left(t\right)}{2}{x}_3^{s+1}-q^{s+1}\left(t\right)\frac{{\mu }_c\left(t\right)}{2}{x}_4^{s+1}+Y\\ \le -M\gamma +\frac{Mr\left(t\right)D_1\left(t\right)}{K\left(t\right)}{x}_1+A_1\\ \le -M\gamma +\frac{M{r}^u{D}_1^u}{K^l}\epsilon +A_1,\end{array}$(13)

其中：

$\begin{array}{c}{A}_1=\underset{\left(x_1,x_2,x_3,x_4\right)\in R_{+}^4}{\sup }\{-\frac{p^{s+1}\left(t\right)r\left(t\right)D_1\left(t\right)}{2K\left(t\right)}{x}_1^{s+2}-\frac{\delta \left(t\right)p^s\left(t\right)D_2^{s+1}\left(t\right)}{4D_1^s\left(t\right)}\left(\frac{2p\left(t\right)}{D_1\left(t\right)}-g\left(t\right)\right)x_2^{s+1}\\ -{\left(\frac{g\left(t\right)\delta \left(t\right)}{2\eta \left(t\right)}\right)}^{s+1}\frac{{\mu }_v\left(t\right)}{2}{x}_3^{s+1}-q^{s+1}\left(t\right)\frac{{\mu }_c\left(t\right)}{2}{x}_4^{s+1}+Y\}.\end{array}$

根据(4)和(5)式，可得$ℒ\widehat{\mathcal{V}}\le -1$ 在${\mathcal{U}}_1$ 上成立。

情形2，若$\left(x_1,x_2,x_3,x_4\right)\in {\mathcal{U}}_2$ ，我们有：

$\begin{array}{c}ℒ\widehat{\mathcal{V}}\le -\frac{\beta \left(t\right)D_1\left(t\right)x_1{x}_3}{D_2\left(t\right)x_2}+\frac{Mr\left(t\right)D_1\left(t\right)}{K\left(t\right)}{x}_1-\frac{p^{s+1}\left(t\right)r\left(t\right)D_1\left(t\right)}{2K\left(t\right)}{x}_1^{s+2}\\ -\frac{\delta \left(t\right)p^s\left(t\right)D_2^{s+1}\left(t\right)}{4D_1^s\left(t\right)}\left(\frac{2p\left(t\right)}{D_1\left(t\right)}-g\left(t\right)\right)x_2^{s+1}\\ -{\left(\frac{g\left(t\right)\delta \left(t\right)}{2\eta \left(t\right)}\right)}^{s+1}\frac{{\mu }_v\left(t\right)}{2}{x}_3^{s+1}-q^{s+1}\left(t\right)\frac{{\mu }_c\left(t\right)}{2}{x}_4^{s+1}+Y\\ \le -\frac{\beta \left(t\right)D_1\left(t\right)x_1{x}_3}{D_2\left(t\right)x_2}+A_2\\ \le -\frac{{\beta }^l{D}_1^l}{D_2^u\epsilon }+A_2,\end{array}$(14)

其中：

$\begin{array}{c}{A}_2=\underset{\left(x_1,x_2,x_3,x_4\right)\in R_{+}^4}{\sup }\{\frac{Mr\left(t\right)D_1\left(t\right)}{K\left(t\right)}{x}_1-\frac{p^{s+1}\left(t\right)r\left(t\right)D_1\left(t\right)}{2K\left(t\right)}{x}_1^{s+2}-\frac{\delta \left(t\right)p^s\left(t\right)D_2^{s+1}\left(t\right)}{4D_1^s\left(t\right)}\left(\frac{2p\left(t\right)}{D_1\left(t\right)}-g\left(t\right)\right)x_2^{s+1}\\ -{\left(\frac{g\left(t\right)\delta \left(t\right)}{2\eta \left(t\right)}\right)}^{s+1}\frac{{\mu }_v\left(t\right)}{2}{x}_3^{s+1}-q^{s+1}\left(t\right)\frac{{\mu }_c\left(t\right)}{2}{x}_4^{s+1}+Y\}.\end{array}$

根据(6)式，可得$ℒ\widehat{\mathcal{V}}\le -1$ 在${\mathcal{U}}_2$ 上成立。

情形3，若$\left(x_1,x_2,x_3,x_4\right)\in {\mathcal{U}}_3$ ，则：

$\begin{array}{c}ℒ\widehat{\mathcal{V}}\le -\frac{\eta \left(t\right)D_2\left(t\right)x_2}{x_3}+\frac{Mr\left(t\right)D_1\left(t\right)}{K\left(t\right)}{x}_1-\frac{p^{s+1}\left(t\right)r\left(t\right)D_1\left(t\right)}{2K\left(t\right)}{x}_1^{s+2}-\frac{\delta \left(t\right)p^s\left(t\right)D_2^{s+1}\left(t\right)}{4D_1^s\left(t\right)}\left(\frac{2p\left(t\right)}{D_1\left(t\right)}-g\left(t\right)\right)x_2^{s+1}\\ -{\left(\frac{g\left(t\right)\delta \left(t\right)}{2\eta \left(t\right)}\right)}^{s+1}\frac{{\mu }_v\left(t\right)}{2}{x}_3^{s+1}-q^{s+1}\left(t\right)\frac{{\mu }_c\left(t\right)}{2}{x}_4^{s+1}+Y\\ \le -\frac{\eta \left(t\right)D_2\left(t\right)x_2}{x_3}+A_2\\ \le -\frac{{\eta }^l{D}_2^l}{\epsilon }+A_2.\end{array}$(15)

根据(7)式，可知$ℒ\widehat{\mathcal{V}}\le -1$ 在${\mathcal{U}}_3$ 上成立。

情形4，若$\left(x_1,x_2,x_3,x_4\right)\in {\mathcal{U}}_4$ ，我们有：

$\begin{array}{c}ℒ\widehat{\mathcal{V}}\le -\frac{\lambda (t)}{x_4}+\frac{Mr\left(t\right)D_1\left(t\right)}{K\left(t\right)}{x}_1-\frac{p^{s+1}\left(t\right)r\left(t\right)D_1\left(t\right)}{2K\left(t\right)}{x}_1^{s+2}-\frac{\delta \left(t\right)p^s\left(t\right)D_2^{s+1}\left(t\right)}{4D_1^s\left(t\right)}\left(\frac{2p\left(t\right)}{D_1\left(t\right)}-g\left(t\right)\right)x_2^{s+1}\\ -{\left(\frac{g\left(t\right)\delta \left(t\right)}{2\eta \left(t\right)}\right)}^{s+1}\frac{{\mu }_v\left(t\right)}{2}{x}_3^{s+1}-q^{s+1}\left(t\right)\frac{{\mu }_c\left(t\right)}{2}{x}_4^{s+1}+Y\\ \le -\frac{\lambda \left(t\right)}{x_4}+A_2\\ \le -\frac{{\lambda }^l}{\epsilon }+A_2.\end{array}$(16)

通过(8)式，有$ℒ\widehat{\mathcal{V}}\le -1$ 在${\mathcal{U}}_4$ 上成立。

情形5，若$\left(x_1,x_2,x_3,x_4\right)\in {\mathcal{U}}_5$ ，我们有：

$\begin{array}{c}ℒ\widehat{\mathcal{V}}\le -\frac{p^{s+1}r{D}_1}{2K}{x}_1^{s+2}+\frac{Mr{D}_1}{K}{x}_1-\frac{\delta p^s{D}_2^{s+1}}{4D_1^s}\left(\frac{2p}{D_1}-g\right)x_2^{s+1}-{\left(\frac{g\delta }{2\eta }\right)}^{s+1}\frac{{\mu }_v}{2}{x}_3^{s+1}-q^{s+1}\frac{{\mu }_c}{2}{x}_4^{s+1}+Y\\ \le -\frac{p^{s+1}r{D}_1}{2K}{x}_1^{s+2}+A_3\\ \le -\frac{p^{s+1}r{D}_1}{2K}\frac{1}{{\epsilon }^{s+2}}+A_3,\end{array}$ (17)

其中：$A_3=\underset{\left(x_1,x_2,x_3,x_4\right)\in R_{+}^4}{\sup }\left\{\frac{Mr{D}_1}{K}{x}_1-\frac{\delta p^s{D}_2^{s+1}}{4D_1^s}\left(\frac{2p}{D_1}-g\right)x_2^{s+1}-{\left(\frac{g\delta }{2\eta }\right)}^{s+1}\frac{{\mu }_v}{2}{x}_3^{s+1}-q^{s+1}\frac{{\mu }_c}{2}{x}_4^{s+1}+Y\right\}$ 。

从(9)式，可以得出$ℒ\widehat{\mathcal{V}}\le -1$ 在${\mathcal{U}}_5$ 上成立。

情形6，若$\left(x_1,x_2,x_3,x_4\right)\in {\mathcal{U}}_6$ ，我们有：

$\begin{array}{c}ℒ\widehat{\mathcal{V}}\le -\frac{\delta p^s{D}_2^{s+1}}{4D_1^s}\left(\frac{2p}{D_1}-g\right)x_2^{s+1}+\frac{Mr{D}_1}{K}{x}_1-\frac{p^{s+1}r{D}_1}{2K}{x}_1^{s+2}-{\left(\frac{g\delta }{2\eta }\right)}^{s+1}\frac{{\mu }_v}{2}{x}_3^{s+1}-q^{s+1}\frac{{\mu }_c}{2}{x}_4^{s+1}+Y\\ \le -\frac{\delta p^s{D}_2^{s+1}}{4D_1^s}\left(\frac{2p}{D_1}-g\right)x_2^{s+1}+A_4\\ \le -\frac{\delta p^s{D}_2^{s+1}}{4D_1^s}\left(\frac{2p}{D_1}-g\right)\frac{1}{{\epsilon }^{3(s+1)}}+A_4,\end{array}$ (18)

其中：$A_4=\underset{\left(x_1,x_2,x_3,x_4\right)\in R_{+}^4}{\sup }\left\{\frac{Mr{D}_1}{K}{x}_1-\frac{p^{s+1}r{D}_1}{2K}{x}_1^{s+2}-{\left(\frac{g\delta }{2\eta }\right)}^{s+1}\frac{{\mu }_v}{2}{x}_3^{s+1}-q^{s+1}\frac{{\mu }_c}{2}{x}_4^{s+1}+Y\right\}$ 。

根据(10)式，可知$ℒ\widehat{\mathcal{V}}\le -1$ 在${\mathcal{U}}_6$ 上成立。

情形7，若$\left(x_1,x_2,x_3,x_4\right)\in {\mathcal{U}}_7$ ，则：

$\begin{array}{c}ℒ\widehat{\mathcal{V}}\le -{\left(\frac{g\delta }{2\eta }\right)}^{s+1}\frac{{\mu }_v}{2}{x}_3^{s+1}+\frac{Mr{D}_1}{K}{x}_1-\frac{p^{s+1}r{D}_1}{2K}{x}_1^{s+2}-\frac{\delta p^s{D}_2^{s+1}}{4D_1^s}\left(\frac{2p}{D_1}-g\right)x_2^{s+1}-q^{s+1}\frac{{\mu }_c}{2}{x}_4^{s+1}+Y\\ \le -{\left(\frac{g\delta }{2\eta }\right)}^{s+1}\frac{{\mu }_v}{2}{x}_3^{s+1}+A_5\\ \le -{\left(\frac{g\delta }{2\eta }\right)}^{s+1}\frac{{\mu }_v}{2}\frac{1}{{\epsilon }^{2(s+1)}}+A_5,\end{array}$ (19)

其中：$A_5=\underset{\left(x_1,x_2,x_3,x_4\right)\in R_{+}^4}{\sup }\left\{\frac{Mr{D}_1}{K}{x}_1-\frac{p^{s+1}r{D}_1}{2K}{x}_1^{s+2}-\frac{\delta p^s{D}_2^{s+1}}{4D_1^s}\left(\frac{2p}{D_1}-g\right)x_2^{s+1}-q^{s+1}\frac{{\mu }_c}{2}{x}_4^{s+1}+Y\right\}$ 。

根据(11)式，可得$ℒ\widehat{\mathcal{V}}\le -1$ 在${\mathcal{U}}_7$ 上成立。

情形8，若$\left(x_1,x_2,x_3,x_4\right)\in {\mathcal{U}}_8$ ，我们有：

$\begin{array}{c}ℒ\widehat{\mathcal{V}}\le -q^{s+1}\frac{{\mu }_c}{2}{x}_4^{s+1}+\frac{Mr{D}_1}{K}{x}_1-\frac{p^{s+1}r{D}_1}{2K}{x}_1^{s+2}-\frac{\delta p^s{D}_2^{s+1}}{4D_1^s}\left(\frac{2p}{D_1}-g\right)x_2^{s+1}-{\left(\frac{g\delta }{2\eta }\right)}^{s+1}\frac{{\mu }_v}{2}{x}_3^{s+1}+Y\\ \le -q^{s+1}\frac{{\mu }_c}{2}{x}_4^{s+1}+A_6\\ \le -q^{s+1}\frac{{\mu }_c}{2}\frac{1}{{\epsilon }^{s+1}}+A_6,\end{array}$ (20)

其中：$A_6=\underset{\left(x_1,x_2,x_3,x_4\right)\in R_{+}^4}{\sup }\left\{\frac{Mr{D}_1}{K}{x}_1-\frac{p^{s+1}r{D}_1}{2K}{x}_1^{s+2}-\frac{\delta p^s{D}_2^{s+1}}{4D_1^s}\left(\frac{2p}{D_1}-g\right)x_2^{s+1}-{\left(\frac{g\delta }{2\eta }\right)}^{s+1}\frac{{\mu }_v}{2}{x}_3^{s+1}+Y\right\}$ 。

根据(12)式，可得$ℒ\widehat{\mathcal{V}}\le -1$ 在${\mathcal{U}}_8$ 上成立。

显然，根据(13)~(20)式可得，对充分小的$\epsilon$ 有：

$ℒ\widehat{\mathcal{V}}\le -1,\left(x_1,x_2,x_3,x_4\right)\in R_{+}^4\backslash {\mathcal{U}}_{\epsilon }.$

因此，系统(2)具有唯一遍历性平稳分布。我们考虑$R_{ik}\left(i=1,2\right)$ 在每个脉冲时刻$t_k$ 取固定值或考虑周期性的脉冲效应，这种周期性尺度变化不会破坏系统的遍历性，那么系统(1)具有唯一遍历性平稳分布。

3. 数值模拟

在下文中，我们将验证本文所得结果的有效性并给出肿瘤脉冲治疗的建议，我们采用Milsteins高阶方法[17]进行数值模拟。

我们模拟了模型(1)的平稳分布情况，见图1，其参数值固定如下：$r\left(t\right)=0.206+0.01\sin t$ ， $K\left(t\right)=1.404\times {10}^8+0.01\sin t$ ，$\beta \left(t\right)=4.7\times {10}^{-3}+0.01\sin t$ ，$\alpha \left(t\right)=0.01+0.01\sin t$ ，${\mu }_u\left(t\right)=0.1+0.01\sin t$ ，$\delta \left(t\right)=0.14+0.01\sin t$ ，$\eta \left(t\right)=140+0.01\sin t$ ，$\omega \left(t\right)=0.1+0.01\sin t$ ，$\theta \left(t\right)=2+0.01\sin t$ ， ${\mu }_v\left(t\right)=0.6+0.01\sin t$ ，$\lambda \left(t\right)=1.65+0.01\sin t$ ，$\xi \left(t\right)=1\times {10}^{-8}+0.01\sin t$ ，$\rho \left(t\right)=0.265+0.01\sin t$ ， ${\mu }_c\left(t\right)=0.6+0.01\sin t$ ，$R_{1k}=R_{2k}=0.1$ ，此外，初始值和白噪声的具体取值见图例所示。

Figure 1. The stationary distribution of Model (1) in deterministic and stochastic models

图1. 确定性模型与随机模型中模型(1)的平稳分布

经过计算可知：当$D_1\left(t\right)>0.246$ 时，$ℛ\approx 4.08>1$ 。当$R_{1k}>0$ 时，根据定理2.2，模型(1)具有唯一遍历平稳分布；当$R_{1k}<0$ 时，我们给出一些有关肿瘤治疗的策略：假设$R_{1k}=-0.1$ ，则：

$D_1\left(t\right)={\displaystyle \prod _{0<t_k<t}\left(1+R_{1k}\right)}={0.9}^n>0.246$ ，其中$n$ 是满足$0<t_k<t$ 中$t_k$ 的个数，求解可得$n<\frac{-0.607}{-0.046}\approx 13.2$ ，所以

$n\le 13$ ，这意味着$n\le 13$ 时$D_1\left(t\right)>0.246$ ，即$ℛ>1$ 成立，模型(1)具有唯一遍历平稳分布。也就是说，在时间段$0<t_k<t$ 内，最多允许进行13次CAR-T细胞注射治疗，以确保$D_1\left(t\right)>0.246$ 。如果治疗次数超过13次，则$D_1\left(t\right)$ 将降到0.246以下，模型(1)未达到平稳分布的条件，表示治疗效果下降或有其他负面影响。因此，这个约束条件可以用来控制治疗的频率或间隔，避免治疗过度而导致效果减弱。这与实际应用中的治疗方案设计有关，通常在CAR-T细胞治疗中，治疗次数和间隔时间需要精心设计，以达到最佳的治疗效果。

4. 结论

本文建立了一种具有脉冲输注和随机扰动的肿瘤免疫动力学模型，理论研究了系统的平稳分布和遍历性，揭示其长期演化的特性和稳定状态。接着通过数值模拟验证理论分析的正确性，进一步地，基于平稳分布条件，推导了脉冲治疗次数的约束条件，探讨脉冲输注策略对肿瘤控制的作用。研究表明过度输注可能影响治疗效果。研究结果为CAR-T细胞联合溶瘤病毒疗法的优化提供了数学依据。未来研究将进一步优化脉冲输注策略，探索不同治疗参数对系统稳定性的影响或可结合临床数据对模型进行验证，推动数学模型在CAR-T细胞治疗中的应用。

基金项目

国家自然科学基金(12271068)。

NOTES

^*通讯作者。