|
|
||||||||
1 Department of Physiology, Kanazawa Medical University, Ishikawa 920-0293; and 2 First Department of Internal Medicine, Tottori University School of Medicine, Yonago 683-0826, Japan
| |
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 (Ist) recently identified in primary pacemaker cells, 2) reformulation of voltage- and Ca2+-dependent inactivation of the L-type Ca2+ channel current (ICa,L), 3) new expressions for activation kinetics of the rapidly activating delayed rectifier K+ channel current (IKr), and 4) incorporation of the subsarcolemmal space as a diffusion barrier for Ca2+. 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 Ca2+ dynamics on pacemaker activity. Our model represents significant improvements over the previous models, because it can 1) simulate whole cell voltage-clamp data for ICa,L, IKr, and Ist; 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 Ca2+ buffers on pacemaker activity more accurately than the previous models.
rabbit sinoatrial node; nonlinear dynamical system; computer simulation; bifurcation diagram
| |
INTRODUCTION |
|---|
|
|
|---|
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 Ca2+ buffering (by calmodulin, troponin, and calsequestrin) and new formulations for Ca2+ 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+/Ca2+ exchanger current (INaCa) and muscarinic K+ channel current (IK,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 Ca2+ channel current (ICa,L) were not based on experimental data. Second, the Ca2+-dependent inactivation of ICa,L, known as an essential property of the L-type Ca2+ channel in the rabbit SA node (10, 67), was not formulated [or not directly dependent on the intracellular Ca2+ concentration ([Ca2+]i), but a function of the inactivation gating variable fL in Wilders et al. (124)]. Although Dokos et al. (24) incorporated a second inactivation gating variable to represent the [Ca2+]i-dependent inactivation process, they assumed an extraordinarily large maximum ICa,L conductance (gCa,L) and a very high affinity for Ca2+ binding to the inactivation site of L-type Ca2+ channels; such a large gCa,L with a high affinity for Ca2+ binding is rather unlikely and has not been demonstrated. Third, these models do not include the slow component of the delayed rectifier K+ current (IKs), 4-aminopyridine (4-AP)-sensitive currents, or sustained inward current (Ist), whereas these currents are known to be present in primary pacemaker cells. Fourth, intracellular Ca2+ 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 Ca2+ 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 IKs 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 Ca2+ during an AP cycle is not zero), whereas intracellular Ca2+ 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 Ist and IK,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 ICa,L or the rapid component of the delayed rectifier K+ current (IKr) 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) Ist, not included in the previous models, has been incorporated; 2) voltage- and Ca2+-dependent inactivation kinetics of ICa,L have been reformulated; 3) expressions for activation kinetics of IKr 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 (INaK) have been reformulated; and 6) the subsarcolemmal space as a diffusion barrier for intracellular Ca2+ 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 Ca2+ dynamics on pacemaker activity. The experimental findings for verification include 1) whole cell voltage-clamp data for ICa,L, IKr, and Ist; 2) waveshapes of spontaneous APs and ionic currents (such as ICa,L and IKr) as observed during AP-clamp recordings; 3) modulations and cessation of pacemaker activity by applications of ICa,L or IKr blockers; 4) modifications of AP waveforms (changes in AP parameters) by blocking the T-type Ca2+ channel current (ICa,T), hyperpolarization-activated current (Ih), or 4-AP-sensitive currents; 5) effects of blocking SR Ca2+ release (by ryanodine) on pacemaker frequency; and 6) negative chronotropic effects of Ca2+ 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 ICa,L or IKr blockers and also the differential effects of BAPTA and EGTA on pacemaker frequency. During inhibition of ICa,L or IKr, 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
| Ib,Ca | Background Ca2+ current |
| Ib,K | Background K+ current |
| Ib,Na | Background Na+ current |
| ICa,L | L-type Ca2+ channel current |
| ICa,T | T-type Ca2+ channel current |
| Ih | Hyperpolarization-activated current |
| IK | Delayed rectifier K+ current |
| IK,ACh | Muscarinic K+ channel current |
| IKr | Rapid component of the delayed rectifier K+ current |
| IKs | Slow component of the delayed rectifier K+ current |
| INa | Fast Na+ channel current |
| INaCa | Na+/Ca2+ exchanger current |
| INaen | Electroneutral Na+ influx current |
| INaK | Na+-K+ pump current |
| INaKmax | Maximum INaK |
| Ip,Ca | Sarcolemmal Ca2+ pump current |
| Ist | Sustained inward current |
| Isus | Sustained component of the 4-AP-sensitive current |
| Ito | Transient component of the 4-AP-sensitive current |
Cell Geometry
| Cm | Cell membrane capacitance |
| Vcell | Cell volume |
| Vi | Myoplasmic volume available for Ca2+ diffusion |
| Vrel | Volume of junctional SR (Ca2+ release store) |
| Vsub | Subspace volume |
| Vup | Volume of network SR (Ca2+ uptake store) |
Ionic Concentrations
| Ca2+i | Myoplasmic Ca2+ concentration |
| [Ca2+]o | Extracellular Ca2+ concentration |
| [Ca2+]rel | Ca2+ concentration in the junctional SR |
| [Ca2+]sub | Subspace Ca2+ concentration |
| [Ca2+]up | Ca2+ concentration in the network SR |
| [K+]i | Intracellular K+ concentration |
| [K+]o | Extracellular K+ concentration |
| [Mg2+]i | Intracellular Mg2+ concentration |
| [Na+]i | Intracellular Na+ concentration |
| [Na+]o | Extracellular Na+ concentration |
Equilibrium (Reversal) Potentials
| ECa,L | Apparent reversal potential of ICa,L |
| ECa,T | Apparent reversal potential of ICa,T |
| EK | Equilibrium (Nernst) potential for K+ |
| EKr | Reversal potential of IKr |
| EKs | Reversal potential of IKs |
| Emh | Reversal potential of INa |
| ENa | Equilibrium (Nernst) potential for Na+ |
| Est | Apparent reversal potential of Ist |
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 |
Ca2+ 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 |
Ca2+ Buffering
| CQtot | Total calsequestrin concentration |
| [CM]tot | Total calmodulin concentration |
| fCMi | Fractional occupancy of calmodulin by Ca2+ in myoplasm |
| fCMs | Fractional occupancy of calmodulin by Ca2+ in subspace |
| fCQ | Fractional occupancy of calsequestrin by Ca2+ |
| fTC | Fractional occupancy of the troponin-Ca site by Ca2+ |
| fTMC | Fractional occupancy of the troponin-Mg site by Ca2+ |
| fTMM | Fractional occupancy of the troponin-Mg site by Mg2+ |
| kbBAPTA | Ca2+ dissociation constant for BAPTA |
| kbCM | Ca2+ dissociation constant for calmodulin |
| kbCQ | Ca2+ dissociation constant for calsequestrin |
| kbEGTA | Ca2+ dissociation constant for EGTA |
| kbTC | Ca2+ dissociation constant for the troponin-Ca site |
| kbTMC | Ca2+ dissociation constant for the troponin-Mg site |
| kbTMM | Mg2+ dissociation constant for the troponin-Mg site |
| kfBAPTA | Ca2+ association constant for BAPTA |
| kfCM | Ca2+ association constant for calmodulin |
| kfCQ | Ca2+ association constant for calsequestrin |
| kfEGTA | Ca2+ association constant for EGTA |
| kfTC | Ca2+ association constant for troponin |
| kfTMC | Ca2+ association constant for the troponin-Mg site |
| kfTMM | Mg2+ 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, Ca2+ handling by the SR, and Ca2+ 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 APPENDIX.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 Ca2+ 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 Ca2+ 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 Ca2+ diffusion to the myoplasm. The subspace volume Vsub 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. 1A). The time constant of the Ca2+ diffusion from subspace to myoplasm (
dif,Ca) was set to 40 µs; this value corresponded to
a diffusion coefficient of
1 × 10
9
cm2/ms if the mean distance between the openings of the
L-type Ca2+ channel or SR Ca2+ 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
[Ca2+]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+/Ca2+ exchanger as well as
to the L-type Ca2+ channel and Ca2+-gated
Ca2+ channel in the SR (see Refs. 51 and
62). Note that it is different from the "diadic" space
between the L-type Ca2+ channel and the ryanodine receptor
(Ca2+-gated Ca2+ channel), where
Ca2+ concentrations increase to several tens of
micromolars; in ventricular myocytes, the ryanodine receptor
and SR Ca2+ release are closely linked to the L-type
Ca2+ channel and its Ca2+-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 Ca2+ channel and SR
Ca2+ 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. 17
and 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
|
(1) |
|
(2) |
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 Q10 of 1.6~3.0.
The functions x
(V) and
x(V) for individual gating
variables are provided in the APPENDIX (Fig. 16).
L-type Ca2+ channel current. The kinetics of ICa,L are described with activation (dL), voltage-dependent inactivation (fL), and Ca2+-dependent inactivation (fCa) gating variables. The inactivation of ICa,L consists of two exponential terms: rapid and slow components. The rapid component is mediated by the intracellular Ca2+-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 Ca2+-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 fL and fCa.
The voltage dependences of ICa,L activation and inactivation (steady-state probabilities and time constants for dL and fL) are shown in Fig. 2, top. For the steady-state activation and inactivation curves (dL,
and fL,
), 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
ICa,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 Q10 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 Ca2+-dependent inactivation
fCa are similar to those used by DiFrancesco and
Noble (22). We included the subsarcolemmal domain and
modeled the Ca2+-dependent inactivation as a function of
the subspace Ca2+ concentration
([Ca2+]sub). The half-maximum
Ca2+ concentration for the Ca2+-dependent
inactivation was set to 0.35 µM, the same value as used by
Courtemanche et al. (15). To reproduce the
Ca2+-dependent inactivation with a time constant of ~10
ms during voltage-clamp pulses, the rate constant of Ca2+
binding was set to 60 ms
1 · mM
1.
|
T-type Ca2+ channel current. Contributions of ICa,T to pacemaker depolarization in the previous models are different: Wilders et al. (124) and Dokos et al. (24) assumed small contribution of ICa,T according to the experimental report of Hagiwara et al. (37), whereas Demir et al. (17) assumed much larger ICa,T based on the data from Doerr et al. (23). With the use of the ICa,T expressions of Demir et al. (17), our standard model could reproduce the relatively large effects of blocking ICa,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 ICa,T expressions of Wilders et al. (124), with CL increasing only by 2.1% on eliminating ICa,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 dT and inactivation gating variable fT. The conductance property of ICa,T was formulated using a linear voltage relation, with the reversal potential ECa,T fixed at +45 mV and the maximum conductance gCa,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 IKr, exhibiting strong inward rectification, and 2) the slowly activating component IKs, exhibiting only weak rectification (2, 34, 64, 66, 100, 101, 123). Whereas IK in the rabbit SA node appears to predominantly reflect only IKr-type behavior (2, 108), the IKs-type component has also been identified in isolated rabbit SA node cells (34, 48, 64). For the standard pacemaker model, therefore, we assumed both IKr- and IKs-type components of IK.
The voltage dependences of IKr activation and inactivation (steady-state probabilities and time constants for the activation gating variable pa and inactivation gating variable pi) are shown in Fig. 3, top, along with those for the previous models. To describe the gating kinetics of IKr, 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 IKr in single rabbit SA node cells and mathematically described the activation kinetics with two gating variables, paF and paS. According to their report, therefore, we described the general activation variable pa as a weighted sum of the fast (paF) and slow (paS) activation variables and used their original expressions for IKr activation kinetics (pa,
,
paF, and
paS). A modified version of the
formulation of Ono and Ito (90) for
IKr activation was also utilized by Zhang et al.
(130). We also employed the expression of Ono and Ito
(90) for steady-state inactivation
(pi,
). No detailed experimental data are
available on the time constant of the voltage-dependent IKr inactivation
(
pi); thus we adopted the
expression of Shibasaki (108) for
pi.
|
60 and
55 mV
(approximately equal to
58 mV) to be achieved. The standard
gKr value of our model is smaller than that
determined by Ono and Ito (90); this difference in
gKr may reflect the regional difference in
IKr 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 IKr single channel conductance by [K+]o,
as reported by Shibasaki (108). Dokos et al.
(24) pointed out that the reversal potential of
IKr (EKr) is positive to
EK by 10~19 mV, thus assumed that
IKr channels are slightly permeable to
Na+. According to recently published reports (e.g., Refs.
90 and 119), however, EKr is nearly
equal to EK; therefore, the
IKr channel was assumed to be highly selective
to K+. The model-generated IKr
during the voltage-clamp pulses simulating the experiment of Ono and
Ito (90) are shown in Fig. 3, bottom.
The kinetics of IKs were described by the
formulation of Zhang et al. (130). The steady-state
activation curve for IKs 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 IKs in rabbit SA node cells. Thus
the time constant of IKs 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 (Ito) and sustained (Isus) 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 Ito and Isus.
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 Q10 of 2.18. However, the use of
Q10 = 2.18 yielded a pronounced phase 1 notch in APs.
Thus we slightly accelerated the gating kinetics by using a
Q10 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
Ito and Isus were
determined by exploring the change in peak overshoot potential (POP)
and variation of [K+]i . The experimentally
measured densities of Ito and
Isus were significantly correlated with
Cm, i.e., cell size, and are larger in cells
with a higher Cm (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
Cm, because the current densities in the central
SA node region would be relatively small (see Refs. 64 and
130). Whereas Zhang et al. (130) set the Isus 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 Ih 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
Ih; however, the contributions of
Ih 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
Ih 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
Q10 = 2.3. According to van Ginneken and Giles
(118), the maximum conductance and reversal potential of
Ih were assumed to be 0.375 nS/pF and
24 mV,
respectively. These values were achieved by setting
gh,K = 7.4 nS and
gh,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 (Ist), is carried by Na+ under physiological conditions, blocked by both organic and inorganic Ca2+ channel blockers as well as by external Ca2+ and Mg2+, and enhanced by isoprenaline or a Ca2+ agonist, BAY K 8644. These biophysical and pharmacological characteristics are compatible with those of the monovalent cation conductance of the L-type Ca2+ channel; thus Guo et al. (32) concluded that Ist is generated by a novel subtype of the L-type Ca2+ channel. Ist has also been recorded in SA node cells of other animal species, such as guinea pigs and rats (31, 77, 78, 109). Because Ist was observed only in spontaneously beating SA node cells but was absent in quiescent cells, Ist may play an essential role in the generation of intrinsic cardiac automaticity (32, 78). Therefore, we incorporated Ist in our standard SA node model for a primary pacemaker cell, whereas previous SA node models did not include Ist.
The voltage dependences of Ist activation and inactivation (steady-state probabilities and time constants for the activation gating variable qa and inactivation gating variable qi) are shown in Fig. 4, top. Because there is no detailed report on the gating kinetics (time constant) of Ist in the rabbit SA node, we modified the formulation of Shinagawa et al. (109) for the rat SA node Ist. 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 Ist, as
reported in Guo et al. (32), appear to be slower than those of the rat Ist reported by Shinagawa et
al. (109). Thus the time constants of inactivation and
recovery of the rabbit Ist were assumed to be
6.65 times larger than those of the rat Ist. The
maximum conductance and reversal potential of
Ist 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 Ist could reproduce the kinetic properties of Ist
reported by Guo et al. (32).
|
Na+ channel current. Although most of the previous SA node models incorporated the fast Na+ channel current (INa), the contribution of INa to model cell dynamics is relatively small: in the models of Wilders et al. (124) and Dokos et al. (24), eliminating INa 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 INa into our model to assess the possible contribution of INa to the pacemaker activity of our model cell. The kinetics of INa 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 INa was described by a weighted sum of two gating variables, hF and hS, with the fraction of slow inactivation (FNa) 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 Q10 of 1.7 was used to correct the experimental data for 37°C. Incorporating INa 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 INa to pacemaking is relatively small in our model cell as well.
INa 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 INa in the SA node decrease with development, with INa in adult rabbit SA node cells being very small or negligible. Thus INa 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 INa. In our simulations for primary pacemaker (central SA node) cells, therefore, INa was assumed to be zero (negligible or completely inactivated).Na+-dependent background current
and muscarinic K+ channel current.
Our model includes two background current components: 1)
Na+-dependent background current
(Ib,Na), reported by Hagiwara et al.
(38), in which Ib,Na was measured
as 0.73 ± 0.21 pA/pF at
50 mV; and 2) muscarinic
K+ channel current (IK,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
IK,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
Ib,Na with an ohmic I-V
relationship and defined IK,ACh as a background
K+ current component exhibiting inward rectification. The
Ib,Na conductance gb,Na
was set to 5.4 pS/pF, yielding a current density of 0.65 pA/pF at
50
mV. The gb,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 IK,ACh, we chose
the formula of Dokos et al. (24) adopted from Egan and
Noble (26). The standard value of
IK,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 INaK were not based on experimental data from SA node cells. The novel INaK 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): Km values for Na+ and K+ binding were set to 14.0 and 1.4 mM, respectively. The voltage dependence of INaK 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 INaK is linked to the influx of Na+ through Ib,Na (and Ist) and efflux of K+ through K+ channels. According to the data of Sakai et al. (99), we set the maximum attainable current (INaKmax) 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+/Ca2+
exchange current.
Most existing models of SA node activity have relied on the
hypothetical INaCa formulation utilized by
DiFrancesco and Noble (22), describing the
Na+/Ca2+ 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 INaCa at large
[Ca2+]i and for negative potential-low
[Ca2+]i conditions. To describe
INaCa kinetics, we adopted the formulation of
Dokos et al. (24), a modified version of the
Matsuoka-Hilgemann "E4" model. The scaling parameter
kNaCa 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
INaCa density of 1 pA/pF at
40 mV and
[Ca2+]i = 0.5 µM. A
kNaCa value of 125 pA/pF yielded an
INaCa of 0.83 pA/pF at
40 mV and
[Ca2+]i = 0.5 µM, giving a diastolic
free Ca2+ level in the presumed physiological range of
<0.3 µM.
SR functions. The SR was modeled as consisting of two compartments: the Ca2+ 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 Ca2+ 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 Ca2+ buffering by calsequestrin. Their formulation is, however, overly complex, given the complete lack of associated data in SA node cells; in their model, Ca2+ 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 Ca2+ 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 Ca2+ uptake and release in SA node cells, the parameter values were determined from previous modeling studies. The formula for the Ca2+ release mechanism was adopted from Dokos et al. (24), with the kinetics of Ca2+ release to subspace expressed as a function of [Ca2+]sub (not [Ca2+]i). The values of Krel used in previous models ranged from 0.8 to 2.0 µM; we chose a medium value of 1.2 µM. The Prel value was determined so as to make the peak [Ca2+]sub and [Ca2+]i nearly maximum with a Krel of 1.2 µM. The formulation for Ca2+ 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).
Ca2+ buffering. Our model includes the dynamics of three Ca2+ buffers: calmodulin, troponin, and calsequestrin. The amounts of calmodulin and troponin available for Ca2+ binding, and the rate constants for Ca2+ 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 Ca2+ binding to calmodulin and troponin were scaled with the use of a Q10 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 Ca2+ binding to calsequestrin were based on the study of Cannell and Allen (12), adjusted to 37°C via a Q10 of 1.6 (see Ref. 68).
We also examined the effects of exogenous Ca2+ buffers, BAPTA and EGTA, on pacemaker activity. The rate constants for Ca2+ 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 [Ca2+]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 Ca2+. 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