為了探究特發性脊柱側凸術后融合器的植入對椎骨細觀生物力學性能與骨組織成骨性的影響，本文建立了術后植入融合器的宏觀尺度有限元模型，利用圣維南原理子模型方法建立骨單元細觀尺度模型。基于上述模型，本文模擬了人體生理工況，對比相同邊界條件下宏觀皮質骨與細觀骨單元生物力學性能差異，分析細觀尺度下融合器的植入對骨組織生長的影響。結果顯示，腰椎細觀結構的應力較宏觀有所增加，此病例中細觀應力為宏觀應力的2.606～5.958倍，融合器上部骨單元應力較下部大；上部椎體端面平均應力按大小排序分別為右側、左側、后部、前部；下部椎體應力排序為左側、后部、右側、前部；旋轉為骨單元應力值最大工況。據此推測，融合器上端面骨組織成骨性較下端面好，上端面骨組織生長快慢順序為右側、左側、后部、前部；下端面為左側、后部、右側、前部；患者術后常做旋轉運動，有利于骨骼生長。通過上述結果，期望本文研究可為特發性脊柱側凸手術方案的設計與融合器優化提供理論依據。
引用本文： 王召耀, 富榮昌, 馬原, 葉鵬. 特發性脊柱側凸骨單元宏細觀生物力學分析. 生物醫學工程學雜志, 2023, 40(2): 303-312. doi: 10.7507/1001-5515.202212053 復制
0 引言
特發性脊柱側凸（idiopathic scoliosis，IS）是指頭側端椎上緣垂線與尾側端椎下緣垂線交角 ≥ 10°[1]，且排除其他先天性、功能性、病理性等原因所產生的脊柱畸形[2]。IS在青少年中發病率為1%～3%[3]，占脊柱疾病的80%[4]，對患者外觀、心理健康皆有影響，會造成脊髓神經功能障礙，合并嚴重的胸廓畸形甚至會導致心肺功能衰竭甚至危及生命，個別患者還需進行截骨與個性化的手術。IS發病機制眾多[5]，同類型IS病例樣本收集難度較大，且對患者進行體外生物實驗涉及諸多倫理問題，幾乎無法收集同類型IS病例進行體外生物力學實驗，因此建立有限元模型并模擬生理工況，成為IS生物力學研究的主要手段。
隨著科技水平的進步，骨骼的多尺度分析已在生物醫學領域得到普遍的應用。但對骨細觀結構的研究大都為承載較為單一的長骨離體實驗，對脊柱等不規則骨的研究相對較少。IS有限元研究以宏觀尺度居多[6]，但術后骨細胞的轉移及骨組織的生長無法直接從宏觀尺度獲得；且人體正常生理活動所產生的宏觀載荷，會對骨組織細觀結構造成直接影響；骨組織細觀結構的力學性能變化，亦會在骨宏觀尺度上有所體現。因此，探究宏觀載荷與細觀結構之間的聯系，對掌握骨骼生物力學性能有著重要意義。
針對當前研究現狀，本文以IS病例為例，以畸形矯正后腰椎為研究對象，基于逆向工程、三維建模方法、有限元建模技術，建立宏觀腰椎與細觀骨單元模型。本文所建模型用以模擬人體生理活動最常見的直立、前屈、后伸、左旋、右旋、左屈與右屈，共７種工況。然后，利用子模型方法，建立宏觀模型與細觀結構的連接，分析宏觀皮質骨與細觀骨單元在相同邊界條件下力學性能的差異，從細觀尺度分析融合器的植入對周圍骨組織生長的影響，以期為IS手術的改進和融合器的優化設計提供理論依據。
1 材料與方法
1.1 建立畸形脊柱模型
模型原始數據于新疆醫科大學第六附屬醫院影像科采集，志愿者基本信息：女，23歲；維吾爾族；身高140 cm；體重57 kg；站立位上胸彎45°，主胸彎154°，仰臥位上胸彎41°，主胸彎159°，左側屈位上胸彎36°，主胸彎167°，右側屈位上胸彎52°，主胸彎144°。志愿者同意此次研究，并簽署知情同意書；本研究通過了新疆醫科大學第六附屬醫院倫理審查，獲得相關數據使用授權。
用計算機斷層掃描（computed tomography，CT）對其胸腰段脊柱進行螺旋掃描，橫向掃描層距0.64 mm，縱向掃描層距1 mm，獲得347張圖像。掃描結果以醫學數字成像與通信（digital imaging and communications in medicine，DICOM）格式保存。患者凹側為左側，屬于極重度脊柱畸形，胸第7～11節共5個椎體形成融合椎，幾乎無任何柔韌度，該病例是一例重度僵硬性脊柱側后凸畸形。
將該患者的DICOM文件導入醫學圖像處理軟件Mimics 19.0（Materialise Inc.，比利時）中，提取胸第1~腰第5椎體幾何輪廓；并導入逆向工程軟件Geomagic Wrap 2017（Geomagic Inc.，美國）中，重建椎體幾何模型輪廓，對復雜曲面進行擬合處理。將重建后的椎體曲面導入三維設計建模軟件SolidWorks 2018（Dassault Inc.，美國）中，構建出椎體間的纖維環、髓核與關節軟骨區域。纖維環與髓核在椎間盤中的體積占比分別為56%與44%[7]。由于本文研究節段較多，為減小計算量，建模時椎間盤僅考慮纖維環基質與髓核部分。將模型導入有限元分析軟件ABAQUS 6.14（Dassault Inc.，美國）中，根據解剖學結構，繪制前縱韌帶、后縱韌帶、橫突韌帶、棘間韌帶、棘上韌帶、黃韌帶、關節囊韌帶等7種韌帶[8]，完成術前患者的脊柱有限元模型建立。
1.2 腰椎矯正及宏觀模型建立
針對本文案例，新疆醫科大學第六附屬醫院預先設計了兩種矯形方案，并使用有限元法分別對兩種方案進行模擬計算。綜合對比各方案術后預期胸腰彎、胸彎、截骨數量、截骨難度以及內固定系統結構等指標，最終確定的矯形方案為：顯露脊柱后方結構，置釘節段為左側胸第2～7節椎體；腰第2～5節椎體；右側胸第2～5節椎體；腰第3～5節椎體。在腰1椎體進行截骨，以部分解除矯正側凸的“限制因素”，該步驟可以矯正局部脊柱后凸畸形和側凸畸形。因胸第7～11節椎體是融合椎，可看成一個僵硬的整體，且從臨床角度對其進行截骨操作風險極高，故未對這一融合體進行截骨操作，安裝凹側連接棒對相鄰螺釘之間進行適當撐開；安裝凸側連接棒，相鄰螺釘進行適當閉合操作，以矯正脊柱側凸畸形。然后，于截骨椎體處植入融合器以促進術后骨組織的生長。宏觀模型的建模過程如圖1所示。
圖1
IS手術建模過程
Figure1.
IS surgical modeling process
1.3 宏觀模型網格劃分與載荷約束
皮質骨采用三角形與矩形線性應變殼單元混合模擬，厚度設為1 mm[9]；松質骨采用四面體與五面體混合一階線性實體單元進行模擬；終板采用矩形線性應變殼單元模擬；纖維環基質、髓核與軟骨皆采用十節點四面體二階線性實體單元模擬；韌帶僅受拉不受壓，設置為超彈性、各向同性材料，使用一維非線性桿單元模擬。融合器采用六面體一階線性實體單元進行模擬，各結構的材料屬性如表1所示[10-15]。
表1
腰椎模型各部分單元選擇與材料參數
Table1.
Unit selection and material parameters for each part of the lumbar spine model
正常活動時，脊柱端面除了受相鄰椎體傳遞的重力法向分量，亦會在不同工況力矩的綜合作用下產生相應運動。以前文中體重57 kg的患者為研究對象，探究融合器與椎骨各工況下的力學性能。根據人體測量方法，如圖2所示，對腰5椎體的下表面施加全約束，限制其6個自由度；沿腰1椎體上端面法向施加380 N的均布載荷，模擬直立時的受力情況；再建立6個力矩，分別對應前屈、后伸、左旋、右旋、左屈與右屈等6個工況；力矩大小為5.7 N·m[16]。
圖2
宏觀模型載荷與約束
Figure2.
Macro model loads and constraints
1.4 建立皮質骨骨單元細觀尺度模型
宏觀尺度為10-2 m量級，所建立細觀模型為10-4 m量級，骨單元為細觀尺度下骨的主要組成部分。骨單元由哈佛氏管以及周圍的同心骨板、骨板介質、粘合線組成，據文獻[17]記載，包絡哈佛氏管的同心骨板為4～20層不等，皮質骨骨板橫截面上骨單元的面密度為11～20 個/mm2。根據統計學樣本分析[18]，設單個骨單元有6層同心骨板，面密度為16 個/mm2，骨單元直徑為200 μm，哈佛氏管直徑為60 μm，同心骨板厚度為10 μm，粘合線的厚度為5 μm，骨板介質厚度為1 μm。
本文選取截骨腰1椎體處8個骨單元進行細觀研究，腰1上半椎體下端面與下半椎體上端面各取4個骨單元，上半椎體骨單元以逆時針順序選取左、后、右、前各一個，給予編號S1～S4；下半椎體骨單元以同樣方法選取，編號為X1～X4。分別繪制骨板、骨板介質、粘合線的三維模型。將模型導入有限元分析軟件ABAQUS 6.14（Dassault Inc.，美國）中，骨單元結構較為規則，同心骨板與粘合線皆設為一階線性實體單元，骨板介質厚度較小，設為矩形線性應變殼單元。骨單元各部分采用共節點連接。皮質骨骨單元細觀尺度建模過程如圖3所示。
圖3
細觀骨單元建模過程
Figure3.
Mesoscopic bone unit modeling process
1.5 細觀模型網格劃分與載荷約束
為了研究皮質骨宏觀與細觀力學性能的差異，應控制其宏觀模型與細觀模型的邊界條件與載荷相一致。宏觀模型的邊界條件施加在腰5椎體處，各工況下載荷施加于腰1椎體上表面。細觀模型在尺寸上較宏觀模型小2個數量級，細觀模型空間位置與宏觀模型一致。但與宏觀模型相比，細觀模型幾何信息中沒有腰5椎體下表面與腰1椎體上表面，無法直接對細觀模型施加與宏觀模型一致的邊界條件與載荷。若將整個腰椎皆以骨單元尺度進行建模分析，其工作量會增加且計算的結果必然存在誤差。例如，將松質骨小梁與骨單元的區域過度、椎間盤終板與骨單元的橋接、骨小梁和終板與骨單元三者接觸處耦合等問題在整體模型中進行綜合考慮，將極大地增加計算的工作量，且并非此次研究的重點。
為了精確且高效地建立宏觀模型與細觀結構的聯系，本文引入了基于圣維南原理的子模型功能。子模型法又稱切割邊界位移法或特定邊界位移法，是獲得模型部分區域更精確解的有限單元技術。將骨單元看作整體腰椎的一部分，視腰椎的宏觀模型為整體模型，骨單元為子模型。切割邊界即為子模型輪廓被整體模型所包絡的部分，整體模型于切割邊界處計算所得位移值即為子模型的初始邊界條件。即，將宏觀腰椎模型的計算結果當作細觀骨單元模型的邊界條件，進行細觀骨單元模型的計算。且子模型是與整體模型相互獨立存在，而并非像片段一樣依存于整體模型。采用子模型法進行分析，可以減小甚至消除整體模型內部復雜載荷的傳遞，使得計算結果更加精確。
骨單元各組織均設為連續均勻且各向同性的線彈性材料，各組織的材料參數如表2所示[18-20]。
表2
骨單元各部分單元選擇與材料參數
Table2.
Unit selection and material parameters for each part of the bone unit
2 結果
2.1 腰椎模型有效性驗證
本文模型為畸形腰椎，國內外暫無同類型IS生物力學加載實驗做對比，無法直接對整體腰椎進行有效性驗證，所以本文通過與分段加載實驗對比，進行模型有效性驗證。取相對健康的腰5、腰4節段進行驗證，對腰5椎體的下表面施加約束，固定其6個自由度；在腰4椎體的上表面分別均勻施加500、1 000、1 500、2 000 N的軸向載荷，觀察其軸向位移[21]。然后，將計算結果與實驗結果相對比，從而驗證模型的有效性。
在500、1 000、1 500、2 000 N的均布軸向壓力下，椎體的軸向位移分別為0.35、0.54、0.81、1.08 mm。如圖4所示，將其軸向壓力—軸向位移曲線與文獻[22-25]相對比。相同條件下，離體實驗與有限元分析二者結果相接近，故模型的可靠性與有效性得以驗證。
圖4
不同均布壓力下椎體軸向位移曲線
Figure4.
Axial displacement curve of vertebral body under different uniform pressures
2.2 宏觀模型有限元結果
宏觀模型整體應力與骨單元宏觀尺度各工況最大應力如圖5所示。
圖5
宏觀模型有限元結果
Figure5.
Macroscopic model finite element results
宏觀骨單元最大應力取于骨單元S1左屈時，最小應力取于骨單元X4直立時。上半骨單元S1～S4平均應力為4.007 MPa；下半骨單元X1～X4平均應力為1.936 MPa。宏觀尺度下，上半骨單元的平均應力較下半骨單元大。
以不同工況角度分析宏觀骨單元整體的平均應力。在直立、前屈、后伸、左旋、右旋、左屈與右屈等7種工況下，骨單元整體所對應的平均應力分別為1.391、3.113、2.994、4.034、3.734、3.263、2.272 MPa，旋轉工況下，宏觀骨單元整體的平均應力最大。
2.3 細觀模型有限元結果
細觀骨單元應力云圖結果如圖6所示，骨單元由同心骨板、骨板介質、粘合線三部分組成，各部分細觀尺度應力結果如圖7所示。
圖6
細觀骨單元應力云圖
Figure6.
Stress clouds of mesoscopic bone unit
圖7
細觀骨單元各結構最大應力
Figure7.
The maximum equivalent stress of each structure of the mesoscopic bone unit
細觀骨單元整體應力較宏觀模型有著較大的提高。細觀骨單元的最大應力值取S3骨單元右屈時，且S1、S2骨單元在大部分工況下皆受到較大應力。S1～S4骨單元平均應力為14.529 MPa，X1～X4骨單元平均應力為7.587 MPa，與宏觀結果相同，細觀尺度下，上半骨單元的平均應力亦較下半骨單元大。
分析不同工況下細觀骨單元整體的平均應力。在直立、前屈、后伸、左旋、右旋、左屈與右屈等7種工況下，所對應的平均應力分別為5.889、11.826、12.078、14.416、12.410、12.279、9.330 MPa，細觀尺度下，骨單元整體應力在旋轉工況下最大，這與宏觀尺度結果一致。
即使處于同一端面的骨單元，因所處空間位置不同，骨單元間的力學性能也有所不同。上半椎體端面的4個骨單元中，S3骨單元各工況下平均應力最大，為18.110 MPa，其次為S1、S2骨單元，S4骨單元平均應力最小，為10.040 MPa。以空間位置排序，腰1上半椎體端面平均應力從大到小分別為右側、左側、后部、前部。下半椎體中，X1骨單元平均應力最大，為11.611 MPa，其次為X2、X3，平均應力最小的骨單元為X4，其值為5.348 MPa。下半椎體端面平均應力由大到小分別為左側、后部、右側、前部。
骨單元細觀與宏觀尺度下的最大應力比如圖8所示。宏觀骨單元材料屬性與皮質骨一致，細觀骨單元由骨板、骨板介質、粘合線組成，其中骨板應力為三者最大，故骨單元細觀與宏觀最大應力比可看作細觀骨板與對應宏觀皮質骨的比值。應力比根據工況與骨單元空間位置不同有所改變，最小應力比為2.606，位于前屈工況下X4骨單元處；最大應力比為5.958，位于左旋工況下X2骨單元處。腰1椎體上半骨單元平均應力比為3.687，下半骨單元應力比為4.113，上半骨單元平均應力比較下半骨單元小。
圖8
細宏觀骨單元最大應力比
Figure8.
Ratio of the maximum equivalent stress of the mesoscopic bone unit to the macroscopic bone unit
3 討論
在IS矯形后7種工況不同尺度的對比分析中，不論宏觀與細觀，旋轉工況整體應力皆為各工況最大。同一工況下，腰椎皮質骨細觀尺度的應力較宏觀尺度有所增加，本文中細觀骨單元應力約為宏觀的2.606～5.958倍。這種應力增大，是由于骨單元各部分材料性質不同所造成的。相較性質均一的宏觀皮質骨，骨板的楊氏模量高于皮質骨，而粘合線與骨板介質部分則低于皮質骨。在相同的邊界條件下，楊氏模量與應力呈正相關。同理，無論何種工況，骨單元的應力最大值均出現在同心骨板處，且大都為靠近哈佛氏管的最內側骨板。
宏觀尺度下，腰1上半椎體骨單元平均應力較下半骨單元大；細觀尺度下，上半骨單元平均應力亦大于下部。從宏觀模型來說，腰1上半椎體上端面承擔各工況的載荷，下端面與椎間融合器直接接觸；腰1下半椎體則相反，上端承受融合器所傳遞載荷，下端面與椎間盤相接觸。融合器的楊氏模量較椎間盤有著數量級上的增加，使上半椎體承受較大載荷。加之在載荷的傳遞過程中部分應變能被融合器吸收，使得宏觀尺度下，上半骨單元平均應力大于下半骨單元。細觀骨單元載荷與邊界條件與宏觀模型相一致，所得結果與宏觀類似。
腰1上半椎體端面平均應力按大小排序分別為右側、左側、后部、前部；下半椎體平均應力排序為左側、后部、右側、前部。這是由于截骨后腰1椎體與融合器接觸端面的輪廓不規則，融合器系統結構偏心所造成的。且腰1椎體受載端面的輪廓不規則，進一步加劇了融合器的受力不均、截面變形不同，使得融合器中部偏右處應力較高。又因腰1上截骨面較下截骨面小，融合器上表面較下表面更易形成較大應力。反映在與融合器接觸的椎體端面處，則有端面四周應力分布不均的現象。據此推測，因受力不同，椎體端面骨組織的生長速度亦有所差異。僅從力學角度推測，融合器上端面骨組織生長快慢順序為右側、左側、后部、前部；下端面為左側、后部、右側、前部。此種生長速度差異，與蔣成明等[26]臨床研究結果相一致。
在細觀骨單元中，應力最大值在骨板處。根據生物細胞學可知，骨單元高應力處將會刺激該區域組織液的流動，從而刺激細胞再生[27]。粘合線與部分骨板有一定程度上的應力集中。在相同邊界條件下，材料屬性分布不均更易造成應力集中。根據伍爾夫定律可知，骨在需要的地方生長，在不需要的地方消亡。骨單元的高應力區域，促進了骨胞元向成骨細胞的轉化，使得骨在此區域的密度以及硬度增加。骨板介質與粘合線處的應力相較骨板有著顯著的降低，這種應力中斷現象，對骨組織的自我保護，起到了積極的作用。
綜上所述，得到如下結論：① 腰椎細觀尺度的應力較宏觀有所增加，本文中細觀骨單元應力約為宏觀的2.606～5.958倍。② 與融合器接觸的椎體端面，上端面平均應力較下端面大。根據伍爾夫定律，上端面骨組織的成骨性較下端面好。③ 各椎體端面四周骨組織的生長順序有所不同，上端面骨組織生長快慢順序為右側、左側、后部、前部；下端面為左側、后部、右側、前部。④ 各工況綜合對比下，旋轉工況整體應力值較大，該患者術后常做旋轉運動，有利于骨骼生長。
目前，IS生物力學特性的研究以宏觀尺度居多，而骨細觀力學性能的研究多以股骨等長骨為研究對象。本研究以不規則骨——椎骨為研究對象，建立其骨單元細觀模型，從細觀尺度探究IS生物力學特性。本研究利用子模型法，獲得截骨椎骨端面四周應力，據此推測骨組織生長速度的差異。據此研究，可使IS患者進行針對性的康復運動，相關研究人員也可據此研究對椎間植入物進行針對性優化。本研究針對皮質骨所采用的子模型細觀建模方法，不僅適用于長骨與不規則骨，對短骨與扁骨也同樣適用。且此模型中邊界條件的選取規則，也可為后續同類研究提供一定的參考依據。
本研究亦存在一些局限性，宏觀尺度下的模型，未考慮肌肉等組織的影響，所建立的微米級細觀骨單元模型，未考慮四周不規則的間骨板與骨陷窩等尺度更小結構。有關骨組織生長機制還需從骨細胞層次進一步剖析，后續研究應將模型繼續細化，并嘗試建立微觀尺度下納米級模型，為宏觀載荷對微觀尺度的影響，提供參考依據。
重要聲明
利益沖突聲明：本文全體作者均聲明不存在利益沖突。
作者貢獻聲明：王召耀進行了三維建模、有限元分析與論文寫作；富榮昌對論文內容進行修改；馬原進行了方案設計；葉鵬負責數據收集與分析。
倫理聲明：本研究通過了新疆醫科大學第六附屬醫院倫理委員會審批（編號：LFYLLSC20170301-02）。
0 引言
特發性脊柱側凸（idiopathic scoliosis，IS）是指頭側端椎上緣垂線與尾側端椎下緣垂線交角 ≥ 10°[1]，且排除其他先天性、功能性、病理性等原因所產生的脊柱畸形[2]。IS在青少年中發病率為1%～3%[3]，占脊柱疾病的80%[4]，對患者外觀、心理健康皆有影響，會造成脊髓神經功能障礙，合并嚴重的胸廓畸形甚至會導致心肺功能衰竭甚至危及生命，個別患者還需進行截骨與個性化的手術。IS發病機制眾多[5]，同類型IS病例樣本收集難度較大，且對患者進行體外生物實驗涉及諸多倫理問題，幾乎無法收集同類型IS病例進行體外生物力學實驗，因此建立有限元模型并模擬生理工況，成為IS生物力學研究的主要手段。
隨著科技水平的進步，骨骼的多尺度分析已在生物醫學領域得到普遍的應用。但對骨細觀結構的研究大都為承載較為單一的長骨離體實驗，對脊柱等不規則骨的研究相對較少。IS有限元研究以宏觀尺度居多[6]，但術后骨細胞的轉移及骨組織的生長無法直接從宏觀尺度獲得；且人體正常生理活動所產生的宏觀載荷，會對骨組織細觀結構造成直接影響；骨組織細觀結構的力學性能變化，亦會在骨宏觀尺度上有所體現。因此，探究宏觀載荷與細觀結構之間的聯系，對掌握骨骼生物力學性能有著重要意義。
針對當前研究現狀，本文以IS病例為例，以畸形矯正后腰椎為研究對象，基于逆向工程、三維建模方法、有限元建模技術，建立宏觀腰椎與細觀骨單元模型。本文所建模型用以模擬人體生理活動最常見的直立、前屈、后伸、左旋、右旋、左屈與右屈，共７種工況。然后，利用子模型方法，建立宏觀模型與細觀結構的連接，分析宏觀皮質骨與細觀骨單元在相同邊界條件下力學性能的差異，從細觀尺度分析融合器的植入對周圍骨組織生長的影響，以期為IS手術的改進和融合器的優化設計提供理論依據。
1 材料與方法
1.1 建立畸形脊柱模型
模型原始數據于新疆醫科大學第六附屬醫院影像科采集，志愿者基本信息：女，23歲；維吾爾族；身高140 cm；體重57 kg；站立位上胸彎45°，主胸彎154°，仰臥位上胸彎41°，主胸彎159°，左側屈位上胸彎36°，主胸彎167°，右側屈位上胸彎52°，主胸彎144°。志愿者同意此次研究，并簽署知情同意書；本研究通過了新疆醫科大學第六附屬醫院倫理審查，獲得相關數據使用授權。
用計算機斷層掃描（computed tomography，CT）對其胸腰段脊柱進行螺旋掃描，橫向掃描層距0.64 mm，縱向掃描層距1 mm，獲得347張圖像。掃描結果以醫學數字成像與通信（digital imaging and communications in medicine，DICOM）格式保存。患者凹側為左側，屬于極重度脊柱畸形，胸第7～11節共5個椎體形成融合椎，幾乎無任何柔韌度，該病例是一例重度僵硬性脊柱側后凸畸形。
將該患者的DICOM文件導入醫學圖像處理軟件Mimics 19.0（Materialise Inc.，比利時）中，提取胸第1~腰第5椎體幾何輪廓；并導入逆向工程軟件Geomagic Wrap 2017（Geomagic Inc.，美國）中，重建椎體幾何模型輪廓，對復雜曲面進行擬合處理。將重建后的椎體曲面導入三維設計建模軟件SolidWorks 2018（Dassault Inc.，美國）中，構建出椎體間的纖維環、髓核與關節軟骨區域。纖維環與髓核在椎間盤中的體積占比分別為56%與44%[7]。由于本文研究節段較多，為減小計算量，建模時椎間盤僅考慮纖維環基質與髓核部分。將模型導入有限元分析軟件ABAQUS 6.14（Dassault Inc.，美國）中，根據解剖學結構，繪制前縱韌帶、后縱韌帶、橫突韌帶、棘間韌帶、棘上韌帶、黃韌帶、關節囊韌帶等7種韌帶[8]，完成術前患者的脊柱有限元模型建立。
1.2 腰椎矯正及宏觀模型建立
針對本文案例，新疆醫科大學第六附屬醫院預先設計了兩種矯形方案，并使用有限元法分別對兩種方案進行模擬計算。綜合對比各方案術后預期胸腰彎、胸彎、截骨數量、截骨難度以及內固定系統結構等指標，最終確定的矯形方案為：顯露脊柱后方結構，置釘節段為左側胸第2～7節椎體；腰第2～5節椎體；右側胸第2～5節椎體；腰第3～5節椎體。在腰1椎體進行截骨，以部分解除矯正側凸的“限制因素”，該步驟可以矯正局部脊柱后凸畸形和側凸畸形。因胸第7～11節椎體是融合椎，可看成一個僵硬的整體，且從臨床角度對其進行截骨操作風險極高，故未對這一融合體進行截骨操作，安裝凹側連接棒對相鄰螺釘之間進行適當撐開；安裝凸側連接棒，相鄰螺釘進行適當閉合操作，以矯正脊柱側凸畸形。然后，于截骨椎體處植入融合器以促進術后骨組織的生長。宏觀模型的建模過程如圖1所示。
圖1
IS手術建模過程
Figure1.
IS surgical modeling process
1.3 宏觀模型網格劃分與載荷約束
皮質骨采用三角形與矩形線性應變殼單元混合模擬，厚度設為1 mm[9]；松質骨采用四面體與五面體混合一階線性實體單元進行模擬；終板采用矩形線性應變殼單元模擬；纖維環基質、髓核與軟骨皆采用十節點四面體二階線性實體單元模擬；韌帶僅受拉不受壓，設置為超彈性、各向同性材料，使用一維非線性桿單元模擬。融合器采用六面體一階線性實體單元進行模擬，各結構的材料屬性如表1所示[10-15]。
表1
腰椎模型各部分單元選擇與材料參數
Table1.
Unit selection and material parameters for each part of the lumbar spine model
正常活動時，脊柱端面除了受相鄰椎體傳遞的重力法向分量，亦會在不同工況力矩的綜合作用下產生相應運動。以前文中體重57 kg的患者為研究對象，探究融合器與椎骨各工況下的力學性能。根據人體測量方法，如圖2所示，對腰5椎體的下表面施加全約束，限制其6個自由度；沿腰1椎體上端面法向施加380 N的均布載荷，模擬直立時的受力情況；再建立6個力矩，分別對應前屈、后伸、左旋、右旋、左屈與右屈等6個工況；力矩大小為5.7 N·m[16]。
圖2
宏觀模型載荷與約束
Figure2.
Macro model loads and constraints
1.4 建立皮質骨骨單元細觀尺度模型
宏觀尺度為10-2 m量級，所建立細觀模型為10-4 m量級，骨單元為細觀尺度下骨的主要組成部分。骨單元由哈佛氏管以及周圍的同心骨板、骨板介質、粘合線組成，據文獻[17]記載，包絡哈佛氏管的同心骨板為4～20層不等，皮質骨骨板橫截面上骨單元的面密度為11～20 個/mm2。根據統計學樣本分析[18]，設單個骨單元有6層同心骨板，面密度為16 個/mm2，骨單元直徑為200 μm，哈佛氏管直徑為60 μm，同心骨板厚度為10 μm，粘合線的厚度為5 μm，骨板介質厚度為1 μm。
本文選取截骨腰1椎體處8個骨單元進行細觀研究，腰1上半椎體下端面與下半椎體上端面各取4個骨單元，上半椎體骨單元以逆時針順序選取左、后、右、前各一個，給予編號S1～S4；下半椎體骨單元以同樣方法選取，編號為X1～X4。分別繪制骨板、骨板介質、粘合線的三維模型。將模型導入有限元分析軟件ABAQUS 6.14（Dassault Inc.，美國）中，骨單元結構較為規則，同心骨板與粘合線皆設為一階線性實體單元，骨板介質厚度較小，設為矩形線性應變殼單元。骨單元各部分采用共節點連接。皮質骨骨單元細觀尺度建模過程如圖3所示。
圖3
細觀骨單元建模過程
Figure3.
Mesoscopic bone unit modeling process
1.5 細觀模型網格劃分與載荷約束
為了研究皮質骨宏觀與細觀力學性能的差異，應控制其宏觀模型與細觀模型的邊界條件與載荷相一致。宏觀模型的邊界條件施加在腰5椎體處，各工況下載荷施加于腰1椎體上表面。細觀模型在尺寸上較宏觀模型小2個數量級，細觀模型空間位置與宏觀模型一致。但與宏觀模型相比，細觀模型幾何信息中沒有腰5椎體下表面與腰1椎體上表面，無法直接對細觀模型施加與宏觀模型一致的邊界條件與載荷。若將整個腰椎皆以骨單元尺度進行建模分析，其工作量會增加且計算的結果必然存在誤差。例如，將松質骨小梁與骨單元的區域過度、椎間盤終板與骨單元的橋接、骨小梁和終板與骨單元三者接觸處耦合等問題在整體模型中進行綜合考慮，將極大地增加計算的工作量，且并非此次研究的重點。
為了精確且高效地建立宏觀模型與細觀結構的聯系，本文引入了基于圣維南原理的子模型功能。子模型法又稱切割邊界位移法或特定邊界位移法，是獲得模型部分區域更精確解的有限單元技術。將骨單元看作整體腰椎的一部分，視腰椎的宏觀模型為整體模型，骨單元為子模型。切割邊界即為子模型輪廓被整體模型所包絡的部分，整體模型于切割邊界處計算所得位移值即為子模型的初始邊界條件。即，將宏觀腰椎模型的計算結果當作細觀骨單元模型的邊界條件，進行細觀骨單元模型的計算。且子模型是與整體模型相互獨立存在，而并非像片段一樣依存于整體模型。采用子模型法進行分析，可以減小甚至消除整體模型內部復雜載荷的傳遞，使得計算結果更加精確。
骨單元各組織均設為連續均勻且各向同性的線彈性材料，各組織的材料參數如表2所示[18-20]。
表2
骨單元各部分單元選擇與材料參數
Table2.
Unit selection and material parameters for each part of the bone unit
2 結果
2.1 腰椎模型有效性驗證
本文模型為畸形腰椎，國內外暫無同類型IS生物力學加載實驗做對比，無法直接對整體腰椎進行有效性驗證，所以本文通過與分段加載實驗對比，進行模型有效性驗證。取相對健康的腰5、腰4節段進行驗證，對腰5椎體的下表面施加約束，固定其6個自由度；在腰4椎體的上表面分別均勻施加500、1 000、1 500、2 000 N的軸向載荷，觀察其軸向位移[21]。然后，將計算結果與實驗結果相對比，從而驗證模型的有效性。
在500、1 000、1 500、2 000 N的均布軸向壓力下，椎體的軸向位移分別為0.35、0.54、0.81、1.08 mm。如圖4所示，將其軸向壓力—軸向位移曲線與文獻[22-25]相對比。相同條件下，離體實驗與有限元分析二者結果相接近，故模型的可靠性與有效性得以驗證。
圖4
不同均布壓力下椎體軸向位移曲線
Figure4.
Axial displacement curve of vertebral body under different uniform pressures
2.2 宏觀模型有限元結果
宏觀模型整體應力與骨單元宏觀尺度各工況最大應力如圖5所示。
圖5
宏觀模型有限元結果
Figure5.
Macroscopic model finite element results
宏觀骨單元最大應力取于骨單元S1左屈時，最小應力取于骨單元X4直立時。上半骨單元S1～S4平均應力為4.007 MPa；下半骨單元X1～X4平均應力為1.936 MPa。宏觀尺度下，上半骨單元的平均應力較下半骨單元大。
以不同工況角度分析宏觀骨單元整體的平均應力。在直立、前屈、后伸、左旋、右旋、左屈與右屈等7種工況下，骨單元整體所對應的平均應力分別為1.391、3.113、2.994、4.034、3.734、3.263、2.272 MPa，旋轉工況下，宏觀骨單元整體的平均應力最大。
2.3 細觀模型有限元結果
細觀骨單元應力云圖結果如圖6所示，骨單元由同心骨板、骨板介質、粘合線三部分組成，各部分細觀尺度應力結果如圖7所示。
圖6
細觀骨單元應力云圖
Figure6.
Stress clouds of mesoscopic bone unit
圖7
細觀骨單元各結構最大應力
Figure7.
The maximum equivalent stress of each structure of the mesoscopic bone unit
細觀骨單元整體應力較宏觀模型有著較大的提高。細觀骨單元的最大應力值取S3骨單元右屈時，且S1、S2骨單元在大部分工況下皆受到較大應力。S1～S4骨單元平均應力為14.529 MPa，X1～X4骨單元平均應力為7.587 MPa，與宏觀結果相同，細觀尺度下，上半骨單元的平均應力亦較下半骨單元大。
分析不同工況下細觀骨單元整體的平均應力。在直立、前屈、后伸、左旋、右旋、左屈與右屈等7種工況下，所對應的平均應力分別為5.889、11.826、12.078、14.416、12.410、12.279、9.330 MPa，細觀尺度下，骨單元整體應力在旋轉工況下最大，這與宏觀尺度結果一致。
即使處于同一端面的骨單元，因所處空間位置不同，骨單元間的力學性能也有所不同。上半椎體端面的4個骨單元中，S3骨單元各工況下平均應力最大，為18.110 MPa，其次為S1、S2骨單元，S4骨單元平均應力最小，為10.040 MPa。以空間位置排序，腰1上半椎體端面平均應力從大到小分別為右側、左側、后部、前部。下半椎體中，X1骨單元平均應力最大，為11.611 MPa，其次為X2、X3，平均應力最小的骨單元為X4，其值為5.348 MPa。下半椎體端面平均應力由大到小分別為左側、后部、右側、前部。
骨單元細觀與宏觀尺度下的最大應力比如圖8所示。宏觀骨單元材料屬性與皮質骨一致，細觀骨單元由骨板、骨板介質、粘合線組成，其中骨板應力為三者最大，故骨單元細觀與宏觀最大應力比可看作細觀骨板與對應宏觀皮質骨的比值。應力比根據工況與骨單元空間位置不同有所改變，最小應力比為2.606，位于前屈工況下X4骨單元處；最大應力比為5.958，位于左旋工況下X2骨單元處。腰1椎體上半骨單元平均應力比為3.687，下半骨單元應力比為4.113，上半骨單元平均應力比較下半骨單元小。
圖8
細宏觀骨單元最大應力比
Figure8.
Ratio of the maximum equivalent stress of the mesoscopic bone unit to the macroscopic bone unit
3 討論
在IS矯形后7種工況不同尺度的對比分析中，不論宏觀與細觀，旋轉工況整體應力皆為各工況最大。同一工況下，腰椎皮質骨細觀尺度的應力較宏觀尺度有所增加，本文中細觀骨單元應力約為宏觀的2.606～5.958倍。這種應力增大，是由于骨單元各部分材料性質不同所造成的。相較性質均一的宏觀皮質骨，骨板的楊氏模量高于皮質骨，而粘合線與骨板介質部分則低于皮質骨。在相同的邊界條件下，楊氏模量與應力呈正相關。同理，無論何種工況，骨單元的應力最大值均出現在同心骨板處，且大都為靠近哈佛氏管的最內側骨板。
宏觀尺度下，腰1上半椎體骨單元平均應力較下半骨單元大；細觀尺度下，上半骨單元平均應力亦大于下部。從宏觀模型來說，腰1上半椎體上端面承擔各工況的載荷，下端面與椎間融合器直接接觸；腰1下半椎體則相反，上端承受融合器所傳遞載荷，下端面與椎間盤相接觸。融合器的楊氏模量較椎間盤有著數量級上的增加，使上半椎體承受較大載荷。加之在載荷的傳遞過程中部分應變能被融合器吸收，使得宏觀尺度下，上半骨單元平均應力大于下半骨單元。細觀骨單元載荷與邊界條件與宏觀模型相一致，所得結果與宏觀類似。
腰1上半椎體端面平均應力按大小排序分別為右側、左側、后部、前部；下半椎體平均應力排序為左側、后部、右側、前部。這是由于截骨后腰1椎體與融合器接觸端面的輪廓不規則，融合器系統結構偏心所造成的。且腰1椎體受載端面的輪廓不規則，進一步加劇了融合器的受力不均、截面變形不同，使得融合器中部偏右處應力較高。又因腰1上截骨面較下截骨面小，融合器上表面較下表面更易形成較大應力。反映在與融合器接觸的椎體端面處，則有端面四周應力分布不均的現象。據此推測，因受力不同，椎體端面骨組織的生長速度亦有所差異。僅從力學角度推測，融合器上端面骨組織生長快慢順序為右側、左側、后部、前部；下端面為左側、后部、右側、前部。此種生長速度差異，與蔣成明等[26]臨床研究結果相一致。
在細觀骨單元中，應力最大值在骨板處。根據生物細胞學可知，骨單元高應力處將會刺激該區域組織液的流動，從而刺激細胞再生[27]。粘合線與部分骨板有一定程度上的應力集中。在相同邊界條件下，材料屬性分布不均更易造成應力集中。根據伍爾夫定律可知，骨在需要的地方生長，在不需要的地方消亡。骨單元的高應力區域，促進了骨胞元向成骨細胞的轉化，使得骨在此區域的密度以及硬度增加。骨板介質與粘合線處的應力相較骨板有著顯著的降低，這種應力中斷現象，對骨組織的自我保護，起到了積極的作用。
綜上所述，得到如下結論：① 腰椎細觀尺度的應力較宏觀有所增加，本文中細觀骨單元應力約為宏觀的2.606～5.958倍。② 與融合器接觸的椎體端面，上端面平均應力較下端面大。根據伍爾夫定律，上端面骨組織的成骨性較下端面好。③ 各椎體端面四周骨組織的生長順序有所不同，上端面骨組織生長快慢順序為右側、左側、后部、前部；下端面為左側、后部、右側、前部。④ 各工況綜合對比下，旋轉工況整體應力值較大，該患者術后常做旋轉運動，有利于骨骼生長。
目前，IS生物力學特性的研究以宏觀尺度居多，而骨細觀力學性能的研究多以股骨等長骨為研究對象。本研究以不規則骨——椎骨為研究對象，建立其骨單元細觀模型，從細觀尺度探究IS生物力學特性。本研究利用子模型法，獲得截骨椎骨端面四周應力，據此推測骨組織生長速度的差異。據此研究，可使IS患者進行針對性的康復運動，相關研究人員也可據此研究對椎間植入物進行針對性優化。本研究針對皮質骨所采用的子模型細觀建模方法，不僅適用于長骨與不規則骨，對短骨與扁骨也同樣適用。且此模型中邊界條件的選取規則，也可為后續同類研究提供一定的參考依據。
本研究亦存在一些局限性，宏觀尺度下的模型，未考慮肌肉等組織的影響，所建立的微米級細觀骨單元模型，未考慮四周不規則的間骨板與骨陷窩等尺度更小結構。有關骨組織生長機制還需從骨細胞層次進一步剖析，后續研究應將模型繼續細化，并嘗試建立微觀尺度下納米級模型，為宏觀載荷對微觀尺度的影響，提供參考依據。
重要聲明
利益沖突聲明：本文全體作者均聲明不存在利益沖突。
作者貢獻聲明：王召耀進行了三維建模、有限元分析與論文寫作；富榮昌對論文內容進行修改；馬原進行了方案設計；葉鵬負責數據收集與分析。
倫理聲明：本研究通過了新疆醫科大學第六附屬醫院倫理委員會審批（編號：LFYLLSC20170301-02）。
