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 cotransport; thus it is difficult to estimate [Na⁺]_i (and [K⁺]_i) accurately from model simulation. In our simulations, net electroneutral Na⁺ and K⁺ fluxes were assumed to be zero; in the free running model with the standard parameter values listed in the APPENDIX, the intracellular free Na⁺ and K⁺ levels averaged 9.44 and 140.0 mM, respectively, during a long period of normal pacemaking. [Na⁺]_o, [K⁺]_o, and [Ca²⁺]_o were set equal to 140, 5.4, and 2 mM, respectively.

Numerical Integration (Dynamic Simulation)

6

The APA as a voltage difference between the local potential minimum and maximum, as well as the CL, were determined for each calculation of a cycle. Numerical integration was continued until the differences in both APA and CL between the newly calculated cycle and the preceding one became <1 × 10³ of the preceding APA and CL values. When periodic behavior was irregular or unstable, model dynamics were computed for 30~60 s; all potential extrema and CL values were then plotted in the diagram. The value of current conductance is expressed as a ratio to the control value unless otherwise stated.

Numerical computations were performed on Power Macintosh G4 computers (Apple Computers; Cupertino, CA) using Matlab 5.2 (MathWorks; Natick, MA).

	RESULTS AND DISCUSSION

TOP ABSTRACT INTRODUCTION THEORY AND METHODS RESULTS AND DISCUSSION APPENDIX REFERENCES

Dynamic Properties of Simulated Pacemaker Activity

Spontaneous APs and sarcolemmal ionic currents. Figure 5, left, shows spontaneous APs and temporal behavior of sarcolemmal ionic currents simulated by the present model with the standard parameter values listed in the APPENDIX. The MDP, POP, and APA are 58.6, +16.6, and 75.2 mV, respectively; the CL and AP duration at 50% repolarization (APD₅₀) are 307.5 and 107.0 ms, respectively. The simulated AP parameters of our model cell are listed in Table 1, along with those of previous SA node models as well as corresponding experimental data for comparison (see also Fig. 6). The AP parameters of our model appear to be in reasonable agreement with the mean experimental values recently determined for central SA node (primary pacemaker) cells and spindle- or spider-shaped cells. The previous models developed by Wilders et al. (124), Demir et al. (17), and Dokos et al. (24) are all likely to reflect the activity of a typical transitional cell. Compared with these previous models for transitional cells, our model for primary pacemaker cells exhibited 1) more positive MDP, 2) relatively small APA, and 3) long APD₅₀. These AP characteristics of our model are comparable to those of primary pacemaker cells in the central region of the rabbit SA node (8, 44, 55, 57), reflecting the regional difference in AP parameters.

View larger version (35K): [in this window] [in a new window]   Fig. 5.   Model-generated behavior of spontaneous APs, underlying transmembrane ionic currents (left) and intracellular Ca2+ dynamics (right) in the steady state. The differential equations were numerically solved for 11 s with an initial condition appropriate to a zero current potential (24.3 mV) and a 1-ms stimulus of 1 pA/pF for triggering an AP. Model cell behavior during the last 1 s (from MDP) of the simulated free-running activity is depicted. Note the differences in the ordinate scales for individual currents. The intracellular Ca2+ dynamics shown on the right include the changes in Ca2+ concentrations (in the subspace, myoplasm, and SR), associated percent changes in the occupancies of intracellular Ca2+ buffers, and Ca2+ fluxes produced by the SR during Ca2+ uptake, release, and transfer. Changes in [Na+]i and [K+]i during pacemaker activity are also shown at the bottom. See the Glossary for definitions of the abbreviations.

View larger version (26K): [in this window] [in a new window]   Fig. 6.   Spontaneous APs and underlying transmembrane time-dependent currents in various SA node models developed by Wilders et al. (124) (A), Demir et al. (17) (B), Dokos et al. (24) (C), Zhang et al. (130) (D; central model), and 5) the present study (E). The differential equations for previous models were numerically solved for 11 s with the initial conditions provided in the original articles. Model cell behaviors during the last 600 ms (from MDP) are depicted. Spontaneous APs recorded in the central SA node (leading pacemaker site) and the currents experimentally determined by the AP clamp method are also shown (F) for validation of the models: AP, Miyamae and Goto (79); ICa,L, Zaza et al. (129); IKr, Ono and Ito (90); ICa,T, Doerr et al. (23); and Ih, Zaza et al. (129). See the Glossary for definitions of the abbreviations.

View larger version (25K): [in this window] [in a new window]   Fig. 7.   Time-independent background and transporter currents during spontaneous pacemaking in various SA node models formulated by Wilders et al. (124) (A), Demir et al. (17) (B), Dokos et al. (24) (C), Zhang et al. (130) (D; central model), and the present study (E). See the Glossary for definitions of the abbreviations.

Intracellular Ca²+ dynamics. Intracellular Ca²⁺ dynamics (concentration changes, SR Ca²⁺ uptake and release, and Ca²⁺ buffering) during normal pacemaking are shown in Fig. 5, right. The free Ca²⁺ concentrations in subspace ([Ca²⁺]_sub) and myoplasm ([Ca²⁺]_i) during the spontaneous APs were 0.18~1.81 and 0.25~0.68 µM, respectively. Although there are no available data on intracellular Ca²⁺ concentrations in the rabbit SA node cell, these values are comparable to those measured in other pacemaker cells (47, 53).

View larger version (38K): [in this window] [in a new window]   Fig. 8.   Intracellular Ca2+ dynamics during spontaneous pacemaking in various SA node models formulated by Wilders et al. (124) (A), Demir et al. (17) (B), Dokos et al. (24) (C), and the present study (D). Note the differences in the ordinate scales for individual models. See the Glossary for definitions of the abbreviations.

Contribution of I_st to pacemaker activity. Figure 9 shows the simulated effects of block of I_st on pacemaker activity of our model cell. When the maximum I_st conductance (g_st) was reduced from the control value to zero, pacemaking slowed, with CL increasing from 307.5 to 473.2 ms, but did not cease. The effects of eliminating I_st on MDP and POP were relatively small, with MDP hyperpolarized from 58.6 to 63.6 mV. These alterations in the AP waveform by the removal of I_st could partly be compensated for by reducing the maximum I_Kr conductance (g_Kr): the I_st-removed model cell with g_Kr reduced to 75% of the control value produced spontaneous APs with MDP = 58.9 mV and CL = 379.8 ms (see Fig. 9, bottom right). Thus I_st is not essential for generating pacemaker activity.

View larger version (27K): [in this window] [in a new window]   Fig. 9.   Contribution of Ist to pacemaker activity of the SA node. Left: dynamics of the model cell during inhibition of Ist. Changes in MDP/POP (top) and CL (bottom) are shown as functions of the maximum Ist conductance (gst). The gst value, shown as a ratio to the control value of 15 pS/pF, was reduced at an interval of 0.01. Right: simulated spontaneous APs of the standard (control) system (top), Ist-removed system (middle), and Ist-removed system with the IKr conductance (gKr) reduced to 75% of the control (bottom). See the Glossary for definitions of the abbreviations.

Effects of Modulating Sarcolemmal Ionic Currents on Pacemaker Activity

Simulated blockade of I_Ca,L. We first examined the effects of blocking I_Ca,L on pacemaker activity of the model cell by decreasing the maximum I_Ca,L conductance (g_Ca,L). Figure 10A shows the dynamic behaviors of the I_st-incorporated [I_st(+)] and I_st-removed [I_st()] systems during g_Ca,L decrease, depicting MDP/POP and CL as functions of g_Ca,L (i.e., bifurcation diagrams for the bifurcation parameter g_Ca,L). As g_Ca,L diminished, APA was reduced and CL increased; POP gradually decreased with reducing g_Ca,L, whereas MDP was relatively stable. Blocking I_Ca,L by 76.6% (reducing g_Ca,L to 23.4% of the control) abolished spontaneous activity of the standard system, with V settling at 29.5 mV; the complete block of I_Ca,L yielded a resting potential of 32.6 mV. In the I_st-removed system, pacemaker activity abruptly ceased at g_Ca,L = 0.298 (V = 37.6 mV) via irregular (chaotic) dynamics between g_Ca,L = 0.439 and 0.304; the complete block of I_Ca,L yielded a resting potential of 42.0 mV. Note that the bifurcation structure as a way to quiescence during inhibition of I_Ca,L depends on whether the model cell includes I_st or not: APA of our standard model gradually and continuously decreased during the g_Ca,L decline, whereas the irregular (chaotic) dynamics occurred when I_st was removed or concomitantly blocked.

View larger version (28K): [in this window] [in a new window]   Fig. 10.   Effects of block of ICa,L on pacemaker activity of the SA node. A: dynamics during ICa,L inhibition of our standard model system including Ist [Ist(+)] and the Ist-removed system with gKr reduced to 75% of the control [Ist()]. Changes in MDP/POP (top) and CL (bottom) are shown as functions of the maximum ICa,L conductance (gCa,L). The gCa,L value, shown as a ratio to the control value of 0.58 nS/pF, was reduced at an interval of 0.001. Note the difference in the scales of CL. B: experimentally observed behaviors of SA node cells in the presence of a Ca2+ antagonist (ICa,L blocker) and simulated dynamics of our model cell with blocked ICa,L. The experimental data are from Kohlhardt et al. (58) (top) and Vinogradova et al. (122) (bottom). Simulated spontaneous APs of the standard (top) or Ist-removed (bottom) system with reduced gCa,L are shown at the bottom of each experimental record. See the Glossary for definitions of the abbreviations.

View larger version (17K): [in this window] [in a new window]   Fig. 11.   Effects of block of ICa,L on pacemaker activities of various SA node models developed by Wilders et al. (124) (A), Demir et al. (17) (B), Dokos et al. (24) (C), Zhang et al. (130) (D; central model), and the present study (E; standard model). Changes in MDP/POP (top) and CL (bottom) during gCa,L decrease are depicted. The gCa,L value, shown as a ratio to the control value for each model, was reduced at an interval of 0.001. Note the differences in the scales of CL. See the Glossary for definitions of the abbreviations.

Simulated blockade of I_Kr. Figure 12 illustrates the simulated effects of blocking I_Kr [decreasing the maximum I_Kr conductance (g_Kr)] on pacemaker activity of the standard and I_st-removed model systems along with experimental data showing the behaviors of central SA node cells during applications of a selective I_Kr blocker, E-4031 (56). In the simulations, reducing g_Kr markedly depolarized MDP, whereas POP was little changed until g_Kr diminished to ~0.5; MDP depolarization during the g_Kr decline was more prominent than that during the g_Ca,L decline (compare Figs. 10 and 12). With reducing g_Kr, CL in the standard system first increased and then decreased, being fairly stable; in contrast, CL in the I_st-removed system monotonously decreased. Whether I_st is incorporated or not, the bifurcation structure of the model system during g_Kr reduction was essentially the same: APA gradually and continuously declined to zero. Pacemaker activities of the standard and I_st-removed systems ceased when g_Kr decreased by 64.0 and 64.9%, respectively, with V settling at 13.8 and 13.5 mV, respectively. The complete block of I_Kr yielded a resting potential of 7.6 mV in the standard system and 8.3 mV in the I_st-removed system.

View larger version (20K): [in this window] [in a new window]   Fig. 12.   Effects of block of IKr on pacemaker activity of the SA node. Left: simulated dynamics during IKr inhibition of our standard model [Ist(+)] and Ist-removed model [Ist()]. Changes in MDP/POP (top) and CL (bottom) are shown as functions of gKr. The gKr value, shown as a ratio to the control value, was reduced at an interval of 0.001. Note that the control gKr value for the Ist-removed system is 0.75 times that for the standard system. Right: experimentally observed behaviors of SA node cells during applications of a selective IKr blocker, E-4031, from Kodama et al. (56). Note the gradual and continuous decline in APA with accentuated depolarization of MDP on the way to cessation of pacemaker activity. See the Glossary for definitions of the abbreviations.

View larger version (16K): [in this window] [in a new window]   Fig. 13.   Effects of block of IKr on pacemaker activities of various SA node models formulated by Wilders et al. (124) (A), Demir et al. (17) (B), Dokos et al. (24) (C), Zhang et al. (130) (D; central model), and the present study (E; standard model). Changes in MDP/POP (top) and CL (bottom) during gKr decrease are depicted. The gKr value, shown as a ratio to the control value for each model, was reduced at an interval of 0.001. See the Glossary for definitions of the abbreviations.

Block of I_Ca,T. The effects of block of I_Ca,T on pacemaker frequency have been largely different in previous experimental reports: Doerr et al. (23) found that the block of I_Ca,T by 40 µM Ni²⁺ exerted a pronounced negative chronotropic effect on pacemaker activity by nearly doubling the basal CL (see also Refs. 84 and 103). In contrast, Wilders et al. (124) reported that the Ni²⁺-induced block of I_Ca,T had a negligible influence on pacemaker frequency. Because of this inconsistency in experimental results, the simulated effects of eliminating I_Ca,T on CL of the previous models were also quite different: the increase in CL was relatively small in the models of Wilders et al. (124) and Dokos et al. (24), whereas it was relatively large in the others (see Table 3). Incorporating the I_Ca,T expressions of Demir et al. (17) based on the AP-clamp experiment of Doerr et al. (23), our model predicted a 17.0~22.4% increase in CL by the complete block of I_Ca,T (Fig. 14, top); this relatively large effect is comparable to the experimental data from Hagiwara et al. (37) and Satoh (103). Owing to the large variation in experimental results, however, it is difficult to validate a model from the contribution of I_Ca,T to diastolic depolarization.

View larger version (23K): [in this window] [in a new window]   Fig. 14.   Effects of block of ICa,T (top), Ih (middle), or 4-AP-sensitive currents (bottom) on spontaneous APs of our model cell. Both the standard system including Ist (left) and the Ist-removed system with gKr reduced to 75% of the control (middle) were tested. Control APs are shown as thin lines; dashed lines indicate zero potentials. Experimentally observed responses of rabbit SA node cells to application of 40 µM Ni2+ (to block ICa,T), 2 mM Cs+ (to block Ih), or 5 mM 4-AP are also shown for comparison (right). The data are from Hagiwara et al. (37) (top), Nikmaram et al. (83) (middle), and Boyett et al. (7) (bottom). See the Glossary for definitions of the abbreviations.

Block of I_h. We also simulated the effect of Cs⁺-induced I_h block on pacemaker activity by setting I_h = 0. As listed in Table 3, the contributions of I_h to pacemaker depolarization of previous models were different: the increase in CL on the removal of I_h was relatively large in the models of Wilders et al. (124) and Dokos et al. (24), whereas it was relatively small in the others. Incorporating the I_h expressions of Wilders et al. (124), our model predicted a 8.4~25.3% increase in CL by the complete block of I_h (Fig. 14, middle). The effect of block of I_h on CL was much greater in the I_st-removed system than in the standard system.

Block of 4-AP-sensitive currents. The effects of block of 4-AP-sensitive currents (I_to and I_sus) on AP configurations are shown in Fig. 14, bottom. Complete block of 4-AP-sensitive currents in the standard and I_st-removed systems prolonged APD₅₀ by 26.7 and 22.3%, respectively, and also caused positive shifts in POP (from +16.6 to +32.9 mV in the standard system). On elimination of both I_to and I_sus, the CL of the standard system increased by 5.1%, whereas that of the I_st-removed system decreased by 7.2%. The simulated increase in APD₅₀ (22.3~26.7%) is in good agreement with the experimental data of Boyett et al. (7): block of 4-AP-sensitive currents by 5 mM 4-AP caused an AP duration prolongation of 25 ± 5% in small balls from central SA node tissues. They also reported a 5 ± 2% decrease in CL by 5 mM 4-AP, consistent with the behavior of our I_st-removed system.

Effects of Modulating Intracellular Ca²+ Dynamics on Pacemaker Activity

Effects of block of SR Ca²+ release. SR Ca²⁺ release is known to affect intracellular Ca²⁺ transients, Ca²⁺-dependent inactivation of I_Ca,L, and activation of I_NaCa, possibly playing an important role in regulating pacemaker activity (39, 67, 104). We therefore assessed the contributions of SR Ca²⁺ release to pacemaker activities of the model cells by reducing the SR volume to 0.1% of the control. As shown in Table 4, the influences of block of SR Ca²⁺ release on pacemaker frequencies of the model cells were different but relatively small. In our standard model for primary pacemaker cells, elimination of SR Ca²⁺ release decreased the peak [Ca²⁺]_sub only by 22.4% (from 1.81 to 1.41 µM) and increased CL by only 3.4% (from 307.5 to 317.6 ms). Thus, in our primary pacemaker model, SR Ca²⁺ release plays only a minor role in generating subsarcolemmal Ca²⁺ transients and regulating pacemaker activity (via modifications of the Ca²⁺-dependent inactivation of I_Ca,L and activation of I_NaCa). This simulated result agrees well with the reports of Janvier and Boyett (52) and Miyamae and Goto (79), in which application of ryanodine to abolish SR Ca²⁺ release resulted in little or no slowing of pacemaker activity in SA node cells. The Ca²⁺-dependent inactivation of L-type Ca²⁺ channels in rat ventricular or atrial myocytes is mediated primarily by Ca²⁺ release from the SR to the microdomain (5, 60). In our SA node model, however, the Ca²⁺-dependent inactivation of L-type Ca²⁺ channels is chiefly mediated by Ca²⁺ influx through the L-type Ca²⁺ channel itself. Consistent with our model behavior, ryanodine little slowed the fast decay of I_Ca,L in SA node cells (67), suggesting that SR Ca²⁺ release does not contribute to the Ca²⁺-dependent inactivation of L-type Ca²⁺ channels in the SA node. Morphological studies indicate that the SR is poorly formed in typical SA node cells and is therefore unlikely to play a major role in modulation of pacemaker activity (see Ref. 17). As suggested by Janvier and Boyett (52), the contributions of SR Ca²⁺ release to intracellular Ca²⁺ transients and to pacemaker activity (via inactivation of I_Ca,L or activation of I_NaCa) would be relatively small in SA node cells because of the poor development of the SR.

Effects of Ca²+ buffers. To buffer intracellular Ca²⁺, Ca²⁺ buffers such as EGTA and BAPTA are routinely used for whole cell patch-clamp experiments. These Ca²⁺ buffers to reduce intracellular Ca²⁺ transients have been shown to affect the pacemaker activity of SA node cells (67, 122). We therefore examined the effects of the Ca²⁺ buffers BAPTA and EGTA on the pacemaker activity of SA node models. Wilders et al. (124) simulated EGTA buffering by fixing [Ca²⁺]_i to 80 nM. However, it is known that the slow buffer EGTA cannot efficiently buffer free Ca²⁺ in the subsarcolemmal space: although both BAPTA and EGTA can efficiently suppress global (myoplasmic) Ca²⁺ transients, only the fast buffer BAPTA is efficient in buffering local (subsarcolemmal) Ca²⁺ transients (see Refs. 5, 106, 122, and 128). Therefore, we calculated intracellular Ca²⁺ dynamics in the presence of BAPTA or EGTA to simulate the differential effects of Ca²⁺ buffers on [Ca²⁺]_sub and on spontaneous APs. As shown in Fig. 15 and Table 5, the fast buffer BAPTA at 10 mM remarkably reduced both the subsarcolemmal and myoplasmic Ca²⁺ transients and dramatically slowed pacemaker activity with CL increasing by 43.4% in our standard model. This negative chronotropic effect is a consequence of the reduction in I_NaCa during the late phase of pacemaker depolarization. In contrast to BAPTA, the slow buffer EGTA at 10 mM could not sufficiently inhibit the subsarcolemmal Ca²⁺ transient, whereas it almost completely suppressed the myoplasmic Ca²⁺ transient to ~0.1 µM; pacemaker frequency was little affected by 10 mM EGTA. These differential effects of BAPTA and EGTA on pacemaker frequency have been experimentally observed (122): only BAPTA significantly reduced the rate of spontaneous APs (by 54%), whereas EGTA did not (see Fig. 15). In contrast to our model, previous models have failed to simulate the differential responses of SA node cells to Ca²⁺ buffers: both BAPTA and EGTA dramatically suppressed intracellular Ca²⁺ transients, exerting negative (or positive) chronotropic effects (see Table 5). Thus only our model could reproduce the differential responses of SA node cells to BAPTA and EGTA. Our model incorporates two compartments for intracellular Ca²⁺, including the subsarcolenmal space as a diffusion barrier for Ca²⁺, whereas previous models assumed only one intracellular compartment. Subspace Ca²⁺ dynamics are indispensable for simulating the differential effects of BAPTA and EGTA on SA node dynamics. A recent experimental report suggests local high Ca²⁺ gradients in the subsarcolemmal microdomain in SA node cells (122).

View larger version (28K): [in this window] [in a new window]   Fig. 15.   A: effects of Ca2+ buffers, BAPTA and EGTA (10 mM), on spontaneous APs (top) and intracellular Ca2+ concentrations (middle and bottom) in our standard model. Simulating the whole cell patch-clamp recording for internally perfused SA node cells, [Na+]i and [K+]i were fixed at 10 and 140 mM, respectively. The equations were numerically integrated for 11 s with an initial condition of a steady state at 60 mV, and membrane potential behavior and intracellular Ca2+ dynamics during the last 1 s (from MDP) were then depicted. Control APs and Ca2+ transients are shown as thin lines. B: experimentally observed responses of rabbit SA node cells to BAPTA and EGTA, the data from Vinogradova et al. (122), are also shown for comparison. See the Glossary for definitions of the abbreviations.

Achievements and Limitations of the Present Model

Despite several improvements over previous models, the present model still has some limitations. Quantitatively, our model predictions (e.g., simulated changes in pacemaker frequency or resting potentials during inhibition of I_Kr) exhibited some inconsistencies with experimental data. These discrepancies may be due to 1) incomplete experimental data on the kinetics (or density) of ionic currents as well as on intracellular Ca²⁺ dynamics (and SR functions) in the SA node; 2) the large heterogeneity of SA node cells (e.g., regional differences in current densities); 3) poor selectivity of the agents used experimentally to block ionic currents or state-dependent kinetics of channel blockade; 4) the existence of some ionic current components not included in our model (e.g., Ca²⁺-activated Cl current; see Ref. 92); 5) the spatial heterogeneity of intracellular distributions of Ca²⁺ and Na⁺ (see Ref. 13); 6) the lack of intracellular or intramembrane regulatory systems in our model; and 7) inappropriate kinetic formulations based on the Hodgkin-Huxley formalism for ionic channel currents such as I_Ca,L and I_Kr, which are essential to pacemaker generation (see Refs. 14, 43, 50, 88, 115, and 127). These points are of great importance in future modeling to develop more sophisticated models.

It should be noted that one limitation in estimating the contribution of a current to pacemaking, as well as in measuring a current by AP clamp, is the need of a pharmacological block of the current (see Ref. 72). Bindings of channel blockers used for experiments are known to be nonspecific and voltage (state) dependent (5, 70, 119, 129). Therefore, the currents recorded during AP clamp may be different from the real pure current, rather reflecting the total current blocked by an agent. As Verheijck et al. (119) suggested, changes in AP parameters induced by a drug are most likely the result of a combination of direct and indirect effects on various ionic currents. An agent used to block a current may directly affect other ionic currents, or, alternatively, secondarily affect ion channel properties via changes in intracellular Ca²⁺ concentrations and/or modifications of second messengers (see Refs. 9 and 129). Thus we should be cautious when simulated current waveforms or blocking effects are compared with experimental data.

There are numerous factors involved in the regulation of ion channel and SR functions in intact cardiac myocytes: intracellular second messengers (e.g., Ca²⁺, cyclic nucleotides, and protein kinases), as well as intramembrane modulators (e.g., receptors, G proteins, and enzymes), have been shown to modify various ionic current systems (for reviews, see Refs. 5, 9, 16, 25, and 91). Therefore, incorporating the dynamics of these modulating factors would be indispensable for predicting SA node behavior more accurately and might at least in part solve the discrepancies between model predictions and experimental data. Tentative models including putative intracellular Ca²⁺-dependent changes in ionic currents or modifications of channel functions by receptors and second messengers have recently been published (9, 16, 25, 75). In this study, however, we did not incorporate these modulating factors, because the aim of this study was not to develop a complete model but to develop an improved model for investigating the dynamical mechanisms of pacemaker generation. The present model could also serve as a base model for the development of a complete model incorporating the intramembrane modulators and intracellular second messengers and further the genetic regulation of ion channels and transporters.

Despite many limitations, our model and simulations can provide the guidance for future developments of more sophisticated single SA node cell models, multicellular (two or three dimensional) models of the intact SA node, and whole heart models.