Simulation of cardiac electrophysiology on next-generation high-performance computers

Rafel Bordas Bruno Carpentieri Giorgio Fotia Fabio Maggio Ross Nobes () Joe Pitt-Francis James Southern Articles on similar topics can be found in the following collections computational biology (52 articles) computer modelling and simulation (82 articles) Receive free email alerts when new articles cite this article - sign up in the box at the top right-hand corner of the article or click here - Email alerting service Simulation of cardiac electrophysiology on next-generation high-performance computers 1Oxford University Computing Laboratory, University of Oxford, Wolfson Building, Parks Road, Oxford OX1 3QD, UK 2CRS4 Bioinformatica, 09010 Pula, Italy 3Fujitsu Laboratories of Europe Ltd, Hayes, Middlesex UB4 8FE, UK Models of cardiac electrophysiology consist of a system of partial differential equations (PDEs) coupled with a system of ordinary differential equations representing cell membrane dynamics. Current software to solve such models does not provide the required computational speed for practical applications. One reason for this is that little use is made of recent developments in adaptive numerical algorithms for solving systems of PDEs. Studies have suggested that a speedup of up to two orders of magnitude is possible by using adaptive methods. The challenge lies in the efficient implementation of adaptive algorithms on massively parallel computers. The finite-element (FE) method is often used in heart simulators as it can encapsulate the complex geometry and small-scale details of the human heart. An alternative is the spectral element (SE) method, a high-order technique that provides the flexibility and accuracy of FE, but with a reduced number of degrees of freedom. The feasibility of implementing a parallel SE algorithm based on fully unstructured all-hexahedra meshes is discussed. A major computational task is solution of the large algebraic system resulting from FE or SE discretization. Choice of linear solver and preconditioner has a substantial effect on efficiency. A fully parallel implementation based on dynamic partitioning that accounts for load balance, communication and data movement costs is required. Each of these methods must be implemented on nextgeneration supercomputers in order to realize the necessary speedup. The problems that this may cause, and some of the techniques that are beginning to be developed to overcome these issues, are described. Keywords: cardiac simulation; high-performance computing; finite elements; spectral elements; adaptive mesh refinement; linear algebra 1. Introduction The bidomain and monodomain equations (Plonsey 1988; Sepulveda et al. 1989; Keener & Sneyd 1998) have been widely used for many years to model the propagation of electrical waves across cardiac tissue. The numerical solution of One contribution of 15 to a Theme Issue The virtual physiological human: tools and applications I. the bidomain equations is usually calculated using standard finite-element (FE), finite-volume (FV) or finite-difference (FD) methods. However, as the discretization of a whole heart with an average nodal spacing of 0.2 mm generates a mesh with many millions of nodes, the linear systems resulting from FD or FE methods are very large. Whole-heart simulation using the bidomain model is therefore a demanding scientific computing problem. Several advanced problem-solving environments for cardiac simulation already exist (Vigmond et al. 2003; Watanabe et al. 2004; Pitt-Francis et al. 2008; http:// www.cmiss.org; http://www.continuity.ucsd.edu; http://www1.pacific.edu/ jeason), but none of these currently provide the required computational speed to allow researchers in, for example, the pharmaceutical industry to run simulations quickly enough to make them viable as an everyday research tool. For example, using CARP, 6.4 hours of CPU time is required on a 64-processor machine to run a 200 ms simulation of the bidomain model on a rabbit ventricular mesh containing 862 515 grid points using a modified BeelerReuter ionic model (Beeler & Reuter 1977; Plank et al. 2007). Extending this to the approximately 30 million grid points that would be required to compute an average (human-sized) 3 heart (250 cm ) at 0.2 mm spatial resolution (Hunter & Borg 2003), it would take more than six weeks to run a 1 s simulation on 64 processors (using an oversimplistic cellular model), even assuming, over-optimistically, linear scaling of the computational time with the number of mesh nodes. Furthermore, it may be the case that an even greater spatial resolution may be required at some critical locations (in areas of complex geometry or high heterogeneity). Clearly, this is not practicable for any real application, which would be likely to require simulations over much longer time scales: minutes, not seconds. There are a number of projects currently underway around the world to develop petascale-class computer systems with tens or hundreds of thousands of CPU cores. At least some of the required performance improvement to make bidomain simulation more tractable could be realized by running heart simulators on such massively parallel systems. However, the required speedup factor is so large that simply porting existing codes to even the most powerful supercomputers would not be sufficient, even assuming linear scaling. (The scaling is, in fact, almost certain to be well below linear.) So, much of the necessary speedup will need to be found from improved numerical algorithms rather than simply implementing existing methods efficiently. In this paper, we describe the current procedure for numerical modelling of the bidomain equations (2) and review recent developments in adaptive numerical algorithms, the development of spectral element (SE) methods as a high-performance alternative to FE, and state-of-the-art parallel linear solvers for large-scale algebraic systems (3). Furthermore, we discuss the problems that are certain to be encountered in designing simulators that exploit these algorithms and are intended to be run on massively parallel computers (4). Some preliminary ideas for overcoming these issues and progressing towards a petascale heart simulation code are also introduced. 2. Numerical modelling of cardiac electrophysiology The propagation of a cardiac action potential is a multi-scale phenomenon, with the electrical activity through ion channels in individual, but electrically connected, cells driving a wavefront over the whole heart. This wave can Next-generation HPC cardiac simulation propagate between neighbouring cells via electrically active gap junctions. However, the cell membrane is also electrically active as it contains ion channels and pumps that selectively allow ions to cross it. Hence, a potential gradient is also seen in the extracellular space and the cardiac action potential can propagate in this also. In order to simulate this behaviour numerically, it is necessary to derive a system of governing equations (...truncated)