Article Information

• PDF (408 kb)
• Full Text with Large Figures
• Cited by in Scopus (7)

PubMed

• Articles by Yasutaka Kurata
• Articles by Toshishige Shibamoto

Related Articles

• Regional Difference in Dynamical Property of Sinoatrial Node...

  Regional Difference in Dynamical Property of Sinoatrial Node Pacemaking: Role of Na Channel Current Biophysical Journal, Volume 95, Issue 2, 15 July 2008, Pages 951-977 Yasutaka Kurata, Hiroyuki Matsuda, Ichiro Hisatome and Toshishige Shibamoto Abstract To elucidate the regional differences in sinoatrial node pacemaking mechanisms, we investigated 1), bifurcation structures during current blocks or hyperpolarization of the central and peripheral cells, 2), ionic bases of regional differences in bifurcation structures, and 3), the role of Na channel current () in peripheral cell pacemaking. Bifurcation analyses were performed for mathematical models of the rabbit sinoatrial node central and peripheral cells; equilibrium points, periodic orbits, and their stability were determined as functions of parameters. Structural stability against applications of acetylcholine or electrotonic modulations of the atrium was also evaluated. Blocking L-type Ca channel current () stabilized equilibrium points and abolished pacemaking in both the center and periphery. Critical acetylcholine concentration and gap junction conductance for pacemaker cessation were higher in the periphery than in the center, being dramatically reduced by blocking . Under hyperpolarized conditions, blocking , but not eliminating , abolished peripheral cell pacemaking. These results suggest that 1), is responsible for basal pacemaking in both the central and peripheral cells, 2), the peripheral cell is more robust in withstanding hyperpolarizing loads than the central cell, 3), improves the structural stability to hyperpolarizing loads, and 4), -dependent pacemaking is possible in hyperpolarized peripheral cells. Abstract | Full Text | PDF (3305 kb)

• Dynamics of the Cell Cycle: Checkpoints, Sizers, and Timers

  Dynamics of the Cell Cycle: Checkpoints, Sizers, and Timers Biophysical Journal, Volume 85, Issue 6, 1 December 2003, Pages 3600-3611 Zhilin Qu, W. Robb MacLellan and James N. Weiss Abstract We have developed a generic mathematical model of a cell cycle signaling network in higher eukaryotes that can be used to simulate both the G1/S and G2/M transitions. In our model, the positive feedback facilitated by CDC25 and wee1 causes bistability in cyclin-dependent kinase activity, whereas the negative feedback facilitated by SKP2 or anaphase-promoting-complex turns this bistable behavior into limit cycle behavior. The cell cycle checkpoint is a Hopf bifurcation point. These behaviors are coordinated by growth and division to maintain normal cell cycle and size homeostasis. This model successfully reproduces sizer, timer, and the restriction point features of the eukaryotic cell cycle, in addition to other experimental findings. Abstract | Full Text | PDF (479 kb)

• Evidence for a Novel Bursting Mechanism in Rodent Trigeminal...

  Evidence for a Novel Bursting Mechanism in Rodent Trigeminal Neurons Biophysical Journal, Volume 75, Issue 1, 1 July 1998, Pages 174-182 Christopher A. Del Negro, Chie-Fang Hsiao, Scott H. Chandler and Alan Garfinkel Abstract We investigated bursting behavior in rodent trigeminal neurons. The essential mechanisms operating in the biological systems were determined based on testable predictions of mathematical models. Bursting activity in trigeminal motoneurons is consistent with a traditional mechanism employing a region of negative slope resistance in the steady-state current-voltage relationship (Smith, T. G. 1975. . 253:450–452). However, the bursting dynamics of trigeminal interneurons is inconsistent with the traditional mechanisms, and is far more effectively explained by a new model of bursting that exploits the unique stability properties associated with spike threshold (Baer, S. M., T. Erneux, and J. Rinzel. 1989. . . . 49:55–71). Abstract | Full Text | PDF (248 kb)

• …more

Copyright © 2005 The Biophysical Society. All rights reserved.
Biophysical Journal, Volume 89, Issue 4, 2865-2887, 1 October 2005

doi:10.1529/biophysj.105.060830

Electrophysiology

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
^† Division of Regenerative Medicine and Therapeutics, Tottori University Graduate School of Medical Science, Yonago 683-0826, Japan

Address reprint requests to Yasutaka Kurata, Dept. of Physiology, Kanazawa Medical University, 1-1 Daigaku, Uchinada-machi, Kahoku-gun, Ishikawa 920-0293, Japan.

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,5,6,7,8,9,10,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.

A mathematical model of cardiac myocytes is described as an n-dimensional nonlinear dynamical system, i.e., a set of ordinary differential equations of the form

(1)

(2)

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,19,20. Modifications of the model are summarized in Table 2; details on major modifications are described below. All expressions (Eqs. (3)) and the standard parameter values used are provided in Appendix 1 and Table 5, respectively.

The kinetics of I_Ca,L was described by Eqs. (3) with the activation (d_L), voltage-dependent inactivation (f_L), and Ca²⁺-dependent inactivation (f_Ca) gating variables. Voltage dependences of I_Ca,L gating kinetics in the current (modified) and original human heart models are shown in Fig. 1 (top), together with temperature-corrected experimental values for comparison. As ten Tusscher et al. 16 pointed out, the steady-state activation and inactivation curves (d_L∞, f_L∞) yielded by the original PB formulas were somewhat unusual for unknown reasons (see Fig. 1, left top), resulting in a very large window current. We therefore renewed the expressions for I_Ca,L. For the steady-state activation and inactivation (d_L∞, f_L∞), we used the data of Bénitah et al. 21. The expression of the activation time constant (τ_dL) was adopted from the ten Tusscher et al. model 16. The formulas for the time constants of the voltage-dependent inactivation and recovery (τ_fL) employed by the original models appeared not to fit the experimental data for single HVMs from Bénitah et al. 21 or other articles (see Fig. 1, right top, and also Figure 2E of ten Tusscher et al. 16). Therefore, we originally formulated τ_fL from the data of Bénitah et al. 21, which were corrected for a temperature of 37°C with a Q₁₀ of 2.2 and [Ca²⁺]_o-dependent factor 16. The τ_fL data were fitted to a function similar to that used by Courtemanche et al. 19 for their human atrial model, using a least square minimization procedure. The formula for the Ca²⁺-dependent inactivation (f_Ca,∞) is also the same as used by Courtemanche et al. 19. Maximum I_Ca,L was formulated as a fully selective Ca²⁺ current, with its reversal potential (E_Ca,L) fixed at a constant value of +52.8mV as reported by Bénitah et al. 21 and the maximum conductance (g_Ca,L) set equal to 249.6 pS/pF at 2mM [Ca²⁺]_o.

Display large version of this figure

Figure 1

Kinetics of I_Ca,L. (Top) Voltage dependence of steady-state probabilities (d_L∞, f_L∞) and time constants for I_Ca,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 d_L∞, f_L∞, 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 I_Ca,L (left) and peak I_Ca,L-V relationship (right). Currents were evoked by 300-ms step pulses from a holding potential of −50mV to test potentials ranging from −30 to +50mV in 10mV increments. Simulating the whole-cell perforated-patch recording, [Ca²⁺]_i was not clamped, whereas [Na⁺]_i and [K⁺]_i were fixed at 5 and 140mM, respectively. The experimental data for the peak I_Ca,L-V relation are from Bénitah et al. 21 (open circles), Pelzmann et al. 22 (squares), and Li et al. 23 (hexagons).

The model-generated I_Ca,L during voltage-clamp pulses and peak I_Ca,L-V relationship are depicted in Fig. 1 (bottom). The simulated peak I_Ca,L-V relation is comparable to the experimental data from Bénitah et al. 21, Pelzmann et al. 22, and Li et al. 23.

The kinetics of I_Kr was described by Eqs. (9) 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,27,28,30,31,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 12pS/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 1Hz 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.

Display large version of this figure

Figure 2

Kinetics of I_Kr (top) and I_Ks (bottom). (Left, middle) Voltage dependence of steady-state probabilities and time constants of I_Kr activation (p_a,∞, τ_pa) and inactivation (p_i,∞) as well as I_Ks 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 p_a,∞, τ_pa, and n_∞ are from Li et al. 24 (circles) and Wang et al. 25 (squares). (Right) Computed voltage-clamp records for I_Kr and I_Ks. Currents were elicited by 3-s step pulses from a holding potential of −60mV to test potentials ranging from −50 to +50mV in 10mV increments.

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). 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 36pS/pF so as to yield the APD at 90% repolarization (APD₉₀) during 1Hz pacing of 336±16ms, as reported by Li et al. 24.

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 2ms 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). 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 1Hz. 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.

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.4mM, 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) (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–2Hz) were relatively small (<5mM), being too small to significantly alter the model cell behaviors. For stability and bifurcation analyses, therefore, [K⁺]_i was fixed at 140mM. The degenerate system with the fixed [K⁺]_i has the finite number of EPs and limit cycles, being suitable for bifurcation analyses.

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×10⁻⁶.

The steady-state AP, sarcolemmal currents, and [Ca²⁺]_i transient of the current model paced at 1Hz with the standard parameter values are shown in Figure 3A. 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 1Hz pacing were −84.9mV, 411V/s, and 338ms, respectively. The diastolic [Ca²⁺]_i and peak [Ca²⁺]_i transient in the model HVM paced at 1Hz 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.41mM (identical) for the resting state and 0.59–2.68 and 2.84–2.92mM, respectively, in a steady state during 1Hz 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.

Display large version of this figure

Figure 3

Simulated dynamics of the model HVM. (A) Steady-state behaviors of the AP, underlying sarcolemmal currents, and [Ca²⁺]_i transient. The model HVM was paced at 1Hz with 1-ms stimuli of 80pA/pF. Differential equations (Eqs. (55)) were numerically solved for 20min with the initial conditions appropriate to a resting state, which were determined by Eqs. (65) with [K⁺]_i fixed at 140mM. Model cell behavior after the last stimulus (during the last AP) is depicted. (B) Rate dependence of APD, I_Ca,L, and [Ca²⁺]_i transients. The model HVM was paced at 0.5, 1, and 2Hz with 1-ms stimuli of 80pA/pF. The differential equations were numerically solved for 20min at each pacing rate; model cell behaviors after the last stimulus are depicted.

We also tested the rate dependence of APD (dynamic restitution), as well as that of the I_Ca,L waveform, [Ca²⁺]_i transients and [Na⁺]_i (Figure 3B). The steady-state values of APD₉₀ at 0.5, 1, and 2Hz were 369, 338, and 298ms, respectively, compatible with the averaged experimental values 23,42. I_Ca,L attenuated with increasing the pacing rate as observed in the AP clamp experiment of Li et al. 23. The values of [Na⁺]_i during stimuli at 0 (resting state), 0.25, 0.5, 1, and 2Hz averaged 6.14, 7.47, 8.26, 9.22, and 9.57mM, respectively, comparable to experimental data 43,44.

The peak [Ca²⁺]_i transient predicted by the current model was a little smaller at 2Hz than at 1Hz. In most experiments, however, the peak [Ca²⁺]_i transient and/or developed tension increased as the pacing rate increased up to 2.5Hz 45,46,47. This inconsistency may result from the lack of intracellular modulating factors such as those involved in the rate-dependent potentiation of I_Ca,L48,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.

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) (see Appendix 2 ). Asymptotic stability of the EP was also determined by computing 14 eigenvalues of a 14×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×10⁻³ 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 60s; 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.

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 B11. 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).

Display large version of this figure

Figure 4

Simulated steady-state BP activities (spontaneous APs, ionic currents, and [Ca²⁺]_i dynamics) in the model HVM at g_K1=0.15 (left) and 0 (right). Differential equations (Eqs. (55)) were numerically solved for 20min at each g_K1 with initial conditions appropriate to an EP and a 1-ms stimulus of 1pA/pF for triggering an AP; model cell behaviors during the last 2.5s starting from MDP are depicted. Note the differences in the ordinate scales for individual currents.

An nth-order “autonomous” continuous-time dynamical system is defined by the state equation dx/dt=f(x), where the vector field f, a smooth function of the state variable x=x(t), does not depend on time, but depends only on the state variable x13.

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 t13.

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.

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

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.

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.

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.

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.

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.

We first examined whether I_K1 downregulation leads to BP generation in the model HVM by decreasing g_K1 at an interval of 0.01 (1%). Spontaneous oscillations (BP activity) abruptly appeared when g_K1 was reduced to 0.15 (15%). Fig. 4 shows the simulated BP activity of the I_K1-downregulated HVM with g_K1 of 0.15 or 0, depicting the membrane potential oscillations, underlying sarcolemmal currents, and [Ca²⁺]_i transients in steady state. The values of MDP, POP, APA, CL, and maximum upstroke velocity were determined for the BP oscillations at each g_K1 value. The simulated AP parameters are listed in Table 4, together with the data from the guinea pig ventricle 3 and rabbit SA node 51 models, as well as from real BP systems 1,52, for comparison. The average [Na⁺]_i during the BP oscillations at g_K1 of 0.15 and 0 were 6.55 and 6.42mM, respectively. Decreasing g_K1 caused the decrease of APA with depolarization of MDP, decrease of CL, increase of [Ca²⁺]_i, and decrease of [Na⁺]_i. Consistent with the result from the guinea pig ventricular model 3, I_NaCa was the dominant inward current during phase 4 depolarization in the HVM pacemaking. In contrast, I_Ca,L was very small during the early phase of the depolarization, although predominantly contributing to phase 0 upstroke. I_Na was also very small due to the voltage-dependent inactivation.

Before investigating the dynamical mechanisms of BP generation via bifurcation analyses of the HVM model, we examined the ionic basis of phase 4 depolarization of the BP activity. The contributions of I_Kr and I_Ks deactivation, as well as individual inward currents (I_Ca,L, I_Na, I_NaCa, I_Na,b, I_Ca,b), to the pacemaker depolarization were assessed by the conventional freezing and blocking methods 53,54,55. Phase 4 depolarization of the I_K1-reduced BP system disappeared on clamping I_Kr deactivation (gating variable p_a) at the time of MDP, but not on clamping I_Ks deactivation (gating variable n). Thus, the g_Kr-decay theory for the natural SA node pacemaker 53,54,55 appears to be applicable to the HVM pacemaker as well. Removal of I_Ca,L at the time of MDP little altered the initial phase of pacemaker depolarization, although reducing the rate of the later phase depolarization with cessation of BP activity. Eliminating I_Na did not significantly affect BP oscillations. Furthermore, in the I_K1-removed system, eliminating I_NaCa (the most dominant pacemaker current) at the MDP abolished phase 4 depolarization, leading to cessation of BP activity, whereas the removal of I_Na,b (the second dominant pacemaker current) or I_Ca,b only reduced the depolarization rate.

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.

Display large version of this figure

Figure 5

Bifurcation structure of the model HVM during g_K1 decreases. (Left) Steady-state I/V relations for g_K1=0−1, depicted at an interval of 0.1 with [K⁺]_i fixed at 140mM. The zero-current crossings, corresponding to EPs, are designated by EP1, EP2, and EP3. (Middle) A bifurcation diagram for g_K1 with the steady-state (EP1–3) and stable periodic (MDP, POP) branches (top), and CL plotted as a function of g_K1 (bottom). EPs were determined by the algebraic equations (Eqs. (65)). The model cell dynamics were computed by numerically solving the differential equations (Eqs. (55)) for 10min at each g_K1, which was reduced at an interval of 0.001, with the initial [Na⁺]_i at g_K1=1 set equal to 6.14mM. The saddle-node bifurcation point at which EP2 and EP3 merge together and disappear is located (labeled SNB at g_K1=0.154). (Right) Steady-state [Ca²⁺]_i (top) and [Na⁺]_i (bottom) at EPs as functions of g_K1. The minimum (Min) and maximum (Max) of [Ca²⁺]_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 g_K1.

The steady-state [Na⁺]_i, 6.14mM in the control resting state (at g_K1=1), reduced as g_K1 decreased, reaching 4.64mM just before the saddle-node bifurcation at g_K1=0.154 (see Fig. 5, right bottom). This reduction in [Na⁺]_i during I_K1 suppression was due to the enhancement of Na⁺ efflux via the Na⁺-K⁺ pump, as well as decrease of Na⁺ influx via the Na⁺/Ca²⁺ exchanger and I_Na,b, with resting potential depolarization. To determine how the [Na⁺]_i decrease affects BP generation, we constructed bifurcation diagrams for g_K1 using the reduced system with [Na⁺]_i fixed at the control value of 6.14mM or critical value of 4.64mM (data not shown). The critical g_K1 at which BP activity emerges via a saddle-node bifurcation was larger in the system with the reduced [Na⁺]_i of 4.64mM (0.153) than with the control [Na⁺]_i of 6.14mM (0.117), suggesting that the [Na⁺]_i reduction facilitates BP generation during I_K1 suppression.

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.