Ionic mechanisms underlying human atrial action potential properties: insights from a mathematical model

Marc Courtemanche, Rafael J. Ramirez, Stanley Nattel


The mechanisms underlying many important properties of the human atrial action potential (AP) are poorly understood. Using specific formulations of the K+, Na+, and Ca2+ currents based on data recorded from human atrial myocytes, along with representations of pump, exchange, and background currents, we developed a mathematical model of the AP. The model AP resembles APs recorded from human atrial samples and responds to rate changes, L-type Ca2+ current blockade, Na+/Ca2+ exchanger inhibition, and variations in transient outward current amplitude in a fashion similar to experimental recordings. Rate-dependent adaptation of AP duration, an important determinant of susceptibility to atrial fibrillation, was attributable to incomplete L-type Ca2+ current recovery from inactivation and incomplete delayed rectifier current deactivation at rapid rates. Experimental observations of variable AP morphology could be accounted for by changes in transient outward current density, as suggested experimentally. We conclude that this mathematical model of the human atrial AP reproduces a variety of observed AP behaviors and provides insights into the mechanisms of clinically important AP properties.

  • action potential morphology
  • action potential rate dependence
  • transient outward current
  • L-type calcium current
  • sodium/calcium exchanger

the human atrial action potential (AP) and ionic currents that underlie its morphology are of critical importance to our understanding of the electrical properties of atrial tissues in normal and pathological conditions. AP shape and underlying ionic current densities are frequency dependent and may be modulated by specific pharmacological agents (33). Their properties are known to change during the course of atrial arrhythmias or as a result of other cardiac pathologies (32, 39). In addition, APs have shown significant variability in morphology, ranging from triangular APs with no sustained plateaus to long APs with definite spike-and-dome morphology, even when recorded from cells isolated from the same atrial sample (4,33, 54). There is a need to integrate the information obtained from measurements of single currents to understand complex current interactions during the AP and their role in controlling AP morphology.

Until recently, limited data were available on the various currents underlying the AP. New findings have revealed important properties and significant interspecies differences. For example, the transient outward current (I to), present in human and rabbit atrial cells (8, 16, 25, 47), has been shown to recover from inactivation at least two orders of magnitude faster in humans than in rabbits (20, 47). A novel delayed rectifier current (I Kur) has been identified in human atrial cells (1, 18, 41, 48, 53) and is different from the sustained current (I sus) observed in rabbit atrial cells (12).I Kur also differs from delayed rectifiers in other species that show kinetically similar characteristics (3, 6, 29). Additional data have been obtained on the classic delayed rectifier currents (I Kr and I Ks) (54,55), the time-independent inward rectifier current (I K1) (30, 31), the fast Na+ current (I Na) (42), the L-type Ca2+ current (I Ca,L) (35) and its inactivation by voltage and intracellular Ca2+ (26, 49), and the Na+/Ca2+ exchanger current (I NaCa) (4).

Models of atrial cells based on animal data only have been published (14, 36, 40, 60). Although these models have provided valuable insights into the mechanisms underlying AP generation in these species, the existence of significant interspecies differences and the amount of human data available indicate a need for an AP model based specifically on direct measurements of human atrial currents. To address this issue, we developed a mathematical model of the AP based on ionic current data obtained directly in human atrial cells. When human data were insufficient to characterize a given atrial current fully, we supplemented with animal data and referred to existing published models of atrial and ventricular APs (14, 36, 37, 40, 60). The most recent published models of mammalian cardiac cells are those of Luo and Rudy (37) (LR2 model) and Linblad et al. (36) (LMCG model), based on measurements from guinea pig ventricular cells and rabbit atrial cells, respectively. Our model builds mostly on the work of Luo and Rudy to develop a working model of the human atrial AP.

Our primary goal was to develop a useful model of the AP from which we could gain insights into experimental observations made on human atrial cells and tissues and make predictions regarding the behavior of these cells under previously untested conditions. Throughout our description of the model, we clearly indicate the sources of data used to derive model parameters, how the parameters were obtained, and whether they were derived directly from data or modified to improve our representation of the AP. We then use the model to investigate the mechanism of experimentally measured AP rate dependence, the changes in AP morphology in the presence of pharmaceutical blockers of theI Ca,L and the Na+/Ca2+exchanger, and the variability in experimentally observed AP morphology.


Transmembrane potential
Gas constant
Faraday constant
Membrane capacitance
Human atrial action potential
Maximal AP upstroke velocity (dV/dtmax)
AP amplitude (peak AP voltage minus diastolic voltage at onset of AP)
AP overshoot (peak AP voltage >0 mV)
AP duration
APD measured from AP onset to 50% of total APA repolarization
APD measured from AP onset to 90% of total APA repolarization
Atrial fibrillation
Temperature adjustment factor, k(T ) =k(T0) Q10(TT0)/10


Sarcoplasmic reticulum


Junctional SR, SR release compartment


Network SR, SR uptake compartment


Cell volume


Intracellular volume


SR uptake compartment volume


SR release compartment volume

[X ]o

Extracellular concentration of ion X

[X ]i

Intracellular concentration of ion X


Calmodulin, sarcoplasmic Ca2+ buffer


Troponin, sarcoplasmic Ca2+ buffer


Calsequestrin, JSR Ca2+ buffer


Ca2+ concentration in release compartment


Ca2+ concentration in uptake compartment


Ca2+-bound calmodulin concentration


Ca2+-bound troponin concentration


Ca2+-bound calsequestrin concentration


Equilibrium potential for ion X

I ion

Total ionic current

I st

Stimulus current


Forward rate constant for gating variable x


Backward rate constant for gating variable x


Time constant for gating variable x


Steady-state relation for gating variable x

I Na

Fast inward Na+ current

g Na

Maximal I Na conductance


Activation gating variable for I Na


Fast inactivation gating variable for I Na


Slow inactivation gating variable for I Na

I K1

Inward rectifier K+ current

g K1

Maximal I K1 conductance

I to

Transient outward K+ current

g to

Maximal I to conductance

K Q10

Q10-based temperature adjustment factor

o a

Activation gating variable for I to

o i

Inactivation gating variable for I to

I Kur

Ultrarapid delayed rectifier K+ current

g Kur

Maximal I Kur conductance

u a

Activation gating variable for I Kur

u i

Inactivation gating variable for I Kur

I Kr

Rapid delayed rectifier K+ current

g Kr

Maximal I Kr conductance

x r

Activation gating variable for I Kur

I Ks

Slow delayed rectifier K+ current

g Ks

Maximal I Ks conductance

x s

Activation gating variable for I Ks

I Ca,L

L-type inward Ca2+ current

g Ca,L

Maximal I Ca,L conductance


Activation gating variable for I Ca,L


Voltage-dependent inactivation gating variable forI Ca,L

f Ca

Ca2+-dependent inactivation gating variable forI Ca,L

I p,Ca

Sarcoplasmic Ca2+ pump current

I p,Ca(max)

Maximal I p,Ca


Na+-K+ pump current

I NaK(max)

Maximal I NaK

f NaK

Voltage-dependence parameter for I NaK


[Na+]o-dependence parameter forI NaK

K m,Na(i)

[Na+]i half-saturation constant forI NaK

K m,K(o)

[K+]o half-saturation constant forI NaK

I NaCa

Na+/Ca2+ exchanger current

I NaCa(max)

I NaCa scaling factor

K m,Na

[Na+]o saturation constant forI NaCa

K m,Ca

[Ca2+]o saturation constant forI NaCa

k sat

Saturation factor for I NaCa


Voltage-dependence parameter for I NaCa

I b,Na

Background Na+ current

g b,Na

Maximal I b,Na conductance

I b,Ca

Background Ca2+ current

g b,Ca

Maximal I b,Ca conductance

I rel

Ca2+ release current from the JSR

k rel

Maximal Ca2+ release rate for I rel


Activation gating variable for I rel


Ca2+ flux-dependent inactivation gating variable forI rel


Voltage-dependent inactivation gating variable forI rel


Sarcoplasmic Ca2+ flux signal for I rel

I up

Ca2+ uptake current into the NSR

I up(max)

Maximal Ca2+ uptake rate forI up


Maximal Ca2+ concentration in NSR

I tr

Ca2+ transfer current from NSR to JSR


Ca2+ transfer time constant

I up,leak

Ca2+ leak current from the NSR


Total calmodulin concentration in myoplasm


Total troponin concentration in myoplasm


Total calsequestrin concentration in JSR

K m,Cmdn

Ca2+ half-saturation constant for calmodulin

K m,Trpn

Ca2+ half-saturation constant for troponin

K m,Csqn

Ca2+ half-saturation constant for calsequestrin

V rest

Resting membrane potential



We model the cell membrane as a capacitor connected in parallel with variable resistances and batteries representing the ionic channels and driving forces. The time derivative of the membrane potential V(with the assumption of an equipotential cell) is given bydVdt=(Iion+Ist)Cm Equation 1where I ion and I stare the total ionic current and stimulus current flowing across the membrane and C m is the total membrane capacitance. The total ionic current is given byIion=INa+IK1+Ito+IKur+IKr+IKs+ICa,L +Ip,Ca+INaK+INaCa+Ib,Na+Ib,Ca Equation 2The ionic and pump currents, along with the handling of intracellular Ca2+ concentration ([Ca2+]i) by the sarcoplasmic reticulum (SR) system, are described in Membrane Currents. The model keeps track of the [Ca2+]i as well as intracellular concentrations of Na+([Na+]i) and K+([K+]i). No extracellular cleft spaces are included in the model (extracellular ion concentrations are fixed). Figure 1 shows a schematic representation of the currents, pumps, and exchangers included in the model.C m is taken to be 100 pF. The length and diameter of the cell are set to 100 and 16 μm, respectively. Unless otherwise noted, physical units are as follows: time (t) is in milliseconds, V is in millivolts, C m is in picofarads, current density is in picoamperes per picofarad, conductance is in nanosiemens per picofarad, and concentrations are in millimoles per liter. Physical parameters for the model cell are given in Table 1. Cell compartment volumes are identical to those used by Luo and Rudy (37).

Fig. 1.

Schematic representation of currents, pumps, and exchangers included in model. Cell includes 3 intracellular compartments: myoplasm, sarcoplasmic reticulum (SR) release compartment [junctional SR (JSR)], and SR uptake compartment [network SR (NSR)]. SeeGlossary for definitions.

View this table:
Table 1.

Model constants

A modified Euler method is used to integrate Eq. 1 with a fixed time step Δt = 0.005 ms (see , Numerical Integration). All simulations were performed on a local workstation (SGI Indigo 2XZ R4400) using double-precision arithmetic. A description of the model current formulations and their general mathematical representations are given below. A detailed listing of all equations and parameters can be found in the .

Membrane Currents

Fast Na+ current.

The model implements I Na as a modification of the widely used Ebihara-Johnson model (15) proposed by Luo and Rudy (37). The current is given byINa=gNam3hj(VENa),gNa=7.8 Equation 3Using the formulation of Luo and Rudy, we find that the model steady-state inactivation and its time constants are shifted 20 mV positive compared with the data of Sakakibara et al. (42). However, the formulation for steady-state inactivation is consistent with more recent human atrial steady-state inactivation data from Schneider et al. (46) (see discussion for a detailed evaluation of available human atrial data). Maximum Na+ conductance (g Na) was adjusted to produce an appropriate value for the AP amplitude (APA) and maximal upstroke velocity ( max). We useg Na = 7.8 nS/pF, which is ∼1.3 times the temperature-adjusted value (with the assumption of Q10 = 3) reported by Sakakibara et al. Figure2 shows the steady-state activation and inactivation curves, as well as time constants, used in the model. Data from Sakakibara et al., shifted by the amounts described above and temperature corrected, are overlaid on the model curves. For activation time constants, temperature-adjusted data from Schneider et al. are shown. The peak current-voltage (I-V ) relationship for I Na is displayed in Fig.3. No experimental I-Vcurve was included, because available data (19, 42, 46), as a result of equal intra- and extracellular Na+ concentrations, display a nonphysiological reversal potential for I Na close to 0 mV.

Fig. 2.

Model representation of parameters describing fast Na+current. Steady-state activation and inactivation curves (A) and gating variable time constants (B) are shown. Experimental data from Sakakibara et al. (42) (m ,h j , τh, τj) and Schneider et al. (46) (τm), voltage shifted (+20 mV) and temperature corrected (Q10 = 3), are included for comparison. Experimental steady-state relationships are lines fitted by Sakakibara et al. to their experimental data.

Fig. 3.

Current-voltage (I-V) relationships of major currents in model with corresponding experimental data. Model curves (solid lines, dashed line without symbols forI Na) represent peak current measured during 3,000-ms voltage steps to various potentials from a holding potential of −80 mV. Experimental data sets (dashed lines with symbols) are scaled so that maximum current for a given experimentalI-V curve matches current amplitude of corresponding model I-V curve at that potential. Data sources and scaling factors are as follows: Li and Nattel (35) forI Ca,L scaled by 0.38, Feng et al. (19) forI Kur scaled by 0.81 and I toscaled by 1.2, Wang et al. (55) for I Ks scaled by 3.0 and I Kr scaled by 0.38. Refer to original studies for exact experimental voltage-clamp protocols used. Experimental I-V curves for I Kurand I to were shifted 10 mV negative before they were scaled, matching voltage shift introduced in model to account for known effects of divalent cations used experimentally to inhibitI Ca,L (13, 28). Experimental I-Vdata for I K (55) were transformed to pA/pF with use of C m = 100 pF before scaling was performed. No experimental data are shown for I Na.

Inward rectifier K+ current.

The model implements I K1 on the basis of available current data and resting membrane resistance measurements. The current is given byIK1=gK1(VEK)1+exp[0.07(V+80)],gK1=0.09 Equation 4We assume no temperature dependence of the inward rectifier current on the basis of preliminary experimental data. The current is assumed to reverse at E K, with a voltage-dependent conductance and no additional dependence on outward K+ concentration ([K+]o). Koumi et al. (30, 31) reported measurements of I K1from human atrial cells. When incorporated in the model, their large conductance value of 0.25 nS/pF at E K yielded a low input membrane resistance (∼40 MΩ for a 100-pF cell). We decreasedI K1 conductance to compensate for this using g K1 = 0.09 nS/pF. When the latter is used, input resistance calculated for a voltage step from −80 to −90 mV is ∼174 MΩ, closer to experimentally recorded values (e.g., ∼150 MΩ in Ref. 52). However, there is a large range of measured values for the input resistance of human atrial cells that may reflect differences in the isolation procedure and its effect on isolation-sensitive currents such as I K1. TheI-V relationship of I K1 is shown in Fig. 4.

Fig. 4.

I-V relationships for I K1 andI NaK obtained from model. Sum of instantaneous time-independent membrane currents (I v) is also shown. For I NaK, the amplitude of which depends on Na+ and K+ concentrations, resting concentrations from Table 2 are used.

Transient outward and ultrarapid rectifier K+currents.

The model implements I to andI Kur on the basis of data from our laboratory (20,53). The currents are given byIto=gtooa3oi(VEK),gto=0.1652 Equation 5 IKur=gKurua3ui(VEK), gKur=0.005+0.051+exp(V15)13 Equation 6Both currents share similar activation kinetics but differ in their inactivation properties. I to has fairly rapid inactivation kinetics; I Kur displays only partial slow inactivation. Currents were modeled at room temperature, and Q10 = 2.2 was used to adjust their kinetics to 37°C (53). Current amplitudes were adjusted to match AP morphology. Figure5 shows steady-state activation, inactivation, and corresponding time constants forI Kur, along with experimental data from our laboratory (53). Figure 5 shows similar results forI to. All model fits shown in Fig. 5 are shifted 10 mV negative for computations to account for the known positive shift (13, 28) induced by divalent cations such as Cd2+, which was used experimentally to block I Ca,L. The peakI-V relationship of both currents, along with scaled experimental I-V curves from Feng et al. (19), are shown in Fig. 3. Figure 6 displays traces for I Kur and I to in response to voltage-clamp steps from a holding potential of −80 mV. The traces ofI to + I Kur show typical rapid inactivation of I to followed by a sustained outward current carried by I Kur. The combined model current traces are overlaid with scaled experimental recordings from Wang et al. (53).

Fig. 5.

Model representation of parameters describing ultrarapid delayed rectifier current (A and B) and transient outward current (C and D). Steady-state activation and inactivation curves (A and C) and gating variable time constants (B and D) are shown. Experimental data from Wang et al. (53), to which model relations are fitted, are included for comparison.

Fig. 6.

I Kur and I to generated by model in response to voltage-clamp steps. Results are shown for combined currents (A), I to alone (B), and I Kur alone (C). From resting state, model is held at −80 mV for 1,000 ms, followed by potential steps to −20, 0, and +20 mV. Experimental recordings ofI Kur + I to, from Wang et al. (55), are included for comparison (dashed lines in A). Experimental recordings were scaled using C m = 100 pF to match peak amplitude of model trace at +20 mV.

Rapid and slow delayed rectifier K+ currents.

The model implements the two components of the classic delayed rectifier, I Kr and I Ks, on the basis of data from our laboratory (54, 55). The currents are given byIKr=gKrxr(VEK)1+expV+1522.4,gKr=0.0294 Equation 7 IKs=gKsxs2(VEK),gKs=0.129 Equation 8Both currents have only a single activation gate, althoughI Ks activation uses a squared (x 2) activation gate, as suggested in previous studies (61). Partial inactivation of I Kr is observed in some experiments (55) but is not included in this model. Figure 7 shows the steady-state activation curves and time constants for both delayed rectifier currents as well as experimentally measured time constants. The steady-state activation parameters are taken directly from fits to experimental data presented by Wang et al. (55). Some experimental evidence (27, 43, 45) suggests that activation and deactivation of I Ks may be multiexponential. To account for these measurements and obtain results consistent with the overall experimental literature, the model computations use time constants for I Ks activation that are one-half the fitted values shown in Fig. 7. In addition, observations of slow deactivation of the delayed rectifier at negative potentials (27, 43) were considered by using asymmetric relationships to represent the dependence of τx(r) and τx(s) on V. Current amplitudes were adjusted to match AP morphology. The I-V relationships of both currents, along with scaled experimental I-Vcurves from Wang et al., are shown in Fig. 3. Figure8 displays traces forI Kr and I Ks in response to 3,000-ms voltage-clamp steps from a holding potential of −80 mV. The model traces are overlaid with scaled experimental recordings from Wang et al. The combined model I K trace shows rapid and slow activation phases linked to activation of I Krand I Ks, respectively.

Fig. 7.

Model representation of parameters describing rapid and slow delayed rectifier currents. A: steady-state activation curves;B: gating variable time constants. Experimental data from Wang et al. (55), to which model relations are fitted, are included for comparison (open symbols in B).

Fig. 8.

I Kr and I Ks generated by model in response to voltage-clamp steps. Results are shown for combined currents (A), I Ks alone (B), and I Kr alone (C). From resting state, model is held at −80 mV for 1,000 ms, followed by potential steps to −20, 0, and +20 mV. Experimental recordings ofI Kr and I Ks, from Wang et al. (55), are included for comparison (dashed lines in B andC). With use of C m = 100 pF to obtain units of pA/pF, experimental recordings were further scaled by a factor of 3.5 (I Ks) and 0.47 (I Kr) to match peak amplitude of model currents at +20 mV.

L-type Ca2+ current.

The model implements the slow inward Ca2+ current on the basis of previous models and data from human atrial cells (35, 37, 40,49). The current is given byICa,L=gCa,LdffCa(V65.0),gCa,L=0.1238 Equation 9Our formulation follows that of Luo and Rudy (37) and includes voltage-dependent activation and inactivation gates, as well as a Ca2+-dependent inactivation gate. The steady-state activation relation was fit to human atrial data from Li and Nattel (35) (Fig. 9) and hence differs slightly from that of Luo and Rudy and Rasmusson et al. (40) by providing a better fit to the slower onset of activation observed in the experimental data. Because of the absence of activation time constant measurements from human atrial myocytes, we retained the activation time constant formulation reported previously (36, 37, 40). Steady-state voltage-dependent inactivation data from Li and Nattel were accurately fit by the relation of Rasmusson et al. (40) (Fig. 9). Hence, we retained their formulation in our model. Our formulation of steady-state Ca2+-dependent inactivation is slightly modified from that of Luo and Rudy to better match the Ca2+transient and reproduce the extent of Ca2+-dependent inactivation of I Ca,L observed experimentally. A fixed intrinsic time dependence (τf,Ca = 2 ms) is incorporated into the f Ca gate, as reported by Friedman et al. (21, 22) (see SR Ca2+ storage and release), to better reproduce the time course of Ca2+-dependent inactivation of I Ca,L.

Fig. 9.

Model representation of parameters describingI Ca,L. A: steady-state activation curves;B: gating variable time constants. Experimental data from Li and Nattel (35), to which steady-state relationships are fitted, are included for comparison (open symbols inA).

In experiments on human atrial cells, Sun et al. (49) showed two fast Ca2+-dependent inactivation time constants forI Ca,L, on the order of 10–40 and 300–400 ms, respectively, at room temperature. In addition, when monovalent cations or Sr2+ was used as charge carrier, a very slow voltage-dependent process was apparent with inactivation time constants in the range of 600–800 ms at room temperature. Our ability to reproduce these results is dependent on our formulation of the current, which implements two inactivation mechanisms forI Ca,L. We have selected the time constant for voltage-dependent inactivation in our model to match qualitatively the measurements of Sun et al. The voltage dependence of the time constant uses the form proposed by Luo and Rudy (37), but the magnitude of the time constant is scaled to better agree with the data of Sun et al. The time constant at 0 mV is ∼240 ms in the model, which is within the experimentally measured range with the assumption of a threefold temperature correction to 37°C. Ca2+-dependent inactivation is mediated directly through the f Cagate. We rely on the Ca2+ transient to drive Ca2+-dependent inactivation, with its time constants being set indirectly through the dynamics of Ca2+ release and reuptake. This provides rapid inactivation of I Ca,Lon a time scale on the order of 10 ms. The model does not account for the existence of a slower Ca2+-dependent inactivation, as identified by Sun et al. (see discussion). TheI-V relationships for I Ca,L, along with a scaled experimental I-V curve from Li and Nattel (35), are shown in Fig. 3. Figure10 A displays model traces forI Ca,L in response to voltage-clamp steps from a holding potential of −80 mV. The current traces exhibit rapid activation, followed by two distinct inactivation time scales. Rapid inactivation mediated by intracellular Ca2+ inactivates a large fraction of the current soon after the onset of the voltage step and is followed by a slower voltage-dependent inactivation. The model response after a voltage-clamp step to 0 mV is compared with a scaled averaged experimental recording from Sun et al. in Fig.10 B.

Fig. 10.

A: I Ca,L generated by model in response to voltage-clamp steps. From resting state, model is held at −80 mV for 1,000 ms, followed by potential steps to −40 to +20 mV.B: model response at a voltage-clamp potential of +20 mV, overlaid with experimental data from Sun et al. (49). Experimental data represent an average experimental trace from 7 voltage-clamp current recordings to +20 mV. Peak amplitude of normalized recording was scaled to match model peak amplitude at +20 mV.

Na+-K+ pump.

Our formulation of the Na+-K+ pump follows the model of Luo and Rudy (37). The amplitude of the pump current [I NaK(max)] is adjusted to maintain stable intracellular ion concentrations at rest (37, 59). The pump current is given byINaK=INaK(max)fNaK11+{Km,Na(i)/[Na+]i}1.5 ·[K+]o[K+]o+Km,K(o),INaK(max)=0.6 Equation 10The voltage dependence of I NaK at resting [Na+]i([K+]o is fixed at 5.4 mM in the model) is shown in Fig. 4.

Na+/Ca2+ exchanger.

The exchanger current follows the revised formulation of Luo and Rudy (37). The exchanger current is given by

Formula Equation 11

As with I NaK, the amplitude is adjusted to maintain stable physiological intracellular ion concentrations at rest (Table 2). Figure11 shows the voltage dependence of the exchanger current at resting [Na+]i for two values of [Ca2+]i: at rest ([Ca2+]i = 0.1 μM) and near peak Ca2+ release concentration ([Ca2+]i = 1 μM). The downward shift in the curve is associated with a shift from Ca2+ entry to Ca2+ extrusion mode at depolarized potentials as Ca2+ concentration increases.

View this table:
Table 2.

State variables of the AP model at rest

Fig. 11.

I-V relationships for I NaCaobtained from model. Current amplitude depends on Na+ and Ca2+ concentrations. Resting concentration shown in Table 2is used for [Na+]i. Two different Ca2+ concentrations are used to illustrate shift from reverse to forward exchanger modes at plateau potentials in response to an increase in intracellular Ca2+.

Background Ca2+ and Na+ currents.

The Ca2+ leak current is adjusted to maintain a stable [Ca2+]i at rest. The Na+ leak current is adjusted along with I NaK to maintain stable [Na+]i and [K+]i at rest. Both approaches have been used extensively in previous AP models (37, 59). The background currents are given byIb,Ca=gb,Ca(VECa),gb,Ca=0.00113 Equation 12 Ib,Na=gb,Na(VENa),gb,Na=0.000674 Equation 13

Ca2+ pump current.

A sarcolemmal Ca2+ pump is included in the model to maintain [Ca2+]i at physiological levels. The pump current formulation (37, 40) is given byIp,Ca=Ip,Ca(max)[Ca2+]i0.0005+[Ca2+]i, Ip,Ca(max)=0.275 Equation 14Figure 4 shows the total time-independentI-V relationship for the model, which includes the background currents, I K1, and the pump and exchanger currents.

SR Ca2+ storage and release.

SR Ca2+ uptake and release are implemented using a two-compartment model, as reported by Luo and Rudy (37). Intracellular Ca2+ is taken up into an SR uptake compartment [network SR (NSR), 5.52% of cell volume] coupled to a release compartment [junctional SR (JSR), 0.48% of cell volume]. Uptake is [Ca2+]i dependent and includes a leak from the uptake compartment back into the intracellular space. The main currents (in mM/ms) involved in SR Ca2+ handling are given byIrel=krelu2vw([Ca2+]rel[Ca2+]i),krel=30 Equation 15 Itr=[Ca2+]up[Ca2+]relτtr,τtr=180 Equation 16 Iup=Iup(max)1+(Kup/[Ca2+]i),Iup(max)=0.005 Equation 17

Triggering of Ca2+ release from the JSR is modified from the original formulation of Luo and Rudy (37) following the work of Friedman et al. (21, 22). Release from the JSR is triggered by Ca2+ flux into the cell with emphasis placed on flux through the closely coupled I Ca,L and JSR release channels. Release is controlled by two Ca2+ flux-dependent activation and inactivation gates, u and v. Steady-state activation and inactivation curves for u andv, as well as their respective time constants, are shown in Fig. 12. This trigger mechanism results in a transient activation of the SR Ca2+ release current in response to a suprathreshold Ca2+ flux. Friedman (21) showed that, for AP simulations, this formulation is equivalent to using a threshold on [Ca2+]i entry, as done by Luo and Rudy (37). The formulation used here has the advantage of being a continuous analytic representation, thereby avoiding possible discontinuities and artifacts in the cellular response near the threshold for Ca2+ release (21). The voltage-dependent inactivation gate w contributes to the decrease in SR release current amplitude at positive potentials, when membrane voltage approaches the I Ca,L reversal potential. The steady-state voltage-dependent inactivation curve for w and its time constant are also shown in Fig. 12. Detailed expressions for the flux computations and gating variables are given in the . In a recent experimental study on human atrial myocytes, Hatem et al. (26) recorded Ca2+ transients in response to 150-ms voltage-clamp pulses to 0 mV applied at frequencies of 0.1 and 1 Hz. Figure 13 shows a comparison of experimental data from Hatem et al. and model-generated Ca2+ transients elicited using the same voltage-clamp protocol. A detailed comparison of model and experimental Ca2+ transient morphologies can be found in thediscussion.

Fig. 12.

Model representation of parameters describing Ca2+-induced Ca2+ release. Steady-state activation and inactivation curves (A) and corresponding time constants (B) for Ca2+ flux-dependent gating variables u andv are shown. C: steady-state inactivation curve and corresponding inactivation time constant for voltage-dependent gating variable w.

Fig. 13.

Comparison of model-generated (solid lines) and experimentally measured (dashed lines) Ca2+ transients in response to 150-ms voltage-clamp pulses to 0 mV applied at a frequency of 0.1 and 1 Hz. Experimental data are based on published measurements of Hatem et al. (26) in human atrial myocytes.

Myoplasmic and SR Ca2+ buffers.

Ca2+ buffering within the cytoplasm is mediated by troponin and calmodulin. Ca2+ buffering in the release compartment is mediated by calsequestrin. In particular, calsequestrin is essential in providing sufficient Ca2+ stores for release within the limited JSR volume. Buffers are considered at equilibrium throughout our simulations, and incoming Ca2+in the myoplasm and JSR compartments is instantaneously equilibrated into free and buffer-bound fractions (see for computational methods). Buffer concentrations and binding constants are given in Table 1 and follow the models of Luo and Rudy (37) and Rasmusson et al. (40).


Resting Properties

The model has a stable resting potential near −81 mV. All intracellular concentrations are stable at rest with [Ca2+]i = 0.1 μM, [Na+]i = 11.2 mM, and [K+]i = 139.0 mM. As stated earlier, the resting input resistance measured as the current change from −80 to −90 mV is ∼174 MΩ. The first 100-ms period of Fig.14 shows the balance of membrane currents in diastole during stimulation at 1,000 ms. During the diastolic phase the steady state involves a balance between the pump and exchanger currents, the background currents, andI K1.

Fig. 14.

Model action potential during stimulation at 1 Hz. Membrane potential and corresponding ionic currents are shown. Model was stimulated from rest at a period of 1,000 ms. Model output from 12th action potential is depicted.

Model AP

Figure 14 shows a model AP generated during stimulation at 1,000 ms in response to 2-ms pulses of 2-nA amplitude (twice diastolic threshold). The time course of membrane currents during the AP is also shown in Fig. 14. AP duration (APD) at 50 and 90% repolarization (APD50 and APD90), APA, AP overshoot (APO), and max are reported in Table 3. The AP generated with the control parameters exhibits a spike-and-dome morphology commonly observed in human atrial AP recordings. This morphology is similar to experimentally recorded spike-and-dome APs reported by Benardeau et al. (4) and Wang et al. (53). Intracellular Ca2+ dynamics during the AP are shown in Fig. 15.

View this table:
Table 3.

AP characteristics and their rate dependence

Fig. 15.

Model action potential and intracellular Ca2+ dynamics during stimulation at 1 Hz. Membrane potential, intracellular Ca2+ flux (trigger for Ca2+ release), intracellular, uptake, and release compartment Ca2+concentrations, and uptake, transfer, and release Ca2+fluxes are shown. Model was stimulated from rest at a frequency of 1 Hz. Model output from 12th action potential is depicted.

AP Rate Dependence

Rate dependence of APD and AP refractoriness is an essential property of atrial cells that is central to our understanding of excitability and propagation patterns in the atria, particularly during reentrant arrhythmias. Table 3 summarizes the main properties of the model APs at various pacing periods. Figure 16displays model APs, after 12 s of pacing, for various stimulation periods ranging from 1,000 to 300 ms. There is no significant increase in APD at >1,000 ms. This is in agreement with data from Fermini et al. (20) that showed no change in AP morphology and APD on changing the stimulation period abruptly from 10 to 1 s. Larger relative rate-dependent changes in APD are measured with APD50than with APD90, as in experimental preparations (20, 32). Figure 16 also shows APs obtained in a tissue preparation at stimulation periods of 1,000, 600, and 300 ms (56). The experimentally observed changes in APD and AP morphology with stimulation rate are qualitatively similar to those observed in the model.

Fig. 16.

Rate dependence of action potential duration and morphology. Top left: model action potentials recorded after 12 s of pacing at various basic cycle lengths (BCLs). Action potential duration decreases as BCL decreases from 1 to 0.3 s. Top right: experimental recordings of action potentials obtained at pacing periods of 1, 0.6, and 0.3 s (data from Ref. 56). Bottom panels: currents obtained during model action potentials at BCLs of 1 and 0.3 s. Currents during 2 consecutive action potentials are shown at BCL of 0.3 s to facilitate comparison to currents at BCL of 1 s.

The mechanism underlying APD changes may be investigated by examining the detailed current traces for APs at stimulation periods of 1,000 and 300 ms depicted in Fig. 16. The traces reveal that an increase in stimulation rate conspicuously reduces I Ca,L. We inspected the time course of gating variables forI Ca,L, I Kr, andI Ks to quantify their role in rate-dependent AP changes. We first consider a comparison of the slow gating variablesf, f Ca, x r, and xs2 at the onset of the AP for stimulation periods of 1,000 and 300 ms. These slow variables are most susceptible to alteration at fast rates. ForI Kr and I Ks, we found that incomplete deactivation at 300 ms resulted in an increase of 0.05 forx r and 0.004 for xs2 compared with their values at 1,000 ms. This compares with reductions of 0.2 for fand 0.16 for f Ca. The maximal conductances ofI Ca,L and I Ks are comparable, whereas the maximal conductance of I Kr is four times smaller. Hence, the changes in the f andf Ca gate result in larger changes in current at the faster rate, pointing to I Ca,L as an important determinant of AP rate dependence. However, the absolute magnitude of current changes with rate must be considered relative to the amplitude of other currents activating at a similar stage during the AP. For this reason, a small rate-dependent change in I Kamplitude may have significant effects during late repolarization, whenI K is the dominant repolarizing current.

To better ascertain the relative importance of I Kand I Ca,L in AP rate adaptation, we use a simulation protocol where the rate-induced changes in a specific current are identified and implemented individually to determine its specific contributions to AP rate dependence. We begin with a simulation of pacing at 1,000 ms that is interrupted just before the onset of a stimulus. Before continuation of the simulation and application of the scheduled stimulus, the values of the gating variables for the current under study are substituted with their values just before a stimulus during pacing at 300 ms. When the simulation is continued and the scheduled stimulus is applied, the resulting AP displays the rate-adaptation properties associated with the current under study, the state of which has been reset to its value during rapid pacing. We have carried out this procedure to examine the rate-dependent effect of I Ca,L (via the dand f gating variables), I Kr (via the x r gating variable), I Ks(via the x s gating variable), and [Ca2+]i (via [Ca2+]i, [Ca2+]up, [Ca2+]rel, and the u, v, and w gating variables) on the AP.

Figure 17 displays the APs after 12 s of pacing at 1,000 and 300 ms, along with the APs arising from isolating the rate-adaptation effects of I Ca,L,I K(I Kr + I Ks), andI Ca,L + I K. We note that when the effect of either I Kr or I Ksalone was examined, the observed change in APD90 was ∼50% of the total reduction seen with both currents combined (I K in Fig. 17). Also, the effect on APD90 of [Ca2+]i andI Ca,L combined was similar to that ofI Ca,L alone, which is why the results for [Ca2+]i were not included in Fig. 17. The limited role of [Ca2+]i at fast pacing rates may reflect the opposing effects of a higher resting [Ca2+]i and a lower Ca2+transient amplitude on I Ca,L inactivation. Our first observation is that the effect of the rate-dependent decrease in available I Ca,L alone is to lower the plateau potential and to significantly accelerate early repolarization. The lower plateau level reduces I K activation, slows terminal repolarization, and results in no net change in total APD. Second, we observe that the effect of the rate-dependent increase in available I K alone is to lower the plateau potential and to compensate for the effect of the reduced plateau level on I K activation, leaving the rate of terminal repolarization unchanged. The effect of I K rate adaptation alone explains less than one-half of the total rate-dependent reduction in APD90. Third, the effect ofI Ca,L and I K combined is synergistic, lowering the plateau potential and accelerating late repolarization relative to the effect of I Ca,Lalone. Changes in these two currents together are sufficient to account for the total rate-dependent AP shortening observed on reducing the pacing period from 1,000 to 300 ms. Hence, we conclude thatI Ca,L and I K act together to produce the observed rate-dependent AP shortening.

Fig. 17.

Role of I Ca,L and I K in producing rate-dependent shortening of action potential on reducing pacing period from 1,000 to 300 ms. Steady-state action potentials during pacing at 1,000 and 300 ms are shown. In addition, action potentials resulting from rate adaptation of I Ca,Lalone, I K(I Kr + I Ks) alone, andI Ca,L + I K are shown. Seeresults for a detailed description of numerical simulation procedure. Combination of I Ca,L andI K rate adaptation is sufficient to explain full extent of rate-dependent action potential shortening on reduction of pacing period from 1,000 to 300 ms.

ICa,L and Rate Dependence

Li and Nattel (35) showed that rate dependence can be abolished in isolated human atrial myocytes by blocking I Ca,L. To compare these experimental data with the behavior of the model, we evaluated the effects of a strong decrease in L-type Ca2+channel conductance on model AP rate dependence. Figure18 A shows model APs after 12 s of pacing with stimulation at 5,000 and 300 ms with a 90% reduction inI Ca,L conductance. Figure 18 B shows corresponding experimental results from Li and Nattel in which nifedipine was used to block I Ca,L. In the model and experiments, block of I Ca,L abolishes AP rate dependence. Although we have shown that I Ca,L andI K play a role in rate adaptation, it appears that when I Ca,L is strongly inhibited, as in Fig. 18, the plateau level is lowered to the point where I Kactivation is greatly reduced and can no longer contribute to rate adaptation.

Fig. 18.

Role of I Ca,L in rate dependence of action potential. A: model action potentials recorded after 12 s of pacing at various BCLs with 90% decrease in I Ca,Lconductance. No rate dependence is observed. B: experimentally recorded action potentials during stimulation at BCL = 5 and 0.3 s in presence of 10 μM nifedipine to blockI Ca,L. No rate dependence is observed.

Block of INaCa and Role in APD

Recently, Benardeau et al. (4) investigated the role of Na+/Ca2+ exchange in AP morphology and APD. In particular, they examined the effect of I NaCa block by Li+ on AP morphology. Their results revealed a shortening of APD and transient hyperpolarization after the AP in the presence of Li+. Figure 19shows control and I NaCa-blocked (90% decrease in maximum exchanger activity compared with control) model APs with stimulation at 1 s. Model APs display both of the changes observed in experimental preparations (see Fig. 7A in Ref. 4), i.e., AP shortening and hyperpolarization. Inspection of the current records shows that the decrease in APD may be related to two mechanisms:1) an increase in Ca2+-dependentI Ca,L inactivation caused by an increase in resting [Ca2+]i and Ca2+ transient amplitude and 2) a reduction in depolarizing exchanger current during the late phase of the AP. The observed membrane hyperpolarization is due to a reduced inward exchanger current at rest, causing a negative shift of the resting membrane potential.

Fig. 19.

Model action potentials during stimulation at 1 Hz under control conditions and with a 90% decrease in I NaCamaximal amplitude. Block of I NaCa is meant to mimic experimental conditions where Li+ is used to block exchanger. Action potential shortening and slight hyperpolarization are observed, as in experimental recordings (see Fig. 7A in Ref.4).

Variability in AP Morphology

Using the model, we attempted to simulate changes in current density that could explain the changes in human atrial AP morphology seen experimentally by many investigators (1, 4, 17, 24, 50, 51, 54). It has been suggested that variations in the relative density of the transient outward current may contribute to some of the AP variability (17, 54). We concentrated on three currents: I to,I Kur, and I Ca,L. Figure20 shows the result of varying the density of each of these three currents, from one-tenth to three times the control conductance value, on AP morphology. Figure 20,AD, shows a continuum of AP morphologies after 12 s of pacing at a stimulation period of 1,000 ms for various modifications of the current densities. ForI Ca,L, we observe an increase in APD90 as I Ca,L conductance increases. There is little change in APD90 as I Kurconductance decreases. For I to, we observe an increase in APD90 as I to conductance decreases for large conductances, then a decrease in APD90as I to conductance decreases for small conductances. Of the three currents we have investigated, it appears that only the I to variations can explain the whole spectrum of AP morphologies observed experimentally (4, 33, 54). This ranges from rectangular APs with a short plateau at the lowestI to conductance values to shorter triangular APs at high I to conductance values, through long spike-and-dome morphologies at intermediate I toconductance values. Both I Kur andI Ca,L alone can explain a change from a spike-and-dome to a triangular AP shape, but a reduction inI to is required in the model to produce the high plateau and absence of notch characteristic of a rectangular AP (54). In addition, the relation between I to and AP morphology in the model parallels the relative changes inI to observed experimentally in the AP classification of Wang et al. (54).

Fig. 20.

Effect of varying I Ca,L (A),I Kur (B), and I to(C), conductances (10, 25, 50, 75, 100, 150, 200, 250, 300% of control) on model action potential morphology during stimulation at 1 Hz. D: experimental action potential types described by Wang et al. (54).


We have developed a mathematical model of the human atrial AP. Whenever possible, we used data directly measured on human atrial cells to derive model parameters. When human data were not available to completely characterize a given current, we relied on existing AP models based mainly on guinea pig ventricular (37) and rabbit atrial (36) data. The model yields AP morphologies that are consistent with a variety of experimental observations and gives potential insights into their underlying ionic mechanisms.

Model Current Formulations Compared with Experimental Data

Fast Na+ current.

I Na is a large, rapid inward current that has proved difficult to measure accurately using voltage-clamp techniques. Data from isolated human atrial cells are available elsewhere (19, 42). However, these data for steady-state activation and especially inactivation of the current are shifted to negative potentials compared with other human (7, 46) and animal (57) data. According to the steady-state inactivation curves in these studies, ∼90% of the available I Na would be inactivated at resting potentials positive to −80 mV. In conjunction with the small maximum conductances measured in these human cells (40 pA/pF maximum current in Ref. 19 compared with 200 pA/pF in Ref. 57 for rabbit atrial cells), there would be insufficient current available to generate the large max (167 V/s in Ref. 56) and nanoampere-range current amplitudes typical of the AP upstroke. In their model, Lindblad et al. (36) use the temperature correction suggested by Colatsky (10) (+3 mV/10°C) to shift the model curves to more positive potential. In the case of the human data reported by Feng et al. (19) and Sakakibara et al. (42), the resulting shift of ∼5–10 mV still leaves an inactivation half-maximal potential (V ½) well negative of the value observed in other mammalian cells and insufficient to generate the required fast inward current during the AP.

Recent human atrial data from Schneider et al. (46) display more realistic values for the V ½ of steady-state inactivation. The I-V curve of Schneider et al. resembles that measured by Sakakibara et al. (42), and both show a profile similar to the model I-V curve of Fig.3, although the data of Sakakibara et al. display a less abrupt activation of the current. The time constants of activation measured by Schneider et al. agree with those of the model; however, rapid inactivation time constants are on the order of 1 ms for a wide range of potentials (between −60 and +90 mV). The model displays fast inactivation time constants in this range for potentials positive to −40 mV, but at more negative potentials (between −80 and −40 mV), fast inactivation time constants increase in the model (to a maximum of 25 ms at about −65 mV). The rapid kinetics of fast inactivation over a wide voltage range reported by Schneider et al. may be due in part to their use of a fitting procedure allowing for overlap of activation and inactivation processes. It is unclear whether fast inactivation slows down in the study of Schneider et al. at potentials negative to −60 mV. They do not provide any fast time constants for recovery from inactivation at potentials between −95 and −60 mV, making it difficult to determine the value of τh at these potentials. On the basis of their recovery from inactivation data, there appear to be three time constants involved, whereas the data of Sakakibara et al. suggest two time constants. Given the difficulties in interpreting the data in these studies within the context of the simple formulation of Eq. 3, we chose to formulate the kinetics ofI Na exactly as in the work of Luo and Rudy (37), which nevertheless matches experimental data quite well, as shown in Fig. 2.

Transient outward current.

I to has been shown to differ substantially in its faster recovery kinetics in humans than in other mammalian species (20,47). On the basis of voltage-clamp current measurements, we found thatI to activation was best fit using a cubed activation gate. Activation and inactivation were fit to data from Wang et al. (55). A single inactivation gate was included in the model. In some current traces, slow inactivation can be seen, but this was attributed to slow inactivation of I Kur and was not included directly in our formulation of I to. The simulated current recordings shown in Fig. 6 indicate that the model reproduces the essential properties of experimental recordings ofI to + I Kur (55).

Ultrarapid delayed rectifier current or sustained outward current.

I Kur has been shown to play an important role in AP current balance at plateau potentials in human atrial myocytes (1, 47,55). On the basis of data from Wang et al. (55), we found thatI Kur had activation kinetics similar to those ofI to but activated at more negative potentials and displayed only partial slow inactivation. Those features are included in our formulation of the model. As with I to, changes in I Kur density can have considerable effects on AP morphology, as shown in Fig. 20.

L-type Ca2+ current.

I Ca,L has been studied extensively in human atrial tissue (9, 32, 35, 39, 49). Its complex inactivation kinetics, mediated in part by intracellular Ca2+, are only partially understood (49). This has made modeling the current more difficult and highly dependent on an accurate formulation of intracellular Ca2+ homeostasis and its control by the SR system (seeSR Ca2+ handling and buffering). Our formulation follows that of Luo and Rudy (37) and includes voltage-dependent activation and inactivation gates as well as a Ca2+-dependent inactivation gate. We have made modifications to the original formulation of Luo and Rudy on the basis of available experimental data (35, 49). However, as discussed by Luo and Rudy, fundamental questions regarding the complex voltage- and Ca2+-driven inactivation of I Ca,Lremain, from an experimental and a modeling standpoint. For example, data from Sun et al. (49) suggest that strictly voltage-dependent inactivation of I Ca,L is incomplete even for very long inactivation pulses. This suggests that steady-state inactivation ( f ) saturates at an inactivation fraction greater than zero (0.45 measured in Ref. 49). As reported by Luo and Rudy, we found that implementing such a relationship gave rise to uncharacteristic results, such as large I Ca,Lwindow currents observed on reuptake of Ca2+ during long APs. Another example can be found in the data of Sun et al., where an apparent second Ca2+-dependent inactivation time constant, on the order of 300–400 ms, could be observed during inactivation of I Ca,L. This is not accounted for in the model, and its origin has not been clarified experimentally. It may be that additional dynamical complexity in I Ca,Linactivation arises as a result of a more complex time course of the Ca2+ transient. This is in agreement with data from Hatem et al. (26) that point to rapid and slow phases during the rise of the Ca2+ transient in human atrial myocytes. These factors are related to limitations in our modeling of the SR Ca2+handling, which represents a simplified formulation of a temporally and spatially complex intracellular process (see SR Ca2+ handling and buffering andLimitations).

Delayed rectifier current.

I K has been characterized by Wang et al. (54) in human atrial cells. As in other mammalian species (44, 45),I K in human myocytes is the sum of a rapid and a slow component. We adopted the current formulation of Zeng et al. (61) fit to human atrial data (54). Steady-state activation and inactivation curves from human atrial cells are similar to those in other species (e.g., rabbit atria in Ref. 45). As pointed out by Zeng et al., there is conflicting evidence regarding the activation and deactivation kinetics of both currents. Wang et al. concluded that activation was well fit with a single time constant for both currents; however, they measured extremely slow activation time constants forI Ks. When included in the model, these slow kinetics limit the role of I Ks during the APO. Recent evidence from experiments on rabbit and guinea pig ventricular myocytes suggests fast and slow time constants for activation and deactivation of I Kr and I Ks(27, 43). Fast time constants for activation of I Ksreached a maximum of 330 ms at +10 mV and 37°C. We account for possible faster kinetics of I Ks in the model by using time constants that are one-half the values measured by Wang et al., putting the maximum time constant for I Ksactivation at ∼500 ms at +10 mV. By use of this formulation,I Ks plays a more significant role in repolarization during the AP but remains much less important thanI Kr at the lower plateau potentials typical of the human atrial AP.

SR Ca2+ handling and buffering.

SR Ca2+ handling and buffering involve uptake, storage, and release of Ca2+ from various cellular compartments (58). We followed the formulation of Luo and Rudy (37) but included a modification to the release mechanism suggested by Friedman et al. (21,22). Release of Ca2+ from the JSR compartment depends on the flux of Ca2+ through the sarcolemmal and SR release channels. This formulation aims to reflect the close coupling of L-type Ca2+ channels in the sarcolemma and JSR Ca2+release channel (34). For closely juxtaposed channels, one can expect that true local Ca2+ concentration in the vicinity of the release channel is more accurately reflected by local Ca2+flux than by bulk cellular Ca2+ concentration. The formulation is also adapted to afford a greater sensitivity of the release mechanism to flux through the L-type Ca2+ channel, as opposed to Ca2+ entry through the Na+/Ca2+ exchanger (22). Release is triggered when the computed Ca2+ flux exceeds a predefined threshold through the opening of an activation gate u. Inactivation is mediated through one flux-dependent (v) and one voltage-dependent (w) inactivation gate. An added advantage of this formulation is that it is completely analytic and does not require identification of qualitative trigger events such as the occurrence of max.

As shown in Fig. 13, Ca2+ transients elicited in the model are qualitatively similar in amplitude and late time course to the experimental transients recorded in response to voltage-clamp steps to 0 mV (26). The discrepancy in the initial rise of the transient has potentially complex causes. The experimental transients were recorded at room temperature, which may in part explain their slower rise. In addition, their rising phase included rapid and slow components. It was postulated by Hatem et al. (26) that the slow rising phase was due to Ca2+ release from sites not closely juxtaposed to sarcolemmal L-type Ca2+ channels. Because of the absence of intracellular Ca2+ diffusion or specific myoplasmic Ca2+ compartments in our model, it is unrealistic to expect to be able to directly reproduce these features of the transient. We may consider the modeled transient as a compromise formulation that reproduces the rapid rise in local Ca2+ concentration expected in the vicinity of the sarcolemma. This ensures our ability to use [Ca2+]i as a direct inactivator ofI Ca,L through the f Ca gate. An updated formulation of intracellular Ca2+ handling should include local myoplasmic compartments in an attempt to capture the essential features of intracellular Ca2+ diffusion.

Inward rectifier K+ current.

I K1 plays a major role in the late repolarization phase of the AP and in determining resting membrane potential and resistance. As stated earlier, there is considerable variability in the measured input resistance of isolated human atrial cells (1, 52), with resting potentials often displaced to positive potentials compared with multicellular tissue preparations (1, 52, 54, 56). Our formulation ofI K1 was adjusted to reproduce the slow late repolarization rate observed in human atrial AP recordings. In doing so, we were able to produce a realistic resting membrane resistance of 180 MΩ (52) and a resting potential characteristic of healthy human atrial multicellular preparations of around −80 mV (56). The modelI-V relationship for I K1 (Fig. 4) has a shape characteristic of direct recordings in human atrial cells (31), but the specific parameters were adjusted to match the AP characteristics described above.

Behavior of the Model AP

Control model AP.

The control model AP exhibits a spike-and-dome morphology. As discussed below, several AP morphologies have been recorded from human atrial cells. In the model the transient outward current produces the initial phase 1 repolarization after the I Na-driven upstroke. On inactivation of I to,I Ca,L, which remains after rapid Ca2+-dependent inactivation, causes a slight depolarization of the membrane and produces the dome of the AP. This net inward current is eventually countered by the slowly activating outwardI K, resulting in slow repolarization. Slow late repolarization is a balance between deactivatingI K, time-independent I K1, and opposing inward current generated by the Na+/Ca2+ exchanger during the late phase of the AP.

Na+/Ca2+ exchanger.

The Na+/Ca2+ exchanger plays an important role in determining APD. As in experiments, we have shown that the exchanger is partly responsible for producing the slow late repolarization typical of human atrial APs. Block of the exchanger resulted in a significant reduction in APD90 (30% in Fig. 19). Because of the slow repolarization rate and small current amplitudes during late repolarization, APD90 is sensitive to small changes in late AP currents such as I K,I K1, and I NaCa. This may explain some of the variability in APD90 measurements from human atrial cells and tissues (5, 23, 24, 32, 33).

Rate dependence of model AP.

Rate dependence of the model AP is described in Fig. 16 and 17 and Table 3. The model is consistent with experimental data in several respects. First, it exhibits little rate dependence at frequencies <1 Hz, consistent with data from Fermini et al. (20) and Le Grand et al. (32). The model predicts greater relative AP shortening at APD50 than at APD90. This is seen in the data of Le Grand et al., even though their APD90 values are somewhat greater than those of the model. Comparing the top panels of Fig. 16, we observe that in experimental recordings and model simulations the AP waveforms at various rates tend to converge at some positive potential before complete recovery. Late repolarization is slowed at faster rates, in part because of a reduced activation ofI K in APs with shorter, more negative plateaus that opposes the rate-dependent increase in I Kavailability. This tends to reduce rate adaptation of the model AP measured using APD90 and is a direct consequence of the sensitivity of late AP repolarization to small changes inI K, I K1, andI NaCa, as discussed above.

Variability in AP morphology.

Variability in AP morphology is often observed in recordings of human atrial APs. Early observations of APs from isolated human atrial tissue demonstrated a considerable variability in their morphology (17, 24,50, 51). Some cells showed little phase 1 repolarization and a positive plateau; others exhibited a prominent phase 1 with a “spike-and-dome” morphology (17, 24, 50). It was suggested that age-related differences may be due to changes inI to (17). Recently, there have been attempts to define various types of APs on the basis of their morphological characteristics and/or underlying ionic current densities (1,4, 54). Wang et al. (54) identified three “types” of APs on the basis of morphology, ranging from type 1, a rectangular AP with a positive plateau, through type 2, a spike-and-dome AP with a plateau at ∼0 mV, to type 3, a triangular AP with little plateau. They associated these changes in morphology with a measured increase in the relative density of I to to I Kfrom type 1 to type 3 APs. Benardeau et al. (4) recently characterized two AP types in their isolated human atrial cell preparations: type A APs are spike-and-dome APs with a high plateau similar to Wang’s type 1 or 2 APs; type B APs were triangular APs similar to Wang’s type 3. In an attempt to identify the source of differences in AP morphology between types A and B, these investigators found no discernible difference in I Ca,L density between the two cell types. All three AP types have been recorded from human atrial multicellular preparations with fine-tipped microelectrodes (38).

The control model AP is typical of the type 1 APs described previously (38, 54). Our investigation of the role of I to,I Kur, and I Ca,L in AP morphology in Fig. 20 reveals that, in our model, variations inI to alone are able to generate many of the variations in AP morphology observed experimentally. AP shapes range from a triangular morphology (type 3 in Refs. 38 and 54 and type A in Ref. 4) for large I to, through a spike-and-dome morphology (type 1 in Refs. 38 and 54) at intermediateI to conductances, to a rectangular morphology (type 2 in Refs. 38 and 54 and type B in Ref. 4) for smallI to amplitudes. In contrast, variations inI Kur and I Ca,L alone cannot reproduce the rectangular morphology observed only in the presence of a small I to. The role of I to as an important determinant of AP morphology is supported by current measurements (54) showing that changes in AP morphology types (1, 2, and 3) were accompanied by changes in the relative amplitude ofI to consistent with the results of Fig. 20. An increase in I to density was also shown to accompany the change from rectangular to spike-and-dome AP morphology observed in young vs. adult human atrial myocytes (17). Although the amplitude ofI to appears to be an important determinant of AP morphology, it is clear that it cannot be considered independently of other membrane currents and that AP morphology will depend on the relative amplitudes of all currents involved.

Clinical Relevance

Atrial fibrillation (AF) is the most common sustained arrhythmia in clinical practice and that for which treatment remains the most problematic (38). The AP model provides insights into the basic mechanisms underlying a variety of important clinical determinants of AF. Impaired atrial APD accommodation to rate has been observed in patients with AF (5) and likely underlies the flat refractory period-heart rate relationship of patients susceptible to AF (2). Our model points to rate-dependent I Ca,L inactivation combined with incomplete I K deactivation as the basic mechanism of APD accommodation to rate. Work by Li and Nattel (35) suggests that inhibition of I Ca,L alone is sufficient to abolish rate adaptation of the human atrial AP, and it has been shown that human myocytes from dilated atria have a substantially reduced I Ca,L (32, 39). In the present study, strong I Ca,L reduction alone was sufficient to abolish rate adaptation of the model AP (Fig. 18). Although I K and I Ca,Lcontribute to model AP rate adaptation, important reductions inI Ca,L lower plateau height to the point whereI K is greatly reduced and AP accommodation is virtually abolished. The model explains how strong reductions inI Ca,L, as caused by atrial disease in humans or by exposure experimentally to I Ca,L blockers, can largely abolish AP rate adaptation and contribute to the susceptibility to AF.

Another property strongly associated with the occurrence of AF is AP heterogeneity (38). The present model indicates that variations in plateau currents (particularly I to) over the experimentally measured range can account for much of the AP variability reported in human atria. These mechanistic insights hold the promise of an understanding of how drugs that alter ionic currents can affect these important AP properties and, thus, may help in the development of improved therapeutic approaches to treating atrial reentrant arrhythmias like AF.

Potential Limitations

Although we have tried to develop a model that is closely based on experimental data, limitations exist because of availability of data and considerable variability in existing experimental data. Limitations of this work are similar to those discussed for recent ionic models (36). The following issues should be considered in evaluating the model.

Data from human atrial cells and tissues typically display considerable variability, with a wide range of current densities having direct consequences on AP shape, rate dependence, and response to pharmacological intervention. The use of averaged data from selected preparations does not take this into consideration and introduces possible complications in comparing model output with experimental AP recordings under various experimental conditions. This is in part why we have considered carefully the effect of varying specific current conductances on AP morphology. In addition, the temperature dependence of current kinetics and amplitude is often poorly characterized and must be approximated, with the help of available data, during the course of model development.

When possible, we have attempted to show experimental data along with corresponding simulated behavior. However, we note that, for voltage-clamp current recordings and experimental I-Vcurves, there may be significant bias in the measured current amplitudes. Selection factors such as the ability to record measurable currents and the viability of cells contribute to bias, especially with currents that are difficult to measure or display significant variability, such as I Ca,L orI K. These amplitude-dependent selection factors are unlikely to bias the measured kinetics of a given current. Because of bias in current amplitude measurements, it is rare that exact current densities reported in experimental studies can be incorporated directly into a model and produce a realistic AP morphology. When necessary, current densities from experimental data were adjusted in the model to improve overall AP properties. The same limitations apply to direct comparisons of current recordings to corresponding model output. Such comparisons can be misleading unless one selects a current trace with an amplitude that exactly matches the amplitude used in the model or unless a scaling factor is used to match the peak current amplitudes. This is why we chose to present scaled experimental data for current recordings and I-V curves, with the scaling factor selected to match the maximum current amplitudes in the experiments to those in the model. Whenever such scaling is carried out, we ensured that the corresponding scaling factor was explicitly given to allow comparison of current magnitudes between model and experimental data.

Limited data are available on intracellular Ca2+ dynamics, Ca2+ storage and release kinetics from the SR, and its role in the inactivation of I Ca,L. We have used a modified version of the intracellular Ca2+-handling scheme proposed by Luo and Rudy (37). This produces a realistic intracellular Ca2+ transient in response to suprathreshold Ca2+ flux into the cell. However, recent experimental results suggest that Ca2+-dependentI Ca,L inactivation may involve complex local phenomena requiring an additional level of model complexity (49), beyond an accurate representation of the average Ca2+transient. This is also important in view of the sensitive dependence of the model AP and its rate dependence on the dynamics of late repolarization, controlled in part by the Ca2+-dependentI NaCa.

Our formulation of the delayed rectifier current represents a compromise between a simple current representation and available experimental evidence that I Kr andI Ks may display complex, multiexponential deactivation kinetics (27, 43). Slow deactivation of these currents may have important effects on late repolarization, refractoriness, and AP rate dependence. Whenever possible, we have tried to select the simplest formulation that captures the essential features of available experimental data. It is possible that the complex kinetics of some currents, such as I K, require more detailed current formulations involving several channel states not readily modeled using a simple Hodgkin-Huxley formalism.

We have chosen not to consider time-varying extracellular ion concentrations in our model at this stage. Certain currents, such asI K1, are known to depend on extracellular ion concentrations. Accumulation of ions in the extracellular cleft space may be an important modulator of ionic currents and AP characteristics at fast pacing rates. These phenomena may be taken into account by modifying appropriate current formulations and including dynamic extracellular cleft space ion concentrations in the model.

Finally, our model was adjusted to produce stable ionic concentrations at rest. In response to periodic stimulation, ionic balance is disturbed and slow changes in intracellular ionic concentrations occur in the model. Initial transients due to kinetic rate adaptation of the currents dissipate over a few seconds. However, small changes in ionic concentrations can take several minutes to develop and result in slow changes in the morphology of the AP. These changes are dependent on the nature and the formulation of the stimulus. The long-term behavior of the model depends on whether the charge applied through the stimulus is carried by an ion, the balance of which is explicitly included in the model. The choice of the ion selected to carry the stimulus charge also affects the outcome. We have chosen to simulate the model without the stimulus having a direct effect on the ion concentrations. The sensitivity of long-term ionic concentrations to stimulus formulation and frequency is a property of ionic models that attempt to preserve a detailed ionic balance (25a, 50a). When discussing results under periodic pacing, we have chosen to present model results after initial rapid transients have dissipated, after 12 s of pacing from rest. The AP morphology at this point is independent of the form and implementation of the stimulus and takes into account the kinetic adaptation of the currents at the specified pacing frequency. However, it does not include the effect of concentration changes occurring on a slower time scale in the model. This is also true of the experimental data on which the model is based. Most frequency-dependent data are collected with the assumption that preparations stabilize within seconds of the onset of stimulation. A detailed investigation of the long-term dynamics of the model and its dependence on the formulation of the stimulus is left for further study.


The authors thank Drs. Jianlin Feng, Normand Leblanc, Gui-Rong Li, Hui Sun, and Zhiguo Wang for help with experimental data and Jean Perrault and Charles Dupont for the use of computer equipment.


  • Address for reprint requests: M. Courtemanche, Research Center, Montreal Heart Institute, 5000 E. Belanger St., Montreal, PQ, Canada H1T 1C8.

  • This work was funded by grants from the Medical Research Council of Canada (S. Nattel) and the Natural Sciences and Engineering Research Council (M. Courtemanche). M. Courtemanche is supported by a Chercheur-Boursier Scholarship from the Fonds de Recherche en Santé du Québec.


Some fractional equations require evaluation of a limit to determine their values at membrane potentials for which their denominator is zero.

Differential Equations

Instantaneous equilibration of Ca2+ with buffers is assumed in all casesdVdt=(Iion+Ist)Cm Equation 18 Iion=INa+IK1+Ito+IKur+IKr+IKs+ICa,L +Ip,Ca+INaK+INaCa+Ib,Na+Ib,Ca Equation 19 dydt=yyτy(whereyis any gating variable) Equation 20 d[Na+]idt=3INaK3INaCaIb,NaINaFVi Equation 21 d[K+]idt=2INaKIK1ItoIKurIKrIKsIb,KFVi Equation 22 d[Ca2+]idt=B1B2 Equation 23 B1=2INaCaIp,CaICa,LIb,Ca2FVi +Vup(Iup,leakIup)+IrelVrelVi Equation 24 B2=1+[Trpn]maxKm,Trpn([Ca2+]i+Km,Trpn)2+[Cmdn]maxKm,Cmdn([Ca2+]i+Km,Cmdn)2 Equation 25 d[Ca2+]updt=IupIup,leakItrVrelVup Equation 26 d[Ca2+]reldt=(ItrIrel)1+[Csqn]maxKm,Csqn([Ca2+]rel+Km,Csqn)21 Equation 27

Equilibrium Potential

Ex=RTzFlog[X]o[X]i,forX=Na+,K+,Ca2+ Equation 28

Fast Na+ Current

INa=gNam3hj(VENa) Equation 29 αm=0.32V+47.131exp[0.1(V+47.13)]3.2,ifV=47.13 βm=0.08expV11 Equation 30 αh=0.135expV+806.80,ifV40 βh=3.56exp(0.079V)+3.1×105exp(0.35V)0.131+expV+10.6611.11,ifV40 Equation 31

Formula Equation 32

βj=0.1212exp(0.01052V)1+exp[0.1378(V+40.14)]0.3exp[2.535×107V]1+exp[0.1(V+32)],ifV40 Equation 33 τφ=(αφ+βφ)1,φ=αφτφ,forφ=m,h,j Equation 34

Time-Independent K+ Current

IK1=gK1(VEK)1+exp[0.07(V+80)] Equation 35

Transient Outward K+ Current

Ito=gtooa3oi(VEK) Equation 36 αo(a)=0.65expV+108.5+expV3059.01 βo(a)=0.652.5+expV+8217.01 Equation 37 τo(a)=[αo(a)+βo(a)]1/KQ10 oa()=1+expV+20.4717.541 Equation 38 αo(i)=18.53+expV+113.710.951 βo(i)=35.56+expV+1.267.441 Equation 39 τo(i)=[αo(i)+βo(i)]1/KQ10 oi()=1+expV+43.15.31 Equation 40

Ultrarapid Delayed Rectifier K+ Current

IKur=gKurua3ui(VEK) Equation 41 gKur=0.005+0.051+expV1513 Equation 42 αu(a)=0.65expV+108.5+expV3059.01 βu(a)=0.652.5+expV+8217.01 Equation 43 τu(a)=[αu(a)+βu(a)]1/KQ10 ua()=1+expV+30.39.61 Equation 44 αu(i)=21+expV185281 βu(i)=expV15816 Equation 45 τu(i)=[αu(i)+βu(i)]1/KQ10 ui()=1+expV99.4527.481 Equation 46

Rapid Delayed Outward Rectifier K+ Current

IKr=gKrxr(VEK)1+expV+1522.4 Equation 47 αx(r)=0.0003V+14.11expV+14.15 βx(r)=7.3898×105V3.3328expV3.33285.12371 Equation 48 τx(r)=[αx(r)+βx(r)]1 xr()=1+expV+14.16.51 Equation 49

Slow Delayed Outward Rectifier K+ Current

IKs=gKsxs2(VEK) Equation 50 αx(s)=4×105V19.91expV19.917 βx(s)=3.5×105V19.9expV19.991 Equation 51 τx(s)=12[αx(s)+βx(s)]1 xs()=1+expV19.912.71/2 Equation 52

L-Type Ca2+ Current

ICa,L=gCa,LdffCa(V65) Equation 53 τd=1expV+106.240.035(V+10)1+expV+106.24 d=1+expV+1081 Equation 54 τf=9{0.0197exp[0.03372(V+10)2]+0.02}1 f=1+expV+286.91 Equation 55 τf(Ca)=2,fCa()=1+[Ca2+]i0.000351 Equation 56

Na+-K+ Pump Current

INaK=INaK(max)fNaK11+{Km,Na(i)/[Na+]i}1.5[K+]o[K+]o+Km,K(o) Equation 57 fNaK=1+0.1245exp0.1FVRT+0.0365ςexpFVRT1 Equation 58 ς=17exp[Na+]o67.31 Equation 59

Na+/Ca2+ Exchanger Current

Formula Equation 60

Background Currents

Ib,Ca=gb,Ca(VECa) Equation 61 Ib,Na=gb,Na(VENa) Equation 62

Ca2+ Pump Current

Ip,Ca=Ip,Ca(max)[Ca2+]i0.0005+[Ca2+]i Equation 63

Ca2+ Release Current From JSR

Irel=krelu2vw([Ca2+]rel[Ca2+]i) Equation 64 τu=8.0,u=1+expFn3.4175×101313.67×10161 Equation 65 τv=1.91+2.091+expFn3.4175×101313.67×10161 v=11+expFn6.835×101413.67×10161 Equation 66 τw=6.01expV7.951+0.3expV7.95(V7.9) w=11+expV40171 Equation 67 Fn=1012VrelIrel5×1013F12ICa,L15INaCa Equation 68

Transfer Current From NSR to JSR

Itr=[Ca2+]up[Ca2+]relτtr Equation 69 τtr=180 Equation 70

Ca2+ Uptake Current by the NSR

Iup=Iup(max)1+(Kup/[Ca2+]i) Equation 71

Ca2+ Leak Current by the NSR

Iup,leak=[Ca2+]up[Ca2+]up(max)Iup(max) Equation 72

Ca2+ Buffers

[Ca2+]Cmdn=[Cmdn]max[Ca2+]i[Ca2+]i+Km,Cmdn Equation 73 [Ca2+]Trpn=[Trpn]max[Ca2+]i[Ca2+]i+Km,Trpn Equation 74 [Ca2+]Csqn=[Csqn]max[Ca2+]rel[Ca2+]rel+Km,Csqn Equation 75

Numerical Integration

At time step p, the updated value of a time-dependent variable is given byχ(p)=χ(p1)+Δtdχdt Equation 76for χ = V, any time-dependent ionic concentration, and

y(p)=y+[y(p1)y]expΔtτy Equation 77 fory=any gating variable.


  1. 1.
  2. 2.
  3. 3.
  4. 4.
  5. 5.
  6. 6.
  7. 7.
  8. 8.
  9. 9.
  10. 10.
  11. 11.
  12. 12.
  13. 13.
  14. 14.
  15. 15.
  16. 16.
  17. 17.
  18. 18.
  19. 19.
  20. 20.
  21. 21.
  22. 22.
  23. 23.
  24. 24.
  25. 25.
  26. 25a.
  27. 26.
  28. 27.
  29. 28.
  30. 29.
  31. 30.
  32. 31.
  33. 32.
  34. 33.
  35. 34.
  36. 35.
  37. 36.
  38. 37.
  39. 38.
  40. 39.
  41. 40.
  42. 41.
  43. 42.
  44. 43.
  45. 44.
  46. 45.
  47. 46.
  48. 47.
  49. 48.
  50. 49.
  51. 50.
  52. 50a.
  53. 51.
  54. 52.
  55. 53.
  56. 54.
  57. 55.
  58. 56.
  59. 57.
  60. 58.
  61. 59.
  62. 60.
  63. 61.
  64. 62.
View Abstract