李 穎，章蘭珠

(華東理工大學(xué)機(jī)械與動(dòng)力工程學(xué)院，上海 200237)

1 引言

近年來，由于交通事故、自身疾病等導(dǎo)致的大面積骨缺損，已成為臨床醫(yī)學(xué)治療的重要課題之一。作為骨組織體外培養(yǎng)的重要載體，仿生骨支架為種子細(xì)胞的黏附和增殖提供生存空間，為細(xì)胞獲取營(yíng)養(yǎng)和新陳代謝提供通道，在成骨階段為組織提供必要的力學(xué)支撐[1]。因此，對(duì)于仿生骨支架的構(gòu)建必須考慮其微觀孔結(jié)構(gòu)、力學(xué)性能、生物活性及打印材料降解速度等因素。

饒嵩[2](2004年)提出了采用分形理論和蒙特卡洛法構(gòu)建支架內(nèi)部微觀孔結(jié)構(gòu)的建模方法；楊楠[3](2006年)提出基于知識(shí)的人工骨三維仿生設(shè)計(jì)，將知識(shí)庫(kù)、分形理論和體視學(xué)相結(jié)合；以上兩種方法所獲得的微觀孔結(jié)構(gòu)為隨機(jī)均勻分布，缺乏梯度變化規(guī)律，且無法控制微孔的位置，可操控性低。Cai S[4]( 2008年)等人提出了基于形狀函數(shù)和共形全六面體網(wǎng)格細(xì)化的骨支架設(shè)計(jì)方法。但該方法仍缺乏梯度變化規(guī)律。蔡升勇(2009年)[5-6]提出了一種利用形函數(shù)進(jìn)行組織工程骨架孔隙建模的方法，該方法雖然得到了具有一定梯度的孔隙模型，但其梯度性仍缺乏理論及實(shí)際數(shù)據(jù)的支持，不具有普適性。尤飛(2010年)[7]基于多約束背包問題模型，利用混合遺傳算法構(gòu)建仿生骨支架模型，該方法構(gòu)建出的微觀孔負(fù)模型仍為隨機(jī)均勻分布，缺乏梯度變化規(guī)律。Gómez， S[8](2016年)等人運(yùn)用Voronoi曲面細(xì)分方法進(jìn)行骨骼三維多孔支架的設(shè)計(jì)。史建平[9](2018年)等人提出基于TPMS的仿生骨組織工程多孔支架建模方法。TPMS方法構(gòu)建的多孔骨支架同時(shí)具備了精確的外輪廓和復(fù)雜的內(nèi)部結(jié)構(gòu)，目前正在成為支架內(nèi)孔精細(xì)化建模的重要選擇[10]，但對(duì)于骨骼梯度化的構(gòu)造過渡較不自然，存在斷層，銜接性有待改進(jìn)。

本文首先對(duì)實(shí)際牛骨進(jìn)行Micro CT掃描，通過CTan軟件分析掃描數(shù)據(jù)，獲得其微觀結(jié)構(gòu)分布規(guī)律后利用Matlab軟件構(gòu)建出相應(yīng)的梯度函數(shù)；然后利用遺傳算法，對(duì)梯度函數(shù)中微孔的位置進(jìn)行適當(dāng)?shù)钠钐幚?，得到既擬合實(shí)際骨截面梯度分布規(guī)律又具有一定隨機(jī)性的微觀孔支架負(fù)模型；最后將實(shí)體模型與微觀孔負(fù)模型進(jìn)行布爾減運(yùn)算，得到具有梯度變化規(guī)律的仿生骨支架模型。

2 基于梯度函數(shù)的仿生骨支架建模

2.1 微觀孔結(jié)構(gòu)的單元體模型

對(duì)于微觀孔結(jié)構(gòu)單元體模型的設(shè)計(jì)，其基本原則是根據(jù)實(shí)際骨骼內(nèi)部微觀孔的結(jié)構(gòu)進(jìn)行仿生設(shè)計(jì)。在對(duì)實(shí)際骨骼進(jìn)行Micro CT掃描后，將得到的微觀孔結(jié)構(gòu)灰度圖像進(jìn)行進(jìn)一步的二值化處理，并采用分形理論對(duì)不規(guī)則曲線進(jìn)行統(tǒng)計(jì)分析[11]，如圖1，最終選擇橢球體作為微觀孔結(jié)構(gòu)的單元體模型，將其模型定義為八元組模型[1]

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

(1)

其中，a——橢球體的長(zhǎng)軸

b——橢球體的短軸

x，y，z——橢球體的中心點(diǎn)坐標(biāo)

θx，θy，θz——橢球體長(zhǎng)軸與三個(gè)坐標(biāo)平面的夾角

圖1 微觀孔結(jié)構(gòu)截面曲線擬合

2.2 Micro CT 實(shí)物掃描及分析

選擇與人體具有較高相似性的牛骨進(jìn)行實(shí)際數(shù)據(jù)的獲取，在對(duì)材料進(jìn)行處理和篩選后，最終確定選用4塊牛脊柱骨作為掃描樣本，如圖2所示。經(jīng)過Micro CT掃描后，利用CTAn軟件對(duì)掃描所得的數(shù)據(jù)圖像集，進(jìn)行感興趣區(qū)域(ROI，Region of Interest)分析，查看和記錄建模所需的定量參數(shù)(如圖3、表1所示)。

圖2 牛骨樣品

圖3 Tb.Sp和Tb.Th變化規(guī)律

表1 ROI數(shù)據(jù)分析

2.3 微觀孔結(jié)構(gòu)的梯度函數(shù)構(gòu)建

為了構(gòu)建符合實(shí)際骨骼內(nèi)部微觀結(jié)構(gòu)的模型，本文提出通過構(gòu)建擬合實(shí)際骨骼內(nèi)部孔隙分布規(guī)律的梯度函數(shù)，控制模型中橢球的大小及位置情況，使橢球的分布規(guī)律與實(shí)際骨骼內(nèi)部孔隙分布規(guī)律相同。梯度函數(shù)是描述橢球中心點(diǎn)位置分布和橢球個(gè)數(shù)及長(zhǎng)軸大小隨中心點(diǎn)位置變化而變化的函數(shù)?，F(xiàn)將梯度函數(shù)定義如下

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

(2)

式中，F(xiàn)(x)——微觀孔結(jié)構(gòu)的總體梯度函數(shù)

fa(x)——橢球長(zhǎng)軸的梯度函數(shù)

fr(x)——橢球位置的梯度函數(shù)

假設(shè)，仿生骨支架為圓柱體包容盒，則微觀孔結(jié)構(gòu)的負(fù)模型構(gòu)建步驟如圖4所示。

圖4 橢球體的位置生成過程

a)將其沿圓柱高度方向劃分為m層，劃分寬度為橢球體短軸大小，為避免同一截面橢球出現(xiàn)中間連通而邊緣不連通的情況，規(guī)定橢球短軸大小為長(zhǎng)軸的最小值，即b=amin。則橢球體Z軸方向坐標(biāo)可表示為：

(3)

式中，zl——第i層橢球的Z軸坐標(biāo)

H——圓柱體包容盒高度

b——橢球體的短軸

m——橢球體層數(shù)

b)再將每一層截面劃分為n個(gè)同心圓(此為不均等分)，每圈同心圓上的橢球體長(zhǎng)軸都相同但橢球體個(gè)數(shù)由同心圓半徑r的大小而確定，其同心圓半徑r由橢球位置的梯度函數(shù)fr(x)確定，fr(x)具體表達(dá)式如下

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)——橢球體間厚度的梯度函數(shù)

fΔ(x)——橢球體間孔徑差

R——ROI分析數(shù)據(jù)中橫截面距骨骼中心截面的距離

Tb.Sp——ROI分析數(shù)據(jù)中骨小梁的分離度值

Tb.Th——ROI分析數(shù)據(jù)中骨小梁的厚度值

n——不同長(zhǎng)軸的橢球體個(gè)數(shù)，即每層截面上同心圓個(gè)數(shù)

式中橢球體間厚度的梯度函數(shù)由ROI分析中的數(shù)據(jù)擬合所得，該微觀結(jié)構(gòu)模型首先從包容盒內(nèi)靠近圓心位置處生成一組橢球體，其長(zhǎng)軸為梯度函數(shù)最大值，再以此橢球體為基準(zhǔn)，其位置半徑下骨小梁對(duì)應(yīng)的厚度為橢球體間間距，再加上前后兩橢球體長(zhǎng)軸的二分之一，這三段距離與上一組橢球位置半徑的和作為下一組橢球的位置半徑r，如圖5所示。其中，下一組橢球的長(zhǎng)軸以上一組橢球減去橢球體間孔徑差計(jì)算而得，存在的誤差忽略不計(jì)。

圖5 橢球位置半徑的確定

c)不同同心圓位置處的長(zhǎng)軸大小由橢球長(zhǎng)軸的梯度函數(shù)fa(x)確定，該函數(shù)由ROI分析中的數(shù)據(jù)擬合所得；

d)最后將每一層截面的每一同心圓一次劃分為k1，k2，…，kn個(gè)區(qū)域(n為同心圓個(gè)數(shù))，每一區(qū)域放置一個(gè)橢球體，其區(qū)域個(gè)數(shù)和各橢球體位置的確定如下(如圖6所示)：

(8)

(9)

式中，kn——第n個(gè)同心圓上橢球體的個(gè)數(shù)

θj——第j個(gè)橢球體的角度位置

圖6 橢球角度位置的確定

由此，某一截面上的橢球體的中心點(diǎn)位置坐標(biāo)可由(ricosθj，risinθj，zl)唯一確定。對(duì)于其姿態(tài)參數(shù)(θx，θy，θz)，由于沒有固定約束，所以隨機(jī)生成即可。由此，現(xiàn)可得圓柱體包容盒內(nèi)橢球體個(gè)數(shù)N為

(10)

基于以上步驟得出的微觀孔結(jié)構(gòu)負(fù)模型，但相比實(shí)際骨骼該模型整體較為規(guī)整，缺乏隨機(jī)性，所以需為每個(gè)橢球體增加隨機(jī)偏移量，在符合梯度分布的基礎(chǔ)上增加隨機(jī)性。

3 基于遺傳算法的孔隙空間位置調(diào)整

在得到具有梯度分布規(guī)律的微觀孔結(jié)構(gòu)負(fù)模型之后，本文進(jìn)一步運(yùn)用遺傳算法對(duì)每個(gè)橢球單元體的中心點(diǎn)坐標(biāo)設(shè)置偏置值進(jìn)行位置調(diào)整，使最終得到的仿生骨支架模型滿足期望的孔隙空間分布。對(duì)于偏置值的設(shè)置流程如圖7所示。

圖7 基于遺傳算法的孔隙空間位置調(diào)整

偏置值的設(shè)置流程具體如下：

步驟1：將建好的具有梯度分布規(guī)律的負(fù)模型作為輸入；

步驟2：對(duì)負(fù)模型進(jìn)行染色體浮點(diǎn)數(shù)編碼，將橢球體位置調(diào)整問題轉(zhuǎn)化為遺傳算法求解最有個(gè)體問題；

步驟3：確定迭代次數(shù)和染色體個(gè)數(shù)，通過種群初始化和染色體的選擇、交叉、變異和擾動(dòng)等遺傳操作對(duì)問題進(jìn)行求解，并用適應(yīng)度函數(shù)評(píng)估產(chǎn)生的多組新個(gè)體；

步驟4：重復(fù)步驟3直至完成迭代次數(shù)，輸出最優(yōu)個(gè)體，即適應(yīng)度值最高的個(gè)體。

為了使微觀孔結(jié)構(gòu)的負(fù)模型在梯度分布的基礎(chǔ)上具有一定的隨機(jī)性，構(gòu)建出更擬合實(shí)際骨骼的模型，本文對(duì)每個(gè)橢球體中心點(diǎn)的位置調(diào)整如下

(11)

(12)

(13)

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

(14)

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

(15)

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

Δθ(xj)——第j個(gè)橢球角度的偏置值

α——同心圓半徑偏置系數(shù)

β——橢球角度偏置系數(shù)

理論上，可以對(duì)橢球體模型的八個(gè)要素均設(shè)置偏置值，但考慮到該模型是根據(jù)實(shí)際骨骼的梯度函數(shù)建立的，為了降低分布調(diào)整對(duì)梯度分布規(guī)律的影響，本文決定只對(duì)橢球體中心點(diǎn)位置進(jìn)行適當(dāng)調(diào)整，橢球體長(zhǎng)短軸仍保持由梯度函數(shù)計(jì)算得出的數(shù)值，姿態(tài)參數(shù)仍保持隨機(jī)生成的原數(shù)值。對(duì)于橢球體中心點(diǎn)位置的調(diào)整范圍，應(yīng)保持在骨小梁厚度范圍內(nèi)，避免大跨度的位置變動(dòng)，如圖8所示。

圖8 孔隙空間位置調(diào)整

4 系統(tǒng)開發(fā)與例證

由于微觀孔結(jié)構(gòu)負(fù)模型是由一定數(shù)量的橢球體構(gòu)成，每個(gè)橢球體包含8個(gè)描述參數(shù)，所以由遺傳算法得到的最優(yōu)個(gè)體的數(shù)據(jù)解集合較為龐大。為了提高建模效率和自動(dòng)化程度，減少人機(jī)交互操作，本文利用Solidworks二次開發(fā)工具，編寫了微觀孔結(jié)構(gòu)負(fù)模型自動(dòng)生成程序。

本文選擇圓柱體作為包容盒外形，首先在包容盒范圍內(nèi)按照ROI分析得到的實(shí)際骨骼的梯度函數(shù)初步生成微觀孔結(jié)構(gòu)負(fù)模型數(shù)據(jù)集，再運(yùn)用遺傳算法對(duì)負(fù)模型中微觀孔的空間位置進(jìn)行適當(dāng)調(diào)整，選擇最優(yōu)個(gè)體輸出數(shù)據(jù)集，將該數(shù)據(jù)集鏈接到微觀孔結(jié)構(gòu)負(fù)模型自動(dòng)生成程序完成負(fù)模型的生成，最后將包容盒實(shí)體模型與微觀孔結(jié)構(gòu)負(fù)模型進(jìn)行布爾減運(yùn)算，得到仿生骨支架模型。具體流程及結(jié)果可視化如圖9所示。

圖9 仿生骨支架建模流程

該仿生骨支架模型以1：2.8的比例還原原試樣，圖10所示為微觀孔負(fù)模型優(yōu)化前后對(duì)比圖，可以看出：1)在由梯度函數(shù)得到的微觀孔結(jié)構(gòu)負(fù)模型中，隨著位置半徑的增大其橢球體長(zhǎng)軸逐漸減小，橢球體間厚度逐漸增大；2)運(yùn)用遺傳算法適當(dāng)調(diào)整橢球體中心點(diǎn)位置后，其微觀孔空間分布在梯度變化的基礎(chǔ)上增加了一定的隨機(jī)性，使整體結(jié)構(gòu)中有序和無序并存。

圖10 微觀孔負(fù)模型優(yōu)化前后對(duì)比

圖11為實(shí)際骨骼與仿生骨支架模型的截面對(duì)比圖，從圖中可以看出本文所建仿生骨支架模型與實(shí)際骨骼存在一定的相似性，分布規(guī)律基本相同。圖12為仿生骨支架模型的微觀結(jié)構(gòu)梯度變化示意圖，從圖中可以看出：1)骨小梁厚度由中心向兩邊逐漸增大；2)隨著位置半徑逐漸增大孔隙逐漸減??；

圖11 實(shí)際骨骼與仿生骨支架截面對(duì)比

圖12 仿生骨支架模型微觀結(jié)構(gòu)梯度變化

基于以上步驟和梯度函數(shù)，建立力學(xué)拉伸實(shí)驗(yàn)?zāi)Ｐ?，?duì)該方法構(gòu)建的仿生骨支架進(jìn)行力學(xué)強(qiáng)度等方面的試驗(yàn)。建模結(jié)果如圖13所示，使用鈦合金材料對(duì)模型進(jìn)行3D打印，實(shí)物如圖14所示。

圖13 實(shí)驗(yàn)?zāi)Ｐ?/p>

圖14 鈦合金3D打印

5 仿生骨支架模型微觀結(jié)構(gòu)的評(píng)價(jià)

仿生骨支架模型的微觀孔結(jié)構(gòu)對(duì)體內(nèi)營(yíng)養(yǎng)的輸送、骨細(xì)胞的黏附、生長(zhǎng)和仿生骨支架的力學(xué)性能息息相關(guān)，而對(duì)于微觀孔結(jié)構(gòu)的影響因素主要包括：支架孔隙率、孔之間的連通性、微孔分布的梯度性和均勻性等。

5.1 孔隙率

孔隙率是描述骨支架內(nèi)部孔隙含量的參數(shù)，孔隙作為人體內(nèi)營(yíng)養(yǎng)輸送和骨細(xì)胞黏附的空間載體，較大時(shí)易形成骨質(zhì)疏松降低力學(xué)性能，較小時(shí)不利于骨組織長(zhǎng)入誘導(dǎo)新骨形成，所以孔隙率是描述支架性能的重要參數(shù)之一。根據(jù)定義可得孔隙率的計(jì)算公式如下：

(17)

式中，δ——支架孔隙率

Ve——支架內(nèi)孔隙體積

Vt——包容盒體積

其中，支架內(nèi)孔隙體積可通過ROI分析計(jì)算得出。

5.2 梯度擬合度

梯度擬合度是本文對(duì)仿生骨支架評(píng)價(jià)中的新增參數(shù)，用來表示所構(gòu)建的支架模型與實(shí)際骨骼梯度變化規(guī)律的擬合程度。擬合值越高說明對(duì)實(shí)際骨骼模型的建模結(jié)果越貼合實(shí)際情況，梯度擬合度的計(jì)算方法與遺傳算法中適應(yīng)度值的計(jì)算方法相同

(18)

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

Ai——實(shí)際骨骼孔隙分布向量

5.3 連通性

連通性是描述骨支架內(nèi)部微觀孔互相連接導(dǎo)通情況的參數(shù)，仿生骨支架必須具有一定的連通性才可以保證骨組織和血管的長(zhǎng)入、營(yíng)養(yǎng)物質(zhì)的輸送和細(xì)胞代謝物的排出。

本文對(duì)連通性的計(jì)算定義如下

(19)

式中，η——支架連通率

N——橢球體總個(gè)數(shù)

ni——每個(gè)橢球周圍連通的其他橢球個(gè)數(shù)

對(duì)于橢球體之間的連通情況，主要分為三種：a)獨(dú)立橢球體(η=0)；b)兩個(gè)橢球體導(dǎo)通(η=1)，由獨(dú)立橢球體生成的盲孔和兩個(gè)橢球體生成的單通孔，這兩種情況較多存在于實(shí)際骨骼中處于邊緣位置的部分；但在構(gòu)建仿生骨負(fù)模型時(shí)要盡量減少盲孔和單通孔的存在；c)三個(gè)或三個(gè)以上橢球?qū)?η>=2)，該種情況為最佳情況，可保證仿生骨支架內(nèi)部的各橢球之間大概率導(dǎo)通。

5.4 指標(biāo)得分

基于多屬性決策模型，本文根據(jù)以上評(píng)價(jià)指標(biāo)提出對(duì)模型優(yōu)劣的打分制度：

步驟1：分別對(duì)基于隨機(jī)均分分布的仿生骨支架模型、基于梯度函數(shù)分布的仿生骨支架和基于梯度函數(shù)分布但橢球空間位置未調(diào)整的仿生骨支架模型進(jìn)行評(píng)價(jià)指標(biāo)的計(jì)算，得到原始屬性值(見表2)，即孔隙率、梯度擬合性和連通性得到原始決策矩陣E。

表2 仿生骨支架評(píng)價(jià)

矩陣的行代表模型類別，列代表屬性；

步驟2：對(duì)原始屬性值進(jìn)行歸一化處理，得到新的決策矩陣E′

步驟3：根據(jù)成對(duì)比較陣標(biāo)度表，構(gòu)建各屬性的成對(duì)比較陣F，并運(yùn)用SPSS軟件計(jì)算各屬性權(quán)重G

G=[0.196， 0.311， 0.493]

即，孔隙率、梯度擬合度和連通性的屬性權(quán)重依次為0.196，0.311，0.493。

步驟4：由所得的屬性權(quán)重G和歸一化處理后的屬性值E′，計(jì)算模型得分Si

(20)

由表3結(jié)果得：1)基于梯度函數(shù)構(gòu)建的仿生骨支架模型比隨即均勻分布的模型在各指標(biāo)上均具有明顯優(yōu)勢(shì)；2)通過遺傳算法進(jìn)行空間位置調(diào)整后的模型在連通性和孔隙率方面均有一定的進(jìn)步性；

綜上，本文所構(gòu)建的基于梯度函數(shù)和遺傳算法的仿生骨支架模型為最優(yōu)模型。

表3 各仿生骨支架得分

6 結(jié)論

1)本文首先基于對(duì)實(shí)際骨骼進(jìn)行的Micro CT掃描和ROI分析，構(gòu)建了擬合實(shí)際骨骼微觀孔結(jié)構(gòu)分布規(guī)律的梯度函數(shù)。以橢球體作為仿生骨支架微觀結(jié)構(gòu)的單元體模型，圓柱體外形為包容盒，根據(jù)得到的梯度函數(shù)控制橢球體的空間位置和長(zhǎng)軸分布，之后再運(yùn)用遺傳算法對(duì)橢球體的中心點(diǎn)位置進(jìn)行適當(dāng)調(diào)整，得到微觀孔結(jié)構(gòu)負(fù)模型，最后通過布爾減運(yùn)算得到仿生骨支架模型；

2)本文為仿生骨支架建模方法提出了一條新的技術(shù)路線和建模方法。該方法對(duì)所有不同形態(tài)和分布的骨骼均適用，且以實(shí)際骨骼數(shù)據(jù)為依據(jù)進(jìn)行建模，具有一定的科學(xué)依據(jù)和參考價(jià)值；

3)以仿生骨支架的孔隙率、梯度擬合度和連通性作為評(píng)價(jià)指標(biāo)，并基于多屬性決策模型提出打分制度，從而建立出仿生骨支架微觀結(jié)構(gòu)的評(píng)價(jià)體系，并提出了各指標(biāo)與打分制度的計(jì)算方法，以實(shí)現(xiàn)對(duì)所建模型的評(píng)價(jià)和優(yōu)化；

4)微觀孔結(jié)構(gòu)負(fù)模型的生成可由程序控制生成，提高了建模效率和自動(dòng)化程度，為實(shí)現(xiàn)計(jì)算機(jī)輔助組織工程奠定了理論和技術(shù)基礎(chǔ)。