載荷作用下松質骨孔隙中的液體流動是刺激骨組織細胞產生生物學響應并調控骨重建的主要因素。因此，闡明牙槽骨內孔隙結構中的液體流動情況對于深入理解力學作用在牙槽骨內的傳導過程以及牙齒發育、正畸牙移動等細胞水平的調控機制具有重要意義。此工作首先進行了大鼠牙齒正畸的動物實驗，并基于微計算機斷層掃描（micro-CT）圖像構建了牙齒-牙周韌帶-牙槽骨有限元模型，分析了咬合力或正畸力作用下牙槽骨中的應變狀態；進而構建了理想模型，應用流固耦合數值模擬方法，分析了動態咬合力加載下無正畸加載、正畸拉伸加載、正畸壓縮加載三種情況下骨內液體的流動情況。模擬結果表明，動態咬合力作用下，沿咬合方向排列的骨小梁表面流體剪應力水平高于非咬合方向排列的骨小梁，正畸力對骨內液體的流動沒有影響。上述結果說明，臨床上通過調整牙齒咬合面形狀等方法改變咬合力的方向，會在牙槽骨表面引起不同水平的流體剪應力，進而刺激骨組織表面的細胞產生響應，最終調控牙槽骨的結構重建。
引用本文： 羅睿, 趙振達, 冷慧杰, 霍波. 大鼠牙槽骨理想模型的流固耦合數值模擬研究. 生物醫學工程學雜志, 2020, 37(1): 87-95. doi: 10.7507/1001-5515.201903019 復制
引言
牙槽骨包括固有牙槽骨、密質骨及松質骨等幾部分，其中松質骨位于密質骨和固有牙槽骨之間，由骨小梁和骨髓組成。牙槽骨的結構變化非常活躍，與牙齒發育、乳牙替換、恒牙移動等生理過程密切相關。這種結構的重建過程可由牙齒受到的載荷（如咬合力或正畸力）調控，主要表現為骨小梁的排列方向與咬合力相適應，以最有效的排列方式來抵抗外加載荷。實際上，骨的結構重建主要由成骨細胞主導的骨形成過程及破骨細胞主導的骨吸收過程組成。在咬合力或正畸力的作用下，骨基質的變形會引起附著于其上的細胞產生相應的應變；松質骨內孔隙體積的改變也會引起孔隙內液體的流動，從而在細胞表面產生流體剪應力[1]。許多研究已經證明，相對于基質的變形，流體刺激會在骨組織細胞中引起更為顯著的生物學響應[2-4]。因此，闡明牙槽骨內孔隙結構中的液體流動情況對于深入理解力學作用在牙槽骨內的傳導過程以及牙齒的發育形成、正畸牙移動等細胞水平的調控機制具有重要意義。
受限于在體測量實驗技術的缺乏，數值模擬是研究骨內液體流動空間分布和動態變化的有效方法。2012 年有學者構建了包含骨小梁孔隙結構的理想松質骨模型，使用流固耦合數值模擬方法計算了其中的液體流動、壓力和流體剪應力等參數的分布[5]。2015 年的另一項研究使用微計算機斷層掃描技術（micro computed tomography，micro-CT）掃描豬的股骨，通過三維重建得到其真實的松質骨幾何模型，也應用流固耦合數值模擬方法研究了骨小梁孔隙結構中的流體剪應力和壓力分布[6]。但目前還沒有關于牙槽骨內液體流動情況的數值模擬研究，因此咬合力或正畸力作用下牙槽骨內骨小梁表面流體剪應力的大小和分布尚不清楚，從而難于判斷牙槽骨中細胞所處的力學微環境。
本文針對大鼠的正畸牙移動模型，利用有限元方法分析了牙槽骨中的應變分布，進而構建了理想模型，應用流固耦合數值模擬方法，得到了動態咬合力加載下，無正畸力、正畸壓縮和正畸拉伸三種加載情況下骨內液體的流動情況，對比分析了咬合力方向、非咬合力方向骨小梁表面流體剪應力的差異，以及拉伸、壓縮對骨小梁表面流體剪應力的影響。
1 材料和方法
1.1 動物實驗
本實驗從北京大學醫學部實驗動物科學部購置一只 6～8 周齡的 SPF 級雌性 SD 大鼠［許可證號：SCXK（京）2016-0010］。將大鼠麻醉后置于手術臺上，以上頜 2 顆切牙作為支抗牙，使用鎳鈦拉簧對其上頜右側第一磨牙施加 0.8 N 正畸力 20 天（見圖 1a）。鎳鈦拉簧的一端通過直徑為 0.2 mm 的結扎絲固定于第一磨牙與第二磨牙之間，另一端連接測力計并測定正畸力，之后將其固定在兩顆切牙頸部 0.5～1.0 mm 的固位槽中。進而將光固化樹脂均勻覆蓋于切牙表面和固位溝處，并用 LED 光固化燈固化樹脂 40 s。將大鼠下頜切牙磨短并用樹脂覆蓋斷面處，以防止其磨斷口腔內的結扎絲。在正畸治療期間，每日檢查彈簧是否正常。大鼠在接受正畸治療前使用 Co60滅菌飼料，正畸治療后，每天被喂食三次動物奶粉。實驗過程中所有操作均得到北京大學實驗動物倫理委員會的批準，并且符合《實驗動物護理和使用指南》的要求。
圖1
牙齒-牙周韌帶-牙槽骨有限元模型的建立
a. 大鼠正畸實驗。b. 牙槽骨 micro-CT 照片，虛線表示牙槽骨受到的位移載荷。牙齒-牙周韌帶-牙槽骨有限元模型：c. 頰側視角；d. 遠中視角。牙槽骨的 von Mises 應變分布：e. 咬合力加載，f. 正畸力加載
Figure1. Finite element analysis on alveolar bonea. orthodontic experiment. b. micro-CT image of alveolar bone, the dash lines indicate the displacement transduced on alveolar bone. Finite element model of tooth-periodontal ligament-alveolar bone: c. buccal view; d. distal view. Strain distribution of alveolar bone: e. occlusal loading; f. orthodontic loading
1.2 牙齒-牙周韌帶-牙槽骨的有限元模型
將大鼠處死后取出其上頜骨，經 micro-CT 掃描后導入 Mimics 軟件對其右側第一磨牙處牙周組織進行三維重建（見圖 1b）。如圖 1c、d 所示，將重建后的牙齒-牙周韌帶-牙槽骨結構導入 ANSYS 中進行有限元計算。牙齒、牙周韌帶、牙槽骨均采用均勻各向同性線彈性材料，彈性模量分別為 20 000、0.7、2 000 MPa，泊松比分別為 0.3、0.45、0.3[7-8]。在牙槽骨的底部及遠中面施加固定約束，不同組織間設定為綁定接觸。模擬大鼠咬合時，對牙齒咬合面施加 4 N 咬合力[9]；模擬正畸實驗時，對牙齒遠中面施加大小為 0.8 N 的線載荷。圖 1e、f 分別為該模型在 4 N 咬合力或 0.8 N 正畸力加載下牙槽骨的應變分布。咬合、正畸加載下牙槽骨前 1% 應變的平均值分別為 650、167 με。此數據將作為理想模型咬合、正畸位移加載的幅值。
1.3 理想模型幾何模型的建立
松質骨是由大量骨小梁連接而成的多孔網狀結構，骨小梁的排列方向主要沿咬合方向。骨小梁的直徑在 40～120 μm 之間[10-11]，骨小梁間隙在 150～700 μm 之間[10, 12-14]。本研究以大鼠上頜骨右側第一磨牙近中牙根遠中附近松質骨為基礎，構建了由一定數目骨小梁組成的總體尺寸為 2 160 μm × 2 160 μm × 2 160 μm 的理想模型，見圖 2a。該理想模型包含 7 × 7 根沿 Z 方向（咬合方向）的骨小梁，4 × 7 根沿X、Y方向（非咬合方向）的骨小梁。骨小梁直徑為 80 μm，在X、Y方向上間隔為 280 μm、Z方向上間隔 560 μm。該理想模型位于充滿骨髓液的大小為 2 660 μm × 2 660 μm × 2 660 μm 的液體域中。
圖2
牙槽骨理想模型
a. 邊界條件；b. 咬合位移加載曲線；c. 正畸位移加載曲線
Figure2. Establishment of ideal model of alveolar bone and boundary conditiona. boundary condition; b. displacement loading along occlusal direction; c. displacement loading along orthodontic direction
1.4 理想模型邊界條件及材料參數
本文中松質骨與骨髓液的材料參數均參考前人的文獻所得。松質骨骨小梁采用均勻各向同性線彈性材料模型，彈性模量為 3 GPa，泊松比為 0.3[7, 15]。骨髓液的粘度系數會隨應變率的升高而降低，但有研究表明在低應變率下，使用 85 Pa?s 的常值粘度系數或實驗獲取的隨應變率變化的粘度系數進行計算，得到的流體剪應力分布無顯著差異。因此本研究將骨髓液假設為密度 1 g/cm3、粘度系數 85 Pa?s 的牛頓流體[6, 16]。本研究中骨髓液會在松質骨變形的驅動下產生流動，反之液體流動也會使骨發生相應變形，屬于雙向流固耦合問題，即在流場與固體結構中分別求解，并在每個時間步長內耦合迭代，直至流體與固體邊界處的解收斂后再繼續推進。理想模型的邊界條件如圖 2a 所示。咬合過程中，該處松質骨的上表面受到牙齒傳來的咬合力，產生沿咬合方向的正弦時變位移，在其底部施加固定約束。考慮正畸加載時，臨床上一般認為牙齒一側的牙槽骨受壓力，發生骨吸收，其對側受拉力，發生骨形成。但是已經有研究表明正畸力作用下牙槽骨不會出現明顯的拉力側或壓力側，牙齒遠中側的牙槽骨既會受到拉力也會受到壓力[17]，因此在固體域的一個側面施加拉伸或壓縮位移載荷。
牙槽骨在咬合力或正畸力作用下最大應變分別為 650 με 和 167 με，該應變值被假設為理想模型中骨小梁在咬合力或正畸力作用下的最大應變，將其乘以骨小梁長度即可求得最大變形量分別為 1.4 μm 和 0.36 μm，以此作為理想模型咬合、正畸位移加載的幅值。咬合位移加載曲線如圖 2b 所示，頻率為 2 Hz，頻率參考大鼠咬合的頻率[18-19]。正畸壓縮位移加載曲線如圖 2c 所示，選擇在 0.24 s 內將正畸位移加至 0.36 μm，之后維持恒定位移加載至 1 s。為避免載荷加載速度突變，在 0.24 s 附近對位移加載速度進行了平滑處理。考慮到被分割區域的松質骨應與其他骨小梁相連的實際情況，在固體域的其余側面施加彈性系數為 3 GPa 的彈性基礎約束，以模擬其周圍區域對中心區域的力學作用。流體域的六個面為允許流體自由進出的開邊界。該模型在 COMSOL 軟件中以四面體單元進行網格劃分，總計 1 412 222 單元。為避免邊界影響，分析數據時選取固體域內部的 2 × 5 × 5 根骨小梁作為研究對象。
1.5 模擬數據的統計分析
本文模擬計算結果的數據顯示為平均值 ± 標準差，并對不同組數據進行了單因素方差分析，P < 0.05 時認為差異有統計學意義。
2 結果
2.1 咬合加載對牙槽骨內不同方向骨小梁表面流體剪應力的影響
圖 3a 顯示了咬合位移加載下 1/4 周期時刻（t = 0.875 s）骨小梁表面流體剪應力分布，圖 3b 為相應時刻液體流速及骨小梁 von Mises 應力分布圖。從上述結果可以看出，沿咬合方向的骨小梁表面有較大的流體剪應力。在咬合加載下，牙槽骨內高水平應力集中在沿咬合方向的骨小梁及其與非咬合方向骨小梁的連接處；越靠近加載端液體流速越大。為了更加直觀地觀察松質骨內部的液體流動，圖 3c、d 展示了沿咬合方向 1/3 位置處截面內的流體剪應力及流速分布。孔隙內的固液交界面上具有較高的流體剪應力與流速。外層液體相較于內層具有更高的流體剪應力。
圖3
咬合加載下牙槽骨內液體流動（t = 0.875 s）
a. 骨小梁表面流體剪應力分布。b. 松質骨應力與流線圖。咬合方向（
a. distribution of fluid shear stress on trabecular surface. b. von Mises stress and streamline. Fluid flow within the section at 1/3 position along occlusal direction（axis
圖 4a 為骨小梁表面流體剪應力按照大小順序排列時四分位點及平均值隨時間的變化。此結果表明，骨小梁表面流體剪應力在每個加載周期的 1/4 時刻附近達到最大值，其均值最高為 2.4 Pa，應該注意的是此時間點處咬合位移加載速度最大。圖 4b、c 分別為咬合方向排列骨小梁（Z方向）及非咬合排列骨小梁（X、Y方向）表面流體剪應力不同分位點及平均值隨時間的變化。咬合方向排列骨小梁表面流體剪應力的 25、50、75 分位點及均值都高于非咬合方向排列的骨小梁。體外骨組織細胞培養實驗已經證實，當施加的流體剪應力在 1～4 Pa 時，成骨細胞將會產生力刺激響應[20]。因此，在圖 4d 中對咬合排列及非咬合排列骨小梁表面流體剪應力 1～4 Pa 的節點數占比進行了統計比較，可以看出沿咬合方向排列骨小梁表面約有 70% 的區域在超過一半周期的時間內受到利于成骨細胞響應的剪應力，而非咬合方向排列骨小梁只有 40% 左右。以上結果表明，咬合過程中沿咬合方向排列的骨小梁表面受到的流體刺激要顯著高于非咬合方向骨小梁，推測力學作用下此處成骨細胞會被激活，從而利于啟動骨形成過程。
圖4
咬合加載下骨小梁表面流體剪應力
骨小梁表面流體剪應力四分位點及平均值隨時間的變化：a. 所有骨小梁；b. 咬合方向排列骨小梁；c. 非咬合方向排列骨小梁。d. 流體剪應力 1～4 Pa 節點占總結點百分比
Figure4. Fluid shear stress on trabecular surface under occlusal loadingThe quartile and average values of fluid shear stress on trabecular surface
2.2 正畸加載對牙槽骨骨小梁表面流體剪應力的影響
圖 5a、c 分別為咬合與正畸壓縮、咬合與正畸拉伸加載下 1/4 周期時刻（t = 0.875 s）骨小梁表面流體剪應力分布，圖 5b、d 為相應時刻液體流速及骨小梁 von Mises 應力分布圖。由上述圖可以看出，正畸壓縮加載和正畸拉伸加載下骨小梁表面流速分布、流體剪應力分布在t = 0.875 s 與咬合加載下基本一致，但是沿正畸方向排列的骨小梁固體應力更大。正畸壓縮和正畸拉伸加載下，骨內液體的流動僅在壓縮或者拉伸恒定前有區別。如圖 5e、f 所示，加載過程中，液體主要沿動態加載的合力方向流動。圖 6a、b 為正畸壓縮、拉伸加載下骨小梁表面流體剪應力四分位點及平均值隨時間的變化。在t = 0.375 s，即正畸位移加載恒定后的第一個周期的幅值處，正畸壓縮加載下骨小梁表面流體剪應力的 25、50 分位點及平均值稍大于正畸拉伸加載下骨小梁，略小于咬合加載下骨小梁。但是在兩個周期后，拉伸或是壓縮狀態下的骨小梁表面流體剪應力都與咬合加載下骨小梁表面流體剪應力沒有差異。在圖 6c、d 中，對t = 0.375 s、t = 0.875 s 兩個時刻咬合、咬合與正畸壓縮、咬合與正畸拉伸加載條件下骨小梁表面的最大流體剪應力均值進行了比較。t = 0.375 s 時，咬合加載下骨小梁表面流體剪應力均值略大于正畸壓縮加載下骨小梁，而正畸拉伸加載下骨小梁表面流體剪應力最低。然而在t = 0.875 s 時，咬合與正畸共同作用下沿X、Y、Z方向排列的骨小梁表面流體剪應力都與無正畸作用時沒有顯著差異。該結果表明，正畸力施加在骨小梁上的瞬間會影響骨小梁表面的流體剪應力，但是在一段時間后，無論是拉伸還是壓縮，對骨小梁表面的流體剪應力都幾乎沒有影響。
圖5
正畸加載下牙槽骨內液體流動
Dstribution of fluid shear stress on trabecular surface at
圖6
不同加載條件下牙槽骨內液體流動
骨小梁表面流體剪應力四分位點及平均值隨時間的變化：a. 正畸壓縮加載；b. 正畸拉伸加載。骨小梁表面的最大流體剪應力均值比較：c.
Quartile and average values of fluid shear stress on trabecular surface
3 討論
松質骨是骨重建最為活躍的區域，目前越來越受到骨生物力學領域研究人員的關注，但較少有研究關注牙槽骨中的液體流動問題。牙槽骨的孔隙結構較為復雜，受限于計算量的限制，很難對整體牙槽骨真實模型進行計算。因此本研究使用多尺度分析的方法，首先建立牙齒-牙周韌帶-牙槽骨的有限元模型，得到咬合或正畸加載下牙槽骨的應變，將其作為理想模型的邊界條件。其次以牙槽骨孔隙結構的幾何參數為基礎，構建了由不同排列方向骨小梁組成的理想模型。該理想模型相較于真實模型更利于分析不同加載條件和骨小梁排列方向對其表面流體剪應力的影響。相較于使用流體動力學模型對骨內液體流動進行計算，本研究使用的流固耦合數值模擬方法能更為準確地反映外力作用下骨組織變形導致的液體流動。
大量研究已經表明，骨組織會優化其結構以適應所處的力學載荷環境。在缺乏力學刺激的情況下會發生骨礦丟失，而運動和力學刺激會促進骨內礦物含量增加。松質骨中骨小梁主要沿主應力方向生長，從而以最有效的結構形式承受外部載荷。流體剪應力作用下，骨組織細胞會發生鈣響應[21- 22]，從而參與調控骨形成[23-24]。除此之外，骨細胞受到液體刺激還會抑制破骨細胞的形成[25-26]。本研究結果顯示，咬合過程中沿咬合方向排列的骨小梁表面有更多的區域在超過一半的時間內受到 1～4 Pa 的流體剪應力，推測其表面的成骨細胞更容易發生生物學響應，從而促進咬合方向骨小梁的形成，這也在一定程度上解釋了牙槽骨內骨小梁主要沿咬合方向排列的原因。
本研究結果表明，正畸力施加在骨小梁上的瞬間會影響骨小梁表面的流體剪應力，這可能是由于壓縮和拉伸完全加載前，骨內液體流動的不同導致的。但是在正畸加載保持恒定一段時間后，無論是拉伸還是壓縮加載，對骨小梁表面的流體剪應力都沒有影響。骨小梁表面的流體剪應力只會受到動態加載的影響，而恒定載荷對其沒有影響。這說明正畸過程中牙槽骨發生的骨重建與其內部孔隙的液體流動無關。牙周韌帶中混雜著不同種類的細胞，如成纖維細胞、成骨細胞以及具有多種分化潛能的間充質干細胞等，而力學刺激會對間充質干細胞的分化產生影響[27-28]。除此之外成纖維細胞也可以感受和響應力學刺激，通過分泌骨保護素和破骨細胞分化因子等參與調控牙槽骨吸收，以影響破骨前體細胞的破骨向分化[29]。因此，正畸力作用下牙周韌帶的調控作用才是牙槽骨結構重建行為的主要原因。
本文以大鼠真實松質骨孔隙結構為基礎，構建了理想模型，應用流固耦合數值模擬方法，得到了動態咬合力加載下，無正畸加載、正畸壓縮加載和正畸拉伸加載三種情況下骨內液體的流動情況。模擬結果表明，沿咬合方向排列的骨小梁表面流體剪應力水平高于非咬合方向排列的骨小梁，正畸力對骨內液體的流動幾乎沒有影響。考慮到臨床上有時會通過調整牙齒咬合面形狀等方法改變咬合力的方向，本文結果表明這會在牙槽骨表面引起不同水平的流體剪應力，進而刺激骨組織表面的細胞產生響應，最終調控牙槽骨的結構重建。
利益沖突聲明：本文全體作者均聲明不存在利益沖突。
引言
牙槽骨包括固有牙槽骨、密質骨及松質骨等幾部分，其中松質骨位于密質骨和固有牙槽骨之間，由骨小梁和骨髓組成。牙槽骨的結構變化非常活躍，與牙齒發育、乳牙替換、恒牙移動等生理過程密切相關。這種結構的重建過程可由牙齒受到的載荷（如咬合力或正畸力）調控，主要表現為骨小梁的排列方向與咬合力相適應，以最有效的排列方式來抵抗外加載荷。實際上，骨的結構重建主要由成骨細胞主導的骨形成過程及破骨細胞主導的骨吸收過程組成。在咬合力或正畸力的作用下，骨基質的變形會引起附著于其上的細胞產生相應的應變；松質骨內孔隙體積的改變也會引起孔隙內液體的流動，從而在細胞表面產生流體剪應力[1]。許多研究已經證明，相對于基質的變形，流體刺激會在骨組織細胞中引起更為顯著的生物學響應[2-4]。因此，闡明牙槽骨內孔隙結構中的液體流動情況對于深入理解力學作用在牙槽骨內的傳導過程以及牙齒的發育形成、正畸牙移動等細胞水平的調控機制具有重要意義。
受限于在體測量實驗技術的缺乏，數值模擬是研究骨內液體流動空間分布和動態變化的有效方法。2012 年有學者構建了包含骨小梁孔隙結構的理想松質骨模型，使用流固耦合數值模擬方法計算了其中的液體流動、壓力和流體剪應力等參數的分布[5]。2015 年的另一項研究使用微計算機斷層掃描技術（micro computed tomography，micro-CT）掃描豬的股骨，通過三維重建得到其真實的松質骨幾何模型，也應用流固耦合數值模擬方法研究了骨小梁孔隙結構中的流體剪應力和壓力分布[6]。但目前還沒有關于牙槽骨內液體流動情況的數值模擬研究，因此咬合力或正畸力作用下牙槽骨內骨小梁表面流體剪應力的大小和分布尚不清楚，從而難于判斷牙槽骨中細胞所處的力學微環境。
本文針對大鼠的正畸牙移動模型，利用有限元方法分析了牙槽骨中的應變分布，進而構建了理想模型，應用流固耦合數值模擬方法，得到了動態咬合力加載下，無正畸力、正畸壓縮和正畸拉伸三種加載情況下骨內液體的流動情況，對比分析了咬合力方向、非咬合力方向骨小梁表面流體剪應力的差異，以及拉伸、壓縮對骨小梁表面流體剪應力的影響。
1 材料和方法
1.1 動物實驗
本實驗從北京大學醫學部實驗動物科學部購置一只 6～8 周齡的 SPF 級雌性 SD 大鼠［許可證號：SCXK（京）2016-0010］。將大鼠麻醉后置于手術臺上，以上頜 2 顆切牙作為支抗牙，使用鎳鈦拉簧對其上頜右側第一磨牙施加 0.8 N 正畸力 20 天（見圖 1a）。鎳鈦拉簧的一端通過直徑為 0.2 mm 的結扎絲固定于第一磨牙與第二磨牙之間，另一端連接測力計并測定正畸力，之后將其固定在兩顆切牙頸部 0.5～1.0 mm 的固位槽中。進而將光固化樹脂均勻覆蓋于切牙表面和固位溝處，并用 LED 光固化燈固化樹脂 40 s。將大鼠下頜切牙磨短并用樹脂覆蓋斷面處，以防止其磨斷口腔內的結扎絲。在正畸治療期間，每日檢查彈簧是否正常。大鼠在接受正畸治療前使用 Co60滅菌飼料，正畸治療后，每天被喂食三次動物奶粉。實驗過程中所有操作均得到北京大學實驗動物倫理委員會的批準，并且符合《實驗動物護理和使用指南》的要求。
圖1
牙齒-牙周韌帶-牙槽骨有限元模型的建立
a. 大鼠正畸實驗。b. 牙槽骨 micro-CT 照片，虛線表示牙槽骨受到的位移載荷。牙齒-牙周韌帶-牙槽骨有限元模型：c. 頰側視角；d. 遠中視角。牙槽骨的 von Mises 應變分布：e. 咬合力加載，f. 正畸力加載
Figure1. Finite element analysis on alveolar bonea. orthodontic experiment. b. micro-CT image of alveolar bone, the dash lines indicate the displacement transduced on alveolar bone. Finite element model of tooth-periodontal ligament-alveolar bone: c. buccal view; d. distal view. Strain distribution of alveolar bone: e. occlusal loading; f. orthodontic loading
1.2 牙齒-牙周韌帶-牙槽骨的有限元模型
將大鼠處死后取出其上頜骨，經 micro-CT 掃描后導入 Mimics 軟件對其右側第一磨牙處牙周組織進行三維重建（見圖 1b）。如圖 1c、d 所示，將重建后的牙齒-牙周韌帶-牙槽骨結構導入 ANSYS 中進行有限元計算。牙齒、牙周韌帶、牙槽骨均采用均勻各向同性線彈性材料，彈性模量分別為 20 000、0.7、2 000 MPa，泊松比分別為 0.3、0.45、0.3[7-8]。在牙槽骨的底部及遠中面施加固定約束，不同組織間設定為綁定接觸。模擬大鼠咬合時，對牙齒咬合面施加 4 N 咬合力[9]；模擬正畸實驗時，對牙齒遠中面施加大小為 0.8 N 的線載荷。圖 1e、f 分別為該模型在 4 N 咬合力或 0.8 N 正畸力加載下牙槽骨的應變分布。咬合、正畸加載下牙槽骨前 1% 應變的平均值分別為 650、167 με。此數據將作為理想模型咬合、正畸位移加載的幅值。
1.3 理想模型幾何模型的建立
松質骨是由大量骨小梁連接而成的多孔網狀結構，骨小梁的排列方向主要沿咬合方向。骨小梁的直徑在 40～120 μm 之間[10-11]，骨小梁間隙在 150～700 μm 之間[10, 12-14]。本研究以大鼠上頜骨右側第一磨牙近中牙根遠中附近松質骨為基礎，構建了由一定數目骨小梁組成的總體尺寸為 2 160 μm × 2 160 μm × 2 160 μm 的理想模型，見圖 2a。該理想模型包含 7 × 7 根沿 Z 方向（咬合方向）的骨小梁，4 × 7 根沿X、Y方向（非咬合方向）的骨小梁。骨小梁直徑為 80 μm，在X、Y方向上間隔為 280 μm、Z方向上間隔 560 μm。該理想模型位于充滿骨髓液的大小為 2 660 μm × 2 660 μm × 2 660 μm 的液體域中。
圖2
牙槽骨理想模型
a. 邊界條件；b. 咬合位移加載曲線；c. 正畸位移加載曲線
Figure2. Establishment of ideal model of alveolar bone and boundary conditiona. boundary condition; b. displacement loading along occlusal direction; c. displacement loading along orthodontic direction
1.4 理想模型邊界條件及材料參數
本文中松質骨與骨髓液的材料參數均參考前人的文獻所得。松質骨骨小梁采用均勻各向同性線彈性材料模型，彈性模量為 3 GPa，泊松比為 0.3[7, 15]。骨髓液的粘度系數會隨應變率的升高而降低，但有研究表明在低應變率下，使用 85 Pa?s 的常值粘度系數或實驗獲取的隨應變率變化的粘度系數進行計算，得到的流體剪應力分布無顯著差異。因此本研究將骨髓液假設為密度 1 g/cm3、粘度系數 85 Pa?s 的牛頓流體[6, 16]。本研究中骨髓液會在松質骨變形的驅動下產生流動，反之液體流動也會使骨發生相應變形，屬于雙向流固耦合問題，即在流場與固體結構中分別求解，并在每個時間步長內耦合迭代，直至流體與固體邊界處的解收斂后再繼續推進。理想模型的邊界條件如圖 2a 所示。咬合過程中，該處松質骨的上表面受到牙齒傳來的咬合力，產生沿咬合方向的正弦時變位移，在其底部施加固定約束。考慮正畸加載時，臨床上一般認為牙齒一側的牙槽骨受壓力，發生骨吸收，其對側受拉力，發生骨形成。但是已經有研究表明正畸力作用下牙槽骨不會出現明顯的拉力側或壓力側，牙齒遠中側的牙槽骨既會受到拉力也會受到壓力[17]，因此在固體域的一個側面施加拉伸或壓縮位移載荷。
牙槽骨在咬合力或正畸力作用下最大應變分別為 650 με 和 167 με，該應變值被假設為理想模型中骨小梁在咬合力或正畸力作用下的最大應變，將其乘以骨小梁長度即可求得最大變形量分別為 1.4 μm 和 0.36 μm，以此作為理想模型咬合、正畸位移加載的幅值。咬合位移加載曲線如圖 2b 所示，頻率為 2 Hz，頻率參考大鼠咬合的頻率[18-19]。正畸壓縮位移加載曲線如圖 2c 所示，選擇在 0.24 s 內將正畸位移加至 0.36 μm，之后維持恒定位移加載至 1 s。為避免載荷加載速度突變，在 0.24 s 附近對位移加載速度進行了平滑處理。考慮到被分割區域的松質骨應與其他骨小梁相連的實際情況，在固體域的其余側面施加彈性系數為 3 GPa 的彈性基礎約束，以模擬其周圍區域對中心區域的力學作用。流體域的六個面為允許流體自由進出的開邊界。該模型在 COMSOL 軟件中以四面體單元進行網格劃分，總計 1 412 222 單元。為避免邊界影響，分析數據時選取固體域內部的 2 × 5 × 5 根骨小梁作為研究對象。
1.5 模擬數據的統計分析
本文模擬計算結果的數據顯示為平均值 ± 標準差，并對不同組數據進行了單因素方差分析，P < 0.05 時認為差異有統計學意義。
2 結果
2.1 咬合加載對牙槽骨內不同方向骨小梁表面流體剪應力的影響
圖 3a 顯示了咬合位移加載下 1/4 周期時刻（t = 0.875 s）骨小梁表面流體剪應力分布，圖 3b 為相應時刻液體流速及骨小梁 von Mises 應力分布圖。從上述結果可以看出，沿咬合方向的骨小梁表面有較大的流體剪應力。在咬合加載下，牙槽骨內高水平應力集中在沿咬合方向的骨小梁及其與非咬合方向骨小梁的連接處；越靠近加載端液體流速越大。為了更加直觀地觀察松質骨內部的液體流動，圖 3c、d 展示了沿咬合方向 1/3 位置處截面內的流體剪應力及流速分布。孔隙內的固液交界面上具有較高的流體剪應力與流速。外層液體相較于內層具有更高的流體剪應力。
圖3
咬合加載下牙槽骨內液體流動（t = 0.875 s）
a. 骨小梁表面流體剪應力分布。b. 松質骨應力與流線圖。咬合方向（
a. distribution of fluid shear stress on trabecular surface. b. von Mises stress and streamline. Fluid flow within the section at 1/3 position along occlusal direction（axis
圖 4a 為骨小梁表面流體剪應力按照大小順序排列時四分位點及平均值隨時間的變化。此結果表明，骨小梁表面流體剪應力在每個加載周期的 1/4 時刻附近達到最大值，其均值最高為 2.4 Pa，應該注意的是此時間點處咬合位移加載速度最大。圖 4b、c 分別為咬合方向排列骨小梁（Z方向）及非咬合排列骨小梁（X、Y方向）表面流體剪應力不同分位點及平均值隨時間的變化。咬合方向排列骨小梁表面流體剪應力的 25、50、75 分位點及均值都高于非咬合方向排列的骨小梁。體外骨組織細胞培養實驗已經證實，當施加的流體剪應力在 1～4 Pa 時，成骨細胞將會產生力刺激響應[20]。因此，在圖 4d 中對咬合排列及非咬合排列骨小梁表面流體剪應力 1～4 Pa 的節點數占比進行了統計比較，可以看出沿咬合方向排列骨小梁表面約有 70% 的區域在超過一半周期的時間內受到利于成骨細胞響應的剪應力，而非咬合方向排列骨小梁只有 40% 左右。以上結果表明，咬合過程中沿咬合方向排列的骨小梁表面受到的流體刺激要顯著高于非咬合方向骨小梁，推測力學作用下此處成骨細胞會被激活，從而利于啟動骨形成過程。
圖4
咬合加載下骨小梁表面流體剪應力
骨小梁表面流體剪應力四分位點及平均值隨時間的變化：a. 所有骨小梁；b. 咬合方向排列骨小梁；c. 非咬合方向排列骨小梁。d. 流體剪應力 1～4 Pa 節點占總結點百分比
Figure4. Fluid shear stress on trabecular surface under occlusal loadingThe quartile and average values of fluid shear stress on trabecular surface
2.2 正畸加載對牙槽骨骨小梁表面流體剪應力的影響
圖 5a、c 分別為咬合與正畸壓縮、咬合與正畸拉伸加載下 1/4 周期時刻（t = 0.875 s）骨小梁表面流體剪應力分布，圖 5b、d 為相應時刻液體流速及骨小梁 von Mises 應力分布圖。由上述圖可以看出，正畸壓縮加載和正畸拉伸加載下骨小梁表面流速分布、流體剪應力分布在t = 0.875 s 與咬合加載下基本一致，但是沿正畸方向排列的骨小梁固體應力更大。正畸壓縮和正畸拉伸加載下，骨內液體的流動僅在壓縮或者拉伸恒定前有區別。如圖 5e、f 所示，加載過程中，液體主要沿動態加載的合力方向流動。圖 6a、b 為正畸壓縮、拉伸加載下骨小梁表面流體剪應力四分位點及平均值隨時間的變化。在t = 0.375 s，即正畸位移加載恒定后的第一個周期的幅值處，正畸壓縮加載下骨小梁表面流體剪應力的 25、50 分位點及平均值稍大于正畸拉伸加載下骨小梁，略小于咬合加載下骨小梁。但是在兩個周期后，拉伸或是壓縮狀態下的骨小梁表面流體剪應力都與咬合加載下骨小梁表面流體剪應力沒有差異。在圖 6c、d 中，對t = 0.375 s、t = 0.875 s 兩個時刻咬合、咬合與正畸壓縮、咬合與正畸拉伸加載條件下骨小梁表面的最大流體剪應力均值進行了比較。t = 0.375 s 時，咬合加載下骨小梁表面流體剪應力均值略大于正畸壓縮加載下骨小梁，而正畸拉伸加載下骨小梁表面流體剪應力最低。然而在t = 0.875 s 時，咬合與正畸共同作用下沿X、Y、Z方向排列的骨小梁表面流體剪應力都與無正畸作用時沒有顯著差異。該結果表明，正畸力施加在骨小梁上的瞬間會影響骨小梁表面的流體剪應力，但是在一段時間后，無論是拉伸還是壓縮，對骨小梁表面的流體剪應力都幾乎沒有影響。
圖5
正畸加載下牙槽骨內液體流動
Dstribution of fluid shear stress on trabecular surface at
圖6
不同加載條件下牙槽骨內液體流動
骨小梁表面流體剪應力四分位點及平均值隨時間的變化：a. 正畸壓縮加載；b. 正畸拉伸加載。骨小梁表面的最大流體剪應力均值比較：c.
Quartile and average values of fluid shear stress on trabecular surface
3 討論
松質骨是骨重建最為活躍的區域，目前越來越受到骨生物力學領域研究人員的關注，但較少有研究關注牙槽骨中的液體流動問題。牙槽骨的孔隙結構較為復雜，受限于計算量的限制，很難對整體牙槽骨真實模型進行計算。因此本研究使用多尺度分析的方法，首先建立牙齒-牙周韌帶-牙槽骨的有限元模型，得到咬合或正畸加載下牙槽骨的應變，將其作為理想模型的邊界條件。其次以牙槽骨孔隙結構的幾何參數為基礎，構建了由不同排列方向骨小梁組成的理想模型。該理想模型相較于真實模型更利于分析不同加載條件和骨小梁排列方向對其表面流體剪應力的影響。相較于使用流體動力學模型對骨內液體流動進行計算，本研究使用的流固耦合數值模擬方法能更為準確地反映外力作用下骨組織變形導致的液體流動。
大量研究已經表明，骨組織會優化其結構以適應所處的力學載荷環境。在缺乏力學刺激的情況下會發生骨礦丟失，而運動和力學刺激會促進骨內礦物含量增加。松質骨中骨小梁主要沿主應力方向生長，從而以最有效的結構形式承受外部載荷。流體剪應力作用下，骨組織細胞會發生鈣響應[21- 22]，從而參與調控骨形成[23-24]。除此之外，骨細胞受到液體刺激還會抑制破骨細胞的形成[25-26]。本研究結果顯示，咬合過程中沿咬合方向排列的骨小梁表面有更多的區域在超過一半的時間內受到 1～4 Pa 的流體剪應力，推測其表面的成骨細胞更容易發生生物學響應，從而促進咬合方向骨小梁的形成，這也在一定程度上解釋了牙槽骨內骨小梁主要沿咬合方向排列的原因。
本研究結果表明，正畸力施加在骨小梁上的瞬間會影響骨小梁表面的流體剪應力，這可能是由于壓縮和拉伸完全加載前，骨內液體流動的不同導致的。但是在正畸加載保持恒定一段時間后，無論是拉伸還是壓縮加載，對骨小梁表面的流體剪應力都沒有影響。骨小梁表面的流體剪應力只會受到動態加載的影響，而恒定載荷對其沒有影響。這說明正畸過程中牙槽骨發生的骨重建與其內部孔隙的液體流動無關。牙周韌帶中混雜著不同種類的細胞，如成纖維細胞、成骨細胞以及具有多種分化潛能的間充質干細胞等，而力學刺激會對間充質干細胞的分化產生影響[27-28]。除此之外成纖維細胞也可以感受和響應力學刺激，通過分泌骨保護素和破骨細胞分化因子等參與調控牙槽骨吸收，以影響破骨前體細胞的破骨向分化[29]。因此，正畸力作用下牙周韌帶的調控作用才是牙槽骨結構重建行為的主要原因。
本文以大鼠真實松質骨孔隙結構為基礎，構建了理想模型，應用流固耦合數值模擬方法，得到了動態咬合力加載下，無正畸加載、正畸壓縮加載和正畸拉伸加載三種情況下骨內液體的流動情況。模擬結果表明，沿咬合方向排列的骨小梁表面流體剪應力水平高于非咬合方向排列的骨小梁，正畸力對骨內液體的流動幾乎沒有影響。考慮到臨床上有時會通過調整牙齒咬合面形狀等方法改變咬合力的方向，本文結果表明這會在牙槽骨表面引起不同水平的流體剪應力，進而刺激骨組織表面的細胞產生響應，最終調控牙槽骨的結構重建。
利益沖突聲明：本文全體作者均聲明不存在利益沖突。
