We have developed a mathematical model of the mouse ventricular myocyte action potential (AP) from voltage-clamp data of the underlying currents and Ca2+ transients. Wherever possible, we used Markov models to represent the molecular structure and function of ion channels. The model includes detailed intracellular Ca2+ dynamics, with simulations of localized events such as sarcoplasmic Ca2+ release into a small intracellular volume bounded by the sarcolemma and sarcoplasmic reticulum. Transporter-mediated Ca2+ fluxes from the bulk cytosol are closely matched to the experimentally reported values and predict stimulation rate-dependent changes in Ca2+ transients. Our model reproduces the properties of cardiac myocytes from two different regions of the heart: the apex and the septum. The septum has a relatively prolonged AP, which reflects a relatively small contribution from the rapid transient outward K+ current in the septum. The attribution of putative molecular bases for several of the component currents enables our mouse model to be used to simulate the behavior of genetically modified transgenic mice.
- cardiac myocytes
- computer modeling
the cardiac action potential (AP) is complex, resulting from the interaction of multiple nonlinear time-dependent currents. Mathematical models of the AP have therefore been an important tool in understanding and exploring the complexities of membrane electrophysiology for more than half a century. The first cardiac models were based on the Hodgkin and Huxley equations of Na+ and K+ currents in the squid giant axon (see Ref. 70). Progress in modeling the AP has mirrored progress in the development of experimental techniques and approaches. For example, despite the considerable limitations and artifacts of the sucrose-gap method, data from these experiments led to the discovery of several new and important mechanisms, such as Ca2+ currents (79) and multiple, distinct K+ currents that are active during repolarization. These currents were incorporated into some early cardiac models (8, 66).
The pioneering model of DiFrancesco and Noble (27) was the first to incorporate active transport and secondary active transport mechanisms. These were included in response to the growing experimental evidence that these electrogenic processes were important in repolarization. Development of the whole cell patch-clamp technique gave rise to a number of animal- and region-specific models, based on detailed descriptions of ionic currents from isolated single-myocyte experiments (41, 55, 62, 63, 76, 77, 91, 98). These models began to address the basis for the diverse AP behavior observed in different species and regions of the heart. Models of bullfrog atrial and cardiac pacemaker cells (76, 77) replicated the lack of a significant intracellular Ca2+ release contribution in these cells. Similarly, models of the guinea pig ventricular myocyte (62, 63, 98) reflected the absence of the transient outward current in these cells. Lindblad et al. (55) developed a model of the rabbit atrial cell with prominent transient outward (Ito) and rapid delayed rectifier (IKr) currents. Human atrial models were developed by Nygren et al. (72) and Courtemanche et al. (25), based partially on measurements from atrial appendage tissue obtained from bypass surgery. Ito had a prominent role in these AP models, reflecting the experimentally obtained data. Recently, a rat ventricular myocyte model was published by Pandit et al. (73). The rat AP has a short AP duration (APD) and relatively dominant transient outward current.
APs of cardiac cells from different animals and different heart regions have different shapes, durations, and sets of ionic currents. These variations reflect corresponding differences in the molecular basis of repolarization and hence pharmacological responses (32, 37, 56). For example, APD at 50% of repolarization (APD50) in guinea pig and canine ventricular myocytes is ∼300 ms (41, 62, 63, 91), whereas APD50 for rat ventricular myocytes from the epicardial region is only ∼13 ms (73). The guinea pig ventricular myocyte AP has an elevated plateau (41, 62, 63, 91), whereas rat and mouse ventricular APs are relatively short, with a “triangular” shape (11, 93). The human atrial AP has a relatively sharp, short peak followed by a slow repolarization phase, which takes up to 250 ms (33, 72). In general, APD seems to be correlated to the size of the heart or the heart rate for a particular species, with smaller, faster-beating hearts having shorter APDs than larger, slower-beating hearts. The shape and underlying currents for the various APs have considerable variability. This variability is less related to heart size than it is to the species and region of the heart. Presumably, the restitution patterns, pharmacology, and arrhythmogenic potential of each waveform reflects the different molecular bases and timing of underlying currents. Understanding the timing of these currents under different situations can be critical. For example, one of the arrhythmogenic defects in human ether-à-go-go-related gene (HERG) (which encodes the rapid delayed rectifier current IKr) may be spontaneously triggered arrhythmias resulting from a complex interaction of IKr with the L-type Ca2+ current and subcellular Ca2+-handling mechanisms during repolarization (67). Clearly, understanding abnormalities in cardiac repolarization requires an understanding of the functioning of the cell as an integrated system. The many complex potential interactions of transmembrane currents and intracellular ions (particularly Ca2+) leave computer modeling as the only method currently available to synthesize and understand these multiple nonlinear interactions.
In general, most models simulate sarcolemmal currents, pumps, and exchangers and combine them with intracellular ionic homeostasis to reproduce an AP. One limitation of many older models is their description of intracellular Ca2+ dynamics (55, 62, 63, 76, 77, 98). These older models did not take into account the importance of spatial localization in Ca2+ cellular dynamics, which results in generation of relatively high-amplitude Ca2+ sparks (see, for example, Ref. 11) and can only roughly approximate the relationship between Ca2+ release and Ca2+ channel inactivation. More recent models (12, 41, 91) have begun to address this issue. Our model builds on the strengths of these previous models and was developed to reproduce a mouse myocyte with a characteristically short AP.
There is a great diversity of anatomic and electrophysiological features in myocytes from different species, e.g., transverse tubules, Ca2+-handling systems, as well as the specific properties of the major ionic currents. Wherever possible, we developed our model of the mouse ventricular myocytes with experimental data from adult mouse ventricular cardiac cells. This has been possible largely because of the increase in detailed examination of currents in various transgenic mice.
Although a majority of the most important currents and parameters are well constrained by this large body of data, there is still uncertainty in the exact kinetic parameters of some current systems. Where there was uncertainty or an absence of mouse data, we used data from less closely related species and heterologously expressed clones and adjusted them to match related mouse data.
Our mouse model is specifically designed to make use of the information from recent advances in molecular biology. The mouse is the most widely used animal in genetic research. There is now an abundance of transgenic mice with genetic defects associated with human diseases. Where possible, the component currents of the model have been assigned putative molecular bases, so the model can be used to test predictions about the consequences of transgenic experiments with relative ease. Our model mouse myocyte has a diversity of K+ channels, not all of which may be present in the same cell, depending on the region from which the cells were isolated. The model has seven distinct K+ currents: a rapid transient outward K+ current (IKto,f), a slow transient outward K+ current (IKto,s), a time-independent K+ current (IK1), an ultrarapidly activating delayed rectifier K+ current (IKur), a noninactivating steady-state K+ current (IKss), a rapid delayed rectifier K+ current (IKr), and a slow delayed rectifier K+ current (IKs).
Molecular biology has provided considerable information on the structure and function of ion channels, e.g., mechanistic information concerning ion channel gating and its modification by either genetic disease or pharmacological intervention. To be able to assess the consequences of these modifications, we used Markov models for several important currents: the Na+ current (INa), the L-type Ca2+ current (ICaL), and IKr. Many genetic manipulations and diseases alter the kinetic properties of a single ion channel and can be related to changes in specific rate constants in the Markov models for these channels. Our mouse model can therefore be used to predict the kinetic consequences of many transgenic or pharmacological modifications on channel properties.
The main goal of this paper was to develop a computer model of the mouse ventricular AP that will provide insight into the ionic basis of AP behavior. We modeled the behavior of myocytes from the apex and the septum regions of the heart (36, 37, 69, 92) to demonstrate that the different regional expression of IKto,f, IKto,s, IKur, and IKss can account for regional differences in myocyte repolarization in the mouse heart. The model contains transmembrane pumps, currents, and exchangers and a detailed intracellular Ca2+-handling system based on data from voltage-clamp experiments on individual currents (e.g., steady-state inactivation, recovery from inactivation). We developed Markov models for the fast Na+ current (INa), the L-type Ca2+ current (ICaL) and IKr. These formalisms allow for future investigation of the cellular mechanisms of arrhythmia generation due to genetic mutations or state-dependent binding that alter functional electrical properties such as restitution.
All animal procedures conformed to the “Guiding Principles for Research Involving Animals and Human Beings” of the American Physiological Society.
Preparation of Left Ventricular Myocytes
Mouse cardiac myocytes were prepared as described previously (48). Briefly, the heart was rapidly excised and submerged in Ca2+-free Tyrode solution containing (in mM) 140 NaCl, 5.4 KCl, 1 MgCl2, 0.33 NaH2PO4, 10 glucose, and 5 HEPES adjusted to pH 7.4. The aorta was cannulated with a blunt-tip needle (20 gauge) placed on a perfusion apparatus. The heart was retrogradely perfused for 3 min with Tyrode solution and then perfused for 18 min with Tyrode solution with 2% calf serum (Sigma-Aldrich) and 75 U/ml each of collagenase type I and type II (Worthington). All solutions were continuously bubbled with 95% O2-5% CO2 at 37°C. Isolated myocytes were stored at room temperature in low-Cl−, high-K+ solution (in mM: 70 KOH, 147 l-glutamic acid, 40 KCl, 20 taurine, 20 KH2PO4, 10 glucose, 10 HEPES, and 0.5 EGTA, adjusted to pH 7.3 with Tris base) before experiments. Calcium-tolerant, rod-shaped ventricular myocytes with clear striations were randomly selected for electrophysiological studies.
Cellular Electrophysiological Studies
Cell-attached patch-clamp macroscopic currents were recorded with an Axopatch 1D amplifier (Axon Instruments). A DigiData 1200 (Axon Instruments) controlled by pCLAMP software (Axon Instruments) was used to generate command pulses and acquire data. Electrodes (Garner Glass) had tip resistances of 1–4 MΩ. For Na+ current measurement, the following extracellular solution was used (in mM): 40 NaCl, 5 CsCl, 20 TEA-Cl, 2.5 MgCl2, 5 CoCl2, 5 4-aminopyridine (4-AP), 10 glucose, 80 sucrose, and 5 HEPES adjusted to pH 7.3. The Na+ current intracellular solution contained (in mM) 110 CsCl, 20 TEA-Cl, 100 aspartic acid, 5 EGTA, 2 MgCl2, 5 HEPES, 2 MgATP, and 0.1 Na3GTP adjusted to pH 7.3. For Ca2+ current measurement, the extracellular solution contained (in mM) 130 TEA-Cl, 2 CaCl2, 1 MgCl2, 10 HEPES, 5 4-AP, and 10 sucrose adjusted to pH 7.3. The intracellular solution was the same as that used for INa measurement. For K+ current measurement, the pipette was filled with a standard internal solution containing (in mM) 100 l-aspartic acid, 110 KOH, 20 KCl, 2 MgCl2, 5 EGTA, 5 HEPES, 2 Mg2ATP, and 0.3 Na3GTP, with pH adjusted to 7.2 with Tris base. Mouse ventricular myocytes were superfused at room temperature (20–22°C) with a HEPES-buffered Tyrode solution containing (in mM) 135 NaCl, 5.4 KCl, 1 MgCl, 2 CaCl2, 5 HEPES, and 10 glucose, with pH adjusted to 7.3 with NaOH. ICaL was blocked by 1 μM nifedipine, and Mg2ATP in the pipettes suppressed the ATP-sensitive K+ current.
Simulation Methods and Simulated Voltage-Clamp Protocols
The model is based on a set of 40 ordinary differential equations solved by a fourth-order Runge-Kutta method, with a time step of 0.0001 ms. A summary of the equations, model parameters, and initial conditions is given in the appendix. The model is nominally adjusted for room temperature of 25°C (298 K). Steady-state initial conditions were obtained by running the model until changes in all variables did not exceed 0.01%. Ryanodine receptor modulation factor PRyR was set to 0 at the beginning of each AP.
We used two types of simulated voltage-clamp protocols: a double-pulse steady-state inactivation protocol and a variable-interval gapped pulse protocol. Usually, the double-pulse steady-state inactivation protocol involved a pulse (P1) to a voltage between the holding potential and + 50 mV (in 10-mV intervals), followed immediately by a second pulse (P2) to a holding potential specific to the current under investigation. The durations of the two pulses were set appropriately for the gating time constants of the given current. Other protocols are described in the text and figures. For ICaL, P1 and P2 were separated by a 2-ms return to the holding potential, to match experimental protocols. APs were simulated for at least 100 cycles until a steady state was established, except for simulation of the negative staircase of Ca2+ transients, which began from a quiescent steady state.
The variable-interval gapped pulse protocol was used to determine the rate of recovery from inactivation. A P1 pulse was applied to activate and inactivate the channel, followed by a second pulse after a variable time interval. The degree of recovery from inactivation was calculated by comparing the peak P2 current to the peak P1 current. The depolarization, duration, and range of interpulse intervals used were specific to the current being studied. Unstimulated ventricular myocytes are quiescent, so we used an 0.5-ms 80 pA/pF stimulus current (Istim) applied at a frequency of 0.5–6.7 Hz to trigger simulated APs.
The electrophysiological data were analyzed by using Clampfit (Axon Instruments), Excel (Microsoft), and Origin (OriginLab). Clampfit was used for the exponential fitting of inactivation kinetics. Other nonlinear curve fittings were performed in Origin. The quality of the fit was evaluated by visual inspection and comparison of χ2-values.
We modeled cellular electrical activity in the mouse by using a nondistributed equivalent electrical circuit with subcellular compartmental spaces. This assumes that there are no electrical gradients within the cell itself, i.e., the membrane potential is spatially homogeneous and all subcellular compartments are uniform or “well stirred.” A schematic representation of the currents, fluxes, and physical compartments of the model is shown in Fig. 1. In generating APs, the model mouse membrane potential (V) was determined by the following differential equation: (1) where Cm is membrane capacitance, INa is the fast Na+ current, ICaL is the L-type Ca2+ current, IKto,f is the rapidly recovering transient outward K+ current, IKto,s is the slowly recovering transient outward K+ current, IKr is the rapid delayed rectifier K+ current (mERG), IKur is the ultrarapidly activating delayed rectifier K+ current, IKss is the noninactivating steady-state voltage-activated K+ current, IK1 is the time-independent inwardly rectifying K+ current, IKs is the slow delayed rectifier K+ current, INaCa is the Na+/Ca2+ exchange current, Ip(Ca) is the Ca2+ pump current, INaK is the Na+/K+ pump current, ICl,Ca is the Ca2+-activated Cl− current, ICab and INab are the background Ca2+ and Na+ currents, and Istim is the externally applied stimulation current. Our model has differential expression of IKto,f, IKto,s, IKur, and IKss in myocytes from the apex and the septum region of the mouse heart.
Derivation and Simulation of Component Currents
Fast Na+ current INa.
The fast transient Na+ current drives rapid depolarization at the upstroke of the ventricular AP. Most older models of INa were based on the Hodgkin and Huxley formalism of the neuronal AP, where inactivation is a single independent and completely absorbing state. INa inactivation is fast and relatively complete after depolarization, so the exact gating mechanisms of this channel were not considered in most previous AP reconstructions, despite evidence that a small noninactivated component of INa might play a role in repolarization (20, 22–24). Na+ channel inactivation has at least two components: fast and slow. Na+ channel mutations in the mouse heart that alter the slow component of inactivation can underlie a form of long QT syndrome (a frequently fatal genetic disease in other species). This suggests an important role for INa in the pathophysiology of arrhythmias (71). Consequently, more advanced Markov-type models are now beginning to be used to characterize INa during the AP (20, 22, 23, 29).
We used a Markov model for INa (22) with three closed states (CNa1, CNa2, and CNa3), an open state (ONa), a fast (IFNa) and two intermediate (I1Na and I2Na) inactivated states, and two closed inactivation states (ICNa2 and ICNa3), as shown in Fig. 2. There is a reversible fast inactivation ONa → IFNa transition, which we have attributed to the initial binding of the inactivation domain to the core domain. The second inactivation transition, IFNa → I1Na, would then represent stabilization of the inactivation domain-receptor complex. Only a few channels enter the I2Na state via slow transitions (22). A direct transition pathway between IFNa, the fast-inactivated state, and the CNa1 closed state was introduced to match inactivation recovery data (29). Inclusion of two closed-inactivation states, ICNa2 and ICNa3, allows for more accurate representation of Na+ channel availability (22). The transition rates were based on the model of Clancy and Rudy (22), with parameters adjusted to match mouse data. The transition rates and dynamic equations for INa are detailed in the appendix. INa is described by the equation: (2) where GNa is the whole cell conductance (mS/μF), ENa is the Na+ reversal potential, and ONa is the probability of the Na+ channel being in the open state.
Figure 3 shows simulated INa voltage-clamp experiments and the equivalent experimental results. We used a conventional double-pulse protocol with a 500-ms P1 pulse to potentials between −130 and +50 mV in 10-mV steps from a holding potential of −140 mV. This was immediately followed by a 180-ms P2 pulse to −20 mV. We also used the voltage-clamp protocols and ionic conditions of References 9, 51, and 86 to simulate their experimental data. Figure 3A shows simulated INa for different P1 amplitudes under control conditions. Only the first 30 ms of the P1 pulse is shown, so that activation and inactivation can be seen clearly. The simulations are comparable to typical INa experimental results (9, 38). Figure 3B shows normalized P1 peak INa as a function of voltage from our model and similar experiments (9, 51, 86). Simulations were run with extracellular Ca2+ concentration ([Ca2+]o) 1.8 mM and three different values for extracellular Na2+ concentration ([Na+]o): 10, 52.5, and 140 mM. These simulations were repeated with [Ca2+]o = 0.5 mM, to mimic the experimental conditions of References 51 and 86. With low [Ca2+]o, the activation threshold for INa decreased by 15 mV (38).
Simulated and experimental steady-state inactivation data are shown in Fig. 3C. The simulated steady-state inactivation function is in agreement with the experimental data (9, 51, 86) within the range of experimental variability. The fast component of inactivation is dominant at most potentials. We therefore fit the simulated traces with a single exponential, to mimic the analysis of available experimental data and derive an apparent time constant of inactivation. This apparent time constant is plotted against the membrane potential in Fig. 3D for both simulated and experimental data (86). Simulations were made with voltage shifts consistent with [Ca2+]o = 1.8 mM and repeated with a 15-mV decrease in activation threshold (Fig. 3D) to mimic the low-[Ca2+]o conditions used in experiments (86).
Recovery from inactivation in the Na+ channel was simulated by using a gapped-pulse protocol with steps to −20 mV with an interstimulus holding potential of either −80 mV or −100 mV. The rate of recovery from inactivation is strongly dependent on the interstimulus holding potential, so two very different interpulse intervals (1–50 ms in 1-ms steps and 10–500 ms in 10-ms steps) were used. Fitting simulated recovery curves with a single exponential function showed that the time constant of recovery increased by a factor of ∼8 when the holding potential was reduced from −80 mV (time constant for recovery from inactivation τrec,Na = 71.5 ms) to −100 mV (τrec,Na = 8.6 ms). This compares well with our data from mouse myocytes (τrec,Na = 66.0 ± 1.9 at −80 mV; n = 5).
One of the advantages of using a Markov model is that transition rates can be modified to simulate the effects of mutations on gating. However, our model of INa gating has several important limitations. The most important of these stems from the size of INa and the consequential difficulty in measuring this current. INa is large and can cause resistive drops of 10 mV or more across the microelectrode tip, resulting in loss of voltage clamp control. Most investigators, therefore, partially block the channel or otherwise reduce conductance to maintain voltage clamp. These manipulations (e.g., reducing [Na+]o or applying partial doses of blocker) are known to alter gating kinetics and equilibrium. Consequently, we also validated our model parameters for INa by comparing the calculated and experimental rate of maximal upstroke velocity (dV/dtmax) due to INa AP upstroke (see Mouse Action Potential below).
L-type Ca2+ current ICaL.
Perturbations in the amplitude and gating kinetics of ICaL have been associated with many arrhythmias and pathological states such as failing, hypertrophic, stunned, ischemic, and hibernating myocardium (47). Activation of the L-type Ca2+ channel is a purely voltage-dependent process, similar to Na+ and K+ channel activation. The molecular basis of L-type Ca2+ channel inactivation is not fully understood. Inactivation of the L-type Ca2+ channel is both a time-dependent (i.e., apparently voltage dependent because of its coupling to activation) and a Ca2+-dependent process (43, 52). Ca2+-induced inactivation of the Ca2+ channel appears to be modulated by calmodulin (74), which may bind to the Ca2+-binding motif (EF hand) on the carboxy tail of the main α1C-subunit (53), thus transducing calmodulin binding into channel inactivation (75).
Describing the complex inactivation behavior of the L-type Ca2+ channel requires assumptions about coupling and inactivation biophysics. Jafri et al. (41) used a mode-switching Markov model for Ca2+ inactivation of the L-type Ca2+ channel to link inactivation to subspace [Ca2+]. The channel is modeled with four independent subunits that can close the channel and whose mode of operation is determined by intracellular Ca2+ binding. Transitions between the two modes are controlled by a Ca2+-dependent factor, γ. The probability of entering the Ca2+ mode is dependent on voltage, conformation, and subspace [Ca2+]. This model makes little use of the emerging structure and function information from voltage-gated channels and fails to model recovery from inactivation.
We developed a new Markov model to describe L-type Ca2+ channel dynamics (Fig. 4), based on homology and extrapolated structure-function relationships of Shaker B voltage-gated K+ channels. Our model has two inactivation states, one of which is a Markov representation of “ball and chain” inactivation in which the ball binding is Ca2+ dependent (74, 75). The model Ca2+ channel has four closed states (C1, C2, C3, and C4), one open state (O), and three inactivated states (I1, I2, and I3). Activation is controlled by the voltage-dependent rate constants, α and β. Transitions from O to I1 correspond to a relatively fast Ca2+-dependent inactivation, controlled by the Ca2+-dependent rate constant γ. Transitions from O to I2 represent relatively slow voltage-dependent inactivation and are controlled by the rate constants α and Kpcf (forward voltage-insensitive rate constant). The I3 inactivated state represents a channel with both Ca2+- and voltage-dependent inactivation. Direct transitions between the closed state C4 and the inactivated states I1, I2, and I3 result in a voltage-dependent recovery from inactivation. All rate constants fulfill the condition of thermodynamic reversibility (for review see Ref. 38), i.e., the product of the rate constants in any clockwise loop is equal to that of the same loop measured in a counterclockwise direction.
The ICaL equation has the form: (3) where GCaL is the whole cell conductance (mS/μF), ECa,L is the L-type Ca2+ channel reversal potential, and O is the channel open probability.
No matter how well the Ca2+ channel is modeled, it cannot be considered in isolation from its surroundings. In ventricular myocytes, L-type Ca2+ channels are found mostly in clusters in the T-tubules directly opposed to the ryanodine receptors on the sarcoplasmic reticulum (SR) (13). Ca2+ that enters the cell as part of ICaL initiates Ca2+-induced Ca2+ release from the SR, which then contributes to Ca2+-dependent inactivation of the Ca2+ channel (10, 30, 31). A more physiologically realistic model of the L-type Ca2+ channel must therefore include an appropriate representation of the SR release channels and the local [Ca2+] near the L-type Ca2+ channel.
Ca2+ fluxes in and around the SR are shown in Fig. 1. Ca2+ entering the cell via through ICaL results in a localized increase in [Ca2+] ([Ca2+]ss) in a restricted subsarcolemmal space. [Ca2+]ss both inactivates the L-type Ca2+ channel and initiates Ca2+-induced Ca2+ release, which further adds to Ca2+ channel inactivation. Ca2+ diffuses from the subsarcolemmal space to the general cytosol, where it binds to contractile proteins and initiates contraction. Ca2+ is removed from the cytosol by translocation across the sarcolemmal membrane by the Na+/Ca2+ exchanger and the sarcolemmal Ca2+ pump and taken back into the network SR by Ca2+-ATPase. Ca2+ diffuses from the network SR to the junctional SR, where it can be released into the subsarcolemmal space once more. The equations governing the time dependence of the ryanodine receptor are based on those of Keizer and Levine (45). We also modified the Jafri et al. (41) model of Ca2+ dynamics to match Ca2+ transients reported for the mouse. A Gaussian distribution function and integral of the time course of ICaL were used to modulate Ca2+ release and simulate graded release of Ca2+ from the intracellular store (Eq. A15, appendix). Incorporation of local effects of Ca2+ release in our mouse model allowed simulation of Ca2+ sparks, the elementary units of Ca2+ signaling, and a precise description of Ca2+ handling and Ca2+ channel inactivation that resembles that of more complex biophysically based models (12).
Simulated ICaL voltage-clamp studies are shown in Fig. 5. A conventional double-pulse protocol was used with a 250-ms P1 pulse stepping from −80 to between −70 and +40 mV followed by a 2-ms return to −80 mV and then a 250-ms P2 pulse to +10 mV. The biexponential current decay is typical for ICaL from mouse myocytes (see, for example, Ref. 88). The half-time of simulated ICaL decay at +10 mV was 4.2 ms, which is smaller than the corresponding experimental value of 8 ms obtained by Wang et al. (90). However, the calculated fast inactivation time constant τ1,CaL = 4.45 ms compares well to the experimental value of 4.7 ± 0.2 ms obtained by Kirchhefer et al. (49), but is somewhat smaller than 12.28 ± 1.29 ms (82) and 7.92 ± 0.74 ms (81) obtained by Santana et al. The simulated slow inactivation time constant τ2,CaL = 16.1 ms is close to the value of 22.47 ± 1.89 ms obtained by Santana et al. (82); however, it is somewhat smaller than the measured values of 56.1 ± 2.6 ms observed by Kirchhefer et al. (49) and 48 ± 1.35 ms recorded by Santana et al. (81). The variability in time constants for the components of inactivation between laboratories probably reflect different experimental conditions (e.g., amount of EGTA, effectiveness of perfusion through electrode tip, isolation procedures, ionic conditions) that can effect both Ca2+- and voltage-dependent inactivation. The simulated traces in Fig. 5A mimic control conditions without EGTA or other nonintrinsic buffers and so are closest to the fastest time constants, which we assumed had the least perturbation from altered Ca2+ handling. When Ca2+-release flux during simulation is decreased by 20-fold to mimic the effect of strong buffering by EGTA, the rate of decay of ICaL is reduced, with fast and slow time constants equal to 8.8 and 57.3 ms, respectively (V = +10 mV).
Normalized peak P1 ICaL as a function of P1 voltage is shown in Fig. 5B. Our model data for normal and high-BAPTA conditions compare well with corresponding experimental data (96). The maximum value of simulated ICaL is 6.9 pA/pF, which is similar to our experimental mouse data (7.7 ± 0.6 pA/pF; n = 17) and within the experimentally recorded range of 5.1–11.6 pA/pF (42, 80, 88). Figure 5C shows the simulated and experimental steady-state inactivation function for ICaL. Figure 5D shows the voltage dependence of the simulated time constants of inactivation.
Figure 6A shows the bulk intracellular [Ca2+] ([Ca2+]i) and subspace volume [Ca2+]ss transients in response to a step depolarization to 0 mV and corresponding experimental data (42). The simulated [Ca2+]i time to peak of 14.8 ms and half-time of decay of 33.6 ms compare well with the experimental values of 17 ms and 34 ms, respectively (90).
Ca2+ spark duration estimated from experiments is ∼10 ms (14), which is comparable to the 9.4-ms simulated time constant for decay of [Ca2+]ss. The bulk [Ca2+]i transient is slower than the [Ca2+]ss transient in both experiment and model, with a tail of hundreds of milliseconds (14, 42). Figure 6B shows the bell-shaped voltage dependence of the simulated [Ca2+]i, which is in qualitative and quantitative agreement with experimental observations (40, 42, 80, 83). The peak ICaL has a similar voltage dependence and shows the dependence of [Ca2+]i on ICaL, mediated by graded Ca2+ release from the SR.
Figure 7 shows the results of a simulated variable-gap double-pulse protocol. Two 250-ms pulses to 0 mV were applied with an interstimulus interval of between 2 and 500 ms for three different interpulse holding potentials. Recovery from inactivation was approximated by a single exponential function with time constants τrec,Ca equal to 19, 39, and 77 ms for three holding potentials, −90, −80, and −70 mV, respectively. These simulations suggest a moderate voltage dependence of recovery. However, no data exist on the voltage-dependent recovery of ICaL in mouse myocytes. The L-type Ca2+ channel in bullfrog sinus venosus and atrium (17) and in rat and rabbit ventricles (97) shows voltage dependence in this range, with recovery becoming faster with hyperpolarization. Voltage-dependent recovery is also consistent with the general properties of recovery from other voltage-gated channels (see Ref. 38 for a review).
Rapidly inactivating transient outward K+ current IKto,f.
There are several rapidly activating voltage-gated K+ currents with differing kinetics in mouse ventricular cells. The most prominent of these is a rapidly activating and inactivating transient outward K+ current IKto,f. The molecular bases of this current are the Kv4.2–Kv4.3 family of channels (60, 92). The kinetic description of this current was derived from the native Ito model for the ferret right ventricle (58) with parameters adjusted to match experimental mouse data (87, 92, 100). The model formulation of IKto,f is: (4) where GKto,f is the maximum whole cell conductance (mS/μF), EK is the K+ reversal potential, and ato,f and ito,f are the activation and inactivation gating variables. Equations for ato,f and ito,f are given in the appendix. This current is equivalent to Ito,f in mouse ventricular myocytes (92).
Mouse ventricular myocytes show substantial regional heterogeneity in repolarization characteristics and can be divided into several types depending on both the anatomic location and the number and type of K+ currents present (37, 92). We modeled ventricular myocytes from two regions: the apex and the septum. The fundamental differences between these cells are the APD and density of K+ currents present (37, 92).
Figure 8 shows the combined depolarization-activated simulated K+ currents (IK,sum) from the apex and septum when depolarized from a holding potential of −80 mV to between −70 and +50 mV in 10-mV steps. Both simulations are in qualitative and quantitative agreement with the experimental results (92). In our model, there are seven types of K+ currents: the rapid transient outward K+ current (IKto,f), the slow transient outward K+ current (IKto,s), the rapid delayed rectifier K+ current (IKr), the ultrarapidly activating delayed rectifier current (IKur), the noninactivating steady-state K+ current (IKss), a very small slow delayed rectifier (IKs), and the time-independent K+ current (IK1).
Figure 9 shows a simulation of apical IKto,f in response to a double-pulse protocol (similar simulations from the septum are not shown). Two 500-ms pulses were applied, one from the holding potential (−80 mV) to between −100 and +50 mV in 10-mV intervals and the second to +50 mV. The peak current-voltage relationships in Fig. 9B are from simulated and experimental data for apical cells (92). In our simulations, IKto,f corresponds to the fast component of Ito in other experimental papers (87, 92). Figure 9C shows the simulated and experimental steady-state inactivation relationships for IKto,f (87, 92). Figure 9D shows simulated and experimental time constants for activation and inactivation as functions of voltage.
In the range of potentials where activation is complete, inactivation of IKto,f shows little or no voltage dependence. This is a common finding for K+ channels that inactivate by either an N-type or a C-type inactivation mechanism (78). However, in the range where activation is not complete, inactivation is generally voltage sensitive (for a review see Ref. 38). This does not imply that inactivation becomes intrinsically voltage dependent in this range, rather that coupling of inactivation to activation must result in a slowing of the rate of inactivation because of a substantial fraction of channels being in a subthreshold state for development of inactivation. Positive to 0 mV, both model and experimental data have voltage-insensitive inactivation. Our model shows voltage dependence below ∼0 mV, which is somewhat inconsistent with the data for mice, as shown in Fig. 9. However, below ∼0 mV, currents are relatively small and inactivation is difficult to measure directly. In the range in which only little or no current is activated, experimenters studying transient outward currents in other preparations have observed voltage-dependent inactivation (18), as we have modeled here. Nonetheless, both Hodgkin-Huxley and Markov models are constrained to have voltage-dependent inactivation in the range in which inactivation is less than complete. The unusual behavior of the native channel may be due to incompletely understood processes related to multiple inactivated states.
A simulated gapped-pulse protocol was used to determine the rate of recovery of IKto,f from inactivation. Two pulses to +50 mV from the holding potential of −70 mV were applied with an interstimulus interval of 10–500 ms. The simulated recovery time constant was τrec,Kto,f = 27.4 ms, which is close to the experimental value of 27 ms recorded in the mouse (92). The time constants of recovery are identical for IKto,f from the apex and the septum (92).
Slow-inactivating transient outward K+ current IKto,s.
The slow-inactivating K+ current, IKto,s, presumably encoded by Kv1.4, is present in cardiac myocytes from the septum region, but largely absent from other regions of the heart (92). The slowly recovering inactivating K+ current might be considered a compensatory current, because it is present in septal cells where there is a considerably smaller IKto,f and in ventricular myocytes from mice in which IKto,f has been transgenically deleted (7, 36, 69). The model formulation of IKto,s is: (5) where GKto,s is the maximum whole cell conductance (mS/μF) and ato,s and ito,s are the activation and inactivation gating variables.
Figure 10A shows simulated current traces in response to a series of step depolarizations from −100 to between −90 and +50 mV in 10-mV increments. The inactivation rate of IKto,s is voltage independent, with a time constant of 270 ms that corresponds to the experimental value of 258 ± 15 ms (92). The simulated and experimental peak current-voltage relationships are shown in Fig. 10B. Figure 10C shows the simulated and experimental steady-state activation and inactivation functions of IKto,s. Our model was developed to simulate experimental currents recorded in the absence of Ca2+ channel blockers (16). Data obtained in the presence of divalent ionic blockers of the Ca2+ channel, such as 2 mM CoCl2, showed a positive shift in activation and inactivation due to the surface charge effect (99). We therefore adjusted our model gating parameters to fit data and account for the shift in activation and inactivation. Figure 10D shows the voltage dependence of time constants of activation for IKto,s. The simulated time constant for recovery from inactivation is 1.306 s, which compares well with the experimental value of 1.298 s (92).
Rapid delayed rectifier K+ current IKr.
The kinetics of IKr are complex and cannot be adequately described by a simple Hodgkin and Huxley formalism (58). Because a detailed kinetic characterization of IKr has not been accomplished in mouse myocytes, we modeled IKr (HERG in humans), using a variation of the Markov model of Wang et al. (89): (6) IKr is described by the following equation: (7) where GKr is the whole cell conductance (mS/μF), OK is the probability of the channel being in the open state, R is the ideal gas constant, F is the Faraday constant, and [K+]o and [K+]i are the K+ concentrations outside and inside the cell, respectively. The driving force for permeation through IKr has some permeability to Na+ to account for the experimentally observed reversal potential in ventricular myocytes (∼10 mV positive to EK). Simulated IKr in response to a double-pulse voltage-clamp protocol is shown in Fig. 11. P1 was a 1-s pulse from the holding potential (−80 mV) to between −70 and +50 mV in 10-mV steps, and P2 was a second 1-s pulse to −40 mV. The simulated currents are comparable to experimental results from IKr (57). The calculated value of maximum IKr amplitude of 0.25 pA/pF was chosen to be close to experimentally observed value in adult mice (57). Our representation of IKr as a linear Markov model does not permit inactivation from closed states. Single-channel experiments suggest that these transitions can occur (46), and they have been included in other IKr models (21). Kiehn et al. (46) described a relatively rapid flicker state, which allows the channel to inactivate without opening when inactivation is sufficiently fast and complete at positive potentials. There is no experimental evidence that the flicker state is part of the activation process or that it has any consequences for drug binding or burst duration. For macroscopic currents, this flicker is readily addressed as a reduced conductance with a mean lifetime equal to burst duration. We have treated it as such in this model, so as not to complicate the model with a rapid time constant.
Ultrarapidly activating delayed rectifier K+ current IKur and noninactivating steady-state K+ current IKss.
In our model, IKur and IKss are described by equations: (8) (9) where GKur and GKss are the maximum whole cell conductances (mS/μF), aur and aKss are activation gates, and iur and iKss are inactivation gates, described in detail in the appendix. IKur and IKss correspond to IK,slow and Iss, respectively (92). We chose to use the IKur terminology to emphasize the rapid activation, as opposed to the slow inactivation that is characteristic of these channels. Use of this terminology also emphasizes the similarity of this current to IKur in human atrium (for review, see Ref. 68). These currents are present in both the apex and the septum, but with different expression levels. The model parameters were chosen to fit mouse experimental data (92).
Figure 12 shows simulated currents in response to a series of step depolarizations from −100 to between −90 and +50 mV in 10-mV intervals. The inactivation rate of IKur is voltage independent and can be fit with a single exponential function (92). The simulated time constant of inactivation is 1,198 ms, which compares well with the experimental value of 1,180 ± 45 ms (92). IKss shows little inactivation in either simulation or experiment (92).
The value of the maximum whole cell conductances, GKur and GKss, were based on voltage-clamp step depolarizations from −70 to +40 mV in myocytes from both the apex and the septum (92). Figure 13, A and B, show the simulated and experimental peak current-voltage relationships, respectively, for IKur and IKss from the apex and septum (92). The simulated and experimental data are in good qualitative and quantitative agreement. Figure 13C shows the simulated and experimental steady-state inactivation functions of IKur (IKss does not inactivate). Our simulations in Fig. 13C mimic experimental data recorded without divalent Ca2+ channel blockers (16). The data in Fig. 13B obtained in the presence of 5 mM CoCl2 showed a positive shift in activation due to the effect of the divalent ions (92). We therefore adjusted our model gating parameters to fit data to account for the shift in activation. Figure 13D shows the voltage dependence of time constants of activation for IKur and IKss. The simulated time constant for recovery from inactivation of IKur (1.032 s) compares well with the experimental value of 1.079 s (92).
Slow delayed-rectifier K+ current IKs.
Expression of IKs in mouse heart is questionable. It may be present in some strains (6, 28) but not in others (19, 65). There is growing evidence that IKs knockout mice are more susceptible to arrhythmias. This indirect evidence suggests an important role for IKs in AP generation (5, 6, 28). Where IKs is detected, it is present in fewer than 10% of myocytes (6). The IKs in our model has a relatively small maximum conductance and slow activation, which is similar to the experimental data (28). We used the formulation of Rasmusson et al. (76) with minor modifications: (10) The equations governing the activation gate, nKs, are shown in the appendix. Simulated voltage-clamp data for IKs are shown in Fig. 14. IKs does not significantly affect AP shape under control conditions; however, it can be of potential importance in experiments where other K+ currents are knocked out.
Ca2+-activated Cl− current ICl,Ca.
ICl,Ca is determined by the equation: (11) where GCl,Ca is the maximum whole cell conductance (mS/μF), OCl,Ca is the probability of the channel being in the open state, Km,Cl is the half-saturation constant, and ECl is the Cl− reversal potential. This small-amplitude current was recently discovered by Xu et al. (94). Experiments (94) and simulations (not shown) suggest that the current does not affect the AP shape under control conditions. However, the current may be important under pathophysiological conditions that result in abnormal cellular Ca2+ loading.
Time-independent K+ current IK1.
The IK1 equation in our model is based on that of DiFrancesco and Noble (27): (12) Figure 15 shows simulated and experimental (7, 50, 88) peak current-voltage relationships in response to a pulse from the holding potential of −80 mV to voltages from −150 to −40 mV. Simulated and experimental results give currents of similar magnitude.
Background currents ICab and INab.
ICab and INab are formulated as linear ohmic currents: (13) (14) where GCab is the background Ca2+ conductance, GNab is the background Na+ conductance, ECaN is the Ca2+ reversal potential, and ENa is the Na+ reversal potential. These currents normally maintain ionic homeostasis.
Na+/Ca2+ exchange current INaCa.
The equation representing INaCa is based on that of Luo and Rudy (62): (15) where Km,Na is the Na+ half-saturation constant for INaCa, Km,Ca is the Ca2+ half-saturation constant, ksat is the saturation factor at very negative potentials, and η = 0.35 is the position of the energy barrier that controls the voltage dependence of INaCa. INaCa plays an important role in mouse myocyte Ca2+ dynamics (84). The absolute value of INaCa does not exceed 1.0 pA/pF within the range −120 to +40 mV. The scaling factor kNaCa was chosen to ensure equilibrium intracellular ionic concentrations of Na+ and Ca2+ within experimental values and to match the experimental ratio of Ca2+ extruded by the SR ATPase to Ca2+ expelled from the cell via INaCa during a myocyte twitch.
Sarcolemmal Ca2+ pump current Ip(Ca).
Ip(Ca) is taken from the Luo and Rudy (62) model: (16) where Ip(Ca)max is the maximum Ca2+ pump current and Km,p(Ca) is the Ca2+ half-saturation constant.
Na+/K+ pump current INaK.
INaK was modeled by using the Luo and Rudy (62) formulation: (17) (18) (19) where INaKmax is the maximum pump current, Km,Nai is the Na+ half-saturation constant, and Km,Ko is the K+ half-saturation constant. This current maintains Na+ and K+ electrochemical gradients across the cell membrane. INaKmax was optimized to reproduce experimental values of quiescent [Na+]i and [K+]i for mouse myocytes.
Quiescent Cellular Properties
To calculate the steady-state resting cellular parameters, the mouse model was allowed to run for at least 1,000 s with no external stimulation. The currents that determine the resting membrane potential are IK1, INaCa, Ip(Ca), INaK, and the background currents ICab and INab. The reversal potential of the time-independent K+ current is the major determinant of the resting potential of the quiescent model myocyte. The resting potential of the cell is only 1.8 mV positive to the reversal potential for IK1. INaK pumps Na+ out of and K+ into the cell, thus maintaining [Na+]i and [K+]i. At rest, INaCa extrudes Ca2+ and brings Na+ into the myocyte, contributing to equilibration of intracellular ion concentrations. The sarcolemmal Ca2+ pump carries Ca2+ out of cell. The magnitudes of these currents were chosen so that the resting potential and equilibrium intracellular ion concentrations were within experimental values (see Table 8 in appendix).
Mouse Action Potential
The mouse APD is region dependent (37), with cells from the septum in general having longer APDs than those from the apex (37). Figure 16 , A and D, shows simulated APs for these two regions obtained with a 1-Hz stimulation rate. APs were obtained after at least 100 s of continuous simulation to ensure that a steady-state solution has been reached. In the apex, the experimental value of APD50 is 4.5 ± 0.3 ms, which is somewhat longer than the simulated value of 4.0 ms. In the septum, the experimental APD50 is 8.5 ± 0.3 ms and the simulated APD50 is 7.6 ms. Regional differences in AP shape are also evident at other potentials. There is close agreement at APD at 25% of repolarization (APD25), with simulated values of 2.0 and 3.4 ms for apex and septum, respectively, compared with the experimentally observed values of 2.3 ± 0.1 and 4.0 ± 0.2 ms (37). There is also close agreement at APD at 75% of repolarization (APD75), with simulated values of 11.3 and 16.1 ms for apex and septum, respectively, compared with the experimental values of 10.0 ± 1.2 and 20.9 ± 1.3 ms (37). These data are presented in Table 1.
Figure 16, B, C, E, and F, shows the major currents underlying the AP. The largest current is INa, which peaks at −150 pA/pF in both the apex and septum. This rapid inward current is responsible for normal initiation of the cardiac AP. The maximum upstroke velocity (dV/dtmax) of the simulated AP for both the apex and septum is 156 V/s, which compares well with the values of 154 ± 10 (4) and 161 ± 44 (85) V/s observed experimentally.
Immediately after depolarization, the membrane potential falls to approximately −20 mV and the shape of the AP is determined by competition between several currents. The major contributors in apical myocytes are IKto,f, IKur, IKss, and ICaL. In septal cells, IKto,s also contributes to repolarization. The outward K+ currents slow the rate of depolarization and initiate repolarization, opposing and dominating the depolarizing inward current, ICaL. IKto,f is relatively large in the apex, where the membrane potential decreases rapidly, resulting in an APD50 of ∼4.0 ms. In the septum, IKto,f is relatively small and comparable to IKur and IKss. This results in the longer septal AP (APD50 ∼7.6 ms). This relatively prolonged depolarization affects other voltage-dependent currents, which in turn affect cellular Ca2+ dynamics and excitation-contraction coupling. The maintained depolarization decreases the amplitude but increases the duration of ICaL. It also increases the duration of Ca2+ influx via INaCa (Fig. 16, C and F).
The final stage of repolarization is relatively slow and is controlled by the slower K+ currents and the Na+/Ca2+ exchanger (IKto,s, IKur, IKss, IK1, and INaCa). IKur and IKss are slightly larger in the apex than in the septum. However, the septum has an additional K+ current, IKto,s (4.2 pA/pF), and a relatively small IKto,f (6.8 pA/pF) compared with the apex (25.9 pA/pF). These differences in K+ current magnitudes result in differing repolarization rates and APDs. The resting potential is maintained by the balance between inward (ICab, INab, and INaCa) and outward [Ip(Ca), IK1, and INaK] currents (Fig. 16, C and F). Our model suggests that the regional difference in shape and duration of the mouse cardiac AP is a result of the relative magnitudes of IKto,f and IKto,s. IKto,f dominates the rapid repolarization of the apex where IKto,s is absent, whereas the smaller IKto,f and additional IKto,s mediate the slower septal repolarization.
Ca2+ Dynamics During AP
Figure 17 shows simulated APs from the apex and septum for different pacing periods, ranging from 150 to 2,000 ms. Simulated data were obtained after 100 s of stimulation. There is little change in AP shape with stimulation frequency, which is consistent with experimental data (36). Rat epicardial myocytes show a similar rate insensitivity (73). The simulated time course of decay of [Ca2+]i after an action potential is 183 and 165 ms for the apex and the septum, respectively. These values are well within the range of experimental observations of 126 ± 11 (34), 188 ± 14 (54), 190 (39), 184 ± 22 (95), ∼210 (15), and 293 ± 19 (64) ms.
In contrast to the consistency of APD, the amplitude of [Ca2+]i does change with pacing frequency. Simulated and experimental data (2, 40) relating basal and peak [Ca2+]i to stimulation frequency are shown in Fig. 18. In our model, diastolic [Ca2+]i increases with increasing pacing frequency, which is consistent with experimentally obtained data (2, 40). Systolic [Ca2+]i shows a minimum at a pacing frequency of 1 Hz and increases at higher or lower frequencies (40). Our model reproduces the behavior of systolic [Ca2+]i, with a minimum value between pacing frequencies of 1 and 2 Hz. The shorter APDs in apical myocytes produce somewhat smaller [Ca2+]i transients. The presence of a minimum in the frequency dependence of [Ca2+]i at low pacing frequencies was observed in some experiments (35, 40), but other studies reported only weak monotonic increases in systolic [Ca2+]i with frequency (2, 34).
Our model of the mouse ventricular myocyte is capable of reproducing one of the key observations in mouse myocytes, the negative staircase of electrically evoked [Ca2+]i transients (39). Figure 19A shows experimental records of [Ca2+]i transients evoked with a stimulation period of ∼1.3 s. The peaks of the transients decline monotonically with time. A similar negative staircase is obtained from the model in both apical and septal cells (Fig. 19, C and D).
Determining Ca2+ influx and efflux, and Ca2+ cycling in general, is central to understanding the regulation of contractile force in cardiac muscle (10, 30, 31). An interesting result from our simulations is a detailed description of Ca2+ dynamics during AP (twitch) that closely resemble experimental data. Figure 20 shows the behavior of major Ca2+ fluxes during apical and septal APs and the integrated Ca2+ fluxes associated with the major Ca2+ transport pathways. Ca2+ release from the SR is initiated by a relatively large but brief Ca2+ influx through L-type Ca2+ channel (JCaL) that peaks at ∼300 and 250 μM/s for apex and septum cells, respectively. This value is close to the experimental estimates of 300 μM/s (11). The simulated Ca2+-release fluxes (Jrel) are larger and longer than the fluxes through L-type Ca2+ channel, with a maximum value of 3.0 and 3.7 mM/s for the apex and the septum, respectively (Fig. 20, A and C). These magnitudes are similar to the experimental observations from the rat, ∼2.6–3.0 mM/s (61), which has similar Ca2+ cycling (11). No experimental data are available for the mouse. Ca2+ extrusion from the cytosol is dominated by the SR pump (Jup) and the Na+/Ca2+ exchanger (JNaCa). There is also a substantial Ca2+ leak flux from the SR (Jleak). Ca2+ uptake to the SR, calculated as the difference Jup − Jleak, achieves a maximum value of 200 and 260 μM/s for the apex and septum, respectively. These values are close to the experimental estimation for the mouse of 338 μM/s (54). The maximum efflux through the Na+/Ca2+ exchanger in the model is ∼26 μM/s for both septal and apical myocytes. This compares well with the experimental data, 30.5 μM/s for the mouse (54).
Figure 20 shows the integrated Ca2+ flux for several Ca2+-handling systems. This demonstrates the relative amount of Ca2+ passing through each system during a typical cycle in the apex (Fig. 20B) and the septum (Fig. 20D). Only ∼1.0 (apex) and 1.4 (septum) μM Ca2+ enters the subspace volume of the cytosol to trigger 33 (apex) and 40 (septum) μM Ca2+ release from the SR. The latter values compare well with the corresponding experimental data of total cytosolic [Ca2+] (Δ[Ca2+]Tot) for mouse measured as 39 ± 3 μM (64). About 90% of Ca2+ that enters the cytosol (from intracellular and extracellular sources) during a simulated twitch is taken up by the SR pump (31.5 and 38.0 μM for the apex and septum, respectively). Only ∼10% of the Ca2+ is extruded by the Na+/Ca2+ exchanger (3.6 and 3.9 μM for the apex and septum, respectively). These proportions compare well with the experimental data of Li et al. (Ref. 54; 90.3% and 9.2% through SR pump and Na+/Ca2+ exchanger, respectively) and Brittsan et al. (Ref. 15; 88% and 12% through SR pump and Na+/Ca2+ exchanger, respectively).
Development of Model
The aim of this study was to develop a model of the isolated mouse ventricular myocyte that reproduces the cellular AP. The model is based on available experimental voltage-clamp data from various membrane currents present in the mouse heart. Wherever possible, the currents have been related to structure and function data, including regional heterogeneity. We tested the model currents against experimental voltage-clamp data and used molecularly based Markov models for the fast Na+ (INa), L-type Ca2+ (ICaL), and rapid delayed-rectifier K+ (IKr) channels. Simulations of AP upstroke dVm/dt, decline of [Ca2+]i transients during twitch, or APD25, APD50, and APD75 served as additional validation of the model. We included a variety of optional K+ channel types to account for the experimentally observed regional diversity in AP shape and duration.
The simulated AP was relatively short, which is characteristic of the mouse myocyte. The kinetic behavior of the currents present in our mouse model are similar to those found in other mammalian cardiac myocytes with much longer APDs, consistent with the molecular similarity of these currents. The main factor in generating the short mouse AP is not differing channel kinetics but the relative magnitude of the currents present. Mouse ventricular myocytes have large outward K+ currents, which ensure rapid repolarization. An additional factor, which has been noted previously (11), is that mice, like other rodents, have extremely fast Ca2+-handling mechanisms. This means that they require significantly less Ca2+ influx via sarcolemmal mechanisms to trigger Ca2+ release from the intracellular stores. The released Ca2+ leads to a rise in systolic Ca2+ and ultimately triggers contraction. In mammalian models with longer APDs (25, 55, 62, 63, 72, 91), a larger percentage of the systolic Ca2+ transient is due to transsarcolemmal sources rather than to Ca2+ release from intracellular stores. In summary, it appears that multiple ionic mechanisms are adapted to the rapid beating required of a mouse heart.
A model of the mouse ventricular myocyte is an important tool in understanding the many type of transgenic mice that have been created to mimic various aspects of abnormal electrical behavior in the heart (7, 42, 48, 60, 69, 93, 100). Gene knockout experiments on ion channels almost invariably lead to compensatory changes and upregulation of other ion channels (7, 37, 69). For example, an α-myosin mutation (associated with human atrial arrhythmias) alters ionic homeostasis when expressed in the mouse ventricle. This illustrates the limitations of the reductionist approach to understanding biological function. A primary defect may be localized to only a few atoms on a single molecule. However, the final phenotype resulting from the mutation is a consequence of both the defect and the response of the organism, which may result in seemingly unrelated pathologies. This pattern of adaptation of multiple subcellular systems to fit the physiological requirements for function carries through to the cellular responses to transgenic manipulation.
The lack of complete knowledge about cellular parameters may be a limitation in some respects, but dealing with this uncertainty has also been one of the strengths of systems biology in general and the study of repolarization in particular. In cardiac muscle, models of repolarization have long played an important role. The models have been developed in parallel with or even in advance of experimental validation. The earliest model of cardiac repolarization was made by postulating the presence of a diminished outward K+ current with slow kinetics, coupled with a sustained inward current (70). Although the exact nature of the K+ currents has changed and the inward current is generally thought to be carried primarily by Ca2+ instead of Na+, the general concepts spurred further investigation. Similarly, subsequent models, such as those of McAllister, Noble, and Tsien (66), Beeler and Reuter (8), DiFrancesco and Noble (27), and Luo and Rudy (62, 63, 98), all had severe limitations in parameters and descriptions of component currents but provided an essential framework in which the behavior of the cell and repolarization as a system could be quantitatively and qualitatively discussed and potential mechanisms could be evaluated and explained.
Our model, and the formulations within it, reflects advances in our understanding of how the mouse ventricular myocyte functions during excitation and repolarization. Equally importantly, the formulations within this model have been developed to parallel advances in our ability to manipulate the cellular system. Thus the main currents and transporters have been associated with specific molecular bases. Although this renders the model more complex in terms of the total number of currents included, it allows the model to be used to interpret data from molecular biology and transgenic experiments.
Incorporation of structure and function information from molecular biophysics is another important component of this cellular model. For many channels, we have moved away from the Hodgkin and Huxley type of formalism in favor of Markov-type models. Markov models have the advantage of being more closely related to the underlying structure and conformation of the ion channel proteins. Clancy and Rudy (20) demonstrated that this is extremely important for relating genetic defects that result in altered channel function to the pathological phenotype. In other cases, e.g., HERG or the L-type Ca2+ channel, conventional Hodgkin and Huxley formalisms simply fail to reproduce channel behavior (58). However, models of individual ion channels and subcellular systems are, by themselves, areas of intense research interest. For example, there is considerable debate as to how the Ca2+ channel and its coupling to the Ca2+-release subsystem should be modeled. We had to make choices between alternative formulations, some of which are clearly controversial. Some of these choices and their relationship to other models and mechanisms of repolarization in other cardiac myocytes are discussed in Mechanisms, Formalisms, and Limitations.
Mechanisms, Formalisms, and Limitations
Different species have APs with differing AP shapes and durations, with various mechanisms of repolarization. The molecular basis for AP initiation in atrial and ventricular myocytes of working myocardium is consistently the fast Na+ current INa (8, 41, 62, 63, 66, 70, 72, 73, 76). Its role in initiating depolarization and sustaining conduction is similar to that originally described for the squid axon AP. On a molecular level, cardiac INa is mediated by a cardiac-specific α-subunit isoform that is relatively insensitive to the neurotoxin TTX. Many of the kinetic properties of INa are dependent on coassembly of the α-subunit with ancillary subunits. These subunits strongly modify the inactivation and recovery properties of the α-subunit of the channel. INa has two types of inactivation with differing molecular bases, termed “fast” and “slow.” Fast inactivation is strongly dependent on a small “flap” of cytoplasmic-facing lipophilic amino acids, and slow inactivation is related to other regions. These two mechanisms have been suggested to be similar to the N- and C-type inactivation mechanisms studied in the distantly related voltage-gated K+ channels (78).
Although early cardiac models hypothesized a major role for a noninactivating component of INa during repolarization, later models considered INa to inactivate completely and contribute little to the plateau phase of repolarization. However, the discovery of a human Na+ channel defect that alters the slow component of inactivation and results in QT prolongation generated new interest in the role of INa and slow inactivation during repolarization. The role of slow inactivation in producing this phenotype was elegantly described by Clancy and Rudy (20) and demonstrated the importance of the use of Markov models that relate to the molecular basis of gating.
Recent advances in the molecular basis of voltage-gated ion channels indicates that activation and inactivation are not independent processes (38). Furthermore, similar mechanisms of inactivation appear to operate across a wide variety of channels. In Ca2+ channels, the Ca2+ dependence of inactivation is mediated by the highly localized increase of [Ca2+] that binds to calmodulin, which is constitutively bound to the Ca2+ channel (74, 75). Our model assumes that this regulation is mediated by an N-type-like inactivation “ball” associated with the channel. The Markov model implementation of this mechanism produces inactivation behavior that cannot be reproduced by a gating-variable type model and is distinct in its molecular basis from the allosteric model of Ca2+ channel inactivation used in other recent formulations.
Our model is able to reproduce experimental observations of Ca2+ fluxes in mouse myocytes. As with other animals, the Ca2+ cycle in the mouse begins with activation of L-type Ca2+ current ICaL, which transports a relatively small amount of Ca2+ across the sarcolemmal membrane into a small subspace volume, where it triggers Ca2+ release from the SR (10, 30). The Ca2+ released from the SR is considerably larger (∼40 μM) than the Ca2+ influx through L-type Ca2+ channels. Additional amounts of Ca2+ also enter the cell through a reverse mode of the Na+/Ca2+ exchanger. Ultimately, the elevation in [Ca2+]i leads to Ca2+ binding to troponin C and activation of myocyte contraction (10, 30). During the diastolic phase, Ca2+ gradually unbinds from troponin and other buffers, resulting in muscle relaxation. Excess Ca2+ is removed from the cytosol mostly through the SR Ca2+ pump and the Na+/Ca2+ exchanger. For mice and rats ∼90% of Ca2+ is pumped back to the SR and ∼10% of Ca2+ is extruded from the cell via the Na+/Ca2+ exchanger (15, 54). In rabbits, cats, dogs, and ferrets, all of which have longer APs, a considerably larger fraction of Ca2+ (30–40%) is extruded from the myocyte by the Na+/Ca2+ exchanger (11). In our model, 90% of the total Ca2+ transported during myocyte relaxation is carried by the SR Ca2+ pump and 10% by the sarcolemmal Na+/Ca2+ exchanger.
Many currents are more accurately described by Markov models rather than Hodgkin-Huxley formulations. For example, the complex voltage-dependent kinetics of the rapid delayed-rectifier current, IKr, can more accurately be reproduced by a Markovian model. IKr is encoded by the ether-à-go-go-related gene, HERG in humans and mERG in mouse. A number of familial and spontaneous mutations of the HERG channel are associated with arrhythmias in humans (44). We used the HERG model developed by Wang et al. (89) to reproduce mouse IKr (57). IKr is not consistently observed in mouse ventricular myocytes, and expression appears to vary with age. It is more prominent in embryonic and neonatal myocytes than in the adult. In adult myocytes, IKr is small relative to other currents (57). Accordingly, in our model of the mouse ventricular myocytes IKr is very small and does not significantly alter the shape and duration of AP under control conditions. However, the HERG current plays an important role in myocyte repolarization in the other species [guinea pig (41, 98), rabbit (55), and human (25, 72)] and may play an important role in species-dependent toxicity to K+ channel blockers. It may also underlie differences in restitution properties.
In the mouse ventricle, as in human and rabbit atrial cells, the rapid transient outward K+ current, IKto,f, dominates repolarization. Elimination of IKto,f can result in a marked increase in mouse ventricular APD and cardiac remodeling (7, 93). IKto,f plays an important role in repolarization of both the apex and the septum in the mouse (37). To evaluate the contribution of IKto,f in these regions, a point mutation (Trp to Phe) was introduced at position 362 in the α-subunit of Kv4.2 channels to produce a nonconducting mutant subunit (7). This “dominant negative” mutation eliminated IKto,f and increased APD90 from 36 ± 10 ms in wild-type mice to 116 ± 14 ms in transgenic mice (37). In our model eliminating IKto,f resulted in significant APD prolongation, but the effect was larger than that observed experimentally. The difference may be due to upregulation of other channels in the mouse.
There is a correlation between the amplitude of IKto,f and APD. In the apex, peak IKto,f is 34.5 ± 2.3 pA/pF and APD50 is 4.5 ± 0.3 ms, whereas in the septum peak IKto,f is only 6.8 ± 0.4 pA/pF and APD50 is 8.5 ± 0.3 ms (37). We used the experimentally determined amplitudes for IKto,f to simulate regional APs. The difference between the simulated APD50 from the septum and the apex was 3.6 ms, which is close to the corresponding experimental value of 4.0 ms (37).
The slow transient outward K+ current, IKto,s, also makes a substantial contribution to the heterogeneity of myocyte repolarization from apical and septal regions. IKto,s is absent in apical cells, but is especially important in the septum, where its amplitude is comparable to that of IKto,f. Most likely, IKto,s plays a compensatory role in repolarization in the septum because of the relatively small amplitude of IKto,f. IKto,s has been tentatively given a molecular basis of Kv1.4 and has been called Ito,s in some cases (36, 37). The inactivation rate of IKto,s is somewhat slower than the heterologously expressed α-subunit of the Kv1.4 channel, but IKto,s has the characteristic slow recovery properties of this channel.
The ultrarapidly activating delayed rectifier K+ current, IKur, and the noninactivating steady-state K+ current, IKss, both play a substantial role in repolarization. IKur has a tentative molecular basis of Kv1.5 (59). The mediator(s) for IKss is controversial (69, 92). Our model allows simulation of AP modification in knockout mice with modulation of IKur and IKss. However, it should be remembered that this is only an approximation of a “knockout” experiment, and several other currents are likely to show compensatory changes in addition to removal or modification of the target current (7, 37, 69).
Our model of the mouse ventricular AP reconciles regional differences in repolarization with the measured differential expression of IKto,f, IKto,s, IKur, and IKss. In the apex IKto,f is dominant, leading to a rapid early repolarization and a relatively short AP. In the septum IKto,s prevails. Differential expression and recovery from inactivation of the transient outward K+ current are also responsible for heterogeneity in rat ventricular AP (73). The rat AP is short compared with other species, with APD50 ≈ 13 and 38 ms for model epicardial and endocardial cells, respectively (73).
Heterogeneity of expression of IKto also contributes to heterogeneity of repolarization in larger animals, e.g., canine ventricle (56), although a clear IKto,f and IKto,s distinction has not been made in larger animals. Regional differences in IKto have also been found in cat, rabbit, and human ventricular myocytes (3). IKto is prominent in both canine epicardial and midmyocardial (M) cells, which have a “notch” in the AP. Conversely, endocardial cells have relatively small IKto and a less clear notch. There is a similar comparative reduction in IKto and notch in left versus right ventricular myocytes in dogs (26). There is also a difference in APD between epicardial and M cells. M cells have smaller IKs expression and a longer APD50 compared with epicardial myocytes (56), a mechanism that is unlikely to be active in adult mouse ventricles.
IKto,f is the dominant outward current in early repolarization or the notch phase of repolarization in mouse, rat, and human cardiac cells. It is absent, however, in bullfrog and guinea pig myocytes (41, 62, 63, 76, 77) and does not play an important role in the early repolarization of rabbit atrial cells (55). During the later phases of repolarization in mouse and rat ventricular and human atrial myocytes, IKur and ICaL act in opposite directions and play a more important role in repolarization. IKur repolarizes membrane potential, whereas ICaL depolarizes the cell to form a plateau (25) or slow down repolarization (72, 73). In rabbit atrial cells, the interplay between IKto,f and ICaL is responsible for the decaying plateau just after early stage of repolarization (55). In the guinea pig ventricular myocyte, the plateau is supported by the interaction of outward K+ currents [the plateau K+ current (IKp), IKr, and IKs; Ref. 98] and total inward Ca2+ current, which includes L-type and T-type Ca2+ currents (62, 63). In bullfrog atrial cells, the middle to late plateau is formed by the activation of the outward delayed-rectifier K+ current and reactivation of the inward Ca2+ current (76).
Mouse and rat APs do not have an elevated plateau, so after inactivation of IKur and ICaL final repolarization is mediated by a balance between IK1, INaCa, INaK, and the background currents. This is similar to the final stage of repolarization in human atrial myocytes, which involves IKr, IKs, and IK1 (25, 72). In rabbit atrial and guinea pig ventricular myocytes, IKr and IKs speed up repolarization after the plateau phase, with the role of IK1 and INaCa increasing at the end of repolarization (55, 62, 63, 98).
Late repolarization in the mouse ventricular AP (often measured as APD90) can be highly variable. It is arguable whether this phase represents a true repolarization “plateau” or an “afterdepolarization.” The high current densities and fast rate of beating of the mouse ventricle put a large demand on the ionic homeostasis mechanisms, especially the Na+/K+ pump and the Na+/Ca2+ exchanger. In this sense, it resembles the DiFrancesco and Noble (27) mechanism of repolarization, in which INaCa contributes significantly to the plateau phase of repolarization in Purkinje fibers. Late repolarization is relatively slow, with small dVm/dt and net membrane current. As a consequence, it is a limitation of the model that this phase may be sensitive to physiological and experimental factors, including accumulation and depletion of ions in restricted extracellular spaces, or buffering of [Ca2+]i by indicator dyes, such as BAPTA or fura-2. Similarly, pump and exchanger current densities may be altered indirectly in cells by activation of electroneutral transport mechanisms not included in the model, such as Cl−/HCO3− exchange or furosemide-sensitive Na+-K+-2Cl− cotransport.
Extracellular and intracellular ion accumulation can be an additional confounding factor. Flux through the Na+ channel is large and localized to the transverse tubule system (13). Flux through IKto,f is also large, further complicating quantitative measurement. Local accumulation and depletion artifacts can occur for both Na+ and K+ in the restricted extracellular space. Changes in ionic concentrations in these extracellular spaces have been reported to be of physiological and experimental significance. We have chosen to model extracellular ions as being in a well-stirred single extracellular space. With the exception of colocalization of Ca2+ channels and ryanodine receptors in narrow diadic spaces, the cytosol is also treated as lumped well-mixed space. This may limit the performance of the model in the face of sustained depolarization or fast repetitive firing.
One limitation of any model of cardiac cellular activity is reconciling disparate experimental data obtained under widely different experimental conditions. All of these data must be “corrected” by the modeler to correspond to some “physiological” reference state. Voltage-clamp data are obtained under experimental conditions different from those used to measure action potentials. For example, IKto,f is measured in the presence of divalent Ca2+ channel blockers. In many cases we have used several strategies to validate the “corrections” made. The magnitude and kinetics of the Na+ current were at least partially validated by simulation for dVm/dt and comparison with the corresponding experimental measurements. Voltage-clamp measurements of ICaL were performed in the presence of Ca2+ buffers (EGTA, BAPTA) that substantially change ICaL gating properties. Therefore, the mechanism of Ca2+ handling, including ICaL and the Ca2+-release system, was validated with the time course of [Ca2+]i transients measured during myocyte twitch. However, this measurement is at least partially affected by the binding of Ca2+ to the Ca2+ indicator. The gating properties of K+ currents responsible for the AP shape and duration at different levels of repolarization were tested with APs from different heart regions (the apex and the septum), where particular K+ channels have differential expression. For example, the gating properties of IKto,f obtained in the presence of divalent ions during voltage-clamp experiments were adjusted in accordance with experimental data (1) and validated by the APD50 from the apical cells, where it is primarily responsible for the early stage of repolarization.
Thus our computer model of the mouse ventricular myocyte successfully reproduces experimental voltage-clamp data for the major ionic currents: a fast Na+ current, an L-type Ca2+ current, a Na+/Ca2+ exchange current, transient outward K+ currents, a rapid delayed rectifier K+ current, an ultrarapid delayed rectifier K+ current, a noninactivating steady-state K+ current, a time-independent K+ current, and a slow delayed rectifier K+ current. It closely describes experimentally observed Ca2+ transients and Ca2+ fluxes during the AP and myocyte relaxation. Our model also reproduces the characteristically short-duration mouse AP and can account for the heterogeneity of AP shape from the apex and the septum. The heterogeneity in our model arises through simulated differential expression of transient outward and ultrarapid delayed-rectifier K+ channels. Our model uses Markov representations that are related to structure-function data, so it will provide an important tool for interpreting the complex changes observed in experiments on transgenic mouse myocytes.
APPENDIX: MOUSE MODEL SUMMARY
(A2) (A3) (A4) (A5) (A6) (A7) (A8)
(A9) (A10) (A11) (A12) (A13) (A14) (A15)
(A18) (A19) (A20) (A21)
l-type calcium current.
(A22) (A23) (A24) (A25) (A26) (A27) (A28) (A29) (A30) (A31) (A32) (A33) (A34)
calcium pump current.
na+/ca2+ exchange current.
calcium background current.
Fast Na+ current.
(A40) (A41) (A42) (A43) (A44) (A45) (A46) (A47) (A48) (A49) (A50) (A51) (A52) (A53) (A54) (A55) (A56) (A57) (A58) (A59) (A60) (A61) (A62) (A63) (A64)
Background Na+ current.
Transient outward K+ current IKto,f.
(A67) (A68) (A69) (A70) (A71) (A72) (A73) (A74)
Transient outward K+ current IKto,s.
(A75) (A76) (A77) (A78) (A79) (A80) (A81)
Time-independent K+ current.
Slow delayed rectifier K+ current.
(A83) (A84) (A85) (A86)
Ultrarapidly activating delayed rectifier K+ current.
(A87) (A88) (A89) (A90) (A91)
Noninactivating steady-state K+ current.
(A92) (A93) (A94) (A95)
Rapid delayed rectifier K+ current (mERG).
(A96) (A97) (A98) (A99) (A100) (A101) (A102) (A103) (A104) (A105) (A106) (A107)
Na+/K+ Pump Current
(A108) (A109) (A110)
Ca2+-Activated Cl− Current
This work was supported in part by the National Heart, Lung, and Blood Institute (R01-HL-59526-01 to R. L. Rasmusson; HL-59139, HL-69020, and HL-33107 to S. F. Vatner; and HL-62442 to S.-J. Kim); an Established Investigator Award to R. L. Rasmusson (9940185N) and a Scientist Development Grant to G. C. L. Bett (0430051N) from the American Heart Association; a National Science Foundation Knowledge and Distributed Intelligence Grant (DBI-9873173); R01-HL-59614; and a Pittsburgh Supercomputer Center Grant to V. E. Bondarenko (MCB010020P).
The authors thank Drs. Jacques Barhanin, Barry London, and Guy Salama for helpful comments and Dr. Jeanne Nerbonne for early access to data in press.
The costs of publication of this article were defrayed in part by the payment of page charges. The article must therefore be hereby marked “advertisement” in accordance with 18 U.S.C. Section 1734 solely to indicate this fact.
- Copyright © 2004 by the American Physiological Society