## Abstract

To elucidate the roles of hyperpolarization-activated current (*I*_{f}) in sinoatrial node (SAN) pacemaking, we theoretically investigated *1*) the effects of *I*_{f} on stability and bifurcation during hyperpolarization of SAN cells; *2*) combined effects of *I*_{f} and the sustained inward current (*I*_{st}) or Na^{+} channel current (*I*_{Na}) on robustness of pacemaking against hyperpolarization; and *3*) whether blocking *I*_{f} abolishes pacemaker activity under certain conditions. Bifurcation analyses were performed for mathematical models of rabbit SAN cells; equilibrium points (EPs), periodic orbits, and their stability were determined as functions of parameters. Unstable steady-state potential region determined with applications of constant bias currents shrunk as *I*_{f} density increased. In the central SAN cell, the critical acetylcholine concentration at which bifurcations, to yield a stable EP and quiescence, occur was increased by smaller *I*_{f}, but decreased by larger *I*_{f}. In contrast, the critical acetylcholine concentration and conductance of gap junctions between SAN and atrial cells at bifurcations progressively increased with enhancing *I*_{f} in the peripheral SAN cell. These effects of *I*_{f} were significantly attenuated by eliminating *I*_{st} or *I*_{Na}, or by accelerating their inactivation. Under hyperpolarized conditions, blocking *I*_{f} abolished SAN pacemaking via bifurcations. These results suggest that *1*) *I*_{f} itself cannot destabilize EPs; *2*) *I*_{f} improves SAN cell robustness against parasympathetic stimulation via preventing bifurcations in the presence of *I*_{st} or *I*_{Na}; *3*) *I*_{f} dramatically enhances peripheral cell robustness against electrotonic loads of the atrium in combination with *I*_{Na}; and *4*) pacemaker activity of hyperpolarized SAN cells could be abolished by blocking *I*_{f}.

- nonlinear dynamics
- computer simulation
- electrotonic modulation
- biological pacemaker

elucidation of the mechanisms of sinoatrial node (SAN) pacemaking is one of the most important subjects in cardiac electrophysiology. A large body of information on the ionic mechanisms of SAN pacemaking has been obtained from numerous experimental and theoretical studies (27, 35, 41, 55), yielding four principal hypotheses that attribute pacemaker depolarization primarily to *1*) activation of the L-type Ca^{2+} channel current (*I*_{Ca,L}) (58); *2*) deactivation of the delayed-rectifier K^{+} current (*I*_{K}) (39); *3*) development of the hyperpolarization-activated cation current (*I*_{f}) (9); and *4*) spontaneous Ca^{2+} release from sarcoplasmic reticulum (SR) and subsequent activation of the Na^{+}/Ca^{2+} exchanger current (*I*_{NaCa}) (28, 30, 33, 34, 45). Nevertheless, there is no dominant pacemaker mechanism, and the roles of individual time-dependent, voltage-gated channel currents in pacemaker generation remained to be clarified. From the theoretical studies using mathematical models for rabbit SAN cells, we have provided significant insights into the dynamic mechanisms of the natural SAN pacemaking and the roles of *I*_{Ca,L}, *I*_{K}, and Na^{+} channel current (*I*_{Na}) in basal pacemaking (22, 25). We have suggested that *1*) the instability of an equilibrium point (EP) is essentially important for the generation of robust pacemaking without annihilation; *2*) *I*_{Ca,L} is responsible for EP instability and generation of pacemaker activity; *3*) *I*_{K} makes the action potential amplitude larger, eases changes in oscillation frequency during hyperpolarization or depolarization, and plays a pivotal role in preventing bifurcation to quiescence; and *4*) *I*_{Na} contributes to EP instability and robust pacemaking against hyperpolarizing loads in the peripheral SAN. However, the roles of *I*_{f} in SAN pacemaking, remaining the subject of controversy (27), have not yet been determined by the theoretical approach.

*I*_{f}, encoded by the hyperpolarization-activated cyclic nucleotide-gated (HCN) channel gene, is known to contribute to pacemaker depolarization as the pacemaker current (9, 27) and is a key current for engineering of biological pacemakers (43, 46, 48, 53, 56). *I*_{f} contributes to prevention of excess hyperpolarization and excessively slow pacemaking against hyperpolarizing loads (16, 27, 35), autonomic regulations of spontaneous activity (7), stabilization of pacemaker frequency (16, 33, 40), and generation of phase 4 depolarization in the periphery of intact SAN, suffering electrotonic influences of the atrium (38). In isolated SAN cells or the center of intact SAN, however, *I*_{f} does not appear to be indispensable for generation of spontaneous oscillations under normal conditions because *1*) *I*_{f} blockers did not abolish spontaneous activity of real SAN cells or the intact SAN, although the block of *I*_{f} might be incomplete at pacemaker potentials (16, 27, 38, 41); *2*) any model SAN cells did not undergo a bifurcation to quiescence during inhibition of *I*_{f}, exhibiting automaticity, even in the absence of *I*_{f} (21, 55); and *3*) in our preliminary study, the effects of reducing *I*_{f} on stability and oscillation dynamics of model cells were much smaller than those of reducing *I*_{Ca,L} or *I*_{K}. Nevertheless, the most prominent electrophysiological property of SAN cells is the possession of relatively large *I*_{f}, which is very small or negligible in atrial or ventricular myocytes; thus the roles of the pacemaker current *I*_{f} in SAN pacemaking must be understood more profoundly. Like *I*_{Ca,L} and *I*_{Na}, *I*_{f} may contribute to EP destabilization as a requisite for robust pacemaking without annihilation and may also enhance the robustness of pacemaker activity against hyperpolarizing loads.

The aim of this study was to provide more profound insights into the roles of *I*_{f} in the natural pacemaking of SAN cells in terms of the nonlinear dynamics and bifurcation theory. In our laboratory's previous studies, bifurcation structures (i.e., ways of changes in the number or stability of equilibrium and periodic states) of ventricular or SAN model systems were investigated for elucidating the dynamic mechanisms of normal and abnormal pacemaker activities (22–25). 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. In this study, therefore, we performed bifurcation analyses for mathematical models of the central and peripheral SAN cells (21, 25, 59). Bifurcation diagrams were constructed by calculating EPs, periodic orbits called limit cycles (LCs), their stability, and bifurcation points as functions of model parameters, such as the maximum conductance of *I*_{f}. We particularly focused on the effects of *I*_{f} on the stability of steady-state membrane potentials and robustness of pacemaker activity against hyperpolarizing loads and thus evaluated stability and oscillation dynamics of the model cells during injections of constant bias current (*I*_{bias}), applications of acetylcholine (ACh), or electrotonic modulations by the atrium. Furthermore, we explored whether and how *I*_{f}-dependent pacemaking, which is defined as the pacemaker activity to be abolished by blocking *I*_{f}, is possible. These theoretical analyses provide significant insights into the contributions of *I*_{f} to EP instability and robust pacemaking of SAN cells.

## THEORY AND METHODS

### Mathematical Formulation

#### Base models for central and peripheral SAN cells.

As base pacemaker models of the rabbit SAN, we used the Kurata et al. central (21) and peripheral (25) cell models, as well as the Zhang et al. (59) models modified by Garny et al. (12) and Kurata et al. (25). A complete list of the equations and standard parameter values for their models, as well as their dynamics, has been provided in the previous articles; these models are implemented in a CellML-based open resource for public access (http://COR.physiol.ox.ac.uk/).

The base models for the normal pacemaker activity include 14 membrane current components. The membrane current system includes *I*_{Ca,L}, *I*_{Na}, *I*_{f}, T-type Ca^{2+} channel current (*I*_{Ca,T}), sustained inward current (*I*_{st}), rapidly activating (*I*_{Kr}) and slowly activating (*I*_{Ks}) components of *I*_{K}, 4-AP-sensitive currents consisting of transient (*I*_{to}) and sustained (*I*_{sus}) components, background currents carried by Na^{+} (*I*_{b,Na}) and K^{+} (*I*_{b,K}), muscarinic K^{+} channel current (*I*_{K,ACh}), Na^{+}-K^{+} pump current (*I*_{NaK}), and *I*_{NaCa}. All of the currents charge the membrane capacitance (*C*_{S}). Then, the time-dependent change in the membrane potential (*V*_{S}) is described by

#### I_{f} formulas.

Three different *I*_{f} schemes and formula sets of the Kurata et al. (21), modified Zhang et al. (59) and Maruoka et al. (37) models were used for the central cell; one of the three *I*_{f} formula sets was incorporated into the *I*_{f}-removed Kurata et al. model used as the base model. For the peripheral cell, the Kurata et al. (25) and Zhang et al. (59) models were tested. The results obtained from the central cell models with different *I*_{f} formula sets were qualitatively the same; those from the two peripheral cell models were essentially the same as well. Thus the results from the Kurata et al. central (21) and peripheral (25) models are shown in the results section. The schemes, formulas, and voltage-dependent kinetics of the three *I*_{f} models are given in online data supplements (Expanded Theory and Methods and the Supplemental Fig. S1; the online version of this article contains supplemental data).

#### Incorporation of ACh effects on ionic currents.

To investigate the bifurcation structure of the model cells during applications of ACh, we incorporated the formulas of Zhang et al. (60) for *I*_{K,ACh}, as well as for the modifications of *I*_{Ca,L} and *I*_{f} by ACh, into the base models. In this study, *I*_{K,ACh} density and sensitivities to ACh of *I*_{K,ACh}, *I*_{Ca,L}, and *I*_{f} were assumed to be the same in central and peripheral cells (for more details, see Refs. 25 and 60). The expressions for these currents with the ACh modifications of Zhang et al. (60) are given in the online data supplements.

#### Formulation of a coupled-cell model.

To investigate the electrotonic influences of atrial myocytes on stability and dynamics of SAN cells, we employed a coupled-cell model, as introduced by Watanabe et al. (54). According to their method, a peripheral SAN cell model was connected to a passive membrane model [resistance-capacitance circuit] for an atrial myocyte via the gap junction conductance (G_{C}) of 0–1,000 nS. Watanabe et al. (54) used the atrial membrane model with a capacitance (*C*_{A}) of 51 pF and a resistance (*R*_{A}) of 220–500 MΩ. However, crista terminalis cells were reported to have relatively large capacitance of 134 ± 22 pF and input resistance of 495 ± 227 MΩ (57). We, therefore, used the *C*_{A} of 134 pF and *R*_{A} of 100–900 MΩ for the membrane model of an atrial myocyte (crista terminalis cell). A peripheral SAN model cell was connected to one to three atrial model cells with an identical G_{C}. For comparison, a central SAN model cell (20–32 pF) was connected to an atrial model cell with the same *C*_{A} and *R*_{A} and tested for electrotonic influences. An offset voltage, corresponding to a resting potential of the atrial myocyte, was set equal to −80 mV.

Time-dependent changes in the membrane potentials of the SAN (*V*_{S}) and atrium (*V*_{A}) were calculated by the equations
*I*_{total} and *I*_{A} represent the sum of sarcolemmal ionic currents in the SAN cell (pA/pF) and the atrial membrane current (pA), respectively. *I*_{C} denotes the gap junction current (pA).

### Bifurcation Analysis

The SAN cell models are 15- or 26-order autonomous continuous-time dynamical systems, each of which is described as a set of 15 or 26 first-order, nonlinear ordinary differential equations. Dynamical properties and bifurcation structures of model systems were determined by analytically and numerically handling the differential equations. For bifurcation analyses, the intracellular K^{+} concentration was fixed at 140 mM to remove degeneracy (for details, see Refs. 20 and 23). 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 (25). Bifurcation parameters chosen in this study include *1*) the maximum conductance of *I*_{f} (g_{f}); *2*) activation time constant of *I*_{f} (τ_{y}); *3*) inactivation time constant of *I*_{st} (τ_{qi}) and *I*_{Na} (τ_{h}); *4*) amplitude of *I*_{bias}; *5*) ACh concentration ([ACh]); and *6*) G_{C}. The maximum conductance and gating time constants of the ionic channel currents are expressed as normalized values, i.e., ratios to the control values. Detailed procedures for *1*) locating EPs and LCs; *2*) constructing one- and two-parameter bifurcation diagrams; and *3*) detecting codimension-one bifurcations (determination of EP and LC stabilities) are provided in our laboratory's previous article (25) and in the online data supplements (Expanded Theory and Methods). Definitions of the specific terms for the nonlinear dynamics and bifurcation theory are also given there. We used *1*) Newton-Raphson algorithm (a nonlinear equation solver, *fsolve*, available in optimization toolbox of MATLAB) to locate EPs, as well as to detect bifurcations of EPs; *2*) brute-force approach using a MATLAB ODE solver, *ode15s*, to calculate stable LCs and arrhythmic (period-*k* periodic, quasi-periodic, or chaotic) dynamics; and *3*) autonomous shooting method (42) using the MATLAB ODE solver or CL_MATCONT, a continuation toolbox for MATLAB (8), to locate unstable LCs and detect bifurcations of LCs. CL_MATCONT is a powerful tool for continuation of EPs and LCs; however, it is not capable of computing arrhythmic dynamics or could not detect the points of bifurcation of LCs in our models for unknown reasons. Thus we used the combined method with multiple tools. The type of bifurcations of LCs and whether oscillation behavior is chaotic or not were determined by calculating the characteristic multiplier and the first Lyapunov exponent, respectively (42, 47).

For evaluation of the stability of membrane potentials, as well as the robustness of pacemaker activity against hyperpolarizing and depolarizing loads, the steady-state membrane potential (*V*_{E}) at EPs corresponding to the zero-current crossings of the steady-state current-voltage (*I*-*V*) curve and its stability were determined with *I*_{bias} being applied at various amplitudes to change *V*_{E}. *I*_{bias}-*V*_{E} curves were constructed by plotting *V*_{E} as functions of the amplitude of *I*_{bias}, with saddle-node (SN) and Hopf (H) bifurcation points located on the curves (see Figs. 1 and S2; for more details, see Ref. 25). To further assess the robustness of pacemaker activity against hyperpolarizations via parasympathetic stimulations or electrotonic modulations, we also determined the critical [ACh] or G_{C} value at which a bifurcation to quiescence or arrhythmic dynamics occurs (for details, see Expanded Theory and Methods in the online data supplements).

## RESULTS

### Influences of *I*_{f} on Bifurcation Structures of SAN Model Cells

As summarized in previous articles (21, 55, 59), eliminating *I*_{f} did not abolish, but only slowed, pacemaker activity of the real and model SAN cells. In our central or peripheral model cells, complete block of *I*_{f}, although slowing pacemaking, did not either stabilize EPs or abolish spontaneous activity. Thus *I*_{f} is not required for generation of basal SAN pacemaking under normal conditions. Nevertheless, this current is expected to destabilize EPs at hyperpolarizing potentials and enhance the robustness of pacemaker activity against hyperpolarizing loads, as well as play a central role in phase 4 depolarization under hyperpolarized conditions, as in pacemaking of Purkinje fibers (11). We, therefore, examined the influences of modulating *I*_{f} on the stability of steady-state potentials (EPs) during *I*_{bias} applications and bifurcation structures during ACh applications or electrotonic modulations of the model cells, with g_{f} increased from zero to the control value.

#### Influences of I_{f} on EP stability and bifurcations during I_{bias} applications.

EP instability is essentially important for generation and robust maintenance of automaticity: if an EP is stable, a system could come to a rest at the stable EP; annihilation or pause of spontaneous activity may occur (e.g., via premature atrial beat), if a system is in a bistable zone over which a stable EP and LC coexist (15, 17, 22, 44). Thus we first examined the influences of *I*_{f} on stability and bifurcation of EPs during *I*_{bias} applications, i.e., whether *I*_{f} enlarges the unstable *V*_{E} (and *I*_{bias}) region in the *I*_{bias}-*V*_{E} curve, for the model cells. This approach would reveal the contributions of *I*_{f} to EP (*V*_{E}) instability, as well as robustness against hyperpolarizing and depolarizing loads (22–25, 52).

Figure 1 shows the effects of increasing g_{f} on stability and bifurcation during *I*_{bias} applications of the Kurata et al. central (21) and peripheral (25) cell models. The results for the Kurata et al. (21) central model incorporating the *I*_{f} formulas of Maruoka et al. (37) and for the Zhang et al. (59) peripheral cell model are also shown in Supplemental Fig. S2, for comparison. Bifurcation points on the *I*_{bias}-*V*_{E} curve were determined with changing g_{f} and plotted as functions of g_{f} on both the potential (*V*_{E}) and current (*I*_{bias}) coordinates. In the central cell, increasing g_{f} did not enlarge but slightly shrunk the unstable *V*_{E} region between two H bifurcation points (H_{1} and H_{2}) by shifting the bifurcation point H_{2} toward more positive potentials. During g_{f} increases, the *I*_{bias} region for unstable *V*_{E} was first enlarged, with the SN bifurcation point shifting toward more negative *I*_{bias} values; however, this effect was limited because further increases of g_{f} caused a H bifurcation at smaller hyperpolarizing *I*_{bias}. As in the central cell model, increasing g_{f} slightly shrunk the unstable *V*_{E} regions above H_{2} (between H_{1} and H_{2}) of the peripheral cell model. Different from the central cell, the *I*_{bias} region for unstable *V*_{E} progressively enlarged with increasing g_{f}.

#### Effects of I_{f} on bifurcations during ACh applications.

The *I*_{f}-induced enlargement of the hyperpolarizing *I*_{bias} region with unstable EPs (Figs. 1 and S2, *bottom*) suggest that *I*_{f} enhances the SAN cell robustness against hyperpolarization; spontaneous activity of the model cells tested was abolished at the points very close to the H or SN bifurcation points (data not shown). To assess the contribution of *I*_{f} to the SAN cell robustness against hyperpolarizing loads, therefore, we examined the effects of *I*_{f} on bifurcation structures of the model cells during ACh applications. One-parameter bifurcation diagrams to illustrate the model cell stability and oscillation dynamics during [ACh] increases were first constructed for different g_{f} values. We further examined the g_{f}-dependent changes in H and SN bifurcation points, i.e., the critical [ACh] value at which a stable EP emerges, and thus constructed two-parameter bifurcation diagrams for g_{f} and [ACh]; the critical [ACh] at which oscillation dynamics became arrhythmic (period-*k* periodic, quasi-periodic, or chaotic) via bifurcations of LCs were also determined as functions of g_{f}.

Figure 2*A* shows one-parameter bifurcation diagrams for [ACh], illustrating bifurcation structures (changes in EPs and LCs, their stability, and bifurcation points) during [ACh] increases of the *I*_{f}-reduced systems, as well as the control system, for the original Kurata et al. (21) central model. During increases of [ACh], EPs on the steady-state branch (*V*_{E}) of the central cells with g_{f} = 0, 0.25, and 1.0 were stabilized via H bifurcations (H_{1}). In the normal cell, LCs were destabilized via a period-doubling bifurcation (PD_{1}); spontaneous oscillations became arrhythmic, abruptly shrunk in amplitude via another period-doubling bifurcation (PD_{2}), and finally vanished at the H bifurcation, as more clearly shown in Supplemental Fig. S3. The PD_{1} bifurcation spawned *period 2* periodic orbits, followed by the emergence of chaotic behavior (a sequence of periodic and chaotic windows), which was demonstrated by calculating the first Lyapunov exponent (see Supplemental Figs. S4 and S5). In the *I*_{f}-removed cell, a LC became unstable via a Neimark-Sacker bifurcation (NS) with emergence of arrhythmic (period-*k* periodic, quasi-periodic, or chaotic) dynamics vanishing in the vicinity of a H bifurcation (H_{2}), following a SN bifurcation (SN_{1}). There was a very steep growth in the cycle length (CL) of LCs just before termination of unstable periodic branches, suggesting the existence of homoclinic bifurcations (*hom* in Fig. 2*A*), which was also reported for a previous *I*_{f}-removed SAN model (15). The *I*_{f}-reduced cell (g_{f} = 0.25), but not the normal cell, exhibited a cascade of PD bifurcations (*period 2*, *4*, and *8* periodic orbits), followed by chaotic behavior (i.e., PD route to chaos).

Bifurcation structures during [ACh] increases were also determined for the normal and *I*_{f}-reduced Kurata et al. (25) peripheral model (Fig. 2*B*). During [ACh] increases, EPs of the peripheral cell with g_{f} = 0, 0.5, and 1.0 were stabilized via H bifurcations. In the normal cell (g_{f} = 1) and with g_{f} = 0.5, LCs became unstable via NS bifurcations, with arrhythmic (period-*k* periodic, quasi-periodic, or chaotic) dynamics emerging at these bifurcations and abruptly vanishing at the H bifurcations (for details, see Supplemental Fig. S6). In the *I*_{f}-removed cell, although bifurcations of LCs could not be detected, there was a very steep growth in the CL of LCs just before cessation of spontaneous oscillations (in the vicinity of the H bifurcation), suggesting the existence of a homoclinic bifurcation. In contrast to the result from the central cell, the critical [ACh] values at the bifurcations to arrhythmic dynamics or quiescence became greater as g_{f} increased. The cells with higher g_{f} exhibited more stable CL (smaller g_{f}-dependent increases in CL) during increases of [ACh].

Figure 3 illustrates two-parameter bifurcation diagrams for g_{f} and [ACh] in which the critical [ACh] values for H and SN bifurcations of EPs, as well as PD and NS bifurcations of LCs, to occur are plotted as functions of g_{f}. The parametric planes were divided into three different areas: *1*) unstable EP and rhythmic (*period 1* periodic) oscillation (UR); *2*) unstable EP and arrhythmic (period-*k* periodic, quasi-periodic, or chaotic) oscillation (UA); and *3*) stable EP and quiescence (no oscillation) (SQ). Note that H and SN bifurcation points on the steady-state branch (locus of EPs) of one-parameter bifurcation diagrams are just or very close to the point at which spontaneous activity ceases. In the central model, increasing g_{f} first enlarged the [ACh] region of unstable EP (UR + UA) via positively shifting the SN bifurcation point (SN_{1}), and then slightly shrunk it via a negative shift of the H bifurcation point (H_{1}). Similarly, the region UR first enlarged but then shrunk as g_{f} increased, whereas the region UA first shrunk and then enlarged with increasing g_{f}. In contrast to the result for the central cell, increasing g_{f} progressively enlarged the [ACh] region of unstable EP of the peripheral cell by increasing the critical [ACh] to cause a H bifurcation. Both the regions UR and UA broadened with increasing g_{f} in a complex manner.

#### Effects of I_{f} on bifurcations during electrotonic modulations.

For further clarifying how *I*_{f} contributes to the robustness of peripheral SAN cells against electrotonic loads of the atrium, we also examined the effects of *I*_{f} on bifurcation structures of the coupled-cell model during G_{C} increases. To first determine the size effect (influence of the number of atrial cells) and influence of *R*_{A} on bifurcations, H bifurcation points were determined as functions of *R*_{A} for a SAN model cell connected to one to three atrial model cells with identical G_{C} (Supplemental Fig. S7A). The experimentally determined *R*_{A} of 495 ± 227 MΩ was too large for a single atrial model cell to significantly affect stability and bifurcation of a peripheral SAN model cell with a relatively large capacitance of 65 pF; at least two atrial model cells are required for yielding H bifurcations. *I*_{f} effects on bifurcations, as well as the electrotonic influence of the atrial cell, strongly depended on *R*_{A}; the g_{f}-dependent increase in bifurcation values was apparently attenuated with smaller *R*_{A} (Supplemental Fig. S7B).

Figure 4 shows two parameter bifurcation diagrams for G_{C} and g_{f}, where the critical G_{C} values at which H bifurcations of EPs and NS bifurcations of LCs occur are plotted against g_{f} for a peripheral SAN model cell connected to two atrial cells of *R*_{A} = 268–495 MΩ. The areas UR, UA, and SQ are determined for the g_{f}-G_{C} plane as for the g_{f}-[ACh] plane (Fig. 3). In the periphery, increasing g_{f} progressively enlarged the G_{C} region of unstable EP (UR + UA) and the region UR by increasing the critical G_{C} for H and NS bifurcations, respectively. In contrast, increasing g_{f} slightly but continuously shrunk the G_{C} region of unstable EP by reducing bifurcation values of G_{C} in the center (data not shown). Although the degree of the *I*_{f} effects depended on the values of total *C*_{A} and *R*_{A}, the results were qualitatively the same for different *C*_{A} and *R*_{A}.

One-parameter bifurcation diagrams for G_{C}, illustrating bifurcation structures during G_{C} increases of the normal and *I*_{f}-reduced SAN model cell connected to two atrial model cells of *R*_{A} = 268 (MΩ), are shown in Fig. 5. During G_{C} increases, EPs of the peripheral cell with g_{f} = 0, 0.3, and 1.0 were stabilized via H bifurcations, with the cell coming to a rest. In both the normal and *I*_{f}-reduced cells, LCs were destabilized via NS bifurcations with the emergence of arrhythmic dynamics between the NS and H bifurcations. As shown in Supplemental Fig. S8, period-*k* periodic and chaotic (or quasi-periodic) windows alternately appeared after the NS bifurcations. In the arrhythmic regions, the model cells exhibited subthreshold prepotentials and skipped-beat runs (Supplemental Fig. S9), which are very similar to the behaviors previously observed in other SAN models (15, 33).

#### Influences of I_{f} gating kinetics on I_{f}-induced changes in bifurcations during hyperpolarization.

It is known that stability and bifurcations of a pacemaker system depend on gating kinetics of ion channels (22). *I*_{f} has relatively slow gating kinetics; the slow gating of *I*_{f} may be critical and advantageous to its beneficial effects on SAN cell robustness against hyperpolarization. Thus we further examined how the time constant for *I*_{f} activation gating variable (τ_{y}) influences *I*_{f}-induced changes in bifurcations during ACh applications or electrotonic modulations. Figure 6 shows two-parameter bifurcation diagrams for g_{f} and [ACh] or G_{C}, which were constructed for various *I*_{f} with different τ_{y} (1–0.01 times the control). In the center, decreasing τ_{y} by a factor of 100 slightly enhanced the shrinkage in the [ACh] region of unstable EP via reducing H bifurcation values at higher g_{f}, while not affecting the SN bifurcation at lower g_{f} (Fig. 6*A*). In the periphery, *I*_{f} with the decreased τ_{y} could not enlarge the [ACh] or G_{C} region of unstable EP as significantly as the normal *I*_{f} (Fig. 6*B*). The influence of *I*_{f} gating was greater on the G_{C} increase-induced bifurcation than on the [ACh] increase-induced one; with the smaller τ_{y} of 0.01, increasing g_{f} did not enlarge but shrunk the G_{C} region of unstable EP.

### Combined Effects of *I*_{f} and Other Currents on Bifurcation Structures of SAN Cells

Inward currents other than *I*_{f} also contribute to SAN pacemaking (35, 41). In particular, *I*_{st} and *I*_{Na} are present at high density in the central and peripheral cells, respectively, with their voltage range of activation being different from, but overlapping with, that of *I*_{f}; thus *I*_{f} may contribute to SAN pacemaking in combination with these currents. Therefore, we further examined the combined effects of *I*_{f} and *I*_{st} or *I*_{Na} on bifurcation structures of the model cells during hyperpolarization.

#### Combined effects of I_{f} and I_{st} on central cell bifurcation during hyperpolarization.

Figure 7*A* shows two-parameter bifurcation diagrams in which the critical [ACh] to yield SN and H bifurcations were plotted as functions of g_{f}. Shown are the loci of bifurcation points during g_{f} increases for the normal and *I*_{st}-removed central cell model, as well as for the model cell with the τ_{qi} being 100–1,000 times smaller than the control values. Decreasing τ_{qi} caused a reduction in the critical [ACh] for the H bifurcation, but not for the SN bifurcation, and thus further limited the *I*_{f}-induced enlargement of the [ACh] region of unstable EP without substantially affecting the *I*_{f} effect on H bifurcations. In the *I*_{st}-removed system, the critical [ACh] to yield a stable EP was relatively low, a SN bifurcation as observed in the *I*_{st}-incorporated system at lower g_{f} did not occur, and increasing g_{f} did not significantly enlarge the [ACh] region of unstable EP.

#### Combined effects of I_{f} and I_{Na} on peripheral cell bifurcation during hyperpolarization.

The *I*_{f}-induced enlargements of the [ACh] and G_{C} region of unstable EP were much greater in the periphery than in the center. *I*_{Na}, which is present only in the periphery, may be responsible for the greater effect of *I*_{f}. We, therefore, determined the critical [ACh] and G_{C} at H bifurcations as functions of g_{f} for the normal and *I*_{Na}-removed peripheral cell model, as well as for the model cell with the τ_{h} being 10–100 times smaller than the control values (Fig. 7*B*). Decreasing τ_{h} caused a reduction in the critical [ACh] and G_{C} for H bifurcations; the increases in bifurcation values during g_{f} increases were much greater in the normal system than in the system with *I*_{Na} of accelerated inactivation. The removal of *I*_{Na} further shrunk the [ACh] and G_{C} region of unstable EP. In the *I*_{Na}-removed model cell and the model cell with accelerated *I*_{Na} inactivation (τ_{h} = 0.1–0.01), the G_{C} region of unstable EP was not significantly enlarged by g_{f} increases.

### Searching for *I*_{f}-dependent Pacemaking under Hyperpolarized Conditions

Figures 3 and 4 imply the existence of the [ACh] and G_{C} region, where blocking *I*_{f} causes a SN or H bifurcation, i.e., where the normal cell but not the *I*_{f}-reduced cell exhibits spontaneous oscillations around an unstable EP. This suggests that *I*_{f}-dependent pacemaking, which is defined as the pacemaker activity to be abolished by blocking *I*_{f}, is possible under hyperpolarized conditions (see Supplemental Fig. S10). Thus we examined whether *I*_{f}-dependent pacemaking occurs in the model cells hyperpolarized by applications of ACh.

#### I_{f} inhibition causes SN bifurcation to quiescence in hyperpolarized central cell.

Figure 8*A* shows EP stability, oscillation dynamics, and bifurcations during *I*_{f} block of the central cell model in the absence and presence of ACh. Under the normal condition (with [ACh] = 0), reducing *I*_{f} did not yield a SN or H bifurcation to create a stable EP, not abolishing LCs. Under the hyperpolarized conditions, however, blocking *I*_{f} led to *1*) de novo creation of EPs at more negative potentials (EP_{2} and EP_{3} in Supplemental Fig. S10A) via SN bifurcations (SNs in Fig. 8*A*); *2*) destabilization of LCs via a PD bifurcation (PD in Fig. 8*A*, *right*) with emergence of *period 2* and *4* periodic orbits, followed by an alternate sequence of periodic and chaotic windows (see Supplemental Fig. S11); and *3*) cessation of spontaneous activity via H bifurcations in the vicinity of the SN bifurcation points, with the cell coming to a rest at the stable EP.

#### I_{f} inhibition causes H bifurcation to quiescence in hyperpolarized peripheral cell.

We also determined EP stability, oscillation dynamics, and bifurcations during *I*_{f} block of the peripheral cell model in the absence and presence of ACh (Fig. 8*B*). In the normal cell, reducing *I*_{f} did not cause stabilization of an EP or abolition of LCs. In the cell hyperpolarized by ACh applications, however, blocking *I*_{f} caused *1*) negative shifts of *V*_{E} (EP in Supplemental Fig. S10B) with its stabilization at hyperpolarizing potentials via H bifurcations (H in Fig. 8*B*); *2*) destabilization of LCs via PD bifurcations, followed by successive PDs and sequences of period-*k* periodic and chaotic dynamics (Supplemental Fig. S12); and *3*) cessation of spontaneous activity at the H bifurcation points with the cell coming to a rest at the stable EP.

## DISCUSSION

### Roles of *I*_{f} in SAN Pacemaking

#### I_{f} itself does not destabilize an EP.

Inward currents, such as *I*_{Ca,L} and *I*_{Na}, have been shown to destabilize an EP of SAN cells, being essentially important for robust pacemaking without annihilation (22, 25). Thus *I*_{f} is also expected to contribute to EP instability of SAN cells. In both the central and peripheral models, however, increasing *I*_{f} did not enlarge, but rather shrunk, the unstable potential region in the *I*_{bias}-*V*_{E} curve determined with *I*_{bias} applications (Figs. 1 and S2). Unlike *I*_{Ca,L} or *I*_{Na}, therefore, *I*_{f} itself cannot destabilize (but may even stabilize) an EP, not contributing to EP instability; excess *I*_{f} may rather counteract the destabilizing effect of *I*_{st} or *I*_{Na} on EPs. EP destabilization may be yielded only by the inward current, with rapid activation and relatively slow inactivation kinetics on depolarization (25).

#### I_{f} enhances SAN cell robustness to hyperpolarizing loads by preventing bifurcations.

Our results (Figs. 2 and 3) suggest that *I*_{f} at lower conductance enhances the robustness of central SAN cells against ACh-induced hyperpolarization (parasympathetic stimulation) by preventing emergence of a stable EP via SN and subsequent H bifurcations. Nevertheless, the *I*_{f} effect on the central cell was relatively small and limited by occurrence of a H bifurcation at more positive steady-state potentials. The central cell robustness to hyperpolarizing loads was not enhanced, but slightly attenuated at higher g_{f}, possibly indicating that *I*_{f} density should be relatively small in the central region of pacemaker systems.

In contrast to the central cell, the peripheral cell exhibited the continuous enlargements of the [ACh] and G_{C} regions of unstable EP and rhythmic oscillation during enhancement of *I*_{f}, although the arrhythmic region might also be enlarged by *I*_{f} (Figs. 3 and 4). Thus *I*_{f} may enhance the robustness of peripheral SAN cells against hyperpolarizing loads, such as parasympathetic stimulation and electrotonic modulations, by preventing stabilization of an EP (H bifurcation), as well as destabilization of LCs (NS bifurcation). Note the regional difference in the way of enhancing SAN cell robustness to parasympathetic stimulation: preventing a SN bifurcation in the center, but preventing a H bifurcation in the periphery (Fig. 3).

Although there is no direct experimental evidence that *I*_{f} enhances SAN robustness against hyperpolarizing loads by preventing bifurcations, previous studies revealed sinus dysrhythmia, recurrent sinus pause, and cessation of spontaneous activity in mice lacking HCN2 or HCN4 (16, 31, 32), possibly reflecting enhancement of the stability of SAN pacemaking by *I*_{f}. The repetitive sinus pause was especially prominent at low heart rates, e.g., under muscarinic stimulation (low sympathetic tone) or in the presence of *I*_{f} blockers (16, 32). These arrhythmic behaviors observed experimentally are very similar to those of the hyperpolarized model cells reproduced in this study (Supplemental Figs. S5, S9, and S11). Maltsev et al. (33) reported that *I*_{f}, enlarging the area of rhythmic firing in the parameter plane of *I*_{Ca,L} conductance and sarcoplasmic reticulum Ca^{2+} pumping rate, enhances the robustness of their SAN model system. Thus *I*_{f} may enhance SAN cell robustness not only to hyperpolarizing loads, but also to other endogenous or exogenous factors.

#### Regional differences in I_{f} effects suggest different roles of I_{f} in center and periphery.

*I*_{f}-induced enhancement of the SAN cell robustness was limited and relatively small in the center, while continuous and relatively large in the periphery (Figs. 3 and 4). *I*_{f} appears to play a pivotal role in preventing a bifurcation to quiescence or arrhythmic dynamics of the peripheral cell, especially that via electrotonic modulations. The greater effect of *I*_{f} on the peripheral robustness to electrotonic loads, as well as higher *I*_{f} density in the periphery (59), are reasonable because in vivo peripheral cells directly suffer the electrotonic load of adjacent atrial myocytes (18, 19, 51, 54) and thus must be more robust to electrotonic modulations than central cells, which are distant from the crista terminalis. These regional differences in the *I*_{f} effects may reflect different roles of *I*_{f} in the center and periphery of the rabbit SAN: *I*_{f} may contribute mainly to the robust pacemaking against electrotonic loads of the atrium in the periphery, but mainly to the sympathetic regulation of pacemaker frequency in the center.

A marked difference in the response to modulating *I*_{f} between the central and peripheral models was that the [ACh] and G_{C} regions of unstable EP, where spontaneous oscillations emerge, were enlarged in the periphery but shrunk in the center by increasing *I*_{f}. We found in the Kurata et al. (21, 25) models that the *I*_{f}-induced increases in H bifurcation values were reversed in the center and limited in the periphery via concomitant increases of the critical [Na^{+}]_{i} (data not shown). Thus the *I*_{f}-induced shrinkage of the unstable region in the central model appears to be due to the [Na^{+}]_{i} increase during enhancement of *I*_{f}, along with the relatively small *I*_{f} effects.

#### Slow kinetics of I_{f} is crucial for enhancement of SAN cell robustness to electrotonic loads.

At least two mechanisms can be assumed for the *I*_{f}-induced enlargements of the [ACh] and G_{C} regions of unstable EP: *1*) supply of inward currents, and *2*) enhancement of EP instability via the relatively slow gating kinetics. The supply of inward currents appears to be essential for the *I*_{f} effects, because *I*_{f} itself did not enhance the instability of *V*_{E} (Figs. 1 and S2). In the center, accelerating *I*_{f} gating little affected the SN bifurcation during ACh applications (Fig. 6*A*); in the periphery, *I*_{f} with accelerated kinetics was still effective in enlarging the [ACh] region of unstable EP (Fig. 6*B*). The current of very rapid kinetics is regarded as a time-independent background current. Thus these findings further support that *I*_{f} improves SAN cell robustness against parasympathetic stimulations chiefly by the supply of inward currents to prevent hyperpolarization of *V*_{E} or emergence of a stable *V*_{E} at more negative potentials. In contrast, slow-gating kinetics of *I*_{f} appear to be crucial for the enhancement of peripheral cell robustness to electrotonic loads of the atrium (Fig. 6*B*); slow deactivation of *I*_{f} on depolarization is probably indispensable for retaining EP instability, although *I*_{f} itself does not destabilize EPs. Further slowing the kinetics of *I*_{f} did not significantly enhance its effects (data not shown), suggesting that the native *I*_{f} kinetics are nearly optimal for enhancing the SAN cell robustness. HCN channel isoforms have been shown to exhibit distinct gating properties (1–3), which are different from those of native *I*_{f} channels, possibly reflecting their distinct roles in SAN pacemaking. These findings may particularly be important in engineering of *I*_{f}-based biological pacemakers by transfer of modified HCN genes for the treatment of bradyarrhythmias: HCN4 exhibits slower activation kinetics (3, 31, 32), and thus may be more suitable for biological pacemakers, than other isoforms.

#### I_{st} and I_{Na} are involved in the mechanisms for I_{f} enhancement of SAN cell robustness.

This study suggest that the *I*_{f} effects strongly depend on *I*_{st} in the center and *I*_{Na} in the periphery (Fig. 7). *I*_{f} appears to enhance the central cell robustness against parasympathetic stimulation and the peripheral cell robustness against hyperpolarizing loads (especially to electrotonic loads of the atrium) in combination with *I*_{st} and *I*_{Na}, respectively. The *I*_{f} effects on the periphery were also shown to be strongly influenced by the inactivation kinetics of *I*_{Na} (Fig. 7*B*). These results suggest that *I*_{Na}-induced destabilization of an EP at hyperpolarizing *V*_{E}, as suggested previously (25), is involved in the mechanisms of *I*_{f}-induced enhancement of the peripheral cell robustness to hyperpolarizing loads, especially to electrotonic modulations by the atrium. Our laboratory's previous report (25) also suggested the contribution of *I*_{Na} to high robustness of the peripheral SAN against electrotonic modulations. Thus the combined effects of *I*_{f} and *I*_{Na} may be indispensable for prevention of EP stabilization and robust maintenance of SAN pacemaking against electrotonic loads of the atrium.

#### I_{f}-dependent pacemaking is possible in hyperpolarized cells.

Our results further suggest that *I*_{f}-dependent pacemaking to be abolished by blocking *I*_{f} is possible under hyperpolarized conditions (Figs. 8 and S10–S12). There are many experimental and theoretical reports suggesting the *I*_{f}-dependent cardiac pacemaker: *1*) Cs^{+}, an *I*_{f} blocker, abolished spontaneous activity of rabbit SAN cells when hyperpolarizing *I*_{bias} was applied (49); *2*) the instantaneous background current in the pacemaker potential range was outward before *I*_{f} activation in rabbit SAN cells (9); *3*) HCN4-deficient mouse SAN cells were quiescent under low-cAMP conditions (16). Previous studies regarding the biological pacemaker engineering also demonstrated that *I*_{f}-dependent biological pacemakers can be created in the atrium and ventricle by HCN gene transfer (43, 48, 56). Thus bifurcations leading to *I*_{f}-dependent pacemaking may actually occur in the SAN or other regions of the heart under hyperpolarized or other nonphysiological conditions.

It may be worthwhile noting that *I*_{f} block yielded bifurcations to quiescence of the central and peripheral models in different ways, i.e., via a SN bifurcation (emergence of a stable EP at more negative potentials) in the center, but via a H bifurcation (stabilization of an EP) in the periphery. Emergence of a SN bifurcation in the central cell is comparable to the report of DiFrancesco (9). The presence of *I*_{Na} is crucial for the generation of *I*_{f}-dependent pacemaking in hyperpolarized peripheral SAN cells, because the G_{C} and [ACh] ranges over which the peripheral cell exhibits *I*_{f}-dependent pacemaking are much broader in the *I*_{Na}-incorporated system than in the *I*_{Na}-removed one (Fig. 7*B*).

### Significance of Nonlinear Dynamical Approach in Cardiac Electrophysiology

#### Construction of general theory for pacemaker mechanisms.

Despite numerous experimental studies, the contribution of *I*_{f} to SAN pacemaking remained the subject of controversy (27); thus the experimental approach has not been able to elucidate the roles of *I*_{f} in pacemaker generation. As suggested by Guevara and Jongsma (15), the theoretical approach–bifurcation analysis–would be necessary for properly appreciating results of experimental studies, as well as for clarifying contributions of individual currents to pacemaker mechanisms. For instance, apparent controversial experimental findings that pacemaker activity was abolished by eliminating *I*_{f} in some studies but not in others may totally be accounted for from the aspect of bifurcation phenomena: a bifurcation to quiescence may occur during *I*_{f} inhibition in some preparations, but not in others, depending on the net background current and other factors.

Our laboratory's previous (22, 25) and present study further indicate that significant insights into the roles of individual time-dependent currents, such as *I*_{Ca,L}, *I*_{Kr}, *I*_{Na}, and *I*_{f}, in natural SAN pacemaking, can be obtained from the nonlinear dynamical study based on bifurcation theory. Moreover, as reported previously, investigating bifurcation structures of SAN or ventricular models is also useful for elucidating the dynamical mechanisms of abnormal automaticity (4, 5, 15, 23, 29, 50, 52) and triggered activity (14). Thus the nonlinear dynamical approach would be necessary for general understanding and systematic descriptions of the mechanisms of normal and abnormal pacemaker activities, possibly leading to the construction of a general theory for pacemaker mechanisms.

#### Application to engineering of biological pacemakers.

Our works also suggest that exploring bifurcation structures of model cells allows us to predict and control dynamics and bifurcation structures of real cells to a certain degree. The mathematical approach may provide a theoretical background for the engineering of biological pacemakers from native quiescent cardiomyocytes in vivo or from human embryonic stem cells ex vivo for gene or cell therapy of bradyarrhythmias (e.g., sick sinus syndromes, atrioventricular block), requiring implantation of electronic pacemakers (for review, see Refs. 13 and 36).

Recent studies regarding gene therapy for bradyarrhythmias have focused on the genetic engineering of *I*_{f}-based biological pacemakers in the atrium or ventricle (43, 46, 48, 56). This study on the roles of *I*_{f} in SAN pacemaking may, therefore, be useful as a theoretical background for the development of *I*_{f}-based biological pacemakers with robust pacemaking and driving, like the SAN system: e.g., our findings suggest that relatively slow gating kinetics of *I*_{f} are favorable for *I*_{f}-based biological pacemakers (Fig. 6), and that *I*_{st} or *I*_{Na} is necessary for *I*_{f} to contribute to robust pacemaking against hyperpolarizing loads (Fig. 7). As already mentioned, distinct kinetic properties of HCN isoforms (1–3) may particularly be important and must be taken into account for the development of *I*_{f}-based biological pacemakers, because *I*_{f} gating kinetics strongly influenced its own effects on bifurcation properties of SAN cells. In future studies, we would like to explore the effects of incorporating various types of native or modified HCN channels on pacemaking and driving ability of pacemaker systems to determine which HCN isoform is most favorable (and how it should be modified) for the development of *I*_{f}-based biological pacemakers.

### Limitations of Study

#### Incompleteness of base SAN and I_{f} models.

As described in our laboratory's preceding article (25), there are many limitations of our theoretical study, including incompleteness of the models and lack of experimental evidence. We assumed that basal pacemaking of SAN cells is attributable to the interplay of sarcolemmal ionic currents, such as *I*_{Ca,L}, *I*_{Na}, *I*_{Kr}, and *I*_{f} (i.e., membrane clock). However, SR Ca^{2+} release and intracellular Ca^{2+} dynamics (i.e., Ca^{2+} clock) may play a central role in pacemaker generation and in balancing robustness and flexibility of pacemaker activity (27, 28, 33, 34). 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 the membrane clock, as well as the Ca^{2+} clock, in SAN pacemaking.

Another shortcoming of the present study may be the use of the classical Hodgkin-Huxley formalism or simplified Markovian state scheme 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 (10) 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 for SAN cells in future investigations.

#### Model dependences of I_{f} effects.

The effects of *I*_{f} block on SAN pacemaking frequency were so much different among model systems, as well as experimental studies (21, 55). One may, therefore, suppose that our conclusions derived on the roles of *I*_{f} are applicable only to the SAN and *I*_{f} models used in this study. The degree of the *I*_{f} effects on bifurcation properties of the SAN model cells was certainly greater, as *I*_{f} activation occurred at more positive potentials and thus affected oscillation frequency more strongly. However, the effects of different *I*_{f} models on bifurcation structures of the model cells were, although quantitatively different, qualitatively the same, regardless of how they affect oscillation frequency or which model was used as a base SAN model. Nevertheless, the influences of *I*_{f} on bifurcation properties should be tested for future novel SAN cell and *I*_{f} models for understanding the roles of *I*_{f} in SAN pacemaking more profoundly.

#### Requirement of more realistic multicellular models.

In the present study, we examined the effects of *I*_{f} on SAN cell robustness to electrotonic loads of the atrium using the passive atrial membrane model, but not more realistic atrial myocyte models. We believe that the coupled-cell model used in this study is useful in exploring electrotonic modulations of SAN pacemaking. However, the influences of *I*_{f} on capability of the SAN to drive the atrium should be examined by using a coupled-cell model incorporating an atrial myocyte (not membrane) model. Furthermore, 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 (12, 18, 19).

## GRANTS

This work was supported in part by Ministry for Education, Science, Sports and Culture of Japan Grant-in-Aid for Scientific Research (C) 17590192 and 20590220 (to Y. Kurata and T. Shibamoto), Kanazawa Medical University Grant for Promoted Research S2008-3 and S2009-1 (to Y. Kurata), and Kanazawa Medical University Grant for Collaborative Research C2009-6 (to Y. Kurata).

## DISCLOSURES

No conflicts of interest, financial or otherwise, are declared by the author(s).

- Copyright © 2010 the American Physiological Society