Simulation of the undiseased human cardiac ventricular action potential: model formulation and experimental validation

Abstract

Cellular electrophysiology experiments, important for understanding cardiac arrhythmia mechanisms, are usually performed with channels expressed in non myocytes, or with non-human myocytes. Differences between cell types and species affect results. Thus, an accurate model for the undiseased human ventricular action potential (AP) which reproduces a broad range of physiological behaviors is needed. Such a model requires extensive experimental data, but essential elements have been unavailable. Here, we develop a human ventricular AP model using new undiseased human ventricular data: Ca(2+) versus voltage dependent inactivation of L-type Ca(2+) current (I(CaL)); kinetics for the transient outward, rapid delayed rectifier (I(Kr)), Na(+)/Ca(2+) exchange (I(NaCa)), and inward rectifier currents; AP recordings at all physiological cycle lengths; and rate dependence and restitution of AP duration (APD) with and without a variety of specific channel blockers. Simulated APs reproduced the experimental AP morphology, APD rate dependence, and restitution. Using undiseased human mRNA and protein data, models for different transmural cell types were developed. Experiments for rate dependence of Ca(2+) (including peak and decay) and intracellular sodium ([Na(+)](i)) in undiseased human myocytes were quantitatively reproduced by the model. Early afterdepolarizations were induced by I(Kr) block during slow pacing, and AP and Ca(2+) alternans appeared at rates >200 bpm, as observed in the nonfailing human ventricle. Ca(2+)/calmodulin-dependent protein kinase II (CaMK) modulated rate dependence of Ca(2+) cycling. I(NaCa) linked Ca(2+) alternation to AP alternans. CaMK suppression or SERCA upregulation eliminated alternans. Steady state APD rate dependence was caused primarily by changes in [Na(+)](i), via its modulation of the electrogenic Na(+)/K(+) ATPase current. At fast pacing rates, late Na(+) current and I(CaL) were also contributors. APD shortening during restitution was primarily dependent on reduced late Na(+) and I(CaL) currents due to inactivation at short diastolic intervals, with additional contribution from elevated I(Kr) due to incomplete deactivation.

Figure 1. Undiseased human I_CaL experiments and model validation.

A) Experiments: I_CaL traces for currents carried by Ca²⁺ (top), Ba²⁺ (middle), and Sr²⁺ (bottom). The voltage protocol is below the Ca²⁺ traces. Ca²⁺ current decay was visibly more rapid than decay for Ba²⁺ or Sr²⁺ currents. Simulations: I_CaL 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 Ba²⁺ 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 Ca²⁺ as charge carrier (CDI and VDI are concurrent). FRC for CDI+VDI was significantly smaller at more hyperpolarized potentials (V_m = −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. (black squares), Fulop et al. (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 (I_CaL, 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. Ca²⁺ ions bind to CaM (on-rate k₁ and off-rate k_-1). With Ca²⁺ ions bound, the Ca²⁺/CaM/channel complex may activate CDI mode (asterisk and black color indicate CDI activation, on-rate k₂ and off-rate k_-2).

Figure 2. Undiseased human I_to, I_NaCa, and I_K1 experiments and model validation.

A) I_to. Experimental data are white circles (N = 8 from 5 hearts for inactivation time constants, N = 10 from 5 hearts for recovery time constants, N = 9 from 6 hearts for steady state inactivation, and N = 23 from 8 hearts for the I–V curve). Simulation results are solid lines. From left to right, top to bottom: steady state activation, steady state inactivation, fast time constant for inactivation, slow time constant for inactivation (insets show fast and slow recovery from inactivation), relative weight of the fast component for inactivation and the I–V curve (normalized). B) I_NaCa. Experimental data are digitally averaged time traces (N = 3 from 2 hearts, white circles, gray is standard error of the mean). Simulation results are the solid line. Top) Voltage clamp protocol. Bottom) I_NaCa in response to the clamp. C) I_K1. Experimental data are previously unpublished (white circles, N = 21 from 12 hearts), from Bailly et al. (black squares) and Konarzewska et al.(black triangles). Simulation results are solid lines (black, gray and dashed black for [K⁺]_o = 4, 8 and 20 mM). Top left) Voltage and [K⁺]_o dependence of steady state rectification. Top right) Voltage and [K⁺]_o dependence of steady state inactivation. Bottom left) Time constant for inactivation. Bottom right) I–V curve, and its [K⁺]_o dependence.

Figure 3. Undiseased human I_Kr and I_Ks experiments and model validation.

A) I_Kr. Experimental data are white circles (N = 10 from 7 hearts for steady state activation, N = 7 from 3 hearts for activation and from 2 hearts for deactivation time constants and weights, and N = 10 from 7 hearts for tail currents). Simulation results are lines. From left to right, top to bottom: steady state activation, time constant for activation (fast (solid) and slow (dashed) time constants converge), fast time constant for deactivation, slow time constant for deactivation, relative weight of the fast component for deactivation, and the I–V curve for normalized tail currents. B) Activation/deactivation profiles in response to the voltage steps shown (−40 mV holding potential to +30 mV steps of various duration, followed by a return to −40 mV, top right inset). Experiments are above. Simulations are below. Activation is rapid, occurring within tens of milliseconds. Deactivation is slow, occurring after several seconds. C) Human AP clamp waveform (top), used to elicit 1 µM E-4031 sensitive current (I_Kr, bottom); experiments are on the left, and comparison to simulations using the same AP clamp is on the right. D) I_Ks. Data are from Virág et al. (black circles). Simulation results are solid lines. From left to right: steady state activation, time constant for activation (much slower than deactivation at depolarized potentials), time constant for deactivation (much faster than activation at hyperpolarized potentials), and the I–V curve, showing normalized tail currents.

Figure 4. Nonfailing human fast I_Na and late I_Na experiments and model validation.

A) Fast I_Na. Experiments are from Sakakibara et al. (black squares) and Nagatomo et al. (black triangles). Simulation results are solid lines. From left to right, top to bottom: steady state activation, time to peak (experiment) and activation time constant (simulation), steady state inactivation, fast time constant for development of inactivation, slow time constant for development of inactivation, time constant for recovery from inactivation, and the I–V curve (solid line simulation and data at 17°C, dashed line simulation at 37°C). In other panels, simulations and data were adjusted to 37°C. Time to peak was fit at 33°C. B) Late I_Na. Experiments are from Maltsev et al. (black squares). Simulation results are solid lines. Top) Steady state activation. Middle) Steady state inactivation. Bottom) I–V curve.

Figure 5. Schematic diagram of human ventricular myocyte model.

Formulations for all currents and fluxes were based either directly (gray) or indirectly (white) on undiseased or nonfailing human experimental data. Model includes four compartments: 1) bulk myoplasm (myo), 2) junctional sarcoplasmic reticulum (JSR), 3) network sarcoplasmic reticulum (NSR), and 4) subspace (SS), representing the space near the T-tubules. Currents into the myoplasm: Na⁺ current (I_Na; representing both fast and late components), transient outward K⁺ current (I_to), rapid delayed rectifier K⁺ current (I_Kr), slow delayed rectifier K⁺ current (I_Ks), inward rectifier K⁺ current (I_K1), 80% of Na⁺/Ca²⁺ exchange current (I_NaCa,i), Na⁺/K⁺ pump current (I_NaK), background currents (I_Nab, I_Cab, and I_Kb), and sarcolemmal Ca²⁺ pump current (I_pCa). Currents into subspace: L-type Ca²⁺ current (I_CaL, with Na⁺ and K⁺ components I_CaNa, I_CaK), and 20% of Na⁺/Ca²⁺ exchange current (I_NaCa,ss). Ionic fluxes: Ca²⁺ through ryanodine receptor (J_rel), NSR to JSR Ca²⁺ translocation (J_tr), Ca²⁺ uptake into NSR via SERCA2a/PLB (J_up; PLB - phospholamban), diffusion fluxes from subspace to myoplasm (J_diff,Na, J_diff,Ca, and J_diff,K). Ca²⁺ Buffers: calmodulin (CMDN), troponin (TRPN), calsequestrin (CSQN), anionic SR binding sites for Ca²⁺ (BSR), anionic sarcolemmal binding sites for Ca²⁺ (BSL). Ca²⁺/calmodulin-dependent protein kinase II (CaMK) and its targets are labeled.

Figure 6. Undiseased human endocardial AP traces from experiments (small tissue preparations) and model simulations.

Simulated APs for a range of pacing frequencies (top) and corresponding examples of experimentally recorded APs at 37°C (below). Arrows indicate cycle length (CL) changes. B) Comparison of simulation (black) and experimentally measured (gray, small tissue preparations) basic AP parameters for a single paced beat from quiescence (37°C, N = 32 from 32 hearts). Shown, from top to bottom, are the resting membrane potential (V_m rest), maximum upstroke potential (V_m max), and maximum upstroke velocity (dV_m/dt max).

Figure 7. Undiseased human endocardial AP response to pacing protocols from experiments (small tissue preparations) and model simulations.

A) Steady state APD rate dependence. B) S1S2 APD restitution (DI – diastolic interval). APD30–90 are labeled at right. Solid lines are simulation results; white squares are experiments at 37°C (N = 140 hearts in panel A, N = 50 hearts in panel B). C) Repolarization rate at CL = 1000 ms. Trajectory of APD30 to APD90 is accurate in the ORd model (white squares are experimental data); less so in other models. D) Dynamic restitution protocol (see Methods). Experiments are from Koller et al., measured in nonfailing human hearts with monophasic AP electrodes (black squares). Simulated results are the black line. At very short diastolic intervals (DI<90 ms), the model shows APD bifurcation (alternans).

Figure 8. Pacing protocols with block of various currents.

Experimental data (small tissue preparations) are white squares. A) Steady state APD90 rate dependence. From left to right, top to bottom: N = 140, 5, 5, 5, 5, 4, and 4 hearts. Shown are control, I_Kr, I_Ks, I_CaL, I_K1, RyR, and late I_Na block. B) APD90 restitution (S1 = 1000 ms). From left to right: N = 50, 3, and 4 hearts. Shown are control, I_Kr, and I_Ks block.

Figure 9. Rate dependence of currents at steady state.

Black arrows indicate CL decrease (rate increase). Top Row) Simulated APs, repeated in each column for timing purposes. Lower Rows (left to right, top to bottom): I_Na, peak I_Na detailed time course, late I_Na, I_to, I_CaL, I_CaL increasing peaks with increasing pacing rate, I_Kr, I_Ks, I_K1, I_NaCa,i, I_NaCa,ss, and I_NaK. Insets show greater detail of late small I_to window current, and early I_Kr spiking at fast rates.

Figure 10. Transmural heterogeneity and validation of transmural cell type models.

A–C) Expression ratio in the model (black bars) compared to experimental data from undiseased human hearts (grayscale bars, labeled). D) Transmural heterogeneity of I_to; simulations are lines, experiments are squares. Results for endo are gray; those for epi are black. E1) Rate dependence of APD90 in endo, M, and epi cell types. Epi and M data were obtained by scaling endo data (white squares) by epi/endo and M/endo APD90 ratios from Drouin et al. (black squares). Simulations are black lines. E2) Same format as panel E1, showing epi APD90 at faster pacing rates. Data are from Glukhov et al., (epi/endo scaling, black triangles). F) Top to bottom: Rate-dependence of endo, M, and epi APs. G) Pseudo-ECG, using a simulated transmural wedge. CL changes are indicated by arrows.

Figure 11. Early afterdepolarizations (EADs).

A) Top left) Experiments in isolated nonfailing human endo myocytes from Guo et al. showed EADs with slow pacing (CL = 4000 ms) in the presence of I_Kr block (0.1 µM dofetilide, ∼85% I_Kr block, reproduce with permission). Top right) Following the experimental protocol of Guo et al. (CL = 4000 ms, 85% I_Kr block) the ORd model accurately showed a single large EAD. Bottom) GB (left) and TP (right) models failed to generate EADs (CL = 4000 ms, even with 100% I_Kr block). B) EAD mechanism. APs are on top. I_CaL (black) and I_CaL recovery gate (gray) are below. Slow pacing alone (CL = 4000 ms) did not cause an EAD (left). Slow pacing plus I_Kr block (85%) caused an EAD (solid lines, right). The EAD was depolarized by I_CaL reactivation during the slowly repolarizing AP plateau (solid lines, solid arrows). When I_CaL recovery was prevented, the EAD was eliminated (dashed lines and dashed arrow).

Figure 12. Rate dependence of intracellular ion concentrations.

Simulation results are solid lines. A) [Na⁺]_i versus pacing frequency (normalized to 0.25 Hz). Experiments are from nonfailing myocytes (Pieske et al., black squares). B) Peak Ca²⁺ transient (normalized to 0.5 Hz). Experiments are from nonfailing myocytes (Schmidt et al., black squares). C) Ca²⁺ transients from experiments (Fura-2-AM) and simulations. Results are normalized to illustrate the time course of decay. The arrow indicates pacing CL changes. D) Frequency dependent acceleration of relaxation. Undiseased human experimental data are white circles. Simulations are the black line.

Figure 13. CaMK and Ca²⁺ cycling.

A) Rate dependence of CaMK active fraction. B) Ca²⁺ cycling under control conditions (left) and without CaMK (right). CL changes are indicated by arrows. Top) [Ca²⁺]_i and diastolic values (inset). Middle) [Ca²⁺]_JSR. Bottom) J_up and J_rel (inset, expanded time scale).

Figure 14. AP and Ca²⁺ 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, J_rel (inset is expanded time scale), [Ca²⁺]_i, [Ca²⁺]_JSR, and J_up. B) Rate dependence of APD (top) and peak [Ca²⁺]_i (bottom) at fast rates (alternans bifurcations disappear without CaMK). C) Same as panel B, but at slower rates (no bifurcations).

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 I_Na, I_to, I_CaL, I_Kr, I_Ks, zoom of plateau I_CaL (dashed box section), I_K1, I_NaCa, and I_NaK are shown.

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 I_Na (dashed gray) or both late I_Na and I_CaL 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 I_Na (dashed gray), I_CaL (dash-dot-dot gray), or late I_Na and I_CaL (solid gray) to their slow rate values (no [Na⁺]_i and [Na⁺]_ss clamping to slow rate values). C) [Na⁺]_i and I_NaK increase with pacing rate under control conditions (left). When [Na⁺]_i and [Na⁺]_ss are clamped to slow rate values, I_NaK 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 I_Na reset – dashed gray, I_CaL reset – dash-dot-dot gray, late I_Na and I_CaL reset – solid gray). Shown in dashed black is resetting late I_Na and I_CaL 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).

Figure 17. I_Kr 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 (xr_fast and xr_slow) 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 I_Kr (bottom) during the S2 beat at different DIs (indicated by arrows). Spiking in I_Kr occurred early during the AP at short DI. Spiking was caused by slow deactivation, increasing availability of I_Kr.

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 [Ca²⁺]_i at steady state for CL = 1000 ms. From left to right, top to bottom: AP (with V_m rest inset at far right), I_Na (inset shows peaks), late I_Na (not present in TP or GB models), I_to (inset shows decay rate), I_CaL (arrow shows ORd peak magnitude; inset shows normalized peaks, which are wide in TP and GB), I_Kr (arrow shows ORd early spike peak magnitude), I_Ks, I_K1, I_NaCa, I_NaK, [Ca²⁺]_i, [Ca²⁺]_i decay rate, and [Na⁺]_i.