基于個(gè)性化三維心臟-軀干模型的心磁正問(wèn)題*

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

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

1 引 言

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

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

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

2 方 法

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

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

本研究使用一位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)成.使用最大類(lèi)間方差法[17]對(duì)心臟和軀干圖像切片依次進(jìn)行圖像分割,結(jié)合手繪處理的方式,分離出心臟中心室和心房的內(nèi)外輪廓和軀干表面輪廓,如圖1所示.

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

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

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

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

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

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

FHN模型存在較多變式,為了更為準(zhǔn)確地表示真實(shí)心臟TMP變化(如去除原始FHN模型中實(shí)際不存在的過(guò)極化現(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)為心臟作為一個(gè)孤立的連續(xù)體,即沒(méi)有有功電流流入或流出心臟表面,所以存在紐曼邊界條件?u/?n=0,其中n表示心臟邊界的法向量.

2.3 體表外心臟磁場(chǎng)模型

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

簡(jiǎn)化(15)式,可以求得

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

3 仿真結(jié)果

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

通過(guò)Tanh法,求得了一維條件下簡(jiǎn)化的FHN方程的15個(gè)行波解,本文以其中一個(gè)解為例

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

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

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

圖5 t=20,40,60時(shí),激發(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)入一個(gè)u=1 的激發(fā)電位,作為FHN方程的初值條件模擬TMP興奮的激發(fā).之后通過(guò)有限元法計(jì)算,可得到心臟TMP隨時(shí)間t變化的動(dòng)態(tài)分布.

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

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

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

圖6 一個(gè)心動(dòng)周期的TMP分布圖 (a),(b),(c) 為去極化過(guò)程 (t=30,60,90); (d),(e),(f)為復(fù)極化過(guò)程(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í)間變化曲線 實(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φ解的相對(duì)誤差RE分布.其中散點(diǎn)表示中軸面上離導(dǎo)線中心各距離所對(duì)應(yīng)的誤差分布.而RE曲線由基于最小二乘法的二次多項(xiàng)式擬合得到.從圖9中可以看出,隨著離直導(dǎo)線中心距離的增加,RE也隨之變大.0.1≤r≤1.0m時(shí),RE 最大為2.85%.對(duì)比心臟到體表處的距離,我們認(rèn)為這個(gè)誤差是可以接受的.經(jīng)過(guò)直導(dǎo)線模型的計(jì)算驗(yàn)證,說(shuō)明構(gòu)建的伽遼金法求解體表外心臟磁場(chǎng)模型是可行的.

圖9 磁場(chǎng)的相對(duì)誤差分布Fig.9.Relative error distribution of magnetic field.

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

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

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

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

圖10 體表外磁感應(yīng)強(qiáng)度 Bz 分布 (a)-(h) 模擬的心磁分布圖(t=40,60,80,100,120,140,160,180)；(i) 文獻(xiàn)中的實(shí)測(cè)心磁圖分布圖[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ù)來(lái)自被試的MRI切片數(shù)據(jù).由于MRI切片數(shù)據(jù)的高精度性以及需求的可定制性,為更好地用于分割重建真實(shí)并且各種分辨率的三維心臟-軀干幾何模型提供了基礎(chǔ).并且為后續(xù)定制化仿真被試的心磁分布與實(shí)測(cè)心磁圖的對(duì)比提供了良好的條件.

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

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

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

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

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