Numerical solution of the bidomain equations

BY S. LINGE 1 2 J. SUNDNES () 0 2 M. HANSLIEN 2 G. T. LINES 0 2 A. TVEITO 0 2 0 Department of Informatics, University of Oslo , PO Box 1080, Blindern , Norway 1 Telemark University College , PO Box 203, 3901, Porgsgrunn , Norway 2 Simula Research Laboratory , PO Box 134, 1325 Lysaker , Norway Knowledge of cardiac electrophysiology is efficiently formulated in terms of mathematical models. However, most of these models are very complex and thus defeat direct mathematical reasoning founded on classical and analytical considerations. This is particularly so for the celebrated bidomain model that was developed almost 40 years ago for the concurrent analysis of extra- and intracellular electrical activity. Numerical simulations based on this model represent an indispensable tool for studying electrophysiology. However, complex mathematical models, steep gradients in the solutions and complicated geometries lead to extremely challenging computational problems. The greatest achievement in scientific computing over the past 50 years has been to enable the solving of linear systems of algebraic equations that arise from discretizations of partial differential equations in an optimal manner, i.e. such that the central processing unit (CPU) effort increases linearly with the number of computational nodes. Over the past decade, such optimal methods have been introduced in the simulation of electrophysiology. This development, together with the development of affordable parallel computers, has enabled the solution of the bidomain model combined with accurate cellular models, on geometries resembling a human heart. However, in spite of recent progress, the full potential of modern computational methods has yet to be exploited for the solution of the bidomain model. This paper reviews the development of numerical methods for solving the bidomain model. However, the field is huge and we thus restrict our focus to developments that have been made since the year 2000. The bidomain equations are considered to represent cardiac electrophysiology accurately and have been applied successfully for a wide range of studies. A particularly large volume of research has been carried out on the study of 1. Introduction re-entrant arrhythmias and their treatment by externally applied electrical stimulus (e.g. Trayanova et al. 2006). In spite of being applied successfully for a number of valuable studies, current research based on the bidomain model is slowed by a number of substantial challenges. Some of these challenges are directly related to the huge complexity and variability of living systems, and the consequent difficulty of developing mathematical models and acquiring the necessary experimental data. However, there are also severe challenges related to solving the equations efficiently. Although great advances have been made in both computer hardware and computational algorithms, computational time and, to some extent, memory usage remain a severe obstacle for biomedical research that is based on the bidomain model. The purpose of this review is to give an overview of recent advances and the current state of research in this area, and to indicate directions and requirements for future research. Owing to the large volume of research in this field, we will limit our discussion to developments after the year 2000. (a ) The mathematical model V$MiVv C V$MiVu Z C Iions; v; Here, equation (1.1) is a system of ODEs that describes electrochemical reactions in the cells, while equations (1.2) and (1.3) describe propagation of the electrical signal through the cardiac tissue. The boundary conditions in equations (1.4) and (1.5) state that the normal components of intra- and extracellular currents are zero on the boundary, which is a standard assumption when considering an insulated heart or tissue sample. The model may easily be extended to include heart tissue surrounded by a conductive bath or a conductive body (e.g. Pullan et al. 2005; Sundnes et al. 2006), but, for the present study, it is sufficient to consider the insulated heart. The nonlinear term Iion(v, s) describes the ionic current across the cell membrane, while the symmetric tensors Mi and Me represent tissue conductivities of the intra- and extracellular spaces, respectively. In the present formulation, the ionic current has been scaled with the cell membrane capacitance, and the conductivities have been scaled with both the membrane capacitance and the cell membrane area to cell volume ratio. See, for instance, Sundnes et al. (2006) for details. In addition to the complexity of the involved equation systems, one may argue that the main source of the computational load in solving the bidomain model is the rapid dynamics of the cellular reactions described by equation (1.1) and the ionic current term in equation (1.2). The fast reactions lead to a narrow wavefront of activation that propagates through the cardiac tissue. Capturing these fast variations in space and time requires high resolution in the discretization, which results in large systems of equations that must be solved for each time step. From a numerical point of view, there are three main approaches to reducing the computational burden: (i) using adaptive and/or high-order discretization techniques in space and time, to reduce the number of unknowns in the resulting discrete systems, (ii) applying state-of-the-art solution methods (i.e. multi-level solvers) to solve the resulting systems in an optimal manner, and (iii) adapting the solvers for parallel hardware to take advantage of the increasingly available supercomputers and computing clusters. As will be demonstrated below, all of these approaches have been investigated to some extent for the case of the bidomain model. (b ) Focus of the present paper Although, to date, the most successful approach for reducing the computational load of bidomain simulations has been to simplify or modify the model in some way, we limit our scope to computational methods for the full model. This choice is justified by the facts that: (i) the full model is superior for a number of important applications, in particular for studying the effects of external stimulus and defibrillation, and (ii) the computational techniques derived for the full model are also relevant for simplified models. The latter is particularly true for the monodomain model and various versions of the cable equation (Hodgkin & Rushton 1946; Roth & Krassowska 1998), and to a lesser extent for models such as the eikonal formulation (e.g. Keener 1991). For research contributions prior to 2000, the reader may wish to consult the review by Henriquez (1993), the more recent review by Lines et al. (2003a) and the papers by Panfilov & Holden (1997) and Zipes & Jalife (2000). In addition, the book by Keener & Sneyd (1998) is highly recommended as a basis for understanding the material presented herein, and a general introduction to numerical methods for the bidomain mo (...truncated)