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

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, e = τ f /τ r << 1, reflects the two order-of-magnitude difference in time constants between the fast excitation, τ f , and slow recovery, τ r , proc-esses, 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][8][9][10][11][12][13][14][15][16][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][19][20][21][22][23][24][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][21][22][23][24][25][26][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][32][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: 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: 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 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 0 , the scale of v is given by σ NA U 0 , 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 σ 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 , defined above, where A = 1.2, δ 1 = 2Δx, δ 2 = 15Δx, and T s = 10 3 Δ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 0 . 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 0 = 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 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. 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 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]. 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 steadystate 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.  Figure 3 the value of v min was 10% greater than the value of .

Restitution curves
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.