A parallel solver for reaction-diffusion systems in computational electrocardiology

In this work, a parallel three-dimensional solver for numerical simulations in computational electrocardiology is introduced and studied. The solver is based on the anisotropic Bidomain %(AB) cardiac model, consisting of a system of two degenerate parabolic reaction-diffusion equations describing the intra and extracellular potentials of the myocardial tissue. This model includes intramural fiber rotation and anisotropic conductivity coefficients that can be fully orthotropic or axially symmetric around the fiber direction. %In case of equal anisotropy ratio, this system reduces to The solver also includes the simpler anisotropic Monodomain model, consisting of only one reaction-diffusion equation. These cardiac models are coupled with a membrane model for the ionic currents, consisting of a system of ordinary differential equations that can vary from the simple FitzHugh-Nagumo (FHN) model to the more complex phase-I Luo-Rudy model (LR1). The solver employs structured isoparametric $Q_1$ finite elements in space and a semi-implicit adaptive method in time. Parallelization and portability are based on the PETSc parallel library. Large-scale computations with up to $O(10^7)$ unknowns have been run on parallel computers, simulating excitation and repolarization phenomena in three-dimensional domains

Modeling Ventricular Excitation : axial and orthotropic anisotropy effects on wavefronts and potentials.

By applying the eikonal approximation to the bidomain model of the cardiac tissue we investigate the influence of the axially isotropic and orthotropic conductivity tensors on the propagation of the excitation wave fronts and on the associated potential distribution and electrograms

Computational modeling of the electromechanical response of a ventricular fiber affected by eccentric hypertrophy

AbstractThe aim of this work is to study the effects of eccentric hypertrophy on the electromechanics of a single myocardial ventricular fiber by means of a one-dimensional finite-element strongly-coupled model. The electrical current ow model is written in the reference configuration and it is characterized by two geometric feedbacks, i.e. the conduction and convection ones, and by the mechanoelectric feedback due to stretchactivated channels. First, the influence of such feedbacks is investigated for both a healthy and a hypertrophic fiber in case of isometric simulations. No relevant discrepancies are found when disregarding one or more feedbacks for both fibers. Then, all feedbacks are taken into account while studying the electromechanical responses of fibers. The results from isometric tests do not point out any notable difference between the healthy and hypertrophic fibers as regards the action potential duration and conduction velocity. The length-tension relationships show increased stretches and reduced peak values for tension instead. The tension-velocity relationships derived from afterloaded isotonic and quick- release tests depict higher values of contraction velocity at smaller afterloads. Moreover, higher maximum shortenings are achieved during the isotonic contraction. In conclusion, our simulation results are innovative in predicting the electromechanical behavior of eccentric hypertrophic fibers

Parallel multilevel solvers for the cardiac electro-mechanical coupling

We develop a parallel solver for the cardiac electro-mechanical coupling. The electric model consists of two non-linear parabolic partial differential equations (PDEs), the so-called Bidomain model, which describes the spread of the electric impulse in the heart muscle. The two PDEs are coupled with a non-linear elastic model, where the myocardium is considered as a nearly-incompressible transversely isotropic hyperelastic material. The discretization of the whole electro-mechanical model is performed by Q1 finite elements in space and a semi-implicit finite difference scheme in time. This approximation strategy yields at each time step the solution of a large scale ill-conditioned linear system deriving from the discretization of the Bidomain model and a non-linear system deriving from the discretization of the finite elasticity model. The parallel solver developed consists of solving the linear system with the Conjugate Gradient method, preconditioned by a Multilevel Schwarz preconditioner, and the non-linear system with a Newton\u2013Krylov-Algebraic Multigrid solver. Three-dimensional parallel numerical tests on a Linux cluster show that the parallel solver proposed is scalable and robust with respect to the domain deformations induced by the cardiac contraction

Role of infarct scar dimensions, border zone repolarization properties and anisotropy in the origin and maintenance of cardiac reentry

Cardiac ventricular tachycardia (VT) is a life-threatening arrhythmia consisting of a well organized structure of reentrant electrical excitation pathways. Understanding the generation and maintenance of the reentrant mechanisms, which lead to the onset of VT induced by premature beats in presence of infarct scar, is one of the most important issues in current electrocardiology. We investigate, by means of numerical simulations, the role of infarct scar dimension, repolarization properties and anisotropic fiber structure of scar tissue border zone (BZ) in the genesis of VT. The simulations are based on the Bidomain model, a reaction-diffusion system of Partial Differential Equations, discretized by finite elements in space and implicit-explicit finite differences in time. The computational domain adopted is an idealized left ventricle affected by an infarct scar extending transmurally. We consider two different scenarios: i) the scar region extends along the entire transmural wall thickness, from endocardium to epicardium, with the exception of a BZ region shaped as a central sub-epicardial channel (CBZ); ii) the scar region extends transmurally along the ventricular wall, from endocardium to a sub-epicardial surface, and is surrounded by a BZ region (EBZ). In CBZ simulations, the results have shown that: i) the scar extent is a crucial element for the genesis of reentry; ii) the repolarization properties of the CBZ, in particular the reduction of IKs and IKr currents, play an important role in the genesis of reentrant VT. In EBZ simulations, since the possible reentrant pathway is not assigned a-priori, we investigate in depth where the entry and exit sites of the cycle of reentry are located and how the functional channel of reentry develops. The results have shown that: i) the interplay between the epicardial anisotropic fiber structure and the EBZ shape strongly affects the propensity that an endocardial premature stimulus generates a cycle of reentry; ii) reentrant pathways always develop along the epicardial fiber direction; iii) very thin EBZs rather than thick EBZs facilitate the onset of cycles of reentry; iv) the sustainability of cycles of reentry depends on the endocardial stimulation site and on the interplay between the epicardial breakthrough site, local fiber direction and BZ rim

A comparison of coupled and uncoupled solvers for the cardiac Bidomain model

The aim of this work is to compare a new uncoupled solver for the cardiac Bidomain model with a usual coupled solver. The Bidomain model describes the bioelectric activity of the cardiac tissue and consists of a system of a non-linear parabolic reaction-diffusion partial differential equation (PDE) and an elliptic linear PDE. This system models at macroscopic level the evolution of the transmembrane and extracellular electric potentials of the anisotropic cardiac tissue. The evolution equation is coupled through the non-linear reaction term with a stiff system of ordinary differential equations (ODEs), the so-called membrane model, describing the ionic currents through the cellular membrane. A novel uncoupled solver for the Bidomain system is here introduced, based on solving twice the parabolic PDE and once the elliptic PDE at each time step, and it is compared with a usual coupled solver. Three-dimensional numerical tests have been performed in order to show that the proposed uncoupled method has the same accuracy of the coupled strategy. Parallel numerical tests on structured meshes have also shown that the uncoupled technique is as scalable as the coupled one. Moreover, the conjugate gradient method preconditioned by Multilevel Hybrid Schwarz preconditioners converges faster for the linear systems deriving from the uncoupled method than from the coupled one. Finally, in all parallel numerical tests considered, the uncoupled technique proposed is always about two or three times faster than the coupled approach

Reduced-order modeling for cardiac electrophysiology. Application to parameter identification

A reduced-order model based on Proper Orthogonal Decomposition (POD) is proposed for the bidomain equations of cardiac electrophysiology. Its accuracy is assessed through electrocardiograms in various configurations, including myocardium infarctions and long-time simulations. We show in particular that a restitution curve can efficiently be approximated by this approach. The reduced-order model is then used in an inverse problem solved by an evolutionary algorithm. Some attempts are presented to identify ionic parameters and infarction locations from synthetic ECGs.Comment: No. RR-7811 (2011

Isogeometric approximation of cardiac electrophysiology models on surfaces: An accuracy study with application to the human left atrium

We consider Isogeometric Analysis in the framework of the Galerkin method for the spatial approximation of cardiac electrophysiology models defined on NURBS surfaces; specifically, we perform a numerical comparison between basis functions of degree p ≥ 1 and globally C k -continuous, with k = 0 or p − 1, to find the most accurate approximation of a propagating front with the minimal number of degrees of freedom. We show that B-spline basis functions of degree p ≥ 1, which are C p−1 -continuous capture accurately the front velocity of the transmembrane potential even with moderately refined meshes; similarly, we show that, for accurate tracking of curved fronts, high-order continuous B-spline basis functions should be used. Finally, we apply Isogeometric Analysis to an idealized human left atrial geometry described by NURBS with physiologically sound fiber directions and anisotropic conductivity tensor to demonstrate that the numerical scheme retains its favorable approximation properties also in a more realistic setting