基于亞區結構的人體二尖瓣的建模與仿真

趙曉佳,宋 曉,黃曉陽*,蘇茂龍,鐘 琪

(1.廈門大學信息科學與技術學院,福建 廈門361005;2.南陽師范學院計算機與信息技術學院,河南南陽473061;3.廈門大學附屬中山醫院超聲醫學科,福建廈門361004)

基于亞區結構的人體二尖瓣的建模與仿真

趙曉佳1,宋 曉2,黃曉陽1*,蘇茂龍3,鐘 琪1

(1.廈門大學信息科學與技術學院,福建 廈門361005;2.南陽師范學院計算機與信息技術學院,河南南陽473061;3.廈門大學附屬中山醫院超聲醫學科,福建廈門361004)

通過建立精細的人體二尖瓣有限元模型,仿真前后瓣葉厚度不均勻的二尖瓣模型在心臟收縮初期的動態閉合過程,與前后瓣葉厚度均勻的二尖瓣模型的仿真結果相比較.仿真結果表明:兩種模型的瓣葉應力分布不均勻,最大應力值都出現在后瓣葉亞區與亞區之間的裂口位置,與臨床二尖瓣撕裂常發生位置一致,表明了仿真方法以及結果的正確性,瓣葉厚度不均勻模型中瓣葉應力水平相對較低,更符合實際.最后根據臨床數據,設置厚度不均勻模型中瓣環與乳頭肌的二維運動,再次進行仿真.仿真結果表明該精細模型瓣葉的閉合度比之前瓣環與乳頭肌固定不動模型的閉合度高,瓣葉應力分布相對較均勻,更接近實際生理情況,為后續的病理等的仿真奠定了基礎.

二尖瓣;有限元模型;應力分布;最大應力

二尖瓣是人體心臟的重要瓣膜組織之一,它位于左心房與左心室之間,由二尖瓣環、兩個瓣葉、腿索、乳頭肌等結構組成[1].二尖瓣的瓣葉分為較大的前瓣葉和較小的后瓣葉,后瓣葉可細化為3個較窄的亞區結構[2]:P1、P2、P3.二尖瓣的瓣葉與瓣環、腿索相連,腿索又與乳頭肌相連.根據Lam等[3]的分類方式將腿索分為3級:1級腿索連接于瓣葉的游離緣,2級腿索連接于瓣葉的粗糙部,3級腿索從室壁連接到后瓣葉的基底部.二尖瓣維持正常的生理職能,需要各個組成結構的完整、乳頭肌的正常收縮及腿索的牽引[4].

二尖瓣的生物力學是一個復雜的問題:結構體的大變形,瓣膜組織的各向異性的非線性彈性行為以及心動周期內的血液動力學的載荷等.建模仿真是分析瓣膜功能有效性的技術方法.隨著科技的發展、仿真技術的不斷完善,人們可以對瓣膜功能有更加精確的了解,可以選擇更適合的病理治療方式.通過對人體二尖瓣進行建模仿真,可以更好地認識二尖瓣的力學機理,對后期分析二尖瓣關閉不全等病理結構以及指導瓣膜修復術具有重要的意義.

目前關于二尖瓣的研究分為結構模型和流固耦合模型兩類.結構模型是指模型中只有二尖瓣的結構組件,不考慮血液;流固耦合模型是指二尖瓣的結構與血液流體域相互作用的模型,它主要研究二尖瓣結構在流體作用下的各種行為以及結構對流體域的影響.其中結構模型是流固耦合模型的基礎,主要研究的學者有Kunzelman團隊[5-7]、Lim團隊[8]等.Kunzelman等[5]建立了最原始的人體二尖瓣模型,模擬了二尖瓣的閉合過程.Lim等[8]建立厚度均勻的膜結構瓣葉,模擬正常二尖瓣工作機理,但該模型中瓣葉自由邊上無尖端,忽略了瓣葉形態的局部細節.Kunzelman等[7]、羅云等[9]將二尖瓣的后瓣葉簡化為1個瓣尖,成功模擬出瓣葉的開閉過程,但并未根據實際生理情況將后瓣葉分為3個亞區.這些研究都沒有模擬瓣環與乳頭肌的二維運動過程.本文研究的目的,就是為了比較更加精細的模型對仿真結果的影響.

本文根據人體二尖瓣的實際生理結構,以二尖瓣的亞區結構特征,建立比較精細的二尖瓣三維有限元模型:后瓣葉分為3個亞區結構,腿索分為游離邊緣腿索,粗糙帶以及基底部腿索.二尖瓣瓣膜厚度分為兩種情況:1)前后瓣葉厚度相同的厚度均勻模型; 2)前后瓣葉厚度不相同的厚度不均勻模型.兩種模型分別模擬其瓣葉在心臟收縮初期心房心室壓力差的作用下的閉合過程,并比較了兩種仿真模型的差異.在瓣葉厚度不均勻的模型基礎上,按照臨床數據及相關描述,設置瓣環與乳頭肌的二維運動過程,更加精細地模擬二尖瓣的閉合過程.

1 方 法

采用ANSXS建立二尖瓣有限元模型,并用LSDXNA求解器進行求解,得到二尖瓣模型的應力隨時間的變化的結果,用LS-PREPOST后處理器生成可視化結果.

1.1 二尖瓣幾何模型

二尖瓣模型由瓣環、瓣葉、腿索和乳頭肌組成.本文根據文獻[3,10-11]建立的模型如圖1所示:xy平面的瓣環呈D字,乳頭肌簡化為2個固定的點,瓣葉從瓣環向z軸正方向延伸,后瓣葉細分為3個亞區結構.與瓣葉游離邊緣相連的腿索有24根,與后瓣葉基底部相連的腿索有2根,與粗糙帶相連的腿索分支有18根.

圖1 人體二尖瓣幾何模型Fig.1 Geometric model of human mitral valve

1.2 物理性質

二尖瓣的所有組織均設置為各向同性的彈性材料.瓣葉密度為1 040 kg/m3[10],泊松比設置為0.45,彈性模量設置為4 MPa[12];腿索的密度為1 040 kg/ m3[10],其應力-應變曲線如圖2所示[13].

1.3 網格劃分

在ANSXS中對二尖瓣模型進行網格劃分,采用自由網格劃分,共劃分為4 396個Belytschko殼單元;腿索橫截面積均設置為0.4 mm2[14],共劃分為760個桁架單元.模型中瓣膜的厚度分2種情況討論:1)瓣葉厚度均勻的模型:前后瓣葉厚度均設置為1.3 mm; 2)前后瓣葉厚度不相同的厚度不均勻模型:根據文獻[15-16],前瓣葉厚度設置為1.32 mm,后瓣葉厚度設置為1.26 mm.兩種情況下的網格劃分方式一樣,僅前后瓣葉厚度不同.

圖2 二尖瓣腿索材料應力-應變曲線Fig.2 Stress-strain behavior of mitral valve chordae tendineae

1.4 邊界與負載條件

基于二尖瓣的面對稱結構,在對稱面上節點x軸方向位移為0.二尖瓣瓣環與乳頭肌在整個心動周期內的位置相對固定,因此假設瓣環和乳頭肌固定不動,瓣環上節點的位移為0.

左心室與左心房在單個心動周期內壓力值的臨床數據[17]如圖3所示,由于二尖瓣在心臟收縮初期壓力斜坡(early systolic pressure ramp,ESPR)區,瓣膜載荷會隨著時間動態增長至最大值,此時的二尖瓣較容易出現應力集中現象,因此著重研究ESPR區瓣膜的應力分布情況.通過測量圖3上同一時刻心室與心房的壓強,然后求得它們之間的差值即左心室房壓強差,這個左心室房壓強差就是施加在二尖瓣上的動態載荷,如圖4所示.

圖3 心室心房壓力曲線Fig.3 Atrial and ventricular pressure profiles

2 結 果

2.1 二尖瓣瓣葉閉合

分別模擬兩種模型在ESPR區二尖瓣瓣葉的閉合情況,仿真結果如圖5所示.

由圖5可知,2個模型中的瓣葉均在左房室壓差以及腿索的共同作用下慢慢對合,最終實現前后瓣葉的接觸.兩個模型中瓣環的形狀沒有發生變化,模型的閉合程度基本一樣,無顯著差異.

圖4 二尖瓣載荷曲線Fig.4 Transvalvular pressure load applied to the mitral valve

圖5 二尖瓣瓣葉閉合情況Fig.5 Configuration of the valve at closure

2.2 二尖瓣瓣葉應力分布

分別對二尖瓣在ESPR區應力分布進行仿真,對應的最大von Mises應力分布如圖6所示.

如圖所示,在ESPR區400 ms時刻,左房室壓差趨于穩定,約16 k Pa,兩個模型的二尖瓣von Mises應力均達到最大,最大應力值均出現在亞區的裂口位置,與臨床上瓣葉容易撕裂的位置一致[18],與鐘琪等[19]的結果一致.表1列出了兩種模型的仿真結果對比情況.兩個模型應力分布都不太均勻,均出現了應力集中現象,兩個模型除了后瓣葉亞區間應力最大外,前后瓣葉以及瓣葉與瓣環連合處的von Mises應力也都比較大,但厚度不均勻模型的總體von Mises應力水平比厚度均勻模型的von Mises應力水平低,更不容易出現撕裂的狀況,與實際生理相符.

表1 2種模型仿真結果對比Tab.1 Comparison of the results between two models

2.3 瓣環和乳頭肌的運動

根據文獻[20-23]對瓣環運動的描述,瓣環在xy平面的移動可以描述為:1)前瓣環長度增加11%,但是瓣環總長度縮短13%;2)瓣環橫向長度(即x軸長度)增加5%.

根據文獻[16]描述,以及整個ESPR區,瓣環與乳頭肌的位置是相對不動,乳頭肌沿x軸朝著前后瓣葉連合處移動0.85 mm,沿y軸朝著前瓣葉方向移動1.10 mm.

在厚度不均勻模型基礎上,仿真瓣環以及乳頭肌在xy平面的移動以及瓣葉的閉合過程,仿真結果如圖7所示.

圖6 400 ms時刻瓣葉最大von Mises應力分布圖Fig.6 Maximum von Mises stress distribution on the valve leaflets at 400 ms

如圖7所示,ESPR區初期,二尖瓣的von Mises應力比較低,左心房室的壓力差約為0 kPa;在300 ms時刻,左心房室的壓差約為10 kPa,二尖瓣的von Mises應力分布在0.042～0.253 MPa之間;在375 ms時刻,左心房室壓差約為16 kPa,二尖瓣的von Mises應力主要在0.047～0.280 MPa,腿索與后瓣葉的連接處及前后瓣葉連合處出現應力集中現象;在425 ms時刻,左心房室壓差約為16 kPa,二尖瓣的von Mises應力分布在0.048～0.291 MPa,基底腿索與后瓣葉的連接處及前后瓣葉連合處出現應力集中現象.總體看,二尖瓣瓣環在整個ESPR區的von Mises應力分布相對比較均勻,最大von Mises應力集中在腿索與瓣葉的連接處以及前后瓣葉之間的連合處,最大應力值為0.485 MPa.在二尖瓣瓣葉的閉合過程中,二尖瓣瓣環的形態發生了很大的變化,由原來的D字形變成了一個類似橢圓形,瓣環在x軸方向長度變長,y軸方向長度縮短.除了形態的變化,瓣葉的閉合度也明顯發生了變化,與原先厚度不均勻模型相比,瓣葉的閉合程度明顯提高,與實際生理情況符合,如表2所示.

圖7 不同時刻二尖瓣瓣葉von Mises應力分布圖Fig.7 von Mises stress distribution on the valve leaflets at different time

表2 瓣環與乳頭肌的運動對仿真結果的影響Tab.2 Impact of annulus and papillary muscles movement to the simulation results

3 結 論

為了研究比較精細模型對仿真結果的影響,本文參照實際生理結構,將二尖瓣的后瓣葉細分為3個亞區,建立了瓣膜厚度分布均勻和分布不均勻2種模型,對二尖瓣瓣葉在ESPR區的閉合過程進行了仿真.仿真結果表明,在整個ESPR區,瓣膜應力最大值出現在400 ms左右,位置在P2亞區與P1、P3亞區的裂口位置,此時心房心室的壓力差最大,約16 k Pa.厚度不均勻模型的應力水平比較低,不易出現撕裂現象,更接近實際的生理情況.在此基礎上,本文根據實際生理情況,設置瓣環與乳頭肌的運動,再次仿真瓣葉的閉合過程,仿真結果顯示前后瓣葉的連合處及腿索與后瓣葉連接處應力較大,仿真的結果與臨床二尖瓣撕裂研究[18]相符合,同時與Kunzelman等[14]的模擬結果接近,前瓣葉的最大應力值略小一些,約為0.300 MPa,后瓣葉的應力值均約為0.200 MPa.瓣葉的閉合度更高,更符合實際生理情況.本文研究結果表明,建立瓣膜厚度不均勻模型,并在仿真過程中設置瓣環與乳頭肌的運動,仿真模型及過程更符合人體真實二尖瓣情況,仿真結果也更接近臨床數據,更適合后續二尖瓣病理情況的仿真與分析.

瓣葉實際材料屬性是非線性各向異性的,本文瓣葉簡化為線性各向同性,以后的研究將細化瓣葉屬性.血液的流動對二尖瓣的行為產生影響,二尖瓣瓣膜在血液影響下的流固耦合運動情況將是后續值得深入研究的一個方向.

[1] 馮敏,黃曉陽,蘇茂龍,等.基于雙彈性腔的左心循環系統的建模與仿真[J].廈門大學學報:自然科學版,2012,51 (4):676-681.

[2] Dominik J,Zacek P.Heart valve surgery:an illustrated guide[M].Berlin:Spring-Verlag,2010:10-11.

[3] Lam J H,Ranganathan N,Wigle E D,et al.Morphology of the human mitral valve:Ⅰ.Chordae tendinese:a new classification[J].Circulation,1970,41(3):449-458.

[4] 黃化清.二尖瓣腿索斷裂[J].臨床心血管病雜志,1985,1 (1/2):80-83.

[5] Kunzelman K S,Cochran R P,Chuong C,et al.Finite element analysis of the mitral valve[J].J Heart Valve Dis, 1993,2:326-340.

[6] Kunzelman K S,Cochran R P,Verrier E D,et al.Anatomic basis for mitral valve modeling[J].J Heart Valve Dis, 1994,3:491-496.

[7] Kunzelman K S,Einstein D R,Cochran R P.Fluid-structure interaction models of the mitral valve:function in normal and pathological states[J].Philosophical Transactions of the Royal Society B:Biological Sciences,2007, 362(1484):1393-1406.

[8] Lim K H,Xeo J H,Duran C.Three-dimensional asymmetrical modeling of the mitral valve:a finite element study with dynamic boundaries[J].J Heart Valve Dis, 2005,14:386-392.

[9] 羅云.二尖瓣有限元建模與生物力學分析[D].廈門:廈門大學,2014.

[10] Lau K D,Diaz V,Scambler P,et al.Mitral valve dynamics in structural and fluid-structure interaction models [J].Medical Engineering&Physics,2010,32(9): 1057-1064.

[11] Ranganathan N,Lam J H C,Wigle E D,et al.Morphology of the human mitral valve:II.The valve leaflets[J]. Circulation,1970,41:459-467.

[12] Pan F D,Donzella G,Fucci C,et al.Structural effects of an innovative surgical technique to repair heart valve defects[J].J Biomech,2005,38:2460-2471.

[13] Kunzelman K S,Cochran K P.Mechanical properties of basal and marginal mitral valve chordae tendineae[J]. ASAIO Trans,1990,36(3):405-408.

[14] Kunzelman K S,Cochran R P,Murphree S S,et al.Differential collagen distribution in the mitral valve and its influence on biomechanical behaviour[J].J Heart Valve Dis,1993,2:236-244.

[15] 鐘琪.人體心臟二尖瓣有限元建模與仿真[D].廈門:廈門大學,2014.

[16] Stevanella M,Votta E,Redaelli A.Mitral valve finite element modeling:implicatioons of tissues′nonlinear response and annular motion[J].Journal of Bionechanical Engineering,2009,131:1-9.

[17] Levick J R.An introduction to cardiovascular physiology [M].London:Hodder Arnold Publication,2003.

[18] 盧曉芳,呂清,謝明星,等.實時三維超聲心動圖評價二尖瓣裂的初步研究[J].中國醫學影像技術,2006,22 (4):490-492.

[19] 鐘琪,曾文華,黃曉陽,等.人體二尖瓣建模及生物力學分析[J].醫用生物力學,2014,29(1):53-59.

[20] Lansac E,Lim K H,Shomura X,et al.Dynamic balance of the aortomitral junction[J].J Thorac Cardiovasc Surg,2002,123(5):911-918.

[21] Ormiston J A,Shah P M,Tei C,et al.Size and motion of the mitral valve annulus in man.I.A two-dimensional echocardiographic method and findings in normal subjects[J].Circulation,1981,64(1):113-120.

[22] Fyrenius A,Engvall J,Janerot-Sjoberg B.Major and minor axes of the normal mitral annulus[J].J Heart Valve Dis,2001,10(2):146-152.

[23] Gorman III J H,Jackson B M,Enomoto X,et al.The effect of regional ischemia on mitral valve annular saddle shape[J].Ann Thorac Surg,2004,77(2):544-548.

Modeling and Simulation Methods for Mitral Valve Based on the Mitral Valve Sub-region

ZHAO Xiao-jia1,SONG Xiao2,HUANG Xiao-yang1*,SU Mao-long3,ZHONG Qi1
(1.School of Information Science and Engineering,Xiamen University,Xiamen 361005,China;2.College of Computer and Information Technology,Nanyang Normal University,Nanyang 473061,China;3 Ultrasound Department, Zhongshan Hospital,Xiamen University,Xiamen 361004,China)

A fine model of the human mitral valve with non-uniform thickness has been established.The valve closure has been simulated using the finite element method during early systolic period,compared with the uniform one.Simulation results show that the stress is unevenly distributed on the leaflets,that the mitral valve with non-uniform has lower stress level than the uniform one.Furthermore,the maximum principal stresses occurs at the clefts between the three scallops in the posterior leaflet,indicating positions of the mitral valve cleft in clinic. The result is consistent with the physiology.It suggests that the method and the results are correct.Finally,based on clinical data,the motion of annulus and papillary muscles has been simulated,showing higher closure and relatively uniform stress distribution,closer to the actual physiological conditons.It laid the foundation for the simulation of subsequent pathology.

the mitral valve;finite element model;stress distribution;maximum principal stress

10.6043/j.issn.0438-0479.2015.02.021

TP 391

A

0438-0479(2015)02-0276-05

2014-08-28 錄用日期:2014-11-19

國家自然科學基金(61102137,61271336,61327001);福建省自然科學基金(2014J01438)

*通信作者:xyhuang@xmu.edu.cn

趙曉佳,宋曉,黃曉陽,等.基于亞區結構的人體二尖瓣的建模與仿真[J].廈門大學學報:自然科學版,2015,54(2): 276-280.

:Zhao Xiaojia,Song Xiao,Huang Xiaoyang,et al.Research on modeling and simulation methods for mitral valve based on the mitral valve sub-region[J].Journal of Xiamen University:Natural Science,2015,54(2):276-280.(in Chinese)