基于梯度函数的仿生骨支架建模与评价

李 颖，章兰珠

(华东理工大学机械与动力工程学院，上海 200237)

1 引言

近年来，由于交通事故、自身疾病等导致的大面积骨缺损，已成为临床医学治疗的重要课题之一。作为骨组织体外培养的重要载体，仿生骨支架为种子细胞的黏附和增殖提供生存空间，为细胞获取营养和新陈代谢提供通道，在成骨阶段为组织提供必要的力学支撑[1]。因此，对于仿生骨支架的构建必须考虑其微观孔结构、力学性能、生物活性及打印材料降解速度等因素。

饶嵩[2](2004年)提出了采用分形理论和蒙特卡洛法构建支架内部微观孔结构的建模方法；杨楠[3](2006年)提出基于知识的人工骨三维仿生设计，将知识库、分形理论和体视学相结合；以上两种方法所获得的微观孔结构为随机均匀分布，缺乏梯度变化规律，且无法控制微孔的位置，可操控性低。Cai S[4]( 2008年)等人提出了基于形状函数和共形全六面体网格细化的骨支架设计方法。但该方法仍缺乏梯度变化规律。蔡升勇(2009年)[5-6]提出了一种利用形函数进行组织工程骨架孔隙建模的方法，该方法虽然得到了具有一定梯度的孔隙模型，但其梯度性仍缺乏理论及实际数据的支持，不具有普适性。尤飞(2010年)[7]基于多约束背包问题模型，利用混合遗传算法构建仿生骨支架模型，该方法构建出的微观孔负模型仍为随机均匀分布，缺乏梯度变化规律。Gómez， S[8](2016年)等人运用Voronoi曲面细分方法进行骨骼三维多孔支架的设计。史建平[9](2018年)等人提出基于TPMS的仿生骨组织工程多孔支架建模方法。TPMS方法构建的多孔骨支架同时具备了精确的外轮廓和复杂的内部结构，目前正在成为支架内孔精细化建模的重要选择[10]，但对于骨骼梯度化的构造过渡较不自然，存在断层，衔接性有待改进。

本文首先对实际牛骨进行Micro CT扫描，通过CTan软件分析扫描数据，获得其微观结构分布规律后利用Matlab软件构建出相应的梯度函数；然后利用遗传算法，对梯度函数中微孔的位置进行适当的偏差处理，得到既拟合实际骨截面梯度分布规律又具有一定随机性的微观孔支架负模型；最后将实体模型与微观孔负模型进行布尔减运算，得到具有梯度变化规律的仿生骨支架模型。

2 基于梯度函数的仿生骨支架建模

2.1 微观孔结构的单元体模型

对于微观孔结构单元体模型的设计，其基本原则是根据实际骨骼内部微观孔的结构进行仿生设计。在对实际骨骼进行Micro CT扫描后，将得到的微观孔结构灰度图像进行进一步的二值化处理，并采用分形理论对不规则曲线进行统计分析[11]，如图1，最终选择椭球体作为微观孔结构的单元体模型，将其模型定义为八元组模型[1]

ELL=(a，b，x，y，z，θx，θy，θz)

(1)

其中，a——椭球体的长轴

b——椭球体的短轴

x，y，z——椭球体的中心点坐标

θx，θy，θz——椭球体长轴与三个坐标平面的夹角

图1 微观孔结构截面曲线拟合

2.2 Micro CT 实物扫描及分析

选择与人体具有较高相似性的牛骨进行实际数据的获取，在对材料进行处理和筛选后，最终确定选用4块牛脊柱骨作为扫描样本，如图2所示。经过Micro CT扫描后，利用CTAn软件对扫描所得的数据图像集，进行感兴趣区域(ROI，Region of Interest)分析，查看和记录建模所需的定量参数(如图3、表1所示)。

图2 牛骨样品

图3 Tb.Sp和Tb.Th变化规律

表1 ROI数据分析

2.3 微观孔结构的梯度函数构建

为了构建符合实际骨骼内部微观结构的模型，本文提出通过构建拟合实际骨骼内部孔隙分布规律的梯度函数，控制模型中椭球的大小及位置情况，使椭球的分布规律与实际骨骼内部孔隙分布规律相同。梯度函数是描述椭球中心点位置分布和椭球个数及长轴大小随中心点位置变化而变化的函数。现将梯度函数定义如下

F(x)=fa(fr(x))

(2)

式中，F(x)——微观孔结构的总体梯度函数

fa(x)——椭球长轴的梯度函数

fr(x)——椭球位置的梯度函数

假设，仿生骨支架为圆柱体包容盒，则微观孔结构的负模型构建步骤如图4所示。

图4 椭球体的位置生成过程

a)将其沿圆柱高度方向划分为m层，划分宽度为椭球体短轴大小，为避免同一截面椭球出现中间连通而边缘不连通的情况，规定椭球短轴大小为长轴的最小值，即b=amin。则椭球体Z轴方向坐标可表示为：

(3)

式中，zl——第i层椭球的Z轴坐标

H——圆柱体包容盒高度

b——椭球体的短轴

m——椭球体层数

b)再将每一层截面划分为n个同心圆(此为不均等分)，每圈同心圆上的椭球体长轴都相同但椭球体个数由同心圆半径r的大小而确定，其同心圆半径r由椭球位置的梯度函数fr(x)确定，fr(x)具体表达式如下

fr(xi)=∑fr(xi-1) +fa(xi-1)+fth(xi)

-fΔ(x)，i=1，2，…，n

(4)

fa(x)=f(Ri，Tb.Spi)，i=0，1，2，…，n-1

(5)

fth(x)=f(Ri，Tb.Thi)，i=0，1，2，…，n-1

(6)

(7)

式中，fth(x)——椭球体间厚度的梯度函数

fΔ(x)——椭球体间孔径差

R——ROI分析数据中横截面距骨骼中心截面的距离

Tb.Sp——ROI分析数据中骨小梁的分离度值

Tb.Th——ROI分析数据中骨小梁的厚度值

n——不同长轴的椭球体个数，即每层截面上同心圆个数

式中椭球体间厚度的梯度函数由ROI分析中的数据拟合所得，该微观结构模型首先从包容盒内靠近圆心位置处生成一组椭球体，其长轴为梯度函数最大值，再以此椭球体为基准，其位置半径下骨小梁对应的厚度为椭球体间间距，再加上前后两椭球体长轴的二分之一，这三段距离与上一组椭球位置半径的和作为下一组椭球的位置半径r，如图5所示。其中，下一组椭球的长轴以上一组椭球减去椭球体间孔径差计算而得，存在的误差忽略不计。

图5 椭球位置半径的确定

c)不同同心圆位置处的长轴大小由椭球长轴的梯度函数fa(x)确定，该函数由ROI分析中的数据拟合所得；

d)最后将每一层截面的每一同心圆一次划分为k1，k2，…，kn个区域(n为同心圆个数)，每一区域放置一个椭球体，其区域个数和各椭球体位置的确定如下(如图6所示)：

(8)

(9)

式中，kn——第n个同心圆上椭球体的个数

θj——第j个椭球体的角度位置

图6 椭球角度位置的确定

由此，某一截面上的椭球体的中心点位置坐标可由(ricosθj，risinθj，zl)唯一确定。对于其姿态参数(θx，θy，θz)，由于没有固定约束，所以随机生成即可。由此，现可得圆柱体包容盒内椭球体个数N为

(10)

基于以上步骤得出的微观孔结构负模型，但相比实际骨骼该模型整体较为规整，缺乏随机性，所以需为每个椭球体增加随机偏移量，在符合梯度分布的基础上增加随机性。

3 基于遗传算法的孔隙空间位置调整

在得到具有梯度分布规律的微观孔结构负模型之后，本文进一步运用遗传算法对每个椭球单元体的中心点坐标设置偏置值进行位置调整，使最终得到的仿生骨支架模型满足期望的孔隙空间分布。对于偏置值的设置流程如图7所示。

图7 基于遗传算法的孔隙空间位置调整

偏置值的设置流程具体如下：

步骤1：将建好的具有梯度分布规律的负模型作为输入；

步骤2：对负模型进行染色体浮点数编码，将椭球体位置调整问题转化为遗传算法求解最有个体问题；

步骤3：确定迭代次数和染色体个数，通过种群初始化和染色体的选择、交叉、变异和扰动等遗传操作对问题进行求解，并用适应度函数评估产生的多组新个体；

步骤4：重复步骤3直至完成迭代次数，输出最优个体，即适应度值最高的个体。

为了使微观孔结构的负模型在梯度分布的基础上具有一定的随机性，构建出更拟合实际骨骼的模型，本文对每个椭球体中心点的位置调整如下

(11)

(12)

(13)

Δr(xi)∈[-α·fth(xi)，α·fth(xi)]，α∈(0，1)

(14)

Δθ(xj)∈[-β·θj，β·θj]，β∈(0，1)

(15)

Δr(xi)——第i圈同心圆半径的偏置值

Δθ(xj)——第j个椭球角度的偏置值

α——同心圆半径偏置系数

β——椭球角度偏置系数

理论上，可以对椭球体模型的八个要素均设置偏置值，但考虑到该模型是根据实际骨骼的梯度函数建立的，为了降低分布调整对梯度分布规律的影响，本文决定只对椭球体中心点位置进行适当调整，椭球体长短轴仍保持由梯度函数计算得出的数值，姿态参数仍保持随机生成的原数值。对于椭球体中心点位置的调整范围，应保持在骨小梁厚度范围内，避免大跨度的位置变动，如图8所示。

图8 孔隙空间位置调整

4 系统开发与例证

由于微观孔结构负模型是由一定数量的椭球体构成，每个椭球体包含8个描述参数，所以由遗传算法得到的最优个体的数据解集合较为庞大。为了提高建模效率和自动化程度，减少人机交互操作，本文利用Solidworks二次开发工具，编写了微观孔结构负模型自动生成程序。

本文选择圆柱体作为包容盒外形，首先在包容盒范围内按照ROI分析得到的实际骨骼的梯度函数初步生成微观孔结构负模型数据集，再运用遗传算法对负模型中微观孔的空间位置进行适当调整，选择最优个体输出数据集，将该数据集链接到微观孔结构负模型自动生成程序完成负模型的生成，最后将包容盒实体模型与微观孔结构负模型进行布尔减运算，得到仿生骨支架模型。具体流程及结果可视化如图9所示。

图9 仿生骨支架建模流程

该仿生骨支架模型以1：2.8的比例还原原试样，图10所示为微观孔负模型优化前后对比图，可以看出：1)在由梯度函数得到的微观孔结构负模型中，随着位置半径的增大其椭球体长轴逐渐减小，椭球体间厚度逐渐增大；2)运用遗传算法适当调整椭球体中心点位置后，其微观孔空间分布在梯度变化的基础上增加了一定的随机性，使整体结构中有序和无序并存。

图10 微观孔负模型优化前后对比

图11为实际骨骼与仿生骨支架模型的截面对比图，从图中可以看出本文所建仿生骨支架模型与实际骨骼存在一定的相似性，分布规律基本相同。图12为仿生骨支架模型的微观结构梯度变化示意图，从图中可以看出：1)骨小梁厚度由中心向两边逐渐增大；2)随着位置半径逐渐增大孔隙逐渐减小；

图11 实际骨骼与仿生骨支架截面对比

图12 仿生骨支架模型微观结构梯度变化

基于以上步骤和梯度函数，建立力学拉伸实验模型，对该方法构建的仿生骨支架进行力学强度等方面的试验。建模结果如图13所示，使用钛合金材料对模型进行3D打印，实物如图14所示。

图13 实验模型

图14 钛合金3D打印

5 仿生骨支架模型微观结构的评价

仿生骨支架模型的微观孔结构对体内营养的输送、骨细胞的黏附、生长和仿生骨支架的力学性能息息相关，而对于微观孔结构的影响因素主要包括：支架孔隙率、孔之间的连通性、微孔分布的梯度性和均匀性等。

5.1 孔隙率

孔隙率是描述骨支架内部孔隙含量的参数，孔隙作为人体内营养输送和骨细胞黏附的空间载体，较大时易形成骨质疏松降低力学性能，较小时不利于骨组织长入诱导新骨形成，所以孔隙率是描述支架性能的重要参数之一。根据定义可得孔隙率的计算公式如下：

(17)

式中，δ——支架孔隙率

Ve——支架内孔隙体积

Vt——包容盒体积

其中，支架内孔隙体积可通过ROI分析计算得出。

5.2 梯度拟合度

梯度拟合度是本文对仿生骨支架评价中的新增参数，用来表示所构建的支架模型与实际骨骼梯度变化规律的拟合程度。拟合值越高说明对实际骨骼模型的建模结果越贴合实际情况，梯度拟合度的计算方法与遗传算法中适应度值的计算方法相同

(18)

式中，φ——支架梯度拟合度

Ai——实际骨骼孔隙分布向量

5.3 连通性

连通性是描述骨支架内部微观孔互相连接导通情况的参数，仿生骨支架必须具有一定的连通性才可以保证骨组织和血管的长入、营养物质的输送和细胞代谢物的排出。

本文对连通性的计算定义如下

(19)

式中，η——支架连通率

N——椭球体总个数

ni——每个椭球周围连通的其他椭球个数

对于椭球体之间的连通情况，主要分为三种：a)独立椭球体(η=0)；b)两个椭球体导通(η=1)，由独立椭球体生成的盲孔和两个椭球体生成的单通孔，这两种情况较多存在于实际骨骼中处于边缘位置的部分；但在构建仿生骨负模型时要尽量减少盲孔和单通孔的存在；c)三个或三个以上椭球导通(η>=2)，该种情况为最佳情况，可保证仿生骨支架内部的各椭球之间大概率导通。

5.4 指标得分

基于多属性决策模型，本文根据以上评价指标提出对模型优劣的打分制度：

步骤1：分别对基于随机均分分布的仿生骨支架模型、基于梯度函数分布的仿生骨支架和基于梯度函数分布但椭球空间位置未调整的仿生骨支架模型进行评价指标的计算，得到原始属性值(见表2)，即孔隙率、梯度拟合性和连通性得到原始决策矩阵E。

表2 仿生骨支架评价

矩阵的行代表模型类别，列代表属性；

步骤2：对原始属性值进行归一化处理，得到新的决策矩阵E′

步骤3：根据成对比较阵标度表，构建各属性的成对比较阵F，并运用SPSS软件计算各属性权重G

G=[0.196， 0.311， 0.493]

即，孔隙率、梯度拟合度和连通性的属性权重依次为0.196，0.311，0.493。

步骤4：由所得的属性权重G和归一化处理后的属性值E′，计算模型得分Si

(20)

由表3结果得：1)基于梯度函数构建的仿生骨支架模型比随即均匀分布的模型在各指标上均具有明显优势；2)通过遗传算法进行空间位置调整后的模型在连通性和孔隙率方面均有一定的进步性；

综上，本文所构建的基于梯度函数和遗传算法的仿生骨支架模型为最优模型。

表3 各仿生骨支架得分

6 结论

1)本文首先基于对实际骨骼进行的Micro CT扫描和ROI分析，构建了拟合实际骨骼微观孔结构分布规律的梯度函数。以椭球体作为仿生骨支架微观结构的单元体模型，圆柱体外形为包容盒，根据得到的梯度函数控制椭球体的空间位置和长轴分布，之后再运用遗传算法对椭球体的中心点位置进行适当调整，得到微观孔结构负模型，最后通过布尔减运算得到仿生骨支架模型；

2)本文为仿生骨支架建模方法提出了一条新的技术路线和建模方法。该方法对所有不同形态和分布的骨骼均适用，且以实际骨骼数据为依据进行建模，具有一定的科学依据和参考价值；

3)以仿生骨支架的孔隙率、梯度拟合度和连通性作为评价指标，并基于多属性决策模型提出打分制度，从而建立出仿生骨支架微观结构的评价体系，并提出了各指标与打分制度的计算方法，以实现对所建模型的评价和优化；

4)微观孔结构负模型的生成可由程序控制生成，提高了建模效率和自动化程度，为实现计算机辅助组织工程奠定了理论和技术基础。