This Article Abstract Full Text (PDF) All Versions of this Article: biophysj.105.060830v1 89/4/2865    most recent Alert me when this article is cited Alert me if a correction is posted Services Similar articles in this journal Similar articles in PubMed Alert me to new issues of the journal Download to citation manager Citing Articles Citing Articles via HighWire Citing Articles via Google Scholar Google Scholar Articles by Kurata, Y. Articles by Shibamoto, T. Search for Related Content PubMed PubMed Citation Articles by Kurata, Y. Articles by Shibamoto, T.

Dynamical Mechanisms of Pacemaker Generation in I_K1-Downregulated Human Ventricular Myocytes: Insights from Bifurcation Analyses of a Mathematical Model

Yasutaka Kurata ^*, Ichiro Hisatome , Hiroyuki Matsuda ^* and Toshishige Shibamoto ^*

^* Department of Physiology, Kanazawa Medical University, Ishikawa 920-0293, Japan; and Division of Regenerative Medicine and Therapeutics, Tottori University Graduate School of Medical Science, Yonago 683-0826, Japan

Correspondence: Address reprint requests to Yasutaka Kurata, Dept. of Physiology, Kanazawa Medical University, 1-1 Daigaku, Uchinada-machi, Kahoku-gun, Ishikawa 920-0293, Japan, E-mail: yasu{at}kanazawa-med.ac.jp.

	ABSTRACT

TOP ABSTRACT INTRODUCTION THEORY AND METHODS RESULTS DISCUSSION APPENDIX 1: MODEL EQUATIONS APPENDIX 2: DETERMINATION OF... ACKNOWLEDGEMENTS REFERENCES

 
Dynamical mechanisms of the biological pacemaker (BP) generation in human ventricular myocytes were investigated by bifurcation analyses of a mathematical model. Equilibrium points (EPs), periodic orbits, stability of EPs, and bifurcation points were determined as functions of bifurcation parameters, such as the maximum conductance of inward-rectifier K⁺ current (I_K1), for constructing bifurcation diagrams. Stable limit cycles (BP activity) abruptly appeared around an unstable EP via a saddle-node bifurcation when I_K1 was suppressed by 84.6%. After the bifurcation at which a stable EP disappears, the I_K1-reduced system has an unstable EP only, which is essentially important for stable pacemaking. To elucidate how individual sarcolemmal currents contribute to EP instability and BP generation, we further explored the bifurcation structures of the system during changes in L-type Ca²⁺ channel current (I_Ca,L), delayed-rectifier K⁺ currents (I_K), or Na⁺/Ca²⁺ exchanger current (I_NaCa). Our results suggest that 1), I_Ca,L is, but I_K or I_NaCa is not, responsible for EP instability as a requisite to stable BP generation; 2), I_K is indispensable for robust pacemaking with large amplitude, high upstroke velocity, and stable frequency; and 3), I_NaCa is the dominant pacemaker current but is not necessarily required for the generation of spontaneous oscillations.

	INTRODUCTION

TOP ABSTRACT INTRODUCTION THEORY AND METHODS RESULTS DISCUSSION APPENDIX 1: MODEL EQUATIONS APPENDIX 2: DETERMINATION OF... ACKNOWLEDGEMENTS REFERENCES

 
The cardiac biological pacemaker (BP) was recently created by genetic suppression of the inward-rectifier K⁺ current (I_K1) in guinea pig ventricular myocytes (1), suggesting possible development of the functional BP as a therapeutic alternative to the electronic pacemaker. A first step for creation of the functional BP would be engineering of single BP cells, which requires deep understanding of the BP mechanisms. Using the Luo-Rudy guinea pig ventricle model (2), Silva and Rudy (3) simulated BP activity by reducing I_K1 conductance (g_K1) and investigated the ionic mechanisms of BP generation in the I_K1-downregulated ventricular myocyte. They reported that BP activity was yielded by 81% suppression of I_K1, concluding that Na⁺/Ca²⁺ exchanger current (I_NaCa) was the dominant pacemaker current. However, the mechanistic difference between ventricular pacemaking and natural sinoatrial (SA) node pacemaking, as well as whether I_K1 downregulation also induces BP activity in human ventricular myocytes (HVMs), remains to be clarified.

The aim of this study was to elucidate the mechanisms of BP generation in I_K1-downregulated HVMs and the roles of individual sarcolemmal currents in HVM pacemaking in terms of the nonlinear dynamics and bifurcation theory. In previous studies (4–11), bifurcation structures of ventricular or SA node models, i.e., ways of changes in the number or stability of equilibrium and periodic states of the model systems, were investigated for elucidating the mechanisms of normal and abnormal pacemaker activities. These theoretical works indicate that the mathematical approach provides a convenient way of understanding how individual currents contribute to pacemaker activities. In this study, therefore, local stability and bifurcation analyses, as well as numerical simulations, were performed for a mathematical model of the HVM. We constructed bifurcation diagrams by calculating equilibrium points (EPs), periodic orbits, stability of the EP, and saddle-node or Hopf bifurcation points as functions of bifurcation parameters (for details, see Theory and Methods). During I_K1 suppression, BP activity abruptly emerged around an unstable EP via a saddle-node bifurcation at which a stable EP corresponding to the resting state disappeared. Our results suggested that BP activity could be developed by reducing I_K1 alone in HVMs as well and that the instability of an EP at depolarized potentials is essentially important for BP generation.

To elucidate the ionic mechanisms of EP destabilization and BP generation in the I_K1-downregulated HVM, we further explored bifurcation structures of the I_K1-reduced BP system during decreases or increases in an L-type Ca²⁺ channel current (I_Ca,L), delayed-rectifier K⁺ currents (I_K), and I_NaCa, which appear to be essentially important for pacemaker generation (3,11). Moreover, we determined stability of the BP system at the steady-state potential (V₀) during applications of the constant bias current (I_bias) and how the unstable V₀ and I_bias regions, where BP oscillations occur, change with decreasing or increasing the sarcolemmal currents. This would reveal the contributions of each current to EP instability and robustness of BP activity to hyperpolarizing or depolarizing loads (11). Our study suggests that the dynamical mechanism of the ventricular pacemaking is essentially the same as that of the natural SA node pacemaking as reported by Kurata et al. (11).

This article clearly shows that bifurcation analyses of a mathematical model allow us to notice the essential importance of EP instability for pacemaker generation and to elucidate the roles of individual ionic currents in pacemaking. Thus, the nonlinear dynamical approach is useful for general understanding of the mechanisms of normal and abnormal pacemaker activities and may also be applicable to engineering of a functional BP as a therapeutic alternative to the electronic pacemaker. Definitions of the terms specific to the nonlinear dynamics and bifurcation theory are provided at the end of the Theory and Methods section to help understand the theory and methods for bifurcation analysis as well as our results and conclusions. Abbreviations and acronyms repeatedly used in this article are listed in Table 1.

	THEORY AND METHODS

TOP ABSTRACT INTRODUCTION THEORY AND METHODS RESULTS DISCUSSION APPENDIX 1: MODEL EQUATIONS APPENDIX 2: DETERMINATION OF... ACKNOWLEDGEMENTS REFERENCES

(1)

This is usually written in the vector form

(2)

where the state variable x is a vector-valued function of time t, and the vector field f is a function of the state variable x (12–14). For a model system to be suitable for bifurcation analysis, the vector field f should be continuous and smooth (i.e., sufficiently differentiable) and should not depend on time or initial conditions but depend only on the state variable x. Such a system is called an "autonomous" system (see also Definitions of Terms at the end of this section).

On the basis of single-cell patch-clamp data from undiseased and failing HVMs, Priebe and Beuckelmann (15) have first developed a single HVM model (referred to as the Priebe-Beuckelmann (PB) model) as a modified version of the Luo-Rudy phase II model for the guinea pig ventricle (2). Recently, ten Tusscher et al. (16) and Iyer et al. (17) developed more elaborate HVM models based on detailed experimental data. Their models appear to be superior to the PB model in reproducing experimental data but less suitable for bifurcation analyses for the following reasons: 1), the ten Tusscher et al. model has the vector field containing many complex functions that are not continuous or smooth, thus they are not always differentiable or yield noncontinuous derivatives; and 2), the Iyer et al. model (n = 67) is much larger than the PB model (n = 15) or the ten Tusscher et al. model (n = 17), making bifurcation analyses practically much harder. We have therefore chosen to modify the PB model on the basis of recent experimental findings as well as other human or animal heart models. The original PB model explicitly contains time t in the formula for the sarcoplasmic reticulum (SR) Ca²⁺ release, making it nonautonomous. Thus, modifications of the PB model were required for converting it into an autonomous system, as well as for improving the capability of reproducing experimental data.

The standard model for the normal activity is described as a nonlinear dynamical system of 15 first-order ordinary differential equations. The membrane current system includes I_Ca,L, the rapid and slow components of I_K (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²⁺ (I_Ca,b) currents, Na⁺-K⁺ pump current (I_NaK), I_NaCa, and Ca²⁺ pump current (I_pCa). The expressions for I_Ca,L and I_Ks as well as SR Ca²⁺ release were reformulated, whereas the formulas for other ionic currents and intracellular Ca²⁺ handling are essentially the same as those in the original PB model or adopted from other existing models (16,18–20). Modifications of the model are summarized in Table 2; details on major modifications are described below. All expressions (Eqs. 3–64) and the standard parameter values used are provided in Appendix 1 and Table 5, respectively.

L

L

L

L

Ca,

View larger version (28K): [in this window] [in a new window]   FIGURE 1  Kinetics of ICa,L. (Top) Voltage dependence of steady-state probabilities (dL, fL) and time constants for ICa,L activation (dL) and inactivation (fL). Equations used for the current (modified) model are shown with the thick lines (labeled "PM"). For comparison, those for existing (original) models are also shown with the thin lines: PB, Priebe and Beuckelmann (15); TNNP, ten Tusscher et al. (16); CRN, Courtemanche et al. (19). The experimental data for dL, fL, and fL are from Bénitah et al. (21) (circles), Pelzmann et al. (22) (squares), and Li et al. (23) (hexagons). (Bottom) Computed voltage-clamp records for ICa,L (left) and peak ICa,L-V relationship (right). Currents were evoked by 300-ms step pulses from a holding potential of –50 mV to test potentials ranging from –30 to +50 mV in 10 mV increments. Simulating the whole-cell perforated-patch recording, [Ca2+]i was not clamped, whereas [Na+]i and [K+]i were fixed at 5 and 140 mM, respectively. The experimental data for the peak ICa,L-V relation are from Bénitah et al. (21) (open circles), Pelzmann et al. (22) (squares), and Li et al. (23) (hexagons).

View larger version (31K): [in this window] [in a new window]   FIGURE 2  Kinetics of IKr (top) and IKs (bottom). (Left, middle) Voltage dependence of steady-state probabilities and time constants of IKr activation (pa,, pa) and inactivation (pi,) as well as IKs activation (n, n). The thick lines are for this study (PM); the thin lines are from the existing models (PB, TNNP, CRN). The experimental data for pa,, pa, and n are from Li et al. (24) (circles) and Wang et al. (25) (squares). (Right) Computed voltage-clamp records for IKr and IKs. Currents were elicited by 3-s step pulses from a holding potential of –60 mV to test potentials ranging from –50 to +50 mV in 10 mV increments.

Formulation of I_Kr
The kinetics of I_Kr was described by Eqs. 9–14 with the activation (p_a) and inactivation (p_i) gating variables. Voltage dependences of I_Kr activation and inactivation in the current and original human heart models are shown in Fig. 2 (top), together with experimental data for comparison; the model-generated I_Kr during voltage-clamp pulses are also depicted. There are limited experimental data available for the gating kinetics of native I_Kr channels in HVMs. Nevertheless, experimentally observed activation of native I_Kr in human cardiac myocytes (24,25) appeared to be faster than the I_Kr activation simulated by the original PB formulas. To formulate I_Kr kinetics, the recently developed HVM models (16,17) employed the experimental data from human ether-a-go-go-related gene (HERG) channels expressed in human embryonic kidney (HEK) cells (26,27), Chinese hamster ovary cells (28), or Xenopus oocytes (29). However, voltage- and time-dependent kinetics of expressed HERG channels appear to depend on cell lines or conditions for expression, as well as conditions for current recordings such as temperatures, and thus may be different from those of native I_Kr channels in human hearts (26–28,30–32). We therefore adopted the data from HVMs of Li et al. (24) for the steady-state activation curve (p_a,) and the formula of Courtemanche et al. (19) based on the data from human atrial myocytes of Wang et al. (25) for the activation time constant (_pa). We also employed the expression of Courtemanche et al. (19) for the steady-state inactivation curve (p_i,). No detailed experimental data are available on the time constant of I_Kr inactivation (_pi) in intact HVMs. Inactivation and recovery of native I_Kr in animal hearts or expressed HERG channels are very rapid (28,33); thus, we assumed that I_Kr inactivation is instantaneous. The maximum I_Kr conductance (g_Kr) was set equal to 12 pS/pF, according to the data of Li et al. (24). With this g_Kr value, the complete block of I_Kr prolonged the action potential duration (APD) of the model HVM paced at 1 Hz by 34.6–48.9%; the simulated APD prolongation by I_Kr block is comparable to the experimentally observed effects of the selective I_Kr blocker E-4031 as reported by Li et al. (24).

Formulation of I_Ks
On the basis of the recent data from Li et al. (24) and Virág et al. (34), the kinetics of I_Ks was reformulated as Eqs. 15–18. Voltage dependences of I_Ks kinetics in the current and original models are shown in Fig. 2 (bottom), together with experimental data for comparison; the model-generated I_Ks during voltage-clamp pulses are also depicted. Steady-state I_Ks activation curve (n) is from Li et al. (24). Virág et al. (34) recently reported that I_Ks in HVMs slowly activates and rapidly deactivates. According to their report, therefore, we reformulated the time constant of I_Ks activation (_n) as Eq. 18. The expression of E_Ks is the same as for the PB model. The maximum conductance (g_Ks) was set equal to 36 pS/pF so as to yield the APD at 90% repolarization (APD₉₀) during 1 Hz pacing of 336 ± 16 ms, as reported by Li et al. (24).

Formulation of SR Ca²⁺ release
In the original PB model, SR Ca²⁺ release kinetics was represented by the formula that explicitly contains time t and also [Ca²⁺]_i,2, which is not a state variable but the sum of net Ca²⁺ influx during the first 2 ms after initiation of the action potential (AP). Thus, the original PB model is a nonautonomous system, the vector field of which depends on time and initial conditions, not suitable for bifurcation analysis. We had therefore to renew the SR Ca²⁺ release formula to convert the model into an autonomous system. Owing to the lack of available data for updating the kinetic formulation of SR Ca²⁺ release in HVMs, we utilized simple expressions, Eqs. 47–51. The formula for conductance of the Ca²⁺ release channel (g_rel) was adopted from Faber and Rudy (20), with the P_rel value reduced to one-third of the original value to obtain the peak [Ca²⁺]_i transient of 1 µM during APs elicited at 1 Hz. For gating behaviors of the Ca²⁺ release channel (d_R, f_R), we used the expressions similar to those for the I_Ca,L gating variables to make the model an autonomous system suitable for bifurcation analyses.

Ion concentration homeostasis
The model also includes material balance expressions to define the temporal variation in [Ca²⁺]_i, [Na⁺]_i, and [K⁺]_i, whereas [Ca²⁺]_o, [Na⁺]_o, and [K⁺]_o were fixed at 2, 140, and 5.4 mM, respectively. As pointed out by Hund et al. (35) and Krogh-Madsen et al. (36), the second-generation models incorporating ion concentration changes have two major problems: 1), drift, with very slow long-term trends in state variables; and 2), degeneracy, with nonuniqueness of steady-state solutions. The current model did not show a long-term drift in the intracellular ion concentrations but reached a steady state during pacing, as well as during BP oscillation, when the principle of charge conservation was taken into account (35).

An n-dimensional fully differential system formulated as a cardiac second-generation model can usually be converted into a differential-algebraic system composed of n – 1 differential equations and one algebraic equation, resulting in n – 1 equations for n unknowns (35,36). In such a system, the Jacobian matrix of which is singular, the EP or limit cycle is not unique but depends on initial conditions; e.g., our full system has a continuum of EPs, because Eq. 65 can be derived from Eqs. 66, 69, and 70 (see Appendix 2). Thus, second-generation models including our full system exhibit degeneracy not suitable for bifurcation analysis to be applicable to isolated equilibria. As suggested by Krogh-Madsen et al. (36), one of the ways to remove degeneracy and thus allow bifurcation analyses of isolated equilibria is to make some ionic concentrations fixed. Variations in [K⁺]_i during changes in parameters or stimulation rates (0.25–2 Hz) were relatively small (<5 mM), being too small to significantly alter the model cell behaviors. For stability and bifurcation analyses, therefore, [K⁺]_i was fixed at 140 mM. The degenerate system with the fixed [K⁺]_i has the finite number of EPs and limit cycles, being suitable for bifurcation analyses.

Determination and validation of model HVM dynamics
Numerical methods for dynamic simulations of HVM behaviors
Dynamic behaviors of the model HVM were determined by solving a system of nonlinear ordinary differential equations numerically. Numerical integration, as well as bifurcation analyses, were performed on Power Macintosh G4 computers (Apple Computer, Cupertino, CA) with MATLAB 5.2 (The MathWorks, Natick, MA). We used the numerical algorithms available as MATLAB ODE solvers: 1), a fourth-order adaptive-step Runge-Kutta algorithm which includes an automatic step-size adjustment based on an error estimate (37), and 2), a variable time-step numerical differentiation approach selected for its suitability to stiff systems (38). The former (named ode45) is the best function for most problems. However, the latter (named ode15s) was much more efficient than ode45 and both solvers usually yielded nearly identical results. Therefore, ode15s was usually used, and ode45 was only sometimes used to confirm the accuracy of calculations. The maximum relative error tolerance for the integration methods was set to 1 x 10^–6.

Dynamic properties of the model HVM: simulated APs, ionic currents, and [Ca²⁺]_i dynamics
The steady-state AP, sarcolemmal currents, and [Ca²⁺]_i transient of the current model paced at 1 Hz with the standard parameter values are shown in Fig. 3 A. According to the report of Hund et al. (35), the stimulus current was assumed to carry K⁺ ions into the cell for charge conservation. During pacing, the model cell dynamics reached a steady state, with no long-term drift in the state variables. Steady-state values of the resting potential, maximum upstroke velocity, and APD₉₀ during 1 Hz pacing were –84.9 mV, 411 V/s, and 338 ms, respectively. The diastolic [Ca²⁺]_i and peak [Ca²⁺]_i transient in the model HVM paced at 1 Hz were 0.169 and 0.961 µM, respectively. The simulated AP parameters and [Ca²⁺]_i dynamics of the current and original HVM models are listed in Table 3, together with corresponding experimental data for comparison. The AP parameters and [Ca²⁺]_i dynamics of the current model appear to be in reasonable agreement with the mean experimental values recently determined for single HVMs, as well as those of the original PB and other HVM models. The values of [Ca²⁺]_rel and [Ca²⁺]_up were 0.41 mM (identical) for the resting state and 0.59–2.68 and 2.84–2.92 mM, respectively, in a steady state during 1 Hz pacing. These values are comparable to the experimental data from rat and rabbit ventricular myocytes (40,41), as well as those in the original PB model, although experimental values for HVMs are unknown.

View larger version (28K): [in this window] [in a new window]   FIGURE 3  Simulated dynamics of the model HVM. (A) Steady-state behaviors of the AP, underlying sarcolemmal currents, and [Ca2+]i transient. The model HVM was paced at 1 Hz with 1-ms stimuli of 80 pA/pF. Differential equations (Eqs. 55–64) were numerically solved for 20 min with the initial conditions appropriate to a resting state, which were determined by Eqs. 65–69 with [K+]i fixed at 140 mM. Model cell behavior after the last stimulus (during the last AP) is depicted. (B) Rate dependence of APD, ICa,L, and [Ca2+]i transients. The model HVM was paced at 0.5, 1, and 2 Hz with 1-ms stimuli of 80pA/pF. The differential equations were numerically solved for 20 min at each pacing rate; model cell behaviors after the last stimulus are depicted.

The peak [Ca²⁺]_i transient predicted by the current model was a little smaller at 2 Hz than at 1 Hz. In most experiments, however, the peak [Ca²⁺]_i transient and/or developed tension increased as the pacing rate increased up to 2.5 Hz (45–47). This inconsistency may result from the lack of intracellular modulating factors such as those involved in the rate-dependent potentiation of I_Ca,L (48,49). Alternatively, the kinetic formulation of I_Ca,L and/or SR Ca²⁺ handling may be inappropriate. We found that the use of the _fL formula from the ten Tusscher et al. model (16) led to the successful reproduction of the rate-dependent increase in the peak [Ca²⁺]_i transient. Nevertheless, their _fL formula was not used, because it appeared not to fit the experimental data as shown in Fig. 1 and because a very large I_Ks (>8 times larger than the experimental values) was required to counteract the large I_Ca,L during phase 2 (16). In the preliminary study, BP dynamics or bifurcation structures of the HVM model were not essentially altered by the use of different formulas for I_Ca,L (_fL) or the change in the rate dependence of [Ca²⁺]_i transients.

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 of the ionic channels (g_K1, g_Ca,L, g_Kr, g_Ks) and amplitude of I_NaCa or I_bias; the maximum conductance and I_NaCa amplitude are expressed as normalized values, i.e., ratios to the control values.

In the HVM model with the fixed [K⁺]_i, 14 state variables define a 14-dimensional state point in the 14-dimensional state space of the system. We calculated EPs and periodic orbits in the state space. An EP was determined as a point at which the vector field vanishes (i.e., f(x) = 0); steady-state values of the state variables were calculated by Eqs. 65–69 (see Appendix 2). Asymptotic stability of the EP was also determined by computing 14 eigenvalues of a 14 x 14 Jacobian matrix derived from the linearization of the nonlinear system around the EP (for more details, see Vinet and Roberge (7)). Periodic orbits were located with the MATLAB ODE solvers. When spontaneous oscillation (BP activity) appeared, the AP amplitude (APA) as a voltage difference between the maximum diastolic potential (MDP) and peak overshoot potential (POP), as well as the cycle length (CL), was 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 x 10^–3 of the preceding APA and CL values. Potential extrema (MDP, POP) 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 in bifurcation diagrams.

Codimension one and two bifurcations to occur in the model system were explored by constructing bifurcation diagrams with one and two parameters, respectively (10,11). For construction of one-parameter bifurcation diagrams in which values of a state variable at EPs or extrema of periodic orbits are shown as a function of one bifurcation parameter, a bifurcation parameter was systematically changed while keeping all other parameters at their standard values. The membrane potential (V) at EPs (steady-state branches) and local potential extrema of periodic orbits (periodic branches) were determined and plotted for each value of the bifurcation parameter. The saddle-node bifurcation point at which two EPs coalesce and disappear was determined as a point at which the steady-state current-voltage (I/V) curve and zero-current axis come in touch with each other. The Hopf bifurcation point at which the stability of an EP reverses was also detected by the stability analysis as described above. For construction of two-parameter bifurcation diagrams in which codimension one bifurcation points are plotted in a parameter plane, the secondary parameter was systematically changed with the primary parameter fixed at various different values. The path of saddle-node and Hopf bifurcation points was traced in the parameter plane; i.e., bifurcation values for the secondary parameter were plotted as a function of the primary parameter.

Definition and evaluation of structural stability for the BP system
We also evaluated the structural stability of the BP system, which is defined as the robustness of pacemaker activity to various interventions or modifications that may cause a bifurcation to quiescence or irregular dynamics (11). Interventions or modifications leading to quiescence or irregular dynamics include injections of I_bias, electrotonic loads of normal HVMs, current leakage via myocardial injury, and intrinsic changes in channel conductance or gating kinetics. In this study, the structural stability of the BP system was tested for hyperpolarizing and depolarizing loads by exploring bifurcation structures during applications of I_bias.

The way of evaluating the structural stability to I_bias is essentially the same as in our previous study for the SA node pacemaker (11). Changes in V₀ and its stability with I_bias applications were depicted as the I_bias-V₀ curve (for more details, see Fig. 1 of Kurata et al. (11)). There are one or two Hopf bifurcation points and two saddle-node bifurcation points corresponding to the current exrema in the I_bias-V₀ curve. In the unstable I_bias region between two Hopf (or a Hopf and saddle-node) bifurcation points, the system has an unstable EP only, usually exhibiting stable limit cycles without annihilation (6,8). When the system removes to the stable I_bias region, it would come to rest at the stable EP via gradual decline of limit cycles, annihilation, or irregular dynamics, as reported by Guevara and Jongsma (8). Note that system A is considered structurally more stable than system B when the amplitude of I_bias required for stabilizing an EP or causing a bifurcation to quiescence is greater in system A than in system B (11). Thus, the larger the unstable I_bias range over which limit cycles continue is, the more structurally stable the system is.

The Hopf and saddle-node bifurcation points, as well as the control V₀ at I_bias = 0, in the I_bias-V₀ curve were determined as functions of bifurcation parameters and plotted on both the potential and current coordinates, as in previous studies (7,11). Exploring how the unstable V₀ and I_bias regions change with decreasing or increasing the sarcolemmal currents would enable us to determine the contribution of each current to the structural stability of the BP system to hyperpolarizing or depolarizing loads, as well as to EP instability itself, which is evaluated by positive real parts of eigenvalues of Jacobian matrices (see Fig. 4 of Kurata et al. (11)).

View larger version (26K): [in this window] [in a new window]   FIGURE 4  Simulated steady-state BP activities (spontaneous APs, ionic currents, and [Ca2+]i dynamics) in the model HVM at gK1 = 0.15 (left) and 0 (right). Differential equations (Eqs. 55–63) were numerically solved for 20 min at each gK1 with initial conditions appropriate to an EP and a 1-ms stimulus of 1 pA/pF for triggering an AP; model cell behaviors during the last 2.5 s starting from MDP are depicted. Note the differences in the ordinate scales for individual currents.

Nonautonomous system
An nth-order "nonautonomous" continuous-time dynamical system is defined by the state equation dx/dt = f(x, t), where the vector field f depends on time, i.e., explicitly contains time t (13).

EP
A time-independent steady-state point at which the vector field vanishes (i.e., dx/dt = 0) in the state space of an autonomous system, constructing the steady-state branch in one-parameter bifurcation diagrams. This state point corresponds to the zero-current crossing in the steady-state I/V curve, i.e., a quiescent state of a cell if it is stable.

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

Limit cycle
A periodic limit set onto which a trajectory is asymptotically attracted in an autonomous system. A stable limit cycle corresponds to an oscillatory state, i.e., pacemaker activity, of a cell.

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

Hopf bifurcation
This 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 for the EP have a single complex conjugate pair and its real part reverses the sign through zero.

Saddle-node bifurcation
This 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 eigenvalues of a Jacobian matrix is zero.

Codimension
The codimension of a bifurcation is the minimum dimension of the parameter space in which the bifurcation may occur in a persistent way (12,50). In other words, the codimension is the number of independent conditions determining the bifurcation (14). A bifurcation with codimension one is a point in one-parameter bifurcation diagrams such as Fig. 6 or a line in two-parameter bifurcation diagrams such as Fig. 7, whereas a bifurcation with codimension two is a point in the two-parameter space.

View larger version (37K): [in this window] [in a new window]   FIGURE 6  EP stability and BP dynamics of the IK1-removed system during inhibition of ICa,L, IKr, or IKs. (Top) Steady-state I/V relations for gCa,L = 0 – 1 (left), gKr = 0 – 1 (middle), and gKs = 0 – 1 (right), depicted at an interval of 0.2. (Middle, bottom) Bifurcation diagrams with the steady-state (EP1) and stable periodic (MDP, POP) branches (middle), as well as CL (bottom), are shown as functions of gCa,L, gKr, or gKs, which was reduced at an interval of 0.001. The differential equations were numerically solved for 2–3 min at each conductance value. A Hopf bifurcation occurred at gCa,L = 0.134 (labeled HB), with the loci of MDP and POP converging at the bifurcation point.

View larger version (31K): [in this window] [in a new window]   FIGURE 7  Changes in the unstable regions of the Ibias-V0 curve during inhibition of ICa,L, IKr, or IKs. Displacements of Hopf (H1, H2) and saddle-node (SN1, SN2) bifurcation points in the Ibias-V0 curve with decreasing gCa,L, gKr, or gKs at an interval of 0.001 are shown for the IK1-removed BP system. Values of V0 and Ibias at the bifurcation points were plotted on the potential (top) and current (bottom) coordinates, respectively, as functions of the conductance. The path of the control V0 at Ibias = 0 (i.e., EP1 in Fig. 6) is also shown on the potential coordinates, intersecting the locus of H2 at gCa,L = 0.134 (labeled HB). Note that the unstable region exists between H1 and H2 (or H1 and SN2), where spontaneous oscillations occur. A codimension two saddle-node bifurcation at which the loci of H1 and H2 merge together, i.e., the unstable region disappears occurred at gCa,L = 0.126 (labeled SNc2).

	RESULTS

TOP ABSTRACT INTRODUCTION THEORY AND METHODS RESULTS DISCUSSION APPENDIX 1: MODEL EQUATIONS APPENDIX 2: DETERMINATION OF... ACKNOWLEDGEMENTS REFERENCES

Bifurcation structure of HVM system during I_K1 suppression
BP activity abruptly appears around unstable EP via saddle-node bifurcation
To elucidate the dynamical mechanisms of BP generation, we explored EP stability, dynamics, and bifurcation structures of the HVM system during I_K1 suppression by constructing bifurcation diagrams for g_K1. As shown in Fig. 5, the zero-current potential in the steady-state I/V curve, as well as steady-state [Ca²⁺]_i and [Na⁺]_i, was plotted as a function of g_K1; when BP activity appeared, AP parameters (MDP, POP, CL), extrema of [Ca²⁺]_i, and average [Na⁺]_i were also plotted. There are three EPs corresponding to the zero-current crossings of the I/V curve in the state space of the system at g_K1 = 1 (denoted as EP1, EP2, and EP3). The stability analysis revealed that EP1 with the most depolarized potential as well as EP2 is unstable, whereas EP3 with the most negative potential corresponding to the resting state is stable. The zero-current potential at EP3 positively shifted with decreasing g_K1; the stable EP3 (and EP2) disappeared via a saddle-node bifurcation when g_K1 reduced to 0.154. After the bifurcation, the system has only one EP at depolarized potentials (EP1), which is always unstable, not stabilized during g_K1 decreases. Limit cycles (BP oscillations) abruptly emerged around the unstable EP1 on occurrence of the saddle-node bifurcation, with APA and CL gradually decreasing during further reductions in g_K1.

View larger version (26K): [in this window] [in a new window]   FIGURE 5  Bifurcation structure of the model HVM during gK1 decreases. (Left) Steady-state I/V relations for gK1 = 0 – 1, depicted at an interval of 0.1 with [K+]i fixed at 140 mM. The zero-current crossings, corresponding to EPs, are designated by EP1, EP2, and EP3. (Middle) A bifurcation diagram for gK1 with the steady-state (EP1–3) and stable periodic (MDP, POP) branches (top), and CL plotted as a function of gK1 (bottom). EPs were determined by the algebraic equations (Eqs. 65–69). The model cell dynamics were computed by numerically solving the differential equations (Eqs. 55–63) for 10 min at each gK1, which was reduced at an interval of 0.001, with the initial [Na+]i at gK1 = 1 set equal to 6.14 mM. The saddle-node bifurcation point at which EP2 and EP3 merge together and disappear is located (labeled SNB at gK1 = 0.154). (Right) Steady-state [Ca2+]i (top) and [Na+]i (bottom) at EPs as functions of gK1. The minimum (Min) and maximum (Max) of [Ca2+]i, as well as average [Na+]i, during the BP oscillations are also plotted. Note that the steady-state [Na+]i at EP3 reduces with decreasing gK1.

Ionic bases of EP instability and BP generation in I_K1-downregulated HVM
Instability of the EP at depolarized potentials (EP1) appears to be essentially important for BP generation. To elucidate the ionic mechanisms of EP instability and BP generation, therefore, we further explored how modulating individual sarcolemmal currents or [Ca²⁺]_i transients affects EP stability, oscillation dynamics, and bifurcation structures of the BP system. We focused on I_Ca,L, I_K, and I_NaCa, which appear to be essentially important for BP generation.

Influences of blocking I_Ca,L or I_K on stability and dynamics of the BP system
Fig. 6 shows the effects of reducing the maximum conductance of I_Ca,L (g_Ca,L), I_Kr (g_Kr), or I_Ks (g_Ks) on EP stability and BP dynamics of the I_K1-removed system. The zero-current potential (EP1) negatively shifted with decreasing g_Ca,L; the EP became stable via a Hopf bifurcation when g_Ca,L was reduced by 86.6%. During the g_Ca,L decreases, a limit cycle (BP oscillation) gradually contracted in size and finally disappeared at the Hopf bifurcation point. In contrast, a Hopf bifurcation to abolish BP activity did not occur during decreases in g_Kr or g_Ks, although APA significantly reduced with decreasing g_Kr. Spontaneous oscillation continued even when I_Kr or I_Ks was completely blocked, whereas it was abolished by eliminating both I_Kr and I_Ks.

As shown in Fig. 7, we also examined the effects of blocking I_Ca,L, I_Kr, or I_Ks on the unstable V₀ and I_bias regions in the I_bias-V₀ curve for further clarifying the contributions of the individual currents to EP instability and structural stability to depolarizing or hyperpolarizing loads of the BP system. Hopf and saddle-node bifurcation points in the I_bias-V₀ curve determined with decreasing conductance values were plotted on the potential and current coordinates to display how blocking each current shrinks the unstable ranges where BP oscillations occur (for more details, see Theory and Methods). The more prominent shrinkage in the unstable V₀ region is, the greater the contribution of the current to EP instability is; the more prominent shrinkage in the unstable I_bias region is, the greater the contribution to the robustness to depolarizing or hyperpolarizing loads is. The unstable regions, shrinking with decreases in g_Ca,L, disappeared via a codimension two saddle-node bifurcation at a 87.4% reduction of g_Ca,L. Pacemaker activity never appeared in the system with g_Ca,L reduced to less than the saddle-node bifurcation value. In contrast, reducing g_Kr or g_Ks little affected the unstable V₀ range, eigenvalues of Jacobian matrices at V₀ being not significantly altered. The unstable I_bias region shrunk with decreasing g_Kr or g_Ks, though it did not disappear even when I_Kr or I_Ks was completely blocked; the robustness to depolarizing I_bias was attenuated by decreasing g_Kr or g_Ks.

Effects of incorporating various K⁺ currents on stability and automaticity of an I_K-removed system
The results shown in Figs. 6 and 7 suggest that I_Kr or I_Ks contributes little to the instability of the EP leading to the generation of spontaneous oscillations. To elucidate the roles of I_Kr and I_Ks in the HVM pacemaking, we tested how the stability and dynamics of the reduced BP system not including either I_Kr or I_Ks are affected by incorporating different K⁺ currents. Both I_Kr and I_Ks were first removed from the standard I_K1-eliminated BP system, and the background K⁺ current of linear I/V relation (I_K,b) or the original I_Kr or I_Ks (one of the three) was then incorporated. Eliminating both I_Kr and I_Ks abolished BP activity, with the model cell coming to a rest at a stable zero-current potential of +16.15 mV. Fig. 8 shows how EP stability and automaticity of the reduced system change with increasing the K⁺ current conductance up to 10–50 pS/pF. As the K⁺ conductance increased, the zero-current potential was negatively shifted and eventually destabilized via a Hopf bifurcation, with limit cycles emerging at the bifurcation. It should be noted that all the K⁺ currents, including I_K,b, could individually create an unstable EP and induce spontaneous oscillations. However, I_K,b, the effect of which was very similar to the electrotonic influence of the adjacent normal HVM in a coupled-cell system (data not shown), produced the oscillation of relatively small amplitude and low upstroke velocity, as well as unstable CL and irregular dynamics. The I_K,b conductance (g_K,b) range where an EP is unstable and spontaneous oscillations occur was very limited, with further increases in g_K,b leading to quiescence via emergence of a stable EP. In contrast, I_Kr could yield the large amplitude oscillation with relatively high upstroke velocity and stable CL; however, I_Kr caused irregular dynamics at higher g_Kr values (>31.5 pS/pF), whereas I_Ks did not.

View larger version (38K): [in this window] [in a new window]   FIGURE 8  Effects of incorporating various K+ currents with different kinetics (IK,b, IKr, or IKs) on EP stability and dynamics of the IK-removed model cell (gK1 = 0). (Top) Steady-state I/V relations for gK,b = 0–10 pS/pF (left), gKr = 0–50 pS/pF (middle), and gKs = 0–50 pS/pF (right), depicted at an interval of 1 or 5 pS/pF. (Middle, bottom) Oscillatory behaviors, as well as zero-current potentials and their stability, of the model system, calculated during increases in gK,b, gKr, or gKs at an interval of 0.01–0.1 pS/pF. Bifurcation diagrams with the steady-state (EP1–3) and periodic (MDP, POP) branches, as well as Hopf bifurcation (HB) and saddle-node bifurcation (SNB) points, were constructed for gK,b, gKr, or gKs (middle). Maximum upstroke velocity (Vmax) and CL of the potential oscillations were also plotted against gK,b, gKr, or gKs (bottom). The vertical lines for IKr and IKs represent the standard values of gKr (12 pS/pF) and gKs (36 pS/pF). The symbol "ID" designates the irregular dynamics.

View larger version (28K): [in this window] [in a new window]   FIGURE 9  EP stability and BP dynamics of the IK1-removed system during INaCa inhibition, determined for the unmodified system (left) and the modified system in which fCa was fixed at 0.412, the value for [Ca2+]i = 0.5 µM (middle), or [Ca2+]i itself was fixed at 0.5 µM (right). The amplitude of INaCa as a ratio to the control value was decreased from 1 (control) to 0.001 (0.1%); practically, the log(INaCa) value was reduced from 0 (INaCa = 1) to –3 (INaCa = 1 x 10–3) at an interval of 0.01. The zero-current potential (EP1) and BP dynamics (MDP, POP, CL), as well as the minimum (Min) and maximum (Max) of [Ca2+]i and minimum values of the inward currents (ICa,L, INaCa, INa,b, ICa,b), are plotted as functions of log(INaCa). Note that a Hopf bifurcation to abolish BP activity occurred only in the unmodified system (HB at INaCa = 0.0056).

We also examined the contributions of I_NaCa or [Ca²⁺]_i dynamics to the EP instability and structural stability to depolarizing or hyperpolarizing loads of the BP system by determining the unstable regions in the I_bias-V₀ curve during decreases in I_NaCa or [Ca²⁺]_i transients. As shown in Fig. 10 (left), inhibition of I_NaCa significantly shrunk the unstable regions, though a codimension two saddle-node bifurcation did not occur. When [Ca²⁺]_i or f_Ca was fixed at the control value, however, the unstable regions did not shrink during I_NaCa reduction (Fig. 10, middle, right). The influences of eliminating [Ca²⁺]_i transients or blocking SR Ca²⁺ release on the unstable regions were very small.

View larger version (30K): [in this window] [in a new window]   FIGURE 10  Changes in the unstable regions of the Ibias-V0 curve during INaCa suppression in the IK1-removed systems. Displacements of the bifurcation points (H1–2, SN1–2), as well as the control V0 at Ibias = 0 (EP1), in the Ibias-V0 curve with decreasing INaCa are shown for the unmodified system (left) and the modified system with fCa fixed at the value for [Ca2+]i = 0.5 µM (middle) or [Ca2+]i itself fixed at 0.5 µM (right). The log(INaCa) value was reduced from 0 to –3 at an interval of 0.01. As in Fig. 7, values of V0 and Ibias at the bifurcation points are plotted on both the potential (top) and current (bottom) coordinates as functions of log(INaCa). In the unmodified system, the path of EP1 intersected the locus of H2 at INaCa = 0.0056 (labeled HB).

	DISCUSSION

TOP ABSTRACT INTRODUCTION THEORY AND METHODS RESULTS DISCUSSION APPENDIX 1: MODEL EQUATIONS APPENDIX 2: DETERMINATION OF... ACKNOWLEDGEMENTS REFERENCES

It should be noted that the EP with the most depolarized potential (EP1) is always unstable, not stabilized during I_K1 suppression (Fig. 5). The saddle-node bifurcation may also be yielded by increasing inward currents such as I_Na,b or I_Ca,b. Nevertheless, increasing the inward currents would not lead to BP generation, because the positively shifted EP1 becomes stable via a Hopf bifurcation, resulting in an arrest at depolarized potentials. Another approach to create a BP cell may be overexpression of the hyperpolarization-activated current (I_h), which has been reported to induce BP activity in canine atrial or Purkinje myocytes (52,56). However, incorporating I_h did not cause a saddle-node bifurcation nor BP generation in the HVM model with normal I_K1, although accelerating BP generation during I_K1 suppression via shifting the saddle-node bifurcation point toward higher g_K1 values (data not shown). This discrepancy is at least in part due to the difference in I_K1 density, which is much higher in ventricular cells than in atrial or Purkinje cells (56