鄭詩(shī)強(qiáng) 張鎮(zhèn)西 李怡勇 郭 赤

基于OpenCMISS的心臟電生理模擬解決方案

鄭詩(shī)強(qiáng)① 張鎮(zhèn)西② 李怡勇① 郭 赤①

目的：心臟電生理模擬對(duì)于深入了解心臟的內(nèi)在工作機(jī)制和發(fā)病機(jī)制具有重要意義，為解決由于心臟電生理模擬涉及的微分方程和偏微分方程中實(shí)際編程問題，探討研究心臟電生理模擬的高效編程方案與設(shè)計(jì)。方法：借助新西蘭奧克蘭生物工程研究所開發(fā)的數(shù)學(xué)建模環(huán)境OpenCMISS，給出了OpenCMISS的一些基本概念、編程原則、流程以及模擬實(shí)例。結(jié)果：OpenCMISS的良好編程環(huán)境能夠有效實(shí)現(xiàn)心臟的電生理模擬。結(jié)論：基于OpenCMISS的心臟電生理模擬解決方案可極大減少程序開發(fā)的難度，為從事電生理模擬領(lǐng)域的研究者提供有效的工具。

電生理模擬；心臟；OpenCMISS；工具

[First-author’s address]Department of Medical Engineering, Wuhan General Hospital of Guangzhou Command, Wuhan 430070, China.

心臟病是威脅人類健康的三大殺手之一。據(jù)2008年《世界衛(wèi)生報(bào)告》稱，全球每年因心血管疾病死亡的人數(shù)約為1730萬(wàn)，占死亡總?cè)藬?shù)的30%，是人類死亡的首要原因。深入理解心臟的內(nèi)在工作機(jī)制和發(fā)病機(jī)制，為預(yù)防和治療心臟疾病提供有效手段是一項(xiàng)重大而艱巨的課題。心臟的仿真模擬相對(duì)于其他研究手段，花費(fèi)的時(shí)間和資金成本低，不會(huì)受到倫理問題的困擾，能夠方便地在定量條件下研究和預(yù)測(cè)心臟的工作機(jī)制和發(fā)病機(jī)制，其模型的基礎(chǔ)來(lái)源于實(shí)驗(yàn)，同時(shí)其仿真結(jié)果又為實(shí)驗(yàn)提供理論指導(dǎo)。

1 心臟電生理模擬概述

心臟的仿真模擬是一項(xiàng)綜合而富有挑戰(zhàn)性的工作，其涉及到心臟的各種生理現(xiàn)象，如血液流動(dòng)、機(jī)械收縮、神經(jīng)調(diào)控和電信號(hào)傳導(dǎo)等，需要整合科學(xué)界各種最新的研究成果，包括細(xì)胞、組織、器官的生物化學(xué)、生物物理和形態(tài)學(xué)等數(shù)據(jù)，以從不同時(shí)間、空間尺度下滿足在分子、蛋白、細(xì)胞、組織、器官和生物體等不同層次的模擬研究[1]。心臟仿真模擬偶聯(lián)的生理過程十分復(fù)雜，研究心臟的單一生理過程可極大簡(jiǎn)化模擬的難度，是全心臟仿真模擬的起點(diǎn)。而心臟的電活動(dòng)是心臟興奮收縮的主要?jiǎng)恿λ冢瑢?duì)心臟的電生理模擬是心臟研究的一大熱點(diǎn)。

心臟電生理模擬分為單細(xì)胞層次和多細(xì)胞層次，單細(xì)胞模型通常基于Hodgkin-Huxley方程(HH方程)[2]。隨著實(shí)驗(yàn)技術(shù)的提高(如單細(xì)胞、單通道膜片鉗技術(shù)的發(fā)展)，越來(lái)越多的離子通道及細(xì)胞過程被發(fā)現(xiàn)，同時(shí)更加精確的數(shù)學(xué)語(yǔ)言描述的引入，模型變得越來(lái)越復(fù)雜[3]。通常單細(xì)胞模型是由數(shù)十個(gè)代數(shù)方程和常微分方程組混合而成。對(duì)于微分方程組的算法，目前已經(jīng)非常成熟，常用的有歐拉法(Euler method)、龍格-庫(kù)塔法(Runge-Kutta method)、以及線性多步法等[4]。在利用這些算法提高收斂速度、計(jì)算時(shí)間時(shí)還可以加入時(shí)間步長(zhǎng)控制來(lái)提高計(jì)算速度。通常電壓變化越大，時(shí)間步長(zhǎng)取得越小。

1.1 電興奮傳導(dǎo)模型—雙域模型與單域模型

由于對(duì)計(jì)算要求過高，多細(xì)胞模型發(fā)展歷史較短。最早的多細(xì)胞模型由Winslow于1993年提出[5]。電纜理論是多細(xì)胞模型的數(shù)學(xué)基礎(chǔ)。Tung[6]基于此提出雙域方程。雙域方程可以進(jìn)一步簡(jiǎn)化為單域方程，計(jì)算速度更快，適用于多數(shù)仿真研究，但有時(shí)必須使用雙域[7]。雙域方程式(公式1)：

其中，Am為細(xì)胞表體比，Vm為跨膜電位，Cm為膜電容，Iion為離子通道電流(外向?yàn)檎?，Di為膜內(nèi)電導(dǎo)率張量，Φe為膜外電勢(shì)，t為時(shí)間。

如果膜外電阻與膜內(nèi)電阻相比可以忽略不計(jì)(De無(wú)窮大)，或者膜外電阻同膜內(nèi)電阻分布場(chǎng)相似(即Di=k·De，k為常數(shù))，雙域模型可以進(jìn)一步簡(jiǎn)化為單域模型(公式2)：

該模型可以視作一個(gè)典型的反應(yīng)擴(kuò)散問題，包括反應(yīng)和擴(kuò)散兩部分，即動(dòng)作電位的產(chǎn)生和傳導(dǎo)。其基本算法為有限差分法(finite difference method，F(xiàn)DM)、有限元法(finite element method，F(xiàn)EM)、有限體積法(finite volume method，F(xiàn)VM)、邊界元法(boundary element method，BEM)等[8-10]。由于所涉方程數(shù)目眾多，常微分方程、偏微分方程和非線性代數(shù)方程悉列其中，加之細(xì)胞類型、節(jié)點(diǎn)數(shù)目和傳導(dǎo)參數(shù)復(fù)雜多樣，造成實(shí)際編程求解的困難。為了更好地解決心臟建模中的算法問題，很多研究中心開發(fā)了免費(fèi)軟件或函數(shù)庫(kù)供研究者使用，如：libMesh、deal.II、FETK、OOFEM、Continuity、CARP、Chaste、OpenCMISS等[11-18]。其中后4個(gè)專門用于心臟建模領(lǐng)域。

1.2 OpenCMISS函數(shù)庫(kù)

OpenCMISS是一個(gè)開源的函數(shù)庫(kù)，其前身是連續(xù)介質(zhì)力學(xué)、圖像分析、信號(hào)處理及系統(tǒng)識(shí)別(Continuum Mechanics，Image analysis，Signal processing and System Identification，CMISS)，由新西蘭奧克蘭生物工程研究所自1980年開始開發(fā)，致力于求解涉及多個(gè)偏微分方程耦聯(lián)的生物工程問題(意味著有多個(gè)生理過程耦聯(lián)，如機(jī)械電耦聯(lián))，目前已經(jīng)發(fā)展得非常成熟，可以解決軟組織力學(xué)、流體力學(xué)、熱量物質(zhì)傳輸?shù)葐栴}及其相互耦聯(lián)。該軟件庫(kù)包含計(jì)算引擎和結(jié)果可視化軟件。其算法不僅包括有限元分析，還包括邊界元法和有限差分法。

根據(jù)新一代CellML、FieldML標(biāo)準(zhǔn)的公布以及開源的需要，其新版本為OpenCMISS[19-20]。該函數(shù)庫(kù)基于Fortran95語(yǔ)言開發(fā)，可在Ubuntu或Fedora版本的Linux系統(tǒng)平臺(tái)下運(yùn)行，實(shí)現(xiàn)了并行和串行計(jì)算。

2 OpenCMISS中的基本對(duì)象

(1)基本函數(shù)：即有限元插值函數(shù)構(gòu)造法，包括Lagrange-Hermite法和Simplex。

(2)坐標(biāo)系統(tǒng)：包括笛卡爾坐標(biāo)系、柱坐標(biāo)系、球坐標(biāo)系、長(zhǎng)球面坐標(biāo)系和扁球面坐標(biāo)系。

(3)區(qū)域：即研究所在位置區(qū)域，如心臟、肺、心肌及心室，提供一個(gè)命名空間和容器以存放各種對(duì)象。具有分層結(jié)構(gòu)，可以含有多個(gè)子區(qū)域。

(4)網(wǎng)格：計(jì)算網(wǎng)格的拓?fù)浣Y(jié)構(gòu)，包含節(jié)點(diǎn)、面元和基本函數(shù)。

(5)域：將網(wǎng)格分解為多個(gè)區(qū)域以進(jìn)行并行計(jì)算。

(6)場(chǎng)：即場(chǎng)變量的空間分布。場(chǎng)的類型有：方程描述場(chǎng)、幾何場(chǎng)、纖維場(chǎng)、材料場(chǎng)、獨(dú)立變量場(chǎng)、分析場(chǎng)、變量場(chǎng)和源場(chǎng)。

(7)方程組：解決特定問題的方程組，與區(qū)域相聯(lián)系。

(8)問題和控制回路：?jiǎn)栴}對(duì)應(yīng)著一個(gè)特定的計(jì)算流程。

(9)算法器和求解方程。

3 OpenCMISS的安裝及一般技巧

OpenCMISS運(yùn)行于Linux平臺(tái)，其大致分為4步：①安裝額外必要的常用軟件或函數(shù)庫(kù)，如subversion、cmake、swig等，使用sudo apt-get install命令進(jìn)行安裝；②下載程序及附帶函數(shù)庫(kù)，使用svn co https：//svn.physiomeproject.org/svn/ opencmissextras/utils/trunk utils命令下載額外庫(kù)，從Physiome Project官網(wǎng)的資源鏈接[21]里下載安裝包；③使用sudo gedit～/.bashrc命令，加入環(huán)境變量OPENCMISS_ROOT、OPENCMISSEXTRAS_ ROOT；④使用make命令，安裝OpenCMISS。其具體安裝步驟可參見相關(guān)網(wǎng)頁(yè)[22]。

OpenCMISS的基本對(duì)象的創(chuàng)建必須經(jīng)過初始化、創(chuàng)建、設(shè)置及創(chuàng)建完成的過程。其相應(yīng)函數(shù)的后綴分別為_Initialise、_CreateStart、Set、_CreateFinish。

對(duì)于單域模型電生理問題，只需在OpenCMISS的examples/Bioelectrics/Monodomain/src文件夾里創(chuàng)建Fortran源文件，編寫相應(yīng)代碼即可，其步驟如下：

(1)對(duì)研究區(qū)域進(jìn)行網(wǎng)格化時(shí)需對(duì)node、mesh對(duì)象進(jìn)行創(chuàng)建，同時(shí)在不同的維度方向上創(chuàng)建mesh elements對(duì)象，使用CMISSMeshElements_NodesSet函數(shù)將面元和節(jié)點(diǎn)相互對(duì)應(yīng)，如圖1所示。

圖1 面元節(jié)點(diǎn)對(duì)應(yīng)關(guān)系示意圖

(2)創(chuàng)建GeometricField，使用CMISSField_ ParameterSetUpdateNode函數(shù)將每個(gè)節(jié)點(diǎn)的坐標(biāo)位置標(biāo)定清楚。

(3)創(chuàng)建MaterialsField，設(shè)定心電傳導(dǎo)參數(shù)分布場(chǎng)，即細(xì)胞表體比Am、膜電容Cm以及細(xì)胞間導(dǎo)電率的分布。

(4)初始化CellML環(huán)境，使用CMISS CellML_ ModelImport函數(shù)將CellML語(yǔ)言編輯的單細(xì)胞模型導(dǎo)入。

此類單細(xì)胞模型可以使用可視化編程工具OpenCell進(jìn)行創(chuàng)建，有許多現(xiàn)成的模型可提供下載[23-24]。該軟件基于XML語(yǔ)言，將細(xì)胞模型視為網(wǎng)狀結(jié)構(gòu)，包括組分(含各個(gè)離子通道電流中的所有變量以及數(shù)學(xué)表達(dá)式)、連接(表示各個(gè)組分中相同變量間的依賴關(guān)系，包括private，public兩大類型)、單位(所有變量的單位)、輸入(可選部分，將其他細(xì)胞模型的成分導(dǎo)入，實(shí)現(xiàn)模型組件的重復(fù)利用)以及群組(表達(dá)各個(gè)組分間包含關(guān)系)5大部分。其軟件界面如圖2所示。

圖2 OpenCell界面

(5)根據(jù)需要將單細(xì)胞模型中的變量設(shè)為“Known”或“Wanted”。

(6)將細(xì)胞膜電壓Vm和DependentField場(chǎng)中電壓V相互映射，并設(shè)定其初值。

(7)創(chuàng)建CellMLModelsField，即細(xì)胞類型分布場(chǎng)，將每個(gè)節(jié)點(diǎn)對(duì)應(yīng)已導(dǎo)入的單細(xì)胞模型。

(8)創(chuàng)建CellMLParametersField，對(duì)標(biāo)為“Known”的變量建立場(chǎng)分布。

(9)創(chuàng)建Problem對(duì)象，設(shè)定其問題類型為“MONODOMAIN”；創(chuàng)建ControlLoop，設(shè)定模擬時(shí)間；創(chuàng)建Solver對(duì)象，設(shè)定ODE、PDE問題的算法；創(chuàng)建和設(shè)定邊界條件，一般為Neumann邊界條件。

(10)使用CMISSProblem_Solve函數(shù)求解；使用CMISSFinalise函數(shù)結(jié)束程序。

(11)結(jié)果處理：將生成的時(shí)間序列文件用CMGUI軟件導(dǎo)入，進(jìn)行相應(yīng)圖像處理，將結(jié)果展示出來(lái)[25]。

當(dāng)然也可以使用Perl語(yǔ)言將結(jié)果提取出來(lái)，再到專門圖像處理軟件中進(jìn)行處理。

4 編程實(shí)例

竇房結(jié)作為心臟的正常起搏位點(diǎn)，是一個(gè)結(jié)構(gòu)復(fù)雜的組織。研究發(fā)現(xiàn)，竇房結(jié)組織由竇房結(jié)細(xì)胞、心房細(xì)胞和纖維原細(xì)胞等多種細(xì)胞組成。Verheijck等[26]提出混合模型，不同區(qū)域細(xì)胞組成比例不同，而這一差異正是造成竇房結(jié)細(xì)胞動(dòng)作電位傳導(dǎo)的重要因素。本研究采取20×20細(xì)胞竇房結(jié)方形組織，含有心房細(xì)胞和竇房結(jié)細(xì)胞，但上下部分心房細(xì)胞含量各不相同，分別為60%，40%。采用偽隨機(jī)方法生成偽竇房結(jié)方形組織，如圖3所示。

圖3 竇房結(jié)混合模型

竇房結(jié)細(xì)胞采用Maltsev模型[27]，細(xì)胞尺寸假定為70 μm×4 μm(長(zhǎng)×半徑)的長(zhǎng)柱體，實(shí)驗(yàn)獲得的細(xì)胞間偶聯(lián)電導(dǎo)率為gc=7.5 nS，電導(dǎo)率張量同細(xì)胞間偶聯(lián)電導(dǎo)率的關(guān)系[28]如公式3：

表體比為公式4：

可求得Am=530 mm-1，D=10.4 nS mm-1。

心房細(xì)胞采用Lindblad等[29]提出的細(xì)胞模型。利用OpenCMISS編程，最終產(chǎn)生的竇房結(jié)組織興奮傳播偽彩圖，可見不同區(qū)域細(xì)胞組成比例可以引起興奮向心房細(xì)胞含量多的區(qū)域傳播(如圖4所示)。

圖4 竇房結(jié)混合模型興奮傳播偽彩圖

5 結(jié)論

深入了解心臟的電生理特性等內(nèi)在工作機(jī)理，對(duì)心臟疾病的預(yù)防和治療有極其重要的意義。心臟的模擬仿真可以避免倫理道德問題的尷尬，其投入時(shí)間和金錢均較其他研究方法低，且可在理想?yún)?shù)下研究和預(yù)測(cè)心臟的發(fā)病機(jī)制，為其他研究方法提供理論指導(dǎo)，是必不可少的研究手段。

OpenCMISS提供了一個(gè)極好的編程環(huán)境，其結(jié)構(gòu)的模塊化可以極大減少心臟電生理模擬開發(fā)、調(diào)試的時(shí)間，使研究者更關(guān)注于模擬結(jié)果的生理闡釋。通過有限元方法設(shè)計(jì)更為實(shí)際的心臟二維或三維結(jié)構(gòu)，收集更加完善的生理參數(shù)，將使心臟電生理模型更加符合實(shí)際，進(jìn)而更好的預(yù)測(cè)和研究心臟電生理活動(dòng)。由于OpenCMISS提供了多個(gè)偏微分方程耦聯(lián)的解決方案，未來(lái)將心臟的其他生理活動(dòng)如血液流動(dòng)、機(jī)械收縮、神經(jīng)調(diào)控和電信號(hào)傳導(dǎo)等耦聯(lián)進(jìn)去，通過并行超大計(jì)算，可實(shí)現(xiàn)一個(gè)鮮活心臟的仿真模擬[30]。將對(duì)心臟疾病的研究、藥物等治療手段的研發(fā)提供有效的理論指導(dǎo)，為人類健康事業(yè)的進(jìn)步提供推動(dòng)力。

[1]Hunter P，Robbins P，Noble D.The IUPS human physiome project[J].Pflugers Arch，2002，445(1)：1-9.

[2]Hodgkin AL，Huxley AF.A quantitative description of membrane current and its application to conduction and excitation in nerve.1952[J].Bull Math Biol，1990，52(1-2)：25-71.

[3]Fink M.Noble ND.Markov models for ion channels：versatility versus identifiability and speed[J].Philos Trans A Math Phys Eng Sci，2009，367(1896)：2161-2179.

[4]李榮華.微分方程數(shù)值解法[M].4版.北京：高等教育出版社，2009.

[5]Winslow RL，Varghese A，Noble D，et al.Generation and propagation of ectopic beats induced by spatially localized Na-K pump inhibition in atrial network models[J].Proc Biol Sci，1993，254(1339)：55-61.

[6]Tung，Leslie.A bidomain model for describing ischemic myocardial DC potentials[D].Cambridge：Massachusetts Institute of Technology，1978.

[7]Eason J，Trayanova N.Phase singularities and termination of spiral wave reentry[J].J Cardiovasc Electrophysiol，2002，13(7)：672-679.

[8]Qu Z Garfinkel A.An advanced algorithm for solving partial differential equation in cardiac conduction[J].IEEE Trans Biomed Eng，1999，46(9)：1166-1168.

[9]Nicotson J，Ca E.A practical method for numerical evaluation of solutions of partial differential equations of the heat-conduction type[J].Proc Camb Phil Soc，1947，43(1)：50-67.

[10]Rogers JM，McCulloch AD.A collocation-Galerkin finite element model of cardiac action potential propagation[J].IEEE Trans Biomed Eng，1994，41(8)：743-757.

[11]Kirk JW，Stogner RH.libMesh：a C++ library for parallel adaptive mesh refinement/coarsening simulations[J].Eng Comput，2006，22(2)：237-254.

[12]Bangerth RH，Kanschat.II—a general-purpose object-oriented finite element library[J].ACM Trans Math So ftware，2007，33(1)：24-27.

[13]Olst M.Adaptive numerical treatment of elliptic systems on manifolds[J].AdvComput Math，2001，15(2)：139-191.

[14]Patzák ZB.Design of object oriented finite element code[J].Adv Eng Software，2001，32(3)：759-767.

[15]Jeff Van Dorn.Continuity 6：A problem solving environment for multi-scale Biology[EB/OL]. (2013-05-01)[2013-09-03]www.continuity. ucsd.edu.

[16]Vigmond EJ，Plank MH，Leon LJ.Computational tools for modeling electrical activity in cardiac tissue[J].J Electrocardiol，2003，36(1)：69-74.

[17]Pitt-Francis JP，Bernabeu MO.Chaste：a testdriven approach to software development for biological modelling[J].Comput Phys Commun，2009，180(16)：2452-2471.

[18]Bradley BA，Britten R.OpenCMISS：a multiphysics and multi-scale computational infrastructure for the VPH/Physiome project[J]. Prog Biophys Mol Biol，2011，107(1)：32-47.

[19]Cuellar AA，Lloyd CM， Nielsen PF，et al.An overview of CellML 1.1，a biological model description language[J].SIMULATION：Trans Soc for Model Simulation Int，2003，79(3)：740-747.

[20]Christie GR，PMN，SA.Blackett，et al.Field-ML：concepts and implementation[J].Philos Transact A Math PhysEngSci，2009，367(10)：1869-1884.

[21]Chris Bradley.OpenCMISS[EB/OL].(2011-05-22) [2013-09-03]https：//github.com/OpenCMISS.

[22]Chris Bradley.OpenCMISS documentation[EB/OL]. (2013-09-03)[2013-09-03]http：//cmiss.bioeng. auckland.ac.nz/OpenCMISS/doc/programmer.

[23]Michael T，Cooling，OpenCell[EB/OL].(2010-02-28)[2013-09-03]http：//www.cellml.org/tools/ opencell.

[24]Cell ML.Model Repository[EB/OL].(2010-05-02) [2013-09-03]http：//models.cellml.org/cellml.

[25]Peter Hunter.CMGUI[EB/OL].(2012-04-16)[2013-09-03]http：//www.cmiss.org/ReleaseCenter/ cmgui.

[26]Verheijck EE，Wessels A，van Ginneken AC，et al.Distribution of atrial and nodal cells within rabbit sinoatrial node：models of sinoatrialtransition[J].Circu lation，1998，97(16)：1623-1631.

[27]Maltsev VA，Lakatta EG.Synergism of coupled subsarcolemmal Ca2+clocks and sarcolemmalvoltage clocks confers robust and flexible pacemaker function in a novel pacemaker cell model[J].Am J Physiol Heart Circ Physiol，2009，296(3)：594-615.

[28] Garny A，Kohl P，Hunter PJ，et al.One-dimensional rabbit sinoatrial node models：benefits and limitations[J].J Cardiovasc Electrophysiol，2003，14(10 Suppl)：S121-132.

[29]Lindblad DS，Murphey CR，Clark JW，et al.A model of the action potential and underlying membrane currents in a rabbit atrial cell[J].Am J Physiol，1996，271(4)：1666-1696.

[30]Santos D，Weber R，Rafael SO，et al.Parallel and adaptive numerical techniques applied to cardiac modeling[J].Main Thematic Areas，2012，21：151.

Research on one novel scheme for cardiac electrophysiological simulation based on OpenCMISS

/ZHENG Shi-qiang, ZHANG Zhen-xi, LI Yi-yong, et al// China Medical Equipment,2014,11(1):42-46.

Objective:The cardiac electrophysiological simulation can further the insight of intrinsic mechanism of heart and its related diseases. However, as too many ordinary differential equations, partial differential equations and cell nodes involved in this simulation, it is very hard for ordinary people to solve the problem programmatically (in code).To solve this problem, we intended to investigate effective programming design for cardiac electrophysiological simulation.Methods:With the aid of OpenCMISS developed by Auckland Bioengineering Institute in New Zealand, some key concepts of OpenCMISS, basic programming principle or strategies and one example are illustrated.Results:The friendly programming environment of OpenCMISS can promote the programming of cardiac elctrophysiological simulation.Conclusion:This novel scheme for cardiac electrophysiological simulation based on OpenCMISS can greatly reduce the difficulty of the application development and provide a powerful tool for researchers in electrophysiological simulation.

Electrophysiological simulation; Heart; OpenCMISS; Tools

1672-8270(2014)01-0042-05

R319

A

10.3969/J.ISSN.1672-8270.2014.01.015

2013-08-26

①?gòu)V州軍區(qū)武漢總醫(yī)院醫(yī)學(xué)工程科 湖北 武漢 430070

②西安交通大學(xué)生命學(xué)院生物醫(yī)學(xué)分析及儀器研究所 陜西 西安 710049

鄭詩(shī)強(qiáng),男,(1989- ),碩士，主管技師。廣州軍區(qū)武漢總醫(yī)院醫(yī)學(xué)工程科，從事醫(yī)院信息化、醫(yī)療設(shè)備開發(fā)與維護(hù)等科研工作。