角膜的剛度降低是導致圓錐角膜的關鍵因素，角膜中膠原纖維剛度和纖維分散程度與角膜生物力學性能密切相關。本文通過建立基于角膜微結構即考慮膠原纖維的準分子激光原位角膜磨鑲術（LASIK）術前和術后人眼角膜有限元模型，模擬了可視化角膜生物力學分析儀的檢測過程，并與臨床實際檢測結果對比，確定屈光手術前后人眼角膜材料的超彈性本構參數及角膜膠原纖維剛度模量。研究發現 LASIK 術后，角膜膠原纖維剛度模量較術前明顯增大，并與中央角膜厚度（CCT）高度相關；術前術后角膜膠原纖維剛度模量與對應中央角膜厚度的預測關系分別為：k1 術前 = exp（9.14 ? 0.009CCT術前），k1 術后 = exp（8.82 ? 0.008CCT術后）。根據本文的研究結果，可由患者的中央角膜厚度值估算其術前術后膠原纖維剛度模量，進而建立與患者實際情況更相符的個性化角膜模型，為更準確地預測角膜安全手術切削量提供理論參考。
引用本文： 孟喬宇, 王曉君, 陳維毅, 李曉娜, 賀瑞. Corvis ST 檢測在體角膜生物力學性能的有限元模擬. 生物醫學工程學雜志, 2020, 37(4): 608-613. doi: 10.7507/1001-5515.201911081 復制
引言
我國是世界上近視發病率最高的國家，發病人數多，近視治療需求量龐大[1]。矯正近視的方法有多種，如角膜接觸鏡、框架眼鏡、角膜屈光術等。角膜屈光術手術時間短，恢復時間快，治療效果明顯，因而備受近視患者青睞。角膜屈光力占全眼球的 70%，角膜屈光術通過切削角膜基質層，改變角膜曲率來調節屈光度[2]。但手術切削破壞了角膜的完整性，角膜受力狀態發生改變，進而引起角膜生物力學性能的變化[3]。
角膜基質層占角膜厚度的 90%，由 300 多層膠原纖維板層構成，是角膜生物力學的主要承載者[4]。角膜復雜的力學行為主要依賴于膠原纖維，膠原纖維對角膜形狀的維持以及角膜生物力學性能的穩定起著至關重要的作用[5]。屈光手術后角膜膠原纖維數量減少、纖維剛度改變，導致角膜穩定性下降，當纖維剛度的改變超過角膜所能承受范圍時，可能引發角膜膨隆甚至圓錐角膜等并發癥[6]。因此，膠原纖維剛度對角膜穩定性具有十分重要的意義。
角膜生物力學性能與角膜的穩定性密切相關，因此角膜生物力學特性的準確測量對于防止手術過矯有重要意義。臨床中常采用可視化角膜生物力學分析儀（corneal visualization Scheimpflug technology，Corvis ST）測量在體角膜生物力學性能[7]。Corvis ST 利用成像技術記錄了角膜在脈沖氣流作用下的整個變形過程，可測得中央角膜厚度（central corneal thickness，CCT）、修正眼內壓（biomechanical corrected intraocular pressure，bIOP）、角膜頂點位移、角膜變形幅度等多個參數，能夠實時動態記錄角膜受壓形變的整個過程，檢測準確度比較高，是目前受到推崇的在體檢測設備。但 Corvis ST 僅可測得角膜宏觀生物力學指標，無法對角膜膠原纖維微結構以及微觀力學性能進行分析檢測。
因此，本文建立考慮膠原纖維的角膜微結構有限元模型，結合有限元模擬軟件和 Corvis ST 檢測儀器對術前術后角膜膠原纖維剛度進行研究，有助于臨床更好地優化屈光手術角膜切削量的計算，降低并發癥的風險，為臨床患者個性化評估切削量從而實現個性化手術提供一定的理論參考。
1 方法
1.1 數據采集
納入標準：① 年齡 18 周歲以上；② 屈光度穩定 2 年以上，且每年屈光度增長不超過 0.5 D；③ 雙眼近視等效球鏡度數均不超過－10.00 D，柱鏡度數均不高于－4.00 D；④ 術前角膜厚度 ≥ 460 μm，且術后基質床剩余厚度 ≥ 250 μm，術前眼壓 ≤ 24 mm Hg；⑤ 軟性角膜接觸鏡停戴 1 周以上，角膜塑形鏡停戴 8 周以上；⑥ 既往無眼部手術史或外傷史，排除圓錐角膜、感染性眼病、青光眼及視網膜疾病等其他眼部疾患；⑦ 排除結締組織病、瘢痕體質、免疫性疾病以及精神過分緊張等不適合手術的患者。本研究經山西省眼科醫院倫理委員會批準通過。
選取于山西省眼科醫院行準分子激光角膜原位磨鑲術（laser assisted in situ keratomileusis，LASIK）的近視患者 34 例（67 眼），其中男 15 例（30 眼）、女 19 例（37 眼）。所有患者術前經 Corvis ST 檢查，術后一個月復查，所有檢查均由同一專業人員進行操作，并獲取術前術后角膜真實頂點位移、CCT 值、bIOP 值等生物力學參數。
1.2 幾何模型的建立
在有限元軟件（ABAQUS 6.14）中建立三維軸對稱角膜模型，采用 8 節點六面體實體單元（C3D8H）。術前術后角膜模型采用統一的前表面曲率半徑 R、高度 H、橫徑 L，如圖 1 中角膜截面圖所示，其中 R = 7.97 mm，H = 2.5 mm，L = 11.49 mm[8]。術前角膜模型采用均勻厚度，根據不同患者術前 CCT 值建立對應的術前模型，如圖 1 中術前角膜幾何模型所示。根據手術切削量在術前模型頂點中心 3.5 mm 半徑范圍內進行切削，得到術后模型，圖 1 展示了術后四分之一角膜模型示意圖。角膜密度設為 1 149 kg/m3[9]。
圖1
角膜模型
Figure1.
Model of cornea
1.3 本構方程的選擇
角膜視為各向異性、超彈性、不可壓縮性材料[4]。Gasser、Holzapfel 和 Ogden 共同提出了一種應變能函數，即 Holzapfel-Gasser-Ogden 超彈性模型（HGO 模型），用以描述角膜本構關系[10]。HGO 模型考慮了膠原纖維分布對應變能函數的影響，對角膜力學行為的描述更為準確可靠。本構方程描述如下：
 |
式中 
為應變能密度函數，第一項表示基質的扭轉應變能，
為基質相關材料參數，代表基質的剪切模量，
是第一格林應變不變量。方程第二項表示角膜的體積應變能，
是和溫度相關的材料參數，
代表彈性體積比，角膜視為不可壓縮材料，所以 
= 1，因此第二項為 0。方程第三項代表膠原纖維的各向異性，表示不同膠原纖維族對方程的影響，
代表膠原纖維的剛度模量，
代表膠原纖維非線性的無量綱常數，
表示纖維族數，此處 
= 2。
 |
為膠原纖維族在纖維分布主方向上的變形，
表示膠原纖維的分散程度，0 ≤ 
< 1/3，
= 0 代表膠原纖維嚴格按膠原纖維主方向分布，
= 1/3 代表膠原纖維隨機分布[9]。其中：
 |
 |
 |
 |
 |
 |
式中 
表示右柯西格林應變張量 
與纖維族的平均纖維角度方向向量 
的偽不變量。
為第一主應變不變量，
為主方向伸長率，
為材料的變形梯度，
為右柯西格林張量，
是纖維束的方向向量。
角膜基質層中富含膠原纖維，LASIK 手術通常在膠原纖維層進行切削[11]，所以我們更多地關注膠原纖維的生物力學性能變化。Ariza-Gracia 等[12]和 Xu 等[10]的研究認為在脈沖氣流檢測過程中 
值對角膜頂點位移的影響程度最大。Ariza-Gracia 等[13]對有限元模型進行敏感性分析后認為 CCT 值、bIOP 值和 
值對角膜頂點位移影響最大，且在生理眼壓范圍內，
值和 
值對頂點位移的響應貢獻幾乎可以忽略，因此參考 Ariza-Gracia 等[12]的材料參數，固定 
值、
值、
值、
值這四個材料參數，通過調節 
值，使模擬計算的角膜頂點位移隨時間變化曲線與 Corvis ST 檢測的曲線吻合度達到最高，反推出不同人眼角膜術前術后的膠原纖維剛度模量值。參數設置見表 1。
表1
材料參數值
Table1.
Material parameter value
1.4 檢測過程的模擬
角膜在 Corvis ST 檢測時同時受眼內壓和脈沖氣流荷載，因此在模型內部施加均布荷載模擬眼內壓，在模型外表面施加脈沖氣流荷載，脈沖氣流荷載采用 Wang 等[14]測得的壓力值。約束角膜緣處的位移，不約束轉角。角膜受眼壓和脈沖氣流作用下歷經壓平、凹陷、第二次壓平、回彈的過程，模擬變形過程如圖 2 所示。對有限元模型的網格密度做了驗證分析，當網格數量達到 32 000 個以上時，模型計算結果趨于穩定。
圖2
角膜變形過程
Figure2.
Corneal deformation process
1.5 統計學分析
使用 SPSS 20.0 軟件對結果進行統計學分析。描述性統計資料采用均數±標準差表示；術前組和術后組數據比較采用配對樣本 t 檢驗分析，術前術后的修正眼內壓、角膜厚度與膠原纖維剛度模量之間的相關性分析采用 Pearson 方法檢驗。P < 0.05 認為差異具有統計學意義。統計中采用雙眼數據，未考慮雙眼間關聯。
2 結果
2.1 術前、術后模型有效性驗證
對術前、術后角膜有限元模型進行計算，對比角膜頂點位移的模擬曲線與 Corvis ST 檢測曲線，分別計算出 67 例術前角膜膠原纖維剛度模量（k1 術前）和 67 例術后角膜膠原纖維剛度模量（k1 術后）。以術前角膜厚度為 560 μm、術后角膜厚度為 446 μm 的模型為例，角膜頂點位移與時間關系的模擬與檢測曲線對比如圖 3 所示，從圖中可以看出模擬曲線與檢測曲線的匹配度較高，同時驗證了本模型的可靠性。術前、術后所有模型的模擬曲線與檢測曲線均類似，不一一列舉。
圖3
模擬與檢測角膜頂點位移－時間曲線對比
Figure3.
Comparison of corneal vertex displacement-time curves between simulation and detection
2.2 LASIK 手術前后參數比較分析
對 LASIK 手術前后的 CCT 值、bIOP 值和 
值進行對比分析，如表 2 所示。術后 CCT 值和 bIOP 值較術前顯著減小，
值較術前顯著增大，差異均具有統計學意義（P < 0.05）。
表2
LASIK 手術前后 CCT 值、bIOP 值、
值的對比
Table2.
Comparison of CCT, bIOP, 
before and after LASIK surgery
LASIK 手術前后各指標間的相關性如表 3 所示，相關系數用 r 表示。k1 術前 與術前 CCT 值（CCT術前）、k1 術后 與術后 CCT 值（CCT術后）的關系圖如圖 4 所示，k1 術前、k1 術后 均隨對應角膜厚度的增加而減小。結合表 3 數據可看出，k1 術前 與 CCT術前、k1 術后 與 CCT術后 高度相關，擬合方程分別為：k1 術前 = exp（9.14 ? 0.009CCT術前），k1 術后 = exp（8.82 ? 0.008CCT術后），方程的擬合優度均大于 0.82。
表3
LASIK 術前后 CCT 值、bIOP 值、
值之間的相關關系
Table3.
The correlation between CCT, bIOP, 
before and after LASIK surgery
圖4
k1 與 CCT 的關系
Figure4.
The relationship between k1 and CCT
k1 術前 與術前 bIOP 值（bIOP術前）、k1 術后 與術后對應 bIOP 值（bIOP術后）的關系圖如圖 5 所示，統計結果表明：k1 術前 與 bIOP術前 不相關（P>0.05）；bIOP術后 與 k1 術后 低度相關（r = ? 0.323，P < 0.01）。
圖5
k1 與 bIOP 的關系
Figure5.
The relationship between k1 and bIOP
3 討論
本文建立了基于角膜微結構即考慮角膜膠原纖維剛度的有限元模型，結合 Corvis ST 檢測數據和有限元分析軟件，將 LASIK 手術前后檢測的角膜頂點位移曲線與有限元計算結果對比，通過調整有限元模型中膠原纖維剛度模量等本構參數，得到 LASIK 屈光手術前后人眼角膜膠原纖維剛度模量及其與對應中央角膜厚度的關系。
研究發現術后膠原纖維剛度模量較術前顯著增大，術后膠原纖維剛度模量與術后中央角膜厚度呈顯著負相關。Zhang 等[15]取兔眼 LASIK 術后一個月角膜條進行單軸拉伸，發現術后角膜彈性模量隨剩余基質層厚度的減小而增大，且術后角膜條帶彈性模量顯著高于術前。本課題組前期研究也發現兔眼 LASIK 術后 28 天角膜彈性模量隨剩余基質床厚度的減小而增大[11]，與本文結果具有一致性。術后角膜經切削，中央組織變薄變平，手術切削引起傷口區域的膠原纖維卷曲以及基質層膠原纖維束和細胞外基質的重塑[16]，角膜受到的拉伸應力大于術前，可能導致膠原纖維發生適應性改變，導致膠原纖維剛度模量增大。
趙科超[17]對 24 例 LASIK 術后患者做整體剛度分析，發現術后整體剛度顯著小于術前整體剛度，這與本文所求角膜材料中術后膠原纖維剛度較術前增大并不沖突，手術切削破壞了角膜完整性，使角膜厚度變薄，雖然膠原纖維在術后傷口愈合過程中剛度增大，但角膜厚度對角膜結構穩定性起的主導作用可能大于膠原纖維的作用，隨著角膜厚度的變薄，角膜整體穩定性和整體剛度還是呈下降趨勢。對角膜進行整體生物力學分析時發現，術前術后角膜整體剛度與對應眼內壓均高度相關[17]，而本文研究樣本范圍內發現術前眼壓和術前纖維剛度模量沒有相關性，術后眼壓和術后纖維剛度模量相關性較低，且考慮到眼壓在一天之內有一定波動，故本文未考慮纖維剛度模量與眼內壓的關系。
本文所選樣本的角膜切削量范圍為 33～173 μm，切削量占術前中央角膜厚度比例為 7%～33%。按切削比例分為 3 組，切削量介于 10%～20% 的為 1 組（共 31 例眼），切削量介于 20%～30% 的為 2 組（共 32 例眼），切削量大于 30% 的為 3 組（共 4 例眼）。本課題組對兔眼 LASIK 術后的研究發現，術后一周時，不同切削量組間彈性模量沒有顯著差異；術后一個月時，切削量為 70% 組的角膜彈性模量顯著高于切削量為 50% 組的彈性模量[11]。本文對術后一個月各組的角膜膠原纖維剛度模量進行計算，發現兩兩組間均沒有顯著差異（1 組 vs 2 組 P = 0.361，2 組 vs 3 組 P = 0.467，1 組 vs 3 組 P = 0.233）。這可能是由于本文所選樣本的切削量主要集中于 10%～30%，1 組和 2 組的樣本數量占了總樣本量的 94%，切削量大于 30% 的眼僅 4 例，若擴大樣本的切削量采集范圍，則可能與較小切削量組之間產生顯著差異。
本文在研究中做了一些假定，首先是在建立角膜幾何模型時未考慮不同角膜幾何形態的差異，比如不同角膜的高度和表面曲率并不完全相同，并且不同角膜在脈沖氣流作用下受到的壓力也可能是不同的，但為簡化計算，使用了統一的角膜幾何模型參數和脈沖氣流荷載，同時對所有角膜邊界條件做了約束角膜緣位移的假定；其次，對模型試算時發現 
值和 
值對角膜頂點位移的影響很小，對匹配檢測曲線的影響不大，因而在賦予本構參數時未對這兩個參數做調整，在不同角膜中均設置為同一固定值，但實際上不同角膜中的 
值和 
值并不完全相同。每個角膜中 
值的求解都同時受邊界條件、載荷以及其他參數等多重因素疊加的影響。本文在計算時為得到術前術后角膜中 
值的變化趨勢，假定角膜中其他因素固定不變，未對其他條件和參數做進一步優化處理，所求出的 
值也并非精確值。
對于個體差異較大的人體組織而言，做精確預測存在一定的困難，需同時考慮的因素較多。同樣地，在臨床中對不同患者使用同樣的手術安全切削量評價標準并不十分合理。本研究探尋了 LASIK 手術前后膠原纖維剛度模量的變化及其與中央角膜厚度的關系，以期在臨床中根據患者的中央角膜厚度值估算其術前術后膠原纖維剛度模量，建立與患者實際情況更相符的個性化角膜有限元模型，為更準確地預測角膜安全手術切削量提供理論參考。
利益沖突聲明：本文全體作者均聲明不存在利益沖突。
引言
我國是世界上近視發病率最高的國家，發病人數多，近視治療需求量龐大[1]。矯正近視的方法有多種，如角膜接觸鏡、框架眼鏡、角膜屈光術等。角膜屈光術手術時間短，恢復時間快，治療效果明顯，因而備受近視患者青睞。角膜屈光力占全眼球的 70%，角膜屈光術通過切削角膜基質層，改變角膜曲率來調節屈光度[2]。但手術切削破壞了角膜的完整性，角膜受力狀態發生改變，進而引起角膜生物力學性能的變化[3]。
角膜基質層占角膜厚度的 90%，由 300 多層膠原纖維板層構成，是角膜生物力學的主要承載者[4]。角膜復雜的力學行為主要依賴于膠原纖維，膠原纖維對角膜形狀的維持以及角膜生物力學性能的穩定起著至關重要的作用[5]。屈光手術后角膜膠原纖維數量減少、纖維剛度改變，導致角膜穩定性下降，當纖維剛度的改變超過角膜所能承受范圍時，可能引發角膜膨隆甚至圓錐角膜等并發癥[6]。因此，膠原纖維剛度對角膜穩定性具有十分重要的意義。
角膜生物力學性能與角膜的穩定性密切相關，因此角膜生物力學特性的準確測量對于防止手術過矯有重要意義。臨床中常采用可視化角膜生物力學分析儀（corneal visualization Scheimpflug technology，Corvis ST）測量在體角膜生物力學性能[7]。Corvis ST 利用成像技術記錄了角膜在脈沖氣流作用下的整個變形過程，可測得中央角膜厚度（central corneal thickness，CCT）、修正眼內壓（biomechanical corrected intraocular pressure，bIOP）、角膜頂點位移、角膜變形幅度等多個參數，能夠實時動態記錄角膜受壓形變的整個過程，檢測準確度比較高，是目前受到推崇的在體檢測設備。但 Corvis ST 僅可測得角膜宏觀生物力學指標，無法對角膜膠原纖維微結構以及微觀力學性能進行分析檢測。
因此，本文建立考慮膠原纖維的角膜微結構有限元模型，結合有限元模擬軟件和 Corvis ST 檢測儀器對術前術后角膜膠原纖維剛度進行研究，有助于臨床更好地優化屈光手術角膜切削量的計算，降低并發癥的風險，為臨床患者個性化評估切削量從而實現個性化手術提供一定的理論參考。
1 方法
1.1 數據采集
納入標準：① 年齡 18 周歲以上；② 屈光度穩定 2 年以上，且每年屈光度增長不超過 0.5 D；③ 雙眼近視等效球鏡度數均不超過－10.00 D，柱鏡度數均不高于－4.00 D；④ 術前角膜厚度 ≥ 460 μm，且術后基質床剩余厚度 ≥ 250 μm，術前眼壓 ≤ 24 mm Hg；⑤ 軟性角膜接觸鏡停戴 1 周以上，角膜塑形鏡停戴 8 周以上；⑥ 既往無眼部手術史或外傷史，排除圓錐角膜、感染性眼病、青光眼及視網膜疾病等其他眼部疾患；⑦ 排除結締組織病、瘢痕體質、免疫性疾病以及精神過分緊張等不適合手術的患者。本研究經山西省眼科醫院倫理委員會批準通過。
選取于山西省眼科醫院行準分子激光角膜原位磨鑲術（laser assisted in situ keratomileusis，LASIK）的近視患者 34 例（67 眼），其中男 15 例（30 眼）、女 19 例（37 眼）。所有患者術前經 Corvis ST 檢查，術后一個月復查，所有檢查均由同一專業人員進行操作，并獲取術前術后角膜真實頂點位移、CCT 值、bIOP 值等生物力學參數。
1.2 幾何模型的建立
在有限元軟件（ABAQUS 6.14）中建立三維軸對稱角膜模型，采用 8 節點六面體實體單元（C3D8H）。術前術后角膜模型采用統一的前表面曲率半徑 R、高度 H、橫徑 L，如圖 1 中角膜截面圖所示，其中 R = 7.97 mm，H = 2.5 mm，L = 11.49 mm[8]。術前角膜模型采用均勻厚度，根據不同患者術前 CCT 值建立對應的術前模型，如圖 1 中術前角膜幾何模型所示。根據手術切削量在術前模型頂點中心 3.5 mm 半徑范圍內進行切削，得到術后模型，圖 1 展示了術后四分之一角膜模型示意圖。角膜密度設為 1 149 kg/m3[9]。
圖1
角膜模型
Figure1.
Model of cornea
1.3 本構方程的選擇
角膜視為各向異性、超彈性、不可壓縮性材料[4]。Gasser、Holzapfel 和 Ogden 共同提出了一種應變能函數，即 Holzapfel-Gasser-Ogden 超彈性模型（HGO 模型），用以描述角膜本構關系[10]。HGO 模型考慮了膠原纖維分布對應變能函數的影響，對角膜力學行為的描述更為準確可靠。本構方程描述如下：
|
式中 為應變能密度函數，第一項表示基質的扭轉應變能， 為基質相關材料參數，代表基質的剪切模量， 是第一格林應變不變量。方程第二項表示角膜的體積應變能， 是和溫度相關的材料參數， 代表彈性體積比，角膜視為不可壓縮材料，所以 = 1，因此第二項為 0。方程第三項代表膠原纖維的各向異性，表示不同膠原纖維族對方程的影響， 代表膠原纖維的剛度模量， 代表膠原纖維非線性的無量綱常數， 表示纖維族數，此處 = 2。
|
為膠原纖維族在纖維分布主方向上的變形， 表示膠原纖維的分散程度，0 ≤ < 1/3， = 0 代表膠原纖維嚴格按膠原纖維主方向分布， = 1/3 代表膠原纖維隨機分布[9]。其中：
|
|
|
|
|
|
式中 表示右柯西格林應變張量 與纖維族的平均纖維角度方向向量 的偽不變量。 為第一主應變不變量， 為主方向伸長率， 為材料的變形梯度， 為右柯西格林張量， 是纖維束的方向向量。
角膜基質層中富含膠原纖維，LASIK 手術通常在膠原纖維層進行切削[11]，所以我們更多地關注膠原纖維的生物力學性能變化。Ariza-Gracia 等[12]和 Xu 等[10]的研究認為在脈沖氣流檢測過程中 值對角膜頂點位移的影響程度最大。Ariza-Gracia 等[13]對有限元模型進行敏感性分析后認為 CCT 值、bIOP 值和 值對角膜頂點位移影響最大，且在生理眼壓范圍內， 值和 值對頂點位移的響應貢獻幾乎可以忽略，因此參考 Ariza-Gracia 等[12]的材料參數，固定 值、 值、 值、 值這四個材料參數，通過調節 值，使模擬計算的角膜頂點位移隨時間變化曲線與 Corvis ST 檢測的曲線吻合度達到最高，反推出不同人眼角膜術前術后的膠原纖維剛度模量值。參數設置見表 1。
表1
材料參數值
Table1.
Material parameter value
1.4 檢測過程的模擬
角膜在 Corvis ST 檢測時同時受眼內壓和脈沖氣流荷載，因此在模型內部施加均布荷載模擬眼內壓，在模型外表面施加脈沖氣流荷載，脈沖氣流荷載采用 Wang 等[14]測得的壓力值。約束角膜緣處的位移，不約束轉角。角膜受眼壓和脈沖氣流作用下歷經壓平、凹陷、第二次壓平、回彈的過程，模擬變形過程如圖 2 所示。對有限元模型的網格密度做了驗證分析，當網格數量達到 32 000 個以上時，模型計算結果趨于穩定。
圖2
角膜變形過程
Figure2.
Corneal deformation process
1.5 統計學分析
使用 SPSS 20.0 軟件對結果進行統計學分析。描述性統計資料采用均數±標準差表示；術前組和術后組數據比較采用配對樣本 t 檢驗分析，術前術后的修正眼內壓、角膜厚度與膠原纖維剛度模量之間的相關性分析采用 Pearson 方法檢驗。P < 0.05 認為差異具有統計學意義。統計中采用雙眼數據，未考慮雙眼間關聯。
2 結果
2.1 術前、術后模型有效性驗證
對術前、術后角膜有限元模型進行計算，對比角膜頂點位移的模擬曲線與 Corvis ST 檢測曲線，分別計算出 67 例術前角膜膠原纖維剛度模量（k1 術前）和 67 例術后角膜膠原纖維剛度模量（k1 術后）。以術前角膜厚度為 560 μm、術后角膜厚度為 446 μm 的模型為例，角膜頂點位移與時間關系的模擬與檢測曲線對比如圖 3 所示，從圖中可以看出模擬曲線與檢測曲線的匹配度較高，同時驗證了本模型的可靠性。術前、術后所有模型的模擬曲線與檢測曲線均類似，不一一列舉。
圖3
模擬與檢測角膜頂點位移－時間曲線對比
Figure3.
Comparison of corneal vertex displacement-time curves between simulation and detection
2.2 LASIK 手術前后參數比較分析
對 LASIK 手術前后的 CCT 值、bIOP 值和 值進行對比分析，如表 2 所示。術后 CCT 值和 bIOP 值較術前顯著減小， 值較術前顯著增大，差異均具有統計學意義（P < 0.05）。
表2
LASIK 手術前后 CCT 值、bIOP 值、 值的對比
Table2.
Comparison of CCT, bIOP, before and after LASIK surgery
LASIK 手術前后各指標間的相關性如表 3 所示，相關系數用 r 表示。k1 術前 與術前 CCT 值（CCT術前）、k1 術后 與術后 CCT 值（CCT術后）的關系圖如圖 4 所示，k1 術前、k1 術后 均隨對應角膜厚度的增加而減小。結合表 3 數據可看出，k1 術前 與 CCT術前、k1 術后 與 CCT術后 高度相關，擬合方程分別為：k1 術前 = exp（9.14 ? 0.009CCT術前），k1 術后 = exp（8.82 ? 0.008CCT術后），方程的擬合優度均大于 0.82。
表3
LASIK 術前后 CCT 值、bIOP 值、 值之間的相關關系
Table3.
The correlation between CCT, bIOP, before and after LASIK surgery
圖4
k1 與 CCT 的關系
Figure4.
The relationship between k1 and CCT
k1 術前 與術前 bIOP 值（bIOP術前）、k1 術后 與術后對應 bIOP 值（bIOP術后）的關系圖如圖 5 所示，統計結果表明：k1 術前 與 bIOP術前 不相關（P>0.05）；bIOP術后 與 k1 術后 低度相關（r = ? 0.323，P < 0.01）。
圖5
k1 與 bIOP 的關系
Figure5.
The relationship between k1 and bIOP
3 討論
本文建立了基于角膜微結構即考慮角膜膠原纖維剛度的有限元模型，結合 Corvis ST 檢測數據和有限元分析軟件，將 LASIK 手術前后檢測的角膜頂點位移曲線與有限元計算結果對比，通過調整有限元模型中膠原纖維剛度模量等本構參數，得到 LASIK 屈光手術前后人眼角膜膠原纖維剛度模量及其與對應中央角膜厚度的關系。
研究發現術后膠原纖維剛度模量較術前顯著增大，術后膠原纖維剛度模量與術后中央角膜厚度呈顯著負相關。Zhang 等[15]取兔眼 LASIK 術后一個月角膜條進行單軸拉伸，發現術后角膜彈性模量隨剩余基質層厚度的減小而增大，且術后角膜條帶彈性模量顯著高于術前。本課題組前期研究也發現兔眼 LASIK 術后 28 天角膜彈性模量隨剩余基質床厚度的減小而增大[11]，與本文結果具有一致性。術后角膜經切削，中央組織變薄變平，手術切削引起傷口區域的膠原纖維卷曲以及基質層膠原纖維束和細胞外基質的重塑[16]，角膜受到的拉伸應力大于術前，可能導致膠原纖維發生適應性改變，導致膠原纖維剛度模量增大。
趙科超[17]對 24 例 LASIK 術后患者做整體剛度分析，發現術后整體剛度顯著小于術前整體剛度，這與本文所求角膜材料中術后膠原纖維剛度較術前增大并不沖突，手術切削破壞了角膜完整性，使角膜厚度變薄，雖然膠原纖維在術后傷口愈合過程中剛度增大，但角膜厚度對角膜結構穩定性起的主導作用可能大于膠原纖維的作用，隨著角膜厚度的變薄，角膜整體穩定性和整體剛度還是呈下降趨勢。對角膜進行整體生物力學分析時發現，術前術后角膜整體剛度與對應眼內壓均高度相關[17]，而本文研究樣本范圍內發現術前眼壓和術前纖維剛度模量沒有相關性，術后眼壓和術后纖維剛度模量相關性較低，且考慮到眼壓在一天之內有一定波動，故本文未考慮纖維剛度模量與眼內壓的關系。
本文所選樣本的角膜切削量范圍為 33～173 μm，切削量占術前中央角膜厚度比例為 7%～33%。按切削比例分為 3 組，切削量介于 10%～20% 的為 1 組（共 31 例眼），切削量介于 20%～30% 的為 2 組（共 32 例眼），切削量大于 30% 的為 3 組（共 4 例眼）。本課題組對兔眼 LASIK 術后的研究發現，術后一周時，不同切削量組間彈性模量沒有顯著差異；術后一個月時，切削量為 70% 組的角膜彈性模量顯著高于切削量為 50% 組的彈性模量[11]。本文對術后一個月各組的角膜膠原纖維剛度模量進行計算，發現兩兩組間均沒有顯著差異（1 組 vs 2 組 P = 0.361，2 組 vs 3 組 P = 0.467，1 組 vs 3 組 P = 0.233）。這可能是由于本文所選樣本的切削量主要集中于 10%～30%，1 組和 2 組的樣本數量占了總樣本量的 94%，切削量大于 30% 的眼僅 4 例，若擴大樣本的切削量采集范圍，則可能與較小切削量組之間產生顯著差異。
本文在研究中做了一些假定，首先是在建立角膜幾何模型時未考慮不同角膜幾何形態的差異，比如不同角膜的高度和表面曲率并不完全相同，并且不同角膜在脈沖氣流作用下受到的壓力也可能是不同的，但為簡化計算，使用了統一的角膜幾何模型參數和脈沖氣流荷載，同時對所有角膜邊界條件做了約束角膜緣位移的假定；其次，對模型試算時發現 值和 值對角膜頂點位移的影響很小，對匹配檢測曲線的影響不大，因而在賦予本構參數時未對這兩個參數做調整，在不同角膜中均設置為同一固定值，但實際上不同角膜中的 值和 值并不完全相同。每個角膜中 值的求解都同時受邊界條件、載荷以及其他參數等多重因素疊加的影響。本文在計算時為得到術前術后角膜中 值的變化趨勢，假定角膜中其他因素固定不變，未對其他條件和參數做進一步優化處理，所求出的 值也并非精確值。
對于個體差異較大的人體組織而言，做精確預測存在一定的困難，需同時考慮的因素較多。同樣地，在臨床中對不同患者使用同樣的手術安全切削量評價標準并不十分合理。本研究探尋了 LASIK 手術前后膠原纖維剛度模量的變化及其與中央角膜厚度的關系，以期在臨床中根據患者的中央角膜厚度值估算其術前術后膠原纖維剛度模量，建立與患者實際情況更相符的個性化角膜有限元模型，為更準確地預測角膜安全手術切削量提供理論參考。
利益沖突聲明：本文全體作者均聲明不存在利益沖突。
