The reductionist approach in biology has provided a wealth of detailed information on the genetic and molecular components from which living systems are composed. This is demonstrated by the fact that there are now several hundred biological databases accessible from the World Wide Web. The importance of this achievement cannot be understated—it has transformed the nature of both biology and medicine. There is, however, growing recognition that the tabulation of genetic and molecular building blocks from which biological systems are composed is not sufficient for understanding the functional properties of these systems. Rather, it is becoming clear that the emergent, integrative behaviors of biological systems result from complex interactions between all system components, and that knowledge of each component, however detailed, is not sufficient by itself to understand such behaviors.

Cardiac electrophysiology is a discipline with a long tradition of integrative modeling. In this paper, we describe our approach to the integrative modeling of cardiac function. This approach spans multiple levels of biological analysis ranging from subcellular to tissue. In developing these models, we have found it necessary to apply diverse analytical methods including: (1) imaging techniques for measurement of anatomic structure and biophysical and biochemical responses of cells and tissue; (2) parallel computing techniques for the numerical solution of large systems of model equations; and (3) interactive visual exploration of model dynamic behavior. We demonstrate how these various techniques contribute to the development of an integrative model of the heart.

Measurement and modeling of biophysical and biochemical properties of single cardiac myocytes

A quantitative understanding of the biophysical and biochemical properties of individual cardiac myocytes is now emerging through application of a range of experimental methodologies including recording of: (1) whole-cell membrane currents from acutely dissociated individual myocytes; (2) membrane currents generated by recombinant channels expressed in cell culture; (3) single-channel gating properties; and (4) time-varying intracellular ion concentration using fluorescent indicators. A number of computational models of single ventricular myocytes based on such data have been developed.^1-3 In the following, we describe a computational model of the canine ventricular myocyte,⁴ a species that is frequently used in the study of heart disease such as myocardial infarction and heart failure.^5,6

The canine myocyte model. The canine myocyte model (see Figure 1) describes voltage-gated membrane currents (labeled I_x in Figure 1, where x denotes the type of current), voltage- or concentration-dependent membrane transporters (labeled I_NaCa, I_p(Ca), I_NaK, and I_SR), mechanisms for sequestration and release of calcium (Ca²⁺) within the cell, and time-varying intracellular concentrations of sodium (Na⁺), Ca²⁺, and potassium (K⁺). Inward voltage-gated membrane currents (the fast inward Na⁺ current I_Na and the L-type [low-threshold] Ca²⁺ current I_Ca,L) generate the rapid upstroke of the cardiac action potential (AP) (see Figure 2A, red line). Voltage-gated outward membrane currents carried primarily by K⁺ determine the repolarizing phase of the AP.

Figure 1 Figure 2

The voltage- and time-dependent behavior is described using systems of nonlinear ordinary differential equations.⁷ In this formulation, current produced by the ensemble behavior of a particular class of ion channels is assumed to be given by the product of a voltage-dependent membrane conductance G_x(v(t), a) and an electrical driving force v(t) – E_x(t), where v(t) is time-varying membrane potential, a is a vector of state variables on which the conductance depends, and E_x(t) is the reversal potential of the permeable ionic species, and is given by the Nernst equation.⁸ As an example, the Hodgkin-Huxley model of the fast inward Na⁺ current I_Na(t) is given by the following system of equations (see Hille⁸):

I_Na(t) = G_Na[v(t), a][v(t) – E_Na(t)]
a = [m[v(t)] h[v(t)]]^T
G_Na[v(t), a] = G_Nam[v(t)]³h[v(t)]
[v(t)] = _m[v(t)]{1 – m[v(t)]} – ß_m[v(t)]m[v(t)]
[v(t)] = _h[v(t)]{1 – h[v(t)]} – ß_h[v(t)]h[v(t)]

ENa(t) =	RT	ln (Nao/Nai)	(1)
	zF

where m[v(t)] and h[v(t)] are state variables describing the activation and inactivation of the current, G_Na is the peak Na⁺ conductance, _m, ß_m, _h and ß_h are voltage-dependent rate constants that are determined experimentally, and Na_o and Na_i are external and internal Na⁺ concentration.

A diversity of voltage-dependent membrane currents are defined in the canine ventricular cell model.^4,9 These include: (1) the fast inward Na⁺ current I_Na(t) responsible for the rapid upstroke of the AP; (2) the L-type Ca²⁺ current I_Ca,L(t) that contributes to the upstroke, and helps set and maintain the AP plateau; (3) the rapidly activating delayed rectifier K⁺ current I_Kr, which contributes to repolarization of the AP; (4) the slowly activating delayed rectifier current I_Ks, which also contributes to repolarization of the AP; (5) the Ca²⁺-independent transient outward K⁺ current I_to1, which is a modulator of AP shape and duration; and (6) the instantaneous inward rectifier current I_K1, which functions in repolarization of the AP at relatively hyperpolarized membrane potentials, and which helps set the resting potential of the myocyte.

Voltage- and concentration-dependent membrane transporters, labeled I_NaCa, I_p(Ca), I_NaK, and I_SR in Figure 1, function to maintain internal ion concentrations near equilibrium values. The Na⁺-Ca²⁺ exchanger (I_NaCa) is the principal means by which Ca²⁺ that enters the cell during an AP is removed from the cell following that AP. The exchanger current depends on Na⁺ and Ca²⁺ concentration both inside and outside the cell, and on membrane potential v(t). The Na⁺-K⁺ pump functions to extrude Na⁺ from the cell that enters during an AP. The sarcolemmal Ca²⁺ pump utilizes energy in the form of adenosine triphosphate (ATP) to extrude Ca²⁺ from the cell, and the sarcoplasmic reticulum (SR) Ca²⁺-ATPase uses ATP as an energy source to pump Ca²⁺ from the cytoplasm into the SR. The current I_k^tr(t) generated by each of these transporters is assumed to be an instantaneous function of the relevant state variables, and is described using algebraic equations:

Ikstr(t) = Fk[Ci(t), Co(t), v(t)]

(2)

where C_i(t) and C_o(t) are the time-varying internal and external concentration of each ionic species. The time rate-of-change of ionic concentrations within the cell is given by:

i(t) = –	ICitot(t)	(3)
	zFVeff

where C_i(t) is the intracellular concentration of the ith ionic species, I_Ci^tot(t) is the total membrane current carried by ion species C_i, z is the valence, F is Faraday's constant, and V_eff is the cell volume.

Finally, membrane potential v(t) is itself a state variable, with time derivative proportional to the sum of ionic currents flowing across the cell membrane:

(t) = –	1		N	Ii(t)	(4)
	C
			i=1

where C is total cell membrane capacitance and I_i(t) is the ith voltage-dependent membrane or transporter current.

Cardiac Ca²⁺ dynamics and excitation-contraction coupling. Cardiac myocytes have evolved intricate mechanisms for the sequestration and release of intracellular calcium and regulation of muscle contraction—a process known as excitation-contraction (EC) coupling. EC coupling involves a close interplay between L-type Ca²⁺ channels in the sarcolemmal membrane, and Ca²⁺-induced Ca²⁺-release channels (referred to here as RyR) in the SR membrane. Modeling of this process is described in Jafri¹⁰ and may be understood with reference to Figure 1. During the initial stages of the action potential, voltage-gated L-type Ca²⁺ channels (also called di-hydropyridine receptors—DHPRs) in the sarcolemmal membrane open, allowing the entry of Ca²⁺ into a restricted subspace between these channels and the SR membrane. As subspace Ca²⁺ concentration increases, Ca²⁺ binds to the RyR, increasing their open probability and leading to Ca²⁺ release from the junctional SR (JSR). The amount of Ca²⁺ released from the JSR is significantly more than the amount of trigger Ca²⁺ entering via L-type Ca²⁺ channels. Diffusion of subspace Ca²⁺ into the cytosol leads to the Ca²⁺ transient that may be imaged using Ca²⁺-sensitive fluorescent indicators. An example of a Ca²⁺ transient measured experimentally using the indicator Indo-1 is shown in Figure 2B (red). This transient was measured simultaneously with the AP shown by the solid line in Figure 2A. Model reconstructions of an AP and Ca²⁺ transient are shown in Figures 2C and 2D. Following the AP, cytosolic Ca²⁺ is either extruded from the cell by the Na⁺-Ca²⁺ exchanger, or pumped back into the SR by the SR Ca²⁺-ATPase. The canine ventricular cell model accurately reproduces not only the shape of the intracellular Ca²⁺ transient but also the rates at which Ca²⁺ is extruded and resequestered by these pumping mechanisms.

Numerical solution of model equations. A system of 34 coupled nonlinear ordinary differential equations, with form similar to those presented previously, define the canine ventricular cell model. A complete description of the original equations defining the model is available in Winslow⁴ (Appendix) and source code is available at http://www.bme.jhu.edu/~rwinslow. Equations are integrated using the Merson modified Runge-Kutta 4th order adaptive step algorithm,¹¹ with a maximum step size of 100 microseconds (µsec) and maximum error tolerance of 10^–6. The error from all variables is normalized to ensure that each contributes equally to the calculation of global error. Initial conditions listed in Winslow⁴ (Appendix) were computed in response to a periodic pulse train of frequency 1 hertz (Hz), and were determined immediately preceding the 11th pulse. Action potentials are initiated using 100 microamperes per microfarad (µA/µF) current injection for 500 µsec. Equations describing the dynamics of the RyR have extremely rapid kinetics, and use of a stiff integrator such as DVODE¹² can improve simulation performance substantially.

Applications of the single cell model. The development and analysis of biophysically and anatomically detailed computational models of physiological systems is providing an important tool for relating changes in gene expression to changes in biological function at the levels of cell and tissue in both health and disease. This approach is illustrated using a recently developed computational model of the failing myocyte.

Heart failure is a significant disease in the United States, with more than two million patients diagnosed with this illness. Patient prognosis remains poor, with over 15 percent dying within one year of initial diagnosis, and a greater than 80 percent six-year mortality rate. Up to 50 percent of these patient deaths result from sudden cardiac death (SCD). Heart failure is in fact the leading cause of SCD in the United States. Remarkable advances have been achieved in understanding the genetic and molecular basis of heart failure. Studies in both animal models (specifically, the canine tachycardia pacing-induced model of heart failure^5,6,13) and human patients have shown significant “down regulation” of the transient outward current I_to1 and the fast inward rectifier current I_K1, respectively, in end-stage heart failure.^14-17 Studies have also shown reduced expression of SERCA2 (the gene encoding the SR Ca²⁺-ATPase) and increased expression of NCX1 (the gene encoding the Na⁺-Ca²⁺ exchanger protein).¹⁸ Equally significant advances have been achieved in characterizing electrophysiological properties of failing ventricular myocytes. It is now known that action potential duration is prolonged in failing myocytes. An example of AP prolongation is shown by the black line in Figure 2A, which was recorded in a left ventricular myocyte isolated from a failing canine heart, in which heart failure was induced using the tachycardia pacing procedure. AP duration is more than double that recorded in normal myocytes (red line, Figure 2A). AP prolongation is known to be arrhythmogenic, increasing the likelihood of severe, life-threatening arrhythmias such as early after-depolarizations (EADs).^18,19 Understanding the mechanism of AP prolongation therefore has clinical relevance to the prevention of SCD in heart failure. It is also known that Ca²⁺ transient amplitude and rate of decay are significantly reduced in heart failure. An example of a Ca²⁺ transient, measured during the AP shown by the black line in Figure 2A, is shown by the black line in Figure 2B. This example of a Ca²⁺ transient measured in a failing myocyte exhibits significant reduction in amplitude and slowed decay. Ca²⁺ transient amplitude reduction contributes to reduced contraction of the heart during heart failure.

The myocyte model has been used to investigate the relationship between altered patterns of gene expression and mechanisms of arrhythmia in heart failure.⁴ These computational models are based on data obtained from the canine tachycardia pacing-induced model of heart failure. Models were used to show that prolongation of AP duration (Figure 2C) and reduction of Ca²⁺ transient amplitude and decay rate (Figure 2D) observed in failing myocytes are accounted for primarily by altered expression of Ca²⁺ handling proteins. Unexpectedly, down regulation of outward K⁺ currents appears to play a secondary role in AP prolongation. By enabling a clear identification of the biophysical mechanisms responsible for these altered responses of failing myocytes, the models have suggested specific therapeutic targets for reducing AP prolongation and risk of arrhythmia in heart failure. Specifically, drugs that inhibit phospholamban (an inhibitor of the SR Ca²⁺-ATPase) or that act directly on the SR Ca²⁺-ATPase to increase its transport rate would be expected to improve contractile strength of the heart and reduce AP duration prolongation in heart failure.

Mapping and modeling of cardiac anatomy

Development of biophysically detailed models of single cardiac myocytes has contributed greatly to our understanding of processes underlying excitation and repolarization in the heart. However, knowledge of cellular events alone is not sufficient for understanding the complex patterns of electrical conduction within the cardiac ventricles. The reason for this is that the detailed spatial orientation of muscle fibers within the heart has a profound influence on the ways in which electrical and mechanical activity is both generated and propagated. The challenge now is to acquire a more detailed understanding of this anatomical structure—knowledge that ultimately will enable the development of biophysically and anatomically accurate models of the heart as a functioning organ.

Histological reconstruction of ventricular anatomy. Cardiac fibers are arranged as counter-wound helices encircling the ventricular cavities, and the orientation of these fibers depends on transmural location.^20-23 Fibers tend to lie in planes parallel to the epicardium, approach a longitudinal orientation on the ventricular surfaces, and rotate toward the horizontal near the midwall. Detailed measurements of fiber orientation within the cardiac ventricles have been obtained by Nielsen,²⁴ who used these data to construct a finite-element model describing myofiber structure. LeGrice et al. have made comprehensive measurements of canine ventricular macro- and micro-structure²⁵ in which they report an additional level of organization of the myofibers—they appear to be arranged in distinct myocardial laminae about 4 myocytes thick, which are separated from adjacent laminae by the extracellular collagen network. This finding is also supported by the recent local reconstruction of myocardial lamina in rabbit heart reported by Costa.²⁶

Estimation of ventricular fiber orientation using diffusion-tensor magnetic resonance imaging. The direct anatomical reconstructions just described are extremely labor-intensive and time-consuming. More rapid techniques for achieving ventricular reconstructions would greatly enhance our ability to characterize cardiac structure in both health and disease. Recently, a number of studies have suggested that diffusion-tensor magnetic resonance imaging (DTMRI) may be used to determine muscle fiber orientation of tissues. DTMRI yields estimates of a voxel-averaged diffusion tensor, the eigenvectors and eigenvalues of which specify the principal directions and rates of water diffusion at each voxel of the tissue image. The eigenvector corresponding to the maximum eigenvalue of the diffusion tensor points in the direction of maximum rate of diffusion. Since muscle fibers are long and thin, this direction has been hypothesized to correspond to the orientation of the long axis of a muscle fiber.

Quantitative histologic verification of this hypothesis has been lacking until recently. Hsu²⁷ provided the first quantitative correlation of DTMRI-based and histological estimates of fiber orientation in an excised portion of the right ventricle. Scollan²⁸ showed that DTMRI estimates of fiber orientation in local tissue regions obtained from perfused, nonbeating rabbit heart accurately reproduce histologic measurements of orientation made at the same locations in the same heart. It was demonstrated subsequently that maps of fiber orientation with much improved spatial resolution can be obtained in fixed myocardium—a preparation that minimizes possible artifactual motion of the heart.²⁹

Reconstruction of ventricular anatomy using DTMRI. DTMRI has recently been used to perform the high spatial resolution reconstruction of the entire cardiac ventricles in rabbit heart.³⁰ The following imaging protocol was used in these studies. First, three-dimensional GRASS (gradient recalled acquisition in the steady-state) intensity images were collected for subsequent use in defining the epicardial and endocardial surfaces. This yielded a set of 128 short-axis slices with in-slice spatial resolution of 156 × 312 micrometers (µm), and a slice separation of 469 µm. Following acquisition of the GRASS intensity data set, DTMRI was used to estimate myocardial fiber orientation. These images were obtained using a slice-selective fast spin-echo diffusion-weighted technique followed by a two-dimensional discrete Fourier transform. Slice thickness was varied depending on longitudinal position. In regions at and near the base of the heart where epicardial surface curvature is smallest, slice thickness was set to 2 millimeters (mm). At more apical regions where curvature is larger, slice thickness was set to 1 mm. Eight basal short-axis sections were obtained with a slice thickness of 2 mm, followed by twelve apical short-axis sections at a slice thickness of 1 mm. This protocol was selected as a trade-off between total imaging time (approximately 20 hours) and spatial resolution, and yields one 256 × 128 × 128 matrix of intensity data, which is used to estimate ventricular geometry, and a coincident 256 × 128 × 20 matrix of points at which estimates of the 3 × 3 diffusion tensor are available.

Reconstruction of ventricular anatomy from these data involves four steps. First, epicardial and endocardial boundaries are estimated using both the short-axis GRASS intensity and diffusion-weighted short-axis slices. This is done using a semiautomated image segmentation tool called HeartWorks.³¹ This yields a series of 128 planar contours defining epicardial and endocardial boundaries. Second, these contours are used to reconstruct the ventricular epicardial and endocardial surfaces for visualization of ventricular structure. Third, diffusion tensor eigenvalues and eigenvectors are computed from the diffusion-weighted short-axis sections that are coincident with a subset of the GRASS imaging short-axis sections, thus providing estimates of fiber orientation within the myocardium in each of these sections; and fourth, fiber orientation is interpolated for myocardium in those short-axis sections for which there are no corresponding diffusion data.

The semiautomated image segmentation algorithms in HeartWorks provide a set of 3 × 128 planar contours that specify epicardial and endocardial surfaces, where each contour is represented by a set of ordered points. Given these contours, ventricular surfaces are reconstructed using a piecewise smooth-surface reconstruction algorithm, developed by Hoppe and DeRose, which determines an optimal triangular tiling between points in adjacent contours.³²

Figure 3A shows reconstruction of a rabbit ventricle based on GRASS image data obtained at a spatial resolution of 156 × 312 × 469 µm. Image segmentation was performed using the method of active contours. Papillary muscles projecting into the ventricles were edited from the images using HeartWorks. Left ventricular (LV) and right ventricular (RV) surfaces are shown in gold and red, respectively. The epicardial surface is rendered as a wire-frame mesh. The high spatial resolution afforded by GRASS imaging is sufficient to reveal the highly detailed structure of both the epicardial and endocardial surfaces. Figure 3B–3E shows the reconstructed ventricle of Figure 3A into which 12 of 20 short-axis DTMRI sections are inserted. Each of these sections shows fiber-inclination angle coded according to the indicated color map. These sections show a transmural variation in fiber angle from the epicardial to endocardial surfaces.

Figure 3

Bringing maps and models together: Simulation of electrical excitation in the heart

The development of rapid methods for the anatomical reconstruction of cardiac ventricular geometry and fiber organization described in the previous section now creates the possibility that anatomical and computational models of an extensive library of normal and diseased hearts may be generated, and that structural and electrical properties of these hearts may be computed and compared. The following section describes methods for computational modeling of the ventricular myocardium.

Governing partial differential equations. The bidomain equations describe the flow of electrical current within the myocardium between the intracellular and extracellular domains.^33,34 This approach treats each domain of the myocardial tissue as a continuum, rather than as being composed of discrete cells connected by gap junctions and surrounded by the extracellular milieu. Thus, quantities such as conductivity and transmembrane voltage represent spatial averages.

The bidomain equations are derived by applying conservation of current between the intra- and extracellular domains. The equations consist of coupled parabolic and elliptic equations that must be satisfied within the myocardium, and an additional elliptic equation that must be satisfied in the bath, or tissue, surrounding the heart. Solution of the full bidomain equations is known to be important when modeling the effects of electric-field stimulation on cardiac electrical responses. However, in circumstances not involving field stimulation, the bidomain equations are usually simplified considerably.³⁵ If the tissue surrounding the body surface is taken to be a good insulator, and the assumption of equal anisotropy is assumed to hold, namely that

Mi(x) =	1	Me(x)	(6)

where is called the anisotropy ratio, then the bidomain equations reduce to the (monodomain) parabolic reaction-diffusion equation

v	(x, t) = –	1		Iion(x, t) + Iapp(x, t) –	1				• (Mi(x)v(x, t))		(7)
t		Cm			ß		+ 1

where x is spatial position within the myocardium, v(x, t) is the transmembrane voltage, C_m is the membrane capacitance per unit area, I_ion(x, t) is the sum of the ionic currents per unit area through the membrane (positive outward), as given by ionic models of the myocyte as described earlier, I_app(x, t) is an applied stimulus current per unit area, ß is the ratio of membrane area to tissue volume, and M_e(x) and M_i(x) are the extracellular and intracellular 3 × 3 conductivity tensors at each point x. Equation 7 is solved subject to a no-flux boundary condition at the cardiac surfaces.

The intracellular conductivity tensors at each point within the heart are specified by fiber orientation and by specific intracellular conductivities in each of the local coordinate directions. The intracellular conductivity tensor in the local coordinate system, G_i(x), is defined as

Gi(x) =

1,i

2,i

3,i

(8)

where _1,i is the longitudinal and _2,i and _3,i are the transverse intracellular conductivities, respectively. This local tensor may be expressed in global coordinates to give the intracellular conductivity tensor of Equation 6 using the transformation

Mi(x) = P(x)G(x)PT(x),

(9)

where P(x) is the transformation matrix from local to global coordinates at each point x. When working from DTMRI data, the columns of P(x) are set equal to the eigenvectors of the diffusion tensor estimated at point x.

Parallel solution of the monodomain equation. Solution of the monodomain equation (Equation 7) throughout the ventricular myocardium is computationally demanding and is typically performed using parallel computation. In this section, we describe a parallelization method implemented on an IBM RISC System/6000* (RS/6000*) SP* POWER3 symmetric multiprocessor (SMP) system running IBM Parallel Environment (PE) for AIX* v2.4 and XL FORTRAN v6.1. Each of the four nodes on this system consists of eight 222 megahertz processors and one gigabyte of shared memory. Each processor has 32 kilobytes of primary instruction cache, 64 kilobytes of primary data cache, and 4 megabytes of direct-mapped secondary cache. The implementation achieves multilevel parallelism by utilizing the Message Passing Interface (MPI) for communication among nodes, and OpenMP directives for parallel computation within each node. The implementation is highly portable, running on either SMP architectures, distributed networks of scalar processors, or hybrid architectures.

The first step of this parallel solution method is to find discrete spatial terms for Equation 7. The divergence operator in the diffusion term generates first, second, and mixed second-order spatial derivatives as shown in the following expansion:

(10)

First derivatives (advective term) are approximated using upwind differencing:

V		V(x + 1, y, z) – V(x, y, z)	(11)
x		hx

The second-order derivatives are approximated by finite differencing using a three-point stencil:

2V		V(x + 1, y, z) – 2V(x, y, z) + V(x – 1, y, z)	(12)
x2		hx2

Mixed second-order derivatives are approximated using a five-point stencil:

2V		V(x+1, y+1, z) – V(x – 1, y+1, z) – V(x+1, y – 1, z) + V(x – 1, y – 1, z)
xy		2hxhy

	=	V(x+1, y+1, z) – V(x, y, z) + V(x, y, z) – V(x – 1, y+1, z)
		2hxhy
	–	V(x+1, y–1, z)+V(x, y, z) – V(x, y, z) + V(x – 1, y – 1, z)
		2hxhy

	=	[V(x+1, y+1, z) – V(x, y, z)] + [V(x, y, z) – V(x – 1, y+1, z)]
		2hxhy
	+	[V(x, y, z) – V(x+1, y – 1, z)] + [V(x – 1, y – 1, z) – V(x, y, z)]	(13)
		2hxhy

The discrete terms are achieved using a regular three-dimensional mesh as shown in Figure 4A. The spatial resolution of the grid affects the numerical accuracy of simulation results. To determine an appropriate grid resolution, a tissue block of one cubic centimeter was constructed for a convergence test. One face of the block was stimulated and parameters such as activation time (the interval between the time at which the stimulus is applied and the time at which the action potential exhibits the greatest rate of increase), conduction velocity (distance traveled by the wave per unit time), peak voltage (maximum value of membrane potential), and upstroke velocity (maximum time rate-of-change of membrane potential) were measured at the same points in space as the grid was successively refined. Figure 5, plotting activation time as a function of grid resolution, shows that there is a large error in activation times when grid resolution is larger than 500 µm, and that a minimum spatial resolution of 250–300 µm is required to model propagation accurately. On the basis of these and similar data, the dimensions of the mesh were chosen to be 101 × 167 × 202 points in the x, y, and z coordinate directions (see Figure 4A) for the activation studies reported here. Corresponding mesh resolution was 235 × 240 × 198 µm.

Figure 4 Figure 5

Of the more than three million lattice points, 906183 were within the myocardium (as opposed to the fluid-filled ventricles). Thus, numerical simulation of electrical activation in the heart model required the solution of more than 30 million coupled nonlinear ordinary differential equations subject to an initial condition (that the heart is in an electrically quiescent state). A number of numerical methods have been used to perform the time integration of this system of equations. These include: (1) the Euler method^36-38; (2) fourth-order nonadaptive and Merson modified adaptive Runge-Kutta methods (RK4)^11,36-38; and (3) an operator-splitting method based on the reacting flow model developed by Knio.^39,40 The latter method uses a backward differentiation algorithm implemented in DVODE (double-precision variable coefficient ordinary differential equation solver) to integrate the reaction term. In the simulations reported here, an RK4 fixed-step method with 2.5 µsec time step is used.

Parallel implementation involves partitioning the computational grid and associated data into a set of N subdomains, and then distributing these subdomains to the SMP nodes that communicate with one another using MPI. To do this, as described earlier, the ventricles are “immersed” in a regular three-dimensional lattice (Figure 4A). The number of subdomains N is set equal to the number of processor nodes. In FORTRAN, data are arranged in the main memory of each node according to the leading dimension of the underlying data arrays. In Figure 4, this is taken to be the z, y, and x dimensions. The myocardium is partitioned into subdomains along the x-axis (long axis) so that data within a subdomain will be localized in shared memory. This in turn optimizes memory and cache accesses. Each subdomain therefore consists of a set of short-axis slices in the y-z plane. Since roughly 40–50 percent of grid points in each short-axis section are within the myocardium, these points (referred to as “active points”) are identified and totaled by use of a mask. To optimize load balancing, the algorithm generates subdomains consisting of an integer number of short-axis sections.

After performing domain decomposition, each node loads and initializes only the parameters and data that pertain to its subdomain. Data required to compute the diffusion term in Equations 10–13 and to compute the time evolution of membrane potential v(t) require nearest-neighbor communication. At most grid points, these data are available in the node's shared memory. For grid points on the boundary of a subdomain, membrane potentials of the neighbors are retrieved from adjacent subdomains before each iteration of the computation. This is accomplished by padding the x-dimension of the subdomain with additional “ghost” cells. Values for these cells are computed by the neighboring subdomains and exchanged between nodes via MPI communication calls at each time step. Figure 4B illustrates a volume partitioned using this algorithm. The x-z plane is shown on the right, with circles representing cells updated by processors local to the subdomain and diamonds representing ghost cells updated by processors in neighboring nodes. All other state variables defining properties of membrane currents (Equation 1), membrane transporter currents (Equation 2), and intracellular concentrations (Equation 3) are local variables in the sense that their time evolution at each grid point may be computed given knowledge of membrane potential at that grid point.

Once all required data are in the memory of the SMP node, all processors within a domain begin computation on a specific short-axis section in the subdomain. Upon finishing, each processor proceeds to the next available slice in its subdomain. The distribution of workload between the two processors is guided by OpenMP directives used to make parallel the loops that compute the currents and integrate the state variables. The scheduling of the parallel loop is static, as the overhead of dynamic scheduling is justified only when there is a high degree of load imbalance.

Whole-heart simulation results. Figure 6 shows snapshots of the electrical activation sequence simulated for a rabbit heart reconstructed using DTMRI as described in the previous section. Intracellular conductivities (Equation 8) were assumed to be transversely isotropic (s_2i = s_3,i). The longitudinal value was set to 2.49 milliseconds/centimeter (ms/cm), and the transverse values to 0.28 ms/cm (8.89:1 ratio), yielding longitudinal and transverse conduction velocities (CVs) of 76 centimeters/second (cm/s) and 26 cm/s, respectively. These conduction velocities are similar to those measured in rabbit myocardium. A current stimulus was injected into a set of cells on the endocardial surface in order to initiate electrical depolarization of the heart. In the first panel, the heart is at its resting potential of –80 millivolts (mV), denoted by the blue color. Early activation is pictured in the second panel, where the endocardial surfaces of the septum have taken on an orange color, corresponding to more positive membrane potentials. The simulated value of total activation time was 30.72 milliseconds—very similar to values reported experimentally. The average CV was 50.54 ± 18.72 cm/s. CV reached a maximum value of 114.31 cm/s at the epicardium due to the imposition of the no-flux boundary condition, and thus reduced electrical loading of cells on the surfaces of the heart.

Figure 6

The three-dimensional heart model has also been used to test the hypothesis that a particular class of cellular arrhythmia known as early after-depolarizations (EADs) may trigger re-entrant arrhythmias in the failing heart. To do this, cell properties within one region of the right ventricle, representing about 5 percent of the total myocardial mass, were adjusted to produce the extreme prolongation of action potential duration (APD) observed during heart failure by reducing I_to1, I_K1, and the SR Ca²⁺-ATPase magnitude, and increasing the Na⁺-Ca²⁺ exchanger magnitude. These changes, through their effects on APD, produced EADs in cells within that region. A wave of excitation was initiated in the model using a single excitatory stimulus applied at endocardial stimulation points. EADs triggered within the spatially localized region of the model heart were seen to initiate re-entrant waves that were completely self-sustaining for the full five seconds of the computation. Images of electrical activity on the model ventricular surface are seen in Figure 7.

Figure 7

Interactive visualization of simulation data

As can be appreciated from Figures 6 and 7, the generation and propagation of electrical activity within the heart is governed by large systems of ordinary differential equations defined over a geometrically complex, nonhomogeneous anatomical structure. Understanding the time evolution of these equations within such a complex geometry is a challenging task. To aid in this understanding, we have developed an interactive statistical visualization tool, called WEAVE (Workbench Environment for Analysis and Visual Exploration), and are applying this tool to analysis of cardiac simulation data.⁴¹

WEAVE supports the rapid development of multiple, customized statistical analyses of the simulation data that are linked to the underlying geometry of the simulation data set. Statistical analyses include, but are not limited to, generation of histograms, pie charts, scatterplots, calculation of covariation between multiple data sets, and clustering. An example of the use of WEAVE is shown in Figure 8. The upper-left panel displays the heart geometry, on which is superimposed by means of a color code the spatial values of a selected state variable or model membrane current at a particular instant of time following electrical activation of the model heart by a current stimulus applied to a local region of the endocardial surface. Any state variable or current may be selected for display—in this instance we show values of the transient outward current I_to1 at every grid point in the model.

Figure 8

The I_to1 current is an outward potassium current that activates following membrane depolarization, but then undergoes rapid inactivation. The lower-left panel displays a histogram of I_to1 values. Inspection of this histogram shows that I_to1 current magnitude is near zero for the majority of grid points in the model, as the majority of cells within the model are in a resting state. The histogram exhibits a second, much smaller peak at an ordinate value of approximately 2.3 picoamperes per picofarad (pA/pF). As shown in this panel, the mouse cursor was used to select three different sets of points in this histogram: (1) a set of points about this second peak (green); (2) a set of points about the major peak positioned near 0 pA/pF (red); and (3) an intermediate set of points between these two peaks (blue). On selection, the location of each set of points is displayed on the heart geometry shown in the upper left panel using the chosen color code. The third panel in the upper-right corner of the display is a scatterplot of I_to1 density (in pA/pF for the ordinate) and membrane potential (in mV for the abscissa). Points in this display are also shown using the color code selected during interaction with the histogram display. These data show that I_to1 current is small for cells with low membrane potential (red), because these cells are at rest. Cells for which voltage is high and I_to1 density takes on intermediate values (blue) define the activation wavefront (depolarization of these cells has just occurred, and I_to1 is beginning to activate). Cells with high voltage and high current density are cells that have just activated.

These data illustrate the ability of WEAVE to perform a wide range of statistical analyses on simulation data, and to enable the user to explore properties of subsets of simulation data with reference to the underlying anatomical model. Such interactive statistical analyses are crucial for interpretation and for understanding of the very large simulation data sets that arise in integrative modeling of the heart.

Discussion

This paper has reviewed modeling research in three broad areas: (1) models of single ventricular myocytes; (2) methods for the reconstruction and modeling of ventricular geometry and microanatomy; and (3) integrative modeling of the cardiac ventricles. We have seen that the level of biophysical detail, and hence the accuracy and predictability of current ventricular myocyte models, is considerable. Nonetheless, much remains to be done. In particular, as we acquire increasingly detailed data regarding the kinetic behavior of specific voltage-gated ion currents and transporter systems in the myocyte, and as we begin to identify the genes and gene systems encoding particular ion currents and transporters, we will need to increase the fidelity with which these mechanisms are represented in the single cell models.

In real cardiac myocytes, there exist a diversity of mechanisms that act to modulate cellular excitability. This includes - and ß-adrenergic signaling pathways acting through G protein⁴²-coupled membrane receptors to modulate properties of the L-type Ca²⁺ channel, various K⁺ channels, and Ca²⁺ transporters such as the SR Ca²⁺-ATPase. Transport of Ca²⁺ into the SR and extrusion of Na⁺ from the cell following the AP are just two examples of vitally important processes utilizing substantial energy resources of the cell in the form of ATP. The addition of these modulatory and energy-requiring homeostatic mechanisms to the cell models remains an important goal for the future.

Diffusion tensor magnetic resonance imaging now offers a relatively rapid way to measure ventricular fiber structure at high spatial resolution. The ability to rapidly acquire fiber orientation data throughout the ventricles in large populations of normal and diseased hearts will enable quantitative statistical comparison of normal and abnormal cardiac structure, and will provide insights into the possible structural basis of arrhythmia in heart disease. Availability of high spatial resolution models of the ventricles will also facilitate three-dimensional computer modeling of ventricular activation and its relationship to the underlying fiber structure of the heart in both health and disease. Unfortunately, a detailed understanding of the spatial heterogeneities within the heart, such as variation of intercellular coupling, regional expression of ionic currents, and Ca²⁺ handling proteins, is still unavailable, although significant progress has certainly been made. Clearly, the development of an anatomically and biophysically accurate model of the cardiac ventricles still remains a challenge for the future.

Acknowledgments

This work is supported by the National Institutes of Health grants RO1HL60133-01 and P50HL52307, the Whitaker Foundation, the Falk Foundation, IBM Corporation, and Physiome Sciences, Inc. D. Scollan was supported by a Medical Scientist Training Program fellowship. We thank our colleagues in experimentation Eduardo Marban, Brian O'Rourke, Gordon Tomaselli, and Brad Nuss for their collaboration on the heart failure modeling project.

*Trademark or registered trademark of International Business Machines Corporation.

Cited references and note

Accepted for publication February 3, 2001.