{'fraction_non_alphanumeric': 0.03419263646765834, 'fraction_numerical': 0.017163626573332435, 'mean_word_length': 4.694902261402836, 'pattern_counts': {'":': 0, '<': 1, '<?xml version=': 0, '>': 1, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Action potential restitution and hysteresis in a reaction-diffusion system with pacing rate dependent excitation threshold

We have demonstrated that rate dependent restitution and action potential duration-refractory period hysteresis can be reproduced in a one-dimensional two-variable Chernyak-Starobin-Cohen reaction-diffusion medium with variable excitation threshold. We show that restitution and hysteresis depend on the relationship between pacing period and steady state excitation threshold and also on the rate of excitation threshold adaptation after an abrupt change in pacing period. It was also observed that the onset of action potential duration alternans is determined by the minimal stable wavefront speed, which could be approximated by the analytical critical speed of a stable solitary pulse. This approximation was suitably accurate regardless of the adaptation constant of excitation threshold, its dependence on pacing interval, or magnitude of the slopes of restitution curves.


I. INTRODUCTION
Repetitive pacing of biological reaction diffusion media by over-threshold stimuli gives rise to periodically propagating excitation waves. In cardiac tissue the duration of excitation referred to as action potential duration, T AP , as well as the speeds of excitation wavefront and waveback, depend on previous stimulation periods and refractory (diastolic) intervals, T DI [1,2,3]. The analysis of such dependences known as restitution curves has been established as an effective method for evaluation of normal functioning of the heart [4,5,6,7,8,9].
Pioneering experimental studies [10,11,12] established two major pacing sequences such as steady state (dynamic) and S1-S2 pacing protocols, which led to two standard T AP restitution dependences. It has been shown that the protocol dependent rates of T AP adaptation were different for stepwise S1-S2 perturbations of cycle length acceleration and deceleration [12,13]. This phenomenon was later experimentally introduced as action potential duration cycle length hysteresis [14,15] and has been recently associated with cardiac ischemia, coronary flow reduction [16,17], and cardiac memory [18,19].
A stepwise change in pacing rate following a long series of conditioning S1 stimuli may result in prolonged adaptation of action potential duration to its new steady-state value.
The process of such adaptation can extend well beyond the T AP response to the first S2 test stimulus [11,12] and give rise to the additional constant BCL (basic cycle length) transient restitution component attributed to cardiac memory [2,20,21,22]. Relationship between this phenomenon and stability of pulse propagation as well as with rate dependent T AP -T DI hysteresis has been recently investigated in experimental [17,20,22,23] and theoretical [24,25,26,27,28] studies.
It has been found that in the presence of restitution transients cardiac dynamics is more complex than predicted by Nolasco and Dahlen restitution criterion [29]. Specifically, it was demonstrated that dynamic restitution curve slopes > 1 may not automatically indicate loss of excitation wave stability and subsequent appearance of alternans [20,24,25,26,27,30].
Computational experiments with Fenton-Karma and Mitchell-Schaeffer ionic models readily identified that such effects can be quantified by fitting T AP and T DI for specific restitution curves using time dependent gate variables [24,25,31,32,33].
In this paper, we analyze excitation wave propagation in a one-dimensional cable based on the approach, which follows from direct experimental observations of the dependence of cardiac muscle resting potential on frequency of external pacing [34,35,36]. We implement a two-variable exactly solvable Chernyak-Starobin-Cohen (CSC) reaction-diffusion model [37,38] and modify it accordingly to incorporate pacing rate driven adjustments of resting potential as a rate dependent excitation threshold. We use an exponential-like evolution of excitation threshold, V r , that takes place over the course of multiple heart beats following stepwise changes in pacing rate [35]. We demonstrate that adaptation of V r significantly affects the stability of pulse propagation and gives rise to rate dependent T AP -T DI hysteresis.
In particular, we find that under given medium parameters, regardless of the slopes of restitution curves, the appearance of alternans is determined by proximity of the wavefront speed to the minimal speed of a stationary solitary pulse determined analytically in [37].

II. METHODS
Basic equations that describe a class of exactly solvable models for excitable media have been defined in [37,38]. 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 in the Appendix.
The pacing function, P (x, t), is defined as a product of two functions, X(x) and T (t).
Each is composed of the Heaviside step function, Θ(x), as follows: The right-hand term B in Eq. 3 responds to the stepwise evolution of pacing period, T m .
Stepwise changes in B result in smooth exponential transition of the excitation threshold from one steady-state plateau to another. A steady-state value of excitation threshold at each m th plateau, V m r ≡ B ± m , was chosen to be linearly dependent on the corresponding pacing interval, T m [34,35].
Here α ± and β ± are positive parameters that determine the amplitude of change of B ± m between two consecutive pacing plateaus for the increasing and decreasing rate, respectively.

III. NUMERICAL SIMULATIONS
The system of Eqs. 1-3 was solved numerically on a short cable of 250 grid points with spatial and temporal grid intervals of ∆x = 0.13 and ∆t = 7.2 × 10 −4 , respectively. The length of the cable was approximately equal to the width of the pulse to reflect the relative dimensions of the heart and a propagating cardiac pulse at moderate heart rates. Periodic wavetrains were produced by stimulating the cable with a square wave at the left end using the function P (x, t) = X(x)T (t), defined above, where A = 10, δ 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 [9]. A typical solution is depicted in Fig. 1, showing the propagation of a single pulse at three instants of time.
In order to quantify the dynamics of the system (Eqs x 0 . In order to analyze a developed pulse, we measured intervals and speeds at x 0 = 20 (except Sec. III D) where the speed of the wavefront had reached a constant steady-state value.

A. Restitution for constant excitation threshold
The system of equations 1-3 was initially studied with Eq. 3 replaced by its asymptotic form ∂Vr ∂t = 0, which is equivalent to τ = ∞. The cable was stimulated periodically for forty consecutive pacing periods over the range T m = 70 to 25 with decrements of 1.5. At the end of each plateau, T AP and T DI had reached steady-state values that were used to compose the steady-state restitution curve. Two values of T m , at 30 and 27, were used to obtain conventional S1-S2 restitution curves [10,11]. The pulse, T n+1 AP , resulting from the test stimulus, S2, following the conditioning sequence and the last diastolic interval, T n DI , from the conditioning S1 plateau composed the S1-S2 restitution curve.
Steady state and S1-S2 restitution curves computed for a pair of constant excitation thresholds and a pair of S1 basic cycle lengths are shown in Fig. 2. We observed that the maximal difference between steady state and S1-S2 restitution curves comprised less than 15% of the corresponding steady state value of T AP , which implied that transient responses to any premature stimulus were on average limited to just a single non-stationary pulse.
The ends of both steady state restitution curves shown in Fig. 2  The upper level insert shows the dispersion curve computed analytically for a steady state solitary pulse [37]. The critical speed and duration of a stable solitary pulse for the model parameters, λ, ζ, and ǫ, described above is c crit = 0.48 and T crit AP = 5.4. As shown in the insert, the difference between the numerical and analytical critical speeds is relatively small and amounts only to 10% of the wavefront speed determined from our numerical model at V r = 0.215. This suggests a criterion for determining the critical speed below which pulse durations start to oscillate. On the contrary, the conventional Nolasco-Dahlen critical slope stability criterion is unsuitable as the maximum slope of the dynamic restitution curve at this V r is 30% greater than one [29].
B. Restitution and hysteresis for rate-dependent V r Unlike the previous section, the evolution of V r according to Eq. 3 resulted in a set of prolonged transients initiated by abrupt changes in stimulation rate. These T AP transients, which constituted the constant BCL restitution (negative slope, Fig. 3A), had a duration of 5-50 stimulation periods depending on the adaptation constant, τ (Fig. 3D). The phase of constant BCL adaptation followed the immediate S1-S2 responses, which were on the contrary positive and aligned with the steady state restitution curves for each steady state excitation threshold (gray lines, Fig. 3A).
Different steady state values of excitation threshold during progressively increasing and decreasing stimulation rates gave rise to T AP -T DI hysteresis. The cable was stimulated periodically for fifty consecutive periods at a series of S1 conditioning plateaus. The series consisted of five plateaus with decreasing T m followed by the same number of plateaus with increasing T m . The excitation threshold evolved according to Eq. 3, and its steady state value at each plateau was related to the pacing period by Eq. 4.
Using this protocol we demonstrated that higher values of T AP -T DI hysteresis corre- sponded to greater differences between B + m and B − m during the decremental and incremental stages of the pacing protocol (Fig. 3B). We also found that T AP -T DI hysteresis was dependent on the adaptation constant, τ , if T AP and T DI were measured before their steady-state values were reached. For instance, for a time lag of seven stimulation intervals (Fig. 3D), the magnitude of T AP -T DI hysteresis increased as the adaptation constant increased from τ = 32 to τ = 216 (Fig. 3C). as computed analytically in [37]. Similar critical values for our short cable were described earlier in Sec. III A for a medium with constant excitation threshold. In this section, we analyze perturbations of T AP and wavefront speed when steady state pulse speeds are below a certain value close to the critical speed of a solitary pulse.
We perturbed both pulse duration and wavefront speed near their analytical critical values, T crit AP and c crit , as indicated by arrows labeled "A,B" and "C,D" (Fig. 4). If the speeds of the wavefronts elicited after perturbation exceeded c crit less than 22%, we observed T AP

alternans (inserts A and B) that developed between two values indicated by open markers at
the end of the arrow "A,B". Alternans occurred sooner when the adaptation rate was higher due to a smaller adaptation constant, τ (insert A). The minimal non-alternating wavefront speeds computed at x 0 = 20 for different B + m were the same regardless of changes of B + m or the four fold increase of the steady-state restitution slope (Fig. 5). On the contrary, when the perturbation resulted in a wavefront speed that exceeded c crit by more than 22% (the end of the arrow "C,D"), similar bifurcations did not occur. Instead, we observed typical exponential adaptation from one steady-state to another (inserts C and D).

D. Influence of propagation on the spatial distribution of hysteresis and alternans
We demonstrate that T AP adaptation to a stepwise change in pacing interval is different for different points of observation along the cable. When the distance between the point alternans. The closest point at which alternans can be observed is located at the midpoint of the cable (Fig. 6C).
We also observed that the magnitude of T AP -T DI hysteresis is larger when measured further from the stimulation site. Figure 7 shows the twofold increase of the magnitude of hysteresis and T AP between two observation points at x 0 = 4 and x 0 = 20.

IV. CONCLUSIONS
We have demonstrated that rate dependent restitution and T AP -T DI interval hysteresis can be reproduced in a one-dimensional two-variable CSC reaction-diffusion medium where excitation threshold adjusts to changes in pacing rate. We show that the rate dependence of restitution and hysteresis are influenced by two major factors. The first one is the adaptation constant, τ , of excitation threshold evolution after an abrupt change in pacing interval. The second is the dependence, B ± m , of the steady-state excitation threshold on the pacing period. We show that steady-state and S1-S2 restitution curves coincide if the excitation threshold is constant, which corroborates with earlier findings for other reaction diffusion models [3,26,28]. On the contrary, the steady-state and S1-S2 restitution curves diverge if the magnitude of excitation threshold varies with changes in pacing rate, leading to prolonged T AP transients following an abrupt change in pacing interval.
We found that larger values of T AP -T DI interval hysteresis were associated with greater hysteresis resulted from increasing adaptation constants, τ , if T AP and T DI were measured before their steady-state values were reached.
Our numerical simulations showed that stimulating the cable with short pacing intervals and high excitation thresholds elicited slower pulses that led to T AP alternans. It was observed that the minimal stable wavefront speed could be approximated by the analytical critical speed of a stable solitary pulse [37,38]. This approximation was suitably accurate regardless of values of B ± m and magnitudes of the slopes of restitution curves. We also found that the onset of alternans occurring after an abrupt change in pacing rate was more delayed for larger values of τ .

APPENDIX:
The scale of u is the maximum steady-state action potential amplitude U 0 , the scale of v is given by σ f U 0 , and the time scale is C m /σ f , where σ f corresponds to the maximum sodium conductance and C m is the membrane capacitance. The characteristic length scale is given by D/σ f , where D is the diffusion coefficient. The small parameter ǫ ≪ 1, is equal