BioMed Research International

Volume 2015, Article ID 473279, 15 pages

http://dx.doi.org/10.1155/2015/473279

## Simulating Cardiac Electrophysiology Using Unstructured All-Hexahedra Spectral Elements

^{1}CRS4, Loc. Pixina Manna, Edificio 1, 09010 Pula, Italy^{2}Fujitsu Laboratories of Europe, Hayes Park Central, Hayes End Road, Hayes, Middlesex UB4 8FE, UK

Received 14 December 2014; Revised 20 March 2015; Accepted 9 April 2015

Academic Editor: Joakim Sundnes

Copyright © 2015 Gianmauro Cuccuru et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

#### Abstract

We discuss the application of the spectral element method to the monodomain and bidomain equations describing propagation of cardiac action potential. Models of cardiac electrophysiology consist of a system of partial differential equations coupled with a system of ordinary differential equations representing cell membrane dynamics. The solution of these equations requires solving multiple length scales due to the ratio of advection to diffusion that varies among the different equations. High order approximation of spectral elements provides greater flexibility in resolving multiple length scales. Furthermore, spectral elements are extremely efficient to model propagation phenomena on complex shapes using fewer degrees of freedom than its finite element equivalent (for the same level of accuracy). We illustrate a fully unstructured all-hexahedra approach implementation of the method and we apply it to the solution of full 3D monodomain and bidomain test cases. We discuss some key elements of the proposed approach on some selected benchmarks and on an anatomically based whole heart human computational model.

#### 1. Introduction

Mechanistic models of heart physiology are emerging as an important paradigm able to complement the insight and exploit the information provided by genetics and biological experiments. Very recent years have seen the first attempts to define models for the numerical simulation of the individual response of real cardiology patients. For instance, in [1] the clinical study of a 65-year-old male with ischemic cardiomyopathy was addressed by means of a workflow consisting of detailed anatomic reconstruction and effective modelling of electrical activity, biomechanics, and hemodynamics. The modelling steps are challenging because they require performing (i) multiple simulations to estimate parameters and test outcomes of different treatment options and (ii) many large scale computer analyses in a clinically useful time-frame. These are strong motivations in favour of effective, anatomy-compliant high performance computer methods.

In particular, current software to handle the forward problem of the electrocardiology, that is, modelling the electrical activity of the heart, largely fails to deliver the performance needed for clinical applications. Typical discretization techniques applied to this problem have an average nodal spacing of 0.2 mm or less, resulting in approximately 30 million grid points for an average-size human heart. This possibly gives rise to large scale computational problems which are difficult to solve in real-time. The most common approach relies on the use of large high performance computer systems featuring to cores [2, 3], a choice often beyond the possibility of typical end-users. Several strategies have been suggested for facing this hindrance, including simplifying/optimizing the cardiac cell models [4], code optimization for run on next-generation hardware devices [5], and adoption of more effective solutions for the approximations of the differential equations at hand. This work is based on the latter approach.

#### 2. Related Work

The bidomain equations [6] are a commonly used model for the propagation of depolarization wavefronts across cardiac tissue. They are often solved numerically using finite elements (FE), finite volumes or finite differences, and a variety of algorithms based on these techniques have been proposed [7–14]. However, each of these methods suffers from a well-known limit affecting the numerical simulation of propagation phenomena: to prevent the onset of nonphysical, spurious effects commonly referred to as numerical dispersion it is necessary to properly fix the accuracy of the method. In cardiac simulation this has most commonly been achieved by using grid refinement, that is, choosing a node spacing that is sufficiently fine that numerical dispersion is not a problem.

Alternatively, one can improve the quality of numerical simulations by using a relatively coarse mesh, increasing the polynomial degree of the basis functions used in the numerical scheme. This is the case of high order methods, examples of which are the -version of the finite element method and the spectral element method. High order methods possess many advantages over standard techniques, such as considerably higher convergence rates, and, in our experience, the possibility of using elements with large aspect ratios without significant deterioration in accuracy. Designed to combine the good accuracy properties of pseudospectral techniques such as Legendre or Chebyshev methods with the geometrical flexibility of classical low order FE methods, spectral element methods were first introduced in 1984 [15], combining spectral methods with domain decomposition approaches. During the last decades, spectral elements have firmly established as an effective tool to treat diverse compute-intensive physical problems on complex geometries in fluid dynamics, continuum mechanics, geophysics, and electromagnetics; see, for example, [16] and references therein.

As for our knowledge, few attempts have been done to exploit spectral elements in electrocardiology. An interesting idea is provided by the so-called spectral smoothed boundary method [17], which tries to merge the excellent convergence property of Fourier-type spectral methods (which, on the other hand, are not able to accommodate irregular geometries) with the phase-field method, an embedding approach capable of dealing with complex domains when no-flux boundary conditions are present. In [18], high order spectral/hp elements on curvilinear triangles are used to model electrical propagation in the heart surface, addressing the monodomain problem and including anisotropic heterogeneous propagation in presence of fibre orientation. In [19], which in the authors’ words “should be seen as a proof-of-concept,” high order finite elements based on triangles have been applied to the solution of the monodomain problem in 1D and 2D reference geometries: furthermore, authors provide rigorous analysis and numerical experiments on theoretical error convergence rates and suggest the possibility to extend their method to spatial adaptivity.

When solving numerically the electrocardiology problem, extremely high accuracy is not a valuable goal. To be realistic, none of the numerous parameters involved in the mono-/bidomain model is known with a precision higher than 10%. On a logical ground, there is no use in solving with outstanding accuracy a problem based on a model which is only a partial approximation of a real patient cardiovascular system. In this frame, the rationale behind the use of spectral elements is the possibility to achieve an a priori fixed, reasonable accuracy using less degrees of freedom (DOF) with respect to low order methods. Matrices resulting from space discretization are more dense but have a smaller overall number of nonzero entries: this is important for decreasing both memory footprints and execution times. On the other hand, a reduced number of DOFs mean smaller ODE systems describing cell models, possibly shrinking the CPU-time for the* in silico* analysis.

Based on our experience in several fields of engineering, we propose to use straightforward, fully unstructured all-hexahedra spectral elements (SEM in the sequel). Differently from high order methods based on triangular elements, the SEM rely on a plain, standard, consolidated, widely accepted mathematical formulation based on tensor product, which brings tangible advantages in practical algorithm implementation. Comparison with triangular spectral elements is a difficult task as the latter came up in many different algorithm implementations and most of the works on theory and numerical experiments concern the 2D case. For instance, in [20] Koornwinder-Dubiner polynomials and Fekete points were chosen as orthogonal basis and approximation points, respectively. This choice brings differentiation matrices greater in size with respect to their SEM counterpart, thus raising the computational effort for matrix-vector products and, consequently, the overall CPU-time. Remarkably, numerical experiments suggest that the condition number of triangular spectral matrices is worse than the SEM ones, with clear advantage for the latter in terms of convergence of iterative methods and error propagation. Furthermore, preconditioners for SEM have been long studied and developed; see, for instance [21] and references therein. Apart from these good mathematical properties, simplicity of SEM results in easy to implement, readable computer codes enjoying good portability to high performance platforms and new hardware devices like GPGPUs (as, for instance, there is no need for complex data-structures). Also, since they share a common philosophy with traditional finite elements, they could benefit from several existing and well-consolidated tools (pre- and postprocessing routines, mesh partitioning, iterative solvers, preconditioners, etc.).

The other side of the picture is that all-hexahedral grids are more difficult to generate than simplicial tetrahedral grids. While automatic (or semiautomatic) hex-only mesh generation for complex geometries is an active subject of research, generation of an all hexahedral mesh for applications in electrophysiology requires nowadays more effort than hybrid or all-tet grids. This is particularly true for those studies where anatomical properties should be resolved at a very fine anatomical detail. Clearly, the extent to which a study requires representation of structural detail is determined by the particular application and there might be some cases where this is not a serious limitation. Finally, the spectral element method implementation we propose is restricted to 8-node hexahedra. Extension to (curved) isoparametric elements is currently being addresses and will allow mapping curved boundaries of irregularly shaped anatomical details and represent geometries more closely.

#### 3. Materials and Methods

##### 3.1. The Bidomain Model

For a 3D domain the bidomain equations in parabolic-parabolic form read [22] where and are the intracellular and extracellular potentials, is the transmembrane potential, is the membrane capacitance per unit area, and are the intracellular and extracellular conductivity tensors modeling the anisotropy of the cardiac tissue, and is the cell surface area to volume ratio. Here, and are arrays of ionic gating variables and ionic concentrations, respectively [22], and and are externally applied intracellular and extracellular stimulus currents. Suitable boundary conditions on either or the current flux , and initial conditions for , , and are set. Functions , and the transmembrane current are determined by parameters and data fitting and represent an electrophysiological cell model. Many of such cells models of variable complexity exist; see, for example, [23] for a recent review. In this work we model the ionic current using the phase-I Luo-Rudy (LR1) cell model [24]. However, the derivation of the SEM formulation of the bidomain model follows analogously for any ordinary differential equation (ODE) based cell model. The LR1 model relies on six gating variables along with the intracellular calcium concentration . The compatibility conditionon and should hold for the system to be solvable. Note that the transmembrane potential is uniquely determined, while the intra- and extracellular potentials and are determined up to the same additive time-dependent constant, whose value is usually obtained by imposing a normalization condition on , for example, zero average on .

We remark that a form alternative to (1a)-(1b), the parabolic-elliptic formulation of the bidomain equations, exists, with unknowns . Such formulation is quite popular, even because it allows dealing with the parabolic and elliptic blocks at different steps by Gauss-Seidel method. Nevertheless, (i) recent studies show that decoupling the bidomain equations can be less efficient than solving them as a global system (for the same level of accuracy [25]), and (ii) evidence exists that formulation provides a significantly more efficient numerical framework [26]. Our SEM implementation addresses the coupled parabolic-parabolic problem, even because, being the SEM mass matrix diagonal, this allows a significant save in memory allocation and, quite likely, a reduced computational effort. This is because the number of nonzero entries of the global matrix (hence the cost for sparse matrix-vector products) is greatly reduced with respect to the parabolic-elliptic formulation.

##### 3.2. Variational Formulation of the Bidomain Model

The variational form of the bidomain equations is as follows: given and fulfilling the compatibility condition (2), for all , find satisfying given initial and boundary conditions such that for all admissible test functions from a suitable test space, where . See [27] for a mathematical analysis of the bidomain model. Existence and uniqueness for a solution of the bidomain problem for a wide class of cell models, including LR1, are discussed in [28].

##### 3.3. Spatial Discretization

Let be a decomposition of the model volume into a family of nonoverlapping hexahedra with typical size . Each element is obtained by a transformation from a reference (or parent) element . On the reference element, we consider the space of polynomial functions with degree less than or equal to in each variable , . Then, for each letbe the space of the functions that are images through of polynomials functions in . Hence, for all in . Finally, we define the spectral element space asIn our implementation ; thus, when the mapping is subparametric, meaning that its degree is lower than the spectral degree. This choice is essentially motivated by practical considerations: mesh generators produce 8-point hexahedra, while -order hexahedra need total points (and produce different grids for analyses with different spectral degree). The subparametric mapping is known to enjoy good mathematical properties (see, for instance, [29]). Nevertheless, if high order hexahedra are desirable (for instance, when dealing with domains with curved boundary which can be described in terms of geometrical primitives), they can be incorporated in the spectral element frame with small additional effort.

##### 3.4. Construction of SEM Basis Functions

First, the Gauss-Lobatto (LGL) points in the reference element are obtained via tensor product of the one-dimensional LGL nodes in . The full spectral grid is then built mapping the LGL nodes over the hexahedra and eliminating duplicated points: as for finite elements, a global numbering is associated with the grid points. In we consider the Lagrange polynomials of degree corresponding to the LGL nodes and, in , the elemental basis functions for obtained by mapping such nodal basis functions according to . Finally, a basis for the whole space is obtained as a patchwork of these elemental functions on each element . More precisely, we choose as a basis for the set , that is, the transformation of polynomials of order from to each , such thatwhere is the Kronecker delta. It is easy to see that the restriction of such spectral shape functions to either coincides with a Lagrange polynomial or vanishes. It is easily verified that the above choice of nodal basis functions assures that the corresponding global bases enjoy as much localization as possible.

##### 3.5. SEM Formulation

Taking the test space equal to , the SEM approximation of (3a)-(3b) consists in finding such that for all and for every where , while and denotes the inner product in . For the numerical realization of the SEM we resort to the so-called GNI approach (see, e.g., [16]), which consists in exploiting numerical integration for computing integrals. With this choice, the semidiscrete form of the bidomain equations becomes the following: find satisfying given initial and boundary conditions and such that for all where the notation indicates that the integrals on are computed using Gauss-Lobatto numerical integration formulas on the corresponding reference element .

###### 3.5.1. Semidiscrete Form of the Bidomain Equation

In terms of the SEM basis, the spectral element solution is where represent the unknown nodal value. If we denote by the diagonal spectral element mass matrix and with the symmetric intra and extracellular stiffness matrices; with elementswe may write the semidiscrete form of the bidomain equations in compact matrix form as follows:orwhere the bidomain mass and stiffness matrices and are defined asIn a similar, but easier way, one can show that the SEM formulation of the monodomain problem reads along with given cell model, boundary and initial conditions.

##### 3.6. Semi-Implicit Time Discretization

We use a mixed time-marching scheme which is implicit for the intracellular concentration variables while it is explicit for and . This is motivated by the fact that the components of in (1c) are fully decoupled and each depends only on itself and , whereas the component of has complex, non-linear dependencies on , on one another and on , making an implicit time-step much more expensive to compute.

*Step 1. *Given the potential at the previous time , solvewhere denotes the time step size.

*Step 2. *Find and by solving the linear systemwhere with and . For the solution of (15) we exploit classical Krylov methods, taking into account the linear system consistency: see, for instance, [30, Theorem 2].

Clearly, the effective solution of (15) requires the choice of an optimal preconditioner, especially for the bidomain case. This is a difficult task which truly deserves a dedicated investigation. In particular, algebraic multigrid preconditioners (AMG) have been applied to the bidomain equations, often outperforming standard methods like the ILU or the block Jacobi preconditioner adopted in this work [31]. Quite often, AMG preconditioners are applied to the parabolic-elliptic formulation of the bidomain equations, while our SEM implementation stems from the coupled parabolic-parabolic problem (see Section 3.1). Analysis of AMG preconditioning for the SEM discretization of the coupled parabolic-parabolic problem (3a)-(3b) will be the object of a future paper, along with a study on the condition numbers of SEM matrices, both unpreconditioned and preconditioned with different methods.

#### 4. Results and Discussion

In order to validate the proposed spectral element scheme we start analyzing its performance for monodomain and bidomain simulations using a well know benchmark on a simplified geometry. Then, to further assess the ability of the method to cope with the additional complexities associated with organ scale cardiac electrophysiology models we present the results of simulating monodomain activity on a realistic human heart geometry. Numerical experiments were performed at the CRS4 High Performance Computing centre comprising more than 3000 processors, grouped into 382 nodes of 8 Intel Xeon (2.8 GHz) cores, connected through a low latency Infiniband network.

##### 4.1. The Monodomain Niederer Benchmark

We use the benchmark proposed in [32]: it is a consensus cardiac tissue electrophysiology model problem that is used for evaluating and verifying cardiac tissue electrophysiology simulators. The problem geometry is defined as a parallelepipedal portion of cardiac tissue (*cuboid* or* slab*) with dimensions of cm (see Figure 1), characterized by transversely isotropic conductivity with fibres aligned to the long axis of the computational domain.