Comparisons of wave dynamics in Hodgkin-Huxley and Markov-state formalisms for the Sodium (Na) channel in some mathematical models for human cardiac tissue

We compare and contrast spiral- and scroll-wave dynamics in five different mathematical models for cardiac tissue. The first is the TP06 model, due to ten Tusscher and Panfilov, which is based on the Hodgkin-Huxley formalism; the remaining four are Markov-state models, MM1 WT and MM2 WT, for the wild-type (WT) Na channel, and MM1 MUT and MM2 MUT, for the mutant Na channel. Our results are based on extensive direct numerical simulations of waves of electrical activation in these models, in two- and three-dimensional (2D and 3D) homogeneous simulation domains and also in domains with localised heterogeneities, either obstacles with randomly distributed inexcitable regions or mutant cells in a wild-type background. Our study brings out the sensitive dependence of spiral- and scroll-wave dynamics on these five models and the parameters that define them. We also explore the control of spiral-wave turbulence in these models.


Introduction
The development of an understanding of the dynamics of waves of electrical activation in cardiac tissue is a problem of central importance in research on life-threatening cardiac arrhythmias, because sudden cardiac death is responsible for roughly half of the deaths because of cardiovascular disease, i.e., 15% of all deaths globally [4]. Approximately 80% of sudden cardiac deaths arise from ventricular arrhythmias [4]. Such arrhythmias are often associated with the formation of spiral or scroll waves of electrical activation; unbroken spirals or scrolls lead to ventricular tachycardia (VT), whereas broken waves, with spiral-or scroll-wave turbulence [5,6,7,8,9,10], are responsible for ventricular fibrillation (VF); VT and VF lead to the malfunctioning of the pumping mechanism of the heart, so, in the absence of medical intervention, VF leads to sudden cardiac death. It is very important, therefore, to study VT and VF by using all means possible, namely, in vivo, in vitro, and in silico investigations, the probability of opening of this channel is different for the TP06, MM1 WT, and MM2 WT models. The peak value of this probability and the time duration of this opening are also dissimilar in MM1 WT and MM2 WT models. These differences alter the action potential (AP) and its morphology. We show that, for the mutant (MUT) Na channels, the failure of inactivation leads to early afterdepolarizations (EADs), in the APs in MM1 MUT and MM2 MUT models [20].
The differences in the WT Na peak amplitude lead to disparate upstroke velocities in these models, which manifest themselves in dissimilar CVRs. Furthermore, the conduction velocities (CVs) in MM1 WT and MM2 WT models turn out to be outside (lower than) the accepted range for CV in the human myocardium; we show that we can obtain CVs in this range if we increase the diffusion constant D in both MM1 WT and MM2 WT models. The differences in our single-cell and cable-level results motivate our study of wave dynamics in mathematical models for cardiac tissue, which use these different models.
We carry out a variety of simulations in 2D homogeneous domains to show that spiral-wave dynamics, in TP06, MM1 WT, MM2 WT, MM1 MUT, and MM2 MUT models, depends sensitively on these models. For example, we demonstrate that, in the MM1 WT (MM2 WT) model, the spiral wave is stable (unstable, meandering spiral). The formation of certain EADs can lead to backward propagation of the wave, and rapid spiral breakup, in the MM2 MUT model; by contrast, in the MM1 MUT model, EADs are somewhat different, so we do not find such backward propagation, the mother rotor is unaffected, and there is only far-field breakup of the spiral. Furthermore, the spatiotemporal evolution of a spiral wave in the MM2 WT model depends sensitively on the time τ S2 between the application of the S1 and S2 impulses, which we use to initiate spiral waves.
In the case of mutant models, because of the different kinds of EADs that we find in MM1 MUT and MM2 MUT, these models display qualitatively different electrical-wave dynamics. Furthermore, in the spirit of the studies of Refs. [6,7,21,22], we investigate the effects of two types of inhomogenieties on spiral-wave dynamics in these models: (a) Two-dimensional (2D), circular or three-dimensional (3D), cylindrical obstacles, with a random distribution of inexcitable regions, to model fibrotic patches in Markov-state WT models; P f , the percentage of inexcitable obstacles, and the radius of the obstacle are important control parameters. (b) A circular patch of mutant cells in an otherwise homogeneous, 2D WT domain; we find that a spiral wave is formed in the MM2 MUT model, but not in the MM1 MUT model, if we pace the tissue at a high frequency.
The elimination of spiral-and scroll-wave turbulence is of central importance in developing lowamplitude defibrillation schemes for the elimination of VT and VF. In the Supplementary Material, we describe one such defibrillation scheme (control of spiral waves) for the models of we study.
We carry out a few illustrative studies of scroll waves in 3D TP06, MM1 WT, and MM2 WT models. We show, in a homogeneous domain, that scroll waves are stable in TP06 and MM1 WT models, but not in the MM2 model. We also investigate when scroll waves are anchored or broken up by cylindrical obstacles, of the type described above.
Finally, we perform a parameter-sensitivity analysis for TP06 and MM1 WT and MM2 WT models, in which we consider three important dependent variables, namely, the APD, V max , and V rest , at the cellular level and CV and APD at the cable level (see Supplementary Material).
The remaining part of this paper is organized as follows. Section 2 is devoted to Methods and Simulations. In Section 3 we report our Results for single-cell studies and tissue-level simulations in 2D square and 3D slab domains for WT models; we also present, for MUT models, single-cell and 2D-simulation results. In Section 4 , Discussion and Conclusions, we end with concluding remarks. 3

Model
The electrical behavior of a single cardiac myocyte is governed by the following ordinary differential Equation (ODE) for the transmembrane potential V m : (1) here, I ion is the sum of all the ionic currents, I i is the current because of the i th ion-channel, and C m is the normalized, transmembrane capacitance. In the parent TP06 model, I ion is the sum of the following 12 ionic currents (Table 1): The spatiotemporal evolution of V m , at the tissue level, is governed by the following reaction-diffusion partial differential Equation (PDE): where D is the diffusion constant; for simplicity we consider the case in which D is a scalar.

Hodgkin-Huxley Model
The TP06 model uses the Hodgkin-Huxley formalism for the WT Na channel. The macroscopic current through this channel is governed by the three gating variables m, h, and j [1]; the first of these is an activation gate and the latter two are inactivation gates; the gating dynamics and the Na current are given by da n dt = a ∞ − a n τ n ; here a n can be m, h, or j; G N a is the maximal sodium-channel conductance, a ∞ is the steady-state value of a n , τ n the time constant of this gating variable, and E N a is the sodium-channel Nernst potential.

Markov-state models
We consider four Markov-state models (MMs): two of these are for the wild-type (WT) and the other two for the mutant (MUT) Na channels. We use the Markov-state formalisms of Ref. [2] for the first WT and MUT Na channels; we refer to these as MM1 WT and MM1 MUT, respectively. We use the Markov models of Ref. [3] for the second WT and MUT Na channels, which we label MM2 WT and MM2 MUT, respectively. We then replace the Na current in the TP06 model by these two different WT and two different MUT models. Finally, we have three different WT models, i.e., the original TP06, MM1 WT, and MM2 WT; and we have two different MUT models, namely, MM1 MUT and MM2 MUT. All the other currents in the original TP06 model are unaltered in our studies below. Schematic diagrams of MM1 WT and MM2 WT models are shown in the top panel of Figure 1. The orange, double-headed arrows indicate transitions between such Markov states; transition rates for the rightward (leftward) transition are given above (below) these arrows, e.g., a 111 (b 111 ) for the IC3 → IC2 (IC2 → IC3) transition in MM1 WT. Similar schematic diagrams for the MM1 MUT and MM2 MUT models are given in the bottom panel of Figure 1. The MM1 MUT model has the same number of states as the MM1 WT model, but the transition rates between the Markov states are different. The MM2 MUT model has 12 states: 8 of these are as in the MM2 WT model; in addition, there are 4 bursting states, namely, BO, BC1, BC2, BC3.
The dynamics of transitions between the states of these Markov models and the Na-channel current I N a are given, respectively, by Equations 7 and 8 below: where l, k label the states O, C1, C2, C3, IC2, IC3, IM1, IM2, IF, IS, BO, BC1, BC2, BC3 (Figure 1); α l and β l are generic labels for forward and backward transition rates, respectively; where P O is the probability of the opening of the Na channel.
We use the values of G N a that are employed either in the original TP06 model [1] or in Ref. [2]; specifically, we use The major differences between the two MM (WT and MUT) models are as follows; • The number of states and the connections between them (see Figure 1 • The MM2 WT model has a Na current with a late component; this is absent in the MM1 WT model.
We show below that these differences can have significant effects on single-cell and spiral-wave properties in these models.

Numerical Simulations
For our single-cell simulations, we use the Rush-Larsen method to solve Equation 5 for the HHM; for the MM Equation 7 we use the implicit trapezoidal method of Ref. [3] and the forward-Euler method for Equation 1. In our 2D tissue simulations for Equation 4, we use a square domain with N × N grid points, a fixed space step of size ∆x = 0.025 cm, and a time step ∆t = 0.02 ms; in 3D we use a slab domain (see below). The accuracy of the numerical scheme is tested and is reported in Supplementary Material.For the homogeneous tissue we consider, D is a scalar; we use D = 0.00154 cm 2 /ms, as in Ref. [1]; this yields a maximum plane-wave conduction velocity 70cm/s, which is in the biophysically reasonable range for human ventricular tissue [23]. For the Laplacian in Equation 4, we use five-point and seven-point stencils, respectively, in 2D and 3D. We impose no-flux boundary conditions. When we study the effects of heterogeneities, we introduce, in a localized region of our simulation domain, a patch of inexcitable cells, which are decoupled from adjoining cells (effectively, D = 0 in this patch).
To obtain the inactivation and activation properties of MM1 WT, MM2 WT, MM1 MUT, and MM2 MUT Na channels, we use the voltage-clamp-simulation protocol of Ref. [2]. In the activation protocol, we clamp the cell with a voltage V c , which ranges from the hyper-polarized regime (below the resting membrane potential of an action potential) to the depolarized regime ( 50 mV); we do this in steps of 5 mV; the clamping is maintained for a clamping time t c = 1s. We then record the peak current I peakN a (V c ) and divide it by the driving force (V − E N a ) to obtain the conductance G(V c ), which we normalize to obtain the activation variable A as follows: V c = −100 mV to 50 mV , t = 1s; Similarly, for the inactivation protocol, we use a holding potential (V h ), ranging between hyperpolarized and the depolarized values, and apply it for 250 ms; we then apply a test potential V t = 0 mV , record the peak Na current, and then define the inactivation variable I as follows: , t > 250 ms; We record the single-cell AP and its morphology, after we have paced the cell with n pulses, each with a constant pacing cycle length (PCL); we use n = 500 pulses. We obtain the static action potential duration restitution (APDR) (s1-s2) as follows: (a) we apply several pulses ( 17) with a fixed PCL (s1); then, once the system reaches the steady state, we change the diastolic interval (DI) (s2) by recording the time at which the cell is 90% repolarized; this is the action potential duration (APD 90 ), or simply the APD. We next obtain the dynamic conduction velocity restitution (CVR) (s1-s1) as follows: We consider a cable of cells (of dimension 832 × 10) and pace it by applying a current stimulus (s1) at one of its ends; we obtain CV from the time that an iso-potential line takes to move between two cells, which are separated by a fixed distance. The CVR is the plot of CV versus DI.
In our 2D studies we use two representative square domains, namely, one with 1024 × 1024 grid points, for our spiral-wave studies, and another with 512 × 512 grid points, when we pace of the simulation domain along an edge. In our studies of scroll-wave dynamics we use a 3D slab domain with 1024 × 1024 × 40 grid points ( 25.6cm × 25.6cm × 1cm). We initiate spiral and scroll waves in such domains by using the following S1-S2 cross-field protocol, with stimuli amplitudes of 150 pA/pF and durations of 3 ms: We allow a plane wave (S1) to propagate in the domain along a particular direction; as it propagates, we start another plane wave (S2), in a direction perpendicular to the S1 wave; this results in a conduction block and, eventually, the formation of a spiral wave (2D) or scroll wave (3D). We vary the time interval τ S2 between the S1 and S2 impulses to study the sensitive dependence of the spiral-wave dynamics on τ S2 .
We have carried out the following two sets of simulations of electrical-wave dynamics with certain localised inhomogeneities in an otherwise homogeneous simulation domain. (a) In the first set of simulations we introduce circular (2D) and cylindrical (3D) regions with a random distribution of inexcitable obstacles to mimic localised fibrotic patches in Markov-state WT models; an important control parameter here is P f , the percentage of inexcitable obstacles in the circular or cylindrical regions. (b) In the second set of simulations, we examine electrical-wave dynamics in the presence of a circular patch of mutant cells in an otherwise homogeneous, 2D WT domain.

Activation and Inactivation
We begin with a comparison of the activation A and the inactivation I in the models MM1 (WT and MUT) and MM2 (WT and MUT) with their TP06 counterparts. Figures 2A and 2B show, for different values of V c , plots of I N a versus time t. These plots illustrate, respectively, the activation and inactivation protocols (see 2.1), whence we obtain Figures 2C and 2D, which depict, respectively, the dependences of A and I on V m for the MM1 WT, MM2 WT, MM1 MUT, MM2 MUT, and TP06 models; for the TP06 model A = m 3 ∞ and I = j ∞ × h ∞ . If we contrast the plots of A in Figure 2C, we see that the curves for both MM1 WT and MM1 MUT models lie to the right of, and are less steep than, their TP06, MM2 WT, and MM2 MUT counterparts; hence, activation occurs most slowly (with respect to V m ) in MM1 WT and MM1 MUT models. In particular, the Na channels in the TP06, MM2 WT, and MM2 MUT models reach near-complete activation at V 2mV , V −14mV , V −8mV , respectively, whereas the MM1 WT and MM1 MUT models do so only for V 60mV and V 60mV , respectively. The potential at which Figure 2D, we find that the plots of I, for the MM1 WT, MM1 MUT, MM2 WT, and MM2 MUT models, are shifted to the right (the depolarized-potential side) compared to I for the TP06 model. Note also that the MM1 MUT model shows inactivation earlier than the MM1 WT model. In the MM2 MUT case, the inactivation occurs earlier than in the MM2 WT model in the Our results for the two Markov-state models are in agreement with those shown in Figure 2 of Refs. [2,19].

Probabilities of the Markov states
Let us examine now the temporal evolution of the probabilities of different Markov states during the course of an action potential. As we have mentioned in Section 2.1, there are three main classes of Markov states for the Na channel, namely, the open states (O, BO), the inactivation states (IF, IS, IM1, IM2, IC2, IC3), and the closed states (C1, C2, C3, BC1, BC2, BC3). The probabilities of these three classes of states are as follows: P O is the open state probability; P I is the sum of probabilities of all the inactivation states; and P C is the sum of probabilities of all the closed states. In the case of MM1 WT and MM2 WT models, they are as follows: The probabilities of these three main classes of states, in MM1 MUT and MM2 MUT models, are as follows: Figure 3 shows plots of P O , P I , and P C versus time t for the Na channel in the course of an action potential for MM1 (MM2) models in the top (bottom) panel; the blue and red curves are for WT and MUT models, respectively. We obtain these plots by pacing a single cell with P CL = 3000 ms. By comparing the blue curves in Figures 3 (A), (B) , (D) and (E)) we find that the duration for which the Na channel is in the inactivation or closed states, i.e., the time interval during which P I = 1 and P C = 0 (inactivation state) or P I = 0 and P C = 1 (closed state), is approximately the same in MM1 WT and MM2 WT models. In contrast, the duration for which P O is significantly greater than 0 differs in MM1 WT and MM2 WT models (compare the blue curves in

Action Potential, APDR, and CVR
We pace a single cell with the following three different values of PCL: high frequency (PCL=300 ms), intermediate frequency (PCL=650 ms), and low frequency (PCL=1000 ms). We present the steady-state AP and the Na current I N a at the top panel of Figure 4(A),(B), and (C) (for PCL = 1000 ms), for the three WT models; and we compare the morphological properties of the APs of these models in Table 2.
Given the differences in the activation profiles in Figures 2(C),(D) and the plots of P O for the MM models in Figures 3(C).(F) , we observe that (a) the times at which the Na channels open are different in all the three WT models; and (b) the amplitude of I N a is comparable in MM2 WT (−312.74 pA/pF ) and TP06 (−300.43 pA/pF ) models, but it is significantly lower in the MM1 WT model (−144 pA/pF ) as we show in Figure 4(B) . These differences in the amplitude of I N a affect the maximum voltage and the upstroke-velocity of the AP ( Table 2). The upstroke velocities for TP06, MM1 WT, and MM2 WT models are markedly different (Table 2). Also, there is the late component of the Na current I N a,L (Figure 4(C)) in the case of MM2 WT; this component is clearly absent in TP06 and MM1 WT models. P O becomes significant at 350ms in the MM2 WT model (Figure 3(F)), so the APD for this model is larger than its counterparts in the TP06 and MM1 WT models (see Figure 4(C) and Table 2).

PCL = 300ms
As we decrease PCL, say to 300ms, we find that both I N amax (the maximal value of −I N a ) and the upstroke velocity in the MM2 WT model increase relative to their counterparts in the TP06 model as we show in Table 2 (contrast this with our results for PCL = 1000ms). These increases occur principally because P O,max is higher in the MM2 WT model than in the TP06 model.

PCL = 3000ms, MUT Na channel
In Figures 4(D),(E) and (F), we show, respectively, plots of V m , I N a f , and I N a L versus time t; we use dashed curves for the MM1 MUT (blue) and MM2 MUT (red) models and the illustrative value PCL = 3000ms. Clearly, the MM1 MUT and MM2 MUT APs in Figure 4(D) show early afterdepolarizations (EADs) [20,24], insofar as their APs are prolonged considerably relative to the the APs for MM1 WT and MM2 WT models, because of the failure of inactivation near the repolarisation region (insets in Figure 3).

APDR and CVR (WT)
For the TP06, MM1 WT, and MM2 WT models, we present plots of the single-cell static APDR ( Figure 5(A)) and the dynamic CVR ( Figure 5(B)), for a one-dimensional cable of cells. The APDR profiles for TP06 and MM1 WT lie close to each other, but the MM2 WT curve lies above these, because of the late current component I N a,L (see above). The slopes of the APDR and CVR profiles are given, respectively, in Figures 5(C) and (D). Note that, in all these three models, the maximal slope of the APDR profile > 1 (it is highest in the MM1 WT model). The Na channel determines the upstroke velocity at the cellular level; therefore, this channel plays an important role in determining CV, in cardiac tissue, and also CVR plots ( Figure 5(B)). From these plots we find that, for TP06, MM1 WT, and MM2 WT models, CV is nearly independent of DI, for large DI; the ranges spanned by

2D results
We have explored differences between the TP06, MM1, and MM2 models at the single-cell and the cable levels. We now compare spiral-and scroll-wave dynamics in these models by carrying out detailed numerical simulations in 2D (Section 3.2) and 3D (Section 3.3) domains.

Wild-type Na channel
We contrast, in the top panel of Figure 6, spiral waves in these three models, with D = 0.00154 cm 2 /ms. We find that spiral waves in TP06 and MM1 WT are stable and they rotate with frequencies ω 4.75 Hz and ω 4.25 Hz, respectively; in particular, the low value of CV (40.41 cm/s), in the MM1 WT model with D = 0.00154 cm 2 /ms, does not alter the spiral-wave dynamics qualitatively. By contrast, in the MM2 WT model, the spiral wave is unstable and exhibits transient breakup; it is not possible to isolate a single cause for this break up, but the late Na current I N a,L (Figure 4(C)) plays an important role in this instability; we have checked that, by increasing β 12 , we can reduce the magnitude of this late current and thus suppress spiral-wave turbulence (the spiral meanders but does not break up into multiple spirals as we show in the Movie (M0) in the Supplementary Material 4).
We We employ the S1-S2 protocol to initiate spiral waves in all these models (Section 2). The pseudocolor plots of V m in Figure 7 show that the spiral-wave activity in the TP06 and MM1 WT models is independent of the time τ S2 , at which the S2 pulse is applied after the S1 pulse (we use 560ms ≤ τ S2 ≤ 620ms). By contrast, in the MM2 WT model, we observe spiral-wave breakup for τ S2 = 560ms and 580ms until the end of our simulation, i.e., 10s; but spiral-wave activity vanishes for τ S2 = 600ms at 5.5s and for τ S2 = 620ms at 6.8s (Figure 7 and Movie (M2) in the Supplementary Material 4).
Spiral-wave dynamics in the MM2 WT model depends on the time τ S2 at which we initiate the S2 pulse. It behooves us, therefore, to examine whether obstacles (or conduction inhomogeneities) affect spiral-wave activity in the MM1 WT and MM2 WT models, for it has been shown, for HH-type models for cardiac tissue, that spiral-wave dynamics depends sensitively on the position, size, and shape of such obstacles [6,7,21,22]. Our obstacles consist of inexcitable points that are distributed randomly within a circular region of radius R; P f is the percentage of the area of the circle that has inexcitable obstacles. Given our experience with studies of spiral-wave dynamics with such obstacles in HH-type models, we expect that, as P f increases, such an obstacle should anchor a spiral wave [25,26]. Therefore, we investigate the dependence of spiral-wave dynamics on P f and R in the MM1 WT and MM2 WT models and compare this with its counterpart in the TP06 model, for different values of τ S2 . Illustrative plots from our simulations are shown in Figure 8.
We find that, for the TP06 and MM1 WT models, the anchoring of the spiral wave depends on R and on P f , but not on τ S2 . The time period T of the anchored spiral increases with R and P f as we show in Figure 9(A); but T decreases for lower percentages (e.g., P f = 30%) in TP06 and MM1 WT models at large values of R (Figure 9(A)). The interaction of the tip of the spiral wave with the obstacle is complicated. In particular, this depends on how much of the region, inside the circular patch, is excitable. For low values of P f , this excitable region forms a tortuous but spanning cluster (in the sense of percolation theory [27]), so the tip of the spiral propagates inside the obstacle, the wave of activation is slightly deformed there, but then it re-emerges into the homogeneous part of the simulation domain. If P f is large, the excitable region can still be tortuous, but it does not form a spanning cluster, so the tip of the spiral rotates around the obstacle, and is anchored to it, but does not propagates inside it. To quantify the effect of our obstacle on the spiral wave we calculate δT ≡ (T − T 0 ), where T is the time period (or inverse of the rotation frequency ω), at a given set of values of P f and R, T 0 is the time period (or inverse of the corresponding frequency ω 0 ) with P f = 100% for the same value of R. Clearly T must depend on P f and R. The plots in Figures 9(B)and (D) show, for TP06 and MM1 WT models, the dependence of δT on R for different values of P f . Given these plots, we identify three regions, namely, (i) δT < 0, i.e., ω > ω 0 , (ii) δT > 0, i.e., ω < ω 0 and (iii) δT = 0 i.e., ω = ω 0 . If δT > 0, then the frequency ω ∼ T −1 , for a given pair (R, P f ), is less than ω 0 ∼ T −1 0 (for R, P f = 100%); this may occur because the spiral core penetrates the obstacle because of a spanning cluster of excitable regions inside the obstacle. In Figures 9(C) and (E) we show different colored regions in the (R, P f ) plane for TP06 and MM1 WT models, respectively: light blue indicates an increase in ω relative to ω 0 (caused by penetration of the spiral core); light green is for a decrease in ω relative to ω 0 (accompanied by penetration of the spiral core); yellow indicates no penetration of the spiral core into the obstacle; dark blue depicts regions in which there is no change in ω relative to ω 0 even though the spiral core penetrates into the obstacle.
For the MM2 WT model, the minimum size R min for spiral anchoring is large, compared to that in TP06 and MM1 WT model, and is R min = 1.875 cm (Movie (M5) in the Supplementary Material 4). The threshold percentage in the MM2 WT case is P f,min 50%. Once we reach the values R min and P f,min required for anchoring, the spiral activity is independent of τ S2 , as we show in the fourth row of Figure 8. The dependence of the spiral rotation time period T on R, for different values of P f , is shown in Figure 9(A). Stability diagrams for the spiral-wave activity, in the presence of localized, inexcitable obstacles distributed within a circular region of radius R, are shown in the (R, τ S2 ) plane, for different values of P f in the MM2 WT model, in Figure 10; brown, green, and blue denote regions with an anchored spiral, spiral breakup, and no activity, respectively.

Mutant Na channel
The mutant Na channel fails to inactivate completely in the MM1 MUT and MM2 MUT models; this leads to prolonged EADs, as we have shown in Sec 3.1 and Figure 4. We find that two of the types of EADs that have been discussed in Ref. [20] occur in both these MUT models: there is a single EAD (of type 2 in the nomenclature of Ref. [20]), in the MM1 MUT model, and an oscillatory EAD (roughly of type 3 in the nomenclature of Ref. [20]), in the MM2 MUT model. These two types of EADs affect the wave dynamics differently, as we demonstrate explicitly by simulating plane-wave propagation in our 2D domain, but with all mutant myocytes. We observe backward propagation of the plane wave in the MM2 MUT model because of the oscillatory EADs; by contrast, there is no such backward propagation in the MM1 MUT model. If we initiate a spiral wave in both these models, then, (a) in the MM1 MUT model, we get almost-instantaneous far-field breakup away from the core, but the mother rotor is unaffected, and (b) in the MM2 MUT model, we obtain almost-instantaneous spiral break-up (bottom panels of Figures 6(D) and (E)). The complete spatiotemporal evolution of such spiral-wave dynamics, for both these cases, is shown in the Movie (M3) in the Supplementary Material 4.
Although there are several studies of the effects of different types of inhomogeneities on spiral-wave dynamics in mathematical models for cardiac tissue (see, e.g., Refs. [6,7,21,22] and references therein), to the best of our knowledge there has been no study, based on Markov-state models, of an inhomogeneity comprising mutant myocytes in a background of wild-type myocytes. Therefore, we present a representative study of spiral-wave dynamics in the presence of a clump of only mutant cells, of radius R = 1.125 cm, embedded in a background of wild-type cells. We then explore the possibility of spiral-wave formation via high-frequency stimulation by pacing the simulation domain from the left boundary (pacing frequency 3.7Hz). We find that, in the MM1 MUT model, no spiral wave forms (Figure 11(A)); by contrast, in the MM2 MUT model a spiral wave forms (Figure 11(B)).
[The spatiotemporal evolution of these waves is shown in the Movie(M10) in the Supplementary Material 4.] This qualitative difference arises because of the different types of EADs (discussed above) in MM1 MUT and MM2 MUT models.

3D results
We end with an illustrative study of scroll waves in 3D TP06, MM1 WT, and MM2 WT models. These waves are shown via color isosurface plots of the transmembrane potential V m in Figures 12 A, B, and C, respectively, for both a homogeneous domain (top panel) and with localized obstacles [P f = 10% (middle panel) and P f = 50% (bottom panel)]. In a homogeneous domain, scroll waves are stable in TP06 and MM1 WT models, but not in the MM2 model. An inexcitable obstacle, with P f = 10%, has no significant impact on scroll waves in TP06 and MM1 WT models. An increase in P f , say to P f = 50%, leads to an anchoring of the scroll waves at the obstacle (as in the study of inexcitable obstacles in Ref. [9]). For the MM2 model, with P f = 10%, scroll-wave break-up is enhanced; but for P f = 50%, the scroll wave gets anchored to the obstacle. The spatiotemporal evolution of these scroll waves is shown in the Movie(M10) in the Supplementary Material 4.

Discussion and Conclusions
Earlier studies of Markov models for cardiac myocytes have focused on the effects of mutations in subunits of N a + , and K + channels in the context of the Brugada and LQT syndromes [28,15,14,29]; these studies have elucidated the effects of changes in the kinetic properties of these ion-channel, and their consequences, such as the prolongation of the APD, which leads, in turn, to EADs. Also, Markov models have been used to asses the importance of a particular functionality of an ion channel, e.g., the role of I Ks on AP repolarisation [30]. Furthermore, Markov models have been used to investigate the activation and inactivation properties of N a + , K + and Ca +2 ion channels as, e.g., in [31,32,33,34]. In addition, some studies have focused on theraupetics and drug-channel interactions [18,3,17], from cellular to the anatomically realistic tissue levels.
We have investigated, from cellular to tissue levels, the differences in kinetic properties of Na channels in TP06, MM1 (WT and MUT), and MM2 (WT and MUT) models. We have shown that Na channels in TP06 and MM2 (WT and MUT) models are activated faster, with respect to V m , than their counterparts in the MM1 (WT and MUT) models; also the inactivation of these channels is faster in the TP06 model than in MM1 (WT and MUT) and MM2 (WT and MUT) models. These differences leads to different times of openings of the Na channel and the amplitudes of I N a,f are completeley determined by the amplitude of P O in the MM models. These changes in the amplitudes of I N a,f , I N a,L , and the activation-inactivation dynamics lead to disparate CVR and maximal CVs in cable simulations. To the best of our knowledge, our study is the first to compare spiral-wave dynamics in different Markov models for the Na (WT and MUT) ion channels in realistic mathematical models for cardiac tissue. We have carried out in silico studies, in both homogeneous simulation domains and domains with inhomogeneities, to compare and contrast spiral-and scroll-wave dynamics in five different models for cardiac tissue (Hodgkin-Huxley type, TP06 model [1], and Markov-state models such as MM1 WT and MM2 WT, for the WT Na channel, and MM1 MUT and MM2 MUT, for the mutant Na channel [2,3]). Our study explores the sensitive dependence of spiral-and scroll-wave dynamics on these five models and the parameters that define them. We also examine the control of spiral-wave turbulence in these models. To the best of our knowledge, such a comparative study of wave dynamics in HHM and Markov-state models has not been carried out hitherto. In our opinion, such a comparison is even more valuable than the comparison of single-cell properties of models for cardiac myocytes. We hope our study will lead to more comparisons of wave dynamics in different mathematical models for cardiac tissue and in in vitro experiments. Furthermore, we have carried out a detailed parameter-sensitivity study, principally for the WT models by using multivariable linear regression (see the Supplementary material).
We mention some of the limitations of our study. We have studied the differences in Na ion-channel modeling, which is important in the context of the LQT syndrome; but we have not carried out such a study for Kr-channel modeling, as mutations in the Kr channel also leads to the LQT syndrome. We have considered a few, illustrative Markov-state models for Na channels, given our computational resources; many more such Markov-state models have been developed for the Na channel [16,35]; a comprehensive comparison of all these models lies beyond the scope of this paper. Also, we have used a single base model, i.e., TP06 model upon which we build the Markov-state models (MM1 WT, MM1 MUT, MM2 WT and MM2 MUT); other base models, e.g., the O'Hara Rudy model [36], can be used; a comprehensive comparison of all these models lies beyond the scope of this paper. We do not use an anatomically realistic simulation domain [37], with information about the orientation of muscle fibers [9,38]; and we use a monodomain description for cardiac tissue. These considerations lie beyond the scope of this paper. We note, though, that the study of Ref. [39,40] has compared results from monodomain and bidomain models and has shown that the differences between them are small. I N a fast inward N a + current I CaL L-type inward Ca ++ current I to Transient outward current I Ks Slow delayed rectifier outward K + current I Kr Rapid delayed rectifier outward K + current I K1 Inward rectifier outward K + current I N aCa N a + /Ca ++ exchanger current I N aK N a + /K + pump current I pCa plateau Ca ++ current I pK plateau K + current I bN a background inward N a + current I bCa background inward Ca + current Table 1: The various ionic currents in the TP06 model Ref. [1]; the symbols used for the currents follow Ref. [1].         , for different P f . If ∆T > 0, then the frequency ω ∼ T −1 , for a given pair (R, P f ), is less than ω 0 ∼ T −1 0 (for R, P f = 100%); this may occur because the spiral core penetrates the obstacle because of a spanning cluster of excitable regions inside the obstacle. In (C) and (E) we show different regions in the (R, P f ) plane for TP06 and MM1 WT models, respectively, with the following the color code: Light Blue: an increase in ω relative to ω 0 (caused by penetration of the spiral core). Light Green: decrease in ω relative to ω 0 (accompanied by penetration of the spiral core); (R, P f ). Yellow: No penetration of the spiral core into the obstacle. Dark Blue: No change in ω relative to ω 0 even though the spiral core penetrates into the obstacle.