肿瘤是指机体内某些组织的细胞发生异常增殖, 它有可能是良性的也可能是恶性的. 肿瘤在其增长过程中一直与宿主的免疫系统相互作用. 宿主的免疫系统是一个复杂的细胞、细胞因子、淋巴组织和器官网络, 可帮助宿主对抗抑制肿瘤细胞的生长和分裂^[1-6]. 另一方面, 肿瘤细胞可以通过其抗原来生化刺激免疫细胞(通常称为效应细胞)的产生, 同时又对活性免疫细胞具有一定灭活作用^{[7, 8]}. 目前, 肿瘤与免疫系统之间的作用机理还未被完全理解^[3]. 因此, 许多数学和医学工作者希望通过建立数学建模, 以揭示免疫系统和肿瘤细胞的作用机制, 为治疗和控制肿瘤和癌细胞的增长提供理论和临床指导^[9-15].

早在1994年, Kuznetsov等^[6]就提出了一类肿瘤-免疫相互作用的数学模型

其中$E = E(t)$和$T = T(t)$分别是效应细胞和肿瘤细胞在时间$t$的数量; $s$是效应细胞流入肿瘤部位的正常速率; 具有饱和作用的Michaelis-Menten形式$\frac{cET}{\varepsilon+T}$表示抗原性; $\beta$是肿瘤细胞使效应细胞失活的速率系数; $\mu$是效应细胞自然死亡率系数; 假设在没有效应细胞的情况下, 肿瘤细胞服从logistic增长模型, 内禀增长率为$r$, 肿瘤细胞的最大承载量为$K$; $\alpha$为效应细胞杀死肿瘤细胞的速率系数. Kuznetsov等^[6]主要通过数据分析研究模型(1.1), 揭示了肿瘤刺激免疫细胞的生长、肿瘤免疫逃逸以及肿瘤休眠等一些临床现象. 但由于模型(1.1)会存在三个正平衡点, 而使得一些理论分析难以实现.

基于模型(1.1), 一些数学模型相继提出, 并得到研究^[16-19], 其中Galach^[20]以双线性形式$\theta ET$来替代模型(1.1)中效应细胞的饱和形式刺激增长率, 得到模型

其中$k = \beta-\theta$, 并对其局部动力学性态进行分析, 仅得到一些充分条件.

随后, 黄佩等^[10]对该模型进行了完整的数学分析, 得到其全局动力学性态. 在模型(1.2)中, 没有明确体现肿瘤对免疫系统的抗原性作用, 即肿瘤细胞的存在会刺激效应细胞的增殖. 而在Kirschner和Panetta^[2]建立的模型中用$cT$刻画了效应细胞在肿瘤抗原性作用下的增殖速率. 因此, 将这一抗原性作用形式引入模型(1.2), 便得到模型

本文定性地分析了此模型的全局动力学性态, 发现该模型会出现鞍结点分支以及强有瘤平衡点和无瘤平衡点双稳定的现象. 同时, 进一步的分析发现了肿瘤细胞抗原性的变化会使该模型出现复杂的动力学性态并通过数值模拟进行了展示.

为了便于数学分析, 对(1.3)式作无量纲变换

并记$\eta = \frac{cK}{s}$, $\delta = \frac{\beta K}{\mu}$, $q = \frac{\alpha s}{\mu^2}$, $r = \frac{\mu r}{\alpha s}$, 同时将$\tau$仍写为$t$, 那么模型(1.3)变为

显然$y = 0$是系统(2.1)的一条解轨线, 于是系统(2.1)满足初值条件$y(0) = y_0>0$的解轨线不会与$y = 0$相交, 即当$t>0$时, $y = y(t)>0$恒成立. 因为当$x = 0$时, $\left.\frac{{\rm d}x}{{\rm d}t}\right|_{x = 0} = 1+\eta y>0$, 所以系统(2.1)满足初值条件$x(0) = x_0>0$的解恒有$x = x(t)>0(t>0)$.

对任意的$x>0$, $y\geq0$, 由系统(2.1)的第二个方程有

则${ }\limsup\limits_{t\to\infty} y(t)\leq 1$. 因此对任意的$\varepsilon>0$, 存在$t_1>0$, 使得当$t>t_1$时, $y(t)<1+\varepsilon$. 于是当$t>t_1$时, 由系统(2.1)的第一式有

所以${ }\limsup\limits_{t\to\infty} x(t)\leq 1+\eta(1+\varepsilon)$. 由$\varepsilon$的任意性可得${ }\limsup\limits_{t\to\infty} x(t)\leq 1+\eta$. 因此系统(2.1)有正不变集

即系统(2.1)是有界的. 以下将在$D$上研究系统(2.1)动力学性态.

令系统(2.1)中的$\frac{{\rm d}x}{{\rm d}t} = 0$, $\frac{{\rm d}y}{{\rm d}t} = 0$, 得

由系统(3.1)的第二式得$y = 0$或

将$y = 0$代入系统(3.1)的第一式可得$x = 1$, 故系统(2.1)存在平衡点$P_0(1, 0)$. 因为$P_0$的$y$坐标等于$0$, 所以称其为无瘤平衡点.

当$y\neq0$时, 将(3.2)式代入系统(3.1)的第一个方程得

由于$y\geq1$时, $f(y)>0$, 所以$f(y) = 0$的正解仅在区间$(0, 1)$内. 对于一元二次函数$f(y)$, $f(0) = 1-r$, $f(1) = \eta+1>0$且$f'(y) = 2\delta ry-(\delta r-r-\eta)$.

记方程$f(y) = 0$的判别式为$\Delta = (\delta r-r-\eta)^2-4\delta r(1-r)$, 则当$r>1$时, 由$f(0) = 1-r<0$知$f(y) = 0$在$(0, 1)$内有唯一正解

当$r = 1$且$\eta<(\delta-1)r$时, $f(y) = 0$在$(0, 1)$内只有一个正解

当$r<1$, $\eta<(\delta-1)r$且$\Delta>0$时, $f(y) = 0$在$(0, 1)$内有两个正解

当$r<1$, $\eta<(\delta-1)r$且$\Delta = 0$时, $f(y) = 0$在$(0, 1)$内有一个二重正解

进一步, 相应地将$y_j(j = 1, \cdots, 5)$代入(3.2)式, 可得

于是, 系统(2.1)在相应的条件下有正平衡点$P_j(x_j, y_j)$ (称其为有瘤平衡点). 对于$P_j$, 根据$y_j$与$f(y)$的关系, 易知

注意到条件

等价于条件

因此, 关于系统(2.1)的平衡点存在性, 有如下结论.

定理3.1  系统(2.1)总存在无瘤平衡点$P_0(1, 0)$. 以下为系统(2.1)有瘤平衡点的存在情况:

当$r>1$时, 系统(2.1)存在唯一有瘤平衡点$P_1(x_1, y_1)$;

当$r = 1$, $\delta>1$且$\eta<\delta-1$时, 系统(2.1)存在唯一有瘤平衡点$P_2(x_2, y_2)$;

当$\frac{4\delta}{(\delta+1)^2}<r<1$, $\delta>1$且$\eta<(\delta-1)r-2\sqrt{\delta r(1-r)}$时, 系统(2.1)存在两个不同的有瘤平衡点$P_3(x_3, y_3)$ (称为弱有瘤平衡点)和$P_4(x_4, y_4)$ (称为强有瘤平衡点);

当$\frac{4\delta}{(\delta+1)^2}<r<1$, $\delta>1$且$\eta = (\delta-1)r-2\sqrt{\delta r(1-r)}$时, 系统(2.1)存在唯一有瘤平衡点$P_5(x_5, y_5)$.

根据函数$f(y)$的表达式和定理3.1, 对于情形$\delta\le 1$, 当$r>1$时系统(2.1)有唯一的有瘤平衡点$P_1$, 当$r\le 1$时系统(2.1)没有有瘤平衡点; 对于情形$\delta> 1$, 系统(2.1)的平衡点存在条件相对复杂些. 为了直观地显示这些条件, 借助$r$的函数$\eta = (\delta-1)r-2\sqrt{\delta r(1-r)}$, 在平面$(r, \eta)$上用图 1显示出相应的平衡点存在区域, 其中

由定理3.1可知: 当$(r, \eta)\in D_0$时, 系统(2.1)不存在有瘤平衡点; 当$(r, \eta)\in D_1$时, 系统(2.1)存在一个有瘤平衡点$P_1$或$P_2$; 当$(r, \eta)\in D_2$时, 系统(2.1)存在两个有瘤平衡点$P_3$和$P_4$; 当$(r, \eta)\in D_3$时, 系统(2.1)存在一个二重有瘤平衡点$P_5$.

同时, 由图 1易知, 系统(2.1)有瘤平衡点的存在性对参数$\eta$ (抗原性)的依赖性. 当$\eta\geq\delta-1$时, 随着$r$的增加并通过$1$时, 系统(2.1)的有瘤平衡点由没有变为一个. 而当$\delta>1$且$\eta<\delta-1$时, 随着$r$的增加, 系统(2.1)的有瘤平衡点由没有, 变为两个, 再变为一个. 此时系统(2.1)在由无瘤平衡点变为两个有瘤平衡点时经历了鞍结点分支.

在本节, 将首先考虑系统(2.1)无瘤平衡点的局部稳定性, 然后考虑有瘤平衡点的局部稳定性, 最后分析稳定有瘤平衡点的类型.

系统(2.1)在无瘤平衡点$P_0$处的Jacobi矩阵为

于是它的两个特征值分别为$\lambda_1 = -1$, $\lambda_2 = q(r-1)$. 因此$P_0$当$r<1$时是局部渐近稳定的, 当$r>1$时是鞍点. 当$r = 1$时, $P_0$是高阶平衡点.

下面分$\eta\ne\delta$和$\eta = \delta$两种情形来分析$r = 1$时系统(2.1)在$P_0$附近的动力学性态.

当$\eta\ne\delta$时, 将平衡点$P_0$平移到原点, 对系统(2.1)作平移变换, $u = x-1$, $v = y$, 得

系统(4.1)在$P_0$处的线性化系统为

其系数矩阵的特征值为$-1$和$0$. 对应于特征值$-1$和$0$的特征向量分别为$(1, 0)^T$和$(\eta-\delta, 1)^T$. 于是对系统(4.1)作非退化线性变换

系统(4.1)变为

其中

根据中心流形定理^[14], 设系统(4.2)在原点处的局部中心流形为

则函数$h(w)$满足

将$G(z, w), H(z, w)$和$h(w)$代入上式, 可得

所以系统(4.2)在原点处的中心流形为

再将其代入(4.2)的第二个方程有

当$\delta-\eta = 1$时, (4.3)式变为

此时$P_0$为局部渐近稳定的结点. 当$\delta-\eta\ne 0, 1$时, $P_0$为鞍结点. 进一步, 根据前面所做的变换, 可知当$\delta-\eta<1$时, $P_0$为$D$上局部渐近稳定的结点; 当$\delta-\eta>1$时, $P_0$为$D$上的鞍点.

当$r = 1$且$\eta = \delta$时, 系统(2.1)变为

定义函数$V = \frac{q(x-1)^2}{8}+y$, 则$V$沿着系统(4.4)解的全导数为

则当$y\ge0$时, $\frac{{\rm d}V}{{\rm d}t}\le 0$. 因此, 当$r = 1$且$\eta = \delta$时, $P_0$是稳定的.

综上所述, 对于系统(2.1)的无瘤平衡点$P_0$的局部稳定性有如下结论.

定理4.1  当$r<1$或$r = 1$且$r\geq\delta-\eta$时, 无瘤平衡点$P_0$局部渐近稳定; 当$r>1$或$r = 1$且$r<\delta-\eta$时, 无瘤平衡点$P_0$不稳定.

对于系统(2.1)的任意有瘤平衡点$P^*(x^*, y^*)$, 系统(2.1)在$P^*$处的Jacobi矩阵为

因为$x^*$, $y^*$满足等式(3.2), 所以系统(2.1)在$P^*$处的Jacobi矩阵可以写为

其迹为tr$J_1(P^*) = -(1+\delta y^*+qry^*)<0$, 行列式为det$J_1(P^*) = qy^*[2\delta ry^*-(\delta r-r-\eta)] = qy^*f'(y^*)$.

由$f'(y_3)<0$, 即det$J_1(P_3)<0$, 所以弱有瘤平衡点$P_3$是鞍点.

对于平衡点$P_i(i = 1, 2, 4)$, $f'(y_i)>0$, 即det$J_1(P_i)>0$, 所以$P_i$是局部渐近稳定的.

又$f'(y_5) = 0$意味着det$J (P_5) = 0$, 故$P_5$是高阶平衡点. 此时, $J(P_5)$的特征值分别为$\lambda_1 = -(1+\delta y_5+qry_5)$和$\lambda_2 = 0$. 其特征值对应的特征向量分别为$(1, A)^T$和$(-r, 1)^T$, 其中$A = \frac{-qry_5}{\eta-\delta x_5} = \frac{q(\delta r-r-\eta)}{\delta(\delta r+r-\eta)}>0$.

为了研究系统(2.1)在有瘤平衡点$P_5$附近的动力学性态, 对系统(2.1)进行平移变换, 令$u_1 = x-x_5$, $v_1 = y-y_5$, 则

其中

再对(4.5)式进行非退化线性变换

则(4.5)式变为

其中

类似于分析$r = 1$且$\eta\neq\delta$时无瘤平衡点$P_0$稳定性的方法, 计算得系统(4.6)在原点的局部中心流形为

将其代入系统(4.6)的第二个式子得

故$P_5$为鞍结点.

通过以上分析得到如下定理.

定理4.2  若系统(2.1)的有瘤平衡点$P_j(j = 1, \cdots, 5)$在$D$内存在, 则$P_3$是鞍点, $P_i(i = 1, 2, 4)$是局部渐近稳定的, $P_5$是鞍结点.

平面系统的渐近稳定平衡点有焦点和结点两种类型. 在这两类平衡点的附近, 系统轨线收敛于该平衡点的方式有所不同. 若该平衡点为焦点, 则当$t$趋于无穷大时, 由该平衡点附近出发的解轨线将以阻尼振荡的形式收敛于该平衡点; 若其为结点, 则在其附近出发的解轨线会最终单调趋于平衡点. 为了区分系统(2.1)在稳定有瘤平衡点附近的变化特征, 在此进一步研究这类平衡点的类型.

定理4.2已经显示有瘤平衡点$P_i(i = 1, 2, 4)$是局部渐近稳定的, 那么$P_i$有可能是焦点或结点. 对于系统(2.1)在有瘤平衡点$P_i(x_i, y_i)$处的$J(P_i)$, 当$\Delta_{i} = {\rm tr}J (P_i)^2-4{\rm det}J (P_i)<0$时, $P_i$是焦点; 当$\Delta_{i}\ge 0$时, $P_i$是结点.

注意到平衡点$P_i(i = 1, 2, 4)$的坐标$x_i$, $y_i$与$q$无关, 所以可将$\Delta_i$表示为$q$的一元二次函数, 即

进一步, 需确定函数$F_i(q)$的判别式为

的符号.

由$y_i (i = 1, 2, 4)$满足的方程(3.3)可知$1 = r+(\delta r-r-\eta)y_i-\delta r y_i^2$, 那么

有$\delta ry_i+\eta-\delta r = \frac{\eta-\delta}{1+\delta y_i}$, 因此$\Delta_i^* = 16y_i^2f'(y_i)\frac{\eta-\delta}{1+\delta y_i}$. 又$f'(y_i)>0$, 所以$\Delta_i^*$的符号由$\eta-\delta$判定.

注意到$\eta<\delta$是$P_2$和$P_4$存在的必要条件, 因此$\Delta_k^*<0$$(k = 2, 4)$, 进而有$\Delta_{k}>0$$(k = 2, 4)$, 所以$P_k$$(k = 2, 4)$为结点, 即平衡点$P_2$和$P_4$只要存在, 就是结点.

对于平衡点$P_1$, 当$\eta\le\delta$时, $\Delta_1^*\le 0$也意味着是结点. 又当$3\delta r y_1\le (2\delta r-r-2\eta)$, 即$y_1\le\frac{(2\delta r-r-2\eta)}{3\delta r}$时, 对于任意的$q$, 都有$F_i(q)>0$, 所以此时$P_1$也是结点.

当$\eta>\delta$时, $\Delta_1^*>0$, 因此$F_i(q) = 0$有两个互异同号实根$q_1$和$q_2(q_1<q_2)$. 进一步, $y_1>\frac{(2\delta r-r-2\eta)}{3\delta r}$意味着$q_1$和$q_2$是正的, 则当$q_1<q<q_2$时, $\Delta_{1}<0$. 此时$P_1$即为焦点. 而当$q\le q_1$或$q\geq_2$时, $P_1$为结点.

根据以上推理和定理4.2可得如下结论.

定理4.3  系统(2.1)的有瘤平衡点$P_2$和$P_4$只要存在, 就是结点.

当条件

成立时, $P_1$存在就是焦点, 否则是结点, 其中$q_1$和$q_2$是方程$F_i(q) = 0$的根.

选择$B = \frac{1}{y}$作为Dulac函数, 记系统(2.1) 右端的函数分别为$M(x, y)$, $N(x, y)$, 则

因此, 根据Bendixson-Dulac定理^[14]可知系统(2.1)不存在闭轨线, 即系统(2.1)没有周期解.

再结合上节所得平衡点的局部稳定性结果, 可得如下关于系统(2.1)的全局性态结论.

定理5.1  对于系统(2.1),

$\rm (1)$若无瘤平衡点存在, 无瘤平衡点$P_0$在区域$D$上是全局渐近稳定的(图 2(a)).

$\rm (2)$若存在唯一一个有瘤平衡点$P_1$或$P_2$时, 该平衡点在区域$D$内是全局渐近稳定的(图 2(b)和(c)分别对应焦点和结点的情形).

$\rm (3)$若存在两个互异的有瘤平衡点$P_3$和$P_4$时, 则存在弱有瘤平衡点$P_3$的稳定流形将可行区域(第一象限)分为两部分(图 2(d)). 若初始值在下方区域, 则肿瘤最终消亡; 若初始值在上方区域, 则肿瘤最终向强有瘤平衡点$P_4$发展.

$\rm (4)$若存在有瘤平衡点$P_5$时, 则其稳定流形将可行区域(第一象限)分为两部分(图 2(e)), 使得当初始值在下方区域时, 肿瘤最终消亡; 当初始值在上方区域时, 肿瘤最终向有瘤平衡点$P_5$发展.

本文在考虑肿瘤抗原性的作用下, 提出了一类具有抗原性的肿瘤-免疫系统相互作用的模型. 为了方便分析, 对所提模型首先进行了无量纲变换, 以简化模型. 在得到模型有瘤平衡点的存在条件基础上, 通过分析所有平衡点的局部动力学性态和构造Dulac函数排除系统周期解的存在性, 得到了模型的全局动力学性态. 分析结果表明, 抗原性的引入会使模型发生鞍结点分支, 以及强有瘤平衡点与无瘤平衡点同时稳定的现象. 这种双稳态的发生意味着肿瘤增长的发展结局依赖于它们的初始状态(如图 2(d)和(e)所示). 另外, 本文进一步给出了有瘤平衡点稳定时是焦点还是结点的判定条件, 为临床判定肿瘤最终发展阶段的性态提供了理论指导依据. 根据所得结果来看, 本文所提模型的动力学性态具有一定的复杂性, 在一定程度上反映了肿瘤细胞和免疫细胞的作用机制, 但还未能充分展示一些更为复杂的动力学现象(如, Hopf分支、周期解的存在性和B-T分支等). 这有待于对所建模型进行进一步完善.