Models of cardiac cell

Post-publication activity

In this article, we briefly describe how transmembrane currents and other cellular ionic processes are modeled to produce action potentials. Then we discuss 45 models of cardiac cells and summarize their properties. For many of these models, we also include action potential plots, movie of reentrant spiral waves, and Java applets in which users can change number of stimuli and pacing rates and observe the time evolution of the voltage, currents, and variables. Although we do not present a complete list of all existing cardiac cell models, we include a broad range of models that represent different species, cell types, and model formulations from simple to biophysically detailed. Both historical and recent models are included, with many recent models incorporating more complex dynamics, such as detailed intracellular Ca²⁺ handling and Markov models of ion channels. A number of the models described have been modified in minor ways, such as modification of parameter values, to represent particular dynamics or to represent pathological states. Such minor modifications in general are not discussed here, except when included in the original manuscript describing the model.

Models of cardiac cells are capable of facilitating insights into the mechanisms underlying cardiac electrical dynamics. A model begins with a mathematical description of electrical events at the cellular level that give rise to cardiac action potentials. In particular, models incorporate formulations of transmembrane ionic currents along with the voltage, ionic concentrations, and ion channel kinetics responsible for the currents. Because neural and cardiac cells have many similarities, much of the mathematics of cardiac cell modeling is drawn from the pioneering work of Hodgkin and Huxley, who formulated a mathematical description of the giant squid axon (1952).

Modeling cardiac action potentials

Creating a model of the cardiac action potential usually begins with developing a representation of the underlying transmembrane currents. The cell membrane potential or voltage \(V\) then follows an ordinary differential equation \(dV/dt = -I_{\rm ion}/C_m,\) where \(I_{\rm ion}\) is the sum of all transmembrane currents and \(C_m\) is the membrane capacitance. The negative sign is necessary because, by convention, inward currents are negative and outward currents are positive, so that a depolarizing inward current is negative but nonetheless increases the membrane potential.

More detailed models also may incorporate descriptions of intracellular ion concentrations and processes. An overview of modeling these different components is given below, followed by a description of other approaches to modeling cardiac cells.

Modeling transmembrane currents

Action potentials are produced as a result of ionic currents that pass across the cell membrane, producing a net depolarization or repolarization of the membrane as time passes and different currents are invoked in response to the voltage changes. The currents are produced by the movement of individual ions across the membrane through ion channels, which are specialized pore-forming proteins that span the membrane and form a pathway for ions to cross the membrane, as shown in schematic form in Figure 1. The channels open and close stochastically in response to various stimuli that regulate the transport of ions across the membrane. Channels open for only a short time, typically several ms, before closing spontaneously to an inactive state from which they cannot open directly. Although each individual channel is not open for long, the large number of channels (thousands) ensures that ion transport occurs. Each type of channel is highly selective for a specific ion, the most important of which are Na⁺, K⁺, Ca²⁺, and Cl^-, and different types of channels sensitive to a given ion also can exist with different kinetics governing their opening and closing.

Figure 1: Simplified example of a cell membrane separating intracellular and extracellular ions. Interactive Java applet showing the generation of an action potential can be found here.

The most general form of a transmembrane current \(I_{\rm Y}\) permeable to ion Y is simply \(I_{\rm {\rm Y}}=g_{\rm Y} (V-E_{\rm Y}),\) where \(g_{\rm Y}\) is a conductance term, \(V\) is the membrane potential or voltage, and \(E_{\rm Y}\) is the Nernst potential for ion species Y. (The Nernst potential is the potential at which the electrical and chemical gradients across the membrane are balanced for that ion species and produce no movement of ions. The equation to compute the Nernst potential is \(E_{\rm Y}=RT/(F z_{\rm Y}) ln({\rm [Y]_o/[Y]_i}),\) where \(R\) is the universal gas constant, \(T\) is the temperature in Kelvin, \(F\) is the Faraday constant, \(z_{\rm Y}\) is the valence of ion Y, and \({\rm [Y]_o}\) and \({\rm[Y]_i}\) are the concentrations of ion Y outside and inside the membrane, respectively.)

In practice, many cardiac ion channels are “gated” in a voltage-sensitive manner, so that the channels open and close in response to the membrane potential. These gates represent movement of the protein forming the ion channel to allow or to block ion movement (see Fig. 1). Often these gates are represented using the Hodgkin-Huxley formalism, in which the conductance term is decomposed into the product of a maximal conductance term and one or more separate normalized variables that represent the probability of finding the channel open. These variables follow their own differential equations. The most common formulation for a gating variable \(s_i\) is \(ds/dt = (s_{\infty}-s)/\tau_s,\) where \(s_{\infty}\) is the voltage-dependent steady-state value of the gate and \(\tau_s\) is the voltage-dependent time constant of the gate. Some ion channels open and close in response to other factors as well, such as Ca²⁺ concentration.

An alternative representation of gating can be achieved through Markov chains. In this case, each ion channel is represented by \(n\) discrete states representing various configurations of the channel, along with the rates at which transitions can occur from one state to another (these rates may be zero if it is not possible for the channel to switch directly between two given states). The variables for the states then represent the probability of finding the channel in that state. For instance, a simple Markov model of an ion channel may include three states to represent the channel in its open, closed, and inactivated forms, along with six transitions rates that describe the likelihood of a transition from any one state to any other. Each state \(s_i\) follows a differential equation of the following form\[ds_i/dt = \sum_{j=1,j\ne i}^n (k_{ji} s_j- k_{ij} s_i),\] where \(k_{ji}\) is the rate at which a transition from state \(s_j\) to state \(s_i\) occurs (and the rates \(k_{ij}\) and \(k_{ji}\) are zero if a direct connection between states \(s_i\) and \(s_j\) is not permitted). Because the sum of all the states representing a channel must be 1, Markov models are inherently overdetermined.

Along with ion channels, pumps and exchangers also transport ions across the membrane, in this case using active processes rather than simple diffusion. One of the most important of these is the Na⁺/Ca²⁺ exchanger, which under normal circumstances operates primarily to extrude Ca²⁺ ions from the cytoplasm by exchanging one Ca²⁺ ion for three Na⁺ ions, but which also operates in the reverse mode to extrude Na⁺ ions regularly following the rapid influx of Na⁺ ions through the fast Na⁺ channels but also abnormally in various pathophysiological states. Also important is the Na⁺/K⁺ pump, which extrudes Na⁺ that enters the cell from the fast Na⁺ channels by transporting three Na⁺ ions outside the cell and two K⁺ ions inside. Formulations for these currents follow similar principles but along with voltage must account for intracellular and extracellular concentrations of the relevant ions.

Modeling intracellular ion concentrations and processes

The most common intracellular ion concentrations tracked in cardiac models are calcium, sodium, potassium, and occasionally chloride. The calcium concentration, in particular, is quite important, as it is the trigger for the cell to contract, and most models developed since 1977 include at least some basic intracellular calcium dynamics. The intracellular sodium and potassium concentrations are more straightforward to represent and are updated simply in proportion to the sum of the transmembrane currents involving each ion (with appropriate weighting for the number of ions involved in pump and exchanger currents).

Intracellular calcium models attempt to reproduce, at some level, the intricate dynamics of calcium within the cell. An increased level of intracellular calcium initiates the process of cell contraction, and cardiac cells whose primary function is to contract (the working ventricular myocytes) have specialized structures that facilitate this process. Calcium is stored internally in a structure called the sarcoplasmic reticulum (SR). A relatively small influx of calcium through L-type calcium channels in the membrane triggers a comparatively large release of calcium from the internal store that increases the concentration of calcium in the cytosol by approximately a factor of ten. The calcium released binds to other compounds in the cell, which initiates the process of contraction.

Detailed models of calcium handling incorporate these different processes by including the channel through which calcium is released from the SR, ion diffusion within the SR, and reuptake of calcium back into the SR. More recently, some models also include a separate intracellular calcium concentration for the region near the L-type calcium channels, as these channels are predominantly located in a specific area with SR release channels located nearby, and it is hypothesized that a locally increased concentration is sensed in this region before the effects of diffusion equalize the calcium concentration throughout the cell.

Other more detailed ion handling can include buffering of different ions and ion diffusion throughout different parts of the cell.

Simplified cardiac models

Figure 2: Two variable FitzHugh-Nagumo model, with the fast variable shown in red and the slow variable in green. In this simplified cell model, the fast variable is associated with the membrane potential. An interactive single cell applet can be found here, a 1D interactive applet here, and a 2D interactive applet here.

Simpler families of cardiac cell models, based on models of other excitable media, also can be used. For instance, the FitzHugh-Nagumo (FHN) model (FitzHugh, 1961), shown in Figure 2, is a generic model for excitable media and can be applied to a variety of systems. The model can be written as

\(\frac{du}{dt}=(a-u)(u-1)u-v, \frac{dv}{dt}=\epsilon(\beta u - \gamma v - \delta)\ .\)

FitzHugh referred to his simplified model as the Bon Hoeffer-van der Pol model and derived it as a simplification of the Hodgkin-Huxley equations. The model adiabatically eliminated the \(m\) and \(h\) gates and retained only a slow variable \(v\) similar to \(n\ .\) In tissue, a diffusion term is added to the \(u\) equation, which was first done by Nagumo. Because of its simple two-variable form and generality, it has been used widely. For application to cardiac dynamics, it can be modified to prevent the hyperpolarization phase at the end of repolarization and thus can represent a cardiac action potential by changing the \( -v \) term in the equation for \(u\) to \(-cuv\ ,\) where \(c\) is a constant. It can also be modified to represent auto-oscillatory behavior, such as occurs in sinoatrial cells, by shifting the nullclines so that they intersect at an unstable fixed point (which can be done by setting \(\delta\) to a non-zero value).

Despite some similarities with cardiac cells, this model has been found to lack many important properties, including separate time scales for depolarization and repolarization and rate dependence of action potential duration and conduction velocity. Modifications of the time scale ratio \(\epsilon \) and the cubic function into piecewise functions of the voltage in the FHN model can produce rate dependence and linear tip trajectories, respectively, that more closely mimic the observed cardiac dynamics. Nevertheless, this model does not account for a non-zero minimum diastolic interval following an action potential before a new action potential can be produced, which is a key characteristic of cardiac cells. For this reason, and because it does not account for any realistic ionic processes and concentrations, the FHN model is seldom used to represent more than the most general properties of excitation and propagation.

Cellular automata also have been used to represent cardiac cells and can be made to be relatively sophisticated, with a number of states and spatial variations in properties to portray different types of cells. However, cellular automata do not actually produce action potentials, but instead only indicate what part of the action potential the cell is currently experiencing; voltage cannot be tracked explicitly. In addition, cellular automata lack the ability to produce electrotonic currents, which arise from differences in the voltages of neighboring cells and can have profound effects tissue-level behavior. Such limitations have caused cellular automata to be largely replaced by differential equation-based models.

Examples of cardiac cell models

Models of a number of different types of cardiac cells from different species have been developed to represent different regions of the heart as well as different states, including various diseases. The main regions include the atria and the ventricles, as well as specialized pacemaking cells in the sinoatrial node and specialized conducting cells in the His-Purkinje system. The continued discovery of new ion channels and ionic processes in cardiac cells, as well as improvements in voltage-clamp techniques and data acquisition, has led to an increase over time in the complexity of mathematical models used to describe cardiac cell electrical dynamics. Below we list a number of important cell models according to their region of the heart and summarize their properties.

Purkinje cell models

Noble model (1962) [4 variables]

Figure 3: Example of action potentials using the Noble model. An interactive single cell Java applet of this model can be found here and a 1D Java applet here.

Figure 4: Alternans in the Noble model. Action potential duration (APD) vs basic cycle length (BCL) or period. As the BCL is decreased, alternans arises at a BCL of around 260 ms and terminates at a BCL of about 120 ms. The gray area shows BCLs where propagation fails in 1D and 2D .

The first model of a cardiac cell is the Noble model (1962) of a generic Purkinje cell, shown in Figure 3. It includes three primary currents, a sodium, a potassium, and a background current assumed to be mediated by chloride ions, and uses a total of four variables, the voltage, sodium current activation and inactivation gates, and a potassium current activation gate. The equations are as follows\[ \frac{\text{d}V}{\text{d}t} = -(I_{Na}+I_{K}+I_{bk})/C_m , \quad \] \(I_K=(g_{k1}+g_{k2})(V_m+V_k) , \quad \) \( I_{Na}=(400m^3h +0.14)(V_m-40) , \quad \) and \( I_{bk} = g_{bk}(V-E_{bk}), \) where \( g_{k1}=1.2exp[(-V_m-90)/150]+0.015 exp[(V_m+90)/60], \quad \) \( g_{k2}= 1.2n^4,\) and \(\frac{\text{d}y}{\text{d}t} = \alpha_y(V_m)(1-y)-\beta_y(V_m)y \) for \(\quad y=[m,h,n] \ .\) (See the Noble model for further details on the parameter values.)

This model is also notable because it predated the discovery of calcium channels in cardiac cells, so that the action potential results from a balance of inward sodium and outward potassium currents. The model is auto-oscillatory with a period of approximately 900 ms and exhibits a period-doubling bifurcation resulting in alternans in action potential duration (APD) when paced rapidly (see Figure 4). In this model, the auto-oscillatory behavior can be eliminated through various means, such as changing the constant sodium current conductance from 0.14 to 0.132 or by slightly increasing the potassium current conductance \(g_k\ .\) One main difference between this and other models of cardiac action potential is the value of the membrane capacitance used, which is much larger than the value used in other models.

McAllister-Noble-Tsien model (1975) [10 variables]

An updated Purkinje cell model was published in 1975. The McAllister-Noble-Tsien (MNT) model includes the inward calcium current (which had been discovered since the 1962 Noble model) along with a sodium current, three potassium currents, one chloride current, and three background currents mediated by K⁺, Na⁺, and Cl^-. A total of 10 variables are used, consisting of the voltage and nine gates.

DiFrancesco-Noble model (1985) [16 variables]

Figure 5: Example of action potentials using the DiFrancesco-Noble model. Interactive applet can be found here.

The most recent Purkinje cell model, shown in Figure 5, was published by DiFrancesco and Noble (1985) as an update of the MNT model.

Figure 6: Action potential duration (APD) as a function of basic cycle length (BCL) or period. As the BCL is decreased, no alternans is observed, unlike the Noble model (shown as dotted lines).

It incorporates a number of important advances and was the first model to include varying intracellular sodium and intracellular and extracellular potassium concentrations, more sophisticated calcium handling, and additional currents like the Na⁺-K⁺ pump and Na⁺-Ca²⁺ exchanger. The model uses 16 variables and represents the voltage; intracellular Ca²⁺, K⁺, and Na⁺ concentrations; extracellular K⁺; uptake and release Ca²⁺ concentrations for the sarcoplasmic reticulum (SR); and nine gating variables. Along with the Na⁺-K⁺ pump and Na⁺-Ca²⁺ exchanger, the model includes Na⁺, Ca²⁺, transient outward, time-dependent and time-independent K⁺, Na⁺ and Ca²⁺ background, and K⁺/Na⁺-mediated pacemaker currents, for a total of 10 currents. The model remains auto-oscillatory but does not exhibit alternans when paced at any cycle length (see Figure 6), in contrast to the Noble model and experimentally observed behavior.

Karma model (1993) [2 variables]

Figure 7: Example of action potentials using the Karma model. Fast variable corresponding to the membrane potential (red) and the slow gate variable (green) are shown. A single cell Java applet can be found here and a 2D Java applet that can produce spiral waves and breakup can be found here.

The Karma model is based on an analysis of the Noble model where the fast gate variables \(m\) and \(h\) have been eliminated adiabatically. The model consists of two variables, a fast variable related to the membrane voltage and a slow gate variable (see Figure 7). This model was first developed to show spiral wave breakup due to alternans (Karma, 1993), as shown in Figure 8. The model equations are as follows:

Figure 8: Spiral breakup due to alternans

\( \frac{dV_m}{dt}=-V_m+[1+4\delta-(w/w_B)^M]h(V_m);\quad\) \( \frac{dw}{dt}=\epsilon[\theta(V_m-1)-w] \ ,\) where \(\theta(V_m-1)\) is a Heaviside step function and \(\epsilon\) is the ratio of time scales between the action potential upstroke and its duration. The action potential generated with this model has a more realistic fast upstroke and slow recovery compared to that of the standard FHN model (see Figure 2). Furthermore, it reproduces the oscillatory pulse dynamics (alternans) observed in experiments with cardiac rings. These dynamics are a known mechanism for spiral wave breakup, as alternans in this case produces changes in the wavelength that can result in propagation block, wave breaks, and multiple spiral waves (Karma 1993). The model, however, is constrained by the fact that the minimum action potential duration and the minimum conduction velocity approach zero in the limit \(\epsilon\) \(\rightarrow 0\ ,\) while experiments and simulations of more detailed ionic models show that both of these quantities are finite in this limit. Thus, this model fails to reproduce certain physiological effects, such as linear trajectories for spiral waves similar to those observed experimentally.

Ventricular models

Following the early models of the Purkinje cell, the next region of the heart to be modeled was the ventricles, with the publication in 1977 of the Beeler-Reuter model of a generic ventricular cell. The Beeler-Reuter model included the relatively recently discovered calcium current, along with an intracellular calcium concentration designed to produce a realistic calcium transient that varied over the course of the action potential.

In the decades after the Beeler-Reuter model was published, models that were more specific to particular species began to emerge. Such differences are necessary because different species express different types of ion channels with different kinetics and can produce vastly different action potentials with different properties.

Generic models

Beeler-Reuter model (1977) [8 variables]

Figure 9: Example of action potentials using the Beeler-Reuter model. A single cell interactive Java applet can be found here and an interactive 1D cable can be found here.

Figure 10: Spiral breakup using the Beeler-Reuter model.

Figure 11: Stable spiral wave using the modified Beeler-Reuter model.

The first ventricular model was published by Beeler and Reuter in 1977 and is shown in Figure 9. This model uses four of the eight different ionic currents known at the time in cardiac muscle. They implemented a fast inward Na⁺ current I_Na, similar to the one used by Hodgkin and Huxley, but they added a second slower inactivation gate \(j\ ,\) the time-dependent outward current I_x1, a time-independent K⁺ outward current I_K1, and a slow inward current I_s carried primarily by Ca²⁺. The total ionic current in the Beeler-Reuter model is given by four currents, and the model uses eight variables.

The Beeler-Reuter model was the first ionic model simulated in 2D, by Courtemanche and Winfree (1991). It not only showed for the first time that ionic models could produce spiral waves, but that in fact they could break up into multiple waves (see Figure 10. In this model, breakup is due to slow recovery fronts (Courtemanche, 1996). The breakup can be suppressed by making the calcium dynamics faster to reduce the large variation in wavelengths. Figure 10 shows a stable spiral wave with a hypermeandering trajectory using what is know as the modified Beeler-Reuter (MBR) model, which is obtained when the calcium dynamics is made faster by a factor of 2 (by dividing the activation and inactivation gates \( d \) and \(f\) by a factor of two) (Courtemanche and Winfree, 1991; Courtemanche, 1996).

Fenton-Karma model (1998) [3 variables]

Figure 12: Example of action potentials in the three variable model. A single cell interactive Java applet can be found here and an interactive one dimensional ring applet here.

The Fenton-Karma model (also known as the 3V-model) as developed in 1998 as a three-variable model of the cardiac action potential. The model uses three transmembrane currents, a fast inward, slow inward, and slow outward, which represent summary Na⁺, Ca²⁺, and K⁺ currents, respectively. Along with the voltage, two gating variables are used to regulate inactivation of the fast inward and slow inward currents. In addition to reducing computational requirements, one of the model’s main contributions is its flexibility in reproducing data from experiments or from other cardiac action potential models. Some of the main action potential features that can be set in a straightforward manner in the model are the threshold for excitation, action potential rate of rise, minimum diastolic interval, minimum and maximum action potential durations, and action potential duration and conduction velocity restitution curves. The original article (Fenton and Karma, 1998) introduces four parameter sets, with values set to reproduce these and other properties of the Beeler-Reuter, modified Beeler-Reuter, and modified Luo-Rudy 1 models, along with guinea pig experimental data. A subsequent paper (Fenton et al., 2002) introduced additional parameter sets designed to possess particular electrophysiological properties, and then used these parameter sets to demonstrate ten different mechanisms of spiral wave breakup.

The Fenton-Karma model and its four-variable counterpart (Cherry and Fenton, 2004; Cherry and Fenton, 2007; Bueno-Orovio et al., 2008) are useful alternatives to more complex models, especially in large-scale simulations. This model preserves important properties of cardiac tissue including action potential rate of rise, different time scales for depolarization and repolarization, action potential duration and conduction velocity restitution curves, nonzero minimum diastolic interval, and action potential shape (for the four-variable version). In comparison with FitzHugh-Nagumo-type models, this model does not have a Maxwell point, and thus its dynamics are comparable to complex ionic models as they have a well-defined minimum diastolic interval and minimum action potential duration.

Guinea pig models

Luo-Rudy 1 model (1991) [8 variables]

Figure 13: Example of action potentials using the Luo-Rudy 1 model. A single cell interactive Java applet can be found here and an interactive 1D cable applet here.

Figure 14: Spiral wave breakup using the modified Luo-Rudy 1 model.

In 1991 Luo and Rudy (1991) published an ionic model (LR1) for the cardiac action potential in guinea pig ventricular cells based on the Beeler-Reuter model (Luo and Rudy, 1991), but updated to include more recent experimental results. The model, shown in Figure 13, reformulated the opening and closing rate coefficients for the sodium current from the BR-model, making it a faster process, and added three new currents, a plateau K⁺ current, a background current with constant conductance and an additional K⁺ current with a gate variable that can be approximated by its steady-state value due to a small time constant. They retained the BR formulation for the slow inward calcium current as well as the time-dependent potassium current. In total, the LR1 model describes six different currents and uses nine variables, one of which is approximated by its steady state and is replaced by a function, so that only eight variables are needed in the calculation.

Like the BR model, the LR1 model exhibits breakup in 2D, but the breakup remains even when the calcium dynamics are made faster by a factor of two in the modified Luo-Rudy 1 model (MLR1), as shown in Figure 14). However, in this case the breakup is produced by a Doppler effect due to a very small minimum diastolic interval and very high excitability (Fenton et al., 2002).

Nordin model (1993) [14 variables]

The Nordin model was derived from the DFN model and whenever possible was fitted using data available from guinea pig ventricular myocytes. A unique feature of the model is the separation of the bulk Ca²⁺ concentration into three regions, superficial, mid, and deep, which was done to support Ca²⁺ intake into the cell and internal release while maintaining an appropriate whole-cell transient. The model includes a total of 11 currents: Na⁺; Ca²⁺; three K⁺ currents; Na⁺-Ca²⁺ exchanger; Na⁺-K⁺ pump; nonspecific Ca²⁺-activated cation current; and background Na⁺, K⁺, and Ca²⁺ currents. Overall, 14 variables were used, including the voltage, six gating variables, intracellular Na⁺ and K⁺ concentrations, three distinct regionally-based intracellular Ca²⁺ concentrations, and Ca²⁺ concentrations in the SR uptake and release compartments.

Luo-Rudy 2 model (1994) [15 variables]

The Luo-Rudy phase 2 (LR2) model was developed as an extension of the earlier LR1 model. It includes a more detailed description of intracellular Ca²⁺ processes, including a two-compartment representation of the SR, pump and exchanger currents, Ca²⁺ buffering in the cytoplasm and SR, and Ca²⁺-induced Ca²⁺ release. Intracellular Na⁺ and K⁺ concentrations also can be tracked. Overall, the model includes 11 currents: the fast Na⁺, L-type Ca²⁺, delayed rectifier K⁺, inward rectifier K⁺, plateau K⁺, Na⁺-Ca²⁺ exchanger, Na⁺-K⁺ pump, sarcolemmal Ca²⁺ pump, background Na⁺ and K⁺, and nonspecific Ca²⁺-activated currents. A total of 15 variables are used: voltage; 6 gating variables; cytoplasmic Na⁺, K⁺, and Ca²⁺ concentrations; Ca²⁺ concentrations in the JSR and NSR; and buffer occupancies for troponin and calmodulin in the cytoplasm and for calsequestrin in the SR.

Luo-Rudy dynamic model (1995 - )

Figure 15: Example of action potentials using the Luo-Rudy dynamic model. A single cell interactive Java applet can be found here and an interactive 1D cable can be found here.

The Luo-Rudy dynamic (LRd) model, shown in Figure 15, essentially consists of any and all updates to the LR2 model. As such, its formulation and components change over time. Among the changes that have been made are the separation of the delayed rectifier \( K^+\) current into rapid and slow components and a different method to compute \(Ca^{2+}\) buffering (Zeng et al., 1995); subsequent reformulation of the slow component along with modification of the Ca²⁺-induced Ca²⁺ release mechanism (Viswanathan et al., 1999); and formulation of a Na⁺-activated K⁺ current, updated Na⁺-Ca²⁺ exchanger current, and modified SR Ca²⁺ release (Faber and Rudy, 2000). Since 1999 the model can account for the properties of three regions (epicardial, midmyocardial, and endocardial) through variations in parameter values. Markov descriptions of some ion channels also have been developed and can be used to reproduce channel dysfunction associated with genetic diseases, such as the Brugada Syndrome (Clancy and Rudy, 2001).

Matsuoka-Sarai-Kuratomi-Ono-Noma model (2003) [45 variables]

The Matsuoka et al. (2003) model (also known as the Kyoto Model) was developed as a “compound model” capable of simulating either guinea pig ventricular or rabbit sinoatrial nodal cell by adjusting various model parameter values. The model incorporates the contraction model of Negroni and Lascano (1996) to simulate sarcomere shortening and contractile force generation. A total of 54 different variables are used, 45 of which are used in the ventricular model, including voltage; 4 states for the fast Na⁺ channel; 4 states for the L-type Ca²⁺ channel; 4 states for the Ca²⁺-dependent inactivation of the L-type Ca²⁺ channel; 3 states for the RyR channel; 15 gating variables; Ca²⁺ concentrations in the cytoplasm, JSR, and NSR; intracellular Na⁺, K⁺, and ATP concentrations; 2 variables to track sarcomere shortening as a result of contraction; 4 variables to track Ca²⁺ buffering by troponin and attachment to cross-bridges; and 2 variables to track buffering of Ca²⁺ by calmodulin in the cytoplasm and by calsequestrin in the SR. The model contains expressions for 17 different transmembrane currents, 14 of which are used for the ventricular model: the fast Na⁺, L-type Ca²⁺, T-type Ca²⁺, transient outward K⁺, rapid and slow delayed rectifier K⁺, inward rectifier K⁺, Na⁺-Ca²⁺ exchanger, Na⁺-K⁺ pump, and five different background currents.

Human models

Priebe-Beuckelmann model (1998) [17 variables]

Figure 16: Example of action potentials using the Priebe-Beuckelmann model. An interactive Java applet can be found here.

Figure 17: Spiral wave dynamics of the Priebe- Beuckelmann model.

The first human ventricular model was published by Priebe and Beuckelmann in 1998 and included specific variations to describe both normal and heart failure-remodeled electrophysiology (see Figure 16). It was based on the LR2 model and retained much of the intracellular calcium handling description from that model. The Priebe-Beuckelmann model includes a total of 10 currents: the L-type Ca²⁺, transient outward K⁺, rapid and slow delayed rectifier K⁺, inward rectifier K⁺, fast Na⁺, Na⁺-Ca²⁺ exchanger, Na⁺-K⁺ pump, and background Na⁺ and Ca²⁺ currents, with human data used to derive the first five of these. A total of 17 variables are used, including the voltage, 9 gate variables, intracellular Na⁺ and K⁺ concentrations, total and free intracellular Ca²⁺ concentrations, total and free junctional SR Ca²⁺ concentrations, and the network SR Ca²⁺ concentration. In 2D the model produces quasi-breakup (Bueno-Orovio et al., 2008) with the same mechanisms of the LR1 model, as shown in Figure 17.

Bernus-Wilders-Zemlin-Verschelde-Panfilov model (2002) [6 variables]

Figure 18: Example of action potentials using the Bernus et al. model. An interactive Java applet can be found here.

Figure 19: Spiral wave dynamics of the Bernus et al. model.

A reduction of the Priebe-Beuckelmann model was published by Bernus et al. in 2002. The same currents were included, except that the delayed rectifier was treated as a single current rather than as rapid and delayed components. The number of variables, however, was greatly reduced, to a total of six (see Figure 18), including the voltage and five gating variables to govern activation of the fast Na⁺ and delayed rectifier K⁺ currents as well as inactivation of the fast Na⁺, transient outward K⁺, and L-type Ca²⁺ currents. The remaining gates from the original Priebe-Beuckelmann model were either combined or adiabatically eliminated, while all ionic concentrations were removed as variables. The model also includes variations to reproduce the action potentials corresponding to epicardial, endocardial, and midmyocardial cells. In 2D, like the original Priebe-Beuckelmann model, the Bernus et al. model also produces quasi-breakup (see Figure 19), where the spiral breaks and recombines continuously (Bueno-Orovio et al., 2008) with a dominant frequency of 286 ms.

Ten Tusscher-Noble-Noble-Panfilov model (2004) [17 variables]

Figure 20: Example of action potentials using the Ten Tusscher et al. model. An interactive Java applet can be found here.

Figure 21: Spiral wave dynamics of the ten Tusscher et al. model.

The Ten Tusscher et al. (2004) model includes descriptions for human epicardial, endocardial, and midmyocardial cells (see Figure 20). In total, it includes 12 transmembrane currents: fast Na⁺, L-type Ca²⁺, rapid and slow components of the delayed rectifier K⁺, inward rectifier K⁺ , transient outward K⁺, plateau K⁺, Na⁺-Ca²⁺ exchanger, Na⁺-K⁺ pump, sarcolemmal Ca²⁺ pump, and background Na⁺ and Ca²⁺ currents. The model uses 17 variables: voltage, 12 gating variables, intracellular Na⁺ and K⁺ concentrations, and Ca²⁺ concentrations in the cytoplasm and SR. The model also has been updated to include a fuller description of calcium dynamics (ten Tusscher and Panfilov, 2006a) and has been simplified to a total of 9 variables (voltage and 8 gates) (ten Tusscher and Panfilov, 2006b). In 2D, the model produces stable spiral waves using initial conditions from pacing at a 1-second cycle length, with a frequency of 259 ms (see Figure 21. However, it can produce spiral breakup if a spiral wave is initiated after pacing at 300ms or faster (Bueno-Orovio et al., 2008).

Iyer-Mazhari-Winslow model (2004) [ 67 variables]

Figure 22: Example of action potentials using the Iyer et al. model. An interactive Java applet can be found here.

Figure 23: Spiral wave dynamics of the Iyer et al. model.

The Iyer-Mazhari-Winslow model (see Figure 22) provides a description of the human epicardial ventricular myocyte and is notable for including Markov chain representations for six transmembrane currents as well as the RyR channel. The model incorporates 13 transmembrane currents, which are the same as those of the Winslow et al. (1999) canine ventricular model except that the plateau K⁺ current has been removed and the transient outward K⁺ current has been separated into fast and slow components. Because of the extensive use of Markov descriptions, a large number of variables is used, 67 in all: voltage; 13 states for the fast Na⁺ channel; 10 states each for the fast and slow transient outward K⁺ channels; 5 states for the rapid component of the delayed rectifier K⁺ channel; 4 states for the slow component of the delayed rectifier K⁺ channel; 11 states for the L-type Ca²⁺ channel; 4 states for the RyR channel; 1 gating variable; Ca²⁺ concentrations in the cytoplasm, JSR, NSR, and restricted subspace; intracellular Na⁺ and K⁺ concentrations; and fraction of Ca²⁺-bound high and low affinity troponin sites. In 2D, the model produces spiral wave breakup (see Figure 23) with relatively small wave lengths and a dominant period of 98 ms.

Bueno-Orovio-Cherry-Fenton model [4 variables]

Figure 24: An example of action potentials using the Bueno-Orovio et al. model. An interactive Java applet can be found here.

The Bueno-Orovio et al. (2008) model (see Figure 24) provides a description of the human ventricular cells. Separate parameter sets describing epicardial, endocardial, and midmyocardial cells are provided, with endocardial action potentials shorter than epicardial and midmyocardial action potentials significantly longer than both epicardial and endocardial. Following a formulation similar to that of the Fenton-Karma model (1998), the model incorporates 3 transmembrane currents, a fast inward, slow inward, and slow outward, which represent summary Na⁺, Ca²⁺, and K⁺ currents. Only four variables are used, including the voltage and three gating variables. In 2D, the model produces a stable spiral wave (see Figure 25) with a dominant cycle length of 295 ms.

When different parameter sets are used, the model also can reproduce the dynamics of the Priebe-Beuckelmann and Ten Tusscher et al. human ventricular models (Bueno-Orovio et al., 2008), including the dynamics of reentrant waves in two-dimensional tissue. In addition, parameter changes to incorporate a conduction velocity restitution curve based on data obtained in humans are provided. These modifications result in a longer wavelength and reduced minimum DI and produce quasi-stable reentrant wave dynamics, similar to the Priebe-Beuckelmann model.

Figure 25: Spiral wave dynamics of the Bueno-Orovio et al. model.

Canine models

Winslow-Rice-Jafri-Marban-O'Rourke model (1999) [33 variables]

The first canine ventricular model was developed by Winslow et al. model and is based on the Jafri et al. (1998) guinea pig ventricular model, with updated current descriptions to represent canine ventricular cells. The model is designed to represent a midmyocardial cell, and modifications to represent the electrophysiological remodeling induced by tachycardia pacing are presented. The 13 transmembrane currents used in the model include the fast Na⁺, rapid and slow components of the delayed rectifier K⁺, transient outward K⁺, inward rectifier K⁺, plateau K⁺, L-type Ca²⁺, K⁺ through the L-type Ca²⁺ channel, Na⁺-Ca²⁺ exchanger, Na⁺-K⁺ pump, sarcolemmal Ca²⁺, and background Na⁺ and Ca²⁺ currents. A total of 33 variables are used, including the voltage; 8 gating variables; 12 states for the L-type Ca²⁺ channel; 4 states for the RyR channel; Ca²⁺ concentrations in the cytoplasm, NSR, JSR, and restricted subspace; intracellular Na⁺ and K⁺ concentrations; and concentrations of Ca²⁺-bound high and low affinity troponin sites in the cytoplasm.

Fox-McHarg-Gilmour model (2002) [13 variables]

Figure 26: Example of action potentials using the Fox et al. model. An interactive Java applet can be found here.

Figure 27: Spiral wave dynamics of the Fox et al. model.

Figure 28: Spiral wave dynamics of the Fox et al. model with I_Kr increased by a factor of two.

The Fox et al. (2002) model (see Figure 26) is based on the Winslow et al. (1999) model, with modifications to ensure that alternans occurred at fast pacing rates. The model incorporates a simplified description of intracellular calcium handling based on Chudin et al. (1999). The same 13 currents as the Winslow et al. model are included, but the model uses only 13 variables, including voltage, cytoplasmic, and sarcoplasmic reticulum Ca²⁺ concentrations, and 10 gates. In 2d, an initiated spiral wave may experience transient breakup but reaches a steady state of a stable spiral wave with a slowly precessing linear trajectory (see Figure 27).

Increasing the contribution of rapid delayed rectifier I_Kr eliminates alternans that occurs in a single cell. However, in a one-dimensional cable, alternans still can be observed. Moreover, in two-dimensional tissue, this modification leads to sustained breakup (see Figure 28), in contrast to the original model. Breakup also occurs when I_Kr is decreased by a factor of two.

Cabo-Boyden model (2003) [16 variables]

Cabo and Boyden (2003) published a model of the canine epicardial cell based closely on the Luo-Rudy 2 model, with variations for both normal and infarcted cells. A total of 13 transmembrane currents are used, with canine-specific descriptions for the fast Na⁺, L-type Ca²⁺, rapid and slow components of the delayed rectifier ⁺ (treated as separate currents, unlike in the LR2 model), and inward rectifier K⁺ currents substituted for the corresponding descriptions in the LR2 model. A transient outward K⁺ following the PB model was added, and the remaining 7 currents were left unchanged from the LR2 model. A total of 16 variables was used, one more than for the LR2 because of the separate treatment of the rapid and slow components of the delayed rectifier K⁺ currents. The IZ model incorporated an additional gate for the L-type Ca²⁺ current to model slower recovery from inactivation under these conditions.

Hund-Rudy model (2004) [29 variables]

Figure 29: Example of action potentials using the Hund-Rudy model. An interactive Java applet can be found here.

Figure 30: Spiral wave dynamics of the Hund-Rudy model.

The Hund-Rudy model (Hund and Rudy, 2004) (see Figure 29) was developed to describe canine ventricular epicardial cells. It incorporates more detailed calcium handling description than many previous models, especially through the addition of autophosphorylation through the Ca²⁺/calmodulin-dependent protein kinase (CaMKII), and includes an intracellular Cl^- ion concentration along with Cl^- and late Na ⁺ currents. A total of 14 transmembrane currents, including fast Na⁺, rapid and slow components of the delayed rectifier K⁺, transient outward K⁺, inward rectifier K⁺, plateau K⁺, L-type Ca²⁺, Na⁺-Ca²⁺ exchanger, Na⁺-K⁺ pump, sarcolemmal Ca²⁺, late Na⁺, Ca²⁺-activated Cl^-, and background Ca²⁺ and Cl^- currents. A total of 29 variables are used: voltage, 17 gating variables for the transmembrane currents; 2 gating variables for the Ca²⁺ release current from the SR; a variable power to which the activation gate for the L-type Ca²⁺ current is raised; fraction of autonomous CaMKII binding sites with trapped calmodulin; Ca²⁺ concentrations in the cytoplasm, restricted subspace, NSR, and JSR; and intracellular Na⁺, K⁺, and Cl^- concentrations. In 2d, spiral waves are hypermeandering and can remain stable or experience breakup depending on the initial conditions used. An example of a stable spiral wave is shown in Figure 30.

Greenstein-Winslow model (2002)

The Greenstein-Winslow model is an updated version of the Winslow et al. canine ventricular model that was designed to incorporate more detailed descriptions of intracellular Ca²⁺ handling. In particular, the model implements explicit local control of Ca²⁺ release by the inclusion of distinct calcium release subunits and through this mechanism intends to mimic the behavior of Ca²⁺ sparks in the junction of the SR and T tubule. Each release subunit consists of four dyadic subspace volumes, and each dyadic susbspace volume, in turn, includes a single L-type Ca²⁺ channel, five RyR channels, and a single Ca²⁺-activated Cl^- channel. At the whole cell level, the Ca²⁺-activated Cl^- current is included, along with Markov chain models of the fast component of the delayed rectifier K⁺ current and the fast and slow transient outward K⁺ currents. An exact number of variables cannot be given because the number of calcium release units is itself variable, generally only a fraction of these release units are simulated, and certain components like the Markov chain representing the RyR channel may contain different numbers of states depending on the Ca²⁺ concentration. Nonetheless, because approximately 80 variables are needed to represent each release unit, simulation of even a moderate number of release units (on the order of 1000) is quite computationally intense.

A simplified version of the model that is more computationally tractable was published by Greenstein et al. (2006). The revision includes a total of 76 variables, with 40 representing intracellular Ca²⁺ states, and the other 36 representing the voltage, gate variables, and ion concentrations.

The simplified model was used as the basis for the Flaim et al. model (2007), which also incorporated variations to represent epicardial, endocardial, and midmyocardial cells. The primary changes to the model formulation include modifications to the transient outward K⁺ current, inclusion of the late Na⁺ current through a Markov representation for I_Na, and transmural variations in the conductance of the slow component of the delayed rectifier I_Ks and in the forward and reverse rates of the SERCA pump. This version includes 87 total variables.

Rabbit models

Puglisi-Bers model (2001) [20 variables]

The first model of the rabbit ventricular myocyte was published by Puglisi and Bers (2001). It was based on the LRd model, with appropriate substitutions made for rabbit ventricular currents, and was one of the first models to include the Ca²⁺-activated Cl^- current. The model included 14 transmembrane currents: fast Na⁺, L-and T-type Ca²⁺, rapid and slow components of the delayed rectifier K⁺, inward rectifier K⁺, plateau K⁺, transient outward K⁺, Ca²⁺-activated Cl^-, Na⁺-Ca²⁺ exchanger, Na⁺-K⁺ pump, sarcolemmal Ca²⁺ pump, and background Na⁺ and K⁺ currents. A total of 20 variables are used: voltage; 11 gating variables; intracellular Na⁺ and K⁺ concentrations; Ca²⁺ concentrations in the cytoplasm, JSR, and NSR; and buffer occupancies for troponin and calmodulin in the cytoplasm and for calsequestrin in the SR.

Shannon-Wang-Puglisi-Weber-Bers model (2004) [45 variables]

Figure 31: Action potentials of the Shannon et al. model at 1s and 200ms cycle lengths.

Figure 32: Spiral wave dynamics of the Shannon et al. model.

The Shannon et al. (2004) model (see Figure 31) includes several novel components, such as separate distributions of ion channels within the sarcolemma and the junctional space along with region-specific ionic concentrations, a modified formulation of the Na⁺-Ca²⁺ exchange current, and updated intracellular Ca²⁺ handling. The model incorporates 14 transmembrane currents, including fast Na⁺; L-type Ca²⁺; rapid and slow components of the delayed rectifier K⁺; inward rectifier K⁺; fast and slow transient outward K⁺; Ca²⁺-activated Cl^-; Na⁺-Ca²⁺ exchanger; Na⁺-K⁺ pump; sarcolemmal Ca²⁺ pump; and background Na⁺, K⁺, and Cl^- currents. A total of 45 variables are used: voltage; 13 gating variables; 4 RyR states; free and bound concentrations of Na⁺ in the sarcolemma and junction; cytosolic Na⁺ concentration; concentrations of Ca²⁺ in the cytosol, SR, sarcolemma, and junction; and 9, 1, 4, and 4 Ca²⁺ buffers in the cytosol, SR, sarcolemma, and junction, respectively. Because the dyes indo-1 and fluo-3 are not present natively in rabbit ventricular myocytes, the corresponding buffers for each dye in the cytosol, sarcolemma, and junction can be removed in many cases, thereby decreasing the total number of variables to 39. In 2d, an initiated spiral wave experiences slow recovery fronts as described by Courtemanche (1996), but otherwise remains stable with a linear core, as shown in Figure 32.

Mahajan-Shiferaw et al. model (2008) [27 variables]

Figure 33: Action potentials of the Mahajan et al. model at 1s and 400ms cycle lengths.

Figure 34: Spiral wave dynamics of the Mahajan-Shiferaw et al. model.

The Mahajan et al. model (Mahajan et al., 2008) was designed to reproduce rabbit action potentials ( Figure 33) over a wider range of pacing rates, including rapid rates relevant to ventricular tachyarrhythmias. It draws largely from the Shannon et al. model, with modifications of the L-type calcium current and intracellular calcium cycling. An important dynamical feature of this model is its ability to generate alternans at fast pacing rates, similar to experimental observations of alternans. Both steep restitution and calcium dynamics can lead to alternans in this model. Nine transmembrane currents are used, including fast Na⁺; L-type Ca²⁺; rapid and slow components of the delayed rectifier K⁺; inward rectifier K⁺; fast and slow transient outward K⁺; Na⁺-Ca²⁺ exchanger; and Na⁺-K⁺ pump currents. A total of 27 variables are used: voltage; 10 gating variables; 7 Markov states for the L-type calcium current; 2 troponin buffers; average concentrations of calcium in the cytosol, submembrane space, SR, and active dyadic clefts; average free calcium available for release in the JSR; the calcium release current; and the intracellular sodium concentration. This models has alternans when paced at fast frequencies, however an initiated spiral wave in 2d is stable and follows long linear cores (larger than those to the Shannon-Wang-Puglisi-Weber-Bers model) as shown in Figure 34.

Rat models

Pandit-Clark-Giles-Demir model (2001) [26 variables]

Figure 35: Action potentials of the Pandit et al. model at 1s and 300ms cycle lengths.

Figure 36: Spiral wave dynamics of the Pandit et al. model.

The Pandit et al. model (see Figure 35) is the first published rat ventricular model. A total of 12 currents are represented, including the fast Na⁺; L-type Ca²⁺; transient outward K⁺; steady state K⁺; inward rectifier K⁺; hyperpolarization-activated; Na⁺-Ca²⁺ exchanger; Na⁺-K⁺ pump; sarcolemmal Ca²⁺ pump; and background Na⁺, K⁺, and Ca²⁺ currents. The model uses 26 variables: voltage; 13 gating variables; intracellular Na⁺ and K⁺ concentrations; Ca²⁺ concentrations in the cytoplasm, restricted subspace, JSR, and NSR; 4 RyR states; and concentrations of Ca²⁺-bound high and low affinity troponin sites in the cytoplasm. It includes variants for epicardial and endocardial cells in the left ventricle, with a smaller density and slower reactivation kinetics of the transient outward K⁺ current in the endocardial cells contributing to shorter action potentials and increased rate dependence in endocardial cells.

At fast pacing rates, the model can exhibit a secondary less polarized fixed point. In two-dimensional tissue sheets, this additional fixed point may contribute to the development of spiral wave meandering and breakup, which arises after several rotations (see Figure 36).

Mouse models

Bondarenko-Szigeti-Bett-Kim-Rasmusson model (2004) [44 variables]

Figure 37: Action potentials in the Bondarenko et al. model at 1s and 50ms cycle lengths for the apex and septum parameter sets.

Figure 38: Spiral wave dynamics of the Bondarenko et al model.

The first model of mouse ventricular myocytes was published by Bondarenko et al. in 2004 and is shown in Figure 37. A total of 15 transmembrane currents are incorporated, including the fast Na⁺; L-type Ca²⁺; rapid, slow, and ultrarapid delayed rectifier K⁺; non-inactivating steady-state K⁺; inward rectifier K⁺; fast and slow components of the transient outward K⁺; calcium-activated transient outward Cl^-; Na⁺-Ca²⁺ exchanger; Na⁺-K⁺ pump; sarcolemmal Ca²⁺ pump; and background Na⁺ and Ca²⁺ currents. The model includes 44 variables: voltage; bulk and dyadic cleft intracellular Ca²⁺, junctional and network SR Ca²⁺, intracellular Na⁺ and K⁺, and myoplasmic troponin low- and high-affinity site Ca²⁺ bound concentrations; 8 L-type Ca²⁺ channel states; 9 fast Na⁺ channel states; 5 rapid delayed rectifier channel states; 4 ryanodine receptor (RyR) channel states; a modulation factor for the RyR Ca²⁺ release current; and 8 gating variables governing activation of the fast and slow transient outward, slow delayed rectifier, ultrarapid delayed rectifier, and steady-state K⁺ currents and inactivation of the fast and slow transient outward and ultrarapid delayed rectifier K⁺ currents. Variations for apical and septal cells are provided, as shown in Figure 37. In 2D, the model exhibits stable spiral waves (see Figure 38.

Atrial cell models

Rabbit models

Hilgemann-Noble model (1987)

The Hilgemann-Noble model (1987) was one of the earliest atrial models. It represented a multicellular preparation and was later scaled down to represent a single cell by Earm and Noble (1990). It focused on investigating a number of Ca²⁺ processes, including buffering, Ca²⁺-induced Ca²⁺release, both Ca²⁺- and voltage-dependent inactivation processes for the Ca²⁺ current, and the SR Ca²⁺ pump, while time- and Ca²⁺-dependence of K⁺ conductances were not simulated. The number of variables depends on the formulation chosen, as the model presents a number of options in terms of which processes are included and which descriptions of the included processes are used.

Lindblad-Murphy-Clark-Giles model (1996) [28 variables]

The Lindblad et al. model (1996) contains a more extensive description of the rabbit atrial cell and includes 12 transmembrane currents: fast Na⁺; L- and T-type Ca²⁺; transient outward K⁺; rapid and slow components of the delayed rectifier K⁺; inward rectifier K⁺; Na⁺-Ca²⁺ exchanger; Na⁺-K⁺ pump; Ca²⁺ sarcolemmal pump; and background Na⁺, K⁺, and Cl^- currents. A total of 28 variables are used including the voltage; 14 gating variables; intracellular Na⁺ and K⁺ concentrations; Ca²⁺ concentrations in the cytoplasm, JSR, and NSR; 3 variables used to determine Ca²⁺ release from the SR; and five variables used to describe Ca²⁺ buffering by troponin and calmodulin in the cytoplasm and by calsequestrin in the SR.

Bullfrog models

Rasmusson-Clark-Giles-Robinson-Clark-Shibata-Campbell model (1990) [16 variables]

One of the earliest atrial models developed was that of Rasmusson et al. (1990) for bullfrog atrial cells. Some noteworthy features of this model include the formulations for buffering in the cytoplasm along with varying extracellular ionic concentrations. This model contains 8 currents: fast Na⁺, Ca²⁺, delayed and inward rectifier K⁺, background, Na⁺-K⁺ pump, Na⁺-Ca²⁺ exchanger, and sarcolemmal Ca²⁺ pump currents. A total of 16 variables are used and represent the voltage; 5 gating variables; intracellular and extracellular Na⁺, K⁺, and Ca²⁺ concentrations; and 4 variables related to calcium buffering by troponin and calmodulin in the cytoplasm.