Hysteretic Dynamics of Multi-Stable Early Afterdepolarisations with Repolarisation Reserve Attenuation: A Potential Dynamical Mechanism for Cardiac Arrhythmias

Some cardiovascular and non-cardiovascular drugs frequently cause excessive prolongation of the cardiac action potential (AP) and lead to the development of early afterdepolarisations (EADs), which trigger lethal ventricular arrhythmias. Combining computer simulations in APs with numerical calculations based on dynamical system theory, we investigated stability changes of APs observed in a paced human ventricular myocyte model by decreasing and/or increasing the rapid (I Kr) and slow (I Ks) components of delayed rectifying K+ current. Upon reducing I Kr, the APs without EADs (no-EAD response) showed gradual prolongation of AP duration (APD), and were annihilated without AP configuration changes due to the occurrence of saddle-node bifurcations. This annihilation caused a transition to an AP with EADs as a new stable steady state. Furthermore, reducing repolarisation currents (repolarisation reserve attenuation) evoked multi-stable states consisting of APs with different APDs, and caused multiple hysteretic dynamics. Depending on initial ion circumstances within ventricular myocytes, these multi-stable AP states might increase the local/global heterogeneity of AP repolarisations in the ventricle. Thus, the EAD-induced arrhythmias with repolarisation reserve attenuation might be attributed to the APD variability caused by multi-stability in cardiac AP dynamics.

where T is the pacing cycle length (i.e. period of the current stimuli) and Td is the duration for which the stimulus current is sustained at the maximum value (Istim,max). The first and second equations in Equation S1 corresponds to the current injection phase and the no-injection phase of the external current stimulus, respectively. During the current injection phase and the noinjection phase, the ventricular myocyte (VM) model becomes an autonomous system with constant parameters Istim,max and zero. Therefore, the paced VM model can be formalised as a composite dynamical system such that the two corresponding autonomous systems are successively switched over time.

Poincaré map of the composite dynamical system
To investigate stability changes in periodic oscillations when a system parameter is changed, we directly assessed the dynamical stability of AP responses using a method involving the Poincaré map. Even though the VM model driven by a discontinuous periodic force (Eq. S1) had a discontinuous nature, the Poincare map could be constructed numerically as successive Consider the following non-autonomous differential equations consisting of Equation 1 and a periodic parameter variation of Equation S1 during the time satisfying t−t0 (mod T) [0,T): where t R denotes time, x R n is the state vector, R m-1 denotes common parameters for f and a R is a parameter specifying f1. The parameters a correspond to Istim,max. We also assumed that f is periodic in T so that f(t + T, x, ) = f(t, x, ), for all t. If we assumed that the solution to Equation S2 is described as a mixed solution of the first and second equations of , then the solution with the initial condition x = x0 at t = t0 is represented by: Let φ1 and φ2 correspond to solutions to the first and second equations of Equation S2 ), ( ), , ; Then, the Poincaré map can be defined as a composite map:

Numerical calculation of the fixed point and detection of bifurcation points
We defined a fixed-point equation based on the fixed point on the Poincaré map (Eq. S4): where x0 R n denotes the initial value at t = t0. If an initial value (x0) satisfies Equation S5, then this point is denoted a fixed point. We can also define the characteristic equation of the fixed point as follows: where In is the n×n identity matrix, and ∂M(x0)/∂x0 denotes the derivative of M with respect to The fixed point equation of Equation S5 cannot be analytically solved; therefore, we used Newton's method, which is a numerical approach to computation. In the following, we assume that the initial value in Equation S2 is given by x0, and x0 (k) is a first-guess of the fixed point.
The recurrent formula for Newton's method is given by: where δ is the correction term, and DH(x0 (k) ) is the Jacobian matrix with respect to the initial value x0, denoted by: Then, the derivative of M with regard to x0 is expressed by: The second equation of Equation S7 must be solved for δ, using a suitable method such as Gauss elimination. H(x0 (k) ) can be obtained from the original Equation S5. Consequently, to obtain each element of the Jacobian matrix, we need to differentiate M with respect to x0.
We defined the solution starting from x0 at t = t0 as: ). , Substituting the solution of this equation into Equation S2, we get: Differentiating this equation by x0 yields: The order of differentiation on the left-hand side is commutative, and the following equation is obtained from the right-hand side: (S13) This equation is of the following form: where Y  ∂/∂x0 is the matrix solution of a variable coefficient linear differential equation, 13 referred to as a variational equation for Equation S2. Then, from Equation S10: (S14) Therefore, we can use Equation S14 as the initial value to calculate the derivative of M with regard to x0, by numerically integrating Equation S13 from t = t0 to t0 + T. The Runge-Kutta method was used, after which Newton's method was performed. The fixed point x0 is accurately located by iteration. However, the derivatives ∂/∂x0 cannot be directly defined due to discontinuities in the derivative of the solution. To avoid the impossibility of the derivative at t = t0 + Td and t0 + T, the first derivative of M with respect to x0 is given by obtaining the derivatives of the sub-maps, successively: In Equation S6, the derivatives of φi, for i = 1, 2 relate to the initial value xk for k = 0, 1, which corresponds to fundamental matrix solutions (i.e. ∂φi /∂xk), and can be obtained by numerically integrating each of the first-order variational equations: and and substituting t = t0 + Td and t0 + T in the respective solutions to Equation S7 and S8.