## Abstract

We developed an improved mathematical model for a single primary pacemaker cell of the rabbit sinoatrial node. Original features of our model include *1*) incorporation of the sustained inward current (*I*
_{st}) recently identified in primary pacemaker cells, *2*) reformulation of voltage- and Ca^{2+}-dependent inactivation of the L-type Ca^{2+}channel current (*I*
_{Ca,L}), *3*) new expressions for activation kinetics of the rapidly activating delayed rectifier K^{+} channel current (*I*
_{Kr}), and *4*) incorporation of the subsarcolemmal space as a diffusion barrier for Ca^{2+}. We compared the simulated dynamics of our model with those of previous models, as well as with experimental data, and examined whether the models could accurately simulate the effects of modulating sarcolemmal ionic currents or intracellular Ca^{2+} dynamics on pacemaker activity. Our model represents significant improvements over the previous models, because it can *1*) simulate whole cell voltage-clamp data for*I*
_{Ca,L}, *I*
_{Kr}, and*I*
_{st}; *2*) reproduce the waveshapes of spontaneous action potentials and ionic currents during action potential clamp recordings; and *3*) mimic the effects of channel blockers or Ca^{2+} buffers on pacemaker activity more accurately than the previous models.

- rabbit sinoatrial node
- nonlinear dynamical system
- computer simulation
- bifurcation diagram

pacemaker activityof 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^{2+} buffering (by calmodulin, troponin, and calsequestrin) and new formulations for Ca^{2+}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^{2+} 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^{2+} channel current (*I*
_{Ca,L}) were not based on experimental data. Second, the Ca^{2+}-dependent inactivation of*I*
_{Ca,L}, known as an essential property of the L-type Ca^{2+} channel in the rabbit SA node (10,67), was not formulated [or not directly dependent on the intracellular Ca^{2+} concentration ([Ca^{2+}]_{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^{2+}]_{i}-dependent inactivation process, they assumed an extraordinarily large maximum *I*
_{Ca,L}conductance (*g*
_{Ca,L}) and a very high affinity for Ca^{2+} binding to the inactivation site of L-type Ca^{2+} channels; such a large *g*
_{Ca,L}with a high affinity for Ca^{2+} 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^{2+} 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. 68and 71), probably too high for SA node cells. Finally, SR volumes or Ca^{2+} 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^{2+} during an AP cycle is not zero), whereas intracellular Ca^{2+} 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^{2+}-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^{2+} 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^{2+} 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^{2+} channel current (*I*
_{Ca,T}), hyperpolarization-activated current (*I*
_{h}), or 4-AP-sensitive currents; *5*) effects of blocking SR Ca^{2+} release (by ryanodine) on pacemaker frequency; and *6*) negative chronotropic effects of Ca^{2+} 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
- APD
_{50} - 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

- Ib,Ca
- Background Ca
^{2+}current *I*_{b,K}- Background K
^{+}current *I*_{b,Na}- Background Na
^{+}current *I*_{Ca,L}- L-type Ca
^{2+}channel current *I*_{Ca,T}- T-type Ca
^{2+}channel current *I*_{h}- Hyperpolarization-activated current
*I*_{K}- Delayed rectifier K
^{+}current *I*_{K,ACh}- Muscarinic K
^{+}channel current *I*_{Kr}- Rapid component of the delayed rectifier K
^{+}current *I*_{Ks}- Slow component of the delayed rectifier K
^{+}current *I*_{Na}- Fast Na
^{+}channel current *I*_{NaCa}- Na
^{+}/Ca^{2+}exchanger current *I*_{Naen}- Electroneutral Na
^{+}influx current *I*_{NaK}- Na
^{+}-K^{+}pump current *I*_{NaKmax}- Maximum
*I*_{NaK} *I*_{p,Ca}- Sarcolemmal Ca
^{2+}pump current *I*_{st}- Sustained inward current
*I*_{sus}- Sustained component of the 4-AP-sensitive current
*I*_{to}- Transient component of the 4-AP-sensitive current

### Cell Geometry

- Cm
- Cell membrane capacitance
- V
_{cell} - Cell volume
- V
_{i} - Myoplasmic volume available for Ca
^{2+}diffusion - V
_{rel} - Volume of junctional SR (Ca
^{2+}release store) - V
_{sub} - Subspace volume
- V
_{up} - Volume of network SR (Ca
^{2+}uptake store)

### Ionic Concentrations

- Ca2+i
- Myoplasmic Ca
^{2+}concentration - [Ca
^{2+}]_{o} - Extracellular Ca
^{2+}concentration - [Ca
^{2+}]_{rel} - Ca
^{2+}concentration in the junctional SR - [Ca
^{2+}]_{sub} - Subspace Ca
^{2+}concentration - [Ca
^{2+}]_{up} - Ca
^{2+}concentration in the network SR - [K
^{+}]_{i} - Intracellular K
^{+}concentration - [K
^{+}]_{o} - Extracellular K
^{+}concentration - [Mg
^{2+}]_{i} - Intracellular Mg
^{2+}concentration - [Na
^{+}]_{i} - Intracellular Na
^{+}concentration - [Na
^{+}]_{o} - Extracellular Na
^{+}concentration

### Equilibrium (Reversal) Potentials

- ECa,L
- Apparent reversal potential of
*I*_{Ca,L} *E*_{Ca,T}- Apparent reversal potential of
*I*_{Ca,T} *E*_{K}- Equilibrium (Nernst) potential for K
^{+} *E*_{Kr}- Reversal potential of
*I*_{Kr} *E*_{Ks}- Reversal potential of
*I*_{Ks} *E*_{mh}- Reversal potential of
*I*_{Na} *E*_{Na}- Equilibrium (Nernst) potential for Na
^{+} *E*_{st}- Apparent reversal potential of
*I*_{st}

### Sarcolemmal Ionic Currents

- L
- Activation gating variable for
*I*_{Ca,L} *d*_{L,∞}- Steady-state
*d*_{L} *d*_{T}- Activation gating variable for
*I*_{Ca,T} *d*_{T,∞}- Steady-state
*d*_{T} *f*_{Ca}- Ca
^{2+}-dependent inactivation gating variable for*I*_{Ca,L} *f*_{Ca,∞}- Steady-state
*f*_{Ca} *f*_{L}- Voltage-dependent inactivation gating variable for
*I*_{Ca,L} *f*_{L,∞}- Steady-state
*f*_{L} *f*_{T}- Inactivation gating variable for
*I*_{Ca,T} *f*_{T,∞}- Steady-state
*f*_{T} *F*_{Na}- Fraction of slow inactivation of
*I*_{Na} *g*_{b,Na}- Background Na
^{+}conductance *g*_{Ca,L}- Maximum
*I*_{Ca,L}conductance *g*_{Ca,T}- Maximum
*I*_{Ca,T}conductance *g*_{h}- Maximum
*I*_{h}conductance *g*_{h,K}- K
^{+}current component of*g*_{h} *g*_{h,Na}- Na
^{+}current component of*g*_{h} *g*_{K,ACh}- Scaling factor for
*I*_{K,ACh} *g*_{Kr}- Maximum
*I*_{Kr}conductance *g*_{Ks}- Maximum
*I*_{Ks}conductance *g*_{Na}- Maximum
*I*_{Na}conductance *g*_{st}- Maximum
*I*_{st}conductance *g*_{sus}- Maximum
*I*_{sus}conductance *g*_{to}- Maximum
*I*_{to}conductance *h*- Inactivation gating variable for
*I*_{Na} *h*_{∞}- Steady-state
*h* *h*_{F}- Fast inactivation gating variable for
*I*_{Na} *h*_{S}- Slow inactivation gating variable for
*I*_{Na} *j*_{Ca,sm}- Net Ca
^{2+}flux through the sarcolemmal membrane *k*_{NaCa}- Scaling factor for
*I*_{NaCa} *K*_{1ni}- Dissociation constant for intracellular Na
^{+}binding to first site on*I*_{NaCa}transporter *K*_{2ni}- Dissociation constant for intracellular Na
^{+}binding to second site on*I*_{NaCa}transporter *K*_{3ni}- Dissociation constant for intracellular Na
^{+}binding to third site on*I*_{NaCa}transporter *K*_{1no}- Dissociation constant for extracellular Na
^{+}binding to first site on*I*_{NaCa}transporter *K*_{2no}- Dissociation constant for extracellular Na
^{+}binding to second site on*I*_{NaCa}transporter *K*_{3no}- Dissociation constant for extracellular Na
^{+}binding to third site on*I*_{NaCa}transporter *K*_{ci}- Dissociation constant for intracellular Ca
^{2+}binding to*I*_{NaCa}transporter *K*_{co}- Dissociation constant for extracellular Ca
^{2+}binding to*I*_{NaCa}transporter *K*_{cni}- Dissociation constant for intracellular Na
^{+}and Ca^{2+}simultaneous binding to*I*_{NaCa}transporter *K*_{mfCa}- Dissociation constant for Ca
^{2+}-dependent inactivation of*I*_{Ca,L} *K*_{mKp}- Half-maximal [K
^{+}]_{o}for*I*_{NaK} *K*_{mNap}- Half-maximal [Na
^{+}]_{i}for*I*_{NaK} *m*- Activation gating variable for
*I*_{Na} *m*_{∞}- Steady-state
*m* *n*- Activation gating variable for
*I*_{Ks} *n*_{∞}- Steady-state
*n* *p*_{a}- Activation gating variable for
*I*_{Kr} *p*_{a,∞}- Steady-state
*p*_{a} *p*_{aF}- Fast activation gating variable for
*I*_{Kr} *p*_{aS}- Slow activation gating variable for
*I*_{Kr} *p*_{i}- Inactivation gating variable for
*I*_{Kr} *p*_{i,∞}- Steady-state
*p*_{i} *q*- Inactivation gating variable for
*I*_{to} *q*_{∞}- Steady-state
*q* *q*_{a}- Activation gating variable for
*I*_{st} *q*_{a,∞}- Steady-state
*q*_{a} *q*_{i}- Inactivation gating variable for
*I*_{st} *q*_{i,∞}- Steady-state
*q*_{i} *Q*_{ci}- Fractional charge movement during intracellular Ca
^{2+}occlusion reaction of*I*_{NaCa}transporter *Q*_{co}- Fractional charge movement during extracellular Ca
^{2+}occlusion reaction of*I*_{NaCa}transporter *Q*_{n}- Fractional charge movement during Na
^{+}occlusion reactions of*I*_{NaCa}transporter *r*- Activation gating variable for
*I*_{to}and*I*_{sus} *r*_{∞}- Steady-state
*r* *x*- Gating variable
*x*_{∞}- Steady-state value of
*x* *y*- Activation gating variable for
*I*_{h} *y*_{∞}- Steady-state value of
*y* - α
_{dL} - Opening rate constant of
*d*_{L} - α
_{fCa} - Ca
^{2+}dissociation rate constant for*I*_{Ca,L} - α
_{qa} - Opening rate constant of
*q*_{a} - α
_{qi} - Opening rate constant of
*q*_{i} - α
_{n} - Opening rate constant of
*n* - β
_{dL} - Closing rate constant of
*d*_{L} - β
_{fCa} - Ca
^{2+}association rate constant for*I*_{Ca,L} - β
_{n} - Closing rate constant of
*n* - β
_{qa} - Closing rate constant of
*q*_{a} - β
_{qi} - Closing rate constant of
*q*_{i} - τ
_{dL} - Time constant of
*d*_{L} - τ
_{dT} - Time constant of
*d*_{T} - τ
_{fCa} - Time constant of
*f*_{Ca} - τ
_{fL} - Time constant of
*f*_{L} - τ
_{fT} - Time constant of
*f*_{T} - τ
_{hF} - Time constant of
*h*_{F} - τ
_{hS} - Time constant of
*h*_{S} - τ
_{m} - Time constant of
*m* - τ
_{n} - Time constant of
*n* - τ
_{paF} - Time constant of
*p*_{aF} - τ
_{paS} - Time constant of
*p*_{aS} - τ
_{q} - Time constant of
*q* - τ
_{qa} - Time constant of
*q*_{a} - τ
_{qi} - Time constant of
*q*_{i} - τ
_{r} - Time constant of
*r* - τ
_{x} - Time constant for a gating variable
*x* - τ
_{y} - Time constant of
*y*

### Ca^{2+}* Diffusion*

- jCa,dif
- Ca
^{2+}diffusion flux from subspace to myoplasm - τ
_{dif,Ca} - Time constant of Ca
^{2+}diffusion from the subspace to myoplasm

### SR Function

- jrel
- Ca
^{2+}release flux from the junctional SR to subspace *j*_{tr}- Ca
^{2+}transfer flux from the network to junctional SR *j*_{up}- Ca
^{2+}uptake flux from the myoplasm to network SR *K*_{rel}- Half-maximal [Ca
^{2+}]_{sub}for Ca^{2+}release from the junctional SR *K*_{up}- Half-maximal [Ca
^{2+}]_{i}for Ca^{2+}uptake by the Ca^{2+}pump in the network SR *P*_{rel}- Rate constant for Ca
^{2+}release from the junctional SR *P*_{up}- Rate constant for Ca
^{2+}uptake by the Ca^{2+}pump in the network SR - τ
_{tr} - Time constant for Ca
^{2+}transfer from the network to junctional SR

### Ca^{2+}* Buffering*

- CQtot
- Total calsequestrin concentration
- [CM]
_{tot} - Total calmodulin concentration
*f*_{CMi}- Fractional occupancy of calmodulin by Ca
^{2+}in myoplasm *f*_{CMs}- Fractional occupancy of calmodulin by Ca
^{2+}in subspace *f*_{CQ}- Fractional occupancy of calsequestrin by Ca
^{2+} *f*_{TC}- Fractional occupancy of the troponin-Ca site by Ca
^{2+} *f*_{TMC}- Fractional occupancy of the troponin-Mg site by Ca
^{2+} *f*_{TMM}- Fractional occupancy of the troponin-Mg site by Mg
^{2+} *k*_{bBAPTA}- Ca
^{2+}dissociation constant for BAPTA *k*_{bCM}- Ca
^{2+}dissociation constant for calmodulin *k*_{bCQ}- Ca
^{2+}dissociation constant for calsequestrin *k*_{bEGTA}- Ca
^{2+}dissociation constant for EGTA *k*_{bTC}- Ca
^{2+}dissociation constant for the troponin-Ca site *k*_{bTMC}- Ca
^{2+}dissociation constant for the troponin-Mg site *k*_{bTMM}- Mg
^{2+}dissociation constant for the troponin-Mg site *k*_{fBAPTA}- Ca
^{2+}association constant for BAPTA *k*_{fCM}- Ca
^{2+}association constant for calmodulin *k*_{fCQ}- Ca
^{2+}association constant for calsequestrin *k*_{fEGTA}- Ca
^{2+}association constant for EGTA *k*_{fTC}- Ca
^{2+}association constant for troponin *k*_{fTMC}- Ca
^{2+}association constant for the troponin-Mg site *k*_{fTMM}- Mg
^{2+}association constant for the troponin-Mg site - [TC]
_{tot} - Total concentration of the troponin-Ca site
- [TMC]
_{tot} - Total concentration of the troponin-Mg site

## THEORY AND METHODS

### Model Development

We formulated a mathematical model describing the dynamic properties of a single primary pacemaker cell of the rabbit SA node under space-clamp conditions (at 37°C). Our model is an extension of previous SA node models (17, 24, 124, 130) that utilized classical Hodgkin-Huxley formalism, including variations in intracellular ion concentrations, Ca^{2+} handling by the SR, and Ca^{2+} buffering. The standard model for normal pacemaker activity is described as a nonlinear dynamical system of 27 simultaneous, first-order, ordinary differential equations. A complete list of the equations and standard parameter values is presented in the
.

#### 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^{2+} 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).

There is now evidence that there exists a small restricted subsarcolemmal domain where Ca^{2+} concentrations may transiently reach higher levels than in the bulk myoplasm (29,47, 62, 86, 88, 122). Therefore, we assumed the subsarcolemmal space as a barrier for Ca^{2+} diffusion to the myoplasm. The subspace volume V_{sub} was tentatively set to 1% of the cell volume (0.035 pl), assuming that the subspace is limited to 20 nm below the sarcolemmal membrane (see Fig.1
*A*). The time constant of the Ca^{2+} diffusion from subspace to myoplasm (τ_{dif,Ca}) was set to 40 μs; this value corresponded to a diffusion coefficient of ≈1 × 10^{−9}cm^{2}/ms if the mean distance between the openings of the L-type Ca^{2+} channel or SR Ca^{2+} release channel and the myoplasmic compartment is assumed to be 200 nm (see Refs.29, 61, and 113). A τ_{dif,Ca} value of 40 μs was chosen to yield a peak [Ca^{2+}]_{sub} of ∼1.8 μM, as reported in the modeling studies of Glukhovsky et al. (29) and Snyder et al. (113). The subspace of our model corresponds to the “fuzzy” space, a functionally restricted intracellular space accessible to the Na^{+}/Ca^{2+} exchanger as well as to the L-type Ca^{2+} channel and Ca^{2+}-gated Ca^{2+} channel in the SR (see Refs. 51 and62). Note that it is different from the “diadic” space between the L-type Ca^{2+} channel and the ryanodine receptor (Ca^{2+}-gated Ca^{2+} channel), where Ca^{2+} concentrations increase to several tens of micromolars; in ventricular myocytes, the ryanodine receptor and SR Ca^{2+} release are closely linked to the L-type Ca^{2+} channel and its Ca^{2+}-dependent inactivation in the diadic space (1, 5, 50, 93, 106, 114,115). In SA node cells, however, such a microdomain or cross-signaling between the L-type Ca^{2+} channel and SR Ca^{2+} release channel has not been demonstrated, being possibly absent because the SR is poorly developed in SA node cells compared with ventricular or atrial myocytes (see Refs. 17and 47). Thus we did not incorporate the diadic space into our model.

#### 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
Equation 1where *I*
_{Ca,L} and*I*
_{Ca,T} represent the L-type and T-type Ca^{2+} 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^{2+}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).

Model equations for channel gating behaviors are essentially the same as those of the previous Hodgkin-Huxley type models. A gating variable,*x*, can be computed as a solution of the first-order differential equation of the form
Equation 2where *x*
_{∞} and τ_{x} are the steady-state *x* value and relaxation time constant, respectively, as functions of *V*. Relaxation time constants have been appropriately scaled for the temperature of 37°C with the use of a Q_{10} of 1.6∼3.0. The functions *x*
_{∞}(*V*) and τ_{x}(*V*) for individual gating variables are provided in the
(Fig. 16).

#### L-type Ca^{2+}* channel current.*

The kinetics of *I*
_{Ca,L} are described with activation (*d*
_{L}), voltage-dependent inactivation (*f*
_{L}), and Ca^{2+}-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^{2+}-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^{2+}-dependent inactivation process is independent of the voltage-dependent one (see Refs. 35 and105); the inactivation process was described by the two Hodgkin-Huxley type gating variables *f*
_{L} and*f*
_{Ca}.

The voltage dependences of *I*
_{Ca,L} activation and inactivation (steady-state probabilities and time constants for*d*
_{L} and *f*
_{L}) are shown in Fig. 2, *top*. For the steady-state activation and inactivation curves (*d*
_{L,∞} and *f*
_{L,∞}), we used the formulas of Demir et al. (17) based on the data from Fermini and Nathan (27). Expressions of the activation time constant τ_{dL} were adopted from Demir et al. (17), who modified the formulation of Nilius (84). To formulate the time constant of the voltage-dependent inactivation/recovery (τ_{fL}), the previous models employed different equations with different voltage dependences; however, these equations did not fit the recovery time course of*I*
_{Ca,L} experimentally observed in single rabbit cardiac myocytes (see Fig. 2). Therefore, we originally formulated the inactivation time constant τ_{fL}from the data of Nakayama et al. (81), Hagiwara et al. (37), and Kawano and Hiraoka (54). According to Demir et al. (17), a Q_{10} factor of 2.3 was applied to scale the gating time constant data for a temperature of 37°C. The inactivation/recovery time constant data were fitted to a function similar to that used by Lindblad et al. (68) for a rabbit atrial model and by Nygren et al. (88) for a human atrial model, using a least-square minimization procedure. Formulas for the Ca^{2+}-dependent inactivation*f*
_{Ca} are similar to those used by DiFrancesco and Noble (22). We included the subsarcolemmal domain and modeled the Ca^{2+}-dependent inactivation as a function of the subspace Ca^{2+} concentration ([Ca^{2+}]_{sub}). The half-maximum Ca^{2+} concentration for the Ca^{2+}-dependent inactivation was set to 0.35 μM, the same value as used by Courtemanche et al. (15). To reproduce the Ca^{2+}-dependent inactivation with a time constant of ∼10 ms during voltage-clamp pulses, the rate constant of Ca^{2+}binding was set to 60 ms^{−1} · mM^{−1}.

Maximum *I*
_{Ca,L} has been formulated as a fully selective Ca^{2+} current, with its reversal potential (*E*
_{Ca,L}) fixed at +45 mV and the maximum conductance *g*
_{Ca,L} set to 0.58 nS/pF at 2 mM [Ca^{2+}]_{o}. As Dokos et al. (24) have pointed out, *I*
_{Ca,L} has a relatively low*E*
_{Ca,L} albeit its high Ca^{2+}selectivity. Wilders et al. (124) adopted the constant-field formulation to describe the conductance property of*I*
_{Ca,L} with a low *E*
_{Ca,L}; as suggested by Dokos et al. (24), however, their formulas introduce large unnecessary fluxes of Na^{+} and K^{+}, leading to an overestimation of*I*
_{NaK}. The *E*
_{Ca,L} has been reported to shift with changing [Ca^{2+}]_{o} in accordance with an ideal Nernstian Ca^{2+} electrode, although displaced negatively below *E*
_{Ca} (see Ref.11); thus Dokos et al. (24) formulated *I*
_{Ca,L} as a fully selective Ca^{2+} current with *E*
_{Ca,L} displaced at a constant 75 mV negative to *E*
_{Ca}. The negative displacement of *E*
_{Ca,L} may be attributable to the existence of a high-concentration Ca^{2+} domain at the intracellular surface of the L-type Ca^{2+} channel pore. Nevertheless, there are no available data on the [Ca^{2+}]_{i} dependence of*E*
_{Ca,L}. It has been suggested that once the L-type Ca^{2+} channel is open, the Ca^{2+}concentration at the inner mouth does not change; the Ca^{2+}concentration in the vicinity of the inner mouth of the channel may be nearly constant (see Refs. 50, 111, and128). In our model, therefore,*E*
_{Ca,L} was set to a constant value of +45 mV, as in the model of Demir et al. (17) utilizing a fully Ca^{2+}-selective formulation of *I*
_{Ca,L}with the *E*
_{Ca,L} fixed at +46.4 mV. It might be difficult to express *E*
_{Ca,L} (or permeation kinetics) as a simple elementary function of ion concentrations, because the mechanisms of ion permeation in the L-type Ca^{2+}channel have been shown to involve complex ion-ion or ion-channel interactions (e.g., Refs. 41, 76, and116).

The model-generated *I*
_{Ca,L} during voltage-clamp pulses and the peak *I*
_{Ca,L}-*V*relationship are depicted in Fig. 2, *bottom*. Our model can simulate the Ca^{2+}-dependent inactivation of*I*
_{Ca,L} experimentally observed in rabbit SA node cells (10, 67); inactivation time courses of the simulated currents are very similar to the experimental data reported by Hagiwara et al. (37), Honjo et al. (44), Verheijck et al. (120), and Vinogradova et al. (122). The simulated peak *I*
_{Ca,L}-*V* relation is comparable to the experimental data as shown by the symbols. The maximum amplitude of *I*
_{Ca,L} measured in the simulated voltage-clamp experiment was 9.44 pA/pF, attained at 0 mV; this value is in good agreement with experimental data (44, 81,119, 120, 122).

#### T-type Ca^{2+}* 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}.

The voltage dependences of *I*
_{Kr} activation and inactivation (steady-state probabilities and time constants for the activation gating variable *p*
_{a} and inactivation gating variable *p*
_{i}) are shown in Fig.3, *top*, along with those for the previous models. To describe the gating kinetics of*I*
_{Kr}, the previous models (17, 24,124) used the equations provided by Shibasaki (108). However, Ono and Ito (90) recently reported the complete quantitative data on the activation kinetics of*I*
_{Kr} in single rabbit SA node cells and mathematically described the activation kinetics with two gating variables, *p*
_{aF} and *p*
_{aS}. According to their report, therefore, we described the general activation variable *p*
_{a} as a weighted sum of the fast (*p*
_{aF}) and slow (*p*
_{aS}) activation variables and used their original expressions for*I*
_{Kr} activation kinetics (*p*
_{a,∞}, τ_{paF}, and τ_{paS}). A modified version of the formulation of Ono and Ito (90) for*I*
_{Kr} activation was also utilized by Zhang et al. (130). We also employed the expression of Ono and Ito (90) for steady-state inactivation (*p*
_{i,∞}). No detailed experimental data are available on the time constant of the voltage-dependent*I*
_{Kr} inactivation (τ_{pi}); thus we adopted the expression of Shibasaki (108) for τ_{pi}.

The conductance of *I*
_{Kr} was chosen as to allow a maximum diastolic potential (MDP) between −60 and −55 mV (approximately equal to −58 mV) to be achieved. The standard*g*
_{Kr} value of our model is smaller than that determined by Ono and Ito (90); this difference in*g*
_{Kr} may reflect the regional difference in*I*
_{Kr} density between the center and periphery of the SA node (see Ref. 56). We also included a term to describe the “square root” activation of *I*
_{Kr}single channel conductance by [K^{+}]_{o}, as reported by Shibasaki (108). Dokos et al. (24) pointed out that the reversal potential of*I*
_{Kr} (*E*
_{Kr}) is positive to*E*
_{K} by 10∼19 mV, thus assumed that*I*
_{Kr} channels are slightly permeable to Na^{+}. According to recently published reports (e.g., Refs.90 and 119), however, *E*
_{Kr} is nearly equal to *E*
_{K}; therefore, the*I*
_{Kr} channel was assumed to be highly selective to K^{+}. The model-generated *I*
_{Kr}during the voltage-clamp pulses simulating the experiment of Ono and Ito (90) are shown in Fig. 3, *bottom*.

The kinetics of *I*
_{Ks} were described by the formulation of Zhang et al. (130). The steady-state activation curve for *I*
_{Ks} used in this study is from Lei and Brown (64). There are limited experimental data available for the time constant of the voltage-dependent activation of *I*
_{Ks} in rabbit SA node cells. Thus the time constant of *I*
_{Ks} activation was described using the expressions of Zhang et al. (130), i.e., the equations of Heath and Terrar (40) based on their data from guinea pig ventricular myocytes.

#### 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}.

The steady-state activation and inactivation curves (*r*
_{∞} and *q*
_{∞}) were based on the experimental data from Honjo et al. (46) and Lei et al. (65). The time constant of the voltage-dependent activation (τ_{r}) was formulated from the data of Giles and van Ginneken (28) for rabbit crista terminalis cells; the inactivation time constant (τ_{q}) was from Giles and van Ginneken (28) and Honjo et al. (46). To correct the data collected at 20.5∼25°C for 37°C, Zhang et al. (130) used a Q_{10} of 2.18. However, the use of Q_{10} = 2.18 yielded a pronounced phase 1 notch in APs. Thus we slightly accelerated the gating kinetics by using a Q_{10} of 3.0. The corrected time constants were comparable to those reported by Honjo et al. (46) and Uese et al. (117).

The values of the scaling parameters (conductances) for*I*
_{to} and *I*
_{sus} were determined by exploring the change in peak overshoot potential (POP) and variation of [K^{+}]_{i} . The experimentally measured densities of *I*
_{to} and*I*
_{sus} were significantly correlated with*C*
_{m}, i.e., cell size, and are larger in cells with a higher *C*
_{m} (46, 65); this probably reflects the regional differences of the current densities (see Ref. 7). We selected relatively small conductance values for our primary pacemaker cell model with relatively small*C*
_{m}, because the current densities in the central SA node region would be relatively small (see Refs. 64 and130). Whereas Zhang et al. (130) set the*I*
_{sus} conductance to a very small value of 3.3 pS/pF for their central model, we used a larger value of 20 pS/pF to accentuate the prolongation of AP duration during the blockade of the 4-AP-sensitive currents, i.e., to reproduce the experimental data of Boyett et al. (7).

#### 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_{10} = 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^{2+} channel blockers as well as by external Ca^{2+} and Mg^{2+}, and enhanced by isoprenaline or a Ca^{2+} agonist, BAY K 8644. These biophysical and pharmacological characteristics are compatible with those of the monovalent cation conductance of the L-type Ca^{2+} channel; thus Guo et al. (32) concluded that*I*
_{st} is generated by a novel subtype of the L-type Ca^{2+} 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}.

The voltage dependences of *I*
_{st} activation and inactivation (steady-state probabilities and time constants for the activation gating variable *q*
_{a} and inactivation gating variable *q*
_{i}) are shown in Fig.4, *top*. Because there is no detailed report on the gating kinetics (time constant) of*I*
_{st} in the rabbit SA node, we modified the formulation of Shinagawa et al. (109) for the rat SA node*I*
_{st}. From the data of Guo et al. (32), the half-activation voltage and slope factor for the activation curve were set to −57.0 and 5.0 mV, respectively. The inactivation and recovery of the rabbit *I*
_{st}, as reported in Guo et al. (32), appear to be slower than those of the rat *I*
_{st} reported by Shinagawa et al. (109). Thus the time constants of inactivation and recovery of the rabbit *I*
_{st} were assumed to be 6.65 times larger than those of the rat *I*
_{st}. The maximum conductance and reversal potential of*I*
_{st} were set to 15.0 pS/pF and +37.4 mV, respectively, according to Guo et al. (32). As shown in Fig. 4, *bottom*, our equations for *I*
_{st}could reproduce the kinetic properties of *I*
_{st}reported by Guo et al. (32).

#### 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_{10} 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.

*I*
_{Na} appears to be completely absent or almost negligible in primary pacemaker cells located in the central region of the SA node (57, 59, 82, 118, 130). Baruscotti et al. (3, 4) reported that both the density and frequency of occurrence of *I*
_{Na} in the SA node decrease with development, with *I*
_{Na} in adult rabbit SA node cells being very small or negligible. Thus *I*
_{Na}would not play an important role in normal pacemaking of primary pacemaker cells. The central model of Zhang et al. (130) for the leading pacemaker cell did not include*I*
_{Na}. In our simulations for primary pacemaker (central SA node) cells, therefore, *I*
_{Na} was assumed to be zero (negligible or completely inactivated).

#### Na^{+}*-dependent background current and muscarinic K*^{+} channel current.

^{+}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).

The background Ca^{2+} current (*I*
_{b,Ca}) has been frequently incorporated in mathematical models of the SA node to balance the Ca^{2+} extrusion via Na^{+}/Ca^{2+} exchange during diastole: the models of Wilders et al. (124), Demir et al. (17), and Zhang et al. (130) included*I*
_{b,Ca} with a conductance of 0.66∼1.25 pS/pF. In our model, however, incorporating an *I*
_{b,Ca} of 0.66∼1.25 pS/pF attenuated pacemaker activity (decreased APA to <70 mV) and induced Ca^{2+} overload (diastolic [Ca^{2+}]_{i} > 0.3 μM) via unnecessary Ca^{2+} influx, unfavorable for improving the model. Although the sarcolemmal Ca^{2+} pump current (*I*
_{p,Ca}) may balance*I*
_{b,Ca}, the contribution of*I*
_{p,Ca} to Ca^{2+} efflux through the SA node cell membrane is unknown. *I*
_{b,Ca} has not yet been recorded directly in SA node cells; the report of Hagiwara et al. (38) suggests that the rabbit SA node membrane is impermeable to Ca^{2+}. In our modeling, therefore,*I*
_{b,Ca} (as well as *I*
_{p,Ca}) was assumed to be negligible.

#### Na^{+}*-K*^{+}pump current.

^{+}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*^{2+}exchange current.

^{2+}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^{2+} 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^{2+}]_{i} and for negative potential-low [Ca^{2+}]_{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^{2+}]_{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^{2+}]_{i} = 0.5 μM, giving a diastolic free Ca^{2+} level in the presumed physiological range of <0.3 μM.

#### SR functions.

The SR was modeled as consisting of two compartments: the Ca^{2+} uptake store (network SR) and release store (junctional SR), as shown in Fig. 1
*B*. Owing to the lack of available data for updating the kinetic formulation of Ca^{2+}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^{2+}buffering by calsequestrin. Their formulation is, however, overly complex, given the complete lack of associated data in SA node cells; in their model, Ca^{2+} concentrations in the SR are unreasonably high compared with those in other models, much higher than experimental values for ventricular myocytes (see Refs. 63and 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^{2+} 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.

Because no experimental data are available to estimate the kinetic parameters for SR Ca^{2+} uptake and release in SA node cells, the parameter values were determined from previous modeling studies. The formula for the Ca^{2+} release mechanism was adopted from Dokos et al. (24), with the kinetics of Ca^{2+}release to subspace expressed as a function of [Ca^{2+}]_{sub} (not [Ca^{2+}]_{i}). The values of*K*
_{rel} used in previous models ranged from 0.8 to 2.0 μM; we chose a medium value of 1.2 μM. The*P*
_{rel} value was determined so as to make the peak [Ca^{2+}]_{sub} and [Ca^{2+}]_{i} nearly maximum with a*K*
_{rel} of 1.2 μM. The formulation for Ca^{2+} uptake was adopted from Luo and Rudy (71). The values of τ_{tr} used in previous SA node models ranged from 6.64 to 400 ms. In our modeling, the τ_{tr} value was set to a medium value of 60 ms, although some reports have suggested a smaller τ_{tr} or single compartment for SR (e.g., Refs. 29 and 113).

#### Ca^{2+}* buffering.*

Our model includes the dynamics of three Ca^{2+} buffers: calmodulin, troponin, and calsequestrin. The amounts of calmodulin and troponin available for Ca^{2+} binding, and the rate constants for Ca^{2+} 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^{2+} binding to calmodulin and troponin were scaled with the use of a Q_{10}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^{2+} binding to calsequestrin were based on the study of Cannell and Allen (12), adjusted to 37°C via a Q_{10} of 1.6 (see Ref. 68).

We also examined the effects of exogenous Ca^{2+} buffers, BAPTA and EGTA, on pacemaker activity. The rate constants for Ca^{2+} binding to and release from these buffers were adopted from Sham (106), Smith et al. (112), and You et al. (128).

#### Ion concentration homeostasis.

Our model also includes material balance expressions to define the temporal variations in [Na^{+}]_{i}, [K^{+}]_{i}, and [Ca^{2+}]_{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^{2+}. 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
, 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^{2+}]_{o}were set equal to 140, 5.4, and 2 mM, respectively.

### Numerical Integration (Dynamic Simulation)

Dynamic behavior of the model cell was determined by solving the simultaneous nonlinear ordinary differential equations numerically. We employed a fourth-order adaptive-step Runge-Kutta algorithm, which includes an automatic step-size adjustment based on an error estimate, or a variable time-step numerical differentiation approach selected for its suitability to stiff systems. The maximum relative error tolerance for our integration methods was set to 1 × 10^{−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^{−3} 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

### 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
. The MDP, POP, and APA are −58.6, +16.6, and 75.2 mV, respectively; the CL and AP duration at 50% repolarization (APD_{50}) 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_{50}. 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.

As shown in Fig. 6, our model-generated *I*
_{Ca,L}waveform during spontaneous APs was very similar to that recorded by Zaza et al. (129) and Doerr et al. (23) during the “AP clamp” of a single rabbit SA node cell. The model-generated *I*
_{Kr} waveform was also very similar to that measured as an E-4031-sensitive current during AP clamp (90, 129). Furthermore, the time course of the simulated*I*
_{Ca,T} was similar to that of*I*
_{Ca,T} measured as a Ni^{2+}-sensitive current by Doerr et al. (23); the*I*
_{h} during pacemaker depolarization was similar to the Cs^{+}-sensitive current recorded by Zaza et al. (129) and also to the *I*
_{h} computed by Maruoka et al. (73). Thus the simulated changes in *I*
_{Ca,L}, *I*
_{Kr},*I*
_{Ca,T}, and *I*
_{h} during pacemaking are comparable to the experimentally observed changes in these currents during spontaneous APs in rabbit SA node cells.

Simulated ionic currents, as well as an AP waveform, in our model were compared with those in the previous models (Figs. 6 and7). The time courses of*I*
_{Ca,L} and *I*
_{Kr} during pacemaker activity in the models were different. In the models of Wilders et al. (124), Demir et al. (17), and Dokos et al. (24), the fast inhibition of*I*
_{Ca,L} was followed by a secondary increase (inward “hump”); in contrast, the inward hump was absent in our simulated *I*
_{Ca,L}. Consistent with our model simulation, the secondary increase in *I*
_{Ca,L} was not observed in the AP-clamp experiments (23, 129). The rapid inactivation of *I*
_{Kr} on AP upstroke (during phase 0) observed in the AP-clamp experiments (90, 129) was reproduced only by our model and the model of Zhang et al. (130). Thus our model is apparently superior to previous models in generating the *I*
_{Ca,L} and*I*
_{Kr} waveforms experimentally observed in single rabbit SA node cells. The amplitudes of *I*
_{b,Na}and *I*
_{K,ACh} (as *I*
_{b,K}) during pacemaker activity of our model are comparable to those of previous models (see Fig. 7). On the other hand, the time courses of*I*
_{NaCa} in the model cells were quite different, reflecting the different intracellular Ca^{2+} dynamics. In our model, *I*
_{NaCa} during phase 4 was relatively large, suggesting the large contribution of*I*
_{NaCa} to pacemaker depolarization.

#### Intracellular Ca^{2+}* dynamics.*

Intracellular Ca^{2+} dynamics (concentration changes, SR Ca^{2+} uptake and release, and Ca^{2+} buffering) during normal pacemaking are shown in Fig. 5, *right*. The free Ca^{2+} concentrations in subspace ([Ca^{2+}]_{sub}) and myoplasm ([Ca^{2+}]_{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^{2+} concentrations in the rabbit SA node cell, these values are comparable to those measured in other pacemaker cells (47, 53).

As shown in Fig. 8, the simulated intracellular Ca^{2+} dynamics in our model cell were quite different from those in previous models, because we assumed a subsarcolemmal space and relatively small SR volumes. The myoplasmic Ca^{2+} transient in our model, the peak value of which is 0.68 μM, was much smaller than that in previous SA node models but is comparable to that experimentally or theoretically determined for atrial or ventricular myocytes (e.g., Refs. 68,71, and 110). In contrast, peak [Ca^{2+}]_{i} values in previous SA node models were >2.5 μM, unreasonably larger than those in atrial or ventricular models. The Ca^{2+} concentrations in the SR ([Ca^{2+}]_{rel} and [Ca^{2+}]_{up}) of our model are similar to those of previous models except for the model of Demir et al. (17), in which the SR Ca^{2+} concentrations (>10 mM) are much higher than experimental values for ventricular myocytes (see Refs. 63 and 107). The Ca^{2+} release from the SR in our model is much smaller than that in the model of Wilders et al. (124) or Demir et al. (17), whereas it is comparable to that in the model of Dokos et al. (24), suggesting a relatively small contribution of SR Ca^{2+} release to intracellular Ca^{2+} transients in primary pacemaker cells.

#### 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.

### Effects of Modulating Sarcolemmal Ionic Currents on Pacemaker Activity

We simulated the effects of inhibition of sarcolemmal ionic currents on pacemaker activity and then compared the simulated behaviors of our model cell with experimental findings as well as with those of previous models. We used two different systems, the*I*
_{st}-incorporated (standard) system and the*I*
_{st}-removed (*g*
_{Kr} = 0.75) system, because *I*
_{st} is not always present in spontaneously beating SA node cells (see Ref. 120).

#### 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}). Figure10
*A* 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.

Dynamic behaviors of real SA node cells, as well as the model cell, during inhibition of *I*
_{Ca,L} are shown in Fig.10
*B*. In experiments using Ca^{2+} antagonists (L-type Ca^{2+} channel blockers), such as nifedipine and verapamil, two ways of abolishing pacemaker activity have been observed: one is characterized by a gradual and continuous decline in APA (with a marked decrease in POP) to quiescence as well as an increase in CL (57, 58) and the other is characterized by irregular (chaotic) dynamics, called “skipped-beat runs,” as observed in the rabbit SA node (122) or cat subsidary pacemaker cells (97). Our model could simulate both of the two distinct bifurcation structures experimentally observed during applications of Ca^{2+} antagonists. In the experiment for small balls of rabbit SA node tissues, nifedipine (2 μM) abolished spontaneous APs in the center of the SA node (i.e., the leading pacemaker site); the resting potential at which SA node cells settled after block of *I*
_{Ca,L} was −40 ± 5 mV (57). The resting potential in this experiment was ∼10 mV more negative than the prediction of our standard model (approximately −29.5 to −32.6 mV), whereas it was close to the value in the *I*
_{st}-removed system (approximately −37.6 to −42.0 mV). This inconsistency may be due in part to a concomitant block of *I*
_{st} by the Ca^{2+} channel blocker, which is known to inhibit *I*
_{st} (see Refs. 32 and 78), as well as to the absence of *I*
_{st} in some preparations; the complete block of both *I*
_{Ca,L} and *I*
_{st} in the standard model yielded a resting potential of −43.7 mV.

The effects of blocking *I*
_{Ca,L} (decreasing*g*
_{Ca,L}) on spontaneous APs of previous SA node models are shown in Fig. 11 for comparison. In all of the models, blocking *I*
_{Ca,L}caused cessation of pacemaker activity; resting potentials of the quiescent model cells at critical *g*
_{Ca,L} values and with complete block of *I*
_{Ca,L} are compared in Table 2. The simulated behaviors of previous models during *g*
_{Ca,L} decrease are inconsistent with experimental observations (30, 57, 58,122) in the following respects: *1*) the model of Dokos et al. (24) did not exhibit the gradual decline in APA (POP) observed in experiments, with pacemaker activity being abruptly abolished via irregular dynamics; *2*) in the central model of Zhang et al. (130), pacemaker activity abruptly ceased with no irregular dynamics; *3*) in the models of Demir et al. (17) and Zhang et al. (130), CL monotonously decreased with reducing *g*
_{Ca,L}, whereas pacemaking slowed during applications of Ca^{2+} antagonists; and *4*) the resting potentials predicted by the models of Wilders et al. (124), Demir et al. (17), and Zhang et al. (130) were 5∼15 mV more positive than those experimentally determined by Kodama et al. (57) and predicted by our *I*
_{st}-removed system. Thus our model is superior to previous models in predicting changes in AP parameters, bifurcation structures, and resting potentials during applications of Ca^{2+} antagonists.

#### 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.

Selective *I*
_{Kr} blockers such as E-4031 have been experimentally found to cause a gradual and continuous decline in APA with marked depolarization in MDP and then quiescence (see Fig. 12,*right*). These dynamic changes during applications of*I*
_{Kr} blockers are in good agreement with the predictions of our model cell. Our model predicted only a slight increase or even decrease in CL during the *g*
_{Kr}decline. In most experiments, however, CL was increased by applications of *I*
_{Kr} blockers (e.g., see Refs.56, 64, and 119), although 0.1 μM E-4031 did not significantly change CL in Ono and Ito (90). This inconsistency may result from the state-dependent kinetics of *I*
_{Kr} block. Alternatively, *I*
_{Kr} blockers may also modulate other ionic currents directly or secondary via changing intracellular ion concentrations. In recent experiments, block of*I*
_{Kr} by 1 μM E-4031 caused a cessation of spontaneous activity; after block of *I*
_{Kr}, the*V* of rabbit SA node cells or tissues settled at −32 ± 2 mV (56), −37.4 ± 2.9 mV (90), or −24.5 ± 1.8 mV (119). These experimentally observed resting potentials after block of *I*
_{Kr} were 10∼20 mV more negative than the predictions of our model. In other reports, however, the resting or zero current potential after block of*I*
_{Kr} was relatively positive: *1*) Yanagihara and Irisawa (126) reported a cessation of pacemaker activity with a resting potential of −10 mV on application of Ba^{2+} to block *I*
_{K} and*2*) in the report of Verheijck et al. (119), 10 μM E-4031 stabilized *V* at −19.6 ± 1.8 mV and yielded a zero current crossing of approximately −10.0 mV in the*I*-*V* curve. These experimental data are in reasonable agreement with the predictions of our model. The relatively deep resting potentials in some experiments may be due to *1*) the regional difference between the central (primary pacemaker) cell and the transitional or peripheral (latent pacemaker) cell (e.g., in the density of background currents); *2*) the difference in recording methods or other experimental conditions (when [Na^{+}]_{i} and [K^{+}]_{i}were fixed at constant values of 10 and 140 mM, respectively, as in the perforated-patch recording for single SA node cells, the resting potential of our standard system was −18.5 mV at*g*
_{Kr} = 0.11 and −17.5 mV with the complete block of *I*
_{Kr}); *3*) a Ca^{2+}-dependent increase of background Cl^{−} or K^{+} conductance in the region of*I*
_{Ca,L} window current (e.g., see Ref.92); or *4*) direct modifications of ionic currents other than *I*
_{Kr} by*I*
_{Kr} blockers.

The effects of block of *I*
_{Kr} (decreasing*g*
_{Kr}) on spontaneous APs of previous SA node models are shown in Fig. 13 for comparison. During the decrease in *g*
_{Kr}, pacemaker activity of previous models abruptly ceased or attenuated with MDP only slightly depolarized; in the model of Dokos et al. (24), irregular dynamics appeared before the cessation of pacemaker activity. In experiments using *I*
_{Kr}blockers such as E-4031, however, application of an*I*
_{Kr} blocker always induced the gradual and continuous decline in APA to quiescence with a marked depolarization in MDP and no irregular dynamics, as simulated by our model (see Fig. 12). Abrupt cessation of pacemaker activity or irregular dynamics as predicted by previous models has never been observed during applications of *I*
_{Kr} blockers. Thus the bifurcation structures during *I*
_{Kr} inhibition of real SA node cells are different from those of previous models, whereas they are essentially the same as those of our model. Block of*I*
_{Kr} caused a cessation of pacemaker activity in all of the models; the resting potentials at critical*g*
_{Kr} values and with the complete block of*I*
_{Kr} for SA node models are compared in Table 2, along with the experimental data. Resting potentials in the models of Wilders et al. (124), Demir et al. (17), and Dokos et al. (24) were more positive than those predicted by our model as well as those observed in the experiments. Taken together, our model appears to be superior to previous models in predicting the behavior of real SA node cells during applications of *I*
_{Kr} blockers.

#### 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^{2+} exerted a pronounced negative chronotropic effect on pacemaker activity by nearly doubling the basal CL (see also Refs. 84 and103). In contrast, Wilders et al. (124) reported that the Ni^{2+}-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.

The discrepancy in the apparent contribution of*I*
_{Ca,T} to pacemaker activity may in part result from the heterogeneity of SA node cells (regional difference) in the density of *I*
_{st} as well as of*I*
_{Na} or *I*
_{Ca,T} itself. In most experiments, Ni^{2+} was used as the selective*I*
_{Ca,T} blocker. However, Ni^{2+} has been found to block *I*
_{st} as well: Ni^{2+} at 40 μM, which blocks *I*
_{Ca,T}completely, also decreases the amplitude of *I*
_{st}to 46∼64% of the control; Ni^{2+} at 1 mM completely abolishes *I*
_{st} (32, 78). As shown in Fig. 9, the pacemaker activity of our model cell was dramatically slowed by reducing the *I*
_{st} conductance: the complete block of *I*
_{st} increased CL by 53.0% (from 307.5 to 470.5 ms). When a SA node cell possesses*I*
_{st} sensitive to Ni^{2+}, therefore, the effect of block of *I*
_{Ca,T} would be overestimated because of a concomitant block of*I*
_{st} by Ni^{2+}. Thus the differences in the Ni^{2+}-induced CL prolongation between experimental studies may reflect the difference in *I*
_{st}density. The Ni^{2+}-sensitive current recorded by Doerr et al. (23) during pacemaker depolarization, considerably larger (∼20 pA/cell) than *I*
_{Ca,T} in our model (∼7 pA/cell), may also include *I*
_{st}.

#### 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.

The simulated effect of elimination of *I*
_{h} in our standard model for the primary pacemaker cell (8.4% increase in CL) is in good agreement with experimental data from central SA node tissues: Nikmaram et al. (83) reported that in central SA node tissues, block of *I*
_{h} by 2 mM Cs^{+}caused a 5∼10% increase in CL (see also Ref. 130). In other previous reports, the effects of block of i_{h} on CL were greater than in the report of Nikmaram et al. (83): Denyer and Brown (18) and van Ginneken and Giles (118) reported increases in CL of 19∼30% and 57%, respectively (see Table 3). Nikmaram et al. (83) also reported a regional difference in the effect of Cs^{+} on CL: the decrease in the spontaneous rate induced by 2 mM Cs^{+}was largest in the periphery (∼19%) and least in the center (∼7%). Because Honjo et al. (44) have shown that the density of *I*
_{h} in the rabbit SA node was significantly greater in larger cells from the periphery than in smaller cells from the center, the different effects of block of*I*
_{h} on CL may reflect the regional difference in*I*
_{h} density. Furthermore, the degree of CL prolongation on block of *I*
_{h} strongly depends on basic CL, MDP, and thus on densities of ionic currents other than*I*
_{h} itself; values of MDP and CL would vary depending on both the experimental condition and area of the SA node from which cells were isolated. As already mentioned, the effect of elimination of *I*
_{h} on CL was much greater in the*I*
_{st}-removed system with longer CL than in the standard system with shorter CL. When *I*
_{st} was eliminated from the standard system (MDP = −63.64 mV), the complete block of *I*
_{h} increased CL by 43.9% (from 470.5 to 677.0 ms). Thus the large differences in the effect of block of *I*
_{h} on CL (different contributions of*I*
_{h} to pacemaker depolarization) reported previously (18, 20, 83, 87, 118) would arise from the regional variations in MDP and basic CL as well as in the density or activation threshold of *I*
_{h}.

#### 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_{50}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_{50} (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^{2+}* Dynamics on Pacemaker Activity*

#### Effects of block of SR Ca^{2+}* release.*

SR Ca^{2+} release is known to affect intracellular Ca^{2+} transients, Ca^{2+}-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^{2+} 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^{2+} 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^{2+} release decreased the peak [Ca^{2+}]_{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^{2+} release plays only a minor role in generating subsarcolemmal Ca^{2+} transients and regulating pacemaker activity (via modifications of the Ca^{2+}-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^{2+} release resulted in little or no slowing of pacemaker activity in SA node cells. The Ca^{2+}-dependent inactivation of L-type Ca^{2+}channels in rat ventricular or atrial myocytes is mediated primarily by Ca^{2+} release from the SR to the microdomain (5,60). In our SA node model, however, the Ca^{2+}-dependent inactivation of L-type Ca^{2+}channels is chiefly mediated by Ca^{2+} influx through the L-type Ca^{2+} 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^{2+} release does not contribute to the Ca^{2+}-dependent inactivation of L-type Ca^{2+}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^{2+} release to intracellular Ca^{2+} 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.

In some experimental reports (6, 39, 67, 104), the block of SR Ca^{2+} release by ryanodine has been shown to exert a negative chronotropic effect on rabbit SA node cells, suggesting the important roles of SR Ca^{2+} release, intracellular Ca^{2+} transients, and subsequent activation of*I*
_{NaCa} in regulating pacemaker activity of the SA node. This inconsistency may be due to the regional difference in the development of the SR between the central and peripheral SA node, i.e., poor development of the SR in primary pacemaker cells (see Refs.47 and 96). Alternatively, ryanodine may directly block sarcolemmal ionic currents, such as*I*
_{Ca,T} and *I*
_{NaCa}, or secondarily modify Ca^{2+}-dependent currents including*I*
_{Ca,L}, *I*
_{Kr}, and*I*
_{h} via affecting intracellular Ca^{2+}dynamics (see Refs. 6, 39, 67,94, and 104). Recent experiments have shown that the *I*
_{Ca,T}-triggered focal Ca^{2+}release from the junctional SR to subspace stimulates*I*
_{NaCa}, thereby playing an important role in regulating pacemaker activity (6, 47). Our model cannot simulate such local events (i.e., Ca^{2+} sparks and waves) due to the spatial heterogeneity of intracellular Ca^{2+}distribution; further compartmentalization for intracellular Ca^{2+} or the use of partial differential equations is required for simulating these phenomena involved in pacemaker regulation.

#### Effects of Ca^{2+}* buffers.*

To buffer intracellular Ca^{2+}, Ca^{2+} buffers such as EGTA and BAPTA are routinely used for whole cell patch-clamp experiments. These Ca^{2+} buffers to reduce intracellular Ca^{2+} transients have been shown to affect the pacemaker activity of SA node cells (67, 122). We therefore examined the effects of the Ca^{2+} buffers BAPTA and EGTA on the pacemaker activity of SA node models. Wilders et al. (124) simulated EGTA buffering by fixing [Ca^{2+}]_{i} to 80 nM. However, it is known that the slow buffer EGTA cannot efficiently buffer free Ca^{2+} in the subsarcolemmal space: although both BAPTA and EGTA can efficiently suppress global (myoplasmic) Ca^{2+} transients, only the fast buffer BAPTA is efficient in buffering local (subsarcolemmal) Ca^{2+}transients (see Refs. 5, 106,122, and 128). Therefore, we calculated intracellular Ca^{2+} dynamics in the presence of BAPTA or EGTA to simulate the differential effects of Ca^{2+} buffers on [Ca^{2+}]_{sub} and on spontaneous APs. As shown in Fig. 15 and Table5, the fast buffer BAPTA at 10 mM remarkably reduced both the subsarcolemmal and myoplasmic Ca^{2+} 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^{2+}transient, whereas it almost completely suppressed the myoplasmic Ca^{2+} 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^{2+} buffers: both BAPTA and EGTA dramatically suppressed intracellular Ca^{2+} 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^{2+}, including the subsarcolenmal space as a diffusion barrier for Ca^{2+}, whereas previous models assumed only one intracellular compartment. Subspace Ca^{2+} dynamics are indispensable for simulating the differential effects of BAPTA and EGTA on SA node dynamics. A recent experimental report suggests local high Ca^{2+} gradients in the subsarcolemmal microdomain in SA node cells (122).

### Achievements and Limitations of the Present Model

On the basis of recently published experimental data, we were able to develop an improved SA node model incorporating *1*) the novel pacemaker current *I*
_{st} not included in previous models, *2*) new formulations for voltage- and Ca^{2+}-dependent inactivation kinetics of the L-type Ca^{2+} channel, *3*) new expressions for activation kinetics of *I*
_{Kr}, *4*) revised kinetic formulas for 4-AP-sensitive currents (*I*
_{to} and*I*
_{sus}), *5*) new formulations for voltage- and concentration-dependent kinetics of*I*
_{NaK}, and *6*) the subsarcolemmal space as a diffusion barrier for Ca^{2+}. The present model provides well-integrated explanations of the electrophysiological behavior of primary pacemaker cells in the rabbit SA node, representing significant improvements over earlier models. As described above, our model can*1*) simulate whole cell voltage-clamp data for*I*
_{Ca,L}, *I*
_{Kr}, and*I*
_{st}; *2*) reproduce the waveshapes of spontaneous APs and sarcolemmal ionic currents (*I*
_{Ca,L}, *I*
_{Kr},*I*
_{Ca,T}, and *I*
_{h}) observed during AP-clamp recordings; *3*) reproduce the bifurcation structures seen during applications of *I*
_{Ca,L} or*I*
_{Kr} blockers; *4*) mimic the reported effects of block of *I*
_{Ca,T},*I*
_{h}, or 4-AP-sensitive currents on pacemaker activity; and *5*) simulate the differential effects of BAPTA and EGTA on pacemaker frequency more accurately than previous SA node models.

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^{2+}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^{2+}-activated Cl^{−} current; see Ref. 92); *5*) the spatial heterogeneity of intracellular distributions of Ca^{2+} 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^{2+} concentrations and/or modifications of second messengers (see Refs. 9 and129). 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^{2+}, 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^{2+}-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.

## Acknowledgments

This work was supported in part by Ministry for Education, Science, Sports and Culture of Japan Grant-in-Aid for Scientific Research (c) 11670717 (to S. Imanishi) and by Kanazawa Medical University Grant for Project Research P99-6 (to Y. Kurata).

## Appendix

The mathematical expressions used in our SA node model are given below. The units used are millivolts, picoamps, nanosiemens, milliseconds, nanofarads, millimolars, and liters. The temperature assumed for the computations was 37°C; the experimental data of gating kinetics obtained at <37°C were corrected for temperature with a Q_{10} = 1.6∼3.0 (see theory and methods). The functions *x*
_{∞}(*V*) and τ_{x}(*V*) for individual gating variables and steady-state *I*-*V* relations for P individual current systems are plotted in Figs. A1 and A2, respectively; the model constants (standard parameter values) are given in Table A1.

### Sarcolemmal Ionic Currents

#### L-type Ca^{2+}* channel current.*

#### T-type Ca^{2+}* channel current.*

#### Rapidly activating delayed rectifier K^{+}* current.*

#### Slowly activating delayed rectifier K^{+}* current.*

#### 4-AP-sensitive currents.

#### Hyperpolarization-activated current.

#### Sustained inward current.

#### Na^{+}* channel current.*

#### Na^{+}*-dependent background current.*

#### Background muscarinic K^{+}* channel current.*

#### Na^{+}*-K*^{+}pump current.

^{+}pump current.

#### Na^{+}*/Ca*^{2+}exchange current.

^{2+}exchange current.

### Intracellular Ca^{2+}* Dynamics*

#### Ca^{2+}* diffusion flux.*

#### Ca^{2+}* handling by the SR.*

### State Variables and Differential Equations

#### Membrane potential.

#### Gating variables.

#### Intracellular ion concentrations.

#### Ca^{2+}* buffering.*

## Footnotes

Address for reprint requests and other correspondence: Y. Kurata, Dept. of Physiology, Kanazawa Medical University, 1-1 Daigaku, Uchinada-machi, Kahoku-gun, Ishikawa 920-0293, Japan (E-mail:yasu{at}kanazawa-med.ac.jp).

The costs of publication of this article were defrayed in part by the payment of page charges. The article must therefore be hereby marked “

*advertisement*” in accordance with 18 U.S.C. Section 1734 solely to indicate this fact.10.1152/ajpheart.00900.2001

- Copyright © 2002 the American Physiological Society