基于精細(xì)解剖結(jié)構(gòu)的左心室心肌缺血仿真

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

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

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

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

1 方法

1.1 心室組織模型

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

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

1.2 數(shù)值求解方法

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

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

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

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

引入輔助域φ來區(qū)分心室組織的內(nèi)外界，同時(shí)建立一個(gè)規(guī)則圖形將心室組織容納進(jìn)去，其中φ=0為在心室外部，φ=1為在心室內(nèi)部.

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

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

下面舉1個(gè)簡單的例子來解釋?duì)针S時(shí)間變化在不規(guī)則圖形上的分布情況.建立一個(gè)100×100的正方形，內(nèi)置一內(nèi)徑20，外徑40的圓環(huán)，在正方形內(nèi)對(duì)該圓環(huán)求解式(3)，當(dāng)t=500 ms時(shí)刻的相值分布如圖1所示.由圖1可見，圓環(huán)區(qū)域內(nèi)相值為1，在圓環(huán)邊界處相值光滑且連續(xù)地從1變?yōu)?，圓環(huán)邊界處外相值穩(wěn)定為0.這樣，利用相場法可以在不破壞計(jì)算精確性的同時(shí)直接使用規(guī)則圖形的邊界求解方法對(duì)復(fù)雜邊界問題進(jìn)行數(shù)值求解.由于使用了統(tǒng)一的控制方程來自動(dòng)處理邊界條件而無需跟蹤，因而大大降低了計(jì)算復(fù)雜度.

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

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

化簡式(5)可變?yōu)?/p>

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

1.3 數(shù)據(jù)預(yù)處理

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

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

2 結(jié)果

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

2.1 正常組織內(nèi)的折返波

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

圖3 正常左心室組織內(nèi)不同時(shí)刻的折返波形態(tài)

2.2 心肌缺血時(shí)的折返波

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

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

本文在左心室組織切片左上角設(shè)置一塊大小約為30×30的不規(guī)則區(qū)域?yàn)槿毖獏^(qū)域.如圖4所示，I區(qū)為缺血區(qū)域，II區(qū)為健康區(qū)域，I區(qū)和II區(qū)之間的連接區(qū)域III區(qū)為過渡區(qū).

實(shí)驗(yàn)發(fā)現(xiàn)由于心內(nèi)膜細(xì)胞的動(dòng)作電位持續(xù)時(shí)間太短，導(dǎo)致T波倒置，為了更接近真實(shí)情況，這里把其IKs通道的電導(dǎo)率GKs按照文獻(xiàn)［11］調(diào)整為0.149 nS/pF，延長其動(dòng)作電位時(shí)程.

圖4 缺血區(qū)域

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

在t=1 350 ms時(shí)，由圖5可以看到，折返波已經(jīng)到達(dá)缺血區(qū)域，波前出現(xiàn)不同程度的缺口，這是由于折返波在缺血區(qū)域內(nèi)的傳播速度減慢所致.隨著［K+］o升高，缺血區(qū)域?qū)φ鄯挡ǖ淖枞絹碓酱?局部輕度缺血如圖5(a)所示的情況下，波前只是出現(xiàn)很小的缺口，缺血程度越嚴(yán)重，所出現(xiàn)的缺口越明顯如圖5(b)，(c)所示.

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

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

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

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

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

對(duì)比上述仿真結(jié)果，不難看出，缺血情況越嚴(yán)重，對(duì)折返波的影響越大.當(dāng)缺血嚴(yán)重到一定程度的時(shí)候，折返波變得很不穩(wěn)定，最終很難穿過缺血區(qū)域，這與文獻(xiàn)［4］的研究結(jié)果一致.

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

3 結(jié)論

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

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

3)實(shí)驗(yàn)結(jié)果表明，心肌缺血情況下，在缺血區(qū)域，折返波的傳播受到影響.缺血程度越嚴(yán)重，折返波越不穩(wěn)定，甚至?xí)l(fā)生斷裂現(xiàn)象.

4)通過分析實(shí)驗(yàn)結(jié)果發(fā)現(xiàn)，即使折返波穿過了缺血區(qū)域，缺血組織的心肌細(xì)胞動(dòng)作電位也發(fā)生改變，導(dǎo)致動(dòng)作電位的空間不一致性，這是誘發(fā)心律失常的本質(zhì)原因.

［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］.生物醫(yī)學(xué)工程學(xué)雜志，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，F(xiàn)ENTON 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］朱浩，尹炳生，朱代謨.基于單細(xì)胞電位計(jì)算心電:若干異常仿真心電圖［J］.生物物理學(xué)報(bào)，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.