針對骨組織工程支架降解問題,研究利用支架的非均勻結構來控制骨缺損修復效果的可行性。在分析了組織工程骨骼修復過程及主要影響因素的基礎上,構建了骨成形控制方程、支架降解控制方程為主體的骨缺損修復預測模型,并與有限元法相結合,預測骨缺損修復過程和修復后的骨結構形態。基于以上方法,分別模擬了均勻及非均勻支架情形下的骨缺損修復效果。預測結果表明,支架結構對修復效果具有較大影響,可以利用支架的非均勻結構來控制骨缺損修復效果。
引用本文： 周根, 劉云峰, 姜獻峰. 可降解支架修復骨缺損的有限元模擬預測. 生物醫學工程學雜志, 2014, 31(3): 601-605. doi: 10.7507/1001-5515.20140113 復制
引言
骨骼作為人體的支架，擔負著保持、承重、造血、貯鈣、代謝等功能[1]。然而，創傷、腫瘤切除和畸形矯正造成的骨缺損在臨床中又常常出現。小量骨的缺損可以通過人體生物組織的再生進行修復，但大段骨缺損無法通過骨骼自行生長來彌補[2]。對此，除了自體移植、異體移植和人工替代材料等傳統骨缺損修復方法外，理想的辦法是應用骨組織工程方法重建其力學系統和生物學系統[3]。
在骨組織工程用于骨骼缺損的修復中，生物可降解骨組織工程支架起著重要作用，它是由多孔微結構組成的復雜框架，為細胞提供結構支撐，引導組織再生，控制組織結構，代替缺損的骨組織提供力學性能[3]。組織工程可降解支架可以創造一種微環境，以利于細胞的黏附、增殖，并且誘導特定的細胞功能，引導和調節細胞間的相互作用，并能夠最終被新生組織替代[4]。
組織工程骨骼修復過程就是支架逐漸降解，同時新骨不斷生成的過程[5]。這是既獨立又相互關聯的兩個過程，且與所處應力環境密切相關。Hollister等[6]通過培養實驗和有限元模擬證實了應力刺激對骨生長的影響。徐健斌等[7]在前人研究的基礎上，提出了應力環境下骨支架降解的模型。Sanz-Herrera等[8-9]在不考慮支架降解情形下，建立了特定類型支架內部骨組織再生數學模型，并針對支架微結構決定內部骨定向生長做了數值研究。Sung等[10]對支架降解率在細胞生長和血管生成方面的作用做了實驗研究。但是這些研究主要以支架的規則幾何結構為基礎，也缺乏對各關鍵因素及相互作用的綜合考慮。
為了定量預測骨缺損修復過程和修復后骨的結構形態，整個過程中與時間空間相關的量要能方便地確定并記錄下來，有限元法就是完成此任務的有效工具。通過支架降解骨生長模型和有限元法相結合，預測骨缺損修復過程和修復后骨的結構形態。
1 方法
骨組織有著復雜的結構和生理機制，骨修復過程實際上是一個相當復雜的生理過程，受到生物學因子、力學因素等多重調控和影響。根據組織工程骨缺損修復過程的特點，重點考慮支架降解、新骨成型、力學環境等關鍵因素，建立相關的控制方程，進而建立完整的有限元模型進行模擬預測。
1.1 骨成形控制方程
盡管支架內的骨生長情形與正常的骨改建有所區別，但同是作為骨細胞的生理行為，兩者應該具有某些共性，因此可以從骨功能適應性理論的基本方程出發[11]，推導適用于支架內部的骨成形控制方程。假設基本的骨成形控制方程滿足如下形式，即
| $\partial \rho \left( \Omega ,t \right)/\partial t=B\left( t \right)\left( S\left( \Omega ,t \right)-K\left( t \right) \right)~,$ |
式中ρ為表征骨生長狀態的特征量，一般取為骨的表觀密度，Ω為空間位置，t為時間，B為靈敏度系數，S為力學激勵，K為生長閾值。
趙文志等[12]結合動物實驗和參數反演研究了大鼠股骨應力與生長生物模型，繪制了B和K關于時間t的曲線，并給出了表達式為：
| $B=1\times {{10}^{-6}}{{e}^{-0.3816t}}$ |
| $K=2294.4{{t}^{0.6216}}$ |
推廣到一般情形下，可假設B和K與t存在如下關系，即
| $B={{b}_{0}}e-\alpha t,$ |
| $K={{k}_{0}}t\beta ,$ |
式中b0、k0、α、β是生物因素相關的系數。
考慮到生物體內生長的“惰性區”現象[13]，得出綜合的骨成形控制方程：
| $\partial \rho \left( \Omega ,t \right)/\partial t=\left\{ \begin{align} & {{b}_{0}}e-\alpha t[S\left( \Omega ,t \right)-(1+\delta ){{k}_{0}}{{t}^{\beta }}];S>(1+\delta ){{k}_{0}}{{t}^{\beta }} \\ & 0;~else \\ & {{b}_{0}}{{e}^{-\alpha t}}[S\left( \Omega ,t \right)-(1-\delta ){{k}_{0}}{{t}^{\beta }}];\text{ }S ＜(1-\delta ){{k}_{0}}t\beta \\ \end{align} \right.,$ |
式中δ為表示“惰性區”范圍的系數。
骨組織在支架內微結構的引導下生長鈣化，初期形成的結構與支架類似，根據微結構幾何特征可以建立?ρ/?t與?h/?t之間的關系，即
| $\partial h\left( \Omega ,t \right)/\partial t=\frac{{{C}_{1}}{{L}^{3}}}{A{{\rho }_{c}}}\partial \rho \left( \Omega ,t \right)/\partial t,$ |
式中h為骨微結構的特征壁厚，C1為幾何相關的系數，L
在骨-支架系統的局部微結構孔隙內，可以簡單地認為骨分布受如下約束：
| $L-2\left( h+H \right)\ge \Phi ~,$ |
式中H為支架微結構特征壁厚，Φ為孔隙容忍度系數。
1.2 支架降解控制方程
可降解支架植入體內后，隨著時間的推移，各部位會發生不同程度的降解，從應力環境、幾何微結構等方面影響骨成型，反之又受骨成型的影響，因此有必要建立一個能反映支架降解主要影響因子的控制方程。
各國研究人員對可降解支架的降解在理論和實驗上都做了很多的研究。可降解支架在體內的降解機制主要可以分為：化學降解、生物降解和物理降解。前期的降解模型主要考慮的是化學降解和生物降解。徐健斌等[7]在前人研究的基礎上，提出了應力環境下骨支架降解的模型。本文對其模型加以修正，并引入周圍骨細胞濃度對支架降解的影響系數η，得到如下支架降解控制方程：
| $\partial V\left( \Omega ,t \right)/\partial t=\eta {{k}_{s}}A\lambda \left( a-\frac{b}{1+exp[(\varepsilon \left( \Omega ,t \right)-c)/d]} \right)\text{ }/{{\rho }_{s}},$ |
式中V為溶解體積，ks為分解率常量，ε為應變，ρs為材料密度，λ、a、b、c、d為與可溶解材料特性相關的常數。
與骨成型類似，可以建立?V/?t與?H/?t之間的關系，根據微結構幾何特征推導得：
| $\partial H/\partial t=\frac{{{C}_{2}}}{A}\partial V/\partial t,$ |
式中C2為幾何相關系數。
設正常體液中的骨細胞濃度系數為η*，則骨-支架系統微結構中的骨細胞濃度系數為D3L3η*，其中D表示孔隙直徑，考慮到孔隙大小的權重影響，周圍骨細胞濃度對支架降解的影響系數η可以由下式估算：
| $\eta =\sum\limits_{i=1}^{n}{\left| \frac{{{D}_{i}}^{3}}{{{L}^{3}}}{{\eta }^{*}}\frac{{{D}_{i}}^{3}}{\sum\limits_{j=1}^{n}{{{D}_{j}}^{3}}} \right|}$ |
1.3 等效材料特性估算
隨著支架在體內慢慢降解吸收，以及新生骨組織逐漸發育，整個骨-支架系統的結構和力學性能會發生連續性變化，通過計算機直接仿真其結構的實時變化比較復雜，計算成本也相當高，因此研究人員借助于各種理論來間接估計骨-支架系統力學性能的動態變化，如多孔材料特性、均勻化理論、有限元方法等[14]，其中利用多孔固體的力學性能直接與孔穴形狀和結構相關這個特性，通過相關理論建立起等效力學性能與結構之間的關系[15]，是一個比較簡單可靠的手段。一般地認為存在如下關系，即
| $\frac{{{E}^{*}}}{{{E}_{s}}}={{c}^{*}}{{\left( \frac{{{\rho }^{*}}}{{{\rho }_{s}}} \right)}^{r}}$ |
式中E*為多孔固體的等效材料常數，Es為實體材料常數，ρ*為多孔材料密度，c*、r為幾何結構相關系數。
為了使估算更加精確，借助prony級數擴展上述表達式：
| $\frac{{{E}^{*}}}{{{E}_{s}}}=\sum\limits_{i=1}^{n}{\left( {{c}_{i}}^{*}(\frac{{{\rho }^{*}}}{{{\rho }_{s}}}){{r}_{i}} \right)}$ |
1.4 有限元模擬
完整的有限元模型包括前處理、模擬計算、后處理三個步驟，先進行幾何模型、材料定義、載荷和邊界條件、網格劃分等前處理，這些可根據具體的應用情形而設置，再嵌入相關的控制方程，就可進行有限元迭代分析并輸出結果，模擬出利用可降解支架修復骨缺損的過程，并預測修復效果。具體的流程如圖 1所示。
圖1
有限元模擬流程圖
Figure1.
Flow chart of finite element simulation
2 實例
參照文獻[16]的邊界載荷條件，利用上述方法對椎骨的缺損修復做模擬預測，主要參數如表 1所示，預測過程及結果如圖 2所示。
表1
主要參數取值
Table1.
Values of main parameters
3 討論
多孔支架內骨組織生長是一個復雜的生理過程，受生物條件、力學刺激等多重因素的調控，本文的骨生長模型以生物條件良好為前提，因此以力學刺激為主導因素，借鑒骨改建基本方程，使其中的常系數改變為時間的函數，推導出骨成形控制方程，并考慮了支架結構對其的約束。支架降解模型則主要考慮為應力環境下受周圍細胞濃度影響的簡單水解。兩者的結合便構成了利用可降解支架修復骨缺損的預測模型。
實例中對椎骨密度的預測結果(見圖 2a)表明骨成形控制方程的有效性，接著模擬在正常椎骨一側出現骨缺損，并分別植入與缺損外形匹配的可降解均勻支架(見圖 2b)和非均勻支架(見圖 2c)，最后預測了支架完全降解后，兩種情形下缺損處的骨生長情況(見圖 2(d、e))。從預測結果來看，缺損處修復后的骨密度分布與原來差別較大，其中一個很重要的原因是植入支架后引起了骨骼內部受力狀況的改變，從而影響了骨生長；再者就是支架微結構對新骨具有約束作用，從而影響了骨成形。相比而言，支架結構非均勻情形下比均勻情形下要好一些，這說明支架結構對修復效果具有較大影響，可以通過調整支架微孔的非均勻分布來獲得較理想的修復效果，使修復后的骨密度更接近正常骨密度分布。
圖2
椎骨缺損修復模擬預測
椎骨密度預測結果如(a)所示；缺損處分別植入可降解均勻(b)和非均勻(c)支架；均勻(d)和非均勻(e)支架的骨缺損修復預測結果
Figure2. Simulation prediction of vertebra defect repairpredicted outcome of vertebra density was showed as (a); implanting biodegradable homogeneous (b) and inhomogeneous (c) scaffold respectively for defect; predicting outcome of bone defect repair with homogeneous (d) and inhomogeneous (e) scaffold
4 結論
通過構建骨成形和支架降解控制方程為主體的骨缺損修復預測模型，并借助有限元的手段實現仿真，趨勢性地預測了骨缺損修復過程和修復后的骨結構形態。通過對預測結果的比較分析，可得出以下結論：
(1)支架降解骨生長模型和有限元法相結合是預測骨缺損修復過程及效果的有效方法。
(2)支架結構對骨缺損修復效果具有較大影響。
(3)可以利用支架的非均勻結構來控制骨缺損修復效果。
本文的模型還略顯粗糙，在整個模擬過程中有很多影響結果的參數，這些參數的力學和生理意義將進一步研究，以期獲得更加符合實際的結果。
引言
骨骼作為人體的支架，擔負著保持、承重、造血、貯鈣、代謝等功能[1]。然而，創傷、腫瘤切除和畸形矯正造成的骨缺損在臨床中又常常出現。小量骨的缺損可以通過人體生物組織的再生進行修復，但大段骨缺損無法通過骨骼自行生長來彌補[2]。對此，除了自體移植、異體移植和人工替代材料等傳統骨缺損修復方法外，理想的辦法是應用骨組織工程方法重建其力學系統和生物學系統[3]。
在骨組織工程用于骨骼缺損的修復中，生物可降解骨組織工程支架起著重要作用，它是由多孔微結構組成的復雜框架，為細胞提供結構支撐，引導組織再生，控制組織結構，代替缺損的骨組織提供力學性能[3]。組織工程可降解支架可以創造一種微環境，以利于細胞的黏附、增殖，并且誘導特定的細胞功能，引導和調節細胞間的相互作用，并能夠最終被新生組織替代[4]。
組織工程骨骼修復過程就是支架逐漸降解，同時新骨不斷生成的過程[5]。這是既獨立又相互關聯的兩個過程，且與所處應力環境密切相關。Hollister等[6]通過培養實驗和有限元模擬證實了應力刺激對骨生長的影響。徐健斌等[7]在前人研究的基礎上，提出了應力環境下骨支架降解的模型。Sanz-Herrera等[8-9]在不考慮支架降解情形下，建立了特定類型支架內部骨組織再生數學模型，并針對支架微結構決定內部骨定向生長做了數值研究。Sung等[10]對支架降解率在細胞生長和血管生成方面的作用做了實驗研究。但是這些研究主要以支架的規則幾何結構為基礎，也缺乏對各關鍵因素及相互作用的綜合考慮。
為了定量預測骨缺損修復過程和修復后骨的結構形態，整個過程中與時間空間相關的量要能方便地確定并記錄下來，有限元法就是完成此任務的有效工具。通過支架降解骨生長模型和有限元法相結合，預測骨缺損修復過程和修復后骨的結構形態。
1 方法
骨組織有著復雜的結構和生理機制，骨修復過程實際上是一個相當復雜的生理過程，受到生物學因子、力學因素等多重調控和影響。根據組織工程骨缺損修復過程的特點，重點考慮支架降解、新骨成型、力學環境等關鍵因素，建立相關的控制方程，進而建立完整的有限元模型進行模擬預測。
1.1 骨成形控制方程
盡管支架內的骨生長情形與正常的骨改建有所區別，但同是作為骨細胞的生理行為，兩者應該具有某些共性，因此可以從骨功能適應性理論的基本方程出發[11]，推導適用于支架內部的骨成形控制方程。假設基本的骨成形控制方程滿足如下形式，即
| $\partial \rho \left( \Omega ,t \right)/\partial t=B\left( t \right)\left( S\left( \Omega ,t \right)-K\left( t \right) \right)~,$ |
式中ρ為表征骨生長狀態的特征量，一般取為骨的表觀密度，Ω為空間位置，t為時間，B為靈敏度系數，S為力學激勵，K為生長閾值。
趙文志等[12]結合動物實驗和參數反演研究了大鼠股骨應力與生長生物模型，繪制了B和K關于時間t的曲線，并給出了表達式為：
| $B=1\times {{10}^{-6}}{{e}^{-0.3816t}}$ |
| $K=2294.4{{t}^{0.6216}}$ |
推廣到一般情形下，可假設B和K與t存在如下關系，即
| $B={{b}_{0}}e-\alpha t,$ |
| $K={{k}_{0}}t\beta ,$ |
式中b0、k0、α、β是生物因素相關的系數。
考慮到生物體內生長的“惰性區”現象[13]，得出綜合的骨成形控制方程：
| $\partial \rho \left( \Omega ,t \right)/\partial t=\left\{ \begin{align} & {{b}_{0}}e-\alpha t[S\left( \Omega ,t \right)-(1+\delta ){{k}_{0}}{{t}^{\beta }}];S>(1+\delta ){{k}_{0}}{{t}^{\beta }} \\ & 0;~else \\ & {{b}_{0}}{{e}^{-\alpha t}}[S\left( \Omega ,t \right)-(1-\delta ){{k}_{0}}{{t}^{\beta }}];\text{ }S ＜(1-\delta ){{k}_{0}}t\beta \\ \end{align} \right.,$ |
式中δ為表示“惰性區”范圍的系數。
骨組織在支架內微結構的引導下生長鈣化，初期形成的結構與支架類似，根據微結構幾何特征可以建立?ρ/?t與?h/?t之間的關系，即
| $\partial h\left( \Omega ,t \right)/\partial t=\frac{{{C}_{1}}{{L}^{3}}}{A{{\rho }_{c}}}\partial \rho \left( \Omega ,t \right)/\partial t,$ |
式中h為骨微結構的特征壁厚，C1為幾何相關的系數，L
在骨-支架系統的局部微結構孔隙內，可以簡單地認為骨分布受如下約束：
| $L-2\left( h+H \right)\ge \Phi ~,$ |
式中H為支架微結構特征壁厚，Φ為孔隙容忍度系數。
1.2 支架降解控制方程
可降解支架植入體內后，隨著時間的推移，各部位會發生不同程度的降解，從應力環境、幾何微結構等方面影響骨成型，反之又受骨成型的影響，因此有必要建立一個能反映支架降解主要影響因子的控制方程。
各國研究人員對可降解支架的降解在理論和實驗上都做了很多的研究。可降解支架在體內的降解機制主要可以分為：化學降解、生物降解和物理降解。前期的降解模型主要考慮的是化學降解和生物降解。徐健斌等[7]在前人研究的基礎上，提出了應力環境下骨支架降解的模型。本文對其模型加以修正，并引入周圍骨細胞濃度對支架降解的影響系數η，得到如下支架降解控制方程：
| $\partial V\left( \Omega ,t \right)/\partial t=\eta {{k}_{s}}A\lambda \left( a-\frac{b}{1+exp[(\varepsilon \left( \Omega ,t \right)-c)/d]} \right)\text{ }/{{\rho }_{s}},$ |
式中V為溶解體積，ks為分解率常量，ε為應變，ρs為材料密度，λ、a、b、c、d為與可溶解材料特性相關的常數。
與骨成型類似，可以建立?V/?t與?H/?t之間的關系，根據微結構幾何特征推導得：
| $\partial H/\partial t=\frac{{{C}_{2}}}{A}\partial V/\partial t,$ |
式中C2為幾何相關系數。
設正常體液中的骨細胞濃度系數為η*，則骨-支架系統微結構中的骨細胞濃度系數為D3L3η*，其中D表示孔隙直徑，考慮到孔隙大小的權重影響，周圍骨細胞濃度對支架降解的影響系數η可以由下式估算：
| $\eta =\sum\limits_{i=1}^{n}{\left| \frac{{{D}_{i}}^{3}}{{{L}^{3}}}{{\eta }^{*}}\frac{{{D}_{i}}^{3}}{\sum\limits_{j=1}^{n}{{{D}_{j}}^{3}}} \right|}$ |
1.3 等效材料特性估算
隨著支架在體內慢慢降解吸收，以及新生骨組織逐漸發育，整個骨-支架系統的結構和力學性能會發生連續性變化，通過計算機直接仿真其結構的實時變化比較復雜，計算成本也相當高，因此研究人員借助于各種理論來間接估計骨-支架系統力學性能的動態變化，如多孔材料特性、均勻化理論、有限元方法等[14]，其中利用多孔固體的力學性能直接與孔穴形狀和結構相關這個特性，通過相關理論建立起等效力學性能與結構之間的關系[15]，是一個比較簡單可靠的手段。一般地認為存在如下關系，即
| $\frac{{{E}^{*}}}{{{E}_{s}}}={{c}^{*}}{{\left( \frac{{{\rho }^{*}}}{{{\rho }_{s}}} \right)}^{r}}$ |
式中E*為多孔固體的等效材料常數，Es為實體材料常數，ρ*為多孔材料密度，c*、r為幾何結構相關系數。
為了使估算更加精確，借助prony級數擴展上述表達式：
| $\frac{{{E}^{*}}}{{{E}_{s}}}=\sum\limits_{i=1}^{n}{\left( {{c}_{i}}^{*}(\frac{{{\rho }^{*}}}{{{\rho }_{s}}}){{r}_{i}} \right)}$ |
1.4 有限元模擬
完整的有限元模型包括前處理、模擬計算、后處理三個步驟，先進行幾何模型、材料定義、載荷和邊界條件、網格劃分等前處理，這些可根據具體的應用情形而設置，再嵌入相關的控制方程，就可進行有限元迭代分析并輸出結果，模擬出利用可降解支架修復骨缺損的過程，并預測修復效果。具體的流程如圖 1所示。
圖1
有限元模擬流程圖
Figure1.
Flow chart of finite element simulation
2 實例
參照文獻[16]的邊界載荷條件，利用上述方法對椎骨的缺損修復做模擬預測，主要參數如表 1所示，預測過程及結果如圖 2所示。
表1
主要參數取值
Table1.
Values of main parameters
3 討論
多孔支架內骨組織生長是一個復雜的生理過程，受生物條件、力學刺激等多重因素的調控，本文的骨生長模型以生物條件良好為前提，因此以力學刺激為主導因素，借鑒骨改建基本方程，使其中的常系數改變為時間的函數，推導出骨成形控制方程，并考慮了支架結構對其的約束。支架降解模型則主要考慮為應力環境下受周圍細胞濃度影響的簡單水解。兩者的結合便構成了利用可降解支架修復骨缺損的預測模型。
實例中對椎骨密度的預測結果(見圖 2a)表明骨成形控制方程的有效性，接著模擬在正常椎骨一側出現骨缺損，并分別植入與缺損外形匹配的可降解均勻支架(見圖 2b)和非均勻支架(見圖 2c)，最后預測了支架完全降解后，兩種情形下缺損處的骨生長情況(見圖 2(d、e))。從預測結果來看，缺損處修復后的骨密度分布與原來差別較大，其中一個很重要的原因是植入支架后引起了骨骼內部受力狀況的改變，從而影響了骨生長；再者就是支架微結構對新骨具有約束作用，從而影響了骨成形。相比而言，支架結構非均勻情形下比均勻情形下要好一些，這說明支架結構對修復效果具有較大影響，可以通過調整支架微孔的非均勻分布來獲得較理想的修復效果，使修復后的骨密度更接近正常骨密度分布。
圖2
椎骨缺損修復模擬預測
椎骨密度預測結果如(a)所示；缺損處分別植入可降解均勻(b)和非均勻(c)支架；均勻(d)和非均勻(e)支架的骨缺損修復預測結果
Figure2. Simulation prediction of vertebra defect repairpredicted outcome of vertebra density was showed as (a); implanting biodegradable homogeneous (b) and inhomogeneous (c) scaffold respectively for defect; predicting outcome of bone defect repair with homogeneous (d) and inhomogeneous (e) scaffold
4 結論
通過構建骨成形和支架降解控制方程為主體的骨缺損修復預測模型，并借助有限元的手段實現仿真，趨勢性地預測了骨缺損修復過程和修復后的骨結構形態。通過對預測結果的比較分析，可得出以下結論：
(1)支架降解骨生長模型和有限元法相結合是預測骨缺損修復過程及效果的有效方法。
(2)支架結構對骨缺損修復效果具有較大影響。
(3)可以利用支架的非均勻結構來控制骨缺損修復效果。
本文的模型還略顯粗糙，在整個模擬過程中有很多影響結果的參數，這些參數的力學和生理意義將進一步研究，以期獲得更加符合實際的結果。
