1. 前言
血管的生成是固体肿瘤生长过程中起到一个重要作用的转折点。血管未形成时，肿瘤生长所需的营养通过周围组织的扩散过来供应，当肿瘤达到一定的体积后，肿瘤中心的细胞由于缺少营养物质而死亡，当肿瘤新生成的细胞数量与死亡的细胞数量达到一个动态平衡时，肿瘤停止生长。因此，它们能长到的最大直径只有几毫米。在缺氧情况下，肿瘤细胞分泌出VEGF (血管内皮生长因子)，刺激周围的血管产生新的毛细血管并向肿瘤方向生长，形成新的血管。血管的生成增加了肿瘤的营养供给，促使肿瘤继续生长，并侵入到邻近的健康组织。不仅如此，少量的肿瘤细胞有可能进入到血管，通过血管循环被运送到身体的其它地方，迁移形成新的肿瘤[1] 。
目前，人们用了各种方法对肿瘤的生长进行建模，如使用常微分方程组在空间进行建模。Owen MR [2] 提出了一种二维多尺度模型模拟血管肿瘤的生长，模型包括了血流量，血管生成，血管重塑，多个细胞的动态以及单位个细胞的周期等。最近，有研究在Owen MR模型的基础上，结合巨噬细胞和缺氧定位，分析了基因疗法与化学疗法的协同抗癌作用[3] 。
本文主要在Owen MR模型的基础上，结合元胞自动机模型，提出一个肿瘤细胞生长的多尺度建模，并研究了相关的数学模型和模拟计算实现机制，并进行了仿真实验，在一个固定边界浓度的健康组织中嵌入一小块肿瘤细胞，随着时间的增加，观察肿瘤细胞的生长变化。
2. 元胞自动机
2.1. 元胞自动机概述
元胞自动机(Cellular Automata)也称细胞自动机，是一种空间和时间变量都是离散值的动力系统理想化模型，在计算机飞速发展的今天受到各领域研究仿真模拟的研究人员的欢迎[4] [5] 。元胞自动机是目前应用最广、研究最多的一种离散模型。元胞自动机有三个特征[6] ：
平行计算：每一个细胞个体同时同步的改变；
局部性：细胞的状态变化只受到周围细胞的影响；
一致性：所有的细胞受到同样的规则支配。
元胞自动机最常见的是二维模型，细胞自动机中每一个格点的变化受到其邻近的格点状态的影响(这些邻近的格点我们称之为邻居)。邻居的选择很大程度上影响模型计算的速度和结果，通常选择与中心格点相接触的几个为邻居[7] 。二维空间邻居模型，如图1所示。
Figure 1. The black grid is the central cell and the gray grids are the neighbor cells [7]
图1. 黑色格为中心细胞，灰色格为邻居细胞[7]
元胞自动机系统是一个微缩的复杂系统。只有两种状态的50个格点的细胞自动机模型有250 (约1.13*1015)种状态。一个元胞自动机在一个有限时间内不可完全预测，在大量研究生物体内细胞时，仿真度更高，模拟结果更加多元性，具有很强的现实意义。
元胞自动机也有自身的不足，如当一个两个细胞相距较远，发生作用时需要一段较长的时间传递，反应速率降低。
2.2. 肿瘤细胞生长元胞自动机模型
我们的肿瘤细胞生长元胞待机模型是由一个L*L个格点组成的二维矩形网格阵列。每个细胞有3种可能的状态：正常肿瘤细胞，休眠肿瘤细胞，死亡肿瘤细胞。细胞与格点的关系定义为[8] ：
每个细胞只能占据一个格点；
每个格点的细胞只受到邻居细胞的影响；
格点内可以没有细胞；
每个细胞只能向邻居格点移动，当邻居格点没有空间时，细胞不发生分裂。
3. 血管肿瘤生长的多尺度建模
我们的计算模型描述了在一个有血管的组织中肿瘤的生长变化，每一个细胞被看作一个独立的个体，每个细胞拥有自己的细胞周期，细胞周围的组织提供细胞所需营养。在亚细胞尺度，我们通过建立偏微分方程描述细胞周期与细胞周围的氧气浓度密切相关性。细胞消耗氧气，在缺氧条件下分泌出VEGF，细胞外的VEGF刺激已有血管生成新的内皮芽尖，并引导芽尖向缺氧区域生长。当两个芽尖相遇或与已有的血管连结时，新的血管形成。血管直径受到周围氧气浓度和流量相关参数的影响，血管向整个组织提供氧气[9] 。各尺度间的相互作用如示意图2。
我们的模型由一规则网格在规定范围内细分为模拟细胞的格子，每个细胞占据一个格子，这些细胞的增殖由细胞周期决定。在系统初始化后，首先更新扩散领域的状态，然后更新细胞状态(包括细胞的分裂和移动)，由于细胞消耗营养物质和氧气，细胞的状态影响扩散层的营养物质和氧气浓度。一直重复这个过程直到模拟结束。
Figure 2. The connections between the different modeling layers [9]
图2. 模型各尺度之间的相互关系 [9]
3.1. 亚细胞尺度
细胞周期(cell cycle)是指细胞从一次分裂完成开始到下一次分裂结束所经历的全过程，一般认为分为四个阶段。G1期是从有丝分裂到DNA复制前的一段时期，主要合成RNA和核糖体，细胞体积显著增大，为下阶段S期的DNA复制作好物质和能量的准备。S期，即DNA合成期。G2期为DNA合成后期，是有丝分裂的准备期。M期即细胞分裂期。胞周期内有两个阶段最为重要：G1到S和G2到M；这两个阶段正处在复杂活跃的分子水平变化的时期，容易受环境条件的影响[10] 。为了方便进行计算机模拟，把细胞周期简单地表示为F [0,1]，当F = 1时，细胞发生分裂，分裂后，F重置为0。细胞当前状态用如下公式表示[11] ：
(1)
C代表当前氧气浓度，Tmin表示细胞完成一个周期需要的最少时间。C1是我们设置的一个常数，可以看出，当C越大，细胞周期的进程的速度越快。
当氧气浓度小于一定值时，细胞进入休眠，此时：
(2)
当休眠时间超过一定的时间，细胞死亡。当细胞周围氧气高于休眠时氧气的阀时，细胞离开休眠状态[12] 。
3.2. 细胞尺度
这部分主要考虑新肿瘤细胞的生成，移动和死亡。当细胞分裂时产生新的细胞，我们称子细胞。
发生分裂的细胞所在邻居有空间时，子细胞占据所在的空间；如果所在位置没有空间，子细胞移动到氧气浓度较高的邻居空间上；如果所在位置和附近邻居都没有空间，则父细胞死亡，细胞不发生分裂。细胞的分裂过程如图3所示。
3.3. 扩散尺度
扩散层主要考虑氧气的浓度，在这里我们使用一个反应扩散方程来模拟氧气浓度的变化[13] 。周围环境作为氧气的来源，细胞消耗氧气，氧气浓度为设为C，每个细胞消耗氧气的速率为K，固定边界条件，氧气浓度扩散方程为：
Figure 3. The cell division process diagram
图3. 细胞分裂进程图
(3)
其中D表示氧气扩散系数,当格点
有细胞时，
,当格点
没有细胞时，
。
由(3)得：
(4)
令
，有：
(5)
即
(6)
当
没有细胞时，
,当
上有细胞时，
。
3.4. 用有限差分法求解反应扩散方程
我们把整个区域被划分为50*50的均匀单元网格，x方向和y方向的步长均为
，各单元格记为
，
，边界浓度为固定值，当
时，
，为了简便，我们记
，
。根据二阶中心差商公式
(7)
(8)
根据(7)，(8)，式(6)可化为：
(9)
省略
，得
(10)
由此可得点
处氧气的浓度值：
(11)
4. 算法描述
模拟的流程如图4所示。
Figure 4. Tumor cell growth simulation flow chart
图4. 肿瘤细胞生长模拟流程图
4.1. 细胞移动
细胞可沿着自己的邻居自由移动，当细胞所在格点邻居被其它细胞占据时，移动的概率为零。不考虑趋化因素的影响下，细胞由x移动到y的概率为：
(12)
其中，D表示细胞自由移动系数，
是x，y之间的距离。
4.2. 细胞分裂
细胞的分裂取决于细胞周期进程和细胞周围的邻居是否有空间。细胞的分裂的算法设计：
1：Input:parent cell //父细胞
2：Output:the new cell created //生成新的细胞
3：Get location of the parent cell; //确定父细胞所在的格点
4：Get neighbouring locations of the parent cell; //获取父细胞邻居的位置
5：Found_available_neighbour ture; //周围邻居无空间
6：For each n in the neighbourhood do //n为父节点的一个邻居
7：If n is available then
8： Found_available_neighbour ture; // 周围有可用空间的邻居
9： Break for;
10：end
11：end
12：If found_available_neighbour then
13：Create a new cell;
14：Place the new cell in n location; //新生成一个细胞并占据邻居的一个格点
15：Else
16：Reset parent cell cycle; //重置父细胞周期为0
17：End
4.3. 细胞的休眠和死亡
当细胞所在位置的氧气浓度低于某个阀值时，肿瘤细胞进入休眠状态，细胞周期进程停止。当肿瘤长时间休眠时，肿瘤细胞死亡。
正常细胞死亡的算法：
1：Foreach normal cell c in the modle do
2：Determine the current oxygen level; //计算氧气浓度
3：If oxygen level of cell c
4:then
5：Update quiescent_time; //更新休眠时间
6：If c is quiescent for a long time then //休眠到过长时肿瘤细胞死亡
7：Kill c;
8：end
9：else
10：if oxygen level of cell c>threshold_quiescence//氧气深度达到细胞一定的阀值时
11：then
12：leave quiescence; //细胞离开休眠状态
13：Reset quiescent_time; //重置休眠时间为0
14：end
15：end
5. 实验模拟的结果
模拟的区域大小为W = [0,50]*[0,50]，组织区域中作为氧气源的位置处，使用固定边界条件，氧气浓度为30。我们在组织中以格点(15，15)为中心，边长为3个格点的正方形区域内插入9个肿瘤细胞，固定边界浓度为30，根据(11)，由计算机计算出氧气浓度见图5。
使用Matlab R2013b，对图4的肿瘤细胞的生长进行数值模拟，模拟的基本步骤为：
1) 参数赋值，离散化计算区域，根据初始值构建数据结构；
2) 根据有限差分方法求解公式(11)，获得区域内每一点氧气浓度的初始分布稳态解；
3) 在每一个时间步长内：
① 获得肿瘤细胞的位置信息；
② 对每一个位置上的肿瘤细胞，根据该位置上营养物质的浓度计算细胞分裂，移动、休眠和坏死的概率；
③ 对每一个肿瘤细胞，随机选择一个行为方式，根据计算好的概率判断能否执行该行为，然后根据肿瘤细胞的行为机制更新肿瘤细胞的分布情况；
④ 根据肿瘤细胞的分布情况，计算氧气浓度的稳态分布；
⑤ 如果有肿瘤细胞到达区域边界，或者计算时间超出了预订时间，则停止迭代；否则，返回步骤①继续执行。
以下是实验都是经过多次模拟得出的结果。观察肿瘤生长随时间的变化可知，在最初的100次迭代中，肿瘤细胞数目增加，没有休眠和坏死的肿瘤细胞，。在第150次迭代时，肿瘤中心出现了部分休眠细胞，整个肿瘤的半径也越来越大。在迭代第250次时，可以明显看到中心部分细胞坏死。图6(d)中可以明显看到肿瘤的三层结构，最中心空白部分为已经坏死的细胞，深蓝色为休眠细胞，浅蓝色为正在生长的细胞。
在各种参数不变的情况下，改变肿瘤细胞周围的氧气的深度，分别将固定边界的浓度设置为15、30、60，比较肿瘤细胞生长的变化，模拟结果如图7所示。
该结果表明，肿瘤细胞的生长速度与细胞周围氧气浓度成正比，即氧气浓度越高，肿瘤细胞生长的速度越快，氧气浓度越低，肿瘤细胞生长的速度越慢。
6. 总结与展望
本文从多尺度的方法，应用元细胞模型研究了氧气浓度对肿瘤细胞生长的影响。研究表明，氧气浓度越高，肿瘤细胞生长速度越快，氧气浓度越低，肿瘤细胞生长速度越慢，甚至停止生长。实际上，肿瘤细胞生长的微环境非常复杂，氧气只是最重要的因素之一。如在不同的环境下，或加入某种药物，细胞的耐氧能力会有改变，细胞进入休眠状态的氧气浓度阀值也会发生改变等。此外，本文并没有从肿瘤细胞如何生成，即从有到无这一过程作出模拟，对血管和细胞外基质的作用也没有加以考虑，与实际中的肿瘤细胞生长。
Figure 5. Initializes the Oxygen concentration distribution
图5. 初始化氧气浓度分布
(a) (b)
(c) (d)
Figure 6. Cell growth variation: shallow blue cells are tumor cells, dark blue cells are quiescent cells. (a) Initial state, time = 0; (b) Cells begin proliferation, time = 100; (c) Quiescent of tumor cells in the hypoxic region, t = 150; (d) Death of tumor cells in the hypoxic region, t = 250
图6. 细胞生长变化图：浅蓝色表示正在生长的肿瘤细胞，深蓝色表示休眠状况下细胞。(a) 初始状态，t = 0；(b) 细胞开始分裂，t = 100；(c) 中心肿瘤细胞由于缺氧，出现部分休眠细胞，t = 150；(d) 出现死亡肿瘤细胞，t = 250
Figure 7. The influence of oxygen concentration on tumor cell growth
图7. 氧气浓度对肿瘤细胞生长的影响
因此，在未来我们可以通过修改网格大小，或者网格内细胞的数量；在扩散层加入各种生长因子如促血管生长因子和各种营养的分布；在血管层考虑血管的生长，研究肿瘤血管的生长和形成；在细胞层研究普通细胞和肿瘤细胞之间的相互影响；在亚细胞层考虑细胞内部的各种蛋白质对细胞周期的作用。这样，通过完善生长模型，使之更加接近现实肿瘤细胞的生长发展，为下一步应用计算机模拟药物对肿瘤细胞的影响打下基础。