INTRODUCTION

Top Abstract Introduction Results Discussion Appendix References

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 Ca²⁺ current (I_Ca,L) (35) and its inactivation by voltage and intracellular Ca²⁺ (26, 49), and the Na⁺/Ca²⁺ 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 the I_Ca,L and the Na⁺/Ca²⁺ exchanger, and the variability in experimentally observed AP morphology.

Glossary

V	Transmembrane potential
R	Gas constant
T	Temperature
F	Faraday constant
Cm	Membrane capacitance
AP	Human atrial action potential
max	Maximal AP upstroke velocity (dV/dtmax)
APA	AP amplitude (peak AP voltage minus diastolic voltage at onset of AP)
APO	AP overshoot (peak AP voltage >0 mV)
APD	AP duration
APD50	APD measured from AP onset to 50% of total APA repolarization
APD90	APD measured from AP onset to 90% of total APA repolarization
AF	Atrial fibrillation
I-V	Current-voltage
Q10	Temperature adjustment factor, k(T ) = k(T0)Q(TT0)/1010

SR

Sarcoplasmic reticulum

JSR

Junctional SR, SR release compartment

NSR

Network SR, SR uptake compartment

V_cell

Cell volume

V_i

Intracellular volume

V_up

SR uptake compartment volume

V_rel

SR release compartment volume

[X ]_o

Extracellular concentration of ion X

[X ]_i

Intracellular concentration of ion X

Cmdn

Calmodulin, sarcoplasmic Ca²⁺ buffer

Trpn

Troponin, sarcoplasmic Ca²⁺ buffer

Csqn

Calsequestrin, JSR Ca²⁺ buffer

[Ca²⁺]_rel

Ca²⁺ concentration in release compartment

[Ca²⁺]_up

Ca²⁺ concentration in uptake compartment

[Ca²⁺]_Cmdn

Ca²⁺-bound calmodulin concentration

[Ca²⁺]_Trpn

Ca²⁺-bound troponin concentration

[Ca²⁺]_Csqn

Ca²⁺-bound calsequestrin concentration

E_X

Equilibrium potential for ion X

I_ion

Total ionic current

I_st

Stimulus current

_x

Forward rate constant for gating variable x

_x

Backward rate constant for gating variable x

_x

Time constant for gating variable x

x

Steady-state relation for gating variable x

I_Na

Fast inward Na⁺ current

g_Na

Maximal I_Na conductance

m

Activation gating variable for I_Na

h

Fast inactivation gating variable for I_Na

j

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

Q₁₀-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 Ca²⁺ current

g_Ca,L

Maximal I_Ca,L conductance

d

Activation gating variable for I_Ca,L

f

Voltage-dependent inactivation gating variable for I_Ca,L

f_Ca

Ca²⁺-dependent inactivation gating variable for I_Ca,L

I_p,Ca

Sarcoplasmic Ca²⁺ pump current

I_p,Ca(max)

Maximal I_p,Ca

I_NaK

Na⁺-K⁺ pump current

I_NaK(max)

Maximal I_NaK

f_NaK

Voltage-dependence parameter for I_NaK

[Na⁺]_o-dependence parameter for I_NaK

K_m,Na(i)

[Na⁺]_i half-saturation constant for I_NaK

K_m,K(o)

[K⁺]_o half-saturation constant for I_NaK

I_NaCa

Na⁺/Ca²⁺ exchanger current

I_NaCa(max)

I_NaCa scaling factor

K_m,Na

[Na⁺]_o saturation constant for I_NaCa

K_m,Ca

[Ca²⁺]_o saturation constant for I_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 Ca²⁺ current

g_b,Ca

Maximal I_b,Ca conductance

I_rel

Ca²⁺ release current from the JSR

k_rel

Maximal Ca²⁺ release rate for I_rel

u

Activation gating variable for I_rel

v

Ca²⁺ flux-dependent inactivation gating variable for I_rel

w

Voltage-dependent inactivation gating variable for I_rel

F_n

Sarcoplasmic Ca²⁺ flux signal for I_rel

I_up

Ca²⁺ uptake current into the NSR

I_up(max)

Maximal Ca²⁺ uptake rate for I_up

[Ca²⁺]_up(max)

Maximal Ca²⁺ concentration in NSR

I_tr

Ca²⁺ transfer current from NSR to JSR

_tr

Ca²⁺ transfer time constant

I_up,leak

Ca²⁺ leak current from the NSR

[Cmdn]_max

Total calmodulin concentration in myoplasm

[Trpn]_max

Total troponin concentration in myoplasm

[Csqn]_max

Total calsequestrin concentration in JSR

K_m,Cmdn

Ca²⁺ half-saturation constant for calmodulin

K_m,Trpn

Ca²⁺ half-saturation constant for troponin

K_m,Csqn

Ca²⁺ half-saturation constant for calsequestrin

V_rest

Resting membrane potential

	MODEL DESCRIPTION

General

(1)

(2)

View larger version (22K): [in this window] [in a new window]   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)]. See Glossary for definitions.

A modified Euler method is used to integrate Eq. 1 with a fixed time step t = 0.005 ms (see APPENDIX, 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 APPENDIX.

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 by

(3)

Using 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 use g_Na = 7.8 nS/pF, which is ~1.3 times the temperature-adjusted value (with the assumption of Q₁₀ = 3) reported by Sakakibara et al. Figure 2 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-V curve 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.

View larger version (17K): [in this window] [in a new window]   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, hj, 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.

View larger version (22K): [in this window] [in a new window]   Fig. 3.   Current-voltage (I-V) relationships of major currents in model with corresponding experimental data. Model curves (solid lines, dashed line without symbols for INa) 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 experimental I-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) for ICa,L scaled by 0.38, Feng et al. (19) for IKur scaled by 0.81 and Ito scaled by 1.2, Wang et al. (55) for IKs scaled by 3.0 and IKr scaled by 0.38. Refer to original studies for exact experimental voltage-clamp protocols used. Experimental I-V curves for IKur and Ito 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 inhibit ICa,L (13, 28). Experimental I-V data for IK (55) were transformed to pA/pF with use of Cm = 100 pF before scaling was performed. No experimental data are shown for INa.

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 by

(4)

We 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_K1 from 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 decreased I_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. The I-V relationship of I_K1 is shown in Fig. 4.

View larger version (16K): [in this window] [in a new window]   Fig. 4.   I-V relationships for IK1 and INaK obtained from model. Sum of instantaneous time-independent membrane currents (Iv) is also shown. For INaK, 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 and I_Kur on the basis of data from our laboratory (20, 53). The currents are given by

(5)

(6)

Both 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 Q₁₀ = 2.2 was used to adjust their kinetics to 37°C (53). Current amplitudes were adjusted to match AP morphology. Figure 5 shows steady-state activation, inactivation, and corresponding time constants for I_Kur, along with experimental data from our laboratory (53). Figure 5 shows similar results for I_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 Cd²⁺, which was used experimentally to block I_Ca,L. The peak I-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 of I_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).

View larger version (28K): [in this window] [in a new window]   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.

View larger version (14K): [in this window] [in a new window]   Fig. 6.   IKur and Ito generated by model in response to voltage-clamp steps. Results are shown for combined currents (A), Ito alone (B), and IKur 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 of IKur + Ito, from Wang et al. (55), are included for comparison (dashed lines in A). Experimental recordings were scaled using Cm = 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 by

(7)

(8)

Both currents have only a single activation gate, although I_Ks activation uses a squared (x²) 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-V curves from Wang et al., are shown in Fig. 3. Figure 8 displays traces for I_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_Kr and I_Ks, respectively.

View larger version (13K): [in this window] [in a new window]   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).

View larger version (11K): [in this window] [in a new window]   Fig. 8.   IKr and IKs generated by model in response to voltage-clamp steps. Results are shown for combined currents (A), IKs alone (B), and IKr 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 of IKr and IKs, from Wang et al. (55), are included for comparison (dashed lines in B and C). With use of Cm = 100 pF to obtain units of pA/pF, experimental recordings were further scaled by a factor of 3.5 (IKs) and 0.47 (IKr) to match peak amplitude of model currents at +20 mV.

L-type Ca²⁺ current. The model implements the slow inward Ca²⁺ current on the basis of previous models and data from human atrial cells (35, 37, 40, 49). The current is given by

(9)

Our formulation follows that of Luo and Rudy (37) and includes voltage-dependent activation and inactivation gates, as well as a Ca²⁺-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 Ca²⁺-dependent inactivation is slightly modified from that of Luo and Rudy to better match the Ca²⁺ transient and reproduce the extent of Ca²⁺-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 Ca²⁺ storage and release), to better reproduce the time course of Ca²⁺-dependent inactivation of I_Ca,L.

View larger version (14K): [in this window] [in a new window]   Fig. 9.   Model representation of parameters describing ICa,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 in A).

View larger version (12K): [in this window] [in a new window]   Fig. 10.   A: ICa,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 by

(10)

The 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⁺/Ca²⁺ exchanger. The exchanger current follows the revised formulation of Luo and Rudy (37). The exchanger current is given by

(11)

View larger version (14K): [in this window] [in a new window]   Fig. 11.   I-V relationships for INaCa obtained from model. Current amplitude depends on Na+ and Ca2+ concentrations. Resting concentration shown in Table 2 is 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 Ca²⁺ and Na⁺ currents. The Ca²⁺ leak current is adjusted to maintain a stable [Ca²⁺]_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 by

(12)

(13)

Ca²⁺ pump current. A sarcolemmal Ca²⁺ pump is included in the model to maintain [Ca²⁺]_i at physiological levels. The pump current formulation (37, 40) is given by

(14)

Figure 4 shows the total time-independent I-V relationship for the model, which includes the background currents, I_K1, and the pump and exchanger currents.

SR Ca²⁺ storage and release. SR Ca²⁺ uptake and release are implemented using a two-compartment model, as reported by Luo and Rudy (37). Intracellular Ca²⁺ 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 [Ca²⁺]_i dependent and includes a leak from the uptake compartment back into the intracellular space. The main currents (in mM/ms) involved in SR Ca²⁺ handling are given by

(15)

(16)

(17)

View larger version (16K): [in this window] [in a new window]   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 and v are shown. C: steady-state inactivation curve and corresponding inactivation time constant for voltage-dependent gating variable w.

View larger version (16K): [in this window] [in a new window]   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 Ca²⁺ buffers. Ca²⁺ buffering within the cytoplasm is mediated by troponin and calmodulin. Ca²⁺ buffering in the release compartment is mediated by calsequestrin. In particular, calsequestrin is essential in providing sufficient Ca²⁺ stores for release within the limited JSR volume. Buffers are considered at equilibrium throughout our simulations, and incoming Ca²⁺ in the myoplasm and JSR compartments is instantaneously equilibrated into free and buffer-bound fractions (see APPENDIX 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).

	RESULTS

Top Abstract Introduction Results Discussion Appendix References

Resting Properties

View larger version (23K): [in this window] [in a new window]   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

View larger version (23K): [in this window] [in a new window]   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

View larger version (39K): [in this window] [in a new window]   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 for I_Ca,L, I_Kr, and I_Ks to quantify their role in rate-dependent AP changes. We first consider a comparison of the slow gating variables f, f_Ca, x_r, and x²_s 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. For I_Kr and I_Ks, we found that incomplete deactivation at 300 ms resulted in an increase of 0.05 for x_r and 0.004 for x²_s compared with their values at 1,000 ms. This compares with reductions of 0.2 for f and 0.16 for f_Ca. The maximal conductances of I_Ca,L and I_Ks are comparable, whereas the maximal conductance of I_Kr is four times smaller. Hence, the changes in the f and f_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_K amplitude may have significant effects during late repolarization, when I_K is the dominant repolarizing current.

To better ascertain the relative importance of I_K and 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 d and f gating variables), I_Kr (via the x_r gating variable), I_Ks (via the x_s gating variable), and [Ca²⁺]_i (via [Ca²⁺]_i, [Ca²⁺]_up, [Ca²⁺]_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), and I_Ca,L + I_K. We note that when the effect of either I_Kr or I_Ks alone was examined, the observed change in APD₉₀ was ~50% of the total reduction seen with both currents combined (I_K in Fig. 17). Also, the effect on APD₉₀ of [Ca²⁺]_i and I_Ca,L combined was similar to that of I_Ca,L alone, which is why the results for [Ca²⁺]_i were not included in Fig. 17. The limited role of [Ca²⁺]_i at fast pacing rates may reflect the opposing effects of a higher resting [Ca²⁺]_i and a lower Ca²⁺ 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 APD₉₀. Third, the effect of I_Ca,L and I_K combined is synergistic, lowering the plateau potential and accelerating late repolarization relative to the effect of I_Ca,L alone. 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 that I_Ca,L and I_K act together to produce the observed rate-dependent AP shortening.

View larger version (22K): [in this window] [in a new window]   Fig. 17.   Role of ICa,L and IK 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 ICa,L alone, IK (IKr + IKs) alone, and ICa,L + IK are shown. See RESULTS for a detailed description of numerical simulation procedure. Combination of ICa,L and IK 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.

I_Ca,L and Rate Dependence

View larger version (12K): [in this window] [in a new window]   Fig. 18.   Role of ICa,L in rate dependence of action potential. A: model action potentials recorded after 12 s of pacing at various BCLs with 90% decrease in ICa,L conductance. 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 block ICa,L. No rate dependence is observed.

Block of I_NaCa and Role in APD

View larger version (21K): [in this window] [in a new window]   Fig. 19.   Model action potentials during stimulation at 1 Hz under control conditions and with a 90% decrease in INaCa maximal amplitude. Block of INaCa 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

View larger version (32K): [in this window] [in a new window]   Fig. 20.   Effect of varying ICa,L (A), IKur (B), and Ito (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).

	DISCUSSION

Top Abstract Introduction Results Discussion Appendix References

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_1/2) well negative of the value observed in other mammalian cells and insufficient to generate the required fast inward current during the AP.

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 that I_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 of I_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 that I_Kur had activation kinetics similar to those of I_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 Ca²⁺ 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 Ca²⁺, are only partially understood (49). This has made modeling the current more difficult and highly dependent on an accurate formulation of intracellular Ca²⁺ homeostasis and its control by the SR system (see SR Ca²⁺ handling and buffering). Our formulation follows that of Luo and Rudy (37) and includes voltage-dependent activation and inactivation gates as well as a Ca²⁺-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 Ca²⁺-driven inactivation of I_Ca,L remain, 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,L window currents observed on reuptake of Ca²⁺ during long APs. Another example can be found in the data of Sun et al., where an apparent second Ca²⁺-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,L inactivation arises as a result of a more complex time course of the Ca²⁺ transient. This is in agreement with data from Hatem et al. (26) that point to rapid and slow phases during the rise of the Ca²⁺ transient in human atrial myocytes. These factors are related to limitations in our modeling of the SR Ca²⁺ handling, which represents a simplified formulation of a temporally and spatially complex intracellular process (see SR Ca²⁺ handling and buffering and Limitations).

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 for I_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_Ks reached 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_Ks activation 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 than I_Kr at the lower plateau potentials typical of the human atrial AP.

SR Ca²⁺ handling and buffering. SR Ca²⁺ handling and buffering involve uptake, storage, and release of Ca²⁺ 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 Ca²⁺ from the JSR compartment depends on the flux of Ca²⁺ through the sarcolemmal and SR release channels. This formulation aims to reflect the close coupling of L-type Ca²⁺ channels in the sarcolemma and JSR Ca²⁺ release channel (34). For closely juxtaposed channels, one can expect that true local Ca²⁺ concentration in the vicinity of the release channel is more accurately reflected by local Ca²⁺ flux than by bulk cellular Ca²⁺ concentration. The formulation is also adapted to afford a greater sensitivity of the release mechanism to flux through the L-type Ca²⁺ channel, as opposed to Ca²⁺ entry through the Na⁺/Ca²⁺ exchanger (22). Release is triggered when the computed Ca²⁺ 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.

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 of I_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 model I-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 Ca²⁺-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 outward I_K, resulting in slow repolarization. Slow late repolarization is a balance between deactivating I_K, time-independent I_K1, and opposing inward current generated by the Na⁺/Ca²⁺ exchanger during the late phase of the AP.

Na⁺/Ca²⁺ exchanger. The Na⁺/Ca²⁺ 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 APD₉₀ (30% in Fig. 19). Because of the slow repolarization rate and small current amplitudes during late repolarization, APD₉₀ 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 APD₉₀ 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 APD₅₀ than at APD₉₀. This is seen in the data of Le Grand et al., even though their APD₉₀ 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 of I_K in APs with shorter, more negative plateaus that opposes the rate-dependent increase in I_K availability. This tends to reduce rate adaptation of the model AP measured using APD₉₀ and is a direct consequence of the sensitivity of late AP repolarization to small changes in I_K, I_K1, and I_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 in I_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_K from 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).

Clinical Relevance

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

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-V curves, 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 or I_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 Ca²⁺ dynamics, Ca²⁺ 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 Ca²⁺-handling scheme proposed by Luo and Rudy (37). This produces a realistic intracellular Ca²⁺ transient in response to suprathreshold Ca²⁺ flux into the cell. However, recent experimental results suggest that Ca²⁺-dependent I_Ca,L inactivation may involve complex local phenomena requiring an additional level of model complexity (49), beyond an accurate representation of the average Ca²⁺ 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 Ca²⁺-dependent I_NaCa.

Our formulation of the delayed rectifier current represents a compromise between a simple current representation and available experimental evidence that I_Kr and I_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 as I_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.