INTRODUCTION

TOP ABSTRACT INTRODUCTION THEORY AND METHODS RESULTS AND DISCUSSION APPENDIX REFERENCES

PACEMAKER ACTIVITY of the sinoatrial (SA) node is well known as the initiation of the spontaneous heart beat. A large body of information on ionic current systems underlying the SA node pacemaker activity has been obtained from recent single cell patch-clamp experiments. On the basis of single cell patch-clamp data, several mathematical models describing the pacemaker activity of a single rabbit SA node cell have been developed in the past decade.

In 1991, Wilders et al. (124) first developed a single SA node cell model based on single cell patch-clamp data to reproduce the electrical behavior of a single rabbit SA node cell quantitatively. Later, Demir et al. (17) proposed a more detailed model for transitional-type cells (rather than for primary pacemaker cells) including intracellular Ca²⁺ buffering (by calmodulin, troponin, and calsequestrin) and new formulations for Ca²⁺ handling by the sarcoplasmic reticulum (SR) based on microanatomic data. Dokos et al. (24) developed a different SA node model incorporating new formulations of the Na⁺/Ca²⁺ exchanger current (I_NaCa) and muscarinic K⁺ channel current (I_K,ACh) as a background K⁺ current. However, these previous models have several difficulties, as follows. First, the formulations of voltage-dependent inactivation and recovery of the L-type Ca²⁺ channel current (I_Ca,L) were not based on experimental data. Second, the Ca²⁺-dependent inactivation of I_Ca,L, known as an essential property of the L-type Ca²⁺ channel in the rabbit SA node (10, 67), was not formulated [or not directly dependent on the intracellular Ca²⁺ concentration ([Ca²⁺]_i), but a function of the inactivation gating variable f_L in Wilders et al. (124)]. Although Dokos et al. (24) incorporated a second inactivation gating variable to represent the [Ca²⁺]_i-dependent inactivation process, they assumed an extraordinarily large maximum I_Ca,L conductance (g_Ca,L) and a very high affinity for Ca²⁺ binding to the inactivation site of L-type Ca²⁺ channels; such a large g_Ca,L with a high affinity for Ca²⁺ binding is rather unlikely and has not been demonstrated. Third, these models do not include the slow component of the delayed rectifier K⁺ current (I_Ks), 4-aminopyridine (4-AP)-sensitive currents, or sustained inward current (I_st), whereas these currents are known to be present in primary pacemaker cells. Fourth, intracellular Ca²⁺ transients in their models (2.5~10 µM at the peak) were much higher than those in atrial or ventricular myocytes (~1 µM) (e.g., Refs. 106 and 110; see also Refs. 68 and 71), probably too high for SA node cells. Finally, SR volumes or Ca²⁺ concentrations in the SR are comparable to or higher than those experimentally determined for ventricular myocytes (see Refs. 17, 63, and 107), probably too large for SA node cells. Recently, Zhang et al. (130) developed separate models for central and peripheral SA node cells based on recent experiments in which the regional differences in action potential (AP) parameters, ionic current densities, and pharmacological responses (i.e., electrophysiological effects of various current blockers) between central and peripheral SA node cells have been studied. Their central and peripheral models provide a theoretical basis for regional differences in AP parameters and are superior to the previous models in that they incorporate novel current systems such as I_Ks and 4-AP-sensitive currents. However, their models have at least two apparent difficulties. First, intracellular ion concentrations were assumed to be constant in their models (the net ion flux of Na⁺, K⁺, or Ca²⁺ during an AP cycle is not zero), whereas intracellular Ca²⁺ dynamics and changes in intracellular Na⁺ ([Na⁺]_i) and K⁺ concentrations ([K⁺]_i) are known to exert substantial effects on pacemaker activity (e.g., see Refs. 67, 98, and 122). Second, their models lack some sarcolemmal currents, such as I_st and I_K,ACh, known to play important roles in regulating the pacemaker activity of rabbit SA node cells. Thus the previous SA node models all have several drawbacks or serious disadvantages; a satisfactory single cell model is not available.

The aim of this study was to develop an improved mathematical model of a single "primary pacemaker cell" of the rabbit SA node that is more suitable than the previous models for investigating the dynamical mechanisms of pacemaker generation. In addition to the drawbacks stated above, the previous single cell models exhibited different bifurcation structures, i.e., different ways of abolishing automaticity (see Ref. 30) from those of real SA node cells during inhibition of I_Ca,L or the rapid component of the delayed rectifier K⁺ current (I_Kr) not suitable for a study on the mechanisms of pacemaker generation (for details, see RESULTS AND DISCUSSION). Thus an improved model cell should have the same bifurcation structures as real SA node cells have as well as the capability of reproducing experimental data more accurately than the previous models. Such a model would also serve as a base model for developing more sophisticated models. On the basis of recent experimental findings, we were able to update the previous models in several ways: 1) I_st, not included in the previous models, has been incorporated; 2) voltage- and Ca²⁺-dependent inactivation kinetics of I_Ca,L have been reformulated; 3) expressions for activation kinetics of I_Kr have renewed; 4) revised kinetic formulas for 4-AP-sensitive currents have been incorporated; 5) voltage- and concentration-dependent kinetics of the Na⁺-K⁺ pump current (I_NaK) have been reformulated; and 6) the subsarcolemmal space as a diffusion barrier for intracellular Ca²⁺ has been incorporated.

To validate our model, we first compared the simulated dynamics of the model (spontaneous APs and ionic currents during pacemaker activity) with experimental data from the rabbit SA node as well as with those of the previous SA node models. We further validated our model by simulating the effects of modulations of sarcolemmal ionic currents or intracellular Ca²⁺ dynamics on pacemaker activity. The experimental findings for verification include 1) whole cell voltage-clamp data for I_Ca,L, I_Kr, and I_st; 2) waveshapes of spontaneous APs and ionic currents (such as I_Ca,L and I_Kr) as observed during AP-clamp recordings; 3) modulations and cessation of pacemaker activity by applications of I_Ca,L or I_Kr blockers; 4) modifications of AP waveforms (changes in AP parameters) by blocking the T-type Ca²⁺ channel current (I_Ca,T), hyperpolarization-activated current (I_h), or 4-AP-sensitive currents; 5) effects of blocking SR Ca²⁺ release (by ryanodine) on pacemaker frequency; and 6) negative chronotropic effects of Ca²⁺ buffers. Our model could reproduce these experimental data more accurately than the previous SA node models. In this study, we particularly focused on the bifurcation structures during applications of I_Ca,L or I_Kr blockers and also the differential effects of BAPTA and EGTA on pacemaker frequency. During inhibition of I_Ca,L or I_Kr, our model exhibited essentially the same bifurcation structures as observed in real SA node cells, whereas previous models did not; only our model could simulate the differential responses of SA node cells to BAPTA and EGTA. Detailed comparisons of our model with the previous models as well as experimental data are addressed in THEORY AND METHODS (on modeling) and RESULTS AND DISCUSSION (on simulated results).

Glossary

General

4-AP	4-Aminopyridine
AP	Action potential
APA	AP amplitude
APD50	AP duration at 50% repolarization
CL	Cycle length
F	Faraday constant
I	Current
MDP	Maximum diastolic potential
POP	Peak overshoot potential
R	Universal gas constant
SA	Sinoatrial
SR	Sarcoplasmic reticulum
t	Time
T	Absolute temperature
V	Membrane potential (in mV)

Ionic Channel Currents

Cell Geometry

Ionic Concentrations

Equilibrium (Reversal) Potentials

Sarcolemmal Ionic Currents

L	Activation gating variable for ICa,L
dL,	Steady-state dL
dT	Activation gating variable for ICa,T
dT,	Steady-state dT
fCa	Ca2+-dependent inactivation gating variable for ICa,L
fCa,	Steady-state fCa
fL	Voltage-dependent inactivation gating variable for ICa,L
fL,	Steady-state fL
fT	Inactivation gating variable for ICa,T
fT,	Steady-state fT
FNa	Fraction of slow inactivation of INa
gb,Na	Background Na+ conductance
gCa,L	Maximum ICa,L conductance
gCa,T	Maximum ICa,T conductance
gh	Maximum Ih conductance
gh,K	K+ current component of gh
gh,Na	Na+ current component of gh
gK,ACh	Scaling factor for IK,ACh
gKr	Maximum IKr conductance
gKs	Maximum IKs conductance
gNa	Maximum INa conductance
gst	Maximum Ist conductance
gsus	Maximum Isus conductance
gto	Maximum Ito conductance
h	Inactivation gating variable for INa
h	Steady-state h
hF	Fast inactivation gating variable for INa
hS	Slow inactivation gating variable for INa
jCa,sm	Net Ca2+ flux through the sarcolemmal membrane
kNaCa	Scaling factor for INaCa
K1ni	Dissociation constant for intracellular Na+ binding to first site on INaCa transporter
K2ni	Dissociation constant for intracellular Na+ binding to second site on INaCa transporter
K3ni	Dissociation constant for intracellular Na+ binding to third site on INaCa transporter
K1no	Dissociation constant for extracellular Na+ binding to first site on INaCa transporter
K2no	Dissociation constant for extracellular Na+ binding to second site on INaCa transporter
K3no	Dissociation constant for extracellular Na+ binding to third site on INaCa transporter
Kci	Dissociation constant for intracellular Ca2+ binding to INaCa transporter
Kco	Dissociation constant for extracellular Ca2+ binding to INaCa transporter
Kcni	Dissociation constant for intracellular Na+ and Ca2+ simultaneous binding to INaCa transporter
KmfCa	Dissociation constant for Ca2+-dependent inactivation of ICa,L
KmKp	Half-maximal [K+]o for INaK
KmNap	Half-maximal [Na+]i for INaK
m	Activation gating variable for INa
m	Steady-state m
n	Activation gating variable for IKs
n	Steady-state n
pa	Activation gating variable for IKr
pa,	Steady-state pa
paF	Fast activation gating variable for IKr
paS	Slow activation gating variable for IKr
pi	Inactivation gating variable for IKr
pi,	Steady-state pi
q	Inactivation gating variable for Ito
q	Steady-state q
qa	Activation gating variable for Ist
qa,	Steady-state qa
qi	Inactivation gating variable for Ist
qi,	Steady-state qi
Qci	Fractional charge movement during intracellular Ca2+ occlusion reaction of INaCa transporter
Qco	Fractional charge movement during extracellular Ca2+ occlusion reaction of INaCa transporter
Qn	Fractional charge movement during Na+ occlusion reactions of INaCa transporter
r	Activation gating variable for Ito and Isus
r	Steady-state r
x	Gating variable
x	Steady-state value of x
y	Activation gating variable for Ih
y	Steady-state value of y
dL	Opening rate constant of dL
fCa	Ca2+ dissociation rate constant for ICa,L
qa	Opening rate constant of qa
qi	Opening rate constant of qi
n	Opening rate constant of n
dL	Closing rate constant of dL
fCa	Ca2+ association rate constant for ICa,L
n	Closing rate constant of n
qa	Closing rate constant of qa
qi	Closing rate constant of qi
dL	Time constant of dL
dT	Time constant of dT
fCa	Time constant of fCa
fL	Time constant of fL
fT	Time constant of fT
hF	Time constant of hF
hS	Time constant of hS
m	Time constant of m
n	Time constant of n
paF	Time constant of paF
paS	Time constant of paS
q	Time constant of q
qa	Time constant of qa
qi	Time constant of qi
r	Time constant of r
x	Time constant for a gating variable x
y	Time constant of y

Ca²+ Diffusion

jCa,dif	Ca2+ diffusion flux from subspace to myoplasm
dif,Ca	Time constant of Ca2+ diffusion from the subspace to myoplasm

SR Function

jrel	Ca2+ release flux from the junctional SR to subspace
jtr	Ca2+ transfer flux from the network to junctional SR
jup	Ca2+ uptake flux from the myoplasm to network SR
Krel	Half-maximal [Ca2+]sub for Ca2+ release from the junctional SR
Kup	Half-maximal [Ca2+]i for Ca2+ uptake by the Ca2+ pump in the network SR
Prel	Rate constant for Ca2+ release from the junctional SR
Pup	Rate constant for Ca2+ uptake by the Ca2+ pump in the network SR
tr	Time constant for Ca2+ transfer from the network to junctional SR

Ca²+ Buffering

	THEORY AND METHODS

TOP ABSTRACT INTRODUCTION THEORY AND METHODS RESULTS AND DISCUSSION APPENDIX REFERENCES

Model Development

Geometrical considerations. We assumed our model cell to be a 70-µm-long by 8-µm-diameter cylinder (i.e., a "spindle-shaped" cell with a length of 70 µm and a mean width of 8 µm) and set the cell volume and membrane capacitance to 3.5 pl and 32 pF, respectively. The cell volume of 3.5 pl was nearly identical to that used by Demir et al. (17) but smaller than the value of 5 pl estimated by Denyer and Brown (19) for spindle-shaped isolated SA node cells and used for the models of Wilders et al. (124) and Dokos et al. (24). Primary pacemaker cells are thought to be smaller than transitional or peripheral cells (see Refs. 44 and 130); thus we chose the smaller value of 3.5 pl for our primary pacemaker cell model. The membrane capacitance of 32 pF is derived from Wilders et al. (124) and Dokos et al. (24), smaller than the value of 55 pF used by Demir et al. (17) for their transitional cell model and larger than the value of 20 pF for the central model of Zhang et al. (130). Microanatomic data such as SR volumes used in this study are from Demir et al. (17). The effective intracellular volume of a SA node cell wherein free Ca²⁺ is available to enter into reaction was set to 46% of the cell volume (1.6 pl), because there are various intracellular structures and organelles, as summarized by Demir et al. (17).

9

View larger version (15K): [in this window] [in a new window]   Fig. 1.   Schematic diagrams depicting a cross section of the model cell (A) and intracellular fluid compartments for Ca2+ (B). R and L represent the radius of the cell (4 µm) and the thickness (depth) of the subsarcolemmal space (20 nm), respectively. Volumes of the cell and each compartment are given in the text. The Ca2+ release from junctional SR (JSR) and Ca2+ uptake to the network SR (NSR) were assumed to be strictly into the restricted subspace and from the bulk myoplasm, respectively. Definitions of the symbols for the Ca2+ concentration in each compartment and Ca2+ flux between compartments are given in the Glossary.

Descriptions of membrane currents. The mathematical expressions used for the membrane current system are essentially the same as those formulated previously (17, 24, 124, 130) with modifications according to newly reported experimental data. The complete model for the normal pacemaking includes 13 membrane current components. The differential equation for the membrane potential (V) is

(1)

where I_Ca,L and I_Ca,T represent the L-type and T-type Ca²⁺ channel currents, respectively. The rapid and slow components of the delayed rectifier K⁺ current are denoted as I_Kr and I_Ks, respectively. The membrane current system also includes the transient (I_to) and sustained (I_sus) components of 4-AP-sensitive currents, hyperpolarization-activated current (I_h), sustained inward current (I_st), Na⁺ channel current (I_Na), background Na⁺ current (I_b,Na), muscarinic K⁺ channel current (I_K,ACh), Na⁺-K⁺ pump current (I_NaK), and Na⁺/Ca²⁺ exchanger current (I_NaCa) charging the membrane capacitance (C_m). For parameter adjustments, the regional differences of current densities in the SA node were taken into account (44, 130).

(2)

L-type Ca²+ channel current. The kinetics of I_Ca,L are described with activation (d_L), voltage-dependent inactivation (f_L), and Ca²⁺-dependent inactivation (f_Ca) gating variables. The inactivation of I_Ca,L consists of two exponential terms: rapid and slow components. The rapid component is mediated by the intracellular Ca²⁺-dependent inactivation with a time constant ranging from 10 to 30 ms for rabbit SA node cells, whereas the slower component reflects the voltage-dependent inactivation with a time constant of 30~70 ms (5, 34, 72, 81, 102). To describe the inactivation process, we adopted a simple model in which the Ca²⁺-dependent inactivation process is independent of the voltage-dependent one (see Refs. 35 and 105); the inactivation process was described by the two Hodgkin-Huxley type gating variables f_L and f_Ca.

L,

L,

1

1

View larger version (29K): [in this window] [in a new window]   Fig. 2.   Kinetics of ICa,L. Top: voltage dependence of steady-state probabilities (dL and fL) and time constants for ICa,L activation (dL) and inactivation (fL). The equations used for the present model are shown as thick lines (K). For comparison, those used for previous models are also shown as thin lines: W, Wilders et al. (124); De, Demir et al. (17); Do, Dokos et al. (24); and Z, Zhang et al. (130). The experimental data for fL are from Kawano and Hiraoka (54) (), Hagiwara et al. (37) (), and Nakayama et al. (81) (). Bottom: computed voltage-clamp records for ICa,L (left) and the peak ICa,L-V relationship (right). Currents were evoked by 100-ms step pulses to test potentials ranging from 30 to +50 mV (in 10-mV increments). The holding potential was 40 mV. Simulating the whole cell perforated-patch recording, [Ca2+]i was not fixed (intracellular Ca2+ was not buffered by EGTA or BAPTA), whereas [Na+]i and [K+]i were fixed at 10 and 140 mM, respectively. The experimental data for the peak ICa,L-V relationship (peak currents normalized to the maximum value at 0 mV) are from Honjo et al. (44) (), Hagiwara et al. (37) (), Vinogradova et al. (122) (), and Verheijck et al. (120) (). See the Glossary for definitions of the abbreviations.

T-type Ca²+ channel current. Contributions of I_Ca,T to pacemaker depolarization in the previous models are different: Wilders et al. (124) and Dokos et al. (24) assumed small contribution of I_Ca,T according to the experimental report of Hagiwara et al. (37), whereas Demir et al. (17) assumed much larger I_Ca,T based on the data from Doerr et al. (23). With the use of the I_Ca,T expressions of Demir et al. (17), our standard model could reproduce the relatively large effects of blocking I_Ca,T on CL as observed experimentally (23, 37, 103); in contrast, the spontaneous AP of our standard model cell was little affected by incorporating the I_Ca,T expressions of Wilders et al. (124), with CL increasing only by 2.1% on eliminating I_Ca,T. Therefore, we adopted the expressions of Demir et al. (17) rather than those of Wilders et al. (124) or Dokos et al. (24) for the steady-state probabilities and time constants of the activation gating variable d_T and inactivation gating variable f_T. The conductance property of I_Ca,T was formulated using a linear voltage relation, with the reversal potential E_Ca,T fixed at +45 mV and the maximum conductance g_Ca,T set to 0.458 nS/pF, according to Demir et al. (17).

Delayed rectifier K+ currents. Recent studies in rabbit, guinea pig, and human cardiac myocytes have identified two types of delayed rectifier K⁺ currents: 1) the rapidly activating component I_Kr, exhibiting strong inward rectification, and 2) the slowly activating component I_Ks, exhibiting only weak rectification (2, 34, 64, 66, 100, 101, 123). Whereas I_K in the rabbit SA node appears to predominantly reflect only I_Kr-type behavior (2, 108), the I_Ks-type component has also been identified in isolated rabbit SA node cells (34, 48, 64). For the standard pacemaker model, therefore, we assumed both I_Kr- and I_Ks-type components of I_K.

a,

i,

View larger version (32K): [in this window] [in a new window]   Fig. 3.   Kinetics of IKr. Top: voltage dependence of steady-state probabilities for IKr activation (pa) and inactivation (pi) and activation time constants (pa). Thick lines represent the present model (K); those for previous models (W, De, Do, and Z) are also shown as thin lines for comparison.  and , Experimental data from Ono and Ito (90). Bottom: computed voltage-clamp records for IKr (left) and amplitudes of IKr measured at the end of test pulses as a function of the test potential (right). Currents were elicited by 1-s step pulses from a holding potential of 60 mV to test potentials ranging from 50 to +40 mV (in 10-mV increments). , Experimental data from Ono and Ito (90). See the Glossary for definitions of the abbreviations.

4-AP-sensitive K+ currents. Recent studies for rabbit SA node cells have identified the transient and sustained outward currents sensitive to 4-AP (46, 65, 119). The models of Zhang et al. (130) incorporated these currents, whereas most previous SA node models did not. Therefore, we incorporated the transient (I_to) and sustained (I_sus) components of 4-AP-sensitive currents into our model. According to Zhang et al. (130), we treated the two 4-AP-sensitive components as separated currents and used the same activation variable r for both I_to and I_sus.

Hyperpolarization-activated current. It is difficult to quantify I_h and its participation in the diastolic depolarization, in part due to the large variation in its threshold of activation ranging from 30 to 70 mV, as well as the variation in current density (18, 21, 73, 81, 118). In previous models, the kinetic data from van Ginneken and Giles (118) were used to formulate the gating kinetics of I_h; however, the contributions of I_h to pacemaker depolarization in the model cells are quite different (refer to Table 3). For our standard model, we adopted the formulation of Wilders et al. (124) to reproduce the relatively large effects of blocking I_h on CL as observed in experiments. The data of van Ginneken and Giles (118) were collected at 30~33°C; following Demir et al. (17), activation time constants were corrected for 37°C by the use of a Q₁₀ = 2.3. According to van Ginneken and Giles (118), the maximum conductance and reversal potential of I_h were assumed to be 0.375 nS/pF and 24 mV, respectively. These values were achieved by setting g_h,K = 7.4 nS and g_h,Na = 4.6 nS for a model cell, as used in Wilders et al. (124).

Sustained inward current. Guo et al. (32) reported a novel pacemaker current activated within the range of pacemaker depolarization in the rabbit SA node (see also Ref. 33). This current, named the sustained inward current (I_st), is carried by Na⁺ under physiological conditions, blocked by both organic and inorganic Ca²⁺ channel blockers as well as by external Ca²⁺ and Mg²⁺, and enhanced by isoprenaline or a Ca²⁺ agonist, BAY K 8644. These biophysical and pharmacological characteristics are compatible with those of the monovalent cation conductance of the L-type Ca²⁺ channel; thus Guo et al. (32) concluded that I_st is generated by a novel subtype of the L-type Ca²⁺ channel. I_st has also been recorded in SA node cells of other animal species, such as guinea pigs and rats (31, 77, 78, 109). Because I_st was observed only in spontaneously beating SA node cells but was absent in quiescent cells, I_st may play an essential role in the generation of intrinsic cardiac automaticity (32, 78). Therefore, we incorporated I_st in our standard SA node model for a primary pacemaker cell, whereas previous SA node models did not include I_st.

View larger version (26K): [in this window] [in a new window]   Fig. 4.   Kinetics of Ist. Top: voltage dependence of steady-state probabilities (qa and qi) and time constants for Ist activation (qa) and inactivation (qi). Thick lines, present model; thin lines, qa and qi formulated by Shinagawa et al. (109) for the rat SA node Ist.  and , Experimental values approximated from the data of Guo et al. (32). Bottom: computed voltage-clamp records for Ist (left) and the peak Ist-V relation (right). Currents were evoked by 500-ms depolarizing test pulses ranging from 70 to +50 mV (in 10-mV increments). The holding potential was 80 mV. See the Glossary for definitions of the abbreviations.

Na+ channel current. Although most of the previous SA node models incorporated the fast Na⁺ channel current (I_Na), the contribution of I_Na to model cell dynamics is relatively small: in the models of Wilders et al. (124) and Dokos et al. (24), eliminating I_Na yielded only a 0.7~0.9% increase in CL, whereas it yielded a 8.4% increase in the model of Demir et al. (17) for transitional cells. In a preliminary study, we incorporated I_Na into our model to assess the possible contribution of I_Na to the pacemaker activity of our model cell. The kinetics of I_Na were reformulated from the recent experimental data: 1) steady-state activation and inactivation curves are based on the data from Muramatsu et al. (80) and Baruscotti et al. (3); 2) for the time constant of activation, we used the formulation of Zhang et al. (130); and 3) the inactivation of I_Na was described by a weighted sum of two gating variables, h_F and h_S, with the fraction of slow inactivation (F_Na) being expressed as a function of V (see Ref. 130). The formulation for the inactivation time constant is based on data from Muramatsu et al. (80). According to Lindblad et al. (68), a Q₁₀ of 1.7 was used to correct the experimental data for 37°C. Incorporating I_Na with a conductance of 1.8 pS/pF, the mean experimental value for transitional or peripheral cells (44), decreased CL only by 5.6% (from 307.4 to 290.2 ms), indicating that the contribution of I_Na to pacemaking is relatively small in our model cell as well.

Na+-dependent background current and muscarinic K⁺ channel current. Our model includes two background current components: 1) Na⁺-dependent background current (I_b,Na), reported by Hagiwara et al. (38), in which I_b,Na was measured as 0.73 ± 0.21 pA/pF at 50 mV; and 2) muscarinic K⁺ channel current (I_K,ACh), reported by Ito et al. (49), who attributed the entire background K⁺ component in rabbit SA node cells to the spontaneous opening of muscarinic K⁺ channels in the absence of an ACh agonist and characterized I_K,ACh as exhibiting inward rectification with an amplitude of 0.33 ± 0.28 pA/pF at 50 mV under physiological conditions. Following their experimental data, we modeled I_b,Na with an ohmic I-V relationship and defined I_K,ACh as a background K⁺ current component exhibiting inward rectification. The I_b,Na conductance g_b,Na was set to 5.4 pS/pF, yielding a current density of 0.65 pA/pF at 50 mV. The g_b,Na value of 5.4 pS/pF is comparable to that in the previous SA node models (2.9~7.5 pS/pF). To model the inward rectifying behavior for I_K,ACh, we chose the formula of Dokos et al. (24) adopted from Egan and Noble (26). The standard value of I_K,ACh amplitude used in this study was 0.23 pA/pF at 50 mV, smaller than the value of 0.46 pA/pF in Dokos et al. (24).

Na+-K⁺ pump current. In previous SA node models, the formulations of I_NaK were not based on experimental data from SA node cells. The novel I_NaK formulation used in this study is based on the recent experimental work of Sakai et al. (99) for rabbit SA node cells. Parameter values were determined from the voltage- and concentration-dependent data provided by Sakai et al. (99): K_m values for Na⁺ and K⁺ binding were set to 14.0 and 1.4 mM, respectively. The voltage dependence of I_NaK in our model is steeper than that in previous models (see Fig. 17), consistent with experimental data from the rabbit SA node (99). The magnitude of I_NaK is linked to the influx of Na⁺ through I_b,Na (and I_st) and efflux of K⁺ through K⁺ channels. According to the data of Sakai et al. (99), we set the maximum attainable current (I_NaKmax) to 3.6 pA/pF, a value high enough to maintain [Na⁺]_i of <10 mM and [K⁺]_i of ~140 mM during long runs of the normal pacemaker activity.

Na+/Ca²⁺ exchange current. Most existing models of SA node activity have relied on the hypothetical I_NaCa formulation utilized by DiFrancesco and Noble (22), describing the Na⁺/Ca²⁺ exchanger as a simultaneous transfer of the ions with an exponential voltage dependency. In contrast, Dokos et al. (24) utilized a more accurate model of consecutive translocation, i.e., the "E4" model originally formulated by Matsuoka and Hilgemann (74), which correctly reproduces the saturation characteristics of I_NaCa at large [Ca²⁺]_i and for negative potential-low [Ca²⁺]_i conditions. To describe I_NaCa kinetics, we adopted the formulation of Dokos et al. (24), a modified version of the Matsuoka-Hilgemann "E4" model. The scaling parameter k_NaCa was set to 125 pA/pF, the same value as determined by Dokos et al. (24) from the study of Hagiwara and Irisawa (36), who reported an I_NaCa density of 1 pA/pF at 40 mV and [Ca²⁺]_i = 0.5 µM. A k_NaCa value of 125 pA/pF yielded an I_NaCa of 0.83 pA/pF at 40 mV and [Ca²⁺]_i = 0.5 µM, giving a diastolic free Ca²⁺ level in the presumed physiological range of <0.3 µM.

SR functions. The SR was modeled as consisting of two compartments: the Ca²⁺ uptake store (network SR) and release store (junctional SR), as shown in Fig. 1B. Owing to the lack of available data for updating the kinetic formulation of Ca²⁺ uptake and release by the SR in SA node cells, we utilized simple formulas similar to those of DiFrancesco and Noble (22), which have been incorporated in previous SA node models (24, 85, 124). The model of Demir et al. (17) utilized more complex SR kinetics based largely on the modeling study by Hilgemann and Noble (42), including the internal SR Ca²⁺ buffering by calsequestrin. Their formulation is, however, overly complex, given the complete lack of associated data in SA node cells; in their model, Ca²⁺ concentrations in the SR are unreasonably high compared with those in other models, much higher than experimental values for ventricular myocytes (see Refs. 63 and 107). We therefore opted for simpler kinetics, which are able to reproduce the essential features of both the uptake and release processes. According to Demir et al. (17), the fractional volumes of Ca²⁺ uptake and release stores in our model cell were taken to be 1.16 and 0.12%, respectively, of the cell volume, corresponding to one-third of the values given for ventricular myocytes.

Ca²+ buffering. Our model includes the dynamics of three Ca²⁺ buffers: calmodulin, troponin, and calsequestrin. The amounts of calmodulin and troponin available for Ca²⁺ binding, and the rate constants for Ca²⁺ binding to and release from these buffers, were adopted from Demir et al. (17) and Lindblad et al. (68), who used data from Robertson et al. (95). The rate constants for Ca²⁺ binding to calmodulin and troponin were scaled with the use of a Q₁₀ of 1.8, according to Lindblad et al. (68). The concentration of calsequestrin within the SR release compartment was set to a relatively low value of 10 mM, the value used by Luo and Rudy (71), Courtemanche et al. (15), and Priebe and Beuckelmann (91). The on and off rates for Ca²⁺ binding to calsequestrin were based on the study of Cannell and Allen (12), adjusted to 37°C via a Q₁₀ of 1.6 (see Ref. 68).

Ion concentration homeostasis. Our model also includes material balance expressions to define the temporal variations in [Na⁺]_i, [K⁺]_i, and [Ca²⁺]_i. [Na⁺]_i and [K⁺]_i were assumed to be homogeneous, because subspace [Na⁺] (or [K⁺]) was nearly equal to myoplasmic [Na⁺] (or [K⁺]) unless the diffusion time constant for Na⁺ (or K⁺) was more than 10 times larger than that for Ca²⁺. As suggested by Nygren et al. (88), there may be electroneutral Na⁺ influx (and K⁺ efflux as well) via electroneutral transport mechanisms such as Na⁺/H⁺ exchange and Na⁺-K⁺-2Cl