Research

Critical scale of propagation influences dynamics of waves in a model of excitable medium

Joseph M Starobin^1,2 , Christopher P Danford² , Vivek Varadarajan² , Andrei J Starobin² and Vladimir N Polotski²

¹  Department of Physics and Astronomy, University of North Carolina at Greensboro, Greensboro, NC, USA

²  Mediwave Star Technology, Inc, Greensboro, NC, USA

author email corresponding author email

Nonlinear Biomedical Physics 2009, 3:4doi:10.1186/1753-4631-3-4

The electronic version of this article is the complete one and can be found online at: http://www.nonlinearbiomedphys.com/content/3/1/4

Received:	23 March 2009
Accepted:	9 July 2009
Published:	9 July 2009

© 2009 Starobin et al; licensee 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.

Abstract

Background

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.

Results

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.

Conclusion

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.

Background

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 [28]. Computational experiments described in [19] 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 [29].

Methods

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:

(1)

(2)

(3)

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-δ₁) - Θ(x-δ₂)] and Y (t) = Θ (t_k) - Θ (t_k + T_s), where A and δ₂ - δ₁ are the amplitude and width of the pulse, respectively, T_s is the pulse duration, and t_k = t_N(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 m^th plateau, , was chosen to be linearly dependent on the corresponding pacing interval, T_m:

(4)

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 U₀, the scale of v is given by σ_NAU₀, 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, δ₁ = 2Δx, δ₂ = 15Δx, and T_s = 10³Δ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 [13]. 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, x₀. 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 x₀ = 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 [29]. 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 [28]. 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).

Conclusion

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 [30]. 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.

Competing interests

JS is Chief Science Officer of Mediwave Star Technology, Inc.

Authors' contributions

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.

Acknowledgements

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.

References

1. Murray J: Mathematical Biology. 2nd edition. New York: Springer; 1993.

2. Mikhailov A: Foundations of Synergetics I. Distributed Active Systems. 2 revised edition. New York: Springer; 1994.

3. 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. PubMed Abstract

4. Noble D: A modification of the Hodgkin-Huxley equations applicable to Purkinje fibre action and pacemaker potentials.

   Journal of Physiology 1962 , 160:317-352. PubMed Abstract | Publisher Full Text | PubMed Central Full Text

5. 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. PubMed Abstract | Publisher Full Text

6. 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. PubMed Abstract | Publisher Full Text

7. 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. PubMed Abstract | Publisher Full Text

8. Karma A: Electrical alternans and spiral wave breakup in cardiac tissue.

   Chaos 1994 , 4:461-472. PubMed Abstract | Publisher Full Text

9. Fenton F, Karma A: Vortex dynamics in three-dimensional continuous myocardium with fiber rotation: filament instability and fibrillation.

   Chaos 1998 , 8:20-47. PubMed Abstract | Publisher Full Text

10. 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. Publisher Full Text

11. 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. PubMed Abstract | Publisher Full Text

12. 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. PubMed Abstract | Publisher Full Text

13. 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. PubMed Abstract

14. 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. PubMed Abstract | Publisher Full Text | PubMed Central Full Text

15. 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. PubMed Abstract | Publisher Full Text

16. 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. PubMed Abstract | Publisher Full Text | PubMed Central Full Text

17. 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. PubMed Abstract | Publisher Full Text

18. 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. Publisher Full Text

19. 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. PubMed Abstract | Publisher Full Text

20. 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. PubMed Abstract | Publisher Full Text

21. 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. PubMed Abstract | Publisher Full Text

22. 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. Publisher Full Text

23. Mitchell C, Schaeffer D: A two-current model for the dynamics of cardiac membrane.

    Bulletin of Mathematical Biology 2003 , 65:767-793. PubMed Abstract | Publisher Full Text

24. 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. PubMed Abstract | Publisher Full Text | PubMed Central Full Text

25. 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. PubMed Abstract | Publisher Full Text

26. Otani N, Gilmour R: A memory model for the electrical properties of local cardiac systems.

    Journal of Theoretical Biology 1997 , 187:409-436. PubMed Abstract | Publisher Full Text

27. 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. PubMed Abstract | Publisher Full Text

28. 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. PubMed Abstract | Publisher Full Text

29. Chernyak Y, Starobin J, Cohen R: Class of exactly solvable models of excitable media.

    Physical Review Letters 1998 , 80:5675-5678. Publisher Full Text

30. 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. Publisher Full Text

31. 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. PubMed Abstract | Publisher Full Text | PubMed Central Full Text

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

33. Ravens U, Wettwer E: Electrophysiological aspects of changes in heart rate.

    Basic Res Cardiol. 1998 , 93 Suppl 1:60-65. PubMed Abstract | Publisher Full Text

34. Rushton WAH: Initiation of the Propagated Disturbance.

    Proceedings of the Royal Society of London 1937 , 124:210-243. Publisher Full Text

35. 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. PubMed Abstract | Publisher Full Text | PubMed Central Full Text

Have something to say? Post a comment on this article!