許煒煒 白明珠 林強(qiáng) 胡正琿?

1) (浙江工業(yè)大學(xué)理學(xué)院,杭州 310023)2) (浙江工業(yè)大學(xué),生物醫(yī)學(xué)物理協(xié)同創(chuàng)新中心,杭州 310023)

1 引 言

心磁正問題是通過構(gòu)建心臟電生理和體表磁場模型,描述心臟電興奮起搏點(diǎn)激發(fā)的心臟電生理活動,并在體表投影磁場的研究過程.這是深入理解心臟體表磁場形成機(jī)制的重要手段,使得心磁信號的生成過程有了理論指導(dǎo).心磁圖(magnetocardiography,MCG)作為一種新型無創(chuàng)檢測技術(shù),通過記錄心臟電生理活動產(chǎn)生體表外的磁場分布變化,來反映心動周期的過程.由于人體組織磁導(dǎo)率接近真空磁導(dǎo)率的特性,MCG不受組織和空間的影響,相比傳統(tǒng)多導(dǎo)聯(lián)心電圖(ECG)可捕獲更微弱的生物信息,對冠心病及心律失常表現(xiàn)出更高的敏感度和準(zhǔn)確度[1,2].而ECG受組織電導(dǎo)率等物理因素制約,抑制了一些不正常的電信號在體表的表現(xiàn)[3].因此我們認(rèn)為研究心臟電生理活動產(chǎn)生的體表外磁場分布的心磁正問題是必要的.另一方面,可以在完成心磁正問題的前提下,從實(shí)測MCG出發(fā),反演重建心臟等效電源或跨膜電位(transmembrane potential,TMP)分布,研究心臟電生理活動過程的心磁逆問題.這或許可以為定位心臟異常信號源頭提供幫助[4].然而心磁信號極為微弱,準(zhǔn)確測量心磁信號比較困難.傳統(tǒng)基于超導(dǎo)量子干涉器件(superconducting quantum interference device,SQUID)的磁力儀制造成本和維護(hù)費(fèi)用高昂,不適于MCG臨床推廣.但是,近年來原子磁力儀技術(shù)發(fā)展迅速,在測量極弱磁場方面展現(xiàn)了優(yōu)異的性能[5],并且原子磁力儀價格大約只有SQUID磁力儀的1/10.何祥等[6]于2017年報道了一種基于非線性磁光旋轉(zhuǎn)效應(yīng)的脈沖泵浦式銣原子磁力儀,并在常溫下清晰地測得了心臟磁場信號,這是國內(nèi)首次用原子磁力儀實(shí)現(xiàn)對心臟磁場信號的探測.隨著基于原子磁力儀技術(shù)的心磁圖儀的發(fā)展和實(shí)用化,MCG將有希望成為臨床上與ECG互為補(bǔ)充的心臟電生理活動檢測技術(shù).心磁正問題的研究可以更好地利用原子磁力儀測得的體表外心磁信號,描述心臟的電生理活動,并為后續(xù)心磁逆問題提供驗(yàn)證.這也為我們研究心磁正問題提供了動力.

在過去幾十年里,隨著生物醫(yī)學(xué)知識和計(jì)算機(jī)技術(shù)的發(fā)展,研究者從細(xì)胞離子通道到組織層面的不同尺度研究并開發(fā)了多種心臟模型.這些模型不僅可以用于更加深入地了解心律失常等疾病的機(jī)理,還可以用于設(shè)計(jì)新的治療方法來治療心臟疾病等.1952年,Hodgkin和Huxley[7]對魷魚軸突的刺激反應(yīng)進(jìn)行了研究,得出了有關(guān)動作電位的復(fù)雜關(guān)系式,這為之后的心臟電生理學(xué)模型的研究奠定了基礎(chǔ)[8-10].最早的三維心臟模型由奧克蘭大學(xué)Nielsen等[11]在1991年基于解剖犬心臟,通過組織學(xué)分析建立,這為后續(xù)研究者創(chuàng)建心臟計(jì)算模型提供了三維幾何基礎(chǔ).之后,各國發(fā)起“可視人計(jì)劃”[12],通過人體解剖切片獲得高精度的人體組織數(shù)據(jù).Xia等[13]在可視人數(shù)據(jù)基礎(chǔ)上,從心臟細(xì)胞模型出發(fā)建立了代表國際先進(jìn)水平的Cardiome-CN虛擬心臟電生理數(shù)學(xué)模型,并將其應(yīng)用于各類心臟疾病的研究.李心雅等[14]基于波面型仿真算法完成了全心的電生理建模和體表電位標(biāo)測圖(body surface potential mapping,BSPM)仿真,降低了仿真所需的計(jì)算量.而醫(yī)學(xué)斷層成像技術(shù)的普及,為研究人員無創(chuàng)獲取個性化的人體心臟幾何數(shù)據(jù)提供了一種便捷的方式.Wang等[4]使用磁共振影像(MRI)數(shù)據(jù),結(jié)合心電正問題和逆問題,構(gòu)建了個性化的三維心肌梗死電生理學(xué)模型.隨著計(jì)算機(jī)技術(shù)的不斷發(fā)展,越來越精細(xì)的心臟電生理學(xué)模型開始建立并用于正常心臟狀態(tài)以及各類心臟疾病的仿真.

本研究將建立一個基于三維個性化幾何的心磁正問題計(jì)算框架,加深對心臟電生理活動演變,以及在軀干表面形成生物磁場的物理過程的理解.本文將對被試的MRI進(jìn)行圖像分割,分離出心臟各腔室與軀干表面,之后對處理好的圖片進(jìn)行三維重建,建立一個個性化的人體心臟和軀干三維幾何模型.在該三維幾何模型的基礎(chǔ)上,用修正后的FitzHugh-Nagumo (FHN)方程建立心臟電生理擴(kuò)散模型,以求得TMP變化和電興奮在心臟內(nèi)的傳播過程.之后,以準(zhǔn)靜態(tài)麥克斯韋方程組為基礎(chǔ)建立體表外心臟磁場模型,模擬由TMP演變產(chǎn)生的心臟磁場在身體內(nèi)和空氣中的傳播過程,求得體表外磁場投影分布.本研究的計(jì)算過程將基于伽遼金法(Galerkin method)的有限元分析(finite element analysis,FEA).最后,將通過簡化模型,對其解析解與數(shù)值解做誤差分析,驗(yàn)證心磁計(jì)算框架的可行性.

2 方 法

2.1 心臟軀干三維幾何模型的構(gòu)建

本研究使用MRI影像數(shù)據(jù),通過圖像分割與三維重建,建立一個個性化心臟軀干三維幾何模型.傳統(tǒng)三維心臟數(shù)據(jù)源來自犬等解剖模型,與真實(shí)人體幾何存在較大差異.另一方面,侵入式的數(shù)據(jù)采集方法不利于獲取人體心臟數(shù)據(jù),因而存在較大局限性.而人體醫(yī)學(xué)斷層影像(如CT,MRI等)技術(shù)可以無創(chuàng)便捷地從被試者獲取醫(yī)學(xué)影像.通過對醫(yī)學(xué)影像進(jìn)行圖像分割與三維重建,即可獲得研究所需的三維模型.該模型可以在更貼近真實(shí)生理狀態(tài)的情況下,模擬心臟TMP傳播以及由此產(chǎn)生的磁場體表外的分布.本研究根據(jù)醫(yī)學(xué)圖像中的灰度值范圍,通過將二值化后的心臟和軀干MRI影像進(jìn)行閾值分割,得到目標(biāo)的基本輪廓,最后將處理好的圖像進(jìn)行三維圖像重建[15,16].在不影響心臟和軀干結(jié)構(gòu)的前提下,對重建完成的三維幾何進(jìn)行適當(dāng)平滑處理,消除圖像分割過程中由于偏差等造成的表面不平滑的狀況.

本研究使用一位26歲健康女性志愿者的平掃M(jìn)RI數(shù)據(jù),構(gòu)建三維幾何模型.其中,心臟圖像由18張像素為 256×256 ,分辨率為1 mm×1 mm斷層分辨率為10 mm的MRI切片構(gòu)成; 軀干數(shù)據(jù)由60張像素為 768×504 ,分辨率為1 mm×1 mm,斷層分辨率為4 mm的MRI切片構(gòu)成.使用最大類間方差法[17]對心臟和軀干圖像切片依次進(jìn)行圖像分割,結(jié)合手繪處理的方式,分離出心臟中心室和心房的內(nèi)外輪廓和軀干表面輪廓,如圖1所示.

對分割完成的圖像進(jìn)行體繪制三維重建[15]并利用移動平均濾波器進(jìn)行平滑處理[12],可以最終得到個性化的心臟和軀干三維結(jié)構(gòu).重建完成后的心臟結(jié)構(gòu)包括左右心房以及左右心室共四個腔室.與仿真體表處心臟電勢分布不同的是,心磁的計(jì)算需要考慮體外區(qū)域的磁場分布.本文構(gòu)建了一個半徑為2 m的球形區(qū)域?qū)⑿呐K-軀干模型包裹,以模擬人體周圍的空氣,如圖2所示.構(gòu)建完成后的幾何模型將被應(yīng)用到心臟電生理活動和體表外心臟磁場分布的仿真計(jì)算.

圖1 MRI切片的圖像分割 (a) 心臟; (b) 軀干Fig.1.Image segmentation of MRI slices: (a) Heart; (b) torso.

圖2 個性化心臟-軀干三維幾何模型Fig.2.Personalized three-dimensional heart-torso model.

2.2 心臟電生理擴(kuò)散模型

根據(jù)心臟TMP演變的特點(diǎn),構(gòu)建能較好模擬興奮傳導(dǎo)且計(jì)算代價較小的電生理擴(kuò)散模型.心臟去極化和復(fù)極化過程中,心肌細(xì)胞膜內(nèi)外離子發(fā)生轉(zhuǎn)移,TMP發(fā)生變化.不同位置的細(xì)胞TMP變化整體上表現(xiàn)為宏觀層面的心臟電興奮傳導(dǎo).FHN模型[18]是一種定性研究心臟電生理學(xué)活動的反應(yīng)擴(kuò)散系統(tǒng)經(jīng)典模型,最早由Hodgkin和Huxley[7]通過簡化多變量的Hodgkin-Huxley細(xì)胞模型得到.與元胞自動機(jī)模型模擬電興奮類似,FHN模型并不直接模擬單個細(xì)胞的興奮特性和細(xì)胞之間的偶聯(lián)關(guān)系.由于不必大規(guī)模計(jì)算心臟細(xì)胞興奮特性,因此可用較小的計(jì)算代價仿真心臟電生理活動的TMP動態(tài)特性.而基于簡單規(guī)則的元胞自動機(jī)模型對解釋電興奮傳導(dǎo)的物理過程不利[19].FHN模型包含描述動作電位產(chǎn)生的細(xì)胞模型和心肌興奮傳導(dǎo)的擴(kuò)散模型,較符合實(shí)際的電興奮傳導(dǎo)特征.FHN模型是由兩個未知變量的非線性偏微分方程組成的反應(yīng)擴(kuò)散系統(tǒng):

其中u表示激發(fā)變量,變化范圍從0到1,ν表示恢復(fù)變量,?·(D?u) 表示擴(kuò)散項(xiàng),D表示擴(kuò)散張量,而f1(u,v) 和f2(u,v) 則定義了動作電位的形狀.

FHN模型存在較多變式,為了更為準(zhǔn)確地表示真實(shí)心臟TMP變化(如去除原始FHN模型中實(shí)際不存在的過極化現(xiàn)象),本研究采用了Aliev等[20]修正的FHN模型

其中參數(shù)a=0.15 ,e=0.01 ,k=8[21],研究初期忽略心肌纖維的各向異性,即假設(shè)擴(kuò)散張量D為單位矩陣.實(shí)際的心臟TMP(Vm)范圍是-80-20 mV,可以表示為

通常認(rèn)為心臟作為一個孤立的連續(xù)體,即沒有有功電流流入或流出心臟表面,所以存在紐曼邊界條件?u/?n=0,其中n表示心臟邊界的法向量.

2.3 體表外心臟磁場模型

根據(jù)電磁學(xué)理論,計(jì)算心臟電興奮產(chǎn)生磁場在體內(nèi)的傳播過程.人體心臟電磁場頻率大約在1-100 Hz之間,對這種頻率的低頻電磁場通常使用準(zhǔn)靜態(tài)電磁學(xué)的知識處理[22].

在宏觀表現(xiàn)為興奮傳導(dǎo)的TMP變化以電流密度的形式,生成磁場在軀干以及體外傳播.心磁信號就是心臟磁場在體外的投影積分.心臟中外電流密度Ji與TMP之間滿足方程:

其中σH表示心臟組織的電導(dǎo)率(下標(biāo)H是heart的縮寫,下同).由于心臟和軀干被視為容積導(dǎo)體,模型中的全電流密度J為外電流密度Ji與無源容積電流密度之和:

其中σ表示電導(dǎo)率,φ表示電勢.電流密度J滿足電流守恒定律[8]

心臟TMP轉(zhuǎn)化為電流密度的變量Ji后,將其應(yīng)用到準(zhǔn)靜態(tài)麥克斯韋方程.由于心磁問題滿足磁準(zhǔn)靜態(tài)場,位移電流JD遠(yuǎn)小于傳導(dǎo)電流,即忽略電場變化率引起的磁場變化磁感應(yīng)強(qiáng)度B滿足

其中μ是磁導(dǎo)率.由于人體組織的無磁性,即相對磁導(dǎo)率μr=1 ,μ=μ0μr=μ0.矢量磁位A滿足B=?×A,并滿足庫侖規(guī)范,得到

由此可得,心臟區(qū)域?H和軀干區(qū)域?T矢量磁位分別滿足:

其中σT表示軀干電導(dǎo)率.由于體外空氣域?A的電導(dǎo)率非常小,忽略體外電流,體表外區(qū)域矢量磁位滿足

根據(jù)電磁學(xué)理論可知,磁場傳播在心臟表面ΓH和體表ΓT的邊界條件都分別滿足,磁感應(yīng)強(qiáng)度法向量方向相同B1n=B2n,磁場切線方向相同H1t=H2t.

2.4 模型的有限元仿真計(jì)算

分別對心臟電生理擴(kuò)散模型和體表外心臟磁場模型進(jìn)行數(shù)值計(jì)算,最終合成完整的心磁正問題計(jì)算框架.求解復(fù)雜邊界下的偏微分方程和物理場問題的數(shù)值解通常包括邊界元法(boundary element method,BEM)[23]、有限元法(finite element method,FEM)[24]、無網(wǎng)格法(meshfree method)[25]等.本研究采用了FEM對心磁正問題進(jìn)行計(jì)算.由于使用等效積分弱形式的方式相對于變分法使用范圍更廣泛,能較好地應(yīng)對偏微分方程和電磁學(xué)問題.另外,伽遼金法在加權(quán)余數(shù)法中擁有更加優(yōu)異的精度[26],因此本研究采用了伽遼金法對心臟電生理擴(kuò)散模型和體表外心臟磁場模型進(jìn)行數(shù)值計(jì)算.

本文采用二階10節(jié)點(diǎn)的四面體單元,對心臟、軀干以及體外空氣域進(jìn)行離散化處理,以便進(jìn)行有限元計(jì)算,如圖3所示.為準(zhǔn)確仿真心臟電生理活動和體表外磁場分布,并考慮計(jì)算成本,采取相對密集的心臟四面體網(wǎng)格與相對稀疏的軀干四面體網(wǎng)格,如圖4所示.心臟部分包括9723個四面體單元,軀干部分包括11698個四面體單元,而空氣域部分包括10762個四面體單元(圖中未顯示空氣域).

圖3 二階10節(jié)點(diǎn)四面體單元Fig.3.The second-order tetrahedral element with 10 nodes.

圖4 離散的心臟-軀干三維模型 (a) 心臟三維模型;(b) 心臟-軀干模型組合Fig.4.Discretized three-dimensional heart-torso model:(a) 3D cardiac model; (b) heart-torso model.

本研究采用的單元形函數(shù)類型為二階拉格朗日單元[27].這種單元類型可以更好模擬彎曲的邊界,對擁有不規(guī)則幾何的心臟和軀干模型具有較好的計(jì)算精度.伽遼金法的特點(diǎn)是選取形函數(shù)作為權(quán)函數(shù),最終獲得含待定系數(shù)的有限元方程,通過計(jì)算有限元方程求得相應(yīng)數(shù)值解.

本研究首先求得心臟電生理模型中不同坐標(biāo)和時間的TMP分布.之后將TMP轉(zhuǎn)化為電流密度參數(shù),應(yīng)用到體表外心臟磁場模型,最終求得在體表外磁場的分布情況.

2.5 心磁計(jì)算框架的解析解驗(yàn)證

本研究分別設(shè)計(jì)了一維狀態(tài)的反應(yīng)擴(kuò)散模型和三維簡化幾何的磁場模型,通過求解上述簡化模型的精確解,并與相應(yīng)FEM數(shù)值解進(jìn)行誤差分析,最終驗(yàn)證心磁正問題計(jì)算框架數(shù)值解的可信度.

本研究對FHN模型的有限元數(shù)值解的準(zhǔn)確性進(jìn)行驗(yàn)證,以保證心臟電生理擴(kuò)散模型的準(zhǔn)確性.為了能簡便求得FHN模型解析解,本研究考慮一維條件下的FHN方程,并與對應(yīng)條件下的伽遼金有限元數(shù)值解進(jìn)行對比.由于FHN方程由兩個相互耦合的非線性偏微分方程構(gòu)成,因此對方程組的求解具有困難.本文考慮適當(dāng)簡化FHN方程組,再對方程進(jìn)行求解.

以原始FHN方程(1)為例,由于在心臟電生理模型中,參數(shù)b滿足 0＜b＜＜1 ,因此可以認(rèn)為因變量v對時間t的偏導(dǎo)為0,即v為常數(shù),進(jìn)一步假設(shè)v為0.FHN方程組可以簡化為一個非線性反應(yīng)擴(kuò)散方程

較多文獻(xiàn)對反應(yīng)擴(kuò)散方程(12)式有研究,如通過李群法[28]、Tanh法[29]、變系數(shù)Bernoulli輔助法[30]等求得反應(yīng)擴(kuò)散方程的多個解析解.

本文參考Tanh法[29],對反應(yīng)擴(kuò)散方程進(jìn)行精確解計(jì)算.假設(shè)方程(12)存在行波解,并設(shè)u(x,t)=u(ξ),其中ξ=k(x-ct) ,將方程轉(zhuǎn)化為u關(guān)于ξ的常微分方程.可以通過構(gòu)造Y=tanh2ξ,變換根據(jù)齊次平衡原則得到關(guān)于Y的齊次方程.根據(jù)方程系數(shù)都為0的特性,可以求得若干行波解,具體過程本文不再詳述.可以將精確解與對應(yīng)數(shù)值解進(jìn)行誤差分析,驗(yàn)證簡化的心臟電生理學(xué)模型的精確性.

另一方面,為了驗(yàn)證體表外磁場數(shù)值解的精確性,本研究建立了一個簡易幾何的測試模型,并在測試模型下求解基于麥克斯韋方程組的磁場解析解.

在BSPM的仿真中,Fischer等[24]采用了同心球模型進(jìn)行解析驗(yàn)證.即用同心內(nèi)球的電勢分布表示心臟TMP分布情況,通過求解外球電勢模擬軀干表面的電勢分布.本文首先假定內(nèi)球中存在均勻并垂直向上的電流密度Ji.模型內(nèi)各點(diǎn)磁感應(yīng)強(qiáng)度已由Geselowitz給出[31]:

其中R為場點(diǎn)到源點(diǎn)的向量,σi表示各區(qū)域的電導(dǎo)率.直接求解內(nèi)球電流產(chǎn)生磁場的分布情況存在困難.由于對稱性,柱坐標(biāo)系下內(nèi)球電流所產(chǎn)生的磁場不存在r和z方向的磁場分量,即磁場方向?yàn)榉轿唤铅辗较?這與有限長直導(dǎo)線產(chǎn)生的磁場方向一致.

假設(shè)長度L,電流為I的直導(dǎo)線等價貢獻(xiàn)了內(nèi)球電流產(chǎn)生的磁場.并假設(shè)外球?qū)щ娐师襎遠(yuǎn)小于1,即忽略(13)式等號右邊的第二項(xiàng),求解磁場在半徑為R的外球中的分布.簡化后的直導(dǎo)線模型產(chǎn)生的磁感應(yīng)強(qiáng)度B為

其中R為場點(diǎn)r=(r,φ,z) 到源點(diǎn)r′=(0,0,z′) 的向量R=r-r′=rer+(z-z′)ez.將(14)式在柱坐標(biāo)系下進(jìn)行積分計(jì)算:

簡化(15)式,可以求得

(16)式表示有限長直導(dǎo)線模型在大球內(nèi)的磁場分布的解析解,之后將與對應(yīng)數(shù)值解進(jìn)行誤差分析.

3 仿真結(jié)果

3.1 心臟電生理擴(kuò)散模型的驗(yàn)證和結(jié)果

通過Tanh法,求得了一維條件下簡化的FHN方程的15個行波解,本文以其中一個解為例

將t=0時(17)式求得的u(x,0) 作為(12)式的初值條件,求解在一維條件下的有限元數(shù)值解.分別求得t=20,40,60三個時刻x在0-40范圍內(nèi)各134個數(shù)據(jù)點(diǎn)的數(shù)值解.通過與相同時刻、坐標(biāo)下的解析解進(jìn)行對比,驗(yàn)證簡化后的FHN方程在一維條件下用伽遼金法求解的精確性.定義相對均方根誤差(relative root mean squared error,RRMSE)

其中an為模型解析解,bn為模型數(shù)值解,N為數(shù)值解的所有時間步長或計(jì)算節(jié)點(diǎn).

圖5顯示了FHN方程的解析解與數(shù)值解的激發(fā)電位曲線圖.求得在t=20,40,60 時,簡化的FHN方程解析解與對應(yīng)數(shù)值解之間的 RRMSE 分別為0.39%,0.53%和0.79%,如表1所列.結(jié)果顯示,在一維條件下用伽遼金法求得簡化后的FHN方程的數(shù)值解與Tanh法求得的解析解的誤差較小,展現(xiàn)了較好的準(zhǔn)確性.可以相信本研究通過伽遼金法求解的心臟電生理擴(kuò)散模型的結(jié)果是有保障的.

圖5 t=20,40,60時,激發(fā)電位u的變化 實(shí)線:FHN方程數(shù)值解; 虛線: FHN方程行波解Fig.5.Evolution of excitation potential u at t=20,40,60.Solid line: numerical solution of FHN equation.Dotted line:analytical solution of FHN equation.

在求解心臟電生理擴(kuò)散模型前,需要假定修改的FHN方程初值條件.由于正常功能的心臟是竇性心律的,即竇房結(jié)作為心臟正常起搏點(diǎn),產(chǎn)生激發(fā)電位,并形成TMP的傳播.因此可以設(shè)定在竇房結(jié)位置導(dǎo)入一個u=1 的激發(fā)電位,作為FHN方程的初值條件模擬TMP興奮的激發(fā).之后通過有限元法計(jì)算,可得到心臟TMP隨時間t變化的動態(tài)分布.

圖6所示為在一個完整的心動周期中,心室的TMP分布的演變情況.其中圖6(a)-圖6(c)表示心室去極化過程,對應(yīng)時刻分別為t=30,60,90.而圖6(d)-圖6(f)表示心室復(fù)極化過程,對應(yīng)時刻分別為t=120,150,180.分別取左右心室的心外膜上一點(diǎn),觀察TMP隨時間的變化,如圖7所示.其中實(shí)線表示右心室心外膜上一點(diǎn)的TMP變化,虛線表示左心室心外膜上一點(diǎn)的TMP變化.這與文獻(xiàn)中正常心率的心室跨膜電位波形一致[32],表明了以修正的FHN方程為基礎(chǔ)的心臟電生理擴(kuò)散模型是符合真實(shí)心臟電生理活動特性的.

3.2 體表外磁場的驗(yàn)證和仿真結(jié)果

如圖8(a)顯示的是在直導(dǎo)線模型中,中軸面上φ方向磁感應(yīng)強(qiáng)度Bφ的分布情況.圖8(b)顯示了Bφ的數(shù)值解與解析解的曲線圖,其中仿真過程中直導(dǎo)線的線徑r=0.1m 無法忽略.使用相對誤差(relative error,RE)求解中軸面上各坐標(biāo)點(diǎn)的磁場數(shù)值解與解析解之間的誤差

圖6 一個心動周期的TMP分布圖 (a),(b),(c) 為去極化過程 (t=30,60,90); (d),(e),(f)為復(fù)極化過程(t=120,150,180)Fig.6.TMP distribution of a cardiac cycle: (a),(b),(c) Depolarization process (t=30,60,90); (d),(e),(f) repolarization process (t=120,150,180).

圖7 TMP隨時間變化曲線 實(shí)線: 右心室TMP變化;虛線: 左心室TMP變化Fig.7.TMP time evolution curve,solid line: right ventricular TMP; dashed line: left ventricular TMP.

圖8 (a) 直導(dǎo)線模型; (b)中軸面磁感應(yīng)強(qiáng)度的分布,實(shí)線: 解析解; 虛線: 數(shù)值解Fig.8.(a) Straight wire model; (b) distribution of magnetic field on the axial plane.Solid line: analytical solution;dotted line: numerical solution.

其中an為模型解析解,bn為模型數(shù)值解.圖9顯示了中軸面上磁感應(yīng)強(qiáng)度Bφ解的相對誤差RE分布.其中散點(diǎn)表示中軸面上離導(dǎo)線中心各距離所對應(yīng)的誤差分布.而RE曲線由基于最小二乘法的二次多項(xiàng)式擬合得到.從圖9中可以看出,隨著離直導(dǎo)線中心距離的增加,RE也隨之變大.0.1≤r≤1.0m時,RE 最大為2.85%.對比心臟到體表處的距離,我們認(rèn)為這個誤差是可以接受的.經(jīng)過直導(dǎo)線模型的計(jì)算驗(yàn)證,說明構(gòu)建的伽遼金法求解體表外心臟磁場模型是可行的.

圖9 磁場的相對誤差分布Fig.9.Relative error distribution of magnetic field.

在心臟電磁場正問題中,為較好地模擬心臟產(chǎn)生的磁場在體表外的分布,通常需要設(shè)置各向異性的心臟電導(dǎo)率,而由于軀干電導(dǎo)率對心磁傳播影響較小,通常設(shè)置為各向同性的軀干電導(dǎo)率.但由于被試的MRI數(shù)據(jù)不包含心肌纖維取向,因此本研究建立的個性化三維心臟-軀干幾何模型前期并未考慮心肌細(xì)胞的纖維取向,并設(shè)定心臟電導(dǎo)率σH=0.48S/m,軀干電導(dǎo)率σT=0.2 S/m[21].

本研究在胸腔前方約5 mm處模擬設(shè)置一個探測傳感器,仿真計(jì)算在該平面的磁場分布情況,即求解體外該平面處的磁場數(shù)值解.通過將心臟電生理擴(kuò)散模型隨時間變化的結(jié)果應(yīng)用到靜態(tài)的體表外心臟磁場模型,即可求得該探測平面的不同時刻的磁感應(yīng)強(qiáng)度分布狀況.

圖10(a)-圖10(h),顯示了本研究仿真的垂直于探測平面體表外磁感應(yīng)強(qiáng)度Bz在一個心動周期內(nèi)的分布狀況.從圖中可以看出,仿真的心臟磁場在垂直于探測面方向存在兩極對稱的形態(tài).在同一個心動周期中,隨著心臟從去極化向復(fù)極化過程的演變,磁場兩極相對位置發(fā)生了變化.而圖10(i)是用9通道心磁圖儀采集人體心磁信號,并通過線性插值得到的健康人體表心磁圖[33].該心磁地圖與心臟功能異常的心磁圖在進(jìn)行了對比,顯示正常心磁地圖在形態(tài)上擁有規(guī)則的兩極形態(tài)特征.

對比圖10(a)-圖10(h)與圖10(i),看到本研究仿真模擬得到的心磁圖與實(shí)測心磁圖在形態(tài)上是類似的,因此可以定性地表明本研究建立的計(jì)算框架是可行的.

圖10 體表外磁感應(yīng)強(qiáng)度 Bz 分布 (a)-(h) 模擬的心磁分布圖(t=40,60,80,100,120,140,160,180)；(i) 文獻(xiàn)中的實(shí)測心磁圖分布圖[26]Fig.10.Extracorporeal magnetic field distribution: (a)-(h) Simulated cardiac magnetic distribution (t=40,60,80,100,120,140,160,180); (i) measured human MCG in the literature.

4 討論和結(jié)論

本研究構(gòu)建的幾何模型數(shù)據(jù)來自被試的MRI切片數(shù)據(jù).由于MRI切片數(shù)據(jù)的高精度性以及需求的可定制性,為更好地用于分割重建真實(shí)并且各種分辨率的三維心臟-軀干幾何模型提供了基礎(chǔ).并且為后續(xù)定制化仿真被試的心磁分布與實(shí)測心磁圖的對比提供了良好的條件.

本研究采用修正的FHN方程組作為心臟電生理擴(kuò)散模型的基礎(chǔ),模擬的心臟跨膜電位傳播過程符合心臟電生理實(shí)際狀況,對深入理解心臟電生理活動過程較為有利.另一方面,反應(yīng)擴(kuò)散模型以較小的計(jì)算代價仿真心臟電生理活動過程,展現(xiàn)了一定的優(yōu)勢.

本研究對心磁正問題計(jì)算框架進(jìn)行了解析解驗(yàn)證.首先,在一維條件下求解簡化的FHN方程解析解,與相同情況方程伽遼金有限元法數(shù)值解對比.其次,簡化體表外心臟磁場模型為直導(dǎo)線模型,求解模型磁感應(yīng)強(qiáng)度解析解與有限元解對比.通過對上述解的誤差分析,證明了該計(jì)算框架計(jì)算心臟電生理活動和體表磁場的可行性.

該計(jì)算框架有較好的擴(kuò)展能力.通過替換所需的心臟-軀干幾何模型,可以個性化仿真不同被試者的心臟電生理活動和體表心磁信號.另一方面,可以通過修改FHN方程組參數(shù),甚至修改FHN方程組形式,模擬不同狀態(tài)或者心律失常等疾病狀況下心臟跨膜電位不同或異常波形產(chǎn)生的異常磁場分布,并加以總結(jié)和分析.此外,通過設(shè)置不同位置的激發(fā)電位模擬心臟異常起搏點(diǎn)位置,模擬非竇性心律情況的心磁分布狀況.也可以通過抑制心臟中某些位置的興奮傳導(dǎo),模擬左束支傳導(dǎo)阻滯(left bundle branch block,LBBB)或右束支傳導(dǎo)阻滯(right bundle branch block,RBBB)等異常傳導(dǎo)狀態(tài)下心臟電生理活動和和心磁信號特征.該計(jì)算框架在這些情況都保留了優(yōu)異的擴(kuò)展能力.

但是,本研究仍然存在一些不足.首先,研究前期的三維心臟幾何模型較為簡陋,雖然采用了人體真實(shí)三維心臟作為模型對象,但是未準(zhǔn)確分割出包括竇房結(jié)具體形態(tài)在內(nèi)的各部分結(jié)構(gòu)細(xì)節(jié),只是在竇房結(jié)對應(yīng)位置設(shè)置了興奮激發(fā)的位置.其次,忽略了重要的心臟內(nèi)外膜與心肌細(xì)胞的纖維取向擴(kuò)散張量,但實(shí)際情況中TMP傳播過程受心臟各向異性結(jié)構(gòu)的影響較大,使得心臟電生理模型與實(shí)際情況存在較大差異.傳統(tǒng)的心臟纖維取向的設(shè)定依靠人為定義,這種方式與心臟實(shí)際的心肌旋向存在差異.而擴(kuò)散張量磁共振成像(DTMRI)等技術(shù)可以映射心臟幾何模型的真實(shí)纖維結(jié)構(gòu)[34],使得心臟電興奮傳播過程能受到纖維各向異性的電導(dǎo)率張量影響.為了心臟三維幾何模型更貼近實(shí)際情況,后期應(yīng)該將心臟纖維數(shù)據(jù)加入模型,以提高心磁正問題計(jì)算框架的精確性.另外,本研究只考慮了靜止空間狀態(tài)下的心磁正問題,并未討論心臟收縮和舒張對心臟電生理活動和體表磁場帶來的影響.因此,本文只定性地將求得的心磁分布圖與實(shí)際的心磁圖進(jìn)行比較,用以驗(yàn)證本研究計(jì)算框架的可行性.但是該心磁計(jì)算框架較好的擴(kuò)展能力為進(jìn)一步優(yōu)化心臟幾何結(jié)構(gòu)等提供了基礎(chǔ).

由于MCG的優(yōu)勢,以及開始在臨床用于互補(bǔ)心電圖的前景,本研究構(gòu)建了一個心磁正向問題計(jì)算框架.該框架包括建立一個個性化三維心臟軀干幾何,在該幾何模型的基礎(chǔ)上,通過構(gòu)建心臟電生理擴(kuò)散模型和體表外心臟磁場模型,研究心臟電生理活動的演變以及由此產(chǎn)生的體表磁場的分布.通過該框架,可以仿真基于真實(shí)心臟-軀干模型的心臟TMP傳播和體表外磁場分布狀況.最后使用簡化的一維FHN方程和直導(dǎo)線模型驗(yàn)證了心磁正問題計(jì)算框架的可行性.這為后續(xù)心磁逆問題打下堅(jiān)實(shí)基礎(chǔ).