IEEE TRANSACTIONS ON BIOMEDICAL ENGINEERING, VOL. 51, NO. 8, AUGUST 2004 1319 Simulation of QRST Integral Maps With a Membrane-Based Computer Heart Model Employing Parallel Processing Marie-Claude Trudel, Bruno Dubé, Mark Potse*, Ramesh M. Gulrajani, Fellow, IEEE, and L. Joshua Leon, Member, IEEE In memory of Dr. Ramesh M. Gulrajani (1944–2004). The authors would like to dedicate this paper to their late colleague. Abstract—The simulation of the propagation of electrical activity in a membrane-based realistic-geometry computer model of the ventricles of the human heart, using the governing monodomain reaction-diffusion equation, is described. Each model point is represented by the phase 1 Luo–Rudy membrane model, modified to represent human action potentials. A separate longer duration action potential was used for the M cells found in the ventricular midwall. Cardiac fiber rotation across the ventricular wall was implemented via an analytic equation, resulting in a spatially varying anisotropic conductivity tensor and, consequently, anisotropic propagation. Since the model comprises approximately 12.5 million points, parallel processing on a multiprocessor computer was used to cut down on simulation time. The simulation of normal activation as well as that of ectopic beats is described. The hypothesis that in situ electrotonic coupling in the myocardium can diminish the gradients of action-potential duration across the ventricular wall was also verified in the model simulations. Finally, the sensitivity of QRST integral maps to local alterations in action-potential duration was investigated. Index Terms—Action-potential duration, cardiac membrane model, computer heart model, parallel processing, QRST integral maps. I. INTRODUCTION C OMPUTER models have long been used for the simula- tion of electrical activity in the heart and the subsequent computation of the electrocardiogram (ECG). Many of these models were of the “cellular automaton” type, with a rule-based Manuscript received July 8, 2002; revised November 11, 2003. This work was supported by the Natural Sciences and Engineering Research Council of Canada. The work of M.-C. Trudel was supported by la Fondation J.A. de Sève and by FCAR, Québec. The work of M. Potse was supported by FCAR, Québec, and by the Netherlands Organization for Scientific Research (NWO). Asterisk indicates corresponding author. M.-C. Trudel was with the Institute of Biomedical Engineering, Université de Montréal, Montréal, QC H3C 3J7, Canada. She is now with EBS Inc., Montréal, QC H2X 3V8, Canada. B. Dubé is with the Institute of Biomedical Engineering, Université de Mon- tréal, Montreal, QC H3C 3J7, Canada (e-mail: gulrajan@igb.umontreal.ca). *M. Potse is with the Institute of Biomedical Engineering, Université de Mon- tréal, Montreal, QC H3C 3J7, Canada (e-mail: mark@potse.nl). R. M. Gulrajani, deceased, was with the Institute of Biomedical Engineering, Université de Montréal, Montreal, QC H3C 3J7, Canada. L. J. Leon was with the Institute of Biomedical Engineering, Université de Montréal, Montréal, QC H3C 3J7, Canada. He is now with the Department of Electrical and Computer Engineering, University of Calgary, Calgary, AB T2N 1N4, Canada. Digital Object Identifier 10.1109/TBME.2004.827934 excitation algorithm based on the known propagation velocities of the action potential along and across the cardiac fibers. An example is the model of Lorange and Gulrajani [1]. Once the excitation reached the individual model cells, a predetermined action-potential waveform was started at each cell. A notable exception to this cellular automaton approach was the model of Leon and Horᡠcek [2] that used subthreshold electrotonic conduction to bring the cells to a threshold membrane potential for excitation, with a predetermined action-potential waveform started at this excitation potential as before. More recently, re- action-diffusion equations have been used to characterize both subthreshold and suprathreshold activity, initially with a simple FitzHugh–Nagumo-type membrane equation characterizing the model cells [3]–[5] and later with a modified Beeler–Reuter membrane equation to represent the model cells by Huiskamp [6]. Here, we present a heart model, similar to Huiskamp’s, that also uses the reaction-diffusion equation but with the modified Beeler–Reuter membrane equation replaced by a modified ver- sion of the more recent Luo–Rudy phase 1 membrane equa- tion [7]. The anatomy of our heart model is the same as that of the earlier one by Lorange and Gulrajani [1]. To enable the numerical integration of the reaction-diffusion equation, how- ever, it was necessary to augment the spatial resolution of the Lorange–Gulrajani model from 1 to 0.25 mm by the addition of intermediate points. Like Huiskamp, we too only simulated excitation of the ventricles, using starting points and starting times at the endocardium taken from the Lorange–Gulrajani model. Since the spatial resolution is much finer than that of Huiskamp’s model (0.25 as opposed to 0.6 mm), simulations with the much larger number of points (approximately 12.5 mil- lion as opposed to 800 000 in Huiskamp’s model) necessitated parallel processing on a multiprocessor Silicon Graphics Origin 2000 computer. Also incorporated into the heart model was a region of “M-cells” in the ventricular midwall, as has been documented in both dogs [8] and in humans [9]. These M-cells have been measured to have an action-potential duration that is longer than those of endocardial or epicardial cells. This sets up a transmural gradient of membrane potential during repolariza- tion that is largely responsible for the T-wave in the ECG being of the same polarity as the QRS complex. Consequently, the integrated area under each ECG lead during the QRS and T complexes is nonzero. A plot of this integrated area versus the two-dimensional (2-D) spatial location of the ECG lead on the 0018-9294/04$20.00 © 2004 IEEE