组织工程软骨培养的细胞自动机算法模型与优化

求冰霞 李 宏* 王钰安 李文康 陶 蒙 陈丰农 陈 洁

1(杭州电子科技大学生命信息与仪器工程学院，杭州 310018)2(上海交通大学医学院第九人民医院，上海 200031)

组织工程软骨培养的细胞自动机算法模型与优化

求冰霞1李 宏1*王钰安1李文康1陶 蒙1陈丰农1陈 洁2

1(杭州电子科技大学生命信息与仪器工程学院，杭州 310018)2(上海交通大学医学院第九人民医院，上海 200031)

建立一个基于细胞自动机算法的软骨细胞和干细胞体外动态培养数学模型，并进行模拟计算与优化，为组织工程软骨的体外培养控制软件的编写提供理论依据。建立细胞自动机算法模型，采用软骨细胞和软骨干细胞作为研究对象进行Matlab仿真实验。实验分为3组，每组各自设置空白对照组、快速扩增组、力学刺激组、扩增与力学刺激复合培养组。空白对照组的结果显示，软骨细胞和干细胞诱导软骨细胞的细胞占空间区域分别为51%和45%；综合考虑细胞生长情况，软骨细胞培养最佳压强施加大小为100 kPa；干细胞诱导培养最佳压强施加大小为15 kPa。运用细胞自动机算法建立的数学模型来动态培养软骨细胞，可以为体外生物反应器控制软件的编写提供理论依据，对于组织工程培养软骨组织及其他器官培养是一种有效的可供参考的方法。

组织工程；仿真；细胞自动机算法；软骨细胞

引言

组织工程技术为体外培养构建工程化软骨带来了希望，而静态培养的组织工程软骨存在力学性能差、形态不完整、质量和安全性难以控制等问题[1-2]；研究发现，适当的力学刺激对于软骨细胞基质分泌、软骨组织塑形具有良好的调控作用[3]，组织工程软骨生物反应器的开发能为软骨组织的体外培养提供生物的和力学的环境，目前，采用软骨生物反应器培养的软骨组织在力学性能和形态上有明显的改善[4]，但是，软骨生物反应器控制软件由于缺乏良好的控制算法而使得软骨生物反应器的培养效果不理想，所以，寻求一种智能优化算法对于完善软骨生物反应器具有重要的意义。

目前已有研究采用数学方法来解决组织工程方面的问题，许多数学算法被运用到这一方面。遗传算法是一种通过模拟自然进化过程搜索最优解的方法，但遗传算法并不能完美模拟细胞的生长方式，只是给出了一种最优细胞培养环境。神经网络算法是由大量的简单基本元件相互联接而成的自适应非线性动态系统算法，具有很强的非线性拟合能力及强大的自学习能力，但存在的问题是难以分析神经网络的各项性能指标，学习时间较长，对实例模型和网络模型之间的矛盾难以解决等。当数据不充分的时候，神经网络算法运行起来就比较困难[5-6]。蚁群算法也是一种进化优化的算法,但由于该算法是一种概率型算法，算法中参数的设定通常依靠实验方法得到，而与人的经验值关系不大，因此运用蚁群算法很难使系统的性能达到最优[7-8]。

利用细胞自动机算法，可以模拟解决各种各样的复杂事物的演化过程。Liu等运用细胞自动机算法建立的模型可以用来处理医学图像[9]。刁元波使用细胞自动机构建了生物序列的离散模型，对细胞信号网络中具有关键作用的跨膜蛋白进行了拓扑结构预测[10]。在生物学中，细胞自动机算法也已有了一定的发展，像肿瘤生长[11]、生物病毒仿真[12]以及菌落的生长[13]等。对于细胞的增长模式，正是在外界环境及自身内部结构的影响下，细胞之间及周围环境相互作用所决定的。因此，采用细胞自动机算法来模拟组织工程软骨的培养就比较合适。

1 材料与方法

1.1 模型建立

假设营养液在注入空间区域后随机流动，充分体现营养液的流动性质。即营养液微粒随机分布在培养区域中，保证每个细胞生长概率一样。

细胞生长的空间区域用一个二维四方格区域L=(Lx×Ly)表示，其中(Lx,Ly∈Z)。一个细胞用区域中的一个元胞空间表示，其位置为(Lx,Ly)，营养微粒和代谢微粒同样用一个元胞空间表示。初始阶段，细胞随机的分布在网格区域，处于自身代谢状态，在吸收营养微粒之后，便开始进行状态的转变——增殖、代谢、死亡。

模拟实验中，在没有吸收到营养时细胞处于自身代谢状态，称之为普通状态细胞。此状态下细胞将会进行周期性的自身代谢，自身代谢消耗存储的营养物质，如果始终得不到营养微粒的补充，那么细胞死亡；如果吸收了营养微粒，细胞状态就会激活，进入增殖状态，这时如果在该细胞周围格子气区域有空位，就增殖为两个细胞，分裂后细胞进入普通状态。此后，若周围格子气区域存在营养微粒，则细胞可以吸收营养进入下次繁殖周期，否则进入潜生阶段(潜生阶段为细胞饥饿状态，在此期间如果有营养微粒到来，可以吸收转变为普通的自身代谢状态，否则死亡)。

定义细胞的4种状态空间分别为F=0,2,3,4，其中0—空，2—自代谢状态，3—增殖状态，4—潜生状态。Ft+1(x,y)表示转换后的细胞状态，Ft(x,y)表示转化前的细胞状态。tf为细胞在吸收了营养微粒后的等待代谢时间。n(x,y)为格子气区域的营养微粒个数。模型的细胞状态转变方程为

Ft+1(x,y)=

在此转变方程的基础上，加入外部受力，外部压强的大小影响到代谢时间tf的长短，进而影响细胞代谢和扩增速率。所以在算法中加入代谢时间tf的影响。细胞的增殖效果，会由于tf的长短而发生改变。

躁狂,是神志失常的一种证侯。多因肝经热盛;或痰火上扰;阳明热盛,热扰心神;或秽浊上干;血蓄下焦,瘀热上冲等所致。在临床治疗的实践中,多采取药物联合的方式予以治疗。根据国内相关研究结果表明,采用利培酮联合碳酸锂治疗躁狂具有很好的疗效,本文结合本院在2015年5月到2017年5月期间收治的32例躁狂患者作为研究对象,对此展开研究,现报道如下:

细胞在转化规则的控制下完成其自身状态的改变，所有细胞的协同行为便形成一定的生物模式，因此利用自动机算法建成了一个软骨细胞培养模型。

1.2 实验方法

采用Matlab仿真软件，通过矩阵和图像的相互转换对细胞生长模型进行实验仿真。在Matlab中，利用矩阵，将初始时刻的自身代谢状态的细胞随机分布在培养区域。然后，在此矩阵基础上建立起细胞自动机模型，并将营养微粒随机分布于矩阵中。之后，细胞就会按照细胞生长模型进行繁殖、生长，最后以图像的形式显示出来。

首先设定软骨组织结构的力学性能和初始培养环境，根据对象的不同实验分为3组，分别为空白对照组、软骨细胞培养组和干细胞诱导培养组。第1组为空白对照组，空白对照组不施加压强，分别对软骨细胞和干细胞进行培养。第2组为软骨细胞培养组，分成A、B、C等3个小组，施加初始压强为55 kPa，之后每次依次增加22.5 kPa，培养方法为：令A组规模最大化增长，不考虑细胞的增长模式是否健康、所繁殖软骨细胞性能是否最佳；对B组施加压强，不考虑软骨细胞的增长速度和规模；C组综合考虑，同时增加压强刺激又考虑增长规模快且多；第3组是干细胞诱导成软骨细胞培养，分成D、E、F等3个小组，初始压强为5 kPa，之后每次依次增加5 kPa，实验条件和培养方法与第2组相同。

对于每组仿真，其培养环境是相同的，仿真的循环周期也是相同的，确保每种仿真细胞生长的时间相同，可充分体现压强的影响。每组仿真反馈了不同的信息，通过图像显示出来。

2 结果

在模拟仿真中，培养区域为100×100的二维区域空间，即有100×100个元胞空间供细胞、营养微粒、代谢微粒占据。图形中红色代表细胞，蓝色代表营养液微粒，深蓝色代表空的元胞空间，黑色代表细胞繁殖过程中的代谢微粒。

2.1 空白对照组

图1 无压强作用下的干细胞培养生长过程。 (a)100 h; (b)300 h; (c)500 hFig.1 Stem cells grow in the case of no mechanical force. (a)100 h; (b)300 h; (c)500 h

图2 无压强作用下的软骨细胞培养生长过程。 (a)100 h; (b)300 h; (c)500 hFig.2 Chondrocytes grow in the case of no mechanical force. (a)100 h; (b)300 h; (c)500 h

2.2 软骨培养组

实验中，对软骨细胞施加55 kPa的压强，可得到明显促进细胞生长的效果。在55 kPa为初始受力条件下，每次依次增加22.5 kPa，随着压强的增加，细胞增殖效果与压强大小呈线性正相关关系，但当压强超过122.5 kPa时，细胞增殖过程开始恶化。仿真结果显示，综合考虑细胞生长时(同时增加压强刺激又考虑增长规模快且大)所施加的压强为100 kPa。

图3显示压强122.5 kPa情况下，软骨细胞生长规模达到最大，体积占整个空间区域的90%；图4显示压强大小为77.5 kPa时，软骨细胞性能达到最佳，体积占整个空间区域的69%；图5显示压强大小为100 kPa情况下，软骨细胞综合生长达到最优，体积占整个空间区域的81%。此外，还可以发现逐渐加大压强对软骨细胞生长的影响，图6、7分别为施加压强77.5和150 kPa情况下的软骨细胞生长情况，细胞最终所占空间区域为62%和66%。

图3 122.5 kPa作用下软骨细胞生长最大化过程。(a)100 h; (b)300 h; (c)500 hFig.3 Chondrocytes expansion process under 122.5 kPa. (a)100 h; (b)300 h; (c)500 h

图4 77.5 kPa作用下软骨细胞性能最优生长过程。(a)100 h; (b)300 h; (c)500 hFig.4 Chondrocytes optimal performance growth process under 177.5 kPa. (a)100 h; (b)300 h; (c)500 h

图5 100 kPa作用下软骨细胞综合生长最优过程。(a)100 h; (b)300 h; (c)500 hFig.5 Chondrocytes comprehensive growth process under 100 kPa. (a)100 h; (b)300 h; (c)500 h

图6 压强77.5 kPa作用下软骨细胞培养生长过程。(a)100 h; (b)300 h; (c)500 hFig.6 Chondrocytes growth process under 77.5 kPa. (a)100 h; (b)300 h; (c)500 h

图7 压强150 kPa作用下软骨细胞培养生长过程。(a)100 h; (b)300 h; (c)500 hFig.7 Chondrocytes growth process under 150 kPa. (a)100 h; (b)300 h; (c)500 h

2.3 干细胞诱导培养组

实验中，对干细胞施加5 kPa的压强，可得到明显促进细胞生长的效果。在5 kPa为初始受力条件下，每次依次增加5 kPa力，可得到如下结果：在压强为20 kPa时细胞生长情况最好，继续增加压强到25 kPa时细胞生长情况转坏。在诱导成软骨细胞后，施加压强122.5 kPa(第2组软骨细胞培养中的最佳压强大小)进行培养。仿真结果显示，综合考虑细胞生长时，干细胞培养阶段压强应为15 kPa，软骨细胞培养阶段压强应为100 kPa。

图8显示压强20 kPa时，软骨干细胞生长规模达到最大，体积占空间区域的86%；图9显示压强10 kPa的情况下，软骨干细胞的细胞性能达到最优，体积占空间区域的67%；图10显示压强大小为15 kPa情况下，软骨干细胞综合生长达到最优，体积占空间区域的75%。此外，还可以发现逐渐加大压强对干细胞生长的影响，图11、12为压强5和25 kPa时的生长情况，细胞最终所占空间区域大小为54%和75%。

图8 20 kPa作用下干细胞生长最大化过程。(a)100 h; (b)300 h; (c)500 hFig.8 Stem cells rapid expansion process under 20 kPa. (a)100 h; (b)300 h; (c)500 h

图9 10 kPa作用下干细胞性能最优过程。(a)100 h; (b)300 h; (c)500 hFig.9 Stem cells optimal performance growth process under 10 kPa. (a)100 h; (b)300 h; (c)500 h

图10 15kPa作用下干细胞综合生长最优过程。(a)100 h; (b)300 h; (c)500 hFig.10 Stem cells comprehensive growth process under 15 kPa. (a)100 h; (b)300 h; (c)500 h

图11 5 kPa作用下干细胞生长情况过程。 (a)100 h; (b)300 h; (c)500 hFig.11 Stem cells growth process under 5 kPa. (a)100 h; (b)300 h; (c)500 h

图12 25 kPa作用下软骨细胞生长情况过程。 (a)100 h; (b)300 h; (c)500 hFig.12 Chondrocytes growth process in 25 kPa. (a)100 h; (b)300 h; (c)500 h

3 讨论和结论

在这一系列仿真中，通过3组的模拟和比较，可以得出组织软骨的最佳培养策略。对于第2组和第3组，可以看到随着压强的增大，组织的生长情况越来越好，但当压强增大到一定数值时，就达到促进组织生长的极限效果，继续增大压强，就会导致组织结构发生破坏，这是因为过大的受力和挤压导致细胞间相互作用被破坏，细胞表面结构损坏，进而导致细胞之间的黏性发生损伤，甚至可能导致细胞的破裂。

第2组和第3组比较，可以得出力学因素诱导干细胞为软骨细胞的过程和影响，为虚拟仿真诱导干细胞培养软骨组织提供控制理论依据，这为组织及其他器官培养提供了虚拟培养方法，利用算法虚拟培养可以节省培养成本，提高效率。

运用细胞自动机算法建立的数学模型来虚拟动态培养软骨细胞，可以很好地控制外界培养条件，实现对关键参数的控制和优化，获得了最佳的培养效果，为完善组织工程软骨生物反应器智能软件的编写提供理论依据，这是一种非常有效的方法。

[1] Matsiko A, Levingstone TJ, O'Brien FJ. Advanced strategies for articular cartilage defect repair[J]. Materials, 2013, 6(2): 637-668.

[2] 杨波,张安,曹峻岭. 组织工程修复关节软骨损伤研究进展[J]. 国外医学(医学地理分册),2007,1:1-6.

[3] 李正章,程应樟,吴险峰,等. 周期性压力培养对兔骨髓间充质干细胞增殖的影响[J]. 中国组织工程研究与临床康复,2009,19:3628-3632.

[4] Couet F, Mantovani D. How to optimize maturation in a bioreactor for vascular tissue engineering: focus on a decision algorithm for experimental planning [J]. Annals of Biomedical Engineering,2010, 38(9): 2877-2884.

[5] 赵其刚，李群湛，彭虎.小波神经网络在NGN流量预测中的应用[J].计算机工程，2006,32(16):20-22.

[6] 柳松青. MATLAB神经网络BP网络研究与应用[J]. 计算机工程与设计,2003,11:81-83,88.

[7] 闵克学，葛宏伟，张毅，等.基于蚁群和粒子群优化的混合算法求解TSP问题[J].吉林大学学报, 2006,24(4):402-405.

[8] 詹士昌，徐婕，吴俊蚁. 群算法中有关算法参数的最优选择[J].科技通报，2003,19(5)：381-386.

[9] 刁元波. 复杂网络理论和细胞自动机在生物信息学中的应用研究[D]. 成都：四川大学,2007.

[10] 黄文高，潘志庚.人工生命在计算机图形学中的应用[J].计算机辅助设计与图形学学报，2005,17(7):1383-1388.

[11] 刘天庆，卜葛丹，程坊，等. 神经干细胞球生长动力学模型[J].中国生物医学工程学报，2006,25(6):708-716.

[12] 冯海林. 生物病毒复杂系统定性仿真研究[D].合肥：中国科学技术大学,2007.

[13] 赵峰, 陶祖莱. 生物模式的元胞自动机模型 (Ⅱ): 菌落的生长模式[J]. 生物医学工程学杂志, 2007, 24(4): 820-823.

The Cellular Automata Model and Optimization Algorithm of Tissue Engineering Cartilage

Qiu Bingxia1Li Hong1*Wang Yu′an1Li Wenkang1Tao Meng1Chen Fengnong1Chen Jie2

1(CollegeofLifeInformationScience&InstrumentEngineering,HangzhouDianziUniversity,Hangzhou310018,China)2(TheNinthPeople′sHospital,CollegeofMedicine,ShanghaiJiaotongUniversity，Shanghai200031,China)

tissue engineering; simulation; cellular automata algorithm; cartilage cells

10.3969/j.issn.0258-8021. 2016. 02.018

2015-06-10， 录用日期:2015-10-20

上海自然科学基金项目 (09ZR1432700)

R318

A

0258-8021(2016) 02-0252-05

*通信作者(Corresponding author)， E-mail: lihong@hdu.edu.cn