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

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: Ca2+ versus voltage dependent inactivation of L-type Ca2+ current (ICaL); kinetics for the transient outward, rapid delayed rectifier (IKr), Na+/Ca2+ exchange (INaCa), 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 Ca2+ (including peak and decay) and intracellular sodium ([Na+]i) in undiseased human myocytes were quantitatively reproduced by the model. Early afterdepolarizations were induced by IKr block during slow pacing, and AP and Ca2+ alternans appeared at rates >200 bpm, as observed in the nonfailing human ventricle. Ca2+/calmodulin-dependent protein kinase II (CaMK) modulated rate dependence of Ca2+ cycling. INaCa linked Ca2+ 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 ICaL were also contributors. APD shortening during restitution was primarily dependent on reduced late Na+ and ICaL currents due to inactivation at short diastolic intervals, with additional contribution from elevated IKr due to incomplete deactivation.

Definitions and Abbreviations reversal potential for ion S (mV) R CG ratio between capacitive and geometric membrane areas (= 2) A geo geometric area (cm 2 ) A cap capacitive area (cm 2 ) v S volume of compartment S (µL) s ∞ steady state value of gate S τ s time constant of gate S (ms) τ s,fast , τ s,slow fast/slow time constant of gate S (ms) A s,fast , A s,slow fraction of channels with gate S undergoing fast/slow process z s valence of ion S I stim stimulus current (µA/µF) I S current through ion channel S (µA/µF) G S ��� maximum conductance of ion channel S (mS/µF) K m,S half-saturation concentration of molecule S (mM) I S maximum current carried through ion channel S (µA/µF) P S permeability to ion S (cm/s) PR Q,S permeability ratio of ion Q to ion S γ S activity coefficient of ion S J S ion flux S (mM/ms) [S] Q concentration of ion S, in sub-cellular compartment Q (mM) [S] maximum concentration of buffer S (mM) myo myoplasmic compartment (also abbreviated by small letter "i") ss subspace compartment (representing submembrane space near t-tubules) SR sarcoplasmic reticulum jsr junctional SR compartment nsr network SR compartment Y compartment Y (e.g. "i", "ss"), as in I NaCa equations

Cell Geometry
Cell geometry was approximated by a cylinder. Cell length (L) was about ten times longer than the radius (Forbes and Sperelakis) 4 .
CaMK active G Na,fast ��������� = 75 mS/μF Slow Delayed Rectifier Potassium Current (I Ks ) x s1,∞ = 1 For, Y ϵ {i, ss} k Na1 = 15 mM, k Na2 = 5 mM, k Na3 = 88.12 mM, k asymm = 12.5 ω Na = 6 • 10 4 Hz, ω Ca = 6 • 10 4 Hz, ω NaCa = 5 • 10 3 Hz k Ca,on = 1.5 • 10 6 mM ms , k Ca,off = 5 • 10 3 Hz q Na = 0.5224, q Ca = 0.1670 Background Currents (I Nab , I Cab , I Kb ) and Sarcolemmal Calcium Pump Current (I pCa ) The formulations for I Nab , I Cab , I Kb , and I pCa were taken from the Hund-Decker-Rudy model 5,6 . I Kb represents small amplitude, rapidly activating K + current observed in the ventricle (I Kp -like 7 or I Kur -like 8 current). The amplitudes of these currents were reduced compared to values used by Decker et al 5 . These choices were made consistent with the following: 1) so that resting [Na + ] i would be similar to values shown in nonfailing human ventricle at 37 °C by Pieske et al. 9 at very slow pacing rates (0.25 Hz), 2) so that the resting [Ca 2+ ] i would be similar to values shown in nonfailing human ventricle at 37° C by Schmidt et al. 10 , and 3) so that the generally lower major current conductances used to match human data in construction of this model would be properly balanced. P Nab = 3.75 • 10 −10 cm/s, z Na = 1 • �I Na + I to + I CaL + I CaNa + I CaK + I Kr + I Ks + I K1 + I NaCa + I NaK + I Nab + I Cab + I Kb

Calcium/Calmodulin-Dependent Protein Kinase (CaMK)
The CaMK model is equivalent to that used in the Hund-Decker-Rudy dog model 5,6 . We assumed that CaMK kinetics are similar in human and dog, in the absence of human ventricle specific measurements.
The time constant for Na + and K + diffusion fluxes are larger than the time constant for Ca 2+ diffusion flux. Physiologically, this amounts to reduced diffusivity for Na + and K + as they exit the subspace.

SR Calcium Release Flux, via Ryanodine Receptor (J rel )
Ca 2+ release channels (ryanodine receptors, RyRs, formulation similar to that in Livshitz et al. 11 ) have been split into two separate populations in this model according to CaMK phosphorylation state, based on observations in dog ventricle 12 . There is a non-phosphorylated release (J rel,NP ) and a CaMK phosphorylated release (J rel,CaMK ). When RyR channels are phosphorylated by CaMK, release amplitude is 1.25 times larger, and the decay time constant is 1.25 times longer. The proportion of the RyR population that behaves in the phosphorylated state is regulated by active CaMK.

ORd Human Model Concentrations and Buffers
In the absence of human ventricle specific measurements, we take Ca 2+ buffering equations and kinetics from the Hund-Decker-Rudy model.

Computational Methodology
Hardware and Software For simulation of the basic human model, we used custom code developed and run using Microsoft Visual C++ 2008 Express Edition on a Windows Vista Dell desktop PC, with an Intel Core 2 Quad processor. Integration was performed as described below (Rapid Integration). We also used custom C++ code run on an array of Dell cluster nodes with 64-bit Intel Xeon processors, running Linux and Sun Microsystems Grid Engine. Execution scripts were written in Python. A fixed time step of 0.01 ms was applied, and the Rush-Larsen Method 16 was used. All simulations were paced to true steady state 17 , unless otherwise noted.
Validation and fitting of individual model components (i.e. time constants, steady state curves) was performed using custom code written in Matlab 2009a running on a Windows Vista Dell desktop PC, with an Intel Core 2 Quad processor. Automated parameter estimation used a sum of least squares objective function, minimized by Matlab functions "fmincon", "ga", and "lsqcurvefit" (interior reflexive Newton's Method for "fmincon" and "lsqcurvefit", genetic algorithm for "ga"). See Matlab documentation for details and references. We used the parallel implementation of "fmincon" and "ga" by opening a matlabpool (size 4). Manual parameter estimation was also used, where minimization was by simple guess and check.

Rapid Integration
The Rush-Larsen Method 16 , applied by Victorri et al. 18 , relies on the assumption that during sufficiently small time intervals, a system of differential algebraic equations becomes effectively uncoupled. One can then readily solve uncoupled differential equations one-by-one to obtain expressions for time evolution of state variables.
Here, identification of sufficiently small time intervals (dt) was determined by comparison to gold standard simulations with fixed dt = 0.005 ms. We match the gold standard when we apply the following rules: 1) dt = 0.005 ms from the start of the stimulus until 25 ms thereafter 2) maximum allowed dt = 1.0 ms 3) dt was adjusted dynamically with changes in membrane voltage, as described in LR1 19 : i. while ∆V ≥ 0.8 mV, dt is reduced tenfold until the condition, ∆V < 0.8 mV, is met (minimum dt = 0.005 ms) Equations for updating gates (e.g. generic gate, s) Equations for updating the n gate, J rel,NP , and J rel,CaMK n = α n • k +2,n k −2,n − �α n • k +2,n k −2,n − n� • exp�−k −2,n • dt� J rel,NP = J rel,NP,∞ − �J rel,NP,∞ − J rel,NP � • exp � −dt τ rel,NP � The Forward Euler Method was applied to update membrane voltage, concentrations, and CaMK trap at each time step.
Using the above method, it took less than one minute of runtime to pace the model to true and accurate steady state at 1 Hz (Microsoft Visual C++ 2008 Express Edition on a Windows Vista Dell desktop PC, with a 2.83 GHz Intel Core 2 Quad processor).

Supplementary Figures
Additional Details for Currents Figure S1. Human I NaCa model faithfully reproduces Kang Figure S3. Human I NaCa model faithfully reproduces Weber et al. 21 intracellular Na + dependence under AP and Ca 2+ clamp conditions. Compare Weber et al. 21 Figure 6A with simulations from our I NaCa model, shown here (same protocols were used   Figure S7. Human I NaK model faithfully reproduces major observations of Nakao and Gadsby 25  Effects of H + , CO 2 and HCO 3 on Na + Handling, I NaK and APD Rate Dependence We kept Clconcentration constant, since the ORd model does not systematically include Clhandling (20 mM intracellular and 100 mM extracellular, as in Decker et al. 5 ). Otherwise, all Crampin and Smith equations were included exactly as described 29 . The I NaK formulation we used, based on Smith and Crampin 23 , includes pH dependence. The simulations below allow I NaK to respond dynamically to pH. Figure S12. H + , CO 2

Effect of KCNE1 Heterogeneity on Transmural I Ks and AP Simulations
Protein forming the β-subunit of I Ks , KCNE1 was measured to be transmurally heterogeneous in the undiseased human ventricle 30 . Western blots showed about two-fold greater intensity for KCNE1 in M-cells compared to epi cells. Considering that KCNE1:KCNQ1 stoichiometry is variable 31 , and that the presence of KCNE1 slows I Ks activation by about five fold and increases I Ks conductance by about five fold, we simulated the effect of KCNE1 transmural heterogeneity on I Ks and the AP. Thus, for heterogeneous KCNE1 simulations in M-cells, I Ks activation was five times slower and conductance was five times greater than in the control M-cell. For epi cells, activation was five times faster, and conductance was five times smaller than in the control epi cell. These conditions are exaggerated (five fold changes compared KNCE1 overabundance to total KCNE1 absence 32 ), showing possible KCNE1 transmural heterogeneity effects in the extreme. As shown, even for the extreme case there was little effect on I Ks or especially on the AP.