正常和早期膝骨關節炎的軟骨生物力學研究1)

林偉健 李俊言 陳瑱賢，2) 靳忠民，，

* (西安交通大學機械制造系統工程國家重點實驗室，西安 710054)

? (西南交通大學機械工程學院，成都 610031)

** (長安大學工程機械學院，西安 710064)

引言

軟骨是人體關節中非常重要的一種承重組織[1-2].它不僅可以減少關節的運動摩擦，而且能通過增大關節的接觸面積來分散載荷、減少應力集中.軟骨由纖維、蛋白聚糖、軟骨細胞和大量的水分組成[1].根據內部纖維結構，軟骨通常分為3 層[1]:表層(占比10%～ 20%)、中層(占比40%～ 60%)和深層(占比約30%).在表層，纖維方向平行于軟骨表面;在中層，纖維分布沒有特定方向;在深層，纖維方向垂直于軟骨表面.由于軟骨復雜的組成成分和結構，其具有特殊的力學性質[1].例如，軟骨是一種固-液雙相材質，其承受瞬間載荷時，軟骨中的液相會承擔大部分載荷，從而減少軟骨固相基質的變形;軟骨內部纖維具有抗拉伸能力，軟骨表層中平行于軟骨表面的纖維的抗拉伸能力可以有效地抵抗軟骨的側向擴張[3]，能有效減小軟骨固相基質的變形，避免軟骨產生過大應變而導致軟骨受損和退變.

膝骨關節炎是老年人中常見的慢性膝關節疾病[4]，其惡化伴隨著軟骨的退變、磨損和最終缺失[4-5].軟骨的退變和磨損通常認為是從表層發展到深層，該過程中軟骨的組成成分、結構和力學性質會發生改變[5]，如蛋白聚糖降解、纖維退變、水含量和滲透率增大以及纖維方向變得無序等等.由于這些改變，軟骨的液相承載能力和纖維的抗拉伸能力等等會大大下降，從而進一步影響了關節的生物力學[4].

早期膝骨關節炎會對軟骨生物力學造成影響[4]，有時也會造成關節疼痛[6].研究表明[4，7]，早期內側膝骨關節炎會引起內側軟骨流體壓力下降以及外側軟骨接觸應力增大等力學變化.此外，與內側早期膝骨關節炎的患者相比，內外側膝骨關節炎的患者通常具有更嚴重滑膜炎[8];而且在日常活動中，內外側早期膝骨關節炎的患者感覺更加疼痛[6].但目前導致正常、內側和內外側早期膝骨關節炎3 種情況下的臨床差異的機理尚未明確.因此，研究正常、內側與和內外側早期膝骨關節炎的軟骨生物力學差異，能從生物力學角度為理解早期膝骨關節炎提供幫助，以便更好地計劃預防和治療方案.

目前通過實驗測定關節軟骨生物力學參數十分困難，因此人們通常應用有限元方法進行關節軟骨生物力學的研究[9-10].由于建模簡單且計算時間短，單相各向同性的軟骨有限元模型使用最廣泛[11-14]，但其提供的計算精度和力學信息有限.近期，張吉超等[15]的綜述報道，國內近5 年膝關節有限元研究中的軟骨模型仍采用單相各向同性模型.然而軟骨的力學行為受其固相和液相的相互作用、纖維網絡和深度相關屬性的影響[10，16-17].因此，為了更精確地反映軟骨內部復雜的纖維結構和組成成份，需要構建更加先進的模型.國外文獻報道中，固-液雙相有限元模型[1]考慮了固-液雙相成份;雙相纖維增強有限元模型[10，18]考慮了軟骨內部的纖維結構;非均質雙相纖維增強有限元模型[17，19]考慮了軟骨深度相關和應變相關的屬性.但截至目前，國內尚未有膝關節有限元研究應用雙相纖維增強的軟骨模型進行力學分析.因此，構建基于雙相纖維增強軟骨模型的膝關節有限元模型對推動國內生物力學的發展有著重要的意義.

前期的膝骨關節炎有限元研究中，國內張震等[12]通過輕度膝骨關節炎的MRI 和CT 數據構建了關節炎早期軟組織的形貌特征，分析了軟骨和半月板的應力分布特征，但是其模型沒有考慮軟骨和半月板的固液雙相屬性，也沒有考慮退變軟骨的屬性變化.國外Dabiri 和Li[4]構建了早期退變軟骨模型來研究膝關節力學的變化，雖然研究考慮了軟骨和半月板的固液雙相纖維增強屬性，但是該模型僅施加的0.1mm 的壓縮載荷，沒有模擬體內膝關節真實的生理載荷.Mononen 等[7]通過構建內側股骨髁關節炎的固液雙相纖維增強軟骨模型，在生理載荷下研究了纖維分布對軟骨承載能力影響，但其未考慮膝關節中真實的固液雙相接觸條件，即接觸區域的流體流動取決于接觸區域界面的流體壓力差，而非接觸區域允許流體自由流動[17，20].綜上所述，目前尚未有研究應用固液雙相纖維增強軟骨模型及真實的雙相接觸條件分析早期膝骨關節炎的軟骨生物力學.

本文構建了基于固液雙相纖維增強軟骨模型的膝關節有限元模型，應用行走姿態下生理載荷和真實的固液雙相接觸條件，對比研究了膝關節正常、內側與內外側早期膝骨關節炎3 種情況下的軟骨生物力學差異，以期為臨床上更好地理解、預防和治療早期膝骨關節炎提供理論支撐.

1 材料和方法

基于生物力學有限元軟件FEBio(版本2.9，www.febio.org)，本文構建了正常膝關節模型、內側早期膝骨關節炎模型和內外側早期膝骨關節炎模型.

1.1 正常膝關節模型

本文采用了來自“open knee project”項目[21]中一名70 歲女性捐獻者(體重77.1 kg)的右膝關節模型，該模型已在之前的研究中采用并驗證[22].通過1.0 Tesla 肢體掃描儀采集膝關節完全伸展位置的磁共振圖像，從圖像中重建出包括股骨、脛骨、軟骨、半月板和韌帶在內的膝關節三維模型(如圖1所示).在Hypermesh (版本 14.0;Altair Engineering，Inc.，Troy，MI，USA)對模型進行網格劃分，骨頭使用四邊形殼單元，其他軟組織使用六面體單元.通過網格收斂性分析，六面體單元數量為95 000 左右時，相對于單元數量為190 000 左右時，接觸應力和流體壓力變化小于5%[10，22]，最終確定六面體單元的尺寸范圍在0.5～ 1.0 mm 之間.

由于骨頭的剛度遠大于其他軟組織，因此將股骨和脛骨設為剛體[12，23].軟骨和半月板設為雙相纖維增強材料[16-17]，由非纖維材料Neo-Hookean、纖維材料和液相組成.參照先前開發的模型[17]，軟骨設為深度相關的材料屬性，包括彈性模量、泊松比、固相比和滲透率[24-25](如表1 所示).軟骨纖維方向在深度方向是變化的，如圖1 所示，纖維在表層平行于軟骨表面，而且呈放射線分布[7];纖維在中間層在3 個方向上均勻分布;纖維在深層垂直于軟骨表面.半月板內部纖維呈周向和徑向分布[22].

表1 本文中軟骨、退變軟骨和半月板的材料屬性(E:彈性模量;v:泊松比;φs:固相比;k:滲透率)Table 1 Material properties of cartilage，degenerated cartilage and meniscus in this study (E:Young's modulus;ν:Poisson's ratio;φs:solid volume fraction;k:permeability)

圖1 (a)膝關節的幾何模型;(b)軟骨、半月板纖維分布模式;(c)軟骨表層纖維的自然放射線分布模式Fig.1 (a) The geometry of the knee joint ;(b) fibrillar distributions pattern of cartilage and meniscus ;(c) the natural split-line patterns in surface zone

根據先前研究[17，22]，材料Neo-Hookean 的應變能函數定義為

式中，J為相對體積(當前體積/初始體積);I1為右柯西-格林變形張量的第一不變量;μ和λ為拉梅常數，通過泊松比v和彈性模量E可以計算獲得

纖維的應變能函數定義為[22，27]

式中，ξ為纖維模量數值的1/4，α為指數參數的系數，β為指數參數的冪，In為纖維伸長量的平方.

當α→0 時，式(3)會產生一個冪法則

纖維只能拉伸，壓縮時不起作用.因此本文中選擇α=0 和β=2，來實現應力從壓縮到拉伸不連續性轉變[17].

軟骨應變相關的滲透率(k)呈指數遞減[28]

式中，k0為初始滲透率，φs為固相比，M為應變相關的指數常數，αe為冪法則指數.根據先前實驗[28]，本研究選擇M=4.638 和αe=0.084 8 來實現應變相關的滲透率.

韌帶為纖維增強材料，由非纖維材料Mooney-Rivlin 和纖維材料組成，包括前交叉韌帶(anterior cruciate ligament，ACL)、后交叉韌帶(posterior cruciate ligament，PCL)、內側副韌帶(medial collateral ligament，MCL)、外側副韌帶(lateral collateral ligament，LCL).材料參數[21]如表2 所示.

表2 本文中韌帶的材料屬性[21]Table 2 Material properties of ligaments in this study[21]

韌帶材料的應變能函數定義為[21]

式中，C1和C2為Mooney-Rivlin 常數，C3至C6為纖維函數常數;和分別為偏右柯西-格林變形張量的第一和第二不變量;K為體積彈性模量;Ei(·)為指數積分函數; λ～ 為沿纖維方向拉伸的偏量; λm為拉直纖維的伸展量.

軟骨的底面和韌帶的末端都固定在骨頭上[12].半月板兩端通過線性彈簧固定在脛骨上，各端彈簧總剛度為2 kN/mm[18].軟骨對軟骨、軟骨對半月板之間的接觸設置為雙相接觸，并且接觸界面為無摩擦[22].在雙相接觸表面，在接觸區域流體的流動取決與接觸界面的流體壓力差，在無接觸區域流體可以自由流動[17，20].該雙相接觸條件可以通過FEBio自帶的“biphasic contact”接觸條件來實現.

1.2 膝骨關節炎模型

通過更改自然膝關節模型的軟骨表層和中間層的材料屬性和纖維分布來模擬早期關節炎.根據研究報道[4，7，26]，相比自然膝關節模型，早期膝骨關節炎模型中的退變軟骨的表層和中間層的非纖維固相基質的彈性模量和纖維的彈性模量分別減少65%和75%;滲透率增加70%;固相比減少5%;纖維方向隨機分布(如表1 和圖2 所示)，其他參數保持不變.至此，如圖2 所示，在本研究中定義正常膝關節為模型1，內側早期膝骨關節炎為模型2，內外側早期膝骨關節炎為模型3.

圖2 不同的膝骨關節炎模型以及它們的纖維分布模式，綠色代表正常組織及正常的纖維方向，紅色代表退變組織及隨機分布的纖維方向Fig.2 Different models of knee OA and their fibrillar distributions pattern.The green denoted the normal tissue and fibrillar direction，and the red denoted the degenerated tissue and randomly distributed fibrillar direction

1.3 載荷和邊界條件

首先，為了驗證膝關節模型的合理性，本研究對模型一的股骨施加1 kN 的垂直載荷，將求解的股骨垂直位移和軟骨接觸應力與文獻結果[29-32]進行對比.脛骨完全固定，股骨屈曲運動被約束.參考膝關節假體磨損試驗的載荷施加ISO 標準(ISO 14243—3:2014)，將1 kN 壓縮載荷施加在股骨的參考點上，該參考點位于關節中心偏內側5 mm 處[22].

在模型驗證之后，在1 s 內采用步態周期內最大載荷時刻(3 倍體重[33]，13%步態周期)和最大屈曲角度時刻(屈曲58°，72%步態周期)的力和運動作為邊界條件，分別對模型1～ 3 進行求解計算，輸出了軟骨的流體壓力、應變和應力等力學信息.

2 結果

如表3 所示，通過與先前實驗測得的股骨垂直位移和關節接觸應力對比，對雙相膝關節建模方法的合理性進行了驗證.在軸向載荷1 kN 下，本文中股骨垂直位移為1.11 mm，文獻[22，29-30]報道范圍在0.70 mm 至1.25 mm 之間;本文中關節峰值接觸應力為3.4 MPa，文獻[30-32] 報道范圍為2.68～5.85 MPa 之間.本文結果均在實驗結果和仿真結果報道范圍之內，雙相膝關節模型的合理性得到了驗證.在許多膝關節有限元研究中[7，11，22]，股骨垂直位移和關節接觸應力分別經常被采用以驗證模型的合理性，它們分別體現了膝關節整體和局部信息.股骨的垂直位移是膝關節軟組織的整體剛度的體現[22，29]，對應的是軟組織在承受載荷時的總變形量，是膝關節整體測得的數據;接觸應力是脛骨軟骨和股骨軟骨兩個相對表面之間載荷傳遞和接觸面積的體現[22，30]，是膝關節局部測得的數據.因此，膝關節整體和局部信息結合可以有效地驗證模型的合理性.

表3 本文中模型1 在1 kN 載荷下的股骨垂直位移和軟骨接觸應力和先前實驗/仿真結果的對比Table 3 Comparison of the femoral vertical displacement and the contact pressure under 1 kN load between the predicted results of Model 1 and published experimental/computational studies

在最大載荷時刻，模型1～ 3 的軟骨流體壓力、固相等效應力和壓縮應變如圖3 所示.與模型1 相比，模型2 的內側軟骨的流體壓力和固相等效應力都明顯減少，分別減少18.1%和13.3%，而外側軟骨的流體壓力和固相等效應力變化小于2%;模型3 的內側軟骨的流體壓力和固相等效應力都明顯減少，分別減少18.8%和13.3%，外側軟骨的流體壓力和固相等效應力也都有明顯減少，分別減少16.7%和14.5%.模型2 的內側軟骨的壓縮應變增大5.0%;模型3 的內外側軟骨的壓縮應變都有明顯增大，分別增大5.0%和13.2%.與模型1 相比，模型2 和模型3 的退變軟骨出現應力與應變相反的變化趨勢是由于相同載荷下軟骨的彈性模量下降所造成的.

圖3 模型1 至模型3 在最大載荷時刻下的預測結果:流體壓力、固相等效應力和壓縮應變Fig.3 The predicted results of Model 1 to 3 under the maximum load situation:fluid pressure;effective solid stress;compressive strain.

在最大載荷時刻，模型1～ 模型3 的軟骨第一主應變和最大剪應變如圖4 所示.從圖中可以得知，與模型1 相比，模型2 的內側軟骨的第一主應變和最大剪應變分別增大了19.8%和20.0%，外側軟骨應變沒有太大變化;模型3 的內外側軟骨的第一主應變和最大剪應變均明顯增大了，增大范圍在19.0%～ 46.0%之間.

圖4 模型1～ 模型3 在最大載荷時刻下的預測結果Fig.4 The predicted results of Model 1 to 3 under the maximum load situation

在最大屈曲角度時刻，模型1～ 模型3 的軟骨流體壓力、固相等效應力和壓縮應變如圖5 所示，結果與最大載荷時刻的結果趨勢一致，但數值相對較小.從圖中可以得知，與模型1 相比，模型2 的內側軟骨的流體壓力和固相等效應力都明顯減少，分別減少25.9%和29.1%，而外側軟骨的流體壓力和固相等效應力變化小于5.5%;模型3 的內側軟骨的流體壓力和固相等效應力都明顯減少，分別減少26.5%和32.1%，外側軟骨的流體壓力和固相等效應力也都有明顯減少，分別減少13.0%和28.3%%.模型2 的內側軟骨的壓縮應變增大26.9%;模型3 的內側和外側軟骨的壓縮應變都有明顯增大，分別增大26.9%和20.8%.

圖5 模型1～ 模型3 在屈曲角度時刻下的預測結果:流體壓力、固相等效應力和壓縮應變Fig.5 The predicted results of Model 1 to 3 under the maximum flexion angle situation:fluid pressure;effective solid stress;compressive strain

在最大屈曲角度時刻，模型1 至模型3 的軟骨第一主應變和最大剪應變如圖6 所示.從圖中可以得知，與模型1 相比，模型2 的內側軟骨的第一主應變和最大剪應變分別增大76.3%和59.6%，外側軟骨應變沒有太大變化;當模型3 的內側軟骨的第一主應變和最大剪應變分別增大73.2%和56.7%，外側軟骨應變沒有太大變化.模型2 和模型3 的第一主應變和最大剪應變結果區別不大，可能是因為最大屈曲角度時刻的載荷較小，因此不同膝骨關節炎模型的第一主應變和最大主應變區別不明顯.

圖6 模型1～ 模型3 在最大屈曲角度時刻下的預測結果Fig.6 The predicted results of Model 1 to 3 under the maximum flexion angle situation

3 討論

本文基于固液雙相纖維增強的軟骨有限元建模方法，構建了正常膝關節雙相有限元模型，再根據正常模型建立了內側和內外側早期膝骨關節炎模型.之后在步態周期內的瞬時承載情況下和真實的雙相接觸條件下，對比研究了正常膝關節、內側和內外側早期膝骨關節炎的軟骨生物力學差異.

早期膝骨關節炎中退變軟骨屬性變化對生物力學特征改變的映射關系是復雜的.總的來說，退變軟骨的非纖維固相基質、纖維、滲透率等的變化導致了軟骨流體壓力下降、固相等效應力下降以及應變增大.首先，滲透率的增加導致退變軟骨的流體壓力下降.滲透率的增加加快了流體流動的速度，從而導致流體壓力下降.流體壓力的下降對軟骨的負載機制造成影響，即導致軟骨液相承載能力下降.其次，非纖維固相基質彈性模量的下降導致非纖維固相基質的強度下降，從而導致退變軟骨的固相等效應力下降，與先前研究的結果一致[34].軟骨非纖維固相基質變軟會導致固相承載能力下降.最后，軟骨表層和中間層纖維的纖維彈性模量下降以及纖維方向變得無序加劇了軟骨側向擴張變形，進而削弱了軟骨的抗拉伸能力.此外，液相承載能力下降、固相承載能力下降以及抗拉伸能力下降都會導致應變增大，而應變增大又導致應變相關的滲透率繼續增加，從而又進一步削弱了軟骨的液相承載能力.因此，膝骨關節炎的退變軟骨屬性變化對軟骨生物力學產生的影響是復雜的、相互關聯的，最終導致軟骨的承載能力和抗拉伸能力下降.

膝關節是一個復雜的系統，退變軟骨的承載能力下降和抗拉伸能力下降會加劇軟骨進一步退變的風險[4]，進而可能加劇膝關節疼痛.軟骨的承載能力下降和抗拉伸能力下降可能導致軟骨變形過大、軟骨細胞死亡和纖維退變[19，35]，甚至可能造成軟骨細胞外基質受損從而釋放出軟骨細胞外基質的內部成分[36-37].這些成分通過與滑膜細胞上的受體結合，從而導致滑膜炎進一步發展[36-37].有研究[38]表明滑膜炎反應與膝骨關節疼痛強烈相關.因此，內外側膝骨退變軟骨的內外側承載能力和抗拉伸能力同時下降以及應變增大可能是疼痛感更明顯的潛在因素.

將雙相纖維增強軟骨模型應用到膝關節雙相有限元模型中是具有挑戰性的.首先，軟骨組織的組成成分和結構十分復雜[1]，因此其力學性能是高度非線性的，例如深度相關性和應變相關性.其成份、結構以及力學參數的測量常常面臨著巨大的挑戰，需要進行多種實驗才能測得[1，27].其次，構建一個層狀雙相纖維增強軟骨有限元模型是具有難度的[17].由于真實軟骨的幾何形狀十分不規則，且厚度不均勻.因此，在有限元網格劃分時，將軟骨劃分成多層網格模型是比較困難的.同時，為了正確設置軟骨內部纖維方向，還需要調整網格分布方式和網格節點順序.再次，膝關節雙相有限元模型收斂難度大于單相模型.雙相有限元模型受載后容易導致網格變形過大，常常造成收斂困難.最后，現有文獻中的膝關節雙相有限元模型[4，12，34]很少考慮和應用了真實的雙相接觸條件[17，20].大部分膝關節雙相有限元研究[4，12，34]都是使用Abaqus 軟件進行分析，而在Abaqus 中需要額外的子程序對接觸面間的液體流動進行修正，才能實現真實的雙相接觸條件;而在FEBio 軟件中，通過在軟件界面設置雙相接觸就可實現真實的雙相接觸條件[20].總而言之，目前膝關節雙相有限元模型的應用仍存在一些挑戰，但其優勢也是十分突出的，例如能比單相模型提供更加深入的力學及摩擦學分析.因此膝關節雙相有限元模型的應用仍有巨大的潛力和很好的發展前景.

本文仍存在一些局限性.首先，當前模型沒有考慮各向異性的滲透率[4，17].據文獻[4]報道，平行于纖維方向的滲透率大于垂直于纖維方向的滲透率，因此，軟骨的滲透率應該是各向異性的.文獻[39]研究表明，各向異性的滲透率對軟骨時間相關行為的影響可以忽略不計，因為在短時間內，軟骨內部的流體還來不及流出.本文主要研究和討論瞬時承載狀態，因此是否采用各向異性的滲透率對本文結果影響不大.此外，目前能精確測得軟骨各向異性的滲透率的裝置設備還有待開發，將各向異性的滲透率應用到膝蓋雙相模型中需要后續工作的不斷補充與完善.其次，當前雙相纖維增強軟骨模型沒有考慮內部的細胞分布[9].軟骨細胞是軟骨中唯一存在的細胞，其生理狀態受應力應變的調控[9，40]，過大的應力應變會導致軟骨細胞的死亡，從而可能進一步導致軟骨退變.因此，通過軟骨的多尺度模型來研究軟骨內部細胞級別的力學信息，能為認識、預防和減緩軟骨退化、關節炎進展提供更有力的幫助.目前已有研究[9]采用非完全耦合的方法實現軟骨的多尺度建模，即將模型分為總模型(軟骨模型)和局部模型(細胞模型).總模型測得的局部應力和局部位移張量轉移到局部模型上，即可以預測軟骨細胞的力學信息.在未來，完全耦合的雙相多尺度軟骨模型還可以耦合到骨肌系統模型中[41]，從而實現個性化邊界條件以及精準生物力學的預測.該方法還有待后續開發和應用.盡管本文建立的模型仍存在局限性，但是現有模型也能為研究內側與內外側膝骨關節炎的軟骨生物力學差異提供一定幫助，從而為臨床上更好的理解、避免和減緩膝骨關節炎發展提供幫助.

4 結論

本文在國內首次構建了基于固液雙相纖維增強軟骨模型的膝關節有限元模型，并采用更真實的在體生理載荷和雙相接觸條件，對比研究了膝關節正常、內側與內外側早期膝骨關節炎的軟骨生物力學差異.與正常膝關節相比，早期的膝骨關節炎導致了關節軟骨承載能力的降低.相比內側早期膝骨關節炎，內外側早期膝骨關節炎導致更大的內外側軟骨流體壓力減少和應變增大.基于固-液雙相纖維增強軟骨模型的膝關節有限元建模方法合理可行，不僅可以用于體內膝骨關節炎模型的仿真，還可以推廣應用于髖、踝和脊柱等其他關節生物力學的研究.本研究為更好的理解早期膝骨關節炎造成的關節疼痛，提供了有力的生物力學依據.