一种基于细胞活动的骨折愈合仿真方法

摘要

一种基于细胞活动的骨折愈合仿真方法，涉及生物医学工程领域。本发明用来预测骨折愈合复杂过程，探寻最优的骨折愈合方案。本发明步骤为：步骤一：建立骨折区域的三维几何模型；步骤二：建立骨折愈合区域有限元模型；步骤三：获取骨折区域单元材料属性和细胞浓度；步骤四：骨折区域有限元分析；步骤五：求解骨折区域单元所受的力学刺激，并根据力调节算法得到新的组织表型；步骤六：分析细胞活动；步骤七：计算新的骨折愈合区域单元材料属性；步骤八：判断程序是否终止。通过本发明可以更加精确的模拟骨折愈合过程，并为医生制定最佳的骨折愈合手术方案提供有益帮助。

说明书

技术领域

本发明涉及生物医学工程领域，特别涉及一种基于细胞活动的骨折愈合仿真方法。

背景技术

骨折是一种常见的创伤，骨折的高发性使得骨折机理及促进愈合的研究尤为迫切。与其他组织损伤修复不同的是，骨折不是靠纤维结缔组织连接，而是骨组织的完全再生。尽管如此，并不是所有的骨折都可以完成愈合，有时会发生延迟愈合甚至是不愈合。骨折延迟愈合或者不愈合会引起患肢疼痛，功能障碍，导致患者失业，由此造成很大的社会经济负担。因此，尽管关于骨折愈合的研究一直备受关注，但仍有5％～10％的骨折因各种原因发生延迟愈合甚至是不愈合。

骨折愈合过程中，包含众多细胞的增殖，分化，凋亡活动。这些细胞包括间充质干细胞，成纤维细胞，软骨细胞和成骨细胞。间充质干细胞具有很高的分化能力，可以分化为成纤维细胞，软骨细胞和成骨细胞。骨折愈合受特定的几何因素、力学因素、生物学因素的影响。其中，力学因素是骨折愈合的总要影响因素，良好的力学刺激会促使间充质细胞分化为不同的细胞，从而形成相应的组织，促使骨折愈合。反之则会导致骨折的延迟愈合甚至是不愈合。目前缺少能精确表达骨折愈合这一复杂过程的计算机仿真方法。现存的骨折愈合仿真方法存在如下缺陷：

1.没有建立专门针对患者的个体化模型；

2.力学因素与细胞分化没有一个确定性关系；

3.仿真方法中只考虑到组织层面，没有考虑细胞的增殖、分化和凋亡等活动；

4.骨折愈合区域的模型和生物学材料设置过于简化。

发明内容

本发明的目的是为了解决骨折愈合仿真中没有考虑细胞的增殖、分化和凋亡活动，力学因和细胞分化没有一个确定性关系和骨折愈合区域的模型和生物力学材料设置过于简化的缺点，而提出的一种基于细胞活动的骨折愈合仿真方法。

本发明的目的通过下述技术方案实现一种基于细胞活动的骨折愈合仿真方法，其特征在于，所述方法的具体步骤为：

步骤一：建立骨折区域的三维几何模型；

步骤二：对步骤一建立的三维几何模型进行网格划分，得到骨折区域的有限元模型；

步骤三：获取骨折区域单元材料属性和细胞浓度；

步骤四：骨折区域有限元分析，求解骨折区域单元所受的八面体剪应变和流速；

步骤五：求解骨折区域单元所受的力学刺激，并根据力调节算法得到新的组织表型；

步骤六：分析细胞活动；

步骤七：计算新的骨折区域单元材料属性；

步骤八：判断程序是否满足终止条件，若不满足，程序则从步骤三开始进入下一个迭代循环；若满足，则程序终止，并记录骨折愈合时间。

本发明的有益效果为：

1.将骨折区域设置为双相多孔弹性模型，相比单相模型，更加符合骨折区域的生物特性，使仿真结果更加精确；

2.通过仿真细胞的扩散、增殖、分化和凋亡来模拟骨折区域的细胞活动，更加符合骨折区域组织分化的实质，是仿真结果更加精确；

3.通过构建骨折愈合仿真方法，可以对医生制定最优的手术方案提供指导，进而提高手术成功率，提高骨折愈合质量，减少骨折骨折不愈合和延迟愈合的情况；

4.通过构建骨折愈合仿真方法，可以对建立的仿真模型进行多次重复实验，减少真实的生物实验，节省时间，提高效率，节省费用，避免人道主义争议。

综上，本发明的仿真方法克服了现有技术的缺点与不足。

附图说明

图1为有限单元模型目标数据获取流程示意图；

图2为基于细胞活动的骨折愈合仿真方法流程图。其中，n＝1代表程序运行的初始步；n>1代表程序运行的后续步。

具体实施方式

具体实施方式1：一种基于细胞活动的骨折愈合仿真方法包括一下步骤：

步骤一：建立骨折区域的三维几何模型；

步骤二：对步骤一建立的三维几何模型进行网格划分，得到骨折区域的有限元模型；

步骤三：获取骨折区域单元材料属性和细胞浓度；

步骤四：骨折区域有限元分析，求解骨折区域单元所受的八面体剪应变和流速；

步骤五：求解骨折区域单元所受的力学刺激，并根据力调节算法得到新的组织表型；

步骤六：分析细胞活动；

步骤七：计算新的骨折区域单元材料属性；

步骤八：判断程序是否满足终止条件，若不满足，程序则从步骤三开始进入下一个迭代循环；若满足，则程序终止，并记录骨折愈合时间。

基于细胞活动的骨折愈合仿真方法流程图如图1所示。

具体实施方式二：本实施方式与具体实施方式一不同的是：所述步骤一中骨折区域的三维几何模型的建立具体为：

采用基于分割的三维医学图像表面重建算法对图像进行三维表面重构，得到三维几何模型；

所述的影像由影像设备CT得到，数据存储格式为DICOM；

由表面模型构建实体模型的过程，就是由表面模型的三角片序列及上下底面构建实体模型的面链表，由表面模型的顶点序列构建实体模型的顶点链表，同时建立起体、环、边、半边链表及各链表中结点的指向关系。以边界模型表达的实体构建过程由一系列欧拉操作实现。基本的欧拉操作包括如下互逆的5对：MVFS，MEV，MEF，MEKR，KFMRH；KVFS，KEV，KEF，KEMR，MFKRH。其中M表示构造，K表示删除，S、E、V、F、R、H分别表示体、边、顶点、面、环、孔。

其他步骤及参数与具体实施方式一相同。

具体实施方式三：本实施方式与具体实施方式一或二不同的是：所述步骤二对三维几何模型采用八面体网格划分单元划分网格，得到骨折区域有限元模型的具体过程为：

将具体实施方式一中建立的三维几何模型导入到网格划分软件中，进行网格划分，得到骨折区域的有限元模型；

经网格划分软件划分网格后会得到很多数据，而本发明中只需要节点坐标和单元编号来表示骨折愈合区域有限元模型，所以需要将得到的数据导入到MATLAB中进行预处理，只提取目标数据，根据目标数据生成后续有限元计算所需要的单元编号和节点坐标两个文件；

所述的节点坐标和单元编号两个文件为txt文本格式的文件；

节点坐标文件包含三列数据，三列数据分别代表每个节点的空间坐标值；

单元标号文件包含四列数据，四列数据分别为每个单元的四个节点的节点序号。

获取目标数据的流程示意图如图2所示；

其他步骤及参数与具体实施方式一或二相同。

具体实施方式四：本实施方式与具体实施方式一至三之一不同的是：所述步骤三获取获取骨折区域单元材料属性和细胞浓度的具体过程为：

骨折区域单元材料属性可分为初始骨折区域单元材料属性赋值和随后的骨折区域单元材料属性赋值；

初始情况下，骨折区域由肉芽组织填充，所以初始骨折区域单元材料属性为肉芽组织的材料属性；

随后的骨折区域单元材料属性由步骤七获得；

骨折区域细胞浓度可由下式表示：

式中，i为细胞类型；nⁱ为细胞i在骨折区域中的浓度；Dⁱ为细胞i的扩散系数；Pⁱ(S)为细胞i的增殖率；Kⁱ(S)为细胞i的凋亡率；S为骨折区域单元的力学刺激；Pⁱ(S)和Kⁱ(S)为力学刺激S的函数；

本发明共涉及四种细胞类型，分别为间充质干细胞，成纤维细胞，软骨细胞和成骨细胞；

扩散系数Dⁱ可由下式求得：

式中，j为组织表型；n_t为组织表型总数量；D_ij为细胞i在组织j中的扩散系数；φ_j为组织j的体积分数；

本发明共涉及四种组织类型，分别为肉芽组织，纤维结缔组织，软骨组织和骨组织；

各组织体积分数有如下关系：

式中，j为组织表型；n_t为组织表型总数量；φ_j为组织j的体积分数。

其他步骤及参数与具体实施方式一至三之一相同。

具体实施方式五：本实施方式与具体实施方式一至四之一不同的是：所述步骤四中利用有限元分析的方法求解骨折区域单元所受的八面体剪应变和流速的具体过程为：

步骤四一：施加外载荷；

在骨折愈合过程中，轴向载荷对骨折愈合起到促进的作用。所以在骨折区域有限元模型的顶端施加轴向外载荷；

步骤四二：设置边界条件；

在骨折区域末端设置完全约束，即骨折区域末端的平移和转动都为0。边界条件的设置还包括间充质细胞的来源的设置。间充质细胞的来源有三个部分：1.骨折区域周围组织；2.骨膜内生层；3.骨髓。具体的间充质细胞来源需要根据骨折区域的CT图像来确定。

步骤四三：骨折区域单元多孔弹性力学模型的建立；

将骨折区域单元看成是多孔弹性材料，则有如下关系式：

a.固体基质、液相和总的应力应变关系如下：

σ^s＝-φ_spI+σ^E>

σ^f＝-φ_fpI>

σ^t＝σ^s+σ^f＝-pI+σ^E>

式中，σ^s、σ^f、σ^t分别为固相、液相和总的应力张量；p为液体压力；φ_s、φ_f分别为固相和液相体积分数；σ^E为有效应力张量；I为单位张量；

线弹性材料的有效应力张量可以表示为：

σ^E＝Cε>

式中，σ^E为有效应力张量；C为刚度张量；ε为总的弹性应变；

刚度张量由下式表示：

式中，E为弹性模量；υ为泊松比；

b.考虑到两相均不可压缩性和各向同性，多孔弹性模型的连续性方程为：

式中，φ_f为液相的体积分数；v^s、v^f分别为固相和液相的速度向量；

c.固相和液相的动量方程如下：

式中，π^s、π^f分别为固相和液相的体力向量；φ^f为液相体积分数；k为渗透率；v^f、v^s分别为固相和液相的速度向量；σ^s、σ^f、σ^E分别为固相、液相和有效的应力张量；p为液体压力；

通过有限元法求解上述方程，可求得弹性应变ε和液相的速度向量v^f。

步骤四四：求解八面体剪应变；

八面体剪应变与弹性应变ε之间的关系如下：

式中，γ₈为八面体剪应变；分别为x方向、y方向、z方向的正应变；分别为xy面、yz面、xz面上的剪应变；

通过上述求解，最终得到骨折区域单元所受的八面体剪应变γ₈和液相的速度向量v^f。

其他步骤及参数与具体实施方式一至四之一相同。

具体实施方式六：本实施方式与具体实施方式一至五之一不同的是：所述步骤五中求解骨折区域单元所受的力学刺激，并根据力调节算法得到新的组织表型的具体过程为：

步骤五一：求解骨折区域单元所受的力学刺激；

骨痂单元的力学刺激与八面体剪应变和液相的速度具有如下关系：

式中，S为骨痂单元的力学刺激；γ₈为八面体剪应变；v^f为液相的速度；a，b分别为经验常数；

液相速度v^f可由下式求得：

式中，v^f为液相的速度向量；v_x、v_y和v_z分别为x方向、y方向和z方向上的速度；

步骤五二：根据力调节算法获得对应的组织表型；

骨折区域中的组织分化是由骨折区域所受的力学刺激决定的，不同的力学刺激会使间充质干细胞向不同的细胞类型分化，进而形成不同的组织表型。本发明中由间充质细胞分化形成的组织表型共有三种：纤维结缔组织、软骨组织和骨组织。其中，骨组织又分为成熟的骨组织和不成熟的骨组织。不同的组织表型所对应的力学刺激边界值如下所示：

其他步骤及参数与具体实施方式一至五之一相同。

具体实施方式七：本实施方式与具体实施方式一至六之一不同的是：所述步骤六中分析细胞活动的具体过程为：

所述的细胞活动是指细胞的扩散、增殖、分化和凋亡。细胞的扩散可由扩散方程表示，表达式如下所示：

式中：i为细胞类型；为细胞i的扩散函数；Dⁱ为细胞i的扩散系数；nⁱ为细胞i的浓度；

细胞的增殖和凋亡都是关于力学刺激的函数，他们之间的关系如下式所示：

式中，Pⁱ(S)为细胞i的增殖率；Kⁱ(S)为细胞i的凋亡率；nⁱ为细胞i的浓度；S₀为八面体剪应变；aⁱ、bⁱ、cⁱ分别为与细胞i相关的系数；

由于在本发明中，细胞类型i共有四种：间充质细胞，成纤维细胞，软骨细胞和成骨细胞。所以上述表达式还可以有如下表达：

式中，a_mesenchymal、b_mesenchymal、c_mesenchymal分别为与间充质干细胞相关的系数；a_fibroblast、b_fibroblast、c_fibroblast分别为与成纤维细胞先关的系数；a_chondrocyte、b_chondrocyte、c_chondrocyte分别为与软骨细胞相关的系数；a_osteoblast、b_osteoblast、c_osteoblast分别为与成骨细胞相关的系数。

其他步骤及参数与具体实施方式一至六之一相同。

具体实施方式八：本实施方式与具体实施方式一至七之一不同的是：所述步骤七中计算新的骨折愈合区域单元材料属性的具体过程为：

在骨折愈合初期，骨折区域有肉芽组织填充。随着骨折愈合的进行，间充质干细胞扩散到骨折区域，并且发生增殖，分化，凋亡过程。随着间充质干细胞的侵入，骨折区域单元的材料属性也随之发生变化。新的骨折区域单元材料属性，即骨折区域单元的弹性模量和泊松比可由下式求得：

式中，E为更新后的骨折区域单元弹性模量；为细胞最大浓度；n_c为细胞总浓度；E_granulation为肉芽组织弹性模量；j为组织类型；E_j为组织j的弹性模量；φ_j为组织j的体积分数；

式中，υ为更新后的骨折区域单元弹性模量；为细胞最大浓度；n_c为细胞总浓度；υ_granulation为肉芽组织弹性模量；j为组织类型；υ_j为组织j的弹性模量；φ_j为组织j的体积分数；

其他步骤及参数与具体实施方式一至七之一相同。

具体实施方式九：本实施方式与具体实施方式一至八之一不同的是：所述步骤八中判断程序是否终止的具体过程如下：

当骨折达到愈合后，骨折区域中只有骨组织。所以程序终止的判断条件为骨折区域所有单元的材料属性与骨的材料属性是否相等。

若不相等，则程序从步骤三开始进入下一个迭代循环；

若相等，则程序终止，并记录骨折愈合时间。

其他步骤及参数与具体实施方式一至八相同。

实施例：

为了说明本仿真的实施方法，下面具体举一个例子来阐述本发明的操作过程。

模拟人体胫骨骨折愈合过程

(1)获取人体胫骨骨折区域的三维几何模型

利用CT成像设备扫描患者胫骨骨折区域，得到胫骨骨折区域的CT三维体DICOM数据，采用基于分割的三维医学表面重建算法对CT图像进行三维表面重构，得到骨折区域三维几何模型；

(2)获取人体胫骨骨折区域有限元模型

将上一步获得的三维几何模型导入到网格划分软件中进行网格划分，由于经过网格划分后会得到很多数据，而本发明的有限元模型可由单元编号和节点坐标来表示，所以将由网格划分得到的数据导入matlab中进行预处理，得到最终需要的单元编号和节点坐标。这两个目标数据为txt文本格式文件，单元编号文件包含四列数据，四列数据分别为每个单元的六个节点的节点序号；节点坐标文件包含三列数据，三列数据分别代表每个节点的空间坐标值。

(3)获取材料属性和细胞浓度

a.材料属性的获取：

各组织的材料属性如表1所示。

表1各组织材料属性

骨折初期，即在初始情况下，骨痂由肉芽组织组成，所以初始情况下，骨痂内单元材料属性为肉芽组织的材料属性。在随后的迭代步中，骨痂内发生了组织分化，所以骨痂内单元材料属性由本发明给出的公式来求得骨痂内单元材料属性。

b.细胞浓度的获取

细胞浓度有下述公式求得：

式中，i为细胞类型；nⁱ为细胞i在骨折区域中的浓度；Dⁱ为细胞i的扩散系数；Pⁱ(S)为细胞i的增殖率；Kⁱ(S)为细胞i的凋亡率；S为骨折区域单元的力学刺激；Pⁱ(S)和Kⁱ(S)为力学刺激S的函数；

(4)由力调节算法获取组织分化后的组织表型

将上述建立的有限元模型，以及初始材料属性和细胞浓度的赋值，对胫骨骨折区域进行有限元分析，求得八面体剪应变和流体速率。由八面体剪应变和流体速率得到骨痂单元所受的力学刺激。通过力调节算法，根据求得的骨痂单元所受的力学刺激，得到相应的组织分化后的组织表型。

(5)判断程序是否终止

计算组织分化后骨痂单元材料属性并于骨组织材料属性进行比较。

若骨痂单元材料属性不等于骨组织材料属性，则进入下一个迭代循环；

若骨痂单元材料属性等于骨组织材料属性，则程序结束并且输出愈合时间。