Simulation of the Effects of Extracellular Calcium Changes Leads to a Novel Computational Model of Human Ventricular Action Potential With a Revised Calcium Handling

The importance of electrolyte concentrations for cardiac function is well established. Electrolyte variations can lead to arrhythmias onset, due to their important role in the action potential (AP) genesis and in maintaining cell homeostasis. However, most of the human AP computer models available in literature were developed with constant electrolyte concentrations, and fail to simulate physiological changes induced by electrolyte variations. This is especially true for Ca²⁺, even in the O’Hara-Rudy model (ORd), one of the most widely used models in cardiac electrophysiology. Therefore, the present work develops a new human ventricular model (BPS2020), based on ORd, able to simulate the inverse dependence of AP duration (APD) on extracellular Ca²⁺ concentration ([Ca²⁺]_o), and APD rate dependence at 4 mM extracellular K⁺. The main changes needed with respect to ORd are: (i) an increased sensitivity of L-type Ca²⁺ current inactivation to [Ca²⁺]_o; (ii) a single compartment description of the sarcoplasmic reticulum; iii) the replacement of Ca²⁺ release. BPS2020 is able to simulate the physiological APD-[Ca²⁺]_o relationship, while also retaining the well-reproduced properties of ORd (APD rate dependence, restitution, accommodation and current block effects). We also used BPS2020 to generate an experimentally-calibrated population of models to investigate: (i) the occurrence of repolarization abnormalities in response to hERG current block; (ii) the rate adaptation variability; (iii) the occurrence of alternans and delayed after-depolarizations at fast pacing. Our results indicate that we successfully developed an improved version of ORd, which can be used to investigate electrophysiological changes and pro-arrhythmic abnormalities induced by electrolyte variations and current block at multiple rates and at the population level.

Introduction

In the last few decades, computational models have been increasingly used to study biological systems, due to the productive synergy between the in silico approach and the collection of experimental data, as well as to the improvement in computational resources. One of the most productive fields of application is cardiac physiology: modeling has provided physiological insights through the prediction of phenomena and mechanisms, later confirmed or disproved experimentally (e.g., Denis and Yoram, 2001; Carusi et al., 2012; Severi et al., 2014; Krogh-Madsen et al., 2016; Ni et al., 2018; Niederer et al., 2019).

One of the main causes of death all over the world is sudden cardiac death, which is correlated with the basic pro-arrhythmic mechanisms at the level of ion currents and single ventricular myocyte action potential (AP). In order to understand these mechanisms, and taking advantage of the increasing availability of specific data on ionic currents gathered from human cardiomyocytes, several mathematical models were developed to describe the biophysical mechanisms underlying the human ventricular AP. Currently, the “gold standard” for in silico human ventricular cellular electrophysiology is the O’Hara-Rudy (ORd) model (O’Hara et al., 2011). In the last few years, ORd was chosen as the consensus in silico model; it was used in multiple studies (e.g., Passini et al., 2015; Dutta et al., 2017b; Moreno et al., 2017) and selected for regulatory purposes (Dutta et al., 2017a). However, since mathematical models are reduced representations of the real biological systems, there is always room for improvements. In fact, the present work aims to propose a new and improved model of the human cardiac ventricular AP, starting from the ORd model.

When a new/updated model is proposed, two major points must be taken into consideration. Firstly, every model has its own validity, and what is the best model always depends on the aim of the investigation for which the model is used. Secondly, in order to claim that a model is actually an “improvement” of a previous one, it should be at least as effective as the previous model at reproducing the experimental data. This means that when a new specific question is addressed, none of the previous capabilities of the model should be lost. Keeping these points in mind, we developed a model oriented toward a specific research application: the investigation of the dependency of ventricular repolarization on extracellular electrolyte concentrations, particularly ionized Ca²⁺ concentration. This is especially meaningful, e.g., when considering new in vitro experiments with cardiomyocytes and ion concentrations set to values different than those considered for developing and tuning the model. At the same time, we sought to reproduce all the experimental protocols reported in the original ORd paper (O’Hara et al., 2011) and its further developments, such as the novel dynamic model for the rapid delayed K⁺ current (I_Kr) that can capture drug-channel dynamic interactions (Grandi et al., 2010; Dutta et al., 2017a). Thus, we can also propose our new model as a general purpose and up-to-date model of human ventricular AP that can be used in a variety of applications, including drug trials.

It is well-known that extracellular Ca²⁺ concentration ([Ca²⁺]_o) affects the cardiac AP: in fact, an increase of [Ca²⁺]_o shortens the AP while a decrease lengthens it. This has been observed in different species (guinea pig, dog, calf, and human) and different cell types (atrial, ventricular, and Purkinje fibers) (Temte and Davis, 1967; Kass, 1976; Leitch, 1996; Bai et al., 2005; Severi et al., 2009; Nagy et al., 2013; Supplementary Figure S1). Since both AP prolongation and shortening may lead to arrhythmia onset, the repolarization dependency on [Ca²⁺]_o may have important implications in all clinical contexts where electrolyte changes occur, such as haemodialysis therapy (Severi et al., 2008), pathological hypo/hypercalcemia, and head-down bed-rest experiments (Passini et al., 2013). Multiple ionic mechanisms are involved in the action potential duration (APD)-[Ca²⁺]_o dependence, which is why the phenomenon is not completely understood. Computational modeling may help elucidating it by analyzing the single ionic currents involved. However, most of the commonly used human ventricular AP models, which were developed using a single [Ca²⁺]_o value, are not able to reproduce the experimentally observed effects of [Ca²⁺]_o changes on APD. In fact, their APDs often display prolongation instead of shortening, and vice versa (as is the case with both the ORd (O’Hara et al., 2011) and the Grandi (Grandi et al., 2010, models).

When simulating the electrical activity of cardiac cells, it is also fundamental to set the same extracellular ionic concentrations used in experimental protocols, in order to properly compare in silico results and in vitro data. On the contrary, it would be incorrect to use the same concentrations when the analysis of in vivo pathophysiological mechanisms is the ultimate aim of simulation (Severi et al., 2009), since the in vivo and in vitro concentrations are different. The advantage of a model is that the extracellular ionic concentrations can be changed according to the study purpose, since they are parameters. However, this can be done only if the model correctly reproduces the main effects of such parameter changes. The non-physiological behavior of the ORd model (in response to changes of some extracellular ion concentrations) is one clear shortcoming of the model and of its usage for replicating specific in vitro and in vivo conditions. For this reason, we have modified the ORd model so that it can be used to simulate conditions where [Ca²⁺]_o changes and consequently also APD.

This work proposes a new human ventricular model (Bartolucci-Passini-Severi, BPS2020) which corrects the APD-[Ca²⁺]_o dependence of the ORd model, while still reproducing the large amount of experimental human data that was used for its development and validation. We paid particular attention in setting in silico the same extracellular concentrations used for in vitro experiments. Moreover, we constructed a population of models (Britton et al., 2013; Muszkiewicz et al., 2015; Paci et al., 2018a) to show that the BPS2020 model is a reliable baseline to reproduce the variability observed in experimental data from human undiseased heart (Britton et al., 2017). The population was also used to evaluate the ability of the model to reproduce pro-arrhythmic phenomena, such as early and delayed afterdepolarizations (EADs and DADs) and alternans, not systematically observed in experiments.

Materials and Methods

Model Development

Summary of the Model Development Strategy

The BPS2020 model was developed starting from the ORd model. Since one of the motivations of this study was to reproduce the physiological APD-[Ca²⁺]_o dependence, we first targeted the components of the model most likely contributing to this relationship. In particular, our earlier studies, based on the Ten Tusscher-Noble-Noble-Panfilov model (ten Tusscher et al., 2004), indicated that the L-Type Ca²⁺ current (I_CaL) was primarily responsible for APD-[Ca²⁺]_o dependence. Two contrasting mechanisms are involved. When [Ca²⁺]_o is higher, it increases the I_CaL driving force which, by itself, would enhance the current density, and prolong the APD. On the other hand, a larger I_CaL also increases its Ca²⁺-dependent inactivation (CDI), which would reduce the current, and shorten the APD. Since the physiologically-observed outcome is APD shortening, we hypothesized that the increase in CDI must play the predominant role. Based on this hypothesis, the first change to the ORd model was a new I_CaL formulation (see section L-type Ca²⁺ Current), which strengthens the sensitivity of CDI to [Ca²⁺]_o compared to the ORd model. With this formulation, the increase in CDI induced by variations in [Ca²⁺]_o overcomes the effect of the increase in driving force, thus achieving the inverse APD-[Ca²⁺]_o relation.

After changing the I_CaL formulation, an overall adjustment of the Ca²⁺ handling was required, in order to maintain a correct APD rate dependence and a physiological Ca²⁺-induced Ca²⁺ release (CICR), as in the ORd model (see section Changes to the Ca²⁺ Handling). This was followed by an automated parameter optimization (similar to what was done in a previous study for a different model (Fabbri et al., 2017). The cost function was defined based on the consistency between the simulation results obtained with the BPS2020 model and the human experimental data presented in O’Hara et al. (2011) across a variety of protocols. We listed in Supplementary Table S3 the parameters that underwent the automatic optimization. At this stage, we also modified the extracellular K⁺ concentration ([K⁺]_o) to match the one used in the experiments (see section Extracellular K⁺ Concentration and Liquid-Junction Potential Correction). More details on the optimization procedure and the cost function are available in Supplementary Material.

We reformulated a few currents/fluxes compared to the ORd model, as described below (see section L-type Ca²⁺ Current to Other Changes), while for others we only tuned the maximal conductances.

L-Type Ca²⁺ Current

The I_CaL was completely revisited: its original Hodgkin-Huxley formulation was replaced by a new Markov model (Figure 1A, top panel), based on the structure proposed by Decker-Rudy for canine epicardial cells (Decker et al., 2009). This new formulation separates voltage dependent inactivation (VDI) and CDI in two loops, each consisting of the same four states: one closed (C), one open (O) and two inactivated (I₁ and I₂) states. Each loop includes four transitions: activation (from C to O), fast inactivation (from O to I₁), slow inactivation (from I₁ to I₂) and recovery (from I₂ to C). Activation and recovery rates (α/β and ψ/ω, respectively) are the same in the two loops; they were mainly derived from the ORd time constant and steady state values of the corresponding gating variables. In particular, the activation rates (α/β) were derived from the formulation of the d activation gate of ORd (see equations in Supplementary Material, section 3.4.3.4). The rates which govern the recovery from inactivation (ψ/ω), were derived from the steady state voltage-dependence of the f inactivation gate of ORd, and from a time constant tuned in order to fit the recovery from inactivation voltage-clamp data (Figure 1E). The γ/δ rates were based on the ORd inactivation gate (f) as well, but their range was limited to 0.2 to take into account the incomplete fast inactivation typical for I_CaL and the time constant was adjusted to reproduce voltage-clamp data. Finally, the rates η/θ modulate the slow inactivation from I₁ to I₂, and the equations were obtained by imposing the microscopic reversibility to the Markov model (see equations in Supplementary Material, section 3.4.3.7). Fast and slow inactivation rates (γ/δ and η/θ, respectively) in the CDI loop are K_CDI times faster than the ones in the VDI loop. This closely reflects the observation by Kim et al. (2004), that CDI works as a faster VDI, activated by elevated Ca²⁺. The CDI and VDI loops are connected by up/down rates (r_up/r_down), modulated by intracellular Ca²⁺ concentration, and controlled by the n gate. The K_CDI factor was initially set to 10, that is CDI was supposed to be 10 times faster than VDI. In fact, no published data are available which directly quantify this factor in the undiseased human ventricle. O’Hara et al. (2011), who determined the CDI time constants in an attempt to represent the shape and magnitude of their fractional remaining current (FRC) measurements, came to a formulation for the CDI fast time constant that is about (the ratio is voltage dependent there) six times faster than the VDI fast time constant. In our model, K_CDI was then automatically optimized, getting the final value of 9.

FIGURE 1

Figure 1. Summary of the design and validation of the new I_CaL model. (A) Schematic representation of the Markov model structure: voltage-dependent inactivation (VDI) and Ca²⁺-dependent inactivation (CDI) are represented as two separate loops, with four states each (C, closed; O, open; I₁ and I₂, inactivated). Inactivation rates in the CDI loop are K_CDI times faster than the corresponding rates in the VDI loop. VDI and CDI loops are interconnected by up/down rates (r_up/r_down), dependent on the n-gate, which directly depends on intracellular Ca²⁺ and its binding to Calmodulin (CaM), as shown in the equation at the bottom of the panel, modified from Decker et al. (2009); O’Hara et al. (2011), and Passini and Severi (2013). (B,C) Comparison of the simulated I_CaL I–V curve (B) and steady state inactivation curve (C) between ORd (light blue), BPS2020 (dark blue), and the experimental data from Magyar et al. (2000) (black squares). (D) Evaluation of the effect of CDI inactivation, by comparing VDI-only and VDI+CDI voltage clamp protocols for I_CaL: both ORd (light blue, left) and BPS2020 (dark blue, right) show a faster inactivation when CDI is included, in agreement with experimental recordings with and without Ba²⁺ from O’Hara et al. (2011) (not shown). (E) Comparison of the recovery from inactivation between ORd (light blue), BPS2020 (dark blue), and the experimental data from Fülöp et al. (2004) (black diamonds), obtained using the P1/P2 protocol.

In the ORd model, the n gate represents the fraction of channels operating in CDI mode. It is the only I_CaL state variable directly dependent on Ca²⁺ concentration, and its formulation is based on the interaction between Ca²⁺ and the Calmodulin (CaM) bound to I_CaL channels (Figure 1A, bottom panel): when four Ca²⁺ ions bind to CaM (k₁/k_–1 rates), the Ca²⁺-CaM complex may activate CDI (k₂/k_–2 rates). In the BPS2020 model, the n gate controls r_up/r_down, thus modulating the fraction of channels in the VDI and CDI loop. Kinetics rates K_mn and k_–2n (see Supplementary Material) were modified compared to ORd, to increase the sensitivity of the n gate to Ca²⁺. All I_CaL equations are included in Supplementary Material.

The new I_CaL model was validated using four different I_CaL voltage-clamp protocols, the same used for the validation of the original ORd model. Simulation results obtained with the BPS2020 model were compared with the ones obtained with the ORd model, as well as with the corresponding experimental data.

The I–V and steady state inactivation curves were compared with data from Magyar et al. (2000) (Figures 1B,C). Results with BPS2020 and ORd are almost overlapped, and in agreement with the data. Simulations in VDI-only conditions (no CDI) are in agreement between BPS2020 and ORd: both models show a faster inactivation when both CDI and VDI are present (Figure 1D). These results are also in agreement with experimental recordings where Ba²⁺ was used as the current carrier, instead of Ca²⁺ (O’Hara et al., 2011). Recovery from inactivation was evaluated using a P₁/P₂ protocol, as in Fülöp et al. (2004): results with BPS2020 are closer to the experimental data than the ones obtained with ORd (Figure 1E). Our model effectively reproduces also the main features of FRC during the activation protocol for I_CaL (Supplementary Figure S4). In particular, FRC for CDI+VDI does not decrease monotonically with increasing voltage; rather, FRC curves appear to reproduce the experimentally observed “U shape” of the current inactivation. Moreover, VDI-alone is relatively weak, whereas CDI causes additional inactivation, whose maximal value is around 80%. On the other hand, some discrepancies with respect to the ORd data and model can be observed: at negative potentials CDI is very weak and at positive potentials, although the amount of CDI after 50 ms is very similar in the data and in the two models, CDI is a bit slower in ours, especially during the first milliseconds (e.g., it is still negligible after 5 ms).

Changes to the Ca²⁺ Handling

In addition to the I_CaL formulation, other changes in Ca²⁺-handling were needed to refine the BPS2020 model. The SERCA pump (Sarco-Endoplasmic Reticulum Ca²⁺ ATPase, J_up) and the background Ca²⁺ currents were increased by factors of 3.13 and 4, respectively. The Ca²⁺ diffusion from subspace to bulk myoplasm was speeded up, by reducing the corresponding time constant. The sarcoplasmic reticulum (SR) was reduced to a single compartment, and its total volume was reduced by 5%, in agreement with other human computational models (Grandi et al., 2010; Paci et al., 2017) (see also Discussion, section APD Changes Induced by Extracellular Ca²⁺ Variations and Implications on Intracellular Ca²⁺ Handling). Consequently, Calsequestrin concentration ([CSQN]) was re-distributed in the whole SR: this reduced its concentration from 10 to 1 mM.

The SR Ca²⁺ release flux (J_rel) via the ryanodine (RyR)-sensitive channels was replaced by the phenomenological formulation used by Paci et al. (2018b), and originally proposed by Koivumäki et al. (2011) for the human atrial myocyte. This formulation employs a Hodgkin-Huxley formalism, and it allows reproducing delayed after-depolarizations (DADs), which could not be reproduced with the original ORd formulation. Minor parameter changes were made to adapt the formulation to the different Ca²⁺ concentrations of BPS2020 (see Supplementary Material, section 3.4.4).

Extracellular K⁺ Concentration and Liquid-Junction Potential Correction

All the in vitro data published by O’Hara et al. (2011) were obtained with a [K⁺]_o = 4 mM, while all the ORd simulations were done with [K⁺]_o = 5.4 mM. In addition, as reported by O’Hara in the reader comments online (O’Hara et al., 2011), the in vitro experiments did not correct for the liquid-junction potential (LJP) at the electrode tip. In fact, in the I_K1 experiments from O’Hara et al. (2011) the apparent reversal potential corresponded to [K⁺]_o = 5.4 mM despite the actual experimental concentration of 4 mM. The Nernst potential difference between the cases of 4 and 5.4 mM extracellular K⁺ is approximately +8 mV. Therefore, we introduced a positive shift equal to +8 mV (approximately corresponding to LJP) in the current equations built using the experimental data with [K⁺]_o = 4 mM. The affected equations are: the steady-state activation/inactivation and time constant of the transient outward K⁺ current (I_to), I_K1 formulation, and the steady state activation/deactivation and time constant of I_Ks. In addition, we set [K⁺]_o = 4 mM in all the simulation protocols which are compared to these experimental data (e.g., rate dependence, restitution, drug block effects). The extracellular Na⁺ concentration ([Na⁺]_o) was also changed from 140 to 144 mM.

Na⁺/Ca²⁺ Exchanger

The Na⁺/Ca²⁺ exchanger (I_NaCa) formulation is the same as in the ORd model. We increased the maximum conductance by a factor of 2.4 upon automatic optimization.

Rapid Delayed Rectifier K⁺ Current (I_Kr)

During the development of the BPS2020 model, a modification of the ORd model was published by Li et al. (2017) and Dutta et al. (2017a) within the Comprehensive in Vitro Proarrhythmia Assay (CiPA) initiative (Sager et al., 2014; Colatsky et al., 2016). This model was optimized for simulations of ion channel block and dynamic drug-channel interactions, and it used a new Markov model for I_Kr (Dutta et al., 2017a), that we decided to incorporate in the BPS2020 model. I_Kr maximum conductance was increased by 20% by the parameter optimization.

Inward Rectifier K⁺ Current (I_K1)

The inward rectifier K⁺ current (I_K1) maximum conductance was decreased by 29% and the steady state rectification slope was increased by 9%, following parameter optimization.

Slow Delayed Rectifier K⁺ Current (I_Ks)

The slow delayed rectifier K⁺ current (I_Ks) conductance was doubled, as in Dutta et al. (2017a). This was done to improve APD rate dependence under I_Kr channel block.

Late Na⁺ Current (I_NaL)

The late Na⁺ current (I_NaL) conductance was increased by a factor of 2.8, as in Dutta et al. (2017a).

Fast Na⁺ Current (I_NaF)

The steady state inactivation and recovery from inactivation gates for I_NaF (h_ss and j_ss) were modified as in Passini et al. (2015) (including also the CaMK phosphorylation pathway), to prevent propagation failure in hyperkalemia conditions. The time constants of inactivation gates were also modified, based on Dutta et al. (2017b). I_NaF conductance was reduced by 73% compared to its original value to preserve the peak current magnitude.

Na⁺/K⁺ ATPase Current (I_NaK)

I_NaK maximum current was increased (twofold) to improve APD rate dependence (see Supplementary Figure S16). We consider this change reasonable, since no direct measurements of the Na⁺/K⁺ ATPase current (I_NaK) in non-failing or healthy human ventricular cells are available, as also reported by O’Hara et al. (2011).

Other Changes

The current stimulus duration was set to 1 ms, with −53 μA/μF amplitude (twice the diastolic threshold), as in Passini et al. (2015). This allows to increase the maximum time step for integration to 1 ms.

Single Cell Simulations

Simulations with the BPS2020 model were run using MATLAB R2018a (Mathworks Inc., Natick, MA, United States) on a Win 10 PC with an Intel Core i7. Numerical integration was performed with the Matlab function ode15s, a variable-step, variable-order solver, based on numerical differentiation formulas (Shampine and Reichelt, 1997). Pacing was maintained until steady state AP was reached (1,000 beats), and APD_X was measured once membrane voltage reached X% of the resting value (as in O’Hara et al., 2011). The code for running single cell simulations is available online at https://www.mcbeng.it/en/downloads/software/16-bps2020-model.html

Population of Models

As in Britton et al. (2013) and Paci et al. (2017), a population of non-diseased AP models accounting for biological variability was constructed, by assuming that variability is mostly due to by cell-to-cell differences in ion channel densities rather than kinetics.

An initial random population of 5,000 human endocardial AP models was generated by sampling 11 parameters in the BPS2020 model, i.e., the maximum conductances, currents and fluxes of I_NaL, I_NaF, I_CaL, I_to, I_Kr, I_Ks, I_K1, I_NaCa, I_NaK, J_rel, and J_up. The parameters were sampled in the range [20–200%] of their original values, using the Latin Hypercube Sampling (McKay et al., 1979).

The random population was then calibrated to select those models whose APs were in agreement with the in vitro non-diseased data (O’Hara et al., 2011) (summarized in Passini et al., 2017). The AP biomarkers used for the experimental calibration were: AP duration at 40, 50, and 90% of repolarization (APD₄₀, APD₅₀, and APD₉₀), AP triangulation (Tri90-40, computed as the difference between APD₉₀ and APD₄₀), maximum upstroke velocity (dV/dt_max), AP peak voltage (V_peak) and resting potential (RMP). Each biomarker calibration range was built considering its minimum and maximum in vitro values. As an additional condition for calibration, we required that the models simulated correctly the inverse APD-[Ca²⁺]_o dependence, as detailed in section APD – [Ca²⁺]_o Dependence. Only models showing all the AP biomarkers within their in vitro ranges and an inverse APD-[Ca²⁺]_o dependence were accepted in the final calibrated population. The other models were discarded. As in Paci et al. (2017), we applied further additional conditions to avoid including models with clearly non-physiological intracellular ion concentrations. Only the models that had [Na⁺]_i in the range [5, 15] mM and the sarcoplasmic Ca²⁺ concentration ([Ca²⁺]_SR) in [0, 10] mM (in steady state) were added to the calibrated population.

To assess the occurrence of delayed afterdepolarizations (DADs) we used the same protocol as in Li and Rudy (2011): starting from steady state each model in the population was paced for 1,500 beats at fast pacing (BCL = 300 ms, 3.3 Hz) and then for 1 long beat (BCL = 10,000 ms). DADs were identified visually.

We also used the population approach to study the occurrence of repolarization abnormalities, such as EADs and repolarization failures (RFs). The occurrence of repolarization abnormalities in the whole population was assessed after the administration of 0.1 μM dofetilide at a cycle length (CL) of 4,000 ms. The same extracellular concentration experimentally used by Guo et al. (2011) ([K⁺]_o = 5 mM, [Ca²⁺]_o = 2 mM, [Na⁺]_o = 137 mM) were used. The drug effect was simulated using the I_Kr drug binding values reported by Dutta et al. (2017a).

Finally, we challenged our in silico population with high pacing rate to assess variability in rate dependence and potential alternans occurrence: each model underwent increasing pacing rates, as in the 1998 work by Koller et al. (2017). The tested CLs were {600, 545, 500, 462, 429, 400, 375, 353, 333, 316, 300, 286, 273, 261, 250, 230} ms. For CLs greater than 300 ms models were paced for 30 s; otherwise for 15 s. The extracellular concentrations were set to [K⁺]_o = 5.4 mM, [Ca²⁺]_o = 1.2 mM and [Na⁺]_o = 140 mM. For this protocol, we performed an additional calibration step to filter the models displaying unphysiological concentrations as consequence of fast pacing. The upper bounds of [Na⁺]_i and [Ca²⁺]_SR were expanded by 50% with respect to the values used to calibrate the population in basal conditions. Therefore, models showing [Na⁺]_i and [Ca²⁺]_SR over 22.5 and 15 mM, respectively, were not considered for further analysis. Models were grouped in three classes: (i) adapting to the pacing rate increment without developing alternans (ADAPT), (ii) failing to adapt (ADAPT FAIL), and (iii) developing alternans (ALT).

To generate the population of models and running the DADs, dofetilide and alternans simulations, we used the Taito supercluster of CSC – IT Center for Science (Finland). We used Matlab R2017b on 50 CPU cores, using 3 GB/core memory. Generating the population took 21 h, while running DADs, dofetilide and alternans took 3, 15, and 8 h, respectively.

Results

APD-[Ca²⁺]_o Dependence

The challenge of reproducing the physiological APD-[Ca²⁺]_o dependence in computational cardiac models was first described by Severi et al. (2009). Most of the published models respond with AP prolongation to increases in [Ca²⁺]_o, i.e., opposite of what observed in vitro and in vivo (Leitch, 1996; Severi et al., 2009; Nagy et al., 2013; Supplementary Figure S1), including ORd. The main goal of this work was to change the ORd model in order to achieve a physiological APD-[Ca²⁺]_o dependence, while preserving all its positive properties. This led to the development of the BPS2020 model.

Figure 2A illustrates the behavior of ORd (light blue) and BPS2020 (dark blue) for three different [Ca²⁺]_o with [K⁺]_o = 5.4 mM. When [Ca²⁺]_o is set to the control value (1.8 mM), simulations with BPS2020 and ORd produce very similar APs: APD₉₀ at 1 Hz is 274 ms and 285 ms, for BPS2020 and ORd, respectively (Figure 2A, solid lines). This is in line with human experimental values of 277 ± 8 ms (O’Hara et al., 2011). However, when simulating different [Ca²⁺]_o, an opposite behavior emerges. In both models, an increase in [Ca²⁺]_o causes an increase in driving force (Ψ_Ca) and CDI (Figure 2A, dashed lines). However, in ORd the first mechanism prevails: I_CaL is enhanced and the consequence is AP prolongation. On the contrary, in BPS2020, the increase in Ψ_Ca is compensated by a much larger increase in CDI: this causes a faster I_CaL inactivation (see also Supplementary Figures S2, S3 for the states occupancy change over the course of action potential in the different conditions), leading to AP shortening. Interestingly, both models show larger intracellular Ca²⁺ concentration peaks, due to the increased I_CaL, in the initial part of the AP (∼100 ms), before CDI starts playing a role. The same mechanisms are involved when considering a decrease in [Ca²⁺]_o, which is causing AP shortening in ORd and prolongation in BPS2020 (Figure 2A, dotted lines). Therefore, the BPS2020 model is able to reproduce the physiological APD-[Ca²⁺]_o dependence (Figure 2B), not well simulated by ORd. An additional quantitative analysis of AP biomarkers dependence on extracellular electrolytes ([K⁺]_o and [Na⁺]_o) was performed (see Supplementary Figure S14). The upper panels show that the model correctly reproduces the expected effects of hypo- and hyperkalemia: strong inverse relationship between APD and [K⁺]_o and strong positive dependence of the resting potential. The middle panels show that only a mild change (prolongation) of APD is caused by increasing [Na⁺]_o over a very large range. Finally, the bottom panels show that the inverse APD dependence on [Ca²⁺]_o is even larger for APD₇₀ and APD₅₀ than for APD₉₀, which is consistent with the progressive increase in CDI with [Ca²⁺]_o and therefore the progressive reduction in I_CaL, which exerts its predominant role in the phase 2 of the AP. On the other hand, the slight increase in AP peak is consistent with the progressive increase in I_CaL peak shown in Figure 2A.

FIGURE 2

Figure 2. Comparison of ORd and BPS2020 behavior for [Ca²⁺]_o variations. (A) Simulated action potential (AP) for the original ORd and the BPS2020 models (left and right panels, respectively) for three different [Ca²⁺]_o. In control conditions ([Ca²⁺]_o = 1.8 mM, solid lines), results with the two models are quite similar. However, when [Ca²⁺]_o increases ([Ca²⁺]_o = 2.7 mM, dashed lines) or decreases ([Ca²⁺]_o = 0.9 mM, dotted lines), they behave in two opposite ways. Only BPS2020 reproduces the inverse APD-[Ca²⁺]_o relationship observed experimentally. (B) APD-[Ca²⁺]_o relationship for ORd (light blue) vs. BPS2020 (dark blue). (C) Changes observed in the APD-[Ca²⁺]_o relationship of BPS2020 when restoring I_CaL (dotted line) or J_rel (dashed line) to the original ORd formulations.

To investigate if the inverse APD-[Ca²⁺]_o dependence is due to the new I_CaL formulation, or rather dependent on the new J_rel formulation, we ran two additional simulations with the BPS2020 model, separately restoring I_CaL and J_rel to the ORd formulation (Figure 2C). When restoring I_CaL, the physiological APD-[Ca²⁺]_o dependence was lost (Figure 2C, dotted line), while it was preserved when restoring J_rel (Figure 2C, dashed line).

These results highlight the major role of I_CaL in controlling this phenomenon, and they support our initial hypothesis that the inverse APD-[Ca²⁺]_o dependence is mainly controlled by a stronger CDI.

APD Rate Dependence

Figure 3A compares simulation results obtained with BPS2020 and ORd for the APD rate dependence (left panel) and restitution protocols (right panel), computed as in O’Hara et al. (2011). APD rate dependence was computed at steady state for different CLs, while APD restitution was computed using the S₁/S₂ protocol, i.e., steady state at CL = 1,000 ms (S₁), followed by a single S2 extra-systolic stimulus, delivered at various diastolic intervals (DIs). APDs at 30, 50, 70, and 90% of repolarization were compared against the same human experimental data used for ORd validation (O’Hara et al., 2011). Simulations were run with [K⁺]_o = 4 mM, as in the experiments. BPS2020 is in qualitative and quantitative agreement with the experimental data, both for APD rate dependence and restitution. ORd model is in qualitative agreement with the experimental data, but simulated AP are overall longer and slightly different from those published in O’Hara et al. (2011), due to the changes in [K⁺]_o.

FIGURE 3

Figure 3. Rate dependence properties of BPS2020 vs. ORd. In all panels, simulation results for BPS2020 and ORd are shown in dark blue and light blue, respectively. Experimental data from O’Hara et al. (2011); Pieske et al. (2002) and Schmidt et al. (1998) are shown as black squares, black circles and black diamonds, respectively. (A) Steady state action potential duration (APD) rate dependence (CL – cycle length) and APD restitution obtained with the S₁S₂ protocol (DI – diastolic interval). APDs computed at 30, 50, 70, and 90% of repolarization are labeled on the right. (B) Steady state APD₉₀ rate dependence changes induced by specific current blocks; stars are the APD₉₀ values in control conditions. (C) APD₉₀ restitution changes induced by specific current blocks. (D) [Na⁺]_i (left) and peak [Ca²⁺]_i (middle and right) vs. pacing frequency.

Figures 3B,C illustrate how APD rate dependence (Figure 3B) and restitution (Figure 3C) vary in presence of specific channel blockers. Again, simulations were run with [K⁺]_o = 4 mM for both BPS2020 and ORd, as in the experiments (O’Hara et al., 2011). Some current block percentages were refined compared to the ones used in O’Hara et al. (2011), based on current literature: for E-4031 1 μM and HMR-1556 1 μM, we simulated 70% I_Kr block and 90% I_Ks block, respectively, as in O’Hara et al. (2011); for nisoldipine 1 μM we simulated 100% I_CaL block, as in Walsh et al. (2007); for BaCl₂ 100 μM we simulated 84% I_K1 block as in Biliczki et al. (2002); for mexiletine 1 μM we simulated its multichannel action, i.e., 54% I_NaL block, 20% I_CaL block and 9% I_Kr block, as in Dutta et al. (2017a); for ryanodine 5 μM we simulated 30% J_rel block. Indeed, according to literature data (Zucchi and Ronca-Testoni, 1997; Xu et al., 1998; Thomas and Williams, 2012), 5 μM ryanodine seems a too low concentration to induce 90% block of J_rel, as assumed in O’Hara et al. (2011). E.g., ryanodine at 1 nM–10 μM, in bilayer studies, causes the RyR to open permanently to a sub conductance level (∼half of the fully open state, Rousseau et al., 1987). Moreover, ryanodine effects also depend on time of exposure. As reported in Zucchi and Ronca-Testoni, 1997), the incubation time required for inhibition was 5–10 min with 500 mM ryanodine and 1 h with 1 mM ryanodine. With 5 μM the time required is likely even longer. However, unfortunately, O’Hara-Rudy did not give details about the time of exposure (O’Hara et al., 2011). Thus, the quantification of RyR block by ryanodine is very tricky, but it is likely that 5 μM produces a block level definitely lower than 90%.

Intracellular Na⁺ and Ca²⁺ Rate Dependence

Figure 3D shows how different pacing frequencies affect [Na⁺]_i (left panel) and peak [Ca²⁺]_i (middle and right panels). We compared simulation results against the human experimental data by Pieske et al. (2002) for [Na⁺]_i and Schmidt et al. (1998) for [Ca²⁺]_i, normalized to 0.25 Hz and 0.5 Hz, respectively. To facilitate a quantitative comparison between BPS2020 and ORd, simulated [Ca²⁺]_i data are also shown as non-normalized (Figure 3D right panel). Extracellular concentrations were set to the experimental values: [Na⁺]_o = 152 mM, [K⁺]_o = 3.6 mM and [Ca²⁺]_o = 1.8 mM, for [Na⁺]_i; [Na⁺]_o = 146 mM, [K⁺]_o = 5.9 mM and [Ca²⁺]_o = 2.5 mM for [Ca²⁺]_i. BPS2020 and ORd were first paced at 0.5 Hz until steady-state, and then paced additional 100 beats for each frequency. Results with BPS2020 qualitatively and quantitatively reproduce the experimental data for [Na⁺]_i, and they qualitatively reproduce the ones for [Ca²⁺]_i, even though the decrease in peak [Ca²⁺]_i is observed at a faster pacing rate (4 Hz) than in vitro (2.5 Hz, simulations vs. experiments). Simulations with the ORd model are different from those shown in O’Hara et al. (2011), possibly due to differences in the extracellular concentrations or simulation protocol, since full simulation details are not provided in O’Hara et al. (2011).

In addition, since CaMK is important for controlling rate dependence of Ca²⁺ cycling, its effects are reported in Supplementary Figure S5, the results were in agreement with ORd. Supplementary Figure S6 shows the Ca²⁺ concentration in the SR of the BPS2020 model compared to the Ca²⁺ concentration in the SR junctional and network compartments in the ORd model.

APD Accommodation

APD accommodation, i.e., the time course of APD response to abrupt changes in pacing rate, was measured in human patients by Franz et al. (1988), and highlighted as a marker of arrhythmic risk in a simulation study by Pueyo et al. (2010), performed with the Ten Tusscher-Panfilov model (ten Tusscher, 2006). Figure 4 shows a comparison between the in vivo experimental recordings (Figure 4A), and the corresponding simulation results for ORd and BPS2020 (Figures 4B,C, respectively). The BPS2020 model shows a very good qualitative agreement with the experiments, and the APD accommodation curve has a more physiological shape compared to ORd.

FIGURE 4

Figure 4. APD₉₀ accommodation. At t = 0 s, the pacing cycle length (CL) is abruptly reduced from 750 to 480 ms (black circles) or 410 ms (white circles). At t = 180 s, the CL is abruptly increased to its original value. (A) Action potential duration (APD) accommodation measured experimentally by Franz et al. (1988). (B) APD accommodation simulated with the ORd model. (C) APD accommodation with the BPS2020 model. (A,B) Are adapted from O’Hara et al. (2011).

Transmural Heterogeneity

To reproduce transmural heterogeneity, we created two additional versions of the BPS2020 model (EPI and M cells), by scaling specific ionic current conductances in the endocardial version (ENDO) of the model described until now. Changes were based on the same mRNA and protein expression data used for the ORd model (O’Hara et al., 2011), and full details are included in Supplementary Table S1. Simulation results for the three different cell types are shown in Supplementary Figure S10, compared with the human experimental data reported by Drouin et al. (1995) and Glukhov et al. (2010). Experimental data have been rescaled with respect to the ORd ENDO APD₉₀ data (O’Hara et al., 2011), to facilitate the comparison in terms of EPI/ENDO and M/ENDO ratios.

Population of Models (EADs, DADs, and Alternans)

The experimental calibration selected a population of 342 models out of the initial 5,000; their APs and Ca²⁺ transients are shown in Figures 5A,B while the AP biomarker distributions are reported in Figure 5C. This population was used to test the occurrence of repolarization abnormalities and alternans.

FIGURE 5

Figure 5. Experimentally calibrated population. (A) Action potentials and (B) Ca²⁺ transients with the BPS2020 model (baseline, white traces), the experimentally calibrated population (blue traces, representing 342 APs) and the rejected models (gray traces). (C) Biomarker distributions in the experimentally calibrated population. The black vertical lines are the experimental boundaries Passini et al. (2017) used to calibrate the population. APD, AP duration at the % specified repolarization; Tri_90–40, APD₉₀-APD₄₀; dV/dt_max, maximum upstroke velocity; Vpeak, AP peak voltage; RMP, resting membrane potential.

Administration of 0.1 μM dofetilide induced early afterdepolarizations (EADs) in 9 models, repolarization failure (RF) in 40 models and simple AP prolongation with no pro-arrhythmic events (REP) in 293 models. Illustrative AP traces for the three groups are reported in Figure 6A, while significant (p < 0.05) differences in the sampled parameters among the REP, EAD, and RF classes are reported in Figure 6B. Interestingly, models belonging to the EAD class showed smaller I_Ks, denoting smaller repolarization reserve, compared to the REP models. The RF models showed smaller I_K1 than the REP ones, in line with the I_K1 role of stabilizing the resting potential. Models developing EADs showed significantly higher I_CaL than the other classes; in fact, it is well-known that one of the mechanisms leading to EADs is actually I_CaL reactivation during phase three of the AP (Zeng and Rudy, 1995; Figure 6B, and expanded version in Supplementary Figure S13). Another outward current affecting AP repolarization is I_NaK, therefore, if it is lower (as in the RF group), it provides less contribution to AP repolarization. However, in Priori and Corr (1990), another mechanism for the occurrence of EADs was highlighted: they concluded that Ca²⁺ release from SR activates an inward current (e.g., I_NaCa) that carries the depolarizing charge for EADs (Figure 6C). These mechanisms underlying EAD generation are in agreement with previous simulation studies (Passini et al., 2015, 2017).

FIGURE 6

Figure 6. Repolarization abnormalities – EADs. (A) Examples of different responses to dofetilide (0.1 μM), CL = 4,000 ms (Guo et al., 2011). Top: repolarizing models (black, REP); middle: models developing early afterdepolarizations (magenta, EAD) and bottom panel failing to repolarize (cyan, RF). (B) Distribution of the scaling factors which show statistically significant differences between the three categories: models repolarizing (REP, black), models developing EADs (EAD, magenta) and models failing to repolarize (RF, cyan) (*p < 0.05). Red crosses represent outliers. (C) EAD induced by I_CaL reactivation (left) and EAD induced by Ca²⁺ release from SR (right).

We also identified models that, in spite their regular AP profiles when paced at CL = 1,000 ms, developed DADs with fast pacing at CL = 300 ms and one long beat at CL = 10,000 ms (Figure 7). Figures 7B,C show two illustrative models developing one DAD and one anticipated AP. In both models the mechanism is the same. Fast pacing induced Ca²⁺ accumulation in the SR (the sarcoplasmic Ca²⁺ concentration in steady state was 1.85 and 1.65 mM, respectively), causing a slow Ca²⁺ efflux from SR before the DAD/anticipated AP by means of the Ca²⁺ leakage flux from SR (J_leak), that slowly increased [Ca²⁺]_i and the Ca²⁺ concentration in the subspace ([Ca²⁺]_SS). Consequently, the RyR-sensitive channels sensed an increased [Ca²⁺]_SS and opened spontaneously (transition 0 to 1 of the RyR_o gating variable), allowing a Ca²⁺ efflux from SR via J_rel. This Ca²⁺ release triggered a remarkable inward current by I_NaCa that depolarized the membrane potential, thus resulting in the DAD/anticipated AP.

FIGURE 7

Figure 7. Repolarization abnormalities – DADs. (A) Illustrative action potential (AP) traces for four models that produced delayed afterdepolarizations (DADs). (B) Illustrative model producing a DAD. (C) Example of DAD degenerating into an anticipated spontaneous AP. In both models the leakage J_leak from the overloaded SR increased the Ca²⁺ concentrations in cytosol ([Ca²⁺]_i) and subspace ([Ca²⁺]_SS). The RyR-sensitive channels sensed the increased [Ca²⁺]_SS and triggered a spontaneous SR Ca²⁺ release (through J_rel) that was translated by the Na⁺/Ca²⁺ exchanger (I_Naca) into the depolarization of the membrane potential, thus determining the DAD and the anticipated AP.

The fast pacing rate protocol we used to trigger alternans resulted in 287 ADAPT models, 28 ADAPT FAIL models and 14 ALT models. 13 models were excluded from this analysis since their intracellular ion concentrations were out of the boundaries defined in the Materials and Methods section. Figure 8 shows three illustrative models for the ADAPT, ADAPT FAIL, and ALT groups. In Figures 8A,B, a representative ADAPT model (the baseline BPS2020) shows a progressive shortening of APD₉₅ in response to the pacing rate increment. Figure 8C,D show a representative ALT model that produces a bifurcation. Finally, Figures 8E,F show a representative ADAPT FAIL model whose AP fails to adapt for CLs shorter than 250 ms. In ADAPT FAIL models, the next stimulus is overlapped more and more to the previous AP. In the 14 models showing alternans, we observed that J_rel failed to recover from inactivation (Supplementary Figure S11), or that the pacing failed to trigger I_Na (and consequently I_CaL, Supplementary Figure S12). In both cases we observed one Ca²⁺ transient every other AP. We also tested whether a J_up increment in the ALT models would have suppressed this phenomenon, since Cutler et al. (2009) showed that 30% upregulation of J_up suppressed alternans. In fact, following a J_up increment, 5 of the 14 ALT models stopped developing alternans, 1 failed to adapt, 5 still developed alternans at higher pacing rate and in 2 no differences were observed. One model was discarded due to out-of-range intracellular concentrations. Therefore, our populations of models confirmed the role of J_up in suppressing/reducing alternans, as also showed by O’Hara et al. (2011). Figures 8C,D (magenta line) illustrate the alternans protocol results obtained with 30% J_up upregulation. Figure 8G shows the distribution of the maximum conductances/currents that are significantly different from the ADAPT class (p < 0.05).

FIGURE 8

Figure 8. Repolarization abnormalities – alternans. Action potentials (APs) at different cycle lengths (CLs) (pacing at 30 s if CL ≥ 300 ms, otherwise 15 s) and APD₉₀-CL relationship for three models from the in silico population. The model (black) in (A,B) belongs to the ADAPT class. The model (green) in (C,D) from the ALT class showed alternans and produced a bifurcation. (E,F) Show an ADAPT FAIL model (red, whose AP fails to adapt for CL shorter than 250 ms). The magenta trace in (C,D) show the alternans suppression due to 30% J_up upregulation. (G) Distribution of the scaling factors showing statistically significant differences (^∗p < 0.05) between the three categories: models adapting to changes in pacing rate (ADAPT, black), models failing to adapt (ADAPT FAIL, red) and models developing alternans (ALT, green). Red crosses represent outliers.

Discussion

In this study we present the BPS2020 model, an updated version of the ORd human ventricular AP model (O’Hara et al., 2011): our primary aim was to develop a model able to simulate the physiological inverse APD-[Ca²⁺]_o dependence, observed both in vitro and in vivo (Severi et al., 2009), and to reproduce the in vitro experiments at the correct [K⁺]_o. The ORd model is currently considered the gold standard among the human ventricular AP model, since it was created using experimental data from over 100 undiseased human hearts, and validated against a plethora of different protocols. However, this model does not reproduce the inverse APD-[Ca²⁺]_o dependence. Therefore, it cannot be used in those cases where the effect of [Ca²⁺]_o variation on the AP is under investigation. The BPS2020 model fills this gap, fixing this mechanism without losing the ability to simulate the in vitro experiments reproduced by the original ORd model (e.g., APD rate dependence, restitution, response to administration of drugs). The major changes in BPS2020 compared to ORd are: (i) a novel Markov formulation for I_CaL; (ii) a single compartment description of SR; (iii) a modified J_rel formulation, as in Paci et al. (2018b); (iv) the dynamic I_Kr formulation from Li et al. (2017); (v) the optimization of the maximum conductance of the other currents to better fit the in vitro data from O’Hara et al. (2011), both in control conditions and under drug block.

APD Changes Induced by Extracellular Ca²⁺ Variations and Implications on Intracellular Ca²⁺ Handling

The APD dependence on [Ca²⁺]_o variations should be considered in all clinical contexts where electrolyte modifications occur, since APD changes are very important triggers for arrhythmia onsets. Unfortunately, most of the published human AP models do not take this dependence into account, and therefore respond in a non-physiological way: [Ca²⁺]_o increases lengthen APD, and vice versa.

In our BPS2020 model, the APD-[Ca²⁺]_o dependence has been corrected: the original L-type Ca²⁺ current has been replaced by a new Markov model, where the CDI sensitivity to intracellular Ca²⁺ has been strengthened. Indeed, we found that the APD-[Ca²⁺]_o inverse relation is mainly mediated by the sensitivity of CDI to [Ca²⁺]_i, in turn modulated by [Ca²⁺]_o variations. Upon changes of [Ca²⁺]_o two counteracting effects are elicited, both affecting I_CaL: an increase in the driving force, which would promote APD prolongation, and an increase in CDI, which would promote a faster I_CaL inactivation and APD shortening. The latter prevails, thus leading to the physiological inverse relation. By quantitatively analyzing CDI using the approach by Limpitikul et al. (2018) we found that: (i) the “CDI strength” (i.e., the maximum amount of current that can be inactivated by CDI) is exactly the same in the ORd and our models (74%); (ii) the CDI sensitivity, i.e., changes in CDI strength around the working point, is definitely different. In the ORd model, the working point is placed into the superior saturation part of the nonlinear CDI curve, making the local sensitivity (slope) almost null. On the contrary, in BPS2020, the local sensitivity curve remains steep for the entire intracellular Ca²⁺ range. These results are shown in Figure 9. Sigmoidal models are quite common to represent physiological mechanisms, and the working point is usually in the high-gain range, rather than in the saturation range. Our results suggest that a poor local sensitivity of the CDI inactivation mechanism could be the main reason why most cardiac models, whose domain of applicability is limited to fixed standard [Ca²⁺]_o = 1.8 mM, fail to reproduce the physiological APD-[Ca²⁺]_o dependence.

FIGURE 9

Figure 9. Ca²⁺-dependence inactivation. Left: framework for isolating the Ca²⁺-dependent inactivation (CDI) (from Limpitikul et al., 2018). Top: voltage step. Middle: [Ca²⁺]_ss example of concentration step. Bottom: I_CaL current in response to voltage activation at baseline [Ca²⁺]_ss (black) and in response to a [Ca²⁺]_ss step (magenta). Right: Steady state CDI (evaluated as a/b) as a function of [Ca²⁺]_ss with ORd (light blue) and BPS2020 model (dark blue), fitted with a Boltzmann curve; CDI quantification are also highlighted for different [Ca²⁺]_o values (magenta for BPS2020 and purple for ORd).

We also tested if some other mechanisms influence the APD-[Ca²⁺]_o dependence, following the same approach used by Grandi et al. (2009), i.e., by removing the Ca²⁺ dependence from the I_Ks formulation, and varying the I_NaCa amplitude (−30% and +100% of its nominal value). In both cases, the correct APD-[Ca²⁺]_o dependence was maintained (see Supplementary Figure S15, left and right panel, respectively), therefore excluding both I_Ks and I_NaCa as relevant factors in modulating APD vs. [Ca²⁺]_o.

By itself, the increased CDI sensitivity in BPS2020 had a negative impact on the restitution curve (S1S2 protocol). This was caused by the relative slow Ca²⁺ diffusion between the two SR compartments (time constant = 100 ms in ORd). When considering short diastolic intervals in the restitution protocol, the junctional SR doesn’t have time to refill, thus leading to a very small Ca²⁺ release. This has not a big effect in ORd, while in BPS2020 – where the sensitivity of CDI to changes in intracellular Ca²⁺ is increased – a smaller Ca²⁺ release leads to a significant decrease of CDI, with a consequent, unphysiological, AP prolongation at short diastolic intervals. A very fast diffusion (as described by a single SR compartment) was sufficient to restore the physiological behavior at all the pacing rates, and this is why in BPS2020 we decided to use a single SR compartment. To enforce this observation we tested the role of Ca²⁺ diffusion within the SR by restoring a SR two-compartments (JNS+NSR) in BPS2020, with the original diffusion time constant used in ORd (100 ms), and with faster diffusion (10 and 1 ms). As shown in Supplementary Figure S9, the physiological restitution curve is not reproduced when considering 100 ms (dashed line), but it is progressively restored when considering smaller time constants (dashed-dotted and solid line, for 10 and 1 ms, respectively). This result supports the new single compartment SR as a reasonable approximation of the diffusion process taking place within the complex SR structure. Indeed, a similar description is used also by other human ventricular models (Grandi et al., 2010; Paci et al., 2017).

At the same time, the assumption of two SR compartments is not biologically justified, since the SR is not anatomically divided into two substructures, and the experimental and modeling results still contain some controversy on the velocity of free Ca²⁺ diffusion within the SR (Bers and Shannon, 2013), with larger diffusion coefficient (quite consistent with a single compartment description) describing better the available data. No detectable [Ca²⁺]_SR gradients between SR release sites and Ca²⁺ uptake sites during normal global Ca²⁺ release have been reported by Bers’ Lab (Picht et al., 2011). Furthermore, potential differences in species have yet to be examined in full. Our simulation results fully support the idea of very rapid availability of SR Ca²⁺ for global release in normal conditions.

In the ORd model, the SR Ca²⁺ release flux J_rel via ryanodine receptor is proportional to the I_CaL current. With such direct dependence of J_rel on I_CaL excitation-contraction coupling models are not able to produce pro-arrhythmic triggers, such as DADs (Fink et al., 2011). For this reason we decided to replace the J_rel equations with a phenomenological formulation (Koivumäki et al., 2011; Paci et al., 2018b). The BPS2020 model is now capable of reproducing complex features of RyR Ca²⁺ release, i.e., the dependence of the release on both [Ca²⁺]_i and [Ca²⁺]_SR and the adaptation of the RyR open probability dependence on [Ca²⁺]_i (Koivumäki et al., 2011). Since the quite unnatural direct dependence of J_rel on I_CaL is usually assumed to confer graded release on phenomenological models of CICR, we tested if it was still possible to reproduce graded Ca²⁺ release with our new formulation of J_rel. This is actually the case for voltages below 0 mV while the model strays from experimental results for positive voltages (see Supplementary Figure S7): this feature of the model deserves further quantitative investigation in the future.

Comparison Against the Experimental Data Used to Validate the ORd Model

We challenged our model by comparing our simulations with all the experimental data that the ORd model simulated in O’Hara et al. (2011). This includes a large variety of protocols: AP morphology in physiological conditions, APD rate dependence and restitution, both in control conditions and with specific channel blockers, transmural heterogeneity, APD accommodation, [Ca²⁺]_i and [Na⁺]_i rate dependence, EADs induced by I_Kr block during slow pacing, and AP and Ca²⁺ alternans induced by high pacing rate.

The BPS2020 model reproduced all these data, with an accuracy at least as good as the ORd model, and – in most cases – better, as shown by comparing the root mean square error to quantify the distance between experimental data and simulations (see Supplementary Table S2).

It is worth noting that many in vitro experiments were actually reproduced significantly better by BPS2020 than ORd, and that our simulations with the ORd model differ from those presented by O’Hara et al. (2011). This is mainly due to an accurate revision of the electrolyte concentrations used in the different simulated conditions. In fact, to achieve a proper comparison between experiments and simulations, it is crucial to simulate the cardiac AP by imposing the same extracellular concentrations used in the corresponding experimental protocol. However, this constraint has often been neglected in the past. All the comparison between experiments and simulations presented here followed this principle. In O’Hara et al. (2011), comparisons with experiments were instead simulated with basal extracellular concentrations ([Na⁺]_o = 140 mM, [Ca²⁺]_o = 1.8 mM and [K⁺]_o = 5.4 mM), although the in vitro data were mainly recorded with [KCl] = 4 mM in the bath solution.

Comparison Against Other Human Models

Although several human ventricle AP models were developed and published in the past two decades, the ORd model is currently considered as the state-of-the-art in the field, based on its extensive validation with human data. Therefore, we have performed a stringent comparison between the two models, while we considered a systematic benchmarking with all the available models beyond the scope of this paper.

Two recent studies presented optimized versions of the ORd model: the first to fit various LQTS profiles (Mann et al., 2016), the second to optimize the ORd model for its use in safety pharmacology, within the CiPA initiative (Dutta et al., 2017a) (ORd CiPA). The scaling factors they proposed for the maximum conductances of the main ionic currents are different from the ones in our model: I_Ks is increased by 5.75 in Mann et al. model (Mann et al., 2016) and by 1.87 in Dutta et al. (2017a). Ours is increased by an intermediate amount (by 2). I_Kr is increased in the BPS2020 model as it is in Dutta et al. (2017a), but the increase in the former is slightly greater. I_NaL is unchanged in the Mann model while it is scaled by 2.661 in the Dutta model; in line with the latter (O’Hara et al., 2011), we incremented its conductance by 2.8. In the present model, I_NaCa and I_NaK were increased (by factors of 2.4 and 2, respectively) as done by Mann et al. (2.95 and 9.12) (Mann et al., 2016). In addition, I_K1 was decreased in the BPS2020 model.

APD rate dependence and APD-[Ca²⁺]_o relationship for ORd, BPS2020, ORd CiPA (Dutta et al., 2017a) and Grandi et al. (2010) models are compared in Supplementary Figure S18. As shown, the ORd CiPA model is very similar to the ORd one, and both ORd CiPA and Grandi cannot reproduce the physiological APD-[Ca²⁺]_o relationship. In addition, Supplementary Figure S17 shows a comparison of the AP and the major ionic currents for BPS2020, ORd, and ORd CiPA; simulations were in single ENDO cells paced to steady state at CL = 1,000 ms.

Another modification of the ORd model has been very recently published by Tomek et al. (2019). The main aim driving the development of that model (simulation of the negative inotropic effect of sodium current block) was different from the one of BPS2020 (APD dependence on [Ca²⁺]_o). As a matter of fact, both models have been used to investigate and integrate into the equations additional knowledge of basic cardiac cell electrophysiology (some of them in the form of model-generated hypotheses), e.g., sensitivity of CDI to intracellular [Ca²⁺], extraction of the activation curve from I–V data, to mention only a couple. However, a rigorous comparison of the two models, as well as the possibility of combining them into a comprehensive AP model, is far beyond the scope of the present work.

Simulation of the Biological Variability of Human Cardiac AP Through a Population of Models

The population of models approach (Britton et al., 2013; Muszkiewicz et al., 2015) has been used for two different purposes: (i) to explore the ability of the model to reproduce the variability observed in human experimental data in control conditions; (ii) to assess the variability in drug response and occurrence of pro-arrhythmic events, following I_Kr block. The first goal was achieved by first varying the maximum conductances of the main ionic currents around their nominal values (Figure 5A), under the assumption that the variability mostly depends on the difference in the ionic channel expressions (Groenendaal et al., 2015). Then, model outputs were calibrated against experimental data in order to select models with AP biomarkers in the physiological range, and discard the others. In this way, each model becomes a representation of a viable cell within the plausible boundaries of the biological variability observed in the population (Britton et al., 2013; Passini et al., 2015). Of note, the population generated from the BPS2020 model by varying maximum conductances in the (20, 200%) range, covers the range of experimentally measured biomarkers, especially those related to the repolarization phase (APD at different level of repolarization and triangulation, Figure 5C). This suggests that (20, 200%) is a good variability range to develop a population of in silico adult cells. This directly links to the second goal. We used our population to assess the capability of the BPS2020 to reproduce phenomena whose occurrence cannot be quantified in terms of mean value (e.g., EADs, DADs, and alternans), and that are not easy to reproduce with the average undiseased cell model. In vivo, EADs, DADs, and alternans are elicited by specific ‘stressfull’ protocols (using drugs, very low or very high pacing rates, etc.). But even in these challenging conditions only few cells show them (for example, as shown in Guo et al., 2011; Coppini et al., 2013). This, of course makes sense, since healthy cells have to be very robust with respect to pro-arrhythmic phenomena occurrence. In the same way, a computational model aiming to reproduce a healthy cardiomyocyte should not necessarily reproduce these arrhythmias and the population of models approach is particularly well suited for these investigations The fact, with our in silico population we were able to obtain three different pro-arrhythmic abnormalities (EADs, DADs, and alternans), providing further validation of the BPS2020 model.

Limitations

The BPS2020 model was extensively validated based on the experimental protocols shown in the original O’Hara et al. paper (O’Hara et al., 2011). However, some limitations can be highlighted, and they should be investigated in more in detail in the future. In the ORd model the I_CaL I–V curve does not perfectly fit the positive potentials and this remains the case for the BPS2020: this shortcoming could lead to I_CaL overestimation at positive potentials. Although the new J_rel formulation, with no direct dependence on I_CaL, was still able to reproduce graded Ca²⁺ release in good agreement for negative voltages, for positive voltages the model strays from experimental data: a detailed investigation of this phenomenon, including a comparison with local control models of Ca²⁺ release, would be necessary but is beyond the scope of the present work. Another limitation is that in our formulation, the RyR opening dependence on Ca²⁺ is very steep (Supplementary Figure S8). Li and Chen (2001) showed that the RyR opening dependence on Ca²⁺ in lipid bilayers is indeed very steep for low Ca²⁺ concentrations, however, our steady-state RyR opening curve is actually even steeper than the experimental one.

The weakness of CDI at negative potentials is a limitation of our model that could be improved in the future. However its impact under AP is likely very limited, since the membrane potential crosses very quickly the voltage interval [−20, −10] mV, firstly during the upstroke and secondly during the late repolarization phase. Besides, the amount of elicited current at −10 mV would be very small (see I–V curve in Figure 1B). The CDI slowness coming from FRC results is another limitation: in the future it could be fixed, further enhancing its role in modulating APD in case of extracellular Ca²⁺ changes.

Due to the relatively poor availability of human cardiomyocytes data, it is reasonable to exploit all those available to the utmost, in order to obtain a model with properties as close as possible to human cells. For this reason, we used all the available human data for model calibration. However, the separation of model calibration and validation (on independent data sets) could be taken into account, once more human data will be available.

Conclusion

The new BPS2020 model correctly simulates the inverse APD-[Ca²⁺]_o dependence, generally not considered by the human ventricular models available in literature, including ORd. BPS2020 accurately simulates a variety of in vitro experiments, and it also takes into consideration the actual electrolyte concentrations used in the experimental setup. These are the experiments used to develop and validate the ORd model, of which BPS2020 represents a modified version. Furthermore, the new Ca²⁺ release formulation enables BPS2020 to produce DADs. In conclusion, the BPS2020 model expands the domain of applicability of the current ORd model, by adding the possibility to explore changes in ventricular electrophysiology induced by electrolyte changes, e.g., effects of hemodialysis or pathological changes in Ca²⁺ concentrations. Therefore, the BPS2020 model can be used to simulate and investigate a variety of conditions, and constitutes what we propose to be deemed an advanced general-purpose model of human ventricular cardiac electrophysiology.

Data Availability Statement

All datasets generated for this study are included in the article/Supplementary Material.

Author Contributions

CB, EP, and SS conceived and designed the study. CB and EP developed and validated the in silico model. MP and JH enabled the access to parallel computing facilities, developed the population of models, and run the in silico experiments on the population. CB, EP, MP, and SS analyzed the in silico data, prepared the figures, and drafted the manuscript. All the authors interpreted the results and revised the manuscript.

Funding

This research did not receive any specific grant from funding agencies in the public, commercial, or not-for-profit sectors. MP was supported by the Academy of Finland (project CardSiPop, Decision No. 307967). EP was supported by an NC3Rs Infrastructure for Impact Award (NC/P001076/1).

Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Acknowledgments

We wish to acknowledge CSC – IT Center for Science, Finland, for generous computational resources. We also thank Dr. Kristina Mayberry for language revision.

Supplementary Material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fphys.2020.00314/full#supplementary-material

References

Bai, C. X., Namekata, I., Kurokawa, J., Tanaka, H., Shigenobu, K., and Furukawa, T. (2005). Role of nitric oxide in Ca²⁺ sensitivity of the slowly activating delayed rectifier K⁺ current in cardiac myocytes. Circ. Res. 96, 64–72.

PubMed Abstract | Google Scholar

Bers, D. M., and Shannon, T. R. (2013). Calcium movements inside the sarcoplasmic reticulum of cardiac myocytes. J. Mol. Cell. Cardiol. 58, 59–66. doi: 10.1016/j.yjmcc.2013.01.002

PubMed Abstract | CrossRef Full Text | Google Scholar

Biliczki, P., Virág, L., Iost, N., Papp, J. G., and Varró, A. (2002). Interaction of different potassium channels in cardiac repolarization in dog ventricular preparations: Role of repolarization reserve. Br. J. Pharmacol. 137, 361–368.

PubMed Abstract | Google Scholar

Britton, O. J., Bueno-Orovio, A., Van Ammel, K., Lu, H. R., Towart, R., Gallacher, D. J., et al. (2013). Experimentally calibrated population of models predicts and explains intersubject variability in cardiac cellular electrophysiology. Proc. Natl. Acad. Sci. U.S.A. 110, E2098–E2105. doi: 10.1073/pnas.1304382110

PubMed Abstract | CrossRef Full Text | Google Scholar

Britton, O. J., Bueno-Orovio, A., Virág, L., Varró, A., and Rodriguez, B. (2017). The electrogenic Na⁺/K⁺ pump is a key determinant of repolarization abnormality susceptibility in human ventricular cardiomyocytes: a population-based simulation study. Front. Physiol. 8:278. doi: 10.3389/fphys.2017.00278

PubMed Abstract | CrossRef Full Text | Google Scholar

Carusi, A., Burrage, K., and Rodríguez, B. (2012). Bridging experiments, models and simulations: an integrative approach to validation in computational cardiac electrophysiology. Am. J. Physiol. Heart Circ. Physiol. 303, H144–H155. doi: 10.1152/ajpheart.01151.2011

PubMed Abstract | CrossRef Full Text | Google Scholar

Colatsky, T., Fermini, B., Gintant, G., Pierson, J. B., Sager, P., Sekino, Y., et al. (2016). The comprehensive in vitro Proarrhythmia Assay (CiPA) initiative — Update on progress. J. Pharmacol. Toxicol. Methods 81, 15–20. doi: 10.1016/j.vascn.2016.06.002

PubMed Abstract | CrossRef Full Text | Google Scholar

Coppini, R., Ferrantini, C., Yao, L., Fan, P., Del Lungo, M., Stillitano, F., et al. (2013). Late sodium current inhibition reverses electromechanical dysfunction in human hypertrophic cardiomyopathy. Circulation 127, 575–584. doi: 10.1161/CIRCULATIONAHA.112.134932

PubMed Abstract | CrossRef Full Text | Google Scholar

Cutler, M. J., Wan, X., Laurita, K. R., Hajjar, R. J., and Rosenbaum, D. S. (2009). Targeted SERCA2a gene expression identifies molecular mechanism and therapeutic target for arrhythmogenic cardiac alternans. Circ. Arrhythmia Electrophysiol. 2, 686–694. doi: 10.1161/CIRCEP.109.863118

PubMed Abstract | CrossRef Full Text | Google Scholar

Decker, K. F., Heijman, J., Silva, J. R., Hund, T. J., and Rudy, Y. (2009). Properties and ionic mechanisms of action potential adaptation, restitution, and accommodation in canine epicardium. Am. J. Physiol. Heart Circ. Physiol. 296, H1017–H1026. doi: 10.1152/ajpheart.01216.2008

PubMed Abstract | CrossRef Full Text | Google Scholar

Denis, N., and Yoram, R. (2001). Models of cardiac ventricular action potentials: iterative interaction between experiment and simulation. Philos. Trans. Math. Phys. Eng. Sci. 359, 1127–1142.

Google Scholar

Drouin, E., Charpentier, F., Gauthier, C., Laurent, K., and Le Marec, H. (1995). Electrophysiologic characteristics of cells spanning the left ventricular wall of human heart: evidence for presence of M cells. J. Am. Coll. Cardiol. 26, 185–192.

PubMed Abstract | Google Scholar

Dutta, S., Chang, K. C., Beattie, K. A., Sheng, J., Tran, P. N., Wu, W. W., et al. (2017a). Optimization of an in silico cardiac cell model for proarrhythmia risk assessment. Front. Physiol. 8:616. doi: 10.3389/fphys.2017.00616

CrossRef Full Text | Google Scholar

Dutta, S., Mincholé, A., Quinn, T. A., and Rodriguez, B. (2017b). Electrophysiological properties of computational human ventricular cell action potential models under acute ischemic conditions. Prog. Biophys. Mol. Biol. 129, 40–52. doi: 10.1016/j.pbiomolbio.2017.02.007

PubMed Abstract | CrossRef Full Text | Google Scholar

Fabbri, A., Fantini, M., Wilders, R., and Severi, S. (2017). Computational analysis of the human sinus node action potential: model development and effects of mutations. J. Physiol. 595, 2365–2396. doi: 10.1113/JP273259

PubMed Abstract | CrossRef Full Text | Google Scholar

Fink, M., Noble, P. J., and Noble, D. (2011). Ca²⁺ -induced delayed afterdepolarizations are triggered by dyadic subspace Ca²⁺ affirming that increasing SERCA reduces aftercontractions. Am. J. Physiol. - Hear. Circ. Physiol. 301, 1–36. doi: 10.1152/ajpheart.01055.2010

PubMed Abstract | CrossRef Full Text | Google Scholar

Franz, M. R., Swerdlow, C. D., Liem, L. B., and Schaefer, J. (1988). Cycle length dependence of human action potential duration in vivo. Effects of single extrastimuli, sudden sustained rate acceleration and deceleration, and different steady-state frequencies. J. Clin. Invest. 82, 972–979.

PubMed Abstract | Google Scholar

Fülöp, L., Bányász, T., Magyar, J., Szentandrássy, N., Varró, A., and Nánási, P. (2004). Reopening of L-type calcium channels in human ventricular myocytes during applied epicardial action potentials. Acta Physiol. Scand. 180, 39–47.

PubMed Abstract | Google Scholar

Glukhov, A. V., Fedorov, V. V., Lou, Q., Ravikumar, V. K., Kalish, P. W., Schuessler, R. B., et al. (2010). Transmural dispersion of repolarization in failing and nonfailing human ventricle. Circ. Res. 106, 981–991. doi: 10.1161/CIRCRESAHA.109.204891

PubMed Abstract | CrossRef Full Text | Google Scholar

Grandi, E., Pasqualini, F. S., and Bers, D. M. (2010). A novel computational model of the human ventricular action potential and Ca transient. J. Mol. Cell. Cardiol. 48, 112–121. doi: 10.1016/j.yjmcc.2009.09.019

PubMed Abstract | CrossRef Full Text | Google Scholar

Grandi, E., Pasqualini, F. S., Pes, C., Corsi, C., Zaza, A., and Severi, S. (2009). Theoretical investigation of action potential duration dependence on extracellular Ca²⁺ in human cardiomyocytes. J. Mol. Cell. Cardiol. 46, 332–342. doi: 10.1016/j.yjmcc.2008.12.002

PubMed Abstract | CrossRef Full Text | Google Scholar

Groenendaal, W., Ortega, F. A., Kherlopian, A. R., Zygmunt, A. C., Krogh-Madsen, T., and Christini, D. J. (2015). Cell-specific cardiac electrophysiology models. PLoS Comput. Biol. 11:e1004242. doi: 10.1371/journal.pcbi.1004242

PubMed Abstract | CrossRef Full Text | Google Scholar

Guo, D., Liu, Q., Liu, T., Elliott, G., Gingras, M., Kowey, P. R., et al. (2011). Electrophysiological properties of HBI-3000: a new antiarrhythmic agent with multiple-channel blocking properties in human ventricular myocytes. J. Cardiovasc. Pharmacol. 57, 79–85. doi: 10.1097/FJC.0b013e3181ffe8b3

PubMed Abstract | CrossRef Full Text | Google Scholar

Kass, R. S. (1976). Control of action potential duration by calcium ions in cardiac Purkinje fibers. J. Gen. Physiol. 67, 599–617.

PubMed Abstract | Google Scholar

Kim, J., Ghosh, S., Nunziato, D. A., and Pitt, G. S. (2004). Identification of the components controlling inactivation of voltage-gated Ca²⁺ channels. Neuron 41, 745–754.

PubMed Abstract | Google Scholar

Koivumäki, J. T., Korhonen, T., and Tavi, P. (2011). Impact of sarcoplasmic reticulum calcium release on calcium dynamics and action potential morphology in human atrial myocytes: a computational study. PLoS Comput. Biol. 7:e1001067. doi: 10.1371/journal.pcbi.1001067

PubMed Abstract | CrossRef Full Text | Google Scholar

Koller, M. L., Riccio, M. L., and Gilmour, R. F. (2017). Dynamic restitution of action potential duration during electrical alternans and ventricular fibrillation. Am. J. Physiol. Circ. Physiol. 275, H1635–H1642.

Google Scholar

Krogh-Madsen, T., Sobie, E. A., and Christini, D. J. (2016). Improving cardiomyocyte model fidelity and utility via dynamic electrophysiology protocols and optimization algorithms. J. Physiol. 594, 2525–2536. doi: 10.1113/JP270618

PubMed Abstract | CrossRef Full Text | Google Scholar

Leitch, S. (1996). Effect of raised extracellular calcium on characteristics of the guinea-pig ventricular action potential. J. Mol. Cell. Cardiol. 28, 541–551.

PubMed Abstract | Google Scholar

Li, P., and Chen, S. R. W. (2001). Molecular basis of Ca²⁺ activation of the mouse cardiac Ca²⁺ release channel (ryanodine receptor). J. Gen. Physiol. 118, 33–44. doi: 10.1371/journal.pone.0184177

PubMed Abstract | CrossRef Full Text | Google Scholar

Li, P., and Rudy, Y. (2011). A model of canine purkinje cell electrophysiology and Ca²⁺ cycling: Rate dependence, triggered activity, and comparison to ventricular myocytes. Circ. Res. 109, 71–79.

Google Scholar

Li, Z., Dutta, S., Sheng, J., Tran, P. N., Wu, W., Chang, K., et al. (2017). Improving the in silico assessment of proarrhythmia risk by combining hERG (Human Ether-à-go-go-Related Gene) channel-drug binding kinetics and multichannel pharmacology. Circ. Arrhythmia Electrophysiol. 10:e004628. doi: 10.1161/CIRCEP.116.004628

PubMed Abstract | CrossRef Full Text | Google Scholar

Limpitikul, W. B., Greenstein, J. L., Yue, D. T., Dick, I. E., and Winslow, R. L. (2018). A bilobal model of Ca²⁺-dependent inactivation to probe the physiology of L-type Ca²⁺ channels. J. Gen. Physiol. 150, 1688–1701. doi: 10.1085/jgp.201812115

PubMed Abstract | CrossRef Full Text | Google Scholar

Magyar, J., Iost, N., Körtvély, A., Bányász, T., Virág, L., Szigligeti, P., et al. (2000). Effects of endothelin-1 on calcium and potassium currents in undiseased human ventricular myocytes. Pflugers Arch. Eur. J. Physiol. 441, 144–149.

PubMed Abstract | Google Scholar

Mann, S. A., Imtiaz, M., Winbo, A., Rydberg, A., Perry, M. D., Couderc, J. P., et al. (2016). Convergence of models of human ventricular myocyte electrophysiology after global optimization to recapitulate clinical long QT phenotypes. J. Mol. Cell. Cardiol. 100, 25–34. doi: 10.1016/j.yjmcc.2016.09.011

PubMed Abstract | CrossRef Full Text | Google Scholar

McKay, M. D., Beckman, R. J., and Conover, W. J. (1979). A comparison of three methods for selecting values of input variables in the analysis of output from a computer code. Technometrics 21, 239–245.

Google Scholar

Moreno, C., Oliveras, A., Bartolucci, C., Peraza, D. A., Gimeno, J. R., Severi, S., et al. (2017). D242N, a KV7.1 LQTS mutation uncovers a key residue for IKs voltage dependence. J. Mol. Cell. Cardiol. 110, 61–69. doi: 10.1016/j.yjmcc.2017.07.009

PubMed Abstract | CrossRef Full Text | Google Scholar

Muszkiewicz, A., Britton, O. J., Gemmell, P., Passini, E., Sánchez, C., Zhou, X., et al. (2015). Variability in cardiac electrophysiology: using experimentally-calibrated populations of models to move beyond the single virtual physiological human paradigm. Prog. Biophys. Mol. Biol. 120, 115–127. doi: 10.1016/j.pbiomolbio.2015.12.002

PubMed Abstract | CrossRef Full Text | Google Scholar

Nagy, N., Acsai, K., Kormos, A., Sebok, Z., Farkas, A. S., Jost, N., et al. (2013). I-induced augmentation of the inward rectifier potassium current (IK1) in canine and human ventricular myocardium. Pflugers Arch. Eur. J. Physiol. 465, 1621–1635. doi: 10.1007/s00424-013-1309-x

PubMed Abstract | CrossRef Full Text | Google Scholar

Ni, H., Morotti, S., and Grandi, E. (2018). A heart for diversity: simulating variability in cardiac arrhythmia research. Front. Physiol. 9:958. doi: 10.3389/fphys.2018.00958

PubMed Abstract | CrossRef Full Text | Google Scholar

Niederer, S. A., Lumens, J., and Trayanova, N. A. (2019). Computational models in cardiology. Nat. Rev. Cardiol. 16, 100–111.

Google Scholar

O’Hara, T., Virág, L., Varró, A., and Rudy, Y. (2011). Simulation of the undiseased human cardiac ventricular action potential: model formulation and experimental validation. PLoS Comput. Biol. 7:e1002061. doi: 10.1371/journal.pcbi.1002061

PubMed Abstract | CrossRef Full Text | Google Scholar

Paci, M., Casini, S., Bellin, M., Hyttinen, J., and Severi, S. (2018a). Large-scale simulation of the phenotypical variability induced by loss-of-function long QT mutations in human induced pluripotent stem cell cardiomyocytes. Int. J. Mol. Sci. 19:3583. doi: 10.3390/ijms19113583

PubMed Abstract | CrossRef Full Text | Google Scholar

Paci, M., Passini, E., Severi, S., Hyttinen, J., and Rodriguez, B. (2017). Phenotypic variability in LQT3 human induced pluripotent stem cell-derived cardiomyocytes and their response to antiarrhythmic pharmacologic therapy: an in silico approach. Heart Rhythm 14, 1704–1712. doi: 10.1016/j.hrthm.2017.07.026

PubMed Abstract | CrossRef Full Text | Google Scholar

Paci, M., Pölönen, R., Cori, D., and Penttinen, K. (2018b). Automatic optimization of an in Silico model of human iPSC derived cardiomyocytes recapitulating calcium handling abnormalities. Front. Physiol. 9:709. doi: 10.3389/fphys.2018.00709

PubMed Abstract | CrossRef Full Text | Google Scholar

Passini, E., Britton, O. J., Lu, H. R., Rohrbacher, J., Hermans, A. N., Gallacher, D. J., et al. (2017). Human in silico drug trials demonstrate higher accuracy than animal models in predicting clinical pro-arrhythmic cardiotoxicity. Front. Physiol. 8:668. doi: 10.3389/fphys.2017.00668

PubMed Abstract | CrossRef Full Text | Google Scholar

Passini, E., Mincholé, A., Coppini, R., Cerbai, E., Rodriguez, B., Severi, S., et al. (2015). Mechanisms of pro-arrhythmic abnormalities in ventricular repolarisation and anti-arrhythmic therapies in human hypertrophic cardiomyopathy. J. Mol. Cell. Cardiol. 96, 72–81. doi: 10.1016/j.yjmcc.2015.09.003

PubMed Abstract | CrossRef Full Text | Google Scholar

Passini, E., Pellegrini, A., Caiani, E., and Severi, S. (2013). Computational analysis of Head-Down Bed Rest effects on cardiac action potential duration. Comput. Cardiol. 40, 357–360.

Google Scholar

Passini, E., and Severi, S. (2013). Extracellular calcium and L-type calcium current inactivation mechanisms: a computational study. Comput. Cardiol. 40, 839–842.

PubMed Abstract | Google Scholar

Picht, E., Zima, A. V., Shannon, T. R., Duncan, A. M., Blatter, L. A., and Bers, D. M. (2011). Dynamic calcium movement inside cardiac sarcoplasmic reticulum during release. Circ. Res. 108, 847–856. doi: 10.1161/CIRCRESAHA.111.240234

PubMed Abstract | CrossRef Full Text | Google Scholar

Pieske, B., Maier, L. S., Piacentino, V., Weisser, J., Hasenfuss, G., and Houser, S. (2002). Rate dependence of [Na⁺]i and contractility in nonfailing and failing human myocardium. Circulation 106, 447–453.

PubMed Abstract | Google Scholar

Priori, S. G., and Corr, P. B. (1990). Mechanisms underlying early and delayed afterdepolarizations induced by catecholamines. Am. J. Physiol. Circ. Physiol. 258, H1796–H1805.

PubMed Abstract | Google Scholar

Pueyo, E., Husti, Z., Hornyik, T., Baczkó, I., Laguna, P., Varró, A., et al. (2010). Mechanisms of ventricular rate adaptation as a predictor of arrhythmic risk. Am. J. Physiol. Circ. Physiol. 298, H1577–H1587. doi: 10.1152/ajpheart.00936.2009

PubMed Abstract | CrossRef Full Text | Google Scholar

Rousseau, E., Smith, J. S., and Meissner, G. (1987). Ryanodine modifies conductance and gating behavior of single Ca²⁺ release channel. Am. J. Physiol. Cell Physiol. 253(3 Pt 1), C364–C368.

PubMed Abstract | Google Scholar

Sager, P. T., Gintant, G., Turner, J. R., and Pettit, S. (2014). Stockbridge, “Rechanneling the cardiac proarrhythmia safety paradigm: a meeting report from the Cardiac Safety Research Consortium. Am. Heart J. 167, 292–300. doi: 10.1016/j.ahj.2013.11.004

PubMed Abstract | CrossRef Full Text | Google Scholar

Schmidt, U., Hajjar, R. J., Helm, P. A., Kim, C. S., Doye, A. A., and Gwathmey, J. K. (1998). Contribution of abnormal sarcoplasmic reticulum ATPase activity to systolic and diastolic dysfunction in human heart failure. J. Mol. Cell. Cardiol. 30, 1929–1937.

PubMed Abstract | Google Scholar

Severi, S., Corsi, C., and Cerbai, E. (2009). From in vivo plasma composition to in vitro cardiac electrophysiology and in silico virtual heart: the extracellular calcium enigma. Philos. Trans. A Math. Phys. Eng. Sci. 367, 2203–2223. doi: 10.1098/rsta.2009.0032

PubMed Abstract | CrossRef Full Text | Google Scholar

Severi, S., Grandi, E., Pes, C., Badiali, F., Grandi, F., and Santoro, A. (2008). Calcium and potassium changes during haemodialysis alter ventricular repolarization duration: in vivo and in silico analysis. Nephrol. Dial. Transplant. 23, 1378–1386.

PubMed Abstract | Google Scholar

Severi, S., Rodriguez, B., and Zaza, A. (2014). Computational cardiac electrophysiology is ready for prime time. Europace 16, 382–383.

Google Scholar

Shampine, L. F., and Reichelt, M. W. (1997). The MATLAB ODE Suite. SIAM J. Sci. Comput. 18, 1–22. doi: 10.1039/b813810a

PubMed Abstract | CrossRef Full Text | Google Scholar

Temte, J. V., and Davis, L. D. (1967). Effect of calcium concentration on the transmembrane potentials of purkinje fibers. Circ. Res. 20, 32–44.

Google Scholar

ten Tusscher, K. H. W. J. (2006). Alternans and spiral breakup in a human ventricular tissue model. AJP Heart Circ. Physiol 291, H1088– H1100.

PubMed Abstract | Google Scholar

ten Tusscher, K. H. W. J., Noble, D., Noble, P. J., and Panfilov, A. V. (2004). A model for human ventricular tissue. Am. J. Physiol. Circ. Physiol. 286, H1573– H1589.

PubMed Abstract | Google Scholar

Thomas, N. L., and Williams, A. J. (2012). Pharmacology of ryanodine receptors and Ca²⁺-induced Ca²⁺ release. Wiley Interdiscip. Rev. Membr. Transp. Signal. 1, 383–397.

Google Scholar

Tomek, J., Bueno-Orovio, A., Passini, E., Zhou, X., Minchole, A., Britton, O., et al. (2019). Development, calibration, and validation of a novel human ventricular myocyte model in health, disease, and drug block. eLife 8:e48890. doi: 10.7554/eLife.48890

PubMed Abstract | CrossRef Full Text | Google Scholar

Walsh, K. B., Zhang, J., Fuseler, J. W., Hilliard, N., and Hockerman, G. H. (2007). Adenoviral-mediated expression of dihydropyridine-insensitive L-type calcium channels in cardiac ventricular myocytes and fibroblasts. Eur. J. Pharmacol. 565, 7–16.

PubMed Abstract | Google Scholar

Xu, L., Tripathy, A., Pasek, D. A., and Meissner, G. (1998). Potential for Pharmacology of Ryanadine Receptor/Calcium Release Channelsa. Ann. N. Y. Acad. Sci. 853, 130–148.

Google Scholar

Zeng, J., and Rudy, Y. (1995). Early afterdepolarizations in cardiac myocytes: mechanism and rate dependence. Biophys. J. 68, 949–964.

PubMed Abstract | Google Scholar

Zucchi, R., and Ronca-Testoni, S. (1997). The sarcoplasmic reticulum Ca2+ channel/ryanodine receptor: modulation by endogenous effectors, drugs and disease states. Pharmacol. Rev. 49, 1–50.

Google Scholar