The O'Hara-Virag-Varro-Rudy (ORd) model (O'Hara et al., 2011) of the human ventricular cardiomyocyte was used as the baseline model for our investigations, as it is one of the most recent, widely used and extensively tested models of the human ventricular cardiomyocyte, and is particularly well-suited for studying human ventricular repolarization, as key currents involved in repolarization and EAD formation (including I_Kr, I_Ks, I_K1, and I_CaL) are parameterized using data exclusively from undiseased human ventricular cardiomyocytes.

We constructed a population of 10,000 human models, each based on the ORd formulation but with differences in the conductance values for 9 key sarcolemmal currents, based on our assumptions that ionic conductances are highly variable, due to both intracellular factors and responses to external stimuli. The conductances varied were those with the most influence on the human ventricular AP: I_Na (fast Na⁺ current); I_NaL (late Na⁺ current); I_CaL (L-type Ca²⁺ current); I_to (transient outward K⁺ current); I_Kr (rapid delayed rectifier K⁺ current); I_Ks (slow delayed rectifier K⁺ current); I_K1 (inward rectifier K⁺ current); I_NCX; (Na⁺/Ca²⁺ exchanger current) and I_NaK (Na⁺/K⁺ pump current). Latin hypercube sampling (McKay et al., 1979) was used to sample sets of conductances uniformly while co-varying all conductances.

As in Britton et al. (2013) we selected a wide sampling range of 0–2 times the baseline value of each conductance to allow models with a wide variety of underlying ionic current configurations, including both increased and reduced channel densities, to be evaluated and potentially accepted into the population if they passed the filtering process. This range is necessarily an assumption as this range cannot be measured in vivo, and voltage clamp measurements of ionic current conductances in isolated cells are affected by the isolation process (Borg and Terracio, 1990). Using a lower sampling limit of 0 also allows the possibility of models with very low, potentially pathological values for some conductances, providing the resulting model can still produce a normal control AP. This choice of limit is intentional, to allow the investigation of abnormal ionic profiles that produce normal APs in control conditions but may have increased susceptibility to repolarization abnormalities under channel block. All conductance values in this study are given as scaling factors of the baseline ORd model's conductances.

Experimental recordings of human AP were obtained as described in O'Hara et al. (2011) from human right ventricular trabeculae and papillary tissue preparations of <2 cm in diameter, dissected from non-failing human hearts, perfused and paced at 1 Hz. Small tissue preparations, rather than isolated cells, were used in this study to avoid the damage to ion channels caused by the isolation process. AP recordings were acquired from these preparations using conventional microelectrode techniques from 62 experiments performed using tissue preparations from 37 hearts (16 females, 21 males, mean age 43 ± 13 years). The recordings used in this study were obtained from a database of recordings from multiple previous studies, all conducted by the same laboratory, which we then further analyzed to calculate biomarker values and experimental ranges. Further details of the preparation process and recording equipment are available in Jost et al. (2005). The preparations used in this study were obtained from hearts donated for research in compliance with the Declaration of Helsinki and approved by the Scientific and Research Ethical Committee of the Medical Scientific Board of the Hungarian Ministry of Health (ETT-TUKEB), under ethical approval No 4991-0/2010-1018EKU (339/PI/010).

The data were analyzed to obtain the experimental ranges of seven biomarkers, which are shown in Figure 1. For each biomarker, experimental ranges were calculated from the minimum and maximum values observed in the data set at 1 Hz control conditions pacing, after excluding three clear outliers, one due to an exceptionally long APD90 (543 ms, data set mean = 293 ± 66 ms) and two due to exceptionally long times to peak (26 and 38 ms, data set mean = 6.7 ± 5.4 ms). The calculated biomarker ranges are displayed in Table S2. Biomarker values were generated for each experiment by calculating the value of each biomarker in the final recorded AP trace of the control portion of each experiment. These experimental ranges were then used to constrain the population of models to only those models that produced normal AP behavior under 1 Hz pacing in control conditions. We refrained from using any additional constraints on the population, to create a population of models with normal APs under normal pacing conditions, but with a wide range of underlying conductances, that we hypothesized would lead to a diverse range of susceptibilities to drug-induced repolarization abnormalities. Performing this filtering led to a total of 568 models of the initial 10,000 being accepted as having normal human ventricular APs; their parameter values are provided in Table S3 in the Supplementary Material.

Simulations were carried out with the 568 models to investigate the occurrence of repolarization abnormalities following block of I_Kr, I_Ks, I_K1, and I_CaL, often considered the key currents that contribute during repolarization (January and Riddle, 1989; Weiss et al., 2010; Varro and Baczko, 2011). We simulated block of all possible two-current pairs of these four currents (6 combinations) and for each pair simulated the 16 possible combinations of 4 block strengths (a 25, 50, 75, or 90% reduction in the appropriate conductance), for a total of 96 simulations performed on each of the 568 models. These simulations represent possible effects of non-selective blockers such as verapamil and dronedarone, both of which block I_Kr and I_CaL (Zhang et al., 1999; Varró et al., 2001), or combinations of selective blockers, e.g., dofetilide and HMR-1556 (Jost et al., 2013). We performed a systematic analysis of block levels, rather than matching specific compounds at set concentrations, to analyze the effects of a wide range of possible block scenarios, and because of the large uncertainties in IC50 values measured for concentration-dependent channel block (Polak et al., 2009).

Repolarization abnormalities (Figure 2) were detected by analyzing the voltage-time gradient of the final two pacing cycles of a simulation. A positive gradient greater than 0.02 mV/ms more than 100 ms after the initial upstroke was used as an automatic detector of afterdepolarizations. We used a 100 ms delay to avoid models with either delayed initial upstroke or a spike and dome configuration from being misclassified as displaying afterdepolarizations.

For each drug block simulation we recorded which models developed repolarization abnormalities. Each model in the population was then categorized by the number of different channel block simulations in which they developed repolarization abnormalities. Models which did not develop repolarization abnormalities in any of the simulations were classified as not susceptible (NS), models which developed repolarization abnormalities in 1–9 of the simulations were classified as moderately susceptible (MS) and models which developed repolarization abnormalities in 10 or more of the simulations were classified as highly susceptible (HS). The threshold separating MS and HS models was selected as it marked a particularly large drop in the number of models displaying repolarization abnormalities in N or more simulations occurred between N = 8 and N = 10 (Figure S1).

We used the Mann-Whitney U-test to determine statistical differences between groups of biomarker values, as it is non-parametric and therefore did not require any assumption of the underlying distributions of biomarker values. Partial correlation coefficients were used to determine correlations between conductances, while controlling for the effects of all other co-varying conductances. Multiple comparisons were corrected for with Bonferroni correction.

A logistic regression model was developed to determine whether models that were highly susceptible to repolarization abnormalities could be classified based on their biomarker values alone. The z-scored values of the seven biomarkers used in this study, obtained from each model in the population, were used as features to train the logistic regression model, which was then used classify models in the test sets as either belonging to the HS category, or not. This classification was then compared to the true categorization. Ten-fold cross validation was used to generate multiple training and testing data sets.

Additional details on simulation protocols and the choice of model used in this study can be found in the Supplementary Material.

Figure 3A shows voltage traces from human models in the population spanning the range seen in experiments, and the baseline ORd model, which has a slightly higher peak membrane potential (45.6 mV) than the highest values seen in our data (39.6 mV), but otherwise would be a viable model in our population. Figure 3B shows that variability in conductances can reproduce the variability in the experimental biomarkers. One possible exception is RMP, which is highly dependent on the extracellular K⁺ concentration, considered constant in our simulations.

Figures 3C–J shows that most of the conductances span the full sampled range (0–2 times the ORd model's baseline value for that conductance), highlighting the robustness of the human AP against conductances variations. However, G_Na, G_CaL, and G_Kr are constrained by the calibration (Figures 3E,G). G_Na is confined to the smallest range, 92% of the models in the population have G_Na values clustered between 0.2 and 1. G_CaL values are densely distributed between 0 and 1 (81% of the population), so that few models have high G_CaL values (e.g., Figure 3E). Low values of G_Kr (below 0.13) are not present in the population (e.g., Figure 3G), suggesting that a minimum amount of I_Kr is necessary in the cell for successful repolarization under normal conditions, but the distribution of G_Kr values spans the rest of the sampled range. All possible pairings of ionic conductances within the population of models are shown in Figure S2. Results demonstrate that a normal AP can be generated by a wide variety of different balances of ionic currents, although these balances may have very different responses to changing conditions.

To test whether the balance of repolarizing currents in the population of models resulted in a realistic range of APD prolongation values under I_Kr block, we compared the APD₉₀ prolongation caused by application of 0.05 μM dofetilide (a selective I_Kr blocker) to the prolongation caused by equivalent I_Kr block on the population of models (Figure 4). Dofetilide was modeled as a single-pore I_Kr inhibitor with a Hill coefficient of 1.2 and an IC50 of 0.03 μM based on data from Kramer et al. (2013), giving a resulting I_Kr block of 65%. The distributions of APD prolongation overlap between experiment and the population. This provides confidence that a realistic balance of I_Kr and other repolarizing currents is captured by the population. There are some outliers in the population that are not seen in experiments, and 7/568 models were excluded from the figure due to developing repolarization abnormalities. This can be explained by the wider range of conductances explored in the population of models compared to the limited set of experimental data in non-diseased hearts.

To investigate our hypothesis that the different conductance profiles in our population would result in different levels of susceptibility to repolarization abnormalities, multiple channel block simulations blocking pairs of the currents I_Kr, I_Ks, I_K1, and I_CaL at a range of block strengths, were conducted as described in Methods. The percentage of models in the population that generated abnormalities in each simulation, are shown in Figures 5A–C for the combinations of currents where abnormalities were detected in at least one simulation. Out of the 96 multichannel block simulations, 29 simulations caused at least one model to develop abnormalities. I_Kr block was the primary cause of repolarization abnormalities. No abnormalities were observed in any simulation where I_Kr was blocked by less than 50%, even though models with very low values of G_Kr were considered in the population. Block of either I_K1 or I_Ks acted to augment I_Kr block and increased the number of models that produced abnormalities (Figures 5A,B), but were not sufficient to produce abnormalities without I_Kr block. I_CaL block attenuated the effects of I_Kr block, reducing the number of abnormalities seen in the population (Figure 5C). At 75% I_Kr block, 50% block of I_CaL was sufficient to abolish all repolarization abnormalities. For 90% I_Kr block, 75% I_CaL block was required to achieve the same effect. In both cases, even 25% block of I_CaL was sufficient to substantially reduce the incidence of repolarization abnormalities in the population, from 4 to 1% of models for 75% I_Kr block, and from 18 to 10% for 90% I_Kr block. Overall, when the population was exposed to both I_Kr and I_CaL block, few models produced abnormalities. However, we found that the APD prolongation caused by I_Kr block was only slightly reduced by I_CaL block (Figure S3).

To analyze how susceptibility to repolarization abnormalities depended on underlying conductance and biomarker values, we classified each model as non-susceptible (NS), moderately susceptible (MS), or highly susceptible (HS) to repolarization abnormalities (see Methods). Out of 568 models, 380 models (67%) were NS, 147 (26%) were MS and 41 (7%) were HS.

We investigated whether the susceptibility of models could be distinguished prior to current block by the control values of any of the seven biomarkers we used to calibrate the population. We found that the median control conditions biomarker values for all biomarkers except for V_m Time were significantly different (Mann-Whitney U-test, p < 0.001) between the sub-population of HS models and the rest of the population (Figure 5D). To determine whether control biomarker values were sufficient to predict which models were highly susceptible to repolarization abnormalities, we constructed a logistic regression model using the biomarkers of the population as predictor variables, to estimate the probability of each parameter set belonging to either the HS or non-HS groups. Over 1,000 iterations of 10-fold cross validation, the logistic regression model only achieved a 61% mean rate of correctly predicting HS models to be in the HS category. This low level of sensitivity suggests that the biomarkers used to calibrate the population are insufficient to predict model susceptibility prior to current block. This is likely due to the high level of overlap in biomarker distributions between the HS and the two other model classifications.

We then tested whether any of the nine conductances that were varied to create the population of models were significantly different between any of the three categories. Median G_Na, G_CaL, and G_NaK values were significantly different between all three categories of model (p < 0.05/27, Figure 5E). G_Na and G_CaL were larger in the higher susceptibility (MS and HS) categories, while G_NaK was smaller. G_NCX was also significantly different between the NS and HS categories, but not between the MS category and either of the other categories. Additionally, the range of G_NCX for all three categories spanned the full sampled range of the conductance, unlike G_Na, G_CaL, and G_NaK (Figure 5E). Therefore, we focused on G_Na, G_CaL, and G_NaK as potential contributors to repolarization abnormality susceptibility.

Using the models in the HS category, we investigated the relative importance of I_Na, I_CaL, and I_NaK, and repolarizing currents I_Kr, I_Ks and I_K1, during early and late periods of repolarization (we define early repolarization as from the AP peak to 40% repolarization, and late repolarization as from 40 to 90% repolarization). Figure 6 displays the average magnitudes of those six currents during each period of repolarization for the 41 HS models, under control conditions and 75% I_Kr block (data for NS and MS models are shown in Figure S4). Under control conditions I_Kr was the largest repolarizing current (mean current density from V_m Peak to APD₉₀: 0.49 ± 0.21 A/F). Unexpectedly, I_NaK was the second largest repolarizing current overall (0.15 ± 0.05 A/F), substantially larger than I_Ks (0.04 ± 0.03 A/F), although I_K1 was larger than I_NaK during late repolarization only (I_K1: 0.18 ± 0.07 A/F, I_NaK: 0.12 ± 0.04 A/F). During 75% I_Kr block, I_NaK became the largest repolarizing current (I_NaK: 0.11 ± 0.04 A/F, I_Kr under 75% block: 0.03 ± 0.02 A/F). I_CaL was a major contributor of inward current during early repolarization under both control conditions and I_Kr block. During late repolarization the magnitude of I_CaL under 75% I_Kr block was increased by over 200% compared to control conditions (control: −0.05 ± 0.03 A/F, I_Kr block: −0.16 ± 0.01 A/F), which supports the importance of I_CaL reactivation in the generation of repolarization abnormalities. I_Na magnitude was not elevated under I_Kr block compared to control conditions (control: −0.007 ± 0.005 A/F, I_Kr block: −0.003 ± 0.004 A/F), suggesting it does not contribute to repolarization abnormality generation directly. However, we found that models with high G_Na correlated with models with low G_NaK (Figure 3C). Therefore, high G_Na may be compensatory for the increased intracellular Na⁺ in models with low G_NaK (mean intracellular Na⁺ was 16.4 ± 6.4 mM in HS models, compared with 8.9 ± 3.9 mM in the rest of the population), and the resulting decrease in inward Na⁺ driving force.

Overall, we identified low G_NaK and high G_CaL as key factors that increased susceptibility to repolarization abnormalities and distinguished highly susceptible models from the rest of the population, and found that I_NaK makes a key contribution to repolarization current, particularly under I_Kr block.

A large majority (38/41) of the HS models in the population had G_NaK values of less than 0.5 times its value in the baseline ORd model, and values of G_CaL greater than 0.78 times the baseline value (Figure 7A). Repolarization abnormality occurrence across the HS models under 75 and 90% I_Kr block was 51 and 95% of models respectively (Figure 7B). After increasing G_NaK to the ORd model's baseline value (i.e., a scaling factor of 1.0) in each HS model, repolarization abnormality occurrence fell to 10 and 71% of models respectively. The reversion of an example HS model to a normal AP phenotype following increase of G_NaK in this way is shown in Figures 7D–G.

We hypothesized that the importance of G_NaK in determining susceptibility to repolarization abnormalities could be due to I_NaK's role in maintaining intracellular Na⁺, the significant contribution of I_NaK to repolarization reserve during I_Kr block (Figure 6), or a combination of both effects. We therefore investigated the effects of intracellular Na⁺ concentration and the electrogenic action of I_NaK individually, to evaluate their contributions to the development of repolarization abnormalities. To investigate how intracellular Na⁺ affected repolarization abnormality occurrence, we clamped its concentration to a physiologically normal value of 7 mM (the baseline concentration in the ORd model), and to an overloaded value of 20 mM, which is within the range reported in heart failure (Pieske et al., 2002) and for acute ischemia. For the 20 mM condition, we also clamped intracellular K⁺ to 145 mM to prevent rundown.

For intracellular Na⁺ clamped at 7 mM, repolarization abnormalities occurred in 56 and 95% of models at 75 and 90% I_Kr block respectively, close to the values for unclamped Na⁺. However, clamping intracellular Na⁺ at 20 mM and K⁺ at 145 mM caused a reduction in repolarization abnormalities to 27 and 46% of models at 75 and 90% I_Kr block respectively (Figure 7C). These results suggested that normal Na⁺ concentration alone was not sufficient to prevent repolarization abnormalities, and that elevated intracellular Na⁺ may make an important contribution to repolarization reserve and decrease occurrence of repolarization abnormalities by increasing the outward current through I_NaK.

To test whether the outward current provided by I_NaK itself was an important contributor to repolarization reserve during I_Kr block, we removed the electrogenic effects of I_NaK on membrane potential, while still allowing ion transport to occur normally. Without the electrogenic contribution of I_NaK, repolarization abnormalities occurred in 93% of low G_NaK HS models following 75% I_Kr block, and in 100% of models following 90% I_Kr block (Figure 7C). This large increase in abnormalities when only the electrogenic component of I_NaK was removed, even in models with a low baseline G_NaK, combined with the large relative magnitude of I_NaK during repolarization shown in Figure 6, provides evidence supporting an important role for I_NaK as a contributor to repolarization reserve in humans during I_Kr block, over a wide range of different conductance profiles.

Hypokalemia and Ca²⁺ loading are both known to increase the likelihood of afterdepolarization generation (Sato et al., 2010; Weiss et al., 2010; Despa and Bers, 2013). We investigated how much each phenomenon altered the incidence of repolarization abnormalities in the population of models under I_Kr block (Figure 8). Mild hypokalemia was simulated by reducing extracellular potassium from 5.4 to 4.0 mM in each model. Each model in the population was simulated under these conditions at five levels of I_Kr block: 0/25/50/75/90%.

Reducing extracellular K⁺ alone was not sufficient to generate repolarization abnormalities, but, combined with I_Kr block it resulted in a small increase in the number of models that displayed abnormalities relative to control K⁺ levels (Figure 8B).

To investigate whether susceptibility was increased by Ca²⁺ loading, we used rapid pacing at 2.5 Hz, which is known to cause Ca²⁺ loading (Maruyama et al., 2010), followed by a reduction in pacing rate to 1 Hz for the final two beats of each simulation, to match the pacing rate to the other simulations in this study. This protocol resulted in 52% of models displaying an average increase in intracellular Ca²⁺ of 50% or more compared to pacing at 1 Hz. The final two beats of each simulation showed no development of afterdepolarizations across the population under control conditions. However, pacing at 2.5 Hz lead to increased repolarization abnormality occurrence following I_Kr block, as shown in Figure 8C. Afterdepolarization occurrence across the population increased from 4 to 14% of models at 75% I_Kr block, and from 18 to 30% at 90% I_Kr block.

In this study, a population of human ventricular cardiomyocytes, calibrated using electrophysiological recordings to yield normal APs under control conditions, provides a quantitative understanding of the relative importance of ionic currents in the occurrence of repolarization abnormalities under ionic drug block. Our main findings include:

I_Kr block was necessary for repolarization abnormalities to occur. No model in the population displayed repolarization abnormalities in any of the simulations in this study where I_Kr was not blocked. This result is supported by the inclusion of simulations with simultaneous 90% block of the repolarizing K⁺ currents I_Ks and I_K1, and the presence of models in our population with very low values of G_Kr (Figure 3C). This finding highlights the particular importance of I_Kr in human ventricular cardiomyocytes, and the unique role of I_Kr block in destabilizing repolarization.

Unexpectedly, we found that I_NaK was an important component of repolarization reserve, specifically, that it made the second largest overall contribution to repolarizing current after I_Kr. While I_Kr inhibition and I_CaL availability were sufficient for repolarization abnormalities to occur in the strongest drug block scenarios (i.e., 90% I_Kr block combined with block of either I_Ks or I_K1), when we analyzed the 41 models in the HS sub-population, we found that G_NaK was less than 50% of its baseline value in 38/41 of these models (Figure 7A). Even in these models, I_NaK was the largest contributor to repolarization current after I_Kr, and of greater magnitude than I_Ks following 75% I_Kr block. A number of disease conditions are associated with reduced expression of the Na⁺/K⁺ pump (Bueno-Orovio et al., 2014), including ischemia (Fuller et al., 2003), heart failure (Xu et al., 1996; Schwinger et al., 1999), and additional comorbidities such as diabetes (Bossuyt et al., 2005). The pump is also selectively blocked by cardiac glycosides, such as digoxin, digitoxin and digitalis (Suhail, 2010). Reports of the effects of reduced G_NaK have focused on the role of the pump in maintaining Na⁺ homeostasis, with elevation of intracellular Na⁺ concentration linked to compromised I_NaK (Wasserstrom and Aistrup, 2005). Our results identified G_NaK as an important determinant of repolarization reserve across a wide range of conductance profiles, and suggest that reduction of repolarization reserve could be a further pro-arrhythmic effect of reduced I_NaK.

As well as being reduced in disease conditions, I_NaK is inhibited by a range of drugs, including multichannel blockers and cardiac glycosides. For example, amiodarone, a commonly used anti-arrhythmic compound that affects multiple channels, also inhibits I_NaK (Gray et al., 1998). Furthermore, there is evidence that cardiac glycosides can reduce trafficking of the hERG channel, which conducts I_Kr in humans, and in particular, that digitoxin could cause a significant reduction in hERG trafficking at therapeutic concentrations (Wang et al., 2007). Another cardiac glycoside, digoxin, has been associated with increased risk of sudden cardiac death (Niemeijer et al., 2015). Therefore, there is a need for further experimental investigation to clarify the role of I_NaK in ventricular repolarization, particularly following I_Kr block. In addition, as the ORd model currently represents our most comprehensive quantitative model of human ventricular cardiomyocyte repolarization, experimental testing of this study's predictions will be important for determining if the ORd model accurately represents the role of I_NaK in repolarization, or if the current model needs to be updated.

Within the population of models a wide variety of conductance profiles produced viable APs (Figure 3A), far from the original conductances of the baseline ORd model. This is consistent with results from previous studies showing the robustness of APs to changes in ionic conductances (Marder and Taylor, 2011; Sarkar and Sobie, 2011; Sánchez et al., 2014; Passini et al., 2016; Zhou et al., 2016). Therefore, all models in the population displayed normal APs with important AP biomarkers within the experimental physiological range. However, across the population there were large differences in models' responses to current block. Even under high levels of I_Kr block, most models did not develop repolarization abnormalities (Figures 5A–C), but almost all models did exhibit substantial APD prolongation (Figure S3B). Therefore, while strong I_Kr block alone is sufficient to cause AP prolongation, it is not sufficient to cause repolarization abnormalities, and also requires an additional mechanism, such as the combination of high I_CaL and low I_NaK. Further, through a logistic regression model, we found that different levels of susceptibility to drug-induced repolarization abnormalities could not be identified prior to current block using standard AP biomarkers and control conditions simulations alone. The logistic regression model could, on average, only correctly assign 61% of the HS models to the HS category. However, combining classification techniques from machine learning and recordings from additional experimental conditions could potentially discover combinations of key features that are more predictive of response to drug block.

In this study, we focused on variability of channel densities, and therefore sarcolemmal conductances, as an important source of variability because of the wide array of biological phenomena known to alter ion channel expression and trafficking, and because conductance alteration through block is the primary mechanism through which many pharmaceutical compounds act on the heart. We have shown that variation in sarcolemmal conductances alone is sufficient to produce drug-induced repolarization abnormalities in models with normal control APs and APDs. However, as ion channel kinetics can vary significantly between individuals, e.g., through genetic mutations (Marban, 2002), future studies could investigate how variation in channel kinetics interacts with conductance variability, e.g., for inherited channelopathies associated with long QT syndrome. For studies investigating phenomena that are known to increase the current through particular channels (e.g., beta-adrenergic stimulation), increasing the range of sampled conductance values for those particular currents could also be useful.

The experimental recordings used in this study, and therefore the sources of AP biomarker variability used to calibrate the population of models, were recorded from endocardial tissue in the right ventricle. Electrophysiological properties vary intramurally and between left and right ventricles (Bueno-Orovio et al., 2012), and these additional sources of electrophysiological variability were not considered in this study. However, we believe the wide ranges of AP biomarker values and conductance profiles considered here mitigate this limitation, particularly relative to using a single model of ventricular electrophysiology.

Although this was not the focus of our investigations, reduced I_NaK has been linked to DAD occurrence (Rosen et al., 1973; Rosen, 1985; Despa and Bers, 2013). DADs were not observed in any model in the population using the stimulation protocols in our simulations, although we did not specifically design simulation protocols to generate them. Therefore, specific simulation protocols and possibly alterations in model structure, particularly in the calcium sub-system, may be required for these investigations.

Even under heavy I_Kr block, only a minority of models in the population displayed abnormalities, and these models occupied a limited parameter regime. It is therefore likely that EAD formation in vivo requires the combined action of multiple destabilizing factors, as in healthy tissue electrotonic coupling suppresses repolarization abnormalities unless they occur simultaneously in a large number of coupled cells (Xie et al., 2010). Electrotonic coupling would also likely prevent the complete repolarization failure that occurred in some models under high I_Kr block from occurring in vivo. However, in diseased or damaged tissue, reduced cellular coupling, combined with increased structural heterogeneity could potentially overcome electrotonic coupling (Shaw and Rudy, 1997; Xie et al., 2010), or even cause electrotonic current to become a generative mechanism for EAD development (Dutta et al., 2016). This tissue-level destabilization could increase the range of ionic profiles that are vulnerable to developing repolarization abnormalities and reduce the level of drug block required to initiate them. Therefore, understanding how tissue heterogeneity interacts with variability in conductances and drug block to influence repolarization abnormality formation is an important topic for future investigations.

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.