Towards an interactive electromechanical model of the heart

Please download to get full document.

View again

of 12
All materials on our website are shared by users. If you have any questions about copyright issues, please report us to resolve them. We are always happy to assist you.
Document Description
Towards an interactive electromechanical model of the heart
Document Share
Document Tags
Document Transcript
  , 20120091, published 21 February 2013 3 2013 Interface Focus    Hervé DelingetteHugo Talbot, Stéphanie Marchesseau, Christian Duriez, Maxime Sermesant, Stéphane Cotin and   Towards an interactive electromechanical model of the heart   References  This article cites 21 articles, 1 of which can be accessed free Subject collections  (12 articles)medical physics  (17 articles)biomechanics  Articles on similar topics can be found in the following collections Email alerting service   here right-hand corner of the article or click Receive free email alerts when new articles cite this article - sign up in the box at the top go to:  Interface Focus  To subscribe to on March 12, 2014rsfs.royalsocietypublishing.orgDownloaded from on March 12, 2014rsfs.royalsocietypublishing.orgDownloaded from Research Cite this article:  Talbot H, Marchesseau S,Duriez C, Sermesant M, Cotin S, Delingette H.2013 Towards an interactive electromechanicalmodel of the heart. Interface Focus 3:20120091. contribution of 25 to a Theme Issue ‘Thevirtual physiological human: integrativeapproaches to computational biomedicine’. Subject Areas: medical physics, biomechanics Keywords: electrophysiology, cardiac mechanics, trainingsimulation, interactive, SOFA framework  Author for correspondence: Hugo Talbote-mail: Towards an interactive electromechanicalmodel of the heart Hugo Talbot 1,2 , Ste´phanie Marchesseau 2 , Christian Duriez 1 ,Maxime Sermesant 2 , Ste´phane Cotin 1 and Herve´ Delingette 2 1 Shacra Team, Inria Lille - North Europe, Lille, France 2 Asclepios Team, Inria Sophia Antipolis - Me´diterrane´e, Sophia Antipolis, France In this work, we develop an interactive framework for rehearsal of andtraining in cardiac catheter ablation, and for planning cardiac resynchroniza-tion therapy. To this end, an interactive and real-time electrophysiologymodel of the heart is developed to fit patient-specific data. The proposedinteractive framework relies on two main contributions. First, an efficientimplementation of cardiac electrophysiology is proposed, using the latestgraphics processing unit computing techniques. Second, a mechanical simu-lation is then coupled to the electrophysiological signals to produce realisticmotion of the heart. We demonstrate that pathological mechanical andelectrophysiological behaviour can be simulated. 1. Introduction Cardiovascular diseases remain one of the most important causes of deathworldwide. Cardiac arrhythmia and heart failure are life-threatening cardiacdiseases that consist of an abnormal electrical and/or mechanical activity of the heart muscle (myocardium). Such pathologies can occur upon changes inthe heart structure following coronary artery disease (heart attack) or aschronic consequences of hypertension, diabetes or cardiomyopathy. Severaltherapies may be pursued depending on the pathology. For ventricular tachy-cardia, radiofrequency ablation of the ventricles may be performed, whereascases of severe heart failure may be treated through cardiac resynchronizationtherapy (CRT).Most of these procedures are minimally invasive and complex. They need to be performed by experienced and highly skilled cardiologists. Therefore, build-ing training systems for rehearsing these endovascular interventions wouldhelp not only junior electrophysiologists to train for those procedures butalso experienced electrophysiologists to rehearse complex procedures. Thereexists a limited number of commercial endovascular simulators [1,2], but they are somewhat limited because they rely on pre-stored electrophysiology data but not on electromechanical modelling (figure 1).The objective of this paper is to propose an innovative framework for theinteractive simulation of a cardiac electromechanical model. With the constantimprovement in computational methods [3], patient-specific cardiac modellingis now being considered as a promising means of assessing therapy and increas-ing pathophysiological understanding. Physiological models not only canreproduce cardiac motion and electrophysiology signals but also can predict in silico  the impact of therapies [4].Our framework is mainly composed of two physiological models. The first isthe modelling of cardiac electrophysiology. The electrical stimulation inside themyocardium is initiated by a natural pacemaker, called the sinoatrial node,located inside the right atrium. This generates a depolarization wave that propa-gates inside the heart wall, and generated by the exchange of ions through themembrane of cardiac cells. The difference in potential between the cardiac cellsand the extracellular matrix is called the action potential (AP) or transmembranepotential. It drives the contraction of the cardiac fibres. Several mathematicalmodels have been investigated by many research groups and they can besorted into three different classes: & 2013 The Author(s) Published by the Royal Society. All rights reserved.  on March 12, 2014rsfs.royalsocietypublishing.orgDownloaded from   — biophysical models [5] are complex models simulatingelectrophysiology at the organ scale and involve a largenumber of parameters;— phenomenological models [6–9] are simplified models (involving fewer parameters) derived from the biophysi-cal models and capturing the AP shape at the organ scale;— Eikonal models [10] are static nonlinear partial differentialequations for the depolarization times obtained from theprevious models; they are based on non-physiologicalparameters and do not take into account all the complexphysiological phenomena, such as re-entries.The second physiological model of our framework rep-resents the mechanical coupling, i.e. how the depolarization of a cardiac cell generates tension in the sarcomeres, the basicunit of contraction, followed by relaxation. Such cellular mech-anisms are highly complex and involve potentially many ionictypes and parameters. Several models [11–15] have been pro- posed at the tissue scale to describe the relationship betweenthe AP and active tension along the cardiac fibres.In this study, we use the Mitchell–Schaeffer model [9] todescribe the cardiac electrophysiology, because it has onlytwo variables and depends on only a few parameters whilecapturing the main properties of restitution of the actionpotential duration (APD) and conduction velocities. Wealso chose the Bestel–Cle´ment–Sorine (BCS) model to simu-late electromechanical coupling owing to its compact natureand its consistency with the law of thermodynamics (balanceof energy). The advantage in having a compact mechanicalmodel, i.e. a model with few variables and parameters, liesin the ability to calibrate its parameters from a limited setof patient-specific data as shown by Marchesseau  et al . [16]for cardiac mechanics. However, despite their relativesimplicity, it must be stressed that both electrical and mechan-ical models lead to partial differential equations that areknown to be time-consuming to solve. For electrophysio-logical models, the steep potential upstroke constrains thetimeresolution,whereas,formechanicalmodels,isovolumetricconstraints make the numerical system highly nonlinear.Our main contribution lies in adapting and optimizingthese models and the solver algorithms in order to reach afully interactive electromechanical simulation of the heart.In the remainder, we present a real-time simulation of theventricular electrophysiology that is compatible with userinteraction. Then, an interactive mechanical model coupledwith electrophysiology is computed to simulate the associatedheart motion.This study is therefore a first step towards realistic trainingsystems and safer learning procedures for the future. It has been developed using SOFA 1 [17], an open source frameworkfor interactive numerical simulations in medicine. 2. Interactive cardiac electrophysiology 2.1. Geometrical and electrophysiology model Inthissection,wefirstfocusonthereconstructionofthebiven-tricular geometry and on the simulation of electrophysiology based on the Mitchell–Schaeffer [9] model. 2.1.1. Geometry acquisition We apply our interactive framework on both steady-state freeprecession (SSFP) and cine-MRI sequences, which are, respect-ively, dedicatedtoanatomical description andmotiontracking.In order to build a patient-specific geometry from thesesequences, a preliminary image-processing stage is necessary.First, we extract the two ventricles from the SSFP sequence,using a semi-automatic segmentation algorithm available inCardioViz3D. 2 This method requires landmarks to beselected inside and outside the myocardium to fit an implicitelliptical surface. Second, the myocardium mask is meshedusing CGAL software 3 . This results in a tetrahedral meshmade up of around 65 500 linear tetrahedra. Both mechanicaland electrophysiology models are highly related to the fibredirections. Realistic cardiac fibres are generated by syntheti-cally varying the elevation angle across the myocardiumwall. Regarding the left ventricle (including the septum as apart of it), the elevation angle is defined from  þ 70 8  on theendocardium, to 0 8  at the middle of the endocardium andto  2 70 8  in the epicardium. Subsequently, the same methodis applied to the right ventricle (figure 2). 2.1.2. Electrophysiology modelling of the heart: theMitchell–Schaeffer model As introduced previously, there exist several models todescribe the ventricular AP. However, only phenomenological Figure 1.  A view of an MRI image coupled with a three-dimensionalrepresentation of the heart. Figure 2.  A view of the heart mesh including the mapped fibres. r     s   f        s    .r     o     y   a  l        s    o   c   i        e   t        y     p   u   b     l       i        s   h     i       n     g   . o  r       g  I       n   t      e  r    f        a   c    e  F      o   c    u   s    3     :    2     0    1    2     0     0     9    1     2  on March 12, 2014rsfs.royalsocietypublishing.orgDownloaded from   models should be considered for our objective of interac-tive simulation. Indeed, biophysical models characterize theelectrophysiology at the cellular scale, which implies a highlevel of complexity and an important computational cost.Eikonal models are efficient algorithms but have difficultysimulating complex electrophysiology phenomena withnon-physiological parameters.We chose the Mitchell–Schaeffer model because of its sixparameters that can be physiologically identified. Moreover,it captures well the restitution properties of the APD com-pared with phenomenological models of similar complexity.The AP shape described by the Mitchell–Schaeffermodel [9] is given in figure 3. The Mitchell–Schaeffer model, derived from the Fenton–Karma model [18], has two variables:  u , the transmembranepotential, and  z , a secondary variable controlling the repolar-ization phase. The temporal evolution of these two variablesis governed by the following equations: @  t u ¼ div ð D r u Þþ  zu 2 ð 1  u Þ t  in   u t  out þ  J  stim ð t Þ and  @  t  z ¼ ð 1   z Þ t  open if   u , u gate   z t  close if   u . u gate : (9>>>>=>>>>; ð 2 : 1 Þ The diffusion term is defined by a 3  3 anisotropic diffu-sion tensor  D ¼ d  diag ð 1 ;  r ;  r Þ , so that the planar conductionvelocity (CV) in the fibre direction is 2.5 times greater than inthe transverse plane ( r ¼ 1/(2.5) 2 ).  d  is the diffusion coeffi-cient. The parameters  t  in  and  t  out  define the repolarizationphase, whereas the constants  t  open  and  t  close  manage thegate opening or closing depending on the changeover voltage u gate . The default values (describing the common AP) of theseparameters are given by Mitchell & Schaeffer [9], and a simu-lation of an electrophysiological wave is shown in figure 4.The term  J  stim ( t ) is the stimulation current applied in thepacing area. Only a set of nodes corresponding to the apexof the endocardial surface of the left and right ventricleswill include this stimulation current to initiate the depolariz-ation wave (in sinus rhythm) in the ventricles. A secondaryarea can also be defined to model the stimulation induced by a catheter. For our simulations, we initiate the sinusstimulus during 0.1ms so that  J  stim  d t ¼ 0 : 2 ,  J  stim  ¼ 2000s  1 :  ð 2 : 2 Þ Based on the work of Coudiere and colleagues [19], extra-cellular potentials can nevertheless be estimated, becausethe authors have shown that it is possible to recover theextracellular potential from the transmembrane potentialunder a reasonable hypothesis. To do so, an additional ellip-tical equation has to be solved on the surface using analgebraic equation or an integral formulation. We can there-fore easily estimate an extracellular potential using themonodomain Mitchell–Schaeffer model. 2.2. Numerical settings We first focus on the spatial discretization method and theintegration schemes under consideration. We subsequentlydetail a comparison study on the numerical settings of oursimulation. Finally, we explain the way we handle thepersonalization of the Mitchell–Schaeffer model to make itpatient specific. 2.2.1. Spatial discretization As explained in §2.1.1, we compute our three-dimensionalvolumetric mesh from MR images. The biventricular geome-try is meshed with 65500 linear tetrahedra, using CGAL.The Mitchell–Schaeffer model and the diffusion term areimplemented, using the finite element method (FEM). Theedge length for cardiac electrophysiology is usually less than0.5mm. However, we use larger tetrahedra (with an edgelength of approx. 2mm) in our simulations. Using largerelements means fewer tetrahedra within the mesh, thus ensur-ing higher computational performances. The influence of theelement size will be studied in more detail in §2.2.3.Recent work considered other spatial discretizationmethods. For instance, Rapaka  et al.  [20] presented theiralgorithm applying the lattice Boltzmann method to cardiacelectrophysiology (LBM-EP). In this work, the authorsexplained that traditional FEMs could not offer fast electro-physiology simulations. In our paper, we propose using theFEM model on a graphics processing unit (GPU) as a goodalternative to compensate for the FEM computational cost. 2.2.2. Integration schemes We now detail the integration schemes we considered.Ethier & Bourgault [21] presented the main schemes used 11   a  c   t   i  o  n   p  o   t  e  n   t   i  a   l 2340APD t   (m   s –1 )500 Figure 3.  Transmembrane potential as described by Mitchell & Schaeffer [9]. ( a )( b )( c )( d  ) Figure 4.  Four steps of a depolarization wave propagating inside a FEMmesh: depolarization starts in ( a ) and the wave then propagates in ( b , c  )until both ventricles are depolarized ( d  ). r     s   f        s    .r     o     y   a  l        s    o   c   i        e   t        y     p   u   b     l       i        s   h     i       n     g   . o  r       g  I       n   t      e  r    f        a   c    e  F      o   c    u   s    3     :    2     0    1    2     0     0     9    1     3  on March 12, 2014rsfs.royalsocietypublishing.orgDownloaded from   for cardiac electrophysiology. We implemented the modifiedCrank–Nicholson/Adams–Bashforth (MCNAB), a second-order,semi-implicitscheme.Thediffusiontermisimplicitlyinte-grated, whereas the ionic current term is explicitly computed.Electrophysiology equations can be generalized as follows: r  @ u @ t  ¼  f  ð u Þþ  g ð u ;  z Þ and  @  z @ t  ¼ h ð u ;  z Þ ; 9>>=>>; ð 2 : 3 Þ where  g  and  h  denote the ionic functions of the Mitchell–Schaeffer model, and  f   is the additional diffusion term (seeequation(2.1)).TheMCNAB schemedescribestheintegration as  M  _ u n þ 1  ¼  916  F ð u n þ 1 Þþ 38  F ð u n Þþ  116  F ð u n  1 Þþ 32 G ð u n ;  z n Þ 12 G ð u n  1 ;  z n  1 Þ and  M  _  z n þ 1  ¼ 32  H  ð u n ;  z n Þ 12  H  ð u n  1 ;  z n  1 Þ ; 9>>>>>>>=>>>>>>>; ð 2 : 4 Þ where  M  represents the operator obtained by integrating theterm of mass density  r  ;  G  and  H   denote the operator obtained by integration of the ionic term; and  F  denotes the operatorobtained by integration of the diffusion term. In this lastsystem, the index number  n  refers to the  n th time step. Usingthe MCNAB, our simulation includes a conjugate gradientcoupled with a Jacobi preconditioner in order to efficientlysolve the linear system ( A x ¼ b ). Because the matrix  A  isdiagonal dominant, the choice of this preconditioneris straightforward. The preconditioned matrix given by the Jacobi can be updated during the simulation in order to takeinto account conductivity changes (owing to thermo-ablation).We also implemented a fully explicit scheme with a back-ward differentiation. In the following, we refer to it as‘explicit-BDF’. It can be solved without a linear solver. Usingthe notation from equation (2.3), it can be written as: u n þ 1  ¼ 4 u n  u n  1 3  þ 23  F ð u n Þþ G ð u n ;  z n Þð Þ  M  1 and  z n þ 1  ¼ 4  z n   z n  1 3  þ 23  H  ð u n ;  z n Þð Þ  M  1 : 9>>=>>; ð 2 : 5 Þ Maps of depolarization times obtained using bothMCNAB and explicit-BDF solvers are shown in figure 5 a and figure 5 b , respectively. 2.2.3. Numerical study To make the best compromise between performance andaccuracy, we study each of the numerical parameters of oursimulation. The influence of the element size regarding theintegration method is considered as well as the influence of the time step. Elementsize. Wefirstfocusonacrucialfeatureregardingtheper-formance context: the element size. Pathmanathan  et al.  [22]studied the relationship between the CV and the element sizedepending on the integration method. We reproduced theone-dimensional planar propagation wave, based on the cellmodel of Mitchell–Schaeffer, on a simplified geometry: athree-dimensional cuboid such as the one described byPathmanathan  et al.  [22]. We implemented the lumped integra-tion and the ionic current integration (ICI) that were described by Pathmanathan  et al.  [22]. The lumping method consists of summing all the coefficients of a line onto the diagonal. TheICIinterpolatesthenodalioniccurrentlinearlyoneachelement.The results using the lumped integration and the ICIare presented in figure 6. We obtain the same trend asPathmanathan  et al.  in their work. The ICI method is amore accurate technique than the lumped integration. How-ever, the ICI integration method is more computationallydemanding than the lumping method. We therefore decidedto rely on a lumped integration of the ionic term of theMitchell–Schaeffer model.As presented in figure 6, the coarser the mesh, the slower theCV (with a constant diffusion coefficient). In our approach, wepropose adapting the nodal diffusion coefficient  d  in order to fitthe patient-specific map of depolarization times. The integrationerror owing to the lumping approximation would therefore benumericallycompensatedwhilebenefitingfromlargerelements,i.e. better performances. This will be carried out during theelectrophysiology personalization detailed in §2.2.4. Time step.  Depending on the integration scheme used, the timestep is limited either owing to the diffusion term or owing tothe ionic term, depending on whether the explicit-BDF or theMCNAB scheme (implicit diffusion) is used. These stability con-ditionshavealreadybeentackledbyEthier&Bourgault [21]and by Coudie`re & Pierre [23]. It can be shown that:— using the semi-implicit MCNAB scheme, the ionic currentdefines the limit of stability: d t  1 = j inf  ð @ R Þj¼ t  in t  out = t  in þ t  out  ¼ 0 : 286ms ;  where  R  denotes the ionic term;— using the explicit-BDF scheme, the diffusion termdetermines the limit of stability. The time step has to bed t  min(d x /2  CV), where d t  is the characteristiclength of the elements, and CV is the CV of the depolar-ization wave inside this tetrahedron. However, the limitowing to the ionic term is still active. This means that d t cannot exceed 0.286 ms whatever the scheme. Usingthe mesh described in §2.1.1, the stability limit using theexplicit-BDF scheme is d t , 0.13 ms. depo la r i z a t i on t i me0. 14 8 4 290. 1 20.080.0 4 0depo la r i z a t i on t i me0. 14 8 4 290. 1 20.080.0 4 0 ( a ) ( b ) Figure 5.  Isochrones computed using ( a ) the MCNAB solver and ( b ) the explicit-BDF solver. Red areas represent scars segmented from MR images. r     s   f        s    .r     o     y   a  l        s    o   c   i        e   t        y     p   u   b     l       i        s   h     i       n     g   . o  r       g  I       n   t      e  r    f        a   c    e  F      o   c    u   s    3     :    2     0    1    2     0     0     9    1     4  on March 12, 2014rsfs.royalsocietypublishing.orgDownloaded from 
Similar documents
View more...
Search Related
We Need Your Support
Thank you for visiting our website and your interest in our free products and services. We are nonprofit website to share and download documents. To the running of this website, we need your help to support us.

Thanks to everyone for your continued support.

No, Thanks