分支理论在肿瘤增长模型中的应用

Application of Bifurcation Theory in Tumor Growth Model

摘要: 本文运用中心流形定理和局部分支理论研究了一个四维肿瘤增长模型的分支动态，理论证明了该模型Transcritial分支和Hopf分支的存在性。本文研究结果表明，若能通过治疗抑制肿瘤细胞的增长率到某个临界值以下，肿瘤细胞数目会在一段时间内较快地衰减至0，人体免疫系统能杀死肿瘤细胞，病人会获得痊愈。

Abstract: The paper does qualitative research on a four-dimensional tumor growth model by using center manifolds theory and local bifurcation theory and proves the existence of Transcritical bifurcation and Hopf bifurcation theoretically. It turns out that the patient will recover from the tumor within some time by self-immune system if the tumor growth rate could drop below the threshold with the treatment.

1. 引言

肿瘤增长威胁到人类的健康甚至生命，所以需要对其进行严格的理论分析。伴随着科技进步，分析手段逐渐多样化，数学理论分析也可应用到肿瘤的研究中，肿瘤增长模型便是分支理论在生物医学上应用的一个典例。

肿瘤增长是一个复杂的过程，受各类细胞之间的相互作用影响 [1]。建立肿瘤增长模型的目的是能找到肿瘤增长过程中的关键影响因素。1994年，Kuznetsov [2] 提出了一个关于肿瘤细胞和效应免疫细胞之间关系的二维模型，为研究肿瘤增长模型奠定了基础。2003年，De pillis和Radunskaya [3] 提出了一个关于肿瘤细胞，健康宿主细胞和效应免疫细胞之间关系的三维模型，该模型选取的参数符合生物特征，即具有实际意义，且适用于一般类型的肿瘤研究。他们发现当轨线进入无肿瘤平衡点的吸引域时化疗就可以停止了。2014年，Louise Viger [4] 在研究肿瘤增长模型时，考虑到了内皮细胞是肿瘤细胞发生转移的关键，因为内皮细胞会促进新的血管生成，给肿瘤细胞的繁殖提供养分，同时也给肿瘤细胞的扩散开辟了路径，所以内皮细胞是肿瘤增长并发生转移的一个重要因素。Viger [4] 基于De pillis和Radunskaya [3] 提出的三维模型，引入内皮细胞的数量，建立了一个新的四维模型。他发现当免疫细胞数量过少或者内皮细胞数量过多，就会导致单肿瘤部位的肿瘤细胞数目过多，进而开始向其他部位转移，即会发生癌细胞转移；同时内皮细胞对免疫细胞数量也具有一定促进作用。

2. 理论分析

2.1. 模型

文献 [4] 中给出的模型为：

$\{\begin{array}{l}\stackrel{\dot{}}{x}={\rho }_{1}x\left(1-x\right)-{\alpha }_{13}xz\\ \stackrel{\dot{}}{y}=\frac{{\rho }_{2}yz}{1+z}-{\alpha }_{23}yz-{\delta }_{2}y+{\alpha }_{24}yw\\ \stackrel{\dot{}}{z}={\rho }_{3}z\left(1-z\right)-{\alpha }_{31}zx-{\alpha }_{32}zy+\frac{{\alpha }_{34}zw}{1+w}\\ \stackrel{\dot{}}{w}=\frac{{\rho }_{4}wz}{1+z}-{\delta }_{4}w\end{array}$ (1)

其中， $x,y,z,\omega$ 分别表示健康宿主细胞的数量，效应免疫细胞的数量，肿瘤细胞的数量和内皮细胞的数量。 ${\rho }_1$ 表示健康宿主细胞的增长率， ${\rho }_2$ 表示效应免疫细胞的增长率， ${\rho }_3$ 表示肿瘤细胞的增长率， ${\rho }_4$ 表示内皮细胞的增长率， ${\alpha }_{13}$ 表示宿主细胞被肿瘤细胞杀死率， ${\alpha }_{23}$ 表示效应免疫细胞被肿瘤细胞杀死率， ${\alpha }_{31}$ 表示肿瘤细胞被宿主细胞杀死率， ${\alpha }_{32}$ 表示肿瘤细胞被效应免疫细胞杀死率， ${\alpha }_{24}$ 表示效应免疫细胞受内皮细胞影响的增长率， ${\alpha }_{34}$ 表示肿瘤细胞受内皮细胞影响的增长率， ${\delta }_2$ 表示效应免疫细胞的自然死亡率， ${\delta }_4$ 表示内皮细胞的自然死亡率。

本文对模型(1)进行研究。

2.2. 平衡点的存在性

首先求解平衡点的参数表达式，令：

$\{\begin{array}{l}{\rho }_1{x}_1\left(1-x_1\right)-{\alpha }_{13}{x}_1{x}_3=0\\ \frac{{\rho }_2{x}_2{x}_3}{1+x_3}-{\alpha }_{23}{x}_2{x}_3-{\delta }_2{x}_2+{\alpha }_{24}{x}_2{x}_4=0\\ {\rho }_3{x}_3\left(1-x_3\right)-{\alpha }_{31}{x}_3{x}_1-{\alpha }_{32}{x}_3{x}_2+\frac{{\alpha }_{34}{x}_3{x}_4}{1+x_4}=0\\ \frac{{\rho }_4{x}_4{x}_3}{1+x_3}-{\delta }_4{x}_4=0\end{array}$ (2)

解得全部平衡点的参数表达式为：

$P_0\left(0,0,0,0\right)$, $P_1\left(0,0,1,0\right)$, $P_2\left(1,0,0,0\right)$, $P_3\left(0,0,b,\frac{{\rho }_3\left(b-1\right)}{{\alpha }_{34}-{\rho }_3\left(b-1\right)}\right)$,

$P_4\left(0,\frac{{\rho }_3}{{\alpha }_{32}}\left(1-c_1\right),c_1,0\right)$, $P_5\left(0,\frac{{\rho }_3}{{\alpha }_{32}}\left(1-c_2\right),c_2,0\right)$, $P_6\left(\frac{{\rho }_3\left({\alpha }_{13}-{\rho }_1\right)}{{\alpha }_{31}{\alpha }_{13}-{\rho }_1{\rho }_3},0,\frac{{\rho }_1\left({\alpha }_{31}-{\rho }_3\right)}{{\alpha }_{31}{\alpha }_{13}-{\rho }_1{\rho }_3},0\right)$,

$P_7\left(0,\frac{1}{{\alpha }_{32}}\left({\rho }_3\left(1-b\right)+\frac{{\alpha }_{34}d}{1+d}\right),b,d\right)$, $P_8\left(1-\frac{{\alpha }_{13}}{{\rho }_1}b,0,b,\frac{{\alpha }_{31}\left(1-\frac{{\alpha }_{13}}{{\rho }_1}b\right)-{\rho }_3\left(1-b\right)}{{\rho }_3\left(1-b\right)-{\alpha }_{31}\left(1-\frac{{\alpha }_{13}}{{\rho }_1}b\right)}\right)$,

$P_9\left(1-\frac{{\alpha }_{13}}{{\rho }_1}{c}_1,\frac{1}{{\alpha }_{32}}\left[{\rho }_3\left(1-c_1\right)-{\alpha }_{31}\left(1-\frac{{\alpha }_{13}}{{\rho }_1}{c}_1\right)\right],c_1,0\right)$,

$P_{10}\left(1-\frac{{\alpha }_{13}}{{\rho }_1}{c}_2,\frac{1}{{\alpha }_{32}}\left[{\rho }_3\left(1-c_2\right)-{\alpha }_{31}\left(1-\frac{{\alpha }_{13}}{{\rho }_1}{c}_2\right)\right],c_2,0\right)$,

$P_{11}\left(1-\frac{{\alpha }_{13}}{{\rho }_1}b,\frac{{\rho }_3}{{\alpha }_{32}}\left(1-b\right)-{\alpha }_{31}\left(1-\frac{{\alpha }_{13}}{{\rho }_1}b\right)+\frac{{\alpha }_{34}d}{1+d},b,d\right)$.

其中： $b=\frac{{\delta }_4}{{\rho }_4-{\delta }_4}$， $d=\frac{1}{{\alpha }_{24}}\left({\alpha }_{23}b+{\delta }_2-\frac{{\rho }_2{\delta }_4}{{\rho }_4}\right)$，

$c_1=\frac{{\rho }_2-{\alpha }_{23}-{\delta }_2-\sqrt{{\left({\rho }_2-{\alpha }_{23}-{\delta }_2\right)}^2-4{\alpha }_{23}{\delta }_2}}{2{\alpha }_{23}}$ ;

$c_2=\frac{{\rho }_2-{\alpha }_{23}-{\delta }_2+\sqrt{{\left({\rho }_2-{\alpha }_{23}-{\delta }_2\right)}^2-4{\alpha }_{23}{\delta }_2}}{2{\alpha }_{23}}$.

文献 [4] 中给出了给定参数下具有实际意义(细胞的数目没有出现负的和全为0的情况)的6个平衡点的参数表达式： $P_0,P_1,P_2,P_4,P_9,P_{11}$。

2.3. 平衡点的个数，类型和稳定性

本文选取文献 [4] 中没有研究过的参数 ${\rho }_3$ 为分支参数，其他参数取值固定，来源于文献 [4]：

$\begin{array}{l}{\rho }_1=0.518,{\rho }_2=4.5,{\rho }_4=0.86,{\delta }_2=0.5,{\delta }_4=1/11,{\alpha }_{13}=1.5,\\ {\alpha }_{23}=0.2,{\alpha }_{24}=0.3,{\alpha }_{31}=1,{\alpha }_{32}=2.5,{\alpha }_{34}=0.75\end{array}$

同时为了方便后面计算和表示，令动态变量 $x,y,z,\omega$ 分别表示为 $x_1,x_2,x_3,x_4$，得到系统(3)

$\{\begin{array}{l}{\stackrel{\dot{}}{x}}_1=0.518x_1\left(1-x_1\right)-1.5x_1{x}_3\\ {\stackrel{\dot{}}{x}}_2=\frac{4.5x_2{x}_3}{1+x_3}-0.2x_2{x}_3-0.5x_2+0.3x_2{x}_4\\ {\stackrel{\dot{}}{x}}_3={\rho }_3{x}_3\left(1-x_3\right)-x_3{x}_1-2.5x_3{x}_2+\frac{0.75x_3{x}_4}{1+x_4}\\ {\stackrel{\dot{}}{x}}_4=\frac{0.86x_4{x}_3}{1+x_3}-\frac{1}{11}{x}_4\end{array}$ (3)

考虑参数的实际意义，选取 ${\rho }_3$ 大于0；考虑生物相关性，排除细胞数出现负的和细胞数全为0的平衡点。随着分支参数 ${\rho }_3$ 的变化，系统(3)的平衡点个数，类型及稳定性如下：

1) 当 ${\rho }_3\in \left(0,0.09104\right)$，系统有8个平衡点，其中1个sink $P_2$ 和7个saddles $P_i\left(i=1,4,6,7,8,9,11\right)$。

2) 当 ${\rho }_3={\rho }_3^1=0.09104$，系统有8个平衡点，其中1个sink $P_2$，6个saddles $P_i\left(i=1,4,6,8,9,11\right)$ 和1个特征值为 $\left\{-0.01076,\pm 0.25794i,0.3407\right\}$ 的非双曲平衡点 $P_7=\left(x_1^1,x_2^1,x_3^1,x_4^1\right)=\left(0,0.07346,0.1182,0.1598\right)$。

3) 当 ${\rho }_3\in \left(0.09104,0.62866\right)$，系统有8个平衡点，其中1个sink $P_2$ 和7个saddles $P_i\left(i=1,4,6,7,8,9,11\right)$。

4) 当 ${\rho }_3={\rho }_3^2=0.62866$，系统有7个平衡点，其中1个sink $P_2$，5个saddles $P_i\left(i=1,4,6,7,9\right)$ 和1个特征值为 $\left\{-0.5781,0.186,-0.023,0\right\}$ 的非双曲平衡点 $P_8=P_{11}=\left(x_1^2,x_2^2,x_3^2,x_4^2\right)=\left(0.6577,0,0.1182,0.1598\right)$。

5) 当 ${\rho }_3\in \left(0.62866,0.71044\right)$，系统有8个平衡点，其中1个sink $P_2$ 和7个saddles $P_i\left(i=1,4,6,7,8,9,11\right)$。

6) 当 ${\rho }_3={\rho }_3^3=0.71044$，系统有7个平衡点，其中1个sink $P_2$，5个saddles $P_i\left(i=1,4,7,8,11\right)$ 和1 个特征值为 $\left\{-0.5743,0.161,0,0.0097\right\}$ 的非双曲平衡点 $P_6=P_9=\left(x_1^3,x_2^3,x_3^3,x_4^3\right)=\left(0.6163,0,0.1325,0\right)$。

7) 当 ${\rho }_3\in \left(0.71044,0.74588\right)$，系统有8个平衡点，其中1个sink $P_2$ 和7个saddles $P_i\left(i=1,4,6,7,8,9,11\right)$。

8) 当 ${\rho }_3={\rho }_3^4=0.74588$，系统有7个平衡点，其中1个sink $P_2$，5个saddles $P_i\left(i=1,4,7,9,11\right)$ 和1个特征值为 $\left\{-0.5785,0.1497,-0.048,0\right\}$ 的非双曲平衡 $P_6=P_8=\left(x_1^4,x_2^4,x_3^4,x_4^4\right)=\left(0.6577,0,0.1182,0\right)$。

9) 当 ${\rho }_3\in \left(0.74588,1\right)$，系统有8个平衡点，其中1个sink $P_2$ 和7个saddles $P_i\left(i=1,4,6,7,8,9,11\right)$。

10) 当 ${\rho }_3={\rho }_3^5=1$，系统有7个平衡点，其中6个saddles $P_i\left(i=1,4,7,8,9,11\right)$ 和1个特征值为 $\left\{-0.518,-0.5,0,-0.0909\right\}$ 的非双曲平衡点 $P_2=P_6=\left(x_1^5,x_2^5,x_3^5,x_4^5\right)=\left(1,0,0,0\right)$。

11) 当 ${\rho }_3\in \left(1,1.24771\right)$，系统有8个平衡点，其中1个sink $P_6$ 和7个saddles $P_i\left(i=1,2,4,7,8,9,11\right)$。

12) 当 ${\rho }_3={\rho }_3^6=1.24771$，系统有8个平衡点，其中1个sink $P_6$，6个saddles $P_i\left(i=1,2,4,7,8,11\right)$ 和一个特征值为 $\left\{-0.4846,\pm 0.3669i,0.0097\right\}$ 的非双曲平衡点 $P_9=\left(x_1^6,x_2^6,x_3^6,x_4^6\right)=\left(0.6163,0.1864,0.1325,0\right)$。

13) 当 ${\rho }_3\in \left(1.24771,1.31131\right)$，系统有8个平衡点，其中1个sink $P_6$，7个saddles $P_i\left(i=1,2,4,7,8,9,11\right)$。

14) 当 ${\rho }_3={\rho }_3^7=1.31131$，系统有8个平衡点，其中1个sink $P_6$，6个saddles $P_i\left(i=1,2,4,7,8,9\right)$ 和1个特征值为 $\left\{-0.4858,\pm 0.4074i,-0.0099\right\}$ 的非双曲平衡点 $P_{11}=\left(x_1^7,x_2^7,x_3^7,x_4^7\right)=\left(0.6577,0.2408,0.1182,0.1598\right)$。

15) 当 ${\rho }_3\in \left(1.31131,+\infty \right)$，系统有8个平衡点，其中2个sinks $P_i\left(i=6,11\right)$ 和6个saddles $P_i\left(i=1,2,4,7,8,9\right)$。

2.4. 平衡点的静态分支

我们注意到当 ${\rho }_3={\rho }_3^2,{\rho }_3^3,{\rho }_3^4,{\rho }_3^5$ 时，系统(3)存在的非双曲平衡点的特征根有单零根，此时对应的平衡点 $P_{11},P_6,P_8,P_2$ 可能会发生静态分支，我们将应用中心流形定理研究当分支参数 ${\rho }_3$ 分别取 ${\rho }_3^2,{\rho }_3^3,{\rho }_3^4,{\rho }_3^5$ 时系统(3)在平衡点 $P_{11},P_6,P_8$ 和 $P_2$ 附近的动态变化。

情形1 ${\rho }_3={\rho }_3^5=1$，平衡点 $P_2=\left(1,0,0,0\right)$，对应的特征值为 $\left(-0.518,-0.5,0,-0.0909\right)$。首先通过坐标变换 $y_1=x_1-x_1^5,y_2=x_2-x_2^5,y_3=x_3-x_3^5,y_4=x_4-x_4^5,\mu ={\rho }_3-{\rho }_3^5$ 将平衡点 $P_2$ 平移到坐标原点，这样模型(3)就变为：

$\{\begin{array}{l}{\stackrel{\dot{}}{y}}_1=-0.518y_1\left(1+y_1\right)-1.5\left(y_1+1\right)y_3\\ {\stackrel{\dot{}}{y}}_2=\frac{4.5y_2{y}_3}{1+y_3}-0.2y_2{y}_3-0.5y_2+0.3y_2{y}_4\\ {\stackrel{\dot{}}{y}}_3=\left(\mu +1\right)y_3\left(1-y_3\right)-y_3\left(y_1+1\right)-2.5y_3{y}_2+\frac{0.75y_3{y}_4}{1+y_4}\\ {\stackrel{\dot{}}{y}}_4=\frac{0.86y_4{y}_3}{1+y_3}-\frac{1}{11}{y}_4\end{array}$ (4)

选用一个线性变换矩阵T (T由系统(4)的平衡点原点的特征值对应的特征向量组成)

$T=\left(\begin{array}{cccc}1& 0& 0& -0.9452\\ 0& 1& 0& 0\\ 0& 0& 0& 0.3264\\ 0& 0& 1& 0\end{array}\right)$，令 $\left(\begin{array}{c}{y}_1\\ y_2\\ y_3\\ y_4\end{array}\right)=T\left(\begin{array}{c}{z}_1\\ z_2\\ z_3\\ z_4\end{array}\right)$，将系统（4）变换为：

$\left(\begin{array}{c}{\stackrel{\dot{}}{z}}_1\\ {\stackrel{\dot{}}{z}}_2\\ {\stackrel{\dot{}}{z}}_3\\ {\stackrel{\dot{}}{z}}_4\end{array}\right)=\left(\begin{array}{cccc}-0.518& 0& 0& 0\\ 0& -0.5& 0& 0\\ 0& 0& -\text{0}\text{.0909}& 0\\ 0& 0& 0& 0\end{array}\right)\left(\begin{array}{c}{z}_1\\ z_2\\ z_3\\ z_4\end{array}\right)+\left(\begin{array}{c}{f}_1\left(z_1,z_2,z_3,z_4,\mu \right)\\ f_2\left(z_1,z_2,z_3,z_4,\mu \right)\\ f_3\left(z_1,z_2,z_3,z_4,\mu \right)\\ f_4\left(z_1,z_2,z_3,z_4,\mu \right)\end{array}\right)$ (5)

其中：

$\begin{array}{c}{f}_1\left(z_1,z_2,z_3,z_4,\mu \right)=-0.518z_1^2-0.456z_1{z}_4-2.363z_2{z}_4+0.585z_4^2\\ +0.945\mu z_4-0.309\mu z_4^2+0.709z_3{z}_4/\left(1+z_3\right)\end{array}$

$f_2\left(z_1,z_2,z_3,z_4,\mu \right)=0.3z_2{z}_3-0.065z_2{z}_4+1.469z_2{z}_4/\left(0.3264z_4+1\right)$

$f_3\left(z_1,z_2,z_3,z_4,\mu \right)=0.2807z_3{z}_4/\left(1+0.3264z_4\right)$

$f_4\left(z_1,z_2,z_3,z_4,\mu \right)=-z_1{z}_4-2.5z_2{z}_4+0.6188z_4^2+z_4\mu -0.3264\mu z_4^2$

我们把参数 $\mu$ 作为新的变量：

$\stackrel{\dot{}}{\mu }=0$ (6)

根据中心流形定理可知系统(5)和(6)在平衡点原点O的小邻域内有一个二维的局部中心流形：

$\begin{array}{l}{W}_{loc}^2\left(O\right)=\{\left(z_1,z_2,z_3,z_4,\mu \right)\in R^3\times R^2|z_1=h_1\left(z_4,\mu \right),z_2=h_2\left(z_4,\mu \right),z_3=h_3\left(z_4,\mu \right),\\ h_i\left(0,0\right)=0,D{h}_i\left(0,0\right)=0,i=1,2,3\},其中\left|z_4\right|,\left|\mu \right|充分小\end{array}$

假定中心流形有如下形式：

$h\left(z_4,\mu \right)=\left(\begin{array}{c}{h}_1\left(z_4,\mu \right)\\ h_2\left(z_4,\mu \right)\\ h_3\left(z_4,\mu \right)\end{array}\right)=\left(\begin{array}{c}{a}_1{z}_4^2+a_2{z}_4\mu +a_3{\mu }^2+\cdots \\ b_1{z}_4^2+b_2{z}_4\mu +b_3{\mu }^2+\cdots \\ c_1{z}_4^2+c_2{z}_4\mu +c_3{\mu }^2+\cdots \end{array}\right)$ (7)

首先我们将 $\left(\begin{array}{c}{z}_1\\ z_2\\ z_3\end{array}\right)=\left(\begin{array}{c}{h}_1\left(z_4,\mu \right)\\ h_2\left(z_4,\mu \right)\\ h_3\left(z_4,\mu \right)\end{array}\right)$ 的两边对时间t求导，得到：

$\left(\begin{array}{c}{\stackrel{\dot{}}{z}}_1\\ {\stackrel{\dot{}}{z}}_2\\ {\stackrel{\dot{}}{z}}_3\end{array}\right)=D_{z_4}h\left(z_4,\mu \right){\stackrel{\dot{}}{z}}_4+D_{\mu }h\left(z_4,\mu \right)\stackrel{\dot{}}{\mu }=\left(\begin{array}{c}{D}_{z_4}{h}_1\left(z_4,\mu \right)\\ D_{z_4}{h}_2\left(z_4,\mu \right)\\ D_{z_4}{h}_3\left(z_4,\mu \right)\end{array}\right){\stackrel{\dot{}}{z}}_4$ (8)

再把(5)和(7)代入(8)，得到：

$\begin{array}{l}\left(\begin{array}{c}{D}_{z_4}{h}_1\left(z_4,\mu \right)\\ D_{z_4}{h}_2\left(z_4,\mu \right)\\ D_{z_4}{h}_3\left(z_4,\mu \right)\end{array}\right)\cdot \left(A\cdot z_4+f_4\left(h_1\left(z_4,\mu \right),h_2\left(z_4,\mu \right),h_3\left(z_4,\mu \right),z_4,\mu \right)\right)\\ -B\left(\begin{array}{c}{h}_1\left(z_4,\mu \right)\\ h_2\left(z_4,\mu \right)\\ h_3\left(z_4,\mu \right)\end{array}\right)-\left(\begin{array}{c}{f}_1\left(h_1\left(z_4,\mu \right),h_2\left(z_4,\mu \right),h_3\left(z_4,\mu \right),z_4,\mu \right)\\ f_2\left(h_1\left(z_4,\mu \right),h_2\left(z_4,\mu \right),h_3\left(z_4,\mu \right),z_4,\mu \right)\\ f_3\left(h_1\left(z_4,\mu \right),h_2\left(z_4,\mu \right),h_3\left(z_4,\mu \right),z_4,\mu \right)\end{array}\right)=0\end{array}$ (9)

其中 $A=0$， $B=\left(\begin{array}{ccc}-0.518& 0& 0\\ 0& -0.5& 0\\ 0& 0& -\text{0}\text{.0909}\end{array}\right)$，令(9)中同类项的系数等于0，求得：

$h\left(z_4,\mu \right)=\left(\begin{array}{c}{h}_1\left(z_4,\mu \right)\\ h_2\left(z_4,\mu \right)\\ h_3\left(z_4,\mu \right)\end{array}\right)=\left(\begin{array}{c}0.2898z_4^2+0.743z_4\mu +\cdots \\ 0\\ 0\end{array}\right)$ (10)

把(10)代入系统(5)和(6)，得到系统(5)和(6)限制在平衡点原点的局部中心流形 $W_{loc}^2\left(O\right)$ 上的系统为：

$\{\begin{array}{l}{\stackrel{\dot{}}{z}}_4=-0.2898z_4^3-1.069\mu z_4^2+0.6188z_4^2+\mu z_4+\cdots \\ \stackrel{\dot{}}{\mu }=0\end{array}$ (11)

令 ${\stackrel{\dot{}}{z}}_4=H\left(z_4,\mu \right)$，由于 $H\left(z_4,\mu \right)$ 满足：

$H\left(0,0\right)=0$， $\frac{\partial H}{\partial z_4}\left(0,0\right)=0$， $\frac{\partial H}{\partial \mu }\left(0,0\right)=0$， $\frac{{\partial }^{2}H}{\partial z_4\partial \mu }\left(0,0\right)=1\ne 0$， $\frac{{\partial }^{2}H}{\partial z_4^2}\left(0,0\right)=0.6188\ne 0$，故系统(11)

的平衡点原点在 $\mu =0$ 处发生了Transcritical分支。因为系统(11)与系统(3)拓扑等价，所以系统(3)的平衡点 $P_2$ 在分支参数 ${\rho }_3={\rho }_3^5=1$ 处发生了Transcritical分支，且在1的左侧，平衡点 $P_2$ 是稳定的，平衡点 $P_6$ 是不稳定的；在1的右侧，平衡点 $P_2$ 是不稳定的，平衡点 $P_6$ 是稳定的。

情形2当 ${\rho }_3={\rho }_3^4=0.74588$，同情形1，系统(3)在 ${\rho }_3=0.74588$ 处发生了Transcritical分支，在0.74588的两侧平衡点 $P_6$ 和平衡点 $P_8$ 的稳定性没有发生改变。

情形3当 ${\rho }_3={\rho }_3^3=0.71044$ 时，同情形2，系统(3)在 ${\rho }_3=0.71044$ 处发生了Transcritical分支，在0.71044的两侧平衡点 $P_6$ 和平衡点 $P_9$ 的稳定性没有发生改变。

情形4当 ${\rho }_3={\rho }_3^2=0.62866$ 时，同情形2，系统(3)在 ${\rho }_3=0.62866$ 处发生了Transcritical分支，在0.62866的两侧平衡点 $P_8$ 和平衡点 $P_{11}$ 的稳定性没有发生改变。

2.5. 平衡点的Hopf分支

当 ${\rho }_3={\rho }_3^1,{\rho }_3^6,{\rho }_3^7$ 时，系统(3)存在的非双曲平衡点的特征根有一对纯虚根，且其他特征根实部非零，此时平衡点可能会发生Hopf分支。

情形1： ${\rho }_3={\rho }_3^7=1.31131$ 时

同2.4情形1，计算得到限制在平衡点原点的局部中心流形上的系统为：

$\{\begin{array}{l}\left(\begin{array}{c}{\stackrel{\dot{}}{z}}_3\\ {\stackrel{\dot{}}{z}}_4\end{array}\right)=\left(\begin{array}{cc}0& -\text{0}\text{.4074}\\ \text{0}\text{.4074}& 0\end{array}\right)\left(\begin{array}{c}{z}_3\\ z_4\end{array}\right)+\left(\begin{array}{c}{f}_3\left(h_1,h_2,z_3,z_4,\mu \right)\\ f_4\left(h_1,h_2,z_3,z_4,\mu \right)\end{array}\right)\\ \stackrel{\dot{}}{\mu }=0\end{array}$ (12)

其中：

$\begin{array}{l}{f}_3\left(h_1\left(z_3,z_4,\mu \right),h_2\left(z_3,z_4,\mu \right),z_3,z_4,\mu \right)\\ =-0.0477\mu z_3-0.0282\mu z_4+0.0443z_3^2-0.1304z_3{z}_4+0.0436z_4^2\\ +0.0008{\mu }^2{z}_3+0.00005{\mu }^2{z}_4+0.01946\mu z_3^2-0.0029\mu z_3{z}_4\\ -0.0059\mu z_4^2+0.005z_3^3+0.001z_3^2{z}_4+0.0553z_3{z}_4^2-0.0546z_4^3+\cdots \end{array}$

$\begin{array}{l}{f}_4\left(h_1\left(z_3,z_4,\mu \right),h_2\left(z_3,z_4,\mu \right),z_3,z_4,\mu \right)\\ =-0.0332\mu z_3+0.0196\mu z_4-0.7992z_3^2-0.1998z_3{z}_4+0.9209z_4^2\\ -0.00056{\mu }^2{z}_3-0.000034{\mu }^2{z}_4-0.01964\mu z_3^2-0.04772\mu z_3{z}_4\\ -0.000254\mu z_4^2+0.0203z_3^3-0.0627z_3^2{z}_4-0.087761z_3{z}_4^2-0.0916z_4^3+\cdots \end{array}$

令 $f=f_3\left(h_1,h_2,z_3,z_4,\mu \right),g=f_4\left(h_1,h_2,z_3,z_4,\mu \right)$

$f\left(0,0,0\right)=0,g\left(0,0,0\right),d={\frac{\text{d}\mathrm{Re}\left(\lambda \left(\mu \right)\right)}{\text{d}\mu }|}_{\mu =0}=-0.146<0$

$\begin{array}{c}a=\frac{1}{16}{\left[f_{z_3{z}_3{z}_3}+f_{z_3{z}_{\text{4}}{z}_4}+g_{z_3{z}_{\text{3}}{z}_4}+g_{z_4{z}_4{z}_4}\right]|}_{\left(0,0,0\right)}\\ +\frac{1}{16\omega }{\left[f_{z_3{z}_4}\left(f_{z_3{z}_3}+f_{z_4{z}_4}\right)-g_{z_3{z}_4}\left(g_{z_3{z}_3}+g_{z_4{z}_4}\right)-f_{z_3{z}_3}{g}_{z_3{z}_3}+f_{z_4{z}_4}{g}_{z_4{z}_4}\right]|}_{\left(0,0,0\right)}\\ =-0.0056<0\end{array}$

根据Poincare-Andronuv-Hopf分支定理 [5]，因为 $a<0,d<0$，所以当分支参数 ${\rho }_3>{\rho }_3^7=1.31131$ 时，平衡点 $P_{11}$ 是稳定的。当 ${\rho }_3$ 逐渐减小到 ${\rho }_3^7$ 时，系统(3)在平衡点 $P_{11}$ 处发生了Hopf分支，随着 ${\rho }_3$ 继续减小，平衡点 $P_{11}$ 变得不稳定，系统在 ${\rho }_3^7$ 的左侧邻域生成了一族稳定的周期轨，发生的Hopf分支类型为supercritical。

情形2： ${\rho }_3={\rho }_3^6=1.24771$ 时，同情形1，因为 $a>0,d>0$，所以当 ${\rho }_3$ 减小到 ${\rho }_3^6$ 时，系统(3)在平衡点 $P_9$ 处发生了Hopf分支，随着 ${\rho }_3$ 继续减小，在 ${\rho }_3^6$ 的左侧邻域内生成了一族不稳定的周期轨，发生的Hopf分支类型为subcritical。

情形3： ${\rho }_3={\rho }_3^1=0.09104$ 时，同情形1，因为 $a>0,d<0$，所以当 ${\rho }_3$ 从 ${\rho }_3^1$ 左侧到 ${\rho }_3^1$ 时，系统(3)在平衡点 $P_7$ 处发生了Hopf分支，随着 ${\rho }_3$ 继续增加，在 ${\rho }_3^1$ 的右侧邻域内生成了一族不稳定的周期轨，发生的Hopf分支类型为subcritical。

3. 数值模拟

Figure 1. The bifurcation diagram

图1. 分支图

图1中 $L_1,L_2,L_4,L_6,L_7,L_8,L_9,L_{11}$ 为对应的平衡点曲线；标记B表示系统(3)在这个点处经历了Transcritical分支；标记H表示在系统这个点处经历了Hopf分支；标记D表示系统在这个点处经历了倍周期分支；实线表示平衡点曲线(或者周期轨)是稳定的，虚线表示平衡点曲线(或者周期轨)是不稳定的。

表1和表2分别给出了随 ${\rho }_3$ 变化平衡点的分支和周期解生成的分支。为方便，将主要出现的几个分支类型简写为：HB (Hopf分支)；TB (Transcritical分支)；PDB (倍周期分支)；

Table 1. Bifurcation occured at equilibrium points

表1. 平衡点的分支

Table 2. Bifurcation generated by periodic solution

表2. 周期解生成的分支

当 ${\rho }_3$ 从1.31131的右侧减小到1.31131时，系统(3)在 $P_{11}$ 处发生了HB (6)；随着 ${\rho }_3$ 继续减小，在1.31131的左侧，平衡点 $P_{11}$ 的邻域内生成一族稳定的周期轨，如图2所示。当 ${\rho }_3$ 减小到1.24771时，系统在 $P_9$ 处发生了HB (5)，随着 ${\rho }_3$ 继续减小，在1.24771的左侧，平点 $P_9$ 的邻域内生成一族不稳定的周期轨。当 ${\rho }_3$ 减小到1.05072时，对应的周期解的4个Floquet乘子为 $\left\{1,-1,0.767665,-4.76211\times {10}^{-6}\right\}$，其中一个Floquet乘子在−1处由内而外穿过单位圆导致了PDB (10)的发生，稳定的周期轨失去了稳定性。当 ${\rho }_3$ 小于1.05072时，我们模拟出了倍周期分支的级联现象，如图3(a)~(d)所示。

Figure 2. The projection of periodic orbit and time series chart of $x_3$ for ${\rho }_3=1.30124$

图2. ${\rho }_3=1.30124$ 处周期轨的相图和 $x_3$ 的时间序列图

(a) (b) (c) (d)

Figure 3. The period-doubling cascade: (a) Period 1 orbit at ${\rho }_3=1.1$ ; (b) Period 2 orbit at ${\rho }_3=1.03934$ ; (c) Period 4 orbit at ${\rho }_3=1.02193$ ; (d) Chaotic attractor at ${\rho }_3=1.00196$

图3. 倍周期分支级联现象(a) ${\rho }_3=1.1$ 处的倍1轨线；(b) ${\rho }_3=1.03934$ 处的倍2轨线；(c) ${\rho }_3=1.02193$ 处的倍4轨线；(d) ${\rho }_3=1.00196$ 处的混沌吸引子

当 ${\rho }_3=1.00196$ 时，以 $\left(0.\text{68414},0.\text{2469}0,0.0\text{262}0,0.23384\right)$ 为初始点的轨线是一个混沌吸引子，其对应的Lyapunov指数为 $\left(0.0\text{14222},0,-0.0\text{14}0\text{86},-0.\text{44887}\right)$。若 ${\rho }_3\in \left(0.947834,1.01265\right)$，系统在 $P_{11}$ 邻域内出发的解都是趋于混沌的。

图4给出了当 ${\rho }_3=0.96$ 时，平衡点 $P_{11}$ 邻域内以 $\left(0.\text{6698},0.2501,0.\text{1073},0.2081\right)$ 为初始点出发的混沌轨道在 $x_1,x_2$ 上的投影和 $P_2$ 邻域内以 $S_0\left(0.\text{85},0.4,0.1,0.1\right)$ 为初始点出发随时间增长收敛到平衡点 $P_2=\left(1,0,0,0\right)$ 的稳定收敛轨道。当 ${\rho }_3\in \left(0.947834,1\right)$ 时两个稳定态始终共存。

当 ${\rho }_3\in \left(0.949534,1\right)$ 时，宿主细胞和免疫细胞的初始数目越多，肿瘤细胞和内皮细胞的初始数目越少，轨线就越容易收敛到无肿瘤平衡点 $P_2$。当 ${\rho }_3<0.949534$ 时，如图5所示，当 ${\rho }_3=0.85$， $P_{11}$ 邻域内以 $\left(0.\text{6577},0.\text{24},0.\text{11}0\text{8},0.\text{16}\right)$ 为初始点出发的解随着时间增加会较快地收敛到平衡点 $P_2$，肿瘤细胞的数量随时间增长趋于0，本文研究的肿瘤细胞的增长率的临界值即为0.949534。模拟结果表明，若肿瘤的增长率通过治疗能降低到临界值以下，肿瘤细胞的数量会在较短时间内衰减至0，其他细胞数量趋于定值。

Figure 4. The projection of coexistence of two steady states for ${\rho }_3=0.96$

图4. ${\rho }_3=0.96$ 处的两个稳定态的共存

Figure 5. The projection of phase diagram for ${\rho }_3=0.85$

图5. ${\rho }_3=0.85$ 处的相图

图1~5的数值模拟结果表明肿瘤细胞的增长率在肿瘤增长模型的分支和动态中有着重要作用。

4. 结论

本文研究了肿瘤细胞的自然增长率 ${\rho }_3$ 这个参数对模型动态的影响，理论分析结合数值模拟的结果表明，若能通过化疗或者其他治疗手段抑制肿瘤细胞的增长率到临界值以下并持续一段时间，在自身免疫系统的作用下，肿瘤细胞的数目随时间增加会较快地衰减至0，就可以依靠病人自身的免疫能力获得痊愈。肿瘤细胞的增长率在临界值之上的一定范围内，病人能否痊愈还和细胞初始数目有关，宿主细胞和免疫细胞的初始数目越多(少)，肿瘤细胞和内皮细胞的初始数目越少(多)，那么病人痊愈的几率就越大(小)；在临界值以下，病人几乎不受初始情况影响，随时间增加依靠免疫系统即可痊愈。研究结果对肿瘤的治疗具备一定的理论指导意义。

致谢

在本文的完成过程中，我真诚地想感谢我的导师，她踏实严谨的治学态度，求真求精的工作作风对我产生了深刻影响。本文从选题到最终完成，几经修改，她皆悉心指导，倾注了很大精力。她教会我的“求真”和“严谨”让我受益无穷，在此我向导师致以深切的敬意与感谢。

同时也要感谢其他老师同学的关心和帮助，感谢文献中的学者们关于肿瘤增长模型的研究成果。