Global Optimization of Ventricular Myocyte Model to Multi-Variable Objective Improves Predictions of Drug-Induced Torsades de Pointes

In silico cardiac myocyte models present powerful tools for drug safety testing and for predicting phenotypical consequences of ion channel mutations, but their accuracy is sometimes limited. For example, several models describing human ventricular electrophysiology perform poorly when simulating effects of long QT mutations. Model optimization represents one way of obtaining models with stronger predictive power. Using a recent human ventricular myocyte model, we demonstrate that model optimization to clinical long QT data, in conjunction with physiologically-based bounds on intracellular calcium and sodium concentrations, better constrains model parameters. To determine if the model optimized to congenital long QT data better predicts risk of drug-induced long QT arrhythmogenesis, in particular Torsades de Pointes risk, we tested the optimized model against a database of known arrhythmogenic and non-arrhythmogenic ion channel blockers. When doing so, the optimized model provided an improved risk assessment. In particular, we demonstrate an elimination of false-positive outcomes generated by the baseline model, in which simulations of non-torsadogenic drugs, in particular verapamil, predict action potential prolongation. Our results underscore the importance of currents beyond those directly impacted by a drug block in determining torsadogenic risk. Our study also highlights the need for rich data in cardiac myocyte model optimization and substantiates such optimization as a method to generate models with higher accuracy of predictions of drug-induced cardiotoxicity.

In silico cardiac myocyte models present powerful tools for drug safety testing and for predicting phenotypical consequences of ion channel mutations, but their accuracy is sometimes limited. For example, several models describing human ventricular electrophysiology perform poorly when simulating effects of long QT mutations. Model optimization represents one way of obtaining models with stronger predictive power. Using a recent human ventricular myocyte model, we demonstrate that model optimization to clinical long QT data, in conjunction with physiologically-based bounds on intracellular calcium and sodium concentrations, better constrains model parameters.
To determine if the model optimized to congenital long QT data better predicts risk of drug-induced long QT arrhythmogenesis, in particular Torsades de Pointes risk, we tested the optimized model against a database of known arrhythmogenic and non-arrhythmogenic ion channel blockers. When doing so, the optimized model provided an improved risk assessment. In particular, we demonstrate an elimination of falsepositive outcomes generated by the baseline model, in which simulations of nontorsadogenic drugs, in particular verapamil, predict action potential prolongation. Our results underscore the importance of currents beyond those directly impacted by a drug block in determining torsadogenic risk. Our study also highlights the need for rich data in cardiac myocyte model optimization and substantiates such optimization as a method to generate models with higher accuracy of predictions of drug-induced cardiotoxicity.

INTRODUCTION
Mathematical models of cardiac electrophysiology are at the cusp of usage in a variety of clinical and pre-clinical applications, including safety pharmacology (Mirams et al., 2012;Zhang et al., 2016). In particular, mathematical modeling forms a central component in the Comprehensive in Vitro Proarrhythmia Assay (CiPA) initiative, a proposed strategy for progressing drug safety testing (Sager et al., 2014;Colatsky et al., 2016;Fermini et al., 2016;Gintant et al., 2016).
In terms of cardiotoxicity, drug safety testing aims to avoid Torsades de Pointes (TdP), a life-threatening ventricular tachycardia. Indeed, occurrences of drug-induced TdP in patients have lead to regulatory bans and market withdrawals of several drugs (Mirams et al., 2011). TdP risk is associated with prolongation of the QT interval on the electrocardiogram, in particular due to block of the hERG channel, which carries the rapid delayed rectifier current (I Kr ; Sanguinetti et al., 1995;Straus et al., 2005;Hoffmann and Warner, 2006). However, multiple other currents and dynamics are of importance to torsadogenesis, and including measured effects of drugs on multiple channels, rather than just hERG, into TdP risk stratification improves risk prediction (Kramer et al., 2013;Mistry et al., 2015). Mechanistically, TdP initiation is linked to early afterdepolarizations (EADs) at the cellular level. Triggering of these EADs may depend directly on multiple different ionic currents, including the L-type calcium current (I CaL ) and the late sodium current (I NaL ) (Lankipalli et al., 2005;Hale et al., 2008), and may also depend on intracellular calcium and sodium dynamics (Terentyev et al., 2014;Kim et al., 2015;Xie et al., 2015;Krogh-Madsen and Christini, 2017), implying that the levels of the ionic transporters that control these concentrations (e.g., the sodium-calcium exchanger and the sodium-potassium pump) are important for torsadogenesis. Indeed, recent in silico work have pointed to the magnitudes of these two transporters as having large impact on TdP risk (Lancaster and Sobie, 2016).
Despite the proposed usage of mathematical models in safety pharmacology, even recent and sophisticated models of human ventricular myocyte electrophysiology perform poorly when simulating each of the most typical congenital long QT (LQT) syndromes (Mann et al., 2016). This naturally raises concerns about the ability of these in silico models to predict drug-induced LQT and TdP. However, using a global optimization strategy, in silico models can optimized to reproduce repolarization delays consistent with those seen clinically in the congenital LQT patient datasets, providing optimism for clinically-related model usage (Mann et al., 2016). A concern remains, however, as to whether these in silico models, optimized in terms of action potential properties, replicate dynamics of intracellular ionic concentrations well enough to reliably predict TdP risk. For example, when optimizing a model in terms of its electrical activity only, it can be difficult to correctly identify parameters that mainly control ionic concentrations (Groenendaal et al., 2015). Indeed, previous modeling studies have shown how identical-looking action potentials, modeled using different combinations of model parameters, can have differing calcium transients (Sarkar and Sobie, 2010).
To investigate this possible limitation, we therefore carried out a multi-variable optimization, using both clinical congenital LQT data and constraints on the concentrations on intracellular Ca 2+ and intracellular Na + ([Ca 2+ ] i and [Na + ] i ). We then asked whether optimized models that better represent the congenital LQT syndromes might allow for more accurate and more reliable modeling of acquired LQT and TdP risk. To this end, we simulated 86 cases of multi-channel drug block with known TdP risk level (Lancaster and Sobie, 2016) and found that the model optimized in terms of both action potential and [Ca 2+ ] i , and [Na + ] i data, better predicts TdP risk.

Cell Model and Drug Simulations
Simulations were performed using the O' Hara-Rudy (ORd;O'Hara et al., 2011) human ventricular ionic model as the baseline model, as this is the model proposed to be used in the CiPA initiative Fermini et al., 2016). We used endocardial myocyte parameter settings except where otherwise noted ( Figure 4A). We used a 1 Hz pacing rate and corresponding steady-state initial conditions (O'Hara et al., 2011). For each perturbation to the model (simulating drug block or LQT syndromes and parameter changes during the optimization; detailed below), the model was simulated for 500 beats prior to collecting data. We quantified action potential duration to 90% (APD 90 ) or 50% (APD 50 ) repolarization, as indicated. Calcium transients were characterized by diastolic level (the minimum [Ca 2+ ] i attained within an action potential cycle cycle) and systolic concentration (the peak [Ca 2+ ] i reached during an action potential). The [Na + ] i varies little within a single action potential and was measured as the maximum value.
For our drug simulations, we used the datasets of Kramer et al. (2013) and Mirams et al. (2011) as curated by Lancaster and Sobie (2016) with an associated yes/no risk of torsadogenesis. For each drug, the dataset gives its estimated effective free therapeutic plasma concentration (EFTPC), along with IC 50 values for block of the channels generating I Kr , I CaL , and the fast sodium current (I Na ). Drug effect on each channel type was modeled as a conductance block based on a Hill equation with a coefficient of 1: where G x,drug is the maximal conductance of channel x in the presence of the drug. The dataset contains 86 entries, with some duplicate drugs modeled differently by the two original sources. There are therefore 68 different compounds in the set, covering a variety of intended clinical use, including anti-arrhythmics, anti-histamines, antipsychotics, hypertension/angina drugs, and others.

Drug Classification
To classify model output generated by these drug simulations we used a Support Vector Machine (SVM; Ben-Hur et al., 2008). We used linear decision boundaries separating two categories of data: TdP risk and no TdP risk. These decision boundaries were computed as the solution to a minimization of an error (E) calculated as the sum of squared distances between the location of miscategorized points and the boundary. Because the two variables used in the classification (APD 50 and diastolic [Ca 2+ ] i ) have very different absolute values, we normalized them to baseline (i.e., no drug) values. In general, the minimum value of E will take on different values when using different cell models to simulate drug effects, indicating that the separation of data points by category is better for some models than others. Therefore, to compare goodness of the classification between models and also to determine sensitivity of the decision boundaries, we calculate regions for which E remains below a threshold value (E * ), which we set to twice the value of the lowest value of E found among the four tested models.

Model Optimization
We optimized the baseline ORd model based on clinical data from LQT patients, following a similar strategy as Mann et al. (2016) and using their QT interval data for control patients and patients with one of the three most prevalent congenital LQT syndromes: LQT1, LQT2, or LQT3. The QT interval data from LQT1 and LQT2 patients came from a patient cohort with heterozygote nonsense mutations only, as that can be mimicked in the model by decreasing the conductances of I Ks and I Kr by 50%, respectively. The LQT3 cohort data is more heterogenous and the subtype was simulated by increasing the conductance of I NaL by a factor that was allowed to vary as part of the optimization process. The amount of QT interval prolongation in these patient groups was 12.2% for LQT1, 16.6% for LQT2, and 16.2% for LQT3. Mapping the delayed repolarization measured clinically as QT interval prolongation directly to APD 90 prolongation in the cell model, the objective data set was 267.97 ms (control), 301.14 ms (LQT1), 312.20 ms (LQT2), and 311.55 (LQT3).
In its simplest setup, the optimization was designed to minimize a sum-of-squares error from the APD 90 objective when subjecting the model to control conditions and each of the LQT subtypes 1, 2, and 3. We refer to this as the "APD LQT " optimization. In other optimizations, we included [Ca 2+ ] i and [Na + ] i information in the objective to improve the optimization outcome. This "multi-variable" optimization was done by adding a hefty error (200 ms squared) if [Ca 2+ ] i and [Na + ] i fell outside a certain range during the control condition. We used a range of 0.05-0.15 µM for diastolic [Ca 2+ ] i , 0.3-0.7 µM for systolic [Ca 2+ ] i , and 7-10 mM for [Na + ] i based on measurements in human ventricular myocytes and recent modeling work (Beuckelmann et al., 1992;Piacentino et al., 2003;Grandi et al., 2010).
For the optimization method, we used a genetic algorithm (GA), which is a global optimization method that has been successful in optimizing cardiac ionic models to experimental and simulated data (Syed et al., 2005;Bot et al., 2012;Kaur et al., 2014;Groenendaal et al., 2015). We used a population size of 200 individual model instantiations and ran each GA for 50 generations. All other settings specific to the GA (detailing selection, crossover, mutation, and elitism) were defined as detailed previously (Bot et al., 2012). Because of the stochasticity inherent to the GA, each optimization was run ten times. We used the run resulting in the lowest error as the optimized model.
The parameters to be determined in the optimization process are scaling factors for the currents I Kr , I CaL , I NaL , the slow delayed rectifier current (I Ks ), the sodium-calcium exchange current (I NCX ), the sodium-potassium pump current (I NaK ), and the extent of I NaL increase with simulated LQT3. All scaling parameters were allowed to vary from 0.1% to 10-fold their values in the baseline model.
Note that for ease and consistency we will refer to current scaling factors as scaling of maximal conductances (and use G Kr , G CaL , G NaL , G Ks , G NCX , and G NaK for the currents defined above), although some currents are technically scaled by a permeability or a maximal charge carried.

Sensitivity of APD, [Ca 2+ ] i , and [Na + ] i in Baseline Model
To help guide our optimization procedure, we first did a sensitivity analysis to the major conductances as parameters with low sensitivity are problematic to determine in an optimization.
In the baseline ORd model, the action potential duration of the ORd model is highly sensitive to changes in I Kr ( Figure 1A). For example, when decreasing G Kr by 50% to simulate LQT2, the response is a prolongation of APD 90 by 117 ms (44%), substantially larger than the QT interval prolongation of 68 ms (16.6%) seen in LQT2 patients with heterozygote nonsense mutations (Mann et al., 2016). The APD has an intermediate sensitivity on G CaL , but shows little sensitivity to variations in G Ks and G NaL , the currents associated with LQT1 and LQT3, respectively. For example, reducing G Ks by 50% to mimic LQT1, gives a modest 8-ms (3%) APD 90 prolongation, much shorter than the 51 ms (12%) QT interval prolongation seen clinically (Mann et al., 2016).
As expected, the calcium transient has a very different parameter sensitivity dependence. It depends strongly on the conductances of I CaL and I NCX , with a 50% increase in G CaL or a 50% reduction of G NCX increasing systolic [Ca 2+ ] i by almost 0.2 µM ( Figure 1B). The calcium dynamics also has a significant dependence on G NaK , which only controls [Ca 2+ ] i indirectly via [Na + ] i changes that regulate I NCX . Indeed, [Na + ] i depends sensitively on G NaK , with a 50% reduction in G NaK resulting in a 1.7 mM increase in [Na + ] i ( Figure 1C). Variations in the remaining key conductances have little influence on [Na + ] i levels.
These results are consistent with those presented previously for ±10 and ±20% parameter variations in the ORd model (O'Hara et al., 2011).

Model Optimization
As it is difficult to estimate parameters to which an output is not sensitive, the above analysis suggests that if optimizing the baseline ORd model to APD data only, it will be problematic to estimate many of the conductance parameters. Including repolarization delay data from LQT types 1, 2, and 3 as additional information to the objective may help determine the scaling of G Ks , G Kr , and G NaL . Further, pinpointing these parameters may narrow down other conductances that correlate with these more directly determined parameters (Groenendaal et al., 2015). The sensitivity analysis also indicates that inclusion of calcium transient data to the optimization objective should help determine G NCX and G NaK scaling, and that incorporation of [Na + ] i may further help determination of G NaK scaling.
We therefore optimized the baseline ORd model to both clinical QT interval data from LQT patients and [Ca 2+ ] i and [Na + ] i as detailed in section 2.3. The model optimized to this multi-variable objective produces APD 90 values that are within 3% of their target values (Figure 2). The optimized parameter  scaling factors are given in Table 1. The optimized model has a much increased G Ks , resulting in a larger response to the simulated LQT1 condition, matching the target data (Figure 2).
Optimizing the baseline model to APD values only (i.e., omitting the [Ca 2+ ] i and [Na + ] i constraints) results in slightly better matching of the objective (Figure 2; errors with 2%). Optimized parameter values are very different, with large increases in scaling of G CaL and G NaK in addition to the enhanced G Ks scaling ( Table 1). Despite the diversity in parameter scalings among the baseline and the optimized models, the action potential morphology is quite similar across these differently parameterized models ( Figure 3A). We also include for comparison the action potential generated by Mann et al. in an optimization to clinical LQT data under both baseline and β-adrenergic conditions ("APD LQT±βAdr " optimization, Mann et al., 2016). The main difference among the action potential waveforms is a depolarization of phase 2 of the action potential, the amount of which correlates with the upscaling of G CaL from the baseline model (about 2-4 in the multivariable and APD LQT±βAdr models, and almost 10 in the APD LQT model).
However, calcium transients and [Na + ] i levels are vastly different across models (Figures 3B,C). When including [Ca 2+ ] i and [Na + ] i constraints in the optimization, the optimized calcium transient and [Na + ] i level are very close to those of the baseline model, despite the allowed ranges being relatively large (0.05-0.15 µM for diastolic [Ca 2+ ] i , 0.3-0.7 µM for systolic [Ca 2+ ] i , and 7-10 mM for [Na + ] i ). In our optimizations, when optimizing to APD only, the calcium transient is significantly enhanced. This is consistent with the much boosted G CaL and the decreased G NCX scaling (relative to the multi-variable optimization) both of which favor a larger calcium transient. In addition, the G NaK scaling is much increased when optimizing to APD only, consistent with the lower [Na + ] i . For the APD LQT±βAdr optimized model (Mann et al., 2016), both G NCX and G NaK were much increased relative to baseline (scaling of 2.95 and 9.12, respectively), resulting in a small-amplitude calcium transient and a low [Na + ] i .

TdP Prediction
To test how well the optimized models predict TdP risk, we used a dataset consisting of drugs blocking I Kr , I CaL , and I Na to varying degrees, and their associated risk of torsadogenesis (TdP positive or TdP negative; Lancaster and Sobie, 2016). As demonstrated previously, while many of the drugs in this dataset that carry a TdP risk do prolong APD, some drug simulations predict an increased APD for TdP negative drugs (Figure 4A, Lancaster and Sobie, 2016). In particular, in the baseline model, simulations of three non-torsadogenic drugs results in action potential prolongation of 15-25 ms (one of these is noted by a black dot in Figure 4A). As demonstrated by Lancaster and Sobie, simulations of those three drugs also led to a decreased diastolic [Ca 2+ ] i , which was not seen in the TdP positive drugs. Therefore, including diastolic [Ca 2+ ] i as a second metric (APD 50 being the first) by which to classify the drugs, allows for a correct TdP risk categorization of these three otherwise false positives ( Figure 4A, Lancaster and Sobie, 2016). Indeed, using APD 50 and diastolic [Ca 2+ ] i in combination correctly classifies drugs in the dataset with high specificity and sensitivity (Figure 4A, Lancaster and Sobie, 2016).
For the APD LQT±βAdr optimized model (Mann et al., 2016), qualitatively similar results are observed ( Figure 4B). The particular values of diastolic [Ca 2+ ] i over which TdP positive drugs are separated from TdP negative drugs are shifted, reflecting baseline differences. The same three TdP negative compounds that resulted in APD prolongation in the baseline model, give increased APD 50 in this optimized model as well.
When simulating drug application in the multi-variable optimized model, predictions are improved ( Figure 4C). In particular, none of the TdP negative drugs result in APD 50 prolongation beyond 5 ms, implying that APD 50 prolongation in itself is a strong predictor of torsadogenic risk in this model. In addition, many of the TdP negative compounds result in more substantial reductions in APD 50 and/or diastolic [Ca 2+ ] i compared to the baseline and the APD LQT±βAdr optimized models. There is therefore an increased flexibility to the positioning of the decision boundary separating the TdP positive from the TdP negative drugs (dashed lines in Figure 4C mark off area within which the categorization error remains less than the threshold value, E * ).
For the APD LQT optimized model, the classification is less successful (Figure 4D). The same three TdP negative drugs that resulted in false positive APD prolongation in the baseline and in the APD LQT±βAdr optimized model do so here. Further, simulation of several TdP positive drugs result in decreased diastolic [Ca 2+ ] i without much change in APD and therefore form false negatives in this categorization. The presence of these false positives and negatives pose a challenge to the classification and prevents the categorization error from getting below the threshold value regardless of the location of the decision boundary.
What are the ionic mechanisms underlying the improvement in predictive ability by the multi-var model? The simulated drugs that give the false positive APD prolongations for the baseline and APD-optimized models are piperacillin and verapamil (note that two independent measurements for verapamil are included in the dataset). These drugs block both I CaL and I Kr . We investigated the ionics of verapamil (black dots in Figure 4) in more detail, as verapamil is a well-known example of an I Kr -blocking agent that does not prolong the QT-interval and is not torsadogenic (Redfern et al., 2003).
A simulation of verapamil application in the baseline model is shown in Figure 5A. Drug-induced reductions in outward I Kr and inward I CaL are seen to be of similar amplitude and  balance each other during the first 200 ms of the action potential. When I CaL inactivates at this time, the loss of outward I Kr is largely unopposed, leading to a decreased rate of repolarization and APD prolongation. In the multi-variable optimized model, the non-drug action potential is generated through a near-balance between a much increased I CaL and a larger I NCX providing inward current against the outward currents I Kr and I Ks , with I Ks now being of similar size as I Kr (Figure 5B). When simulating verapamil application in the multi-var optimized model, there is a loss of inward current by the direct effect of I CaL conductance block and because of a reduction of I NCX due to the decreased calcium transient. As both I CaL and I NCX are increased in the multivariable optimized model relative to the baseline model, the loss of inward current with verapamil application is amplified, preventing repolarization delays. Further, the increased I Ks in the multi-var model helps maintain repolarization under verapamil application. Thus, factors beyond the scaling of the directly blocked currents I Kr and I CaL contribute to the drug-induced response.

DISCUSSION
We investigated optimization of conductance parameters in a human ventricular myocyte model to match clinical data from LQT patients using constraints on the concentrations of intracellular calcium and sodium ions. Without these constraints, parameter optimization can lead to models with unphysiological calcium transient and [Na + ] i . To test the hypothesis that the optimization would allow the model to make improved predictions of drug-induced arrhythmogenesis, we investigated the ability of the model to determine TdP risk in a large set of known drugs. We found that using the optimized model improves TdP prediction in two complementary manners. First, simulations of three TdP negative drugs that result in APD prolongation using the baseline model result in no or minimal APD prolongation when using the optimized model. Second, when using both diastolic [Ca 2+ ] i and APD 50 for the modelbased drug classification, the optimized model gives an improved separation between the TdP positive and negative drugs, measured as an increased flexibility in the positioning of the decision boundary. Based on these findings, our main conclusion is that intracellular ionic concentrations are important for safety pharmacology modeling.

In Silico TdP Prediction
Other studies have investigated ionic-model-based TdP prediction using different approaches. It is clear from these studies that a range of strategies can be applied to improve TdP prediction. First, there is improvement to the baseline model, which in itself can involve a number of approaches. One is to optimize a model to clinical LQT data, as done here or previously (Mann et al., 2016). Conceptually, fitting a cellular model to clinical ECG data rather than to experimental cellular-level data may appear counter-intuitive, but it makes sense given that the model is used to predict an organ-level, rather than a cellularlevel, arrhythmia risk. Another model optimization approach is to tune the model to experimental data obtained with ion channel blockers. Such data can be additional to the data used originally to build the baseline model, as in the example of a canine model that upon optimization delivered improved prediction of test drug data . In another article in this Research Topic, Dutta et al. (2017) used drug data presented in the original ORd model paper to reparameterize the ORd model. This optimization was done in conjunction with another model improvement strategy: re-casting the I Kr description as a Markov model with state-and voltage-dependent drug block (Li et al., 2017). An alternative strategy to optimizing a model is to generate a population of models to represent inter-individual and/or inter-cell variability, potentially recapitulating variability in drug response across a heterogeneous population (Lancaster and Sobie, 2016;Britton et al., 2017). Another contribution to this Research Topic demonstrates that such population models predict TdP risk better than using a single baseline model (Passini et al., 2017). Population models may also be used to gain mechanistic insights into arrhythmogenesis. For example, Passini et al. determined that different sub-populations of models had different propensities to repolarization abnormalities, with low conductances for the outward currents I Kr and I NaK and increased levels of I CaL and I NCX making models more prone to repolarization abnormalities, emphasizing that currents other than I Kr are important in this aspect. Second, TdP prediction may be improved by selection of better risk measures. While repolarization delay (APD or QT interval prolongation) has high sensitivity to TdP positive drugs, its specificity is more limited, with some drugs prolonging QT interval, yet carrying only low TdP risk (e.g., amiodarone; Sager et al., 2014). Measures that may be useful in risk stratification include diastolic [Ca 2+ ] i (Lancaster and Sobie, 2016) as employed here. While this measure was selected, in combination with APD 50 , from a range of action potential and calcium transient biomarkers through a machine learning process, there is a mechanistic basis as to why [Ca 2+ ] i levels may be associated with TdP risk, as abnormal intracellular calcium dynamics and spontaneous calcium release is associated with EAD formation, a cellular-level trigger of TdP (Lancaster and Sobie, 2016;Němec et al., 2016). Another risk measure proposed from in silico work is the net charge carried during the action potential by six major ionic currents Li et al., 2017). This measure may also be mechanistically linked to TdP arrhythmogenesis, as it is indicative of robustness against EAD generation under a G Kr -reduction challenge . Use of repolarization abnormality occurrence (i.e., EADs or incomplete repolarization) in simulations as a metric for TdP risk may also present a viable stratification pathway (Passini et al., 2017). Given the direct link to arrhythmogenesis, this seems like a promising risk marker, but a possible limitation lies in its use of highly elevated drug concentrations to trigger the repolarization abnormalities, which may lead to an overestimation of the number of false negatives. Use of this metric rather than APD prolongation improves TdP prediction in a population of models, but not in the baseline ORd model (Passini et al., 2017).
In summary, it is clear that in silico cell models can be improved to better predict TdP risk and that measures beyond APD prolongation are helpful to this end, but it also apparent that significant uncertainties remain as to how to best carry out the modeling and the arrhythmia risk prediction.

Kr/Ks Balance
Our optimization resulted in significant rescaling of many parameters, in particular G Ks , which was increased approximately eight-fold. This is comparable to the scaling of 5.75 found in Mann et al. (2016). In the baseline ORd model, I Ks is relatively small under control conditions, its peak value during an action potential being roughly 10 times smaller than peak I Kr . Significant upscaling of this current is therefore necessary to recapitulate the clinical LQT1 phenotype showing substantial QT interval prolongation with loss of I Ks . Likewise, the Grandi-Bers model, another recent human ventricular myocyte model (Grandi et al., 2010) that has little reliance on I Ks under control conditions, requires sizeable upscaling of G Ks (about 25-fold) to reproduce the LQT1 clinical data (Mann et al., 2016). In contrast, the ten Tusscher-Panfilov human ventricular myocyte model (ten Tusscher and Panfilov, 2006) which has similarly sized I Ks and I Kr , requires increased G Kr (2.65-fold) and decreased G Ks (0.41-fold) to reproduce the clinical LQT dataset (Mann et al., 2016).
These substantial increases in G Ks required for the ORd and the Grandi-Bers models to reproduce the clinical LQT data are at odds with the I Ks ranges recorded experimentally. Factors that may contribute to this disagreement include: (1) Differences between levels of β-adrenergic-dependent kinases and phosphorylation which regulate I Ks and exacerbate LQT1 (Wu et al., 2016); (2) Transmural or other intra-heart heterogeneity with some regions having especially delayed repolarization; (3) Methodological experimental limitations with I Ks rundown and/or damage of the I Ks channel due to enzymatic digestion-however, the recordings that formed the basis of G Ks in the ORd model were done in small tissue preparations using microelectrodes for the express purpose of mitigating these complications (O'Hara et al., 2011).
The I Ks conductance was also increased (by 87%) in the recent optimization of the ORd model with the Markov I Kr formulation to the original O'Hara et al. data . While this approximately doubled I Ks at baseline, I Ks remained much smaller than I Kr (by about five-fold) and would not be expected to be able to reproduce the LQT1 phenotype. Because the particular balance between I Kr and I Ks can be important for action potential stability and EAD generation (Devenyi et al., 2017), one may expect a model with large G Ks to behave differently from a model with smaller G Ks model in terms of arrhythmogenesis. Because the reasons for the discrepancy between the experimental and the clinical I Ks data are not known, it will likely be useful to the field to have both a model with large G Ks , replicating the clinical LQT data, and a model with smaller G Ks , replicating the experimental data.
The mismatch between the experimental and the clinicallybased estimations of G Ks also raises broader questions regarding how to best handle inconsistent data in model development. Our approach here has been to use data of perceived highest relevance to the particular type of predictions made, i.e., to use clinicallybased parameter estimations to predict clinical arrhythmia risk. This approach is in line with the general strategy of using data specific to a particular system (e.g., a cell or a patient) to generate a model specific to that system. However, the best way forward may be to couple rigorously uncertainty in model parameters to uncertainty in model predictions using uncertainty quantification tools (Johnstone et al., 2016).

Effect of Verapamil on Action Potential Duration
The balance between different currents is also important for determining a model's response to simulated drug block. The anti-hypertension and anti-angina drug verapamil blocks I CaL and I Kr , does not prolong the QT interval, and does not prolong APD in recordings from human trabeculae (Redfern et al., 2003;Britton et al., 2017). However, different human in silico models give different responses to simulated verapamil application. The Grandi-Bers and the ten Tusscher-Panfilov models predict action potential shortening in response to verapamil (Mirams et al., 2011(Mirams et al., , 2012. In variations of the ORd model, verapamil almost always prolongs the APD, but the response varies depending on drug concentration, on how block is modeled, and on whether a Markov model is used for I Kr Dutta et al., 2017;Passini et al., 2017).
One hypothesis as to why verapamil does not prolong APD is that its block of I Kr is compensated for by block of I CaL . Using our multi-var optimized model, we show here that in addition to the block of I CaL , a secondary reduction in I NCX (due to the decreased calcium transient) is important in offsetting the I Kr block by verapamil. The size of the I Ks current is also important in determining APD under I Kr block conditions as I Ks provides a repolarization reserve. However, I Ks level in itself is not predictive of APD shortening with verapamil since in the APD LQT optimized model, which has a much increased repolarization reserve in I Ks , verapamil leads to APD prolongation.

Limitations
There are several limitations to our modeling and optimization approach. We allowed large ranges of the scaling (0.1% to 10-fold) of the parameters to be estimated in the optimization. Consequentially, the conductance scalings may be unphysiologically large, with, e.g., G Ks becoming larger than estimated experimentally. However, we are explicitly not attempting to make the best model of a single cell or small tissue, but, rather, a model capable of making clinically relevant predictions. We did not include the clinical data from control and LQT types 1, 2, and 3 patients during β-adrenergic stimulation (Mann et al., 2016) in the optimization objective as preliminary optimizations with this addition resulted in adrenergically stimulated action potentials having unsmooth repolarization profiles, characterized by slow late repolarization. Due to experimental difficulties in determining absolute values of [Ca 2+ ] i and [Na + ] i , we based the allowed ranges of these mainly on modeling work, particularly the ORd and the Grandi-Bers models. While experimental measurements of [Ca 2+ ] i in human ventricular myocytes are consistent with the simulated values (Beuckelmann et al., 1992;Piacentino et al., 2003), reported measurements of [Na + ] i are much higher (∼20 mM), but may be overestimated (Pieske et al., 2002;Grandi et al., 2010). While the inclusion of bounds on [Ca 2+ ] i and [Na + ] i provided additional information to constrain conductance parameters, it is likely that inclusion of more data into the objective would help constrain parameters further. Such data could include more repolarization markers, further calcium transient features, and drug block data.
We modeled the drug application using a simple conductance block, although some drugs are known to block in a statedependent manner (Mirams et al., 2011;Di Veroli et al., 2014;Britton et al., 2017;Dutta et al., 2017). However, use of this simpler approach allowed us to simulate a larger drug data set. We used a single model for the drug simulations. It might be valuable in future work to generate a population of models around the optimized model to potentially improve predictions and to give insights into ionic mechanisms underlying population heterogeneity in drug responses.