article: Elucidating the underlying mechanisms of fatal cardiac arrhythmias requires a tight integration of electrophysiological experiments, models, and theory. Existing models of transmembrane action potential (AP) are complex (resulting in over parameterization) and varied (leading to dissimilar predictions). Thus, simpler models are needed to elucidate the “minimal physiological requirements” to reproduce significant observable phenomena using as few parameters as possible. Moreover, models have been derived from experimental studies from a variety of species under a range of environmental conditions (for example, all existing rabbit AP models incorporate a formulation of the rapid sodium current, INa, based on 30 year old data from chick embryo cell aggregates). Here we develop a simple “parsimonious” rabbit AP model that is mathematically identifiable (i. e. , not over parameterized) by combining a novel Hodgkin-Huxley formulation of INa with a phenomenological model of repolarization similar to the voltage dependent, time-independent rectifying outward potassium current (IK). The model was calibrated using the following experimental data sets measured from the same species (rabbit) under physiological conditions: dynamic current-voltage (I-V) relationships during the AP upstroke; rapid recovery of AP excitability during the relative refractory period; and steady-state INa inactivation via voltage clamp. Simulations reproduced several important “emergent” phenomena including cellular alternans at rates > 250 bpm as observed in rabbit myocytes, reentrant spiral waves as observed on the surface of the rabbit heart, and spiral wave breakup. Model variants were studied which elucidated the minimal requirements for alternans and spiral wave break up, namely the kinetics of INa inactivation and the non-linear rectification of IK. The simplicity of the model, and the fact that its parameters have physiological meaning, make it ideal for engendering generalizable mechanistic insight and should provide a solid “building-block” to generate more detailed ionic models to represent complex rabbit electrophysiology.
  The simulation of action potentials (APs) and their propagation has been an integral part of the field of electrophysiology for over half a century. [1] Over a hundred cellular AP models from many regions of the heart encompassing a variety of species have been published. These “ionic” models have been derived almost exclusively from voltage clamp data from a variety of species, often recorded under conditions that are not physiological. Ionic models are increasingly complex, and are comprised of numerous sub-models (mostly transmembrane currents), and contain tens of variables and hundreds of parameters. For example, a recent meta-analysis found that 50% of the data used in the development of the ten-Tusscher-Panfilov [2] and Iyer [3] human ventricular models were of non-human origin [4]. While their complexity has many advantages, there are a variety of limitations including: over-parametrization[5]; non-uniqueness/multi-stability [6,7]; dissimilar predictions [8–12]; limited validation outside their development domain [13] (e. g. using a cell model derived from voltage clamp data to study reentrant arrhythmias); difficulty in representing natural “physiological redundancy” (e. g. repolarization reserve [14]); interpreting the relative roles of model sub-components on model predictions; and generalizing insight from “model-specific” simulations results. An alternative, and complementary approach compared to these ionic models, is to use “phenomenological” models which are designed to represent specific, often macroscopic phenomena (e. g. , rate dependence of action potential duration) [15–19]. Phenomenological models are designed to reproduce one (or two) specific phenomena (on) very well, and are: simplistic, computationally efficient, and sometimes amenable to analytical approaches. However, unlike ionic models, phenomenological models do not provide a direct link to physiologically meaningful model parameters derived experimentally. Therefore phenomenological models provide only limited mechanistic insight, and are not amenable for reproducing numerous phenomena, nor can they be easily extended. The two most important membrane currents for cardiac excitability are the rapid sodium current (INa) and the rectifying potassium current (IK1), which are both highly non-linear functions of transmembrane potential (Vm). INa is responsible for the rapid all-or-none depolarization of the AP upstroke, and IK1 is responsible for maintaining the resting potential near -85 mV. In fact, it has been shown recently that transfecting cells with genes encoding for cardiac INa, IK1, and gap junctions Cx43 (only) into unexcitable somatic cells can transform monolayers into excitable tissue capable of supporting propagating “cardiac” action potentials including reentrant “spiral” waves! [20] Due to the ethical and practical limitations involved in studying human physiology, animal experiments are necessary to understand and develop clinical strategies aimed at treating the pathological mechanisms underlying heart disease. The isolated rabbit heart has been shown to share certain electrophysiological characteristics as the human, including repolarization reserve [21] and ventricular fibrillation dynamics. [22] as well as the effects of many pharmaceuticals. [23] All existing ionic models of the rabbit ventricular AP [24–27] incorporate the INa equations from the LR1 guinea pig model, [28] albeit with different conductance values. The LR1 INa equations incorporate the gating kinetics of fast activation, fast inactivation and slow inactivation. The fast activation and fast inactivation equations were derived from voltage clamp (VC) data from spherical clusters of eleven day-old embryonic chick heart cells, [29] and slow inactivation parameters were based on VC results from sheep and pig ventricular trabeculae. [30] Obviously, there is a need to update rabbit AP models to incorporate an INa model based on rabbit data, but this is not a simple task. To date, INa activation parameters have been derived exclusively from VC protocols; however, since the INa current in the adult myocyte is too large to allow adequate voltage control under physiological conditions, [31] these VC protocols are carried out at low temperatures with low extracellular sodium concentrations or in HEK cells or oocyctes transfected with the SCN5A gene which encodes the cardiac INa channel (Nav 1. 5). The first aim of this paper is to develop a “parsimonious” model for the rabbit myocyte AP based only on data recorded from the rabbit under physiological conditions. A model is considered parsimonious if it accomplishes a desired level of explanation or prediction with as few parameters as possible (although, as far as we are aware, there is no definitive and scientifically rigorous definition of parsimonious). For this paper, the desired level of prediction is the following well-known and important electro-physiological phenomena, measured from rabbit ventricular myocytes/tissue under nearly identical and physiological conditions: 1) steady-state inactivation as determined from voltage clamp experiments [31]; 2) action potential depolarization in single cells; 3) recovery of AP excitability in single cells[32]; and 4) action potential depolarization dynamics during propagation in the whole heart. [33] For calibration and evaluation of the model we coupled this INa model with various models of IK1, and also explored a range of cellular and tissue behavior in an AP model incorporating a previously published phenomenological two parameter model of repolarization. [33] The overall model will be referred to as the parsimonious rabbit (PR) model. The second aim of this paper is to demonstrate that this simple and parsimonious model is sufficient for reproducing important “emergent” complex phenomena such as alternans and spiral wave breakup, which suggests a novel physiological mechanism for arrhythmia instabilities.
  The equations of INa are of the form pioneered by Hodgkin-Huxley (HH) [1] INa=gNam3h (Vm−ENa), (1) where Vm is the transmembrane potential, gNa is the maximal conductance of INa, m (fast activation) are gating variables, and h (fast inactivation), and ENa is the Nernst equilibrium potential for sodium. Beeler and Reuter [34] introduced a slow inactivation gate (j) INa=gNam3hj (Vm−ENa), (2) which has been incorporated into all of the most recent cardiac HH INa models. [24–27]. The HH equations for the gating variables are of the form: dydt=αy (1−y) −βyy≡y∞−yτy, (3) where y represents each gating variable, αy (Vm) and βy (Vm) represent the on and off rate constants for gating (respectively) which are voltage dependent, y∞ (Vm) and τy (Vm) are the steady state fraction of activation or inactivation and time constant (respectively) and are functions of Vm. Hodgkin and Huxley [1] fit their voltage clamp data to empirically determine the functions αy (Vm) and βy (Vm) using only 9 parameters; in contrast the LR1 INa equations contain 31 parameters. Since 1952 there has been significant advancement in the derivation of gating equations based on first principles. [35,36] The simplest functional form for the rate equations of voltage dependent gating based on thermodynamics and chemical reaction rates are: [37] αy=αy0exp (−δyVmky), βy=βy0exp ( (1−δy) Vmky), (4) where the constants αy0, βy0, and 0 ≤ δy ≤ 1 are positive, while ky < 0 for activation and ky > 0 for inactivation. Using Eq 3 we can derive the equations for y∞ (Vm) and τy (Vm) from Eq 4: y∞ (Vm) =11+exp[ (Vm−Ey) ky], (5) and τy (Vm) =2τy0exp[δy (Vm−Ey) ky]1+exp[ (Vm−Ey) ky]. (6) There are 4 parameters per gate (Ey, ky, τy0, δy) in this formulization. The main limitation of this formulization is that τm (Vm) →0 away from Em which is not realistic, because gating transformations cannot be instantaneous. For example, the LR1 equations predict τm (−85mV) = 0. 0055 ms and τm (+30mV) = 0. 041 ms at 37 C, however these values are below the resolution of voltage clamp measurements. We address this problem by removing the Vm dependence of τm (making it a constant, thus eliminating parameter δm). To summarize, we define the novel PR INa submodel using Eqs 1 and 3, with m∞ and h∞ given by Eq 5, τh given by Eq 6, and τm equal to a constant. All model parameters are provided in Table 1. Values for Eh and kh where taken directly from Table 1 in [31]. We set ENa = 65mV, and all other parameters were determined via manual fitting to the experimental data described below using no more than 2 significant digits; τhmax and δh were computed via “inverting” Eq 6 given the two value pairs τh (−80mV) = 4. 0ms and τh (0mV) = 0. 45ms. Simulations required for the calibration were carried out in single cells and a one-dimensional cable by integrating the three ordinary differential equations using exact integration for the gating variables and forward Euler integration for Vm with dt = 0. 001 ms which is actually a very conservative choice made because of the low computational cost. The fact that τm is constant provides for a well-defined minimum time constant of the model which allows for the use of larger values of dt and computational efficiency compared to other models. [38] Cable simulations were performed using a central difference approximation of the Laplacian in a 4 cm long cable using Iion = INa + IK1L (in order to establish a steady-state resting membrane potential) with diffusion coefficient (D) equal to 0. 001 cm2/ms and dx = 0. 01 cm, where IK1L is the formulation of IK1 provided by Livshitz and Rudy [39]. INa and Vm signals were recorded at the center of the cable to construct the I-V curves for calibration to the experimental data as described below. The current during the AP upstroke is predominately INa, therefore the precise time course of Vm (t) contains a wealth of information regarding INa kinetics. [40–42] In isolated myocytes the time course of INa can be estimated as Iion=−CmdVmdt (after the stimulus current is turned off), where Cm is the specific membrane capacitance which we assume is equal to 1 μF/cm2. During propagation, however, the total transmembrane current is not proportional to dVmdt; but we recently developed a methodology and justified the computation of INa as Cm (D (CV) 2d2Vmdt2−dVmdt) during the AP upstroke during planar wave propagation, where CV is the conduction velocity [33]. These relationships were used to calibrate our model to dynamic I-V curves from glass microelectrode Vm measurements during AP upstrokes from isolated rabbit ventricular myocytes[43] and from the epicardial surface during planar propagation on the surface of the isolated Langendorrff-perfused heart. [33] In isolated myocytes Vm traces from APs stimulated with near threshold stimuli with a latency of 2. 5 ms (range of 1. 2 to 3. 7 ms) were chosen from a previous study[43], because this latency provided datum during the early decrease INa of near –40 mV. Similar to Joyner et al. [32] we found a slight decrease (5. 2 ± 0. 5 mV/ms per ms) of (dVmdt) max as a function of latency which corresponds to only a 3. 4% decrease compared to the traditional latency value of 1 ms. It is well known that (dVmdt) max is decreased during the relative refractory period of the AP, and this is thought to be a cellular phenomenon [32]. Here we assume that this phenomena (“recovery of AP excitability”) is controlled only by the voltage and time dependent INa current and the AP shape. To reproduce the experimental results of recovery of AP excitability in isolated myocytes, [32] we developed a novel “virtual protocol”. This protocol involves simulating an action potential clamp of INa using APs measured during pacing from six isolated myoctes recorded in a previous study. [43] During the 11th beat the simulations were switched to current-clamp mode and stimulated with a 2 ms stimulus at various coupling intervals; the amplitude was adjusted such that the latency was 1 ms. The resulting peak INa was normalized to that recorded during the 10th beat and presented as a function of coupling interval. Since the INa equations in all previous rabbit models are identical to LR1, except for scaling (conductances), we present previous model predictions using LR1 with two conductance values (23 mS/μF and 8 mS/μF) to cover the entire range (grey lines in figures). In addition we present results from the Ebihara-Johnson (EJ) model [29] (dashed grey lines) which is identical to LR1, except it does not contain slow inactivation. To create our parsimonious rabbit (PR) cell model, we combined our model of INa with a previously published phenomenological model of repolarization (IK) designed to reproduce rabbit AP repolarization during propagation: IK=gKe−b (Vm−EK) (Vm−EK), (7) where gK is the maximal conductance of IK, b is a parameter controlling AP shape, and EK is the reversal potential for potassium (set equal to the average resting potential of the six cells in Fig 1 of -83 mV); the nominal parameters values for IK are provided in Table 1, as previously determined. [33] We define our PR cell model via a total ionic current of INa + IK in Eqs 1 and 5–7 with parameters in Table 1.
  The experimental traces of Vm during AP depolarization of isolated rabbit ventricular myocytes are shown in Fig 1A; the traces are aligned to (dVmdt) max. The average dynamic I-V curve during the AP upstroke for these traces are plotted as symbols (● representing the mean ± standard error of the mean) in Fig 1B. The calibrated PR results are shown as a thick black line. Previous model predictions (grey lines) bracket the range of peak INa but do not accurately capture the initial inward current deflection of the experimental I-V relationship (see box). The experimental dynamic I-V curves during the AP upstroke during propagation from [33] are shown as symbols (●) in Fig 2A. Similar to above, previous model predictions (grey lines) span the range of peak INa but do not accurately capture the initial inward current deflection of the experimental I-V curve (see box). Although the PR model fits the initial decrease of INa and maximum Vm well, it does not match the experimental dynamic I-V curve during the second half of the upstroke. It is important to note that the experimental dynamic I-V curves for isolated myocytes (Fig 1B) is different than that obtained during propagation (Fig 2A). To quantify these differences we fit both curves to fifth order polynomials and identified the regions in which the 95% confidence intervals of these fits did not overlap (see Fig A in S1 Text). The inward current during the upstroke for myocytes was more negative compared to during propagation for -41 mV < Vm < -22 mV and for Vm > +4. 5 mV. The normalized magnitude of peak INa as a function of recovery time are shown in Fig 3B; experimental results from [32] are shown as symbols (●), PR model results are shown as a thick black line, and previous model predictions shown in grey. Simulation results (i. e. , lines) include standard error bars because the virtual protocol includes six different APs and thus include the effect of normal variation of AP shape. Unlike previous models with slow inactivation (Eq 2), the PR and EJ80 models which include only fast inactivation (Eq 1) are consistent with the experiments. The voltage dependence of steady state INa activation and inactivation as well as the corresponding time constants are shown in Fig 4 for our PR model, as well as LR1, and EJ models. Our PR model exhibits a steeper activation slope at slightly less depolarized voltages compared to LR1 and EJ (panel A). As described above, τm in the PR model does not exhibit voltage dependence but is similar to the mean value of LR1 and EJ (panel B). The steady state INa inactivation for PR is more negative compared to LR1 and EJ (panel C) and the maximum value of τh for PR is less than that for τh for LR1 and EJ80 and is considerably less than τj for LR1 (panel D). The PR results in Figs 1 and 2 represent the results from for our novel INa sub-model only; the results in Fig 3 include coupling to the IK1L sub-model. To assess the effect of repolarization current on the PR results shown in Figs 1–3 we reran all simulations with cell models incorporating IK1L or the phenomenological model of repolarization Eq 7 (PR cell model). The results from these simulations were nearly identical to the PR INa sub-model only; in fact the results are shown as thick dashed (IK1L) and dotted (IK) lines in Figs 1–3 but cannot be seen because they superimpose upon the solid black line. Simulations using our PR cell model resulted in APs with the following characteristics: INamin=−233μA/cm2, Vmmax=34mV and APD = 197ms for single cells; and INamin=−289μA/cm2, Vmmax=24mV, CV = 55 cm / s and APD = 142ms for propagation in a cable. With the exception of APD, these values were not significantly altered (< 2%) compared to simulations using IK1L [39] or IK [33] indicating that the predictions of our PR INa sub-model are quite insensitive to the specific choice of repolarization current. We performed parameter sensitivity analysis by quantifying the effect of a 1% variation of the eight INa gating parameters on eight quantities of interest; results are provided in the Table A in S1 Text as an 8x8 matrix. We also characterized the effect of AP shape on cell and tissue-level behavior for our new PR model (INa + IK) by varying repolarization parameters gK and b over the entire “physiological range” for all mammals: reported values of gK range from 0. 1 to 0. 5 mS/μF and the range of b values chosen, 0. 03 ≤ b ≤ 0. 05 mV-1, correspond to APDs ranging from 23 ms to 516 ms. The results are displayed in the Figs B and C in S1 Text. Over this range of gK and b, INamin varied 1. 6% and Vmmax varied 1. 2% in single cells, while during propagation the following quantities exhibited a fairly small dependence on IK parameters: INamin: 1. 8%; Vmmax: 7. 5%; CV: 4. 5%. The fact that action potential depolarization characteristics are fairly insensitive to repolarization in single cells and during propagation provides support for the robustness of our INa model. To confirm that the data used in the model calibration is representative, we compared several model results to experimental data from previously published studies. Table B in S1 Text demonstrates that the following model results are within the range found in the literature: resting membrane potential, action potential amplitude, (dVm/dt) max, and CV. In addition, we compared the predictions of spiral waves in the PR and Mahajan (M08) models (5 cm x 5 cm) to those recorded in the whole rabbit heart by Schalij et al. [44] Schalij et al found that only fibrillation was observed in the intact heart but that, after freezing the inside of the heart, stable reentry occurred around a linear line of conduction block of ~ 2 cm with a cycle length of 130 ± 11 ms. [44] Simulation results for PR and M08 are shown in Fig 5. The cycle length (~ 150 ms) and line of block (~ 3. 2 cm) for PR using the nominal parameter values were larger than the experiments. Incidentally, decreasing parameter b by only 4% decrease improved the correspondence considerably (cycle length: ~125 ms; line of block: 2. 5 cm as shown in the middle, right panel of Fig 6 below). This change is justified because the PR model does not include any repolarization kinetics and is not developed to reproduce APD shortening with rate. In contrast to the experimental results of Schalij et al. , [44] the line of block in our simulations was not static, but rotated which is typical of spiral wave simulations. [2] Simulations with M08 resulted in non-sustained unstable reentry with average cycle length of length (~ 150 ms); the tip trajectory traversed the majority of the region, and eventually terminated by hitting the top boundary. There was evidence of wave break (see green and blue trajectories in Fig 5B consistent with previous results indicating SWB for this model in a 7. 5 cm x 7. 5 cm region. The level of validation required for a model depends on its “context of use” and the consequences of incorrect model predictions[13]. Hence further validation is expected to be required for each specific context. The preliminary validation performed here is appropriate for using the model to gain insight into underlying physiological mechanisms as demonstrated below. During “functional” spiral wave reentry the depolarization and repolarization processes are interdependent. As shown in Fig 6, the dynamics of spiral waves in our model varied considerably as a function of repolarization (i. e. , gK and b). Spiral waves were not stable (i. e. , rotationally symmetric) for any parameter pairs (gK, b) in the range studied. For 8 of 9 parameter pairs, the spiral wave tip followed typical non-circular trajectories, while spiral wave breakup (SWB) occurred for one parameter set (gK = 0. 5, b = 0. 035). The frequency of activity for single spiral waves ranged from ~3 Hz to ~18 Hz and the frequency of activity during SWB was ~30 Hz. Since SWB was observed (albeit for parameters corresponding to APD = 23ms which is not consistent with normal rabbit physiology) we investigated whether cellular alternans occurred in our PR model. We paced the PR cell model with a 2 ms stimuli of -20 pA/pF and found that alternans occurred for cycle lengths between 202 and 210 ms, as shown in Fig 7, which is within the range (190–240 ms) reported previously for isolated rabbit myocytes. [26]. Since the PR model contains no repolarization kinetics, we hypothesized that APD alternans occurred as a result of alternations of peak INa, which resulted in alternans in Vmmax, which caused alternans in APD. Since INa is negligible except during depolarization, APD is a function of only Vmmax (and can be computed by integrating Eq 7 with initial condition Vm=Vmmax). Since alternations in the variable h were more pronounced than m (see Fig 7B), we hypothesized that only INa inactivation kinetics are necessary for cellular alternans. A “reduced” two variable (Vm, h) model, in which variable m3 was replaced with function m∞3 (Vm) and gNa was reduced to 5. 8 mS/μF (to maintain CV = 55cm / s), confirmed that INa inactivation by itself was capable of generating cellular alternans in our PR model.
  Traditionally, cardiac AP models have been developed from equations and parameters derived from voltage clamp experiments. Many currents have voltage-dependent gating variables whose characterization via voltage clamp provides intuitive and meaningful parameters. However, in order to obtain such parameters it is necessary to isolate individual currents either chemically, functionally, or genetically. In addition, the relationship of voltage-clamp derived parameters to AP characteristics in single cells and during propagation is complex, non-intuitive, and often involves the interaction of multiple variables. Here we provide an alternative approach in which we calibrate our model in such a way as to reproduce certain AP behavior under physiological conditions while including all available voltage clamp data. In this way, our PR model is both physiologically meaningful and parsimonious. It is physiological because it incorporates HH-like equations and parameters with only one phenomenological parameter b which controls AP shape. It is parsimonious because it was constructed was constructed to reproduce the specific set of experimental data sets (Figs 1–3) that capture important electrophysiological behavior with as few parameters as possible. In addition, we have recently shown that this model is mathematically identifiable (i. e. , not over parameterized) [45]. Incidentally, as described by Biktashev et al. [46] FitzHugh–Nagumo models fail to reproduce some features of cardiac ionic models such as: slow repolarization, slow sub-threshold response, non-Tikhonov asymptotic properties of excitability [47], wave front dissipation [48], and different action potential amplitude in single cells versus propagation [9] which our model captures (see Fig 2). Physiological models are developed using a set of experimental calibration data, which for our model is for rabbit tissue and myocytes at physiological temperature and heart rate. Our approach should be relatively straight-forward to implement for other species, although this will require the high-fidelity measurement of the dynamic I-V curve during propagation [33] and in single cells as well as the recovery of excitability in isolated myocytes [32]. As far as we are aware, our PR model is the first to include dynamic I-V curves during propagation into its development (i. e. , calibration) which we believe is superior to simply adjusting gNa to match CV. We believe that proper extensions of our PR model will include additional constraints obtained using appropriate data sets during recalibration. One reason we think our model is useful is because it has only one inward current which is non-zero only during depolarization, and only one outward current. Unlike traditional models in which multiple combinations of currents can give rise to the same AP shape [49], in our model repolarization is characterized only by a time-independent, voltage dependent current (IK), and therefore, IK parameters completely determines AP shape [33]. In addition, each current typically incorporates one or two variables (often gates). Despite the non-uniqueness of the relationship between model parameters and AP shape, model developers have largely ignored the interaction of the various currents and variables (components) in simulations. For a model with n > 1 components, the total number of possible interactions between two or more components, excluding the full model of all n components, is 2n − n − 2. Traditional models have at least 10 components, which results in over 1000 interactions! Characterizing these interactions appropriately both physiologically and mathematically continues to be a significant challenge. [2,31,50,51]The simplicity of our model allows unique physiological insight regarding the interaction of depolarization and repolarization processes. We found that IK conductance (gK) is related to the cycle length of reentry (see Fig 6) consistent with the experimental findings of Warren et al. in guinea pig [52]. In addition, our PR model exhibited the expected relationship between APD and spiral wave behavior [53]; specifically, the frequency of activity for single spiral was inversely related to APD, and also inversely related to the region covered by the spiral wave tip. Surprisingly, our PR model (comprised of only a HH INa current with only two gating variables and a time-independent, voltage-dependent IK current) exhibited important “emergent” behavior including cellular alternans (Fig 7), unstable spiral waves (Figs 5 and 6), and even spiral wave breakup (Fig 6, bottom left panel). We initially thought that alternans and spiral wave breakup would require dynamic repolarization variables with time constants with values close to the APD capable of generating APD alternans. We only discovered spiral breakup serendipitously while performing a sensitivity analysis of the spiral wave meandering behavior on the IK repolarization parameters. It is important to note, however, that the parameter set for which breakup occurred (gK = 0. 5, b = 0. 035) resulted in a rate (~30 Hz) which is twice that observed in experiments. [54,55] While our PR model reproduces important physiologically meaningful behavior, it does have limitations. It should be noted that the alternans resulting from the interaction of INa and IK are the ‘minimal’ requirements, other currents such as Ito and ICaL may either enhance or suppress this behavior. Our PR model does not include intracellular calcium dynamics, nor the majority of repolarization membrane currents. It should be noted that the Tyrode’s solution used to perfuse the hearts from which the propagation I-V experimental data were obtained (Fig 2) in the presence of 15 mM diacetyl monoxime, however the longitudinal propagation velocity was 59 ± 3 cm/s which is within the range reported for the rabbit without any uncoupling agents (see Table B in S1 Text). The INa equations are of HH type, Markov models are required to replicate certain features of voltage clamp experiments including drug binding kinetics. [50] Importantly, the APD in our model is dependent on Vmmax which depends on the stimulus current, which is different in single cells compared to in tissue during propagation. When the PR model is stimulated to unphysiologically large Vm the APD is unphysiologically long; to avoid this in tissue simulations we initiated propagation in the cable by holding Vm at the end (1 mm) at 0 mV for 2 ms. One significant difference between our INa sub-model and previous rabbit INa sub-models is that we did not include a slow inactivation gate. While slow inactivation of INa can be justified from voltage clamp experiments [31], we have found that it can induce unphysiological post repolarization refractoriness as shown in Fig 3 and [51]. In addition, the INa sub-model in the existing rabbit models are the same as LR1 in which the slow inactivation kinetics were derived from sheep and calf trabeculae muscles in which spatial homogeneity of Vm and VC control are lacking. We believe that much more research is needed to appropriately incorporate slow inactivation of INa (as well as late INa) into AP models, especially those used to simulate propagation. Incidentally, ten Tusscher et al. interpolated experimental time constants of fast and slow activation as a function of voltage over the range -80 to -40 mV in which data is not available to values more than four times values outside of this range (see Fig 1 D & E in [2]). Our novel PR model provides unique and physiological generalizable insight into the relationship among gating kinetics and AP behavior in single cells and during propagation. We show unequivocally that only two currents (INa, IK) and three variables (Vm, m, h) are required to reproduce cellular alternans and spiral wave breakup. In fact, only two variables (Vm, h) are required to produce alternans and implicate a destabilizing effect resulting from the interaction of the non-linear rectification of IK [56] and the voltage and time dependence of h. We believe that our PR model is an ideal tool to quantify excitability and propagation related to a variety of important physiological and pathophysiological issues. For example our model can be used to study the effects of drugs (our model predicts that an 83% decrease of gNa is required for propagation failure), genetic mutations of SCN5A such as LQT3, and the differences between atrial and ventricular tissue. In addition we believe that it provides a solid “building-block” to develop more realistic rabbit AP models by including additional sub-models (e. g. , ICaL, Ito, IKr and intracellular calcium handling). All experimental data and PR simulation results presented in this paper (with the exception of the videos) will be made available online at Research Gate. The model will be uploaded to CellML.
summary: Understanding and preventing life-threatening irregular electrical heart rhythms includes basic experimental, numerical, and theoretical research. Computer models of the electrical dynamics of cardiac cells and impulse propagation throughout the heart are essential tools of this research. For example, simulations of electrical activity such as rotating ‘spiral waves’ have been used to understand how irregular heart rhythms are maintained. However, existing models are derived exclusively from sub-cellular data under a variety of environmental conditions and species. These models tend to be exceedingly complex including hundreds of variables and parameters which make them difficult to validate and analyze. Using experimental data from one species (rabbit) under nearly identical physiological conditions, we developed a simple model of the electrical activity of a cardiac cell derived from sub-cellular, cellular, and tissue experimental data. This model reproduces cellular excitability and its recovery as well as several important “emergent” phenomena including beat-to-beat cellular alterations and unstable spiral waves. Under some conditions, unstable spiral waves in this model give rise to continuous formation of new spiral waves (i. e. , “spiral wave breakup”) which is thought to be the underlying cause of cardiac fibrillation and sudden cardiac death.
keywords: action potentials
  medicine and health sciences
  muscle tissue
  depolarization
  membrane potential
  vertebrates
  rabbits
  electrophysiology
  neuroscience
  animals
  mammals
  reproductive physiology
  simulation and modeling
  animal models
  model organisms
  physiological parameters
  research and analysis methods
  animal cells
  biological tissue
  muscle cells
  cell biology
  anatomy
  physiology
  biology and life sciences
  cellular types
  amniotes
  neurophysiology
  organisms
title: A Parsimonious Model of the Rabbit Action Potential Elucidates the Minimal Physiological Requirements for Alternans and Spiral Wave Breakup