## Abstract

A cardiac biological pacemaker (BP) has been created by suppression of the inward rectifier K^{+} current (*I*_{K1}) or overexpression of the hyperpolarization-activated current (*I*_{h}). We theoretically investigated the effects of incorporating *I*_{h}, T-type Ca^{2+} current (*I*_{Ca,T}), sustained inward current (*I*_{st}), and/or low-voltage-activated L-type Ca^{2+} channel current (*I*_{Ca,LD}) on *1*) creation of BP cells, *2*) robustness of BP activity to electrotonic loads of nonpacemaking (NP) cells, and *3*) BP cell ability to drive NP cells. We used a single-cell model for human ventricular myocytes (HVMs) and also coupled-cell models composed of BP and NP cells. Bifurcation structures of the model cells were explored during changes in conductance of the currents and gap junction. Incorporating the pacemaker currents did not yield BP activity in HVM with normal *I*_{K1} but increased the critical *I*_{K1} conductance for BP activity to emerge. Expressing *I*_{h} appeared to be most helpful in facilitating creation of BP cells via *I*_{K1} suppression. In the coupled-cell model, *I*_{st} significantly enlarged the gap conductance (*G*_{C}) region where stable BP cell pacemaking and NP cell driving occur, reducing the number of BP cells required for robust pacemaking and driving. In contrast, *I*_{h} enlarged the *G*_{C} region of pacemaking and driving only when *I*_{K1} of the NP cell was relatively low. *I*_{Ca,T} or *I*_{Ca,LD} exerted effects similar to those of *I*_{st} but caused shrinkage or irregularity of BP oscillations. These findings suggest that expressing *I*_{st} most effectively improves the structural stability of BPs to electrotonic loads and the BP ability to drive the ventricle.

- mathematical model
- bifurcation analysis
- computer simulation

a cardiac biological pacemaker (BP) has recently been created by genetic suppression of the inward rectifier K^{+} current (*I*_{K1}) in guinea pig ventricular myocytes (28) or overexpression of the hyperpolarization-activated current (*I*_{h}) in canine atrial or Purkinje myocytes (32, 36), suggesting possible development of the functional BP as a therapeutic alternative to the electronic pacemaker (8, 37).

A first step for creation of the functional BP would be engineering of single BP cells, which requires deep understanding of the BP mechanisms. By bifurcation analyses of a mathematical model for human ventricular myocytes (HVMs), we have elucidated (22) the dynamical mechanisms of BP generation in *I*_{K1}-downregulated HVMs and the roles of individual sarcolemmal currents in HVM pacemaking. We have suggested that *1*) BP activity can be developed by reducing *I*_{K1} alone in HVMs as in guinea pig ventricular myocytes (41), *2*) the instability of an equilibrium point (EP) with depolarized potentials is essentially important for BP generation, and *3*) the dynamical mechanism of ventricular pacemaking is essentially the same as that of natural sinoatrial (SA) node pacemaking as reported by Kurata et al. (21). In the previous study, however, whether BP cells can be developed from HVMs by incorporating “pacemaker currents,” which contribute to phase 4 depolarization in SA node pacemaking, was not clearly shown, with the most efficient way to create BP cells remaining unknown; the pacemaker currents include *I*_{h}, T-type Ca^{2+} channel current (*I*_{Ca,T}), sustained inward current (*I*_{st}), and low-voltage-activated L-type Ca^{2+} channel current (*I*_{Ca,LD}). Therefore, we first investigated whether incorporating these pacemaker currents can yield BP activity and how this affects BP generation during *I*_{K1} suppression in the model HVM. Our findings suggested that *I*_{h} facilitates BP generation during *I*_{K1} suppression in combination with the other pacemaker currents, although *I*_{K1} downregulation is necessary for constructing BP cells from HVMs.

Our preliminary study also revealed some drawbacks of the *I*_{K1}-downregulated HVM pacemaker, compared with the SA node pacemaker, to prevent the creation of functional BPs. They include *1*) relatively low robustness to hyperpolarizing (electrotonic) loads and *2*) low ability to drive adjacent nonpacemaker cells, which may be due to the lack of pacemaker currents. Using a coupled-cell model composed of *I*_{K1}-deleted BP cells and a nonpacemaking HVM [referred to as “nonpacemaking (NP) cell”], therefore, we further explored the effects of incorporating pacemaker currents on the structural stability of BP cells (robustness of BP activity) to electrotonic modulations and BP cell ability to drive the adjacent NP cell. Our results suggested that *I*_{Ca,T}, *I*_{st}, and *I*_{Ca,LD} can, but *I*_{h} cannot, significantly improve the structural stability of BPs to electrotonic loads of HVMs and the BP ability to drive the ventricle.

This study indicates that the theoretical approach based on nonlinear dynamics and bifurcation theory may allow us to predict the effects of current modulations on pacemaker cell dynamics accurately and find out how to control the dynamical properties of real myocytes properly. Exploring bifurcation structures of the mathematical model would provide a theoretical background for engineering of BP cells and functional BPs from native cardiomyocytes or embryonic stem (ES) cells, i.e., for gene or cell therapy of bradyarrhythmias (8, 37). This study also provides novel insights into the roles of the pacemaker currents in natural SA node pacemaking. Definitions of terms specific to nonlinear dynamics and bifurcation theory are given at the end of theory and methods (see also Refs. 21–23, 30).

## THEORY AND METHODS

### Base Mathematical Model for Human Ventricular Myocytes

We used our HVM model (22) as a modified version of the Priebe-Beuckelmann model (33). More elaborate HVM models were recently developed by Iyer et al. (12) and ten Tusscher et al. (42), which are probably superior to the original or modified Priebe-Beuckelmann model in reproducing experimental data. Nevertheless, these models are much more complex or include vector field functions that are not continuous or smooth, being less suitable for bifurcation analyses (for more details, see Ref. 22). We have therefore chosen to use our model, which is more suitable for bifurcation analyses.

The standard model for the normal activity of single HVMs is described as a nonlinear dynamical system of 15 first-order ordinary differential equations. The membrane current system includes the L-type Ca^{2+} channel current (*I*_{Ca,L}), rapid and slow components of delayed rectifier K^{+} currents (denoted *I*_{Kr} and *I*_{Ks}, respectively), 4-aminopyridine-sensitive transient outward current (*I*_{to}), Na^{+} channel current (*I*_{Na}), *I*_{K1}, background Na^{+} (*I*_{Na,b}) and Ca^{2+} (*I*_{Ca,b}) currents, Na^{+}-K^{+} pump current (*I*_{NaK}), Na^{+}/Ca^{2+} exchanger current (*I*_{NaCa}), and Ca^{2+} pump current (*I*_{pCa}). Time-dependent changes in the membrane potential (*V*) are described by the equation (1) where *I*_{stim} represents the stimulus current (in pA/pF), being set equal to zero for simulations of BP activity. Details on modifications and expressions of the HVM model are described in our previous article (22).

Our full system includes material balance expressions to define the temporal variation in intracellular Ca^{2+}, Na^{+}, and K^{+} concentrations ([Ca^{2+}]_{i}, [Na^{+}]_{i}, and [K^{+}]_{i}), whereas extracellular Ca^{2+}, Na^{+}, and K^{+} concentrations ([Ca^{2+}]_{o}, [Na^{+}]_{o}, and [K^{+}]_{o}) were fixed at 2, 140 and 5.4 mM, respectively. As pointed out by Hund et al. (11) and Krogh-Madsen et al. (18), the second-generation model incorporating ion concentration changes exhibits degeneracy (i.e., nonuniqueness of steady-state solutions), not suitable for bifurcation analysis to be applicable to isolated equilibria. One of the ways to remove degeneracy and thus allow bifurcation analyses of isolated equilibria is to make some ionic concentrations fixed (18). Variations in [K^{+}]_{i} during changes in bifurcation parameters were relatively small, not significantly altering the model cell behaviors. For stability and bifurcation analyses, therefore, [K^{+}]_{i} was fixed at 140 mM.

### Incorporation of Pacemaker Currents

To investigate how pacemaker currents contributing to phase 4 depolarization in SA node pacemaking affect bifurcation structures of HVM model cells, we incorporated the formulas of *I*_{h}, *I*_{Ca,T}, and *I*_{st}, which were used for our SA node model (20) or experimentally determined for human cardiomyocytes. The low voltage-activated L-type Ca^{2+} channel (D-LTCC) was recently reported to activate at pacemaker potential regions, i.e., at potentials 10–20 mV more negative than the high-voltage-activated Ca^{2+} channel, contributing to mouse SA node pacemaking (16, 27, 31, 51). Therefore, we also tested the effect of the current mediated by D-LTCC (denoted as *I*_{Ca,LD}). The formulas for individual pacemaker currents are given in appendix b. The origins of the expressions and conductance values tested for each current, as well as those chosen for the original SA node models as mean experimental values, are summarized in Table 1 (for details, see below).

For comparison of the voltage-dependent gating behaviors of the currents tested, steady-state probabilities and time constants of gating variables, as well as the steady-state (window) currents calculated by the steady-state equations for the gating variables, are shown in Fig. 1. To determine the conductance range to be tested for individual currents, we first examined their effects on stability and dynamics of a single BP cell (Fig. 2). Higher conductance of *I*_{h} caused a saddle-node bifurcation and cessation of BP activity with drastic increases in [Na^{+}]_{i} and [Ca^{2+}]_{i}; higher conductance of *I*_{Ca,T}, *I*_{st}, or *I*_{Ca,LD} caused shrinkage of BP oscillations with decrease in upstroke velocity and then cessation of BP activity via Hopf bifurcations. The maximum conductance value to be tested was limited for each current so as not to significantly shrink BP oscillations or not to cause Na^{+}- or Ca^{2+}-overload conditions (for more details, see below).

#### I_{h}.

The expressions for *I*_{h} were adopted from the rabbit SA node model (20, 49), as well as from experimental studies for human *I*_{h} (3, 29). Although all the *I*_{h} formulas exerted qualitatively the same effects, the effect of the rabbit SA node *I*_{h} was most dramatic, which is probably due to the most positive activation threshold. Thus we show data only for the rabbit SA node *I*_{h}. The maximum conductance (*g*) value for *I*_{h} (*g*_{h}), 0.375 nS/pF in the original SA node model (20, 49), was limited to 0.12 nS/pF (for the effects on BP generation) or 0.5 nS/pF (for the effects on pacemaking and driving ability of BP cells); higher *g*_{h} caused inaccurate calculations of steady states with drastic increases of [Na^{+}]_{i} and [Ca^{2+}]_{i} to unreasonable values.

#### I_{Ca,T}.

The expressions for *I*_{Ca,T} were adopted from the rabbit SA node (5, 20) and ventricular (34) models. The effects of these *I*_{Ca,T} currents were qualitatively the same, with the effect of the SA node *I*_{Ca,T} being more remarkable; thus the data for the rabbit SA node *I*_{Ca,T} are shown. Relatively high *g*_{Ca,T} was required to affect dynamical properties of the model cell, probably because the window region (i.e., overlap of activation and inactivation curves) of the model *I*_{Ca,T} is relatively small (see Fig. 1). Therefore, the *g*_{Ca,T} value was increased up to 3 nS/pF, although it was 0.458 nS/pF in the original SA node model (5, 20). Further increases in *g*_{Ca,T} caused shrinkage of BP oscillations and stabilization of an EP via a Hopf bifurcation (see Fig. 2).

#### I_{st}.

The formulas for *I*_{st} were adopted from the rabbit SA node model (20) and the experimental report for the rat SA node (40). The effects of the two *I*_{st} currents were essentially the same; thus only data for the rabbit SA node *I*_{st} are shown. The *g*_{st} value, 0.015 nS/pF in the original SA node model (20), was limited to 0.02 nS/pF; further increases in *g*_{st} yielded extremely large window currents leading to Na^{+}-overload conditions, because the window region of the model *I*_{st} is relatively large (see Fig. 1).

#### I_{Ca,LD}.

The kinetics of *I*_{Ca,LD} was formulated on the basis of experimental data from mouse SA node cells (27, 31, 51) and human embryonic kidney cells expressing D-LTCC (16). These previous studies suggested that *1*) the voltage dependence of the activation and inactivation of *I*_{Ca,LD} is 10–20 mV more negative than that of the high-voltage-activated Ca^{2+} channel current (*I*_{Ca,LC}), *2*) the activation of *I*_{Ca,LD} is faster than that of *I*_{Ca,LC}, and 3) the inactivation of *I*_{Ca,LD} is 1.5–2 times slower than that of *I*_{Ca,LC}. According to these data, therefore, we formulated the kinetics of *I*_{Ca,LD} as provided in appendix b (*Eqs. A38*–*A42*). Because the effects on model cell dynamics of the model *I*_{Ca,LD} were relatively small, the *g*_{Ca,LD} value was increased up to 1 nS/pF. Further increases in *g*_{Ca,LD} caused Ca^{2+}-overload conditions as well as shrinkage of BP oscillations (see Fig. 2).

### Formulation of Coupled-Cell Model

To investigate the electrotonic influences of adjacent HVMs on stability and dynamics of BP cells as well as BP cell ability to drive adjacent HVMs, we used a coupled-cell model (Fig. 3), in which *I*_{K1}-deleted BP cells were connected to a NP cell with normal or reduced *I*_{K1} via the gap junction conductance (*G*_{C}) of 0–20 nS. One NP cell was connected with one to seven BP cells; more than one BP cells were assumed to be well-coupled enough to synchronize completely and act as a cluster of isopotential cells (or assumed to connect with the NP cell via an identical *G*_{C}). Pacemaker currents were incorporated into the BP cells only (not into the NP cell).

Time-dependent changes in the membrane potentials of the BP (*V*_{BP}) and NP (*V*_{NP}) cells were calculated by the equations (2) (3) where *C*_{BP} and *C*_{NP} represent the membrane capacitance of the BP and NP cells, respectively. *I*_{total(BP)} and *I*_{total(NP)} are the sum of sarcolemmal ionic currents in the BP and NP cells, with *I*_{GJ} denoting the gap junction current. *I*_{stim} is the stimulus current applied to the BP cells, usually set equal to zero. The units for *V*, *C*, *I*_{total} (*I*_{stim}), and *I*_{GJ} are millivolts, picofarads, picoamperes per picofarad, and picoamperes, respectively. Time-dependent changes in all the other state variables of the BP and NP cells were computed independently, as described for the single HVM in our previous article (22).

### Numerical Methods for Dynamic Simulations and Bifurcation Analyses

Dynamic behaviors of the model cells were determined by solving a system of nonlinear ordinary differential equations numerically. Numerical integration as well as bifurcation analyses were performed on Workstation CELSIUS X630 (Fujitsu, Tokyo, Japan) with MATLAB 7 (MathWorks, Natick, MA). We used the numerical algorithms available as MATLAB ODE solvers: *1*) a fourth-order adaptive-step Runge-Kutta algorithm that includes an automatic step-size adjustment based on an error estimate (*ode45*) (6) and *2*) a variable time-step numerical differentiation approach selected for its suitability to stiff systems (*ode15s*) (39). Both solvers usually yielded nearly identical results. The latter, much more efficient than the former, was usually used; the former as the best function for most problems was only sometimes used to confirm the accuracy of calculations. The maximum relative error tolerance for the integration methods was set to 1 × 10^{−6}.

### Stability and Bifurcation Analyses

We examined how the stability and dynamics of the model cells 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 of *I*_{K1} (*g*_{K1}) and pacemaker currents (*g*_{h}, *g*_{Ca,T}, *g*_{st}, *g*_{Ca,LD}) as well as *G*_{C}; *g*_{K1} is expressed as a normalized value, i.e., in ratio to the control value of 3.9 nS/pF, unless otherwise stated. The methods for bifurcation analyses of the single BP or NP cell are described in our previous article (22).

Although the methods of bifurcation analyses for the coupled-cell system are essentially the same as for the single cell, the number of state variables increases twice: in the coupled-cell system with fixed [K^{+}]_{i}, 28 state variables define a 28-dimensional state point in the 28-dimensional state space of the system. We calculated EPs, i.e., steady-state values of the variables, and periodic orbits in the state space of the coupled-cell system. Steady-state values of the variables were calculated by *Eqs. A53*–*A58* in appendix c. Asymptotic stability of the EP was also determined by computing 28 eigenvalues of a 28 × 28 Jacobian matrix derived from the linearization of the nonlinear system around the EP (for details, see Ref. 45). Periodic orbits were located with the MATLAB ODE solvers. When spontaneous oscillation (BP activity) appeared, the action potential amplitude (APA) as a voltage difference between the maximum diastolic potential (*V*_{min}) and peak overshoot potential (*V*_{max}) 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. Potential extrema (*V*_{min}, *V*_{max}) and CL of the steady-state oscillation were plotted against bifurcation parameters. When periodic behavior was irregular or unstable, model dynamics were computed for 60 s; all potential extrema and CL values were then plotted.

We constructed one-parameter bifurcation diagrams for *G*_{C} and two-parameter bifurcation diagrams for *G*_{C} and pacemaker current conductance. For construction of one-parameter bifurcation diagrams, the membrane potential at EPs (steady-state branches) and local potential extrema (*V*_{min}, *V*_{max}) of periodic orbits (periodic branches) were determined and plotted as a function of *G*_{C}, which was systematically changed while keeping all other parameters at their standard values. The saddle-node bifurcation point at which two EPs coalesce and disappear with emergence of robust BP activity was determined; the Hopf bifurcation point at which an unstable EP is stabilized with cessation of BP activity was also located. For construction of two-parameter bifurcation diagrams, the critical *G*_{C} values at bifurcation points were determined as functions of the secondary parameters such as the pacemaker current conductance, the number of BP cells, and NP cell *g*_{K1}; the primary parameter (*G*_{C}) was systematically changed, with the secondary parameter fixed at various different values. The path of Hopf and saddle-node bifurcation points was traced in the parameter plane, i.e., bifurcation values for the primary parameter were plotted as a function of the secondary parameter.

We also evaluated the “structural stability” of BP cells, which is defined as the robustness of BP activity to various interventions or modifications that may cause a bifurcation to quiescence or irregular dynamics (21, 22). In this study, we tested the structural stability of BP cells to electrotonic (hyperpolarizing) loads of an adjacent NP cell, using the coupled-cell system, which is also useful in exploring BP cell ability to drive the adjacent NP cell. The method of evaluating the structural stability to electrotonic loads is essentially the same as that for constant bias currents as described in our previous articles (21, 22). In the unstable *G*_{C} region where the system has no stable EP, the system usually exhibited robust pacemaking and driving without annihilation (9, 24). When the system removes to the stable *G*_{C} region with a stable EP, it will come to a rest at the stable EP via gradual decline of limit cycles, annihilation, or irregular dynamics, as reported by Guevara and Jongsma (9). In the coupled-cell system, the larger the *G*_{C} region over which all EPs are unstable (there is no stable EP) and thus BP activity appears is, the more structurally stable the BP system is (for more detail, see Ref. 22).

### Definitions of Terms Specific to Nonlinear Dynamics and Bifurcation Theory

#### Equilibrium point.

The equilibrium point (EP) is a time-independent steady-state point at which the vector field vanishes in the state space of a dynamical system, constructing the steady-state branch in one-parameter bifurcation diagrams. This state point corresponds to the zero-current crossing in the steady-state current-voltage (*I-V*) curve, i.e., a quiescent (resting) state of a cell if it is stable.

#### Periodic orbit.

A periodic orbit is a closed trajectory in the state space of a system, constructing the periodic branch in one-parameter bifurcation diagrams.

#### Limit cycle.

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.

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 a quiescent state. Bifurcation phenomena we can see in cardiac myocytes include generation or cessation of pacemaker activity and occurrence of irregular dynamics such as skipped-beat runs and early afterdepolarizations.

#### Saddle-node bifurcation.

Saddle-node bifurcation is a bifurcation at which two EPs (steady-state branches) or two periodic solutions (periodic branches) emerge or disappear. The saddle-node bifurcation of EPs occurs when one of the eigenvalues of a Jacobian matrix for the EP is zero.

#### Hopf bifurcation.

The Hopf bifurcation is a bifurcation at which the stability of an EP reverses with emergence or disappearance of a limit cycle, occurring when eigenvalues of a Jacobian matrix have a single complex conjugate pair and its real part reverses the sign through zero.

## RESULTS

### Effects of Pacemaker Currents on BP Generation in Single HVM Model

Our preceding study (22) suggests that *I*_{K1} downregulation, ensuring both the occurrence of a saddle-node bifurcation and instability of the EP with depolarized potentials, yields spontaneous oscillations (BP activity) in HVMs without incorporating any pacemaker current. However, a saddle-node bifurcation to induce BP activity may also be yielded by incorporating pacemaker currents such as *I*_{h}, *I*_{Ca,T}, *I*_{st}, and *I*_{Ca,LD}, which are abundant in SA node primary pacemaker cells but absent or very small in native HVMs. In the preliminary study, we found that these pacemaker currents could not yield BP generation in the model HVM with normal *I*_{K1} (*g*_{K1} = 1): BP activity did not appear when *g*_{h}, *g*_{Ca,T}, *g*_{st}, and *g*_{Ca,LD} were increased up to the limited values of 0.12, 3, 0.02 and 1 nS/pF, respectively (see theory and methods) or even when *g*_{Ca,T}, *g*_{st}, or *g*_{Ca,LD} was further increased to the critical value above which an EP at depolarized potentials would be stable and thus BP oscillations could not occur (see Fig. 2). Nevertheless, incorporating the pacemaker currents is expected to accelerate BP generation during *I*_{K1} suppression by increasing the critical *g*_{K1} value at which a saddle-node bifurcation occurs and thus BP oscillations emerge. We therefore investigated how these currents affect the saddle-node bifurcation and BP generation during *I*_{K1} suppression.

#### Effects of incorporating each pacemaker current on BP generation during I_{K1} suppression.

As clearly shown in Fig. 4, the saddle-node bifurcation point below which BP oscillations occur shifted toward higher *g*_{K1} values with increasing *g*_{h}, suggesting that expression of *I*_{h} accelerates BP generation during *I*_{K1} suppression by increasing the critical *g*_{K1} value to induce BP activity. To further examine the effects of incorporating the pacemaker currents on BP generation in the *I*_{K1}-reduced HVM, we determined the critical *g*_{K1} value for BP generation (i.e., *g*_{K1} value at a saddle-node bifurcation) as functions of *g*_{h}, *g*_{Ca,T}, *g*_{st}, or *g*_{Ca,LD}, as well as exploring the influences of the pacemaker currents on the steady-state *I-V* relation at *g*_{K1} = 1 (Fig. 5). The inward pacemaker currents counteracted the outward component of *I*_{K1}, inwardly shifting the steady-state *I-V* curve. The critical *g*_{K1} value increased with increasing *g*_{h}; however, it did not reach the control *g*_{K1} value but reached the maximum of 0.287 (28.7%) at *g*_{h} = 0.1039 nS/pF. *I*_{Ca,T} at higher densities (*g*_{Ca,T} ≤ 3 nS/pF) also increased the critical *g*_{K1}, with its increase limited to <0.266 (26.6%). The effects of *I*_{st} and *I*_{Ca,LD} were relatively small within the limited conductance ranges of ≤0.02 and ≤1 nS/pF, respectively.

#### Effects of coexpression of I_{h} and other pacemaker currents on BP generation.

The critical *g*_{K1} value did not reach the control *g*_{K1} value when the conductance of individual currents was increased up to the limited values, or even when *g*_{Ca,T}, *g*_{st}, or *g*_{Ca,LD} was increased to more than the Hopf bifurcation point at which BP oscillation would disappear (see Fig. 2). However, the critical *g*_{K1} may further be increased by coexpression of *I*_{h} and the other pacemaker currents that have a different voltage range of activation. As shown in Fig. 6, therefore, we also examined the effects of coexpressing *I*_{h} and *I*_{Ca,T}, *I*_{st}, or *I*_{Ca,LD} on BP generation during *I*_{K1} suppression. *I*_{Ca,T}, *I*_{st}, or *I*_{Ca,LD} coexpressed with *I*_{h} further increased the critical *g*_{K1} value, i.e., enhanced the accelerating effect of *I*_{h}. Nevertheless, the increase in the critical *g*_{K1} was limited even when *I*_{h} and the other pacemaker currents were coexpressed. Drastic increases of [Na^{+}]_{i} and [Ca^{2+}]_{i} at the bifurcation points were caused by coexpression of *I*_{h} and the other currents at higher conductance.

### Structural Stability and Driving Ability of HVM Pacemaker

One of the drawbacks of the *I*_{K1}-downregulated HVM pacemaker may be the relatively low structural stability to hyperpolarizing loads such as electrotonic modulations of surrounding nonpacemaker cells, because it lacks the pacemaker currents that are abundant in SA node primary pacemaker cells and contribute to pacemaker depolarization. We therefore examined the structural stability of BP cells to electrotonic loads of the adjacent NP cell, using a coupled-cell model in which BP cells (*g*_{K1} = 0) were connected to a NP cell (*g*_{K1} = 0.2–1) via the *G*_{C} ranging from 0 to 20 nS. The coupled-cell system is also useful in exploring BP cell ability to drive adjacent NP cells.

#### G_{C}-dependent behaviors of the coupled-cell system.

In the system with one BP cell and one normal HVM (*g*_{K1} = 1), an EP was stabilized via a Hopf bifurcation at *G*_{C} = 0.689 nS, with the BP cell not exhibiting pacemaker activity at *G*_{C} values higher than the bifurcation value. Thus the coupling of one BP cell to one normal HVM produced either *1*) BP cell pacemaking at slower rates without driving the NP cell for lower *G*_{C} values or *2*) inhibition of BP cell pacemaking for higher *G*_{C} values. We could not find a *G*_{C} value for which one BP cell could exhibit pacemaking itself and drive the normal HVM but found that three or more BP cells, assumed to be well coupled to form a cluster, were required for driving the normal HVM.

Figure 7 shows the simulated membrane potentials and *I*_{GJ} at various *G*_{C} values for the coupled-cell system composed of five BP cells and one normal HVM. With increasing *G*_{C}, the BP rate decreased, with the oscillation amplitude getting larger. The BP cells could induce only subthreshold responses of the HVM at *G*_{C} = 2.0 nS, whereas they induced driven action potentials at *G*_{C} = 2.5–3.5 nS. Further increase in *G*_{C} to 5.0 nS abolished BP activity, which was probably due to hyperpolarizing loads of the NP cell as indicated by positive *I*_{GJ}.

#### Influences of size factor and NP cell g_{K1} on G_{C}-dependent bifurcation structures.

The structural stability and driving ability of pacemaker systems are known to depend on the size factor (i.e., the number of BP cells relative to that of NP cells) and *g*_{K1} of adjacent NP cells (13, 14, 19, 47), which are practically important for engineering of functional BPs in vivo. Therefore, we more minutely explored the effects of the size factor and NP cell *g*_{K1} on *G*_{C}-dependent bifurcation structures of coupled-cell systems.

Stability and dynamics during *G*_{C} increases of the coupled-cell system composed of one to seven BP cells and one normal HVM are shown in Fig. 8*A*. EPs and potential extrema were determined for both the BP and NP cells, plotted as functions of *G*_{C} values. Hopf and saddle-node bifurcation points were also determined and located in the bifurcation diagrams. In the systems with one to five BP cells, hyperpolarizing loads of the NP cell caused *1*) prolongation of CL with an increase in oscillation amplitude, *2*) irregular dynamics, and then *3*) a bifurcation to quiescence (Hopf bifurcation), with increasing *G*_{C}. These behaviors are essentially the same as those reported for interactions of the SA node cell and the atrial or ventricular myocyte (13, 44, 47, 48). The critical *G*_{C} value at the Hopf bifurcation to cause stabilization of EPs and cessation of BP activity increased with increasing number of BP cells. One or two BP cells could not drive the normal HVM. With three to five BP cells, the increase in *G*_{C} yielded a BP cell ability to drive the normal HVM; however, further increases in *G*_{C} led to the cessation of pacemaking and driving via the Hopf bifurcation. The *G*_{C} region of driving the normal HVM was very limited with three to five BP cells, whereas it dramatically enlarged when the number of BP cells increased to more than five. With six or more BP cells, a saddle-node bifurcation at which two EPs (EP2, EP3) disappear occurred to ensure the robust pacemaking and driving without annihilation (see Fig. 8*A*, *far right*, for 7 BP cells).

The relatively low ability to drive the normal HVM and low structural stability to electrotonic loads of the BP cell appeared to be due to the relatively high *I*_{K1} density in the normal HVM (14). Therefore, we also explored stability and dynamics during *G*_{C} increases of a coupled-cell system consisting of one BP cell and one NP cell with reduced *I*_{K1}, as shown in Fig. 8*B* for the normalized NP cell *g*_{K1} of 0.35, 0.30, and 0.25. With the reduced g_{K1} of 0.30 (1.17 nS/pF) or less, increasing *G*_{C} did not cause a Hopf bifurcation to quiescence but yielded a saddle-node bifurcation to ensure robust pacemaking and driving.

### Effects of Pacemaker Currents on Structural Stability and Driving Ability of HVM Pacemaker

Using the coupled-cell model, we examined how incorporating the pacemaker currents affects the structural stability of BP cells to electrotonic loads of the adjacent NP cell, as well as BP cell ability to drive the NP cell. Figure 8 suggests that a Hopf bifurcation point is just at or close to the point at which BP activity is abolished by hyperpolarizing loads of the NP cell, whereas a saddle-node bifurcation point is the point above which robust pacemaking and driving occur. To assess the effects of the pacemaker currents on the structural stability and driving ability of BP cells, therefore, we focused on the shift in the Hopf bifurcation point, as well as emergence of a saddle-node bifurcation.

#### Effects of incorporating pacemaker currents on G_{C}-dependent bifurcation structures.

Bifurcation structures during *G*_{C} increases were first determined at various conductances of each pacemaker current. We constructed two-parameter bifurcation diagrams in which the critical *G*_{C} value at Hopf or saddle-node bifurcation points was plotted as functions of the maximum current conductance. Figure 9 shows the effects of incorporating *I*_{h}, *I*_{Ca,T}, *I*_{st}, or *I*_{Ca,LD} on the bifurcation structure during *G*_{C} increases of the coupled-cell system in which three BP cells were connected to one normal HVM or one BP cell was connected to one NP cell with reduced *I*_{K1} (*g*_{K1} = 0.35). Stability and dynamics during *G*_{C} increases of the coupled-cell systems incorporating each pacemaker current are illustrated in Fig. 10. The results suggest that individual pacemaker currents affect the *G*_{C}-dependent bifurcation structure of the system in different ways.

##### I_{H} EFFECT.

In the coupled-cell system with three BP cells and one normal HVM, incorporation of *I*_{h} (∼0.5 nS/pF) did not significantly change the critical *G*_{C} value at Hopf bifurcations where an EP is stabilized and BP activity is abolished, not enlarging the region of pacemaking labeled “P^{+}D^{±}” (Fig. 9, *top*). As shown in Fig. 10*A*, *I*_{h} (0.2 nS/pF) did not enlarge the *G*_{C} region where all EPs are unstable and BP activity appears. When *I*_{K1} of the NP cell was relatively small (*g*_{K1} = 0.35), however, *I*_{h} enlarged the pacemaking and driving region labeled P^{+}D^{+} by producing a saddle-node bifurcation (Fig. 9, *bottom*; Fig. 10*B*). Overexpression of *I*_{h} (*g*_{h} > 0.2043 nS/pF) caused cessation of BP activity at higher *G*_{C} values, with a stable resting potential dramatically hyperpolarized probably via drastic increases in [Na^{+}]_{i} and [Ca^{2+}]_{i}.

##### I_{CA,T} EFFECT.

Similar to *I*_{h}, *I*_{Ca,T} (∼3 nS/pF) did not significantly change the critical *G*_{C} value at Hopf bifurcations, not enlarging the P^{+}D^{±} region, in the coupled-cell system with three BP cells and one normal HVM (Fig. 9, *top*). However, higher densities of *I*_{Ca,T} yielded a saddle-node bifurcation and thus enlarged the P^{+}D^{+} region; *I*_{Ca,T} at 2 nS/pF enlarged the *G*_{C} region where all EPs are unstable and BP activity appears (Fig. 10*A*). When the NP cell *I*_{K1} was relatively small, *I*_{Ca,T} dramatically enlarged the P^{+}D^{+} region by producing a saddle-node bifurcation (Fig. 9, *bottom*; Fig. 10*B*). With relatively high conductance of *I*_{Ca,T}, the amplitudes of BP oscillations and driven action potentials of the NP cell were relatively small.

##### I_{ST} EFFECT.

In contrast to *I*_{h} or *I*_{Ca,T}, incorporation of *I*_{st} (∼0.02 nS/pF) dramatically shifted Hopf bifurcation points toward higher *G*_{C} values and further produced a saddle-node bifurcation to ensure robust pacemaking and driving, i.e., enlarged the P^{+}D^{+} region (Fig. 9). *I*_{st} at 0.01 nS/pF dramatically enlarged the *G*_{C} region where BP oscillations occur, with CL of BP oscillations being relatively stable during increases in *G*_{C} (Fig. 10).

##### I_{CA,LD} EFFECT.

The effects of *I*_{Ca,LD} were similar to those of *I*_{st}, although much larger conductance was required for *I*_{Ca,LD} to exert the effects as dramatically as *I*_{st}. Incorporation of *I*_{Ca,LD} (∼1 nS/pF) increased the Hopf bifurcation value and produced a saddle-node bifurcation to enlarge the P^{+}D^{+} region (Fig. 9). As shown in Fig. 10, *I*_{Ca,LD} (0.5 nS/pF) enlarged the *G*_{C} region where BP oscillations occur. However, *I*_{Ca,LD} at high densities of 0.5–1 nS/pF caused unstable CL and irregular dynamics during increases in *G*_{C}.

#### Influences of altering pacemaker current kinetics on their modulation of bifurcation structures.

Individual pacemaker currents affected the *G*_{C}-dependent bifurcation structure in different ways, as well as with different potencies. The effects of *I*_{h} were completely different from those of the other pacemaker currents, which probably resulted from the difference in the voltage range of activation. Furthermore, the effects of *I*_{Ca,T}, *I*_{st}, and *I*_{Ca,LD} were apparently different in the system with three BP cells and one normal HVM: the locus of Hopf bifurcation points at lower *G*_{C} was separated from that of saddle-node bifurcation points with *I*_{Ca,T}, but not with *I*_{st} or *I*_{Ca,LD}. The differences in their effects may result from the large differences in their inactivation kinetics, whereas the different potencies are probably due to the different degrees of the overlap of activation and inactivation curves (see Fig. 1). Thus we examined the effects of incorporating the modified *I*_{Ca,T} or *I*_{st} with slowed or accelerated inactivation on *G*_{C}-dependent bifurcation structures of the system with three BP cells and one normal HVM (Fig. 11). Changing the inactivation kinetics altered the loci of Hopf bifurcations: the modified *I*_{Ca,T} with the inactivation time constant (τ_{fT}) 10 times higher than the control yielded a bifurcation diagram similar to that for the original *I*_{st}; in contrast, the modified *I*_{st} with the inactivation time constant (τ_{qi}) decreased to 0.01 of the control yielded a locus of Hopf bifurcations separated from that of saddle-node bifurcations, like the original *I*_{Ca,T}. Thus the different effects of the pacemaker currents are at least in part ascribable to the difference in their inactivation kinetics, which appeared to affect the Hopf bifurcation, i.e., stability of EPs, but not the saddle-node bifurcation.

#### Effects of pacemaker currents on size factor or NP cell g_{K1}-dependent bifurcation structures.

To further examine how the pacemaker currents affect the pacemaking and driving of BP cells, bifurcation structures during *G*_{C} increase were determined with change in the size factor (the number of BP cells) or NP cell *g*_{K1} for the modified systems with BP cells incorporating *I*_{h}, *I*_{Ca,T}, *I*_{st}, or *I*_{Ca,LD}, as well as for the standard system not including a pacemaker current (control). We constructed two-parameter bifurcation diagrams in which the critical *G*_{C} value at saddle-node or Hopf bifurcation points was plotted as functions of the size factor or normalized NP cell *g*_{K1}; one to six BP cells were connected to one normal HVM (Fig. 12) or one BP cell was connected to one NP cell with reduced *g*_{K1} of 0.2–0.7 (Fig. 13).

Incorporating *I*_{h} (0.1–0.2 nS/pF) exerted relatively small effects on the size factor-dependent bifurcation structure of the coupled-cell system with a normal HVM, enlarging the P^{+}D^{+} area only in the region of relatively large size factors and high *G*_{C} values (Fig. 12). Decreasing *g*_{K1} of the NP cell, *I*_{h} enlarged the *P*^{+}*D*^{+} area by producing a saddle-node bifurcation at relatively low *g*_{K1} values (Fig. 13). In contrast to *I*_{h}, *I*_{st} (0.01–0.02 nS/pF) dramatically enlarged the P^{+}D^{+} area by facilitating the emergence of saddle-node bifurcations at smaller size factors and higher *g*_{K1} values. The effects of *I*_{Ca,T} (1–2 nS/pF) or *I*_{Ca,LD} (0.5–1 nS/pF) were similar to those of *I*_{st}, although much larger conductance was required for *I*_{Ca,T} or *I*_{Ca,LD} to exert effects as dramatically as *I*_{st}.

## DISCUSSION

### Possible Roles of Pacemaker Currents in Creation of BP Cells

We first investigated whether incorporating pacemaker currents exerts beneficial effects on creation of BP cells by constructing two-parameter bifurcation diagrams in which the critical *g*_{K1} value (saddle-node bifurcation point) to produce BP activity was plotted as functions of the maximum current conductance.

#### Pacemaker currents facilitate BP generation during I_{K1} suppression.

Our results suggest that the pacemaker currents, especially *I*_{h}, can facilitate BP generation during *I*_{K1} suppression, although their effects were limited (Fig. 5). The mechanism underlying the facilitation of BP generation appears to be relatively simple: the inward pacemaker currents counteract the outward component of *I*_{K1}, inwardly shift the steady-state *I-V* curve, and thus lead to a saddle-node bifurcation at a higher *g*_{K1} value (see Fig. 5, *top*). Nevertheless, the critical *g*_{K1} value never reached the control value of 1, suggesting that any pacemaker currents cannot yield BP activity in the HVM with normal *I*_{K1}.

#### I_{K1} downregulation is a requisite for BP generation in HVM.

The limited effects of the pacemaker currents appear to be due to the induction of intracellular Na^{+} and Ca^{2+} overloads: as shown in Fig. 6, overexpressions of the pacemaker currents cause increases in [Na^{+}]_{i} and [Ca^{2+}]_{i}, which enhance the outward Na^{+}-K^{+} pump current, attenuate the inward background Na^{+} and Ca^{2+} currents (by decreasing the driving force for Na^{+} and Ca^{2+}), and thus counteract the inward shift of *I-V* curves that leads to a saddle-node bifurcation. When [Na^{+}]_{i}, which is a state variable in the standard system, was fixed at the control value of 6.14 mM (for *g*_{K1} = 1) to prevent the pacemaker currents causing Na^{+} and Ca^{2+} overloads, increasing the pacemaker currents led to BP generation via a saddle-node bifurcation without inhibition of *I*_{K1}, albeit with the requirement of very high conductance (e.g., *g*_{h} > 0.996 nS/pF) (data not shown). Thus Na^{+} and Ca^{2+} overloads would prevent a bifurcation to BP generation during enhancements of the pacemaker currents. In conclusion, incorporating *I*_{h} and/or the other pacemaker currents would not yield BP generation in HVMs with normal *I*_{K1}, and thus *I*_{K1} downregulation would be necessary for constructing a BP cell system from HVMs. Inconsistent with our result, Viswanathan et al. (46) reported that the rabbit SA node *I*_{h} could induce BP activity in the Luo-Rudy guinea pig ventricular model without reducing *I*_{K1}. We are not sure about the reason for this inconsistency. Further studies are required to clarify whether the effects of *I*_{h} and other pacemaker currents depend on species (human vs. guinea pig) or on the property of a model itself as well as whether expressing *I*_{h} can actually induce BP activity in real NP myocytes of the ventricle.

Overexpression of *I*_{h} has been reported to induce BP activity (or increase BP rates) in canine atrial or Purkinje myocytes (32, 36). Figure 5 indicates that increasing *I*_{h} alone can induce BP activity via a saddle-node bifurcation if *g*_{K1} is smaller than the control value of 3.9 nS/pF, e.g., when a *g*_{K1} value is 25% of the control value (0.975 nS/pF), increasing *g*_{h} to 0.0403 nS/pF or more yields BP generation. Koumi et al. (17) reported that a slope conductance at a reversal (K^{+} equilibrium) potential of the human atrial *I*_{K1} (0.134 nS/pF) is only 18.3% of that of the human ventricular *I*_{K1} (0.732 nS/pF). As suggested by previous reports (14, 32, 36), expression of *I*_{h} alone would induce pacemaker activity in myocytes with relatively low density of *I*_{K1}, such as atrial and Purkinje myocytes.

### Possible Roles of Pacemaker Currents in Robust Pacemaking and Driving of Functional BP

Using the coupled-cell model, we examined whether incorporation of pacemaker currents exerts beneficial effects on BP functions, i.e., whether it improves the structural stability of BPs to electrotonic loads of adjacent NP cells and the BP ability to drive the ventricle. The results suggest that individual pacemaker currents affect the BP functions in different ways.

#### I_{h} effect is relatively small and NP cell g_{K1} dependent.

The effects of *I*_{h} on *G*_{C}-dependent bifurcation structures were quite different from those of the other pacemaker currents, which is probably due to the difference in the voltage range of activation. *I*_{h} did not significantly improve the structural stability of BP cells to electrotonic loads of the normal HVM or the BP cell ability to drive the normal HVM (Figs. 9 and 10*A*). When NP cell *I*_{K1} was relatively small, however, *I*_{h} improved the structural stability and driving ability of BP cells by producing a saddle-node bifurcation (Figs. 9 and 10*B*). Moreover, the effects of I_{h} were shown to depend on I_{K1} density of the NP cell, as well as the size factor (Figs. 12 and 13). These results suggest that *I*_{h} can support the robust pacemaking and driving of BPs only when *I*_{K1} density in adjacent NP cells is relatively small. *I*_{K1} density is much larger in the human ventricle than in the human atrium (17); thus *I*_{h} may support the robust pacemaking and driving of BPs (as well as the SA node) in the atrium with relatively small *I*_{K1}, but not in the ventricle with relatively large *I*_{K1}. Possible roles of *I*_{h} in the working of ventricular BPs could be to modulate pacemaker frequency in response to autonomic regulations, as well as to prevent excess hyperpolarization, rather than to improve the structural stability and driving ability by destabilizing EPs or yielding saddle-node bifurcations.

#### I_{Ca,T} effect is NP cell g_{K1} dependent and characterized by shrinkage of voltage oscillations.

Similar to *I*_{h}, *I*_{Ca,T} did not improve the structural stability of BP cells to electrotonic loads of the normal HVM by shifting the Hopf bifurcation point at lower *G*_{C} values (Fig. 9, *top*); the low efficacy of *I*_{Ca,T} appeared to be due to the rapid inactivation kinetics, because the modified *I*_{Ca,T} with slowed inactivation significantly increased the Hopf bifurcation value (Fig. 11). Nevertheless, higher densities of *I*_{Ca,T} improved the structural stability and driving ability of BP cells by producing a saddle-node bifurcation. The effects of *I*_{Ca,T} also depend on the NP cell *I*_{K1}, as well as the size factor (Figs. 9, 12, and 13); *I*_{Ca,T} may effectively improve the functions of BPs in the atrium, whereas much larger *I*_{Ca,T} is required for the ventricle. However, higher densities of *I*_{Ca,T} caused shrinkage in BP oscillations and driven action potentials of the NP cell (Fig. 10), which is probably due to the rapid kinetics of *I*_{Ca,T} inactivation and is not desirable for functional BPs.

#### I_{st} most dramatically improves structural stability and driving ability of BP cells.

In contrast to *I*_{h} or *I*_{Ca,T}, *I*_{st} at relatively low conductance dramatically improved the structural stability and driving ability of BP cells by shifting Hopf bifurcation points toward higher *G*_{C} values and yielding a saddle-node bifurcation (Fig. 9), producing BP oscillations with stable frequency over the broad *G*_{C} range (Fig. 10). Furthermore, *I*_{st}, as well as *I*_{Ca,T} and *I*_{Ca,LD}, may contribute to reducing the critical size of a BP system (i.e., the number of BP cells) to be required for driving the ventricle: with *I*_{st} at 0.01 nS/pF, four BP cells are enough to yield robust pacemaking and driving of a normal HVM via creation of a saddle-node bifurcation, whereas more than five BP cells are required without a pacemaker current (see Fig. 12). Our results suggest that *I*_{st}, exhibiting relatively slow inactivation kinetics, most dramatically improves the structural stability and driving ability of the HVM pacemaker. *I*_{st} may be indispensable for robust pacemaking and driving of ventricular BPs.

#### I_{Ca,LD} effect is similar to I_{st} effect but characterized by generation of irregular dynamics.

The effects of *I*_{Ca,LD} were qualitatively similar to those of *I*_{st}, although its potency was much smaller (Figs. 9, 12, and 13). However, *I*_{Ca,LD} at high densities caused unstable CL and irregular dynamics during increases in *G*_{C}, which is possibly a consequence of producing Ca^{2+} overload conditions, not desirable for functional BPs. We suppose that the role of *I*_{Ca,LD} in the working of pacemaker systems is essentially the same as that of *I*_{Ca,LC}: *I*_{Ca,LD} may play a major role in producing the action potential upstroke as well as excitation-contraction coupling, rather than in supporting robust pacemaking and driving.

### Implications and Significance of Study

In the present study, we used stability and bifurcation analyses to investigate the effects of pacemaker currents on creation and modulations of the HVM pacemaker. Previous reports indicate that exploring bifurcation structures of SA node or ventricular models is useful for general understanding and systematic descriptions of the dynamical mechanisms of normal and abnormal pacemaker activities, e.g., the roles of individual currents in pacemaking (4, 9, 21, 22, 24, 45). Pacemaker currents are not indispensable for basal pacemaking of SA node cells (20, 21) or BP cells (22, 46). As shown in this study for BP cells, however, the pacemaker currents may be necessary for constructing a pacemaker system with high structural stability and high driving ability, as well as for finer regulations of pacemaking rates.

Bifurcation analyses of model cell systems may also allow us to accurately predict and properly control the dynamics and bifurcation structures of real cells, possibly applicable to cell system engineering to develop functional BPs from native cardiomyocytes or ES cells for gene or cell therapy of bradyarrhythmias (8, 10, 15, 25, 28, 32, 36, 50). Pacemaker currents are not necessarily required for creating BP cells (22, 46). Nevertheless, stable transfection of pacemaker channel genes or transplantation of BP cells expressing pacemaker channels (e.g., ES cell-derived SA node cells) to host tissues may be helpful or even necessary for creating functional BPs with robust pacemaking and driving, as well as responsiveness to autonomic regulations. This study suggests that *1*) expression of *I*_{h} (or the other pacemaker currents) is not indispensable but is helpful for creating BP cells from HVMs, *2*) expression of *I*_{h} improves structural stability and driving ability of BPs when *I*_{K1} density in adjacent NP cells is relatively small, as in the atrium, and *3*) expression of *I*_{st}, i.e., a pacemaker current with relatively slow inactivation kinetics, may especially be effective in yielding robust BP activity with stable frequency.

### Limitations and Perspectives of Study

As mentioned in our previous article (22), the limitations of our study include incompleteness of the HVM model due to the lack of experimental data from HVMs as well as the lack of experimental evidence for bifurcation structures of real myocytes. More elaborate HVM models, such as recently published ones (12, 42), with refinements of formulas to be suitable for bifurcation analyses would have to be developed for more detailed investigations using the theoretical approach. In addition, the pacemaker current formulas based on experiments for human cardiomyocytes should be used in future studies.

As suggested by Silva and Rudy (41), a potential advantage of the BP over the electronic pacemaker is responsiveness to β-adrenergic stimulation (β-AS). BP responsiveness to β-AS was very limited in the guinea pig ventricular model (41). We also found in the simulations according to their method that the HVM pacemaker did not exhibit a significant rate increase in response to β-AS, suggesting that the ventricular pacemaker is less sensitive to β-AS than the SA node pacemaker. Such low sensitivity to β-AS of the ventricular pacemaker would be due at least in part to the lack of the pacemaker currents that are abundant in SA node primary pacemaker cells (20, 21, 27, 31, 51) and are known to be enhanced by β-AS (26, 35, 43). Nevertheless, we did not show the data on this issue, because it is very important but beyond the aim of this study and because our model lacks the intracellular modulating factors such as cAMP and protein kinases. Further sophisticated models incorporating the modulating factors, as developed by Saucerman et al. (38), would be required.

We used the simple coupled-cell system to investigate the structural stability and driving ability of BP cells, believing that this study is valuable as a first step toward development of a theoretical background for engineering of functional BPs. Nevertheless, a real pacemaker system such as the intact SA node has much more complex architectures to facilitate optimization of the electrical loading by surrounding atrial or ventricular tissues (1, 2, 14, 44). Thus more elaborate multicellular models are required for investigation of the structural stability to electrotonic loads and ability to drive the heart of a BP system in vivo and of how to create BP systems with robust pacemaking and driving. In the coupled-cell model, the critical size factor (the critical number of BP cells) for stable driving of a normal HVM was quite large, e.g., at least six BP cells were required in the absence of a pacemaker current (see Figs. 8*A* and 12). A similar size effect, a size factor of 5, was previously reported for other coupled-cell systems composed of SA node cells and ventricular myocytes (19, 47). However, the coupled-cell model cannot predict the critical size of a BP to be required for stable driving of the ventricle. Using a two-dimensional system in which real atrioventricular node or model SA node cells were coupled to ventricular or atrial model cells, Joyner et al. (14) determined the critical size of the pacemaker cells for driving the surrounding quiescent cells and suggested that the critical size depends on the scale of the systems, *I*_{K1} density in surrounding NP cells, *G*_{C} values, and many other factors. Moreover, the roles of pacemaker currents in pacemaking and driving of BPs may also change depending on the dimension of the systems, as well as the number of BP cells, *I*_{K1} of adjacent NP cells, and *G*_{C} values (see Figs. 12 and 13). Further studies using multicellular (2 or 3 dimensional) models are required to determine the critical BP size for driving the ventricle and to determine how incorporating pacemaker currents affects critical BP size and BP functions.

## APPENDIX A: GLOSSARY

### General

*F*- Faraday constant
*V*- Membrane potential (mV)

### Cell Geometry

*C*_{m}- Cell membrane capacitance
*C*_{BP}*C*_{m}of BP cells*C*_{NP}*C*_{m}of NP cells- V
_{i} - Myoplasmic volume available for Ca
^{2+}diffusion - V
_{rel} - Volume of junctional SR (Ca
^{2+}release store) - V
_{up} - Volume of network SR (Ca
^{2+}uptake store)

### Ionic Concentrations

- [Na
^{+}]_{o} - Extracellular Na
^{+}concentration - [K
^{+}]_{o} - Extracellular K
^{+}concentration - [Ca
^{2+}]_{o} - Extracellular Ca
^{2+}concentration - [Na
^{+}]_{i} - Intracellular Na
^{+}concentration - [K
^{+}]_{i} - Intracellular K
^{+}concentration - [Ca
^{2+}]_{i} - Myoplasmic Ca
^{2+}concentration - [Ca
^{2+}]_{rel} - Ca
^{2+}concentration in junctional SR - [Ca
^{2+}]_{up} - Ca
^{2+}concentration in network SR

### Equilibrium Potentials

*E*_{Na}- Equilibrium (Nernst) potential for Na
^{+} *E*_{K}- Equilibrium (Nernst) potential for K
^{+}

### Sarcolemmal Ionic Currents

*g*_{h}- Maximum
*I*_{h}conductance *g*_{Ca,T}- Maximum
*I*_{Ca,T}conductance *g*_{st}- Maximum
*I*_{st}conductance *g*_{Ca,LD}- Maximum
*I*_{Ca,LD}conductance *d*_{L}- Activation gating variable for
*I*_{Ca,L} *f*_{L}- Voltage-dependent inactivation gating variable for
*I*_{Ca,L} *d*_{T}- Activation gating variable for
*I*_{Ca,T} *f*_{T}- Inactivation gating variable for
*I*_{Ca,T} *p*_{a}- Activation gating variable for
*I*_{Kr} *n*- Activation gating variable for
*I*_{Ks} *q*- Inactivation gating variable for
*I*_{to} *y*- Activation gating variable for
*I*_{h} *q*_{a}- Activation gating variable for
*I*_{st} *q*_{i}- Inactivation gating variable for
*I*_{st} - α
_{qa} - Opening rate constant for
*q*_{a} - β
_{qa} - Closing rate constant for
*q*_{a} - α
_{qi} - Opening rate constant for
*q*_{i} - β
_{qi} - Closing rate constant for
*q*_{i} *d*_{LD}- Activation gating variable for
*I*_{Ca,LD} *f*_{LD}- Voltage-dependent inactivation gating variable for
*I*_{Ca,LD} *f*_{Ca,∞}- Ca
^{2+}-dependent inactivation gating variable for*I*_{Ca,LD} *h*- Inactivation gating variable for
*I*_{Na} *x*_{∞}- Steady-state value of a gating variable
*x* - τ
_{x} - Time constant for a gating variable
*x* *J*_{Ca,net}- Net Ca
^{2+}flux through the sarcolemmal membrane *J*_{Na,net}- Net Na
^{+}flux through the sarcolemmal membrane

### Intracellular Ca^{2+} Dynamics (SR Function and Ca^{2+} Buffering)

*J*_{rel}- Ca
^{2+}release flux from junctional SR to subspace *J*_{up}- Ca
^{2+}uptake flux from myoplasm to network SR *J*_{tr}- Ca
^{2+}transfer flux from network SR to junctional SR *J*_{leak}- Ca
^{2+}leak flux from network SR to myoplasm - [CM]
_{tot} - Total calmodulin concentration
- [CQ]
_{tot} - Total calsequestrin concentration
- [TC]
_{tot} - Total concentration of troponin-Ca complex
*f*_{TC}- Fractional occupancy of troponin-Ca complex by Ca
^{2+} *f*_{CM}- Fractional occupancy of calmodulin by Ca
^{2+}in myoplasm *f*_{CQ}- Fractional occupancy of calsequestrin by Ca
^{2+} *K*_{dCM}- Ca
^{2+}-binding constant for calmodulin *K*_{dCQ}- Ca
^{2+}-binding constant for calsequestrin - B
_{CM} - Scaling factor for fast Ca
^{2+}buffering by calmodulin - B
_{CQ} - Scaling factor for fast Ca
^{2+}buffering by calsequestrin

## APPENDIX B: MODEL EQUATIONS

The mathematical expressions used for the pacemaker currents (*I*_{h}, *I*_{Ca,T}, *I*_{st}, *I*_{Ca,LD}), as well as differential equations for state variables of the coupled-cell system, are given below. Units are millivolts, picoamperes, nanosiemens, picofarads, milliseconds, millimolar, and liters. The temperature assumed for the model is 37°C. The symbols used and their definitions are the same as those in our rabbit SA node model (20). Expressions for other currents and dynamics of SR Ca^{2+} handling, as well as standard parameter values and initial conditions for computations, are given in our previous article (22). See appendix a for symbol definitions.

### Hyperpolarization-Activated Current (I_{h})

#### Rabbit SA node I_{h} from Kurata et al. (20) model [adopted from Wilders et al. (49) model].

(A1) (A2) (A3) (A4) (A5)

#### Human ventricle I_{h} from Cerbai et al. (3).

(A6) (A7) (A8) (A9) (A10)

#### Human HCN2 current from Moroni et al. (29).

(A11) (A12) (A13) (A14) (A15)

### T-type Ca^{2+} Channel Current (I_{Ca,T})

#### Rabbit SA node I_{Ca,T} from Kurata et al. (20) model [adopted from Demir et al. (5) model].

(A16) (A17) (A18) (A19)

#### Rabbit ventricle I_{Ca,T} from Pugligi-Bers (34) model.

(A20) (A21) (A22) (A23)

### Sustained Inward Current (I_{st})

#### Rabbit SA node I_{st} from Kurata et al. (20) model.

(A24) (A25) (A26) (A27) (A28) (A29)

#### Rat SA node I_{st} from Shinagawa et al. (40).

(A30) (A31) (A32) (A33) (A34) (A35) (A36) (A37)

### Low-Voltage-Activated L-Type Ca^{2+} Channel Current (I_{Ca,LD})

(A38) (A39) (A40) (A41) (A42)

### Differential Equations for State Variables

#### Membrane potential (V).

(A43) (A44) (A45)

#### Gating variables.

(A46)

#### Intracellular ion concentrations.

(A47) (A48) (A49) (A50) (A51) (A52) (A53) (A54)

## APPENDIX C: DETERMINATION OF EPS AS INITIAL CONDITIONS

The methods for determination of EPs as initial conditions for the single cell were described previously (22). To determine an EP of the coupled-cell system, steady-state values of the state variables *V*, [Ca^{2+}]_{i}, [Na^{+}]_{i}, [Ca^{2+}]_{rel}, and [Ca^{2+}]_{up} for BP and NP cells were calculated numerically by the algebraic method with a nonlinear equation solver available in MATLAB 7. For bifurcation analyses and simulations using the [K^{+}]_{i}-fixed system, steady-state values of the state variables for BP and NP cells were computed by the following algebraic equations derived from the differential equations. (A55) (A56) (A57) (A58) (A59) (A60) *Equations A57*–*A60* are for either the BP or the NP cell; in total, 10 equations are required for calculating the steady state of the coupled-cell system.

## 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 (to Y. Kurata and T. Shibamoto), Kanazawa Medical University Grant for Collaborative Research C2005-1 (to Y. Kurata and T. Shibamoto), and Kanazawa Medical University Grant for Promoted Research S2006-6 (to Y. Kurata).

## Acknowledgments

Present address of H. Matsuda: Group in Leading Project for Biostimulation, Kyoto Univ., Kyoto, Japan.

## 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 © 2007 by the American Physiological Society