Hyperpolarization-Activated Current, I f, in Mathematical Models of Rabbit Sinoatrial Node Pacemaker Cells

A typical feature of sinoatrial (SA) node pacemaker cells is the presence of an ionic current that activates upon hyperpolarization. The role of this hyperpolarization-activated current, I f, which is also known as the “funny current” or “pacemaker current,” in the spontaneous pacemaker activity of SA nodal cells remains a matter of intense debate. Whereas some conclude that I f plays a fundamental role in the generation of pacemaker activity and its rate control, others conclude that the role of I f is limited to a modest contribution to rate control. The ongoing debate is often accompanied with arguments from computer simulations, either to support one's personal view or to invalidate that of the antagonist. In the present paper, we review the various mathematical descriptions of I f that have been used in computer simulations and compare their strikingly different characteristics with our experimental data. We identify caveats and propose a novel model for I f based on our experimental data.


A Funny Current
The sinoatrial node (SA node) is the normal pacemaker of the mammalian heart and generates the electrical impulse for the regular, rhythmic contraction of the heart. The intrinsic pacemaker activity, or spontaneous electrical activity, of an SA nodal pacemaker cell is based on the spontaneous diastolic depolarization that depolarizes the cell towards the action potential threshold. During this diastolic depolarization phase, there is a tiny net inward current across the cell membrane of no more than a few picoamps in amplitude. Animal studies, almost exclusively carried out on cells isolated from rabbit heart, have learned that this net inward current is the result of a complex interaction of multiple inward and outward ion currents, including a hyperpolarization-activated current of mixed ionic nature, known as "funny current, " , as reviewed elsewhere [1][2][3][4][5][6][7][8][9]. In line with its activation at hyperpolarized membrane potentials [10], thus generating an inward current during diastole, its enhancement by direct binding of cyclic AMP (cAMP) [11], and its principal presence in primary [12] and secondary pacemakers [13,14], is traditionally also named "pacemaker current. " Of note, the pore-forming subunits of the channel are encoded by members of the HCN gene family, with members HCN1-4 (see [8] and primary references cited therein). HCN channels are not only expressed in the heart but also in the brain. In neuroscience, the HCN current is usually designated ℎ .
In recent years, has regained interest in several fields. First, has become a target for pharmacological reduction of heart rate, which may be beneficial for heart failure patients. This reduction is achieved through the specific blocker ivabradine, which has become available for clinical use, and represents a new approach in selective heart rate reduction [15]. Second, mutations in the HCN4 gene, encoding the major HCN isoform of the human SA node [16], have been associated with hereditary SA nodal dysfunction in several families [17][18][19][20][21][22][23]. Third, HCN channels are used in engineering a biological pacemaker, as summarized in numerous review papers, for example, [24][25][26][27][28][29]. In all of these fields, an appropriate quantitative model of the electrical activity of is a desirable tool.
It should be noted that depends on intracellular calcium levels, and conversely, through the mutual interactions between amplitude, spontaneous firing frequency and intracellular Ca 2+ cycling [30,31]. Thus, a blockade 2 BioMed Research International of by Cs + [30] or ivabradine [31] does not only affect firing frequency, but also intracellular Ca 2+ cycling. In the interactions between the "membrane clock" (composed of voltage-dependent sarcolemmal currents and also designated "M clock, " "voltage clock, " or "ion channel clock") and the "calcium clock" or "Ca 2+ clock" (composed of tightly coupled sarcoplasmic reticulum Ca 2+ cycling molecules together with the electrogenic sodium-calcium exchanger), cAMP plays an important role [9]. It is, therefore, important to make measurements on with the amphotericin-perforated patch configuration of the patch-clamp technique to avoid the dialyzing effects of the common whole-cell configuration. Instead of rupturing the membrane as in the whole-cell patch clamp configuration, amphotericin B is used to gain electrical access to the cell [32]. This substance makes minute holes in the membrane that allow the passage of small monovalent ions, thus leaving the cytosolic composition intact. In particular, the intracellular Ca 2+ and cAMP levels are preserved.
Differences in recording conditions like the aforementioned patch clamp configurations may readily explain part of the variability in experimental data on between laboratories. A further source of variability is observed in measurements on HCN channels in heterologous expression systems like HEK-293 cells [8]. Here, part of the variability can be attributed to differences in the expression level of members of the MinK family of single transmembrane spanning proteins, which are encoded by the KCNE gene family and can act as -subunits for the HCN family of pore-forming -subunits [33][34][35]. In particular, the MinK-related peptide 1 (MiRP1, encoded by KCNE2), with high mRNA levels in the rabbit SA node [33], may associate with HCN1, HCN2, and HCN4 and modulate the HCN channel expression and kinetics.
Originally, was termed "funny current" because of its atypical characteristics, including its slow activation upon hyperpolarization rather than depolarization [10], its direct activation by cAMP [11], and its highly selective permeability to both Na + and K + ions [36] with a Na : K permeability ratio of approximately 1 : 4 under physiological conditions [37]. As a result of its mixed ionic nature, exhibits a reversal potential of ≈ −30 mV if corrected for the liquid junction potential [38][39][40]. Thus, is an inward current carried by sodium ions at diastolic membrane potentials. However, this Na + current is critically dependent on the presence of extracellular K + ions. It increases with increasing extracellular K + concentration, as does the Na : K ratio, which saturates near the physiological extracellular K + concentration of ≈5 mM [41]. More recently, it was shown that channels are also permeable, albeit to it small extent, to Ca 2+ ions [42,43], with the Ca 2+ flux contributing to ≈0.5% of the current produced by [42]. However, the functional relevance of this permeability to Ca 2+ remains unclear [44].
Another characteristic feature of HCN channels, and thus possibly of , is their ability to undergo a "mode shift" in their voltage gating. This mode shift or "voltage hysteresis" has been studied for HCN1, HCN2, and HCN4 channels that were heterologously expressed in Xenopus oocytes or mammalian COS-7 or HEK-293 cells [45][46][47][48][49]. The voltage hysteresis is clearly present in HCN1 channels [45,46,48] but less pronounced in HNC2 channels [46,47]. Whereas both Azene et al. [46] and Elinder et al. [47] concluded that voltage hysteresis is absent or almost absent in HNC4 channels, Xiao et al. [49] reported a clear hysteresis. Given that in most species, including rabbit, HCN4 is by far the most abundant HCN isoform in the SA node [8], it remains to be elucidated whether voltage hysteresis of plays a functional role in cardiac pacemaker activity. However, voltage hysteresis may prove important in fine-tuning the firing frequency of geneand cell-based biological pacemakers, in particular if these make use of HCN1 or HCN2 [46].
Despite the numerous experimental studies, the contribution of to SA nodal pacemaker activity has been and still is a matter of, often intense, debate, particularly in relation to the calcium clock [50][51][52][53][54][55][56][57][58][59][60][61][62][63][64]. A complicating factor in this ongoing debate is created by the slow activation kinetics and negative activation profile of relative to the time scale and voltage range of SA nodal diastolic depolarization, which makes only a small number of channels activated during diastolic depolarization. In a total of 12 mathematical models of the pacemaker activity of rabbit SA nodal cells published between 1982 and 2003 quantitatively widely different mechanisms underlie spontaneous diastolic depolarization, as demonstrated by the 0.9-30% increase in cycle length upon block of [65]. It is, therefore, not surprising that computer simulations, albeit with "updated" models, have not only been used to support a limited role for [66] but also to underscore that is "the major inward diastolic ionic current" [67]. Figure 1 shows experimental data on in isolated rabbit SA nodal myocytes obtained at physiological temperature using the amphotericin-perforated patch-clamp technique under voltage clamp conditions. To minimize contamination by K + and Ca 2+ currents, data were recorded in the presence of 1 mM BaCl 2 , 5 M E4031, 5 M chromanol 293b, 0.5 mM 4-aminopyridine, 5 M nifedipine, and 40 M NiCl 2 in the bath solution. The concentrations of K + and Na + in the bath solution were 5.4 and 140 mM, respectively, whereas those in the pipette solution were 145 and 10 mM, respectively. Of note, the data of Figure 1 are consistent with the data that we collected in previous studies under similar conditions [40,68].

Experimental Data
As an example of a voltage clamp trace, Figure 1(a) shows the membrane current in response to a voltage clamp step to −130 mV from a holding potential of −40 mV. In this example, the recorded current was not corrected for the capacitive transient that accompanies a voltage clamp step (Figure 1(a), arrows). At the end of the 2-s step to −130 mV, is fully activated. From the "tail current" that is observed upon stepping back from −130 mV to potentials ranging between −110 and 0 mV (−40 mV in case of Figure 1(a)), one can derive the fully activated current-voltage relationship of Figure 1(b), in which the recorded current is normalized to membrane capacitance and thus expressed in pA/pF. From the linear fit to the data of Figure 1(b), assuming an ohmic current-voltage relationship, one obtains a fully activated conductance of 0.224 nS/pF and an reversal potential of −34.8 mV. Together with the sodium and potassium concentrations in the pipette and the bath, this reversal potential yields a ,Na : ,K ratio of 0.491 and a Na : K permeability ratio, according to the Goldman equation, of 1 : 4.3, in line with aforementioned permeability ratio of 1 : 4 under physiological conditions [37].
If the tail current associated with test potentials between −120 and −30 mV is normalized to the tail current associated with the step to −130 mV, at which is fully activated, one determines the voltage dependence of activation, that is, the steady-state activation of at a series of membrane potentials 4 BioMed Research International ( Figure 1(c)). The red curve in Figure 1(c) results from the Boltzmann equation: in which ∞ is the degree of steady-state activation, is the membrane potential, 0.5 is the half-activation voltage, and is the slope factor. The fitting procedure yielded 0.5 and values of −73.0 and 9.0 mV, respectively.
Ignoring the variable initial delay in (de)activation, the current traces closely resemble a monoexponential process, both during activation and deactivation (cf. Figure 1(a)), which at a given membrane potential can be well fitted using a single time constant . The resulting time constant data are shown in Figure 1(d). The red bell-shaped curve in Figure 1 in which is the time constant of (de)activation in and is the membrane potential in mV. It is important to note that the fitted curve levels off at 50 ms for membrane potentials > −10 mV, in accordance with the experimental observation that deactivation is fast but certainly not instantaneous at depolarized potentials [39,[69][70][71]. It should be noted that widely different values have been reported for the rate of deactivation near 0 mV, approximately ranging from 10-50 s −1 [39,[69][70][71], which translates into a time constant of deactivation of 20-100 ms. Thus, our value of 50 ms is in line with these experimental data but somewhat uncertain.

State Diagrams of the Channel
Several kinetic schemes have been proposed to describe the behavior of the channel, varying from a simple two-state scheme as diagrammed in Figure 2(a) to a complex scheme with as many as five open and three closed configurations [72]. In mathematical models of in rabbit SA nodal pacemaker cells-either as a model of per se or as part of a model of the pacemaker activity of SA nodal myocytestwo-, three-, and five-state kinetic schemes have been used in relation to .

Two-State Model.
In the two-state model (Figure 2(a)), the channel flips between its open (conducting) state O and its closed (nonconducting) state C, controlled by a Hodgkin and Huxley [73] type gate with voltage-dependent rate constants and . Accordingly, is given by in which is the fully activated conductance, is the reversal potential, and the gating variable , with 0 ≤ ≤ 1, obeys the first-order differential equation: or, equivalently, with ∞ = ( + ) , In case of a voltage clamp step at time zero, the analytical solution to (4) becomes The gating variable thus attains a new steady-state value ∞ in a monoexponential fashion with time constant . The advantage of the two-state model is that it allows a direct translation of experimental data on , which are commonly acquired under voltage clamp conditions and presented in terms of a Boltzmann equation like that of Figure 1(c) and time constants of (de)activation, into a mathematical description through The two-state model has been used by DiFrancesco and Noble [70], Dokos et al. [74], and Zhang et al. [75].

Three-State Model.
Experimentally, an initial delay or sigmoidal onset may be observed in both activation and deactivation of under voltage clamp conditions. Following the approach by Hodgkin and Huxley [73] to describe the sigmoid course of activation of the potassium current in their nerve fibers, van Ginneken and Giles [39] introduced a second gate in their mathematical description of to account for the observed delay "semiquantitatively. " Accordingly, is now given by The experimentally determined Boltzmann curve then corresponds with 2 ∞ rather than ∞ , while determination of the rate constants and requires detailed analysis of the voltage clamp traces, as carried out by van Ginneken and Giles [39]. Unfortunately, they erroneously doubled their experimentally observed decay rate to obtain the decay rate of , whereas they should have halved it, as set out in detail by Dokos et al. [74].
With two independent gates, there are four different states of the channel (Figure 2(b)). However, because the two gates are kinetically identical, the four-state scheme of Figure 2   The three-state model of Figure 2(c), with two independent, kinetically identical gates, has been employed by Demir et al. [76], Kurata et al. [77], Maltsev and Lakatta [66], and Severi et al. [67].
Unlike the two-state model of Figure 2(a), the three-state model accounts for the sigmoidal onset of activation that may be observed in voltage clamp traces. It should, however, be realized that this delay in activation shows a remarkable variability and may even be absent [38,39]. The threestate model accounts for the initial delay in activation, but not for any delay in deactivation, but such delay in deactivation is not observed in case of moderate and short hyperpolarizations as take place during normal SA nodal pacemaker activity [39].

Five-State Model.
In many mathematical models of in rabbit SA nodal pacemaker cells, a two-or three-state kinetic scheme is used [39,66,67,70,[74][75][76][77]. In the mathematical model of an SA nodal pacemaker cell by Sarai et al. [78], however, a five-state scheme is used with two closed states (C 1 and C 2 ) and three open states (O 1 , O 2 , and O 3 ), as diagrammed in Figure 2(d). This five-state scheme was introduced by Maruoka et al. [71] to describe under voltage clamp conditions and included in the "Kyoto model, " which provides a common set of equations for ventricular myocytes and SA nodal pacemaker cells, by Matsuoka et al. [79] to account for the "sigmoidal onset of activation on hyperpolarization, and delayed removal of activation on depolarization beyond the reversal potential".

Mathematical Models of
Apart from differences in the kinetic schemes, the aforementioned models of in rabbit SA nodal pacemaker cells differ in the associated rate constants as well as the conductance and reversal potential of , which together determine the course of during an action potential, that is, under current clamp conditions. In this section, we present and discuss, in chronological order, the various models of that we mentioned in the previous sections.  Table 1 and then move on to Section 5.
In several cases, the fully activated conductance and/or reversal potential of were not explicitly stated in the model description and had to be determined from the relative sodium and potassium conductance of , in combination with the intracellular and extracellular sodium and potassium concentrations of the model cell through where , K and , Na are the sodium and potassium conductance of , respectively, and K and Na are the Nernst potentials for sodium and potassium, respectively. Conversely, (11) can be used to estimate the , Na : , K ratio if and the ion concentrations are known.

Model of DiFrancesco and Noble.
The aim of DiFrancesco and Noble [70] was to provide a simple description of "relevant to the reconstruction of the diastolic phase of the spontaneous action potential. " They noted that, for reconstruction purposes, it would be sufficient to describe the kinetics of with a simple first-order Hodgkin and Huxley [73] type model and that, in this context, a more complex model was not justified. Thus, they adopted the two-state model of Figure 2  Membrane potential (mV) as equations for ∞ and (in s), respectively, from which and can be derived through (8) and (9). The black curves in Figures 3(a) and 3(b), are constructed from (12) and (13), respectively. Figure 3(a) shows the steady-state value of , that is, ∞ , for each of the two-state models discussed in Section 3, whereas Figure 3  for each of the models discussed in Section 3. In addition, the curve that we determined experimentally (Figure 1(c)) is shown as a dark gray dashed curve. In case of the simple two-state model of DiFrancesco and Noble [70], the black curve of Figure 4 Figure 4(a) illustrates that the steady-state activation curve of DiFrancesco and Noble [70] is similar in shape to the curve that we presented in Section 2-Boltzmann curves with slope factors of 9.2 and 9.0, respectively-but that the half-activation voltage is more positive (−64 versus −73 mV, as listed in Table 1).
In their description of , DiFrancesco and Noble [70] used a reversal potential of −10.3 mV and a fully activated conductance of 9909 pS, that is, 0.3303 nS/pF when normalized to their membrane capacitance of 30 pF. Thus, (3) reads with expressed in mV and in pA/pF. This yields the fully activated current that appears as a black line in Figure 4(b) and, on expanded current and membrane potential scales, in Figure 4(d). As can be appreciated from Figures 4(b) and 4(d), the fully activated conductance of 0.3303 nS/pF is ≈50% larger and the reversal potential of −10.3 mV is ≈25 mV more positive than in our experimental data.

Model of van Ginneken and Giles. van Ginneken and
Giles [39] carried out a comprehensive study on in isolated rabbit SA nodal pacemaker cells. They analyzed their : fully activated conductance; : reversal potential; ,Na : ,K : ratio of sodium and potassium conductance; 0.5 : half-activation voltage; −10 and +20 : time constant of deactivation at −10 and +20 mV, respectively. Experimental data are mean ± SEM. Data of DiFrancesco and Noble [70] and van Ginneken and Giles [39] are not corrected for liquid junction potential. experimental data in terms of the three-state kinetic scheme of Figure 2 for the associated rate constants, both expressed in ms −1 , which can be turned into ∞ and for display in Figures 3(c) and 3(d) (green curves), through (6). The steady-state activation curve of can be obtained by squaring ∞ and is shown in Figures 4(a) and 4(c). In the physiological membrane potential range, it is highly similar to our Boltzmann fit of Figure 1(c), which appears in Figures 4(a) and 4(c) as a dark gray dashed curve.
In their experiments, van Ginneken and Giles [39] observed mean values of −24 mV and 12.0 nS for the reversal potential and fully activated conductance of , respectively. When normalized to their mean membrane capacitance of 55 pF, the fully activated conductance becomes 0.2182 nS/pF, which is remarkably similar to the value of 0.224 nS/pF that we determined from our experiments (Figure 1). With these values, (10) reads with expressed in mV and in pA/pF. This yields the fully activated current that appears as a green line in Figures 4(b)  and 4(d). The reversal potential of −24 mV is more positive than we observed, but this may reflect differences in the bath and pipette solutions and the extent to which the data were corrected for the liquid junction potential.
Unfortunately, as already mentioned in Section 3.2, van Ginneken and Giles [39] erroneously doubled their experimentally observed decay rate in their analysis instead of halving it. This affects and , but not ∞ . Furthermore, it should be noted that all experiments were carried out at a fixed temperature between 30 and 33 ∘ C (±1 ∘ C), which may have led to an underestimation of the rates of activation and deactivation.

Model of Demir et al.
In their mathematical model of an SA nodal pacemaker cell, Demir et al. [76] based their equations for on the data by van Ginneken and Giles [39]. They also used the three-state model of Figure 2(c), but reanalyzed the raw data of van Ginneken and Giles [39], validating their equations in a comparison of modelgenerated voltage clamp traces with those reported by van Ginneken and Giles [39]. Furthermore, they converted the thus obtained time constant values to a temperature of 37 ∘ C, assuming a Q 10 of 2.3. This reanalysis led to  [39]. This also holds for the steady-state activation curve of obtained by squaring ∞ , which is shown in Figures 4(a) and 4(c).
In their model, Demir et al. [76] selected an reversal potential of −30 mV. For the fully activated conductance, they chose of value of 19.63 nS, which turns into a value of 0.3569 nS/pF when normalized to the model's membrane capacitance of 55 pF, which Demir et al. [76] based on the mean value reported by van Ginneken and Giles [39]. With these values, (10) reads with expressed in mV and in pA/pF. The associated fully activated current is shown in Figures 4(b) and 4(d).

Model of Dokos et al.
In 1996, two years after Demir et al. [76], Dokos et al. [74] also published a mathematical model of an SA nodal pacemaker cell. However, they used the simpler first-order Hodgkin and Huxley type model of Figure 2 for the rate constants and , both expressed in s −1 , which have been turned into ∞ and for display in Figures 3(a) and 3(b) (blue curves), using (6). The steady-state activation curve is shown in Figures 4(a) and 4(c) and differs not only from that of van Ginneken and Giles [39] but also from that of Demir et al. [76].
Dokos et al. [74] selected values of the sodium and potassium conductance of "to produce a reversal potential of −25 mV and a cycle length prolongation of ≈27% in the free-running model when was abolished. " The resulting reversal potential is −24.97 mV, and the resulting fully activated conductance is 5.102 nS, which turns into a value of 0.1595 nS/pF when normalized to the model's membrane capacitance of 32 pF, which Dokos et al. [74] adopted from the early model by Wilders et al. [80]. With these values, (3) reads = × 0.1595 × ( + 24.97) , with expressed in mV and in pA/pF. The associated fully activated current is shown in Figures 4(b) and 4(d).

Model of Zhang et al.
Zhang et al. [75] also used a singlegate description of , as in Figure 2(a), in the mathematical model of an SA nodal pacemaker cell that they published in 2000. They created two versions, one for a small cell, with a membrane capacitance of 20 pF, presumably originating from the center of the SA node, and one for a large cell, with a membrane capacitance of 65 pF, presumably originating from the periphery of the SA node. They introduced rate constants and , both in s −1 , given by which they validated in a comparison with the experimental data for ∞ and of van Ginneken and Giles [39] and Liu et al. [81]. The associated ∞ and are shown in Figures 3(a) and 3(b) (red curves) and the steady-state activation curve in Figures 4(a) and 4(c). These hold for both versions of the model. Zhang et al. [75] assumed identical sodium and potassium conductance values for . This results in a reversal potential of −5.25 mV, which differs from our experimentally observed value (Figure 1(b)) by 30 mV. The fully activated conductance of was validated against the current-voltage relationships reported by Honjo et al. [82], who found that the fully activated current density of (in pA/pF) of SA nodal cells increases with the membrane capacitance of the cells. Zhang et al. [75] for the central cell model and for the peripheral cell model, both with expressed in mV and in pA/pF. The associated fully activated current is shown as a red solid line for the central cell model and a dashed line for the peripheral cell model in Figures 4(b) and 4(d). Kurata et al. In 2002, Kurata et al. [77] published "an improved mathematical model for a single pacemaker cell of the rabbit SA node. " This primary cell model has a membrane capacitance of 32 pF, in line with the earlier models by Wilders et al. [80] and Dokos et al. [74]. The formulation for was adopted from the model of Wilders et al. [80], who had in turn used the equations and parameters of van Ginneken and Giles [39] and arrived on potassium and sodium conductances of 7.4 and 4.6 nS, respectively, for , based on the observed reversal potential of −24 mV. Thus, ∞ is given by (15), which explains that the ∞ curve in Figure 3(c) and the steady-state activation curve in Figures  4(a) and 4(c), of Kurata et al. [77] and van Ginneken and Giles [39], coincide.

Model of
The equation for the time constant was derived from the rate constants of (16) but corrected for 37 ∘ C by the use of a Q 10 factor of 2.3 following Demir et al. [76], thus arriving at in which is in ms and 0.71665 is the correction factor for the temperature of 30-33 ∘ C in the experiments of van Ginneken and Giles [39]. Thus, the Kurata et al. [77] curve in Figure 3(d) is similar in shape but smaller in magnitude than the van Ginneken and Giles [39] curve. For the fully activated conductance of , Kurata et al. [77] used the value of 12 nS of van Ginneken and Giles [39], which, with their membrane capacitance of 32 pF, results in a normalized value of 0.375 nS/pF in the equation for : in which is again expressed in pA/pF and in mV. The reversal potential of −26.02 mV deviates from the value of −24 mV of the Wilders et al. [80] model, because there are some minor differences in sodium and potassium concentrations between the models. The fully activated current is shown in Figures 4(b) and 4(d).
where and are both in ms −1 and in mV. The remaining three transitions (see Figure 2(d)) are all controlled by rate constants and given by  Figure 3(e). This yields the curve shown in Figures 4(a) and 4(c), which is strikingly steep when compared to each of the other curves.
Unlike the other models, the fully activated current ,full is not ohmic and thus linearly dependent on the membrane potential but determined by ,full = 1.821 × CF Na + 7.7286 × CF K , (30) in which 1.821 and 7.7286 are the permeabilities for Na + and K + in pA/mM, and CF Na and CF K are given by the "constantfield equations"

Model of Maltsev and Lakatta. In 2009, Maltsev and
Lakatta [66] published a mathematical model of an SA nodal pacemaker cell that is based on the model by Kurata et al. [77] but incorporates a submembrane "calcium clock" that interacts with the "membrane clock" in producing the pacemaker activity of the model cell [6,9]. Maltsev and Lakatta [66] adopted the kinetics of Kurata et al. [77]. This explains why the Maltsev and Lakatta [66] and Kurata et al. [77]

Model of Severi et al.
The most recent mathematical model of a rabbit SA nodal pacemaker cell is that of Severi et al. [67], which they presented as "an updated computational model of rabbit sinoatrial action potential to investigate the mechanisms of heart rate modulation. " In this model, the kinetic and conductive properties of are largely based on the work of DiFrancesco and Noble from the 1980s. The kinetics are adopted from the early SA nodal cell model by Noble et al. [83], who used a Hodgkin and Huxley type model with two identical gates, as in Figures 2(b) and 2(c). However, Severi et al. [67] shifted the associated ∞ and curves to more depolarized potentials by ≈11 mV, based on experimental data from Altomare et al. [84] and Barbuti et al. [85], producing where is in mV and in s. The ∞ and curves are shown in Figures 3(c) and 3(d) and the steady-state activation curve in Figures 4(a) and 4(c).
According to the original description by DiFrancesco and Noble [86] and in line with Noble et al. [83], Severi et al. [67] assumed identical conductance values for the sodium and potassium components of , with a total conductance of 6.429 nS. Because, as in the models by Wilders et al. [80], Dokos et al. [74], and Maltsev and Lakatta [66], the model cell has a membrane capacitance of 32 pF, the normalized conductance amounts to 0.2009 nS/pF (Table 1). With the model concentrations for sodium and potassium ions, the reversal potential amounts to −4.39 mV, which is 7.9 mV more positive than in the model by Noble et al. [83], due to differences in ion concentrations between the two models. Thus, is given by in which is expressed in pA/pF and in mV and which yields the fully activated current shown in Figures 4(b) and 4(d).

Novel Model.
Apart from the ten models detailed in Sections 4.1-4.9, Figures 3 and 4 show curves labeled "present paper. " These curves represent a novel model for based on the experimental data that we presented in Section 2. We use the kinetic scheme of Figure 2(a) and describe by where is in pA/pF, is in mV, and in s. As can be appreciated from Figure 4(c), the steady-state activation curve of our model closely matches with that of van Ginneken and Giles [39], Kurata et al. [77], and Maltsev and Lakatta [66] in the physiological membrane potential range, whereas there are significant discrepancies with those of other models, in particular the models by DiFrancesco and Noble [70], Sarai et al. [78], and Severi et al. [67]. The conductance of 0.224 nS/pF, on the other hand, closely matches the values of 0.2182, 0.2123, and 0.2009 nS/pF of the models by van Ginneken and Giles [39], Zhang et al. [75] (peripheral cell), and Severi et al. [67], respectively, as can be appreciated from the slope of the lines in Figure 4(d). However, Figure 4(d) also illustrates that the reversal potential of the latter two models differs from that of our model by as much as 30 mV, which creates an almost twofold difference in driving force near the maximum diastolic potential of an SA nodal action potential.
Our model does not have an explicit cAMP dependence. However, autonomic modulation of can be incorporated through a shift of the steady-state activation curve along the voltage axis, by up to ≈ 10 mV in the positive direction for adrenergic modulation and up to ≈10 mV in the negative direction for cholinergic modulation. Such shift has been observed experimentally [11,87] and has been incorporated in several models, for example, [88][89][90][91][92], to reflect the autonomic modulation of through acetylcholine and (nor) adrenalin. In addition, a similar voltage shift should be applied to the time constant curve to account for the experimentally observed cAMP dependence of this curve [93]. The latter shift has been ignored in most models, but not in the recent model of Severi et al. [67].

Reconstruction of
In the previous section, we have identified models of in terms of characteristics derived from and related to voltage clamp experiments, including rate constants, steady-state activation, fully activated conductance, and reversal potential. In the present section, we show how these characteristics determine the course of during an action potential, that is, under current clamp conditions.

Steady-State Current.
Before reconstructing during an action potential in Section 5.2, we use the data of Figure 4 to compute the steady-state current at −60 mV to get an impression of the amplitude of that one might expect for each of the models. Figure 5(a) shows the steady-state activation at −60 mV in each of the 11 models. This steadystate activation ranges from 0.042 in the Demir et al. [76]  model to 0.486, that is, almost 50% activation of at −60 mV, in the Severi et al. [67] model, which constitutes an almost 12fold difference. Most models, including the one based on our experimental data, predict a value near 19% for the steadystate activation at −60 mV (Figure 5(a)).
As for the steady-state activation, there is a wide difference in the fully activated current amplitude at −60 mV, which is not only determined by the conductance but also by the reversal potential as a determinant of the driving force. With a value of 16.4 pA/pF, the model by DiFrancesco and Noble [70] has the largest amplitude, whereas the central cell model by Zhang et al. [75] shows the smallest amplitude with a value of 3.0 pA/pF, a 5.5-fold difference.
In combination, the steady-state activation of Figure 5(a) and the fully activated current amplitude of Figure 5(b) determine the amplitude of that can be activated at −60 mV. We multiplied the fully activated current amplitude of Figure 5 Figure 5 shows that the amount of that can be activated at −60 mV varies widely between the models, but this does not imply that this is also the case during the course of action potential. In the latter case, the rate at which activates and deactivates plays an important role. Therefore, we subjected each of the models to an "action potential clamp": we reconstructed during the experimentally recorded action potentials of Figure 6(a). These action potentials were applied as part of a sufficiently long train, and was computed according to the equations listed in Section 4. The resulting traces are shown in Figure 6(c), together with the net membrane current, net , which was computed from the time derivative ( / ) of the membrane potential trace of Figure 6(a), as shown in Figure 6(b). From the current traces shown in Figure 6(c), we computed the diastolic current amplitude at −60 mV (Figure 7(a)) as well as the maximum current amplitude during diastole (Figure 7(b)). Also, we computed the charge carried by during the 200 ms, 25 mV diastolic depolarization from the maximum diastolic potential of −63 mV to −38 mV (Figure 7(c)). Both Figure 6(c) and Figure 7(a) demonstrate that only a fraction of the steady-state current of Figure 5(c) is activated during an action potential. This fraction varies from 0.3% for the Sarai et al. model [78] to 42% for the Demir et al. model [76]. Figure 7(c) shows that, in the absence of other inward or outward membrane currents, the charge carried by during the 200 ms depolarization would be sufficient or almost sufficient to depolarize the membrane by the observed 25 mV.

Dynamics of .
For example, our model based on the experimental data of Section 2 predicts a charge carried by of 0.024 pC/pF, which is equivalent to a depolarization of 24 mV. Notably, the models of Maltsev and Lakatta [66] and Severi et al. [67] predict depolarizations of 12 and 67 mV, respectively.
With a peak inward current of only 0.027 pA/pF (Figures 6(c) and 7(b)), the smallest is obtained with the Sarai et al. model [78], although a relatively large can be activated under steady-state conditions ( Figure 5(c)). This emphasizes the important role of the rate at which activates and deactivates, but in this particular case the exceptional steepness of the steady-activation curve (Figures 4(a) and 4(c)) also plays an important role. The importance of the rate of (de)activation is perhaps better illustrated with the Zhang et al. [75] peripheral and Kurata et al. [77] models that show a similar amplitude in Figure 5(c) but clearly different peak inward currents of 0.68 pA/pF for the Zhang et al. [75] peripheral cell model and threefold less, 0.23 pA/pF, for the Kurata et al. [77] model (Figures 6(c) and 7(b)). Of note, this cannot be explained by the use of a single-gate kinetic scheme by Zhang et al. [75] versus a double-gate scheme by Kurata et al. [77], as the effects of the selection of a single or double gate on the reconstructed are minimal [58,83]. Rather, it may reflect the erroneously overestimated deactivation rate of (see Sections 3.2 and 4.2) in the Kurata et al. [77] model, which also affects the traces obtained with the related models by van Ginneken and Giles [39] and Maltsev and Lakatta [66]. The latter two traces almost coincide despite the smaller conductance in the Maltsev and Lakatta [66] model. Here, the faster kinetics of the Maltsev and Lakatta [66] compensate for this smaller conductance.
In most models, deactivates almost instantaneously at depolarized potentials, resulting in almost complete deactivation of near the overshoot of the action potential and a slowly developing during the subsequent diastolic phase. However, as set out in Section 2, deactivation is fast but certainly not instantaneous at depolarized potentials. If this is taken into account, as in our novel model, is available early in diastole and of relatively constant amplitude during diastolic depolarization. Interestingly, Zaza et al. [94] already noted that the presence of measurable inward 2 mM Cs + sensitive current almost immediately after repolarization in their action potential clamp experiments on rabbit SA nodal cells is apparently at odds with the slow kinetics of activation at diastolic potentials, that this suggests that may not deactivate completely during repetitive activity, and that this would also increase the amount of available during diastolic depolarization.

Conclusion
The various mathematical descriptions of that have been used in computer simulations show strikingly different characteristics when reconstructing the course of during an action potential. This explains-at least to some extent-that some successfully use computer simulations to support their view that plays a fundamental role in the generation of pacemaker activity and its rate control, while others provide computer simulation results in favor of their view that the role of is limited to a modest contribution to rate control. We have identified some important caveats regarding the reconstruction of the course of during an action potential through computer simulations. An obvious one is the use of appropriate activation kinetics and an appropriate value for the conductance. The half-activation voltage and fully activated conductance of most models, with values of −80 to −60 mV and ≈0.2 nS/pF, respectively, are in line with the experimental data. However, the steep steady-state activation curve of the complex five-state model of Sarai et al. [78] is clearly at odds with the experimental data. A somewhat less obvious caveat is the selection of an reversal potential that is in line with the experimental data. This reversal potential should be around −30 mV, but this is definitely not the case for the models of DiFrancesco and Noble [70], Zhang et al. [75], and Severi et al. [67]. The final and most important caveat to be taken into account is that deactivation is not almost instantaneous at depolarized potentials. In most models, including the recent models by Maltsev and Lakatta [66] and Severi et al. [67], this deactivation is much faster than can be deduced from the scarce experimental data, which likely results in an underestimation of during diastolic depolarization. Our novel model for the reconstruction of in mathematical models of SA nodal pacemaker cells is simple and straightforward but takes care of all of these caveats.