基于梯度函數的仿生骨支架建模與評價

李 穎，章蘭珠

(華東理工大學機械與動力工程學院，上海 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將梯度函數定義如下

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)微觀孔結構負模型的生成可由程序控制生成，提高了建模效率和自動化程度，為實現計算機輔助組織工程奠定了理論和技術基礎。