[发明专利]一种基于力学刺激和生物联合刺激的骨折愈合仿真系统

权利要求书

1.一种基于力学刺激和生物联合刺激的骨折愈合仿真系统，其特征在于，所述系统包括：

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

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

所述的骨折区域分析模型建立模块(1)实现其功能的具体过程为：

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

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

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

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

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

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

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

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

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

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

所述的骨折区域力学调控建模模块(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) (6)

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为液体压力；

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

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

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

所述的骨折区域力-生物联合调控建模模块(4)实现其功能的具体过程为：

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

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

程序终止判断模块(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)程序结束

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

2.根据权利要求1所述的一种基于力学刺激和生物联合刺激的骨折愈合仿真系统，其特征在于：所述的骨折区域生物调控建模模块(3)实现其功能的具体过程为：

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

1)相关细胞浓度建模

间充质干细胞浓度建模：

式中，c_mes为间充质干细胞浓度，t为时间，D_mes为间充质干细胞随机移动扩散系数，C_mesCT为间充质干细胞化学响应系数，g_bone为骨生长因子浓度，g_vessel为血管生成生长因子浓度，C_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为血管生成细胞增殖函数，d_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₂为与软骨生长因子相关的最低响应浓度，g_cart为软骨生长因子浓度；

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

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

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

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

软骨生长因子浓度建模：

式中，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血管生成细胞生长因子饱和常数，K_gvessel为血管生成细胞生长因子胞外基质增加常数；

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

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

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

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

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

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

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

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

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

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

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