一种基于力学刺激和生物联合刺激的骨折愈合仿真系统

摘要

一种基于力学和生物学联合调控的骨折愈合仿真系统，涉及生物医学工程领域。所述系统包括骨折区域分析模型建立模块，骨折区域力学调控建模模块，骨折区域生物调控建模模块，骨折区域力‑生物联合调控建模模块和程序终止判断模块。骨折区域分析模型建立模块用于建立骨折区域几何模型和有限元模型。骨折区域力学调控建模模块用于求解骨折区域单元所受的力学刺激。骨折区域生物调控建模模块用于相关细胞浓度、生长因子浓度和胞外基质浓度。骨折区域力‑生物联合调控建模模块用于实现力学刺激和生物刺激的联合调控。程序终止判断模块用于判断程序是否终止。本发明通过对骨折愈合过程进行仿真，为探寻最佳的骨折愈合方案提供有益帮助。

说明书

技术领域

本发明涉及生物医学工程领域，特别涉及一种基于力学刺激和生物刺激联合调控的骨折愈合仿真系统。

背景技术

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

骨折愈合是一个非常复杂的生物过程，涉及到力学刺激和生物刺激。骨折愈合过程包括众多细胞的增殖，分化，凋亡活动。这些细胞包括间充质干细胞，成纤维细胞，软骨细胞和成骨细胞。在特定几何因素，力学因素和生物学因素的刺激下，具有很高分化能力的间充质干细胞会向软骨细胞和骨细胞分化，从而形成相应的骨组织，促使骨折愈合。反之，若几何因素，力学因素或生物学因素不满足骨折愈合发生的条件，则会导致骨折延迟愈合甚至是不愈合。目前缺少能精确表达由力学刺激和生物刺激联合调控的骨折愈合这一复杂过程的计算机仿真系统。现存的骨折愈合仿真系统存在如下缺陷：

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

2.没有建立力学因素和生物学因素的联合调控；

3.力学因素与细胞的增殖，分化以及凋亡没有一个确定性关系；

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

发明内容

本发明的目的是为了解决骨折仿真愈合中没有将力学刺激和生物刺激联合起来调控骨折愈合过程，力学刺激与细胞的增殖，分化以及凋亡没有一个确定性关系和骨折愈合区域的模型和材料属性设置过于简化的缺点，而提出的一种基于力学刺激和生物刺激联合调控的骨折愈合仿真系统。

本发明的目的通过下述技术方案实现一种基于力学刺激和生物刺激联合调控的骨折愈合仿真系统，其特征在于，所述系统包括：

骨折区域分析模型建立模块、骨折区域力学调控建模模块、骨折区域生物调控建模模块、骨折区域力-生物联合调控建模模块和程序终止判断模块。

骨折区域分析模型建立模块用于建立骨折区域几何模型和有限元模型；

骨折区域力学调控建模模块用于对骨折区域的力学特性进行建模，求解骨折区域单元所受的力学刺激；

骨折区域生物调控建模模块用于对骨折区域的生物特性进行建模，求解与骨折愈合相关的细胞浓度、生长因子浓度和胞外基质浓度；

所述的与骨折愈合相关的细胞浓度、生长因子浓度和胞外基质浓度分别为：间充质干细胞浓度，成纤维细胞浓度，软骨细胞浓度，骨细胞浓度，血管生成细胞浓度，软骨生长因子浓度，骨生长因子浓度，血管生成生长因子浓度，成纤维细胞胞外基质浓度，软骨细胞胞外基质浓度，骨细胞胞外基质浓度，血管生成细胞胞外基质浓度。

骨折区域力-生物联合调控建模模块用于建立力学刺激与由力学刺激引起的相关生理活动的函数关系，从而实现力学刺激和生物刺激的联合调控；

程序终止判断模块用于判断程序是否终止，如不满足终止条件，则程序进入下一迭代步；若满足终止条件，则程序结束，记录骨折愈合时间，绘制相关细胞浓度、相关生长因子浓度及相关细胞胞外基质浓度随时间变化的云图。

本发明的有益效果为：

1.本发明提出的一种基于力学刺激和生物刺激联合调控的骨折愈合仿真系统是基于 windows开发语言平台来开发软件，通过自主编程实现骨折愈合过程的动态模拟，基于对话框的形式，易于操作，培训周期短；

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

3.采用梯形函数实现力学刺激和生物刺激的联合调控；

4.采用耦合偏微分方程组来对骨折愈合过程进行仿真，可以得到相关细胞浓度、生长因子浓度、细胞胞外基质浓度随时间的变化关系；

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

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

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

附图说明

图1为基于力刺激和生物刺激联合调控的骨折愈合仿真系统流程图；

图2为骨折区域几何模型建立流程图；

图3为骨折区域有限元模型建立流程图；

图4为骨折愈合过程中相关细胞的迁移、增殖、分化及分泌过程示意图；

图5力学刺激调控骨折愈合过程中相关生理过程梯形函数图。

具体实施方式

具体实施方式一：如图1所示，本实施方式所述的一种基于力学刺激和生物刺激联合调控的骨折愈合仿真系统包括：

骨折区域分析模型建立模块1、骨折区域力学调控建模模块2、骨折区域生物调控建模模块3、骨折区域力-生物联合调控建模模块4和程序终止判断模块5；

骨折区域分析模型1建立模块用于建立骨折区域几何模型和有限元模型；

骨折区域力学调控建模模块2用于对骨折区域的力学特性进行建模，求解骨折区域单元所受的力学刺激；

骨折区域生物调控建模模块3用于对骨折区域的生物特性进行建模，求解与骨折愈合相关的细胞浓度、生长因子浓度和胞外基质浓度；

所述的与骨折愈合相关的细胞浓度、生长因子浓度和胞外基质浓度分别为：间充质干细胞浓度，成纤维细胞浓度，软骨细胞浓度，骨细胞浓度，血管生成细胞浓度，软骨生长因子浓度，骨生长因子浓度，血管生成生长因子浓度，成纤维细胞胞外基质浓度，软骨细胞胞外基质浓度，骨细胞胞外基质浓度，血管生成细胞胞外基质浓度。

骨折区域力-生物联合调控建模模块4用于建立力学刺激与由力学刺激引起的相关生理活动的函数关系，从而实现力学刺激和生物刺激的联合调控；

程序终止判断模块5用于判断程序是否终止，如不满足终止条件，则程序进入下一迭代步；若满足终止条件，则程序结束，记录骨折愈合时间，绘制相关细胞浓度、相关生长因子浓度及相关细胞胞外基质浓度随时间变化的云图。

具体实施方式二：如图1-5所示，本实施方式中，所述的骨折区域分析模型建立模块1 实现其功能的具体过程为：

1)骨折区域三维表面几何模型的建立；

采用基于分割的三维医学影像表面重建算法对图像进行表面重构，通过阈值筛选、交互式分割和三维重建过程得到三维表面几何模型；

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

2)骨折区域有限元模型的建立；

将骨折区域三维表面几何模型进行网格划分，使连续的几何模型离散化，得到骨折区域有限元模型；

所述的网格划分包括面网格划分和体网格划分两个步骤：面网格划分过程用于将三维表面模型进行优化，包括：表面模型优化、平滑处理、修补漏洞；表面模型的优化通过减小表面模型的三角面片来实现，该过程只需将相邻的两个顶点合并到一个新的顶点上，并延续原有的拓扑关系；平滑处理的过程中，对三维的面网格模型进行去噪；修补漏洞的过程中，通过将模型当中的空洞提取成空间多边形，然后对空洞多边形进行三角化的方法实现；体网格划分的过程是将面网格模型进行拉伸、旋转步骤实现的；

通过网格划分得到的骨折区域有限元模型包括单元编号和节点坐标两部分；

单元编号包含四列数据，四列数据分别为每个单元的四个节点的节点序号；

单元坐标包含三列数据，三列数据分别为每个节点的空间坐标值；

其他组成及连接与具体实施方式一相同。

具体实施方式三：如图1-5所示，本实施方式中，所述的骨折区域力学调控建模模块2 实现其功能的具体过程为：

1)在骨折区域有限元模型上施加外加载荷，并设置边界和初始条件；

所述的外加载荷的大小由骨所承受的力的大小决定，实验对象不同，所受的外加载荷也不同；

所述的边界和初始条件由CT图像信息获得；

2)将骨折区域看作双相多孔弹性模型，由多孔弹性理论得到骨折区域单元的本构方程，平衡方程和几何方程，并通过有限单元法计算骨折区域单元收到的畸变应变D和流体流速V，具体过程为：

a.本构方程

式中，σ_rr，σ_θθ，σ_zz为正应力，τ_rθ，τ_θz，τ_rz为剪应力；ε_rr，ε_θθ，ε_zz为正应变，γ_rθ，γ_θz，γ_rz为主应变；α，α'分别为各向同性弹性面的Biot系数和轴向Biot系数；p为骨痂单元中的流体压力；M₁₁，M₁₂，M₁₃，M₃₃，M₄₄，M₅₅分别为脱水的弹性模量矩阵分量；

其中，M₁₁，M₁₂，M₁₃，M₃₃，M₄₄，M₅₅脱水的弹性模量矩阵分量表达式如下所示：

M₄₄＝E_r/2(1+ν_r)>

M₅₅＝G'(7)

式中，E_r，ν_r分别是各同性弹性层的弹性模量和泊松比；E_z，ν_z分别是轴向弹性模量和泊松比；G'为剪切模量；

b.平衡方程

式中，σ_rr，σ_θθ，σ_zz为正应力，τ_rθ，τ_θz，τ_rz为剪应力；r为径向半径；

c.几何方程

式中，ε_rr，ε_θθ，ε_zz为正应变，γ_rθ，γ_θz，γ_rz为主应变；u_r，u_θ，u_z分别为三个方向上的位移；r为径向半径；

通过上述方程的求解得到骨痂单元的正应变ε_rr，ε_θθ，ε_zz，由正应变可得到骨痂单元受到的畸变应变：

式中，D为骨痂单元受到的畸变应变；ε_rr，ε_θθ，ε_zz分别为各个方向上的正应变；

骨痂单元中液体的流速V为：

其中，k为骨痂中液体的达西渗透系数；u液体粘度；p为液体压力；

其他组成及连接与具体实施方式一或二之一相同。

具体实施方式四：如图1-5所示，本实施方式中，所述的骨折区域生物调控建模模块3 实现其功能的具体过程为：

采用耦合偏微分方程组的形式对骨折区域生物调控模型进行建模，分别包括：相关细胞浓度建模、相关生长因子浓度建模和相关细胞胞外基质浓度建模；

1)相关细胞浓度建模

间充质干细胞浓度建模：

式中，c_mes为间充质干细胞浓度，t为时间，D_mes为间充质干细胞随机移动扩散系数，C_mesCT为间充质干细胞化学响应系数，g_bone为骨生长因子浓度，g_vessel为血管生成生长因子浓度，>mesHT为间充质干细胞运动系数，m为胞外基质总和，A_mes为间充质细胞增殖函数，F₁为由骨生长因子和血管生成生长因子调控的间充质干细胞向骨细胞分化的分化函数，F₂为由软骨生长因子调控的间充质干细胞向软骨细胞分化的分化函数，F₄为间充质干细胞向成纤维细胞分化的分化函数；

成纤维细胞浓度建模：

式中，c_fibro为成纤维细胞浓度，D_fibro为成纤维细胞随机扩散系数，C_fibro为成纤维细胞化学响应系数，A_fibro为成纤维细胞增殖函数，F₃为软骨骨化函数，d_fibro为成纤维细胞衰退常数；

软骨细胞浓度建模：

式中，c_cart为软骨细胞浓度，A_cart为软骨细胞增殖函数；

骨细胞浓度建模：

式中，c_bone为骨细胞浓度，C_bone为骨细胞化学响应系数，A_bone为骨细胞增殖函数，d_bone为骨细胞衰退常数；

血管生成细胞浓度建模：

式中，c_vessel为血管生成细胞浓度，D_vessel为血管生成细胞随机扩散系数，C_vesselCT为血管生成细胞化学响应系数，C_vesselHT为血管生成细胞运动系数，A_vessel为血管生成细胞增殖函数，>vessel为血管生成细胞衰退常数；

式中，D_hmes为间充质细胞最大扩散系数，K_hmes为间充质细胞最大扩散速率；

式中，C_kCTmes为间充质干细胞最大化学响应系数，K_kCTmes为间充质干细胞最大化学响应速率；

式中，C_kHTmes为间充质干细胞最大运动系数，K_kHTmes为间充质干细胞最大运动速率；

其中，i＝mes,fibro,cart,bone,vessel；

式中，A_i0为相应细胞i最大增殖系数，K_i为相应细胞i最大增殖速率；

式中，Y₁₁为与骨生长因子相关的转化常数，Y₁₂为与血管生成生长因子相关的转化常数，H₁₁为与骨生长因子相关的最低响应浓度，H₁₂为与血管生成生长因子相关的最低响应浓度；

式中，Y₂为与软骨生长因子相关的转化常数，H₂为与软骨生长因子相关的最低响应浓度；

式中，m_cart为软骨细胞胞外基质浓度，m_vessel为血管声称细胞胞外基质浓度，Y₃为软骨向骨转化的转化常数，B_ecart为软骨骨化常数，B_vessel为血管生成成熟，H₃为软骨向骨转化的最低响应浓度；

m＝m_fibro+m_cart+m_bone+m_vessel(31)

式中，m_fibro为成纤维细胞胞外基质浓度，m_bone为骨细胞胞外基质浓度；

2)相关生长因子浓度建模

软骨生长因子浓度建模：

式中，g_cart为软骨生长因子浓度，D_gcart为软骨生长因子扩散系数，E_gcart为软骨生长因子生成系数，d_gcart软骨生长因子衰退常数；

骨生长因子浓度建模：

式中，D_gbone为骨生长因子扩散系数，E_gbone为骨生长因子生成系数，d_gbone为骨生长因子衰退常数；

血管生成生长因子浓度建模：

式中，D_gvessel为血管生成生长因子扩散系数，E_gvesselbone为由骨细胞分泌血管生成生长因子系数，E_gvesselcart为由软骨细胞分泌血管生成生长因子系数，d_gvessel为血管生成生长因子衰退常数，d_gvesselcart为软骨细胞分泌血管生成生长因子衰退常数；

式中，G_gcart为软骨生长因子生成常数，H_gcart为软骨生长因子饱和常数，K_gcart为软骨细胞胞外基质增加常数；

式中，G_gbone为骨生长因子生成常数，H_gbone为骨生长因子饱和常数；

式中，G_gvesselbone为与骨相关的血管生长因子生成常数，H_gvessel血管生长因子饱和常数，>gvessel为血管生长因子胞外基质增加常数；

式中，G_gvesselcart为与软骨相关的血管生长因子生成常数；

3)相关细胞胞外基质浓度建模

成纤维细胞胞外基质浓度建模：

式中，P_fibros为成纤维细胞胞外基质合成系数，Q_fibro为成纤维细胞胞外基质降解系数；

软骨细胞胞外基质浓度建模：

式中，P_cart为软骨细胞胞外基质合成系数，Q_cart为软骨细胞胞外基质降解系数；

骨细胞胞外基质浓度建模：

式中，P_bone为骨细胞胞外基质合成系数；

血管生成细胞胞外基质浓度建模：

式中，P_vessel血管生成细胞胞外基质合成系数；

联立方程(19)-(42)构成耦合偏微分方程组，采用有限元方法求解，得到间充质干细胞浓度，成纤维细胞浓度，软骨细胞浓度，骨细胞浓度，血管生成细胞浓度，软骨生长因子浓度，骨生长因子浓度，血管生成生长因子浓度，成纤维细胞胞外基质浓度，软骨细胞胞外基质浓度，骨细胞胞外基质浓度，血管生成细胞胞外基质浓度随时间的变化关系。

其他组成及连接关系与具体实施方式一至三之一相同。

具体实施方式五：如图1-5所示，本实施方式中，所述的骨折区域力-生物联合调控建模模块4实现其功能的具体过程为：

骨折愈合过程中，与力学刺激相关的过程分别有：血管生成细胞增殖过程，血管生成细胞胞外基质合成过程，由软骨细胞分泌血管生长因子过程，骨细胞增殖过程，间充质干细胞分化为骨细胞过程，骨细胞胞外基质合成过程，软骨骨化过程，血管生成细胞凋亡；

采用梯形函数分别建立畸变应变、流体流速与上述相应过程之间的关系；定义非刺激值和刺激值，当畸变应变和流体流速未达到相应过程的刺激值时，抑制相应过程的发生；当畸变应变和流体流速达到相应过程的刺激值时，促进相应过程的发生；

其他组成及连接关系与具体实施方式一至四之一相同。

具体实施方式六：如图1-5所示，本实施方式中，所述的程序终止判断模块5实现其功能的具体过程为：

1)判断骨折区域单元材料属性

判断当前骨折区域单元材料属性是否与骨的材料属性相同，若不相等，则程序执行2) 步骤；若相等，程序执行3)步骤；

2)更新骨折区域单元材料属性

若骨折区域单元材料属性不等于骨的材料属性，则骨折区域单元材料属性进行更新，并进入下一个迭代步，骨折区域材料属性更新公式如下：

式中，E_ele,i为单元i弹性模量，m_fibro,i为单元i成纤维细胞胞外基质浓度，为成纤维细胞胞外基质参考浓度，E_fibro为成纤维细胞胞外基质弹性模量，m_cart,i为单元i软骨细胞胞外基质浓度，为软骨细胞胞外基质参考浓度，E_cart为软骨细胞胞外基质弹性模量，m_bone,i为单元i骨细胞胞外基质浓度，为骨细胞胞外基质参考浓度，E_bone为骨细胞胞外基质弹性模量；

式中，ν_ele,i为单元i泊松比，ν_fibro为成纤维细胞胞外基质泊松比，ν_cart为软骨细胞胞外基质泊松比，ν_bone为骨细胞胞外基质泊松比；

式中，μ_ele,i为单元i渗透率，μ_fibro为成纤维细胞胞外基质渗透率，μ_cart为软骨细胞胞外基质渗透率，μ_bone为骨细胞胞外基质渗透率；

3)程序结束

若所有骨折区域单元材料属性等于骨的材料属性，则程序结束，记录骨折愈合时间，绘制相关细胞浓度、相关生长因子浓度及相关细胞胞外基质浓度随时间变化的云图。

其他组成及连接关系与具体实施方式一至五之一相同。

实施例：

为了说明本系统的使用方法，下面具体举一个例子阐述本发明的操作过程。

模拟羊跖骨骨折愈合过程

1.建立羊跖骨分析模型

1)建立羊跖骨几何模型

将CT图像数据利用三维医学图像表面重建算法对图像进行三维表面重构，得到羊跖骨骨折区域三维几何模型。

2)建立羊跖骨有限元分析模型

将上述三维几何模型导入网格划分软件中进行网格划分，将得到的有限元模型导入到 matlab中进行预处理，只提取目标数据，根据目标数据生成后续有限元计算需要的单元编号和节点坐标文件。单元编号和节点坐标两个文件为txt文本格式的文件。单元编号文件包含四列数据，四列数据分别为每个单元的四个节点序号。节点坐标文件包含三列数据，三列数据分别为每个节点的空间坐标值。

2.羊跖骨力学模型的建立

1)施加边界条件及材料属性

施加边界条件：对羊跖骨有限元分析模型底部施加完全固定约束。

施加材料属性：羊跖骨骨折区域材料属性赋值包括皮质骨材料属性赋值和骨痂单元材料属性赋值。将骨折区域看成是基于混合物的双相多孔弹性模型，皮质骨和骨痂单元材料属性如表1所示。

2)骨折区域力学刺激有限元分析

对建立好的羊跖骨有限元模型施加载荷。由于羊在正常行走时，跖骨所受力的大小为 500N，所以在骨折区域施加大小为500N的载荷。

对羊跖骨骨折区域进行有限元分析，得到骨折区域单元受到的畸变应变和流体流速。

表1各组织材料属性

3.羊跖骨生物学模型建模

采用偏微分方程组建立羊跖骨生物学模型，经有限元方法得到相关细胞浓度，相关生长因子浓度和相关细胞胞外基质浓度随时间的变化关系。

4.羊跖骨力-生物联合调控建模

采用梯形函数建立畸变应变和流体流速与与力相关的生物过程之间的关系。当力未达到生物过程刺激值时，该过程抑制；当力达到刺激值时，该过程发生。具体梯形图如图5所示，具体关系如下表所示：

表二生物过程参数的非刺激值与刺激值

5.判断程序终止条件

判断当前骨折区域单元是否等于骨的材料属性，若不等，则更新当前材料属性，程序进入下一个迭代步；若相等，则程序结束，记录愈合时间，绘制相关细胞浓度，相关生长因子浓度和相关细胞胞外基质浓度随时间的变化云图。