5
views
0
recommends
+1 Recommend
0 collections
    0
    shares
      • Record: found
      • Abstract: found
      • Article: found
      Is Open Access

      Dissecting the roles of calcium cycling and its coupling with voltage in the genesis of early afterdepolarizations in cardiac myocyte models

      research-article
      1 , 2 , 1 , * ,
      PLOS Computational Biology
      Public Library of Science

      Read this article at

      Bookmark
          There is no author summary for this article yet. Authors can add summaries to their articles on ScienceOpen to make them more accessible to a non-specialist audience.

          Abstract

          Early afterdepolarizations (EADs) are abnormal depolarizations during the plateau phase of the action potential, which are known to be associated with lethal arrhythmias in the heart. There are two major hypotheses for EAD genesis based on experimental observations, i.e., the voltage ( V m )-driven and intracellular calcium (Ca)-driven mechanisms. In ventricular myocytes, Ca and V m are bidirectionally coupled, which can affect each other’s dynamics and result in new dynamics, however, the roles of Ca cycling and its coupling with V m in the genesis of EADs have not been well understood. In this study, we use an action potential model that is capable of independent V m and Ca oscillations to investigate the roles of V m and Ca coupling in EAD genesis. Four different mechanisms of EADs are identified, which are either driven by V m oscillations or Ca oscillations alone, or oscillations caused by their interactions. We also use 5 other ventricular action potential models to assess these EAD mechanisms and show that EADs in these models are mainly V m -driven. These mechanistic insights from our simulations provide a theoretical base for understanding experimentally observed EADs and EAD-related arrhythmogenesis.

          Author summary

          Early afterdepolarizations (EADs) are dangerous abnormal electrical activities in the heart, which may cause lethal arrhythmias. Although EADs have been widely investigated in both experimental and computer simulation studies, their mechanisms remain incompletely understood. In the present work, we carry out computer simulations using action potential models with detailed formulations of ionic currents and calcium cycling to investigate the genesis of EADs. Different mechanisms and causes are identified in our simulations, which agree with experimental observations. The mechanistic insights from our work provide a theoretical base for understanding the mechanisms of EADs and EAD-related arrhythmogenesis.

          Related collections

          Most cited references58

          • Record: found
          • Abstract: found
          • Article: not found

          hERG potassium channels and cardiac arrhythmia.

          hERG potassium channels are essential for normal electrical activity in the heart. Inherited mutations in the HERG gene cause long QT syndrome, a disorder that predisposes individuals to life-threatening arrhythmias. Arrhythmia can also be induced by a blockage of hERG channels by a surprisingly diverse group of drugs. This side effect is a common reason for drug failure in preclinical safety trials. Insights gained from the crystal structures of other potassium channels have helped our understanding of the block of hERG channels and the mechanisms of gating.
            Bookmark
            • Record: found
            • Abstract: found
            • Article: found
            Is Open Access

            Simulation of the Undiseased Human Cardiac Ventricular Action Potential: Model Formulation and Experimental Validation

            Introduction The first step toward preventing sudden cardiac death is understanding the basic mechanisms of ventricular arrhythmias at the level of ion channel currents and the single myocyte action potential (AP), using both experiments[1] and theoretical models[2]. Obtaining ventricular myocytes from human hearts for the study of arrhythmia mechanisms is both rare and technically challenging. Consequently, these mechanisms are usually studied with human channels expressed in non myocytes, or with non human (rodent or other mammalian) myocytes. However, these approaches have limitations, because functionally important accessory subunits and anchoring proteins native to ventricular myocytes[3] are absent in expression systems, and even among mammalian ventricular myocytes, ion channel kinetics[4], [5], [6] and consequently arrhythmia mechanisms are strongly species dependent. These issues limit the applicability of results from animal studies to human cardiac electrophysiology and clinical arrhythmia[7]. Measurements from undiseased human ventricular myocytes are a requisite for understanding human cell electrophysiology. Here, we present data from over 100 undiseased human hearts for steady state rate dependence, and restitution of the ventricular AP. Importantly, we also obtained essential new measurements for the L-type Ca2+ current, K+ currents, and Na+/Ca2+ exchange current from undiseased human ventricle. These previously unavailable data are critically important for correct formulation of mathematical models for simulation of electrophysiology and cellular arrhythmia mechanisms[8]. Using the new data together with previously published experiments, a detailed mathematical model of undiseased human ventricular myocyte electrophysiology and Ca2+ cycling was developed and thoroughly validated over the entire range of physiological frequencies. This model is referred to as the ORd (O'Hara-Rudy dynamic) model throughout the text. Model comparisons are conducted with the ten Tusscher-Panfilov (TP) model[9], and the Grandi-Bers (GB) model[10]. The ORd model was used to describe cellular electrophysiology mechanisms specific to human ventricular myocytes. Underlying mechanisms of AP duration (APD) rate dependence and APD restitution were investigated. The effects of Ca2+/calmodulin-dependent protein kinase II (CaMK) on known ionic current and Ca2+ cycling targets were incorporated and studied. Early afterdepolarizations (EADs) and alternans were reproduced by the model. These are important arrhythmogenic phenomena that must be reproduced in order to study the mechanisms of cardiac arrhythmias in human and simulate clinical interventions such as drugs. Results Throughout, new undiseased human ventricle experimental data are represented by white circles or white squares for isolated myocyte or small tissue preparation measurements, respectively. Previously published nonfailing human ventricle experimental data are shown with black symbols. Other data classification schemes are provided case by case in figure legends. Formulation, Validation and Properties of Simulated Currents: L-type Ca2+ Current (ICaL) Data for ICaL in the undiseased human ventricle are from Magyar et al.[11] and Fulop et al.[12] (both at 37°C, model validation in Figure 1C). Magyar et al. measured steady state activation, steady state inactivation, and the current voltage (I–V) curve. Fulop et al. measured recovery from inactivation. However, neither study separated Ca2+ dependent inactivation (CDI) from voltage dependent inactivation (VDI). In fact, no published data are available which separate CDI and VDI in the undiseased or nonfailing human ventricle. For this measurement, we made new recordings in undiseased human ventricular myocytes at 37°C (Figure 1, current traces and white circles). 10.1371/journal.pcbi.1002061.g001 Figure 1 Undiseased human ICaL experiments and model validation. A) Experiments: ICaL traces for currents carried by Ca2+ (top), Ba2+ (middle), and Sr2+ (bottom). The voltage protocol is below the Ca2+ traces. Ca2+ current decay was visibly more rapid than decay for Ba2+ or Sr2+ currents. Simulations: ICaL in response to the same voltage protocol with CDI (CDI+VDI, top), and without CDI (VDI-only, bottom). B) Experimental data are on the left (white circles, N = 5, from 3 hearts). Simulation results are on the right (solid lines). FRC is fractional remaining current. Times after peak current shown are from 5 to 55 ms, in 5 ms steps (indicated by arrow). Top left) Experiments showing the voltage and time dependence of FRC with Ba2+ as charge carrier (VDI only). Top right) Simulations of FRC, with n−gate = 0 (representing VDI only; see text and panel E). Bottom left) Experiments showing FRC with Ca2+ as charge carrier (CDI and VDI are concurrent). FRC for CDI+VDI was significantly smaller at more hyperpolarized potentials (Vm = −20 to 0 mV, dashed box) than FRC for VDI-alone. Bottom right) Simulations of FRC with free running n gate, allowing both CDI and VDI to occur. C) Data are from Magyar et al.[11] (black squares), Fulop et al.[12] (black diamonds), and previously unpublished (white circles, N = 5, from 3 hearts). Simulation results are solid lines. From left to right, top to bottom: steady state activation, steady state inactivation, fast time constant for VDI, slow time constant for VDI, relative weight of the fast component for VDI, I–V curve, experiments showing recovery from inactivation, and corresponding simulations. D) Human AP clamp waveform, used to elicit 1 µM nisoldipine sensitive current (ICaL, experiments, left) and comparison to simulations using the same AP clamp (right). E) Schematic diagram for the n gate, representing the fraction of L-type channels undergoing CDI. Calmodulin (CaM) is constitutively attached to the L-type channel. Ca2+ ions bind to CaM (on-rate k1 and off-rate k-1). With Ca2+ ions bound, the Ca2+/CaM/channel complex may activate CDI mode (asterisk and black color indicate CDI activation, on-rate k2 and off-rate k-2). Measurements were carried out with Ca2+ as charge carrier, allowing both CDI and VDI, or with Ba2+ as charge carrier, allowing only VDI. Results for Sr2+ were similar to those for Ba2+. From holding potential of −60 mV, 75 ms steps were to potentials ranging from −40 to +50 mV (10 mV increments, 3 second interpulse interval, Figure 1A). 75 ms was sufficient for comparison of CDI and VDI, since it is in the early phase of decay in which CDI effects are most prominent[13]. Simulated current traces for CDI+VDI and for VDI–alone were similar to the experiments. Fractional remaining current (FRC, at time t and voltage Vm, FRC(t,Vm) = I(t,Vm)/Ipeak(Vm)) quantified the voltage and time dependence of inactivation for comparison between charge carriers. Figure 1B compares FRC for Ba2+ (experiments top left, simulations right), and Ca2+ (experiments bottom left, simulations right). With Ba2+ as the charge carrier, FRC monotonically decreased with increasing voltage at all times after peak current. This finding is consistent with dependence of inactivation on voltage alone. In contrast, for Ca2+ currents, FRC did not decrease monotonically with increasing voltage. Rather, Ca2+ current FRC curves appear to be effectively voltage independent. FRC for CDI+VDI was statistically different from FRC for VDI-alone at the more hyperpolarized potentials (−20 to 0 mV, unpaired two-tailed t-test, p >τx2 at hyperpolarized potentials, where deactivation dominated, and τx2 200 bpm). The amplitude of APD alternans was ∼10 ms. These findings were reproduced by the human model (APD alternans of 11 ms at CL = 250 ms, Figure 14). Pacing at rates faster than 230 ms in the model caused 2 to 1 block (i.e. failed APs every other beat), because APD began to encroach upon the pacing cycle length, leading to enhanced refractoriness of Na+ current due to incomplete repolarization. 10.1371/journal.pcbi.1002061.g014 Figure 14 AP and Ca2+ alternans at fast pacing. Black lines are control. Gray lines are without CaMK. The two consecutive beats are labeled 1 and 2. A) Pacing at CL = 250 ms. From left to right, top to bottom: AP, expanded time scale showing AP repolarization, Jrel (inset is expanded time scale), [Ca2+]i, [Ca2+]JSR, and Jup. B) Rate dependence of APD (top) and peak [Ca2+]i (bottom) at fast rates (alternans bifurcations disappear without CaMK). C) Same as panel B, but at slower rates (no bifurcations). Since Koller measurements were performed in intact hearts, electrotonic coupling effects would have played a role. Therefore, simulations in a strand of 100 coupled endo cells were conducted to test whether alternans occurred in coupled tissue as well. Indeed, during CL = 280 ms steady state pacing, alternans developed in the multicellular fiber (results shown in Supplement Figure S10 in Text S1). As in Livshitz et al.[77], beat to beat alternans in the Ca2+ subsystem were the cause of the APD alternans in the model. Longer APs coincided with larger Ca2+ transients. For steady state pacing at 250 ms pacing cycle length (shown in Figure 14A), we found that clamping the subspace Ca2+ concentration to either the odd or even beat waveforms eliminated alternans, but clamping of the voltage, myoplasmic Ca2+, ICaL, or INaCa did not eliminate alternans (odd or even beat clamp, not shown). Cutler et al.[78] demonstrated that 30% SERCA upregulation eliminated alternans. Similarly, in our human model, a 20% increase in Jup magnitude eliminated alternans (shown in Supplement Figure S11 in Text S1). CaMK suppression also eliminated alternans in the model (Figure 14A and 14B, gray traces). At slower pacing rates, APD was minimally affected by CaMK suppression. However, the peak Ca2+ concentration was markedly reduced, especially at faster rates (Figure 14C). Currents Participating in Steady State APD Rate Dependence and APD Restitution In order to describe the mechanisms underlying steady state rate dependence and restitution of the APD in the model, it is instructive to first systematically determine which currents participate in these phenomena. In Figure 15, currents were plotted versus Vm during steady state and S1S2 restitution pacing for a variety of CLs and DIs, respectively. If I–V curves are CL or DI independent (i.e. curves overlap), then that current did not participate in steady state rate dependence or restitution, respectively. Conversely, if I–V curves depended greatly on CL or DI, then that current played at least some role in these phenomena. 10.1371/journal.pcbi.1002061.g015 Figure 15 I–V curves during steady state rate dependent pacing at various CLs (panel A), and S1S2 restitution at various DIs (panel B). Arrows indicate decreasing CL or DI. From left to right, top to bottom, results for late INa, Ito, ICaL, IKr, IKs, zoom of plateau ICaL (dashed box section), IK1, INaCa, and INaK are shown. As CL or DI decreased, fast INa, responsible for the maximum AP upstroke velocity and maximum Vm, was reduced (see Figure 9, and principles detailed in Luo and Rudy[79]). This is because shortened time at resting potential between beats prevents complete recovery from inactivation. Thus, at fast pacing rates, and short DIs, maximum Vm and upstroke velocity were reduced, explaining some of what follows. During steady state pacing, IKs was strongly rate dependent (Figure 15A). The I–V curves were dramatically different at different pacing CLs. However, IKs was a relatively small contributor to the rate dependence of APD because IKs density in human ventricle is small under basal conditions (no β-adrenergic stimulation), and changes relative to slow rate values produced minimal additional outward current at fast rates. Late INa, ICaL, INaCa and INaK also showed CL dependent changes during steady state pacing (Figure 15A). INaK became more outward at fast rates. The changes in INaK were dramatic, and the current density was relatively large. Thus, INaK was an important contributor to APD shortening at fast pacing rates. Late INa became dramatically less inward at fast rates, making it a secondary contributor to APD shortening at fast rates. Changes in ICaL and INaCa opposed APD shortening at fast rates; these currents became more inward at short CLs. INaCa increased to match the increased Ca2+ extrusion burden. Importantly, ICaL increased despite reduced channel availability. ICaL inactivation gates recovered less between beats as pacing rate increased (∼20% less at CL = 300 ms compared to CL = 2000 ms). The same mechanism caused reduced late INa at fast rates (availability at CL = 300 ms was ∼1/3 that at CL = 2000 ms). However, influences of increased CaMK facilitation combined with increased driving force (reduced maximum Vm) actually caused ICaL to become larger at fast rates. If Na+ is clamped to small values associated with slow pacing ([Na+]i and [Na+]ss = 6.2 mM at CL = 2000 ms), preventing its accumulation at fast rates, INaK remains small and CL independent (this mechanism is described later in detail), causing plateau voltages to become relatively CL independent. Thus, with Na+ clamp, ICaL changes with pacing rate are different than under control conditions. CL independent plateau voltages confer CL independence to the driving force for plateau ICaL. Na+ clamping reduced Ca2+ (via INaCa) which reduced activated CaMK and thus ICaL facilitation. An interesting consequence is that with Na+ clamp, ICaL changes with CL help to cause APD shortening at fast rates, whereas in control (i.e. no Na+ clamp), ICaL changes with CL oppose APD shortening. During restitution, late INa, Ito, ICaL, IKs and INaCa showed DI dependent changes (Figure 15B). Dramatically less inward late INa at short DIs helped shorten the APD. The mechanism was reduced availability due to residual inactivation at the start of the S2 beat. ICaL was reduced for the same reason. This was evident during the plateau. CaMK facilitation did not depend on DI because Ca2+ accumulation (necessary for CaMK activation) is slow, occurring only after long term pacing to steady state. Similarly, Na+ did not accumulate at short DIs, which kept INaK constant. Therefore, plateau potentials and ICaL driving force during the plateau were relatively DI independent. Just as in the case of Na+ clamp, these properties combined to allow reduced availability of ICaL at short DI to dominate the behavior. However, reduced maximum Vm increased the driving force during the time of peak ICaL, which caused peak current to generally increase at short DIs. At extreme DI of 5 ms, the slow AP upstroke (i.e. reduced dVm/dt) caused mild ICaL inactivation coincident with activation, so the peak current was reduced compared to DI = 10 ms. Changes in other currents (Ito, IKs and INaCa), though nonzero, were relatively minor due to timing. DI dependent changes that increased or reduced current during phase-1 of the AP had little effect on final repolarization time. The exception is IKr. IKr is large enough that early spiking helped shorten APD at very short DIs (detailed simulations follow). Ionic Basis for APD Rate Dependence and Restitution Steady state rate dependence of the APD was largely caused by accumulation of intracellular Na+ at fast rates. This is illustrated in Figure 16A. When [Na+]i and [Na+]ss were clamped to values from steady state pacing at CL = 2000 ms, APD lost much of its sensitivity to pacing rate and remained relatively long. Conversely, when the clamp was to [Na+]i and [Na+]ss from steady state pacing at CL = 300 ms, the APD remained relatively short at all rates. Pacing rate dependent [Na+]i and [Na+]ss changes are linked to the AP via INaK, which responds to [Na+]i levels. INaK increased with [Na+]i at fast rate. However it did not increase, regardless of the pacing rate, when [Na+]i and [Na+]ss were kept low (Na+ at CL = 2000 ms; Figure 16C, right). Moreover, APD remained long at all CLs when INaK was clamped to its slow rate waveform (not shown). 10.1371/journal.pcbi.1002061.g016 Figure 16 Major causes of steady state APD rate dependence and S1S2 APD restitution. A) APD rate dependence in control (solid black), and with [Na+]i and [Na+]ss clamped to slow rate (solid gray) or fast rate (dashed black) values. When late INa (dashed gray) or both late INa and ICaL inactivation gates were reset to their slow rate values (dash-dot-dot gray) in addition to [Na+]i and [Na+]ss slow rate clamp, APD lost almost all rate dependence. Note that slow rate [Na+]i and [Na+]ss clamp alone left residual APD rate dependence, especially at fast rates (CL = 300 to 700 ms, shaded box). B) APD rate dependence (control, solid black) was largely unaffected by resetting inactivation gates for late INa (dashed gray), ICaL (dash-dot-dot gray), or late INa and ICaL (solid gray) to their slow rate values (no [Na+]i and [Na+]ss clamping to slow rate values). C) [Na+]i and INaK increase with pacing rate under control conditions (left). When [Na+]i and [Na+]ss are clamped to slow rate values, INaK is small and rate independent (right). D) APD restitution in control (solid black), and when inactivation gates were reset to S1 values upon S2 delivery (late INa reset – dashed gray, ICaL reset – dash-dot-dot gray, late INa and ICaL reset – solid gray). Shown in dashed black is resetting late INa and ICaL inactivation to the DI = 5 ms value. E) [Na+]i and [Na+]ss clamp to S1 values (dashed gray) did not affect APD restitution (control, solid black). Steady state APD rate dependence was not completely eliminated by Na+ clamp alone. That is, clamping [Na+]i and [Na+]ss to slow rate values did not cause APD curves to become absolutely flat with respect to CL, especially at fast pacing rates (Figure 16A, shaded box CL = 300 to 700 ms, solid gray line). This signifies that other mechanisms are involved. When in addition to clamping [Na+]i and [Na+]ss to their slow rate values, we also reset the inactivation gates for late INa, and especially for both late INa and ICaL to their CL = 2000 ms values at the start of each beat, the APD curve flattened further at fast rates (Figure 16A, dashed gray and dashed-dot-dot gray lines, respectively). Importantly, resetting these inactivation gates alone, without also clamping Na+, had little effect on APD rate dependence (Figure 16B). As described previously, without Na+ clamp, fast pacing caused late INa reduction and ICaL increase; the former helped while the latter opposed APD shortening. However, with Na+ clamp, both currents became less inward with fast pacing. Thus, resetting ICaL inactivation gates to slow rate values had different effects with, versus without Na+ clamping. Na+ clamp prolonged the APD. The prolongation and changed ICaL behavior after Na+ clamp rendered late INa and ICaL gate resetting more potent effectors of further AP prolongation; especially at fast rates where residual inactivation between beats was substantial. Rate dependent Na+ changes only occurred with the steady state pacing protocol due to slow ion accumulation after lengthy pacing regimes. For APD restitution, clamping [Na+]i and [Na+]ss to values from S1 pacing during the S2 beat did not affect APD (Figure 16E). However, restitution was dramatically affected by resetting inactivation gates for late INa and/or ICaL to their S1 starting values at the start of the S2 beat (Figure 16D). APD remained long for all DIs. Conversely, when late INa and/or ICaL inactivation gates were reset to S2 starting values for DI = 5 ms, APD remained short for all DIs. Again, resetting these inactivation gates to their slow rate values had only minor effects on steady state APD rate dependence (Figure 16B). At very short DIs, IKr played an important role in APD restitution. In Figure 17A, the fast and slow time dependent deactivation gates (xrfast and xrslow, respectively) were reset to their value at DI = S1 = 1000 ms (dashed gray line, compare to control solid black line). Deactivation of IKr is slow (Figure 3B). For DI = S1, deactivation was complete between beats. At short DIs, it was incomplete at the start of the S2 beat, enhancing IKr availability (early IKr spiking, Figure 17B, bottom) and outward current that contributes to APD shortening. The enhanced availability only mattered at very short DIs, because at these DIs APD was short enough that increased outward current during phase-1 of the AP affected final repolarization time. Changes to currents during later AP phases 2 and 3 (during the plateau and early repolarization, e.g. late INa and ICaL), generally have greater impact on the APD. Early IKr spiking reduced maximum Vm, which affected all other currents, including late INa and ICaL. Comparison with Other Human Ventricular AP Models Several important differences exist between the ORd model presented here and other human models (e.g. TP[9] and GB[10] models). Notably, model differences in the rate of repolarization and EAD formation were examined in direct comparison with experiments (Figures 7C, and Figure 11A, respectively). Readers wishing to simulate the human ventricular AP have a choice of models. To help further differentiate the models, additional comparisons are shown in Figure 18. 10.1371/journal.pcbi.1002061.g017 Figure 17 IKr deactivation is important for APD restitution at very short DIs. A) APD restitution in control (solid black), and when the fast and slow deactivation gates (xrfast and xrslow) were reset to the DI = S1 = 1000 ms value at the start of the S2 beat (dashed gray). Bottom) Zoom in to more clearly show the consequence of deactivation resetting at short DIs (section outlined by dashed box above). B) Traces for the AP (top) and IKr (bottom) during the S2 beat at different DIs (indicated by arrows). Spiking in IKr occurred early during the AP at short DI. Spiking was caused by slow deactivation, increasing availability of IKr. 10.1371/journal.pcbi.1002061.g018 Figure 18 Comparison with other human ventricular AP models. Single endo cell simulations from ORd, TP, and GB models are solid black, gray, and dashed black lines, respectively. Experimental results (small tissue preparations) are white squares. A) APD rate dependence. Results for APD90, 70, 50 and 30 are shown top left, top right, bottom left, and bottom right, respectively. B) The AP, major currents, [Na+]i, and [Ca2+]i at steady state for CL = 1000 ms. From left to right, top to bottom: AP (with Vm rest inset at far right), INa (inset shows peaks), late INa (not present in TP or GB models), Ito (inset shows decay rate), ICaL (arrow shows ORd peak magnitude; inset shows normalized peaks, which are wide in TP and GB), IKr (arrow shows ORd early spike peak magnitude), IKs, IK1, INaCa, INaK, [Ca2+]i, [Ca2+]i decay rate, and [Na+]i. Undiseased human ventricular measurements of steady state rate dependence of APD90, 70, 50 and 30 were accurately reproduced by the ORd model (Figure 18A, same data as in Figure 7A). Rate dependence of APD90 is fairly accurate in the TP model. However, rate dependence of APD70, 50 and especially APD30 are not accurate. The GB model repolarization rate is more accurate, but divergence from the measurements is large for APD30. At fast pacing rates, GB model APD90 is accurate. Slow pacing APD90 is long compared with experiments (at CL = 2000 ms, APD90 is ∼40 ms longer than in experiments). In addition, APD rate dependence does not plateau at CL = 2000 ms. In Figure 18B, the AP, major currents, and [Na+]i and [Ca2+]i were compared between models. Simulations were in single endo cells paced to steady state at CL = 1000 ms. Of note, the TP and GB models do not include late INa. The width of the ICaL peak and the morphology were model dependent. It was “cigar shaped” in the TP model. In the GB model, the ICaL peak was broad and poorly defined. The ORd model ICaL peak was sharp, as seen in undiseased human ventricle experiments (AP clamp, Figure 1D). IKr was relatively small in the GB model, but shared a similar morphology with the ORd model. The TP IKr morphology is characterized by an early spike and a wider late spike. The IKs density in the TP model was much larger than in the other models (∼10 fold larger). Density and morphology of INaCa was model dependent. INaCa was smallest in the ORd model (based on undiseased human measurements, Figure 2B). INaK was roughly 1.5 and 2 fold greater in GB and TP models, respectively, compared with ORd. The Ca2+ transient peak was much larger in the TP model than in the other models, which were similar to each other. The decay rate of [Ca2+]i was somewhat slower in the ORd model (accurate to undiseased human measurements; Figure 12 panels C and D). Model [Na+]i was 7.2, 8.2, and 9.7 mM in ORd, GB, and TP models, respectively. Discussion Though the available undiseased human ventricle dataset has been missing essential elements, several human ventricle AP models have been developed and published. The Priebe and Beuckelmann model[35] lacks human specific data for reformulation of major currents, and so was based in large part on its guinea pig predecessor[80]. The TP model[81] and updated version[9] is easy to use, includes many reformulated currents, and simulates physiological restitution and alternans. However, both the TP and GB[10] models lack sufficient ICaL data for validation, and cannot produce EADs. The GB model includes K+ current reformulations using undiseased human data for validation, but does not produce AP or Ca2+ transient alternans. EADs and alternans are both important mechanisms of arrhythmogenesis and should be reproduced in simulation studies of human arrhythmias. The Iyer et al. model[34] is based almost entirely on data from human channels expressed in non myocytes. Though the expressed channels are human, native myocyte ion channels in the ventricle are composed of a variety of protein isoform combinations, auxiliary subunits, cytoskeletal elements, and membrane lipid composition, all of which may influence channel behavior. Anchoring and other regulatory proteins present in native cells also define the local environment for ICaL in particular[82], but are not present in expression systems. Fink et al. modified the TP model[36] to include updated IKr and IK1 (with [K+]o dependence) formulations, based on undiseased human ventricular measurements. The rate of AP repolarization in this modified scheme is more accurate compared with the original TP model. For these advantages, the model sacrifices runtime speed (Markov formulation is used for IKr). Other core issues of the TP model carry over to this modified version (incorrect ICaL, non-physiologically large IKs, and no EAD generation under appropriate conditions). We believe that the new undiseased human data presented here are essential, and substantially increase human specific model accuracy. Due to extensive validation using these new data, our model reproduces all of the following important physiological behaviors: 1) CDI versus VDI inactivation of ICaL; 2) reformulated, detailed and accurate kinetics (using weighted time constants) for Ito, INaCa, IK1, IKr, IKs, fast INa, and late INa; 3) AP repolarization rate from 30% to 90% repolarization; 4) APD at all physiological pacing rates with/without block of major currents, 5) APD restitution with/without block of delayed rectifier currents; 6) transmural heterogeneity causing upright pseudo-ECG T-wave; 7) early afterdepolarizations (EADs); 8) effects of CaMK; and 9) AP and Ca2+ transient alternans. EADs and Repolarization Rate One of the most important aspects of the model is its close correspondence to experimental measurements of not only APD90, but also to APD30, 50 and 70 at all physiologically relevant pacing rates and for S1S2 restitution. This large pool of data has previously been unavailable. Accurate repolarization rate (i.e. time between APD30 and 90) for the restitution protocol is crucial for simulating any phenomenon related to reentrant arrhythmia, where head-tail interactions determine refractoriness and vulnerability[83]. Use of new undiseased data for currents that are active during the plateau and phase-3 of the AP (ICaL, INaCa, IKr and IKs) contributed to the correct repolarization rate. The rate of repolarization and its effects on ICaL control EAD formation in this model, as in canonical EAD explanations [75], [84]. Failure of the TP and GB models to reproduce EADs may be due in part to their accelerated repolarization rates (Figure 7C). It may also be caused by inaccurate formulation of ICaL inactivation, developed in absence of the essential undiseased human data presented here. Steady State APD Rate Dependence Due to the small amplitude and rapid deactivation kinetics of IKs in the human ventricle in absence of β-adrenergic stimulation, it does not play a major role in determining APD, APD rate dependence, or APD restitution under basal conditions[85] (Figure 8). This is in contrast to guinea pig ventricle, where slower deactivation and larger amplitude IKs make it the most important current for steady state APD rate dependence (simulations[86] and experiments[87]). Phosphorylation by PKA in the case of β-adrenergic stimulation greatly enhances both the activation rate and amplitude of IKs [88]. With β-adrenergic stimulation, IKs plays an important role in steady state APD rate dependence[89]. Clearly, IKs is important under various circumstances – the AP repolarizes in human ventricle experiments even when IKr is blocked[85], and clinical long QT syndrome type 1 is caused by IKs loss of function[90]. Typically, isolated myocyte patch clamp experiments underestimate IKs due to enzymatic degradation[42]. In ORd, the role of IKs was validated using small tissue preparations, where selective IKs block prolonged APD, but only very modestly under basal conditions (no β-adrenergic stimulation, <15 milliseconds in experiments and simulations at CL = 1000 ms, Figure 8). Block of IKr caused the most severe changes to the human AP (rate dependence and restitution, Figure 8). However, Supplement Figure S5 in Text S1, and Figure 15A show that IKr is rate independent, as in experiments[39] and therefore was not responsible for causing APD changes with pacing rate. Rather, our simulations identified rate dependent changes in INaK secondary to [Na+]i accumulation as a primary cause of APD rate dependence (Figure 16A, 16C). This finding is not new. Simulations in dog ventricle[19], human atrium[91], and in the GB human ventricle[10] models all led to this conclusion. However, findings from the Iyer human model[34] differ, at least in part, regarding this mechanism. In the Iyer model, [Na+]i affected APD rate dependence via INaCa, which is primarily outward at fast rates. Rate dependence in the TP model[9] is less [Na+]i dependent because, as Grandi discussed[10], IKs is exaggerated. Experiments by Pieske et al.[57] investigated [Na+]i in heart failure, versus nonfailing human ventricular myocytes. Pieske experiments demonstrate that rate dependent [Na+]i accumulation is an important phenomenon in health and disease. However, additional experiments are needed to determine whether and how [Na+]i affects INaK and APD in human ventricle. In addition to INaK and INaCa (both included in the ORd model), intracellular Na+ is also mediated by fluxes related to H+, CO2, and HCO3 - homeostasis. Exchangers and cotransporters move Na+ ions down the electrochemical gradient in order to offset the cost of H+, CO2, and HCO3 - pumping. Na+ rate dependent handling and consequently INaK should be affected by these processes, which were not explicitly included in the ORd model. In the absence of H+, CO2, and HCO3 - fluxes, it is possible that the role of INaK might have been over estimated. It is important to address this because INaK and its response to Na+ accumulation was a major cause of APD rate dependence in the model. Thus, we performed simulations where H+, CO2, and HCO3 - effects on Na+ were explicitly included, using Crampin and Smith equations[92] (Supplement Figure S12 in Text S1). Quantitative details of Na+ handling, INaK and APD rate dependence were affected when we included H+, CO2, and HCO3 - handling processes. However, the qualitative outcomes were not affected. INaK increase with fast pacing, secondary to Na+ accumulation, was still the primary determinant of APD rate dependence during steady state pacing. Removal of the effects of Na+ accumulation on steady state APD rate dependence by clamping [Na+]i and [Na+]ss did not completely eliminate APD rate dependence. Especially at fast rates (Figure 16A, shaded box CL = 300 to 700 ms, solid gray line), APD was not absolutely flat with respect to CL. APD rate dependence was largely unaffected by resetting inactivation gates for late INa, and/or ICaL to their slow rate values at the start of each beat (Figure 16B). Interestingly, if these gates were reset while also clamping Na+ to slow rate values, the APD-CL curve became almost completely flat, even at fast rates (Figure 16A, dashed gray and dashed-dot-dot gray lines, respectively). Thus, accumulation of Na+ and consequent effects on INaK is a major cause of APD rate dependence, however, not the only cause. Other currents also participate at fast pacing rates. Though the GB model[10] demonstrated the Na+/INaK/APD rate dependence mechanism, it did not show the additional effects of late INa and ICaL. The GB model cannot show these multi-factorial causes of APD rate dependence because it does not include late Na+ current (Figure 18), and because ICaL kinetics are inaccurate due to lack of experimental data. Due to charge conservation, accumulation of [Na+]i is associated with an equal reduction in [K+]i and a volume converted [K+]o increase in tissue clefts and interstitial spaces[93]. This can affect behavior by increasing IK1 (its [K+]o sensitivity is included in this model), which depolarizes resting voltage and reduces excitability. However, our simulations represent experiments in an isolated myocyte in a large bath, where [K+]o is constant. Even in vivo, [K+]o is tightly controlled via regulation by the lymphatic system and kidneys. APD Restitution We showed that in contrast to steady state rate dependence, [Na+]i had no effect on APD restitution. Rather, restitution was primarily caused by the time course of recovery from inactivation of late INa and ICaL; processes which had little effect on steady state-rate dependence of APD (absent Na+ clamp). At very short DIs, slow deactivation of IKr caused increased availability and spiking, which helped shorten the APD. APD rate dependence was caused primarily by concentration changes, while APD restitution was caused by gating kinetics. Previous studies have not made this important distinction between steady state rate dependence and restitution mechanisms in human. The role of ICaL and its inactivation kinetics in APD restitution reiterates the primacy of ICaL in determining basic physiological behaviors, highlighting the importance of the new ICaL experimental data, presented here, to model development and validation. A role for late INa in restitution could not have been hypothesized using TP or GB models, which have no late INa. The density of late INa was constrained in the ORd model by experiments from nonfailing human ventricular myocyte measurements by Maltsev et al.[51], where the late current was measured 200 ms after the start of the voltage clamp step (∼0.35 µA/µF I–V curve maximum). The maximum late INa during the free running AP model was much smaller (∼0.15 µA/µF, about half the I–V curve maximum) even at slow pacing rates, where late INa was largest. Late current is difficult to measure directly, and it is possible that the current density was overestimated due to selection bias. That is, late INa is small, and not all cells produced measureable late current (2 of 3 myocytes[51]). However, we consider the model density of late INa to be accurate based on model reproduction of experiments which consistently showed substantial APD90 shortening following application of 10 µM mexiletine in undiseased human myocardium (90% late INa block in simulations, Figure 8A). Ca2+ Cycling, CaMK and Alternans Previously published human ventricle AP models have not incorporated the CaMK signaling pathway. Our human simulations show, as in dog simulations[31], [77], that CaMK plays an important role in determining frequency dependence of Ca2+ cycling (Figure 13). The model also shows that the integrated electrophysiological consequence of CaMK effects on target channels is minimal. That is, CaMK suppression had only minor effects on APD rate dependence and AP morphology. At very fast pacing (CLs <300 ms), the Ca2+ cycling consequences of CaMK phosphorylation were central to alternans formation. Suppression of CaMK eliminated alternans. CaMK related findings are in agreement with simulations using other models developed by our group[77], models from other groups[94], and experiments[95]. However, experiments showing the effects of pharmacological suppression of CaMK on rate dependent behaviors (e.g. by Wehrens et al.[96] with KN-93 in rabbit) should be performed in human ventricular myocytes to validate model predictions. Transmural Heterogeneity The method used for implementation of the transmural cell types (M and epi cell), based on the thoroughly validated endo cell framework, was simplistic. That is, we considered that channel conductance was proportional to transmural gradients in mRNA or protein expression for alpha subunits of ion channels. Only in the case of Ito were functional current measurement data available[70]. Staying within error bars for mRNA or protein data[67], [68], [69], channel conductances were modulated so that the simulated transmural AP differences were consistent with experiments[50], [71]. The effect of transmural heterogeneity of accessory β-subunits was not considered in the transmural cell type definitions. However, in the case of IKs, the KCNE1 β-subunit is transmurally heterogeneous. KCNE1 protein was about two times greater in M-cells compared to epi cells[69]. The presence of KCNE1 carries two important functional consequences 1) ∼5 fold slower activation and 2) ∼5 fold larger conductance[97]. Therefore, theoretically, twice as much KCNE1 in M-cells may increase the variable stoichiometry ratio of KCNE1 to alpha subunit KCNQ1[98], slowing activation and increasing conductance. We conducted simulations to evaluate the influence of KCNE1 heterogeneity on IKs and the AP (Supplement Figure S13 and related description in Text S1). Due to the small amplitude of human IKs in the absence of β-adrenergic stimulation, implementation of KCNE1 heterogeneity did not appreciably affect the AP (Supplement Figures S13 and S19 in Text S1, where transmural APDs are shown to be relatively insensitive to changes in IKs conductance). Interestingly, the simulated effects of KCNE1 on activation rate and conductance offset one another. That is, slowed activation and larger conductance in M-cells yielded IKs current that was remarkably close to the control case. Similar results were found for epi cell simulations: the effects of faster activation and reduced conductance were offsetting such that their combined effect was minimal. APD Accommodation Steady state rate dependence of APD and APD restitution were the focus of the simulations and experiments in this study. However, the time course of APD response to abrupt changes in pacing rate has been shown in human by Franz et al.[99], and simulated in the TP model by Pueyo et al.[100] as a marker for arrhythmia risk. Simulations of APD accommodation in our model compare favorably to Franz experiments (same pacing protocols used in experiments were used in the simulations, Supplement Figure S14 in Text S1). Single exponential curves were fit to the time dependence of APD changes. For abrupt CL reduction from 750 to 480 ms, the time constant was 250 and 284 seconds in experiments and simulations, respectively. Time constants were 300 and 299 seconds in experiments and simulations, respectively, when CL was abruptly returned to 750 ms. When the CL reduction was more severe, from CL = 750 to 410 ms, the time constants were 252 and 165 seconds, in experiments and simulations, respectively. For return to CL = 750 ms, the time constants were 350 and 289 seconds, respectively. Pueyo used time to 90% accommodation to compare model with experiments demonstrating similar accuracy. Both simulation studies also show initial overshoot, or “notching”, as observed and described by Franz. Parameter Sensitivity As in Romero et al.[101], we performed a sensitivity analysis to determine factors participating in important model outputs, including 1) steady state APD90 rate dependence (Supplement Figure S15 in Text S1), 2) S1S2 APD90 restitution (Supplement Figure S16 in Text S1), 3) rate dependence of maximum (systolic) [Ca2+]i (Supplement Figure S17 in Text S1), 4) rate dependence of [Na+]i (Supplement Figure S18 in Text S1), and 5) transmural cell type APD90 at steady state (Supplement Figure S19 in Text S1). The findings from our analysis were similar to those shown by Romero et al. using the TP human AP model[101]. That is, in ORd and TP models, IKr and ICaL affect APD90 while ICaL, INaCa, and INaK affect peak [Ca2+]i. One important difference is the role for IKs. A much larger role was played by IKs in the TP model (∼10 fold larger density than in other human models, Figure 18B). In the TP model, IKs is responsible for steady state rate dependence of the APD (shown by Grandi et al.[10]). IKr conductance changes affect APD90 substantially in our model. This was expected, since IKr is the largest outward current (also in experiments, Figure 8, and in Romero's analysis using the TP model). Though IKr affects APD, it is not responsible for its rate dependence. Conductance changes in INaK did not substantially affect APD90 because INaK is a relatively small current. Yet, rate dependent changes in INaK (secondary to Na+ accumulation at fast rate) were the primary determinant of APD rate dependence. [Na+]i at different pacing rates, and thus its relative changes with rate, was by far most sensitive to INaK conductance (Supplement Figure S18 in Text S1). This supports our strategy for setting INaK conductance to reproduce rate dependence of [Na+]i in nonfailing human myocytes[57]. Computational Tractability and Model Stability To keep the ORd model computationally efficient and parameters well constrained, the Hodgkin-Huxley formalism was used in formulating current equations. This choice was made as a design principal with the thought that interested users can modularly replace any current or flux with more detailed Markov formulations of mutation or drug effects as desired (e.g.[53], [66]). Similarly, intracellular Ca2+ handling can be modified (e.g. more spatial detail, Markov ryanodine receptor implementation), or various signaling pathways and related effects on ion channels can be added (e.g.[31], [89], [102]). The basic ORd model has 41 state variables. In the absence of CaMK and its effects on target currents and fluxes, the number of state variables is 31. Exclusion of Markov models increases parameter constraint. It also prevents the system of differential-algebraic equations from being overly stiff. This promotes model stability and computational tractability. Using the rapid integration technique described in Supplement Text S1 (Computational Methodology section), the model arrives at true and accurate steady state in under one minute of runtime (∼1000 beats are needed, depending on the CL, Visual C++ running on a desktop PC; more details in Supplement Text S1). ORd equations are all smoothly varying functions, free of singularities and “if” conditionals. Thus, the model can readily be implemented in any of a variety of automated numerical integrators, such as Matlab (The MathWorks Inc.), CellML (http://www.cellml.org/), CHASTE[103], or CARP (CardioSolv LLC.). Limitations Direct measurement of INaK in the undiseased or nonfailing human ventricular myocyte is lacking. Therefore, INaK was validated by reproduction of important biophysical properties (Supplement Figure S7 in Text S1), and by reproduction of [Na+]i rate dependence measured in nonfailing human ventricular myocytes (Pieske et al.[57], Figure 12A). However, independent and direct experimental measurement of INaK in undiseased or nonfailing human ventricular myocytes would provide additional support for the mechanistic conclusion that INaK changes secondary to Na+ accumulation at fast pacing rates is a major determinant of steady state APD rate dependence. This conclusion is consistent with several other modeling studies which proposed the same mechanism (dog ventricle[19], human atrium[91], and human ventricle[10]). The relationship between [Na+]i, INaK and steady state APD rate dependence was robust. It was not disrupted by including the effects of Na+/H+ and Na+/HCO3 - exchange fluxes on Na+ handling (Crampin and Smith equations[92], Supplement Figure S12 in Text S1). Na+ accumulation and INaK response were not the only cause of APD rate dependence in the ORd model. At fast pacing rates (CL = 300 to 700 ms), late INa and ICaL were also involved (Figure 16A, and related discussion). Measurements of undiseased human endo APs were performed in small tissue preparations (1–3 gram pieces). This was to avoid possible enzymatic degradation of K+ channel proteins[40], [42], affecting currents and the AP. However, electrical loading in tissue subtly affects behavior[19]. We performed simulations using a multicellular fiber model to include loading effects, which had only minor consequences (Figure S8). APD was ∼275 ms in our human endo preparations at CL = 1000 ms, well matched by the model (273 ms). In vivo noninvasive electrocardiographic imaging of the activation-recovery interval, an indicator of the cellular epi APD, was ∼260 ms in healthy adults[104]. Human monophasic AP measurements are also in this range[58]. Measurements from Drouin et al. showed longer APDs (∼350 ms in endo cells on the cut transmural face at CL = 1000 ms). Having validated the endo model based on more than 100 of our own endo AP measurements, we thought it reasonable to use Drouin transmural APD ratios, rather than the uniformly longer APDs themselves, for validation of the transmural cell type models. The presence of M cell APs in the nonfailing human heart was observed by Drouin et al.[50], and more recently by Glukhov et al.[71]. However, there is controversy regarding the M cell definition and its role in human. Our M cell model was based on data where the M cell was defined by its transmural location. The resulting simulated M cell AP corresponded with the “max cell” observed by Glukhov. Recently, Sarkar and Sobie developed a method for quantitative analysis of parameter constraint and relationships between parameters and target outputs in AP models[105]. We did not apply this analysis during model development. However, the extensive validation of channel kinetics and the emergent response of the AP to a variety of dynamic pacing protocols, used in development and validation of the model, ensures sufficient parameter constraint. The parameter sensitivity tests we performed were instructive, though relatively limited (conductance changes only). Application of Sarkar and Sobie's analysis to our model is beyond the scope of this paper, but should provide worthwhile insights regarding inter-relatedness of processes in the human ventricle, in addition to formally testing parameter constraint. Materials and Methods Characteristics of Human Tissue During the last 15 years, undiseased hearts were donated for research in compliance with the Declaration of Helsinki and were approved by the Scientific and Research Ethical Committee of the Medical Scientific Board of the Hungarian Ministry of Health (ETT-TUKEB), under ethical approval No 4991-0/2010-1018EKU (339/PI/010). Data from 140 hearts were used in this study. Of these, 78 were from male donors (56%). The average donor age was 41 years old with standard deviation of 12 years. Tissue Preparation Tissue transport and ventricular endo preparations were performed as previously described[85]. Tissue was carefully pinned and placed in a modified Tyrodes superfusate (in mM: NaCl 115, KCl 4, CaCl2 1.8, MgCl2 1, NaHCO3 20, and glucose 11, pH 7.35, 37°C), and point stimulation was via bipolar platinum electrodes. Drug solutions were made fresh on the day of use. Drugs included in this study were, in µM: E-4031 1, HMR-1556 1, nisoldipine 1, BaCl2 100, ryanodine 5, mexiletine 10. Simulated application of these drugs was 70% IKr [59], and 90% IKs [60], ICaL [61], IK1 [62], RyR[63], and late INa [64] block, respectively. Myocyte Isolation Tissue transport and myocyte isolation for the undiseased donor hearts were as previously described[85]. Myocyte isolation commenced immediately upon arrival in the laboratory, using the perfusion disaggregation procedure, previously described[85]. Electrophysiology Data were obtained using conventional whole cell patch-clamp techniques. Micropipette fabrication and data acquisition were as described previously for undiseased donor heart[85]. Axopatch 200 amplifiers, Digidata 1200 converters, and pClamp software were used (Axon Instruments/Molecular Devices). Experiments were performed at 37°C. The standard bath solution contained, in mM: NaCl 144, NaH2PO4 0.33, KCl 4.0, CaCl2 1.8, MgCl2 0.53, Glucose 5.5, and HEPES 5.0 at pH of 7.4, and pipette solutions contained K-aspartate 100, KCl 25, K2ATP 5, MgCl2 1, EGTA 10 and HEPES 5. The pH was adjusted to 7.2 by KOH (+15−20 mM K+). For L-type Ca2+ current measurement, the bath solution contained in mM: tetraethylammonium chloride (TEA-Cl) 157, MgCl2 0.5, HEPES 10, and 1 mM CaCl2, or BaCl2, or SrCl2 (pH 7.4 with CsOH). The pipette solution contained (in mM) CsCl 125, TEA-Cl 20, MgATP 5, creatine phosphate 3.6, EGTA 10, and HEPES 10 (pH 7.2 with CsOH). For Na+/Ca2+ exchange current measurement, the bath solution contained, (in mM): NaCl 135, CsCl 10, CaCl2 1, MgCl2 1, BaCl2 0.2, NaH2PO4 0.33, TEACl 10, HEPES 10, glucose 10 and (in µM) ouabain 20, nisoldipine 1, lidocaine 50, pH 7.4. The pipette solution contained (in mM): CsOH 140, aspartic acid 75, TEACl 20, MgATP 5, HEPES 10, NaCl 20, EGTA 20, CaCl2 10 (pH 7.2 with CsOH). Ca2+ Transient Florescence Isolated myocytes from the undiseased donor hearts were used to measure the Ca2+ transient during point stimulation via bipolar platinum electrodes, indicated by Fura-2-AM, as was described previously[106]. Bath temperature was 37°C. Pacing Protocols For both experiments and simulations, we determined APD at 30, 50, 70 and 90% of complete repolarization (APD30–90, in ms). The start of the AP was the time of maximum dVm/dt. The time of APDX occurred once membrane voltage was X% of the resting value. Resting voltage was measured immediately prior to each paced beat. For APD rate dependence, pacing was to steady state. For APD restitution (S1S2, or static restitution), S1 pacing was at cycle length (CL) = 1000 ms. The S2 beat was delivered at variable diastolic intervals (DIs), measured relative to APD90. The dynamic restitution protocol was simulated as in experiments by Koller et al.[58]. Pacing was at a variety of rates (30 seconds at CLs from 230 to 1000 ms, no S2 beats). APD95 was plotted against DI (where DI = CL – APD95). Unlike static S1S2 restitution, the dynamic restitution protocol allows for more than one APD to be associated with a given DI. This is significant because bifurcation in the dynamic restitution curve is believed to be arrhythmogenic[107]. Population Based CaMK effects For all channels affected by CaMK, we created separate models for the fully CaMK phosphorylated channels, and the totally non phosphorylated channels. Then, based on the degree of CaMK activation (CaMKactive), we determined the proportion of channels affected by CaMK. To calculate the CaMK affected current (or flux), we added the weighted sum of fully affected and totally unaffected channels, using the proportionality. The model employed for CaMK activity was validated previously[31], [77]. Relative Weights in a Two Time Constant Scheme When measurements called for a gating process to be represented by both a fast and a slow process, we included separate fast and slow gates. However, we did not simply multiply fast and slow gates to modulate conductance as others have done previously. To do so allows the fast process alone to completely control deactivation/inactivation, and the slow process alone to completely control activation/recovery. Rather, since measurements of bi-exponential behaviors provide the relative weight of fast/slow processes, we modeled the measurements accordingly, and used the weighted sum of fast and slow processes. Transmural Wedge Simulation We computed the pseudo-ECG using a 1-dimensional model of the transmural wedge preparation[108], [109]. In brief, the spatially weighted sum of the voltage gradient was determined at a point 2 cm from the epi end of a heterogeneous multicellular fiber, along the fiber axis. Cells 1–60 were endo, 61–105 were M, and 106–165 were epi. The stimulus was delivered to cell 1. Cells 15 from both ends of the fiber were excluded from the gradient measurement due to confounding edge effects. Pacing was for 100 beats using steady state initial conditions from paced single cells. Equations, Computers, and Software All model equations, hardware and software used, and rapid integration methods are provided in Supplement Text S1. Model code can be downloaded from the research section of our website: http://rudylab.wustl.edu. Supporting Information Text S1 Supplementary materials. (PDF) Click here for additional data file.
              Bookmark
              • Record: found
              • Abstract: found
              • Article: not found

              A model for human ventricular tissue.

              The experimental and clinical possibilities for studying cardiac arrhythmias in human ventricular myocardium are very limited. Therefore, the use of alternative methods such as computer simulations is of great importance. In this article we introduce a mathematical model of the action potential of human ventricular cells that, while including a high level of electrophysiological detail, is computationally cost-effective enough to be applied in large-scale spatial simulations for the study of reentrant arrhythmias. The model is based on recent experimental data on most of the major ionic currents: the fast sodium, L-type calcium, transient outward, rapid and slow delayed rectifier, and inward rectifier currents. The model includes a basic calcium dynamics, allowing for the realistic modeling of calcium transients, calcium current inactivation, and the contraction staircase. We are able to reproduce human epicardial, endocardial, and M cell action potentials and show that differences can be explained by differences in the transient outward and slow delayed rectifier currents. Our model reproduces the experimentally observed data on action potential duration restitution, which is an important characteristic for reentrant arrhythmias. The conduction velocity restitution of our model is broader than in other models and agrees better with available data. Finally, we model the dynamics of spiral wave rotation in a two-dimensional sheet of human ventricular tissue and show that the spiral wave follows a complex meandering pattern and has a period of 265 ms. We conclude that the proposed model reproduces a variety of electrophysiological behaviors and provides a basis for studies of reentrant arrhythmias in human ventricular tissue.
                Bookmark

                Author and article information

                Contributors
                Role: Formal analysisRole: InvestigationRole: MethodologyRole: Writing – review & editing
                Role: ConceptualizationRole: Data curationRole: Funding acquisitionRole: InvestigationRole: Writing – original draftRole: Writing – review & editing
                Role: ConceptualizationRole: Data curationRole: Formal analysisRole: Funding acquisitionRole: InvestigationRole: SupervisionRole: Writing – original draftRole: Writing – review & editing
                Role: Editor
                Journal
                PLoS Comput Biol
                PLoS Comput Biol
                plos
                PLOS Computational Biology
                Public Library of Science (San Francisco, CA USA )
                1553-734X
                1553-7358
                February 2024
                28 February 2024
                : 20
                : 2
                : e1011930
                Affiliations
                [1 ] Department of Physics, South China University of Technology, Guangzhou, China
                [2 ] Department of Medicine, David Geffen School of Medicine, University of California, Los Angeles, Los Angeles, California, United States of America
                University of Virginia, UNITED STATES
                Author notes

                The authors have declared that no competing interests exist.

                Author information
                https://orcid.org/0000-0001-5217-6022
                https://orcid.org/0000-0003-1524-5919
                Article
                PCOMPBIOL-D-23-01792
                10.1371/journal.pcbi.1011930
                10927084
                38416778
                17f6b672-2f1e-425c-9141-10112b9ddda5
                © 2024 Wang et al

                This is an open access article distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited.

                History
                : 4 November 2023
                : 19 February 2024
                Page count
                Figures: 6, Tables: 1, Pages: 17
                Funding
                Funded by: Guangdong Basic and Applied Basic Research Foundation
                Award ID: 2021A1515010500
                Award Recipient :
                Funded by: funder-id http://dx.doi.org/10.13039/100000002, National Institutes of Health;
                Award ID: R01 HL134709
                Award Recipient :
                Funded by: funder-id http://dx.doi.org/10.13039/100000002, National Institutes of Health;
                Award ID: R01 HL139829
                Award Recipient :
                Funded by: funder-id http://dx.doi.org/10.13039/100000002, National Institutes of Health;
                Award ID: R01 HL157116
                Award Recipient :
                Funded by: funder-id http://dx.doi.org/10.13039/100000002, National Institutes of Health;
                Award ID: P01 HL164311
                Award Recipient :
                This work is supported by the Guangdong Basic and Applied Basic Research Foundation under Grant No. 2021A1515010500 (XH); National Institutes of Health grants R01 HL134709, R01 HL139829, R01 HL157116, and P01 HL164311 (ZQ). The funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript.
                Categories
                Research Article
                Computer and Information Sciences
                Data Management
                Data Visualization
                Phase Diagrams
                Computer and Information Sciences
                Computer Applications
                Computer Modeling
                Biology and Life Sciences
                Biochemistry
                Bioenergetics
                Ionic Current
                Biology and Life Sciences
                Physiology
                Electrophysiology
                Membrane Potential
                Depolarization
                Research and Analysis Methods
                Bioassays and Physiological Analysis
                Electrophysiological Techniques
                Membrane Electrophysiology
                Patch Clamp Techniques
                Computer and Information Sciences
                Systems Science
                Nonlinear Dynamics
                Physical Sciences
                Mathematics
                Systems Science
                Nonlinear Dynamics
                Biology and Life Sciences
                Cell Biology
                Cellular Types
                Animal Cells
                Muscle Cells
                Biology and Life Sciences
                Anatomy
                Biological Tissue
                Muscle Tissue
                Muscle Cells
                Medicine and Health Sciences
                Anatomy
                Biological Tissue
                Muscle Tissue
                Muscle Cells
                Biology and Life Sciences
                Physiology
                Electrophysiology
                Membrane Potential
                Action Potentials
                Biology and Life Sciences
                Physiology
                Electrophysiology
                Neurophysiology
                Action Potentials
                Biology and Life Sciences
                Neuroscience
                Neurophysiology
                Action Potentials
                Custom metadata
                vor-update-to-uncorrected-proof
                2024-03-11
                All relevant data are within the manuscript and its Supporting information files.

                Quantitative & Systems biology
                Quantitative & Systems biology

                Comments

                Comment on this article