F.B. Sachse

G.T. Lines

J. Grice Le

J. Rogers

L. Gerardo-Giorda

L. Mirabella

L.F. Pavarino

M. Pennacchio

M. Potse

M.J. Gander

N. Trayanova

P. Colli Franzone

R.H. Clayton

English

Luca Gerardo-Giorda

Lucia Mirabella

Mauro Perego

Alessandro Veneziani

Crossref

Optimized Schwarz Methods and Model Adaptivity in Electrocardiology Simulations

Gerardo-Giorda, L.

Mirabella, L.

Perego, M.

Veneziani, A.

BCAM's Institutional Repository Data

Optimized Schwarz Methods and modeladaptivity in electrocardiology simulationsLuca Gerardo-Giorda, Lucia Mirabella, Mauro Perego and Alessandro VenezianiAbstract The Bidomain model is nowadays one of the most accurate mathematicaldescriptions of the action potential propagation in the heart. However, its numericalapproximation is in general fairly expensive as a consequence of the mathemati-cal features of this system, and several works have been devoted to devise effec-tive solvers and preconditioners, [10, 11, 15] among others. A simplification of thismodel, called Monodomain problem is often adopted in order to reduce computa-tional costs of the numerical solution of the cardiac potential. The latter model ishowever less accurate. A possible trade-off between accuracy and cost is a modeladaptive strategy. The computational domain is subdivided into regions, coupledthrough an Optimized Schwarz Method, in which either the Bidomain or the Mon-odomain problem is solved, according to an a-posteriori model error estimator fol-lowing the spatio-temporal evolution of the action potential propagation. Here wepresent a possible implementation of this approach, following up previous works onthe error estimation and the Optimized Schwarz coupling.1 The Numerical ModelAt the macroscopic level, the myocardial tissue is regarded as the superposition oftwo continuous and anisotropic media, the intra-cellular and the extra-cellular one.Luca Gerardo-GiordaBCAM - Basque Center for Applied Mathematics, Bilbao, Spain, e-mail: lgerardo@bcamath.orgLucia Mirabella“W.H. Coulter” Dept. Biomed. Engrg., GaTech, Atlanta, GA, USA e-mail: lu-cia.mirabella@bme.gatech.eduMauro PeregoDept. of Sci. Comp., Florida State University, Thallahassee, FL, USA e-mail: mperego@fsu.eduAlessandro VenezianiDept. of Math. and CS, Emory University, Atlanta, GA, USA e-mail: ale@mathcs.emory.edu12 Luca Gerardo-Giorda, Lucia Mirabella, Mauro Perego and Alessandro VenezianiThey coexist at each point and are connected by a cell membrane. The tissue con-ductivity depends upon its cells orientation, and in the most general case the associ-ated tensor is anisotropic [7, 14]. In any point x ∈Ω , where Ω is the spatial domainunder consideration, it is possible to identify an orthonormal triplet of directions,al(x), at(x), an(x), with al(x) parallel to the fibers direction, and we denote by σ li,e,σ ti,e, and σni,e the corresponding intra and extracellular conductivity coefficients. Theconductivity tensors are given byDτ(x) = σ lτ(x)al(x)aTl (x)+σtτ(x)at(x)aTt (x)+σnτ (x)an(x)aTn (x), τ = i,e. (1)We assume that Di,e fulfill in Ω a uniform elliptic condition.The Bidomain model. The Bidomain model is a nonlinear reaction-diffusion sys-tem of parabolic type describing the spatio-temporal dynamics of the intra and ex-tracellular potentials, denoted by ui and ue, while the cell membrane is regarded asdislocated in the domain [2]. We rely in this paper on a non-symmetric formulationin terms of the transmembrane potential u = ui− ue, and the extracellular one [4].We denote by u = (u,ue)T the unknown, and by lettingD = σleDiσ li +σleσ leDi−σ li Deσ li +σleDi Di +De E1 = [1 00 0]e1 =[10]the Bidomain system readsχCm E1∂u∂ t−∇ ·D∇u+χ Iion(u)e1 = Iapp (2)where Iion(u) is a nonlinear function of the transmembrane potential u, specified bya ionic model, and where Iapp represent the applied current stimuli. Several ionicmodels are available in literature, from more phenomenological to more accurateones, but the choice of the nonlinear term Iion(u) does not have any influence on theprocedure highlighted in what follows. The problem is completed by suitable initialconditions, and by homogeneous Neumann boundary conditions on ∂Ω , modelingan insulated myocardium. The transmembrane potential u is then uniquely deter-mined from (2), while the extracellular potential ue is determined up to a functionof time, whose value is usually obtained by imposing that∫Ωue = 0 .The Monodomain model. The Monodomain model is a simplified model for thepropagation of the electrical stimulus, based upon a proportionality assumption be-tween the intracellular and the extracellular conductivity tensors, namely assumingDe = λDi, where λ is a constant to be properly chosen. We assume here λ = σ le/σ li[6], and the Monodomain model readsχCm∂u∂ t−∇ · σleDiσ li +σle∇u+χIion(u) = Iapp, (3)Optimized Schwarz Methods and model adaptivity in electrocardiology 3still coupled with suitable initial conditions, and homogeneous Neumann boundaryconditions on ∂Ω . Differently from the Bidomain, the Monodomain model featuresa unique solution and is way cheaper to solve numerically. In absence of appliedcurrents, the Monodomain model is accurate enough to catch the desired dynamicsand effects of the action potential propagation [12]. However, the Bidomain modelbecomes necessary when current stimuli are applied in the extracellular space. Also,the Monodomain is inadequate to simulate defibrillation [16].1.1 Numerical approximationTime integration. For simplicity in presentation, we consider a fixed time step ∆ t,and we denote with superscript n the variables computed at time tn = n∆ t. Boththe Bidomain (2) and the Monodomain equations (3) are advanced in time by asemi-implicit scheme, where the nonlinear term (the ionic current) is evaluated atthe previous time steps. More precisely, moving from tn to tn+1 we solve in ΩχCm E1un+1−un∆ t−∇ ·D∇un+1 = Iapp−χ Iion(un)e1 (4)for the Bidomain system, andχCmun+1−un∆ t−∇ · σleDiσ li +σle∇un+1 = Iapp−χ Iion(un) (5)for the Monodomain one.Space discretization. Both Bidomain (4) and Monodomain (5) models are dis-cretized in space by finite elements. When solving the Bidomain system, the un-knowns of the fully discrete problem are represented by the vector(uh,ue,h)T , stor-ing the nodal values of the transmembrane and extracellular potentials. The matrixassociated with the discrete Bidomain models is given byB =[Buu BueBeu Bee]= χCm∆ t M+ σleσ li +σleKi σleσ li +σleKi−σ liσ li +σleKeKi Ki +Ke , (6)where M is the mass matrix and Ki,e are the stiffness matrices associated with thechosen finite elements space.When solving the Monodomain system, the unknown of the fully discrete problemis uh, and the associated matrix is simply block Buu of the Bidomain matrix.4 Luca Gerardo-Giorda, Lucia Mirabella, Mauro Perego and Alessandro Veneziani2 A model adaptive strategyStemming from the observation that the Monodomain error is concentrated aroundthe wavefront (particularly crosswind with respect to the fibers, see Figure 1), weconsider here a model adaptive strategy that solves the Bidomain only where ac-tually needed. In a first implementation of this approach [9] an intermediate Hy-bridomain model based on a suitable a posteriori model estimator was introduced.The Hybridomain assembles the block Bue only in correspondence with the nodesidentified as Bidomain ones by the model estimator, while the second equation staysuntouched. This simplifies significantly the implementation of the adaptive proce-dure, however the computational advantage is limited, since also in the Monodomainregions an extended problem with the same size of the Bidomain one is solved. Analternative procedure consists of a genuine heterogenous coupling by splitting thedomain. This coupling guarantees a significant reduction of the computational costs,but raises non trivial issues when matching the two models, featuring a differentsize. This has been considered in [6], where the Optimized Schwarz method hasbeen advocated for the heterogeneous coupling.From the practical standpoint, an immediate idea for setting up a model adap-tive solver would be to use the a posteriori error estimator for detecting the regionswhere to solve the Bidomain problem and then to couple these subdomains with theMonodomain regions, so to have the Bidomain model capturing the wavefron prop-agation. However, this approach is barely doable. As a matter of fact. the Robin-type interface conditions in the Optimized Schwarz setting require the assembly ofmass matrices on the interfaces. As a consequence, every time the Bidomain regionchanges, one should identify the new interfaces and then recompute the matrices,with an additional computational cost that is anticipated to reduce the advantage ofthe Optimized Schwarz coupling. The model adaptive strategy we propose here re-lies instead on a a priori subdivision of Ω into smaller subdomains Ω j. The modelerror estimator will associate runtime each subdomain with either the Bidomain orthe Monodomain problem. In this way, the interfaces matrices needed for the cou-pling can be computed once at the beginning of the time loop. The non-symmetricformulation of the Bidomain system ensures that by assembling the blocks of everysubdomain matrix, also the Monodomain ones are available. Details on the Opti-mized Schwarz coupling are given in the next Section.Fig. 1 Error in the propagation of the membrane potential u between Bidomain and Monodomainsimulation (from [6]): fibers oriented along the x axis.Optimized Schwarz Methods and model adaptivity in electrocardiology 52.1 Coupling conditions and Optimized Schwarz MethodsWe outline here the coupling conditions for the three different types of interfaces. Ifthe subdomains involved have the same characteristic (Bido/Bido and Mono/Mono)the corresponding solutions are labeled by subscript 1 and 2, while if the subdo-mains have different characteristics (Bido/Mono) the corresponding solutions arelabeled with subscript B and M.Bidomain/Bidomain interface. The coupling conditions on the Bidomain/Bidomaininterface have been introduced in [5], and are given bynT1 D∇u1 +α1 Σ u1 = nT1 D∇u2 +α1 Σ u2nT2 D∇u2 +α2 Σ u2 = nT2 D∇u1 +α2 Σ u1, where Σ = σ leσ li +σ le 01 σli +σleσ li . (7)The convergence of the Optimized Schwarz Algorithm based on the interface con-ditions (7) was analyzed in [5], where also optimal parameters have been identifiedby means of Fourier analysis.Bidomain/Monodomain interface. Due to a dimensional mismatch between thetwo models, two interface conditions are needed on the Bidomain side of the inter-face, and one on the Monodomain side [6]. Possible coupling conditions arenTBσ leDiσ li +σle(∇uB +∇ue,B)−nTBσ li Deσ li +σle∇ue,B +σ leασ li +σleuB = nTBσ leDiσ li +σle∇uM +σ leασ li +σleuBnTBDi(∇uB +∇ue,B)+nTBDe∇ue,B +αuB +σ li +σleσ liαue,B = αurest(8)for the Bidomain subproblem, andnTMσ leDiσ li +σle∇uM +σ leασ li +σleuM = nTMσ leDiσ li +σle(∇uB +∇ue,B)−nTMσ li Deσ li +σle∇ue,B +σ leασ li +σleuB(9)for the Monodomain one. To cope with the mismatch, the second condition in (8)is a transparent boundary condition, designed to avoid spurious reflexions off theinterface for the extracellular potential wave. The convergence of the OptimizedSchwarz Algorithm based on the interface conditions (8)-(9) was analyzed in [6],where also optimal parameters has been identified by means of Fourier analysis.Monodomain/Monodomain interface. The Optimized Schwarz coupling is waysimpler on the interface between two Monodomain regions. The Monodomainmodel is a parabolic equation equipped with a nonlinear reaction term. The IMEXtemporal integration scheme reduces it, within the single time step, to the solution of6 Luca Gerardo-Giorda, Lucia Mirabella, Mauro Perego and Alessandro Veneziania linear steady reaction-diffusion problem, whose solution by means of OptimizedSchwarz Methods has been widely studied, and an optimal parameter has been iden-tified [3]. The coupling on the interface is thus given bynT1σ leDiσ li +σle∇u1 +αoptu1 = nT1σ leDiσ li +σle∇up2 +αoptup2nT2σ leDiσ li +σle∇u2 +αoptu2 = nT2σ leDiσ li +σle∇u1 +αoptu1.(10)2.2 The model error estimatorFor choosing between a Bidomain or Monodomain simulation within the singlesubdomain we need a reliable model error estimator. An efficient a posteriori errorestimator was introduced in [9], based on the extracellular potential computed froman extended Monodomain solution. To further reduce the computational cost, wecompute only a Monodomain solution uM over the whole domain.We let De = σleσ liDi+Dε , and the model estimator is computed at the subdomain levelasζ2j =∫ tn+1tn∫Ω j∇uMσ li Dεσ li +σle(D−1i +D−1e) σ li Dεσ li +σle∇uM dxdt.The value ζ 2j is an upper bound for the error between the Monodomain and theBidomain model within the Ω j [9]. The Bidomain model is then activated in Ω jwhenever ζ 2j exceeds a given threshold τ j, depending on the size of the subdomain.2.3 The model adaptive algorithmWe summarize here the resulting model adaptive strategy, by highlighting the pre-processing and the runtime phases.Preprocessing1. Split the computational domain into non-overlapping subregions Ω j ( j = 1, ..,N).2. Identify the interfaces Γi j between subdomains Ωi and Ω j.3. Assemble the local matrices B juu, B jue, B jeu, B jee4. Assemble the interface mass matrices MΓi j5. Compute the incomplete ILU factorization of the local B juu and B jee matricesRuntime (time step tn→ tn+1)1. Run a Monodomain simulation at time tn+1 over the whole domain Ω2. Evaluate the a-posteriori model estimator and compute the local indicatorOptimized Schwarz Methods and model adaptivity in electrocardiology 73. For all subdomains where the local indicator is above a given threshold, activateBidomain4. Run a few iteration of the Optimized Schwarz Algorithm using the solution com-puted in Step 1. as initial guess5. Advance to time tn+1A dynamical allocation of tasks is under investigation to properly balance, in realproblems, the load of each processor in the parallel solver.3 Preliminary numerical resultsNumerical results in this section have the purpose to show the effectiveness of themodel adaptive method: for this reason we consider here only 2D simulations. Thenumerical tests are run in Matlab R© 7.5. The Bidomain problems are solved withby a flexible GMRES (f-GMRES) right preconditioned by the Block-triangular pre-conditioner introduced in [4], while the Monodomain problems are solved by a CGpreconditioned by an ILU factorization.We consider the strip Ω = [0,3]× [0,1] subdivided into the three noverlapping sub-domains Ω1 = [0,1]× [0,1], Ω2 = [1,2]× [0,1], and Ω3 = [2,3]× [0,1]. The fibersare oriented with the principal direction perpendicular to the interfaces, and we im-pose a stimulus on the whole left boundary of Ω1. The ionic model is the Rogers-McCulloch one [13]. We plot in Figure 2 the wavefront position at different times(top row), and the activated subdomains (bottom row).A more detailed presentation of the method will be the subject of a forthcomingwork. Further work needs to be done to identify the proper trade-off between thenumber of subdomains, and the size of the Bidomain region surrounding the wave-front, and to properly handle the processors load balance in a parallel architecture.References1. R. H. Clayton, O. M. Bernus, E. M. Cherry, H. Dierckx, F. H. Fenton, L. Mirabella, A. V.Panfilov, F. B. Sachse, G. Seemann, H. Zhang, Models of cardiac tissue electrophysiology:Progress, challenges and open questions, Progress in biophysics and molecular biology 104,pp. 22-48, 2011.2. P. Colli Franzone, L. Pavarino, G. Savaré. Computational electrocardiology: mathematicaland numerical modeling, in Complex Systems in Biomedicine - A. Quarteroni, L. Formaggia,and A. Veneziani, editors. Springer, Milan, 2006.3. M.J. Gander. Optimized Schwarz methods. SIAM J. Num. Anal., 44(2), pp. 699–731, 2006.4. L. Gerardo-Giorda, L. Mirabella, F. Nobile, M. Perego, and A. Veneziani. A model-basedblock-triangular preconditioner for the Bidomain system in electrocardiology. J. Comp. Phys.,228, pp. 3625–3639, 2009.5. L. Gerardo-Giorda and M. Perego, Optimized Schwarz Methods for the Bidomain system inelectrocardiology M2AN, September 2012, DOI: http://dx.doi.org/10.1051/m2an/2012040.8 Luca Gerardo-Giorda, Lucia Mirabella, Mauro Perego and Alessandro VenezianiFig. 2 Propagation of the membrane potential (in red the excited region, top row), the subdomainsΩ1, Ω2, and Ω3 (middle row), and the activated Bidomain region (in red, bottom row).6. L. Gerardo-Giorda, M. Perego, and A. Veneziani. Optimized Schwarz coupling of Bidomainand Monodomain models in electrocardiology. M2AN, Vol. 45 (2), pp. 309-334, 2011.7. J. Le Grice, B.H. Smaill, and P.J. Hunter. Laminar structure of the heart: a mathematicalmodel. Am. J. Physiol., 272 (Heart Circ. Physiol.)(41):H2466–H2476, 1995.8. G.T. Lines, M.L. Buist, P. Grottum, A.J. Pullan, J. Sundnes, and A. Tveito. Mathematicalmodels and numerical methods for the forward problem in cardiac electrophysiology Comput.Visual. Sci,5:215-239, 2003.9. L. Mirabella, F. Nobile, and A. Veneziani. An a posteriori error estimator for model adaptivityin electrocardiology. Comp. Meth. Appl. Mech. Engrg., (doi:10.1016/j.cma.2010.03.009),2010 (in press).10. L. F. Pavarino and S. Scacchi. Multilevel additive Schwarz preconditioners for the Bidomainreaction-diffusion system. SIAM J. Sci. Comp., 31(1):420–443, 2008.11. M. Pennacchio and V. Simoncini. Algebraic multigrid preconditioners for the bidomain reac-tiondiffusion system. Appl. Num. Math., 59(12):3033–3050, 2009.12. M. Potse, B. Dubé, J. Richer, and A. Vinet. A comparison of Monodomain and BidomainReaction-Diffusion models for Action Potential Propagation in the Human Heart. IEEETrans. Biomed. Eng., 53(12):2425–2435, 2006.13. J. Rogers and A. McCulloch. A collocation-Galerkin finite element model of cardiac actionpotential propagation. IEEE Transactions on Biomedical Engineering, 41:743–757, 1994.14. F. B. Sachse. Computational Cardiology. Springer, Berlin, 2004.15. E.J. Vigmond, R. Weber dos Santos, A.J. Prassl, M. Deo, and G. Plank. Solvers for the caridacbidomain equations. Progress in Biophysics and Molecular Biology, 96(1-3):3–18, 2008.16. N. Trayanova. Defibrillation of the heart: insights into mechanisms from modelling studies.Experimental Physiology, 91: 323–337, 2006

Optimized schwarz methods and model adaptivity in electrocardiology simulations

https://bird.bcamath.org/bitstream/20.500.11824/272/1/GerardoGiorda_MS19.pdf

Optimized schwarz methods and model adaptivity in electrocardiology simulations

Abstract

Similar works

Full text

Available Versions

Crossref

BCAM's Institutional Repository Data