基于精細解剖結構的左心室心肌缺血仿真

呂偉剛，王寬全，左旺孟，黎 捷，張恒貴

(1.哈爾濱工業大學計算機科學與技術學院，150001哈爾濱，lwg-2001@163.com; 2.曼徹斯特大學物理與天文學院，曼徹斯特，英國)

近半個世紀以來，心臟病已成為威脅人類健康最嚴重的疾病之一.根據中華人民共和國衛生部最新統計，2008年我國城市居民心臟病的死亡率已由2007年的第3位(僅次于惡性腫瘤和腦血管病)上升至第2位(僅次于惡性腫瘤)，而且仍保持上升趨勢.其中心臟猝死造成的后果最為嚴重.大多數情況下，心臟猝死是由心律失常所導致［1］.室速和室顫是心律失常的兩大類癥狀，往往發生在心肌缺血期間［2］.

TNNP模型是Tusscher等［3］于2004年提出的人類心室肌細胞模型，該模型基于新的實驗數據，包括了所有主要的離子電流，在電生理特性和仿真性能上更接近人體的真實情況，因而更適用于復雜心室組織電生理活動的仿真研究.袁永峰等［4］研究了規則圖形上的心肌缺血仿真研究問題，但是該模型存在2點問題:1)基于規則模型的仿真結果可能難以準確反映折返波在真實人體心室組織的傳播過程;2)由于心室幾何形狀上的復雜性，因此規則圖形上的邊界求解過程不能直接應用在心室組織上.基于上述問題，本文引入精細的人體心臟解剖數據，同時考慮了心室肌細胞的電異質性，構建了一個整合的人體心室組織模型.在處理無通量邊界問題上，使用相場法來自動計算復雜的邊界條件［5-7］.

1 方法

1.1 心室組織模型

心肌組織是由許多單細胞互相電耦合而組成的一個電耦合胞體.在由心肌細胞組成的傳導系統中，動作電位是通過偏微分方程來描述的.為了描述心室組織的電興奮傳導過程，使用TNNP模型，表述為

式中:Vm為跨膜電位，t為時間，Cm為跨膜電容，σ為電導率，Sv為面積體積比，也可寫為， D為擴散系數，▽為梯度算子，Iion為跨膜離子流總和，Istim為所施加的刺激電流.

1.2 數值求解方法

心肌組織的電生理模型都需要滿足無通量邊界條件［8］，即n·D▽Vm=0，其中n為垂直于表面的法相矢量.其物理意義為沒有電流流入或流出外部環境.

對于規則圖形而言，在不考慮纖維走向各向異性的前提下，此時的邊界條件問題可以簡化為

式中n⊥i，n⊥j.由式(2)可以得到邊界點的跨膜電位Vm(x，y)=Vm(x-1，y)或Vm(x，y)= Vm(x，y-1)，即規則圖形邊界點的跨膜電位等于其內部沿x或y方向最近鄰點的跨膜電位值.

但是對于精細的人體心室組織來說，由于其幾何上的復雜性，計算邊界點的法相矢量n并求解邊界條件非常復雜且計算量大，在保證計算精確性的前提下就需要一個高效的方法來處理邊界問題.這里使用相場法來自動計算復雜的邊界條件.相場法是以金茲堡 -朗道相變理論為基礎，引入新變量——相場φ而得名.

引入輔助域φ來區分心室組織的內外界，同時建立一個規則圖形將心室組織容納進去，其中φ=0為在心室外部，φ=1為在心室內部.

通過求解式(3)來確定φ值為

式中ξ為控制界面寬度的參數.當ξ值足夠小時，無通量邊界條件被滿足，證明過程詳見文獻［7］. G(φ)是任意滿足在φ=0和φ=1有最小值的雙阱函數.這里，選擇G(φ)為

下面舉1個簡單的例子來解釋φ隨時間變化在不規則圖形上的分布情況.建立一個100×100的正方形，內置一內徑20，外徑40的圓環，在正方形內對該圓環求解式(3)，當t=500 ms時刻的相值分布如圖1所示.由圖1可見，圓環區域內相值為1，在圓環邊界處相值光滑且連續地從1變為0，圓環邊界處外相值穩定為0.這樣，利用相場法可以在不破壞計算精確性的同時直接使用規則圖形的邊界求解方法對復雜邊界問題進行數值求解.由于使用了統一的控制方程來自動處理邊界條件而無需跟蹤，因而大大降低了計算復雜度.

圖1 相值φ在圓環上的分布(ξ=0.25 mm，Δx=Δy= 0.25 mm，Δt=0.015 ms，t=500 ms)

同理，由于在心室組織之外，跨膜電容Cm，跨膜電流(Iion+Istim)以及電導率σ并不存在，所以分別修正為 φCm，φ(Iion+Istim)以及 φσ，代入式(1)，可得

化簡式(5)可變為

這里使用前向歐拉法求解式(6)，具體過程略.本文選擇空間步長Δx=Δy=0.33 mm，時間步長Δt=0.02 ms，界面控制參數ξ=0.33 mm.

1.3 數據預處理

本文以美國“可視人計劃”(Visible Human Project)解剖結構數據集為基礎，提取了女性切片數據中的左心室三維數據，沿Z軸方向選擇信息量較多的1片左心室組織.考慮到心室肌細胞的電異質性，按照距離內外邊界的比例從內向外將心室壁分為3層結構:心內膜層、中間層和心外膜層.其細胞比例為4∶3∶4，這樣就建立了人體左心室切片組織(如圖2).

圖2 分層后的人體左心室切片組織

2 結果

要研究心肌缺血下折返波的傳播狀況，首先就要激發出折返波.這里使用S1-S2協議，其中，S2刺激加在S1刺激不應期的尾部，S2和S1的強度2倍于閾值.

2.1 正常組織內的折返波

圖3為正常的左心室組織內不同時刻折返波的傳播過程.此時折返波沒有遇到任何障礙，勻速平衡的傳播.

圖3 正常左心室組織內不同時刻的折返波形態

2.2 心肌缺血時的折返波

當心臟處于異常狀態如心肌缺血的情況下，電興奮的傳播過程將受到影響.心肌缺血是心臟最常見的局部病理改變，其取決于缺血面積、部位、嚴重性等因素.當心肌細胞缺血時，發生的生理和病理變化有［9］:1)局部缺氧，ATP水平下降，酸性代謝產物積累，細胞內外環境酸化;2)因能量代謝障礙，膜離子泵功能下降，細胞外鉀離子濃度(［K+］o)增高;3)胞內累積的脂質代謝物堆積到細胞膜上，影響離子通道蛋白與離子載體蛋白的活性;4)細胞內離子濃度失衡，Ca2+和Mg2+升高.

本文僅從缺血嚴重程度上探討局部心肌缺血對折返波的影響，通過在缺血區域內增加［K+］o來仿真心肌缺血環境［10］.模擬的3種缺血情況分別為局部輕度缺血、局部中度缺血和局部嚴重缺血，對應的［K+］o值分別為7.5、9.0和11.0 mM，而正常情況下［K+］o為5.4 mM.

本文在左心室組織切片左上角設置一塊大小約為30×30的不規則區域為缺血區域.如圖4所示，I區為缺血區域，II區為健康區域，I區和II區之間的連接區域III區為過渡區.

實驗發現由于心內膜細胞的動作電位持續時間太短，導致T波倒置，為了更接近真實情況，這里把其IKs通道的電導率GKs按照文獻［11］調整為0.149 nS/pF，延長其動作電位時程.

圖4 缺血區域

圖5～圖7分別展示了在3種典型的時刻不同缺血情況下折返波的形態.

在t=1 350 ms時，由圖5可以看到，折返波已經到達缺血區域，波前出現不同程度的缺口，這是由于折返波在缺血區域內的傳播速度減慢所致.隨著［K+］o升高，缺血區域對折返波的阻塞越來越大.局部輕度缺血如圖5(a)所示的情況下，波前只是出現很小的缺口，缺血程度越嚴重，所出現的缺口越明顯如圖5(b)，(c)所示.

圖5 不同局部缺血情況下同一時刻折返波形態(t=1 350 ms)

在t=1 400 ms時，由圖6可以看到，局部輕度缺血對折返波影響不大，只是傳播速度略微減慢，之后仍能穿過缺血區域繼續向上傳播如圖6(a)所示.局部中度缺血的情況下，缺血區域內的波速明顯慢于周圍正常組織，但是折返波試圖穿過缺血區域繼續傳播如圖6(b)所示.局部嚴重缺血時，折返波已經很難穿過缺血區域，然而其周圍的正常組織內折返波繼續向上傳播，這樣折返波的傳播發生很大程度的偏轉如圖6(c)所示.

在t=1 600 ms時，由圖7可以看到，局部輕度缺血和中度缺血時，折返波已經完全穿過了缺血區域繼續傳播如圖7(a)，(b)所示.局部嚴重缺血時，折返波在缺血區域處發生斷裂，波前繞過缺血區域，匯合之后繼續傳播如圖7(c)所示.

圖6 不同局部缺血情況下同一時刻折返波形態(t=1 400 ms)

圖7 不同局部缺血情況下同一時刻折返波形態(t=1 600 ms)

對比上述仿真結果，不難看出，缺血情況越嚴重，對折返波的影響越大.當缺血嚴重到一定程度的時候，折返波變得很不穩定，最終很難穿過缺血區域，這與文獻［4］的研究結果一致.

同時，通過對正常組織和缺血組織的心肌細胞分析發現，缺血時期［K+］o的提高可以明顯地影響缺血組織心肌細胞動作電位時程(APD).缺血越嚴重，缺血組織心肌細胞APD縮短程度越大.出現這種現象主要是由于［K+］o的提高在某種程度上增大了內部校正和延遲校正通道的傳導速率，從而導致了APD的減小.即使折返波穿過了缺血區域，缺血組織心肌細胞的動作電位(AP)也已經發生了改變，這種改變使得動作電位空間不一致，引起各種離子以及代謝物的擴散.從電生理角度看，這些變化會改變心肌細胞動作電位的形態、興奮性、傳導速率等，進而導致心臟組織傳導阻塞，這極大的不利于折返波的穩定傳播，從而誘發心律失常.

3 結論

1)基于TNNP模型，考慮了心室肌細胞的電異質性，構建了人體左心室組織電生理模型，在處理無通量邊界問題上，引入了一種高效的算法——相場法來求解.

2)通過設置缺血區域，同時在該區域內增加［K+］o仿真了局部心肌缺血情況下折返波在心室組織的傳播過程.

3)實驗結果表明，心肌缺血情況下，在缺血區域，折返波的傳播受到影響.缺血程度越嚴重，折返波越不穩定，甚至會發生斷裂現象.

4)通過分析實驗結果發現，即使折返波穿過了缺血區域，缺血組織的心肌細胞動作電位也發生改變，導致動作電位的空間不一致性，這是誘發心律失常的本質原因.

［1］Ten TUSSCHER K H，BERNUS O，HREN R，et al. Comparison of electrophysiological models for human ventricular cells and tissues［J］.Progress in Biophysics and Molecular Biology，2006，90(1/3):326-345.

［2］HEIDENREICH E A，RODRIGUEZ J F，DOBLARE M，et al.Electrical propagation patterns in a 3D regionally ischemic human heart.a simulation study［C］// Proceedings of Computers in Cardiology 2009.Park City:IEEE，2009:665-668.

［3］Ten TUSSCHER K H，NOBLE D，NOBLE P J，et al.A model for human ventricular tissue［J］.American Journal of Physiology-Heart and Circulatory Physiology，2004，286(4):1573-1589.

［4］袁永峰，王寬全，田慧麗.二維人體心室心肌缺血模型中的折返波仿真研究［J］.生物醫學工程學雜志，2009，26(6):1329-1334.

［5］KARMA A，RAPPEL W J.Quantitative phase-field modeling of dendritic growth in two and three dimensions［J］.Physical Review E，1998，57(4):4323-4349.

［6］FENTON F H，CHERRY E M，KARMA A，et al.Modeling wave propagation in realistic heart geometries using the phase-field method［J］.Chaos，2005，15(1):013502.

［7］BUENO-OROVIO A，PEREZ-GARCIA V M，FENTON F H.Spectral methods for partial differential equations in irregular domains:the spectral smoothed boundary method［J］.SIAM Journal on Scientific Computing，2006，28(3):886-900.

［8］CLAYTON R H，PANFILOV A V.A guide to modeling cardiac electrical activity in anatomically detailed ventricles［J］.Progress in Biophysics and Molecular Biology，2008，96(1/3):19-43.

［9］朱浩，尹炳生，朱代謨.基于單細胞電位計算心電:若干異常仿真心電圖［J］.生物物理學報，2001，17(1):123-134.

［10］XU A，GUEVARA M R.Two forms of spiral-wave reentry in an ionic model of ischemic ventricular myocardium［J］.Chaos，1998，8(1):157-174.

［11］ZHANG H，HANCOX J C.In silico study of action potential and QT interval shortening due to loss of inactivation of the cardiac rapid delayed rectifier potassium current［J］.Biochemical and Biophysical Research Communications，2004，322(2):693-699.