Short-term action potential memory and electrical restitution: A cellular computational study on the stability of cardiac repolarization under dynamic pacing

Electrical restitution (ER) is a major determinant of repolarization stability and, under fast pacing rate, it reveals memory properties of the cardiac action potential (AP), whose dynamics have never been fully elucidated, nor their ionic mechanisms. Previous studies have looked at ER mainly in terms of changes in AP duration (APD) when the preceding diastolic interval (DI) changes and described dynamic conditions where this relationship shows hysteresis which, in turn, has been proposed as a marker of short-term AP memory and repolarization stability. By means of numerical simulations of a non-propagated human ventricular AP, we show here that measuring ER as APD versus the preceding cycle length (CL) provides additional information on repolarization dynamics which is not contained in the companion formulation. We focus particularly on fast pacing rate conditions with a beat-to-beat variable CL, where memory properties emerge from APD vs CL and not from APD vs DI and should thus be stored in APD and not in DI. We provide an ion-currents characterization of such conditions under periodic and random CL variability, and show that the memory stored in APD plays a stabilizing role on AP repolarization under pacing rate perturbations. The gating kinetics of L-type calcium current seems to be the main determinant of this safety mechanism. We also show that, at fast pacing rate and under otherwise identical pacing conditions, a periodically beat-to-beat changing CL is more effective than a random one in stabilizing repolarization. In summary, we propose a novel view of short-term AP memory, differentially stored between systole and diastole, which opens a number of methodological and theoretical implications for the understanding of arrhythmia development.


Introduction
The longer duration of cardiac action potential (AP), as compared for example with that of skeletal muscle and neuronal cells, has the functional significance of controlling and a1111111111 a1111111111 a1111111111 a1111111111 a1111111111 response for sudden changes in constant CL values [2,5]. However, it has never been explored in dynamic conditions and at high pacing rate.
By systematically studying the behavior of a well-established human cardiac ventricular AP model [34] under dynamic pacing conditions, we found that hysteresis in dER CL also develops and carries relevant information on AP dynamics that is not contained in dER DI data, and ought to be taken into account to understand repolarization stability. In doing so, we also found that a given extent of beat-to-beat CL variability, at a high pacing rate, makes ventricular repolarization more stable. In addition to this previously made observation [35], we found that periodic, rather than stochastic, beat-to-beat CL changes are more effective in stabilizing repolarization. What emerges from our study is a novel picture of cardiac AP dynamics under variable pacing, where short-term memory is a rate-dependent property, differentially shared by systole (APD) and diastole (DI) depending on the dynamics of pacing variability.

Methods
All simulations reported in this study were performed by means of the ten Tusscher et al. 2006 [34] human ventricular AP model, which will be referred to in this manuscript as TP06. The CellML format of the model [36] was recompiled in its Matlab version by means of COR facility at http://www.cor.physiol.ox.ac.uk. The 'ode15s' solver built into the R2016a version of Matlab (The Math-Works, Inc., USA) was used to integrate the model equations. All simulations were run on a PC with Intel Core i5, 2.5 GHz CPU. APs were elicited by simulating 3 ms-long current injections with an amplitude 1.5 times the current threshold. AP duration was measured as APD -60mV , i.e. the time between the maximum first derivative of membrane potential (V m ) during the initial fast depolarization phase and the time when V m had fallen to a value of -60 mV. The rate dependence (RD) of APD was measured as follows: at different CLs, from 300 ms up to 1400 ms, step 20 ms, a 1000 beats pacing train was simulated to allow the AP waveform to reach a steady-state configuration. The average values of the last 20 beats of each sequence was taken for each CL and used to draw the RD curve. Classic electrical restitution (ER) was measured by conditioning training the membrane at a given cycle length (CL Ã ), delivering an extra-stimulus delayed to within CL Ã ± 20 ms, and reporting the resulting APD versus the CL or DI value of the cycle containing the last conditioning beat. ER was also measured during dynamic, random or periodic pacing by saving the 19×1 vector of TP06 variables at the end of the (n-1) th cycle preceding the last conditioning beat (n th beat), and applying the classic ER protocol described above while making the n th CL varying within a given range; thus ER CL is the function f described by APD n+1 = f(CL n ) and ER DI by APD n+1 = f(DI n ). In some simulations the time-dependence of the gating variables for given ion currents (I CaL , I Kr , and I Ks ) was removed using the following procedure. A simulation was run at a given CL (after a conditioning train of 1000 beats) and the current-voltage relationship of the ion current of interest flowing during the AP was plotted and fitted with variable order polynomials. The purely voltage-dependent function was then used to update ion current values during AP simulation, instead of solving the corresponding gating equations. The substitution with the voltage-dependent fitted values was not applied within the first 9 ms of the AP potential due to the fast kinetics of currents in that phase and to the fact that our main interest was in late repolarization dynamics (see also the S3 Fig). Fig 1A and 1B shows the dependence of the TP06 AP model from pacing rate by reporting RD (red) and ER (blue) curves, the latter derived for 3 values of conditioning CL (320, 420, and 520 ms), measured as ER CL (panel A) and ER DI (panel B). At a given steady-state, if CL suddenly changes to another constant value, APD will instantly change according to ER CL and, after a given number of beats (N b ), reach a new steady-state given by RD CL ; in other words, it takes a different N b for each ER curve to rotate (clockwise for ER CL and anti-clockwise for ER DI ) into the RD curve (panels C and D). N b (and the corresponding time, not shown) decreases mono-exponentially with the conditioning CL (panel E).

Steady and dynamic AP states
N b measures the evolution from the sudden (ER) to the stationary (RD) rate dependence of APD for a membrane in its steady-state. Indeed, the membrane characterized in both the RD The inset in panel A shows AP waveforms measured at 320 and 520 ms. Panels C and D show schematically the state of the membrane at the n th beat after a constant pacing train (blue dot), with the angle between ER and RD. The number of beats (N b ) taken for ER to rotate clockwise (panel C) and anti-clockwise (panel D) into RD is shown in panel E as a function of conditioning CL. and the ER protocols described above is in steady-state, i.e. there are no time-dependent electrogenic processes occurring at the time of the measurement. From now on we will call the average value of conditioning CL, "CL Ã " (consequently, this definition can hold either for constant or variable pacing), and define the state of the membrane for a given n th beat with one point in the (CL, APD)-or (DI, APD)-space with coordinates (CL n-1 , APD n ) or (DI n-1 , APD n ); all the information needed to assign the state of the n+1 th beat is in the ER CL (or ER DI ) obtained after the n th beat (which we call its ER curve). In steady-state conditions (constant conditioning CL) each state belongs to its ER curve (Fig 1C and 1D); it generally does not do so in dynamic conditions (beat-to-beat variable conditioning CL), regardless of whether CL n is equal or different from CL n-1 (Fig 2), whether we are considering ER CL or ER DI .
Random changes in high pacing rate. An example of dynamic state is that in which the pacing CL changes randomly within a given percentage of a central value. In this case, the time law for CL is: where rand(-1,1) is a random number between -1 and 1, generated for each beat, and clv the half range of CL variability. The time at which beats are elicited is: Twenty consecutive superimposed APs (from a sequence of 1000) simulated under random pacing (CL Ã = 350 ms, clv = 35 ms) are shown in the left-hand panel of Fig 3A. There is no fixed relationship between consecutive CLs, since they can assume any value within CL Ã ± 35 ms. Given that the states do not belong to their ER curves, the system moves beat-to-beat on a family of ER curves (Fig 3A, second column) along a random path. We will call the collection of the states occupied by the system the "space of states". The space of states for random pacing takes a trapezoidal shape (third column), limited superiorly by the highest ER curve of the family and inferiorly by the lowest one. We note that, for the same conditions, the space of states in the APD vs DI representation is a single ER DI function (fourth column). Moreover, the vertical width of the family of ER CL curves decreases in parallel with the pacing rate and, for CL > 500 ms, the space of states essentially coincides with a single ER CL curve.
Periodic changes in high pacing rate. Another example of a dynamic state is that in which the pacing CL oscillates within a given range around the central CL Ã value. This can be simulated by making CL change like: where ɷ is the frequency of oscillation. Thus, the time at which beats are elicited is: and 2π/ɷ is the number of beats required for a complete oscillation. A sequence of consecutive APs simulated according to Eq 3 (CL Ã = 350 ms, clv = 35 ms, ɷ = 2.4 4) is shown in Fig 3B. Consecutive beats belong to a family of ER CL curves, although in this case their path is forced by the regular oscillations of CL to move clockwise along a closed hysteresis loop, which represents the space of states for this type of pacing variability. Again, for the same conditions the space of states in the APD vs DI representation is a single ER DI function (fourth column). In addition, as pacing frequency decreases, for the same pacing parameters, the hysteresis loop collapses into a single ER CL curve. Spaces of states under beat-to-beat CL-varying pacing describe what we have called dynamic ER (see Methods section), which we will refer to as dER. A third type of dynamic state is that in which the pacing CL alternates around the central value CL Ã , such as CL(N) = CL Ã + (-1) N clv [37]. In this case, consecutive beats belong to only two ER CL curves, which are the upper and the lower ones of the family of ER CL curves obtained by the pacing protocol described by Eqs (1) and (3). Results from this alternating protocol were reported at the 42 nd EWGCCE Meeting of the European Society of Cardiology [38], but have not been included in the present study.

Short-term memory
The three pacing protocols described above reveal an intrinsic property of AP repolarization that is associated with short-term AP memory via the transition from a single to a family of ER curves, and only emerges in certain conditions. Thus, we quantify AP memory as the vertical width of the space of states, i.e. the distance between the lower and upper ER curves, which we will call, depending on the type of representation adopted, Shift CL and Shift DI , since they measure the vertical shift of the family of ER curves introduced during dynamic pacing (see Fig  3B). A necessary condition for such a shift to develop is a high pacing frequency, and a key factor is the size of the beat-to-beat CL changes which, under random pacing is uniquely determined by clv, whereas under periodic pacing it is determined by clv and ω.  Cardiac action potential memory is shared between systole and diastole together with the RD curve (red). For small ɷ values, the states, as in the case of steady-state, belong to their ER curves, at the intersection with the RD curve (panel B, top). As ɷ increases, the average (through consecutive beats) distance between states and ER curves increases. Shift CL is negligible when ɷ is low (< 0.2) but becomes significant for higher ɷ values, leading to hysteresis (top panels in Fig 4C). The opposite is true for ER DI (bottom panels in Fig 4C), where hysteresis is absent for high ɷ values and start to appear for low ɷ values.
The results of a closer inspection for hysteresis in dER CL and dER DI under periodically changing pacing and for different values of CL Ã (300 to 400 ms, step 10 ms, clv = 70 ms) and ɷ (0.0 to 2.5 step 0.1) are shown in a color code (note the different scale) in Fig 5. The dependence of both Shift CL and Shift DI from ɷ for a CL Ã = 350 ms (see also vertical broken line in panel A and B) is shown in panel C. When pacing CL varies periodically within ± 70 ms around a value of 350 ms, hysteresis in dER DI appears only within a narrow range of small ɷ values, whereas dER CL hysteresis is present and much larger for most of the ɷ values and increases with ɷ. Thus, hysteresis in dER develops at high pacing rate under periodically varying CL, where this memory effect is shared between APDvsCL and APDvsDI representations, depending on the rate of beat-to-beat CL changes.
A summary of the dER configurations associated with constant and dynamic pacing is schematically shown in Fig 6 in the case of dER CL representations, and can be extended to dER DI . Three additional AP ventricular models, the human Priebe and Beuckelmann [39], the O'Hara and Rudy [40] and the rabbit Mahajan et al [41], were tested for the same properties and showed qualitatively similar results, the main difference being related to the CL range (at high pacing rate) where the hysteresis effect appears. This range was always chosen to avoid the CL range for period-doubling bifurcation where APD alternans occurs, when constant pacing CL is lowered below 210 ms (O'Hara and Rudy model), 345 ms (Priebe and Beuckelmann), 239 ms (Mahajan model), and 250 ms (TP06 model). The same range was recently explored in a canine model of ventricular AP by McIntyre and colleagues [42] who found that beat-to-beat CL variability (0-6%, rather than the 10-20% considered here), promotes alternans by shifting the upper limit of its occurrence towards higher values.

Ion currents and dER during periodic pacing
So far, we have shown that the space of AP states is one point at constant-, a closed curve in periodic-, and an entire surface in random pacing-conditions (Fig 6). Here we analyze the three ion currents mainly responsible for AP repolarization in the human ventricle because of their role in determining and modulating the size and shape of the space of states under periodic pacing (Fig 7). Role of G max . The steady-state APD of the TP06 model was prolonged by the same amount (12.5%) by decreasing, in turn, G max of I Kr or I Ks , and increasing G max of I CaL (-90%, -60%, and +90% respectively). An opposite maneuver (+100%, +98%, and -40%) led to a 12.5% decrease in APD in the three cases ( Fig 7A). The dER CL was then measured by simulating a pacing train with a CL varying periodically according to Eq 3 (ɷ = 2.4, CL Ã = 350 ms, clv = 35 ms), first in control conditions (black), and then with modified G max values. The greatest increase in the hysteresis loop was obtained under a decrease in G Ks (blue). Results obtained with randomly varying CLs led to comparable results where spaces of states were trapezoidal surfaces instead of a closed loop (not shown). When, from the same data, dER DI curves were derived for the same pacing conditions, they resulted as single unimodal curves in all instances; the decrease in G Kr led to an increase in the dER DI slope from 0.49 (control) to 0.71, a decrease in G Ks to 0.83, and an increase in G CaL to 0.62. Role of time-dependency. A time-independent version of each of the three ion currents was obtained by recording the current during a simulated AP in steady-state control conditions (CL = 350 ms), fitting its current-voltage relationship with a polynomial, and replacing the corresponding TP06 gating equations with the purely voltage-dependent fit. The TP06 model was then paced as in the previous section and dER CL measured in control (Fig 7C,  black), and in conditions where, in turn, I Kr , I Ks , and I CaL were purely voltage-dependent during AP repolarization phase (green, blue, and red respectively). Whereas the time-dependence of I Kr does not seem to play a role in dER CL hysteresis, the removal of the time-dependence of I Ks increased Shift CL while that of I CaL almost abolished it under periodic pacing. The timecourse of the three ion currents after the fitting procedure can be seen in the S3 Fig.

Perturbing AP states
Monitoring dER under different pacing conditions helps to study AP repolarization dynamics when a given relatively regular pacing is suddenly perturbed and immediately restored. Perturbation may be in the form of a pre-or post-mature stimulus or, as in our case, of a single missing beat, which, when seen from the ventricle, can be thought of as a pause in the SA nodal firing or as a single failure in AV nodal conduction. Fig 8 shows results obtained by simulating pacing trains where CL was, in turn, kept constant, made to change periodically (ɷ = 2.4, clv = 35 ms), or randomly (clv = 35 ms). At the 10 th beat of each sequence the current stimulus was set to zero (star), so that no AP was elicited, and immediately restored in the following cycle (panel A). The number of beats (N b ) required for each AP train to resume the un-perturbed APD sequence (with a ΔAPD < 2 ms) after the pause was measured for the 3 protocols.
For the sake of the internal clarity of this discussion, and without claiming generality on this delicate theoretical issue, we use the term 'stability' here as the ability of repolarization to resume its unperturbed state after a pacing perturbation, which, in the present setting, takes the form of a missing beat. Beat-to-beat variable pacing is always associated with greater Nb (smaller stability) than constant pacing (Fig 8B), even if Nb value is progressively decreasing (and stability increasing) with the rate and amplitude (ɷ and clv) of changes and is always smaller under periodic than random ones (see S1 and S2 Figs).
Since N b can change slightly in dynamic control conditions, depending on the position of the missing beat within the sequence, we will consider its average over several beats here; it increases from 2.6 during constant pacing, to about 6.0 during periodic pacing, and 9.2 during random pacing (Fig 8B). When APD was prolonged by increasing G CaL or decreasing G Ks , a stable 1:1 alternant APD oscillation (within the simulated 500 beats) arose each time a beat was missing during constant pacing at CL Ã = 350 ms, whereas APD prolongation obtained by decreasing G Kr only led to an increased N b (with respect to the control). Consecutive AP waveforms and the time course of corresponding APDs are shown in Fig 8 for the three protocols under control conditions (C), and in the case when APD was prolonged by decreasing G Ks (D). Classic ER DI and ER CL curves measured under constant pacing as described in Methods, are shown in Fig 9A and 9C for the two conditions, control (black), and diminished G Ks (red), adopted in Fig 8C and 8D. The numbers in square brackets show the range of the ER DI and ER CL slopes when, as in the present case, the CL of the restitution protocol was made to change within ±35 ms around CL Ã = 350 ms. When measured at the constant conditioning value of CL (350 ms) or DI (119.1 ms in control, and 85.9 ms under treatment), the 60% decrease in G Ks led to an increase of the ER DI and ER CL slopes from 0.5 to 0.8. In the same figure the dER DI and dER CL representations of periodic pacing simulated in Fig 8C and 8D are also shown in panels B and D respectively. The dER DI slope reaches values up to 1.6 under G Ks decrease, whereas the slope of the major axis of the hysteresis loop of dER CL goes from 0.7 to 1.1 in the same conditions. For a given clv value, repolarization stability depends on the frequency of CL changes (ω). This is shown in Fig 10, where the TP06 model, APD-prolonged by a 60%-decrease of G Ks and paced with a constant CL Ã = 350 ms, always developed APD alternans after a missing beat (star). When paced with a CL periodically varying around 350 ms (clv = 35 ms) and with a high enough angular frequency (ω = 2.4), dER CL developed hysteresis, and alternans was stopped after only a few beats (top panels A and B). As ω decreases, the vertical width of the hysteresis loop (Shift CL ) decreases as well (Fig 10C, green), while the number Nb of post- . The difference between APDs following the missing beat and those of the unperturbed sequence were measured and the number of beats (Nb) required for such a difference to reach a value < 2 ms measured. (B) The histogram shows Nb measured during constant (black), periodic (orange), and random (grey) pacing (CL Ã = 350 ms, clv = 35 ms, ω = 2.4) in control conditions and for each AP prolonging treatment (see Fig 7). In the case of an increase in G CaL and a decrease in G Ks the missing beat triggered stable APD alternans. A transient alternans between short APD values was also present under a decrease in G Ks , but progressively disappeared within 30 beats. (C) Example of data generating histogram in B. Consecutive AP waveforms (left-hand column) and APD time course (on the right) are shown for 10 s simulations in control conditions under constant, periodic, and random pacing and with a missing beat (star) in each sequence. (D) Same for AP prolonged by decrease in G Ks .
https://doi.org/10.1371/journal.pone.0193416.g008 perturbation beats falling out of the hysteresis loop increases (Fig 10C, blue). It should be noted that, the anti-arrhythmic effect (ability of preventing alternans) obtained for ω > 2.0 is achieved despite the expected pro-arrhythmic increase in the dER DI slope (single unimodal curve at this ω values) up to values greater than 1 (Fig 9B). Thus, if on the one hand pacing variability within a given range (clv) tends, per se, to make repolarization less stable (Fig 8B), on the other the rate of beat-to-beat CL changes (ω and, with it, Shift CL ) counteracts this effect. This is particularly evident under G Ks downregulation, when a single missing beat during fast pacing leads to stable APD alternans, which is prevented by large beat-to-beat oscillations of CL (Fig 8).
Not only can a randomly or periodically changing pacing frequency prevent persistent APD alternans, which develops at a constant high pacing rate; it can also stop it when already established, bringing the system back to its unperturbed space of states. This is shown in the simulations in Fig 11: the AP model was first paced for about 17 s at a constant CL Ã = 350 ms,  Fig 8C and 8D, in control conditions (black) and in diminished G Ks conditions. In brackets, the range of slopes when CL was made to change ±35 ms around 350 ms. In the case of dER CL (panel D) the slope of hysteresis was approximated by measuring the slope of the major axis (dotted lines). https://doi.org/10.1371/journal.pone.0193416.g009 Cardiac action potential memory is shared between systole and diastole which was suddenly switched to a periodically (ω = 2.4, clv = 35 ms) or randomly (clv = 35 ms) varying CL sequence. A single missing beat (star) during constant pacing made APD alternate. Alternans was stopped within a few beats (double arrow) after the switching. If beat-to-beat periodic CL changes were slowed (ω 2.0), the switch to variable pacing did not stop the alternans. A comparison of the space of dER CL and dER DI states for periodic pacing at ω = 2.4 and ω = 2.0 is reported in Fig 11C; whereas Shift CL decreases from 73 ms to 44 ms (double arrows in figure), the shape and the slope of dER DI remain unmodified.
When, as shown in Fig 12,  Cardiac action potential memory is shared between systole and diastole space of states (see Fig 7C). Such ability was preserved when the same intervention was performed on I Kr and I Ks . Thus, given a certain range of beat-to-beat CL variability, a mechanism emerges that counteracts other instability sources. This can be unmasked by dER CL but not by dER DI representations, appears to be linked to the time-dependence of I CaL , and can make repolarization more stable.

Discussion
Measurement of the electrical restitution of cardiac ventricular AP has been used previously to study memory and stability of repolarization, which have also been related to hysteresis in APD vs DI that develop in dynamic conditions [25,29,43]. By means of simulations carried out on a human ventricular AP model [34], we show that hysteresis affects both APD vs DI and APD vs CL representations of dER, and that the latter is a better marker of memory and stability. We present a generalization of the mechanism underlying hysteresis in different pacing conditions, investigate the role of 3 plateau ion currents in modulating this property, and analyze the link between hysteresis and AP stability after transient perturbations of beating rate. We have found that, at high pacing rate, by increasing the frequency of periodic beat-tobeat CL changes, the space of dER CL states changes from the typical unimodal curve to an hysteresis loop; dER DI does the same, though for much slower beat-to-beat CL changes and to a smaller extent. A randomly varying pacing frequency leads to a scattered dER CL space of states covering an area whose shape and size depend on the pacing parameters. The geometry of the space of states for a given pacing protocol is mainly determined, in the AP model under study, by the interplay of voltage and time-dependency of I CaL and I Ks . At high pacing frequency, it is the size of the dER CL space (and not that of the dER DI ) that correlates with repolarization stability after a single missing beat, while periodic more than random pacing is effective in damping sudden APD oscillations.

ER CL vs ER DI
APD vs DI representations (dER DI ) of sequences of cardiac APs have long been adopted for measuring dynamic electrical restitution [44], as opposed to APD vs CL (dER CL ). The main reason for this choice is the meaning of the ER DI slope for AP dynamics under constant pacing CL, with a slope > 1 amplifying, and < 1 suppressing APD oscillations after CL disturbances [45], even though the general applicability of such a concept remains controversial [3,14,34,46]. Moreover, the concept holds only as long as CL is constant, whereas most of the time we are interested in dynamic conditions where CL is the independent variable driving both APD and DI. Wu and Patwardhan developed an experimental protocol to make DI the independent variable [25] by timing each consecutive current stimulus with a programmed DI, which can follow any time-law, and is therefore independent from the preceding APD. More recently, Elizabeth Cherry and coworkers from one side and Sharon Zlochiver and coworkers from another applied a constant-DI pacing protocol, which appears to be effective in differentiating and controlling alternans mechanisms [47,48]. However, time-clamping DI can hardly be associated with any in vivo condition and, if this does not make it less appropriate for investigating membrane properties, it does make this approach less clear from a practical point of view. In addition, the hysteresis in dER DI representations is less prominent in the TP06 model when compared to that in dER CL . Finally, the conditions that bring about hysteresis in dER DI , which have been described earlier [25], are different from those in dER CL (Fig 5), and, since in our simulations mainly the latter seem to be critical for AP stability, we are devoting most of our work to exploring AP sequences as dER CL .

Hysteresis and memory
The use of restitution hysteresis as a marker and measure of short-term AP memory has been studied elsewhere for dER DI [15,25,26,29,49] and several models have been proposed to investigate it at different CL values [21,43,32]. Short-term AP memory has been seen as a function that accumulates information during systole and dissipates it during diastole [14]. Our own simulations confirm that the information concerning the recent past and future of an AP can be measured by the hysteresis of dER, but they also suggest that such information is shared in a diverse way among systole (APD) and diastole (DI) in different dynamic pacing conditions, with a prevalence of the former in the model under study. Others have found, by looking at dER DI hysteresis, much larger Shift DI values at comparable pacing rates in porcine ventricular tissue [26] and it is likely that a different balance between Shift DI and Shift CL would be found in different AP models and in different species. In fact, heterogeneity of restitution properties has been documented even in different cell types from the same heart [49], and the two components of the delayed rectifier potassium current, differently contributing to dER, have also been shown to be heterogeneously distributed along the base-apex direction [50].
If we pace the TP06 AP model with a periodically changing CL Ã centered at 350 ms within a ± 20% variability range, and regardless of whether we drive CL or DI to achieve such changes (not shown), we observe measurable hysteresis in dER CL for ω values greater than 0.2, whereas in dER DI hysteresis is only measurable for ω values lower than 1.0, with a small ω interval where both are present (Fig 5C). Hysteresis in dER CL and not in dER DI means that information is stored in APD, since at any CL (n th beat), the APD choice for the "future" (n+1 th beat) does not depend on the previous DI (no choice, monotonic dER curve) but on the previous CL (double choice depending on previous history, i.e. hysteresis), and therefore on APD (= CL-DI). The information concerning the recent pacing history is stored in CL and not in DI, thus we can conceive it as accumulated in APD (see for example Fig 4C left, hysteresis in dER CL and not in dER DI ). Previous DI does not distinguish whether CL is increasing or decreasing whereas previous APD (and therefore CL) does. In other words, at high pacing frequency (CL Ã < 400 ms), if CL changes periodically from beat-to-beat at a fast rate (ω > 1, see Fig 5C), AP memory accumulates mainly on APD. For slower CL changes, the opposite happens, since what it is that discriminates between the two possible future states of the system now depends on DI (whether it increases or decreases) and, in fact, can be seen in dER DI (Fig 4C right). Thus, to characterize short-term AP memory, both representations should be explored, since each contains information that cannot be derived from the other, and this suggests that short-term AP memory is shared between systole (APD) and diastole (DI) depending on pacing conditions.

Pacing perturbations
Although the link between hysteresis and repolarization stability has been previously described in the case of dER DI [25], our results show that dER CL can be used to understand and predict, for instance, conditions when otherwise silent pathologies of ion currents might become proarrhythmic, and also how pacing rate variability can control the transition to arrhythmias. An example of this can be found, for instance, in the defective I Ks due to KCNQ1 mutation of Type 1 long QT syndrome [51], which we simulated in the experiment in Figs 11 and 12. The inherited APD prolongation observed in this pathology becomes malignant under adrenergic stimulation when, it has been proposed [52], part of the repolarization reserve (mainly I Kr ) is rate-dependently removed, and a pacing perturbation can initiate alternans [53] and, in turn, trigger arrhythmias. In these conditions, however, pacing variability introduces a protecting mechanism from arrhythmias, which can be monitored in the shape of the corresponding space of dER CL states (Fig 7B). In fact, although under these circumstances (periodic pacing, CL Ã = 350 ms, ɷ = 2.4, clv = 35 ms, see Fig 11A) the increase in dER DI slope (Fig 9B) would predict a decreased repolarization stability, the large beat-to-beat oscillations of CL bring about an increase in Shift CL and, with this, in APD stability too (Figs 10C and 11A). Accordingly, for values of ω = 2 (Shift CL diminished by 40%, Shift DI unmodified) and lower (see Fig  11C), the switch from constant to periodic pacing did not stop alternans. The 68% increase in the slope of classical ER DI (Fig 9A) does predict an increased instability under G Ks decrease, but is measured under constant pacing.
The ability of periodic pacing to prevent alternans is better described by dER, which captures electrical restitution properties under dynamic pacing, like that shown in Fig 8C and 8D (middle panels), and dER CL rather than dER DI , with its hysteresis behavior, provides a link between this ability and short term AP memory.

Role of active membrane properties
Ion currents dynamics is known to determine and affect cardiac electrical restitution [26,54], hence our interest in studying its role in shaping the space of states in general, and particularly that of dER CL . The difficulty with the role of repolarization ion currents and ER is the fact that their up/down-regulation brings about APD changes which, per se, can affect restitution properties [55]. We therefore consider changes in the maximum conductance of different ion currents that lead to the same shortening/prolonging effect on APD (Fig 7A and 7B). Taken together, our data suggest a major involvement of both the amplitude and time-dependency of I CaL in determining size, thickness and slope of the hysteresis loop in dER CL . Accordingly, when G CaL is increased/decreased, Shift CL also increases/decreases, and when the time-dependence of I CaL is removed, Shift CL is virtually abolished. This latter finding is particularly interesting when compared with the results on rabbit ventricular myocytes of Mahajan and coauthors [56], who proposed that flattening ER DI by inhibition of I CaL inactivation rather than by current block could be viewed as a way to improve repolarization stability without depressing contractility. Thus, our own results also confirm that targeting kinetics of I CaL recovery rather than its amplitude may be a promising anti-arrhythmic strategy.
The repolarizing I Ks partially counterbalance I CaL all over the AP and particularly in the late repolarization phase. Thus, when its G max is reduced or the time-dependence of I Ks is removed, the positive effect of I CaL on Shift CL , partially freed from this restraint, can develop even further (Fig 7C, central column). The fast kinetics of I Kr (compared to that of I Ks ) perhaps explain why removal of time-dependence from this current does not modify Shift CL .
By demonstrating the same effect with a verapamil-induced reduction of I CaL in isolated ventricular tissue from pigs, Guzman [26] noted the contrasting effect of slope and memory on stability, suggesting a careful consideration of both effects when testing antiarrhythmic properties of I CaL blockers. We have extended this notion to our data on dER CL and note the fact that in the simulations in Fig 7B, the corresponding dER DI curves are all unimodal (not shown), and only much smaller ω (~0.2) values made them assume hysteresis behavior. Thus, both phenomena pertain to AP memory, though only the one at higher ω values (dER CL ), as pointed out above, can be extrapolated to comparable random pacing conditions, and explain the protective effect on repolarization after pacing perturbations.

Limitations
The present work has been performed on a specific model of human ventricular AP. Although preliminary results obtained with additional models (2 humans and 1 rabbit), which we briefly mention above, seem to confirm our conclusions, only further work could completely rule out a possible model-dependency in our findings. This is particularly true for our discussion about the role of inactivation of I CaL in modulating memory and stability, where models with a more detailed description of calcium current kinetics and intracellular calcium handling could better prove our conclusions.

Conclusions
Our data show that, at high pacing frequency, information concerning the recent pacing history is stored differentially in systole (APD) and diastole (DI). This kind of AP memory, particularly that measured in dER CL , turns on in extreme pacing variability conditions, where it counterbalance the pro-arrhythmic effect of abrupt changes in CL. Any partial or complete loss of this safety mechanism, as in the case of changes in amplitude or kinetics of repolarization ion currents, exposes ventricular repolarization to an increased risk of arrhythmias, which can be monitored in dER CL and not in dER DI . Screening ventricular AP repolarization for APD and DI memory, both in vivo and in the growing family of cardiac human ventricular AP models, may lead to a refinement of clinical criteria for electrical instability and may benefit many aspects of cardiac physio-pathology, such as an understanding of the latent instability of certain inherited repolarization pathologies, the design of new antiarrhythmic drugs, or the challenge of intelligent sensing in artificial pacemakers. Of further relevance in this regard is our finding of the greater effectiveness of periodic rather than random pacing in dampening sudden APD oscillations, for which we have provided a statistical explanation (see S1 and S2 Figs), and which again is more effectively revealed in the dER CL , rather than dER DI representations. One possible application of this finding might be in designing pacing protocols for artificial pacemakers, where a predefined periodic variability can be added to the pacing rate to reduce susceptibility to arrhythmias.
Supporting information S1 Fig. APD displacement after a missing beat. This figure shows how the data in Fig 8B were obtained. AP sequences were simulated under constant, periodic, and random pacing conditions (CL Ã = 350 ms, ω = 2.4, clv = 35 ms), where a given n th beat was missing and the difference between the 50 APDs following the missing beat and the those of the unperturbed sequence was measured. The same was done for sequences where the missing beat was the (n+1) th , and then (n+2) th , until (n+19) th . Panel a reports the averages, for each pacing conditions, of the 20 traces described above. The same was done for AP sequences where APD was 12.5% prolonged by modifying maximum conductance of ion channels (see Results).