Computers in Biology and Medicine 37 (2007) 8 – 20 www.intl.elsevierhealth.com/journals/cobm
Phase response characteristics of sinoatrial node cells D.G. Tsalikakisa, d , H.G. Zhangb , D.I. Fotiadisa, c, d,∗ , G.P. Kremmydasa , Ł.K. Michalisd, e a Unit of Medical Technology and Intelligent Information Systems, Department of Computer Science, University of Ioannina, GR 45110 Ioannina, Greece b University of Manchester, Manchester, UK c Biomedical Research Institute-FORTH, GR 45110 Ioannina, Greece d Michaelidion Cardiology Center, GR 45110 Ioannina, Greece e Department of Cardiology, Medical School, University of Ioannina, GR 45110 Ioannina, Greece
Received 29 December 2004; accepted 20 September 2005
Abstract In this work, the dynamic response of the sinoatrial node (SAN), the natural pacemaker of the heart, to short external stimuli is investigated using the Zhang et al. model. The model equations are solved twice for the central cell and for the peripheral cell. A short current pulse is applied to reset the spontaneous rhythmic activity of the single sinoatrial node cell. Depending on the stimulus timing either a delay or an advance in the occurrence of next action potential is produced. This resetting behavior is quantified in terms of phase transition curves (PTCs) for short electrical current pulses of varying amplitude which span the whole period. For low stimulus amplitudes the transition from advance to delay is smooth, while at higher amplitudes abrupt changes and discontinuities are observed in PTCs. Such discontinuities reveal critical stimuli, the application of which can result in annihilation of activity in central SAN cells. The detailed analysis of the ionic mechanisms involved in its resetting behavior of sinoatrial node cell models provides new insight into the dynamics and physiology of excitation of the sinoatrial node of the heart. 䉷 2005 Elsevier Ltd. All rights reserved. Keywords: Heart; Sinoatrial node; Mathematical models; Cardiac models; Phase resetting; Three-dimensional phase transition curves; Regional differences
1. Introduction Auto-rhythmic excitation in the cardiac muscle is driven by the sinoatrial node (SAN), the natural pacemaker of the heart. Under physiological or pathological conditions during clinical interventions, such us defibrillation of cardiac muscle or accidental situations like electrocution, the pacemaker of the heart is subjected to external stimulation [1]. Experimental [2,3] and numerical [4–6] studies indicate that the application of short current pulses to sinoatrial pacemaker cells result in changes in the cell’s cycle length which depend on both the phase and the amplitude of the stimulus [3,7–31]. In most cases, the applied perturbation does not affect the amplitude or the configuration of the pacemaker cell action potential. The magnitude and direction of the phase shift resulting ∗ Corresponding author. Unit of Medical Technology and Intelligent Information Systems, Department of Computer Science, University of Ioannina, GR 45110 Ioannina, Greece. Tel.: +30 6510 98893; fax: +30 6510 98889. E-mail address:
[email protected] (D.I. Fotiadis).
0010-4825/$ - see front matter 䉷 2005 Elsevier Ltd. All rights reserved. doi:10.1016/j.compbiomed.2005.09.011
from the alteration of the beat-to-beat interval depends on the timing as well as on the strength of the stimulus. In terms of the resetting effect, the total charge injected into the cell seems to be the crucial factor. Thus, any combinations of various pulse amplitudes and durations, which result in the injection of the same amount of current, are equivalent [5]. The response of biological oscillators to external perturbations has been extensively studied using both theoretical and experimental approaches (reviewed in [32]). Using topological arguments, Winfree [33] greatly contributed to the qualitative understanding of phase resetting responses in biological oscillators and pointed out that, in a spatially extended system like cardiac tissue, such responses play a crucial role in the initiation of arrhythmias. Stimuli of critical amplitude and phase might lead to annihilation of normal rhythmic activity and the initiation of complex spiral wave arrhythmias. Mathematical models which describe the electrical activity of SAN have been presented by various research groups [1]. The first models produced by Yanagihara et al. [34] and Noble and Noble [35], where the basis for models presented later [36,37].
D.G. Tsalikakis et al. / Computers in Biology and Medicine 37 (2007) 8 – 20
The analysis of phase response characteristics provides new insight into the dynamics of the mathematical models, of the electrical activities of SAN, and an experimentally verifiable test for the accurate reconstruction of SAN tissue dynamics. Simulation of the spatial variation in the electrical properties of SAN cells and their phase response profiles is necessary for reconstruction of the complex dynamics of the SAN tissue. Zhang et al. [38] recently presented one-dimensional model which describes all of the above. In this work, we investigate and characterize the effects of external stimulation on central and peripheral SAN cells using the Zhang et al. model [38]. The equation of the model along with the appropriate initial conditions are solved using the Runge–Kutta method. Short (0.5 ms) depolarizing and hyperpolarizing electrical current pulses of varying amplitude and of timing spanning the whole period are applied and the phase transition curves (PTC) are obtained [22]. Three-dimensional PTC [4], relating old phase and stimulus amplitude to the new phase after stimulation, are generated in order to locate critical stimuli and singularities in the models. Numerical simulations for the electrical activity of the transitional cells lying between the center and the periphery of the SAN [38] are carried out. Simulations are conducted when specific ionic currents vary. The application of a critical depolarizing stimulus (about 0.4 nA) during the late repolarization phase of action potential resulted in annihilation of the activity in central SAN cells, revealing the existence of a stable singularity in the corresponding model configuration [22]. Detailed analysis of the phase response characteristics of the peripheral cells failed to show a similar singularity and annihilation of normal activity. This difference in phase response behavior of central and peripheral cells is portrayed in the structure of the corresponding PTC.
9
where V is the membrane potential, t is the time, C is the membrane capacitance of SAN and Itot = (INa + ICa,L + ICa,T + Ito + Isus + IK,r + IK,s + If + Ib,Na + Ib,Ca + Ib,K + INaCa + Ip ). The transition from central to peripheral SAN activity was modelled as in [38], using a scaling factor which depends on the distance of each cell from the center of the SAN. The set of differential Eqs. (1)–(16) given in Tables 2–9 is solved using a fourth order Runge–Kutta method with time step 0.1 ms. Further details for the formulation of the equations described in Tables 2–9 and the experimental background behind the equations are sited in [38]. The set of differential Eqs. (1)–(16) are solved twice: first the set of equations is solved for the central cell with initial conditions shown in the left column of Table 9 and second for the peripheral cell with initial conditions shown in the right column in Table 9. The model equations for central, peripheral and transitional cell configurations are integrated until a stable limit cycle is obtained. The minimum diastolic potential is determined on the stable limit cycle and the vector of dynamic variables at this point is recorded and used as the initial condition in the subsequent simulations. For each of the configurations examined, a single 0.5 ms depolarizing or hyperpolarizing current stimulus of a given amplitude is applied at several time steps during the whole period of integration. The starting reference point is selected at the membrane potential level of (−10 mV) during the depolarization phase of the action potential and at each subsequent run the stimulus timing increases by 1 ms. The stimulation current amplitude varies from 0 to 7.5 nA. For a given stimulus amplitude, the corresponding PTC are obtained [22], for the whole range of stimulus amplitudes, three-dimensional PTC are also obtained [4].
2. Methods 3. Results We assume that we apply a short pulse at the SAN and we obtain the PTC using the model proposed by Zhang et al. [38]. The spatial effects are ignored and the problem is reduced to a time dependent problem which is described by the set of ordinary differential equations. 2.1. Model equations The electrical activity of central and peripheral SAN cells is simulated using the model proposed by Zhang et al. [38] which was developed using experimental data. The model includes formulations for the ionic currents (INa , ICa,L , ICa,T , Ito ), sensitive sustained current (Isus , IK,r , IK,s , If ) and background currents (Ib,Na , Ib,Ca , Ib,K , Ip , INaCa ). Table 1 shows the parameters entering the model. The total ionic current Itot appearing in the differential equation which models the membrane potential V is comprised of the currents: INa , ICa,L , Ia,T , Ito , Isus , IK,r , IK,s , If , Ib,Na , Ib,Ca , Ib,K , INaCa and Ip . The total membrane potential is modelled as dV −1 = Itot , dt C
(1)
The normal electrical activity of central, peripheral, and transitional SAN cells is shown in Fig. 1. The membrane potential V of central (dark line—in Fig. 1(a)) and peripheral (dark line—in Fig. 1(c)) SAN cells is plotted vs time. The gray lines indicate the electrical activity of the transitional cells. The duration of the simulation is 1 s. The transition between central and peripheral activity is modelled using a scaling factor (ranging from 0 to 1) which depends on the distance of the cell from the center of the SAN (see Eqs. (80) and (81) in Ref. [38] and Ref. [39]). The scaling factor is used in order to calculate the capacitance and the ionic conductances of the transitional cells given their corresponding values at the center and the periphery of the SAN. For the transitional cell activity shown in Fig. 1 the scaling factor values (0.04 in Fig. 1(a) and 0.13 in plot Fig. 1(c)) were selected to indicate the point where a sharp change in membrane potential configuration is evident. Thus, in Fig. 1(c) the transitional cell appears to have the characteristic spike-anddome configuration of peripheral cells. The total ionic current Itot is plotted vs membrane potential V for each of the corresponding simulations Fig. 1(b) and (d).
25
0.05
0
0.00
-25
-0.05
-50
-0.10
-75
-0.15 0.0
0.2
0.4
(a)
0.6
0.8
1.0
t
-75 (b)
-50
-25 V
0
25
30
0.35
0
0.00
-30
-0.35
-60
-0.70
Itot
V
Itot
D.G. Tsalikakis et al. / Computers in Biology and Medicine 37 (2007) 8 – 20
V
10
-1.05
-90 0.0
0.2
0.4
0.6
0.8
1.0
t
(c)
-90
-60
(d)
-30
0
30
V
Fig. 1. Normal electrical activity in central (a,b), peripheral (c,d) and transitional sinoatrial node cells (gray lines in a–d). In figures (b) and (d) the total ionic current Itot is plotted vs the membrane potential V shown in figures (a) and (c). Figures (a) and (c) indicate the variations of membrane potential V with time t.
0.05
40
T0
T 0.00
V
Itot
0
-0.05
-40
-0.10
-80 0.0 (a)
0.4
0.8 t
1.2
1.6
-75 (b)
-50
-25 V
0
25
Fig. 2. (a) Single pulse perturbation protocol for the central sinoatrial node cell. A depolarizing stimulus of 0.45 nA is applied with a time delay after the detection (at V = −10 mV during the depolarization phase) of the second action potential. The gray trace corresponds to the unperturbed activity. (b) The total ionic current Itot vs the membrane potential V .
The difference in Itot amplitude between central and peripheral cells is mainly due to the contribution of INa which is not present in central cells. The Itot –V plots provide a means to visualize the effects of external perturbation on the SAN cell activity (see Fig. 2) and can be considered as a representation of the complete multi-dimensional phase space of the model. The single pulse perturbation protocol used is shown in Fig. 2. The gray solid line represents the normal unperturbed electrical activity of central SAN cell with a period T0 =334 ms. The onset of the action potential is marked at V = −10 mV (depolarization phase). At a time interval after the onset of the action potential a short (0.5 ms) depolarizing current pulse is
applied resulting in early depolarization of the cell membrane (dark line). The phase of the stimulus is /T0 while the new phase of the perturbed oscillator is = + 1 − T /T0 where T is the period of the perturbed cycle [22]. The electrical activity resumes immediately after the first perturbed cycle as is evident from the Itot –V plot (Fig. 2(b)) in which both simulations are shown (gray line corresponds to the unperturbed limit cycle oscillation of the central SAN cell and the black line corresponds to the limit cycle oscillation of the SAN). The individual ionic currents comprising Itot in both the central and peripheral SAN cell models are shown in Fig. 3(a), (b). Current traces of 1 s of spontaneous activity are plotted.
D.G. Tsalikakis et al. / Computers in Biology and Medicine 37 (2007) 8 – 20
1 Itot
Itot
0.04
-5 0.0 ICa.L
ICa.L
-0.10 0.02 -0.10
-1.4 0.00 ICa.T
-0.004
-0.08
0.001
0.012 IK.s
IK.s
ICa.T
0.000
-0.002
0.004
0.004 INaCa
INaCa
0.000
-0.002
-0.010
INaK
0.026
0.005
0.012
0.00
0.00 Ib
Ib
INaCa
0.008
-0.04
0.001
0.03 If
If
-0.01
-0.002 0.0 (a)
11
0.2
0.4
0.6
0.8
-0.05 0.0
1.0
t
(b)
0.2
0.4
0.6
0.8
1.0
t
Fig. 3. (a) Ionic currents for central sinoatrial node cell model; (b) Ionic currents for peripheral sinoatrial node cell model.
A comparison of the two plots reveals the contribution of each ionic current in Itot . The sodium current INa is the principal component of Itot in peripheral while it is absent in central SAN cells. In the central cell model, the maximum upstroke velocity of the action potential is 6 V s−1 , the minimum diastolic potential is −68 mV and the maximum systolic potential is 22 mV, while in the peripheral cell model the maximum upstroke velocity of the action potential is 86 V s−1 , the minimum diastolic potential is −81 mV and the maximum systolic potential is 25 mV. For each one of the experiments shown in Fig. 1(a) and (c), we repeat those with from 0 to T0 (correspondingly 0 < < 1) and current stimulus amplitudes from 0 to 7.5 nA. The results are shown in Fig. 4 where representative PTC are depicted. Fig. 4(a)–(d) correspond to low amplitude stimuli (−0.3 nA) while Fig. 4(e)–(g),(j) traces correspond to high amplitude stimuli. Fig. 4(a)–(d), (j) correspond to central and peripheral cells while the Fig. 4(b), (c), (f), (g) correspond to the transitional cell configurations shown in Fig. 1 (gray lines in Fig. 1(a) and (c), respectively). The advance or delay in the occurrence of the first action potential after perturbation for a set of experiments for which the stimulus phase varies from 0 to 1 can be visualized by plotting the ratio of the corresponding period T1 (during which
the stimulus was applied) to the period T0 of the unperturbed cycle over the stimulus phase . Such plots were obtained for the central and peripheral SAN (Fig. 5). Low amplitude stimulation shown in Fig. 5(a), (b), in both the central and peripheral SAN cells, produces only minimal advances (T1 /T0 < 1) or delays (T1 /T0 > 1) for most of the cycle excluding a small region around = 0.4 where a significant advance is observed. At this region the membrane is in the absolute refractory period and the stimulus amplitude is sufficient to initiate an action potential and thus to shorten T1 . Increasing the stimulation amplitude extends the region where an advance is observed. In central cells the delay in the occurrence of the action potential at intermediate stimulus levels observed close to = 0.6 is associated with the approximation of a singularity in the phase space. Around this singularity T1 becomes close to 2 × T0 (Fig. 5(c)) while for 0.4 < < 0.6 approaching T1 /T0 = 2, there are missing points which correspond to T1 /T0 → ∞. The PTC shown in Fig. 4 indicate that, in central, peripheral, and transitional SAN cells, both Type 0 (Fig. 4(a)–(d)) and Type 1 (Fig. 4(e)–(j)) phase resetting is obtained. Low amplitude stimuli produce continuous PTC since the effect of the external perturbation on the phase of the limit cycle oscillator is less pronounced. The discontinuities of the PTC for stronger
12
D.G. Tsalikakis et al. / Computers in Biology and Medicine 37 (2007) 8 – 20
φ'
1
φ'
1
0
0 0
(a)
1
φ
0 (b)
1
1
φ'
φ'
1
0
0 1
0
φ
(c)
1
0
φ
(d)
1
φ'
1
φ' 0
0
0
1
φ
(e)
0
1
φ
(f) 1
φ'
φ'
1
0 (g)
φ
0
φ
0
1 (h)
0
φ
1
Fig. 4. Phase transition curves for central ((a), (e)), peripheral ((d), (j)), and transitional cells ((b), (c), (f), (g)).
stimuli (Type 1 phase resetting) are indicative of the existence of a singularity in the limit cycle. Each of the PTC in Fig. 4 represents a set of experiments like the one shown in Fig. 2 where the stimulus amplitude is constant while the interval varies. A set of such curves for various stimulus amplitudes produces a three-dimensional phase transition surface [4] relating both and stimulus amplitude to the new phase of the perturbed oscillator. Such three-
dimensional phase transition plots allow further investigation of the phase resetting behavior of an oscillator. Moreover, such PTC can help in finding the exact stimulus phase and amplitude combination which brings the oscillator close to the singularity range. The three-dimensional phase transition plots for the central, transitional and peripheral SAN cells are shown in Fig. 6 as contour plots with the new phase is color-coded
T1 T0
D.G. Tsalikakis et al. / Computers in Biology and Medicine 37 (2007) 8 – 20 2.0
2.0
1.8
1.8
1.6
1.6
1.4
1.4
1.2
1.2
1.0
1.0
0.8
0.8
0.6
0.6
0.4
0.4 0.2
0.2 0.0
0.2
0.4
T1 T0
0.6
0.8
1.0
(a)
0.0
0.2
0.4
2.0
2.0
1.8
1.8
1.6
1.6
1.4
1.4
1.2
1.2
1.0
1.0
0.8
0.8
0.6
0.6
0.4
0.4
0.6
0.8
1.0
0.6
0.8
1.0
0.6
0.8
1.0
(b)
0.2
0.2 0.0
0.2
0.4
0.6
0.8
1.0
(c)
T1 T0
13
0.0
0.2
0.4
(d)
2.0
2.0
1.8
1.8
1.6
1.6
1.4
1.4
1.2
1.2
1.0
1.0
0.8
0.8
0.6
0.6
0.4
0.4
0.2
0.2 0.0
0.2
0.4
0.6
0.8
1.0
(e)
0.0
0.2
0.4
(f) 2.0 1.8 1.6 1.4
T1 T0
1.2 1.0 0.8 0.6 0.4 0.2 0.0
0.2
(g)
0.4
0.6
0.8
1.0
Fig. 5. Alteration of the timing of the first action potential after stimulation is plotted T1 /T0 vs for various stimulation amplitudes. (a) Peripheral cell at −0.5 nA, (c) peripheral cell at −2.5 nA, (e) peripheral cell at −5 nA, (g) peripheral cell at −6.5 nA and (b) central cell at −0.2 nA, (d) central cell at −0.4 nA, (f) central cell at −0.6 nA (T0,central = 320 ms, T0,peripheral = 160 ms).
on a gray scale from 0 to 1. The ordering of the plots is similar to that of Fig. 4 and the corresponding PTC of Fig. 4 can be considered as horizontal “slices” of the data at the corresponding current stimulus amplitude (Istim ). The difference in the stimulus amplitude for the corresponding plots is due
to the difference in the amplitude of the corresponding total membrane current (a stronger stimulus is required to produce Type 1 phase resetting behavior of Fig. 1). In the three-dimensional phase transition plots shown in Fig. 6, there exist regions where the new phase is very dense
14
D.G. Tsalikakis et al. / Computers in Biology and Medicine 37 (2007) 8 – 20
Central
Peripheral 0.0 0.1
-0.3
-1.5
0.2 0.3
-0.6
-3.0
Istim
0.4 0.5
-0.9
-4.5
-1.2
-6.0
-1.5
-7.5
0.6 0.7 0.8 0.9
0.0
0.2
0.4
0.6
0.8
1.0
(a)
0.0
0.2
0.4
0.6
0.8
1.0
1.0
(b)
Transitional
Transitional 0.0 0.1
-0.4 -1
0.2
-0.8
0.3
Istim
-2
0.4
-1.2
0.5 0.6
-3
-1.6
0.7 -2.0
0.8
-4
0.9 -2.4
1.0
-5 0.0 (c)
0.2
0.4
0.6
0.8
1.0
0.0
0.2
0.4
0.6
0.8
1.0
(d)
Fig. 6. Three-dimensional contour plots of the new phase vs old phase and amplitude: for (a) central, (b) peripheral, and (c), (d) transitional SAN cell.
60
V
30
0
-30
-60 0.0
0.5
1.0 t
1.5
2.0
Fig. 7. A critical stimulus of −0.4 nA applied at t = 0.53 s resulting in permanent annihilation of sinoatrial node repetitive activity.
(e.g. for = 0.55 and Istim = −0.4 nA). Such regions indicate the existence of a singularity and should be further investigated by fine-tuning of both interval and Istim amplitude.
Further experimentation with finer values for and Istim within those regions revealed the existence of a stable singularity for the central SAN cell. In the simulation shown in Fig. 7 a
D.G. Tsalikakis et al. / Computers in Biology and Medicine 37 (2007) 8 – 20
critical stimulus of −0.4 nA is applied at t = 0.53 s resulting in permanent annihilation of the normal repetitive activity. Transitional and the peripheral cells do not exhibit similar behavior. The application of stimulus of appropriate amplitude and phase (e.g. = 0.4 and Istim = −5.4) results in a series of low amplitude oscillations of the membrane potential with normal limit cycle behavior resuming after 2–3 cycles. 4. Discussion Investigation of the dynamic behavior of biological oscillators and their response to external perturbation is of great importance in biological research since biological oscillations are involved in many vital processes in living systems [32,33]. Understanding the dynamic response of the SAN to external perturbations is important in elucidating its behavior under normal and pathological conditions [33]. Studies in the past have concentrated on the elucidation of the phase resetting dynamics of the SAN using biophysical ionic models [15,29,5] but did not address the issue of regional differences between central and peripheral node cells (the corresponding ionic models described central cell behavior). The SAN is a complex structure and regional differences in the electrical activity of regional cells seem to be important in the behavior of the natural pacemaker of the heart [39,40,30]. The Zhang et al. model [38] used in this study accounts for regional differences and can be used to model the electrical activity of central, peripheral, and transitional cells. In this work, the analysis of phase resetting characteristics of each type of SAN cell by means of three-dimensional phase transition plots, revealed important differences in their response to external electrical stimuli (see Fig. 6). Although all types of cells display both Type 0 and Type 1 phase resetting behavior, the stimulus amplitude at which Type 1 phase resetting behavior first appears varies. This is a result of the regional differences in Itot amplitude (see Fig. 1): Itot is higher in the periphery of SAN (mainly due to the dominance of INa ) and thus a higher amplitude of external stimulus is required in order to significantly affect the normal limit cycle behavior. We have produced Itot –V plots to provide a means to visualize the effects of external perturbation on the SAN cell electrical activity (Fig. 2). Our simulations show the existence of a critical stimulus amplitude-phase combination which is capable of annihilating the electrical activity of the central SAN cells (Fig. 7). So far, there is no reference in the literature describing and proofing the existence of such a combination in Zhang model. It is believed that such critical stimuli are involved in the generation of certain cardiac arrhythmias [33] and gives new insights into the investigation of cardiac models. In a complex structure like the SAN, regional differences and spatiotemporal interactions could play an important role in its pathophysiological response to external perturbations and should be further investigated. This must be included in our next set of experiments targeting one-dimensional cells. However, the use of the finite element method might help in the analysis of two- and three-dimensional structures. The comparison of phase resetting characteristics of other models will be included too.
15
Investigation of the effects of increases of various concentration on annihilation of SAN activity must be further studied.
5. Summary In the present work, we use the model proposed from Zhang et al. [38] and we examine the dynamic response of the SAN to external stimulus. We apply a short current to reset the spontaneous rhythmic activity of the single SAN. This application of the external stimulus results in the prolongation of the normal activity and a new phase is obtained. This resetting behavior is quantified in terms of phase transition curves (PTCs) for short electrical current pulses of varying amplitude which span the whole period. Discontinuities, as we show in the figures, reveal critical stimulus. The application of this critical stimulus can result in annihilation of activity on SAN cells. Detailed analysis of the ionic mechanisms which involved in the resetting behavior could be useful and provide new insight into the physiology of the excitation in SAN cell. Our work, is a step forward to a future work on higher dimensionality model were we believe the resetting technic can be more significant.
Acknowledgements This work is partially funded by the Greek Secretariat for Research and Technology PENED-01E511 to D.G. Tsalikakis.
Appendix The parameters used in the model are given in Table 1. The set of differential equations (1)–(16) given in Tables 2–9 was solved using a fourth order Runge–Kutta method with time step 0.1 ms. Table 1 Parameters of the model 4-AP AM APD CL Cm a (x), C s (x) Cm m
dL , dT dNa,Ca dV /dtmax Da , Ds
EK,s ENa , ECa , EK ECa,L , ECa,T F FK,r FNa
4-Aminopyridine Atrial muscle Action potential duration Spontaneous cycle length Cell capacitance Capacitance of atrial muscle cell or SA node cell in one-dimensional model of intact SA node at distance x from center of SA node Activation variables for ICa,L and ICa,T Denominator constant for INa,Ca Maximum upstroke velocity of action potential Diffusion coefficient between atrial muscle cells or SA node cells in one-dimensional model of the intact SA node Reversal potential for IK,s Equilibrium potentials for Na+ , Ca2+ , and K+ Reversal potentials for ICa,L and ICa,T Faraday’s constant Fraction of activation of IK,r that occurs slowly Fraction of inactivation of INa that occurs slowly
16
D.G. Tsalikakis et al. / Computers in Biology and Medicine 37 (2007) 8 – 20 Table 1 (continued)
Table 1 (continued) fL , fT gp , gc g a (x), g s (x)
gNa gCa,L , gCa,T gto , gsus gK,r , gK,s gf,Na , gf,K gb,Na , gb,Ca , gb,K h1 , h2 h INa ICa,L , iCa,T Ito , Isus IK,r , IK,s IK If If,Na , If,K Ib,Na , Ib,Ca , Ib,K INaCa Ip Ip Itot a (x), I s (x) Itot tot
Ist IK,ACh IK,AP kNa,Ca Km,Na , Km,K L
Ls m MDP n∞ pa pa,f , pa,s pi Q10 r R q SA node, SAN t T TOP V V a, V s
V a (x), V s (x)
Inactivation variables for ICa,L and ICa,T Conductance of a current in peripheral or central SA node cell models Conductance of a current in atrial muscle cell or SA node cell in one-dimensional model of intact SA node at distance x from center of SA node Conductance of INa Conductance of ICa,L and gCa,T Conductance of Ito and Isus Conductance of IK,r and IK,s Conductance of Na+ and K+ components of if Conductance of Ib,Na , Ib,Ca and Ib,K Fast and slow inactivation variables for INa Net fractional availability of INa TTX-sensitive Na+ current L- and T-type Ca2+ currents Transient and sustained components of 4-APsensitive current Rapid and slow delayed rectifying K+ currents Sum of IK,r and IK,s Hyperpolarization-activated current Na+ and K+ components of If Background Na+ , Ca2+ , and K+ currents Na+ /Ca2+ exchanger current Na+ –K+ pump current Maximum Ip Total ionic current in a cell Total ionic current in atrial muscle cell or SA node cell in one-dimensional model of intact SA node at distance x from center of SA node Sustained current ACh-activated K+ current ATP-sensitive K+ current Scaling factor for INa,Ca Dissociation constants for Na+ and K+ activation of Ip Length of string of SA node and atrial tissue in one-dimensional model of intact SA node Length of string of SA node tissue in one-dimensional model of intact SA node Activation variable for iNa Maximum diastolic potential Steady-state value of n General activation variable for IK,r Fast and slow activation variables forIK,r Inactivation variable for IK,r Fractional change in a variable with a 10 ◦ C increase in temperature Activation variable for Ito Universal gas constant Inactivation variable for Ito Sinoatrial node Time Absolute temperature Takeoff potential Membrane potential Membrane potential of atrial muscle cell or SA node cell in one-dimensional model of intact SA node Membrane potential of atrial muscle cell or SA node cell in one-dimensional model of intact SA node at distance x from center of SA node
x xs y z [Na+ ]i , [Ca2+ ]i an bn
NaCa n
Distance from center of SA node in one-dimensional model of intact SA node Activation variable for IK,s Activation variable of If Valency of ion Intracellular Na+ , Ca2+ , concentrations Voltage-dependent opening rate constant of n Voltage-dependent closing rate constant of n Position of Erying rate theory energy barrier controlling voltage dependence of INaCa Time constant of n Space constant
Table 2 Formulation of the equations of Na+ current (INa ) The sodium current INa is given by the equation INa = gNa (m)3 h[Na+ ]o
(F )2.0 (e(V −ENa )F /RT − 1.0) V RT (eVF/RT − 1.0)
Sodiumcurrent m gate dm (m∞ − m) = dt m m∞ =
1.0 (1.0 + e−V /5.46 )
(2) 1.0/3.0
m =
0.0006247 (0.832e−0.335(V +56.7) + 0.627e0.082(V +65.01) ) +0.00004)
Sodium current h gate h = ((1.0 − FNa )h1 + FNa h2 ) (h1,∞ − h1 ) dh1 = dt h1
(3)
(h2,∞ − h2 ) dh2 = dt h 2
(4)
h1,∞ =
1.0 (1.0 + e(V +66.1)/6.4 )
h2,∞ = h1,∞ 0.000003171e−0.2815(V +17.11) h1 = + 0.0005977 (1.0 + 0.003732e−0.3426(V +37.76) )
h2 =
0.00000003186e−0.6219(V +18.8) + 0.003556 (1.0 + 0.00007189e−0.6683(V +34.07) )
D.G. Tsalikakis et al. / Computers in Biology and Medicine 37 (2007) 8 – 20 Table 3 Formulation of the equations of L-type Ca2+ current, ICa,L
Table 4 (continued)
L-type Ca channel ICa,L 0.006 ICa,L = gCa,L fL dL + (V − ECa,L ) −(V +14.1)/6.0 (1.0 + e )
d,T =
17
1.0 (d,T + d,T )
dT,∞ =
1.0 (1.0 + e−(V +37.0)/6.8 )
L-type Ca channel d gate (dL,∞ − dL ) ddL = dt d,L d,L = −14.19
(V + 35.0) 42.45V − −0.208V (e − 1.0) (e−(V +35.0)/2.5 − 1.0)
(5)
(fT,∞ − fT ) dfT = dt f,T
(8)
f,T = 15.3e−(V +71.7)/83.3
5.71(V − 5.0) d,L = 0.4(V −5.0) (e − 1.0)
f,T = 15.0e(V +71.7)/15.38
1.0 d,L = (d,L + d,L ) dL,∞ =
T-type Ca channel f gate
1.0 (f,T + f,T )
f,T =
1.0 (1.0 + e−(V +23.1)/6.0 )
fT,∞ =
1.0 (1.0 + e(V +71.0)/9.0 )
L-type Ca channel f gate dfL (fL,∞ − fL ) = dt f,L
f,L =
(6)
Table 5 Formulation of the equations for sensitive currents (Ito and Isus ) and delayed rectifying current (IK,r )
3.12(V + 28.0) (e(V +28.0)/4.0 − 1.0)
Ito
f,L = f,L =
25.0
Ito = gto qr(V − EK )
(1.0 + e−(V +28.0)/4.0 )
Isus Isus = gsus r(V − EK )
1.0 (f,L + f,L )
fL,∞ =
q gate (q∞ − q) dq = dt q
1.0 (1.0 + e(V +45.0)/5.0 )
q∞ =
(9)
1.0 (1.0 + e(V +59.37)/13.1 )
0.06517 + 0.000024e0.1(V +50.93) 0.57e−0.08(V +49.0)
Table 4 Formulation of the equations of T-type Ca2+ current (ICa,T )
q = 0.0101 +
T-type Ca channel ICa,T
4-AP sensitive currents r gate
ICa,T = gCa,T dT fT (V − ECa,T )
(r∞ − r) dr = dt r
T-type Ca channel d gate (dT,∞ − dT ) ddT = dt d,T
d,T = 1068.0e(V +26.3)/30.0
r∞ = (7)
(10)
1.0 (1.0 + e−(V −10.93)/19.7 )
r = 0.00298 +
d,T = 1068.0e−(V +26.3)/30.0
References
0.01559 + 0.369e−0.12(V +23.84) )
(1.037e0.09(V +30.61)
18
D.G. Tsalikakis et al. / Computers in Biology and Medicine 37 (2007) 8 – 20 Table 6 (continued)
Table 5 (continued) Rapid delayed rectifying potassium current IK,r
14.0 (1.0 + e−(V −40.0)/9.0 )
x s =
IK,r = gK,r pa pi (V − EK ) pa = ((1.0 − FK,r )pa,f + FK,r pa,s )
xs = e−V /45.0 Rapid delayed rectifying potassium current pa,f gate (pa,f,∞ − pa,f ) dpa,f = dt pa,f
pa,f,∞ =
pa,f =
(11) xs∞ =
1.0 (1.0 + e−(V +14.2)/10.6 )
1.0 (xs + xs )
xs =
1.0 (37.2e(V −9.0)/15.9 + 0.96e−(V −9.0)/22.5 )
xs (xs + xs )
Hyperpolarization activated current If If = (If,Na + If,K )
Rapid delayed rectifying potassium current pa,s gate If,Na
(pa,s,∞ − pa,s ) dpa,s = dt pa,s
(12)
pa,s,∞ = pa,f,∞
(13)
pa,s =
1.0 (4.2e(V −9.0)/17.0 + 0.15e−(V −9.0)/21.6 )
If,Na = gf,Na y(V − ENa )
If,K If,K = gf,K y(V − EK )
Hyperpolarization activated current y gate (y∞ − y) dy = dt y
Rapid delayed rectifying potassium pi gate dpi (pi,∞ − pi ) = dt p i pi,∞ =
(14)
1.0 (1.0 + e(V +18.6)/10.1 )
y = e−(V +78.91)/26.62
y = e(V +75.13)/21.25
y∞ =
Table 6 Formulation of the equations for slow delayed rectifying current (IK,s ) and hyperpolarization current (If ) Slow delayed rectifying potassium current IK,s IK,s = gK,s (xs )2.0 (V − EK,s ) Slowdelayed rectifying potassium current xs gate dxs (xs∞ − xs ) = dt xs
(15)
y =
y (y + y )
1.0 (y + y )
(16)
D.G. Tsalikakis et al. / Computers in Biology and Medicine 37 (2007) 8 – 20
Table 7 Formulation of the equations for background, pump and exchanger currents Sodium background current Ib,Na Ib,Na = gb,Na (V − ENa ) Potassium background current Ib,K Ib,K = gb,K (V − EK ) Calcium background current Ib,Ca Ib,Ca = gb,Ca (V − ECa ) Sodium calcium exchanger INaCa INaCa = KNaCa (([Na+ ]i )3.0 [Ca2+ ]o e0.03743V NaCa −([Na+ ]o )3.0 [Ca2+ ]i e0.0374V (NaCa −1.0) ) (1.0+dNaCa ([Ca2+ ]i ([Na+ ]o )3.0 +[Ca2+ ]o ([Na+ ]i )3.0 ))
Sodium potassium pump Ip 2.0 3.0 [Na+ ]i [K+ ]o Ip =Ip,max (Km,K + [K+ ]o ) (Km,Na + [Na+ ]i ) 1.6 × (1.5 + e−(V +60.0)/40.0 )
Table 8 Formulation of the equations for reversal and equilibrium potentials Reversal and equilibrium potentials RT [Na+ ]o ln F [Na+ ]i
ENa = EK =
RT [K+ ]o ln + F [K ]i
ECa =
[Ca2+ ]o RT ln 2.0F [Ca2+ ]i
EK,s =
RT ([K+ ]o + 0.12[Na+ ]o ) ln F ([K+ ]i + 0.12[Na+ ]i )
Table 9 Initial conditions for central and peripheral cell Variable
V m h1 h2 dL fL dT fT y r q xs pa,f pa,s pii
Value Central
Peripheral
−55.4434 0.158763 0.0365616 0.0515223 0.00328396 0.853578 0.0622561 0.129637 0.0230307 0.0334583 0.288777 0.137793 0.29045 0.561665 0.973981
−78.4021 0.0531058 0.341077 0.0207796 0.000109311 0.97761 0.00226369 0.278358 0.0215392 0.0110545 0.435506 0.0950833 0.313397 0.3783 0.991546
19
[1] in: D. Zipes, J. Jalife (Eds.), Cardiac Electrophysiology—From Cell to Bedside, second ed., W.B. Saunders, Philadelphia, 1995. [2] W. Zeng, L. Glass, A. Shrier, The topology of phase response curves induced by single and paired stimuli in spontaneously oscillating chick heart cell aggregates, J. Biol. Rhythms. 7 (2) (1992) 89–104. [3] J.M. Shumaker, J.W. Clark, W.R. Giles, A model of the phase-sensitivity of the pacemaking cell in the bullfrog heart, J. Theor. Biol. 151 (2) (1991) 193–230. [4] G.P. Kremmydas, A.V. Holden, A. Bezerianos, T. Bountis, Representation of sino-atrial node dynamics by circle maps, Internat. J. Bifurcation Chaos 6 (10) (1996) 1799–1805. [5] A.C. Coster, B.G. Celler, Phase response of model sinoatrial node cells, Ann. Biomed. Eng. 31 (3) (2003) 271–283. [6] S. Abramovich-Sivan, S. Akselrod, A prc based model of a pacemaker cell: effect of vagal activity and investigation of the respiratory sinus arrhythmia, J. Theor. Biol. 192 (2) (1998) 219–234. [7] E.E. Verheijck, R. Wilders, R.W. Joyner, D.A. Golod, R. Kumar, H.J. Jongsma, L.N. Bouman, A.C.G. van Ginneken, Pacemaker synchronization of electrically coupled rabbit sinoatrial node cells, J. Gen. Physiol. 111 (1) (1998) 95–112. [8] S.S. Demir, J.W. Clark, W.R. Giles, Parasympathetic modulation of sinoatrial node pacemaker activity in rabbit heart: a unifying model, Am. J. Physiol. Heart Circ. Physiol. 276 (6) (1999) H2221–H2244. [9] R.F. Gilmour, Phase-resetting of circus movement reentry in isolated canine purkinje-muscle preparations, Circulation 76 (4) (1987) 112. [10] V.S. Reiner, C. Antzelevitch, Phase resetting and annihilation in a mathematical-model of sinus node, Am. J. Physiol. 249 (6) (1985) 1143 –1153. [11] T.R. Chay, Y.S. Lee, Phase resetting and bifurcation in the ventricular myocardium, Biophys. J. 47 (5) (1985) 641–651. [12] J.M.B. Anumonwo, M. Delmar, A. Vinet, D.C. Michaels, J. Jalife, Phase resetting and entrainment of pacemaker activity in single sinus nodal cells, Circ. Res. 68 (4) (1991) 1138–1153. [13] A.M. Kunysz, A.A. Munk, A. Shrier, Phase resetting and dynamics in isolated atrioventricular nodal cell clusters, Chaos 5 (1) (1995) 184–192. [14] M.R. Guevara, A. Shrier, Phase resetting in a model of cardiac purkinjefiber, Biophys. J. 52 (2) (1987) 165–175. [15] M.R. Guevara, H.J. Jongsma, Phase resetting in a model of sinoatrial nodal membrane—ionic and topological aspects, Am. J. Physiol. 258 (3) (1990) H734–H747. [16] J.R. Clay, R.M. Brochu, A. Shrier, Phase resetting of embryonic chick atrial heart cell aggregates—experiment and theory, Biophys. J. 58 (3) (1990) 609–621. [17] M.R. Guevara, A. Shrier, L. Glass, Phase resetting of spontaneously beating embryonic ventricular heart cell aggregates, Am. J. Physiol. 251 (6) (1986) H1298–H1305. [18] J.R. Clay, M.R. Guevara, A. Shrier, Phase resetting of the rhythmic activity of embryonic heart cell aggregates—experiment and theory, Biophys. J. 45 (4) (1984) 699–714. [19] R.M. Luceri, A. Castellanos, D.S. Kayden, F. Moleiro, R.G. Trohman, J.T. Horgan, R.J. Myerburg, Phase resetting of ventricular parasystolic rhythms produced by extrinsic perturbations, Pace Pacing Clin. Electrophys. 7 (3) (1984) 461. [20] J. Jalife, A.J. Hamilton, V.R. Lamanna, G.K. Moe, Effects of current flow on pacemaker activity of the isolated kitten sinoatrial node, Am. J. Physiol. 238 (3) (1980) H307–H316. [21] J. Jalife, V.A.J. Slenter, J.J. Salata, D.C. Michaels, Dynamic vagal control of pacemaker activity in the mammalian sinoatrial node, Circ. Res. 52 (6) (1983) 642–656. [22] L. Glass, A.T. Winfree, Discontinuities in phase-resetting experiments, Am. J. Physiol. 246 (2) (1984) R251–R258. [23] D.C. Michaels, D.R. Chialvo, E.P. Matyas, J. Jalife, Chaotic activity in a mathematical model of the vagally driven sinoatrial node, Circ. Res. 65 (5) (1989) 1350–1360. [24] V.I. Arnold, Cardiac arrhythmias and circle mappings, Chaos 1 (1) (1991) 20–24. [25] L. Glass, Cardiac arrhythmias and circle maps—a classical problem, Chaos 1 (1) (1991) 13–19.
20
D.G. Tsalikakis et al. / Computers in Biology and Medicine 37 (2007) 8 – 20
[26] A.T. Winfree, Biological rhythms and the behavior of populations of coupled oscillators, J. Theor. Biol. 16 (1) (1967) 15–42. [27] W.P.M. Vanmeerwijk, G. Debruin, A.C.G. Vanginneken, J. Vanhartevelt, H.J. Jongsma, E.W. Kruyt, S.S. Scott, D.L. Ypey, Phase resetting properties of cardiac-pacemaker cells, J. Gen. Physiol. 83 (4) (1984) 613–629. [28] M. Kawato, R. Suzuki, Biological oscillators can be stopped–topological study of a phase response care, Biol. Cybern. 30 (4) (1978) 241–248. [29] D.C. Michaels, E.P. Matyas, J. Jalife, A mathematical model of the effects of acetylcholine pulses on sinoatrial pacemaker activity, Circ. Res. 55 (1) (1984) 89–101. [30] M.R. Boyett, H. Honjo, I. Kodama, The sinoatrial node, a heterogeneous pacemaker structure, Cardiovasc. Res. 47 (4) (2000) 658–687. [31] M.R. Boyett, H. Zhang, A. Garny, A.V. Holden, Control of the pacemaker activity of the sinoatrial node by intracellular Ca2+ . experiments and modelling, Philos. Trans. Soc. London Ser. A Math. Phys. Eng. Sci. 359 (1783) (2001) 1091–1110. [32] A.T. Winfree, The Geometry of Biological Time, Springer, New York, 1980. [33] A.T. Winfree, When Time Breaks Down, Princeton University Press, Princeton, 1987. [34] K. Yanagihara, A. Noma, H. Irisawa, Reconstruction of sinoatrial node pacemaker potential based on the voltage-clamp experiments, Jpn. J. Physiol. 30 (1980) 841–857. [35] D. Noble, S. Noble, A model of sino-atrial node electrical activity based on a modification of the difrancesco-noble equations, Proc. R. Soc. London B Biol. Sci. 222 (1984) 295–304. [36] S. Demir, J. Clark, C. Murphey, W. Giles, A mathematical model of a rabbit sinoartial node cell, Am. J. Physiol. Cell Physiol. 266 (3) (1994) 832–852. [37] S. Dokos, B. Celler, N. Lovell, Modification of difrancesco-noble equations to simulate the effects of vagal stimulation on vivo mammalian sinoatrial node electrical activity, Ann. Biomed. Eng. 21 (1999) 321–335. [38] H. Zhang, A.V. Holden, I. Kodama, H. Honjo, M. Lei, T. Varghese, M.R. Boyett, Mathematical models of action potentials in the periphery and center of the rabbit sinoatrial node, Am. J. Physiol. Heart Circ. Physiol. 279 (1) (2000) H397–H588. [39] A. Garny, P. Kohl, P.J. Hunter, M.R. Boyett, D. Noble, One-dimensional rabbit sinoatrial node models: benefits and limitations, J. Cardiovasc. Electrophysiol. 14 (s10) (2003) S121. [40] H. Zhang, A.V. Holden, M.R. Boyett, Gradient model versus mosaic model of the sinoatrial node, Circulation 103 (4) (2001) 584–588.
Dimitrios G. Tsalikakis was born in Thessaloniki, Greece, in 1977. He received the Diploma degree in Mathematics from the University of Ioannina and he is a PhD candidate in the Department of Cardiology since 2002. His research interests include cardiac modeling, monophasic action potential analysis, automated systems analysis, virtual instrumentation, and scientific computing. He is member of the Unit of Medical Technology and Intelligent Information Systems. Henggui Zhang, born in 1964 in China, educated in Physics (BSc, 1985), Computer Science (MSc, 1985–1988), Non-linear Science (PhD, 1988–1990), and Biomedical Science (PhD, 1991–1994). Since 1991, he has been working as a research assistant (1991–1994, Leeds, funding from MRC), postdoctoral research fellow (1994–1995, JHU, funding from NSF, USA; 1995–1996, Leeds; funding from Welcome Trust; 1996–2000, Leeds; funding from the British Heart Foundation), and then a senior research fellow (2000–2001, Leeds, funding from British Heart Foundation). In October 2001, he moved to UMIST to take up a lectureship. Currently he is a senior lecturer in Biological Physics Group, School of Physics and Astronomy, the University of Manchester. Dimitrios I. Fotiadis was born in Ioannina, Greece, in 1961. He received the Diploma degree in Chemical Engineering from the National Technical University of Athens and the PhD degree in Chemical Engineering from the University of Minnesota, USA. Since 1995, he has been with the Department of Computer Science, University of Ioannina, Greece, where he is currently an Associate Professor. He is the director of the Unit of Medical Technology and Intelligent Information Systems. His research interests include biomedical technology, biomechanics, scientific computing and intelligent information systems. Georgios P. Kremmydas, born in 1972 in Argostoli, Kefalonia Island, Greece. He received his BSc (Hon) Degree in Physiology in 1994, University of Leeds, UK and a PhD in 2000, Department of Biomedical Sciences, University of Leeds, UK. He is currently the Director of Research Technology and Development Department of L.U.M.C. of Kefalonia and Ithaki and a visiting researcher of the Computer Science Department, University of Ioannina, Greece. Lampros K. Michalis was born in Arta, Greece, in 1960. He received the MD degree with distinction from the Medical School, University of Athens, Greece, in 1984. Since 1995, he has been with the Medical School, University of Ioannina, Greece, where he is currently a Professor in Cardiology. He is in charge of the coronary care unit and the catheter laboratory of the University Hospital of Medical School, Ioannina, Greece. His research interests focus on bioengineering, interventional cardiology, and electrophysiology.