于緯倫 武曉剛,2) 李朝鑫 孫玉琴 張美珍 陳維毅,3)

?(太原理工大學生物醫(yī)學工程學院，太原 030024)

?(太原理工大學體育學院，太原 030024)

引言

骨組織是一種非常智能的材料，它能夠根據日常所受的生理載荷調整其骨量和微觀結構.這種自適應的行為被認為是由骨細胞主導的.骨細胞是骨內所有細胞中對力學刺激最敏感的細胞[1]，且負責將力學刺激信號轉化為成骨細胞和破骨細胞的調節(jié)信號，從而使骨適應環(huán)境變化[2].實驗研究指出負荷引起的流體流動在骨細胞膜上產生剪切應力是骨重建過程中生化反應被激活的主要的刺激[3-5]，且與骨的力傳導機制有關[6].因此，研究骨細胞在骨陷窩--骨小管系統內經歷的流體流動行為的特征是非常重要的.

由于骨細胞嵌入礦化的細胞外基質，直接進行實驗研究具有難度，而體外實驗并不能真實呈現骨細胞在體內所受的力學環(huán)境[7].研究人員試圖通過建立骨骼間質流體流動的數學模型來克服這個問題[8].滲透率是影響間質流體通過骨內孔隙流動的一個重要特性.在利用數值模擬方法來探索外力引起的骨內液體流動時，研究人員通過多孔彈性介質模型計算出骨內的液體流速、壓力梯度、流體剪切力等參數，但這些模型都假定骨的滲透性是各向同性的(骨細胞為球型)[7-9]，不能準確反映骨內具體微觀結構的(如骨細胞--骨小管系統的形狀，方向和密度等)液體流動信息.

近年來，人們對骨細胞的空間特性(包括密度和形態(tài))以及這些特性與疾病和衰老的潛在關系越來越感興趣.體外實驗研究表明，骨細胞幾何形狀影響其應變響應[10].Carter 等發(fā)現股骨的前側，后側，內側和外側的骨細胞的密度，形狀和方向都有較大差別，其可能原因是與負載的局部變化有關[11].骨質疏松癥已經成為世界范圍內的一個嚴重問題，最近有研究表明患有骨質疏松的病人的骨細胞形狀更不規(guī)則[12-13]，骨小管更彎曲[13]，而其骨細胞膜上感受的流體剪應力和流體速度也發(fā)生較大的變化[14].年齡也是影響骨細胞形狀的一個重要因素，年齡增大會導致骨細胞面積減小且形狀偏扁平[15-16].以前的研究一般都認為骨陷窩--骨細胞的方向與骨中膠原纖維的方向一致[9,17-19]，即滲透率主方向與坐標軸方向一致.事實上,二者方向會有不一致的情況[11,20-21]，這會使得骨內液體流動主方向(即滲透率主方向)與膠原纖維的方向成一定的夾角.骨陷窩--骨細胞--骨小管的這種形態(tài)和方向變化會導致骨組織內滲透率的各向異性行為.在進行有限元多孔彈性力學分析時，需要準確量化骨陷窩--小管孔隙的滲透性來捕捉骨的各向異性流體流動行為，才能更真實地預測孔隙流體壓力和流速響應.

因此，為了更精確地闡明骨內液體流動的具體形式，本研究開發(fā)一種基于骨組織微觀結構的多孔彈性力學有限元模型.首先根據骨陷窩的形狀，方向和其他微觀結構生理參數，計算得到骨的三維滲透率.然后，基于多孔彈性力學理論建立骨單元的有限元模型，計算在軸向載荷作用下骨單元內液體的壓力和流速，以期為骨的力傳導與骨的功能適應性機制提供一個更深入的理解.

1 基于骨單元微觀結構計算其滲透率

1.1 骨單元滲透率的計算

圖1(a)是整段骨組織，圖1(b)是圖1(a)中的骨單元.假設圖1(b)中的骨陷窩--小管系統的排列規(guī)則且骨小管分布均勻，則骨單元可以看成是由一個個骨基質包裹著的骨陷窩--骨細胞--小管系統的正方體的周期單元(CPUC)組成(圖1(c)).圖1(d)是骨小管的微觀結構，其中rc是骨小管半徑，ro是骨細胞突觸半徑，a0是骨細胞突觸周圍纖維基質的半徑.

圖1 骨組織分層結構示意圖Fig.1 The hierarchical structure for bone tissue

根據骨陷窩--骨細胞--骨小管的微觀結構得出其滲透率klcp[22-23]

其中，q是rc(0.23μm)和ro(0.1μm)之間的無量綱比(q=rc/ro)[7]，γ 是rc與纖維基質填充的單個骨小管的滲透性(kp)平方根的無量綱比γ=，其中kp=(?/a0)2.377[17,22-23]，a0的半徑是5 nm，?是纖維與纖維的間距(7 nm)[23-24].

A1和B1由以下方程得出

其中，I0，K0，I1和K1是第1 類和第2 類的修正貝塞爾函數.

L表示兩個骨陷窩之間的距離，也是CPUC 的邊長，可由以下公式得出

式中，VL為單位體積，NLac為單位體積內骨陷窩的數量.NLac的范圍在2.6 ～9.0×104(NLac/mm3)[5,11,25-27]，本文選取單位體積的NLac值為3.7×104[28]，得到L的值為30μm.根據文獻報告的平均測量值，每個骨陷窩周圍的骨小管總數N=62[28].通過對幾種不同生物的骨細胞的組織形態(tài)進行考察發(fā)現骨陷窩的形狀類似于橢球形狀[18].用橢球的標準方程來表示骨陷窩

其中，a,b,c分別為骨陷窩長半軸、中半軸和短半軸.則骨陷窩--小管的孔隙率? 為

式中，Lc表示骨小管的平均長度.值得注意的是由于骨細胞與骨陷窩空間內的液體會進行液體交換，且與骨陷窩空間的滲透率相差不大[7,29]，所以我們把整個骨陷窩空間和骨細胞體都算作孔隙空間.ni(i=x,y,z) 是骨小管穿過CPUC 某個面的數目.如圖2(a)所示，為了確定x，y和z方向的nx，ny和nz，我們假設骨小管基于骨陷窩橢圓體投影表面積三維分布，并可以用骨陷窩的投影面積比來測量不同方向骨小管的數量[17]

其 中，nx，ny和nz分別是平行于x，y和z軸穿過CUPC 各面的小管數；Sx，Sy和Sz分別是x，y和z方向上骨陷窩的投影表面積，N是一個骨陷窩周圍的骨小管總數.如圖2(b)所示，當骨陷窩繞x軸偏轉θ 角時，骨陷窩在，和軸上的投影面積和骨小管的三維分布同樣可由式(6)求出.

圖2 (a)骨陷窩沿著x 軸，y 軸和z 軸的投影面積(S x，S y 和S z).(b)骨陷窩繞x 軸旋轉θ 角Fig.2 (a)Projected surface areas of the lacuna-osteocyte in the x,y,and z directions(S x,S y and S z)were used to determine the distribution of canaliculi;(b)The lacuna-osteocyte rotates around x-axis by θ

1.2 骨陷窩形狀

為了更清晰描述骨陷窩形狀，首先規(guī)定了3 個特征值EV1，EV2 和EV3(EV1 是長半軸的平方，EV2 是中半軸的平方和EV3 是短半軸的平方).然后，從顆粒形狀的研究中修改了3 種指標，平衡度(degree of equancy)(EV3:EV1)，伸長度(degree of elongation)(1–EV2:EV1)和扁平度(flatness)(1–EV3:EV2)[11,30]，以定義形狀的差異程度.如表1 所示，取幾組不同形狀的骨陷窩來觀察骨陷窩形狀的差異對其內部液體流動的影響.圖3 展示了這3 種指標對骨陷窩形狀影響的重要性.

圖3 骨陷窩形狀的代表性示例.骨陷窩在不同平面的投影形狀(a) x-y 面，(b)y-z 面，(c)x-z 面Fig.3 Representative cases of lacunae shapes.The project of lacunae shapes are shown schematically in(a)the x-y plane,(b)the y-z plane,and(c)the x-z plane

當3 個特征值相等(EV1=EV2=EV3)時，骨陷窩的形狀即為球形(reference).當前兩個特征值相等，第3 個特征值較小(EV1=EV2>EV3)時，骨陷窩的形狀呈扁平形(Case 1 和Case 2)；第1 個特征值較大，后兩個特征值相等(EV1>EV2=EV3)時，骨陷窩的形狀呈伸長形(Case 3 和Case 4)；3 個特征值都不相同(EV1>EV2>EV3)，形狀偏向取決于伸長度和扁平度的大小(Case 5 和Case 6).

表1 骨陷窩代表性示例的幾何參數和指標Fig.1 Geometrical and degree of representative cases

根據表1 中骨陷窩形狀參數計算得到Case1-Case 6 的滲透率分別為K1，K2，K3，K4，K5和K6

參考模型(reference)模型的滲透率是各向同性的其值計算得：1.05×10?20m2.

1.3 骨陷窩方向

在以前的研究中通常假設骨組織孔隙流體的滲流速度方向與坐標軸方向相同.事實上,由于骨陷窩方向的不同會帶動滲透率方向的偏轉，如果還按原來的方法，則計算結果會產生較大誤差.骨陷窩的軸長一般在幾微米到十幾微米之間[2,7,11,14,31]，當如圖2(b)所示的骨陷窩(a=15μm，b=9μm，c=3.9μm)發(fā)生θ 角的偏轉時，骨組織滲透率可以用全方向的滲透率張量表示

當θ 角分別為0?，30?，45?，60?和90?時,對應于o-xyz坐標下滲透率分別為

2 骨單元控制方程及有限元模型的建立

2.1 控制方程

本文將骨單元看成由CUPC 單元組成的固--液耦合多孔彈性材料，其本構方程可由以下方程來描述[33-35]

其中，σ 是整個多孔介質對應的應力張量，ε 為應變張量.M是多孔介質脫水的4 階彈性矩陣，α 和M分別是Biot 有效應力系數張量和模量.p為多孔介質孔隙內液體壓力，ξ 為單位體積內液體含量的改變量.tr()表示矩陣的跡算子.

在不考慮體積力的情況下的平衡方程[8]

液體質量守恒方程

達西定律可寫為[36]

式中，V是流速矢量，k和μ分別是滲透張量和液體的黏度.將式(8)代入式(10)，式(9)和式(12)代入式(11)得到骨單元的控制方程

2.2 有限元模型的建立

在這項工作中，使用COMSOL Multiphysics 軟件探討軸向壓縮載荷作用下骨單元內的流體與固體相互作用的多孔彈性行為.如圖4(a)，數值模擬過程中骨單元模型被定義為一個由CPUC 組成的橫觀各向同性的中空的柱體(忽略了哈弗氏管)，其部分材料和幾何參數設置如表2 所示.

表2 幾何和材料常數Table 2 Geometrical and material constants used in osteon model

考慮到載荷頻率較低，哈弗氏管像蓄水池一樣起到儲蓄的作用，能夠維持流體正常的流進流出，其內的壓力可設置為0，用作參考壓力.哈弗氏管的孔隙尺寸遠大于骨陷窩--小管孔隙，骨單元壁的孔隙壓力可以由此釋放.哈弗氏管表面的應力可以忽略不計，這也導致在哈弗氏管處應力為0 的邊界條件[37].骨單元外壁設置剛性約束(位移為0)，表示骨單元受周圍骨間質約束不能移動[33].如圖4(a)和方程(14)，軸向對稱的循環(huán)載荷代表縱向壓縮的生理運動，其諧波幅值(w)和頻率(f)分別為0.5μm 和1 Hz.在一個周期內的最大軸向應變(ε)發(fā)生在0.5 s 時，其幅度為0.001[33].

圖4 骨單元模型的建立(a)和網格劃分(b)Fig.4 Establishment of osteon model(a)and mesh generation(b)

如圖4(b)所示，借助掃掠的網格剖分方法來生成精確有效的有限元網格，網格包含20 880 個單元.最小單元質量為0.601 5；平均單元質量為0.847 4.

3 結果

3.1 骨陷窩形狀的影響

根據本文的加載形式，當t=0.25 s 時，應變率達到最大值0.001 s?1，骨單元內孔隙液體壓力和流速的變化規(guī)律有最大響應[12,33].如圖5 繪制了t=0.25 s時考慮不同形狀骨陷窩CUPC 的骨單元孔隙壓力及流速的分布圖.在軸向的對稱載荷下，參考模型呈現的是對稱分布與之前的研究結果吻合[12,33,36].由于內部骨陷窩形狀的不同，Case 1 ～Case 6 的液體的壓力和流速(靠近哈弗氏管部位較為明顯)的大小與分布規(guī)律均與參考模型有較大差別，且最大的孔隙壓力值均大于參考模型.這說明在骨陷窩體積相同的情況下，形狀的差異程度會影響壓力和流速的大小和分布規(guī)律.Case 1，Case 2 和Case 5 的骨陷窩形狀都是扁平形的，Case 3,Case 4 和Case 6 的骨陷窩形狀都是伸長形的.Case 1，Case 2 和Case 5 的壓力和流速在分布上與參考模型有些類似，沿著y?z面和x?z面的分布差異不大.Case 3，Case 4 和Case 6 的模型的壓力和流速在y?z面和x?z面的分布差異較大.

由圖6 可以看出，在t=0.25 s 時，Case 1--參考模型的壓力和流速分別為3.21 × 105Pa，4.29 ×105Pa，3.9×105Pa，3.43×105Pa，3.57×105Pa，3.84× 105Pa，2.51 × 105Pa 和6.34 × 10?8m/s 和6.57 ×10?8m/s,6.5 × 10?8m/s，6.81 × 10?8m/s，6.74 ×10?8m/s，6.4 × 10?8m/s，6.72 × 10?8m/s 和5.85 ×10?8m/s.在不同Case 模型中，Case 2 和Case 4 模型沿y軸的壓力差距最大，增大了86%，Case 5 和參考模型沿y軸流速差距最大，增大了18%.在單個Case 模型上，Case 1，Case 2，Case 5 和參考模型沿x軸和y軸的最大壓力和流速差別不大，而Case 3，Case 4 和Case 6 沿x軸和y軸的壓力有明顯區(qū)別，分別增大了62%，47%和58%，而流速差別不大.

圖5 骨陷窩形狀影響下的骨單元孔隙壓力(p)及流速(v)的分布圖Fig.5 The pore pressure(p)and flow velocity(v)distribution of osteon due to the influence of lacunae shape

圖6 骨陷窩形狀對骨單元孔隙壓力和流速沿著不同方向的影響Fig.6 Effect of the lacunae shape on the pore pressure and flow velocity along the different directions of osteon

3.2 骨陷窩方向的影響

圖7 繪制了t=0.25 s 時考慮具有不同方向骨陷窩CUPC 的骨單元孔隙壓力及流速的分布圖.在我們所研究參數范圍內，0?～90?模型的液體的壓力和流速的大小與分布規(guī)律均有較大差別.隨著骨陷窩旋轉角度的增加其沿著y?z面和x?z面的壓力分布差距越大.30?，45?，60?模型的流速在哈弗氏管處出現流速局部增大的情況.

圖7 骨陷窩方向影響下的骨單元孔隙壓力(p)及流速(v)的分布圖Fig.7 The pore pressure(p)and flow velocity(v)distribution of osteon due to the influence of lacunae direction

由圖8 可以看出，骨陷窩在0?～90?的旋轉過程中導致骨單元的模型的壓力值逐漸減小.在y軸方向上30?，45?和60?的模型流速出現明顯的不均勻分布.滲透率的主值方向與坐標軸方向出現偏轉，即液體流動主方向(,和方向)與骨膠原纖維方向(x,y和z) 成θ 角時會導致骨組織內壓力和流速大小和分布規(guī)律的變化.根據所選擇的骨陷窩(a=15μm，b=9μm，c=3.9μm)及其偏轉的角度(繞x軸偏轉0?～90?)，不同模型之間的最大的壓力和流速有著較大差異，其中0?沿y軸的壓力比90?增大了125%；30?沿y軸的流速比90?增大了56%.同一模型中由于偏轉角度的不同，沿不同坐標軸方向最大壓力和流速也存在這較大差異，其中60?模型x軸上的最大壓力比y軸上增大了58%;30?模型沿y軸上的流速比x軸上增大了50%.可見，骨陷窩角度的不同，可明顯使得局部的壓力和流速發(fā)生改變.

圖8 骨陷窩方向對骨單元孔隙壓力和流速沿著不同方向的影響Fig.8 Effects of the lacunae direction on the pore pressure and flow velocity along the different directions of osteon

圖8 骨陷窩方向對骨單元孔隙壓力和流速沿著不同方向的影響(續(xù))Fig.8 Effects of the lacunae direction on the pore pressure and flow velocity along the different directions of osteon(continued)

4 討論

骨組織的滲透率是骨陷窩--小管流體流動空間的滲透率，不能脫離流體流動通道獨立存在，也不可簡單矢量合成.本研究首先通過估計不同形狀和方向的骨陷窩中骨小管的數量，其次利用這些骨小管的三維分布以及骨組織其他微結構參數(骨陷窩密度，骨小管半徑和骨細胞突觸的半徑等)來計算骨組織的滲透率和孔隙率等參數，最后根據計算所得的參數建立骨單元的多孔彈性有限元模型.結果表明，由于骨陷窩形狀和方向的不同對骨單元內液體壓力和流速的大小和分布規(guī)律有著明顯的影響.

孔隙流體壓力是骨陷窩--小管系統中的一種重要負荷誘導現象，它在很大程度上干預骨細胞生長，分化和化學物質運輸[7,38].骨陷窩的形狀和方向發(fā)生變化時，孔隙壓力行為有著明顯的改變.具體來說，對于形狀改變的模型，Case 3，Case 4 和Case 6 在x和y方向上的壓力差距最大，而Case1 Case2，Case5 和參考模型差距不大.這種情況和骨陷窩形狀的改變引起的滲透率的各向異性有關，可以發(fā)現，Case 3，Case 4 和Case 6 模型在x和y方向的滲透率差了一個數量級，而Case1 Case2，Case5 和參考模型在x和y方向的滲透率在一個數量級且差距不大.而在z方向上所有的模型的壓力都沒有明顯的變化，這是因為負載誘導產生流體壓力的機制使得液體通過骨陷窩--小管流入哈弗氏管，使孔隙壓力得到釋放，即骨單元內液體的主要流動是在骨單元外壁和哈弗氏管之間，而在z方向上幾乎沒有流動.

對于角度不同的模型，設定的是骨陷窩繞x軸旋轉θ 角.實際上，兩個坐標系中x軸與x軸是重合的，因此骨單元在x方向上的滲透率是相同的.如圖8(a)和圖8(c)可以看到，沿x軸上的壓力變化遠小于沿y軸的變化.在y方向的滲透率變化大約在一個數量級左右.如圖8(c)，在y軸方向上0?～90?的模型的滲透率是逐漸增大的，其壓力是逐漸減小的.在30?，45?和60?時，骨中膠原纖維的方向與骨陷窩的方向不在一個方向，會導致y-z方向的滲流行為，引起流速局部增大的情況.總體來看，壓力大小的變化與滲透率有很大相關性，隨著滲透率的減小，壓力有增大的趨勢.

雖然骨陷窩方向的改變，使得30?，45?和60?模型局部的流速發(fā)生較大的變化，但就整體來看骨陷窩形狀和方向的改變對流速的影響都不超過一個數量級.這種現象可以解釋為，滲透性由1.0×10?20m2降低到1.0×10?21m2的過程中抑制了流體流動，但孔隙壓力和壓力梯度相應的增加維持了流體速度的大小，從而維持了骨陷窩--小管內正常的對流運輸，因此骨內達西流速變化不大.

在典型的人類日常生理活動的載荷(1.0×103με)和頻率(1 Hz)下[39]，分別計算出了各向異性骨陷窩--小管孔隙的最大孔隙壓力和流速響應，壓力范圍從2.6×105Pa 到4.8×105Pa，流速范圍從4.7×10?8m/s 到9.8×10?8m/s.壓力的大小足以驅動流體跨皮質流動(大于6.0×103Pa)[7].流速大小也足以直接刺激骨細胞(大于2.0×10?8m/s)[36].每一個骨單元(哈弗氏系統)都可以充當其自身的排水系統，并通過哈弗氏管釋放因骨組織變形產生的流體壓力.盡管單個骨單元在整體骨組織中受到骨外膜到骨內膜產生的跨皮質壓力梯度的影響，但一個哈弗氏系統產生的局部的壓力梯度完全可以控制本身骨陷窩--骨小管流體的流動.因此，本文未考慮相鄰骨單元之間的耦合程度，選用不透水的粘合線計算流體壓力和速度大小是符合生理值要求的.

我們研究的一個局限性是骨小管被理想化為直管.實際上，骨細胞突觸是通過彎曲的小管從骨陷窩延伸到CUPC 表面[40].這項研究沒有考慮到骨小管彎曲的影響，因為精確計算彎曲的骨小管非常困難.這個限制可以在模型的下一步改進中解決.另一個局限性是本文把骨單元看成是由一個個相同的CUPC 組成，但是實際上每個CUPC 內的骨陷窩形狀和方向會有區(qū)別，這將引起局部的液體流動的變化.雖然理論上模型需要確定每個CUPC 內骨陷窩的形狀和方向，但據實驗觀察了解到，在骨組織的一定區(qū)域內骨陷窩的形狀是類似的[11]，且骨中膠原纖維的方向與骨陷窩的方向有一定的相關性[18-19].這樣的區(qū)域大到包含一個或者幾個骨單元大小[11]，因此本研究的結果可用于建立骨單元的各向異性多孔彈性有限元模型.因此，只要確定所分析的特定區(qū)域內骨陷窩形狀和方向，就可以把本研究中計算的各向異性滲透率應用到有限元模型中，以求解載荷誘導的流體壓力和速度場.

5 結論

本研究提出了一種根據骨陷窩形狀和方向估計骨單元內骨陷窩--骨小管孔隙各向異性滲透性的方法，并通過確定的骨陷窩--小管滲透率來描述骨單元內液體各向異性的流動.結果表明，骨陷窩的形狀和方向以及骨小管的三維分布情況是決定骨陷窩--小管孔隙內液體流動的各向異性程度的重要參數.形狀不同的模型最大壓力和流速的差異分別可達86%和18%；方向不同的模型最大壓力和流速的差異分別可達125%和56%.本文的研究對進一步理解骨力傳導機制和一些像骨質疏松之類的骨質疾病具有重要參考價值.