一类具有Michaelis-Menten型抑制项的肿瘤-免疫模型动力学分析

2020-08-04张改梅王晓琴李建全

科学技术与工程 2020年18期

张改梅, 王晓琴, 李建全

(陕西科技大学文理学院，西安 710021)

癌症已成为严重危害人类健康的重大公共卫生问题。近期世界卫生组织发表了最新的全球癌症统计数据，2018年全球约有1.81×107例新的癌症病例和9.6×106例癌症死亡病例。由于人口的增长和人口老龄化，全球癌症发病率和死亡率快速增长[1-2]。研究表明免疫系统能够识别和消除恶性肿瘤[3-5]，数学模型模拟肿瘤免疫动力学不仅有助于理解免疫细胞和肿瘤细胞是如何相互作用的，而且也能提供一个有用的工具来预测免疫治疗的结果和改善治疗策略。

作为肿瘤免疫效应阶段的执行场所，肿瘤微环境内存在许多因素参与肿瘤组织与免疫系统的相互作用。由于肿瘤微环境的复杂性，不仅要考虑肿瘤细胞和效应细胞自身的增殖和死亡、免疫效应细胞对肿瘤细胞的抑制作用，也需要考虑肿瘤对效应细胞的作用。肿瘤的出现会刺激效应细胞的增长，同时，肿瘤对效应细胞也有抑制作用，如肿瘤可通过降低微环境中葡萄糖的浓度来抑制免疫效应细胞的活性；肿瘤微环境常呈现酸性，在酸性条件下，免疫效应细胞的活性明显降低，而且肿瘤细胞本身也参与肿瘤微环境中的免疫抑制[6]。

关于肿瘤免疫相互作用的模型已有许多学者研究[7-10]。1994年，Kuznetsov等[11]提出了二维肿瘤免疫反应的数学模型:

(1)

Delisi等[12]、Adam[13]研究了肿瘤细胞与效应细胞相互作用的常微分方程模型，这两个模型均采用Michaelis-Menten型抑制函数来表示效应细胞对肿瘤细胞的抑制作用。研究结果均显示一定范围内，效应细胞的增长会增大肿瘤细胞的存活率。此外，他们还给出了肿瘤生长不受免疫系统控制变为恶性肿瘤的阈值条件。Kirschner等[14]在文献[11-12]的基础上，假设肿瘤对效应细胞的刺激项为线性项，首次提出了肿瘤细胞、免疫效应细胞和IL-2之间的相互作用的数学模型。研究发现肿瘤的抗原性对该模型的动力学性质有着非常大的影响，并解释了肿瘤复发的原因。2015年，Yang等[15]在文献[14]基础上加入了脉冲免疫治疗，建立了一个与肿瘤-免疫及脉冲免疫治疗相关的数学模型，研究发现效应细胞的初始密度、效应细胞与肿瘤细胞的比例、免疫治疗周期等对癌症的治疗至关重要[15]。Parthasakha等[16]提出了具有Monod-Haldane动力学反应的肿瘤免疫与IL-2之间相互作用的时滞数学模型，以揭示相关细胞间现象的动力学机制。确定了模型解的正性、有界性和一致持久性，发现模型在无肿瘤平衡态下存在跨临界分支。Duan等[17]对高斯白噪声驱动的肿瘤免疫系统的稳定性进行了研究。研究发现在这个系统中，有几个稳态。当且仅当噪声较弱时，一个稳态是全局渐近稳定的，而当噪声较强或较弱时，另一个稳态总是不稳定的。Li等[18]提出并研究了随机切换环境下肿瘤细胞与免疫系统之间的竞争模型。得到了肿瘤细胞灭绝和持续存在的理论阈值。此外，还进行了随机模拟来证实和解释这些结论。

由于目前尚未有模型用Michaelis-Menten形式来表示肿瘤对效应细胞的抑制作用，因此在假设肿瘤细胞对免疫效应细胞具有线性刺激率和Michaelis-Menten型抑制函数的基础上，建立了如式(2)所示的肿瘤细胞与效应细胞相互作用的动力学模型:

(2)

分析了式(2)平衡点的存在性、局部渐近稳定性和全局动力学行为，以及模型产生鞍结点分支的条件，并进一步考虑肿瘤细胞对效应细胞的抑制率系数和效应细胞对肿瘤细胞的抑制率系数对式(2)动力学行为的影响。最后对所得结果进行生物学意义的解释和数值模拟来验证所得结果。

1 解的非负性和有界性

关于式(2)解的非负性和有界性有如下结论。

定理1设E(t)、T(t)为式(2)满足非负初始条件的解, 则E(t)、T(t)是非负的且最终有界。

证明首先，易知T=0为式(2)第二个方程的解，根据解的存在唯一性可知，在非负初始条件下T(t)≥0成立。当E(0)=0时，E′(0)>0，所以存在时刻t1，当t∈(0,t1)时，有E(t)>0；当E(0)>0时，假设存在t2>0使得E(t2)=0且当t∈(0,t2)时，有E(t)>0，所以E′(t2)≤0。而由式(2)的第一个式子可得E′(t2)=s+cT(t2)>0，这与E′(t2)≤0矛盾，所以，式(2)在非负初始条件下的解保持非负性。

其次，由式(2)的第二个方程得：

(3)

s+c(K+ε1)-μE

(4)

2 模型分析

2.1 平衡点的存在性

(5)

(6)

由式(6)的第二个方程可以得到：

(7)

当T≠K时，将式(7)代入式(6)的第一个方程得：

(8)

证明由于式(2)有正不变集D，所以可通过判定方程H(T)=0在区间(0,K)上根的存在性来确定的式(2)平衡点的存在性。直接计算可得:

(9)

式(9)中：H′(T)为H(T)关于T求一阶导数；H″(T)为H(T)关于T求二阶导数。

由于式(2)中各参数都为正数，所以在区间(0,K)内有H″(T)>0，即函数H(T)在区间(0,K)的图形是上凹的。

当方程H(T)=0在区间(0,K)内存在一个单根即T*时，其满足H′(T*)>0。

2.2 式(2)平衡点的局部稳定性

(10)

(11)

利用变换：

(12)

将式(11)化为

(13)

式(13)中：l=cαε-pr-pμ。

(14)

将u=av2+o(v2)代入式(14)，并比较两端同次幂系数可得:

(15)

此时:

(16)

进一步，将式(16)代入式(13)的第二个方程中，得到中心流形上的解满足：

(17)

综上所述，对式(2)边界平衡点的局部稳定性有如下结论。

(18)

因此

(19)

由式(6)的第一个方程可得：

(20)

因此：

(21)

进一步，由式(7)可得：

(22)

(23)

式(23)在原点附近可表示为

(24)

(25)

于是式(24)可写为

(26)

(27)

根据中心流形定理[19]，求得式(27)在原点的局部中心流形为

(28)

进一步，将式(28)代入式(27)的第一个方程，可得：

(29)

综上所述，对式(2)正平衡点的局部稳定性有如下结论。

2.3 式(2)平衡点的全局分析

对式(2)，记：

(30)

(31)

因此式(2)在区域D内无闭轨线。根据式(2)平衡点的局部稳定性，对于式(2)平衡点的全局稳定性有如下定理。

P0为无肿瘤平衡点图1 式(2)不存在正平衡点时的情形Fig.1 The case of formula (2) without the positive equilibrium

P*为模型的唯一有瘤平衡点图2 式(2)存在唯一正平衡点时的情形Fig.2 The case of formula (2) with unique positive equilibrium

为模型的有瘤平衡点图3 式(2)存在正平衡点P**时的情形Fig.3 The case of formula (2) with positive equilibrium P**

蓝色线表示P*的稳定流形图4 式(2)存在两个正平衡点时的情形Fig.4 The case of formula (2) with two positive equilibrium

3 数值模拟

为了直观地展示理论分析的正确性，对式(2)进行数值模拟。选取如下参数值：s=1.3×104，c=6.7×10-9，μ=0.0412，ε=1×107，r=1.5×10-1。

(2)取参数值p=5，α=1.041×10-7，K=570 000。此时R0=4.57>1，式(2)有唯一正平衡点P*(4.3×104,5.5×105)，且该点在区域D内是全局渐近稳定的(图2)。

将分析肿瘤细胞对效应细胞的抑制率系数p和效应细胞对肿瘤细胞的抑制率系数α对模型动力学行为的影响。

首先,当s=1.3×104，c=6.7×10-9，μ=0.041 2，ε=1×107，r=1.5×10-1，α=1.041×10-7，K=570 000时，取p分别为0、4、8、12、16。由图5可以看到，当p=0时，肿瘤细胞的数量处于最低水平，而效应细胞的数量处于最高水平。当p增大时，肿瘤细胞的数量会随着肿瘤对效应细胞的抑制率系数p的增大而增加，而效应细胞的数量会随着p的增大而减少。其次，当s=1.3×104，c=6.7×10-9，μ=0.041 2，ε=1×107，r=1.5×10-1，p=5，K=570 000时，由图6可以看到，随着效应细胞对肿瘤细胞的抑制率系数α的增大，肿瘤细胞的数量逐渐减少，而效应细胞的数量逐渐增加，当α=15×10-7时，肿瘤细胞的数量减少到0，也即此时肿瘤会灭绝。