## Abstract

To elucidate the effects of hyperpolarization-activated current *I*_{f} on robustness of sinoatrial node (SAN) pacemaking in connection with intracellular Na^{+} concentration (*Na*_{i}) changes, we theoretically investigated *1*) the impacts of *I*_{f} on dynamical properties of SAN model cells during inhibition of L-type Ca^{2+} channel currents (*I*_{CaL}) or hyperpolarizing loads and *2*) *I*_{f}-dependent changes in *Na*_{i} and their effects on dynamical properties of model cells. Bifurcation analyses were performed for *Na*_{i}-variable and *Na*_{i}-fixed versions of mathematical models for rabbit SAN cells; equilibrium points (EPs), limit cycles (LCs), and their stability were determined as functions of model parameters. Increasing *I*_{f} conductance (*g*_{f}) shrank *I*_{CaL} conductance (*g*_{CaL}) regions of unstable EPs and stable LCs (rhythmic firings) in the *Na*_{i}-variable system but slightly broadened that of rhythmic firings at lower *g*_{f} in the *Na*_{i}-fixed system. In the *Na*_{i}-variable system, increased *g*_{f} yielded elevations in *Na*_{i} at EPs and during spontaneous oscillations, which caused EP stabilization and shrinkage in the parameter regions of unstable EPs and rhythmic firings. As *g*_{f} increased, parameter regions of unstable EPs and stable LCs determined for hyperpolarizing loads shrank in the *Na*_{i}-variable system but were enlarged in the *Na*_{i}-fixed system. These findings suggest that *1*) *I*_{f} does not enhance but rather attenuates robustness of rabbit SAN cells via facilitating EP stabilization and LC destabilization even in physiological *g*_{f} ranges; and *2*) the enhancing effect of *I*_{f} on robustness of pacemaker activity, which could be observed at lower g_{f} when *Na*_{i} was fixed, is actually reversed by *I*_{f}-dependent changes in *Na*_{i}.

- nonlinear dynamics
- bifurcation analysis
- mathematical model
- computer simulation
- hyperpolarizing load

the hyperpolarization-activated cation current (*I*_{f}), encoded by the hyperpolarization-activated cyclic nucleotide-gated (HCN) channel gene, is well known to play a central role in pacemaker depolarization as a pacemaker current in sinoatrial node (SAN) cells (2, 7, 28) and is a key current for engineering of biological pacemakers (44, 45). *I*_{f} contributes to prevention of excess hyperpolarization and excessively slow pacemaking against hyperpolarizing loads (15, 28, 34), autonomic regulations of spontaneous activity (5), stabilization of pacemaker frequency (15, 33, 39), and generation of diastolic depolarization in the periphery of intact SAN suffering electrotonic influences of the atrium (38). These experimental reports suggest enhancement of the robustness of SAN cell pacemaking against hyperpolarizing loads, as well as initiation of pacemaker depolarization, by *I*_{f}. Maltsev and Lakatta (33) have also concluded, using their novel SAN cell model for rabbit central SAN cells, that *I*_{f} enhances robustness of central SAN cell pacemaking against attenuation of the L-type Ca^{2+} channel current (*I*_{CaL}) or sarcoplasmic reticulum (SR) Ca^{2+} cycling, because larger *I*_{f} broadened the area of spontaneous rhythmic firings on the parametric plane of the maximum *I*_{CaL} conductance (*g*_{CaL}) and SR Ca^{2+} pumping rate (*P*_{up}). Thus *I*_{f} may enhance SAN cell robustness against intrinsic and extrinsic disturbances including hyperpolarizing loads and attenuation of the membrane and SR Ca^{2+} clocks.

In our previous theoretical study using our mathematical models for rabbit SAN cells (25), however, larger *I*_{f} did not necessarily enhance but rather attenuated the robustness of central SAN cell pacemaking against hyperpolarizing loads via parasympathetic stimulation while enhancing the robustness of peripheral SAN cell pacemaking against electrotonic loads of the atrium in combination with the Na^{+} channel current (*I*_{Na}): in the central SAN cell model, the parameter regions of unstable stationary states and rhythmic firings shrank at higher maximum *I*_{f} conductance (*g*_{f}) in physiological *g*_{f} ranges. Thus whether *I*_{f} enhances robustness of central SAN cell pacemaking is still controversial.

What are the reasons for this inconsistency? One of possible reasons is that the effect of *I*_{f} is model dependent, and another is the influence of intracellular Na^{+} concentration (*Na*_{i}) changes: the Maltsev-Lakatta model (33) is based on the Kurata et al. model (20) and thus possesses essentially the same *I*_{f} formula as the Kurata et al. model. However, *Na*_{i} was fixed at a constant value of 10 mM (handled as a constant parameter) in the Maltsev-Lakatta model (33), whereas it was handled as a state variable in our model (25). Because *Na*_{i} has been shown to increase during enhancement of *I*_{f} and significantly affect pacemaker activity (4), this discrepancy may be at least in part due to the concomitant change in *Na*_{i} during enhancement of *I*_{f}; *Na*_{i} is likely to vary during changes in *g*_{f} or other bifurcation parameters and thereby influence dynamical properties of SAN cells.

The aim of this study was to reevaluate the roles of *I*_{f} in SAN pacemaking, with particular attention to *I*_{f}-dependent changes in *Na*_{i} that may affect its impacts on robustness of SAN cells. In our previous studies (20–26), bifurcation structures (i.e., ways of changes in the number or stability of equilibrium and periodic states) of SAN and ventricular cell model systems were investigated for elucidating the dynamical mechanisms of normal and abnormal pacemaker activities, which is called “bifurcation analysis.” These theoretical works indicate that the initiation and cessation of pacemaker activity are considered as bifurcation phenomena and that the mathematical approach (bifurcation analysis) provides a convenient way of understanding how individual currents contribute to pacemaker activities (see also Ref. 15). Therefore, we again applied bifurcation analysis to the mathematical models for rabbit SAN cells. Bifurcation diagrams were constructed by calculating equilibrium points (EPs) as stationary states, periodic orbits called limit cycles (LCs), their stability, and bifurcation points as functions of model parameters such as *g*_{f}. We particularly focused on the effects of *I*_{f} on the stability of stationary state membrane potentials (*V*_{m}) and robustness of pacemaker activity against inhibition of *I*_{CaL} or hyperpolarizing loads. The robustness against hyperpolarizing loads was evaluated by determinations of stability and oscillation dynamics of the model cells during injections of hyperpolarizing bias currents (*I*_{bias}), enhancement of the background K^{+} conductance (*g*_{bK}), and applications of acetylcholine (ACh). The *I*_{f} effects were tested for both the *Na*_{i}-fixed system, including *Na*_{i} as a parameter, and the *Na*_{i}-variable system (like a real cell), in which *Na*_{i} is a state variable and thus varies during parameter changes. We further determined *Na*_{i}-dependent bifurcation structures of the *Na*_{i}-fixed system. These theoretical analyses provide significant insights into the contribution of *I*_{f} to robustness of SAN cell pacemaking as well as how *Na*_{i} influences *I*_{f} effects and other bifurcation properties of SAN cells. This study, suggesting that excess *I*_{f} attenuates robustness of pacemaker activity, would be particularly important as a theoretical background for genetic engineering of *I*_{f}-based biological pacemakers with ability of robust pacemaking and driving for gene or cell therapy of bradyarrhythmias (for review, see Refs. 14, 35, 44).

## THEORY AND METHODS

### Mathematical Formulation (Base Models for SAN Cells)

As base pacemaker cell models of the rabbit SAN, we used the Maltsev-Lakatta model (33) as well as other previously developed single-cell models (6, 10, 20, 24, 53) for central or transitional SAN cells. Results and conclusions obtained from these models were essentially the same (model independent); in this article, therefore, the results from the Maltsev-Lakatta model are mainly shown (unless otherwise stated) along with limited comparisons with other previous models. A complete list of the equations and standard parameter values for their models, as well as their dynamics, has been provided in the original articles and is also given on the CellML website (http://www.cellml.org). These models are implemented in a CellML-based open resource for public access (http://www.opencor.ws/).

The base Maltsev-Lakatta model for normal pacemaker activity includes 13 membrane current components. The membrane current system includes *I*_{CaL}, T-type Ca^{2+} channel current (*I*_{CaT}), *I*_{f}, sustained inward current (*I*_{st}), rapidly activating (*I*_{Kr}) and slowly activating (*I*_{Ks}) components of delayed-rectifier K^{+} channel currents, 4-aminopyridine (4-AP)-sensitive currents consisting of transient (*I*_{to}) and sustained (*I*_{sus}) components, background currents carried by Na^{+} (*I*_{bNa}) and Ca^{2+} (*I*_{bCa}), muscarinic K^{+} channel current (*I*_{KACh}), Na^{+}-K^{+} pump current (*I*_{NaK}), and Na^{+}/Ca^{2+} exchanger current (*I*_{NCX}). All the currents charge the membrane capacitance (*C*_{m}). The time-dependent change in the membrane potential *V*_{m} is then described by

### Bifurcation Analysis

The Maltsev-Lakatta model is a 29-order autonomous continuous-time dynamical system, which is described as a set of 29 first-order, nonlinear ordinary differential equations. Dynamical properties and bifurcation structures of the model cell system were determined by analytically and numerically handling the differential equations. For bifurcation analyses, the intracellular K^{+} concentration (*K*_{i}) was fixed at 140 mM to remove degeneracy (for details, see Refs. 19 and 22). *Na*_{i} in the model cell was fixed at 10 mM as in Maltsev and Lakatta (33) or not fixed as in Kurata et al. (20, 25); i.e., bifurcation structures were determined for both the *Na*_{i}-fixed and *Na*_{i}-variable systems for comparisons. Simulated behaviors of *V*_{m}, ionic currents, and intracellular Ca^{2+} concentrations during normal pacemaking of the *Na*_{i}-variable model cell were very similar (nearly identical) to those of the *Na*_{i}-fixed model cell, as well as to those shown in the original article (33), with the mean *Na*_{i} during normal pacemaking being 9.93 mM. Numerical computations were performed with MATLAB 7.5 (The MathWorks, Natick, MA) on Workstation HP xw9400 and Z800 (Hewlett-Packard, Tokyo, Japan).

Theory and methods for bifurcation analyses are essentially the same as described previously (24). Bifurcation parameters chosen in this study include *1*) the maximum ion channel conductances *g*_{CaL} and *g*_{f}; *2*) *Na*_{i}; *3*) amplitude of hyperpolarizing *I*_{bias}; *4*) *g*_{bK}; and *5*) ACh concentration ([ACh]). The maximum conductance *g*_{CaL} is given in picosiemens per picofarads (pS/pF) unless otherwise stated, whereas *g*_{f} is expressed as a normalized value, i.e., ratio to the control value. Detailed procedures for *1*) locating EPs and LCs, *2*) constructing one- and two-parameter bifurcation diagrams, and *3*) detecting bifurcations (determination of EP and LC stabilities) are provided in our previous article (24). Definitions of the specific terms for the nonlinear dynamics and bifurcation theory are also given in our previous article (24) as well as in textbooks (e.g., 27, 40, 46).

The critical parameter value at which a bifurcation occurs (bifurcation point) was determined by evaluating the existence and asymptotic stability of EPs and LCs at each value of a relevant parameter. The Hopf bifurcation point at which the stability of an EP reverses was located by stability analysis (46, 51). Period-doubling and Neimark-Sacker bifurcation points where an LC becomes unstable with emergence of arrhythmic dynamics were also detected (for details, see Refs. 12, 15, and 25).

To assess the robustness of pacemaker activity against hyperpolarizing loads, we determined the critical values of hyperpolarizing *I*_{bias}, *g*_{bK}, or [ACh] at which a bifurcation to quiescence or arrhythmic dynamics occurs. Here, the mathematical term “robustness” for a pacemaker system is defined as its ability to keep rhythmic firings against changes in key parameters, i.e., perturbations (interventions) that may cause bifurcations to cessation of pacemaker activity or arrhythmic dynamics. Hopf bifurcations to yield stable EPs and quiescence, as well as period-doubling or Neimark-Sacker bifurcations of LCs to destabilize LCs with spawning of arrhythmic oscillations, were detected; if the change of a bifurcation parameter required for yielding a bifurcation is greater in *system A* than in *system B*, one can say that *system A* is more robust than *system B* against the relevant perturbation (for more details, see Ref. 21).

## RESULTS

### Impacts of I_{f} on g_{CaL}-Dependent Bifurcation Structures of SAN Model Cells

We first examined the influences of changing *g*_{f} on robustness of SAN cell pacemaking by exploring bifurcation structures of the *Na*_{i}-variable and *Na*_{i}-fixed model cell systems during inhibition of *I*_{CaL} as done in Maltsev and Lakatta (33).

#### Enhancing I_{f} shrank g_{CaL} regions of unstable EPs and stable LCs in the Na_{i}-variable system.

Figure 1 shows one-parameter bifurcation diagrams for *g*_{CaL}, where stationary-state *V*_{m} at EPs, *V*_{m} extrema of LC oscillations (spontaneous firings), and bifurcation points were plotted as functions of *g*_{CaL} for the *Na*_{i}-variable and *Na*_{i}-fixed model cells with different *g*_{f} (0, 1, or 2 times the control value of 150 pS/pF). In the *Na*_{i}-variable system, the *g*_{CaL} region of stable LCs (rhythmic firings) as well as unstable EPs significantly shrank with increasing *g*_{f}. In contrast, in the system with *Na*_{i} fixed at 10 mM (or various values ranging over 8–12 mM), increasing *g*_{f} slightly enlarged the *g*_{CaL} region of spontaneous firings, as shown by Maltsev and Lakatta (33), whereas it slightly shrank the region of unstable EPs. As illustrated in Fig. 2*A*, spontaneous firings in the *Na*_{i}-variable system with larger *g*_{f}, as well as those in the *Na*_{i}-fixed system with smaller *g*_{f}, were abolished at higher *g*_{CaL}. Thus the *I*_{f} effects on stability and bifurcations of LCs were opposite in the *Na*_{i}-variable and *Na*_{i}-fixed systems. The *g*_{CaL} regions of stable EPs became broader with increasing *g*_{f} in both the *Na*_{i}-variable and *Na*_{i}-fixed systems: the systems with larger *g*_{f} did not generate spontaneous oscillations but were quiescent (at a stable EP) at higher *g*_{CaL} (Fig. 2*B*).

For more clearly demonstrating the effect of *I*_{f}, particularly that of its overexpression as achieved in engineering of *I*_{f}-based biological pacemakers (41), the *g*_{CaL}-dependent bifurcation structures of the model systems were further tested over broader ranges of *g*_{f}, as summarized by two-parameter bifurcation diagrams for *g*_{CaL} and *g*_{f} where the loci of bifurcations of EPs and LCs are plotted on the *g*_{f}-*g*_{CaL} parametric plane (Fig. 3). The parametric planes were divided into three different areas: *1*) unstable EP and rhythmic oscillation (UR), *2*) stable EP and quiescence (SQ), and *3*) unstable EP and arrhythmic oscillation (UA) in the *Na*_{i}-variable system or coexistence of stable EP and rhythmic oscillation (SR), i.e., bistable zone (BS), in the *Na*_{i}-fixed system. In the bistable zone of the *Na*_{i}-fixed system (e.g., *g*_{CaL} = 279–340 pS/pF at *g*_{f} = 1 as shown in Fig. 2), the transitions between two states (a stable resting state and a rhythmic firing state) called annihilation and single-pulse triggering occur (for more details, see Refs. 15 and 29). It should be noted that the *I*_{f} effects on *g*_{CaL}-dependent bifurcation structures were remarkably different in the *Na*_{i}-variable and *Na*_{i}-fixed systems, and that the *g*_{CaL} regions of unstable EPs and stable LCs (rhythmic firings) were not broadened but shrank at higher *g*_{f}, particularly in the *Na*_{i}-variable system. In the *Na*_{i}-variable system, increasing *g*_{f} did not broaden but monotonically shrank the *g*_{CaL} regions of unstable EPs and rhythmic firings and also broadened the region of arrhythmic dynamics by shifting LC bifurcations toward larger *g*_{CaL}, as shown in Figs. 1*A* (*right*) and 3*A*. In agreement with the report of Maltsev and Lakatta (33), the *Na*_{i}-fixed system certainly exhibited slight enlargement of the *g*_{CaL} region of rhythmic firings with increasing *g*_{f} (up to 3-fold increase). Nevertheless, *g*_{f} increases of more than threefold led to the shrinkage in the *g*_{CaL} region of rhythmic firings in the *Na*_{i}-fixed system, as well.

#### Enhanced I_{f} facilitated EP stabilization via its voltage-dependent kinetics.

The shrinkage in the parameter region of unstable EPs at higher *g*_{f} is at least in part ascribed to the stabilization of stationary-state *V*_{m} at an EP via voltage-dependent *I*_{f} kinetics. As exemplified in Fig. 4 for the *Na*_{i}-variable system, small depolarization from the stationary-state *V*_{m} caused instantaneous reduction of inward *I*_{f}, which is greater at higher *g*_{f}; resultant outward shift of the total sarcolemmal membrane current (*I*_{total}) prevents further depolarization and facilitates repolarization to the stationary-state *V*_{m}. The greater *g*_{f}-dependent shrinkage of the unstable EP regions in the *Na*_{i}-variable system is related to the *I*_{f}-induced increase in the stationary-state *Na*_{i} at EPs (Fig. 4, *bottom left*), as discussed later.

#### Enhancing I_{f} caused increases in Na_{i} at EPs and during LC oscillations.

Why do the *I*_{f} effects on the parameter-dependent bifurcation structure differ in the *Na*_{i}-variable and *Na*_{i}-fixed systems? Because *Na*_{i} has been demonstrated to increase as *I*_{f} is enhanced (4), the inconsistency may be at least in part due to the concomitant increase in *Na*_{i} during enhancement of *I*_{f}. Thus we further examined the *g*_{f}-dependent changes in *Na*_{i} at EPs and during spontaneous oscillations in the *Na*_{i}-variable system. Figure 5 shows a contour plot on the *g*_{f}-*g*_{CaL} parametric plane of stationary-state values of *Na*_{i} at EPs (*A*) and one-parameter bifurcation diagrams in which *Na*_{i}, as well as *V*_{m}, *I*_{f}, and *I*_{NaK}, at stationary states (EPs) and its extrema during spontaneous firings are plotted as functions of *g*_{f} (*B*). The values of *Na*_{i} became higher with increasing *g*_{f}. As illustrated in Fig. 6 in connection with Fig. 1*A*, the *Na*_{i}-variable system with larger *g*_{f} exhibited higher *Na*_{i} both at EPs and during spontaneous firings at a given g_{CaL} along with the bifurcations at larger *g*_{CaL}. The *Na*_{i} reached during stable rhythmic oscillations were much higher than those in stationary states at EPs (Figs. 5*B* and 6).

### Influences of Na_{i} on Bifurcation Structure of SAN Model Cells

Figures 5 and 6 strongly suggest that *g*_{f}-dependent changes in *Na*_{i} at least in part account for the difference in the *I*_{f} effect on the parameter-dependent bifurcation structures of the *Na*_{i}-variable and *Na*_{i}-fixed systems. Because *I*_{NCX}, as well as *I*_{NaK} and other Na^{+} fluxes, depends on *Na*_{i}, stability and bifurcations of SAN cells during parameter changes may be affected by concomitant variations in *Na*_{i}. We therefore investigated how the parameter *Na*_{i} affects stability and bifurcations of EPs and LCs in the *Na*_{i}-fixed system.

#### Increasing Na_{i} shrank parameter regions of unstable EPs and stable LCs.

Figure 7*A* shows a two-parameter bifurcation diagram for *Na*_{i} and *g*_{CaL} where the loci of bifurcations of EPs and LCs are plotted on the *Na*_{i}-*g*_{CaL} parametric plane for the *Na*_{i}-fixed system, incorporating *Na*_{i} as a parameter. The parametric plane was divided into three different areas as in Fig. 3*A*: *1*) unstable EP and rhythmic oscillation (UR), *2*) stable EP and quiescence (SQ), and *3*) coexistence of stable EP and rhythmic oscillation (SR). The *g*_{CaL} regions of unstable EPs (UR) and rhythmic oscillations (UR + SR) dramatically shrank with increasing *Na*_{i}.

#### Increased Na_{i} yielded EP stabilization and quiescence by affecting ionic currents.

As shown in Fig. 7*B*, increased *Na*_{i} enhanced outward *I*_{NaK} and diminished inward *I*_{NCX} (Ca^{2+} efflux via the Na^{+}/Ca^{2+} exchanger) at EPs as well as during spontaneous firings, which resulted in hyperpolarization of *V*_{m} and increase in total intracellular Ca^{2+} content (*Ca*_{total}) [and subspace Ca^{2+} concentration (*Ca*_{sub}), as well] at EPs. With increasing *Na*_{i}, stationary-state *I*_{CaL} was reduced via voltage-dependent deactivation at hyperpolarized *V*_{m} and Ca-dependent inactivation at higher *Ca*_{sub}; thus the combined reductions in both *I*_{CaL} and inward *I*_{NCX} at higher *Na*_{i} resulted in smaller activation of the inward currents (i.e., outward shift of *I*_{total}) on depolarization from stationary-state *V*_{m}, leading to EP stabilization (Fig. 8).

### Impacts of I_{f} on Robustness of SAN Model Cells Against Hyperpolarizing Loads

*I*_{f} is usually expected to enhance the robustness of SAN pacemaking against hyperpolarizing loads by supplying inward currents and thereby preventing excess hyperpolarization of SAN cells. However, as for the *g*_{CaL}-dependent bifurcations, the influences of *I*_{f} on bifurcation structures of the model cells suffering hyperpolarizing loads may be different in the *Na*_{i}-variable and *Na*_{i}-fixed systems: i.e., even if *I*_{f} enhances SAN cell robustness in the *Na*_{i}-fixed system at lower *g*_{f} as reported previously (33), *I*_{f}-dependent changes in *Na*_{i} may eliminate and reverse the enhancing effect of *I*_{f} in the *Na*_{i}-variable system like a real cell. To determine whether *I*_{f} enhances robustness against hyperpolarizing loads of SAN cell pacemaking, therefore, we further examined the influences of modifying *I*_{f} on stability of EPs (stationary states) and LCs (spontaneous firings), as well as their bifurcations, in the *Na*_{i}-variable and *Na*_{i}-fixed model cells during injections of hyperpolarizing *I*_{bias}, increments of *g*_{bK}, and applications of ACh. Essentially the same results were obtained for all these hyperpolarizing loads.

#### Enhancing I_{f} shrank I_{bias} regions of unstable EPs and stable LCs in the Na_{i}-variable system.

Figure 9 shows one-parameter bifurcation diagrams in which stationary-state *V*_{m} at EPs, *V*_{m} extrema of spontaneous firings, and bifurcation points are plotted as functions of the amplitude of hyperpolarizing *I*_{bias} for the *Na*_{i}-variable and *Na*_{i}-fixed systems. In the *Na*_{i}-variable system, the *I*_{bias} regions of unstable EPs and stable LCs (rhythmic firings) shrank with increasing *g*_{f}; spontaneous firings became unstable and arrhythmic (chaotic) via a bifurcation (destabilization) of LCs and vanished via a bifurcation (stabilization) of EPs. In contrast, the unstable EP and stable LC regions of the *Na*_{i}-fixed system were enlarged by *I*_{f} at the relatively small *g*_{f}.

Bifurcation structures during hyperpolarizing loads of the model systems were further tested for broader ranges of *g*_{f}, as illustrated by two-parameter bifurcation diagrams for hyperpolarizing *I*_{bias} and *g*_{f} where the loci of bifurcations of EPs and LCs are plotted on the *g*_{f}-*I*_{bias} parametric plane (Fig. 10). The parametric planes were divided into three different areas as in Fig. 3: *1*) unstable EP and rhythmic firing (UR), *2*) stable EP and quiescence (SQ), and *3*) unstable EP and arrhythmic firing (UA) in the *Na*_{i}-variable system or coexistence of stable EP and rhythmic firing (SR), i.e., bistable zone (BS), in the *Na*_{i}-fixed system. The *I*_{bias} regions of unstable EPs and stable LCs (rhythmic firings) as well as unstable LCs (arrhythmic firings) in the *Na*_{i}-variable system dramatically shrank with increasing *g*_{f} (Fig. 10*A*). In contrast, in the *Na*_{i}-fixed system, the *I*_{bias} region of unstable EPs was first enlarged by *I*_{f} at relatively small *g*_{f} and then shrunk by greater increases in *g*_{f}, whereas that of stable LCs (rhythmic firings) was broadened with increasing *g*_{f} (Fig. 10*B*). The critical *I*_{bias} values to yield bifurcations to quiescence at lower *g*_{f} (control value or less) were larger, but those at higher *g*_{f} were smaller, in the *Na*_{i}-variable system than in the *Na*_{i}-fixed system.

#### Augmented I_{f} caused dramatic increase in Na_{i} during enhancement of hyperpolarizing loads.

As shown in Fig. 11 in connection with Fig. 9*A*, the *Na*_{i}-variable system with larger *g*_{f} exhibited higher *Na*_{i}, as well as higher amplitudes of *I*_{f}, at EPs and during spontaneous firings, along with the bifurcations at smaller hyperpolarizing *I*_{bias}; during increases in amplitudes of *I*_{bias}, *I*_{f} and *Na*_{i} increased more dramatically at higher *g*_{f}. In the *Na*_{i}-fixed system, increasing the parameter *Na*_{i} shrank the *I*_{bias} regions of unstable EPs and stable LCs (rhythmic firings) as for the *g*_{CaL} regions (data not shown); as in the *g*_{CaL}-dependent bifurcation, the increased *Na*_{i} facilitated bifurcations during *I*_{bias} applications and thus significantly shrank the *I*_{bias} regions of unstable EPs and rhythmic firings at larger *g*_{f}.

### Comparisons Among Previously Developed SAN Cell Models

To determine whether the *I*_{f} effects on SAN cell robustness are model dependent or not, we tested and compared *I*_{f}-dependent bifurcation structures during *I*_{CaL} inhibitions or hyperpolarizing loads (*g*_{bK} increments) of previously developed single SAN cell models, including Wilders et al. (53), Demir et al. (6), and Dokos et al. (10) models for transitional SAN cells, as well as Kurata et al. (20, 24) models for central SAN cells. Figure 12 shows two-parameter bifurcation diagrams where bifurcations of EPs are located on the *g*_{f}-*g*_{CaL} and *g*_{f}-*g*_{bK} parametric planes for the *Na*_{i}-variable and *Na*_{i}-fixed versions of each SAN cell model. The *g*_{CaL} regions of unstable EPs shrank with increasing *g*_{f} in all of the *Na*_{i}-variable systems, whereas those in the *Na*_{i}-fixed systems shrank or enlarged depending on the models. As for *g*_{bK}-dependent bifurcations, enhancing *I*_{f} slightly enlarged the *g*_{bK} region of unstable EPs at *g*_{f} lower than the control value in the Wilders et al. and Dokos et al. models but shrank it at *g*_{f} higher than the control value in all the *Na*_{i}-variable systems. The *g*_{f}-dependent changes in bifurcation points were very small in the Demir et al. model, which has a much more negative threshold potential of *I*_{f} activation than other models; positive shifts in the threshold potential yielded *g*_{f}-dependent shrinkage in the region of unstable EPs (data not shown). In contrast, all the *Na*_{i}-fixed versions exhibited *g*_{f}-dependent enlargements of the *g*_{bK} regions of unstable EPs. The critical parameter values to yield the bifurcation of LCs (to arrhythmic dynamics or quiescence) were very close to or varied in parallel with those at the bifurcation of EPs; the parameter regions of stable LCs (rhythmic firings) were also broader at the control *g*_{f} (*g*_{f} = 1) than at higher *g*_{f} (*g*_{f} ≥ 2) in all the *Na*_{i}-variable systems (data not shown).

## DISCUSSION

In this study, the effects of *I*_{f} on robustness of rabbit SAN pacemaker cells were theoretically investigated in connection with the influence of changes in *Na*_{i}. We found that *I*_{f} effects on bifurcation structures of SAN model cells strongly depend on *Na*_{i} and were different in the *Na*_{i}-variable and *Na*_{i}-fixed systems. The preventing effect of *I*_{f} at lower *g*_{f} (around the control value) on cessation of pacemaker activity via *g*_{CaL} decreases or hyperpolarizing loads as observed in the *Na*_{i}-fixed system (33) was eliminated and reversed by *I*_{f}-dependent changes in *Na*_{i}. These findings suggest that *I*_{f} does not necessarily enhance but may attenuate SAN cell robustness against hyperpolarizing loads, enervation of *I*_{CaL}, or other endogenous and exogenous factors, even in physiological *g*_{f} ranges.

### Impacts of I_{f} on Robustness of SAN Cell Pacemaking

#### I_{f} itself does not necessarily enhance but may attenuate robustness of SAN cells.

Inward currents such as *I*_{CaL} and *I*_{Na} have been shown to destabilize an EP of SAN cells, being essentially important for robust pacemaking without annihilation (21, 24). Thus *I*_{f} is also expected to contribute to EP destabilization in SAN cells. As shown in our previous study (25) and in this study (Figs. 3 and 10), however, increasing *I*_{f} did not enlarge but rather shrank the parameter regions of unstable EPs (stationary states), particularly in the *Na*_{i}-variable system. Unlike *I*_{CaL} or *I*_{Na}, therefore, *I*_{f} itself does not prevent but facilitates EP stabilization, not contributing to EP instability; excess *I*_{f} may rather counteract the destabilizing effect of *I*_{CaL} or *I*_{Na} on EPs. *I*_{f} could contribute to EP stabilization mainly via voltage-dependent kinetics, i.e., a rapid decrease (or increase) in inward *I*_{f} on depolarization (or hyperpolarization) from stationary-state *V*_{m} (Fig. 4).

Maltsev and Lakatta (33) reported for their *Na*_{i}-fixed model cell that *I*_{f} enlarged the area of rhythmic firings on the parametric plane of *g*_{CaL} and *P*_{up}, and thus enhanced SAN cell robustness against inhibitions of the membrane and/or SR Ca^{2+} clock. In our study, however, larger *I*_{f} did not enlarge but rather shrank the *g*_{CaL} (and *P*_{up}), hyperpolarizing *I*_{bias}, and *g*_{bK} (and [ACh]) regions of stable LCs (rhythmic firings) via facilitating bifurcations (destabilization) of LCs, particularly in the *Na*_{i}-variable system (Figs. 3 and 10). In addition, *I*_{f} did not significantly strengthen the enhancing effects of SR Ca^{2+} clock on SAN cell robustness (data not shown). Thus this study suggests that *I*_{f} does not necessarily enhance the robustness of pacemaker activity in central or transitional SAN cells of the rabbit. As suggested previously (25), *I*_{f} may contribute mainly to the sympathetic regulation of pacemaker frequency in the central or transitional SAN cells, while contributing to the robust pacemaking against electrotonic loads of the atrium in peripheral SAN cells.

Although we could not find experimental evidence that larger *I*_{f} attenuates robustness of SAN pacemaker cells, it has been reported that overexpression of HCN-encoded pacemaker current silences biological pacemakers derived from guinea pig atrial myocytes as a cautionary note for development of *I*_{f}-based biological pacemakers (30). This observation may reflect that excess *I*_{f} expression yields EP stabilization and thus attenuation of pacemaker cell robustness, supporting our results.

We tested *I*_{f} effects over broad ranges of *g*_{f}, including those much higher than (4–12 times) the control values, to demonstrate the *I*_{f} effects more clearly as well as to illustrate the effects of *I*_{f} overexpression; however, it should be noted that in the *Na*_{i}-variable system, as in a real cell, enhancing *I*_{f} attenuates the robustness of SAN cell pacemaking regardless of whether *g*_{f} is relatively small (in physiological ranges) or larger. Larger *I*_{f} far beyond the physiological density dramatically shrank the parameter regions of spontaneous rhythmic oscillations (Figs. 3 and 10), which is of particular importance for engineering of *I*_{f}-based biological pacemakers where *I*_{f} could be overexpressed to several times the density of native currents, as demonstrated previously (41, 45), given that *I*_{f} overexpression has experimentally been shown to cause cessation of biological pacemaker activity (30). In the *Na*_{i}-variable systems with physiological *g*_{f}, the critical *g*_{CaL} values to yield bifurcations were much smaller than the normal control *g*_{CaL}; nevertheless, *I*_{f}-dependent bifurcation to quiescence or arrhythmic dynamics may actually occur in real SAN cells with such smaller *g*_{CaL}, considering large cell-to-cell variations in *I*_{CaL} density (37) and possible inhibition of *I*_{CaL} by Ca antagonists (*I*_{CaL} blockers) in patients with cardiovascular diseases.

#### I_{f} effects on stability and bifurcation of SAN cells depend on concomitant changes in Na_{i}.

The effects of *I*_{f} on SAN cell robustness, especially that against hyperpolarizing loads, were quantitatively and qualitatively different in the *Na*_{i}-variable and *Na*_{i}-fixed systems, as summarized by the two-parameter bifurcation diagrams (Figs. 3 and 10). As mentioned above, the differences in the *I*_{f} effects between the two systems come from *I*_{f}-dependent changes of *Na*_{i} in the *Na*_{i}-variable system (Figs. 5 and 6); in the *Na*_{i}-fixed model cell, the parameter *Na*_{i} was shown to exert substantial influences on stability and bifurcations of the model cell during changes in other parameters via modulating *I*_{CaL}, *I*_{NaK}, and *I*_{NCX} (Figs. 7 and 8).

Marked differences between the *Na*_{i}-variable and *Na*_{i}-fixed systems are that *1*) the *g*_{CaL}, *I*_{bias}, and *g*_{bK} regions of unstable EPs are broader in the *Na*_{i}-variable system than in the *Na*_{i}-fixed system at lower *g*_{f} but shrink dramatically and become narrower in the *Na*_{i}-variable system at higher *g*_{f} (Figs. 3 and 10); *2*) with increasing *g*_{f}, the *g*_{CaL} regions of stable LCs (rhythmic firings) shrank much more prominently in the *Na*_{i}-variable system (Figs. 1 and 3) but enlarged at lower *g*_{f} and then shrank only slightly in the *Na*_{i}-fixed system (33); and *3*) as *g*_{f} increased, the *I*_{bias} and *g*_{bK} regions of stable LCs (rhythmic firings) shrank in the *Na*_{i}-variable system but enlarged in the *Na*_{i}-fixed system (Fig. 10). The larger parameter region of unstable EPs in the *Na*_{i}-variable system at the smaller *g*_{f} is associated with lower stationary-state *Na*_{i} at EPs, which further decreases during *g*_{CaL} reductions but increases with enhancing hyperpolarizing *I*_{bias} (Figs. 5*A*, 6 and 11); in the *Na*_{i}-fixed system, decreasing the parameter *Na*_{i} enlarged the parameter regions of unstable EPs and spontaneous firings via increasing *I*_{CaL} and inward *I*_{NCX}, as well as decreasing outward *I*_{NaK}, at EPs (Figs. 7 and 8). At lower *g*_{f} (in physiological *g*_{f} ranges), therefore, the decreased *Na*_{i} at EPs may further contribute to enhancement of SAN cell robustness. On the other hand, the greater *I*_{f}-dependent shrinkage of the unstable regions in the *Na*_{i}-variable system is at least in part due to the *I*_{f}-dependent increase in *Na*_{i} (Na^{+} influx) as predicted theoretically (Figs. 5, 6, 7 and 11; see also Ref. 4) and also observed experimentally (4). In addition, the parameter region of unstable EPs was much narrower than that of stable LCs in the *Na*_{i}-fixed system but not (both are comparable) in the *Na*_{i}-variable system (Fig. 3), which is also ascribable to the difference in *Na*_{i} at EPs and during LC oscillations in the *Na*_{i}-variable system: *Na*_{i} at EPs were much lower than those during LC oscillations (Figs. 5*B* and 6), yielding enlargement of the parameter region of unstable EPs (Fig. 7*A*). These results suggest that changes in *Na*_{i} strongly affect stability and bifurcations of SAN cells and thus must be taken into account in experimental and theoretical studies.

#### I_{f} may enhance SAN cell robustness in combination with other inward currents.

The present study suggests that *I*_{f} does not necessarily enhance but may attenuate robustness of SAN cell pacemaking. However, there are many experimental and theoretical reports suggesting the enhancement of SAN cell robustness by *I*_{f} and the *I*_{f}-dependent cardiac pacemaker defined as pacemaking to be abolished by blocking *I*_{f}: *1*) an *I*_{f} blocker abolished spontaneous activity of rabbit SAN cells when hyperpolarizing *I*_{bias} was applied (50); *2*) the background current in the pacemaker potential range was outward before *I*_{f} activation in rabbit SAN cells (7); and *3*) mouse SAN cells lacking HCN4 or HCN2 exhibited sinus dysrhythmia, recurrent sinus pause, and quiescence, especially at low heart rates, e.g., under muscarinic stimulation or in the presence of *I*_{f} blockers (16, 31, 32). The *I*_{f}-based regular pacemaking was also observed in mouse embryonic cardiomyocytes (48, 56), mouse embryonic stem (ES) cell-derived cardiomyocytes (43), and neonatal rat ventricular myocytes (11). Furthermore, previous studies regarding biological pacemaker engineering have demonstrated that *I*_{f}-based biological pacemakers can be created in the atrium and ventricle by HCN gene transfer (3, 42, 49, 54). These experimental findings, suggesting the enhancement of SAN cell robustness by *I*_{f} and the *I*_{f}-dependent pacemaking, are apparently inconsistent with our results in this study. Nevertheless, our previous study (25) suggests that the *I*_{f}-dependent pacemaking is possible under hyperpolarized conditions in SAN cells and that *I*_{f} may enhance the robustness of central and peripheral SAN cells against hyperpolarizing loads in combination with *I*_{st} and *I*_{Na}, respectively. Thus bifurcations leading to *I*_{f}-dependent pacemaking and *I*_{f}-based enhancement of SAN cell robustness may actually occur in SAN (or other regions of the heart) under hyperpolarized or other conditions via the combined actions of *I*_{f} and other inward currents.

As suggested previously (25), whether *I*_{f} enhances robustness of SAN cells strongly depends on the presence of other inward currents (e.g., *I*_{st}, *I*_{Na}) that destabilize EPs (via Hopf bifurcations) and/or eliminate stable EPs (via saddle-node bifurcations). Such inward currents also include the class D *I*_{CaL} (*I*_{CaLD}) and *I*_{CaT} in the central SAN of mouse or other species and embryonic or fetus hearts. For instance, *I*_{CaLD}, which activates in the pacemaker potential range but is not present in the rabbit SAN, may contribute to enhancement of SAN cell robustness and emergence of *I*_{f}-dependent pacemaking in the mouse SAN; indeed, *I*_{f} at lower conductance (physiological densities) enhanced robustness of the *Na*_{i}-variable system modified by incorporating *I*_{CaLD} to simulate mouse SAN cell pacemaking (data not shown). Requirement of the combined effects of *I*_{f} and other inward currents such as *I*_{st} and *I*_{CaLD} was also suggested for enhancement of biological pacemaker functions (23). In fact, the combined action of *I*_{f} and *I*_{CaLD} or *I*_{CaT} as a requisite to generation of spontaneous rhythmic firings has been suggested for mouse atrioventricular cells (36) and ES cell-derived cardiomyocytes (55).

#### Roles of I_{f} in pacemaker generation may be different depending on regions and species.

It is worthwhile noting that *I*_{f} can enhance robustness of pacemaker cells in two different ways (25): *1*) by preventing a saddle-node bifurcation (emergence of a stable EP at more negative *V*_{m}) as in the rabbit central SAN cell with relatively large *I*_{st}, and *2*) by preventing a Hopf bifurcation (stabilization of an EP) as in the rabbit peripheral SAN cell with relatively large *I*_{Na}. Thus *I*_{f} block can yield bifurcations to quiescence of pacemaker cells via a saddle-node or Hopf bifurcation (25). In the Maltsev-Lakatta model, increasing *I*_{f} did not enhance the robustness of pacemaker activity, because hyperpolarizing loads did not cause a saddle-node bifurcation of EPs in this model cell with relatively small *I*_{st}. However, our previous study (25) suggests that *I*_{f} at lower conductance (*g*_{f} smaller than the control value) can enhance the robustness against hyperpolarizing loads of rabbit central SAN cells by preventing a saddle-node bifurcation, which occurs in the presence of relatively large *I*_{st}. The emergence of a saddle-node bifurcation during *I*_{f} inhibition or hyperpolarizing loads is comparable to the report of DiFrancesco (7) for SAN cells, and is also applicable to pacemaker activity of Purkinje fibers, which possess the inward-rectifier K^{+} channel current (*I*_{K1}) at relatively low density (9, 47). Thus *I*_{f} can improve pacemaker cell robustness against hyperpolarizing loads chiefly by supplying inward currents to prevent emergence of a stable stationary state at more negative *V*_{m}. On the other hand, *I*_{f} may enhance the robustness against hyperpolarizing loads of rabbit peripheral SAN cells, as well as SAN cells of other species and spontaneously firing cells of embryonic hearts, by preventing a Hopf bifurcation (EP stabilization), given that they possess relatively large inward currents such as *I*_{Na}, *I*_{CaLD}, and *I*_{CaT}, which are capable of destabilizing EPs (25).

In our previous study on biological pacemakers (23), *I*_{f} facilitated the generation of pacemaker activity during downregulation of *I*_{K1} in human ventricular myocyte model by yielding a saddle-node bifurcation (elimination of a stable EP with the most negative stationary-state *V*_{m}) at higher *I*_{K1} conductance (*g*_{K1}). Using a coupled-cell model, we also showed that *I*_{f} enhanced the robustness (capability of pacemaking and driving nonpacemaker cells) of biological pacemakers against electrotonic loads of nonpacemaker cells by facilitating a saddle-node bifurcation to eliminate a stable EP with the most negative stationary-state *V*_{m} in the nonpacemaker cells during increases in the gap junction conductance, when *g*_{K1} of the surrounding nonpacemaker cells is relatively small (23). However, larger *I*_{f} did not enhance robustness or driving ability of biological pacemaker cells, suggesting again the undesirability of excess expression of *I*_{f} and the requirement of coexpression of *I*_{f} and other inward currents such as *I*_{st} and *I*_{CaLD} for the development of *I*_{f}-based biological pacemakers (30).

### Limitations of Study

As described in our previous articles (24, 25), there are many limitations of our theoretical study, including incompleteness of the models and lack of experimental evidence. The intracellular factors are known to exert substantial effects on ion channel functions and pacemaker dynamics, possibly affecting bifurcation structures, as well. More elaborate models incorporating detailed descriptions of the intracellular Ca^{2+} handling and many other modulating factors have to be developed and tested for future investigations on the roles of *I*_{f} and other components in SAN pacemaking. Another shortcoming of the present study may be the use of classical Hodgkin-Huxley formalism for the kinetic formulation of *I*_{f} and other time-dependent voltage-gated channel currents. More general and realistic schemes have been developed for native *I*_{f} channels (7) and HCN channel isoforms (1) as well as for other individual channels. Thus further sophisticated models based on more general and realistic Markovian state schemes need to be developed and tested in future investigations.

We showed that our conclusions regarding the roles of *I*_{f} could commonly be derived from all the SAN cell models tested in this study (Fig. 12); the effects of *I*_{f} on robustness of SAN cell pacemaking were qualitatively the same regardless of which model was used as a base cell model or which *I*_{f} formula was incorporated (data not shown; see Ref. 25). Nevertheless, the influences of *I*_{f} on bifurcation properties should be tested for future novel (more elaborate) SAN cell and *I*_{f} models for understanding the roles of *I*_{f} in SAN pacemaking more profoundly.

In the present study, we examined the effects of *I*_{f} on robustness of single SAN cells. However, we need investigations of the *I*_{f} effects on pacemaking and driving ability of the SAN in vivo, which require one- to three-dimensional multicellular models, like those used in previous studies (13, 17, 18).

## GRANTS

This work was supported in part by Ministry for Education, Science, Sports and Culture of Japan Grants-in-Aid for Scientific Research (C) 20590220 and 23590266 (to Y. Kurata and T. Shibamoto) and Kanazawa Medical University Grants for Promoted Research S2011-3 and S2012-2 (to Y. Kurata).

## DISCLOSURES

No conflicts of interest, financial or otherwise, are declared by the authors.

## AUTHOR CONTRIBUTIONS

Y.K. conception and design of research; Y.K. performed experiments; Y.K., I.H., M.T., and T.S. analyzed data; Y.K., I.H., and T.S. interpreted results of experiments; Y.K. and M.T. prepared figures; Y.K. drafted manuscript; Y.K., I.H., M.T., and T.S. approved final version of manuscript; I.H. and T.S. edited and revised manuscript.

- Copyright © 2013 the American Physiological Society