Critical scale of propagation influences dynamics of waves in a model of excitable medium
© Starobin et al; licensee BioMed Central Ltd. 2009
Received: 23 March 2009
Accepted: 09 July 2009
Published: 09 July 2009
Duration and speed of propagation of the pulse are essential factors for stability of excitation waves. We explore the propagation of excitation waves resulting from periodic stimulation of an excitable cable to determine the minimal stable pulse duration in a rate-dependent modification of a Chernyak-Starobin-Cohen reaction-diffusion model.
Various pacing rate dependent features of wave propagation were studied computationally and analytically. We demonstrated that the complexity of responses to stimulation and evolution of these responses from stable propagation to propagation block and alternans was determined by the proximity between the minimal level of the recovery variable and the critical excitation threshold for a stable solitary pulse.
These results suggest that critical propagation of excitation waves determines conditions for transition to unstable rhythms in a way similar to unstable cardiac rhythms. Established conditions were suitably accurate regardless of rate dependent features and the magnitude of the slopes of restitution curves.
Studies of unstable waves in reaction-difusion systems are the subject of significant theoretical and practical importance, particularly with respect to the analysis of biological excitable media such as nerve and cardiac tissue [1, 2]. Pioneering studies suggested that propagation of waves in these media may be governed by a few fundamental parameters which influence formation of the excitation wavefront [3, 4]. One such parameter, ϵ = τ f /τ r ≪ 1, reflects the two order-of-magnitude difference in time constants between the fast excitation, τ f , and slow recovery, τ r , processes, while another one, excitation threshold, v r , is linked to the critical level of excitation necessary to initiate the propagation of an excitation pulse through the medium.
Although examination of instabilities based on this approach would possess an obvious advantage of being theoretically analyzable, attention has primarily focused on experimentally motivated methods. It has been shown in tissue preparations that the duration of excitation, known as action potential duration, and the speeds of excitation wavefronts and wavebacks depend on previous stimulation periods and refractory (diastolic) intervals [5, 6]. The analysis of such dependences, known as restitution curves, has been introduced as a method for evaluation of stability of excitation waves in biological excitable media [7–17].
Numerical simulations also confirmed these findings and demonstrated that T h and T f can be fitted to specific restitution curves (called a restitution portrait) using time dependent gate variables in various reaction (ionic) models [18–25]. Several recent studies implemented an approach based on a projection of different ionic models to discrete maps with sufficient rate-dependence (memory) and multi-component equivalent restitution portraits [20–27]. However, it was later confirmed that none of the stability criteria derived from these complex discrete maps predicted the onset of alternans in tissue studies . Computational experiments described in  suggest that such deficiency is more likely to result from lack of analysis of propagation in the medium as a whole rather than from the lack of particular details of the ionic models of individual excitable cells.
To address this deficiency we introduce an alternative approach based on the analysis of a modified Chernyak-Starobin-Cohen (CSC) reaction-diffusion model which accounts for effects of wave propagation. In contrast to the complexity of the majority of existing reaction-diffusion models, the CSC model employs just a few fundamental parameters and is analytically solvable. It offers a rigorous criterion for determining the stability of excitation wave propagation in an excitable cable [29, 30].
In this paper, we study the stability of excitation waves by examining their propagation in a one-dimensional CSC medium, and modify the model accordingly to incorporate the memory linked pacing rate driven adjustments of excitation threshold, v r . This modification was motivated by direct experimental measurements of v r that demonstrated that exponential-like evolution of excitation threshold takes place over the course of multiple pacing cycles following stepwise changes in pacing rate in guinea pig and human ventricular muscle [31–33]. Such an approach allows us to establish the condition for stability of propagation of excitation waves using analytical and numerical solutions of the reaction-diffusion model while incorporating into the model experimentally observed action potential restitution dynamics.
We demonstrate that in a wide range of excitation thresholds, v r , and pacing rates, T m , regardless of the particular values of slopes of restitution curves and rates of adaptation of v r (more or less memory), the loss of stability of waves in an excitable cable of finite length is determined by proximity of the minimal level of the recovery variable, v, at the foot of the action potential and the critical excitation threshold of the solitary pulse computed analytically in .
u(x, t) and v(x, t) are a membrane potential and slow recovery current, respectively. λ, ϵ, ζ, and τ are the model parameters, where τ-1 < ϵ. The scaling of the system is described below.
The pacing function, P (x, t), is defined as a product of two functions, X(x) and Y (t). Each is composed of the Heaviside step function, Θ (x), as follows: X(x) = A [Θ(x-δ1) - Θ(x-δ2)] and Y (t) = Θ (t k ) - Θ (t k + T s ), where A and δ2 - δ1 are the amplitude and width of the pulse, respectively, T s is the pulse duration, and t k = tN(m-1) + [k - N (m - 1)]T m are the instants of time when stimuli are delivered. T m is the pacing period, N represents the number of stimuli at each pacing interval plateau, the index m denotes each pacing plateau, m = 1, ..., M, and k is an integer in the range k = 1, ..., N M . The number of stimuli N is the same for all plateaus. Overall during the course of the protocol, T m progressively decreases to the minimal value when stable propagation of the wave is still possible.
Here α and β are positive parameters that determine the amplitude of change of B m between two consecutive pacing plateaus.
The scale of u is the maximum steady-state action potential amplitude U0, the scale of v is given by σ NA U0, and the time scale is C m /σ NA , where σ NA corresponds to the maximum sodium conductance and C m is the membrane capacitance. The characteristic length scale is given by , where D is the diffusion coefficient. The small parameter ϵ ≪ 1, is equal to C m /(t K σ NA ) and ζ = σ K /σ NA where σ K and τ K correspond to the maximum potassium conductance and potassium current time constant, respectively.
Results and discussion
In order to quantify the dynamics of the system in response to perturbations in T m (Eqs. (1)–(3)), we computed the action potential duration (T h ), the diastolic interval (T f ) at each pacing period, and the minimal value of the recovery variable at the foot of the propagating pulse, v min (Figure 1). The values of v min are always greater than v r due to incomplete medium recovery from repeated pacing, and v min increases monotonically with decreasing T m . The action potential duration was defined as the interval of time when u > v at a specified node, x0. Accordingly, the refractory period, T f , was defined as the interval of time when u <v. In order to analyze steady-state duration of the developed pulse, we measured T h at the center of the cable x0 = 10.
Stability of propagation for constant excitation threshold
Insert A portrays the analytical dependence of T h on excitation threshold v r for a stable solitary pulse . It also shows the critical level of excitation threshold ( = 0.359, vertical dashed line) beyond which no stable solitary pulse exists. For the repeatedly paced cable, the numerically determined minimal steady state values, , at the ends of the two steady-state restitution curves are depicted by the two circles (open and filled) in the lower right corner of insert A, which are to the right of the vertical dashed line v = (Figure 2). The location of both circles demonstrates that the loss of stability for the sequence of pulses occurs when v min is close to the critical excitation threshold for the stable solitary pulse . The character of the loss of stability at the ends of restitution curves depends on the margin between v min and which decreases for greater values of excitation threshold v r (Figure 2, Insert A, filled circle).
It was observed that for = 0.1 complexity of block type behavior for higher values of v r evolved from 2:1 to 3:2 responses, followed by T h alternans when the difference between v min and was less than 3%. The evolution from stable propagation to propagation block and T h alternans is illustrated by inserts B, C, and D, respectively. Insert B characterizes stable adaptation of T h from one steady-state to another. The adaptation is rapid, and as a result the steady-state and S1–S2 curves nearly coincide. At the end of the restitution curve with lower values of v r stable action potential propagation transforms into a block type N:M (N>M) response shown in insert C (upper curve, v r = 0.25). For higher values v r when v min is closer to , stable propagation transforms to T h alternans shown in insert D (lower curve, v r = 0.31).
Similar trends were observed for higher and lower values of ϵ. For ϵ = 0.06 2:1 conduction blocks transformed directly to alternans regardless of the difference between and v min . For ϵ = 0.14 however, there was a progression from 2:1 responses to more complex 3:2 patterns which culminated in alternans as the difference between and v min decreased from approximately sixteen to one percent.
Obtained results confirm that irrespective of the restitution slopes, alternans do not appear until, at the end of the restitution curve, the value of v min exceeds by some small margin (Figure 2). Under these conditions alternans appear for slopes smaller than one (S = 0.44, Figure 2 lower curve), and conversely, stable propagation occurs for slopes greater than one (S = 1.94, Figure 2 upper curve).
Finally, although very useful for stability analysis, the original CSC model with constant excitation threshold does not allow reconstruction of transitional restitution behavior near a change in pacing rate. Appropriate model enhancements are described in the next section.
Stability of propagation with rate-dependent v r
In order to reproduce experimentally observed restitution curves, the system of equations (1)–(2) was upgraded with Eq. (3), which allowed us to incorporate effects of memory . We will show that in this system the criterion for the appearance of action potential duration alternans will be the same as in the memory-less case detailed in the previous section describing the cable with constant excitation threshold.
The pacing protocol described previously was used again to construct the steady-state restitution curve and determine the minimal stable action potential durations. In addition, a similar stepwise protocol with larger 5% decrements in T m was used. A conventional S1–S2 restitution protocol, in which the cable was stimulated for forty consecutive periods with a given T m (S1) followed by delivery of a premature S2 stimulus, was used to compute S1–S2 curves around given steady-state values of T h .
A gradual phase of constant BCL adaptation followed the immediate S1–S2 response to an abrupt change in pacing rate, and the T h transients had a duration of 5–50 stimulation periods depending on the adaptation constant, τ . Pictured in the restitution domain, this combination of S1–S2 and constant BCL responses formed the triangles typical of experimental and computational studies of ventricular action potential [21, 22, 24].
Evolution of Wave Dynamics.
ϵ = 0.06, = 0.40
ϵ = 0.10, = 0.36
ϵ = 0.14, = 0.33
Similar to simulations with constant v r , for rate dependent excitation threshold the complexity of transition from blocks to alternans also increased for higher values of parameter ϵ. For ϵ = 0.06 transformation from 2:1 blocks to alternans was almost abrupt. However, for higher values of ϵ as the difference between v min and further decreased, the transition to alternans was more complex and evolved through stages with higher 3:2 and 4:3 complexity of propagation blocks (Table 1).
Previous attempts to determine conditions for action potential duration alternans have focused on the analysis of the slopes of restitution curves in various theoretical and experimental models. However, these efforts did not result in a consistent theoretical criterion for prediction of action potential duration instabilities [22, 24, 27, 28]. In part this happened because stability was examined solely in regard to the magnitudes of particular restitution slopes rather than with the analysis of the conditions for a critically propagating pulse [13, 34, 35].
In this paper, we have demonstrated that stable propagation of excitation waves in a paced cable occurs until the minimal level of the recovery variable in front of the rising action potential wave reaches a value that is greater than the critical excitation threshold for the stable solitary pulse by some small margin. This margin decreases for higher v r , which agrees with previous analytical findings for the CSC model approximating the critical speeds of an infinite wavetrain . At the limit of the critically stable pulses, our numerical system revealed two types of unstable behavior. When the minimal value of recovery variable at the foot of the propagating wave approached the critical level of excitation threshold beyond which no stable solitary pulse was able to propagate, we observed conduction blocks of increasing complexity followed by alternans. Alternans developed slower for larger values of the adaptation constant, τ.
It was demonstrated that these conditions were valid with or without rate-dependent features and regardless of the magnitude of the slopes of restitution curves.
This research was funded by Mediwave Star Technology, Inc. and was partially supported by the University of North Carolina at Greensboro. We are grateful to Lanty L. Smith and Thomas R. Sloan for their continuous support. We would also like to thank David Schaeffer, Wanda Krassowska, and Daniel Gauthier for helpful discussions and critical reviews.
- Murray J: Mathematical Biology. 1993, New York: Springer, 2View ArticleMATHGoogle Scholar
- Mikhailov A: Foundations of Synergetics I. Distributed Active Systems. 1994, New York: Springer, 2 revisedView ArticleMATHGoogle Scholar
- Krinsky V, Kokoz Y: Analysis of equations of excitable membranes – I. Reduction of the Hodgkin-Huxley equations to a second order system. Biofizika. 1973, 18: 506-511.Google Scholar
- Noble D: A modification of the Hodgkin-Huxley equations applicable to Purkinje fibre action and pacemaker potentials. Journal of Physiology. 1962, 160: 317-352.View ArticleGoogle Scholar
- Chialvo D, Michaels D, Jalife J: Supernormal excitability as a mechanism of chaotic dynamics of activation in cardiac Purkinje fibers. Circulation Research. 1990, 66: 525-545.View ArticleGoogle Scholar
- Gilmour R, Otani N, Watanabe M: Memory and complex dynamics in canine cardiac Purkinje fibers. Am J Physiol. 1997, 272 (4 Pt 2): H1826-H1832.Google Scholar
- Cain J, Tolkacheva E, Schaeffer D, Gauthier D: Rate-dependent propagation of cardiac action potentials in a one-dimensional fiber. Phys Rev E Stat Nonlin Soft Matter Phys. 2004, 70 (6 Pt 1): 061906-10.1103/PhysRevE.70.061906.View ArticleGoogle Scholar
- Karma A: Electrical alternans and spiral wave breakup in cardiac tissue. Chaos. 1994, 4: 461-472. 10.1063/1.166024.View ArticleADSGoogle Scholar
- Fenton F, Karma A: Vortex dynamics in three-dimensional continuous myocardium with fiber rotation: filament instability and fibrillation. Chaos. 1998, 8: 20-47. 10.1063/1.166311.View ArticleADSMATHGoogle Scholar
- Garfinkel A, Kim Y, Voroshilovsky O, Qu Z, Kil J, Lee M, Karagueuzian H, Weiss J, Chen P: Preventing ventricular fibrillation by flattening cardiac restitution. Proceedings of the National Academy of Sciences. 2000, 97: 6061-6066. 10.1073/pnas.090492697.View ArticleADSGoogle Scholar
- Qu Z, Garfinkel A, Chen P, Weiss J: Mechanisms of discordant alternans and induction of reentry in simulated cardiac tissue. Circulation. 2000, 102: 1664-1670.View ArticleGoogle Scholar
- Schwartz I, Triandaf I, Starobin J, Chernyak Y: Origin of quasiperiodic dynamics in excitable media. Phys Rev E Stat Phys Plasmas Fluids Relat Interdiscip Topics. 2000, 61 (6 Pt B): 7208-7211. 10.1103/PhysRevE.61.7208.Google Scholar
- Chernyak Y, Starobin J: Characteristic and critical excitation length scales in 1D and 2D simulations of reentrant cardiac arrhythmias using simple two-variable models. Critical Reviews in Biomedical Engineering. 1999, 27 (3–5): 359-414.Google Scholar
- Boyett M, Jewell B: A study of the factors responsible for rate-dependent shortening of the action potential in mammalian ventricular muscle. Journal of Physiology. 1978, 285: 359-380.View ArticleGoogle Scholar
- Elharrar V, Surawicz B: Cycle length effect on restitution of action potential duration in dog cardiac fibers. Am J Physiol. 1983, 244 (6): H782-H792.Google Scholar
- Franz M, Swerdlow C, Liem L, Schaefer J: Cycle length dependence of human action potential duration in vivo. Journal of Clinical Investigation. 1988, 82: 972-979. 10.1172/JCI113706.View ArticleGoogle Scholar
- Arnold L, Page J, Attwell D, Cannell M, Eisner D: The dependence on heart rate of the human ventricular action potential duration. Cardiovascular Research. 1982, 16: 547-551. 10.1093/cvr/16.10.547.View ArticleGoogle Scholar
- Fenton F, Evans S, Hastings H: Memory in an excitable medium: a mechanism for spiral wave breakup in the low excitable limit. Physical Review Letters. 1999, 83: 3964-3967. 10.1103/PhysRevLett.83.3964.View ArticleADSGoogle Scholar
- Cherry E, Fenton F: Suppression of alternans and conduction blocks despite steep APD restitution: electrotonic, memory, and conduction velocity restitution effects. Am J Physiol Heart Circ Physiol. 2004, 286 (6): H2332-H2341. 10.1152/ajpheart.00747.2003.View ArticleGoogle Scholar
- Tolkacheva E, Schaeffer D, Gauthier D, Mitchell C: Analysis of the Fenton-Karma model through an approximation by a one-dimensional map. Chaos. 2002, 12: 1034-1042. 10.1063/1.1515170.View ArticleADSMathSciNetMATHGoogle Scholar
- Tolkacheva E, Schaeffer D, Gauthier D, Krassowska W: Condition for alternans and stability of the 1:1 response pattern in a 'memory' model of paced cardiac dynamics. Phys Rev E Stat Nonlin Soft Matter Phys. 2003, 67 (3 Pt 1): 031904-10.1103/PhysRevE.67.031904.View ArticleADSGoogle Scholar
- Kalb S, Tolkacheva E, Schaeffer D, Gauthier D, Krassowska W: Restitution in mapping models with an arbitrary amount of memory. Chaos. 2005, 15 (2): 023701-10.1063/1.1876912.View ArticleADSMathSciNetGoogle Scholar
- Mitchell C, Schaeffer D: A two-current model for the dynamics of cardiac membrane. Bulletin of Mathematical Biology. 2003, 65: 767-793. 10.1016/S0092-8240(03)00041-7.View ArticleGoogle Scholar
- Schaeffer D, Cain J, Gauthier D, Kalb S, Oliver R, Tolkacheva E, Ying W, Krassowska W: An ionically based mapping model with memory for cardiac restitution. Bulletin of Mathematical Biology. 2007, 69: 459-482. 10.1007/s11538-006-9116-6.View ArticleMathSciNetMATHGoogle Scholar
- Shiferaw Y, Sato D, Karma A: Coupled dynamics of voltage and calcium in paced cardiac cells. Phys Rev E Stat Nonlin Soft Matter Phys. 2005, 71 (2 Pt 1): 021903-10.1103/PhysRevE.71.021903.View ArticleADSGoogle Scholar
- Otani N, Gilmour R: A memory model for the electrical properties of local cardiac systems. Journal of Theoretical Biology. 1997, 187: 409-436. 10.1006/jtbi.1997.0447.View ArticleGoogle Scholar
- Banville I, Gray R: Effect of action potential duration and conduction velocity restitution and their spatial dispersion on alternans and the stability of arrhythmias. Journal of Cardiovascular Electrophysiology. 2002, 13: 1141-1149. 10.1046/j.1540-8167.2002.01141.x.View ArticleGoogle Scholar
- Kalb S, Dobrovolny H, Tolkacheva E, Idriss S, Krassowska W: The restitution portrait: A new method for investigating rate-dependent restitution. Journal of Cardiovascular Electrophysiology. 2004, 15: 698-709. 10.1046/j.1540-8167.2004.03550.x.View ArticleGoogle Scholar
- Chernyak Y, Starobin J, Cohen R: Class of exactly solvable models of excitable media. Physical Review Letters. 1998, 80: 5675-5678. 10.1103/PhysRevLett.80.5675.View ArticleADSGoogle Scholar
- Chernyak Y, Starobin J, Cohen R: Where do dispersion curves end? A basic question in theory of excitable media. Physical Review E. 1998, 58: R4108-R4111. 10.1103/PhysRevE.58.R4108.View ArticleADSMathSciNetGoogle Scholar
- Attwell D, Cohen I, Eisner D: The effects of heart rate on the action potential of guinea-pig and human ventricular muscle. Journal of Physiology. 1981, 313: 439-461.View ArticleGoogle Scholar
- Davidenko J, Levi R, Maid G, Elizari M, Rosenbaum M: Rate dependence and supernormality in excitability of guinea pig papillary muscle. American Journal of Physiology – Heart and Circulatory Physiology. 1990, 259: H290-H299.Google Scholar
- Ravens U, Wettwer E: Electrophysiological aspects of changes in heart rate. Basic Res Cardiol. 1998, 93 Suppl 1: 60-65. 10.1007/s003950050220.View ArticleGoogle Scholar
- Rushton WAH: Initiation of the Propagated Disturbance. Proceedings of the Royal Society of London. 1937, 124: 210-243. 10.1098/rspb.1937.0083.View ArticleADSGoogle Scholar
- Noble D: The relation of Rushton's 'liminal length' for excitation to the resting and active conductances of excitable cells. Journal of Physiology. 1972, 226: 573-591.View ArticleGoogle Scholar
This article is published under license to BioMed Central Ltd. This is an Open Access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/2.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.