372 IEEE TRANSACTIONS ON BIOMEDICAL ENGINEERING, VOL. 45, NO. 3, MARCH 1998 Efficient Integration of a Realistic Two-Dimensional Cardiac Tissue Model by Domain Decomposition Weilun Quan,* Member, IEEE, Steven J. Evans, and Harold M. Hastings, Member, IEEE Abstract—The size of realistic cardiac tissue models has been limited by their high computational demands. In particular, the Luo–Rudy phase II membrane model, used to simulate a thin sheet of ventricular tissue with arrays of coupled ventricular myocytes, is usually limited to 100 100 arrays. We introduce a new numerical method based on domain decomposition and a priority queue integration scheme which reduces the compu- tational cost by a factor of 3–17. In the standard algorithm all the nodes advance with the same time step , whose size is limited by the time scale of activation. However, at any given time, many regions may be inactive and do not require the same small and consequent extensive computations. Hence, adjusting locally is a key factor in improving computational efficiency, since most of the computing time is spent calculating ionic currents. This paper proposes an efficient adaptive numerical scheme for integrating a two-dimensional (2-D) propagation model, by incorporating local adjustments of . In this method, alternating direction Cooley–Dodge and Rush–Larsen methods were used for numerical integration. Between consecutive integrations over the whole domain using an implicit method, the model was spatially decomposed into many subdomains, and adjusted locally. The Euler method was used for numerical integration in the subdomains. Local boundary values were determined from the boundary mesh elements of the neighboring subdomains using linear interpolation. Because was defined locally, a priority queue was used to store and order next update times for each subdomain. The subdomain with the earliest update time was given the highest priority and advanced first. This new method yielded stable solutions with relative errors less than 1% and reduced computation time by a factor of 3–17 and will allow much larger (e.g., 500 500) models based on realistic membrane kinetics and realistic dimensions to simulate reentry, triggered activity, and their interactions. Index Terms—Cardiac electrophysiology, computer simulation, finite difference method, realistic membrane model. I. INTRODUCTION R EALISTIC and complex cardiac membrane models have been developed quickly in recent years using data ac- cumulated from single-cell and single-channel recordings. With new membrane models [1]–[5], computer simulations of more complicated electrophysiological phenomena, such as early and delayed afterdepolarizations, effects of antiarrythmic drugs, as well as dynamically tracking ionic environments, Manuscript received September 1, 1995; revised November 7, 1997. This work was supported in part by an award from the Whitaker Foundation and a faculty grant from the Long Island Jewish Medical Center. Asterisk indicates corresponding author. *W. Quan is with the Cardiology Department, Winthrop-University Hospi- tal, 222 Station Plaza North, Room 607, Mineola, NY 11501 USA (e-mail: wquan@winthrop.org). S. J. Evans is with the Cardiology Department, Long Island Jewish Medical Center, New Hyde Park, NY 11042 USA. H. M. Hastings is with the Mathematics Department, Hofstra University, Hempstead, NY 11549 USA. Publisher Item Identifier S 0018-9294(98)01612-7. ischemia, and other pathophysiologic abnormalities are possi- ble. Propagation models based on newly developed membrane models should yield more realistic and quantitative results. These new membrane models require highly intensive compu- tation, however, which limits their implementation into the propagation models. The numerical solution of a realistic two-dimensional (2-D) model poses a difficult problem for nu- merical analysis and, as yet, no completely reliable algorithm exists which is computationally efficient. Barr and Plonsey [6], and Henriquez [7], as well as Pollard et al. [29], used an explicit scheme of dynamically tracking the active region to reduce the total number of nodes in calcu- lations at any time step by assuming that the regions outside activation are either at rest or at a constant potential. This method offers a good estimation of the maximum rate of rise of the upstroke, , and the conduction velocity, but suffers from distortion in the foot of the action potential. In addition, confidence in the accuracy of computed action potentials for a long time periods is not warranted for an explicit numerical scheme due to the possible accumulation of round-off errors. Furthermore, dynamic tracking of active regions, based on a predetermined “threshold for propagation,” is extremely difficult to perform when the resting potential is dynamically changing during reentrant propagation [24], due to strong head-tail interactions or other physiological/pathological con- ditions (e.g., dynamically varying ionic concentrations causing varying resting potentials). Lesh et al. [8] used a numerical scheme which uncouples local membrane integration of ordinary differential equations (ODE’s) from the global step-size control for partial dif- ferential equation (PDE’s). The contribution of total axial current for each membrane element was obtained from the previous global time step of the PDE. It is not known how much error is introduced by using values of the surrounding elements at previous global steps for local updating, espe- cially when the global time step is relatively large compared with the local time steps of membrane potential integration. Leon and Roberge [9], [10] used a parallel cable network to simulate a thin sheet of cardiac tissue, offering a large amount of computational saving in updating the cable equation when lateral connections are spaced far apart. However, in a propagation model, most of the computation time is spent in calculating ionic currents. Computational saving in solving the cable equation will not decrease the overall computational requirement to a great extent. The key to computational efficiency for a propagation model is a good numerical scheme for updating ionic currents. For reconstruction of nonpropagated membrane potentials, 0018–9294/98$10.00  1998 IEEE