- Open Access
Critical scale of propagation influences dynamics of waves in a model of excitable medium
Nonlinear Biomedical Physicsvolume 3, Article number: 4 (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 .
Basic equations that describe a class of exactly solvable models for excitable media have been defined in [29, 30]. Here we introduce a modification of this analytical model by adjusting the excitation threshold, v r , in response to changes in frequency of external pacing. We will consider the model in dimensionless form:
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.
The right-hand term B in Eq. (3) responds to a stepwise evolution of pacing period, T m . Similar to the experimental findings described in [31, 32], stepwise changes in B result in smooth exponential transition of the excitation threshold from one steady-state plateau to another. For the sake of simplicity, a steady-state value of excitation threshold at each mth plateau, , was chosen to be linearly dependent on the corresponding pacing interval, T m :
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
The system of Eqs. (1)–(3) was solved numerically on a short cable of 150 grid points with spatial and temporal grid intervals of Δx = 0.13 and Δt = 7.2 × 10-4, respectively. Periodic wavetrains were produced by stimulating the cable with a square wave at the left end using the function P (x, t) = X(x)Y (t), defined above, where A = 1.2, δ1 = 2Δx, δ2 = 15Δx, and T s = 103Δt. The model parameters λ, ϵ, and ζ were equal to 0.4, 0.1, and 1.2, respectively, for all simulations. Numerical solutions were computed using a second-order explicit-difference scheme with no-flux boundary conditions . A typical solution as shown in Figure 1 depicts the propagation of a single pulse between two successive pacing stimuli. The length of the cable is approximately equal to the width of the fully developed pulse to reflect the physiologically relevant relative dimensions of the heart and a propagating excitation wave in a wide range of heart rates and excitation thresholds.
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
The system of equations was initially studied with Eq. (3) replaced by its asymptotic form, which is equivalent to τ = ∞. The cable was stimulated periodically in a stepwise manner with one percent decrements between each pacing plateau (40 consecutive pacing periods) over the range T m = 61 to 28. At the end of each plateau, T h , T f , and v min had all reached steady state values that were used to compose the steady state restitution curve and determine the minimal stable action potential duration, , for a given excitation threshold. Steady state restitution curves computed for two different excitation thresholds are shown in Figure 2.
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].
Figure 3 shows a superposition of the steady state restitution curve with five transitions between different steady-state values of T h computed with 5% changes in T m starting with T m = 37.7 (β = 5 × 10-3, α = 0.375). At longer BCLs, the S1–S2 curves are shallower than the steady-state restitution curves, and all of the transient constant BCL responses after the first response to the change in stimulation rate lie on a straight line with negative slope, which is referred to as the "BCL-line" [21, 22, 28]. At shorter BCLs, the slope of the S1–S2 curve increases. When the slope of the S1–S2 curve (insert A) is greater than the slope of the steady-state curve, the response (dots near the end of steady-state curve) to a decrease in stimulation interval evolves as a series of oscillations that damp to a new steady-state value.
Figure 4 summarizes some typical T h responses to an abrupt change in pacing rate. In inserts A and B, the response is a gradual adaptation to a subsequent steady-state during which time the BCL is constant. The duration of this adaptation is directly controlled by the constant, τ. When τ is small, the bifurcation occurs almost immediately after five stimulation periods (Figure 4, insert C). When τ is large, the bifurcation is correspondingly delayed as shown in the insert D. This might be understood as "sliding" into the unstable area after an abrupt change in pacing rate due to the gradual increase of v r , rather than simply jumping into it as for the cable with constant excitation threshold.
Our simulations showed that complexity of wave dynamics for the system (1)–(3) was also determined by the difference between the values of v min and critical excitation threshold for a stable solitary pulse . For ϵ = 0.1 at the end of the restitution curve for the case shown in Figure 3 the value of v min was 10% greater than the value of . Under these circumstances, we found that damping oscillations preceded conduction block type 3:2 responses (Figure 3, insert B). When the difference between v min and decreased to 3% of , perturbations of T m transformed into persistent T h alternans (Table 1).
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.
Murray J: Mathematical Biology. 1993, New York: Springer, 2
Mikhailov A: Foundations of Synergetics I. Distributed Active Systems. 1994, New York: Springer, 2 revised
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.
Noble D: A modification of the Hodgkin-Huxley equations applicable to Purkinje fibre action and pacemaker potentials. Journal of Physiology. 1962, 160: 317-352.
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.
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.
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.
Karma A: Electrical alternans and spiral wave breakup in cardiac tissue. Chaos. 1994, 4: 461-472. 10.1063/1.166024.
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.
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.
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.
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.
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.
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.
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.
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.
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.
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.
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.
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.
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.
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.
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.
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.
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.
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.
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.
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.
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.
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.
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.
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.
Ravens U, Wettwer E: Electrophysiological aspects of changes in heart rate. Basic Res Cardiol. 1998, 93 Suppl 1: 60-65. 10.1007/s003950050220.
Rushton WAH: Initiation of the Propagated Disturbance. Proceedings of the Royal Society of London. 1937, 124: 210-243. 10.1098/rspb.1937.0083.
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.
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.
JS is Chief Science Officer of Mediwave Star Technology, Inc.
JS, CD, VV, AS, and VP conceived and designed the theory and numerical studies. CD, VV, and AS performed the numerical simulations. JS, CD, VV, AS, and VP analyzed the results. JS, CD, VV, AS, and VP compiled and the manuscript.