## Abstract

To elucidate the dynamical mechanisms of the sinoatrial (SA) node pacemaker activity, we investigated the roles of L-type Ca^{2+} (*I*_{Ca,L}) and delayed-rectifier K^{+} (*I*_{Kr}) currents in pacemaking by stability and bifurcation analyses of our rabbit SA node model (Kurata Y, Hisatome I, Imanishi S, and Shibamoto T. *Am J Physiol Heart Circ Physiol* 283: H2074–H2101, 2002). Equilibrium points (EPs), periodic orbits, stability of EPs, and Hopf bifurcation points were calculated as functions of conductance or gating time constants of the currents for constructing bifurcation diagrams. Structural stability (robustness) of the system was also evaluated by computing stability and dynamics during applications of constant bias currents (*I*_{bias}). Blocking *I*_{Ca,L} or *I*_{Kr} caused stabilization of an EP and cessation of pacemaking via a Hopf bifurcation. The unstable zero-current potential region determined with *I*_{bias} applications, where spontaneous oscillations appear, shrunk and finally disappeared as *I*_{Ca,L} diminished, but shrunk little when *I*_{Kr} was eliminated. The reduced system, including no time-dependent current except *I*_{Ca,L}, exhibited pacemaker activity. These results suggest that *I*_{Ca,L} is responsible for EP instability and pacemaker generation, whereas *I*_{Kr} is not necessarily required for constructing a pacemaker cell system. We further explored the effects of various K^{+} currents with different kinetics on stability and dynamics of the model cell. The original *I*_{Kr} of delayed activation and inward rectification appeared to be most favorable for generating large-amplitude oscillations with stable frequency, suggesting that *I*_{Kr} acts as an oscillation amplifier and frequency stabilizer. *I*_{Kr} may also play an important role in preventing bifurcation to quiescence of the system.

- rabbit
- primary pacemaker cell
- nonlinear dynamics
- bifurcation diagram
- computer simulation

elucidation of the mechanisms of the initiation of the spontaneous heartbeat, i.e., pacemaker activity of the sinoatrial (SA) node, is one of the most important subjects for cardiac electrophysiology. A large body of information on the ionic mechanisms of SA node pacemaking has been obtained from numerous experimental and theoretical studies. So far, three principal hypotheses for generation of pacemaker activity have been proposed, attributing the pacemaker depolarization primarily to *1*) activation of the L-type Ca^{2+} channel (slow inward) current (69), *2*) deactivation of the delayed-rectifier K^{+} current [referred to as the “K^{+} current conductance (*g*_{K})-decay hypothesis”; Refs. 7, 16, 46], and *3*) development of the hyperpolarization-activated current (14). However, the roles of individual ionic components in the SA node pacemaking remain the subject of controversy, with no consensus on the identity of a single dominant pacemaker current (19, 32, 54). From theoretical analysis using a mathematical model for single SA node cells, Dokos et al. (16) concluded that the combined effect of a number of sarcolemmal currents is responsible for phase 4 depolarization and that the time-independent background Na^{+} current (*I*_{b,Na}) is the dominant inward (pacemaker) current during this phase, with its removal leading to cessation of rhythmic activity. We also obtained the same conclusions in our simulations using an improved SA node model (37): all the current systems directly or indirectly contribute to pacemaker depolarization, with *I*_{b,Na} being the most dominant inward current during phase 4. However, a system not including a time-dependent current never exhibits pacemaker activity; thus the roles of individual time-dependent currents in pacemaker generation must be considered.

The aim of this study was to elucidate the role of each time-dependent current in the normal (natural) pacemaking of SA node primary pacemaker cells in terms of the nonlinear dynamics and bifurcation theory. Local stability and bifurcation analyses, as well as numerical simulations, were performed for an updated mathematical model of the primary pacemaker cell (37). The L-type Ca^{2+} channel current (*I*_{Ca,L}) and rapidly activating delayed-rectifier K^{+} current (*I*_{Kr}) appear to be essentially important for pacemaker generation in primary pacemaker cells of the rabbit SA node, because blocking one of these currents abolishes spontaneous activity (31, 34, 37, 59, 64). In contrast, other time-dependent currents [T-type Ca^{2+} channel current, hyperpolarization-activated current, sustained inward current, 4-aminopyridine (4-AP)-sensitive currents, and slowly activating delayed-rectifier K^{+} current] would be inessential to pacemaker generation, only modulating oscillation frequency or action potential waveforms: in experimental or simulation studies, spontaneous activity did not cease even when these currents were removed or not incorporated (19, 37); in our preliminary study, reducing the maximum conductance or changing gating time constants of these currents did not cause a bifurcation to quiescence, with their effects on the stability, dynamics, and bifurcation structure of the model cell much smaller than those of *I*_{Ca,L} or *I*_{Kr}. In this study, therefore, we focused on the roles of *I*_{Ca,L} and *I*_{Kr} in normal pacemaking. We constructed bifurcation diagrams by calculating equilibrium points (EPs), periodic orbits, stability of the EP, and Hopf bifurcation points as functions of the bifurcation parameters to characterize the dynamic properties of *I*_{Ca,L} and *I*_{Kr} (i.e., the maximum conductance and gating time constants). The “structural stability” of the system (robustness of pacemaker activity) was also evaluated by computing the stability and dynamics of the model cell during applications of constant bias currents (*I*_{bias}). Furthermore, to clarify the roles of *I*_{Kr} in regulating pacemaker activity, we tested the effects of various types of K^{+} currents with different kinetics on the stability and dynamics of the model system.

A salient point of this study is that the theory of bifurcations of nonlinear dynamical systems was used to obtain significant insights into the mechanisms of the SA node pacemaker activity. As suggested by Guevara and Jongsma (25), the bifurcation theory would be necessary to properly appreciate the results of numerical simulations or experimental findings on the pacemaker mechanisms. Stability and bifurcation analyses of the mathematical model would allow us to define the contribution of each current to the SA node pacemaking and to identify key variables or parameters responsible for pacemaker generation. We believe that this paper provides significant insights into the roles of *I*_{Ca,L} and *I*_{Kr} in the normal pacemaking of primary pacemaker cells.

## THEORY AND METHODS

### Mathematical Model

We used our model for the single primary pacemaker cell of the rabbit SA node (37), which can reproduce the dynamic properties and bifurcation structures of real SA node cells more accurately than other existing models (for details, see Ref. 37). A complete list of the equations and standard parameter values was provided previously (37).

The complete model for the normal pacemaking includes 12 membrane current components. The differential equation for the membrane potential (*V*) is (1) where *I*_{Ca,L} and *I*_{Ca,T} represent the L-type and T-type Ca^{2+} channel currents, respectively, and the rapid and slow components of the delayed-rectifier K^{+} current (*I*_{K}) are denoted as *I*_{Kr} and *I*_{Ks}, respectively. The membrane current system also includes the transient (*I*_{to}) and sustained (*I*_{sus}) components of 4-AP-sensitive currents, hyperpolarization-activated current (*I*_{h}), sustained inward current (*I*_{st}), background Na^{+} current (*I*_{b,Na}), muscarinic K^{+} channel current (*I*_{K,ACh}), Na^{+}-K^{+} pump current (*I*_{NaK}), and Na^{+}/Ca^{2+} exchanger current (*I*_{NaCa}), charging the membrane capacitance (*C*_{m}). Because *I*_{st} is not necessarily present in spontaneously beating SA node cells (60), we tested both *I*_{st}-incorporated and *I*_{st}-removed systems (see Ref. 37). The results obtained from the two systems were qualitatively the same; therefore, only those from the *I*_{st}-incorporated system are shown in this article.

Intracellular concentrations of Na^{+} ([Na^{+}]_{i}) and K^{+} ([K^{+}]_{i}) were fixed at 10 and 140 mM, respectively, simulating the whole cell perforated-patch recording for internally perfused SA node cells. Extracellular concentrations of Na^{+}, K^{+}, and Ca^{2+} were set equal to 140, 5.4, and 2 mM, respectively. To reduce the number of variables and speed up computations, the rapid buffering approximation was used for describing the buffering effects of the fast Ca^{2+} buffers calmodulin and calsequestrin (29, 65).

### Stability and Bifurcation Analyses

*Constructing bifurcation diagrams.* We examined how the stability and dynamics of the model cell alter with changes in bifurcation parameters and constructed bifurcation diagrams for one or two parameters. Bifurcation parameters chosen in this study were the maximum conductance [conductances of L-type Ca^{2+} current (*I*_{Ca,L}) (*g*_{Ca,L}) and *I*_{Kr} (*g*_{Kr})] or gating time constants [for activation gating variable (*d*_{L}) for *I*_{Ca,L} (τ_{dL}), voltage-dependent inactivation gating variable (*f*_{L}) for *I*_{Ca,L} (τ_{fL}), activation gating variable (*p*_{a}) for *I*_{Kr} (τ_{pa}), inactivation gating variable (*p*_{i}) for *I*_{Kr} (τ_{pi})] of the two channels expressed as ratios to the original (control) values, amplitude of *I*_{bias}, and conductance of leak currents (*g*_{leak}). In our normal pacemaker model with fixed [Na^{+}]_{i} and [K^{+}]_{i}, 22 variables define a 22-dimensional state point in the 22-dimensional state space of the system. We calculated EPs and periodic orbits in the state space and also determined the asymptotic stability of an EP by computing 22 eigenvalues of a 22 × 22 Jacobian matrix derived from the linearization of the nonlinear system around the EP (for more details, see Ref. 62). For the codimension one-bifurcation analysis to construct the one-parameter bifurcation diagram, a bifurcation parameter was systematically changed while keeping all the other parameters at their standard values; EPs (steady-state branches) and local extrema of periodic orbits (periodic branches) were plotted against the bifurcation parameter. Hopf bifurcation points at which the stability of an EP reverses were detected by the stability analysis as described above. In the codimension two-bifurcation analysis to construct the two-parameter bifurcation diagram, the path of Hopf bifurcation points was traced in the parameter plane.

Specific terms for the nonlinear dynamics and bifurcation theory are defined as follows (see also Refs. 21, 24, 39, 57). Equilibrium point (EP) is a time-independent steady-state point in the phase space of a system at which the system permanently stays unless any perturbations are applied, constructing the steady-state branch in a bifurcation diagram. This state point corresponds to the zero-current crossing in the steady-state current-voltage (*I-V*) curve, i.e., a quiescent state of a cell if it is stable.

The periodic orbit is a closed trajectory in the phase space of a system, constructing the periodic branch in a bifurcation diagram. The limit cycle is a periodic limit set onto which a trajectory is asymptotically attracted. A stable limit cycle corresponds to an oscillatory state, i.e., pacemaker activity, of a cell.

Bifurcation is a qualitative change in a solution of differential equations caused by altering parameters, e.g., a change in the number of EPs or periodic orbits, a change in the stability of an EP or periodic orbit, and a transition from a periodic to quiescent state. The bifurcation phenomena we can see in cardiac myocytes include cessation or generation of pacemaker activity and occurrence of abnormal (irregular) dynamics such as skipped-beat runs and early afterdepolarizations (see Ref. 25). Hopf bifurcation is a bifurcation at which the stability of an EP reverses with emergence or disappearance of a periodic solution (limit cycle), occurring when the eigenvalues of a Jacobian matrix for the EP have a single complex conjugate pair and the sign of its real part changes. Saddle-node bifurcation (SNB) is a bifurcation at which two EPs (steady-state branches) or two periodic solutions (periodic branches) coalesce and disappear (or a bifurcation yielding de novo creation of two new EPs or periodic orbits). The SNB of EPs occurs when one of eigenvalues of a Jacobian matrix is zero.

According to Parker and Chua (52), the structural stability of a dynamical system is determined by whether an infinitesimal perturbation of the system alters the qualitative behavior of the system's solutions: a dynamical system is structurally stable if the dynamics never alter qualitatively with an infinitesimal change in a parameter (for more a rigorous definition, refer to Refs. 23, 36, 67). Typical examples of a qualitative change in the system dynamics are a change in the stability of a limit set (EP or periodic orbit) and the creation or disappearance of a limit set.

Several textbooks are available for further details on the concepts of bifurcation theory and practical numerical methods for bifurcation analyses (see, e.g., Refs. 23, 36, 52, 67).

*Definition and evaluation of structural stability for SA node cell system.* The structural stability of a dynamical system is determined by whether infinitesimal perturbation qualitatively alters the system dynamics: the system is structurally stable if an infinitesimal change in a parameter never causes a qualitative change in the dynamics. According to this definition, biological systems including SA node cells are always structurally stable under normal (physiological) conditions. In this article, however, we redefine the term structural stability to be more convenient for cardiac electrophysiology. We assumed that the structural stability of the SA node cell system could be defined as the robustness of pacemaker activity to various interventions or modifications that may cause a bifurcation to quiescence or irregular dynamics and that *system A* is structurally more stable than *system B* when the change of a parameter required for causing a bifurcation is greater in *system A* than in *system B*. Interventions or modifications leading to quiescence or irregular dynamics include injections of *I*_{bias} (1), electrotonic loads of the atrium (58, 61, 66), current leakage via myocardial injury (49), and intrinsic changes in channel conductance or gating kinetics. In this study, the structural stability of our SA node model was tested for application of *I*_{bias}, which may mimic the electrotonic interaction with atrial myocytes or the change in the net background current. Bifurcation structures of model systems during *I*_{bias} applications have often been explored in previous theoretical studies on the mechanisms of abnormal automaticity to develop in ventricular myocytes or SA node cells (10, 25, 62). How the structural stability to *I*_{bias} was determined is illustrated in Fig. 1, in which the changes in the zero-current potential (*V*_{0}) as *V* at an EP and its stability during applications of *I*_{bias} are depicted as the *I*_{bias}-*V*_{0} curve corresponding to the steady-state membrane *I-V* relation. We compared the standard system and the modified system with *I*_{Kr} activation 100 times faster than the control; these two systems are clearly different in structural stability to *I*_{bias}. The “control” EP (*V*_{0}), defined as the EP (*V*_{0}) at *I*_{bias} = 0, is unstable because two eigenvalues of a Jacobian matrix for the EP are positive real numbers. There are two Hopf bifurcation points (referred to as H1 and H2), where real parts of a pair of the complex conjugate eigenvalues go through zero and reverse their sign. In the whole interval between H1 and H2, the system has no stable EP but has only one unstable EP, generating limit cycle oscillations. It should be noted that the bifurcation value of *I*_{bias} at which a Hopf bifurcation or quiescence occurs is greater in the standard system than in the modified system with accelerated *I*_{Kr}; one can say that the standard system is structurally more stable than the modified system.

### Numerical Integration (Dynamic Simulation)

Dynamic behavior of the model cell was determined by solving the simultaneous nonlinear ordinary differential equations numerically. We used a fourth-order adaptive-step Runge-Kutta algorithm that includes an automatic step size adjustment based on an error estimate or a variable time-step numerical differentiation approach selected for its suitability for stiff systems. The maximum relative error tolerance for our integration methods was set to 1 × 10^{–6}. Numerical computations were performed on Power Macintosh G4 computers (Apple Computers; Cupertino, CA) with MAT-LAB 5.2 (MathWorks; Natick, MA).

The action potential amplitude (APA) as a voltage difference between the maximum diastolic potential (MDP) and peak overshoot potential (POP), as well as the cycle length (CL), were determined for each calculation of a cycle. Numerical integration was continued until the differences in both APA and CL between the newly calculated cycle and the preceding one became <1 × 10^{–3} of the preceding APA and CL values. When periodic behavior was irregular or unstable, model dynamics were computed for 30 s; all potential extrema and CL values were then plotted in the diagram.

## RESULTS

### Effects of Conductance of I_{Ca,L} and I_{Kr} on Stability and Dynamics of SA Node Model

We first examined the effects of reducing the maximum conductance of *I*_{Ca,L} (*g*_{Ca,L}) or *I*_{Kr} (*g*_{Kr}) on EP stability and oscillation dynamics of the model cell system. Figure 2 shows the bifurcation diagrams constructed for *g*_{Ca,L} or *g*_{Kr}, depicting the *V*_{0} at an EP (steady-state branch) and local potential extrema of MDP and POP (periodic branches), as well as CL, as functions of *g*_{Ca,L} or *g*_{Kr}. *V*_{0} negatively or positively shifted with decreasing *g*_{Ca,L} or *g*_{Kr}; the EP became stable via a Hopf bifurcation when *g*_{Ca,L} and *g*_{Kr} reduced to 0.3757 and 0.1192, respectively. During the decreases in *g*_{Ca,L} or *g*_{Kr}, a limit cycle oscillation of the model cell gradually contracted in size and finally disappeared at the Hopf bifurcation points.

### Effects of Conductance of I_{Ca,L} and I_{Kr} on Structural Stability of SA Node Model

We also examined the effects of reducing *g*_{Ca,L} or *g*_{Kr} on the structural stability of the model system to applications of *I*_{bias}. Figure 3 shows the displacement of H1 and H2 in the *I*_{bias}-*V*_{0} curve as well as the shift in the control *V*_{0}, with decreases in *g*_{Ca,L} or *g*_{Kr}. H1 and H2 were plotted as functions of *g*_{Ca,L} or *g*_{Kr} on both the potential (*V*_{0}) and current (*I*_{bias}) coordinates. This approach is essentially the same as that introduced by Vinet and Roberge (62) for bifurcation analyses of their ventricular model (see Fig. 12 in Ref. 62). The unstable *V*_{0} and *I*_{bias} regions between H1 and H2, shrinking with decreasing *g*_{Ca,L}, disappeared via a codimension two-saddle node bifurcation at *g*_{Ca,L} = 0.3376. In contrast, reducing *g*_{Kr} only slightly shrunk the unstable *V*_{0} range. The unstable regions did not disappear even when *I*_{Kr} was completely blocked. How EP (*V*_{0}) stability and oscillation dynamics of the *I*_{Ca,L}-reduced or *I*_{Kr}-removed system change during applications of I_{bias} is shown in Fig. 4. As *g*_{Ca,L} diminished, the positive real parts of eigenvalues [Re(λ)] and limit cycle oscillations dramatically attenuated, finally vanishing at the SNB point. Pacemaker activity never appeared in the *I*_{Ca,L}-removed system. In contrast, eliminating *I*_{Kr} did not significantly decrease the positive Re(λ) values. The *I*_{Kr}-removed system, quiescent at *I*_{bias} = 0, resumed spontaneous activity when an EP was destabilized at H1 by applications of hyperpolarizing *I*_{bias}.

### Effects of Gating Kinetics of I_{Ca,L} and I_{Kr} on Stability and Dynamics of SA Node Model

To identify the key variables or parameters responsible for generation and regulation of pacemaker activity, we further explored the effects of changing gating kinetics (time constants of the voltage-dependent activation and inactivation) of *I*_{Ca,L} (τ_{dL}, τ_{fL}) or *I*_{Kr} (τ_{pa}, τ_{pi}) on EP stability and oscillation dynamics of the model cell. For inactivation of *I*_{Ca,L}, we focused only on the voltage-dependent inactivation (*f*_{L}); in the preliminary study, modifying or removing the Ca^{2+}-dependent inactivation (*f*_{Ca}) did not cause a bifurcation or substantially change the stability, dynamics, or bifurcation structure of the model system. In Fig. 5, the model cell dynamics characterized by MDP, POP, CL, and the maximum upstroke velocity (MUV), as well as the stability of an EP (*V*_{0} = –25.22 mV), are shown as functions of the time constants. Increasing τ_{dL} or decreasing τ_{fL} led to stabilization of the EP (Hopf bifurcation) with a cessation of limit cycle oscillations. In contrast, changing τ_{pa} or τ_{pi} did not stabilize the EP or abolish spontaneous activity. Decreasing τ_{pa} or increasing τ_{pi}, however, caused the decrease in APA and MUV. In the vicinity of the control conditions, CL was especially sensitive to τ_{fL} but relatively insensitive to the others.

### Effects of Gating Kinetics of I_{Ca,L} and I_{Kr} on Structural Stability of SA Node Model

Figure 6 shows how changes in the gating time constants of *I*_{Ca,L} or *I*_{Kr} influence the structural stability to *I*_{bias} applications of the model cell. With increasing τ_{dL}, the unstable regions in the *I*_{bias}-*V*_{0} curve, where stable pacemaking is possible, shrank and finally disappeared via a codimension two-saddle node bifurcation. Similarly, decreasing τ_{fL} markedly shrank the unstable regions, although not causing a SNB. In contrast, changing τ_{pa} or τ_{pi} exerted relatively small effects on the unstable *V*_{0} and *I*_{bias} ranges: decreasing τ_{pa} or increasing τ_{pi} shrank the unstable regions but did not lead to a SNB or a Hopf bifurcation of the control EP.

### Stability and Automaticity of a Reduced System Not Including I_{K}

As shown in Fig. 2, pacemaker activity was abolished by blocking *I*_{Kr}; thus *I*_{Kr} appears to be necessary for pacemaker generation. Does this, however, indicate that *I*_{Kr} is indispensable for constructing a pacemaker cell system? We tested whether the reduced system not including *I*_{Kr} (or *I*_{Ks}) can exhibit spontaneous activity. The *I*_{K}-eliminated system resumed pacemaker activity in response to increasing *I*_{K,ACh}, incorporating the background K^{+} current of the linear *I-V* relation (*I*_{b,K}), exposure to hyperpolarizing *I*_{bias}, or electrotonic interaction with the atrium (data not shown here, but that *I*_{b,K} or hyperpolarizing I_{bias} can destabilize an EP and induce spontaneous activity in the *I*_{K}-removed system is illustrated in Figs. 8 and 10). Furthermore, as shown in Fig. 7, we explored the stability and oscillatory behaviors of a model system reduced by eliminating all the time-dependent currents except *I*_{Ca,L}. In control with [Na^{+}]_{i} = 5 mM, the reduced model cell was quiescent at a stable *V*_{0} of –5.89 mV. When *I*_{b,K} was increasingly added, however, *V*_{0} negatively shifted with increasing *I*_{b,K} conductance (*g*_{b,K}) and was eventually destabilized via a Hopf bifurcation at *g*_{b,K} = 7.073 pS/pF. Stable limit cycles similar to the natural pacemaker activity emerged when an EP became unstable.

### Influences of K^{+} Current Kinetics on Dynamic Properties of SA Node Model

Figure 7 suggests that the time-dependent *I*_{K} is not necessarily required for constructing a pacemaker cell system. What, then, is the role of *I*_{K} in SA node pacemaking? To answer this question, we studied the stability and dynamics of various modified systems incorporating different K^{+} current formulas. Both the original *I*_{Kr} and *I*_{Ks} were first removed from the standard system, and a time-dependent *I*_{K} with different gating kinetics or a time-independent background K^{+} current was then incorporated (for details, see Table 1). Figure 8 shows how EP stability and automaticity of the model systems change with increasing K^{+} conductance (*g*_{K}) up to 100 pS/pF (see also Table 2 for comparison). All the K^{+} currents, except the ultra-rapidly-activating delayed-rectifier K^{+} current (*I*_{Kur})-type current, could create an unstable EP and induce limit cycle oscillations via a Hopf bifurcation. The effects of the various K^{+} currents on dynamic properties of the model system are summarized as follows. *1*) The original *I*_{Kr} formulas (*I*_{Kr,F}, *I*_{Kr,S}) yielded large APA and high MUV with stable CL. *2*) The accelerated *I*_{Kr} with decreased τ_{pa} and *I*_{Ks}-type current formulated by eliminating *I*_{Kr} inactivation produced the oscillation of relatively small APA and low MUV. *3*) The slowed *I*_{K} caused unstable CL and irregular dynamics. *4*) The time-independent *I*_{Kr,b} or *I*_{b,K} caused unstable CL and irregular dynamics, with a shrinkage in the *g*_{K} region where an EP is unstable and spontaneous oscillations appear.

### Influences of K^{+} Current Kinetics on Structural Stability of SA Node Model

As shown in Fig. 9, we further examined how the various K^{+} currents affect the structural stability of the model cell to applications of *I*_{bias} or depolarizing loads of leak currents. Hopf bifurcation points in the *I*_{bias}-*V*_{0} or *g*_{leak}-*V*_{0} curve were determined for individual K^{+} currents with *g*_{K} increased up to 100 pS/pF. The effects of the various K^{+} currents on the structural stability to *I*_{bias} or *g*_{leak} are summarized as follows (see also Table 2). *1*) The original *I*_{Kr} could continuously enlarge the unstable *V*_{0}, *I*_{bias}, and *g*_{leak} regions during *g*_{K} increases, sufficiently improving the structural stability to *I*_{bias} and *g*_{leak}. *2*) The accelerated *I*_{Kr}- or *I*_{Ks}-type current, leading to a shrinkage in the unstable *V*_{0} region, could not enlarge the unstable *I*_{bias} or *g*_{leak} region as sufficiently as the original *I*_{Kr} could. *3*) With the *I*_{Kur}-type current, the unstable *V*_{0} and *I*_{bias} regions shrank and eventually disappeared via a SNB; this current could not create the unstable *g*_{leak} region. *4*) The time-independent *I*_{Kr,b} or *I*_{b,K} could not enlarge the *I*_{bias} or *g*_{leak} region of instability; with increasing *g*_{b,K}, the unstable regions shrank and finally disappeared via a SNB. Figure 10 shows the stability and dynamics during *I*_{bias} applications of the model systems incorporating the different K^{+} currents. When the original *I*_{Kr}- or *I*_{Ks}-type current was incorporated, stable limit cycles appeared with relatively stable CL over the large *I*_{bias} range. In contrast, the system with accelerated *I*_{Kr} or time-independent K^{+} currents exhibited unstable CL and irregular dynamics in the relatively small *I*_{bias} region of instability.

## DISCUSSION

### Nonlinear Dynamical Aspects of Pacemaker Mechanisms

To elucidate the dynamical mechanisms of SA node pacemaker activity, we theoretically investigated the roles of *I*_{Ca,L} and *I*_{Kr} in pacemaking by stability and bifurcation analyses of a primary pacemaker cell model. The term “mechanisms of pacemaker activity” is usually interpreted as meaning how individual currents drive or contribute to phase 4 depolarization (7, 46, 47, 54). From the point of view of the nonlinear dynamics and bifurcation theory, however, pacemaker activity could be defined as a stable limit cycle around an unstable EP, and thus the most important dynamical property of the SA node cell system for stable pacemaking seems to be the “instability” of an EP: if an EP is stable, the system will be quiescent at the stable EP or in a bistable zone where stable steady and periodic states coexist and thus “annihilation” occurs (see Ref. 25). When considering the pacemaker mechanisms, therefore, one should first understand why the *V*_{0} of SA node pacemaker cells is unstable, unlike the stable resting potential of atrial or ventricular myocytes. Our study further suggests that initiation of normal pacemaking (or pacemaker arrest) is considered as a bifurcation phenomenon caused via a Hopf (or saddle node) bifurcation. On the basis of these theoretical aspects, we considered the roles of *I*_{Ca,L} and I_{Kr} in SA node pacemaking by exploring how these currents affect the stability and bifurcation structure of the model cell system.

### Roles of I_{Ca,L} in SA Node Pacemaking

*I _{Ca,L} is responsible for EP instability and pacemaker generation in primary pacemaker cells.* In the dynamic simulations of normal SA node pacemaking, the density of

*I*

_{Ca,L}during the early phase of diastolic depolarization was very small; phase 4 depolarization is not driven by the recovery from the voltage-dependent inactivation of

*I*

_{Ca,L}but by the deactivation of

*I*

_{Kr}(16, 37). Nevertheless,

*I*

_{Ca,L}appears to be indispensable for pacemaker generation in primary pacemaker cells in that it is responsible for the instability of an EP: the

*I*

_{Ca,L}-removed system, an EP of which is always stable, never exhibits pacemaker activity, whereas eliminating

*I*

_{Kr}and/or other time-dependent currents did not significantly attenuate the instability of an EP (

*V*

_{0}) (Figs. 3 and 4). Furthermore,

*I*

_{Ca,L}may be the only time-dependent current necessary for destabilizing an EP and thereby generating pacemaker activity, because spontaneous oscillations emerged in the reduced system including no time-dependent current except

*I*

_{Ca,L}when an EP was destabilized via a Hopf bifurcation (Fig. 7). The normal pacemaking of SA node cells and also the depolarization-induced abnormal automaticity of ventricular myocytes were reported to be the results of the voltage- and time-dependent interaction between

*I*

_{Ca,L}and

*I*

_{K}and thus require both of these two currents (19, 46, 62). However, our results suggest that

*I*

_{Ca,L}is, but

*I*

_{K}is not, responsible for instability of an EP as a requisite to stable pacemaking in primary pacemaker cells and that the time-dependent current required for constructing a pacemaker cell system is

*I*

_{Ca,L}alone. One should not confuse the mechanisms of “destabilizing an EP” with those of “driving phase 4 depolarization”: the early phase of diastolic depolarization in normal pacemaking is driven mainly by

*I*

_{Kr}deactivation (not by

*I*

_{Ca,L}recovery,

*I*

_{h}activation, or

*I*

_{Ca,T}activation).

*Gating kinetics of I _{Ca,L} is major determinant of instability and oscillation dynamics of SA node cells.* Instability of an EP and pacemaker generation in our SA node model appear to depend on the large time scale difference between rapid activation (

*d*

_{L}) and slow inactivation (

*f*

_{L}) of

*I*

_{Ca,L}(Figs. 5 and 6). Vinet and Roberge (62) reported that the rapid activation and slow inactivation of

*I*

_{Ca,L}underlie the development of the depolarization-induced abnormal automaticity in the model of ventricular myocytes; thus the dynamical mechanisms of normal (natural) and abnormal (ectopic) pacemaking may be essentially the same. EP instability and oscillation dynamics of our model system were substantially affected by the gating kinetics of

*I*

_{Ca,L}, especially τ

_{fL}, but not by those of other time-dependent currents, suggesting that the voltage-dependent inactivation of

*I*

_{Ca,L}is the key process responsible for regulation of the stability and dynamics of primary pacemaker cells.

### Roles of I_{Kr} in SA Node Pacemaking

*I _{K} is not necessarily required for constructing a pacemaker cell system.* At least two experimental and theoretical findings appear to indicate that

*I*

_{Kr}is required for pacemaker generation in the rabbit SA node:

*1*) spontaneous activity was abolished by blocking

*I*

_{Kr}in all of the single-cell models as well as in experiments (Fig. 2; see also Refs. 31, 37), and

*2*) deactivation of

*I*

_{Kr}is of critical importance in allowing phase 4 depolarization to occur (7, 16, 46, 75). Consistent with previous reports, our study showed that the decay of

*g*

_{Kr}(decline of

*p*

_{a}) drives the early phase of pacemaker depolarization in our model cell, supporting the

*g*

_{K}decay theory for normal pacemaking (data not shown). Nevertheless, we suggest that

*I*

_{Kr}(or

*I*

_{Ks}) is inessential to pacemaker generation, not required for constructing a pacemaker cell system, for the following reasons.

*1*) The instability of an EP (

*V*

_{0}) did not significantly attenuate even when

*I*

_{Kr}(and

*I*

_{Ks}) was removed (Figs. 3 and 4).

*2*) The

*I*

_{K}-eliminated system restored pacemaker activity when an EP (

*V*

_{0}) was destabilized by incorporating

*I*

_{b,K}, enhancing

*I*

_{K,ACh}, or getting hyperpolarizing

*I*

_{bias}(Figs. 7, 8, and 10).

*3*) Previous reports demonstrated that after automaticity was abolished by a complete block of

*I*

_{K}the rabbit SA node resumed spontaneous activity on application of hyperpolarizing

*I*

_{bias}(68, 69) and that the rabbit SA node retained spontaneous activity in the presence of an

*I*

_{Kr}blocker via electrotonic modulation from the atrium (61).

*I _{Kr} serves as an oscillation amplifier and frequency stabilizer.* Although I

_{K}is not necessary for pacemaker cell construction, the gating kinetics of the K

^{+}current appear to play important roles in regulations of the oscillation dynamics of SA node cells. Accelerating

*I*

_{K}activation or slowing (or eliminating)

*I*

_{Kr}inactivation caused marked decreases in APA and MUV (Figs. 5 and 8). Substituting a time-independent background K

^{+}current for

*I*

_{K}caused oscillations of unstable CL (and irregular dynamics) during

*g*

_{K}increases or

*I*

_{bias}applications (Figs. 8 and 10). Of the K

^{+}current formulas tested, the original

*I*

_{Kr}of delayed activation and very rapid inactivation (inward rectification) appeared to be most suitable for yielding stable oscillations of large amplitude (high conduction velocity) and stable frequency (Table 2). Thus our study suggests that

*I*

_{Kr}acts as an oscillation amplifier and a frequency stabilizer.

*I _{Kr} may contribute to high structural stability of SA node cell system.* This study also suggests that K

^{+}current kinetics affect the structural stability of SA node cells. The structural stability to

*I*

_{bias}or

*g*

_{leak}declined when

*I*

_{Kr}activation was accelerated, when

*I*

_{Kr}inactivation was slowed down or removed, and when

*I*

_{K}was replaced by a time-independent K

^{+}current (Figs. 6 and 9). Substituting a time-independent K

^{+}current for

*I*

_{K}also caused a decline in the structural stability to

*g*

_{K}changes (Fig. 8). Of the K

^{+}currents tested, the original

*I*

_{Kr}appeared to be most suitable for creating the structurally stable system with large unstable

*V*

_{0},

*I*

_{bias}, and

*g*

_{leak}regions (Fig. 9, Table 2). Thus

*I*

_{Kr}may also play a pivotal role in robust maintenance of pacemaker activity, i.e., in preventing a bifurcation to quiescence (or irregular dynamics), which may be caused by interventions or modifications such as electrotonic loads of the atrium, current leakage via myocardial injury, and changes in conductance or gating kinetics of ion channels.

### Significance of Applying Nonlinear Dynamics and Bifurcation Theory to Cardiac Electrophysiology

*General understanding and systematic descriptions of pacemaker activities and arrhythmias.* In this study, we used stability and bifurcation analyses to investigate the roles of *I*_{Ca,L} and *I*_{Kr} in SA node pacemaking. Nonlinear dynamical approaches based on the bifurcation theory have been applied to electrophysiology of excitable cells: bifurcation analyses of Hodgkin-Huxley equations for neurons, skeletal muscles, and cardiac myocytes were performed to elucidate the mechanisms underlying various diseases such as epilepsy, myotonia, and cardiac arrhythmias (3, 8, 10–12, 15, 20, 53, 55, 57). Previous studies indicated that investigating bifurcation structures of ventricular or SA node models is useful for general understanding and systematic descriptions of the mechanisms of abnormal pacemaker activities (10, 11, 22, 39, 62) as well as reentrant arrhythmias or conduction block (9, 12, 38, 63). However, the mechanisms of normal SA node pacemaking have never been explored in terms of the nonlinear dynamics and bifurcation theory. Our work shows that the mathematical approach provides a convenient way of understanding how individual currents contribute to SA node pacemaking. Further investigations of bifurcation structures of our model may also serve to clarify the mechanisms of abnormal SA node activities such as skipped-beat runs, annihilations, and early afterdepolarizations (11, 25, 39, 41, 48).

*Prediction and regulation of cell system dynamics with application to cell system engineering.* Bifurcation analyses of a model cell system may allow us to accurately predict and properly control the dynamics of real cells. One of our goals is to predict the onset of cardiac arrhythmias by exploring the bifurcation structures of a whole heart model and to find an appropriate treatment (i.e., the best way of controlling the stability and dynamics of the heart) with antiarrhythmic agents. Bifurcation theory may also be applicable to cell system engineering to develop pacemaker cells for the therapy of cardiac diseases (43). The theoretical approach would reveal the conditions and modes for oscillatory behaviors to emerge or disappear in cardiac cells, possibly leading to appropriate design of a structurally stable cell system with stable pacemaker activity. Our study suggests that *1*) pacemaker activity could be developed by expressing the L-type Ca^{2+} channel alone provided that *V*_{0} becomes unstable and that *2*) structurally stable cell systems with robust pacemaking could be established by expressing the delayed-rectifier K^{+} channel of *I*_{Kr}-type behavior.

### Limitations and Perspectives of Study

*Possible disadvantage of using a simplified model not including intracellular modulators.* Our model system was constructed as a base model to reproduce normal pacemaker activity under basal conditions without autonomic or hormonal regulations. For simplicity, activities of intracellular factors known to regulate ion channel functions and pacemaker activity (e.g., cAMP, protein kinases) were assumed to be constant at their mean standard values (their effects are considered to be lumped into the parameters of our model). One may wonder whether this simplification leads to misunderstanding of the roles of *I*_{Ca,L} or *I*_{Kr} in pacemaking, because these currents are modulated by the intracellular factors (2, 6, 13, 26, 45, 64, 70). In the nonlinear dynamical aspect, a modulating factor may substantially affect the stability and dynamics of the system if the factor rapidly and largely changes its activity during a pacemaker cycle (voltage changes). However, well-known modulators, like cAMP or protein kinases, appear not to be such fast kinetic factors but to be slow kinetic factors whose activities are nearly constant under basal conditions. We therefore believe that our simplified model is valuable and highly appropriate for exploring the essential mechanisms of normal pacemaking. These factors would have to be incorporated when investigating the mechanisms of pacemaker modulations via autonomic or hormonal mediators (4, 5, 17, 27, 73).

*Possible regional differences in pacemaker mechanisms.* The pacemaking mechanisms addressed in this study are for primary pacemaker cells in the central SA node, i.e., the dominant pacemaker mechanisms. We expect that the same mechanisms underlie the pacemaking of transitional or peripheral SA node cells. However, blocking *I*_{Ca,L} abolished spontaneous activity in central SA node tissues but not in peripheral tissues (33), suggesting a regional difference in pacemaker mechanisms. In the peripheral SA node, the fast Na^{+} channel current not present in the central SA node might contribute to EP instability and pacemaker generation. Previous reports also suggested that *I*_{h} plays a predominant role in the subsidiary pacemaker mechanisms (35, 74, 75). Whether or not the pacemaker mechanisms for the central and peripheral SA node are essentially the same remains to be settled by bifurcation analyses of peripheral SA node models as well as further experiments.

*Possible differences in pacemaker mechanisms between species.* We used the rabbit SA node model, hoping that our conclusions for the rabbit would also be applicable to other species including humans. Drouin (18) reported that the electrical behaviors of the human SA node resemble those of animal SA nodes. However, the current components in SA node cells are known to be different between species, e.g., the major component of *I*_{K} is *I*_{Ks} in the guinea pig (42) and pig (51) and *I*_{Kr} in the rabbit (28, 40, 50, 56). Thus the roles of each current in the SA node pacemaking might be different between species. Bifurcation analyses of SA node models, as well as further experiments, for other species are required to see whether the pacemaker mechanisms are essentially the same in different species.

*Pacemaker mechanisms of intact SA node possibly involving electrotonic interactions.* We investigated the mechanisms of single-cell pacemaking, which could also be applicable to pacemaking of the intact SA node tissue. However, blocking *I*_{Kr} abolished spontaneous activity in the rabbit SA node tissue separated from the atrium as well as in single cells but not in the intact tissue (61, 74), suggesting that electrotonic interaction is involved in the pacemaker mechanisms of the intact SA node. Electrotonic influences of the atrium on pacemaker activity of the SA node have been studied experimentally and theoretically (30, 38, 58, 61, 66): the hyperpolarizing loads of the atrial myocyte model caused *1*) a prolongation of CL, *2*) irregular dynamics (skipped-beat runs), and then *3*) a bifurcation to quiescence (pacemaker arrest) with increase in the gap junction conductance (66) and also yielded *4*) resumption of spontaneous activity in the *I*_{K}-blocked system (61). According to the method of Watanabe et al. (66), we explored the electrotonic effects of the atrial model on the stability and dynamics of our primary pacemaker model; the results were essentially the same as reported previously (data not shown). Nevertheless, multicellular models including central and peripheral SA node cells and also atrial myocytes are required for investigation of the roles of electrotonic interactions in the pacemaking mechanisms of the intact SA node. Modeling the SA node tissue would also be useful in exploring the mechanisms of pacemaker shifts, pacemaker synchronization (44), or abnormal rhythms such as the exit block and intranodal reentry.

*Lack of experimental evidence.* The conclusions in this study are, of course, the predictions from the model cell system; thus they must be verified and supported by experimental work using real SA node cells. Abolition of pacemaker activity by blocking *I*_{Ca,L} or *I*_{Kr}, as well as other fundamental properties of the model cell behaviors, has been verified in our previous paper (37). Restoration of spontaneous activity in the *I*_{K}-removed system by hyperpolarizing currents or electrotonic loads has also been observed experimentally (61, 68). We could not find experimental evidence that substituting a time-independent K^{+} current for *I*_{K} causes unstable CL and irregular oscillations. This should be verified by future experiments, e.g., by exploring the effects of ACh applications to enhance *I*_{K,ACh} on the stability and dynamics of SA node cells in the presence of *I*_{K} blockers (see Ref. 71). It may also be of interest to examine whether the incidence of sinus arrhythmias is higher in long-QT syndromes with a mutation in *I*_{K} channels. To the best of our knowledge, there is no systematic study on the structural stability of SA node cells or tissues. We hope that this study provides a stimulus to further experimentation on this issue.

## DISCLOSURES

This work was supported in part by Ministry for Education, Science and Culture of Japan Grant-in-Aid for Scientific Research (c)15590195 (to Y. Kurata), Kanazawa Medical University Grant for Promoted Research S2001–14 (to Y. Kurata), and Kanazawa Medical University Grant for Collaborative Research C2003–1 (to Y. Kurata and T. Shibamoto).

## Footnotes

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

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

- Copyright © 2003 by the American Physiological Society